Efficient Born Iterative Method for Inverse Scattering Based on Modified Forward-Solver

The traditional Born iterative method (BIM) and distorted Born iterative method (DBIM) are both effective for electromagnetic scattering inversion. But they are inefficient because of the computation of forward scattering problem in each iteration. In this paper, a modified forward-solver scheme based on Born series is proposed for improving the efficiency. By using this scheme, the traditional time-consuming forward scattering procedure, matrix inversion, is replaced by a simple multiplication with low computation complexity in each iteration. Numerical and experimental two-dimensional electromagnetic scattering models are used to validate the proposed method. The imaging results and CPU times indicate that the presented methods have good accuracy and can improve the efficiency greatly.


I. INTRODUCTION
Electromagnetic inverse scattering methods have wide applications in nondestructive testing, biomedical imaging, target identification, etc. Their goal is to solve the inherent non-linear relationship between known and unknown information. Usually, integral equations are used to describe this relationship due to the accuracy and simplicity. Therefore, based on it, the Born iterative method (BIM) was proposed to reconstruct the permittivity in the region of interest (ROI) [1]. In this method, the Green's function remains unchanged in each iteration so it needs many iterative steps and much CPU time to get accurate results. Distorted Born iterative method (DBIM) updates the Green's function by the new computed permittivity of ROI in each iteration [2]. Therefore, it can accelerate the convergence speed and save the computation time. But DBIM has worse robustness than BIM.
Both BIM and DBIM have attracted much attention because they are effective and accurate approaches to solve the nonlinearity of inverse scattering. In time-domain, Chew et al. applied BIM to inverse the permittivity profiles of electrically large scatterers [3] and utilized DBIM to image The associate editor coordinating the review of this manuscript and approving it for publication was Abhishek K. Jha . buried dielectric cylinder [4]. By taking the advantages of BIM and DBIM, a hybrid iterative method is presented in [5], where DBIM is used to accelerate the convergence speed and BIM is used to stabilize the final solution, respectively. For the strong scattering medium, multiple frequency DBIM is developed in [6], where low frequency DBIM solution is used to initialize the procedure at high frequency. However, the forward scattering calculation in iterations is still timeconsuming. Subsequently, for the acceleration of the forward solver, conjugate and biconjugate gradient method with fast Fourier transform (BCGS-FFT) are used in [7] and [8]. In [9], the multilevel fast multipole algorithm is utilized to solve the forward scattering with DBIM to ameliorate the computation burden. In [10] and [11], subspace-based DBIM and variational-BIM are proposed to achieve high efficiency and accurate imaging. In [12], a new scheme of controlling the balance factor to improve the convergence speed is proposed based on BIM. Additionally, BIM and modified BIM have implemented in medical electromagnetic imaging to retrieve more accurate images in [13] and [14].
Although the forward procedure of BIM can be accelerated by some techniques under some specific situation, it is still complicated and time-consuming to calculate the inverse of matrix. As we know, in the traditional BIM and DBIM, the forward-solver calculates the accurate total fields with the updated permittivity distribution in each iteration. In fact, it is unnecessary. The updated permittivity distribution is not the final solution if the iteration continues. Therefore, it can be concluded that the high precision forward scattering calculation in the process can be avoided. In this paper, a fast forward-solver scheme based on Born series is proposed to avoid the traditional time-consuming forward scattering calculation of BIM and DBIM. In the new scheme, the total field in each iteration can be calculated quickly by the former value and the updated permittivity distribution. It leads to the acceleration of iteration. The proposed methods are validated by numerical and experimental examples, which show that they have good accuracy and high efficiency.

II. FORWARD PROBLEM A. SCATTERING MODEL
A two-dimensional (2-D) medium imaging model is considered in Fig. 1, which consists of a bounded homogeneous background domain D with known relative permittivity ε 0 . The domain D contains an inhomogeneous medium with unknown permittivity distribution in xoy-plane. In this paper, the ROI is illuminated by transverse magnetic (TM) plane wave from different directions. For each incidence, the scattered field can be obtained by the receivers, R, which are located around the domain D. For the total field inside the domain D, the following scalar integral equation can be gotten as: where E t z (r) and E i z (r) denote the total and incident electric field, respectively; k 0 is the wavenumber of homogeneous medium; is the 2-D Green's function in homogeneous medium where j represents the imaginary unit and H (1) 0 is the zero-order Hankel functions of first kind; and is the unknown relative permittivity distribution profile where ε r (r) is the real relative permittivity of inhomogeneous medium. It is the key parameter which needs to be reconstructed.
Discretizing the domain D into N discrete squares D 1 , D 2 , . . . , D N , Eq. (1) can be expressed in matrix form as where E t D and E i D are the total and incident field vector within domain D, respectively. G D is the 2-D Green's function matrix counting for internal total field and X is the relative permittivity vector of domain D.
If the relative permittivity distribution of the domain D is known, it means X is known. Therefore, the total field vector E t D can be gotten according to Eq. (4) by where I is the identity matrix. The total field can be calculated by Eq. (5) with known incident field and relative permittivity distribution X. While it can be accelerated by some scheme such as BCGS-FFT method in [8] under some specific situation, the computational cost of the matrix inverse is still large especially for the large complicated cases. If the total field inside the domain D is known, the scattering field at the receivers can be expressed from Eq. (1) as: Therefore, in matrix form, Eq. (6) can be described as: where E s R is the scattered field vector at the receivers, G S is the 2-D free space Green's function matrix counting for the scattered field, and E t D is the total field vector based on the relative permittivity distribution X. With the distorted Green's function, Eq. (7) can be transformed as where n is the number of iteration, E s,n R and G n S are the scattered field vector and the Green's function matrix based on relative permittivity distribution X n , respectively.

B. BORN APPROXIMATION AND BORN SERIES
Since the total field is the sum of the incident and scattered field, Eq. (4) can be written as where E s D is the scattered field vector within the domain D.
If the scattered field is small compared to the total field, the first-order Born approximation can be gotten by using the incident field to substitute the total field as where the Roman numeral, 1, represents the order of Born approximation. It can be seen that only the first-order scattering component is considered, which will limit the scope of the method. Therefore, high-order approximation is considered, and the total field can be expressed as According to Born series expansion, the total field has another expression as: where E s,0 From it, in can be found that the infinite n makes the result in Eq. (11) converge to the actual total field. Obvious, there is a certain convergence range for above Born series and the range is determined by whether G D · X ≤ 1.

III. INVERSE PROBLEM A. BORN AND DISTORTED BORN ITERATIVE METHOD
The purpose of inverse scattering problem is to reconstruct the relative permittivity matrix X of ROI with the known data E s R and E i D . If the total field E t D is known, the inverse scattering problem can be solved directly by inversing Eq. (7) with the regularization scheme in [1].
However, in practice, the total field within the domain D is unknown so that the inverse problem is nonlinear. By BIM or DBIM, the solutions can approach to the real solution by the iterations. The iteration scheme of the two methods is shown in Fig. 2.
It can be found that in Step 2 the total field is updated by the forward solver as Based on Born series, Eq. (13) can be expanded as Obviously, this procedure is time consuming due to the calculation of matrix inverse. Additionally, in Fig. 2, σ is the threshold value. The δ n is the relative residual sum of squares between the calculated and real scattered field at the receivers, which is defined as

B. EFFICIENT ITERATIONS BASED ON NEW FORWARD-SOLVER SCHEME
According to the iteration scheme in Fig. 2, it can be found that the traditional BIM and DBIM contains two kinds of iteration. The first one is calculating the accurate total field within Step 2 by inversing the matrix, which is equal to the infinite-order Born series iteration in Eq. (14). The second one is getting the accurate permittivity distribution by the updated total field through Step 2 to 4. Obviously, as mentioned before, the first kind of iteration is unnecessary. As long as the final total field converges to the accurate solution through the second kind of iteration, the permittivity distribution can be reconstructed correctly. Therefore, a new forward solver based on high-order Born approximation is proposed as Substituting Eq. (16) for Eq. (13) in Step 2 of BIM and DBIM, the total field can be expanded as VOLUME 8, 2020 Similar to Eq. (11), Eq. (17) can also be considered as Born series expansion as where E

C. DISCUSSION OF THE NEW FORWARD-SOLVER
It can be found that E t,n D in Eq. (17) will be close to the value in Eq. (14) when n is infinite. Hence, Eq. (16) can also make the final total field converge to real value by the second kind of iteration with the increasing iterative steps n because the high-order Born series is considered.
For the new forward-solver, it can be found that the values of X 0 , X 1 , . . . , X n−1 are used in Eq. (17) while only the last X n−1 is used in Eq. (14). But it does not mean that they are all stored in the memory. According to Eq. (16), only the last E t,n−1 D is used for the calculation of E t,n D . All previous values are contained in E t,n−1 D . Therefore, there is no extra memory cost.
In terms of computation complexity, it is clear that, by the new forward-solver of Eq. (16), an approximate total field is gotten in Step 2 quickly because the previous calculated total field is reused. This forward procedure avoids to dealing with the large matrix inverse and results in low computation complexity. If the matrix size of G D · X n−1 is M × M , the complexity can be reduced from O(M 3 ) to O(M 2 ).
In addition, in terms of limitation, as mentioned in Section II, the convergence is determined by whether G D · X ≤ 1. Therefore, for the proposed methods, they are limited by the value of ||G D · X n−1 ||.

A. NUMERICAL EXAMPLE
In this section, numerical and experimental examples are given to validate the proposed methods. Since the proposed forward-solver is implemented within the framework of BIM and DBIM to improve the efficiency, the proposed methods are mainly compared with them. A typical 2-D numerical example with ''Austria'' permittivity distribution is considered, firstly. As shown in Fig. 3, the ''Austria profile'' consists of two discs and one ring. More details of the ''Austria profile'' can be referred in [10]. The operating frequency is 0.15 GHz, which makes the domain D is 1.5λ × 1.5λ (λ is the wavelength). To avoid inverse crime, the forward scattering data is generated by finite-difference method in frequency-domain (FDFD) where the mesh size is λ/60. The inversion mesh size is λ/20. 36 receivers are set around the domain D and the incident wave is illuminated from 18 different directions. In the inversion, since the initial scattered field data, generated by FDFD, has numerical residual itself, the threshold value σ is set as 0.05% based on the initial numerical residual. And priori information that relative permittivity parameter is real number and greater than 1 is used.  By the inversion procedures, the results by four methods are shown in Fig. 4, respectively. It can be found that both BIM and DBIM can properly reconstruct parameter distribution and yield good results, where the two discs and one ring are well distinguished. And the results by the modified BIM and DBIM are almost the same with the results by BIM and DBIM, respectively. Therefore, good agreement is achieved in this paper.
Further, for evaluating the quality of results, the norm error ε is defined as where ε simulated r and ε real r are the simulated and real relative permittivity. Table 1 gives a comparison of four methods in terms of CPU time, iteration steps, and ε, respectively. It can be found that BIM method completes the whole operation with the best quality and maximum CPU time. Compared with BIM, aided by the proposed fast forward-solver, the modified BIM can provide comparative result with much less CPU time. Although the modified BIM needs more iteration steps to converge, the total computation efficiency is still very high due to the low computation complexity. As for DBIM-based methods, similar conclusion can be gotten that the modified DBIM can yield almost the same result with DBIM but with less CPU time.
In order to further evaluate the convergence of the modified methods, the convergence curves (scattered field difference δ n and defined norm error ε as functions of number of iterations steps) of four methods are given in Fig. 5. It can be observed that, compared to BIM and DBIM, the modified BIM and DBIM take a few more steps to converge, respectively. After the convergence is completed, the convergence curves of modified BIM and DBIM are almost consistent with  the convergence curves of BIM and DBIM, respectively. The above phenomenon proves that proposed fast forward-solver can maintain the same precision level meanwhile reducing a lot of CPU time.

B. EXPERIMENTAL EXAMPLE
The second example is based on experimental data with real noise, which is collected by Institute Fresnel in database [15]. The configuration of experiment is shown in Fig. 6.  In experiment, 49 receivers and 36 transmitters are placed on the receiver rail and transmitter rail, respectively. A dielectric cylinder with 1.5 cm radius is placed 3 cm below the center of circle rail. The relative permittivity of cylinder is 3 ± 0.3 and the data at 1GHz is used. The mesh size is 0.3 cm in both directions, which results in 61 × 61 pixels in the ROI domain S. More details of the experiment can be referred in [15].
The supplied data only contains the received incident field and total field at receivers. The incident field in ROI, which is not contained in the supplied experimental data, is an essential input information to start the inversion and it usually needs to be estimated [15]. There are several methods [15]- [17] which can approximate the incident field. In this paper, a complex coefficient factor is used to produce the raw data as [16].
The reconstructed results by four methods are shown in Fig. 7, respectively. It can be found that the object is well reconstructed. The accuracies by the BIM-and DBIMbased methods are almost same with each other, respectively. The experimental example proves that the proposed fast forward-solver can also handle real noisy data. A comparison of four methods is also given in Table 2. It can be found that the CPU time of the modified-BIM and DBIM is much less than those by the conventional methods and the errors are almost same.

V. CONCLUSION
In this paper, a new forward-solver scheme based on Born series is proposed and applied in BIM and DBIM. Compared with the traditional calculation of the total field, the proposed forward-solver yields an approximate solution to improve the efficiency. Therefore, the total field can finally converge to the accurate solution by considering high-order components, which guarantees the veracity of the whole inversion operation. The proposed method is well validated by two typical examples based on numerical and experimental data, respectively. The reconstructed results and convergence curves of numerical example show that the modified methods can fit in the framework of BIM and DBIM well and improve the efficiency greatly. Meanwhile, the experimental example shows that modified methods are also capable of dealing with realistic case. In general, the fast computation speed of the proposed method is very helpful for imaging systems.