Distorted Born iterative method with back-propagation improves permittivity reconstruction

Abstract

The distorted Born iterative method (DBIM) is a widely used quantitative reconstruction algorithm to solve the inverse scattering problems of microwave imaging. The major mathematical challenges in solving such problems are non-uniqueness, non-linearity, and ill-posedness. Due to these issues, the optimization algorithm converges to a local minimum. This drawback can be overcome by selecting the correct initial guess solution, which helps to escape local minima and thus guides the inversion algorithm to a satisfactory result. This study uses a back-propagation algorithm to calculate the initial estimate, which significantly accelerates the rate of convergence and improves the accuracy of the standard DBIM approach. The results of this method are compared with zero initialization and Born-approximated initialization. For comparison, weak as well as strong scattering profiles of synthetic and experimental dataset are considered. The results suggest that the proposed method provides a significant improvement in terms of computing cost and efficiency. Furthermore, the proposed technique has the potential to successfully push the limits of reconstructible contrast.

Introduction

Microwave imaging (MWI) technique uses electromagnetic waves to reconstruct the objects located in the region of investigation. It is an evolving technology with many applications in the fields of medical imaging [1], tree health monitoring [2], hidden weapon detection [3], remote sensing [4], and more. MWI is an ill-posed, non-linear inverse problem. The new integral method [5] is one of the most recent methods for reducing the non-linearity of the problem, but the reconstruction in this work is highly dependent on a hyperparameter that changes with the object being reconstructed. In an inverse scattering problem, a total electric field is required to find the exact solution, which is impossible to obtain. Therefore, no exact solution has yet been formulated [6]. However, we can get an approximate inverse solution by using different inversion strategies. They are mainly divided into two categories: linear and non-linear approaches. For strong scattering targets, the linear approaches produce inaccurate reconstructions. As a result, non-linear solutions are important to solve such problems. The two types of non-linear approaches are deterministic and stochastic inversion. Genetic algorithm [7], ant colony optimization [8], simulated annealing [1], differential evolution [9], and particle swarm optimization [10] are examples of stochastic inversion methods that search for the optimal solution of the objective function without requiring a priori knowledge of the target. These approaches, however, have a very high computational cost.

In particular, the contrast source inversion [11], inexact Newton method [12], modified gradient method [13], Born iterative method [14], subspace optimization method [15], and the distorted Born iterative method (DBIM) [16, 17] are the most popular deterministic algorithms. Among them, DBIM is a widely used regularized iterative algorithm because of its ability to reconstruct high-resolution images with a faster rate of convergence. It requires three main procedures to solve the problem: the initial guess distribution, the estimation of the next iterate, and the stopping condition. It is advantageous to begin the inversion procedure with an initial guess solution that incorporates some a priori data [18]. The calculation of a good initial approximation leads to a better solution in terms of convergence rate and accuracy [19]. A frequency hopping technique is utilized in [20], which makes an excellent initial prediction from lower frequencies. This prevents an algorithm from becoming stuck in local minima. However, this increases the system complexity and data acquisition time. The algebraic reconstruction technique [19, 21] also acts as a good initial guess. However, it starts to become unstable as permittivity value increases. In [22], the initial guess is chosen as the objects of arbitrary numbers and shapes. This approach, however, is unable to reconstruct all structures.

To our knowledge, previous research has not reported the effect of the best initial guess on the reconstructed solution of DBIM. Therefore, to address this issue, we conducted a detailed study to find a better initial estimate for DBIM. Several methods are available in the literature to find an appropriate initial estimate. This paper reports on two widely used non-iterative inversion methods of initial estimation: Born approximation and back-propagation. The performance of these methods is compared using relative error (RE). Also, the effect of the regularization parameter on the reconstructed solution is studied. To investigate the reconstruction performance of the proposed approach, numerical examples of synthetic and experimental data [23, 24] are considered. In addition, we have tested the proposed method on the standard Austria profile for various contrasts. The results suggest that the proposed method increases the convergence speed and range of validity of the standard DBIM.

This paper is structured as follows. Section “Formulation” formulates the inverse scattering problem, which explains the mathematical basis of inversion methods, and a detailed description of a proposed algorithm. Numerical simulations and results are presented in Section “Simulation results.” The discussion and brief conclusion of the work is drawn at the end in Section “Conclusion.”

Formulation

The geometrical setup for this problem is shown in Fig. 1. In this configuration, the observation (measurement) domain is denoted by Π, and the investigation domain is denoted by Θ. The unknown scatterer is non-magnetic and dielectric with permittivity of ɛ _r and permeability of μ ₀. The background medium is homogeneous with dielectric permittivity of ɛ ₀ and permeability of μ ₀. The scatterer under investigation is illuminated with N _p incidences using the time-harmonic electromagnetic waves. The incident electric field can be represented as $u_{p}^{i}( r) ,\; \, \ p = 1,\; \, 2,\; \, \ldots ,\; \, N_{p},\; \, \ r \in \Theta$. The scattered electric field $u_{p}^{s}( r) ,\; \, \ r \in \Pi$ is measured in the circular observation domain for all illuminations.

Fig. 1. Basic geometrical configuration of MWI.

In MWI, the forward problem mainly concerns the relation between the incident field and the scattered field. It consists of the following two integral equations [4]:

(1)$$u^{t}\left(r \right) = u^{i}\left(r \right) + k_{0}^{2}\int_{\Theta} g\left(r,\; {r}' \right)J\left({r}' \right)d{r}',\; \quad r \in \Theta$$

(2)$$u^{s}\left(r \right) = k_{0}^{2}\int_{\Pi} g\left(r,\; {r}' \right)J\left({r}' \right)d{r}' ,\; \quad r \in \Pi$$

where J(r) = [ɛ _r(r) − 1] u ^t(r) is a contrast current density, and $k_{0} = \omega \sqrt {\varepsilon _0 \mu _0}$ denotes the wavenumber of the background medium. For simplicity, we introduce the scalar Green's function g(r,  r′) as

(3)$$k_{0}^{2}\int_{\Theta}g\left(r,\; {r}' \right)J\left({r}' \right)d{r}' = \left\{\matrix{G_{\Pi}\left(J \right),\; \quad r \in \Pi \cr G_{\Theta}\left(J \right),\; \quad r \in \Theta }\right.$$

where G _Π and G _Θ are referred as external and internal radiation operators, respectively. Consequently, equations (1) and (2) can be discretized and written in the form of fields as

(4)$$u^{t}\left(r \right) = u^{i}\left(r \right) + G_{\Theta} u^{t} \xi ,\; \quad r \in \Theta$$

(5)$$u^{s}\left(r \right) = G_{\Pi} u^{t} \xi,\; \quad r \in \Pi$$

Here, ξ denotes the contrast function, which can be defined as ξ = (ɛ _r/ɛ _bac − 1), where ɛ _bac is the background permittivity. This is a non-linear equation because contrast function (ξ) and total electric field (u ^t) are both unknown variables. In addition, the operator matrix is severely ill-posed in this case. As a result, such problems are often referred to as non-linear, ill-posed problems.

In the inverse scattering problem, the scattered fields $u_{p}^{s}( r) ,\; \, \ r \in \Pi$ for all incidences N _p are measured in a region of measurement Π (usually disjoint from Θ), and the aim is then to reconstruct ξ(r),   r ∈ Θ. These problems can be effectively solved using the DBIM [16]. This algorithm is an extension of the Born iterative method (BIM). Here, the optimization problem consists of evaluating the difference in contrast function δξ. The cost function is selected as the Euclidean norm of mismatch of the collected and estimated scattered field. To stabilize the optimization, a regularization term is included. Consequently, the cost function can be expressed as

(6)$$O\left(\delta \xi \right) = \sum_{\,p = 1}^{N_p}\left \Vert u_p^s - G_{\Pi} \ \xi \ u_p^{bac} - G^{bs} \ \delta\xi \ u_p^{bac} \right \Vert ^2 + \alpha \left \Vert \delta\xi \right \Vert ^2$$

where ‖ · ‖ represents the Euclidean vector norm, and α denotes the non-negative regularization parameter. This cost function is minimized to estimate δξ. Thereafter the contrast function is updated as

(7)$$\xi_{n + 1} = \xi_n + \delta\xi_n$$

The initial guess distribution (ξ _n=0) is required to start this iterative scheme. A proper selection of an initial guess plays a vital role in obtaining the convergent solution. Generally, it is chosen as zero or the results obtained by non-iterative inversion algorithms [4]. In this work, two such initial guess methods have studied, namely, Born approximation and back-propagation. These methods are easy to implement and produce good reconstruction under certain validity conditions.

After initial guess solution, the total electric field and the system matrix (kernel) are calculated at each iteration. The reconstruction profile of the last step acts as the inhomogeneous background medium, and Green's function for this background is updated as

(8)$$G^{bs}_n = G^s \left(I - \xi_n \ G_{\Theta} \right)^{-1}$$

where G ^bs denotes the inhomogeneous background Green's function operator. Also, the total electric field inside Θ is approximated by the incident field in the presence of a background, which is called the secondary incident field:

(9)$$u_{\,p, n}^{bac} = \left(I-G_{\Theta} \ \xi_n \right)^{-1} \ u_p^{i}$$

After this, the cost function O is estimated and the procedure is repeated until the convergence of the solution is reached.

Initial guess estimation

We solve an under-determined set of equations to obtain a sequence of updates in inverse microwave imaging. Because there are many false minima, convergence to the correct solution is not assured. We employ prior information of contrast, which is available in many real scenarios, to improve the chances of calculating the true contrast. By doing so, we force the solution to have certain predefined properties. In this context, an initial guess solution can be extremely useful in escaping the local minima. As a result, we used two commonly used non-iterative techniques to start the DBIM inversion process in this study.

Born approximation

This model-based approximation is applicable when the scatterer under investigation is weak compared to the background medium [1]. In this case, the scattered field is much weaker than the incident field. Therefore the total field can be replaced by the incident field only, which results in

(10)$$u^{s}\left(r\right)\cong \left(G_{\Pi} u^{i}\right)\xi ,\; \quad r \in \Pi$$

Now the problem of retrieving ξ is ill-posed and linear. Since the scattered field is contaminated by noise, direct inversion produces incorrect results [4]. Therefore, to obtain stable solutions, a regularization strategy must be introduced. If the Tikhonov regularization is applied, the actual linear problem transforms into an optimization problem:

(11)$$f\left(\xi \right) = \sum_{\,p = 1}^{N_{\,p}} \left \Vert \left(G_{\Pi} u^{i}\right)\xi - u_{\,p}^{s} \right \Vert _{\Pi}^{2} + \lambda \left \Vert \xi \right \Vert _{\Theta}^{2}$$

This cost function can be minimized with the help of the singular value decomposition. It decomposes the system matrix as (G _Π u ⁱ) = LSR*, where L and R represent the left and right singular matrices, and S denotes the singular value matrix. After simplification, we obtain the contrast function vector as [25]

(12)$$\xi = \sum_{k = 1}^{\rho} {\sigma_k\over \sigma_k^2 + \lambda} \left \langle l_k,\; u^s \right \rangle r_k$$

where 〈 · , · 〉 represents the inner product, ρ is the rank of the system matrix, σ _k denotes the kth singular value, l _k is the kth left singular vector, and r _k is the k ^th right singular vector.

Back-propagation

It is a fast non-iterative inversion method, which can provide quantitative information [26]. In Born approximation, the original non-linear problem becomes a linear problem, which results in limited applicability. However, in back-propagation, the measured scattered field depends non-linearly on unknowns, and the inverse problem decomposes into numerous linear equations, each of which is solved without iteration [4]. This inverse scheme consists of three steps. The first step involves the determination of the induced current J, which can be expressed as [4]

(13)$$J = \kappa \ G_{\Pi}^{{\dagger}} \ u^{s}$$

where ${\dagger}$ denotes the adjoint operator, and κ is the complex parameter used to minimize the objective function

(14)$$F\left(\kappa\right) = \left \Vert u^{s} - G_{\Pi} \ \kappa \ G_{\Pi}^{{\dagger}} \ u^{s} \right \Vert ^{2}$$

This function can be minimized by taking derivative with respect to κ, which results in [4]

(15)$$\kappa = {\left \langle u^{s},\; G_{\Pi}\left(G_{\Pi}^{{\dagger}} u^{s} \right)\right \rangle\over \left \Vert G_{\Pi} \left(G_{\Pi}^{{\dagger}} u^{s} \right)\right \Vert }$$

From this value of κ, the induced current can be obtained. The next step is to compute the total field u ^t in Θ,

(16)$$u^{t} = u^{i} + G_{\Theta} J$$

Finally, the solution ξ(r) can be calculated as [4]

(17)$$\xi\left(r \right) = {\sum_{\,p = 1}^{N_{\,p}} J_{\,p}\left(r \right)\left(u_{\,p}^{t}\left(r \right)\right)^{\ast}\over \left \vert u_{\,p}^{t}\left(r \right)\right \vert ^{2}}$$

In this work, a hybrid method is proposed to produce fast and accurate reconstructions by combining the effectiveness of the back-propagation method in escaping from local minima with the ability of standard DBIM in iteratively finding a solution. Also, this method has the potential to successfully push the limits of reconstructible contrast. The simulation results for this method are produced using various examples, and the corresponding reconstruction results are reported in the following section.

Simulation results

In this section, numerical simulations are performed to evaluate the performance of the proposed algorithm. Here, Θ is a square of a side of 20 cm. The background is a free space with ɛ _bac = 1. The simulations are performed at a frequency of 3 GHz. The algorithms with an initial guess as zero, BA, and BP are referred to as O-DBIM, BA-DBIM, and BP-DBIM, respectively. Here, O-DBIM symbolize ordinary DBIM, with ɛ(r) = 0,  r ∈ Θ as the initial guess.

To quantify the efficiency of the stated algorithms, RE is calculated, which is defined as [4]

(18)$$RE = {\left \Vert \varepsilon_r - \hat{\varepsilon_r} \right \Vert \over \left \Vert \varepsilon_r \right \Vert }$$

where ɛ _r and $\hat {\varepsilon _{r}}$ represent the actual and the reconstructed permittivity, respectively. For describing convergence, the logarithmic scale is more useful. As a result, the plotting of this error parameter is done using this scale.

Synthetic data

In this section, we present some synthetic data (Austria and MNIST test image) results. Here, the investigation domain is surrounded by 36 receivers, located equidistantly in a circular fashion, at a distance of 60 cm from the center. The imaging domain is discretized to 961 (31 × 31) square cells.

Austria profile

In the first example, the Austria profile is considered, which is one of the commonly used and challenging configurations. It consists of one ring structure and two disks, as shown in Fig. 2. The thickness of the ring is 3 cm, with the outer and inner radii equal 6 and 3 cm, respectively. It is centered at (0, − 2) cm. The disks of radius 2 cm are centered at (− 3, 6) cm and (3, 6) cm. The permittivity distribution on this target is uniform.

Fig. 2. Permittivity distribution for Austria profile.

To test the effectiveness of the proposed method, it is evaluated on the Austria profile for various contrasts (ɛ _r = 1.5,  2,  2.5,  2.6,  3). The permittivity value of 1.5 represents the weak scatterer, and the object becomes the strong scatterer as the permittivity value increases. The simulation results after 30 iterations are reported in Fig. 3. Also, the behavior of RE with respect to iteration number is reported. According to the results, all of the approaches can efficiently reconstruct weak and moderate scattering objects. However, as the contrast increases, the O-DBIM and BA-DBIM techniques begin to converge more slowly. They also do not acquire convergence for highly contrasting objects. The BP-DBIM, on the other hand, can converge faster than the other approaches under all conditions. As a result, it can be concluded that the proposed technique can effectively push the boundaries of reconstructible contrast.

Fig. 3. Permittivity reconstruction for Austria profile with different permittivities: ɛ _r = 1.5 (first row), ɛ _r = 2 (second row), ɛ _r = 2.5 (third row), ɛ _r = 2.6 (fourth row), ɛ _r = 3 (fifth row).

The performance of the BP-DBIM is further evaluated by checking the effect of regularization parameter (α) on the reconstruction. The trajectories of RE for different α values are reported in Fig. 4. It has been observed that a reasonable variation in the value of α does not affect the solution. As iterations progress, RE converges at the same point for all the values of α. This avoids the selection of optimum regularization parameter, making BP-DBIM the robust reconstruction algorithm.

Fig. 4. Effect of α on microwave image reconstruction.

MNIST test image

In this example, the unknown scatterer is modeled using a test image from the MNIST database [27], which contains handwritten digits. As shown in Fig. 5, the permittivity distribution of the object is homogeneous with ɛ _r = 2.75 (ξ = 1.75). As a result, it can be considered a strong scattering profile.

Fig. 5. Permittivity distribution for MNIST test image (ɛ _r = 2.75).

Simulation results using the above-mentioned methods are displayed in Fig. 6. The results show a close match between the original and reconstructed profiles for the BP-DBIM method. For this case, the RE values for O-DBIM, BA-DBIM, and BP-DBIM after 30 iterations are 0.4211, 0.3748, and 0.1946, respectively.

Fig. 6. Permittivity reconstruction for MNIST test image after 30 iterations: (a) O-DBIM, (b) BA-DBIM, (c) BP-DBIM, and (d) RE.

Experimental data

To examine the performance of the BP-DBIM algorithm under realistic conditions, experimental data (made available by the Institute of Fresnel, France) are also considered. The measurement setup consists of linearly polarized horn antennas. Here, the simulation is performed on two examples. The first example consists of a simple homogeneous, circular cylinder from the Fresnel dataset 2001 [23]. Similarly, the second example consists of two circular dielectric cylinders from Fresnel dataset 2005 [24].

Single dielectric cylinder

The scatterer is illuminated at an angle of 0, 10, …, 350°, and the scattered fields are collected at an angle of 60, 65, …, 300° with respect to the corresponding emitter. The scatterer profile consists of a single circular, homogeneous, dielectric cylinder of permittivity 3 ± 0.3. The cylinder has a radius of 1.5 cm and is kept at a distance of 3 cm from the origin. The reconstructed permittivity distributions at a frequency of 3 GHz are presented in Fig. 7. It can be observed that O-DBIM and BA-DBIM provide less accurate results, whereas the quality of the reconstruction obtained using BP-DBIM is quite good.

Fig. 7. Permittivity reconstruction for experimental data: (a) original profile of single dielectric cylinder from the Fresnel dataset 2001, (b) O-DBIM, (c) BA-DBIM, (d) BP-DBIM, (e) 1D plot along y-axis (x = 0), and (f) RE.

Foam and plastic cylinder

In this example, the scatterer is illuminated at an angle of 0, 45, …, 315°, and the scattered fields are collected at an angle of 60, 61, …, 300° with respect to the corresponding emitter. As shown in Fig. 8(a), the scatterer profile consists of two dielectric cylinders. The first one consists of a centered foam cylinder with a radius of 4 cm and permittivity of 1.45 ± 0.15. Another target consists of a plastic cylinder of radius 1.55 cm, permittivity 3 ± 0.3, placed at a distance of 5.55 cm from the origin. The reconstructed permittivity distributions at a frequency of 3 GHz are displayed in Fig. 8. Here again, results indicate that the BP-DBIM provides superior results as compared to O-DBIM and BA-DBIM.

Fig. 8. Permittivity reconstruction for experimental data: (a) original profile of FoamDielExtTM from the Fresnel dataset 2005, (b) O-DBIM, (c) BA-DBIM, (d) BP-DBIM, (e) 1D plot along x-axis (y = 0), and (f) RE.

Conclusion

In this work, a back-propagation-based initial guess estimation method is reported for the non-linear model of the DBIM. Initially, this method (BP-DBIM) is tested on the synthetic data, which consists of complicated scattering objects that are both weakly and strongly contrasting. The reconstruction results are compared with zero initial guess (O-DBIM) and Born-approximated initial guess (BA-DBIM). In addition, the effect of the regularization parameter (α) on the RE is also studied. It has been observed that a reasonable variation in the value of α does not affect the solution. Results show improved reconstruction in terms of accuracy, and convergence rate. Thereafter, the algorithm is tested on the realistic data provided by the Institute of Fresnel, France. Here also, the superior reconstruction results are obtained.