Direct reconstruction for acoustic inverse Born scattering
June 04, 2026
Abstract
We consider the inverse medium scattering problem for the Helmholtz equation in two dimensions, i.e., the task to recover a compactly supported penetrable two-dimensional scatterer from full knowledge of the associated far field data or, equivalently, the far field operator. Although this problem is uniquely solvable, it is severely ill-posed and nonlinear. In the regime of weak scattering, the Born approximation yields a linearized relation between the contrast and the far field data, thus overcoming the second difficulty. This linear setting allows to build on recent work on linearized electrical impedance tomography, which relies on triangular Zernike decompositions, to derive an explicit reconstruction formula that expresses the expansion coefficients of the contrast in terms of those of the far field data. By choosing the expansion functions appropriately, the resulting system matrix decouples into separate (infinite) triangular systems for the spatial angular frequencies in the contrast. Consequently, each of these systems can be solved independently by performing forward substitutions. Our numerical experiments indicate that this approach, combined with an adequate regularization method, remains effective even when applied to full nonlinear far field data beyond the Born regime.
Mathematics subject classifications (MSC2010): 35R30, 78A46, 47A52
Keywords: Inverse medium scattering, Helmholtz equation, far field data, Born approximation, direct reconstruction, QR factorization
Short title: Inverse Born scattering
1 Introduction↩︎
In inverse scattering, one seeks for an unknown material parameter, described by a contrast function, from associated scattered near or far field data, with applications in, e.g., seismology, medical imaging, radar technology, geophysical exploration and nondestructive testing. Comprehensive surveys of this active research area are provided in [1], [2]. In this article, we focus on the inverse medium scattering problem for time-harmonic acoustic waves, a setting governed by the Helmholtz equation. We assume to observe far field data for plane wave incident fields along all possible illumination directions at a fixed frequency and present in two dimensions a direct, fast reconstruction algorithm in the framework of the Born approximation, following the ideas in [3], [4] for electrical impedance tomography (EIT).
It is well-known that full knowledge of the far field data uniquely determines the contrast. However, inverse scattering is ill-posed as small perturbations in the observed data may result in huge errors in the reconstructed contrast. Moreover, the data depend nonlinearly on the contrast due to multiple scattering effects.
There are two classes of reconstruction algorithms for this problem, either aiming to reconstruct (1) the entire contrast function or (2) the shape or boundary of its compact support. Methods focusing on solving problem (2) are referred to as qualitative methods. Prominent examples include the Linear Sampling Method, the Factorization Method, and monotonicity-based techniques (see, e.g., [5]–[9]). Quantitative methods addressing problem (1) typically either solve a computationally expensive regularized optimization problem or seek to invert a linearized model in the weak scattering regime, with the hope that the resulting inversion procedure remains effective even when the measured data contain nonlinear effects (see, e.g., [10]–[12]). We particularly emphasize the recent low-rank method [12] in the latter category as it is conceptually close to the approach proposed in this paper.
We assume that the contrast is supported in an a priori known \({\boldsymbol{c}}\)-centered ball \(B_R({\boldsymbol{c}})\) of radius \(R\). By interpreting the available far field data in the standard way as the kernel of the far field operator on \({L^2}({S^1})\), the linearized forward map can be presented as a linear operator from \({L^2}(B_R({\boldsymbol{c}}))\) to the space of Hilbert–Schmidt operators \(\mathrm{HS}({L^2}({S^1}))\), which allows its parametrization as an infinite matrix after introducing bases for these separable Hilbert spaces. In particular, expanding the angular dependence of the unknown contrast and the observed Born far field data in Fourier bases reveals a decoupling of the system in the angular direction: a certain angular frequency in the contrast only affects a corresponding diagonal in the Born far field operator expanded with respect to a (modulated) Fourier basis. This observation motivates the application of concepts in numerical linear algebra to solve each of the decoupled infinite linear systems that connect the radial behavior of the contrast for a certain angular frequency to a diagonal in the data matrix. Our proposed algorithm essentially applies a QR factorization in an offline stage to the angular subsystems, allowing one to only solve (small) triangular systems when the data become available.
Another key component is the intrinsic low-rank structure of the measured far field data caused by the super-exponential decay of Bessel functions with respect to their order, providing a natural justification for truncating the infinite-dimensional systems while essentially maintaining the retrievable information. This is particularly important for our algorithm since the QR factorization step applies a Gram–Schmidt orthogonalization in an appropriately weighted \(L^2\) topology on a finite interval to certain products of Bessel functions, a process that is guaranteed to become unstable if run too long due to the aforementioned super-exponential decay. The sparsity of far field data has been previously investigated for related inverse problems in [13], [14], and foundational results on this research direction can be found in [15]–[17].
The low-rank structure of the inverse medium scattering problem has also been analyzed from a different perspective in [12], [18]. In these works, the far field data are expanded in terms of eigenfunctions of a restricted Fourier integral operator, known as the generalized prolate spheroidal wave functions.
The remainder of this article structures as follows. In Section 2, we provide theoretical background on acoustic inhomogeneous medium scattering as well as on its linearization given by the Born approximation. Section 3 then derives a representation for the linearized system, decoupling it into triangular systems along the radial direction over the angular frequencies in the contrast. Finally, in Section 4, the essential support of the expansion coefficients of the far field data is utilized to justify a truncation-based regularization for inverting the infinite-dimensional triangular systems. At this stage, we also discuss numerical instability in the construction of the radial basis functions. Numerical results are presented in Section 5, where we test our method using (noisy) linearized Born and nonlinear full far field data and benchmark it against the low-rank method [12] and the MATLAB’s built-in nonuniform fast Fourier transform (NUFFT). We close the article with some conclusions.
2 Inhomogeneous medium scattering↩︎
We consider scattering of time-harmonic acoustic waves by a compactly supported penetrable object lying in a homogeneous background medium. Let \(\kappa>0\) denote the wave number and assume the incident field to be a plane wave \[\tag{1} \begin{equation} u^i({\boldsymbol{x}},{\boldsymbol{d}}) \,:=\, \mathrm{e}^{\mathrm{i}\kappa{\boldsymbol{x}}\cdot{\boldsymbol{d}}}\,, \qquad {\boldsymbol{x}}\in{\mathbb{R}}^2\,, \end{equation} that propagates along an illumination direction {\boldsymbol{d}}\in{S^1}:=\{{\boldsymbol{x}}\in{\mathbb{R}}^2\,:\,|{\boldsymbol{x}}|=1\}. In the following, all dependencies on {\boldsymbol{d}} are marked by a second argument. The incident field u^i(\,\cdot\,,{\boldsymbol{d}}) hits and interacts with a compactly supported penetrable scatterer \Omega\subset {\mathbb{R}}^2, modeled by the exterior support of a real-valued contrast function q\in{L^{\infty}}({\mathbb{R}}^2) satisfying q>-1 a.e.~on \Omega and q=0 a.e.~on {\mathbb{R}}^2\setminus\overline{\Omega}. This interaction generates a total field u(\,\cdot\,,{\boldsymbol{d}})\in H^1_{{\mathrm{loc}}}({\mathbb{R}}^2) that solves \begin{equation} \tag{2} \Delta u(\,\cdot\,,{\boldsymbol{d}}) + \kappa^2(1+q)u(\,\cdot\,,{\boldsymbol{d}}) \,=\, 0 \qquad \text{in } {\mathbb{R}}^2\,, \end{equation} such that the scattered field u^s(\,\cdot\,,{\boldsymbol{d}}):=u(\,\cdot\,,{\boldsymbol{d}})-u^i(\,\cdot\,,{\boldsymbol{d}}) satisfies the Sommerfeld radiation condition \begin{equation} \tag{3} \lim_{r\to\infty} \sqrt{r} \Bigl( \frac{\partial u^s}{\partial r}({\boldsymbol{x}},{\boldsymbol{d}}) - \mathrm{i}\kappa u^s({\boldsymbol{x}},{\boldsymbol{d}}) \Bigr) \,=\, 0 \,, \qquad r=|{\boldsymbol{x}}|\rightarrow\infty \,, \end{equation}\] uniformly with respect to the direction \(\widehat{{\boldsymbol{x}}}= {\boldsymbol{x}}/|{\boldsymbol{x}}|\in{S^1}\).
It is known that the unique weak solution \(u(\,\cdot\,,{\boldsymbol{d}})\in H^1_{{\mathrm{loc}}}({{\mathbb{R}}^2})\) of 1 (see, e.g., [19]) satisfies the Lippmann–Schwinger integral equation \[\label{eq:LippmannSchwinger} u(\,\cdot\,,{\boldsymbol{d}}) \,=\, u^i(\,\cdot\,,{\boldsymbol{d}}) + \kappa^2 \int_\Omega q({\boldsymbol{y}})u({\boldsymbol{y}},{\boldsymbol{d}}) \Phi(\,\cdot\,-{\boldsymbol{y}}) \, \mathop{\mathrm{d\!}}{\boldsymbol{y}} \,=:\, u^i(\,\cdot\,,{\boldsymbol{d}}) + L_q u(\,\cdot\,,{\boldsymbol{d}}) \qquad \text{in } \Omega \,,\tag{4}\] with \(L_q:{L^2}(\Omega)\rightarrow{L^2}(\Omega)\) and \(\Phi({\boldsymbol{x}}):=\mathrm{i}/4\, H^{(1)}_0(\kappa|{\boldsymbol{x}}|)\), \({\boldsymbol{x}}\not={\boldsymbol{0}}\), denoting the fundamental solution to the Helmholtz equation in free space at wave number \(\kappa\) (see, e.g., [19]). The function \(H^{(1)}_0\) is the Hankel function of the first kind and order zero. Due to the asymptotic behavior of \(H^{(1)}_0\) for a large argument, the scattered field \(u^s(\,\cdot\,,{\boldsymbol{d}})\) fulfills the asymptotic far field expansion \[\label{eq:far95field95expansion} u^s({\boldsymbol{x}},{\boldsymbol{d}}) \,=\, \frac{\mathrm{e}^{\mathrm{i}\pi/4}}{\sqrt{8\pi}} \frac{\mathrm{e}^{\mathrm{i}\kappa r}}{\sqrt{\kappa r}} u^\infty(\widehat{{\boldsymbol{x}}},{\boldsymbol{d}}) + O\bigl(r^{-\frac{3}{2}}\bigr) \,, \qquad r=|{\boldsymbol{x}}| \to \infty \,,\tag{5}\] uniformly with respect to the observation direction \(\widehat{{\boldsymbol{x}}}={\boldsymbol{x}}/|{\boldsymbol{x}}|\in{S^1}\). Here, the far field pattern \({u^\infty\in {L^2}({S^1}\times{S^1})}\) is given by \[\label{eq:full95far95field} u^\infty(\widehat{{\boldsymbol{x}}},{\boldsymbol{d}}) \,=\, \kappa^2 \int_\Omega q({\boldsymbol{y}}) u({\boldsymbol{y}},{\boldsymbol{d}}) \mathrm{e}^{-\mathrm{i}\kappa \widehat{{\boldsymbol{x}}}\cdot {\boldsymbol{y}}} \, \mathop{\mathrm{d\!}}{\boldsymbol{y}}\,, \qquad \widehat{{\boldsymbol{x}}}\in{S^1}\,,\tag{6}\] (see, e.g., [19]), and it determines the associated far field operator \[\label{eq:FarfieldOperator} F: {L^2}({S^1}) \rightarrow{L^2}({S^1}) \,, \quad \bigl(Fg\bigr)(\widehat{{\boldsymbol{x}}}) \,:=\, \int_{{S^1}} u^\infty(\widehat{{\boldsymbol{x}}},{\boldsymbol{d}}) g({\boldsymbol{d}})\, \mathop{\mathrm{d\!}}s({\boldsymbol{d}})\,,\tag{7}\] which maps superpositions of plane wave incident fields to the far field patterns of the associated scattered fields. This operator is well-known to be compact, normal and of trace class (see, e.g., [20]). In particular, it is Hilbert–Schmidt on \({L^2}({S^1})\), i.e., \(F \in \mathrm{HS}({L^2}({S^1}))\), and it can thus be viewed as an infinite-dimensional matrix after fixing an orthonormal system for \({L^2}({S^1})\).
In this work, we are interested in the inverse medium scattering problem of recovering the contrast \(q\) from the knowledge of the associated far field operator \(F\), or equivalently, from the knowledge of the far field data \(u^\infty(\widehat{{\boldsymbol{x}}},{\boldsymbol{d}})\) for all \(\widehat{{\boldsymbol{x}}},{\boldsymbol{d}}\in{S^1}\). This broadly studied problem is known to be uniquely solvable (see, e.g., [19]) but severely ill-posed and nonlinear as multiple scattering effects have to be taken into account.
2.1 Linearization by considering the Born approximation↩︎
If \(L_q\) from 4 satisfies \(\|L_q\|_{\mathscr{L}({L^2}(\Omega))}\ll1\) the total field \(u(\,\cdot\,,{\boldsymbol{d}})\) can be accurately approximated by its Born approximation \(u_B(\,\cdot\,,{\boldsymbol{d}}):= (I+L_q)u^i(\,\cdot\,,{\boldsymbol{d}})\) for all \({\boldsymbol{d}}\in{S^1}\). This follows by only accounting for the first two terms in the Neumann series of \((I - L_q)^{-1}\) when viewing the Lippmann–Schwinger equation 4 as a fixed-point equation for \(u(\,\cdot\,,{\boldsymbol{d}})\). The field \(u_B(\,\cdot\,,{\boldsymbol{d}})\) solves a source problem for the Helmholtz equation, \[\Delta u_B(\,\cdot\,,{\boldsymbol{d}}) + \kappa^2u_B(\,\cdot\,,{\boldsymbol{d}}) \,=\, -\kappa^2 qu^i(\,\cdot\,,{\boldsymbol{d}})\qquad \text{in } {\mathbb{R}}^2\,,\] with the associated scattered field \(u^s_B(\,\cdot\,,{\boldsymbol{d}}):=u_B(\,\cdot\,,{\boldsymbol{d}})-u^i(\,\cdot\,,{\boldsymbol{d}})\) fulfilling the Sommerfeld radiation condition 3 . The Born far field pattern \(u_B^\infty(\,\cdot\,,{\boldsymbol{d}})\) is defined by replacing the scattered field \(u^s\) by its Born approximation \(u^s_B\) in the asymptotic expansion 5 , which leads to the representation \[\label{eq:BornFF} u_B^\infty(\widehat{{\boldsymbol{x}}},{\boldsymbol{d}}) \,:=\, \kappa^2 \int_{\Omega} q({\boldsymbol{y}})\mathrm{e}^{-\mathrm{i}\kappa (\widehat{{\boldsymbol{x}}}-{\boldsymbol{d}})\cdot{\boldsymbol{y}}}\, \mathop{\mathrm{d\!}}{\boldsymbol{y}}\,, \qquad \widehat{{\boldsymbol{x}}},{\boldsymbol{d}}\in{S^1}\,.\tag{8}\] The associated Born far field operator \[\label{eq:BornffO} F_B:L^2({S^1}) \rightarrow L^2({S^1})\,, \qquad (F_Bg)(\widehat{{\boldsymbol{x}}}) \,:=\, \int_{{S^1}} u_B^\infty(\widehat{{\boldsymbol{x}}},{\boldsymbol{d}}) g({\boldsymbol{d}}) \, \mathop{\mathrm{d\!}}s({\boldsymbol{d}})\,\tag{9}\] is also a Hilbert–Schmidt operator. We refer to inverse Born scattering as the task to recover the contrast \(q\) from the knowledge of the Born far field operator \(F_B\).
Formula 8 indicates that knowing \(u_B^\infty(\widehat{{\boldsymbol{x}}},{\boldsymbol{d}})\) for all \(\widehat{{\boldsymbol{x}}},{\boldsymbol{d}}\in{S^1}\) is equivalent to knowing the two-dimensional Fourier transform of \(q\) on the origin-centered disk of radius \(2 \kappa\) \[B_{2\kappa}({\boldsymbol{0}})\,=\,\{ \kappa(\widehat{{\boldsymbol{x}}}-{\boldsymbol{d}})\;:\; \widehat{{\boldsymbol{x}}},{\boldsymbol{d}}\in{S^1}\} \subset {\mathbb{R}}^2.\] However, sampling \(\widehat{{\boldsymbol{x}}}\) and \({\boldsymbol{d}}\) uniformly on \({S^1}\) demonstrates that the natural sampling pattern for the Fourier transform of \(q\) over \(B_{2\kappa}({\boldsymbol{0}})\) is non-uniform, as visualized in Figure 1. Since \(q\) has compact support, its Fourier transform is an analytic function by virtue of the Paley–Wiener theorem, and the unique continuation principle thus reveals that the contrast \(q\) is uniquely determined by the knowledge of the Born far field operator \(F_B\). However, it follows from [21] that the eigenvalues of \(F_B\) decay at the same rate as those of \(F\), implying that inverse Born scattering is ill-posed to the same extend as the original inverse medium scattering problem.
Inspired by the treatment of the linearized continuum model of EIT in [3], [22], we adopt formulas 8 –9 as the starting point for deriving a direct reconstruction formula for inverse Born scattering. Since we plan to treat the contrast \(q\) as an element of \(L^2(\Omega)\) and \(F_B\) as a Hilbert–Schmidt operator on \(L^2({S^1})\), i.e., as an element of \(\mathrm{HS}({L^2}({S^1}))\), we complete this section by noting that the forward map of Born scattering is indeed bounded between these spaces.
Proposition 1. The linear forward map \[T_B: q \, \mapsto \, F_B\] defined by 8 –9 is bounded from \(L^2(\Omega)\) to \(\mathrm{HS}({L^2}({S^1}))\) with \[\| T_B \|_{\mathscr{L}(L^2(\Omega),\mathrm{HS}({L^2}({S^1}))} \, \leq \, 2 \pi \kappa^2 \sqrt{|\Omega|}\,.\]
Proof. Since \[(F_B g)(\widehat{{\boldsymbol{x}}}) \, = \, \kappa^2 \int_{{S^1}} \int_{\Omega} q({\boldsymbol{y}})\mathrm{e}^{-\mathrm{i}\kappa (\widehat{{\boldsymbol{x}}}-{\boldsymbol{d}})\cdot{\boldsymbol{y}}}\, \mathop{\mathrm{d\!}}{\boldsymbol{y}}\, g({\boldsymbol{d}}) \, \mathop{\mathrm{d\!}}s({\boldsymbol{d}})\,, \qquad \widehat{{\boldsymbol{x}}}\in{S^1}\,,\] the basic theory on Hilbert–Schmidt integral operators yields (see, e.g., [23]) \[\label{eq:forward95kernel} \| T_B q \|_{\mathrm{HS}({L^2}({S^1}))}^2 \leq \kappa^4 \int_{{S^1}} \int_{{S^1}} \bigg | \int_{\Omega} q({\boldsymbol{y}})\mathrm{e}^{-\mathrm{i}\kappa (\widehat{{\boldsymbol{x}}}-{\boldsymbol{d}})\cdot{\boldsymbol{y}}}\, \mathop{\mathrm{d\!}}{\boldsymbol{y}}\bigg|^2 \, \mathop{\mathrm{d\!}}s({\boldsymbol{d}}) \, \mathop{\mathrm{d\!}}s(\widehat{{\boldsymbol{x}}}) \, .\tag{10}\] As \[\bigg | \int_{\Omega} q({\boldsymbol{y}})\mathrm{e}^{-\mathrm{i}\kappa (\widehat{{\boldsymbol{x}}}-{\boldsymbol{d}})\cdot{\boldsymbol{y}}}\, \mathop{\mathrm{d\!}}{\boldsymbol{y}}\bigg|^2 \leq |\Omega| \| q \|_{L^2(\Omega)}^2 \,\] by the Cauchy–Schwarz inequality, the assertion follows by integrating twice over \({S^1}\) in 10 and taking the square root. ◻
3 Angularly decoupled triangular systems for Born scattering↩︎
We start by formulating Born scattering in a matrix form between orthonormal bases of \(L^2(B_1({\boldsymbol{0}}))\) and \(\mathrm{HS}({L^2}({S^1}))\), then consider a specific choice for the radial parts of the basis for \(L^2(B_1({\boldsymbol{0}}))\), and finally show that our choices lead to angularly decoupled triangular systems for determining the expansion coefficients of the contrast. In the following, we assume prior knowledge of a discoidal region of interest (ROI) \(B_R({\boldsymbol{c}})\subset{\mathbb{R}}^2\) containing \(\Omega\).
3.1 Matrix representation↩︎
We begin by deriving the expansion coefficients for the observed far field data with respect to a modulated Fourier basis following [13]. By linear substitution, we rewrite the Born far field pattern 8 as \[u_B^\infty(\widehat{{\boldsymbol{x}}},{\boldsymbol{d}}) \,=\, \kappa^2 \int_{B_R({\boldsymbol{0}})} q({\boldsymbol{y}}+{\boldsymbol{c}})\mathrm{e}^{-\mathrm{i}\kappa{\boldsymbol{c}}\cdot(\widehat{{\boldsymbol{x}}}-{\boldsymbol{d}})}\mathrm{e}^{-\mathrm{i}\kappa (\widehat{{\boldsymbol{x}}}-{\boldsymbol{d}})\cdot{\boldsymbol{y}}}\, \mathop{\mathrm{d\!}}{\boldsymbol{y}}\,, \qquad \widehat{{\boldsymbol{x}}},{\boldsymbol{d}}\in{S^1}\,.\] Due to the Jacobi–Anger expansion (see, e.g., [2]), the plane wave term in this formula can be expanded as \[\mathrm{e}^{-\mathrm{i}\kappa (\widehat{{\boldsymbol{x}}}-{\boldsymbol{d}})\cdot{\boldsymbol{y}}} \,=\, 2\pi \sum_{m,n\in{\mathbb{Z}}} \mathrm{i}^{n-m} \mathrm{e}^{\mathrm{i}(n-m)\arg{\boldsymbol{y}}}J_m(\kappa|{\boldsymbol{y}}|)J_n(\kappa|{\boldsymbol{y}}|){\boldsymbol{e}}_m(\widehat{{\boldsymbol{x}}})\overline{{\boldsymbol{e}}_{n}({\boldsymbol{d}})}\,,\] where \(({\boldsymbol{e}}_m)_{m\in{\mathbb{Z}}}:=(\exp(\mathrm{i}m\arg(\,\cdot\,))/\sqrt{2\pi})_{m\in{\mathbb{Z}}}\) is the standard Fourier basis of \({L^2({S^1})}\) and \(J_m\) denotes the Bessel function of the first kind and order \(m\). By introducing the modulated Fourier system \(({\boldsymbol{e}}^{\boldsymbol{c}}_m)_{m\in{\mathbb{Z}}}:=(\exp(-\mathrm{i}\kappa{\boldsymbol{c}}\cdot(\,\cdot\,)){\boldsymbol{e}}_m)_{m\in{\mathbb{Z}}}\), which also forms an orthonormal basis for \({L^2}({S^1})\), we conclude that \[\mathrm{e}^{-\mathrm{i}\kappa{\boldsymbol{c}}\cdot(\widehat{{\boldsymbol{x}}}-{\boldsymbol{d}})}\mathrm{e}^{-\mathrm{i}\kappa (\widehat{{\boldsymbol{x}}}-{\boldsymbol{d}})\cdot{\boldsymbol{y}}} \,=\, 2\pi \sum_{m,n\in{\mathbb{Z}}} \mathrm{i}^{n-m} \mathrm{e}^{\mathrm{i}(n-m)\arg{\boldsymbol{y}}}J_m(\kappa|{\boldsymbol{y}}|)J_n(\kappa|{\boldsymbol{y}}|){\boldsymbol{e}}^{\boldsymbol{c}}_m(\widehat{{\boldsymbol{x}}})\overline{{\boldsymbol{e}}^{\boldsymbol{c}}_{n}({\boldsymbol{d}})}\,.\] This immediately yields a representation for the Born far field operator corresponding to \(q\) in the orthonormal basis \((\langle\,\cdot\,,{\boldsymbol{e}}^{\boldsymbol{c}}_n\rangle_{{L^2}({S^1})}{\boldsymbol{e}}^{\boldsymbol{c}}_m)_{m,n\in{\mathbb{Z}}}\) of \(\mathrm{HS}({L^2}({S^1}))\) with expansion coefficients \[\begin{align} a_{m,n} &\,:=\, \langle F_B{\boldsymbol{e}}^{\boldsymbol{c}}_n,{\boldsymbol{e}}^{\boldsymbol{c}}_m\rangle_{L^2({S^1})} \notag\\ &\,=\, 2\pi \kappa^2\mathrm{i}^{n-m} \int_{B_R({\boldsymbol{0}})} q({\boldsymbol{y}}+{\boldsymbol{c}})\mathrm{e}^{\mathrm{i}(n-m)\arg{\boldsymbol{y}}}J_m(\kappa|{\boldsymbol{y}}|)J_n(\kappa|{\boldsymbol{y}}|)\, \mathop{\mathrm{d\!}}{\boldsymbol{y}}\notag\\ &\,=\, (2\pi)^{3/2} \kappa^2\mathrm{i}^{n-m} \int_{B_R({\boldsymbol{0}})} q({\boldsymbol{y}}+{\boldsymbol{c}})\overline{{\boldsymbol{e}}_{m-n}(\widehat{{\boldsymbol{y}}})}J_m(\kappa|{\boldsymbol{y}}|)J_n(\kappa|{\boldsymbol{y}}|)\, \mathop{\mathrm{d\!}}{\boldsymbol{y}}\notag\\ &\,=\, (2\pi)^{3/2} (\kappa R)^2\mathrm{i}^{n-m} \int_{B_1({\boldsymbol{0}})} q(R{\boldsymbol{y}}+{\boldsymbol{c}})\overline{{\boldsymbol{e}}_{m-n}(\widehat{{\boldsymbol{y}}})}J_m(\kappa R|{\boldsymbol{y}}|)J_n(\kappa R|{\boldsymbol{y}}|)\, \mathop{\mathrm{d\!}}{\boldsymbol{y}}\,, \label{eq:exp95coeff95FF} \end{align}\tag{11}\] where \(\widehat{{\boldsymbol{y}}}= {\boldsymbol{y}}/|{\boldsymbol{y}}|\in{S^1}\).
We continue by introducing an orthonormal basis \((\Psi_{j,k})_{j\in{\mathbb{Z}},k\in{\mathbb{N}}_0}\) of \(L^2(B_1({\boldsymbol{0}}))\) for expanding the shifted and rescaled contrast \(q(R(\,\cdot\,)+{\boldsymbol{c}})\) to enable representing the forward map (cf. Proposition 1) \[\label{eq:forward} L^2(B_1({\boldsymbol{0}})) \ni q(R(\,\cdot\,)+{\boldsymbol{c}}) \mapsto F_B \in \mathrm{HS}({L^2}({S^1}))\tag{12}\] between orthonormal bases. Observing that the kernel of the integral transform 11 , mapping \(q(R(\,\cdot\,)+{\boldsymbol{c}})\) to \((a_{m,n})_{m,n \in {\mathbb{Z}}}\), separates into a radial \(|{\boldsymbol{y}}|\)-dependent part and an angular \(\widehat{{\boldsymbol{y}}}\)-dependent factor that is given by the standard Fourier basis, it is natural to search for the basis of \(L^2(B_1({\boldsymbol{0}}))\) in the form \[\label{eq:def95Psi} \Psi_{j,k}({\boldsymbol{y}}):={\boldsymbol{e}}_j(\widehat{{\boldsymbol{y}}})R_k^{|j|}(|{\boldsymbol{y}}|), \qquad {\boldsymbol{y}}\in B_1({\boldsymbol{0}})\,.\tag{13}\] Here, the radial functions \((R_k^{|j|})_{j\in{\mathbb{Z}},k\in{\mathbb{N}}_0}\) are chosen such that \((\Psi_{j,k})_{j\in{\mathbb{Z}},k\in{\mathbb{N}}_0}\) forms an orthonormal basis for \({L^2}(B_1({\boldsymbol{0}}))\); due to the orthonormality of the standard Fourier basis on \({L^2}({S^1})\), it straightforwardly follows that the necessary and sufficient condition is that \((R_k^{|j|})_{k\in{\mathbb{N}}_0}\) is an orthonormal basis for the weighted \(L^2\)-space \[L^2_r(0,1) = \bigg\{ f: (0,1) \to {\mathbb{C}}\text{ measurable} \; : \; \int_0^1 |f(r)|^2 r \, \mathop{\mathrm{d\!}}r< \infty \bigg\}\] for each \(j \in {\mathbb{Z}}\). We may thus expand \[\label{eq:exp95qc} q(R(\,\cdot\,)+{\boldsymbol{c}}) \,=\, \sum_{j\in{\mathbb{Z}}} q_j \,=\, \sum_{j\in{\mathbb{Z}}}\sum_{k=0}^\infty c_{j,k} \Psi_{j,k} \qquad \text{for } c_{j,k}\,:=\, \big \langle q(R(\,\cdot\,)+{\boldsymbol{c}}), \Psi_{j,k} \big \rangle_{{L^2}(B_1({\boldsymbol{0}}))}\,,\tag{14}\] where we call \(q_j\) the \(j\)-th angular frequency in \(q(R(\,\cdot\,)+{\boldsymbol{c}})\). Note that \(q_j\), \((a_{m,n})_{m,n\in{\mathbb{Z}}}\) and \((c_{j,k})_{j\in{\mathbb{Z}},k\in{\mathbb{N}}_0}\) all depend on \({\boldsymbol{c}}\) and \(R\), but we suppress this dependence to improve readability. The specific choice of the radial basis functions \((R_k^{|j|})_{j\in{\mathbb{Z}},k\in{\mathbb{N}}_0}\) will be considered in Subsection 3.2 below.
We address the inverse medium scattering problem by recovering the coefficients \((c_{j,k})_{j\in{\mathbb{Z}},k\in{\mathbb{N}}_0}\) of the unknown contrast from the knowledge of the coefficients \((a_{m,n})_{m,n\in{\mathbb{Z}}}\) of the observed Born far field operator. Inserting 14 into 11 and exploiting the orthonormality of the Fourier basis reveals the underlying infinite-dimensional system matrix \((b_{m,n}^{j,k})_{m,n \in {\mathbb{Z}}, j\in{\mathbb{Z}},k\in{\mathbb{N}}_0}\), that characterizes the Born forward map 12 with respect to the chosen orthonormal basis: \[\label{eq:LS} a_{m,n} \,=\, \sum_{j\in{\mathbb{Z}}}\sum_{k=0}^\infty b_{m,n}^{j,k} c_{j,k}\,, \qquad m,n\in{\mathbb{Z}}\,,\tag{15}\] for \[\label{eq:matrix95coeff} b^{j,k}_{m,n} \,:=\, \begin{cases} (2\pi)^{3/2}(\kappa R)^2(-\mathrm{i})^{j} \big\langle R^{|j|}_{k}, J_m(\kappa R\,\cdot\,)J_{m-j}(\kappa R\,\cdot\,) \big\rangle_{L^2_r(0,1)} & \text{if } n=m-j\,, \\[1mm] 0 & \text{else}. \end{cases}\tag{16}\] Having a closer look at 15 –16 , we make the following central structural observations.
Remark 2.
As in [3], [22] for the case of EIT, 16 uncovers a decoupling of angular frequencies, which enables considering the diagonals \((a_{m,m-j})_{m\in{\mathbb{Z}}}\) of the data matrix separately for \(j\in{\mathbb{Z}}\). Indeed, all available information on the \(j\)-th angular frequency \(q_j\) in \(q(R(\,\cdot\,)+{\boldsymbol{c}})\) is included in the \((-j)\)-th diagonal of the data matrix, which allows to solve for the associated expansion coefficients \(( c_{j,k} )_{k \in {\mathbb{N}}_0}\) from \[\label{eq:LSj} a_{m,m-j} \,=\, \sum_{k=0}^\infty b_{m,m-j}^{j,k} c_{j,k}\,, \qquad m\in{\mathbb{Z}}\, ,\tag{17}\] for each \(j \in {\mathbb{Z}}\).
For each \(j\in{\mathbb{Z}}\), a half of the equations in 17 are redundant due to a symmetry in \(m \in {\mathbb{Z}}\), which means that only a half of the equations needs to be considered and the remaining data can be used for noise filtering. Indeed, the reciprocity relation in the Born far field data \[u^\infty_B(\widehat{{\boldsymbol{x}}},{\boldsymbol{d}}) \,=\, u^\infty_B(-{\boldsymbol{d}},-\widehat{{\boldsymbol{x}}}) \qquad \text{for all } \widehat{{\boldsymbol{x}}},{\boldsymbol{d}}\in{S^1}\,,\] yields [13] \[\label{eq:a95symmetry} a_{-(m-j),-m} \,=\, (-1)^j\, a_{m,m-j} \qquad \text{for all } m, j \in{\mathbb{Z}}\,.\tag{18}\] The same can alternatively be deduced from (see, e.g., [24]) \[J_{-m}J_{-(m-j)} \,=\, (-1)^jJ_mJ_{m-j} \qquad \text{for all } m, j \in{\mathbb{Z}},\] which also gives \[b^{j,k}_{-(m-j),-m} \,=\, (-1)^j\,b^{j,k}_{m,m-j} \qquad \text{for all }m, j \in{\mathbb{Z}}\text{ and } k \in {\mathbb{N}}_0,\] by virtue of 16 . Consequently, the \((-j)\)-th diagonal of the data matrix \((a_{m,n})_{m,n\in{\mathbb{Z}}}\) is symmetric up to the factor \((-1)^j\) with respect to the (possibly virtual) element \(a_{j/2,-j/2}\), and the analogous conclusion also holds for the system matrix \((b_{m,n}^{j,k})_{m,n \in {\mathbb{Z}}}\) with fixed \(j\) and \(k\). In consequence, we do not discard any unique equations in 17 if we only consider \[\label{eq:LSj95reciprocity} a_{m,m-j} \,=\, \sum_{k=0}^\infty b_{m,m-j}^{j,k} c_{j,k}\,, \qquad m\geq \tfrac j 2\,,\tag{19}\] where \(a_{m,m-j}\) could be replaced by the averaged data \[\label{eq:noise95filt} \widetilde{a}_{m,m-j}:=(a_{m,m-j}+(-1)^ja_{-m+j,-m})/2\,,\tag{20}\] without altering the equations, to improve the signal-to-noise ratio. \(\lozenge\)
It remains to construct orthonormal bases \((R_k^{|j|})_{k\in{\mathbb{N}}_0}\), \(j \in {\mathbb{Z}}\), for \(L^2_r(0,1)\) so that inverting 15 , or equivalently 19 , becomes straightforward. Following [3], [22], we aim for a choice that makes the system matrix in 19 lower triangular for every \(j \in {\mathbb{Z}}\). The radial Zernike bases employed in [3], [22] are unsuitable for this purpose, but it turns out they can be replaced by bases generated through a Gram–Schmidt orthonormalization process of the products of Bessel functions appearing in 16 .
3.2 Choice of the radial bases↩︎
For \(j,m\in{\mathbb{Z}}\), we define \[\label{eq:def95Pjm} P_{m}^{j}(r) \,:=\, J_{m}(\kappa R r) J_{m-j}(\kappa R r) \,,\qquad r\in(0,1)\,,\tag{21}\] which appear in 16 and will act as our initial, i.e., non-orthonormalized, radial basis functions. To ease the notation, we set \[{\mathbb{Z}}_{\geq c} := \{ k \in {\mathbb{Z}}\; : \; k \geq c \}\] for \(c \in {\mathbb{R}}\).
Proposition 3. For each \(j \in {\mathbb{Z}}\), the functions \((P_{m}^{j})_{m \in {\mathbb{Z}}_{\geq j/2}}\) are linearly independent and their linear span is dense in \(L^2_r(0,1)\).
Proof. We start by proving that \((P_{m}^{j})_{m \in {\mathbb{Z}}_{\geq j/2}}\) are linearly independent. For \(m \geq j/2\) and any \(j \in {\mathbb{Z}}\), the lowest order term in the converging origin-centered power series representation of \(P_{m}^{j}(r)\) behaves as \(r^{2m-j}\) (e.g., [25]). Hence, each function in \((P_{m}^{j})_{m \in {\mathbb{Z}}_{\geq j/2}}\) has its own distinct polynomial behavior close to the origin, which proves that they are linearly independent.
Let \(\rho \in L^2_r(0,1)\). We prove the assertion on the density by showing that \(\rho\) can be orthogonal in \(L^2_r(0,1)\) to all functions in the set \((P_{m}^{j})_{m \in {\mathbb{Z}}_{\geq j/2}}\) only if it vanishes. Define a shifted and scaled contrast as \(q(R {\boldsymbol{y}}+ {\boldsymbol{c}}) = {\boldsymbol{e}}_j(\widehat{{\boldsymbol{y}}}) \rho(|{\boldsymbol{y}}|)\), \({\boldsymbol{y}}\in B_1({\boldsymbol{0}})\), for a fixed but arbitrary \(j \in {\mathbb{Z}}\), which according to Remark 2 means that the data matrix \((a_{m,n})_{m,n\in{\mathbb{Z}}}\) only has nonzero elements on its \((-j)\)-th diagonal. By inserting our choice of \(q(R (\,\cdot\,) + {\boldsymbol{c}})\) into 11 and integrating over \({S^1}\), we get the representation (cf. 16 ) \[a_{m,m-j} = (2\pi)^{3/2}(\kappa R)^2(-\mathrm{i})^{j} \langle \rho, P_{m}^{j} \rangle_{L^2_r(0,1)}, \qquad m \in {\mathbb{Z}},\] for the elements on the \((-j)\)-th diagonal. If \(\rho\) is orthogonal to all functions in \((P_{m}^{j})_{m \in {\mathbb{Z}}_{\geq j/2}}\), the \((-j)\)-the diagonal of \((a_{m,n})_{m,n\in{\mathbb{Z}}}\) is thus empty due to the symmetry 18 , meaning that the data matrix altogether vanishes. This means that our \(q(R (\,\cdot\,) + {\boldsymbol{c}})\) is in the nullspace of the forward operator 12 , and thus \(q(R (\,\cdot\,) + {\boldsymbol{c}})\) is zero almost everywhere by the unique solvability of the considered inverse Born scattering problem. This completes the proof. ◻
For any fixed \(m, j\in{\mathbb{Z}}\), \[\label{eq:P95jm} P_{m-j}^{-j} \,=\, J_{m-j}J_{(m-j)-(-j)} \,=\, J_mJ_{m-j} \,=\, P_m^{j}\,,\tag{22}\] from which it follows that for any \(j \in {\mathbb{Z}}\), \[\label{eq:lin95ependence95Pjm} (P_{m}^{j} )_{m \in {\mathbb{Z}}_{\geq j/2}} = (P_{m}^{-j} )_{m \in {\mathbb{Z}}_{\geq -j/2}} = (P_{m}^{|j|})_{m \in {\mathbb{Z}}_{\geq |j|/2}},\tag{23}\] with the functions in these sets given in the same ordering with respect to increasing \(m\). Hence, we only need to consider the radial bases \((P_{m}^{|j|})_{m \in {\mathbb{Z}}_{\geq |j|/2}}\), \(j \in {\mathbb{Z}}\), in what follows.
Now we are ready to properly introduce our basis \((\Psi_{j,k})_{j\in{\mathbb{Z}},k\in{\mathbb{N}}_0}\) for \(L^2(B_1({\boldsymbol{0}}))\) by defining the orthonormal basis \((R_k^{|j|})_{k \in {\mathbb{N}}_0}\) of \(L^2_r(0,1)\) in 13 for each \(j \in {\mathbb{Z}}\).
Definition 1. The radial basis functions \((R_k^{|j|})_{k \in {\mathbb{N}}_0}\), \(j \in {\mathbb{Z}}\), in 13 are defined by applying the Gram–Schmidt ortogonalization process with respect to the inner product of \(L^2_r(0,1)\) to the functions \((P_{m}^{|j|})_{m \in {\mathbb{Z}}_{\geq |j|/2}}\). That is, \[\label{eq:def95radial95basis} R_k^{|j|} \,:=\, \frac{\widetilde{R}_k^{|j|}}{\big \|\widetilde{R}_k^{|j|} \big \|_{L^2_r(0,1)}} \qquad \text{with} \; \widetilde{R}_k^{|j|} \,:=\, P_{k+\lceil|j|/2\rceil}^{|j|} - \sum_{m=0}^{k-1} \big \langle P_{k+\lceil|j|/2\rceil}^{|j|},R_{m}^{|j|} \big \rangle_{L^2_r(0,1)}R_{m}^{|j|}\tag{24}\] for \(k \in {\mathbb{N}}_0\). \(\lozenge\)
Since by Proposition 3 the functions \((P_{m}^{|j|})_{m \in {\mathbb{Z}}_{\geq |j|/2}}\) are linearly independent and their linear span is dense in \(L^2_r(0,1)\), the set \((R_k^{|j|})_{k \in {\mathbb{N}}_0}\) is a well-defined orthonormal basis of \(L^2_r(0,1)\) for each \(j \in {\mathbb{Z}}\), which is precisely the requirement for the construction leading to 19 to be valid. On the negative side, the procedure in 24 is numerically unstable, which we will address in Section 4 below.
3.3 Angularly decoupled triangular systems↩︎
Let us then examine how our specific choice for the radial bases simplifies the system 19 . For \(j \in {\mathbb{N}}_0\), \(m \in {\mathbb{Z}}_{\geq j/2}\) and \(k \in {\mathbb{N}}_0\), formula 16 gives \[\label{eq:reduced95b} b_{m,m-j}^{j,k} = (2\pi)^{3/2}(\kappa R)^2(-\mathrm{i})^{j} \big\langle R^{|j|}_{k}, P_m^{j} \big\rangle_{L^2_r(0,1)} = (2\pi)^{3/2}(\kappa R)^2(-\mathrm{i})^{j} \big\langle R^{|j|}_{k}, P_m^{|j|} \big\rangle_{L^2_r(0,1)}.\tag{25}\] The right side of 25 vanishes if \(k > m-\lceil |j|/2\rceil = m-\lceil j/2\rceil\) since \(R^{|j|}_{k}\) is orthogonal to the subspace of \(L^2_r(0,1)\) spanned by the first \(k\) functions in \((P_{m}^{|j|})_{m \in {\mathbb{Z}}_{\geq |j|/2}}\) due to the Gram–Schmidt process 24 . On the other hand, for \(-j \in {\mathbb{N}}_0\) and \(m \in {\mathbb{Z}}_{\geq j/2}\), it follows from 22 that \[\begin{align} b_{m,m-j}^{j,k} &= (2\pi)^{3/2}(\kappa R)^2(-\mathrm{i})^{j} \big\langle R^{|j|}_{k}, P_m^{j} \big\rangle_{L^2_r(0,1)} \nonumber \\[1mm ] &= (-\mathrm{i})^{2j} (2\pi)^{3/2}(\kappa R)^2 (-\mathrm{i})^{-j} \big\langle R^{|j|}_{k}, P_{m-j}^{-j} \big\rangle_{L^2_r(0,1)} = (-1)^{j} b_{m-j,m}^{-j,k}. \label{eq:reduced95b95minus} \end{align}\tag{26}\] Since in this case \(P_{m-j}^{-j} = P_{m+|j|}^{|j|}\), the orthogonalization process 24 dictates that \(b_{m,m-j}^{j,k} = (-1)^{j} b_{m-j,m}^{-j,k}\) vanishes when \[k > m + |j| - \lceil |j|/2\rceil = m + \lfloor |j| /2\rfloor = m - \lceil j/2\rceil.\] Together with 19 , these conclusions yield angularly decoupled triangular systems for determining the expansion coefficients \((c_{j,k})_{j \in {\mathbb{Z}}, k \in {\mathbb{N}}_0}\) from the knowledge of the data matrix \((a_{m,n})_{m,n\in {\mathbb{Z}}}\): \[\label{eq:triangular0} a_{m,m-j} \,=\, \sum_{k=0}^{m-\lceil j/2\rceil} b_{m,m-j}^{j,k}c_{j,k}, \qquad m \in {\mathbb{Z}}_{\geq j/2} \, ,\tag{27}\] for \(j \in {\mathbb{Z}}\). In order to recast 27 for a fixed \(j\in{\mathbb{Z}}\) in a matrix-vector form, we introduce the infinite-dimensional vectors \[\label{eq:full95vectors} {\boldsymbol{c}}^j \,:=\, \big[c_{j,k-1} \big]_{k=1}^{\infty} \qquad\text{and}\qquad {\boldsymbol{a}}^j \,:=\, \big[ a_{(m-1)+\lceil j/2\rceil,(m-1)-\lfloor j/2\rfloor} \big]_{m=1}^{\infty},\tag{28}\] corresponding to the \(j\)-th angular frequency in \(q(R(\,\cdot\,)+{\boldsymbol{c}})\) and a half of the \((-j)\)-th diagonal in the data matrix, respectively, as well as the lower triangular infinite-dimensional system matrix \({\boldsymbol{F}}^j:=[F^j_{m,k}]_{m,k=1}^\infty\) given componentwise as \[\label{eq:F95comp} F^{j}_{m,k} \,=\, \begin{cases} (2\pi)^{3/2}(-\mathrm{i})^j(\kappa R)^2 \big \langle R_{k-1}^{|j|}, P_{(m-1)+\lceil |j|/2\rceil}^{|j|} \big\rangle_{L^2_r(0,1)} & \text{if } k\leq m \,,\\[1mm] 0 & \text{else}. \end{cases}\tag{29}\]
For each \(j\in{\mathbb{Z}}\), the resulting system \[\label{eq:infinite95system} {\boldsymbol{F}}^{j} {\boldsymbol{c}}^j\,=\, {\boldsymbol{a}}^j\tag{30}\] corresponds to 27 when \(m\) runs from \(\lceil j/2\rceil\) to infinity and can be solved through a forward substitution (not accounting for instability). Note that according to 25 and 26 , the latter term in the inner product in 29 should, in fact, be \(P_{(m-1)+\lceil j/2\rceil}^{j}\), but employing 22 for \(j < 0\), one deduces \[P_{(m-1)+\lceil j/2\rceil}^{j} = P_{(m-1)+\lceil j/2\rceil - j }^{-j} = P_{(m-1) + \lceil |j|/2\rceil}^{|j|},\] which allows the presented form.
These considerations lead to the following reconstruction formula that is our main result.
Theorem 4. Denote by \((\Psi_{j,k})_{j\in{\mathbb{Z}},k\in{\mathbb{N}}_0}\) the orthonormal system for \({L^2}(B_1({\boldsymbol{0}}))\) defined by 13 and 24 , and let \((a_{m,n})_{m,n\in{\mathbb{Z}}}\) be the given expansion coefficients of the Born far field operator as in 11 . Then, the expansion coefficients of the shifted and scaled contrast \(q(R(\,\cdot\,)+{\boldsymbol{c}})\) with respect to \((\Psi_{j,k})_{j\in{\mathbb{Z}},k\in{\mathbb{N}}_0}\), as given in 14 , can be computed separately for each \(j \in {\mathbb{Z}}\) via a recursion with respect to \(k\): \[\label{eq:forward95substitution} c_{j,k} \,=\, \frac{1}{(2\pi)^{3/2}(\kappa R)^2(-\mathrm{i})^j \big\| \widetilde{R}_{k}^{|j|} \big \|_{L^2_r(0,1)}} a_{k+\lceil j/2\rceil, k-\lfloor j/2\rfloor} - \sum_{i=0}^{k-1} \frac{\big \langle R_{i}^{|j|}, P_{k+\lceil |j|/2\rceil}^{|j|} \big \rangle_{L^2_r(0,1)}}{\big\| \widetilde{R}_{k}^{|j|} \big\|_{L^2_r(0,1)} } \, c_{j, i}\tag{31}\] for \(k \in {\mathbb{N}}_0\). Here, \((\widetilde{R}_{k}^{|j|})_{k \in {\mathbb{N}}_0}\) are the unnormalized orthogonal basis functions from 24 and the sum in 31 is viewed to be empty for \(k=0\).
Proof. Fix \(j \in {\mathbb{Z}}\). Solving the lower triangular system 30 via forward substitution directly gives \[c_{j,k} \,=\, \frac{1}{(2\pi)^{3/2}(\kappa R)^2(-\mathrm{i})^j \big \langle R_{k}^{|j|}, P_{k+\lceil |j|/2\rceil}^{|j|} \big \rangle_{L^2_r(0,1)}} a_{k+\lceil j/2\rceil, k-\lfloor j/2\rfloor} - \sum_{i=0}^{k-1} \frac{\big \langle R_{i}^{|j|}, P_{k+\lceil |j|/2\rceil}^{|j|} \big \rangle_{L^2_r(0,1)}}{\big \langle R_{k}^{|j|}, P_{k+\lceil |j|/2\rceil}^{|j|}\big \rangle_{L^2_r(0,1)}} c_{j, i} \,\] for \(k \in {\mathbb{N}}_0\). Representing \(P_{k+\lceil |j|/2\rceil}^{|j|}\) as in 24 and utilizing the orthogonality of \((R_{k}^{|j|})_{k \in {\mathbb{N}}_0}\) allows to replace the inner product in the denominator by the norm of \(\widetilde{R}_{k}^{|j|}\), which is nonzero by construction. This completes the proof. ◻
Observe that all inner products and norms appearing in 31 have already been computed (and stored) when running the Gram–Schmidt process in 24 . Hence, applying the reconstruction formula 31 is essentially for free, assuming the Gram–Schmidt process has been run offline prior to having the data in hand. It is also worth noting that our approach to solving the system 19 is essentially an infinite-dimensional QR factorization, with 31 corresponding to solving the triangular system defined by the “R” matrix.
Example 1. As an example we consider \(q=\chi_{B_r({\boldsymbol{0}})}\) for some \(0<r<1\), with \(\chi_{(\, \cdot \, )}\) denoting the characteristic function of a given set. From 11 we conclude that the expansion coefficients of the associated Born far field operator are given by \[a_{m,n} \,=\, \begin{cases} 2\pi\kappa ^2\|J_m(\kappa|\,\cdot\,|)\|^2_{L^2(B_r({\boldsymbol{0}}))} \,=\, 2\pi\|J_m(|\,\cdot\,|)\|^2_{L^2(B_{\kappa r}({\boldsymbol{0}}))} & \text{if } m=n\,, \\ 0 & \text{else}, \end{cases}\] i.e., by an infinite-dimensional diagonal matrix. By [24], we can rewrite \[a_{m,n} \,=\, \begin{cases} 2\pi^2(\kappa r)^2 \left( J^2_m(\kappa r) - J_{m-1}(\kappa r)J_{m+1}(\kappa r)\right) & \text{if } m=n\,, \\ 0 & \text{else}. \end{cases}\] Thus, for this specific choice of \(q\) we have an explicit formula for the data vector in 28 , namely \[{\boldsymbol{a}}^j \,=\, \begin{cases} \big [2\pi^2(\kappa r)^2 \left( J^2_{m-1}(\kappa r) - J_{m-2}(\kappa r)J_{m}(\kappa r)\right) \big]_{m=1}^{\infty} & \text{if } j=0\,, \\[1mm] {\boldsymbol{0}}& \text{else}. \end{cases}\] Consequently, only one infinite-dimensional triangular system for \(j=0\) has to be solved. \(\lozenge\)
Example 2. We visualize some elements of the orthonormal basis 13 of \({L^2}(B_1({\boldsymbol{0}}))\) for \(\kappa R=2.5\) and \(\kappa R=10\), with the radial components computed numerically via the Gram–Schmidt process 24 . We restrict the angular frequency to \(j \in \{0, \dots, 2 N \}\) and the radial index to \(k \in \{ 0, \dots, N - \lceil j/2 \rceil \}\) with \(N = 3\). According to \(\eqref{eq:forward95substitution}\), the associated expansion coefficients of a scaled and shifted contrast, together with those for the corresponding negative frequencies \(j \in \{-2 N, \dots, -1 \}\), can be determined from the knowledge of the truncated data matrix \(( a_{m,n} )_{m,n=-N}^{N}\) (or, more precisely, from the knowledge of slightly more than a half of it). We use a Gauss–Legendre quadrature with \(N_r=100\) nodes for evaluating the integrals involved in the Gram–Schmidt process 24 . The resulting \((N+1)^2 = 16\) basis functions are shown in Figure 2 for \(\kappa R=2.5\) and in Figure 3 for \(\kappa R=10\). The functions in Figure 2 are qualitatively similar to the corresponding Zernike polynomials employed as the spatial basis functions for EIT in [3], [22], but the ones in Figure 3 systematically have more oscillations in the radial direction and more finer details close to the center of the unit disk compared to the Zernike polynomials. To the best of our knowledge, no families of functions obtained via the orthonormalization of such products of Bessel functions have previously been documented in the literature.
Figure 2: The real parts of the basis functions \((\Psi_{j,k})_{j\in\{0,\ldots,2N\}, k\in\{0,\ldots,N-\lceil j/2\rceil\}}\) for \(\kappa R=2.5\) and \(N=3\). First row: \(\Psi_{0,0}\), \(\Psi_{0,1}\), \(\Psi_{0,2}\) , \(\Psi_{0,3}\). Second row: \(\Psi_{1,0}\), \(\Psi_{1,1}\), \(\Psi_{1,2}\) , \(\Psi_{2,0}\). Third row: \(\Psi_{2,1}\), \(\Psi_{2,2}\), \(\Psi_{3,0}\) , \(\Psi_{3,1}\). Fourth row: \(\Psi_{4,0}\), \(\Psi_{4,1}\), \(\Psi_{5,0}\) , \(\Psi_{6,0}\). The color scale is the same in all subfigures..
Figure 3: The real parts of the basis functions \((\Psi_{j,k})_{j\in\{0,\ldots,2N\}, k\in\{0,\ldots,N-\lceil j/2\rceil\}}\) for \(\kappa R=10\) and \(N=3\). First row: \(\Psi_{0,0}\), \(\Psi_{0,1}\), \(\Psi_{0,2}\) , \(\Psi_{0,3}\). Second row: \(\Psi_{1,0}\), \(\Psi_{1,1}\), \(\Psi_{1,2}\) , \(\Psi_{2,0}\). Third row: \(\Psi_{2,1}\), \(\Psi_{2,2}\), \(\Psi_{3,0}\) , \(\Psi_{3,1}\). Fourth row: \(\Psi_{4,0}\), \(\Psi_{4,1}\), \(\Psi_{5,0}\) , \(\Psi_{6,0}\). The color scale is the same in all subfigures..
\(\lozenge\)
4 Regularization of the infinite-dimensional systems↩︎
As elaborated in [13], almost all expansion coefficients \((a_{m,n})_{m,n \in {\mathbb{Z}}}\) of the given far field data are close to zero and thus negligible. More precisely, it is reasonable to estimate that \[\label{eq:truncation} a_{m,n} \,=\, 0 \qquad \text{for } |m|,|n|>N\tag{32}\] with \(N\) chosen slightly larger than \(\kappa R\). The aforementioned work also provides explicit bounds for the approximation error introduced at this stage. The nonzero structure of the expansion coefficients for the observed far field data under such a truncation is visualized in Figure 4 (left), with only the elements needed for 31 included (i.e., without employing the noise filtering step in 20 ).
Figure 4: Left: Expansion coefficients \((a_{m,n})_{|m|,|n|\leq N}\) that are taken into account in the inversion (cf. 31 ) after introducing the truncation index \(N\) (blue), with a representative diagonal \((a_{m,m-j})_{|m|,|m-j|\leq N}\) corresponding to the angular index \(j=-2\) highlighted (red). Right: Orthonormalization error 34 in the Gram–Schmidt process for different values of \(\kappa R\) as functions of the truncation index \(N\)..
To make the effect on the individual diagonals explicit, we rewrite this assumption as \[a_{m+\lceil j/2\rceil, m-\lfloor j/2\rfloor} \,=\, 0 \qquad \begin{cases} \text{for } |j|>2N \,, \\ \text{for } -2N \leq j < 0 \,,\, m>N+\lfloor j/2\rfloor=N-\lceil|j|/2\rceil\,, \\ \text{for } 0 \leq j \leq 2N \,,\, m>N-\lceil j/2\rceil=N-\lceil|j|/2\rceil\, \end{cases}\] for \(m\in{\mathbb{N}}_0\). Consequently, after introducing a truncation index \(N\), only \(4N+1\) systems for \(j\in\{-2N,\ldots,2N\}\) are to be solved instead of the infinite number of systems in 30 , and for a particular angular index \(j\), the corresponding system becomes \((N-\lceil|j|/2\rceil+1)\)-dimensional, cf. Example 2.
For \(j\in\{-2N,\ldots,2N\}\), we define the truncated vectors (cf. 28 ) \[{\boldsymbol{c}}^{j,N} \,:=\, \big[c_{j,k-1}\big]_{k=1}^{N+1-\lceil|j|/2\rceil}\,,\; {\boldsymbol{a}}^{j,N} \,:=\, \big[a_{(m-1)+\lceil j/2\rceil,(m-1)-\lfloor j/2\rfloor} \big]_{m=1}^{N+1-\lceil|j|/2\rceil} \,\in\, {\mathbb{C}}^{N+1-\lceil|j|/2\rceil}\] and the lower triangular finite-dimensional system matrix \[{\boldsymbol{F}}^{j,N} \, = \, \big[F^j_{m,k} \big]_{m,k=1}^{N+1-\lceil|j|/2\rceil} \,\in\, {\mathbb{C}}^{(N+1-\lceil|j|/2\rceil)\times(N+1-\lceil|j|/2\rceil)} \, .\] Solving the resulting finite-dimensional systems \[\label{eq:truncated95system} {\boldsymbol{F}}^{j,N}{\boldsymbol{c}}^{j,N} \,=\, {\boldsymbol{a}}^{j,N}, \qquad j\in\{-2N,\ldots,2N\},\tag{33}\] in place of the infinitely many infinite-dimensional systems of type 30 constitutes the first step of our regularization scheme. If no further regularization is introduced, the forward substitution formula 31 still applies to solving the individual systems in 33 for \(k \in \{0, \dots, N-\lceil|j|/2\rceil \}\).
Remark 5. The choice of the truncation index \(N\) involves a trade-off between two competing criteria. It should be large enough to retain the relevant information in the observed far field data, yet not so large that the Gram–Schmidt orthonormalization becomes numerically unstable. Indeed, due to the super-exponential decay of the Bessel function \(J_n(t)\) for fixed \(t>0\) as \({n\rightarrow\infty}\), choosing \(N\) too large leads to divisions by values that are numerically indistinguishable from zero during the normalization step. To illustrate this phenomenon, we construct for fixed \(\kappa R>0\) and \(N\in{\mathbb{N}}\) the corresponding \((N+1)^2\) Gram–Schmidt basis functions \({(R_k^{j})_{j=0,\ldots,2N, k=0,\ldots,N-\lceil j/2\rceil}}\) as defined by 24 . For the numerical integration, we use \(N_r=100\) Gauss–Legendre quadrature nodes and weights \((r_i,\omega_i)_{i=0,\ldots,N_r-1}\). By arranging the involved function evaluations for each \(j\) as elements of a matrix, we obtain \[{\boldsymbol{Q}}^{j,N} \,:=\, \begin{pmatrix} R_0^{j}(r_0) & R_1^{j}(r_0) & \cdots & R_{N-\lceil j/2\rceil}^{j}(r_0) \\ R_0^{j}(r_1) & R_1^{j}(r_1) & \cdots & R_{N-\lceil j/2\rceil}^{j}(r_1) \\ \vdots & \vdots & \ddots & \vdots \\ R_0^{j}(r_{N_r-1}) & R_1^{j}(r_{N_r-1}) & \ldots & R_{N-\lceil j/2\rceil}^{j}(r_{N_r-1}) \end{pmatrix} \in {\mathbb{R}}^{N_r\times (N+1-\lceil j/2\rceil)}\,,\] whose columns are ideally orthonormal with respect to the inner product defined by the diagonal weight matrix \({\boldsymbol{W}}=\mathop{\mathrm{diag}}(\omega_{0}r_0,\ldots,\omega_{N_r-1}r_{N_r-1})\in{\mathbb{R}}^{N_r\times N_r}\). Consequently, we may quantify how much the constructed basis functions deviate on average from forming orthonormal systems via evaluating the mean error \[\label{eq:GSO95error} \varepsilon_{\text{GSO}} \,:=\, \frac{1}{N+1} \sqrt{\sum_{j=0}^{2N}\big\|({\boldsymbol{Q}}^{j,N})^\top {\boldsymbol{W}}{\boldsymbol{Q}}^{j,N} -{\boldsymbol{I}}_{N+1-\lceil j/2\rceil} \big\|^2_{\mathrm F}}\tag{34}\] where \({\boldsymbol{I}}_{N+1-\lceil j/2\rceil}\in{\mathbb{R}}^{(N+1-\lceil j/2\rceil)\times (N+1-\lceil j/2\rceil)}\) is the identity matrix and \(\|\,\cdot\,\|_{\mathrm F}\) denotes the Frobenius norm. This error is shown in Figure 4 (right) for \(N\in\{1,\ldots,40\}\) and \(\kappa R\in\{5,10,15,20\}\) on a semi-logarithmic scale. All error curves exhibit a regime of very low error for small \(N\), followed by a sharp increase near \(N \approx \kappa R\), after which they plateau. This behavior suggests that choosing \(N\) significantly larger than \(\kappa R\) is not reasonable. The rule \(N:=\lceil \mathrm{e}\kappa R/2\rceil\) used in [13] turns out to be too large for our purposes, but \(N:=\lceil\kappa R\rceil\) proves to be a reliable choice in our numerical experiments. \(\lozenge\)
The truncated systems 33 can be rearranged into a single block-diagonal system \[\label{eq:block95system} {\boldsymbol{F}}^N {\boldsymbol{c}}^N \,=\, {\boldsymbol{a}}^N\tag{35}\] by introducing \[\begin{align} {\boldsymbol{F}}^N \,&:=\, \mathop{\mathrm{diag}}\!\big({\boldsymbol{F}}^{-2N,N},\ldots,{\boldsymbol{F}}^{2N,N} \big) \,\in\, {\mathbb{C}}^{M\times M}\,, \\[1mm] {\boldsymbol{a}}^N \,&:=\, \big[{\boldsymbol{a}}^{-2N,N}; \ldots;{\boldsymbol{a}}^{2N,N}\big] \,\in\, {\mathbb{C}}^M\,, \\[1mm] {\boldsymbol{c}}^N \,&:=\, \big[{\boldsymbol{c}}^{-2N,N};\ldots ;{\boldsymbol{c}}^{2N,N} \big] \,\in\, {\mathbb{C}}^M\,, \end{align}\] where semicolon denotes vertical concatenation and \[M \,:=\, \sum_{j=-2N}^{2N} \left(N+1-\left\lceil\tfrac{|j|} 2\right\rceil\right) \,=\, (N+1)+4 \sum_{j=1}^N j \,=\, (N+1)(2N+1)\,.\] In our numerical experiments documented in Section 5 below, we further employ a truncated singular value decomposition (SVD) for the system 35 to handle full nonlinear far field data and/or additive noise. The block-diagonal structure of the system matrix \({\boldsymbol{F}}^N\) allows its SVD to be completely characterized by SVDs of the individual small blocks \({\boldsymbol{F}}^{j,N}\), \(j=\{-2N,\ldots,2N\}\), thereby avoiding the need to compute a high dimensional SVD. See [3] for more detailed analysis on applying a truncated SVD to a similar block-diagonal system in the framework of EIT.
5 Numerical examples↩︎
In this section, we demonstrate the functionality of our proposed method by numerical examples for both Born far field data and full far field data, with and without additive noise. Moreover, we compare the performance of our method to the recently introduced low-rank method for solving the inverse Born scattering problem [12] and with the MATLAB’s built-in NUFFT, as it is described in [26].
We assume the ability to sample the the Born far field data at \(2L\in{\mathbb{N}}\) equiangular illumination and measurement directions. That is, we assume the availability of the (noisy) matrix \[\label{eq:observation} \frac{\pi}{L} \left[ u_B^\infty(\widehat{{\boldsymbol{x}}}_m,{\boldsymbol{d}}_n) \right]_{1\leq m,n\leq 2L} \in{\mathbb{C}}^{2L\times 2L}\tag{36}\] with \[\widehat{{\boldsymbol{x}}}_l \,=\, {\boldsymbol{d}}_l \,=\, (\cos(\varphi_l),\sin(\varphi_l))^\top \,,\qquad \varphi_l \,=\, \frac{\pi (l-1)}{L} \,,\qquad l=1,\ldots,2L\,.\] Here, \(2L\in{\mathbb{N}}\) is chosen large enough to resolve all relevant information in the far field data, i.e., \(L\) has to be larger than \(\kappa\) times the radius of the smallest origin-centered ball containing the whole scatterer \(\Omega\) (cf. 32 and [13]). The two-dimensional fast Fourier transform of the matrix \[\label{eq:meas95a} \frac{\pi}{L}\left[ \mathrm{e}^{-\mathrm{i}\kappa{\boldsymbol{c}}\cdot({\boldsymbol{d}}_n-\widehat{{\boldsymbol{x}}}_m)} u_B^\infty(\widehat{{\boldsymbol{x}}}_m,{\boldsymbol{d}}_n)\right]_{1\leq m,n\leq 2L}\in{\mathbb{C}}^{2L\times 2L}\tag{37}\] then yields an approximation for the expansion coefficients \((a_{m,n})_{-L\leq m,n\leq L-1}\) as defined in 11 . If we consider reconstruction from full far field data, then \(u_B^\infty(\widehat{{\boldsymbol{x}}}_m,{\boldsymbol{d}}_n)\) in 37 is replaced by \(u^\infty(\widehat{{\boldsymbol{x}}}_m,{\boldsymbol{d}}_n)\) from 6 .
Example 3. As the first example, we consider the piecewise constant contrast \[q=\chi_{B_{r_1}({\boldsymbol{c}}_1)}-0.25\chi_{B_{r_2}({\boldsymbol{c}}_2)}+0.5\chi_{B_{r_3}({\boldsymbol{c}}_3)}\,,\] where \({\boldsymbol{c}}_1=(-0.35,0.4)^\top\), \({\boldsymbol{c}}_2=(-0.1,-0.45)^\top\), \({\boldsymbol{c}}_3=(0.45,0.1)^\top\), \(r_1=r_2=0.3\) and \(r_3=0.2\). We assume the prior knowledge that the ROI \(B_1({\boldsymbol{0}})\) encloses the support of the contrast and consider the wave number \(\kappa = 30\). The exact contrast restricted to the ROI is shown in Figure 5 (left).
Figure 5: Exact contrast in Example 3 (left) and in Example 4 (right)..
In this setting, the Born far field data can be expressed analytically: by separately considering the expansion coefficients of the Born far field operator deduced in Example 1 for the three discoidal inclusions, it follows that \[u_B^\infty(\widehat{{\boldsymbol{x}}}_m,{\boldsymbol{d}}_n) \,=\, \sum_{i=1}^3 \sum_{l=1}^\infty 2\pi^2(\kappa r_i)^2 \left( J_{l-1}^2(\kappa r_i)-J_{l-2}(\kappa r_i)J_l(\kappa r_i)\right){\boldsymbol{e}}^{{\boldsymbol{c}}_i}_l(\widehat{{\boldsymbol{x}}}_m-{\boldsymbol{d}}_n)\] for \(m,n=1,\ldots,2L\). We choose \(L=125\) and truncate the series at \(250\) terms. For expanding the contrast, we use \(N_r=250\) Gauss–Legendre nodes in the radial direction and \(N_\varphi=250\) equiangular nodes in the angular direction.
Figure 6: Example 3. Reconstructed contrast from the exact Born far field data for different truncation indices. Top left: \(N=5\) (too small). Top right: \(N=10\) (too small). Middle left: \(N=15\) (too small). Middle right: \(N=20\) (slightly too small). Bottom left: \(N=29\) (optimal). Bottom right: \(N=31\) (too large)..
We first study the dependence of the reconstruction quality in the noise-free case on the truncation index \(N \in \{1,\dots,35\}\) without employing truncated SVD for further regularization. Selected reconstructions are shown in Figure 6, and the corresponding relative \(L^2\) reconstruction error over the ROI is presented in Figure 7 (left).
Figure 7: Example 3. Left: Relative \(L^2\) reconstruction error \(\varepsilon_{\mathrm{rel}}\) as a function of the truncation index \(N\) with exact Born far field data. The optimal choice \(N=29\) is marked by a vertical line. Right: The best and worst relative \(L^2\) reconstruction errors \(\varepsilon_{\mathrm{rel}}\) as functions of the absolute noise level \(p\) over \(20\) runs with different realizations of noise..
The error curve illustrates the trade-off as described in Remark 5: choosing the truncation index \(N\) too small leads to loss of relevant information, whereas choosing it too large results in dominance of the orthonormalization error introduced by the Gram–Schmidt process. Interestingly, the region in which the contrast is reconstructed accurately expands gradually outward from the center when increasing \(N\) until the minimal relative reconstruction error is achieved for \(N=29\approx\kappa R\). Figures 6 and 7 reveal that our method cannot achieve arbitrarily accurate reconstructions since higher angular and radial modes are excluded from the computation by construction; the relative \(L^2\) error remains above \(20\%\) for all truncation indices. In particular, discontinuities in the contrast, i.e., the jumps at the boundaries of the discs, cannot be reconstructed exactly. However, the reconstructions for \(N = 20\) and \(N = 29\) in Figure 6 can be considered visually satisfactory. According to the reconstruction for \(N=31\) in Figure 6 and the error plot in Figure 7 (left), even a small increase in the truncation index \(N\) beyond its optimal value leads to poor reconstructions. Hence, without further regularization, the quality of the reconstruction is strongly influenced by the quality of the a priori known ROI.
Figure 8: Example 3. Worst reconstructed contrasts in the sense of relative \(L^2\) error from noisy Born far field data over \(20\) runs at the noise levels \(p=20\) (left) and \(p=80\) (right)..
Next, we fix the truncation index \(N=30\) and study the quality of our reconstructions when the observed far field data are contaminated by additional noise. To this end, we add to each element of the exact data matrix in 36 complex noise whose real and imaginary parts are independently drawn from a zero-mean uniform distribution, with the standard deviation scaled a posteriori so that the Frobenius norm of the added noise matrix is \(p\in\{0,5,\ldots,95,100\}\) per cent of the Frobenius norm of the exact data matrix. In what follows, we refer to this noise model by simply saying that the data contain \(p\)% of noise. We use the truncated SVD together with the Morozov discrepancy principle with respect to the Euclidean norm as an additional regularization strategy, with the target vector carrying the noisy truncated coefficients of the Born far field operator (cf. 35 and 37 ) and the employed noise level scaled appropriately to account for the amount of noise no longer present in the data after the truncation. We generate \(20\) independent noise realizations for each \(p\) and plot the relative \(L^2\) errors of the resulting best and worst reconstructions in Figure 7 (right). The related worst-case reconstructions for \(p=20\) and \(p=80\) are presented in Figure 8. In the shown worst cases, \(45\%\) and \(24\%\) of the singular components are taken into account for \(p = 20\) and \(p=80\), respectively. Interestingly, our method turns out to be very robust to this form of additive noise since it is distributed across both high and low modes in the data, and thus a large fraction of the noise is filtered out by the introduction of the truncation index \(N\). In the next example, we will see that this does not apply when the model discrepancy originates from multiple scattering effects included in full far field data. \(\lozenge\)
Example 4. Next we study a smooth contrast of a similar geometric structure as in the previous example, with the aim to compare reconstructions obtained using Born far field data and full far field data as the input for our reconstruction algorithm. We examine how reconstructions from full far field data deteriorate with increasing multiple scattering effects, i.e., as the Born far field data become an increasingly inaccurate approximation of the full far field data observed in practice. This is achieved by keeping the contrast \(q\) fixed as shown in Figure 5 (right) while gradually increasing the wavenumber from \(\kappa=1\) to \(\kappa=56\), which also leads to an increase in the truncation index from \(N=1\) to \(N=56\). When considering full far field data, we accompany the data truncation with truncated SVD with the noise level for the Morozov principle chosen (unrealistically) to be the Euclidean norm of the discrepancy in the data vector in 35 between Born and full far field data; see 37 and the comment succeeding it. For evaluating both Born and full far field data in \(2L=250\) equiangular observation and illumination directions, we use the fast Lippmann–Schwinger solver proposed by Vainikko [27]. For expanding the contrast, we choose \(N_r=N_\varphi=310\). The corresponding relative \(L^2\) reconstruction errors along with the relative Frobenius-based data error (cf. 36 ) induced by approximating the Born far field data with the full far field data are shown in Figure 9. The reconstructed contrasts for \(\kappa=11\) and for \(\kappa=46\) are illustrated in Figure 10. According to Figures 9 and 10, the best reconstruction using full far field data is obtained at \(\kappa = 11\), with practically no visual difference in quality to the corresponding reconstruction based on Born far field data. Although multiple scattering effects dominate at \(\kappa = 46\), amounting to about 80% of the magnitude of the linearized component, the reconstruction from full far field data still provides a clear picture of the positions and sizes of the three components of the scatterer, with the exterior shapes of two of the components still recognizable. However, the reconstruction does not capture the dynamic range of the target contrast, and it includes holes inside the scattering components.
Figure 9: Example 4. Relative \(L^2\) reconstruction error \(\varepsilon_{\mathrm{rel}}\) as a function of the wave number \(\kappa\) with the exact Born far field data and the exact full far field data as the inputs for the proposed method. The relative Frobenius-based approximation error induced by replacing the Born far field data with the full far field data is also shown for reference (cf. 37 )..
Figure 10: Example 4. Reconstructed contrast from the exact Born far field data (left) and the exact full far field data (right). Top: \(\kappa=11\) (the best reconstruction from the full far field data, cf. Figure 9). Bottom: \(\kappa=46\) (multiple scattering effects dominate linearized data)..
For \(\kappa < 11\), accurate reconstructions cannot be obtained even from Born far field data, as \(N\) is chosen so small that the amount of data is insufficient for reconstructing high spatial frequencies in the contrast \(q\). For higher values of \(\kappa\), enough information is available, with the relative error curve for Born far field data plateauing at around \(7\%\) in Figure 10.
According to [13], the expansion coefficients of full far field data have the same essential support as those of Born far field data (cf. 32 ). Hence, it can be argued that the truncation of the infinite-dimensional systems at index \(N\) does not remove retrievable information, which partially explains why structural information about the contrast geometry can be reconstructed reasonably even when multiple scattering effects are dominant. On the other hand, unlike in the case of artificially generated noise in the previous example, no noise filtering of the model error can be expected via truncation of the data matrix. \(\lozenge\)
Example 5. In this final example we compare our method with the low-rank method proposed in [12] and MATLAB’s built-in NUFFT as described in [26].
The method presented in [12] follows a similar strategy to ours: based on the Born approximation, a simple invertible system matrix is derived by expanding both the contrast and the observed far field data in terms of suitably orthonormal systems for \(L^2(B_1({\boldsymbol{0}}))\), namely the prolate spheroidal wave functions. Obtaining a diagonal system matrix, instead of a block-wise triangular one as in our case, comes at the cost of additional modeling error due to an unnatural sampling structure for the far field data.
Regarding the NUFFT, we refer to the state of the art implementation fiNUFFT by the Flatiron Institute (cf. [28], [29]), which outperforms MATLAB’s built-in option in terms of computational time for nonuniform sampling points. Since such efficiency aspect is not relevant for our numerical examples, we nevertheless employed MATLAB’s built-in implementation in the presented examples.
To generate the data, we use the Matlab file “ship2D.m” from the toolbox IPscatt (see [30]), which was also used for the numerical tests in [12]. The real and imaginary parts of the exact contrast are shown in Figure 11, where we assume the prior information that the support of the contrast is contained in \(B_1({\boldsymbol{0}})\). We simulate full far field data using the fast Lippmann–Schwinger solver [27] and \(2L=200\), and subsequently add \(20\%\) of noise to the data (as described in Example 3). We run the three methods for the wave numbers \(\kappa=30\) and \(\kappa=60\). This leads to the truncation indices \(N=30\) and \(N=60\) for our method, and in addition, we employ truncated SVD with the Morozov discrepancy principle as in Example 4, assuming unrealistically the knowledge of the truncated Born far field data for choosing the spectral cut-off. For our method and the low-rank method, we set \(N_r=N_\varphi=400\) for expanding the contrast. For NUFFT, we use a cartesian grid enclosing the ROI and consisting of \(100\) equidistant nodes in both dimensions. As the spectral cut-off in the low rank method, we use \(90\%\) of the prolate eigenvalue of largest magnitude as suggested in [12].
Figure 11: Example 5. Real part (left) and imaginary part (right) of the exact contrast function..
Figure 12: Example 5. Real part (left) and imaginary part (right) of reconstructed contrast function for wave number \(\kappa=30\) from full far field data with \(20\%\) of additive noise. Top: our proposed method. Middle: the low-rank method [12]. Bottom: MATLAB’s built-in NUFFT..
Figure 13: Example 5. Real part (left) and imaginary part (right) of reconstructed contrast function for wave number \(\kappa=60\) from full far field data with \(20\%\) additive noise. Top: our proposed method. Middle: the low-rank method [12]. Bottom: MATLAB’s built-in NUFFT..
The reconstructions for \(\kappa=30\) are shown in Figure 12 and those for \(\kappa=60\) in Figure 13. Our proposed method and the low-rank method from [12] produce reconstructions of comparable quality, whereas NUFFT yields a slightly blurred reconstructions of the contrast. Based on this single example, the performance of our method seems to be on par with other direct reconstruction algorithms that are based on the Born approximation, but drawing more precise conclusions is not possible without more detailed testing. \(\lozenge\)
Conclusions↩︎
Following the ideas in [3], [22] for EIT, we introduced a direct reconstruction method for inverse medium Born scattering for the Helmholtz equation in two spatial dimensions. Choosing appropriate basis for representing the far field operator and the contrast function, the proposed method reduces the inverse problem to solving decoupled triangular systems that correspond to different angular frequencies in the contrast. The bases for representing the far field operator and the angular behavior of the contrast function are of Fourier type and can be given explicitly, but introducing the needed radial bases for the contrast requires numerically orthogonalizing certain products of Bessel functions, which adds an additional unstable step to the algorithm. On the positive side, this orthogonalization can be performed offline before the data is available and can also be stabilized by truncating the to-be-inverted system based on ideas in [13]. Due to the achieved angular decoupling and the triangular structure of the subsystems, the proposed method allows an efficient numerical implementation, as well as an explicit recursive reconstruction formula (Theorem 4) if instability issues are not considered.
The presented numerical experiments demonstrate that our method produces reconstruction of good quality from (noisy) Born far field data, and it can also be applied to full nonlinear far field data to obtain reconstructions ranging in quality from good to reasonable depending on the extent of multiple scattering effects. According to our numerical experiments, the proposed method compares favorably to other algorithms designed for solving the inverse medium Born scattering problem, in particular, producing reconstructions comparable to those by the low-rank method introduced recently in [12].
A natural direction for future work is to derive explicit representations for the radial basis functions; this would also be a first step toward extending the stability results of [22] from EIT to inverse medium scattering. Another promising avenue is to generalize our reconstruction method to three dimensions, building on the ideas in [4].
Acknowledgments↩︎
This work was supported by the Research Council of Finland (Flagship of Advanced Mathematics for Sensing Imaging and Modelling grant 359181).
References↩︎
IPscatt – A
MATLAB toolbox for the inverse medium problem in scattering,” ACM Trans. Math. Software, vol. 45, no. 4, pp. Art. 45, 20, 2019, doi: 10.1145/3328525.Department of Mathematics and System Analysis, Aalto University, 00076 Helsinki, Finland (
nuutti.hyvonen@aalto.fi, lisa.schatzle@aalto.fi).↩︎