Contraction Integral Equation for Three-Dimensional Electromagnetic Inverse Scattering Problems

Inverse scattering problems (ISPs) stand at the center of many important imaging applications, such as geophysical explorations, industrial non-destructive testing, bio-medical imaging, etc. Recently, a new type of contraction integral equation for inversion (CIE-I) has been proposed to tackle the two-dimensional electromagnetic ISPs, in which the usually employed Lippmann–Schwinger integral equation (LSIE) is transformed into a new form with a modified medium contrast via a contraction mapping. With the CIE-I, the multiple scattering effects, i.e., the physical reason for the nonlinearity in the ISPs, is substantially suppressed in estimating the modified contrast, without compromising physical modeling. In this paper, we firstly propose to implement this new CIE-I for the three-dimensional ISPs. With the help of the FFT type twofold subspace-based optimization method (TSOM), when handling the highly nonlinear problems with strong scatterers, those with higher contrast and/or larger dimensions (in terms of wavelengths), the performance of the inversions with CIE-I is much better than the ones with the LSIE, wherein inversions usually converge to local minima that may be far away from the solution. In addition, when handling the moderate scatterers (those the LSIE modeling can still handle), the convergence speed of the proposed method with CIE-I is much faster than the one with the LSIE. Secondly, we propose to relax the contraction mapping condition, i.e., different contraction mappings are used in updating contrast sources and contrast, and we find that the convergence can be further accelerated. Several numerical tests illustrate the aforementioned interests.


Introduction
Inverse scattering problems (ISPs) in electromagnetics and acoustics are of great interest in industries due to important imaging applications in various areas, such as geophysical survey, non-destructive testing, ground-penetrating radar, bio-medical imaging, etc., as solving the ISPs provides rich information about the unknown targets, such as locations, shapes, and material distributions within some structures [1][2][3][4].For instance, microwave imaging has been used to inspect the abnormalities in human bodies like bleeding in the head and tumours in breasts.On the other hand, they are also quite important due to the representative difficulties in solving a large group of inverse problems concerning waves and fields, and thus researchers in mathematics, physics, and engineering societies have devoted great efforts to improving efficiencies and accuracies of the numerical solvers [5].As shown in Figure 1, solving ISPs is to determine the unknown scatterers within a domain, when the domain is illuminated by several different incidences and we can measure the scattered fields outside the domain (usually contaminated by noise) for each incidence.It is well known that the main difficulties in such ISPs are their two intrinsic properties, i.e., ill-posedness and nonlinearity [1].As most of the practical problems are three-dimensional (3-D) ones, which require much more computational resources than those in two-dimensional cases, they are often more difficult to solve due to the data deficiency (measurement aperture usually covers a small solid angle), and therefore also stronger ill-posedness and nonlinearity.Great efforts have been paid to tackle these demanding problems, for instance [6][7][8][9].
Lippmann-Schwinger Eq.Methods for solving ISPs can be catalogued into two types of optimization approaches: deterministic and stochastic.The deterministic type of the inversion methods has been developed for decades, such as the contrast source inversion (CSI) method [6,10,11], the Born iterative method and distorted Born iterative method [12,13], the level set method [14,15], the subspace-based optimization method (SOM) [8,[16][17][18], etc.The second type, the stochastic type of the inversion methods, usually employs a group of initial guesses and uses the stochastic optimization scheme to minimize the objective function, such as the genetic algorithm and the evolutionary optimization.The stochastic type methods increase the possibilities of finding the global minimum rather than being trapped in a local minimum as the deterministic optimization techniques [19,20].Both techniques have also been applied to solve 3-D inverse scattering problems [6,[21][22][23][24].Other than the quantitative methods mentioned above, some qualitative methods, such as the linear sampling method [25] and other methods, like [26], are proposed to retrieve the geometric supports of the unknown scatterers.
To tackle the nonlinearity, the major difficulty in efficiently solving the ISPs, different types of integral equations have been proposed, such as in [27][28][29].In particular, motivated by the contraction integral equation (CIE) in solving the direct scattering problems with highly conductive background media [30], a new type of contraction integral equation for inversion (CIE-I) has been proposed [29] to tackle the two-dimensional highly nonlinear ISPs by transforming the usually employed Lippmann-Schwinger integral equation (LSIE) into a new form with a modified medium contrast via a contraction mapping.With the CIE-I, the global multiple scattering effects (MSE) are substantially suppressed in estimating the modified contrast in the CSI type methods.With the FFT type twofold SOM regularization scheme [8], the inversion solver with CIE-I is capable of effectively alleviating the nonlinearity of the two-dimensional ISPs, especially those highly nonlinear ones with strong scatterers with large contrast and/or large dimensions (in terms of wavelengths).In this paper, this new CIE-I will be implemented in tackling the 3-D ISPs.
In summary, the contributions of the paper include: 1.
The CIE-I is firstly implemented to tackle the computationally costly full 3-D ISPs, to address the highly nonlinear 3-D ISPs, and to accelerate the convergence of the inversions.

2.
A relaxed type of inversion scheme based on CIE-I is proposed, with different auxiliary parameters β (the parameter in CIE-I to control the portion of MSE in estimating the contrast) in updating the contrast sources and in updating the contrast.This means to further accelerate the convergence of the inversions.

3.
Several numerical tests are provided with details, for the sake of further algorithmic studies.
After the Introduction, the proposed 3-D inversion method is detailed in Section 2. In Section 3, three numerical examples are given to validate the proposed method.Finally, conclusions are drawn.

Inversion with CIE-I
In this paper, the domains of interest (DoI) are chosen to be rectangular cuboid in order to implement the conjugate gradient fast Fourier transform (CG-FFT) scheme and apply the Fourier basis in TSOM.For the convenience of reading, in this paper, we denote the one-dimensional tensor as a, two-dimensional tensor as a, three-dimensional tensor as â, and four-dimensional tensor as â.Unless otherwise specified, the subscript of the tensors denotes the index of the element, such as a m,n denotes the element in a with index {m, n}.We use bold symbols to denote vectorial physical quantities, such as the positions r and the electric fields E in 3-D cases.

3-D Modeling
In 3-D cases, there are N i incident waves from different angles impinging onto the rectangular cuboid DoI D (D ⊂ R 3 , the background 3-D homogeneous medium with permittivity 0 and permeability µ 0 ), where nonmagnetic scatterers are located, and these incident waves are expressed as E inc l (r), l = 1, 2, . . ., N i , r ∈ D. For each incidence, the scattered fields are collected by N r receivers located at r j , j = 1, 2, . . ., N r .With all such information, including every incident field inside the domain of interest and the corresponding scattered fields at the positions of all detectors, we aim at determining the dielectric profile (r), r ∈ D.
The scattering models are governed by the following electric field volume integral equation based on the well-known Lippmann-Schwinger integral equation (LSIE).For the l th incidence, the field equation in the domain D is expressed as where χ(r) = ( r (r) − 1) is the contrast, r (r), I l (r) = χ(r)E tot l (r) and E inc l (r) are the relative permittivity, the contrast source and the incident electric field at r, respectively, and (G 3D D I l )(r) is an integral operator with the dyadic Green's function as the integral kernel, which can be written as [8,31], with k 0 = ω √ 0 µ 0 as the wave number of the background medium, and g(r, r ) as the 3-D Green's function for the background homogeneous medium.Notice that the I l (r) is the contrast source in this paper, while, in [8], it is the physically induced current that includes a multiplicative factor −iω 0 compared to the contrast source.Otherwise, as we know, the dyadic Green's function is composed of nine scalar elements, which represent the mapping from the three components of the contrast source to the three components of the scattered fields inside the DoI.
The nonlinearity of ISPs comes from the MSE, as shown in Equation (1).To alleviate the nonlinearity of the model, in [29], by some mathematical manipulation on Equation (1), another new-type integral equation, which is denoted as CIE-I herein, can be obtained where is the modified contrast function, and β(r) is a chosen auxiliary parameter to control the portion of the MSE in estimating the contrast [29], which can be a constant or a variable at the different position in the DoI.For the convenience of discretizing the equation, we rewrite the dyadic Green's function as where u, v = 1, 2, 3 represent the x-, yand z-components of a vector, respectively, and ι 1 = x, ι 2 = y, and ι 3 = z.Following the conventions in [8], we can rewrite the vectorial Equation ( 3) into three coupled scalar equations in the discrete forms.Thus, we first discretize the rectangular domain of interest into small cuboid subdomains, whose dimensions are much smaller than the wavelength and the center of which are located at r m,n,p , with m, n, p integers and m ∈ Here, M 1 , M 2 , and M 3 are the total number of subdomains along x-, y-, and z-directions, respectively, and we let M = M 1 × M 2 × M 3 be the total number of the subdomains.With such discretization, we have βm,n,p Îl;u;m,n,p = Rm,n,p βm,n,p Îl;u;m,n,p + Rm,n,p Êinc l;u;m,n,p where Rm,n,p = βm,n,p χm,n,p βm,n,p χm,n,p + 1 , χm,n,p is the contrast at r m,n,p , whereas Êinc l;u;m,n,p and Îl;u;m,n,p are the incident electric field and the induced current at r m,n,p , respectively.Subscript u = 1, 2, 3 denotes the x, y, and z components of a vector.Note that the convolution-type operators Ĝ3D D;uv are obtained via the Equation ( 4).For further details of the discretization of Equation ( 1) and the finite difference scheme to generate (G 3D D I l )(r), please refer to the Appendix of [6].Similarly, the integral operator relating the contrast sources and the scattered fields could also be expressed as the summation of the contribution from all the subdomains, where 2, 3 denotes the x, y, z component of the corresponding vector, respectively), I l is a 3M dimensional vector obtained by I l = vec Îl .In a 3-D scenario, the vectorization operation vec {•} is defined to vectorize a four-dimensional tensor into a vector, i.e., if I l = vec Îl , we have In Equation ( 6), the scattering operator is defined as a 3N r × 3M matrix, with G S;uv , a N r × M matrix, the mapping from the v component of the induced current to the u component of scattered fields (the subscripts u, v = 1, 2, 3 are not indexed for tensor elements).The explicit expression of G S;uv is where δ(y) is 1 when y = 0 and is 0 otherwise . ., M 2 , and p = 1, 2, . . ., M 3 .Here, (r) u denotes the u component of r.

Objective Function for Inversions
In this subsection, we build the objective function used in the proposed inversion method, in which we will use the FFT-TSOM [8] as the regularization to stabilize the inversion, as done in [29].Emphasize that the twofold subspace constraints have different regularization effects.The first fold, the original SOM, balances the two mismatches in the objective function, whereas the second fold, the TSOM, is the key to stabilize the inversions with CIE-I.Details are given in below.
For the original SOM part, as given in [16], with the spectral information of G 3D S (the singular value decomposition-SVD of r M), assuming that the singular values σ S j is a non-increasing sequence), the contrast sources can be decomposed into two parts, deterministic part of the contrast sources (DPCS) and ambiguous part of the contrast sources (APCS), the former being obtained as where . ., L, and the superscript * denotes the Hermitian operation while superscript + refers to the dominant current subspace, the subspace corresponding to the dominant singular values.The value of L is chosen according to the noise level [16].Later, we will see how this might work after introducing the objective function for the inversion.
For the second fold subspace constraint, according to the FFT-TSOM [8], the APCS can be written as where γl = γx;l ; γy;l ; γz;l , I and u = x, y, and z, IDFT {} is the inverse discrete Fourier transform operator, the vec{•} is the vectorization operator.Note that the inverse discrete Fourier transform (IDFT) is performed by the 3-D FFT algorithm, the computational complexity of which is By using the low-frequency Fourier components, we are able to constrain the APCS within a low-dimensional subspace.The reason [17] is to use only the contrast sources components in this subspace that is influential to the scattered fields within the DoI, such that the stability of the inversions is substantially increased.To reduce the computational costs, such a subspace can be approximately spanned by low-frequency Fourier bases [8].Here, if we use all Fourier bases, the construction of the APCS becomes the one in the original SOM.Otherwise, we can use the low-frequency components, and set the coefficients for those high-frequency components as zeros.This can be achieved via using a mask with eight corners being 1, with size M F (1 ≤ M F ≤ M 1 /2, M 2 /2, M 3 /2), and other positions being 0. The details can be found in [8].
Having expressed the contrast sources as aforementioned, it is straightforward to define the objective function.Firstly, it is natural to give the mismatch of the scattered fields by where I d l and I a l are as in Equations ( 9) and (10), respectively, and • denotes the L 2 norm of a tensor.The current equation in Equation ( 5) is another key equation to satisfy.Using the APCS construction Equation (10), we define an operator as With this definition, we could write the mismatch of Equation ( 5) as where Γ3D The inversion is to minimize this objective function.As in [8], the conjugate gradient (CG) type algorithm that is used in contrast source inversion (CSI) method is adopted to minimize this nonlinear problem by alternatively updating the γl and R at every iteration of the optimization.

Sketch of the Inversion Method
Following the inversion method in [8,29], we summarize it as follows: 1.
Set the background medium and null APCS as the initial guesses and choose an L value such that the corresponding first L singular values of G 3D S are larger than the noise level (assuming that the noise is a white Gaussian one).

2.
Set proper values for β in CIE-I modeling and proper value for M F to control the number of Fourier bases being used.

3.
Carry out the CG type optimization algorithm to alternatively update the two types of variables, where the APCS is updated with a one-step Polak-Ribière CG scheme and the contrast is updated with the least squares method.

4.
Stop the optimization if a termination condition is met, which can be a maximum number of iterations or a pre-defined relative change of APCS coefficients.

5.
If the maximum number of rounds of inversion is met, go to Step 6.Otherwise, the obtained contrast and APCS will be used as the initial guesses for the next round of optimization with smaller β and larger M F in Step 3.

6.
Output the obtained contrast.
As mentioned in [29], the first round inversion with a large β and with a low-dimensional subspace enables a fast convergence to a meaningful coarse result that could be used as an initial guess in the second round.By increasing the dimension of the APCS subspace and including more MSE in estimating the contrast, the inversion gives us a result with better resolution.Since in the second round a (supposedly) good initial guess is given, the convergence is also very fast.For some difficult problems, the inversion could be carried out with multi-round optimizations by gradually increasing the dimensions of the APCS subspace in each round while decreasing the values of β in estimating the contrast.

Updating Contrast and Contrast Sources with Different β
As mentioned above, there are two updates at each iteration for the two types of unknowns.In [29], the same β is used in each round of optimization for both updates.In the first round optimization in inversions, to tackle the nonlinearity, a large β is needed in updating the modified contrast, and therefore we use the same β for the update of APCS.In the subsequent rounds (instead of the first round) of optimization, for the sake of stability, a small β is used in estimating the modified contrast to cope with a high-dimensional APCS subspace, as discussed in [29].However, there is not such a need to use a small β in updating the APCS.Consequently, in the second round of optimization, we can use a CIE-I model with larger value β 1 for the update of the APCS and another with smaller value β 2 for the update of the contrast at the same iteration.The same applies to the third round optimization, if there is such a need.The purpose of doing so is to further accelerate the convergence of these computationally burdensome 3-D inversions after the first round optimization.
As for the value range for β 2 , we carry out a similar calculation as in [29], which reflects the norm of the scattering operator that maps the contrast sources to the scattered fields within DoI.The results shown in Figure 2 indicate that, if we want to suppress the MSE in estimating the contrast when using the least squares method, we might need to choose a value for β 2 that is larger than 3.5, and this is the guideline for the first round optimization.For the second or subsequent rounds, as good initial guesses are provided, we can include more MSE to retrieve the fine features of the unknown scatterers.

Numerical Simulations
In this section, we will test the proposed inversion methods in three examples.We will use the same physical setup of sources and receivers for the 3-D example in [8], as shown in Figure 3.The DoI are all the same, a box with size 3λ × 3λ × 3λ.The box is illuminated by 60 electric dipole antennas at 300 MHz (λ = 1 m in air), located at three circles (with 20 dipole antennas evenly distributed on each) with the same radius 3 m.The three circles are in x − y, y − z and x − z planes, and their centers are at (0.2, 0, -0.1), (0.1, 0, -0.15), and (-0.05, 0.1, 0), respectively.The direction of the electric dipole sources in the y − z plane are in the x-direction, while those in the x − z and x − y planes are in the yand z-directions, respectively.Scattered fields are measured by 60 receivers, located at the same positions as the 60 dipole sources.We collect all three components of the fields.Consequently, we have a 60 × 180 synthetic data matrix.In all three tests, the synthetic data are calculated with a 60 × 60 × 60 mesh, while, in the inversions, a 30 × 30 × 30 mesh is used.All synthetic data are contaminated with additive Gaussian white noise with level of 10%.The termination condition in each round of optimization is to reach a pre-defined relative changing rate of the APCS coefficients, given as where γ(h) l is the APCS Fourier coefficient for the l th incidence at h th iteration.We set δ 3D < 10 −3 as the termination condition.
The reconstruction results are quantitatively evaluated by the following error: When running on a workstation with eight threads and 32 Gb RAM, every iteration of the proposed inversion method costs about 25 seconds CPU time in MATLAB, including one update on the contrast sources and one update on the contrast.Due to the independence between different incidences, we use 8 MATLAB workers to update the contrast sources in parallel.

Example 1
We reuse the 3-D test in [8], where a coated cube is employed.The cube is with an outer layer ( r2 = 1.5 + i0.3) and an inner layer ( r1 = 2 + i0.8),where the outer edge length is b = 2λ and inner edge length is a = 1λ, as shown in Figure 4a.In [8], we observe that the inversion with LSIE and FFT-TSOM is able to satisfactorily retrieve the cube while the LSIE with SOM fails.Here, we illustrate how the proposed inversion method with CIE-I and FFT-TSOM performs.Similarly, we carry out two rounds of optimization for this inversion as done in [8] with M F = 6 and 10, but with CIE-I, we use β = 6 in the first round and β = 1 for the second round.The reconstruction results are shown in Figure 5, with { ˆ est r } as the real part of the reconstructed relative permittivity and { ˆ est r } as the imaginary part, where we see after two rounds of optimization that the coated cube is successfully found.Now, we compare the convergence speeds of inversions with LSIE and CIE-I.We plot the reconstruction errors for both cases, as shown in Figure 6, in which the inversion with CIE-I appears to be converging much faster than the one with LSIE, with the same reconstruction quality.This confirms again the speeding convergence of inversions in two-dimensional cases shown in [29].We would like to emphasize that this fast convergence property is extremely important in handling the 3-D ISPs, since the computational costs in 3-D problems are much higher than in two-dimensional cases.In the four subfigures, the 30 slices of DoI with each at z = z q plane, where z q , q = 1, . . ., 30, are the grid points along z-direction, and z p < z q if p < q.The displaying sequence is with the convention of left to right, and top to down.For instance, the top left corner one is with z 1 , and the top row second column one is with z 2 .The same applies to the figures hereafter.In the first round, β 1 = β 2 = 6.In the second round, With this example, we further test the relaxed scheme with different β for updating the APCS and the contrast.The termination condition remains the same, δ 3D < 10 −3 .We apply this scheme in the second round optimization, with β 1 = 6 for the APCS update and β 2 = 1 for the contrast update.For comparison, the recorded reconstruction error at each iteration is plotted in Figure 6 as well, which shows that the convergence is further accelerated compared with the case of using β 1 = 1.

Example 2
In this test, we use a profile with four cubes with increasing contrasts.Having edge length 0.8λ, they are centered 0.7, −0.7)λ, (0.7, −0.7, −0.7)λ, (−0.7, −0.7, 0.7)λ, and (0.7, 0.7, 0.7)λ, with r = 2, 3, 4 and 5, respectively, as shown in Figure 7.  Two inversions are carried out with LSIE and CIE-I, each with two rounds (M F = 6 and 10).For CIE-I, the choices of β are 6 and 1 in the two rounds.The reconstruction results and errors are shown in Figures 8-10.From these results, we clearly observe that the inversion with LSIE is only able to find the bottom two cubes with lower contrasts, whereas the inversion with CIE-I succeeds in finding all four of four cubes.Reconstruction errors displayed in Figure 9 confirm it, inversion with LSIE diverging while converging with CIE-I.In addition, we test the inversion scheme with different β for updating the APCS and for updating the contrast.Firstly, we choose β 1 = 6 for the update of APCS and β 2 = 1 for the update of the contrast in the second round of optimization with M F = 10, which however immaturely converges.This might be due to the reason of using a large β 1 .Then, we choose β 1 = 3, and the optimization converges to a good solution with a convergence speed faster than when using β 1 = 1.
From this example, the inversion with CIE-I is able to handle the highly nonlinear problems consisting of strong scatterers, those with either large contrasts or large dimensions (in terms of wavelength), while the inversion with LSIE may fail.

Example 3
We will use another highly nonlinear example to confirm the large difference between the performances of inversion methods with CIE-I and LSIE.In this example, we use a profile that is similar to the famous two-dimensional "Austria" that has an annular and two separated disks.Here, we use a coated cube with a hollow inside and two rods outside the cube, as shown in Figure 11, all of which are with r = 2.5.We see that the inversion with LSIE fails to find the profile while the one with CIE-I succeeds.This can also be seen from the errors of the reconstructions in Figure 13a, where the error by the inversion with LSIE diverges.In the reconstruction results in Figure 14, we see that there is an artefact within the coated cube, where the medium is supposed to be air.This artefact appears in the reconstructed result after the first round optimization with M F = 6, meaning that the inversion, though capable of retrieving most of the main features of the scatterers, is still not stable enough.Therefore, we carry out another inversion with three rounds of optimization, with M F = 4, 8, and 10 in each round, and β = 6, 1, and 0.5, respectively.The final reconstruction result is given in Figure 15, and we see that there is no more such an artefact.

Figure 3 .
Figure 3.The physical setup for the numerical tests.
Real part of the true profile of the coated cube.
Imaginary part of the true profile of the coated cube.
{ ˆ est r } after the first round optimization.{ ˆ est r } after the first round optimization.
{ ˆ est r } after the second round optimization.

Figure 5 .
Figure 5. Reconstruction results of inversions with CIE-I.In the four subfigures, the 30 slices of DoI with each at z = z q plane, where z q , q = 1, . . ., 30, are the grid points along z-direction, and z p < z q if p < q.The displaying sequence is with the convention of left to right, and top to down.For instance, the top left corner one is with z 1 , and the top row second column one is with z 2 .The same applies to the figures hereafter.In the first round, β 1 = β 2 = 6.In the second round, β 1 = β 2 = 1.

1 Figure 6 .
Figure 6.Errors of reconstructions obtained by different inversion methods for Example 1.
We still carry out inversions with LSIE and CIE-I with two rounds of optimization, the first with M F = 6 and the second with M F = 10.For the CIE-I, β 1 = β 2 = 6 is used in the first round, andβ 1 = β 2 = 1 inthe second round.The reconstructions and errors are shown in Figures 12, 13a and 14.

Figure 12 .
Figure 12.Reconstruction results of inversions with LSIE for Example 3.
Errors of reconstructions in inversions with M F = 6 and 10.
Errors of reconstructions in inversions with M F = 4, 8, and 10.

Figure 13 .
Figure 13.Errors of reconstructions obtained by different inversion methods for Example 3.