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 . 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  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  for the one-dimensional wave equation. The extension to higher dimensions was obtained in  and was analyzed in . 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 . This gives the single scattering (Born) forward map used in conventional imaging in radar [12, 11], seismic inversion  and elsewhere. (2) Obtain a fast converging, iterative inverse scattering method for the acoustic impedance in a medium with known and smooth wave speed . (3) Develop a non-iterative “backprojection” imaging method that is free of multiple scattering artifacts .
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  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 . 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  and with the conventional, reversed time migration method . 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 transform111If 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  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
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  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
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 likeand satisfying , and by
the eigenfunctions, which form an orthonormal basis ofwith the appropriate boundary conditions, then we have
The pulse is band-limited, so the sum is for , where
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
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
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
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 .
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
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 .
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 .
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