2218
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
Research Article
Multi-resolution subspace-based optimization method for solving three-dimensional inverse scattering problems XIUZHU YE,1 LORENZO POLI,2 GIACOMO OLIVERI,2 YU ZHONG,3 KRISHNA AGARWAL,4 ANDREA MASSA,2,5 AND XUDONG CHEN6,* 1
Department of Electrical and Computer Engineering, Beihang University, Beijing 100191, China ELEDIA Research Center@DISI, University of Trento, via Sommarive 5, 38123 Trento, Italy 3 A*STAR, Institution of High Performance Computing, Singapore 138632, Singapore 4 Singapore-MIT Alliance for Research and Technology (SMART) Centre, CREATE Tower, Singapore 138602, Singapore 5 Laboratoire des Signaux et Systèmes UMR8506 (CNRS-Univ. Paris Sud-CentraleSupélec), 3 rue Joliot Curie, 91192 Gif-sur-Yvette, France 6 Department of Electrical and Computer Engineering, National University of Singapore, Singapore 117576, Singapore *Corresponding author: [email protected] 2
Received 10 June 2015; revised 30 September 2015; accepted 30 September 2015; posted 30 September 2015 (Doc. ID 242704); published 28 October 2015
An innovative methodology is proposed to solve quantitative three-dimensional microwave imaging problems formulated within the contrast source framework. The introduced technique is based on the combination of an efficient iterative multiscaling strategy aimed at mitigating local minimum issue arising in inverse scattering problems, and a local search algorithm based on the subspace-based optimization method (SOM) devoted to effectively retrieving both the “deterministic” and the “ambiguous” parts of the unknown contrast currents. To achieve this goal, a nested iteration process is adopted in which the outer loop iteratively refines the region of interest (ROI) where the scatterers are detected, while the inner loop retrieves the dielectric properties of the scatterers within the ROIs. Selected numerical examples are also given to show the validity and robustness of the proposed algorithm in comparison with state-of-the-art techniques. © 2015 Optical Society of America OCIS codes: (290.3200) Inverse scattering; (180.6900) Three-dimensional microscopy. http://dx.doi.org/10.1364/JOSAA.32.002218
1. INTRODUCTION Active imaging methodologies [1–3] are aimed at retrieving the electromagnetic properties of the unknown scatterers located inside inaccessible domains probed by known incident waves. The importance of these techniques is motivated by their very wide domain of applications, which comprises geological exploration [4], nondestructive evaluation [5–7], biomedical [8,9], and through-wall imaging [10,11], etc. Unfortunately, solving a general inverse scattering problem is challenging due to the large number of unknowns, high nonlinearity, and ill-posedness [12–16]. Several inversion methods exploiting linear and nonlinear formulations of the imaging problem have been developed in the literature to cope with these theoretical issues. However, linear inversion methods, such as those based on the Born approximation [16], are often limited to weak scatterers. On the other hand, fully nonlinear algorithms preserve all the multiple scattering effect and are more accurate when dealing with strong scatterers [13], but they require the solution of a much more challenging inverse problem. So far, most 1084-7529/15/112218-09$15/0$15.00 © 2015 Optical Society of America
inversion methods are developed under the two-dimensional (2D) scenarios and only few nonlinear algorithms have been developed for solving three-dimensional (3D) inverse scattering problems [16]. Apart from the aforementioned difficulties, in comparison with the 2D scenario, 3D ones yield to (i) a larger amount of data and a rapid growth of the computational cost with the size of the investigation domain, (ii) a much more difficult model to be generated and visualized, and (iii) a vectorial formulation of the problem equations rather than a scalar one. However, the development of 3D imaging algorithms is of fundamental importance because of the 3D nature of most real-world problems as well as the increased information (e.g., concerning the field polarization) that can be exploited as compared to simpler 2D methodologies [17]. In this framework, a powerful nonlinear algorithm, namely the subspace-based optimization method (SOM) [18], has been proposed to image 3D scenarios within the contrast source formulation of the inverse scattering problem. This approach, which has also been widely deployed in the 2D case,
Research Article
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
decomposes the unknown contrast source into two orthogonal parts: the deterministic part, which is uniquely determined by the measured scattered field and the first several leading singular values of the mapping from domains of interest to receivers, and the ambiguous part, which is found by solving a nonlinear optimization problem through a local search strategy [19]. The singular value decomposition (SVD) procedure is an effective tool widely used in the inverse problem society and the role of diversity on the singular values of linear scattering operators has been thoroughly studied in [20]. SOM has effectively extended it to the decomposition of induced current. Such a procedure has twofold advantages. First, the deterministic part serves as a good initial guess to the optimization, which speeds up the convergence process. Second, the optimization is carried out in a smaller current “subspace,” which actually mitigates the nonlinearity of the problem [19]. Nevertheless, the SOM method still suffers from local minimums [21] and computational complexity issues when large domains or high resolutions are required, since (a) it yields to a deterministic optimization strategy in a high-dimensional space, and (b) it requires the computation of the singular value decomposition (SVD) of the Green’s operator for the evaluation of the deterministic currents. As a consequence, the integration of SOM with an effective multiresolution strategy, namely, the iterative multiscaling methods (IMSA) [22], has been proposed in [23]. Such a choice has been motivated by the fact that the IMSA can mitigate local minima issues [21], since it is able to keep the ratio between the number of unknowns and the information content of the field as close to unity as possible at each iteration. Moreover, each multizooming step requires the solution of a relatively small problem, thus reducing the computational complexity of the associated SVD [22]. However, the arising IMSA-SOM method has been developed and validated only in the 2D case [22]. Accordingly, this paper aims at generalizing the multiresolution SOM in order to yield an efficient and reliable imaging method to solve the general 3D inverse scattering problems. To validate this introduced technique, the proposed IMSA-SOM algorithm is tested in a selected set of numerical examples dealing with scatterers with different levels of complexity and/or different noise levels. The paper is organized as follows. In Section 2, the formulation of the inversion problem and the detailed IMSA-SOM procedure designed for the 3D case are introduced. Numerical examples are given in Section 3 to indicate the validity and efficiency of the proposed algorithm. Conclusions are finally drawn in Section 4. 2. PROBLEM FORMULATION AND IMSA-SOM SOLUTION APPROACH A. Formulation of the Problem
Consider a 3D setup, as shown in Fig. 1, in which an unknown scatterer (which is an isotropic one) is embedded inside a domain of interest D with lossless and nonmagnetic background of permittivity ε0 and permeability μ0, respectively. The domain D is successively illuminated by V monochromatic plane waves E¯inc p r, p 1; …; V with time dependence expjωt (omitted hereinafter for the sake of conciseness).
2219
Fig. 1. Problem geometry.
The amplitudes for the two orthogonal components of the incident electric field with respect to the polar angle and the azimuth angle are Aθ and Aϕ , respectively. For each illumination, the incident electric field E¯inc p r and the total electric field E¯tot p r are then collected by M receivers uniformly distributed in the xy, yz, xz planes around the circular observation domain, respectively, with positions rq0 , q 1; 2; …; M . In the above assumptions, the following electric field integral equation holds true: Z inc tot 2 ¯ ¯ E r E r − k0 ∇∇· gr; r 0 τr 0 E¯tot r 0 dr 0 ; D
(1) where k0 denotes the wave number in the homogenous background medium, τr εr r − j σr ωε0 − 1 is the contrast −jk jr−r 0 j
0 function of the object, and gr; r 0 e4πjr−r is the scalar 0j Green’s function. It is worth remarking that in the 3D case the operator ∇∇· is not a scalar quantity, unlike in the 2D case. In the following, Eq. (1) will be discretized by the method of moment (MoM) assuming a standard voxel basis and point matching strategy [15].
B. Proposed Inversion Methodology
In order to apply the IMSA-SOM, the guidelines presented in [23] are followed and suitably generalized. More in detail, a nested procedure is adopted that comprises the following two loops: (a) An external multiresolution loop that is devoted to adaptively identifying and zooming the region of interest (ROI) where the scatterers are detected. More specifically, at the sth step of the IMSA iteration procedure (s 1; …; S), the ROIs are meshed into N s cubic voxels, the centers of which s are located at rs n , n 1; 2; …; N , and the contrast currents s s ¯ I n , n 1; 2; …; N are assumed constant within each voxel (the N s is determined by the degrees of freedom (DOF) of the scattered field [24]). At the end of each sth step (s 1; …; S), a filtering and clustering procedure is employed to update the ROIs adaptively. (b) Internal SOM loop that comprises the computation of the deterministic currents within the ROI by means of a truncated SVD procedure aimed at calculating the deterministic current, and the evaluation of the ambiguous part of induced contrast currents through the subspace optimization strategy [18,26].
2220
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
According to the above guidelines, Eq. (1) is expressed in the observation domain in the following matrix form for the sth external iteration step: p 1; …; V ; (2) E¯sca G¯ s I¯s ; p
s
p
where E¯sca p is a 3M -dimensional vector (M being the number of field probes) containing the scattered field values for the pth incidence, and the superscript “sca” stands for the scattered s field. Moreover, I¯s p is a 3N -dimensional vector containing the unknown current coefficients. Furthermore, the discretized version of Eq. (1) can be expressed inside the investigation domain as (3) I¯s ζ¯s fE¯inc gs G¯ s I¯s ; p
p
D
p
¯s ¯s where ζ¯s diagζ¯s 1 ; …; ζ N s , ζ n is a 3 × 3 dimensional diagonal matrix with the entries on the diagonal as τn [18] s (isotropic dielectrics are assumed), and fE¯inc p g denotes the incident field vector for the pth incidence in the sth outer loop. The internal Green’s operator G¯ s D maps the contrast current to the scattered electric field inside ROI, while the external Green’s operator G¯ s s maps the induced current inside ROI to the scattered field on the receivers [26]. The detailed entries of the two mapping operators can be found in [18]. s s ¯ s s By applying the SVD to G¯ s s one yields G s · v¯m χ m u¯m , s s where the χ s m , v¯m and u¯m , m 1; …; 3M denote the singular values, right and left singular vectors, respectively. The singular values of G¯ s s can be ordered according to their magnitude (from larger to smaller) and the induced current can thus be split into two orthogonal parts, namely the deterministic and the ambiguous part, as follows: ¯det s ¯amb s I¯s p fI p g fI p g ;
(4)
s is calculated without where the deterministic part fI¯det p g resourcing to any prior information [18]: s
s fI¯det p g
LA sca X u¯s m · E¯ m1
χ s m
v¯s m:
(5)
The value of the singular value truncation numbers Ls A is selected adaptively by resourcing to a user-defined control threshold α that is obtained by the fraction of the total spectrum energy: PLs s s m1 χ m s 1; …; S; LA arg min P3M s − α ; Ls m1 χ m PLs s χ m − α > 0: (6) subject to Pm1 s 3M m1 χ m
In Eq. (6), χ s m , m 1; …; 3M denotes the singular values of the mapping G¯ s s . The ambiguous part is then constructed by employing the remaining 3N s − Ls A current basis as s fI¯amb p g
3N s −Ls A
X
m1
v¯s
mLs A
¯ s · β¯s : βs m V
(7)
Indeed, the contribution of Eq. (7) to the scattered field on the receivers is negligible [18]. It should be noted that the computational complexity of the SVD has been reduced
Research Article
in [18] since only a thin SVD is needed. This ambiguous part of the currents is then retrieved by minimizing the following cost function [18]: F s ζ¯s ; β¯s ; …; β¯s V
1
0
1 ¯s ¯sca 2 ‖G¯ s s Ip − Ep ‖ 2 V B C X ‖E¯sca p ‖ B C B s s C; s ¯s 2 A ¯ ¯ inc s @ ¯ ¯ ‖I p − ζ fE p g G D I p ‖ p1
(8)
s 2 ‖fI det p g ‖
where ‖ ‖ is the Euclidean norm through a local search strategy as detailed in [18]. The arising IMSA-SOM iteration process for the 3D scenario can be summarized as follows (s being the number of outer iteration steps and t being the number of inner iteration steps): Step 1: Outer loop: initialization s 1. s
s
s
Step 1.1: Calculate G D , G s and the SVD of G s . Obtain s fI¯det p g , p 1; 2; …; V by using Eq. (5). Step 1.2: Initialization: t 0, set a proper initial guess to the contrast τ¯s 0 0 [25], and to accelerate the convergence process we obtain β¯s 0 by back propagation, which is different from [25]. Step 1.3: Two step conjugate-gradient type optimization [27]: t t 1. Step 1.3.1: Update β¯s t . Step 1.3.2: Update τ¯s t . Step 1.4: Stop iteration when a user-defined maximum number of iterations T s is reached or a user-defined threshold is reached. Step 2: Stop iteration if a user-defined outer loop number S is reached, otherwise identify/update the ROIs: s s 1 and go to step 1. The details of the two step conjugate-gradient method are omitted for concision and we refer the readers to [22] and [27] for details. 3. NUMERICAL EXAMPLES The numerical analysis presented in this section has a twofold objective. First, it aims to analyze the sensitivity of the proposed approach on its control parameters (i.e., the fractional parameter α and the number of IMSA steps S) and to provide a guideline on their choice. Second, it is devoted to the assessment of the reconstruction accuracy and efficiency as well as the robustness of the IMSA-SOM when dealing with various scatterer and noise levels. Comparisons between the IMSA-SOM with the BARE-SOM as well as to state-of-the-art multiresolution techniques are also provided to show its advantages and drawbacks. Toward this end, a benchmark cubic investigation domain D with side length 1.5λ is assumed to be illuminated by V 30 plane waves with incident angles uniformly distributed around a circle in the y–z plane. The amplitudes for each component of the incident electric field are Aθ 0.7 and Aϕ 0.7, respectively. For each incidence angle, M 30 probes are evenly distributed on circles of radius 3.2λ in the xy, yz, xz planes, respectively, to measure the field data.
Research Article
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
2221
Synthetic data are generated by dividing the domain of interest D into 20 × 20 × 20 voxels. To avoid inverse crime issues, the domain is partitioned into N s 10 × 10 × 10 subdomains for the inversion procedure. A. Sensitivity Analysis
In order to analyze the behavior of the proposed inversion algorithm versus the choice of its configuration parameters, the behavior of the reconstruction error defined as Z 1 jτr ˜ − τrj ξtot dr; (9) QD D jτr 1j (τr, ˜ τr being the actual and retrieved contrast respectively, and QD the number of cells of the investigation domain) has been computed for different signal-to-noise ratios (SNRs). More specifically, a cubic scatterer with side length 0.3λ placed at the center of D and characterized by τ 1.0 has been taken as a first benchmark test case. The IMSA-SOM method has been configured assuming a maximum number of external iterations equal to S 6 and a maximum number of internal iterations equal to T s 500, s 1; …; 6. The behavior of the total error versus α under different SNR levels is plotted in Fig. 2(a). As can be observed from the above formulation, the value of α balances the relative residues in the state equation and the data equation. Accordingly, the larger the value of α the smaller the residue in the data equation. Therefore, if α is too large, the data equation residue is smaller than the noise level, and thus the relative residue in the state equation will be large, yielding to a deteriorated accuracy. On the other hand, small α values cause the residue in the data equation to be large, thus resulting in possible convergence issues. Accordingly, and as expected from SOM theory [19], the optimal α value should depend on the SNR [Fig. 2(a)]. However, since the noise level cannot be considered as a known parameter, a tradeoff choice is mandatory. The plots of the total
Fig. 3. Sensitivity analysis (centered cube, l obj 0.3λ, τ 1.0). Behavior of total error versus S when (a) SNR 30 dB, (b) SNR 20 dB, (c) SNR 15 dB, and (d) SNR 10 dB.
error versus the SNR for different values of the control threshold [Fig. 2(b)] yield the conclusion that the optimal tradeoff value for α is approximately 0.5, since such value provides the best average accuracy whatever the SNR. To complete the IMSA-SOM sensitivity analysis, its performance versus S is investigated next. The behavior of the total error versus the number of IMSA steps under different noise conditions and α configurations [Figs. 3(a)–3(d)] show that S 5 represents the best tradeoff as the total error does not change significantly after the 5th outer loop. Furthermore, the results in Fig. 3 also confirm that α 0.5 represents an optimal tradeoff choice for the SNR values (from 10 to 30 dB) considered here. B. Numerical Assessment
To assess the performance of proposed IMSA-SOM strategy in dealing with different inverse scattering scenarios, a set of benchmark targets has been imaged in the next numerical experiments. 1. “Off-centered Cube” Profile
In the first example, an off-centered cubic scatterer with side length 0.3λ and contrast τ 0.5 (Fig. 4) is assumed to be
Fig. 2. Sensitivity analysis (centered cube, l obj 0.3λ, τ 1.0). Behavior of total error (a) versus α and (b) versus the SNR.
Fig. 4. Numerical assessment (off-centered cube, l obj 0.3λ, τ 0.5). Actual scatterer.
2222
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
Research Article
Fig. 6. Numerical assessment (off-centered cube, l obj 0.3λ, τ 0.5, SNR 20 dB). Cross section of the IMSA-SOM reconstructed profiles when (a)–(c) s 4, (d)–(f ) s 5, (g)–(i) s 6, for (a), (d), and (g) x −0.3λ, (b), (e), and (h) y 0, and (c), (f ), and (i) z 0. Fig. 5. Numerical assessment (off-centered cube, l obj 0.3λ, τ 0.5, SNR 20 dB). Cross section of the actual object (a)–(c), and IMSA-SOM reconstructed profiles when (d)–(f) s 1, (g)–(i) s 2, and (j)–(l) s 3 for (a), (d), (g), and (j) x −0.3λ, (b), (e), (h), and (k) y 0, and (c), (f), (i), and (l) z 0.
y
The side lengths of each bar are l xobj 0.15λ, l obj 0.45λ, l zobj 0.15λ and the centers are located at −0.375λ; −0.075λ; 0.225λ and −0.075λ; −0.075λ; 0.225λ, respectively.
imaged with an SNR equal to 20 dB, both on the incident field and on the scattered field. The reconstruction process of the unknown scatterer is shown in Figs. 5 and 6. More in detail, the reconstruction quality enhances with the IMSA iterations from a coarse profile [S 1 − 3, Figs. 5(d)–5(f)] to an accurate inversion [S 4 − 6, Figs. 6(d)– 6(f)], where the contrast as well as the approximate shape and location of the scatterer are all clearly reconstructed. The reconstructions do not significantly improve after the 5th outer loop, as expected from the previous sensitivity analysis. However, the result is quite satisfactory considering the 20 dB noise level. To support the above considerations, a wider set of experiments has been carried out by varying the contrast in the range τ ∈ 0.5; 4.0 (Fig. 7). The plots of total errors versus contrast indicate that the IMSA-SOM is able to provide stable reconstruction for scatterers with high contrast even in high noise conditions (Fig. 7). For completeness, the evolution of the IMSA-SOM cost function is reported in Fig. 8. Note that the method rapidly converges at each IMSA step, as in the 2D case [23]. 2. “Two Parallel Bars” Profile
The previous experiments have dealt with the retrieval of single (isolated) targets. This session aims at testing the effectiveness of the IMSA-SOM in reconstructing multiple scatterers, namely, two parallel bars with τ 1 (Fig. 9).
Fig. 7. Numerical assessment (off-centered cube, l obj 0.3λ). Behavior of the total error vs. the IMSA step when (a) SNR 30 dB, (b) SNR 20 dB, (c) SNR 15 dB, (d) SNR 10 dB, and (e) SNR 5 dB.
Research Article
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
2223
Fig. 10. Numerical assessment (“two bars,” τ 1.0, SNR 20 dB). Cross section of (a)–(c) the actual object and (d)–(f) IMSA-SOM reconstructed profiles when (a) and (d) x −0.375λ, (b) and (e) x −0.075λ, and (c) and (f ) z −0.225λ.
Fig. 8. Numerical assessment (off-centered cube, l obj 0.3λ). Behavior of the cost function vs. the iteration index when (a) τ 0.5, (b) τ 1.0, (c) τ 2.0, (d) τ 3.0, and (e) τ 4.0.
function τx; y; z jxj jyj jzj 1. Therefore, the maximum value of the contrast is 2.95 while the minimum is 1.45. The plots of the retrieved two cubes when SNR 20 dB and S 3 are shown in Figs. 11 and 12, respectively. The results show that the IMSA-SOM is able to retrieve the complicated two separated inhomogeneous cubes, with the gradual changing trend of the profile clearly seen. 4. “Off-centered Ring” Profile
Fig. 9. Numerical assessment (“two bars,” τ 1.0). Actual dielectric profile.
A more complicated and challenging example is considered next (Fig. 13). More in detail, a ring-shaped scatterer centered at −0.2λ; 0.0λ; 0.0λ and with side lengths l xobj 0.3λ, y l obj 0.6λ, l zobj 0.6λ, respectively, has been assumed. The central hole has a dimension of 0.3λ × 0.3λ × 0.3λ. Such a target is much more complicated to retrieve because of the discontinuity in the permittivity profile at the center of the scatterer (i.e., the “hole”), thus yielding to high spatial frequency components in the contrast profile that are difficult to retrieve unless the probes are located very close to scatterers [24]. By comparing the results obtained by the IMSA-SOM with those of the state-of-the-art BARE-SOM [19] under the same
The plot of the retrieved dielectric profiles when SNR 20 dB (Fig. 10) shows that the presence as well as the relative permittivity of the two bars are correctly identified by the proposed methodology. Considering that the actual spacing between the two bars is only 0.15λ, it is possible to deduce that the IMSA-SOM is able to achieve super resolution even under 20 dB noise. 3. “Two Off-centered Inhomogeneous Cubes” Profile
This example aims at testing the ability of IMSA-SOM to retrieve inhomogeneous multiple objects. The domain of interest is of the size 1.5λ × 1.5λ. The scatterers to be reconstructed are two inhomogeneous cubes of side length l obj 0.5λ, with the centers located in 0.4λ; −0.4λ; 0.4λ and −0.4λ; 0.4λ; −0.4λ, respectively. The contrast of the two cubes are denoted by the
Fig. 11. Two inhomogeneous off-centered cubes of side l obj 0.5λ, SNR 20 dB—(a)–(c), Actual object #1 and (d)– (f ) IMSA-SOM reconstructed object #1 when (a) and (d) x 0.4λ, (b) and (e) y −0.4λ, and (c) and (f) z 0.4λ.
2224
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
Research Article
Fig. 12. Two inhomogeneous off-centered cubes of side l obj 0.5λ, SNR 20 dB—(a)–(c) Actual object #2 and (d)–(f ) IMSA-SOM reconstructed object #2 when (a) and (d) x 0.4λ, (b) and (e) y −0.4λ, and (c) and (f ) z 0.4λ.
Fig. 14. Numerical assessment (“off-centered ring,” τ 1.0, SNR 20 dB). Cross section of (a)–(c) the actual object and (d)–(f) of the BARE-SOM, and (g)–(i) IMSA-SOM reconstructed profiles when (a), (d), and (g) x −0.2λ, (b), (e), and (h) y 0.0, and (c), (f ), and (i) z 0.0.
Fig. 13. Numerical assessment (“off-centered ring,” τ 1.0). Actual dielectric profile.
assumptions, it turns out that the multifocusing strategy enables enhanced reconstruction accuracy [e.g., Fig. 14(d) versus Fig. 14(g)]. This result is actually expected from the theoretical viewpoint because of the fundamental features of the multiscaling scheme, which include (i) the capability to adaptively filter out the undesired artifacts with the refinement of the ROIs, and (ii) its effectiveness in balancing the available data and problem unknowns, which yields a better resolution, a stronger mitigation of the local minima issues [21], and a faster convergence. To test the performance of IMSA-SOM versus that of BARE-SOM in a more general scenario, the error figures are also evaluated for different contrast values (i.e., τ ∈ 0.5; 2— Fig. 15). As can be noticed, the total error values of the IMSA-SOM for a given SNR are significantly smaller than the corresponding BARE-SOM ones, whatever the noise level (Fig. 15). The average computational time used for one iteration is listed out in Table 1. The numerical simulations are done on a personal computer with one core of a standard I7 CPU at 2.0 G. As can be seen, the IMSA-SOM took significantly less time for a single iteration than the BARE-SOM. In addition, because of the larger number of unknowns involved,
Fig. 15. Numerical assessment (“off-centered ring”). Behavior of the total error versus the scatterer contrast.
Table 1. Comparison of Average Computational Time for Single Iteration Average Single Iteration Time(s) BARE-SOM IMSA-SOM
9.5 2.4
it takes BARE-SOM more than twice the number of iterations to converge compared to IMSA-SOM. This result further supports the previous conclusions concerning the effectiveness and reliability enabled by the proposed nested scheme, which is able to combine the advantages of the multifocusing paradigm and the subspace optimization method.
Research Article
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
2225
Future work, beyond the scope of the current paper, will be aimed at assessing the effectiveness of the proposed IMSASOM technique when dealing with anisotropic media as well as nondestructive testing/nondestructive evaluation applications [7]. Moreover, its validation against experimental data is currently under investigation. Further work will be also devoted to validating and/or extending the methodology to the case of larger-sized investigation domains, for instance considering the integration of IMSA-SOM with time-frequency hopping schemes.
Fig. 16. Numerical assessment (“off-centered ring,” τ 1.0, SNR 10 dB). Cross section of (a)–(c) IMSA-SOM and (d)–(f ) IMSA-CG reconstructed profiles when (a) and (d) x −0.2λ, (b) and (e) y 0.0, and (c) and (f) z 0.0.
4. COMPARISONS WITH STATE-OF-THE-ART MULTI-RESOLUTION TECHNIQUES The final numerical experiments are aimed at comparing the features of the proposed IMSA-SOM scheme with those of a state-of-the-art multiresolution technique [23] exploiting a local search strategy based on the conjugate gradient method (CG) for the solution of the inverse scattering equations [27]. Toward this end, the scattering scenario considered in Fig. 13 has been imaged assuming the same setup and enforcing SNR 10 dB. The plots of the retrieved dielectric profiles (Fig. 16) show that the IMSA-SOM yields a better accuracy than the IMSA-CG technique [e.g., Fig. 16(b) versus Fig. 16(e)] even in a very low SNR scenario. Such a result shows the robustness to noise of the chosen inversion procedure, which is related to the high regularization capabilities of the SOM combined with the reduction of the local minimums issues enabled by the IMSA. 5. CONCLUSIONS AND FINAL REMARKS An innovative multifocusing subspace-optimization technique is introduced to solve the 3D inverse problem by suitably generalizing, from the theoretical and algorithmic viewpoint, the 2D IMSA-SOM strategy [23]. To achieve this goal, a nested process is considered that combines a multiresolution process (which is able to iteratively refine the detected ROIs) and a subspace-optimization technique that is devoted to the retrieval of the deterministic and the ambiguous part of the induced currents. The presented numerical validation has shown that (i) suitable tradeoffs for the configuration of the IMSA-SOM control parameters can be deduced that provide satisfactory performance in many scattering scenarios and noise conditions (Section 3.A); (ii) the proposed technique outperforms its “bare” counterpart [18] in terms of efficiency and reliability thanks to the IMSA capability to mitigate local minima issues (Fig. 15); and (iii) the regularization features of SOM enable the proposed strategy to overcome state-of-the-art multiresolution techniques based on local search techniques (Fig. 16).
Funding. DIGITEO (France) under the “Call for Chairs 2014” (SIRENA project, 2014–2017); Singapore MINDEFNUS Joint Applied R&D Cooperation Programme (MINDEF/NUS/JPP/13/02/01). REFERENCES 1. Q. Liu, Z. Zhang, T. Wang, J. Bryan, G. Ybarra, L. Nolte, and W. Joines, “Active microwave imaging I-2-D forward and inverse scattering-methods,” IEEE Trans. Microw. Theory Tech. 50, 123–133 (2002). 2. A. Liseno and R. Pierri, “Imaging perfectly conducting objects as support of induced currents: Kirchhoff approximation and frequency diversity,” J. Opt. Soc. Am. A 19, 1308–1318 (2002). 3. A. Brancaccio, G. Leone, and R. Pierri, “Information content of Born scattered fields: Results in the circular cylindrical case,” J. Opt. Soc. Am. A 15, 1909–1917 (1998). 4. D. Lesselier and W. C. Chew, “Special section on electromagnetic characterization of buried obstacles,” in Inverse Problems (Institute of Physics, 2004). 5. G. Oliveri, P. P. Ding, and L. Poli, “3-D crack detection in anisotropic layered media through a sparseness-regularized solver,” IEEE Antennas Wireless Propag. Lett. 14, 1031–1034 (2015). 6. S. Caorsi, A. Massa, M. Pastorino, and M. Donelli, “Improved microwave imaging procedure for nondestructive evaluations of two dimensional structures,” IEEE Trans. Antennas Propag. 52, 1386–1397 (2004). 7. S. Caorsi, A. Massa, M. Pastorino, A. Randazzo, and A. Rosani, “A reconstruction procedure for microwave nondestructive evaluation based on a numerically computed Green’s function,” in 20th IEEE Instrumentation and Measurement Technology Conference, Vail, Colorado, 2003, pp. 987–992. 8. A. Abubakar, P. M. van den Berg, and J. J. Mallorqui, “Imaging of biomedical data using a multiplicative regularized contrast source inversion method,” IEEE Trans. Microw. Theory Techn. 50, 1761–1771 (2002). 9. W. Bangerth, A. Joshi, and E. Sevick, “Inverse biomedical imaging using separately adapted meshes for parameters and forward model variables,” in 4th IEEE International Symposium on Biomedical Imaging—Macro to Nano, Arlington, VA (2007), pp. 1368–1371. 10. A. Abubakar, W. Hu, P. M. Van den Berg, and T. M. Habashy, “A finitedifference contrast source inversion method,” Inverse Probl. 24, 065004 (2008). 11. C. Yu, L. Song, and Q. Liu, “Through-wall imaging(TWI) by radar: 2D tomographic results and analyses,” IEEE Trans. Geosci. Remote Sens. 43, 2793–2798 (2005). 12. W. Zhang and Q. Liu, “Three-dimensional scattering and inverse scattering from objects with simultaneous permittivity and permeability contrasts,” IEEE Trans. Geosci. Remote Sens. 53, 429–439 (2015). 13. T. Cui, Y. Qin, G. L. Wang, and W. Chew, “Low-frequency detection of two-dimensional buried objects using high-order extended Born approximations,” Inverse Probl. 20, S41–S62 (2004). 14. A. Devaney, “Inversion formula for inverse scattering within the Born approximation,” Opt. Lett. 7, 111–112 (1982). 15. A. Peterson, S. Ray, and R. Mittra, Computational Methods for Electromagnetics (IEEE, 1998).
2226
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
16. D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory (Springer, 1998). 17. Z. Zhang and Q. Liu, “Three-dimensional nonlinear image reconstruction for microwave biomedical Imaging,” IEEE Trans. Biomed. Eng. 51, 544–548 (2004). 18. Y. Zhong, X. Chen, and K. Agarwal, “An improved subspace-based optimization method and its implementation in solving threedimensional inverse problems,” IEEE Trans. Geosci. Remote Sens. 48, 3763–3768 (2010). 19. X. Chen, “Application of signal-subspace and optimization methods in reconstructing extended scatterers,” J. Opt. Soc. Am. A 26, 1022–1026 (2009). 20. R. Solimene, M. Maisto, and R. Pierri, “Role of diversity on the singular values of linear scattering operators: the case of strip objects,” J. Opt. Soc. Am. A 30, 2266–2272 (2013). 21. T. Isernia, V. Pascazio, and R. Pierri, “On the local minima in a tomographic imaging technique,” IEEE Trans. Geosci. Remote Sens. 39, 1596–1607 (2001).
Research Article
22. S. Caorsi, M. Donelli, D. Franceschini, and A. Massa, “A new methodology based on an iterative multiscaling for microwave imaging,” IEEE Trans. Microw. Theory Techn. 51, 1162–1173 (2003). 23. G. Oliveri, Y. Zhong, X. Chen, and A. Massa, “Multiresolution subspace-based optimization method for inverse scattering problems,” J. Opt. Soc. Am. A 28, 2057–2069 (2011). 24. O. M. Bucci and T. Isernia, “Electromagnetic inverse scattering: retrievable information and measurement strategies,” Radio Sci. 32, 2123–2137 (1997). 25. X. Chen, “Subspace-based optimization method for solving inverse scattering problems,” IEEE Trans. Geosci. Remote Sens. 48, 42–49 (2010). 26. Y. Zhong and X. Chen, “Twofold subspace-based optimization method for solving inverse scattering problems,” Inverse Probl. 25, 085003 (2009). 27. H. Harada, D. Wall, T. Takenaka, and T. Tanaka, “Conjugate gradient method applied to inverse scattering problems,” IEEE Trans. Antennas Propag. 43, 784–792 (1995).
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
Research Article
Multi-resolution subspace-based optimization method for solving three-dimensional inverse scattering problems XIUZHU YE,1 LORENZO POLI,2 GIACOMO OLIVERI,2 YU ZHONG,3 KRISHNA AGARWAL,4 ANDREA MASSA,2,5 AND XUDONG CHEN6,* 1
Department of Electrical and Computer Engineering, Beihang University, Beijing 100191, China ELEDIA Research Center@DISI, University of Trento, via Sommarive 5, 38123 Trento, Italy 3 A*STAR, Institution of High Performance Computing, Singapore 138632, Singapore 4 Singapore-MIT Alliance for Research and Technology (SMART) Centre, CREATE Tower, Singapore 138602, Singapore 5 Laboratoire des Signaux et Systèmes UMR8506 (CNRS-Univ. Paris Sud-CentraleSupélec), 3 rue Joliot Curie, 91192 Gif-sur-Yvette, France 6 Department of Electrical and Computer Engineering, National University of Singapore, Singapore 117576, Singapore *Corresponding author: [email protected] 2
Received 10 June 2015; revised 30 September 2015; accepted 30 September 2015; posted 30 September 2015 (Doc. ID 242704); published 28 October 2015
An innovative methodology is proposed to solve quantitative three-dimensional microwave imaging problems formulated within the contrast source framework. The introduced technique is based on the combination of an efficient iterative multiscaling strategy aimed at mitigating local minimum issue arising in inverse scattering problems, and a local search algorithm based on the subspace-based optimization method (SOM) devoted to effectively retrieving both the “deterministic” and the “ambiguous” parts of the unknown contrast currents. To achieve this goal, a nested iteration process is adopted in which the outer loop iteratively refines the region of interest (ROI) where the scatterers are detected, while the inner loop retrieves the dielectric properties of the scatterers within the ROIs. Selected numerical examples are also given to show the validity and robustness of the proposed algorithm in comparison with state-of-the-art techniques. © 2015 Optical Society of America OCIS codes: (290.3200) Inverse scattering; (180.6900) Three-dimensional microscopy. http://dx.doi.org/10.1364/JOSAA.32.002218
1. INTRODUCTION Active imaging methodologies [1–3] are aimed at retrieving the electromagnetic properties of the unknown scatterers located inside inaccessible domains probed by known incident waves. The importance of these techniques is motivated by their very wide domain of applications, which comprises geological exploration [4], nondestructive evaluation [5–7], biomedical [8,9], and through-wall imaging [10,11], etc. Unfortunately, solving a general inverse scattering problem is challenging due to the large number of unknowns, high nonlinearity, and ill-posedness [12–16]. Several inversion methods exploiting linear and nonlinear formulations of the imaging problem have been developed in the literature to cope with these theoretical issues. However, linear inversion methods, such as those based on the Born approximation [16], are often limited to weak scatterers. On the other hand, fully nonlinear algorithms preserve all the multiple scattering effect and are more accurate when dealing with strong scatterers [13], but they require the solution of a much more challenging inverse problem. So far, most 1084-7529/15/112218-09$15/0$15.00 © 2015 Optical Society of America
inversion methods are developed under the two-dimensional (2D) scenarios and only few nonlinear algorithms have been developed for solving three-dimensional (3D) inverse scattering problems [16]. Apart from the aforementioned difficulties, in comparison with the 2D scenario, 3D ones yield to (i) a larger amount of data and a rapid growth of the computational cost with the size of the investigation domain, (ii) a much more difficult model to be generated and visualized, and (iii) a vectorial formulation of the problem equations rather than a scalar one. However, the development of 3D imaging algorithms is of fundamental importance because of the 3D nature of most real-world problems as well as the increased information (e.g., concerning the field polarization) that can be exploited as compared to simpler 2D methodologies [17]. In this framework, a powerful nonlinear algorithm, namely the subspace-based optimization method (SOM) [18], has been proposed to image 3D scenarios within the contrast source formulation of the inverse scattering problem. This approach, which has also been widely deployed in the 2D case,
Research Article
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
decomposes the unknown contrast source into two orthogonal parts: the deterministic part, which is uniquely determined by the measured scattered field and the first several leading singular values of the mapping from domains of interest to receivers, and the ambiguous part, which is found by solving a nonlinear optimization problem through a local search strategy [19]. The singular value decomposition (SVD) procedure is an effective tool widely used in the inverse problem society and the role of diversity on the singular values of linear scattering operators has been thoroughly studied in [20]. SOM has effectively extended it to the decomposition of induced current. Such a procedure has twofold advantages. First, the deterministic part serves as a good initial guess to the optimization, which speeds up the convergence process. Second, the optimization is carried out in a smaller current “subspace,” which actually mitigates the nonlinearity of the problem [19]. Nevertheless, the SOM method still suffers from local minimums [21] and computational complexity issues when large domains or high resolutions are required, since (a) it yields to a deterministic optimization strategy in a high-dimensional space, and (b) it requires the computation of the singular value decomposition (SVD) of the Green’s operator for the evaluation of the deterministic currents. As a consequence, the integration of SOM with an effective multiresolution strategy, namely, the iterative multiscaling methods (IMSA) [22], has been proposed in [23]. Such a choice has been motivated by the fact that the IMSA can mitigate local minima issues [21], since it is able to keep the ratio between the number of unknowns and the information content of the field as close to unity as possible at each iteration. Moreover, each multizooming step requires the solution of a relatively small problem, thus reducing the computational complexity of the associated SVD [22]. However, the arising IMSA-SOM method has been developed and validated only in the 2D case [22]. Accordingly, this paper aims at generalizing the multiresolution SOM in order to yield an efficient and reliable imaging method to solve the general 3D inverse scattering problems. To validate this introduced technique, the proposed IMSA-SOM algorithm is tested in a selected set of numerical examples dealing with scatterers with different levels of complexity and/or different noise levels. The paper is organized as follows. In Section 2, the formulation of the inversion problem and the detailed IMSA-SOM procedure designed for the 3D case are introduced. Numerical examples are given in Section 3 to indicate the validity and efficiency of the proposed algorithm. Conclusions are finally drawn in Section 4. 2. PROBLEM FORMULATION AND IMSA-SOM SOLUTION APPROACH A. Formulation of the Problem
Consider a 3D setup, as shown in Fig. 1, in which an unknown scatterer (which is an isotropic one) is embedded inside a domain of interest D with lossless and nonmagnetic background of permittivity ε0 and permeability μ0, respectively. The domain D is successively illuminated by V monochromatic plane waves E¯inc p r, p 1; …; V with time dependence expjωt (omitted hereinafter for the sake of conciseness).
2219
Fig. 1. Problem geometry.
The amplitudes for the two orthogonal components of the incident electric field with respect to the polar angle and the azimuth angle are Aθ and Aϕ , respectively. For each illumination, the incident electric field E¯inc p r and the total electric field E¯tot p r are then collected by M receivers uniformly distributed in the xy, yz, xz planes around the circular observation domain, respectively, with positions rq0 , q 1; 2; …; M . In the above assumptions, the following electric field integral equation holds true: Z inc tot 2 ¯ ¯ E r E r − k0 ∇∇· gr; r 0 τr 0 E¯tot r 0 dr 0 ; D
(1) where k0 denotes the wave number in the homogenous background medium, τr εr r − j σr ωε0 − 1 is the contrast −jk jr−r 0 j
0 function of the object, and gr; r 0 e4πjr−r is the scalar 0j Green’s function. It is worth remarking that in the 3D case the operator ∇∇· is not a scalar quantity, unlike in the 2D case. In the following, Eq. (1) will be discretized by the method of moment (MoM) assuming a standard voxel basis and point matching strategy [15].
B. Proposed Inversion Methodology
In order to apply the IMSA-SOM, the guidelines presented in [23] are followed and suitably generalized. More in detail, a nested procedure is adopted that comprises the following two loops: (a) An external multiresolution loop that is devoted to adaptively identifying and zooming the region of interest (ROI) where the scatterers are detected. More specifically, at the sth step of the IMSA iteration procedure (s 1; …; S), the ROIs are meshed into N s cubic voxels, the centers of which s are located at rs n , n 1; 2; …; N , and the contrast currents s s ¯ I n , n 1; 2; …; N are assumed constant within each voxel (the N s is determined by the degrees of freedom (DOF) of the scattered field [24]). At the end of each sth step (s 1; …; S), a filtering and clustering procedure is employed to update the ROIs adaptively. (b) Internal SOM loop that comprises the computation of the deterministic currents within the ROI by means of a truncated SVD procedure aimed at calculating the deterministic current, and the evaluation of the ambiguous part of induced contrast currents through the subspace optimization strategy [18,26].
2220
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
According to the above guidelines, Eq. (1) is expressed in the observation domain in the following matrix form for the sth external iteration step: p 1; …; V ; (2) E¯sca G¯ s I¯s ; p
s
p
where E¯sca p is a 3M -dimensional vector (M being the number of field probes) containing the scattered field values for the pth incidence, and the superscript “sca” stands for the scattered s field. Moreover, I¯s p is a 3N -dimensional vector containing the unknown current coefficients. Furthermore, the discretized version of Eq. (1) can be expressed inside the investigation domain as (3) I¯s ζ¯s fE¯inc gs G¯ s I¯s ; p
p
D
p
¯s ¯s where ζ¯s diagζ¯s 1 ; …; ζ N s , ζ n is a 3 × 3 dimensional diagonal matrix with the entries on the diagonal as τn [18] s (isotropic dielectrics are assumed), and fE¯inc p g denotes the incident field vector for the pth incidence in the sth outer loop. The internal Green’s operator G¯ s D maps the contrast current to the scattered electric field inside ROI, while the external Green’s operator G¯ s s maps the induced current inside ROI to the scattered field on the receivers [26]. The detailed entries of the two mapping operators can be found in [18]. s s ¯ s s By applying the SVD to G¯ s s one yields G s · v¯m χ m u¯m , s s where the χ s m , v¯m and u¯m , m 1; …; 3M denote the singular values, right and left singular vectors, respectively. The singular values of G¯ s s can be ordered according to their magnitude (from larger to smaller) and the induced current can thus be split into two orthogonal parts, namely the deterministic and the ambiguous part, as follows: ¯det s ¯amb s I¯s p fI p g fI p g ;
(4)
s is calculated without where the deterministic part fI¯det p g resourcing to any prior information [18]: s
s fI¯det p g
LA sca X u¯s m · E¯ m1
χ s m
v¯s m:
(5)
The value of the singular value truncation numbers Ls A is selected adaptively by resourcing to a user-defined control threshold α that is obtained by the fraction of the total spectrum energy: PLs s s m1 χ m s 1; …; S; LA arg min P3M s − α ; Ls m1 χ m PLs s χ m − α > 0: (6) subject to Pm1 s 3M m1 χ m
In Eq. (6), χ s m , m 1; …; 3M denotes the singular values of the mapping G¯ s s . The ambiguous part is then constructed by employing the remaining 3N s − Ls A current basis as s fI¯amb p g
3N s −Ls A
X
m1
v¯s
mLs A
¯ s · β¯s : βs m V
(7)
Indeed, the contribution of Eq. (7) to the scattered field on the receivers is negligible [18]. It should be noted that the computational complexity of the SVD has been reduced
Research Article
in [18] since only a thin SVD is needed. This ambiguous part of the currents is then retrieved by minimizing the following cost function [18]: F s ζ¯s ; β¯s ; …; β¯s V
1
0
1 ¯s ¯sca 2 ‖G¯ s s Ip − Ep ‖ 2 V B C X ‖E¯sca p ‖ B C B s s C; s ¯s 2 A ¯ ¯ inc s @ ¯ ¯ ‖I p − ζ fE p g G D I p ‖ p1
(8)
s 2 ‖fI det p g ‖
where ‖ ‖ is the Euclidean norm through a local search strategy as detailed in [18]. The arising IMSA-SOM iteration process for the 3D scenario can be summarized as follows (s being the number of outer iteration steps and t being the number of inner iteration steps): Step 1: Outer loop: initialization s 1. s
s
s
Step 1.1: Calculate G D , G s and the SVD of G s . Obtain s fI¯det p g , p 1; 2; …; V by using Eq. (5). Step 1.2: Initialization: t 0, set a proper initial guess to the contrast τ¯s 0 0 [25], and to accelerate the convergence process we obtain β¯s 0 by back propagation, which is different from [25]. Step 1.3: Two step conjugate-gradient type optimization [27]: t t 1. Step 1.3.1: Update β¯s t . Step 1.3.2: Update τ¯s t . Step 1.4: Stop iteration when a user-defined maximum number of iterations T s is reached or a user-defined threshold is reached. Step 2: Stop iteration if a user-defined outer loop number S is reached, otherwise identify/update the ROIs: s s 1 and go to step 1. The details of the two step conjugate-gradient method are omitted for concision and we refer the readers to [22] and [27] for details. 3. NUMERICAL EXAMPLES The numerical analysis presented in this section has a twofold objective. First, it aims to analyze the sensitivity of the proposed approach on its control parameters (i.e., the fractional parameter α and the number of IMSA steps S) and to provide a guideline on their choice. Second, it is devoted to the assessment of the reconstruction accuracy and efficiency as well as the robustness of the IMSA-SOM when dealing with various scatterer and noise levels. Comparisons between the IMSA-SOM with the BARE-SOM as well as to state-of-the-art multiresolution techniques are also provided to show its advantages and drawbacks. Toward this end, a benchmark cubic investigation domain D with side length 1.5λ is assumed to be illuminated by V 30 plane waves with incident angles uniformly distributed around a circle in the y–z plane. The amplitudes for each component of the incident electric field are Aθ 0.7 and Aϕ 0.7, respectively. For each incidence angle, M 30 probes are evenly distributed on circles of radius 3.2λ in the xy, yz, xz planes, respectively, to measure the field data.
Research Article
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
2221
Synthetic data are generated by dividing the domain of interest D into 20 × 20 × 20 voxels. To avoid inverse crime issues, the domain is partitioned into N s 10 × 10 × 10 subdomains for the inversion procedure. A. Sensitivity Analysis
In order to analyze the behavior of the proposed inversion algorithm versus the choice of its configuration parameters, the behavior of the reconstruction error defined as Z 1 jτr ˜ − τrj ξtot dr; (9) QD D jτr 1j (τr, ˜ τr being the actual and retrieved contrast respectively, and QD the number of cells of the investigation domain) has been computed for different signal-to-noise ratios (SNRs). More specifically, a cubic scatterer with side length 0.3λ placed at the center of D and characterized by τ 1.0 has been taken as a first benchmark test case. The IMSA-SOM method has been configured assuming a maximum number of external iterations equal to S 6 and a maximum number of internal iterations equal to T s 500, s 1; …; 6. The behavior of the total error versus α under different SNR levels is plotted in Fig. 2(a). As can be observed from the above formulation, the value of α balances the relative residues in the state equation and the data equation. Accordingly, the larger the value of α the smaller the residue in the data equation. Therefore, if α is too large, the data equation residue is smaller than the noise level, and thus the relative residue in the state equation will be large, yielding to a deteriorated accuracy. On the other hand, small α values cause the residue in the data equation to be large, thus resulting in possible convergence issues. Accordingly, and as expected from SOM theory [19], the optimal α value should depend on the SNR [Fig. 2(a)]. However, since the noise level cannot be considered as a known parameter, a tradeoff choice is mandatory. The plots of the total
Fig. 3. Sensitivity analysis (centered cube, l obj 0.3λ, τ 1.0). Behavior of total error versus S when (a) SNR 30 dB, (b) SNR 20 dB, (c) SNR 15 dB, and (d) SNR 10 dB.
error versus the SNR for different values of the control threshold [Fig. 2(b)] yield the conclusion that the optimal tradeoff value for α is approximately 0.5, since such value provides the best average accuracy whatever the SNR. To complete the IMSA-SOM sensitivity analysis, its performance versus S is investigated next. The behavior of the total error versus the number of IMSA steps under different noise conditions and α configurations [Figs. 3(a)–3(d)] show that S 5 represents the best tradeoff as the total error does not change significantly after the 5th outer loop. Furthermore, the results in Fig. 3 also confirm that α 0.5 represents an optimal tradeoff choice for the SNR values (from 10 to 30 dB) considered here. B. Numerical Assessment
To assess the performance of proposed IMSA-SOM strategy in dealing with different inverse scattering scenarios, a set of benchmark targets has been imaged in the next numerical experiments. 1. “Off-centered Cube” Profile
In the first example, an off-centered cubic scatterer with side length 0.3λ and contrast τ 0.5 (Fig. 4) is assumed to be
Fig. 2. Sensitivity analysis (centered cube, l obj 0.3λ, τ 1.0). Behavior of total error (a) versus α and (b) versus the SNR.
Fig. 4. Numerical assessment (off-centered cube, l obj 0.3λ, τ 0.5). Actual scatterer.
2222
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
Research Article
Fig. 6. Numerical assessment (off-centered cube, l obj 0.3λ, τ 0.5, SNR 20 dB). Cross section of the IMSA-SOM reconstructed profiles when (a)–(c) s 4, (d)–(f ) s 5, (g)–(i) s 6, for (a), (d), and (g) x −0.3λ, (b), (e), and (h) y 0, and (c), (f ), and (i) z 0. Fig. 5. Numerical assessment (off-centered cube, l obj 0.3λ, τ 0.5, SNR 20 dB). Cross section of the actual object (a)–(c), and IMSA-SOM reconstructed profiles when (d)–(f) s 1, (g)–(i) s 2, and (j)–(l) s 3 for (a), (d), (g), and (j) x −0.3λ, (b), (e), (h), and (k) y 0, and (c), (f), (i), and (l) z 0.
y
The side lengths of each bar are l xobj 0.15λ, l obj 0.45λ, l zobj 0.15λ and the centers are located at −0.375λ; −0.075λ; 0.225λ and −0.075λ; −0.075λ; 0.225λ, respectively.
imaged with an SNR equal to 20 dB, both on the incident field and on the scattered field. The reconstruction process of the unknown scatterer is shown in Figs. 5 and 6. More in detail, the reconstruction quality enhances with the IMSA iterations from a coarse profile [S 1 − 3, Figs. 5(d)–5(f)] to an accurate inversion [S 4 − 6, Figs. 6(d)– 6(f)], where the contrast as well as the approximate shape and location of the scatterer are all clearly reconstructed. The reconstructions do not significantly improve after the 5th outer loop, as expected from the previous sensitivity analysis. However, the result is quite satisfactory considering the 20 dB noise level. To support the above considerations, a wider set of experiments has been carried out by varying the contrast in the range τ ∈ 0.5; 4.0 (Fig. 7). The plots of total errors versus contrast indicate that the IMSA-SOM is able to provide stable reconstruction for scatterers with high contrast even in high noise conditions (Fig. 7). For completeness, the evolution of the IMSA-SOM cost function is reported in Fig. 8. Note that the method rapidly converges at each IMSA step, as in the 2D case [23]. 2. “Two Parallel Bars” Profile
The previous experiments have dealt with the retrieval of single (isolated) targets. This session aims at testing the effectiveness of the IMSA-SOM in reconstructing multiple scatterers, namely, two parallel bars with τ 1 (Fig. 9).
Fig. 7. Numerical assessment (off-centered cube, l obj 0.3λ). Behavior of the total error vs. the IMSA step when (a) SNR 30 dB, (b) SNR 20 dB, (c) SNR 15 dB, (d) SNR 10 dB, and (e) SNR 5 dB.
Research Article
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
2223
Fig. 10. Numerical assessment (“two bars,” τ 1.0, SNR 20 dB). Cross section of (a)–(c) the actual object and (d)–(f) IMSA-SOM reconstructed profiles when (a) and (d) x −0.375λ, (b) and (e) x −0.075λ, and (c) and (f ) z −0.225λ.
Fig. 8. Numerical assessment (off-centered cube, l obj 0.3λ). Behavior of the cost function vs. the iteration index when (a) τ 0.5, (b) τ 1.0, (c) τ 2.0, (d) τ 3.0, and (e) τ 4.0.
function τx; y; z jxj jyj jzj 1. Therefore, the maximum value of the contrast is 2.95 while the minimum is 1.45. The plots of the retrieved two cubes when SNR 20 dB and S 3 are shown in Figs. 11 and 12, respectively. The results show that the IMSA-SOM is able to retrieve the complicated two separated inhomogeneous cubes, with the gradual changing trend of the profile clearly seen. 4. “Off-centered Ring” Profile
Fig. 9. Numerical assessment (“two bars,” τ 1.0). Actual dielectric profile.
A more complicated and challenging example is considered next (Fig. 13). More in detail, a ring-shaped scatterer centered at −0.2λ; 0.0λ; 0.0λ and with side lengths l xobj 0.3λ, y l obj 0.6λ, l zobj 0.6λ, respectively, has been assumed. The central hole has a dimension of 0.3λ × 0.3λ × 0.3λ. Such a target is much more complicated to retrieve because of the discontinuity in the permittivity profile at the center of the scatterer (i.e., the “hole”), thus yielding to high spatial frequency components in the contrast profile that are difficult to retrieve unless the probes are located very close to scatterers [24]. By comparing the results obtained by the IMSA-SOM with those of the state-of-the-art BARE-SOM [19] under the same
The plot of the retrieved dielectric profiles when SNR 20 dB (Fig. 10) shows that the presence as well as the relative permittivity of the two bars are correctly identified by the proposed methodology. Considering that the actual spacing between the two bars is only 0.15λ, it is possible to deduce that the IMSA-SOM is able to achieve super resolution even under 20 dB noise. 3. “Two Off-centered Inhomogeneous Cubes” Profile
This example aims at testing the ability of IMSA-SOM to retrieve inhomogeneous multiple objects. The domain of interest is of the size 1.5λ × 1.5λ. The scatterers to be reconstructed are two inhomogeneous cubes of side length l obj 0.5λ, with the centers located in 0.4λ; −0.4λ; 0.4λ and −0.4λ; 0.4λ; −0.4λ, respectively. The contrast of the two cubes are denoted by the
Fig. 11. Two inhomogeneous off-centered cubes of side l obj 0.5λ, SNR 20 dB—(a)–(c), Actual object #1 and (d)– (f ) IMSA-SOM reconstructed object #1 when (a) and (d) x 0.4λ, (b) and (e) y −0.4λ, and (c) and (f) z 0.4λ.
2224
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
Research Article
Fig. 12. Two inhomogeneous off-centered cubes of side l obj 0.5λ, SNR 20 dB—(a)–(c) Actual object #2 and (d)–(f ) IMSA-SOM reconstructed object #2 when (a) and (d) x 0.4λ, (b) and (e) y −0.4λ, and (c) and (f ) z 0.4λ.
Fig. 14. Numerical assessment (“off-centered ring,” τ 1.0, SNR 20 dB). Cross section of (a)–(c) the actual object and (d)–(f) of the BARE-SOM, and (g)–(i) IMSA-SOM reconstructed profiles when (a), (d), and (g) x −0.2λ, (b), (e), and (h) y 0.0, and (c), (f ), and (i) z 0.0.
Fig. 13. Numerical assessment (“off-centered ring,” τ 1.0). Actual dielectric profile.
assumptions, it turns out that the multifocusing strategy enables enhanced reconstruction accuracy [e.g., Fig. 14(d) versus Fig. 14(g)]. This result is actually expected from the theoretical viewpoint because of the fundamental features of the multiscaling scheme, which include (i) the capability to adaptively filter out the undesired artifacts with the refinement of the ROIs, and (ii) its effectiveness in balancing the available data and problem unknowns, which yields a better resolution, a stronger mitigation of the local minima issues [21], and a faster convergence. To test the performance of IMSA-SOM versus that of BARE-SOM in a more general scenario, the error figures are also evaluated for different contrast values (i.e., τ ∈ 0.5; 2— Fig. 15). As can be noticed, the total error values of the IMSA-SOM for a given SNR are significantly smaller than the corresponding BARE-SOM ones, whatever the noise level (Fig. 15). The average computational time used for one iteration is listed out in Table 1. The numerical simulations are done on a personal computer with one core of a standard I7 CPU at 2.0 G. As can be seen, the IMSA-SOM took significantly less time for a single iteration than the BARE-SOM. In addition, because of the larger number of unknowns involved,
Fig. 15. Numerical assessment (“off-centered ring”). Behavior of the total error versus the scatterer contrast.
Table 1. Comparison of Average Computational Time for Single Iteration Average Single Iteration Time(s) BARE-SOM IMSA-SOM
9.5 2.4
it takes BARE-SOM more than twice the number of iterations to converge compared to IMSA-SOM. This result further supports the previous conclusions concerning the effectiveness and reliability enabled by the proposed nested scheme, which is able to combine the advantages of the multifocusing paradigm and the subspace optimization method.
Research Article
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
2225
Future work, beyond the scope of the current paper, will be aimed at assessing the effectiveness of the proposed IMSASOM technique when dealing with anisotropic media as well as nondestructive testing/nondestructive evaluation applications [7]. Moreover, its validation against experimental data is currently under investigation. Further work will be also devoted to validating and/or extending the methodology to the case of larger-sized investigation domains, for instance considering the integration of IMSA-SOM with time-frequency hopping schemes.
Fig. 16. Numerical assessment (“off-centered ring,” τ 1.0, SNR 10 dB). Cross section of (a)–(c) IMSA-SOM and (d)–(f ) IMSA-CG reconstructed profiles when (a) and (d) x −0.2λ, (b) and (e) y 0.0, and (c) and (f) z 0.0.
4. COMPARISONS WITH STATE-OF-THE-ART MULTI-RESOLUTION TECHNIQUES The final numerical experiments are aimed at comparing the features of the proposed IMSA-SOM scheme with those of a state-of-the-art multiresolution technique [23] exploiting a local search strategy based on the conjugate gradient method (CG) for the solution of the inverse scattering equations [27]. Toward this end, the scattering scenario considered in Fig. 13 has been imaged assuming the same setup and enforcing SNR 10 dB. The plots of the retrieved dielectric profiles (Fig. 16) show that the IMSA-SOM yields a better accuracy than the IMSA-CG technique [e.g., Fig. 16(b) versus Fig. 16(e)] even in a very low SNR scenario. Such a result shows the robustness to noise of the chosen inversion procedure, which is related to the high regularization capabilities of the SOM combined with the reduction of the local minimums issues enabled by the IMSA. 5. CONCLUSIONS AND FINAL REMARKS An innovative multifocusing subspace-optimization technique is introduced to solve the 3D inverse problem by suitably generalizing, from the theoretical and algorithmic viewpoint, the 2D IMSA-SOM strategy [23]. To achieve this goal, a nested process is considered that combines a multiresolution process (which is able to iteratively refine the detected ROIs) and a subspace-optimization technique that is devoted to the retrieval of the deterministic and the ambiguous part of the induced currents. The presented numerical validation has shown that (i) suitable tradeoffs for the configuration of the IMSA-SOM control parameters can be deduced that provide satisfactory performance in many scattering scenarios and noise conditions (Section 3.A); (ii) the proposed technique outperforms its “bare” counterpart [18] in terms of efficiency and reliability thanks to the IMSA capability to mitigate local minima issues (Fig. 15); and (iii) the regularization features of SOM enable the proposed strategy to overcome state-of-the-art multiresolution techniques based on local search techniques (Fig. 16).
Funding. DIGITEO (France) under the “Call for Chairs 2014” (SIRENA project, 2014–2017); Singapore MINDEFNUS Joint Applied R&D Cooperation Programme (MINDEF/NUS/JPP/13/02/01). REFERENCES 1. Q. Liu, Z. Zhang, T. Wang, J. Bryan, G. Ybarra, L. Nolte, and W. Joines, “Active microwave imaging I-2-D forward and inverse scattering-methods,” IEEE Trans. Microw. Theory Tech. 50, 123–133 (2002). 2. A. Liseno and R. Pierri, “Imaging perfectly conducting objects as support of induced currents: Kirchhoff approximation and frequency diversity,” J. Opt. Soc. Am. A 19, 1308–1318 (2002). 3. A. Brancaccio, G. Leone, and R. Pierri, “Information content of Born scattered fields: Results in the circular cylindrical case,” J. Opt. Soc. Am. A 15, 1909–1917 (1998). 4. D. Lesselier and W. C. Chew, “Special section on electromagnetic characterization of buried obstacles,” in Inverse Problems (Institute of Physics, 2004). 5. G. Oliveri, P. P. Ding, and L. Poli, “3-D crack detection in anisotropic layered media through a sparseness-regularized solver,” IEEE Antennas Wireless Propag. Lett. 14, 1031–1034 (2015). 6. S. Caorsi, A. Massa, M. Pastorino, and M. Donelli, “Improved microwave imaging procedure for nondestructive evaluations of two dimensional structures,” IEEE Trans. Antennas Propag. 52, 1386–1397 (2004). 7. S. Caorsi, A. Massa, M. Pastorino, A. Randazzo, and A. Rosani, “A reconstruction procedure for microwave nondestructive evaluation based on a numerically computed Green’s function,” in 20th IEEE Instrumentation and Measurement Technology Conference, Vail, Colorado, 2003, pp. 987–992. 8. A. Abubakar, P. M. van den Berg, and J. J. Mallorqui, “Imaging of biomedical data using a multiplicative regularized contrast source inversion method,” IEEE Trans. Microw. Theory Techn. 50, 1761–1771 (2002). 9. W. Bangerth, A. Joshi, and E. Sevick, “Inverse biomedical imaging using separately adapted meshes for parameters and forward model variables,” in 4th IEEE International Symposium on Biomedical Imaging—Macro to Nano, Arlington, VA (2007), pp. 1368–1371. 10. A. Abubakar, W. Hu, P. M. Van den Berg, and T. M. Habashy, “A finitedifference contrast source inversion method,” Inverse Probl. 24, 065004 (2008). 11. C. Yu, L. Song, and Q. Liu, “Through-wall imaging(TWI) by radar: 2D tomographic results and analyses,” IEEE Trans. Geosci. Remote Sens. 43, 2793–2798 (2005). 12. W. Zhang and Q. Liu, “Three-dimensional scattering and inverse scattering from objects with simultaneous permittivity and permeability contrasts,” IEEE Trans. Geosci. Remote Sens. 53, 429–439 (2015). 13. T. Cui, Y. Qin, G. L. Wang, and W. Chew, “Low-frequency detection of two-dimensional buried objects using high-order extended Born approximations,” Inverse Probl. 20, S41–S62 (2004). 14. A. Devaney, “Inversion formula for inverse scattering within the Born approximation,” Opt. Lett. 7, 111–112 (1982). 15. A. Peterson, S. Ray, and R. Mittra, Computational Methods for Electromagnetics (IEEE, 1998).
2226
Vol. 32, No. 11 / November 2015 / Journal of the Optical Society of America A
16. D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory (Springer, 1998). 17. Z. Zhang and Q. Liu, “Three-dimensional nonlinear image reconstruction for microwave biomedical Imaging,” IEEE Trans. Biomed. Eng. 51, 544–548 (2004). 18. Y. Zhong, X. Chen, and K. Agarwal, “An improved subspace-based optimization method and its implementation in solving threedimensional inverse problems,” IEEE Trans. Geosci. Remote Sens. 48, 3763–3768 (2010). 19. X. Chen, “Application of signal-subspace and optimization methods in reconstructing extended scatterers,” J. Opt. Soc. Am. A 26, 1022–1026 (2009). 20. R. Solimene, M. Maisto, and R. Pierri, “Role of diversity on the singular values of linear scattering operators: the case of strip objects,” J. Opt. Soc. Am. A 30, 2266–2272 (2013). 21. T. Isernia, V. Pascazio, and R. Pierri, “On the local minima in a tomographic imaging technique,” IEEE Trans. Geosci. Remote Sens. 39, 1596–1607 (2001).
Research Article
22. S. Caorsi, M. Donelli, D. Franceschini, and A. Massa, “A new methodology based on an iterative multiscaling for microwave imaging,” IEEE Trans. Microw. Theory Techn. 51, 1162–1173 (2003). 23. G. Oliveri, Y. Zhong, X. Chen, and A. Massa, “Multiresolution subspace-based optimization method for inverse scattering problems,” J. Opt. Soc. Am. A 28, 2057–2069 (2011). 24. O. M. Bucci and T. Isernia, “Electromagnetic inverse scattering: retrievable information and measurement strategies,” Radio Sci. 32, 2123–2137 (1997). 25. X. Chen, “Subspace-based optimization method for solving inverse scattering problems,” IEEE Trans. Geosci. Remote Sens. 48, 42–49 (2010). 26. Y. Zhong and X. Chen, “Twofold subspace-based optimization method for solving inverse scattering problems,” Inverse Probl. 25, 085003 (2009). 27. H. Harada, D. Wall, T. Takenaka, and T. Tanaka, “Conjugate gradient method applied to inverse scattering problems,” IEEE Trans. Antennas Propag. 43, 784–792 (1995).
