Multiresolution Virtual Experiments for Microwave Imaging of Complex Scenarios

Abstract

In this paper, a multiresolution approach for the quantitative microwave imaging of complex scenarios is introduced. The proposed strategy takes advantage of the combined use of a recently introduced iterative method known as distorted iterated virtual experiments (DIVE), based on the paradigm of “virtual experiments”, and a wavelet-based projection scheme. This strategy allows the unknown profiles to be represented at different resolution scales and, as such, it is particularly suitable for the imaging of highly heterogeneous targets. Moreover, the developed algorithm blends together the intrinsic multiresolution feature of the wavelet projection with the one gained by means of a frequency hopping technique. The method was tested against realistic heterogeneous scenarios of practical interest, such as breast and tree trunk phantoms, which are of interest in non-invasive medical diagnostics and the health monitoring of standing trees.

1. Introduction

The reconstruction of the electrical properties (EPs) of unknown objects by means of microwave imaging (MWI) can represent a promising tool for non-invasive diagnostics. In fact, microwaves can interact with matter and provide useful information about it in a non-destructive way, without using ionizing radiation. However, MWI is affected by a limited spatial resolution, and involves the challenging task of solving an electromagnetic inverse scattering problem, which is non-linear and ill-posed [1]. Notably, these difficulties are worsened when the unknown profile is highly heterogeneous. Many strategies have been proposed in the literature to tackle both non-linearity and ill-posedness of inverse scattering problems, especially in the case of complex scenarios, such as the cases dealing with biomedical imaging or underground inspection [2,3,4,5,6,7,8,9,10].

As far as non-linearity is concerned, we can distinguish among approximation-based methods, which, while having a limited range of applicability, actually solve a linear problem [11,12]; and non-linear methods, which instead aim at solving the full non-linear problem. Among the latter, it is worth recalling global optimization approaches [13] and local optimization procedures [14,15,16,17]. The former have the advantage of being free from false solution issues, but the huge dimension of real case problems prevents their actual applicability in practice. For this reason, local optimization procedures are mostly adopted in inverse scattering problems.

As for the ill-posedness of the inverse scattering problem, regularization strategies are necessary in order to get reliable results [18,19]. Finally, it is worth mentioning the importance of exploiting multi-frequency data in inverse scattering, which also plays a beneficial role in the ill-posedness of the problem, allowing the accuracy and stability of the final results to be improved [20].

With respect to this frame, this paper proposes an innovative approach for quantitative MWI that is particularly suitable for the reconstruction of complex scenarios. It is based on a local optimization approach, and combines a recently introduced iterative procedure based on the virtual experiments (VE) framework, called distorted iterated virtual experiments (DIVE) [21], and a multiresolution strategy based on regularization by projection, wherein the wavelet basis is conveniently exploited to represent the unknown EPs of the target [22].

In [21], the DIVE method has been presented and tested in canonical scenarios, whose EPs were retrieved by exploiting the truncated singular value decomposition (TSVD) [23] or the total variation–compressive sensing (TV-CS) approach [24]. Both approaches represent a regularized inversion scheme for the linear problem involved at each DIVE iteration. In the case of TSVD, the regularization is performed by truncating the SVD of the scattering operator to a certain index, which leads to a smooth reconstruction of the actual profile. The TV-CS approach is instead based on the sparsity concept and the minimization of the ${\ell }_1$-norm of the gradient of profiles, which are stepwise constant. These different regularization approaches are adopted in order to maximize the performance of the imaging algorithm, according to the class of profiles one is aiming at retrieving (on the basis of the available a priori information on the scenario at hand). In the case of highly heterogeneous and complex scenarios, both TSVD and TV-CS approaches have turned out to be neither effective nor reliable. In this respect, the aim of this paper is to introduce an enhanced version of DIVE that fits well to complex scenarios by acting on the regularization scheme.

The proposed strategy, named Multiresolution-DIVE (M-DIVE), solves the problem by spanning different physical scales within a stepwise refinement framework, which includes both frequency hopping [25] and multiscale/multiresolution techniques [26]. In more detail, M-DIVE takes advantage of the intrinsic capability of wavelet transform to represent heterogeneous targets at different resolution scales, thus reducing the ill-posedness of the linear problem to be solved at each DIVE iteration. Moreover, for the first time, DIVE is used in conjunction with multi-frequency inversion techniques in order to further improve the accuracy and reliability of the final reconstruction. Accordingly, M-DIVE starts by retrieving a low-resolution reconstruction by processing only the lower-frequency data, and this reconstruction is used as a starting estimation for higher-frequency reconstruction.

Starting from these considerations, and from the preliminary results obtained in [27], this paper further investigates the potential of this approach with respect to complex scenarios. In particular, two test cases of practical interest are considered. The first one concerns medical MWI, which is acquiring a great deal of interest in medical diagnosis because it can offer significant advantages over other medical imaging techniques—namely, the use of non-ionizing radiation and the possibility of cheap and portable devices. In particular, it has gained increasing interest in the field of breast diagnostics [28,29]. Breast profiles are characterized by a very high degree of inhomogeneity, both in terms of geometry and EPs, hence, they allow an opportunity to fully exploit the multiresolution capability of the proposed approach. The adoption of a wavelet representation-based approach for biomedical imaging is not new in the literature [30,31]. It was previously proposed in conjunction with the distorted Born iterative method (DBIM) scheme [31]. However, since the DIVE procedure has shown an extended validity and better performance with respect to DBIM [21], the proposed M-DIVE is expected to provide better results.

The second test case concerns the inspection of wood materials, around which interest has been growing in recent years. The reasons for this are related to the possibility of achieving useful information about the healthy state of standing trees in order to prevent falls and optimize cutting strategies [32,33,34], but also to assess the quality of the wood product in industrial processes [32,33,34]. Moreover, the use of microwaves for these kinds of inspections seems to be attractive with respect to ultrasound or X-ray-based techniques, especially in terms of cost and data acquisition time (see [32,33,34,35] and references therein).

The paper is organized as follows. In Section 2.1, the mathematical formulation of the problem underlying the MWI and DIVE scheme is given. In Section 2.2, the proposed multiresolution approach is introduced and described, while in Section 3, an extensive numerical analysis against realistic breast and tree trunk phantoms is reported. Finally, the discussion of the results and conclusions follow.

2. Materials and Methods

2.1. Mathematical Formulation and DIVE Scheme

Let us consider an unknown non-magnetic object embedded in a homogeneous medium of known electromagnetic features, whose cross section (with compact support) $\Sigma$ is hosted in the imaging domain $\Omega$. With reference to the canonical 2D scalar case, the region being tested is probed by a set of T incident fields (TM polarized) transmitted by an array of antennas located on a closed curve $\Gamma$ around $\Omega$. The resulting scattered fields are measured by N receiver antennas, also located on $\Gamma$. By assuming and dropping the time harmonic factor $exp\left\{j\omega t\right\}$, the equation describing the scattering phenomenon, at the frequency $f$, is [1]:

$${E_s}^{(v)}(\mathit{R},f)={\mathcal{A}}_e[\chi ({\mathit{r}}^{\prime },f)E^{(v)}({\mathit{r}}^{\prime },f)],$$

(1)

where the subscript $v=1,\dots ,T$ identifies the generic incident field; $\mathit{R}\in \Gamma$; ${\mathit{r}}^{\prime }\in \Omega$; $E$ and $E_s$ are the total field in $\Omega$ and the scattered field on $\Gamma$, respectively; and ${\mathcal{A}}_e$ is a short notation for the integral radiation operator. Finally, the EPs of the unknown object are encoded in the contrast function $\chi ({\mathit{r}}^{\prime })=\epsilon ({\mathit{r}}^{\prime })/{\epsilon }_b-1$, wherein $\epsilon$ and ${\epsilon }_b$ are the complex permittivity of target and background, respectively.

The solution of Equation (1) for the estimation of $\chi$ from the measured scattered fields $E_s$ implies the overcoming of both non-linearity, as the total field also depends on the unknown contrast, and the ill-posedness, due to the properties of ${\mathcal{A}}_e$ [1].

A possible effective approach to tackle the non-linearity has been introduced in recent years within the VE framework [21]. The VE framework was derived from the circumstance that the linear superposition of different incident fields gives rise to the same linear superposition of the corresponding scattered fields. In particular, the use of these new scattering experiments has allowed a new field approximation to be introduced, to linearize the problem in (1). Finally, in order to enlarge the range of validity of such VE-based approximation, the iterative DIVE method [21] has recently been proposed and tested in cases of canonical targets. This scheme is intrinsically different from DBIM, as it takes into account the nature of the scatterers from the first step. Notably, its performance has been shown to be remarkably better [21].

The DIVE procedure can be summarized in five steps, as shown in Figure 1. In the initialization step, a first estimate ${\chi }^0$ is achieved by using the VE-based approximation, or other convenient starting guesses, depending on the available a priori information. In the update step, the forward scattering problem pertaining to the current estimate ${\chi }^{k-1}$ is solved to update the Green’s function and compute the new data—that is, the anomalous field $\Delta {E_{\mathrm{s}}}^k$. Then, a convergence control is performed by defining the relative residual error at the kth iteration: if the stopping rule is not satisfied, possible corrections $\Delta {\chi }^k$ with respect to ${\chi }^{k-1}$ are estimated by considering a distorted formulation of the scattering problem. Thereafter, the distorted linear sampling method [36] is exploited to localize $\Delta {\chi }^k$ and design a new VE set. The last step deals with the solution of the relevant distorted linear inverse problem, which is recast as:

$$\Delta {{ℰ}_s}^k(\mathit{R},f)={{\mathcal{A}}_e}^k\left[\Delta {\chi }^k({\mathit{r}}^{\prime },f){ℰ}^k({\mathit{r}}^{\prime },f)\right]=ℒ(\Delta {\chi }^k({\mathit{r}}^{\prime },f)),$$

(2)

where $\Delta {{ℰ}_s}^k$ and ${ℰ}^k$ are the anomalous and approximated total fields arising in the VE, respectively, and ${{\mathcal{A}}_e}^k$ is the external radiation operator pertaining to the reference scenario at the kth iteration, while $ℒ$ indicates the linear operator which applies to $\Delta {\chi }^k$. Note that, unlike DBIM, the total field is approximated by taking into account the contribution of the anomaly, thanks to the VE-based approximation. Once Equation (2) is solved by adopting a convenient regularization technique, the contrast update is pursued by adding the reconstructed perturbations to the current reference scenario, namely ${\chi }^{k+1}={\chi }^k+\Delta {\chi }^k$. Note that, depending on the a priori information about the target, different regularization techniques can be adopted to solve Equation (2). Finally, the procedure continues until the stopping criterion is fulfilled.

More details about the VE framework and basics of DIVE can be found in Reference [21] and references therein.

2.2. Multiresolution DIVE

In this paper, the DIVE approach is tailored with respect to complex and heterogeneous scenarios. To this end, it is integrated into a stepwise refinement framework in order to combine both the intrinsic multiresolution feature of the regularization by projection on wavelet basis [26] and the possibility offered by a frequency-hopping-based strategy of increasing the resolution scale by simply changing the frequency of interest. The processing of multi-frequency data allows an increase in the data-unknown ratio. In particular, the frequency hopping technique plays a very important role, as low frequencies make it possible to locate the objects and to roughly reconstruct them, while higher frequencies allow finer details to be retrieved [20].

On the other hand, the wavelet basis offers the possibility to decompose the unknown profile into two sets of coefficients, namely, coarse and detail coefficients, and to accommodate the trade-off between the efficiency and accuracy of the representation [22]. The thus-obtained M-DIVE scheme is able to zoom in on details progressively, by starting with a coarse estimate of the unknown and then further improving it.

Consequently, the unknown $\Delta {\chi }^k$ in Equation (2) is represented as the superposition of wavelet bases ${\mathcal{W}}_n(\mathit{r})$:

$$\Delta {\chi }^k({\mathit{r}}^{\prime },f)={\displaystyle \sum _n}\Delta {{\mathcal{W}}_n}^k(f){ℳ}_n(f){\mathcal{W}}_n({\mathit{r}}^{\prime }),$$

(3)

where ${\mathcal{W}}_n$ are the wavelet bases, $\Delta {{\mathcal{W}}_n}^k$ are the wavelet coefficients (which are the actual unknowns of the problem), and ${ℳ}_n$ is a binary mask which selects only the coarse coefficients at each frequency and decomposition level. Then, the solution of the linear problem (2) is achieved by minimizing the following cost functional:

$$\Phi (\Delta {{\mathcal{W}}_n}^k)=\Vert \Delta {{ℰ}_s}^k(\mathit{R},f)-ℒ{(\Delta {\chi }^k({\mathit{r}}^{\prime },f))\Vert }_2^2$$

(4)

through an “inner” iterative procedure, ${\Vert ·\Vert }_2$ being the ${\ell }_2$-norm. In particular, the inversion of Equation (4) is performed by integrating the representation (3) into a conjugate gradient least square technique. More details can be found in [31].

Note that, during the overall M-DIVE procedure, the decomposition level of the wavelet representation is changed throughout the different frequencies. Accordingly, M-DIVE starts from a high level of decomposition (or, equivalently, a high-order coarse representation), and then the level is progressively reduced to retrieve finer details in the image, thus improving image resolution.

3. Results

In the following text, two relevant examples are presented in order to give a proof of the performance achievable by means of M-DIVE, both dealing with heterogeneous profiles coming from two real applications involving complex scenarios—that is, breast phantom imaging [29,30,31] and tree trunk inspection [32,33,34].

In both examples, the scattered field data are simulated using an in-house developed forward solver, based on the method of moments, while at each iteration of the M-DIVE scheme, the pertaining model to build the VE is attained by means of full-wave simulations.

3.1. Breast Phantom Imaging

Two-dimensional examples concerned with quantitative reconstructions of the EPs of a slice of two anthropomorphic breast phantoms are presented. To build the 2-D profile, the 92nd slice from phantom ID 062204 (representing a heterogeneously dense breast) and the 65th slice from phantom ID 012304 (representing a very dense breast), both belonging to the University of Wisconsin repository [37] (available online), were extracted and resized in such a way that the actual cell dimension was 0.5 mm. The considered background medium was lossless and had relative permittivity equal to 18 [38]. Moreover, the adopted multi-view multi-static measurement configuration considers T = N = 20 and T = N = 24 sources, modelled as line sources, surrounding the breast on a circumference of radius 0.0824 m and 0.1 m for phantom IDs 062204 and 012304, respectively. This number of antennas permits collection of all the available information at the lowest frequency of 1 GHz in a non-redundant way [39]. The processed data were then corrupted by an additive Gaussian noise, with signal-to-noise ratio (SNR) equal to 30 dB.

The range of the considered frequencies was 1–4 GHz, with a frequency step of 250 MHz [31]. The mother wavelet adopted during the procedure was the Daubechies20. In particular, at the first three frequencies, a 3rd-level decomposition was performed, while in the frequency range 1.75–2.25 GHz, a 2nd decomposition level was considered. Finally, the 1st level of decomposition was adopted in the range 2.5–3 GHz, and 0th-level in the range 3.25–4 GHz.

The obtained reconstruction, and some intermediate results, are shown in Figure 2 and Figure 3. Note that at 1 GHz, the initial guess ${\chi }^0$ was obtained by assuming the actual shape of the breast to be known, and by filling it with homogeneous fatty tissue. Moreover, in the solution of problem (2), constraints were enforced on the admissible EPs values of the investigated region being tested in order to get permittivity values higher than 1 and conductivity values higher than 0.

As it can be seen in Figure 2 and Figure 3, the morphology and the EPs of the imaged breast profiles were accurately retrieved. The accuracy of the obtained reconstructions was quantified by means of the normalized mean square reconstruction error (NMSE), computed separately for the complex permittivity function $\epsilon$, and for the real (${\epsilon }^{\prime }$) and imaginary (${\epsilon }^{″})$ parts of the latter. Note that the NMSEs reported in

Table 1. NMSE of heterogeneously dense breast (ID 062204).

Frequency	NMSE on ε	NMSE on ε′	NMSE on ε″	NMSE on ε from [31]
1 GHz	0.29	0.28	0.69	0.41
2 GHz	0.19	0.18	0.47	Not provided
3 GHz	0.13	0.12	0.45	0.28
4 GHz	0.11	0.10	0.42	-

and

Table 2. NMSE of very dense breast (ID 012304).

Frequency	NMSE on ε	NMSE on ε′	NMSE on ε″	NMSE on ε from [31]
1 GHz	0.22	0.20	0.66	0.39
2 GHz	0.15	0.14	0.47	Not provided
3 GHz	0.10	0.09	0.41	0.29
4 GHz	0.08	0.07	0.40	-

were lower with respect to the ones in [31], reported in the last columns of the tables, thus confirming the better capability of the M-DIVE approach in imaging highly heterogeneous profiles. Note also that in [31], in contrast to the M-DIVE examples, the highest considered frequency was 3 GHz, possibly because no significant improvements were observed as frequency increased.

3.2. Tree Trunk Inspection

The tree being tested was a four-layered oak phantom with dielectric properties fixed by following References [32,40]. In particular, the inner heartwood layer had $\epsilon =9$ and $\sigma =0.005$ S/m; the sapwood layer was characterized by $\epsilon =16$, $\sigma =0.01$ S/m; next, the thicker layer was the phloem, with $\epsilon =22$, $\sigma =0.05$ S/m; finally, the bark had $\epsilon =6.2$, $\sigma =0.02$ S/m. The outer diameter of the trunk was about 50 cm and a lossless matching medium with permittivity equal to 4.7 was chosen according to Reference [40]. In order to mimic a more realistic environment, we added a random variation of ±5% around the average EPs values of the trunk. The investigation domain was discretized into 64 × 64 square cells, while the multi-view multi-static data, corrupted by 30 dB of noise, were collected by T = N = 21 antennas located at a distance of about 12.5 cm from the trunk.

For the sake of monitoring the healthy state of a tree, we simulated some cracks inside the trunk, identified by the dielectric parameters of the void ($\epsilon =1$, $\sigma =0$). When an exact knowledge of the tree structure is available, such anomalies can be detected by exploiting distorted or differential imaging techniques [29,41]. However, in most cases, this kind of a priori information is not available and a full imaging of the trunk under investigation has to be pursued.

In this example, the considered frequency range was 0.1–1 GHz, with a step of 150 MHz. Note that with respect to the previous medical examples, different frequency ranges were adopted due to the different size of the target being tested (50 cm versus 8 cm). Moreover, the wavelet decomposition level was equal to 3 for frequencies 100 and 250 MHz, 2 for 400 and 550 MHz, and 1 for 700 and 850 MHz. Finally, a 0th decomposition level was adopted at the last frequency.

The outcome of the inversion is reported in Figure 4. Note that the starting point of the iterative inversion algorithm was set as a homogeneous trunk having the actual dimension of the phantom and EPs of the phloem layer. As it can be seen, the M-DIVE scheme was able to identify the different layers and the presence of cracks inside the trunk. The accuracy of the reconstruction was also confirmed by the achieved low values of NMSE (

Table 3. NMSE of oak tree trunk.

Frequency	NMSE on ε	NMSE on ε′	NMSE on ε″
100 MHz	0.15	0.14	0.19
400 MHz	0.05	0.05	0.28
700 MHz	0.03	0.03	0.63
1 GHz	0.02	0.02	0.97

).

4. Conclusions

In this contribution, a novel technique for quantitative MWI of complex scenarios was introduced and tested. The proposed strategy exploits the DIVE scheme, which is an iterative method recently introduced and assessed in relation to canonical targets, based on the paradigm of “virtual experiments”. In order to introduce an enhanced version that fits well in complex scenarios, DIVE was herein exploited for the first time in conjunction with a frequency hopping technique and a wavelet-based multi-resolution representation of the unknown. By doing so, a multiresolution version of the DIVE (M-DIVE) scheme was introduced.

The developed imaging technique was tested with respect to two different applications involving heterogeneous profiles, namely, realistic anthropomorphic breast and tree trunk phantoms. The reconstruction images and the synthetic errors confirmed the capability of M-DIVE to deal with complex and heterogeneous scenarios.

Note that while the two applications are different and, as such, consider different measurement setups and settings, the general methodology (i.e., the way in which the data coming from a generic measurement setup are processed) does not change. Of course, when moving to real applications, different settings (e.g., type and polarization of the antennas, frequency range, etc.) must be properly chosen and the data-processing procedure tuned accordingly, for example, by means of a suitable calibration of the measured data.