A multipole-expansion based linear sampling method for solving inverse scattering problems

Krishna Agarwal; Xudong Chen; Yu Zhong

doi:10.1364/OE.18.006366

1. Introduction

Non-destructive evaluation of inaccessible regions forms an important class of problems of practical interest. A detailed quantitative analysis of the region often requires solving severely ill-posed problems and may not be necessary in some applications. Most often, knowledge of the presence, locations and shapes of scatterers is sufficient. Further, for problems requiring complete inverse solution, the information mentioned above may serve as a good initial guess and reduce the complexity of inverse problems considerably [1–6]. Several methods have been proposed for such qualitative imaging of the domain. Some examples include multiple signal classification (MUSIC) [2,3,7,8], decomposition of time reversal operator (DORT) [9–11], synthetic aperture focusing technique (SAFT) [12,13], factorization [14–16], and linear sampling method (LSM) [17–25].

The focus of the present work is to propose a variant of the linear sampling method based upon the multipole expansion of the scattered field. LSM has been proposed and understood primarily as a mathematical approach [17,18]. Following [17,18], some important work in developing LSM from the mathematical perspective has been done [19–23]. Only, recently, some research papers have attempted to provide an insight into LSM from the perspective of physics [24–26]. These papers interpret LSM as a focusing problem and study the distribution of induced current on the scatterer. Despite different interpretations, conventionally some form of mathematical regularization is employed for LSM. As far as we know, Tikhonov regularization is used almost ubiquitously in LSM. Though other mathematical regularization techniques may be used, this is not the scope of the paper.

This paper presents a modified linear sampling method that uses a physical regularization for reconstruction of the scatterer support. To this end, we study the multipole expansion of the scattered field instead of the distribution of the current. Out of the infinite number of multipoles possible, the proposed method uses only the monopole and dipole terms for reconstruction and truncates all the higher order multipoles. Thus, in essence, the proposed method uses a physical approximation instead of the conventionally used mathematical regularization. The truncation of the higher order multipoles is justified using various arguments, including the effect of noise on such regularization. It is also shown that as compared to the conventional LSM that is sensitive to the choice of regularization parameter and the threshold used for deciding the scatterer support, the proposed method is more robust to such variations. Various examples are presented to validate the performance of the proposed method.

The original contributions of the paper are in the following aspects: First, the paper proposes a novel and interesting physical explanation to the LSM. Since LSM was primarily developed by mathematicians, understandably they resorted to mathematical regularization methods rather than physically grounded regularizations. Given that LSM is being accepted very well in the inverse scattering community, it is indeed beneficial that LSM be understood and used in terms of scattering physics. Second, the proposed modified version of LSM (MLSM) is simple to implement and obtains better imaging results than the traditional LSM for scatterers that are not simply connected. The regularization parameter in MLSM involves the choice of multipoles, which is more intuitive and simple than choosing a mathematical constant for the regularization parameter and more computation efficient than applying the general discrepancy rule for generating the regularization parameter as in traditional LSM. Third, the proposed work not only investigates the LSM in a new perspective, it also provides the linkage of the LSM to other imaging methods. For example, the suppression of higher multipoles in MLSM is in spirit similar to the suppression of secondary sources induced at other point-like scatterers in the MUSIC imaging [8]. The proposed work hopefully opens a new research direction that finds intimate links between the LSM and other imaging methods.

2. Background: Linear sampling method and multipole expansion

2.1. Linear sampling method

The two-dimensional scalar electromagnetic scattering problem (i.e., infinitely extending dielectric cylinders under TM illumination) is considered here. We assume that the support of dielectric scatterers is ∑, and its boundary is denoted by ∂∑. Let E _∞(φ, θ) be the scattered far-field pattern measured on Γ in the direction of φ, when a unit amplitude plane wave impinges from the direction θ, where Γ is a circle centered at the origin and has sufficiently large radius compared with the wavelength λ. The two-dimensional Green’s function is denoted by Φ(r⃗,r⃗́) = (i/4)H₀ ⁽¹⁾(k|r⃗-r⃗́|), where k is the wave number. Its far field expression is denoted as Φ_∞(φ, r⃗́), where φ is the angle of point r⃗.

For a generic point r⃗_p, LSM consists of solving the far-field integral equation for the unknown g(θ,r⃗_p) [17–25]:

\int_{Γ} E_{\infty} (φ, θ) g (θ, {\vec{r}}_{p}) dθ = Φ_{\infty} (φ, {\vec{r}}_{p}) .

It is obvious that the right hand side represents a circularly symmetric scattered field with respect to r⃗_p, which, for later convenience, is referred to as the monopole radiation pattern. The value of g(θ,r⃗_p) becomes unbounded if the sampling point r⃗_p does not belong to the scatterer support. While reconstructing, (1) is cast into a matrix equation and the matrix equation is solved for g(θ,r⃗_p) using Tikhonov regularization. The sampling points r⃗_p with corresponding g(θ,r⃗_p) significantly large are concluded to belong to the background. The support of scatterer is identified as those of r⃗_p with corresponding g(θ,r⃗_p) significantly small.

2.2. Linear sampling method as a focusing approach

In the LSM, the right-hand-side of (1) represents a circularly symmetric far field pattern. However, this pattern is not necessarily achieved by an elementary point source or a circularly symmetric current around the sampling point. The reason for this argument is discussed below. It is highlighted that this argument motivates us to investigate the multipole expansion of the far field (see section 2.4) instead of the current distribution, which eventually leads to the modified LSM that is to be discussed in section 3.

First, we derive a current distribution in the free space that is able to produce the far field radiation pattern Φ_∞(φ,r⃗_p). Let Ω be any arbitrary smooth domain containing the sampling point r⃗_p, i.e., r⃗_p ∈ Ω. An arbitrary function u(r⃗) ∈ C ²(Ω¯) satisfying the following boundary conditions is chosen:

\begin{matrix} u (\vec{r}) = Φ (\vec{r}, {\vec{r}}_{p}) \\ \frac{\partial u (\vec{r})}{\partial n} = \frac{\partial Φ (\vec{r}, {\vec{r}}_{p})}{\partial n} \end{matrix}} for \vec{r} \in \partial Ω,

where n denotes the normal to the boundary pointing outward. We set the current distribution in Ω as J(r⃗) = −∇² u(r⃗) − k ² u(r⃗). For r⃗ ∈ Ω, using the Green’s formula [27],

u (\vec{r}) = \int_{\partial Ω} [Φ (\vec{r}, {\vec{r}}^{'}) \frac{\partial u}{\partial n} ({\vec{r}}^{'}) - u ({\vec{r}}^{'}) \frac{\partial Φ (\vec{r}, {\vec{r}}^{'})}{\partial n ({\vec{r}}^{'})}] ds ({\vec{r}}^{'}) - \iint_{Ω} [\nabla^{2} u ({\vec{r}}^{'}) + k^{2} u ({\vec{r}}^{'})] Φ (\vec{r}, {\vec{r}}^{'}) d {\vec{r}}^{'},

and substituting (2):

u (\vec{r}) = \int_{\partial Ω} [Φ (\vec{r}, {\vec{r}}^{'}) \frac{\partial Φ ({\vec{r}}^{'}, {\vec{r}}_{p})}{\partial n ({\vec{r}}^{'})} - Φ ({\vec{r}}^{'}, {\vec{r}}_{p}) \frac{\partial Φ (\vec{r}, {\vec{r}}^{'})}{\partial n ({\vec{r}}^{'})}] ds ({\vec{r}}^{'}) + \iint_{Ω} J ({\vec{r}}^{'}) Φ (\vec{r}, {\vec{r}}^{'}) d {\vec{r}}^{'} .

It is easy to prove that the first integral vanishes by applying Green’s second theorem [27] outside of Ω,

\int_{\partial Ω + Γ_{\infty}} [Φ (\vec{r}, {\vec{r}}^{'}) \frac{\partial Φ ({\vec{r}}^{'}, {\vec{r}}_{p})}{\partial n ({\vec{r}}^{'})} - Φ ({\vec{r}}^{'}, {\vec{r}}_{p}) \frac{\partial Φ (\vec{r}, {\vec{r}}^{'})}{\partial n ({\vec{r}}^{'})}] ds ({\vec{r}}^{'}),

considering the fact that the first integral over the infinitely far surface Γ_∞ is zero. Thus, we have:

u (\vec{r}) = \iint_{Ω} J ({\vec{r}}^{'}) Φ (\vec{r}, {\vec{r}}^{'}) d {\vec{r}}^{'}, \vec{r} \in Ω .

This definition of u(r⃗) can be extended to the domain outside Ω, i.e., for r⃗ ∈ ℝ²\Ω. Thus, u(r⃗) satisfies the radiation condition, coincides with u(r⃗) = Φ(r⃗,r⃗_p) on r⃗ ∈ ∂Ω, and satisfies the Helmholtz equation in ℝ²\Ω. Considering the uniqueness of the exterior Dirichlet problem [27], it can be concluded that u(r⃗) = Φ(r⃗,r⃗_p) in ℝ²\Ω. Consequently, the current distribution J(r⃗) is able to produce the far field pattern Φ_∞(φ, r⃗_p). Due to the arbitrariness of both the domain Ω and the function u(r⃗) ∈ C ²(Ω¯), so long as u(r⃗) satisfies (2), the current distribution J(r⃗) that produces the far field pattern Φ_∞(φ, r⃗_p) is not necessarily a circularly symmetric source centered in r⃗_p or an elementary source located at r⃗_p.

Next, we consider a current distribution in the scatterer (r⃗ ∈ ∑) that is able to produce the far field radiation pattern Φ_∞(φ, r⃗_p). Let Ω in the previous paragraph be ∑, and the sought current distribution be J(r⃗) = −∇² u(r⃗)−k ² u(r⃗). Then following the above argument, the induced current distribution J(r⃗ ∈ ∑) need not necessarily be a circularly symmetric source centered in r⃗_p or an elementary source located at r⃗_p to produce Φ_∞(φ, r⃗_p).

2.3. Example

Let us consider the scatterer shown in Fig. 1(a). In numerical simulations, 13 detectors and 13 sources, distributed uniformly along Γ with radius 10λ, have been used for reconstruction and the LSM reconstruction result thus obtained is shown in Fig. 1(b). Here, −log₁₀∥g(θ,r⃗_p)∥ has been plotted for various sampling points. Considering the sampling point r⃗_p =(0,0), which LSM detects as belonging to the scatterer, the induced current distribution on the scatterer is plotted in Fig. 1(c). It is noticeable that neither there is a significantly strong elementary source at r⃗_p, nor is the current distribution circularly symmetric.

Fig. 1. Illustration of the current distribution for the sampling point r⃗_p = (0,0) and the proposed multipole-based interpretation.

Download Full Size | PDF

2.4. Linear sampling method: expansion of scattered far field in terms of multipole radiation

As shown in the previous sub-sections, the monopole radiation pattern is not necessarily produced by a point source or circularly symmetric distributed source. This indicates that instead of investigating current distribution, we should investigate the radiation field itself. The scattered fields received at the receivers can be decomposed into various independent terms corresponding to the multipole radiation from a sampling point r⃗_p. In fact, to solve (1) is to find a superposition, g(θ,r⃗_p), of E _∞(φ,θ) such that among all multipoles, the monopole radiation is the only dominant contributor in the resultant total radiation.

The scattered field received at a point r⃗_φ, E(r⃗_φ,θ), whose far field notation is E_φ(φ,θ), is radiated by the induced currents on the scatterer, i.e.,

E ({\vec{r}}_{φ}, θ) = \int_{\sum} J (\vec{r}, θ) Φ ({\vec{r}}_{φ}, \vec{r}) d \vec{r},

where, J(r⃗,θ) is the current induced on r⃗. Using addition theorem [28] on Φ(r⃗_φ,r⃗), the expression of E(r⃗_φ,θ) can be rewritten in terms of various multipoles corresponding to a sampling point r⃗_p as below:

E ({\vec{r}}_{φ}, θ) = \sum_{n = - \infty}^{\infty} α^{(n)} ({\vec{r}}_{p}, θ) Φ^{(n)} ({\vec{r}}_{φ}, {\vec{r}}_{p}),

where,

α^{(n)} ({\vec{r}}_{p}, θ) = \int_{\sum} J (\vec{r}, θ) J_{n} (k ∣ {\vec{r}}_{p} - \vec{r} ∣) e^{- in \arg (\vec{r} - {\vec{r}}_{p})} d \vec{r},

Φ^{(n)} ({\vec{r}}_{φ}, {\vec{r}}_{p}) = \frac{i}{4} H_{n}^{(1)} (k ∣ {\vec{r}}_{φ} - {\vec{r}}_{p} ∣) e^{in \arg ({\vec{r}}_{φ} - {\vec{r}}_{p})} .

It is evident that α ⁽ⁿ⁾(r⃗_p, θ) represents the n th effective multipole current at r⃗_p, such that the sum of the radiated fields from all such multipoles is equal to the measured scattered field. From (8), the fundamental equation of LSM (1) is equivalent to:

\int_{Γ} α^{(n)} ({\vec{r}}_{p}, θ) g (θ, {\vec{r}}_{p}) dθ = {\begin{matrix} 1 n = 0 \\ 0 otherwise \end{matrix} .

For convenience of further use, we define:

υ_{n} = \int_{Γ} α^{(n)} ({\vec{r}}_{p}, θ) g (θ, {\vec{r}}_{p}) d θ .

Physically, υ_n can be understood as the multipole currents induced by the source distribution g(θ,r⃗_p). It is shown that (11) holds (meaning that monopole is the single dominant contributor of radiation from the induced sources) even in an arbitrarily shaped support. This means that as long as the multipole expansion of the induced current distribution at a sampling point is such that the monopole is the only prominent contributor, the sampling point will be detected as a scatterer.

For the example in section 2.3, we use 41 multipoles and consider the same sampling point as before, i.e., r⃗_p = (0,0). The values of α ⁽ⁿ⁾(r⃗_p, θ) are determined analytically using (9). The error in the multipole expansion $\frac{∥ E_{\infty} (φ, θ) \sum_{n = - 20}^{20} α^{(n)} ({\vec{r}}_{p}, θ) Φ^{(n)} (φ, {\vec{r}}_{p}) ∥}{∥ E_{\infty} (φ, θ) ∥}$ is approximately 10^-14. Thus, the multipole expansion is reasonably correct. The values obtained for υ_n (see (12)), by substituting α ⁽ⁿ⁾(r⃗_p,θ) and g(θ,r⃗_p), are shown in Fig. 1(d). It is evident that (11) is approximately satisfied, i.e., the monopole radiation dominates the scattered field for the sampling point that is indeed within the envelope detected by the LSM.

3. Proposed reconstruction approach: Modified linear sampling method

3.1. Reconstruction of the scatterer support based on multipole field

Now, we present a reconstruction approach based on the above argument. For ease of reference, especially in figures, we shall call the proposed method multipole-based linear sampling method (MLSM). Although one needs to consider infinite number of multipoles to fully account for the scattered field in (8), usually, a sufficiently large finite number of multipoles is enough to approximate E _∞(φ, θ), especially so in the presence of noise. Considering (2N + 1) number of multipoles, the expression for the far field (8) can be rewritten as:

E_{\infty} (φ, θ) \approx \sum_{n = - N}^{N} α^{(n)} ({\vec{r}}_{p}, θ) Φ^{(n)} (φ, {\vec{r}}_{p}) .

The Eq. (13) suggests that the multipole radiation functions Φ⁽ⁿ⁾(φ, r⃗_p) can be understood as a mapping from the effective multipole current at a sampling point r⃗_p to the measured scattered electric field. Consider a system of N_θ transmitters and N_φ receivers. For each incidence θ, (13) can be written as

\bar{E} = \overset{̿}{Φ} \bar{A},

where, E¯, an N_φ dimensional vector, consists of the various receiver measurements, φ̿, a matrix of dimension N_φ × (2N + 1), consists of the multipole radiation terms Φ⁽ⁿ⁾(φ, r⃗_p), and A¯, a (2N + 1) dimensional vector, contains the effective multipole currents α ⁽ⁿ⁾(r⃗_p, θ).

The value of A¯ can be solved uniquely using least squares pseudoinverse. After obtaining A¯ for each of the N_θ incidences, a discretized version of (11) is written:

\overset{̿}{A} \bar{g} = \bar{D},

where A̿, of dimension (2N + 1)×N_θ, contains the vectors A¯ for various incidences, g¯ is an N_θ dimensional vector that needs to be determined, and D¯ is a vector with all elements except the (N + 1) th element as zero. The (N +1) th element, which corresponds to a monopole, is 1. The value of g(θ,r⃗_p) can be determined from (15) by the least squares pseudoinverse.

3.2. Comparison of the conventional LSM and the proposed method with regard to g(θ,r⃗_p)

For the sake of convenience of reference, in this sub-section, we shall refer to the vector g¯ computed using LSM as g¯_LSM and the vector g¯ computed using MLSM as g¯_MLSM. In the framework of LSM, g¯_LSM is computed using the following equation:

{\bar{g}}_{LSM} = {\overset{̿}{E}}^{+} {\bar{Φ}}_{\infty},

where the matrix E̿ contains the vectors E¯ (defined after (14)) for all the incidences, Φ¯_∞ is the vector containing Φ_∞(φ, r⃗_p) for all the receivers, and the superscript ( + ) denotes pseudoinverse.

In LSM, typically the right hand side of (16) is calculated using Tikhonov regularization. In terms of the singular value decomposition (SVD) [29], the matrix E̿ is given as E̿∙v¯_s = σ_su¯_s, where u¯_s and v¯_s denote the s th left and right singular vectors, and σ_s denotes the s th singular value of the matrix E̿. Then the vector g¯_LSM can be computed using Tikhonov regularization as below:

{\bar{g}}_{LSM} = \sum_{\forall s} \frac{σ_{s}}{σ_{s}^{2} + α^{2}} {\bar{v}}_{s} {\bar{u}}_{s}^{*} {\bar{Φ}}_{\infty},

where α is the Tikhonov regularization parameter and the superscript * denotes the Hermitian operator. Generally, the regularization parameter α is chosen using the general discrepancy principle and is evaluated for each sampling point individually. However, in the current work, inspired from [24], we have used a constant value of the regularization parameter. It should be noted that Φ¯_∞ can be written as Φ¯_∞ = Φ̿D¯. Thus, effectively (14) and (15) can be combined to generate a form similar to (16).

However, for computing g¯_MLSM, first vectors A¯ for various incidences are computed using (14) as below:

\bar{A} = {\overset{̿}{Φ}}^{+} \bar{E} .

In the proposed method, we have used least squares based pseudoinverse. This is followed by the computation of g¯_MLSM using:

{\bar{g}}_{MLSM} = {({\overset{̿}{Φ}}^{+} ∙ \overset{̿}{E})}^{+} \bar{D} .

Once again, we have used least squares based pseudoinverse. Due to the intermediate step of computation of A¯, which depends upon the number of multipoles considered, g¯_MLSM cannot be written in a form similar to (16). Specifically, g¯_MLSM depends upon the number of multipoles considered through the presence of Φ̿⁺. However, when the value of N is very large, in the noise-free case, the solution of g¯_LSM is close to g¯_MLSM.

3.3. Choice of number of multipoles

The reconstruction method presented above is similar to LSM analytically if we consider a large N, such that the multipoles of order higher than the considered (2N +1) multipoles have very less and diminishing contribution to the scattered field E_φ(φ, θ). However, to solve (14) uniquely, at least (2N + 1) receivers are needed. Further, using large number of multipoles is expected to increase the computational complexity of the problem. This problem is further accentuated as the procedure involves two steps of inversion for each sampling point as compared to only one inversion per sampling point for LSM. Thus, it is interesting to study the effect of choosing small number of multipoles.

To this end, we propose to use N = 1 for imaging. The use of N = 1 implies the following. Now, solving (14) would mean that we seek an optimal combination of the monopole and dipole currents such that the resultant radiation fields match the scattered fields as closely as possible. Solving (15) would result in a g(θ,r⃗_p), such that the contribution from the determined dipole current is very small. Thus, effectively the requirement on g(θ,r⃗_p) has reduced, which would otherwise require to suppress the contribution from all higher multipoles. It is worth highlighting that although the monopole and dipole terms are dominant in most scattering problems, there indeed exist cases where other higher order multipoles may be dominant. For such cases, a suitable choice of N, such that the dominant multipole is included, can reconstruct the scatterer. This situation shall be considered in our future research work.

For the sake of illustration and further discussion, we present the reconstruction results for the geometry shown in Fig. 1(a). The measurement setup used is as before. We plot −log₁₀∥g(θ,r⃗_p)∥, where g(θ,r⃗_p) has been computed by the MLSM. Figure 2(a), 2(b) present the reconstruction results for N = 20 and N = 1 respectively. It is evident that N = 20 gives a reconstruction similar to LSM. For N = 1, it is noticeable that the background has a higher value of −log₁₀∥g(θ,r⃗_p)∥ as compared to LSM. This is expected because when N = 1, the strict requirement on g(θ,r⃗_p) to suppress all the higher order multipoles is now significantly eased. Now, g(θ,r⃗_p) requires suppressing only the dipole term to zero. This is also the reason that −log₁₀∥g(θ,r⃗_p)∥ is higher than LSM at the scatterers’ locations as well.

An estimate of the scatterer support $\sum^{´}$ can be determined as below:

\sum^{'} = {{\vec{r}}_{p} : - \log_{10} ∥ g (θ, {\vec{r}}_{p}) ∥ > Min + β (Max - Min)}

where, Min = min(−log₁₀∥g(θ,r⃗_p)∥ : ∀r⃗_p, Max = max(−log₁₀∥g(θ,r⃗_p)∥ : ∀r⃗_p), and β is the threshold used for estimating the scatterer support. The effect of threshold shall be discussed later. Presently we use β = 0.8. The scatterer’s supports estimated for N = 20 and the proposed method (N = 1) are presented in Fig. 2(c), 2(d). The above results presented are for the noise-free scenario. Results for the noisy situation (10% additive white Gaussian noise added to the measured electric field) are presented in Fig. 2(e)–2(h). It is evident that proposed method shows a better estimation of the scatterer’s shape than the LSM.

Fig. 2. Comparison of MLSM (N = 20) and MLSM (N = 1) for noise-free and noisy (10% additive Gaussian noise) scenarios.

Download Full Size | PDF

The truncation of all the higher order multipoles except the monopole and dipole terms is reasonable on account of the basic fact proven in section 2 that for a point belonging to the scatterer support, only the monopole radiation is the primary contributor to the radiation pattern. Based on this fact, two arguments are presented in support of our choice.

3.4. Sufficiency of N = 1: Effect of reduction of multipoles

For every sampling point, we compare the values of A¯ g¯ for N = 20 and N = 1 respectively, where A̿ is obtained for these two values of N, and g(θ,r⃗_p) is obtained by the conventional LSM. The rows of A̿ g¯ correspond to υ_n, n = −N to N, defined in (12). The absolute value of the difference between the effective multipole currents υ_n computed considering N = 20 and N = 1 for n = −1 to 1 is plotted in Fig. 3. It can be seen that the difference is very small (of the order 10^-4) for the three multipoles over the complete domain. Thus, the choice of N = 1 is reasonable.

Fig. 3. Effect of reduction of multipoles. After obtaining the values of υ_n by substituting the source distribution computed using conventional LSM and the matrix A̿ computed using Eq. (14) for N = 20 and N = 1 respectively, the absolute value of difference in υ_n corresponding to these cases for n = − 1 to 1 is computed. First, second, and third columns show this absolute value of difference for n equals to −1, 0, and 1 respectively.

Download Full Size | PDF

3.5. Sufficiency of N = 1: Effect of noise on various multipoles

Here, we study the contribution of noise to the error in the computation of the various multipoles. For this purpose, we consider two sampling points (0,0) and (−0.9λ,−0.9λ), belonging to the scatterer support and the background respectively and study the effect of noise on these sampling points in Fig. 4(a) and Fig. 4(b) respectively. In Fig. 4, the first column shows the values of υ_n, n = −20 to 20 obtained by the proposed method in the noise-free scenario; the second column is similar to the first column with the only difference being the noisy scenario (10% additive Gaussian noise); the third column displays the absolute value of the difference between the data plotted in first and second columns; the fourth column shows the difference relative to the magnitude of υ_n averaged over the noisy and noise-free cases. All the results in Fig. 4 have been averaged over 100 simulations. It can be seen that for the first sampling point, few multipoles are affected by noise. However, the effect of noise is very small (less than 0.05) for any multipole. Though the monopole is corrupted by noise, the relative effect of corruption on the monopole is very less as seen in the fourth column of Fig. 4(a). For the second point, the first and second columns suggest that though the monopole is the largest contributor, it is not the only prominent contributor (as expected). The third column suggests that more number of multipoles are corrupted by noise (as compared to the first sampling point). The level of corruption is also higher. The last column suggests that monopole is least affected by noise in this case as well. Thus, it can be concluded that the monopole shall not become the only prominent contributor even in the presence of noise. Further, by choosing N = 1, in either case, the nature of monopole will remain the same, irrespective of the presence or absence of noise. Thus, considering N = 1 is reasonable.

Fig. 4. Effect of noise on the multipoles. First column shows the values of υ_n obtained by the proposed method in the noise-free scenario. Second column is similar to the first column with the only difference being the noisy scenario (10% additive Gaussian noise). Third column displays the absolute value of the difference between the data plotted in first and second columns. Fourth column shows the difference relative to the magnitude of υ_n averaged over the noisy and noise-free cases.

Download Full Size | PDF

3.6. Nature of approximation: Comparison with Tikhonov regularization

In the conventional LSM, the use of Tikhonov regularization parameter effectively results into the approximation of the monopole radiation. However, the nature of approximation in the proposed method is entirely different. Tikhonov regularization serves to approximate the monopole radiation in terms of the strengths given to the various spectral vectors that span the space of the scattered fields. Thus, it is a mathematical approximation which can be tuned or changed by changing the value of the Tikhonov regularization parameter α. On the other hand, the proposed method employs the approximation in terms of the truncation of higher order multipoles. This approximation is justified from the perspective of physics because the contribution of higher order multipoles to the scattered field is indeed small. Further, this approximation is a definitive approximation, i.e., only monopole and dipole radiation are kept, which is different from the Tikhonov regularization that is sensitive to regularization parameter α.

To illustrate the fact, we apply conventional LSM for reconstruction of the example in Fig. 1(a) using various values of the regularization parameter α. We have chosen α = aσ ₁, where a takes values {0.001, 0.01, 0.1,1} and σ ₁ is the largest singular value of the matrix E̿. The reconstruction results are presented in Fig. 5(a). It is evident that the performance of the conventional LSM is sensitive to α. It is also evident that the conventional LSM cannot be used to obtain reconstruction results similar to the proposed method by choosing a suitable α. In other words, approximation in the sense of Tikhonov regularization is different from the approximation N = 1. It is worth mentioning here that conventionally the value of α is not kept constant. It is rather computed using general discrepancy principle for each pixel in the domain. However, our simulations indicate that in the presence of 10% noise, the value of α on the scatterer support and region surrounding it is close to 0.1. Since our intention is to compare the nature of Tikhonov regularization and the proposed method, choosing a constant value of α does not affect the general conclusions.

At this point, we define an error parameter that can be used to evaluate the quality of reconstruction. The error measure is defined on the convex hull of the scatterer. The threshold β classifies the sampling points as scatterers or non-scatterers. If the number of sampling points inside the convex hull that get classified wrongly is M_err(β), then the error parameter can be defined as below:

Error (β) = \frac{M_{err} (β)}{M},

where M is the total number of sampling points inside the convex hull. This error measure is useful to evaluate the accuracy of the reconstruction method in reconstructing the scatterer support with complicated boundary or boundary that is not simply connected.

For the examples presented in this paper, we vary threshold β in the range [0.6,0.9] and plot the error measure Error (β) for conventional LSM (for various values of α) and the proposed method in Fig. 5(b). It should be noted that in practice, a suitable value of β is not known and β is chosen heuristically in most cases. It is evident that the proposed method has lower value of Error (β) for most values of β.

Fig. 5. Comparison of LSM and MLSM. In (b), ‘Error’ refers to the error measure defined in (21).The above results are in the presence of 10% noise.

Download Full Size | PDF

4. More numerical examples

In order to validate the effectiveness of proposed method, we consider various examples and compare the performance of LSM and the proposed method. Each of the examples considered here is difficult to reconstruct qualitatively using conventional LSM on various accounts. These include the presence of multiple scatterers in close proximity and scatterers with boundaries that are not simply connected. Further details for each example are discussed in the corresponding paragraphs. All the examples consider a square region of size 2×2 m², where the scatterers are placed in free space. The frequency of the incident wave is 300 MHz (wavelength λ = 1 m). The measurement setup is same as described in section 2.3. The measured data is corrupted by 10% additive Gaussian noise. Three examples of dielectric scatterers (Fig. 7) and two examples of perfectly conducting scatterers (Fig. 8) are considered. In Fig. 7 and Fig. 8 the first column shows the scatterer profile (relative permittivity for dielectric cylinders, the contours for perfectly conducting cylinders), the second column shows the reconstruction using conventional LSM and the third column shows the reconstruction using the proposed method. The regularization parameter used in LSM has been computed using the general discrepancy principle.

The error measures for the three examples of dielectric scatterers are plotted in Fig. 6. It should be noted that this error measure is not defined for the perfectly conducting cylinders. This is because in the case of perfectly conducting cylinders, the scatterer profile and quality of reconstruction are described in terms of the contour of the cylinder.

Fig. 6. Plot of error measure defined in (21) for the examples of dielectric scatterers presented in Fig. 7. The presented results are in the presence of 10% noise.

Download Full Size | PDF

4.1. Dielectric scatterers

Example 1: We consider a profile that consists of two circular cylinders and one annular cylinder, each of relative permittivity 2. The circular cylinders are of radius 0.2 m each and are centered at (0.3, 0.6) m and (−0.3, 0.6) m. The ring has an exterior radius of 0.6 m and an inner radius of 0.3 m, and is centered at (0, −0.2) m. This profile is commonly known as the Austria profile.

In the Austria profile [see Fig. 7(a)], it is difficult to detect the inner hollow of the annulus using LSM. The error measure, plotted in Fig. 6(a) and the reconstruction results demonstrate that MLSM performs better than LSM. Using MLSM, not only the inner hollow of the annular structure is visibly detected, the reconstructed profile is also closer to the actual profile.

Example 2: We consider an example where three circular cylinders are placed linearly and are obstructed by two bars on either side of the arrangement. The circular cylinders are of radius 0.2 m each, and are centered at (0, −0.6) m, (0, 0) m, and (0, −0.6) m. The bars are of size 0.2 m × 1.5 m, placed parallel to the y axis, and centered at (−0.6, 0) m and (0.6, 0) m respectively. The relative permittivity of all the scatterers is 2.

This example is particularly difficult for two reasons. First, the circular structures are obstructed by the bars on two sides (left and right). On the sides where the circular cylinders are not obstructed by the bars (top and bottom), they are placed one behind another. Second, all the cylinders are in close proximity with one another. Due to this they are difficult to resolve. Similar to the first example, MLSM demonstrates better performance [see Fig. 6(b) and Fig. 7(b)] and all the scatterers can be resolved using MLSM.

Example 3: We consider an example where two circular scatterers are enclosed within an annular square scatterer (a setting similar to through wall imaging). The annular square scatterer is of external dimension 1.5 m and internal dimension 1.4 m. Its relative permittivity is 2 and it is centered at (0, 0.1) m. The circular cylinders are of radii 0.15 m and 0.2 m respectively and centered at (−0.3, −0.2) m and (0.2, 0.2) m respectively. The circular cylinders have relative permittivity 2.5. See Fig. 7(c).

In this example, though it is expected that the square enclosure is reconstructed well, the reconstruction of the inner circular cylinders is difficult due to the presence of the enclosure. The difference in the sizes of the circular cylinders may also make the reconstruction more challenging. For this example, Fig. 6(c) shows a comparable performance of MLSM and LSM. However, the reconstructed results shown in Fig. 7(c) clearly show that the quality of reconstruction using MLSM is significantly better. MLSM is able to better detect the square enclosure as well as the two circular scatterers as compared to LSM.

Fig. 7. Examples of reconstruction of dielectric cylinders. The first column shows the scatterer profile (relative permittivity), the second column shows the reconstruction using conventional LSM and the third column shows the reconstruction using the proposed method. The presented results are in the presence of 10% noise.

Download Full Size | PDF

4.2. Perfectly conducting scatterers

Example 1: The first example in this category consists of three cylinders placed close together [Fig. 8(a)]. The radii of the cylinders are 0.6 m, 0.3 m, and 0.3 m respectively. They are centered at (0, −0.3) m, (−0.4 0.6) m, and (0.4, 0.6) m respectively.

This example is difficult to reconstruct due to two reasons. First the cylinders are in close proximity with each other. Second, the larger cylinder (being perfectly conducting) is obstructing the view of the smaller cylinders and limiting the aspect available for these cylinders. Due to this, it is generally expected that only the outer contour of the group of cylinders is reconstructed easily. The reconstruction results show that LSM is able to detect only parts of their contours, while MLSM is able to show the presence of three cylinders very well.

Example 2: The second example consists of a two dimensional array of nine circular cylinders [Fig. 8(b)]. Each cylinder is of radius 0.2 m and their centers are on (±0.6, ± 0.6) m, (±0.6,0) m, (0, ±0.6) m, and (0, 0) m.

This example is difficult to reconstruct due to the presence of many cylinders in close proximity. Being perfectly conducting, the outer cylinders limit the aperture available to other cylinders. Due to this, the cylinder in the centre of the arrangement is least likely to be detected. For the other cylinders as well, the outer contour of the group of cylinders is the most likely detection. The reconstruction results show that the conventional LSM is able to detect only the external boundary of the collection of the cylinders, while MLSM is able to detect the internal contours as well.

Fig. 8. Examples of reconstruction of perfectly conducting cylinders. The first column shows the scatterer profile (contours of cylinders), the second column shows the reconstruction using conventional LSM and the third column shows the reconstruction using the proposed method. The presented results are in the presence of 10% noise.

Download Full Size | PDF

5. Conclusion

A modified linear sampling method has been proposed for the reconstruction of the scatterers’ support. The proposed method uses a physical approximation of the LSM formulation rather than the conventionally used mathematical regularization. Instead of studying the circular symmetry of the induced current distribution, the multipole expansion of the scattered field is studied here.

The proposed reconstruction approach truncates all higher order multipoles except the dipole and monopole terms. The error obtained due to the truncation of higher order multipoles and the effect of noise on monopole show that the truncation of higher order multipoles is justified. It is also shown that the physical approximation (obtained by truncating higher order multipoles) is different from Tikhonov regularization and is expected to perform better. The proposed method demonstrates good performance for various complicated scatterer supports.

Current work is a first step towards the physical explanation of LSM. For this reason, the current work considers the simplest scalar electromagnetic scattering problem (two-dimensional scatterers under TM illumination). The next step in this direction would be to develop MLSM for non-scalar inverse scattering problems like three-dimensional imaging and two-dimensional imaging under TE illumination. This extension can be done with the help of some theorems in linear algebra. It is worth mentioning that care should be taken during this extension since the fundamental radiating source is dipole in these two cases, instead of monopole. We are hopeful that a physics based understanding of LSM would eventually lead to wider applicability of LSM for more complicated problems, for example involving special materials, such as anisotropic scatterers. The enhanced imaging ability of the proposed method suggests it can be useful to the through-wall-imaging, as shown in Fig. 7(c), which obviously finds wide real-life applications.

Acknowledgements

The authors wish to acknowledge the financial support from the Singapore Ministry of Education under the grant R263000485112.

References and links

1. K. Agarwal and X. Chen, “Application of differential evolution in 2-dimensional electromagnetic inverse problems,” in IEEE Congress on Evolutionary Computation, (2007), 4305–4312. [CrossRef]

2. K. Agarwal and X. Chen, “Applicability of MUSIC-Type Imaging in Two-Dimensional Electromagnetic Inverse Problems,” IEEE Trans. Antenn. Propag. 56(10), 3217–3223 (2008). [CrossRef]

3. Y. Zhong and X. Chen, “MUSIC imaging and electromagnetic inverse scattering of multiple-scattering small anisotropic spheres,” IEEE Trans. Antenn. Propag. 55(12), 3542–3549 (2007). [CrossRef]

4. M. Brignone, G. Bozza, A. Randazzo, M. Piana, and M. Pastorino, “A Hybrid Approach to 3D Microwave Imaging by Using Linear Sampling and ACO,” IEEE Trans. Antenn. Propag. 56(10), 3224–3232 (2008). [CrossRef]

5. I. Catapano, L. Crocco, M. D’Urso, and T. Isernia, “On the effect of support estimation and of a new model in 2D inverse scattering problems,” IEEE Trans. Antenn. Propag. 55(6), 1895–1899 (2007). [CrossRef]

6. X. Chen and Y. Zhong, “A robust noniterative method for obtaining scattering strengths of multiply scattering point targets,” J. Acoust. Soc. Am. 122(3), 1325–1327 (2007). [CrossRef] [PubMed]

7. H. Ammari, E. Iakovleva, D. Lesselier, and G. Perrusson, “Music-type electromagnetic imaging of a collection of small three-dimensional inclusions,” SIAM J. Sci. Comput. 29(2), 674–709 (2007). [CrossRef]

8. X. Chen and Y. Zhong, “MUSIC electromagnetic imaging with enhanced resolution for small inclusions,” Inverse Probl. 25(1), 015008 (2009). [CrossRef]

9. M. Fink and C. Prada, “Acoustic time-reversal mirrors,” Inverse Probl. 17(1), 201 (2001). [CrossRef]

10. A. J. Devaney, “Time reversal imaging of obscured targets from multistatic data,” IEEE Trans. Antenn. Propag. 53(5), 1600–1610 (2005). [CrossRef]

11. T. Rao and X. Chen, “Analysis of the time-reversal operator for a single cylinder under two-dimensional settings,” J. Electromagn. Waves Appl. 20(15), 2153–2165 (2006). [CrossRef]

12. K. Mayer, R. Marklein, K. J. Langenberg, and T. Kreutter, “Three-dimensional imaging system based on Fourier transform synthetic aperture focusing technique,” Ultrasonics 28(4), 241–255 (1990). [CrossRef]

13. C. K. Liao, M. L. Li, and P. C. Li, “Optoacoustic imaging with synthetic aperture focusing and coherence weighting,” Opt. Lett. 29(21), 2506–2508 (2004). [CrossRef] [PubMed]

14. T. Arens and A. Kirsch, “The factorization method in inverse scattering from periodic structures,” Inverse Probl. 19(5), 1195–1211 (2003). [CrossRef]

15. A. Kirsch, “The factorization method for Maxwell’s equations,” Inverse Probl. 20(6), S117–S134 (2004). [CrossRef]

16. A. Kirsch, “Characterization of the shape of a scattering obstacle using the spectral data of the far field operator,” Inverse Probl. 14(6), 1489–1512 (1998). [CrossRef]

17. A. Kirsch and S. Ritter, “A linear sampling method for inverse scattering from an open arc,” Inverse Probl. 16(1), 89–105 (2000). [CrossRef]

18. D. Colton, J. Coyle, and P. Monk, “Recent developments in inverse acoustic scattering theory,” SIAM Rev. 42(3), 369–414 (2000). [CrossRef]

19. F. Cakoni, D. Colton, and H. Haddar, “The linear sampling method for anisotropic media,” J. Comput. Appl. Math. 146(2), 285–299 (2002). [CrossRef]

20. D. Colton, H. Haddar, and M. Piana, “The linear sampling method in inverse electromagnetic scattering theory,” Inverse Probl. 19(6), S105–S137 (2003). [CrossRef]

21. D. Colton and R. Kress, “Using fundamental solutions in inverse scattering,” Inverse Probl. 22(3), R49–R66 (2006). [CrossRef]

22. F. Cakoni and D. Colton, Qualitative methods in inverse scattering theory: An introduction (Springer, Berlin, 2006).

23. A. Kirsch and N. Grinberg, The factorization method for inverse problems (Lecture series in mathematics and its applications, Oxford University Press, 2008).

24. I. Catapano, L. Crocco, and T. Isernia, “On simple methods for shape reconstruction of unknown scatterers,” IEEE Trans. Antenn. Propag. 55(5), 1431–1436 (2007). [CrossRef]

25. N. Shelton and K. F. Warnick, “Behavior of the regularized sampling inverse scattering method at internal resonance frequencies,” J. Electromagn. Waves Appl. 17(3), 487–488 (2003). [CrossRef]

26. I. Catapano, L. Crocco, and T. Isernia, “Improved Sampling Methods for Shape Reconstruction of 3-D Buried Targets,” IEEE Trans. Geosci. Rem. Sens. 46(10), 3265–3273 (2008). [CrossRef]

27. D. L. Colton and R. Kress, Inverse acoustic and electromagnetic scattering theory, 2nd ed. (Springer, New York, 1998).

28. M. Abramowitz and I. A. Stegun, Handbook of mathematical functions: with formulas, graphs, and mathematical tables, 7th ed. (Dover Publications, New York, 1972).

29. R. A. Horn and C. R. Johnson, Matrix analysis (Cambridge University Press, Cambridge; New York, 1985).

A multipole-expansion based linear sampling method for solving inverse scattering problems

Abstract

1. Introduction