A nonlinear method for imaging with acoustic waves via reduced order model backprojection

TL;DR

This paper presents a nonlinear, data-driven imaging method for acoustic waves that improves resolution and reduces artifacts by using reduced order models and a novel orthogonalization process based on measured data.

Contribution

It introduces a nonlinear imaging approach using reduced order models derived from data, enhancing resolution and artifact suppression over traditional linear methods.

Findings

Improves imaging resolution compared to RTM.

Effectively removes multiple reflection artifacts.

Achieves accurate imaging with imperfect velocity models.

Abstract

We introduce a novel nonlinear imaging method for the acoustic wave equation based on data-driven model order reduction. The objective is to image the discontinuities of the acoustic velocity, a coefficient of the scalar wave equation from the discretely sampled time domain data measured at an array of transducers that can act as both sources and receivers. We treat the wave equation along with transducer functionals as a dynamical system. A reduced order model (ROM) for the propagator of such system can be computed so that it interpolates exactly the measured time domain data. The resulting ROM is an orthogonal projection of the propagator on the subspace of the snapshots of solutions of the acoustic wave equation. While the wavefield snapshots are unknown, the projection ROM can be computed entirely from the measured data, thus we refer to such ROM as data-driven. The image is obtained by backprojecting the ROM. Since the basis functions for the projection subspace are not known, we replace them with the ones computed for a known smooth kinematic velocity model. A crucial step of ROM construction is an implicit orthogonalization of solution snapshots. It is a nonlinear procedure that differentiates our approach from the conventional linear imaging methods (Kirchhoff migration and reverse time migration - RTM). It resolves all dynamical behavior captured by the data, so the error from the imperfect knowledge of the velocity model is purely kinematic. This allows for almost complete removal of multiple reflection artifacts, while simultaneously improving the resolution in the range direction compared to conventional RTM.

Figures

Equations

$$p_{tt}=c(x{)}^2\Delta p+q_t(t)\varphi (x)e(x),-\infty <t<\infty ,x\in \Omega ,$$

$$\text{mbox}suppe\subset {\mathcal{B}}_A,$$

$${\int }_{{\mathbb{R}}^d}f(x)e(x)dx={\int }_{\Omega }f(x)e(x)dx={\int }_{\text{mbox}suppe}f(y)d{\Sigma }_y,$$

$$\underset{t\to -\infty }{\lim }p(x,t)=0.$$

$$\mathbf{n}(x)\cdot \nabla p(x,t)=0,x\in {\mathcal{B}}_A,$$

$$p(x,t)=0,x\in {\mathcal{B}}_I,$$

$$\tilde{q}(s^2)={\int }_0^{\infty }2\cos (ts)q(t)dt.$$

$${\mathcal{M}}_p(\varphi e)={p(x,t)|}_{x\in {\mathcal{B}}_A},-\infty <t<\infty ,$$

$$u(x,t)=p(x,t)+p(x,-t),t>0.$$

$$u_{tt}=c(x{)}^2\Delta u,t>0,x\in \Omega ,$$

$$u(x,0)=[\tilde{q}(-A)\varphi e](x),u_t(x,0)=0,x\in \Omega ,$$

$$A=c(x{)}^2\Delta .$$

$${\left[\mathbf{n}(x)\cdot \nabla p(x,t)\right]|}_{x\in {\mathcal{B}}_A}=0,{p(x,t)|}_{x\in {\mathcal{B}}_I}=0.$$

$${\mathcal{M}}_u(\varphi e)={u(x,t)|}_{x\in {\mathcal{B}}_A}=[{\mathcal{M}}_p(\varphi e)](x,t)+[{\mathcal{M}}_p(\varphi e)](x,-t),t>0,x\in {\mathcal{B}}_A.$$

$$\text{mbox}supp{e}_j\cap \text{mbox}supp{e}_i=\varnothing ,j\ne i.$$

$${\mathbf{D}}_{i,j}(t)={\int }_{\text{mbox}supp{e}_i}\psi (y)[{\mathcal{M}}_u(\varphi e_j)](y,t)d{\Sigma }_y,t>0,$$

$${\mathbf{D}}_{i,j}(t)={\int }_{\Omega }\psi (x)e_i(x)u(x,t)dx,t>0,$$

$${\mathbf{D}}_{i,j}^k={\mathbf{D}}_{i,j}(t_k),i,j=1,\dots ,m,k=0,\dots ,2n-1.$$

$$u_j(x,0)=[\tilde{q}(-A)\varphi e_j](x),\frac{\partial u_j}{\partial t}(x,0)=0,x\in \Omega ,j=1,\dots ,m,$$

$$\mathbf{u}(x,t)=\left[u_1(x,t),\dots ,u_m(x,t)\right],$$

$${\mathbf{u}}_{tt}=A\mathbf{u},t>0,x\in \Omega ,$$

$$\mathbf{u}(x,0)=[\tilde{q}(-A)\varphi \mathbf{e}](x),{\mathbf{u}}_t(x,0)=0,x\in \Omega ,$$

$$\mathbf{e}=[e_1(x),\dots ,e_m(x)].$$

$$\mathbf{u}(\cdot ,t)=\cos \left(t\sqrt{-A}\right)\tilde{q}(-A)\varphi \mathbf{e}.$$

$$[{\mathbf{v}}^{\ast }\mathbf{w}](t)=\left[\begin{array}{cccc}⟨v_1,w_1⟩(t)& ⟨v_1,w_2⟩(t)& \dots & ⟨v_1,w_m⟩(t)\\ ⟨v_2,w_1⟩(t)& ⟨v_2,w_2⟩(t)& \dots & ⟨v_2,w_m⟩(t)\\ ⋮& ⋮& \ddots & ⋮\\ ⟨v_m,w_1⟩(t)& ⟨v_m,w_2⟩(t)& \dots & ⟨v_m,w_m⟩(t)\end{array}\right],$$

$$⟨f,g⟩(t)={\int }_{\Omega }\overline{f(x,t)}g(x,t)dx.$$

$$\mathbf{D}(t)=[(\mathbf{e}\psi {)}^{\ast }\mathbf{u}](t)={\mathbf{e}}^{\ast }\psi \cos \left(t\sqrt{-A}\right)\tilde{q}(-A)\varphi \mathbf{e}.$$

$$\hat{A}=c(x)\Delta c(x),$$

$$\mathbf{D}(t)={\mathbf{e}}^{\ast }\psi {\tilde{q}}^{1/2}(-c\hat{A}{c}^{-1})\cos \left(t\sqrt{-c\hat{A}{c}^{-1}}\right){\tilde{q}}^{1/2}(-c\hat{A}{c}^{-1})\varphi \mathbf{e}.$$

$$\mathbf{D}(t)={\mathbf{e}}^{\ast }\psi c{\tilde{q}}^{1/2}(-\hat{A})\cos \left(t\sqrt{-\hat{A}}\right){\tilde{q}}^{1/2}(-\hat{A})c^{-1}\varphi \mathbf{e}.$$

Full text

1 Introduction

We consider here a problem of estimating a coefficient of a scalar wave equation, an acoustic velocity, in some domain of interest, from the measurements of its solutions, typically referred to as wavefields, on an array of transducers located on the boundary of said domain. Two aspects of this problem are typically considered. First, a quantitative reconstruction of the coefficient is the subject of the so-called full waveform inversion (FWI) [16, 13, 32]. Second, one may seek a qualitative estimate of the acoustic velocity, say its discontinuities. The discontinuities in the acoustic velocity lead to reflections of waves and thus may be referred to as reflectors or scatterers. The problem of estimating the reflectors (scatterers) from the measurements of the reflected (scattered) wavefields is often referred to as the imaging problem [23, 15, 21]. It is the imaging of reflectors with acoustic waves that is of interest to us here.

The dependency of wavefields on the acoustic velocity is highly nonlinear, which makes both the quantitative and qualitative estimation of the velocity very challenging. Most state-of-the-art imaging methods employed in seismic geophysical exploration and medical ultrasound diagnostics are linear in the recorded data. When more than one reflector is present, first reflections (primaries) from them are typically mixed together with multiply-reflected wavefields (multiples). This creates additional events in the data compared to linearized case with primaries only. Linear imaging methods may interpret these events as extra reflectors that are not present in the actual medium. This creates image artifacts known as multiples in geophysics [31, 33] or reverberation (“comet tail”) artifacts in medical ultrasound [25, 35]. Consequently, complex (and often computationally expensive) data pre-processing techniques are applied to filter the multiple reflection events out of the data before feeding it to imaging algorithms.

Many conventional methods for imaging with waves rely on the time-reversibility of the wave equation. In one form or another they back-propagate the data from the transducers using certain assumptions on the kinematics of the unknown medium. Such approaches are commonly known as migration [28, 18, 3, 14] or time-reversal [17, 11] methods. Here we propose a completely different approach. Instead of propagating the data, we use the theory of model order reduction [2, 1, 22] to construct an approximate wave equation propagator. The approximant is derived to satisfy certain data interpolation conditions, as explained later. If the data is discretely sampled in time, then the reduced order model (ROM) interpolating the data is an orthogonal projection of the propagator on the (Krylov) subspace spanned by the wavefield snapshots taken precisely at the data sampling instants. Computing the projection requires an orthogonalization of wavefiled snapshots to be performed. Since the wavefields are not known inside the domain of interest, the orthogonalization is performed implicitly. It is a nonlinear operation that is crucial to our approach, as it allows to suppress the effects of multiple reflections and probe the medium of interest with localized wavefields.

The use of ROMs in inversion and imaging is not an entirely new idea. Resistor network ROMs were used successfully to obtain quantitative reconstructions of conductivity in two-dimensional electrical impedance tomography (EIT) both directly and iteratively [4, 6, 7, 9, 5]. They were also employed in [10] to reconstruct Schrödinger potentials. While the resistor networks interpolate the EIT or Schrödinger data, they are not ROMs of projection type, unlike the ones considered here. Projection-based ROMs were first studied in the context of inversion for the coefficient of a parabolic equation in [20, 8]. Finally, in [19] reduced models interpolating the time domain wave equation data were introduced to obtain quantitative reconstructions of acoustic velocity. Here we consider ROMs of the same type as in [19], and their construction follows a similar procedure. However, the use of ROMs is different. While in [19] the information about the acoustic velocity is extracted directly from the matrix entries of the ROM, here we perform a backprojection operation on the ROM to go from the reduced model space back to the “physical” space of the domain of interest. This method was first introduced by the authors in a conference abstract [27], and we present here a thorough exposition of the approach briefly outlined therein.

The paper is organized as follows. We begin in section 2 with a forward data model and a formulation of inversion and imaging problems. To simplify the subsequent ROM construction, we symmetrize the data model in section 3. Then we define the time-domain interpolatory ROM and discuss some of its properties in section 4. An algorithm for computing the ROM from the data is given in section 5. After the ROM is computed, we introduce its backprojection and the corresponding imaging functional. Both questions are addressed in section 6. An interpretation of the imaging functional that explicitly separates the kinematic and reflective behavior is given in section 7. A thorough numerical study is performed in section 8. In particular, the orthogonalized wavefield snapshots that span the projection subspace are discussed in section 8.2. Approximations of delta functions by the outer products of orthogonalized snapshots are considered in section 8.3. We compare the proposed ROM backprojection imaging method to conventional reverse time migration (RTM) in section 8.4. Stability of ROM computation in the presence of noise and a regularization algorithm are addressed in section 8.5. The numerical study is concluded in section 8.6 with two large scale synthetic examples from geophysical exploration and medical ultrasound diagnostics. Finally, we conclude in section 9 and discuss the three main directions of future research.

2 Problem formulation

Consider a scalar wave equation for the acoustic pressure $p(x,t)$

[TABLE]

where $\Omega$ is the domain of interest in $\mathbb{R}^{d}$, $d=2,3$ . Domain $\Omega$ can be either infinite or finite. Here we assume $\Omega$ is finite with boundary ${\mathcal{B}}=\partial\Omega$ split into the accessible ${\mathcal{B}}_{A}$ and inaccessible ${\mathcal{B}}_{I}={\mathcal{B}}\setminus{\mathcal{B}}_{A}$ parts respectively. We also assume that the acoustic velocity $c(x)$ is piecewise smooth in $\Omega$.

Wave equation (1) is driven by the source term comprised of a wavelet $q(t)$, source weighting function $\phi(x)$ and a distribution $e(x)$. Source distribution $e(x)$ is supported on the accessible boundary

[TABLE]

and is defined so that

[TABLE]

for any smooth test function $f$, where $d\Sigma_{y}$ is the surface element (arc length in the case $d=2$) at point $y\in{\mathcal{B}}$.

Since (1) is defined for all $t\in(-\infty,\infty)$, we enforce a decay condition for large negative times

[TABLE]

We prescribe reflective boundary conditions on the accessible boundary

[TABLE]

where ${\bf n}(x)$ is the unit outward normal at $x\in{\mathcal{B}}_{A}$. We also prescribe zero Dirichlet conditions on the inaccessible boundary

[TABLE]

although zero Neumann conditions can be used instead.

Following [19] we assume that the source wavelet $q$ is an even, sufficiently smooth approximation of $\delta(t)$ with a non-negative Fourier transform

[TABLE]

To simplify the exposition we assume that different sources correspond to varying $\phi(x)$ and $e(x)$, while the wavelet $q(t)$ remains the same, i.e. each source emits the same waveform.

All possible sources and receivers on ${\mathcal{B}}_{A}$ provide the knowledge of the boundary measurement operator

[TABLE]

for all $\phi$ and $e$ satisfying (2).

To simplify the exposition it is beneficial to move the source term to the initial condition using a form of Duhamel’s principle, as explained in [19]. Using (7) and the fact that $q(t)$ is even, we can work instead of the pressure $p$ with its even part

[TABLE]

It satisfies a homogeneous wave equation

[TABLE]

with initial conditions

[TABLE]

where we denote by $A$ the spatial part of the PDE operator in (10):

[TABLE]

The homogeneous boundary conditions (5)–(6) remain unchanged

[TABLE]

The knowledge of boundary measurement operator ${\mathcal{M}}_{p}$ implies the knowledge of another operator

[TABLE]

In practice, we do not have the full knowledge of the boundary measurement operator ${\mathcal{M}}_{u}$. Instead, we take a finite sampling of it in both space and time.

First, consider the spatial measurement sampling. We place $m$ transducers on the accessible boundary that act simultaneously as sources and receivers. The corresponding source-receiver distributions are $e_{j}(x)$, $j=1,2,\ldots,m,$ with disjoint supports

[TABLE]

Since the supports are disjoint, we fix single source and receiver weighting functions $\phi(x)$ and $\psi(x)$ respectively. Their particular choice is discussed in section 3. Then the continuum time data that we can in principle measure with such transducer set-up is an $m\times m$ matrix-valued function of time

[TABLE]

which can also be written using (3) and (14) as

[TABLE]

where $u(x,t)$ satisfies initial conditions (11) with $e=e_{j}(x)$. Hereafter we denote matrices and matrix-valued functions by bold letters. Physically, the entry ${\bf D}_{i,j}(t)$ is the signal recorded at receiver $i$ at time $t$ while the medium is excited with source $j$.

Second, in practice the temporal measurements are also sampled discretely at a certain rate. As explained later, our method relies on the temporal sampling being uniform. Hence, we introduce time instants $t_{k}=\tau k$, $k=0,1,2,\ldots,2n-1$, where $\tau$ denotes the sampling interval and the total number of temporal samples is $2n$.

Once the measurements are sampled in both space and time, we can define the sampled data as

[TABLE]

We are now in position to state the following two problems.

Inversion problem. Given measured data (18) estimate the acoustic velocity $c(x)$ in $\Omega$.

Imaging problem. Given measured data (18) and a smooth kinematic velocity model $c_{o}(x)$ estimate the (location of) discontinuities of the acoustic velocity $c(x)$.

In [19] a method to solve the inversion problem in one dimension and its generalization to higher dimensions was considered. The method relies on the techniques of model order reduction. Here we use similar techniques to address the imaging problem. Unlike the quantitative inversion problem, where no a priori knowledge of the acoustic velocity is available, the imaging can be viewed as a problem of qualitative estimation of the velocity discontinuities on top of a known kinematic background. By referring to $c_{o}(x)$ as the kinematic model or background, we assume that it captures sufficiently well the kinematics of the wave equation (10), i.e. the wave travel times between pairs of points in $\Omega$ are approximately equal for $c(x)$ and $c_{o}(x)$.

3 Symmetrized data model

Our imaging approach relies on the construction of a reduced order model of the sampled data (18). In order to understand the form which such a model has to take, it is beneficial to transform (17) and (18) further.

Denote by $u_{j}(x,t)$ the solution of (10) with initial conditions

[TABLE]

i.e. the wavefield induced by source $j$. In what follows it is convenient to work simultaneously with the wavefields induced by all sources. Thus, we define a row-vector-valued function of time ${\bf u}:\Omega\times(0,\infty)\to\mathbb{R}^{1\times m}$ by

[TABLE]

which solves

[TABLE]

where the action of operators on row-vector-valued functions is understood component-wise.

The initial conditions for (21) are

[TABLE]

with

[TABLE]

Using operator functions we can write the solution of (21) with initial conditions (22) as

[TABLE]

In the derivations below it is convenient to use the following notation for products of row-vector-valued functions. Suppose ${\bf v},{\bf w}:\Omega\times(0,\infty)\to\mathbb{R}^{1\times m}$ then by ${\bf v}^{*}{\bf w}:(0,\infty)\to\mathbb{R}^{m\times m}$ we denote

[TABLE]

where

[TABLE]

Then using (25) we can write (17) in matrix form as

[TABLE]

We notice that in (27) there is a certain asymmetry between the source and receiver terms. Since the sources and receivers are collocated, it is possible to transform (27) to a symmetric form which makes the subsequent derivation of the reduced model much simpler. To achieve this we introduce the symmetrized PDE operator

[TABLE]

a similarity transform of $A$, since $A=c(x)\widehat{A}c(x)^{-1}$. Using the facts that operator functions of $A$ commute and that the Fourier transform $\widetilde{q}$ of the source wavelet is non-negative (see section 2) we can rewrite (27) as

[TABLE]

To transform further we observe that similarity transforms commute with analytic operator functions, thus

[TABLE]

To write (30) in a fully symmetric form we need to make the following assumption on the source and receiver functions, namely

[TABLE]

for some known function $\theta$. Since the source and receiver functions are always multiplied by distributions $e_{j}(x)$ supported on the accessible boundary, assumption above only requires the knowledge of the acoustic velocity $c(x)$ on ${\mathcal{B}}_{A}$ which is not restrictive. Indeed, if one has access to ${\mathcal{B}}_{A}$ to place the transducers, it should be possible to determine $c(x)$ there as well.

Under assumption (31) we obtain the expression for the continuum time data in a fully symmetric form

[TABLE]

where

[TABLE]

is the row-vector-valued transducer function ${\bf b}:\Omega\to\mathbb{R}^{1\times m}$.

Note that $\widehat{A}$ is self-adjoint with respect to the standard $L^{2}$ inner product (26), and so are its operator functions including the one in (32). Since the transducer function ${\bf b}$ enters (32) in a symmetric fashion, the continuum time data is also self-adjoint

[TABLE]

which is known as reciprocity, i.e. the data measured at receiver $i$ induced by source $j$ is the same as data measured at receiver $j$ induced by source $i$.

4 Time-domain interpolatory reduced order model

At the core of our imaging approach is the construction of a reduced order model for the sampled data ${\bf D}^{k}$, $k=0,\ldots,2n-1$. To that end we need to transform the symmetrized data (32) using the fact that the data is sampled uniformly in time. Specifically, since

[TABLE]

we can write

[TABLE]

where $T_{k}$ are Chebyshev polynomials of the first kind of degree $k$ and $P$ is the propagator

[TABLE]

We seek the reduced order model in the form (36), namely

[TABLE]

where $\widetilde{\bf P}\in\mathbb{R}^{mn\times mn}$ is the reduced order propagator and $\widetilde{\bf B}\in\mathbb{R}^{mn\times m}$ is the reduced order transducer matrix.

To compute the reduced order propagator and transducer matrices we solve the following interpolation problem

[TABLE]

Note that since there are $m$ transducers, reduced order matrices $\widetilde{\bf P}$ and $\widetilde{\bf B}$ consist of $m\times m$ blocks. In particular, $\widetilde{\bf P}$ consists of $n\times n$ blocks, while $\widetilde{\bf B}$ is a block column matrix $n$ blocks tall. This allows to use block versions of results in [19] on the existence and construction of the solution to (39). This requires defining the following quantities.

First, following the procedure of section 3 we introduce the row-vector-valued function of symmetrized wavefields $\widehat{\bf u}:\Omega\times(0,\infty)\to\mathbb{R}^{1\times m}$ as solutions of

[TABLE]

with initial conditions

[TABLE]

Second, since the data ${\bf D}^{k}$ is sampled uniformly in time, we also sample the symmetrized wavefields to obtain the symmetrized snapshots

[TABLE]

Obviously, the sampled data ${\bf D}^{k}$ in terms of the symmetrized snapshots is given by

[TABLE]

Third, we introduce the following block Krylov subspace

[TABLE]

where by “block” we refer to the fact that each $\widehat{\bf u}^{k}$, $k=0,1,\ldots,n-1$ is a row-vector-valued function of $x$ and thus the dimension of ${\mathcal{K}}_{n}(P,{\bf b})$ is $mn$. We can consider an orthonormal basis for (44) that we denote by ${\bf V}:\Omega\to\mathbb{R}^{1\times mn}$, which we write as

[TABLE]

The basis ${\bf V}$ is orthonormal in the sense of the following “matrix” product. Let ${\bf Y},{\bf Z}:\Omega\to\mathbb{R}^{1\times mn}$ be

[TABLE]

then

[TABLE]

where each $({\bf y}^{k})^{*}{\bf z}^{l}\in\mathbb{R}^{m\times m}$ is defined by (25) with the distinction that ${\bf y}^{k}(x)$, ${\bf z}^{l}(x)$ do not depend on $t$. The role of time is played by the indices $k$ and $l$ once the temporal sampling (42) is introduced. Then we can write the orthonormality conditions for ${\bf V}$ as

[TABLE]

where ${\bf I}_{mn}$ and ${\bf I}_{m}$ are identity matrices of sizes $mn$ and $m$ respectively.

Now we can address solving interpolation problem (39) to obtain the reduced model $\widetilde{\bf P}$, $\widetilde{\bf B}$. The solution is characterized by the following lemma.

Lemma 1**.**

Reduced order propagator $\widetilde{\bf P}\in\mathbb{R}^{mn\times mn}$ and transducer matrix $\widetilde{\bf B}\in\mathbb{R}^{mn\times m}$ that solve the interpolation problem

[TABLE]

are orthogonal projections

[TABLE]

of the true propagator (37) and transducer function (33) on Krylov subspace (44), where ${\bf V}$ is an orthonormal basis for ${\mathcal{K}}_{n}(P,{\bf b})$ in the sense of (48), and the products in (50) are understood in the sense of (46)–(47), in particular

[TABLE]

The proof of the lemma is in Appendix A. We make the following observations regarding Lemma 1 the implications of which are addressed in detail in the subsequent sections.

Not surprisingly, in order to match the data at time instants $t_{k}=\tau k$, $k=0,\ldots,2n-1,$ we project on the subspace of wavefield snapshots $\widehat{\bf u}^{k}$ sampled at the same instants $t_{k}$. However, the Krylov subspace ${\mathcal{K}}_{n}(P,{\bf b})$ is spanned by only $n$ snapshots, while the data is sampled at $2n$ instants. This question is addressed in the proof in Appendix A.

We observe that $\widetilde{\bf P}$ and $\widetilde{\bf B}$ are not unique as defined by (50). This ambiguity comes from the choice of an orthonormal basis in ${\mathcal{K}}_{n}(P,{\bf b})$. If $\widetilde{\bf P}$ and $\widetilde{\bf B}$ are solutions of (49) then ${\bf Q}^{*}\widetilde{\bf P}{\bf Q}$ and ${\bf Q}^{*}\widetilde{\bf B}$ are also solutions for any unitary ${\bf Q}\in\mathbb{R}^{mn\times mn}$. We make use of this fact in the proof in Appendix A.

However, not all orthonormal basis choices are equally useful if one wants to use the reduced model for imaging (or inversion, see [19]). Inversion and imaging applications require the basis functions ${\bf v}^{k}$ to be localized in $\Omega$ which implies a certain structure of $\widetilde{\bf P}$ and $\widetilde{\bf B}$. The exact structure of $\widetilde{\bf P}$ and $\widetilde{\bf B}$ is established in Lemma 4 in section 5.

While (50) characterizes the solution of (49), it does not provide a practical way of computing $\widetilde{\bf P}$ and $\widetilde{\bf B}$. Indeed, it requires the knowledge of ${\bf V}$ which comes from orthogonalization of symmetrized wavefield snapshots $\widehat{\bf u}^{k}$, $k=0,\ldots,n-1$. Obviously, when solving inversion and imaging problems we do not have the knowledge of $\widehat{\bf u}^{k}(x)$ for $x\in\Omega$. The only measurements we make are those of ${\bf D}^{k}$, $k=0,\ldots,2n-1$. Thus, we must find a way to compute projections (50) from the knowledge of the sampled data ${\bf D}^{k}$ only. Algorithm 2 achieving that is derived in section 5.

5 Reduced order model computation

Let us introduce the following row-vector-valued function ${\bf U}:\Omega\to\mathbb{R}^{1\times mn}$ of symmetrized snapshots

[TABLE]

so that ${\mathcal{K}}_{n}(P,{\bf b})=\mbox{span}({\bf U})$.

If we had the knowledge of ${\bf U}$, then one way to compute an orthonormal basis ${\bf V}$ for ${\mathcal{K}}_{n}(P,{\bf b})$ is to apply block Gram-Schmidt orthogonalization procedure to ${\bf U}$. In matrix notation this can be written as a block QR decomposition

[TABLE]

where ${\bf L}\in\mathbb{R}^{mn\times mn}$ is a block lower triangular matrix

[TABLE]

with blocks ${\bf L}_{k,l}\in\mathbb{R}^{m\times m}$. Note that multiplication of ${\bf V}$ by a matrix from the right follows the standard rules of matrix multiplication.

While we do not have the knowledge of ${\bf U}$, relation (53) implies that

[TABLE]

which is a block Cholesky decomposition of the mass matrix

[TABLE]

The advantage of (56) over (53) is that it does not require the knowledge of wavefield snapshots in $\Omega$, but only the inner products between them.

Recall from (42) that the snapshots are expressed in terms of Chebyshev polynomials as $\widehat{\bf u}^{k}=T_{k}(P){\bf b}$. On the other hand, the products of Chebyshev polynomials satisfy

[TABLE]

and so the inner products between the snapshots are given by

[TABLE]

Combining (42) and (43) the data is ${\bf D}^{k}={\bf b}^{*}T_{k}(P){\bf b}$, hence the blocks of the mass matrix can be expressed just in terms of the sampled data

[TABLE]

where ${\bf M}_{k,l}\in\mathbb{R}^{m\times m}$, $k,l=0,\ldots,n-1$. Thus, in order to compute the block Cholesky factor ${\bf L}$ of ${\bf M}$ only the knowledge of the sampled data is required.

Note that the mass matrix notation for the Gramian ${\bf U}^{*}{\bf U}={\bf M}$ that we use is not a coincidence. Indeed, if we want to approximate (42) on the Krylov subspace (44) with a Galerkin method, ${\bf M}$ is the mass matrix and we can also define the stiffness matrix

[TABLE]

which we can also express in terms of the sampled data.

Let us consider the action of the propagator $P$ on the snapshot $\widehat{\bf u}^{k}$:

[TABLE]

i.e. $P$ takes the average of the symmetrized wavefield propagated by time $\tau$ forward and backward. Now we combine (61) with (42) and (57) to express the blocks of the stiffness matrix in terms of the sampled data

[TABLE]

where ${\bf S}_{k,l}\in\mathbb{R}^{m\times m}$, $k,l=0,\ldots,n-1$.

Formulas (59)–(62) enable us to compute the reduced model $\widetilde{\bf P}$, $\widetilde{\bf B}$ from the knowledge of the sampled data only. Indeed, for the reduced order propagator from (50) and (53) we obtain

[TABLE]

Similarly, for the reduced order transducer matrix

[TABLE]

We summarize the derivations above in the following algorithm.

Algorithm 2** (Reduced order model computation from sampled data).**

Given the sampled data ${\bf D}^{k}\in\mathbb{R}^{m\times m}$, $k=0,\ldots,2n-1$, to compute the reduced order propagator $\widetilde{\bf P}\in\mathbb{R}^{mn\times mn}$ and transducer matrix $\widetilde{\bf B}\in\mathbb{R}^{mn\times m}$ that solve the interpolation problem

[TABLE]

perform the following steps:

Form the mass matrix ${\bf M}\in\mathbb{R}^{mn\times mn}$ consisting of blocks ${\bf M}_{k,l}\in\mathbb{R}^{m\times m}$ for $k,l=0,\ldots,n-1,$ according to

[TABLE] 2. 2.

Form the stiffness matrix ${\bf S}\in\mathbb{R}^{mn\times mn}$ consisting of blocks ${\bf S}_{k,l}\in\mathbb{R}^{m\times m}$ for $k,l=0,\ldots,n-1,$ according to

[TABLE] 3. 3.

Perform block Cholesky decomposition of the mass matrix

[TABLE]

to compute the block lower triangular Cholesky factor ${\bf L}\in\mathbb{R}^{mn\times mn}$. 4. 4.

Compute $\widetilde{\bf P}$ and $\widetilde{\bf B}$ from

[TABLE]

Note that while the regular Cholesky factorization is defined uniquely, there is an ambiguity in defining block Cholesky factorization that comes from the computation of diagonal blocks ${\bf L}_{k,k}$. The particular version of block Cholesky algorithm that we use to compute (68) is given below.

Algorithm 3** (Block Cholesky decomposition).**

*To compute block Cholesky decomposition of a matrix ${\bf M}\in\mathbb{R}^{mn\times mn}$ with $m\times m$ blocks, perform the following steps:

For $k=0,1,\ldots,n-1:$*

[TABLE]

For $l=k+1,\ldots,n-1:$

[TABLE]

Other choices can be used instead of (70) to compute the diagonal blocks ${\bf L}_{k,k}$, as long as they satisfy

[TABLE]

For example, one may use eigendecomposition or regular Cholesky factorization of the right hand side of (72) to obtain ${\bf L}_{k,k}$. No matter the choice, the reduced order propagator $\widetilde{\bf P}$ obtained by Algorithm 2 has the structure established in the following lemma.

Lemma 4**.**

Reduced order propagator $\widetilde{\bf P}$ computed by Algorithm 2 has a block tridiagonal structure

[TABLE]

while the reduced order transducer matrix $\widetilde{\bf B}$ has the structure

[TABLE]

where $\boldsymbol{\alpha}^{k},\boldsymbol{\beta}^{k}\in\mathbb{R}^{m\times m}$, $k=0,\ldots,n-1$.

The proof is given in Appendix B. Here we note that since $\widetilde{\bf P}$ defined by (63) is self-adjoint, so are its diagonal blocks $\boldsymbol{\alpha}^{k}=(\boldsymbol{\alpha}^{k})^{*}$, $k=0,\ldots,n-1$.

6 Backprojection imaging functional

Once the reduced order propagator $\widetilde{\bf P}$ is computed from the sampled data with Algorithm 2, we can use it for inversion and imaging. In order to construct an imaging functional from $\widetilde{\bf P}$, we consider first the Green’s function $g(x,y,t)$ for the symmetrized problem (40)–(41). Then $g(x,y,t)$ satisfies

[TABLE]

with initial conditions

[TABLE]

where we denote by $\delta_{y}$ a delta function concentrated at $y$.

Then similarly to (24) we can formally write the solution to (75)–(76) at time $t=\tau$ as

[TABLE]

in terms of the inner product (26). Alternatively, if we treat $\delta_{x}$ and $\delta_{y}$ as row-vector-valued functions with a single column, then we can use the notation of (25) to write

[TABLE]

If we could probe the propagator $P$ with delta functions, this would allow us to image the Green’s function (78). This is not possible if we only have access to the sampled data ${\bf D}^{k}$ and thus to $\widetilde{\bf P}$. However, we can probe $P$ with other functions. Let us consider

[TABLE]

a function of two variables

[TABLE]

As we will see later, $[{\bf V}{\bf V}^{*}](x,y)$ approximates $\delta(x-y)$ in a certain sense. Thus, we introduce an approximate Green’s function

[TABLE]

by replacing the delta functions in (78) with ${\bf V}{\bf V}^{*}$. Then we can transform

[TABLE]

or explicitly

[TABLE]

where the standard rules for matrix-vector multiplication apply.

While $\widetilde{\bf P}$ can be computed from the sampled data, it is not possible to obtain the orthonormal basis functions ${\bf V}$ from the knowledge of ${\bf D}^{k}$ only. Thus, we need to make another approximation in (82)–(83). Recall that when solving an imaging problem we are given the approximate knowledge of the travel times between pairs of points in $\Omega$ in the form of a smooth kinematic model $c_{o}(x)$. For $c_{o}(x)$ we can compute the symmetrized wavefield snapshots $\widehat{\bf u}_{(o)}^{k}(x)$ in the whole domain $\Omega$, as well as their orthogonalization

[TABLE]

and the projected propagator

[TABLE]

Hereafter we use the subscript $(o)$ to refer to all quantities associated with $c_{o}(x)$ including the Green’s function $g_{(o)}(x,y,\tau)$ and its approximation $\widetilde{g}_{(o)}(x,y,\tau)$ defined in (81).

If $c_{o}(x)$ approximates the kinematics of the problem with unknown desired $c(x)$ sufficiently well, we may expect that the orthogonalized snapshots for the two problems are also close

[TABLE]

Particular reasons for expecting such behavior are discussed later. Then we can make another approximation in (81)–(82) to obtain

[TABLE]

where all terms can be computed either from the knowledge of the sampled data or the kinematic model.

Finally, we have all the components to define the ROM backprojection imaging functional

[TABLE]

that approximates the difference between the Green’s functions for the unknown medium $c(x)$ and the kinematic model $c_{(o)}(x)$, that is

[TABLE]

If in the vicinity of a point $x\in\Omega$ the unknown medium $c(x)$ contains a reflector, then the difference (89) will be large. Otherwise, $c(x)\approx c_{(o)}(x)$ and ${\mathcal{I}}_{BP}(x)\approx 0$.

The use of term “backprojection” should not be confused with the classical method in computed tomography [26]. The (filtered) backprojection formula for inversion of Radon transform in computed tomography is linear in the input data, while the computation of ${\mathcal{I}}_{BP}$ from the data ${\bf D}^{k}$, $k=0,\ldots,2n-1,$ in Algorithm 2 is a highly nonlinear procedure due to block Cholesky factorization (68) and inversion (69).

7 Interpretation in terms of Schrödinger wave equation

Here we present an interpretation of the imaging formula (88) that goes beyond the approximate Green’s function perturbation (89). Note that the imaging procedure in section 6 makes a distinction between the kinematic part of $c(x)$ that is assumed to be well approximated by $c_{o}(x)$, and the reflective part of $c(x)$ that is imaged by ${\mathcal{I}}_{BP}$. Meanwhile, the operator (12) has a single coefficient that does not allow for a clear separation between the kinematic and reflective parts. Thus, instead of the wave equation (10) with operator (12), in a slight abuse of notation we consider the wave equation

[TABLE]

with operator

[TABLE]

where $c(x)$ is the wavespeed and $\sigma(x)$ is the acoustic impedance.

The operator (91) is symmetric with respect to the $1/[\sigma(x)c(x)]$-inner product and nonnegative definite (positive definite for a bounded $\Omega$). Thus, following section 3 we use $u(x,t)=\sqrt{\sigma(x)c(x)}\widehat{u}(x,t)$ to transform (90) into $\widehat{u}_{tt}=\widehat{A}\widehat{u},$ where the symmetrized operator $\widehat{A}$ is given by

[TABLE]

Furthermore, we can transform $\widehat{A}$ to Schrödinger form

[TABLE]

with

[TABLE]

and the Schrödinger potential

[TABLE]

where we omit the dependency on $x$ of $c(x)$, $\sigma(x)$ and $q(x)$ to simplify the notation.

To obtain the separation between the kinematic and reflective parts we assume that the wavespeed is known exactly $c(x)\equiv c_{o}(x)$ and we only image the perturbations of the Schrödinger potential

[TABLE]

where $\sigma_{o}(x)$ is the known background impedance. Note that the knowledge of $c(x)\equiv c_{o}(x)$ is a rather reasonable assumption. Discontinuities of $c(x)$ do not produce significant reflections observable for moderate apertures [29, Chapter 6], so typically the inexactness of $c(x)$ does not introduce significant topological errors in the image, and may only result in a slight image deformation.

Similarly to (75) we introduce the Green’s functions $g(x,y,t)$ and $g_{(o)}(x,y,t)$ for the unknown and reference media respectively. They solve

[TABLE]

where ${\mathcal{A}}_{o}={\mathcal{C}}_{o}-q_{o}c_{o}{\mathcal{I}}$, and {\mathcal{C}}_{o}=\sqrt{c_{o}}\nabla\cdot\Big{(}c_{o}\nabla\sqrt{c_{o}}\Big{)}.

According to (89) the functional ${\mathcal{I}}_{BP}(x)$ is related to the perturbation

[TABLE]

It follows from (97)–(98) that it satisfies approximately

[TABLE]

if we assume

[TABLE]

If we also introduce the Green’s function $G_{(o)}(x,y,t)$ solving

[TABLE]

we can express the solution of (100) as

[TABLE]

and so

[TABLE]

Thus, due to causality of (98) and (102), ${\mathcal{I}}_{BP}(x)$ yields $\delta q$ averaged in the travel time ball of radius $\tau/2$ centered at $x$. Note that since (104) only relies on the solutions of (100) for $t<\tau$, the assumption (101) is valid for $\tau$ small enough compared to the scale of inhomogeneities of $q$ (and thus $\sigma$).

8 Numerical study and discussion

The proposed imaging method differs significantly from existing approaches. In this section we illustrate numerically and discuss its various aspects that distinguish it from conventional methods.

8.1 Numerical model

To study the performance of the proposed imaging method we consider the following numerical model. Let $\Omega=[0,3]\;km\times[0,3]\;km$. The inaccessible boundary ${\mathcal{B}}_{I}$ consists of the left, right and bottom sides of the square, while ${\mathcal{B}}_{A}$ is the top side. The accessible boundary ${\mathcal{B}}_{A}$ supports $m=32$ point-like transducers, i.e.

[TABLE]

where $x_{j}\in{\mathcal{B}}_{A}$, $j=1,\ldots,m,$ are the transducer locations.

We use an acoustic velocity model $c(x)$ consisting of a smooth background and two extended reflectors one of which is branching. The smooth part of $c(x)$ is used as the kinematic model $c_{o}(x)$. Both $c(x)$ and $c_{o}(x)$ are shown in Figure 1. Reflectors are relatively high contrast with velocity of $1\;km/s$ inside the reflecting inclusion, while the surrounding background velocities are about $2-3\;km/s$ depending on the location.

The wave equation is discretized with second order finite differences on a uniform tensor product grid with $300\times 300$ nodes spaced $10\;m$ apart. Each of the two reflectors has thickness of $20\;m$, two grid nodes.

The source wavelet is chosen so that its Fourier transform (7) is

[TABLE]

where $\sigma$ determines the characteristic duration of the source wavelet $q(t)$ in time. As discussed in [19], it is essential for the stability of ROM computation for the temporal sampling interval $\tau$ to be consistent with the Nyquist sampling limit for the source wavelet. Since the Gaussian (106) does not have a sharp frequency cutoff, its Nyquist limit is not precisely defined, but the choice $\tau\approx\sigma$ works well in practice. In particular, in all numerical examples in this section we take $\tau=\sqrt{3}\sigma/2\approx 0.87\sigma=1.5\cdot 10^{-2}\;s$. While it is mostly pointless to talk about the “frequency” of the broadband (infinite-bandwidth) source wavelet (106), we can characterize the source in terms of the sampling frequency $\omega_{\tau}=1/\tau$. Here $\omega_{\tau}=66.66\;Hz$. The data is measured at $2n=130$ time instants $t_{k}=k\tau$, $k=0,\ldots,2n-1$ with a terminal time of $1.715\;s$.

8.2 Orthogonalized wavefield snapshots

At the heart of our imaging approach is the implicit orthogonalization of the symmetrized wavefield snapshots $\widehat{\bf u}^{k}$, $k=0,1,\ldots,n-1$. This is an essentially nonlinear operation (due to the nonlinearity of block Cholesky factorization in (68) and Cholesky factor inversion in (69)) that makes possible accounting for such nonlinear phenomena as multiple reflections, as discussed later. Wavefield snapshots for the true medium $c(x)$ and the kinematic model $c_{o}(x)$ as well as their orthogonalized counterparts are shown in Figure 2. The following features can be observed.

First, when comparing orthogonalized snapshots ${\bf v}^{k}$ to $\widehat{\bf u}^{k}$ (and also ${\bf v}^{k}_{(o)}$ to $\widehat{\bf u}^{k}_{(o)}$) we observe the focusing in both the range and cross-range. Hereafter we refer to the direction orthogonal to the transducer array as the range and to the direction along the array as the range. The orthogonalized snapshots have a well pronounced peak below the corresponding source (up to the lateral bending of the wavefront due to $c_{o}(x)$ variations). This peak becomes wider as the propagation time increases, i.e. we observe cross-range defocusing which is the consequence of the decrease of effective aperture as the wave travels away from the transducer array.

Second, a phenomenon closely related to focusing that can only be observed in the reflecting medium $c(x)$ is the suppression of reflections in the orthogonalized snapshots ${\bf v}^{k}$. While the energy is spread out over reflections in $\widehat{\bf u}^{k}$, after orthogonalization the wavefield is concentrated around the peak of ${\bf v}^{k}$. This is an essential feature of the approach, since it allows to suppress the artifacts caused by multiple reflections that conventional linear imaging methods are unable to handle. No matter how complex the reflections in the true medium $c(x)$ are, the propagator $P$ is probed in (63) with focused orthonormal basis functions ${\bf v}^{k}$ that are localized around the same locations as ${\bf v}_{(o)}^{k}$. Note that the presence of reflections can actually improve the focusing, as we can see comparing ${\bf v}^{k}$ to ${\bf v}^{k}_{(o)}$ in Figure 2.

8.3 Approximations of delta functions

The imaging functional (88) can be written as

[TABLE]

which according to (77) and (89) should approximate

[TABLE]

Thus, in order for the imaging method to work, the following approximations of delta functions must hold

[TABLE]

Given the focused nature of orthogonalized snapshots ${\bf v}^{k}$ observed in section 8.2, we expect that their “outer products” in (109) will be well localized, as expected from good approximations of delta functions. This conjecture is explored in Figure 3, where we plot $[{\bf V}_{(o)}{\bf V}^{*}_{(o)}](x_{\alpha},y)$ and $[{\bf V}_{(o)}{\bf V}^{*}](x_{\alpha},y)$ for a few fixed points $x_{\alpha}\in\Omega$ as functions of $y\in\Omega$. We observe indeed that both $[{\bf V}_{(o)}{\bf V}^{*}_{(o)}](x_{\alpha},y)$ and $[{\bf V}_{(o)}{\bf V}^{*}](x_{\alpha},y)$ are well localized and have a clearly pronounced peak for $y\approx x_{\alpha}$ while being close to zero elsewhere. Moreover, we note that the focusing of the outer products ${\bf V}_{(o)}{\bf V}^{*}_{(o)}$ and ${\bf V}_{(o)}{\bf V}^{*}$ is much tighter than that of the orthogonalized snapshots ${\bf v}^{k}$.

The outer products can be seen as nonlinear analogues of point spread functions quantifying the resolution of the imaging functional ${\mathcal{I}}_{BP}$. The nonlinearity comes in as dependency of ${\bf V}$ in ${\bf V}_{(o)}{\bf V}^{*}$ on the true medium $c(x)$. Note though that the dependency is kinematic in nature, that is if $c_{(o)}(x)$ captures the kinematics of the true medium $c(x)$, we expect

[TABLE]

as observed in Figure 3.

8.4 Images and comparison to reverse time migration

Once the focusing and localization properties of orthogonalized snapshots ${\bf v}^{k}$ and their outer product are established, we can finally assess the performance of the backprojection imaging functional ${\mathcal{I}}_{BP}$ on the model described in section 8.1 (see Figure 1). We compare our approach with reverse time migration (RTM) [3], a linear imaging approach that is the standard method for high quality imaging with waves in geophysics [14, 30, 34]. Here we use a particular version of RTM known as pre-stack migration, which is more accurate than the computationally cheaper post-stack migration. The pre-stack RTM imaging functional is hereafter referred to as ${\mathcal{I}}_{RTM}$.

Note that both ${\mathcal{I}}_{BP}$ and ${\mathcal{I}}_{RTM}$ imaging functionals produce weaker images of reflectors that are further away from the array. This is particularly pronounced when the wave has to pass through reflectors closer to the array before reaching the far ones, thus losing energy in the process, as it gets reflected back to the array at each interface. To account for this behavior, in what follows we modify both imaging functionals by scaling them with a multiplier that depends linearly on the distance from the imaging point to the transducer array. In a slight abuse of notation we still refer to such modified functionals as ${\mathcal{I}}_{BP}$ and ${\mathcal{I}}_{RTM}$.

In Figure 4 we compare the images obtained with ROM backprojection formula (88) and with conventional pre-stack RTM for the model described in section 8.1 with two kinematic models: a kinematic model $c_{(o)}(x)$ from Figure 1 consisting of the true smooth background without the reflectors; and a constant kinematic model with $c_{(o)}\equiv 2.5\;km/s$. We observe that in both cases our approach is clearly superior compared to its conventional linear counterpart. The main advantage of backprojection imaging is the automatic removal of multiple reflection artifacts. These artifacts appear inevitably in any linear imaging method due to the interaction of reflectors with reflective boundaries (here mainly with the top part ${\mathcal{B}}_{A}$ of the boundary) and with each other. Each such interaction produces additional events in the data that linear migration methods image as extra reflectors that are not actually present in the true medium. We observe in Figure 4 that the RTM image has many such artifacts with magnitudes comparable to the actual reflectors. Moreover, the artifacts change their location and magnitudes depending on the kinematic model used. Meanwhile, ROM backprojection image is entirely free of such artifacts with the only features present being the actual reflectors.

We also observe in Figure 4 that reflector images produced by ${\mathcal{I}}_{BP}$ are thinner than those in RTM images. We conjecture that ROM backprojection has higher range resolution than RTM, but the detailed resolution analysis of ROM backprojection is outside the scope of this paper and remains a topic of future research.

8.5 Noisy data and regularization

As mentioned above, construction of the reduced model $\widetilde{\bf P}$ with Algorithm 2 is a procedure nonlinear in the data ${\bf D}^{k}$. While such approach allows for automatic removal of multiple reflection artifacts and improved resolution, as demonstrated in section 8.4, it also has downsides compared to conventional linear migration.

One particular issue that we address here is that of stability. The first nonlinear operation in Algorithm 2 is block Cholesky factorization (68) of the mass matrix ${\bf M}$. Unfortunately, the mass matrix is ill-conditioned. In Figure 5 we display the condition number of ${\bf M}$ for the numerical example from section 8.1 as the half-number of data time samples $n$ increases from $5$ to $65$, the value used to compute the images in section 8.4. The condition number increases rapidly from $10^{4}$ to $10^{10}$ for $n=15$ and the subsequent increase is rather slow.

Since $\mbox{cond}({\bf M})$ is so large, the smallest eigenvalues of the mass matrix are very close to zero. Note that ${\bf M}$ depends linearly on ${\bf D}^{k}$ in (66), and even a small amount of noise in the data can push the eigenvalues of the perturbed mass matrix into the negative values. Thus, a regularization scheme is required for the imaging method to work with noisy data. Here we propose a simple regularization approach that is an analogue of Tikhonov regularization for linear inverse problems. The idea is to push the eigenvalues of the mass matrix ${\bf M}_{(\epsilon)}$ computed from the noisy data ${\bf D}_{(\epsilon)}^{k}$ back into the positive values. However, this has to be done consistently with the algebraic structure of the problem. This is achieved by the following simple algorithm.

Algorithm 5** (Data regularization).**

Given noisy data ${\bf D}_{(\epsilon)}^{k}$, $k=0,\ldots,2n-1$, and an initial regularization parameter $\mu>1$, to compute the regularized data ${\bf D}_{\mu}^{k}$, $k=0,\ldots,2n-1$, perform the following steps:

Compute the regularized data

[TABLE] 2. 2.

Assemble the regularized mass matrix

[TABLE]

and find its smallest eigenvalue $\lambda_{min}({\bf M}_{(\mu)})$. 3. 3.

If $\lambda_{min}({\bf M}_{(\mu)})<0$, increase $\mu$ and go to step 1, otherwise quit.

We observe that the regularized data ${\bf D}_{(\mu)}^{k}$ differs from the input noisy data ${\bf D}^{k}_{(\epsilon)}$ at the initial time instant $k=0$ only. Note that ${\bf D}^{0}={\bf b}^{*}{\bf b}\in\mathbb{R}^{m\times m}$ is a positive definite matrix, moreover, for the sources (105)–(106) it is diagonally dominant. Since ${\bf D}^{0}$ only enters the diagonal blocks of the mass matrix according to (112), for sufficiently large $\mu$ the regularized mass matrix ${\bf M}_{(\mu)}$ also becomes positive definite. Obviously, the size of regularization parameter $\mu$ needed to achieve $\lambda_{min}({\bf M}_{(\mu)})>0$ depends on the amount of noise present in ${\bf D}_{(\epsilon)}^{k}$.

It is important that Algorithm 5 modifies the data itself instead of just the mass matrix. When the regularized data is fed to ROM computation Algorithm 2 it will not only affect the mass matrix (66), but the stiffness matrix (67) as well. This makes the regularization scheme more consistent with the algebraic structure of the problem.

We illustrate the performance of Algorithm 5 in Figure 6, where we compare ${\mathcal{I}}_{BP}$ computed from noiseless data to ${\mathcal{I}}_{BP}^{(\mu)}$ computed from the noisy data that is regularized with Algorithm 5. The noisy data ${\bf D}^{k}_{(\epsilon)}$ is obtained from the simulated noiseless data ${\bf D}^{k}$ using a multiplicative noise model from [8] with a noise level parameter $\epsilon$. Such model ensures that the signal-to-noise ratio is approximately $100\%/\epsilon$. Here we take $\epsilon=10\%$ which results in the value of regularization parameter $\mu=1.55$, as the minimum value for which $\lambda_{min}({\bf M}_{(\mu)})>0$ robustly for many realizations of the noise.

We observe in Figure 6 a certain deterioration of the image obtained from the regularized noisy data. First, the reflectors become somewhat thicker which indicates a slight loss of resolution. Second, ${\mathcal{I}}_{BP}^{(\mu)}$ loses some of the contrast at the slanted parts of reflectors thus emphasizing the parts of reflectors that are nearly parallel to the transducer array. However, ROM backprojection still suppresses the multiple reflections very well. While we notice one weak artifact appearing above the bottom reflector, the image is still almost entirely artifact-free compared to RTM images in Figure 4.

8.6 Large scale numerical examples

We conclude the numerical study with two large-scale examples inspired by the applications of imaging with acoustic waves in geophysics and medical ultrasound tomography.

The first numerical example is the famous Marmousi model [12] that has become a standard synthetic test model in geophysics used to evaluate the performance of inversion and imaging algorithms. The model is discretized on a $10\;m$ uniform grid with $N=1360\times 280=$ $380,800$ nodes. The temporal sampling interval $\tau=1.8\cdot 10^{-2}\;s$ with the corresponding sampling frequency $\omega_{\tau}=55.55\;Hz$. The data is measured for $2n=130$ time instants at the array with $m=102$ transducers spaced uniformly every $110\;m$ over the accessible boundary ${\mathcal{B}}_{A}=\{(x,0)\;|\;1.2\;km\leq x\leq 12.3\;km\}$. The kinematic model $c_{o}(x)$ is obtained by convolving the true velocity $c(x)$ with a Gaussian kernel of width $400\;m$ and height of $280\;m$, see Figure 7.

A large scale example like Marmousi allows us to illustrate a way to make ROM backprojection imaging more computationally tractable for large problems. Instead of computing the ROM from the whole data set ${\bf D}^{k}$, $k=0,\ldots,2n-1,$ at once, we split the transducer array $\{1,\ldots,m\}$ into $s$ (possibly overlapping) sub-arrays $S_{i}=\{m_{i}^{min},\ldots,m_{i}^{max}\}$, $i=1,\ldots,s$. For each sub-array $i$ we compute the imaging functional ${\mathcal{I}}_{BP,i}$ from the restriction of the data to that sub-array, i.e. from ${\bf D}^{k}_{S_{i},S_{i}}$, $k=0,\ldots,2n-1$. Then the composite backprojection image is simply a weighted sum of individual images

[TABLE]

with some positive weights $\eta_{i}$, $i=1,\ldots,s$. Note that such approach essentially discards the large offset data, since ${\bf D}^{k}_{S_{i},S_{i}}$ is a block of the full data matrix ${\bf D}^{k}$ containing the entries at most $m_{i}^{max}-m_{i}^{min}$ positions away from the main diagonal. Nevertheless, if the sub-arrays are sufficiently large, the quality of such a composite image is still high, while the computational cost is substantially lower.

In Figure 7 we display two composite images for the Marmousi model. The first is computed for $s=5$ overlapping sub-arrays of $m_{i}^{max}-m_{i}^{min}+1=34$ transducers each, while the second one is for $s=11$ overlapping sub-arrays of $m_{i}^{max}-m_{i}^{min}+1=17$ transducers each. We observe no visible deterioration of the image when using smaller sub-arrays (which discards more far offset data), in fact one may argue that the image for $s=11$ is somewhat sharper than that for $s=5$.

The second large-scale numerical example is a model for medical ultrasound diagnostics of breast cancer. The model consists of a circular phantom of the acoustic velocity distribution $c(x)$ in a cross-section of a human breast, see Figure 8. The phantom was designed and provided to the authors by the Computational Bioimaging Laboratory at the Washington University in St. Louis, Department of Biomedical Engineering, in collaboration with TomoWave Laboratories, Inc. (Houston, TX).

The overall diameter of the phantom is $17\;cm$. A homogeneous background corresponds to healthy breast tissue. Adjacent to the reflective boundary (skin surface) is a fatty layer with rough inner surface. Embedded into the background are six inclusions that represent both benign and malignant lesions. While four round-shaped inclusions represent benign lesions, two star-shaped inclusions (top left and bottom-most, see Figure 8) correspond to malignant tumors. The shape, dimensions and acoustic velocity of inclusions are chosen according to human physiology.

The model is discretized on a $4.34\;mm$ uniform grid with $N=392\times 392=153,664$ nodes and a $784\times 784$ fine grid ($N=614,656$). The temporal sampling interval $\tau=8.6\cdot 10^{-7}\;s$ with the corresponding sampling frequency $\omega_{\tau}=1.16\;MHz$. The data is measured for $2n=80$ time instants at the array with $m=192$ transducers spaced uniformly over the whole boundary. A composite image is formed with the transducer array split into $s=16$ overlapping sub-arrays of $m_{i}^{max}-m_{i}^{min}+1=24$ transducers each.

We observe in Figure 8 that the backprojection image gives a reasonably good characterization of the shape of inclusions. The star-shaped (malignant) lesions are imaged to have much rougher boundaries than their round-shaped (benign) counterparts. This feature of makes ROM backprojection imaging promising for medical ultrasound diagnostics. By comparison, the shapes of the inclusions in the RTM image are substantially more blurred, which makes the distinction between star-shaped and round-shaped lesions difficult. Partially, this is a consequence of lower resolution of RTM compared to ROM backprojection for the same source frequency. Note that the medical ultrasound imaging problem (for soft tissues) exhibits substantially less nonlinearity than geophysical imaging problem. While the contrast in acoustic velocity in geophysical applications can be as high as $300\%$, a typical velocity contrast in human soft tissues is within $5\%$.

9 Conclusions and future work

We presented a novel method for imaging the discontinuities of the acoustic velocity $c(x)$, a coefficient of the wave equation (1), from discretely sampled time domain data ${\bf D}^{k}$ measured at an array of transducers. Unlike conventional migration methods, including time reversal approaches, the proposed ROM backprojection imaging method is nonlinear in the data. The nonlinearity comes from an implicit orthogonalization of time-domain wavefield snapshots that form a basis for the subspace on which the propagator is projected. This allows the method to handle the nonlinear interactions between reflectors that often lead to artifacts in linear migration images. In particular, ROM backprojection images are almost entirely free of multiple reflection artifacts that require special handling in conventional time reversal approaches. While ROM computation can become unstable in the presence of noise in the data, we propose a computationally inexpensive regularization procedure similar to Tikhonov regularization for linear inverse problems.

Although the numerical study in section 8 demonstrate the potential of ROM backprojection imaging approach and its advantages over conventional linear migration methods, there is a number of questions that must be investigated further for ROM backprojection imaging to become a viable alternative to migration in real world applications. We identify the following directions for future research:

• •

Resolution analysis. The numerical results in section 8.4 imply that ROM backprojection imaging may achieve higher resolution in the range direction compared to conventional RTM. This question requires a careful separate consideration. One may consider first the one-dimensional case from [19], where it might be possible to explicitly construct the orthogonalized snapshots ${\bf V}_{(o)}$ and their outer products ${\bf V}_{(o)}{\bf V}_{(o)}^{*}$ which characterize the resolution.

• •

Advanced regularization techniques. While the simple regularization procedure in Algorithm 5 provides reasonably good images from noisy data, more elaborate approaches can be studied. If one views the approach of Algorithm 5 as an analogue of Tikhonov regularization for linear inverse problems, then it is natural to also consider an analogue of truncated SVD regularization [24]. In the context of ROM computation this means projecting both mass ${\bf M}$ and stiffness ${\bf S}$ matrices on the appropriate subspace of singular vectors of the mass matrix. The subspace should be chosen so that the projected mass matrix is positive definite. Note that since the backprojection imaging formula (88) involves the difference between the projected propagators for the unknown medium and the kinematic model, there must be consistency in the choice of singular vector subspaces for the unknown medium and the kinematic model.

• •

Non-collocated sources and receivers. We assumed that the data is measured at an array of transducers that can act both as sources and receivers, thus allowing us to write the data model ${\bf D}(t)$ in the symmetrized form (32). Consequently, this lead to both left and right projection subspaces in interpolation relation (49) to be the same block Krylov subspace (44), which simplifies the method. While an assumption of collocated sources and receivers is not unreasonable for medical ultrasound imaging, it is often not valid in geophysical exploration, where the sources (shots) and receivers (geo-/hydro-phones) are almost always distinct, with the number of shots being significantly smaller than the number of receivers. This may require considering different left and right projection subspaces, which leads to implicit bi-orthogonalization of wavefield snapshots instead of implicit orthogonalization (48) used currently.

Acknowledgments

This material is based upon research supported in part by the U.S. Office of Naval Research under award number N00014-17-1-2057 to L. Borcea and A.V. Mamonov. The work of A.V. Mamonov was also partially supported by the National Science Foundation grant DMS-1619821 and by the New Faculty Research Program (2015-2016) of the University of Houston.

Appendix A Proof of Lemma 1

As mentioned in section 4, the choice of an orthonormal basis ${\bf V}$ for ${\mathcal{K}}_{n}(P,{\bf b})$ does not affect the interpolation condition (49). Thus, we can work with any basis for ${\mathcal{K}}_{n}(P,{\bf b})$. In particular, for the purposes of this proof it is convenient to take

[TABLE]

which is obviously orthonormal in the sense of (47)–(48). In what follows we use extensively

[TABLE]

where ${\bf E}^{k}\in\mathbb{R}^{mn\times m}$ has ${\bf I}_{m}$ block at the $k^{\mbox{th}}$ block position and zeros elsewhere.

We prove first that

[TABLE]

The proof is by induction with two base cases being $k=0,1$.

If we notice that

[TABLE]

then the base case $k=0$ is established by

[TABLE]

For $k=1$ we have

[TABLE]

Now we apply the three term recurrence relation for Chebyshev polynomials

[TABLE]

to prove the general case using (61) and the induction hypothesis

[TABLE]

This establishes (116) which immediately implies the interpolation condition

[TABLE]

for $k=0,\ldots,n-1$.

To prove the interpolation condition (49) for $k=n,\ldots,2n-1$ we consider the following. From the product formula for Chebyshev polynomials (57) it follows that

[TABLE]

thus

[TABLE]

and

[TABLE]

Now suppose that $1\leq l\leq n-1$, then by induction proof above, (116) holds for $\widehat{\bf u}^{l}$ and also $|n-1-l|\leq n-1$, so (122) holds for ${\bf D}^{|n-1-l|}$. Substituting both into (125) we find

[TABLE]

by (123), which proves the interpolation condition (49) for $k=n-1+l=n,\ldots,2n-2$. So we have established (49) for all $k$ except $k=2n-1$.

The case $k=2n-1$ follows from (125) and (123) because

[TABLE]

Appendix B Proof of Lemma 4

Since $\widetilde{\bf P}$ is self-adjoint, it is enough to consider

[TABLE]

to show that $\widetilde{\bf P}$ is block tridiagonal.

Let us denote ${\bf R}={\bf L}^{-*}$, then from block QR decomposition (53) we have

[TABLE]

Now consider

[TABLE]

Using the action of the propagator on snapshots (61) we can rewrite (130) as

[TABLE]

But ${\bf V}^{*}{\bf U}={\bf L}^{*}$, hence $({\bf v}^{k+l})^{*}\widehat{\bf u}^{i}=({\bf L}^{*})_{k+l,i}$, thus

[TABLE]

The blocks away from the main block diagonal and block sub-diagonal correspond to $l\geq 2$. Since $i\leq k$ we have $k+l>i+1$ and $k+l>|i-1|$, but ${\bf L}^{*}$ is block upper-triangular, thus $({\bf L}^{*})_{k+l,i+1}=({\bf L}^{*})_{k+l,|i-1|}=\boldsymbol{0}$. So we conclude

[TABLE]

i.e. $\widetilde{\bf P}$ is block tridiagonal.

For $\widetilde{\bf B}$ we simply have

[TABLE]

which shows (74) and concludes the proof.