Comparisons between Isotropic and Anisotropic TV Regularizations in Inverse Acoustic Scattering

Ping Liu

School of Mathematics and Statistics, Shandong Normal University, Jinan, China.

DOI: 10.4236/eng.2023.152009                

Abstract

This article compares the isotropic and anisotropic TV regularizations used in inverse acoustic scattering. It is observed that compared with the traditional Tikhonov regularization, isotropic and anisotropic TV regularizations perform better in the sense of edge preserving. While anisotropic TV regularization will cause distortions along axes. To minimize the energy function with isotropic and anisotropic regularization terms, we use split Bregman scheme. We do several 2D numerical experiments to validate the above arguments.

Keywords

Inverse Acoustic Scattering Problem, Regularization, Isotropic TV, Anisotropic TV, Ill-Posedness

1. Introduction

Inverse acoustic scattering problem is to get the refractive index from measurements of the scattered field or its far field pattern data. It has wide applications in many aspects such as radar, sonar, geophysical exploration, medical CT imaging and nondestructive testing [1] .

Inverse acoustic scattering problem is typically non-linear and ill-posedness. To be precise, the small disturbance of the measurement data could cause severe error in the inversion results. For nonlinearity, the most common way is to adopt linear or optimization method. The regularization method can approximate the ill-posed problem to the well-posed problem to estimate the refractive index distribution. Therefore, iterative regularization method is a typical way to deal with such problems. For example, simplified Newton method [2] , modified gradient method [3] , quasi-Newton method [4] , Gauss-Newton method [5] [6] . In addition, other research results on this issue include: Dual-space method [7] , linear sampling method (LSM) [8] [9] structure recognition function $g_z$ , which is the solution of the far-field equation $F_{g_z}={\Phi }_{\infty }\left(\cdot ,z\right)$ , factorization method

[10] is solved by ${\left(F\ast F\right)}^{\frac{1}{4}}{g}_z={\Phi }_{\infty }\left(\cdot \mathrm{,}z\right)$ in place of the far field equations, multiple

signal classification (MUSIC) [11] [12] constructs a non-iterative solution function, level set method [13] is to determine the regional boundary method.

In this article, we use the isotropic and anisotropic TV regularizations to overcome the ill-posed of inverse acoustic scattering, and the split Bregman algorithm to overcome the non-differentiability of the reconstruction problem and accelerate the inversion process. Finally, the numerical experiment is carried out to verify.

2. Forward Problem and Inverse Scattering Problem

For a given time harmonic incident wave $u^i={\text{e}}^{i\omega t}$ , the mathematical model of its propagation in inhomogeneous medium can be given by the following Helmholtz equation

$\Delta u^t+k^2{n}^2{u}^t=0\text{in}{R}^d\mathrm{.}$ (1)

Here, $k>0$ is the wavenumber, the refractive index function n

$n=\{\begin{array}{ll}1\hfill & R^d\backslash D\mathrm{,}\hfill \\ \text{otherwise}\hfill & D\mathrm{,}\hfill \end{array}$

$u^t=u^i+u^s$ is the total field where $u^s$ represents the scatter field. To ensure the uniqueness of the solution, we require Sommerfeld radiation condition

$\underset{\left|x\right|\to \infty }{lim}{\left|x\right|}^{\frac{d-1}{2}}\left(\frac{\partial }{\partial \left|x\right|}-ik\right)u^s\left(x\right)=0$ (2)

uniformly in all directions $\hat{x}=\frac{x}{\left|x\right|}$ . The forward problem is to derive a solution using the Equations ((1), (2)). The forward problem is governed by the following Lippmann-Schwinger equations [1] .

$u^t\left(x\right)=u^i\left(x\right)+k^2{\displaystyle {\int }_D}\Phi \left(x\mathrm{,}y\right)\left(n^2-1\right)u^t\left(y\right)\text{d}y\mathrm{,}y\in R^d\backslash D\mathrm{.}$ (3)

Here, $\Phi$ is the fundamental solution of the Helmholtz equation i.e.

$\Phi \left(x\mathrm{,}y\right)=\{\begin{array}{l}\frac{i}{4}{H}_0^1\left(k\left|x-y\right|\right)\mathrm{,}\text{if}x\mathrm{,}y\in R^{2}x\ne y\mathrm{,}\\ \frac{1}{4\pi }\frac{{\text{e}}^{ik\left|x-y\right|}}{\left|x-y\right|}\mathrm{,\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{if}x\mathrm{,}y\in R^{3}x\ne y\mathrm{,}\end{array}$

where $H_0^{\left(1\right)}$ denotes the zero-th order Hankel function of the first kind. Note that for the sake of solving, we assume that the contrast $q=n^2-1$ is satisfied $\Re \left(q\right)>-1$ and $\Im \left(q\right)\ge 0$ in this article and $\text{D}\in R^d\left(d=\mathrm{2,3}\right)$ is a bound domain where contain the support of the contrast q. For simplicity, we assume that the volume potential V is given by

$\left(\text{V}u\right)\left(x\right)=k^2{\displaystyle {\int }_D}\Phi \left(x\mathrm{,}y\right)u\left(y\right)\text{d}y\mathrm{,}$ (4)

so (3) can be described as the following equation

$u^t=u^i+\text{V}\left(q\cdot u^t\right)\mathrm{.}$ (5)

Hence, ${u^t|}_{\Gamma }={\left(Id-\text{V}\left(q\cdot \right)\right)}^{-1}{u}^i$ and ${u^s|}_{\Gamma }=\text{V}\left(q\left(Id-\text{V}\right){\left(q\cdot \right)}^{-1}{u}^i\right)$ .

The forward operator F is defined via:

$\text{F}\left(q\right)\mathrm{<mo>\; :\; </mo>}=\text{V}\left(q{\left(Id-\text{V}\left(q\cdot \right)\right)}^{-1}{u}^i\right)\mathrm{.}$ (6)

The problem of Inverse acoustic scattering can be recast contrast q by the measured data of (6). But, it is well known that (6) is nonlinear about q, which may cause additional complexity in reconstruction methods, so we consider linearizing the above equation. Given the known fixed value q₀, we linearize the forward operator F with Taylor polynomial at q₀ as the following equation

$\text{F}\left(q\right)\approx \text{F}\left(q_0\right)+\text{<msup> F <mo> ′ </mo> </msup>}\left(q_0\right)\left[q-q_0\right]\mathrm{,}$ (7)

where $\text{F}\left(q_0\right)$ is measured scattered fields at q₀, $\text{<msup> F <mo> ′ </mo> </msup>}\left(q_0\right)$ is Fréhet derivative of $\text{F}\left(q_0\right)$ . For the purpose of computerized reconstruction, we discretize the domain D as $D={\displaystyle {\cup }_{i=1}^{N_D}}{T}_i$ , where $T_i$ is a rectangle element, then ${q|}_{T_i}$ can be considered as a constant $q_i$ , q can be approximated by the vector $q={\left(q_1\mathrm{,}{q}_2\mathrm{,}\cdots \mathrm{,}{q}_{N_D}\right)}^{\text{T}}$ . We consider minimization of the discretized functional

$\delta {\text{F}}_{meas}\approx G\delta q\mathrm{.}$ (8)

Here, $G\in R^{N_t{N}_r\times N_D}$ is the jacobian matrix of $\text{F}\left(q\right)$ , $\delta q=q-q_0\in R^{N_D}$ , $\delta {\text{F}}_{meas}=\text{F}\left(q\right)-\text{F}\left(q_0\right)\in R^{N_t{N}_r}$ ,where $N_t$ , $N_r$ represents the number of transmitters, receivers on $\Gamma$ . Noting that the measurement data $\delta {\text{F}}_{meas}$ may have some error due to perturbed measurements or errors, Intuitively, the problem of inverse acoustic scattering can be recast into the following least square error problem to better approximations of the exact solution.

$\delta q^{\ast }=\arg \underset{\delta q}{\min }\frac{1}{2}{\Vert G\delta q-\delta F_{meas}\Vert }_{{\mathcal{l}}_2}^2\mathrm{,}$ (9)

where $\delta q^{\ast }\in C^{N_D}$ is the searched-for exact contrast, $G\in C^{N_t{N}_r\times N_D}$ , $F_{meas}\in C^{N_t{N}_r}$ . However, the solution (9) is a complex valued function, so it is difficult to use regularization methods. Thus, we convert a complex-valued functional (9) to a real-valued functional (10) as following

$\delta q^{\ast }=\arg \underset{\delta q}{\min }\frac{1}{2}{\Vert G\delta q-\delta F_{meas}\Vert }_{{\mathcal{l}}_2}^2\mathrm{.}$ (10)

Here, $G=\left(\begin{array}{cc}\Re G& -\Im G\\ \Im G& \Re G\end{array}\right)\in R^{2N_t{N}_r\times 2N_D}$ , $\delta q=\left(\begin{array}{c}\Re \delta q\\ \Im \delta q\end{array}\right)\in R^{2N_D}$ ,

$F_{meas}^{\delta }=\left(\begin{array}{c}\Re F_{meas}^{\delta }\\ \Im F_{meas}^{\delta }\end{array}\right)\in R^{2N_t{N}_r}$ . However, (10) solely minimizing the discrepancy

causes numerical instabilities. Thus, we incorporate a-priori information about the solution to avoid instable results. The most commonly used regularization method is Tikhonov regularization (TR), which is to solve

$\delta q^{\ast }=\arg \underset{\delta q}{\min }\frac{1}{2}{\Vert G\delta q-{\text{F}}_{meas}^{\delta }\Vert }_{{\mathcal{l}}_2}^2+\lambda {\Vert \delta q\Vert }_{{\mathcal{l}}_1}$ (11)

with a positive regularization parameter $\lambda$ . However, it has an excessively smooth effect on the solution, which will blur the edge of the reconstructed image. To preserve the shape, we use the isotropic TV

$\delta q^{\ast }=\arg \underset{\delta q}{\min }\frac{1}{2}{\Vert G\delta q-{\text{F}}_{meas}^{\delta }\Vert }_{{\mathcal{l}}_2}^2+\lambda {\Vert \sqrt{{\left|D_x\delta q\right|}^2+{\left|D_y\delta q\right|}^2}\Vert }_{{\mathcal{l}}_1}$ (12)

and anisotropic TV method

$\delta q^{\ast }=\arg \underset{\delta q}{\min }\frac{1}{2}{\Vert G\delta q-F_{meas}^{\delta }\Vert }_{{\mathcal{l}}_2}^2+\lambda {\Vert D_x\delta q\Vert }_{{\mathcal{l}}_1}+\lambda {\Vert D_y\delta q\Vert }_{{\mathcal{l}}_1}$ (13)

with respectively the first-order difference operators along the x and y directions of $D_x,D_y\in R^{2N_D\times 2N_D}$ . However, (12) is non-differentiability. So, we use the Split Bregman algorithm to solve the non-differentiability problem, i.e. by introducing new variable $u_x\mathrm{,}{u}_y\in R^{2N}$ , (12) needs to be converted to minimization problem:

$\{\begin{array}{l}\delta q^{\ast }=\arg \underset{\delta q}{\min }\frac{1}{2}{\Vert G\delta q-F_{meas}^{\delta }\Vert }_{{\mathcal{l}}_2}^2+\lambda {\Vert \sqrt{u_x^2+u_y^2}\Vert }_{{\mathcal{l}}_1}\\ \text{s}\text{.t}\mathrm{.}{D}_x\delta q=u_x{D}_y\delta q=u_y\mathrm{.}\end{array}$ (14)

The minimization of subproblems in (14) can be separated into the minimization of three simple subproblems $dq$ and $u_x\mathrm{,}{u}_y$ . Therefore, the following split bregman iteration stages are proposed:

$\{\begin{array}{l}\delta q^{k+1}=\arg \underset{\delta q}{\min }\frac{1}{2}{\Vert G\delta q-F_{meas}^{\delta }\Vert }_{{\mathcal{l}}_2}^2+\frac{\alpha }{2}{\Vert D_x\delta q-u_x-b_x\Vert }_{{\mathcal{l}}_2}^2+\frac{\alpha }{2}{\Vert D_y\delta q-u_y-b_y\Vert }_{{\mathcal{l}}_2}^2\\ \left(u_x^{k+1}\mathrm{,}{u}_y^{k+1}\right)=\arg \underset{u}{\min }\lambda {\Vert \sqrt{u_x^2+u_y^2}\Vert }_{{\mathcal{l}}_1}+\frac{\alpha }{2}{\Vert D_x\delta q-u_x-b_x\Vert }_{{\mathcal{l}}_2}^2+\frac{\alpha }{2}{\Vert D_y\delta q-u_y-b_y\Vert }_{{\mathcal{l}}_2}^2\\ b_x^{k+1}=b_x^k+\alpha \left(D_x\delta q-u_x\right)\\ b_y^{k+1}=b_y^k+\alpha \left(D_y\delta q-u_y\right)\mathrm{.}\end{array}$ (15)

The (13) is non-differentiability. Similarly, (13) needs to be converted to

$\{\begin{array}{l}\delta q^{\ast }=\arg \underset{\delta q}{\min }\frac{1}{2}{\Vert G\delta q-F_{meas}^{\delta }\Vert }_{{\mathcal{l}}_2}^2+\lambda {\Vert D_x\delta q\Vert }_{{\mathcal{l}}_1}+\lambda {\Vert D_y\delta q\Vert }_{{\mathcal{l}}_1}\\ \text{s}\text{.t}\mathrm{.}{D}_x\delta q=u_x{D}_y\delta q=u_y\end{array}$

$\begin{array}{c}\delta q^{k+1}=\arg \underset{\delta q}{\min }\{\frac{1}{2}{\Vert G\delta q-F_{meas}^{\delta }\Vert }_{{\mathcal{l}}_2}^2+\frac{\alpha }{2}{\Vert D_x\delta q-u_x-b_x\Vert }_{{\mathcal{l}}_2}^2\\ +\frac{\alpha }{2}{\Vert D_y\delta q-u_y-b_y\Vert }_{{\mathcal{l}}_2}^2\}\mathrm{<mo>\; ;\; </mo>}\end{array}$ (16)

$u_x^{k+1}=\arg \underset{u}{\min }\lambda {\Vert u_x+\frac{\alpha }{2}{\Vert D_x\delta q-u_x-b_x\Vert }_{{\mathcal{l}}_2}^2\Vert }_{{\mathcal{l}}_1}\mathrm{<mo>\; ;\; </mo>}$ (17)

$u_y^{k+1}=\arg \underset{u}{\min }\left\{\lambda {\Vert u_y\Vert }_{{\mathcal{l}}_1}+\frac{\alpha }{2}{\Vert D_y\delta q-u_y-b_y\Vert }_{{\mathcal{l}}_2}^2\right\}\mathrm{<mo>\; ;\; </mo>}$ (18)

$b_x^{k+1}=b_x^k+\alpha \left(D_x\delta q-u_x\right)\mathrm{<mo>\; ;\; </mo>}$ (19)

$b_y^{k+1}=b_y^k+\alpha \left(D_y\delta q-u_y\right)\mathrm{.}$ (20)

The first-order optimality condition of (16) is

$\left(G^{\text{T}}G+\alpha D_x^{\text{T}}{D}_x+\alpha D_y^{\text{T}}{D}_y\right)\delta q^{k+1}=\alpha \left(u_x^k-b_x^k\right)+\alpha \left(u_y^k-b_y^k\right)+G^{\text{T}}{F}_{meas}^{\delta }$ (21)

The minimization of ((17), (18)) can be found by the soft threshold formula

$\begin{array}{l}{u}_x^{k+1}\left[k\right]=h_{\frac{\lambda }{\rho }}\left({\left(D\underset{\_}{\delta q}\right)}^{k+1}\left[k\right]+\frac{b_x^{k+1}\left[k\right]}{\rho }\right)k=\mathrm{1,2,}\cdots \mathrm{,2}{N}_D\mathrm{.}\\ u_y^{k+1}\left[k\right]=h_{\frac{\lambda }{\rho }}\left({\left(D\underset{\_}{\delta q}\right)}^{k+1}\left[k\right]+\frac{b_y^{k+1}\left[k\right]}{\rho }\right)k=\mathrm{1,2,}\cdots \mathrm{,2}{N}_D\mathrm{.}\end{array}$ (22)

Here, $b_x^{k+1}\left[k\right]\mathrm{,}{b}_y^{k+1}\left[k\right]$ and $u_x^{k+1}\left[k\right]\mathrm{,}{u}_y^{k+1}\left[k\right]$ are respectively the k-th element of $b_x\mathrm{,}{b}_y$ and $u_x\mathrm{,}{u}_y$ , $h_g(\cdot )$ is the soft-thresholding operator defined as

$h_g(\cdot )=\{\begin{array}{l}\cdot -g\mathrm{sgn}(\cdot )\mathrm{,}\left|\cdot \right|>g\\ \mathrm{0,}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{otherwise}\mathrm{,}\end{array}$ (23)

where $\mathrm{sgn}(\cdot )$ is the sign function. We summarize the above procedures as the following algorithm in the form of pseudocode.

3. Numerical Experiments

To validate the proposed method, we demonstrate the reconstruction quality of isotropic TV, anisotropic TV, and Tikhonov regularization methods using some synthetic data. Firstly, we solve the forward problem (3) using the Integral equation to obtain the near field of $u^s$ as measured data.

All tests were conducted on an Intel Core i7 3.40 GHz CPU and 16 GB of RAM. To avoid “inverse crime,” the discrete computational region of the grid in the forward problem is finer than that in the inverse problem. We show the test set for synthetic data in Figures 1(a)-(f). Here, we choose the domain D as $\left(\mathrm{0,1}\right)\times \left(\mathrm{0,1}\right)$ , wavenumber $k=200$ and 35 transmitters, receivers on $\Gamma$ . we add 1% random noise to measured data for test the stability of method.

In Figure 1, we depict the test set of synthetic data. We apply Tikhonov

regularization, isotropic TV and anisotropic TV methods, respectively, to carry out the numerical simulation. The parameters of the three methods are set to be optimal empirically. It can be seen from Figure 2 that the location of isotropic TV and anisotropic TV reconstruction is more accurate, the shape is closer, and the reconstruction effect is better than the usual Tikhonov regularization (TR) reconstruction.

Figure 1. The 2D numerical model.

Figure 2. The first row (a)-(f) shows that the inversion result of Tikhonov regularization method; the second row (g)-(l) shows the that the inversion result of isotropic TV method; the Third row (m)-(r) shows the that the inversion result of anisotropic TV method.

There are several observations from the two figures. Firstly, all three regularization methods can well capture the main feature of the inner object, including position and shape. From Figures 2(a)-(f) shows that using SB to solve the Tikhonov regularization problem; from Figures 2(g)-(l) show that the images reconstructed by isotropic TV have an obvious ladder effect; from Figures 2(m)-(r) shows that using split Bregman to solve the anisotropic TV regularization problem can obtain an accurate image. The edges of the reconstructed images using the anisotropic TV distort along the coordinate axes.

The numerical experiment shows that the effect of method isotropic and anisotropic TV regularization in inverse acoustic scattering reconstructions is better sharpen the edges and is more robust against data noise than that of method Tikhonov regularization. But, in the anisotropic numerical simulation, the inversion image is distorted along the edge of the coordinate axis alignment, which indicates that there is a large geometric distortion. In order to simplify the operation process and save the running time, we apply the split Bregman algorithm to inverse acoustic scattering problem.

4. Conclusion and Future Work

In this article, we use isotropic TV and anisotropic TV regularization by the split Bergeman method to solve the inverse acoustic scattering problem. These two methods were compared with the Tikhonov method. The simulation results of the measured data show that the isotropic TV regularization can cause a staircase effect; the anisotropic regularization can cause geometric distortion along the coordinate axis. However, these two methods can well preserve boundaries more than Tikhonov regularization. In future work, we will focus on a method that can avoid the distortions along the coordinate axis and does not depend on the selection of regularization parameters.

Acknowledgements

The author would like to thank editor and referees for their valuable advice for the improvement of this article.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Colton, D. and Kress, R. (1998) Inverse Acoustic and Electromagnetic Scattering Theory. Springer, Berlin. https://doi.org/10.1007/978-3-662-03537-5
[2]	Author (2002) A Simplified Newton Method in Scattering Theory: The TE Case. IEEE Antennas and Propagation Society International Symposium (IEEE Cat. No. 02CH37313), 2, 102-105.
[3]	Colton, D. and Kress, R. (1983) Integral Equation Methods in Scattering Theory. John Willey and Sons, Inc., New York.
[4]	Gutman, S. and Klibanov, M. (1993) Regularized Quasi-Newton Method for Inverse Scattering Problems. Mathematical and Computer Modelling of Dynamical Systems, 18, 5-31. https://doi.org/10.1016/0895-7177(93)90076-B
[5]	Hohage, T. and Langer, S. (2010) Acceleration Techniques for Regularized Newton Methods Applied to Electromagnetic Inverse Medium Scattering Problems. Inverse Problems, 26, 15-28. https://doi.org/10.1088/0266-5611/26/7/074011
[6]	Zaeytijd, J.D., Franchois, A., Eyraud, C. and Geffrin, J.-M. (2007) Full-Wave Three-Dimensional Microwave Imaging with a Regularized Gauss-Newton Method—Theory and Experiment. IEEE Transactions on Antennas and Propagation, 2 November 2007, 3279-3292. https://doi.org/10.1109/TAP.2007.908824
[7]	Colton, D. and Monk, P. (1985) A Novel Method for Solving the Inverse Scattering Problem for Time-Harmonic Waves in the Resonance Region. SIAM Journal on Applied Mathematics, 45, 1039C1053. https://doi.org/10.1137/0145064
[8]	Ito, K., Jin, B. and Zou, J. (2013) A Direct Sampling Method for Inverse Electromagnetic Medium Scattering. Inverse Problems, 29, 095018. https://doi.org/10.1088/0266-5611/29/9/095018
[9]	Cakoni, F., Colton, D. and Monk, P. (2011) The Linear Sampling Method in Inverse Electromagnetic Scattering. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898719406
[10]	Kirsch, A. and Grinberg, N.I. (2008) The Factorization Method for Inverse Problems. Oxford University Press, Oxford. https://doi.org/10.1093/acprof:oso/9780199213535.001.0001
[11]	Kirsch, A. (2002) The MUSIC Algorithm and the Factorization Method in Inverse Scattering Theory for Inhomogeneous Media. Inverse Problems, 18, 1025-1040. https://doi.org/10.1088/0266-5611/18/4/306
[12]	Ammari, H., Iakovleva, E., Lesselier, D. and Perrusson, G. (2007) MUSIC-Type Electromagnetic Imaging of a Collection of Small Three-Dimensional Inclusions. SIAM Journal on Scientific Computing, 29, 674-709. https://doi.org/10.1137/050640655
[13]	Dorn, O. and Lesselier, D. (2009) Level Set Methods for Inverse Scattering: Some Recent Developments. Inverse Problems, 25, 125001. https://doi.org/10.1088/0266-5611/25/12/125001
[14]	Goldstein, T, and Osher, S. (2009) The Split Bregman Method for L1 Regularized Problems. SIAM Journal on Imaging Sciences, 2, 323-343. https://doi.org/10.1137/080725891