A Reduced Order Model Approach to Inverse Scattering in Lossy Layered Media

We introduce a reduced order model (ROM) methodology for inverse electromagnetic wave scattering in layered lossy media, using data gathered by an antenna which generates a probing wave and measures the time resolved reflected wave. We recast the wave propagation problem as a passive infinite-dimensional dynamical system, whose transfer function is expressed in terms of the measurements at the antenna. The ROM is a low-dimensional dynamical system that approximates this transfer function. While there are many possible ROM realizations, we are interested in one that preserves passivity and in addition is: (1) data driven (i.e., is constructed only from the measurements) and (2) it consists of a matrix with special sparse algebraic structure, whose entries contain spatially localized information about the unknown dielectric permittivity and electrical conductivity of the layered medium. Localized means in the intervals of a special finite difference grid. The main result of the paper is to show with analysis and numerical simulations that these unknowns can be extracted efficiently from the ROM.

Data sharing not applicable to this article as no datasets were generated or analysed during the current study. The results of this article are fully reproducible by following the implementation presented in the appendix.

1. By \({\widehat{u}}(0,s)\) we mean \(\displaystyle \lim _{T \searrow 0} {\widehat{u}}(T,s)\).

2. Passive means that the dynamical system does not generate energy internally.

3. These “truncated spectral measure” measurements are used for convenience, but ROM’s obtained from other matching conditions, such as \(D^{\text {ROM}}_n(s_j) = D(s_j)\) for 2n properly chosen \((s_j)_{j=1}^{2n}\) can be used as well.

4. The regularity assumptions on r(T) and \(\zeta (T)\) can possibly be relaxed, but since we draw conclusions from [9], we use the assumptions made in that study.

5. The spectrum is defined as the set of \(s \in {\mathbb {C}}\) such that the operator is not boundedly invertible.

6. In [3] the grid was called “optimal", but spectrally matched is a more appropriate name.

7. The Weyl function is denoted by M in [9] and the Laplace frequency by \(\rho \).

This research is supported in part by the AFOSR awards FA9550-18-1-0131 and FA9550-20-1-0079, and by ONR award N00014-17-1-2057. The last two authors of this article are working with Rob Remis and Murthy Guddati on the related Project “Krein’s embedding of dissipative data-driven reduced order models”. We would like to thank them for stimulating discussions. We also thank Alex Mamonov and Mikhail Zaslavsky for stimulating discussions on the topic of this paper.

Authors and Affiliations

1. Department of Mathematics, University of Michigan, Ann Arbor, MI, 48109-1043, USA

   Liliana Borcea & Jörn Zimmerling

2. Department of Mathematics, Worcester Polytechnic Institute, Worcester, MA, 01609, USA

   Vladimir Druskin

Corresponding author

Correspondence to Jörn Zimmerling.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

From the Weyl Function to the Transfer Function

The Weyl function \({{\mathcal {W}}}(s)\) is defined^{Footnote 7} in [9] as

where \(\phi (T,s)\) is the solution of

Let us introduce the following, pairwise linearly independent solutions associated with the operator pencil (33): \(\psi (T,s)\), \(\xi (T,s)\) and \(\eta (T,s)\). The first one satisfies

the second one satisfies

and the third one satisfies

It is easy to check that the Wronskian

is constant in T, so we can define

where we used the boundary condition in (101). Similarly, the Wronskian

is constant in T, so we can define

where we used the boundary condition in (102).

Now it follows that the solution of (99) can be written as

and the solution w(T, s) of the Schrödinger problem (29)–(30) is

The latter is because \({\mathscr {L}}_{q,r}(s) w(T,s) = 0\) for \(T \in (0, T_L)\) by construction, and at the boundary we have

Now using (30) we obtain the jump condition

which corresponds to the Dirac delta forcing \(-s \zeta _0 \delta (T)\) in (29).

Solving for \(\psi (T,s)\) in (107) we get

and since \(\phi (0,s) = 1\), the transfer function has the expression

Moreover, taking the T derivative in (111) at \(T = 0\) and using the definition (98) of the Weyl function and the boundary condition (109), we obtain that

This proves Eq. (32).

The poles of the transfer function are the zeroes of the Weyl function and therefore of the Wronskian (106). They correspond to the eigenvalues of the quadratic operator pencil \({\mathscr {L}}_{q,r}(s)\) with domain \({{\mathcal {S}}}_N\) defined in (34). The zeroes of the transfer function are \(s = 0\) and the set of poles of the Weyl function, which are the eigenvalues of the quadratic operator pencil \({\mathscr {L}}_{q,r}(s)\) with domain \({{\mathcal {S}}}_D\) defined in (35).

The Transfer Function for Small Variations of the Loss Function

To use the perturbation theory in [20], consider the first order system formulation (26)–(27). We made the assumption (24) to simplify the boundary conditions satisfied by w(T, s). Similarly, we shall assume in this appendix that

which implies that \({\widehat{w}}(T,s)\) satisfies a homogeneous Neumann boundary condition at \(T = T_L\). The loss function has small variations, so we model it as

Using (115) in (26) we obtain

where \({{\mathcal {L}}}\) is the differential operator (16) and Q(T) is the skew-symmetric multiplication operator defined in terms of \(\zeta (T)\) in (28). The loss function (115) is in the multiplication operator

Note that \(i \big [{{\mathcal {L}}}+ Q(T)\big ]\) acting on the space of functions satisfying the boundary conditions (27) is a self-adjoint (with respect to the Euclidian inner product) indefinite differential operator on a bounded interval and thus has a countable set of real valued eigenvalues with no finite accumulation point [11, Chapter 7]. The squares of these eigenvalues \((-\theta _j^2)_{j \ge 1}\) are the same as the eigenvalues in the Sturm-Liouville problems

and

where

The Sturm-Liouville theory gives that the eigenvalues are simple. Assuming that the eigenfunctions \(\varphi _j(T)\) and \({\widehat{\varphi }}_j\) in (118)–(119) are normalized as

then the eigenfunctions of the operator \(i \big [{{\mathcal {L}}}+ Q(T)\big ]\) are

and the eigenvalues are \(\pm \theta _j\). Equivalently, the eigenvalues of \({{\mathcal {L}}}+ Q(T)\) are \(\pm i \theta _j\).

Now, if we consider the operator \({{\mathcal {L}}}+ Q(T) + R_0\), for the constant loss \(r_0\), the eigenfunctions \(\varvec{\varphi }_j^\pm (T)\) are still orthonormal and are determined by the components of (123) as follows

where \(s = \lambda _j\) is the root of

Indeed, it is easy to check that with \( \varvec{\varphi }_j^+(T) \) given in (124),

is equivalent to

which leads to (125).

Finally, if we add the perturbation \({R}_\alpha (T)-R_0\), which is clearly a bounded operator with \(O(\alpha )\) norm, we can use the analytic perturbation theory in [20] to obtain that the eigenvalues and eigenprojections are analytic for \(\alpha \) in some vicinity of 0. Thus, we can use these eigenprojections to express the wave w(T, s) as a series and obtain the expression (46) of the transfer function.

Passivity and pH Structure

A dynamical system is called passive if it does not generate energy internally [28, 30]. As stated in [28, Section 2], this property is realized if the transfer function satisfies

1. 1.

   \(D({\overline{s}}) = \overline{D(s)}\), for all \(s \in {\mathbb {C}}\).

2. 2.

   D(s) is analytic for \(\text{ Re }(s) > 0\).

3. 3.

   \(D(s) + \overline{D(s)} \ge 0\) for \(\text{ Re }(s) > 0\).

It is obvious from the expression (18) of the transfer function that it satisfies the first condition. The second condition says that the dynamical system is stable i.e., the poles of D(s) are in the left half complex plane. We know from Sect. 2 that these poles are \(\{\lambda _j, \overline{\lambda _j}, ~j \ge 1 \}\), where \(-\lambda _j\) are the eigenvalues of the operator \({{\mathcal {L}}}+ {{\mathcal {Q}}}(T) + R(T)\) acting on the space of functions satisfying the boundary conditions (27). If we let \({{\varvec{V}}}_j(T)\) be the eigenfunctions, then we have

Taking the real part of the inner product with \({{\varvec{V}}}_j(T)\) and using that \({{\mathcal {L}}}+ {{\mathcal {Q}}}(T)\) are skew-symmetric, we get

That \(\text{ Re }(\lambda _j) \le 0\) follows from this equation and the fact that the diagonal multiplication operator R(T) is positive semidefinite.

It remains to verify the third condition, which has the following physical interpretation: Since \(-u(T,s) \mathbf{e}_{x_2}\) is the electric field and \({\widehat{u}}(T,s) \mathbf{e}_{x_1}\) is the magnetic field, the Poynting vector at \(T = 0\), which determines the power flow, is

Thus, the third condition is equivalent to saying that the power flow is into the medium i.e., in the positive range direction.

To check this condition, we use the first order system formulation (26) to write

where in the second line we took the adjoint of the inverse and recalled that \({{\mathcal {L}}}\) and \({{\mathcal {Q}}}(T)\) are skew-symmetric. Now we can write

where the operator

can be factorized as

Substituting in (130) and using that

we get that for \( \text{ Re }(s) > 0\),

where the inequality is because R(T) is positive semidefinite.

This proves that our dynamical system is passive. Checking that it has a pH structure is a straightforward verification of [1, Definition 3].

Derivation of the ROM Transfer Function

Multiplying Eq. (47) by \(\zeta _j\) and (48) by \(-{\widehat{\zeta }}_j\), we get the following linear system

where \({\mathbb {\varvec{T}}}\) is the tridiagonal matrix

We can symmetrize this matrix using the diagonal matrix \( {\varvec{\varGamma }}= \mathrm{diag}({\widehat{\gamma }}_1, \gamma _1, {\widehat{\gamma }}_2, \ldots , \gamma _n), \) so we rewrite (132) as

where

Finally, we can factor out the square root of the diagonal matrix in (132) to get

where

is the matrix given in (54)–(52).

The ROM transfer function is

as stated in Eq. (53).

Proof of Lemma 9

Let us introduce the weighted inner product

The linear operator \(i {{\mathcal {L}}}_\zeta \) defined in (79) with boundary conditions (80) is self-adjoint with respect to this inner product and thus has a countable set of real eigenvalues with no finite accumulation point [11, Chapter 7]. This implies that \({{\mathcal {L}}}_\zeta \) has purely imaginary eigenvalues. In fact, \({{\mathcal {L}}}_\zeta \) is related via a similarity transformation to the operator \({{\mathcal {L}}}+ Q(T)\) studied in Appendix B

so the eigenvalues are the same \(\{\pm i \theta _j, ~~ j \ge 1\}\). The eigenfunctions of \({{\mathcal {L}}}+ Q(T)\) are the vector valued functions (123), which are orthonormal in the Euclidian inner product. These define the eigenfunctions of \({{\mathcal {L}}}_\zeta \)

and the orthonormality relations (83) follow from (122).

The same discussion applies to the linear operator \({{\mathcal {L}}}_\zeta \) with boundary conditions (81). \(~ \Box \)

Proof of Proposition 10

Recall from Sect. 2.1 that the poles of the transfer function D(s) are the eigenvalues \(\{\lambda _j, \overline{\lambda _j}, ~ j\ge 1\}\) of the operator pencil \({\mathscr {L}}_{q,r}(s)\) defined in (33) acting on the space \({{\mathcal {S}}}_N\) defined in (34), whereas the zeroes of D(s) are the eigenvalues \(\{\mu _j, \overline{\mu _j}, ~ j\ge 1\}\) of the operator (33) acting on the space \({{\mathcal {S}}}_D\) defined in (35). Moreover, as explained in Sect. 2 (recall Eqs. (26) and (29)), \({\mathscr {L}}_{q,r}(s)\) is connected to the first order pencil

with \(R_\alpha (T)\) defined in (117) and \({{\mathcal {L}}}_\zeta \) defined in (79) as follows: First, we have the similarity transformation

Second, the first order system

can be written as the second order equation

with potential q(T) defined in (23) and the other way around. This implies that \({\mathscr {L}}_{q,r}(s)\) with domain \({{\mathcal {S}}}_N\) has the same eigenvalues \(\{\lambda _j, \overline{\lambda _j}, ~ j\ge 1\}\) as \({\mathscr {P}}_{\zeta ,r}(s) \) with boundary conditions (80). Similarly, \({\mathscr {L}}_{q,r}(s)\) with domain \({{\mathcal {S}}}_D\) has the same eigenvalues \(\{\mu _j, \overline{\mu _j}, ~ j\ge 1\}\) as \({\mathscr {P}}_{\zeta ,r}(s) \) with boundary conditions (81).

We know that the ROM based estimated functions \(\zeta ^{(n)}(T)\), \({\mathfrak {r}}^{(n)}(T)\) and \(\widehat{{\mathfrak {r}}}^{(n)}(T)\) define an operator pencil

with the following properties:

1. 1.

   The pencils (145) and (141) with boundary conditions (80) have the same first n eigenvalues \(\lambda _j\), for \(j = 1, \ldots , n\). Moreover, the residues \(y_j\), which equal the jumps of the spectral measures, are also the same, for \(j = 1, \ldots , n.\)

2. 2.

   In the case of boundary conditions (81), the eigenvalues of (145) are approximately equal to the zeroes \(\mu _j\) of the transfer function. We denote the error in their approximation by o(1) in the limit \(n \rightarrow \infty \).

Now, using the assumption (86) on the loss function we can write

which is an \(O(\alpha )\) perturbation of \({\mathscr {P}}_{\zeta ,r_0}(s)\). Since the poles and residues of (145) and (146) match to all orders of \(\alpha \), we can use Proposition 8 to obtain from the O(1) matching that the O(1) primary loss must be \(r_0\) and the dual loss is zero. Moreover, the impedance \(\zeta ^{(n)}(T)\) approximates \(\zeta (T)\) as \(n \rightarrow \infty \). Therefore, we can write pointwise in \((0,T_L)\) that

with functions \(\rho ^{(n)}(T)\) and \({\widehat{\rho }}^{(n)}(T)\) independent of \(\alpha \).

Because the same constant \(r_0\) appears in both the pencils (145) and (146), we can subtract it from both problems. This results in a transformation of the eigenvalues, but since these eigenvalues match, they will be transformed the same way by the subtraction of \(r_0\) i.e., they will still match. The new pencils are

and their eigenvalues are

where \(i \theta _j\) and \(i \vartheta _j\) are the purely imaginary eigenvalues of the lossless problem (see Lemma 9). The \(O(\alpha )\) perturbations of these eigenvalues are

and similarly

where \(\varvec{\varPhi }_j^{+}\) and \(\varvec{\varPsi }^+_j,\) for \(j \ge 1\) are the orthonormal eigenfunctions in Lemma 9. Substituting their expressions in these equations gives (89)–(90).

Finally, note that if \(\zeta ^{(n)}(T)\) had an \(O(\alpha )\) error term, then that would reflect in an imaginary perturbation of the eigenvalues, because \({{\mathcal {L}}}_{\zeta }\) and \({{\mathcal {L}}}_{\zeta ^{(n)}}\) are skew-symmetric operators with respect to the weighted inner product \(\left\langle \cdot , \cdot \right\rangle _{\zeta ^{-1},\zeta }\). However, the perturbations (152)–(152) are real valued, so the error in the impedance must be of higher order in \(\alpha \). \(\square \)

Setup for the Numerical Simulations

In this appendix we describe how we generate the data used in the inversion results in Figs. 3, 4, 5 and 6. We also give details on the calculation of the Jacobian used in the Gauss-Newton iteration for solving the optimization problem (96).

To generate the data, we solve the system (5)–(6), with boundary conditions (7), (13), using finite differences on a staggered grid with constant step size \(\tau \), except for the first dual step size, as illustrated in the following sketch:

The primary wave u(T, s) is discretized on the primary grid, at the nodes illustrated in blue and the dual wave \({\widehat{u}}(T,s)\) is discretized on the dual grid, at the nodes illustrated in red. We suppress the s dependence of the discretized waves in the illustration. The boundary conditions are imposed at the pale colored nodes. The travel time domain is normalized at \(T_L=1\) and we use \(N=3000\) grid steps, so that the discretization never falls below 30 points per wavelength. The derivatives are calculated with the standard, two point forward differentiation rule. The truncated measure data can be calculated using the spectral decomposition of the finite differences matrix. However, to better emulate the measurement process, we evaluate \(D(s) = u_1(s) \) in the interval \(s \in [-i\omega _{\mathrm{max}},i\omega _{\mathrm{max}}]\), discretized at 10000 equidistant points, and then use the vectorfit algorithm [15, 16] to extract the poles and residues of D(s). The value of \(\omega _{\mathrm{max}}\) depends on n and it is given in the following table:

The vectorfit algorithm alone does not give a good estimate of the truncated spectral measure transfer function \(D_n^\text {ROM}(s)\) from D(s), because the poles and residues outside the spectral interval of interest have a large contribution, especially when the mean loss \(r_0\) is large. Thus, we proceed as follows: First, we estimate the mean loss \(r_0\) by fitting D(s) at points with \(|s| \gg 1\) with the transfer function

with poles and residues given by the asymptotes of the spectrum:

Once we estimate \(r_0\), we use the vectorfit algorithm to estimate the poles and residues of \(D(s) - D^{\mathrm{as}}(s;r_0)\), for \(s \in [-i\omega _{\mathrm{max}},i\omega _{\mathrm{max}}]\). These then define the truncated spectral measure transfer function \(D_n^\text {ROM}(s)\) of the ROM, according to Eq. (42).

The Jacobian used in the Gauss-Newton iteration for solving problem (96) is approximated numerically using finite differences. We perturb the Fourier coefficients of \(\zeta ^S(T)\) and \(r^S(T)\) by \(\varDelta _{\mathrm{num}}=0.01\) and compute the resulting ROM parameters \({\mathfrak {r}}^{\mathrm{pert}}_j, \widehat{{\mathfrak {r}}} ^{\mathrm{pert}}_j,\zeta ^{\mathrm{pert}}_j\) and \({\widehat{\zeta }}_j^{\mathrm{pert}}\). The entries of the Jacobian corresponding to \(\zeta _j^S\) are \((\zeta _j^S- \zeta ^{\mathrm{pert}}_j ) \varDelta _{\mathrm{num}}^{-1}\), and similarly for the coefficients corresponding to the loss function.

Reprints and permissions