An approach to reconstruct buried objects is proposed. It is based on the integral equations of the electromagnetic inverse scattering problem, written in terms of the Green’s function for half-space geometries. The full nonlinearity of the problem is exploited in order to inspect strong scatterers. After discretization of the continuous model, the resulting equations are solved in a regularization sense by means of a two-step inexact Newton algorithm. The capabilities and limitations of the method are evaluated by means of some numerical simulations.

1. Introduction

Subsurface imaging is an important subject in several applicative areas, including seismic and geophysical prospecting, nondestructive, testing and medical diagnostics [139]. Electromagnetic techniques are widely applied to face this problem. Ground penetrating radar [40] is the basic instrumentation able to retrieve discontinuities in the lower subspace (the ground). Processing of GPR data can allow an improvement in the detection and localization capabilities of the system. In some cases, neural network-based approaches have been adopted [41]. Moreover, inverse-scattering-based techniques have been proposed, too. By using these methods, the measured values of the field scattered by the discontinuities present in the lower half space are “inverted” in order to obtain the spatial distributions of dielectric parameters (e.g., the dielectric permittivity and the electric conductivity) inside a fixed inspection domain. Approaches of this kind have been proposed, for example, for the detection of water leaking from pipes [42]. If an accurate model of the buried scatterers is available, it would be possible to retrieve information about the buried objects in terms of their dielectric parameters.

The present authors developed in [43] an inverse-scattering-based method for buried object imaging under the second-order Born approximation (SOBA) [21, 4448]. By using the SOBA, the reconstruction of buried objects has been obtained with a better accuracy with respect to linearized approximations. In particular, due to the nonlinear nature of the scattering equations, in which the electric field inside the inspected region is an unknown quantity as well as the spatial distributions of the dielectric properties of the medium, the SOBA allows to approximately retrieve the field distribution, too.

However, in order to detect strong discontinuities, the exact nonlinear inverse scattering formulation must be considered. This approach has been faced in [49, 50], in which strong scatterers have been reconstructed in free space under tomographic imaging conditions. The approach proposed in [49], which is based on the use an inexact Newton method, has been experimentally tested in [51]. The same method has been successfully applied in [52] within the contrast source formulation (under free-space conditions).

In the present paper, the approach developed in [43] for buried object detection under the SOBA is extended to treat strong scatterers. In this case, the exact nonlinear inverse scattering formulation, in terms of the half-space Green’s function, is taken into account. In particular, starting from the integral equations for the scattered field for cylindrical objects, slides of the spatial distributions of the dielectric parameters are reconstructed by numerically solving these equations in a regularization sense.

The paper is organized as follows. In Section 2 the mathematical formulation of the approach is reported. In Section 3, the application of the inexact Newton method is described. Numerical results are reported in Section 4, and conclusions are drawn in Section 5.

2. Mathematical Formulation

Let us consider a cylindrical object buried in a homogeneous-half-space medium (Figure 1). A 2D transverse-magnetic (TM) configuration is assumed. The dielectric properties of the buried object are assumed to be independent of the 𝑧 coordinate. The transmitting and receiving antennas are located in the transverse plane and generate a 𝑧-polarized electric field in both the half spaces (this field is independent of the axial coordinate, too).

The dielectric characteristics of the investigation domain S, in the lower half space, are described by the so-called object function, which is defined as 𝜏𝜖(𝐫)=𝑗𝜔(𝐫)𝜖𝑏,(1) where 𝜖(𝐫) is the complex dielectric permittivity at point r, and 𝜖𝑏 denotes the dielectric permittivity of the soil.

The investigation area is illuminated by 𝑉 line-current sources, with angular frequency ω (in the following, the 𝑒𝑗𝜔𝑡 time dependence is assumed and omitted). A mixed illumination configuration (partially borehole and partially from the top), which allows direct comparisons with the results reported in [43], is exploited. Under such assumptions, also the scattered electric field is 𝑧-polarized and independent of 𝑧. Consequently, only the transverse plane can be considered for the imaging procedure.

The scattered electric field, 𝑢(𝑣)scatt, 𝑣=1,,𝑉, is collected for any incident wave in 𝑀 measurement points 𝐬𝑚(𝑣), which constitute the so-called observation domain 𝑂. The dielectric properties of 𝑆 can be related to the scattered field by means of the following set of nonlinear integral equations: 𝑢(𝑣)scatt(𝐫)=𝐺upperlower𝜏𝑢(𝑣)tot𝑢(𝐫),𝐫𝑂,(𝑣)tot(𝐫)=𝑢(𝑣)inc(𝐫)+𝐺lowerlower𝜏𝑢(𝑣)tot(𝐫),𝐫𝑆,(2) where 𝐺upperlower()(𝐫)=𝑗𝜔𝜇0𝑆()𝑔HS(𝐫,𝐫)𝑑𝐫, 𝐫𝑂 and 𝐺lowerlower()(𝐫)=𝑗𝜔𝜇0𝑆()𝑔HS(𝐫,𝐫)𝑑𝐫, 𝐫𝑆 are the data and state operators, where 𝑔HS is the half-space Green’s function [43] (superscripts in the operators recall where they are calculated, whereas subscripts recall where they are integrated).

Equations (2) can be rewritten in a compact operator form as HS(𝐱)=𝐲,(3) where 𝐱=[𝜏𝑢(1)tot𝑢(𝑉)tot]𝑡 is an array containing the values of the object function and of the total electric field inside the investigation area for any illumination. In (3), 𝐲=[𝑢(1)scatt𝑢(𝑉)scatt𝑢(1)inc𝑢(𝑉)inc]𝑡 is an array containing the values of the scattered electric field acquired by the measurement probes in the upper half space and the values of the incident electric field inside the investigation area. The half-space operator HS can be expressed, according to (2), as 𝐺HS(𝐱)=upperlower𝜏𝑢(1)tot𝐺upperlower𝜏𝑢(𝑉)tot𝑢(1)tot𝐺lowerlower𝜏𝑢(1)tot𝑢(𝑉)tot𝐺lowerlower𝜏𝑢(𝑉)tot.(4) As it is well known, the inverse problem (4) is non-linear and ill-posed. Consequently, nonlinear regularizing methods must be used.

3. Regularized Solution of (3)

Due to the mentioned ill posedness, (3) must be solved in a regularized sense. The inexact Newton regularizing scheme developed in [49, 50] is applied here. This method includes two nested loops, in which a linearization (with respect to the current solution guess 𝐱𝑘) of the operator HS is constructed (outer loop) and solved by means of the truncated Landweber method (inner loop). Exactly as in [49], the algorithm performs the following steps.(1) For 𝑘=0 (k denotes the outer iteration index), a starting guess 𝐱0 is chosen (initialization phase).(2) The following linearized version of equation is constructed HS𝐱𝑘𝐡=𝐲𝑘,(5)where 𝐡=[𝜏𝑢(1)𝑢(𝑉)]𝑡, 𝐲𝑘=𝐲HS(𝐱𝑘), and HS𝐱𝑘 is the Fréchet derivative of HS at point 𝐱𝑘. A regularized solution ̂𝐡𝑘 to (5) is obtained by applying an iterative truncated Landweber algorithm, for which we have 𝐡𝑘,0𝐡=0,𝑘,𝑙+1=𝐡𝑘,𝑙𝛾HS𝐱𝑘𝐻𝑆𝐱𝑘𝐡𝑘,𝑙𝐲𝑘,(6)where HS𝐱𝑘 is the adjoint operator of HS𝐱𝑘.(3) The current solution is updated by using the following relation: 𝐱𝑘+1=𝐱𝑘+̂𝐡𝑘.(7)(4) A predefined stopping rule is finally applied. If this criterion is not satisfied, a fixed maximum number of outer iterations, 𝑘max, is performed. In both cases, at any iteration, the procedure can be stopped, otherwise steps (2)(4) are repeated.

It should be noted that step (2) is actually the inner loop of the algorithm, whereas steps (1)(3) constitute the outer loop.

4. Numerical Results

The developed method has been tested by means of numerical simulations. The working frequency has been set equal to 300 MHz (i.e., similar to that used in ground penetrating radars). A homogeneous lossy soil characterized by a relative dielectric permittivity 𝜖𝑟𝑏=4 and an electric conductivity 𝜎𝑏=0.01S/m has been considered (in order to simulate a dry sandy soil [53]). Similar to [43], a mixed measurement configuration is assumed, where reflection and cross-borehole operating modes are mixed together in order to acquire as much information as possible. Correspondingly, the antennas are located at points (Cartesian coordinates) 𝐫ant𝑛=23𝜆04,3𝜆0+(𝑛1)2𝜆0221,𝑛=1,,15,3𝜆0+(𝑛15)2𝜆0221,0,𝑛=16,,28,3𝜆0,(𝑛29)2𝜆021,𝑛=29,,43,(8) where 𝜆0 is the free-space wavelength in the upper medium. A subset of these antennas act both as transmitter (TX) and receiver (RX). When an antenna operates in TX mode, all the remaining ones collect the scattered electric field (i.e., they work in RX mode). Consequently, for every illumination, 𝑀=42 field samples are collected. In particular, in the cases reported in the following, 𝑉=7 antennas are used. They are located at positions denoted by the indexes given by 𝑛=7(𝑣1)+1,𝑣=1,,𝑉.

The investigation area is a square domain of dimensions 𝜆0×𝜆0, whose upper side coincides with the air-ground interface. In order to numerically solve the equations involved in the inverse-scattering problem, such domain has been discretized into 𝑁=20×20 square subdomains.

The scattered field values used as input data for the considered simulations have been numerically obtained by using a computational code based on the method of moments [54]. In order to avoid the so-called inverse crime, a finer discretization (such as the ratio of the cell dimensions in the direct and inverse problems is not an integer number) has been used for the computation of the synthetic data. Moreover, in order to simulate a more realistic scenario, a Gaussian noise with zero mean value and variance corresponding to a signal-to-noise ratio SNR=25dB has also been added to the computed data.

In all the simulations reported in this section, the initialization phase has been always performed by assuming an empty investigation domain in the lower half space. Consequently, in this phase, the object function is equal to zero, and the total electric field is assumed equal to the one transmitted in the lower half space.

The first reported case concerns the reconstruction of a single circular cylinder of radius 𝑎=0.2𝜆0 centered at 𝐫obj=(0.2𝜆0,0.5𝜆0). The dielectric properties of the cylinder are 𝜖𝑟obj=5.6 and 𝜎obj=0.01S/m. The parameters of the inexact Newton algorithm are the following: maximum number of inner iterations, 𝑘LWmax=30; maximum number of outer iterations, 𝑘INmax=100. Figure 2 reports the behavior of two error parameters for different values of the number of outer iterations of the iterative procedure. In particular, the two parameters are defined as 𝑟data=𝑢actual𝑢calc2𝑢actual2,𝑟solution=𝜏actual𝜏rec2𝜏actual2,(9) with 𝜏actual and 𝜏rec being the actual and reconstructed permittivity distributions (at the current iteration) and 𝑢actual and 𝑢calc the measured and computed values of the electric field data (at the current iteration, too). As can be seen from Figure 2, the plot of the “residual” on the data (the error parameter 𝑟data) decreases monotonically as expected. Moreover, the plot of the solution “residual” (the error parameter 𝑟solution) shows a typical semiconvergence behavior, that is, after a certain value of the iteration number the reconstruction quality degrades [55, 56].

The final image of the distribution of the relative dielectric permittivity inside the investigation region (corresponding to the optimum iteration number) is shown in Figure 3. As can be seen, the method is able to correctly retrieve the dielectric profile of the investigation area. In particular, both the shape of the buried inclusion and the value of its dielectric permittivity are correctly identified.

For the present case, an analysis of the effects of the noise has also been performed. In particular, the mean relative error on the reconstruction, defined as the 𝐿1-norm 1𝑒=𝑁𝑁𝑛=1||𝜖actual𝑛𝜖rec𝑛||||𝜖actual𝑛||(10) has been calculated. In (10), 𝜖actual𝑛 and 𝜖rec𝑛 denote the values of the original and reconstructed complex permittivity inside the 𝑛th discretization subdomain of the investigation region. Moreover, the optimal value of the outer iteration have been evaluated, too. The obtained results are reported in Table 1. As expected, as the noise level increases, the reconstruction “quality” decreases. However, even in presence of strong noise levels, the method is able to yield rather good reconstructions (e.g., the error parameter 𝑒 is equal to about 8% error for SNR=15dB).

In the second set of simulations, the reconstruction of a more complex configuration has been considered. In particular, the presence of a double-layer circular cylinder has been assumed. The outer layer has radius 𝑎1=0.2𝜆0 and is characterized by a relative dielectric permittivity 𝜖𝑟1=5.6 and an electric conductivity 𝜎1=0.01 S/m, whereas the inner one has a radius 𝑎2=0.1𝜆0 and is characterized by a relative dielectric permittivity 𝜖𝑟2=8 and an electric conductivity 𝜎2=0.01 S/m. The object is centered at 𝐫obj=(0.2𝜆0,0.5𝜆0). The parameters of the inexact Newton method are the same as the ones adopted in the first case. The reconstructed distribution of the relative dielectric permittivity (final image) is shown in Figure 4. As can be seen, the object is correctly located, and the values of the dielectric properties are identified with a quite good accuracy.

Finally, the reconstruction of two separate circular cylinders have been considered. The radii of the two cylinders are 𝑎1=0.2𝜆0 and 𝑎2=0.15𝜆0. They are buried with centers located at points 𝐫obj1=(0.2𝜆0,0.5𝜆0) and 𝐫obj2=(0.25𝜆0,0.75𝜆0). The dielectric properties of the two scatterers are the following: 𝜖𝑟1=5.6, 𝜖𝑟2=8, and 𝜎1=𝜎2=0.01 S/m. The signal-to-noise ratio is SNR=25𝑑𝐵. The parameters of the iterative inversion algorithm are the same as in the previous cases. The reconstructed distribution of the relative dielectric permittivity is shown in Figure 5. This image corresponds to the optimal solution. As can be seen, also in this quite complex case, the method is able to correctly reconstruct the investigated scene. In particular, the two objects are clearly identified and separated. Moreover, both shapes and dielectric properties of the two targets are estimated with quite good accuracies (although for the second cylinder the dielectric permittivity is slightly underestimated). In order to better assess the reconstruction accuracy, Figure 6 gives the profiles of the retrieved distribution of Figure 5 along some cuts parallel to the coordinate axes.

5. Conclusion

In this paper, a previously proposed approach for the imaging of buried objects has been extended to inspect strong scatterers, for which the exact nonlinear equations of the inverse scattering problem must be considered. The developed method is based on a two-step inexact Newton method, which allows us to obtain a regularized solution of the (discretized) integral equations relating the samples of scattered electric field collected in the upper half space to the distributions of the dielectric parameters inside the test domain located in the ground. The reported preliminary results of numerical simulations seem to indicate the validity of the proposed approach, which needs further evaluation assessments, including an experimental validation, which will be the subject of future research activity in this area.