Many important problems in remote sensing, where measurements are not available in the domain of interest (radar imaging, seismic exploration, medical array ultrasound, etc.) lead to inverse scattering problems which can be strongly nonlinear in case of large perturbations of the unknown PDE coefficients. The model based nonlinear optimization which is the method of choice for the solution of such problems can be unreliable and prohibitively expensive. Data driven nonlinear transforms emerged as an attractive alternative, however it was recently shown that the most common ReLU neural networks are intractable for reliable solution of non-parametric inverse problems. Data driven ROMs recently emerged as a feasible option for such problems.

The key of this approach is data-driven computation of the state solution in the domain of interest not available for direct measurements for a black-box model via a nonlinear transform. It implicitly was used for different imaging applications with data-driven ROMs, e.g., see, [2, 3, 1]. Here we give its new explicit formulation allowing simple analysis and clear connection with preceding work.

We assume the following wave model problem for a domain \(\Omega \subset {R}^d\)

\begin{equation} \label {waveeq} u_{tt} -\Delta u + q(x) u= 0 \ \mbox {in} \ \Omega \times [ 0,\infty ) \end{equation}

with initial conditions

\begin{eqnarray} \label {waveeqic} u ( t=0) &=& g \ \mbox {in}\ \Omega \\ u_t ( t=0) &=& 0\ \mbox {in}\ \Omega \label {waveeqid} \end{eqnarray}

where \(g(x)\) is an initial condition representing a localized source near the boundary, and we assume homogeneous Neumann boundary conditions on the spatial boundary \(\partial \Omega \). We assume \(q(x)\ge 0\) is our unknown potential, not necessarily small, but with compact support. The exact forward solution to (1-3) is

\begin{equation} \label {exactsolution} u(x,t) = \cos {(\sqrt {A}t)} g(x), \end{equation}

where

\[A={-\Delta +q(x)},\]

with the square root and cosine are defined via the spectral theorem. This solution is assumed to be unknown, except near the receivers.

Assume we measure the SISO transfer function at the receiver collocated with the at the \(2n-1\) evenly spaced time steps \(t= k\tau \) for \(k=0,\dots , 2n-2\), modeled by

\begin{eqnarray} \label {datadef} F(k\tau ) &= &\int _\Omega g(x) u(x,k\tau ) dx\nonumber \\ &=& \int _\Omega g(x) \cos {(\sqrt {A }k\tau )} g(x) dx . \end{eqnarray}

The inverse problem we consider is as follows: Given

\[ \{ F(k \tau ) \} \ \mbox { for } \ k=0,\dots , 2n-2,\]

reconstruct \(q\). This generally nonlinear problem becomes a simple linear problem if \(u_k = u(x,k\tau )\) is known in the entire domain, that may not be accessible to the direct measurement.

So our objective here is to compute, from the measured data \(F(k\tau )\) only, approximations of the internal snapshots \(u_k = u(x,k\tau )\) for \(k=0,\ldots n-1\) assuming that \(q(x)\) is unknown.

We introduce Gramian matrix \(M\) with with elements

\begin{equation} \label {massmatrixj} M_{kl}= \int _\Omega u_k u_l dx \end{equation}

for \(k,l =0,\ldots ,n-1\), can be written as

\begin{equation} \label {massmatrixj2} M_{kl} = \int _\Omega g(x) \cos {(\sqrt {A}k\tau )} \cos {(\sqrt {A }l\tau )} g(x) dx. \end{equation}

thanks to the formula (4). Then, from (5), (7), and the cosine angle sum formula, one has

\begin{equation} \label {massfromdataj} M_{kl} = \frac {1}{2} \left ( F((k-l)\tau ) + F((k+l)\tau ) \right ), \end{equation}

so Gramian \(M\) can obtained directly from the data [2] . Formula (7) is the Chebyshev moment problem yielding \(M\) given by 8 as the sum of Toeplitz and Hankel matrices.

Now, let

\[U = [ u_1(x), \ldots , u_n(x) ] \]

be a row vector of the true snapshots, so that we can write

\begin{equation} \label {truemass} M = \int _\Omega U^\top U\in R^{n\times n}. \end{equation}

Consider also the background field \(u^0(x,t)\), the solution to (1-3) with \(q(x)=0\), which we assume that we know. Let

\[U_0 = [ u_1^0(x), \ldots , u_n^0(x) ] \]

be a row vector of the background snapshots \(u_k^0 = u^0(x,k\tau )\), and let

\[ M_0 = \int _\Omega U_0^\top U_0 \in R^{n\times n}\]

be the background mass matrix. We want to obtain approximation

\[U\approx {\bf U} = U_0 T \]

satisfying (9) via projection, i.,e., condition

\begin{equation} \label {datamass} \int _\Omega {\bf U}^\top {\bf U} = M \end{equation}

that yields

\begin{equation} \label {ROMmass} M = T^\top M_0 T. \end{equation}

Equation (11) was inspired by the celebrated Marchenko-Gelfand-Levitan-Krein (MGLK) Volterra equation, e.g.,see [5] and references wherein, with \(T\) being so-called transmutation matrix. Its discrete analogy first appeared in study of connection between the discrete MGLK and Lanczos algorithms [6]. A critical observation is that the waves in background (\(q=0\)) and true (unknown \(q\) ) media travel with the same speed, so thanks to the causality principle, \(T\) is upper triangular matrix. This restriction leads to its uniqueness, that can be shown by direct calculations.

The Cholesky factorization of \(M\) constitutes the nonlinear part of the data transform. The internal solution generated via SISO data was successfully used for radar imaging applications, however it had significant artifacts due to luck of aperture [4].

For seismic exploration and medical array ultrasound the SISO data (5) can be replaced by the square MIMO transfer function

\begin{eqnarray} \label {datadefMIMO} F(k\tau ) &= &\int _\Omega G(x) u(x,k\tau ) dx\nonumber \\ &=& \int _\Omega g(x) \cos {(\sqrt {A }k\tau )} G(x) dx \in R^{m\times m}, \end{eqnarray}

where \(G(x)=[g_1(x),\ldots , g_m(x)]\) is the row vector of \(m\) transmitters collocated with receivers. The Proposition 1 can be replaced by its block analogy for the data given by (13), which was implicitly used in e.g., [3].

Another important application, the synthetic aperture radars (SAR) used for imaging from airborne platforms can only access \(\hbox {diag }F\), and for computing block-transmutation matrix they require data-completion.

Finally we outline the list of computational linear algebra bottlenecks in the proposal framework:

• • Fast Cholesky factorization and spectral decomposition of sum of block Hankel and Toeplitz matrices

• • Lifting partial data matrices to full square MIMO data, e.g., from diagonal as in the SAR framework

• • Efficient truncation or correction of spurious non-Hamiltonian modes (with negative eigenvalues, appeared due to measurement errors or inaccuracies of the above mentioned data-completion or lifting) of data-driven Gramians

• • Estimation of non-strictly triangular data-driven transmutation matrices

• • Extension to problems with dissipation

References

• [1]  Liliana Borcea, Josslin Garnier Alexander V. Mamonov, and Jörn Zimmerling, When data driven reduced order modeling meets full waveform inversion SIAM Review , 66 (3) (216): 501-532.

• [2]  Druskin, V.r, A. Mamonov, A. Thaler, and M. Zaslavsky. ”Direct, nonlinear inversion algorithm for hyperbolic problems via projection-based model reduction.” SIAM Journal on Imaging Sciences 9, no. 2 (2016): 684-747.

• [3]  Druskin, V., Mamonov, A. V., Zaslavsky, M. (2018). A nonlinear method for imaging with acoustic waves via reduced order model backprojection. SIAM Journal on Imaging Sciences, 11(1), 164-196.

• [4]  Druskin, V., Moskow, S. and Zaslavsky, M., 2024. Reduced Order Modeling Inversion of Monostatic Data in a Multi-scattering Environment. SIAM Journal on Imaging Sciences, 17(1), pp.334-350.

• [5]  T. Habashy. A generalized Gelfand–Levitan–Marchenko integral equation. Inverse Problems, 7, pp.703–711, 1991.

• [6]  F. Natterer. A discrete Gelfand–Levitan theory. Electronic, from the author’s web page.