Reduced order model approach for imaging with waves

08/03/2021 ∙ by Liliana Borcea, et al. ∙ University of Michigan University of Houston Ecole Polytechnique 0

We introduce a novel, computationally inexpensive approach for imaging with an active array of sensors, which probe an unknown medium with a pulse and measure the resulting waves. The imaging function uses a data driven estimate of the "internal wave" originating from the vicinity of the imaging point and propagating to the sensors through the unknown medium. We explain how this estimate can be obtained using a reduced order model (ROM) for the wave propagation. We analyze the imaging function, connect it to the time reversal process and describe how its resolution depends on the aperture of the array, the bandwidth of the probing pulse and the medium through which the waves propagate. We also show how the internal wave can be used for selective focusing of waves at points in the imaging region. This can be implemented experimentally and can be used for pixel scanning imaging. We assess the performance of the imaging methods with numerical simulations and compare them to the conventional reverse-time migration method and the "backprojection" method introduced recently as an application of the same ROM.



There are no comments yet.


page 17

page 22

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

This paper is concerned with an application of reduced order modeling to imaging reflective structures in a known, non-scattering host medium, from data gathered by an active array of sensors that emit probing pulses and measure the resulting waves.

Model reduction is an important topic in computational science, which traditionally has been concerned with finding a low-dimensional reduced order model (ROM) that approximates the response (observables) of a given dynamical system for a set of inputs [2, 3]. In wave-based imaging the inputs and observables are controlled and measured by the sensors in the array, but the dynamical system is not given, as it is governed by the wave equation with unknown coefficients like the wave speed. Thus, we need a data-driven ROM.

Data driven reduced order modeling is a growing field which combines ideas from traditional model reduction and learning [10]. Much of it is concerned with studying dynamical systems using the Koopman operator theory [17, Chapter 1]. Dynamical system identification based on the Koopman operator has been proposed for instance in [9] and [17, Chapter 13]. However, these approaches are difficult to use in imaging because they assume knowing the full state of the dynamical system, aka the snapshot of the wave, at a finite set of time instances. Since we only know the wave at the sensor locations, which are far from the imaging region, a new way of learning the wave propagation from the data is needed.

To our knowledge, the first sensor array data driven ROM for wave propagation was introduced in [13] for the one-dimensional wave equation. The extension to higher dimensions was obtained in [5] and was analyzed in [7]. The latter study showed that wave propagation can be viewed as a discrete time dynamical system governed by a “propagator operator”, where the time step is the data sampling interval. This propagator operator maps the wave from the states at instants and to the future state at time , for any . The ROM in [5, 7] is an algebraic analogue of the dynamical system. Its evolution is controlled by an propagator matrix, given by the Galerkin projection of the propagator operator on the function space spanned by the first snapshots of the wave, assuming that the array records for the duration . What distinguishes the ROM construction from the many other Galerkin projections in the literature, is that it is obtained only from the measurements of the snapshots at the sensors in the array i.e., it does not require knowing the approximation space.

The ROM introduced in [5, 7] has been used so far to: (1) Approximate the Fréchet derivative of the reflectivity to sensor array data map [8]. This gives the single scattering (Born) forward map used in conventional imaging in radar [12, 11], seismic inversion [4] and elsewhere. (2) Obtain a fast converging, iterative inverse scattering method for the acoustic impedance in a medium with known and smooth wave speed [7]. (3) Develop a non-iterative “backprojection” imaging method that is free of multiple scattering artifacts [14].

We propose yet another application of the ROM: Estimate the “internal” acoustic wave that originates from the vicinity of the imaging point and propagates through the unknown medium to the array of sensors. This idea has been tried before for Schrödinger’s equation in the spectral (frequency) domain

[6, 15], where the solutions are smooth functions that are easier to approximate than the internal waves in this paper. We use the internal waves for two novel imaging methods: The first is a computationally inexpensive approach designed to sense rapid changes of the wave speed in the vicinity of the imaging point. Its imaging function is connected to the point spread function of the time reversal process, and we explain how the aperture of the array, the bandwidth of the probing pulse and the medium through which the waves propagate affect the resolution. The second method can be implemented experimentally. It controls the excitation from the array in order to focus waves at the imaging points, and then uses a matched field approach to image with the resulting backscattered wave in a pixel scanning manner.

The paper is organized as follows: We begin in section 2 with the mathematical formulation of the imaging problem and review briefly from [7] the relevant facts about the ROM, needed in the next sections. The estimation of the internal wave is described in section 3. The first imaging method based on this internal wave is introduced and analyzed in section 4. We also give there a comparison with the backprojection imaging method introduced in [14]. The second, pixel scanning imaging method is described in section 5. We use numerical simulations in section 6 to assess the performance of the imaging methods and to compare them with the backprojection approach [14] and with the conventional, reversed time migration method [4]. We end with a summary in section 7.

2 Formulation of the imaging problem and the ROM

We are interested in imaging reflective structures in a non-scattering and known host medium occupying the bounded domain , using data gathered by an active array of sensors located at , for . We suppose that the aperture of the array is planar in three-dimensions or linear in two-dimensions, and call “range” the spatial coordinate in the direction orthogonal to it. The coordinates in the plane (line in two-dimensions) parallel to the aperture are called “cross-range”.

The sensor probes the medium with a pulse and thus generates the wave , the solution of the wave equation


with quiescent initial condition


We assume henceforth a real valued pulse , that is an even function supported in the short interval

and has non-negative Fourier transform

111If the source emits an arbitrary pulse , we can convolve the echos received at the array with . Mathematically, this is equivalent to having the even pulse , with Fourier transform .


which is negligible in the complement of the set , where is the center (carrier) frequency and is the bandwidth.

The reflective structures are modeled in (2.1) by rough changes (jumps) of the wave speed , with respect to the known and smooth reference wave speed of the host medium. These changes are supported in the imaging domain , which is a subset of lying at large distance from the array. The unknown appears as a coefficient in the self-adjoint, second order elliptic operator


with homogeneous boundary conditions at . The self-adjointness of is convenient for the operator calculus in [7] and the next sections, but we note that (2.1) can be written in the standard wave equation form for the acoustic pressure Since is known and equal to at the sensor locations, the measurements of the pressure, sampled in time at interval , define the array data

Figure 1: Illustration of the setup: An array of sensors (indicated with triangles) lying near the accessible boundary probes a medium with incident waves and measures the backscattered waves. The inaccessible boundary is drawn with the dashed line. The sought after reflectors are supported in the remote subdomain .

The domain may be physical or the truncation of an infinite domain, justified by hyperbolicity and the finite duration of the measurements. In either case, we divide the boundary in two parts: The “accessible” boundary , named so because it lies in the immediate vicinity of the array, and the “inaccessible” boundary . The accessible boundary is useful for the ROM construction because the waves propagate only on one side of the array222If such a boundary does not exist, the medium should be known and homogeneous on the other side of the array, so that the waves there can be removed with some additional processing., as illustrated in Fig. 1. We model it as sound hard, using the homogeneous Neumann boundary condition


where denotes the normal derivative. The inaccessible boundary is modeled as sound soft,


and if it is due to the truncation of an infinite domain, it is sufficiently far away from the sensors to affect the waves over the duration of the measurements.

The imaging problem is to estimate the support of the large and localized variations of the wave speed, from the data (2.5) collected by the array.

2.1 Review of the ROM for wave propagation

Here we review briefly from [7] the relevant facts about the ROM, needed to state the new results.

2.1.1 The dynamical system for wave propagation.

We work with the even in time wave


which satisfies for , due to causality and the initial condition (2.2). During the short duration of the pulse, the wave senses only the vicinity of the sensor location , where the wave speed equals the known . Thus, the second term in (2.8) can be calculated and we can work with the data matrices


Because is even, it is easy to obtain from (2.1) and (2.6)–(2.7) that satisfies


with initial conditions derived in [7, Appendix A]


We define throughout functions of the operator in the standard way, using its spectral decomposition deduced from [18, Theorem 4.12]. Specifically, if we denote by

the eigenvalues, ordered like

and satisfying , and by

the eigenfunctions, which form an orthonormal basis of

with the appropriate boundary conditions, then we have


The pulse is band-limited, so the sum is for , where

Note that the solution of (2.10)–(2.13) is


and that the data matrices (2.9) can be written in symmetric inner product form as follows


Here we used that functions of commute, and denoted by

the inner product. We also introduced the “sensor functions”


and used the assumption (2.3) to define the square root of .

The notation in (2.17) reminds us that is a pulse dependent, blurry version of the Dirac . Indeed, comparing (2.17) with (2.13), we note that is the initial state of the solution of (2.1)–(2.2) and (2.6)–(2.7), when the sensor emits the pulse


which is real valued, even and satisfies Because of causality and the finite wave speed, is supported in a ball centered at , with radius of order and it can be computed using the operator in the host medium


Let us group all the sensor functions in the

dimensional row vector field


and define the “snapshots” as the dimensional row vector fields


These are the states of the discrete time dynamical system governed by the “propagator operator”


Indeed, using a trigonometric identity of the cosine, we get that the states evolve like


starting from


2.1.2 Data driven ROM construction.

The ROM is the algebraic analogue of the dynamical system (2.23),


with propagator matrix and states . It corresponds to the Galerkin projection of (2.23) on the dimensional function space spanned by the first snapshots (2.21). Using linear algebra notation, we write this space as


and note that is a dimensional row vector field and that it is unknown. Nevertheless, the construction in [7, Section 2] shows that it is possible to get the ROM (2.25) from what we know: the initial snapshot in (2.20) and the data matrices


where denotes the transpose and we introduced the notation for the integral of the outer product of dimensional row vector functions.

Key to the ROM construction are the data driven, symmetric, positive definite “mass” matrix with blocks


and the “stiffness” matrix with blocks


Here we used definitions (2.21) and (2.27), the self-adjointness of and a trigonometric identity for the cosine.

The ROM propagator is defined in [7, Section 2.2.1] as follows: Let be the block upper triangular matrix obtained from the block Cholesky factorization of the mass matrix


This matrix can be used to write the Gram-Schmidt orthogonalization of


which defines the orthonormal, causal basis of the approximation space (2.26),


with dimensional row vector field components , , called the “orthonormal snapshots”. Then, we have


where the index denotes the inverse and transpose. The first equality in this equation is used to compute from the data, and the second equality shows that it is a projection of the operator .

The first ROM states satisfy


and we note how the algebraic structure of captures the causal wave propagation: The column blocks of are indexed using the time instants , for , while the row blocks of are indexed according to the range locations reached by the wavefront at these instants. The first column block of , which equals , has all but the first block equal to zero, because the true snapshot is supported near the array. The second column block of , which equals , has an additional nonzero block because the true snapshot reaches some range in the medium, and so on.

3 The internal wave

We now use the data driven ROM reviewed above to estimate an internal wave that originates from the vicinity of an arbitrary point and propagates to the array through the true, unknown medium.

The best estimate that we could hope for would be


with the initial state


given by the approximation


of in the space (2.26), blurred a little by the pulse dependent operator . The same reasoning used for the sensor functions (2.17) applies to (3.2) and shows that it is supported in the ball centered at , with radius. We are interested in evaluating the wave (3.1) at and the sensor locations , for . However, we cannot get exactly , because we do not know the approximation space (2.26). The next proposition shows that we can compute instead


at , for where the approximation (3.3) of is replaced by the “ROM point spread function”


calculated with the orthonormal snapshots in the reference medium with known wave speed


We will see in the next section that the tighter the focus of around , the better the imaging using the internal wave . So when can we expect such a result? The answer lies in how well we can approximate the snapshots in in the reference space calculated for the known ,


If it is true that the approximation error is small, then the Gram-Schmidt orthogonalization, which is a stable procedure, gives


and the ROM point spread function (3.5) is an approximation of (3.2).

We discuss in A two setups where we can analyze explicitly the approximation (3.8): in a layered medium and in a waveguide. The error in these two cases is controlled by the time step , the separation between the sensors and the array aperture size. In more general settings we only have numerical evidence that if and the sensor separation are small enough and the aperture is large enough, then the ROM point spread function is peaked at .

Proposition 3.1

Let be the block upper triangular Cholesky factor of the data driven mass matrix , with block entries given by (2.28). The estimated internal wave (3.4) evaluated at the sensor locations and at the time instants , for , is given by the dimensional row vector field


where is the column block of the identity matrix

Proof: Let us begin with the auxiliary dimensional column vectors


and note that


Here the first equality is by the definition (2.26) of , in the third equality we used that the sum over equals the identity matrix and the last equality is due to the Gram-Schmidt orthogonalization (2.31). Applying the operator to both sides of (3.11) and using definition (2.21) we get


Furthermore, using definition (2.19) and the self-adjointness of , we have


for all . Gathering these results in an dimensional column vector and recalling definition (2.20) and the expression (2.27) of the data matrices, we obtain


where the last equality is by (2.28). Finally, we substitute (3.10) in (3.14) and use the Cholesky factorization of the mass matrix to get the result


4 Imaging with the internal wave

In this section we introduce a novel imaging approach based on the internal wave estimated in Proposition 3.1. The imaging function is strikingly simple: it is the squared norm of this wave evaluated at the sensors


Moreover, is easy to compute and it does not even require the full ROM. It just uses the Cholesky factor of the data driven mass matrix with block entries (2.28), and the orthonormal snapshots (3.6) calculated in the reference medium with known wave speed .

Our goal in this section is to analyze (4.1) and show that it gives an estimate of the location of reflective structures embedded in the host medium. The analysis is based on the continuum time approximation, where the sum over is replaced by an integral over time, and assumes a long enough duration of the measurements. We also suppose, as is typical in applications, that the wave speed is constant near the sensors and therefore the accessible boundary,


We begin in section 4.1 with the connection between the internal wave and the Green’s function of the acoustic wave equation. This is useful for obtaining the main result in section 4.2, where we relate to the time reversal process and we discuss its resolution. We end in section 4.3 with a comparison of and the imaging function of the backprojection approach introduced in [14].

4.1 The Green’s function and its connection to the internal wave

The next lemma connects the internal wave (3.4) with the Green’s function of the acoustic wave equation, satisfying


where is an arbitrary point in and is the Laplace operator with respect to .

Lemma 4.1

Let be any point in the imaging domain , which supports the sought after reflective structures. The internal wave (3.4) evaluated at the sensor locations satisfies