Abstract
This work concerns an inverse time-dependent electromagnetic scattering problem of imaging internal defects in a homogeneous isotropic medium. The position and cross section of the defects are detected by transient electromagnetic pulses in the case of TE polarization. We apply the Kirchhoff migration scheme to locate the position of small objects from both synthetic and experimental data. The multiple-input-multiple-out scheme is used to recover extended scatterers from the data generated by the software GprMax. Numerical experiments show that the Kirchhoff migration method is not only efficient but also robust with respect to polluted data at high noise levels. Experimental results show good quantitative agreement with numerical simulations.
1. Introduction
Inverse scattering is to recover physical and geometrical information of inaccessible objects from scattered fields measured outside. It has been one of the most challenging problems with considerable practical applications in many areas of technology such as nondestructive evaluation, subsurface and ground-penetrating radar, geophysical remote sensing, medical imaging, seismology, and target identification; see [1–5]. One of the most prominent features of inverse scattering problem is its noninvasiveness, along with the affordability due to cheap nonionizing sensors. However, solving inverse scattering problems is difficult due to the inherent ill-posedness and nonlinearity. Small variation in the measured data can lead to large errors in the reconstruction of the scatterer, unless regularization methods are used. Extensive investigation has been carried out and a variety of inversion algorithms have been proposed.
In the time-harmonic regime, Kirsch and Kress [6, 7], Colton and Monk [8, 9], and Angrll, Kleiman, and Roach [10] theoretically separate the nonlinearity and the ill-posedness of inverse obstacle problems, giving rise to the decomposition method. Muhammad [11], Liang Yang [12], and Hui H et al. [13] employed iterative methods to solve inverse obstacle problems. In these optimization-based iterative approaches, an efficient forward solver is needed for each iteration, and good a priori information might be required in order to choose an initial guess that ensures numerical convergence. Noniterative sampling methods were intensively studied over the last twenty years, for instance, linear sampling method [14], factorization method [15], enclosure method [16], and singular point source method [17]. The key ingredient of the sampling approach is to design an appropriate indicator function from measured data for characterizing the region occupied by the scatterer. Forward solvers are not needed in the process of inversion. The above-mentioned approaches are mostly applicable to measured data irradiated at a fixed energy with many incoming directions. A recursive linearization method with multifrequency data was investigated in [18–22].
There exist a number of reconstruction algorithms in the time domain with dynamic measurement data. Here we mention the time reversal techniques [23, 24], the reversed time migration [25], and the boundary control method [26]. For time-dependent sampling type methods, we refer to the point source method [27], the enclosure method [28], and the total focusing method (TFM) arising from nondestructive evaluation [29, 30]. In one of the author's previous works [31], the TFM (which is also known as the Kirchhoff migration approach used in geophysics; see, e.g., [32]) was examined for imaging acoustically extended scatterers in two dimensions and the indicator behavior was mathematically analyzed. This work concerns imaging internal defects in a homogeneous isotropic medium. It is shown that the TFM is not only efficient but also robust to measurement noise. The aim of this paper is to test the TFM for inverse electromagnetic scattering problem of imaging perfectly conducting cylinders buried in a half-space homogeneous isotropic background medium. Such kind of inverse problems has many applications in searching for internal defects with radar techniques like ground-penetrating radar (GPR); see, e.g., [5, 33, 34]. We note that the light speed is much faster than the sound speed, leading to additional difficulties in computational simulation [35]. In this paper, we make use of the open software GprMax [35, 36] to generate the forward scattering data. The GprMax is an electromagnetic simulation tool based on the Finite-Difference Time-Domain (FDTD) approach coupled with the perfectly matched layer (PML) technique. In our experiments, the real data are gained by GPR.
In this paper, we apply the TFM (or Kirchhoff migration scheme) to inverse electromagnetic scattering problems. In comparison to other algorithms, the implementation of TFM is quite simple, since the image is formed through a superposition of the scattered signals irradiated by each transducer. In fact, our imaging functions are explicit and involve only integral calculations on the measurement surface. Hence, our inversion algorithm is totally “direct”. This also explains why this method is very robust with respect to measurement noise at high levels. Although the inverse scattering problem is difficult due to the inherent ill-posedness and nonlinearity, we do not need to solve the ill-posed and nonlinear problems, even without approximation and iteration. In this work, experiments with synthetic and experimental data are conducted to demonstrate the feasibility and applicability of this approach in the TE polarization case of inverse electromagnetic scattering problems.
The remaining part of this paper is organized as follows. In Section 2 we rigorously formulate the forward and inverse electromagnetic scattering problems. Section 3 is devoted to the description and implementation of the TFM for imaging small scatterers from both synthetic and real data. Extended scatterers will be reconstructed in Section 4 by using GprMax simulated data. Conclusion follows in Section 5.
2. Problem Setting
Consider the electromagnetic wave propagation in a homogeneous isotropic medium. This background medium can be characterized by the electric permittivity , magnetic permeability , and electric conductivity , all of which are assumed to be constant. The propagation of the electric and magnetic fields is governed by the Maxwell system.Assume that is an infinitely long perfectly conducting cylinder buried in the lower half-space , which remains invariant in -direction; see Figure 1(a). The medium in the exterior of D is assumed to be dielectric, which means that . The cross section of in the -plane is denoted by . Incident cylindrical source waves emitted from the surface will be utilized for the purpose of detecting the position and shape of . Throughout the paper, we assume that the incident wave with , , is a filamentary current generated at the source location with a temporal pulse signal . The temporal function is supposed to have a compact support for some . Moreover, we consider the Transverse Electric (TE) case of the monopole with the polarization direction . Then, the incident electric and magnetic waves are governed by the systemfor all , , and , where is the Dirac distribution. Denote by the scattered fields and write the total field as , . Then, we havetogether with the zero initial condition at in . Here is the unit outward normal to the boundary . Eliminating the magnetic field, we arrive atIn the TM polarization case, the electric field takes the form . Hence, the perfectly conducting boundary condition can be written asimplying that on . Usingwe deduce the reduced wave equation from the Maxwell system (5) as follows:In the particular case that is a -periodic function, where denotes the frequency, the incident field can be explicitly represented asHere is the Hankel function of first kind of order zero; see [37, Chapter 3.4].
(a) Geometrical model
(b) Sample grid subdivision
In this paper we shall consider the following inverse problem: (IP): determine the position and shape of from knowledge of irradiated by incident dipoles located at .
Throughout this paper, we use the temporal function of the formwhich is known as the Ricker pulse. Here is the center frequency and 0 is the amplitude. The incident pulse function with GHz is depicted in Figure 2.
We shall consider two cases:
(i) is a small object (compared to the incident wavelength), is a finite line segment, and the data are given by for all .
(ii) is an extended scatterer, for some such that , and the data are given by for all irradiated by many dipoles emitted from .
Note that in case (i), the positions of the receivers and transmitters are identical, leading to the so-called inverse back scattering problem, whereas the case (ii) corresponds to a multiple-input-multiple-output (MIMO) radar system. The aim of this paper is to apply the total focusing method to the above inverse electromagnetic scattering problems with synthetic and experimental data. The TFM is sometimes described as the “gold standard” in the classical beamforming. Both TFM and the Kirchhoff migration form the image with the superposition of the scattered signals irradiated by each transmitter.
3. Imaging Small Scatterers
In this section we suppose that a single metal stick is buried in sand below the ground plane . The cross section of the metal is so small compared to the wavelength such that can be regarded as a point-like object. The electric dipoles are generated on part of the ground plane for some . The Kirchhoff migration imaging scheme for inverse back scattering consists of the following three steps:
Firstly, we take the searching area as a rectangular domain in the lower half-space and divide it into an grid. The grid points, as shown in Figure 1(b), will be referred to as sampling points . We shall design an imaging function defined on these grid points.
Secondly, the filamentary current is transmitted from the source position and the scattered signals are also recorded at the same position , which we denote by for . Note that stands for the scattered wave field. Calculate the flight time of the signal emitted from to a sampling point and then returned back to the receiver , given byHere denotes the light speed in the background medium. In particular, we have in terms of the light speed in the vacuum.
Thirdly, collect all signals for and plot the indicator functionHere we have assumed that the terminal is large enough. By the mathematical analysis performed in [21], the local maximizers of can represent the location of the small scatterer.
3.1. Imaging with Synthetic Data
To implement the Kirchhoff migration scheme, we assume that the interface between the air and the sand is given by with . The dielectric properties of the sand are considered as 4. We choose receivers and transmitters uniformly lying between and . The cross section of the metal stick is supposed to be the circle of radius 0.009 centered at ; see Figure 3 for the locations of the metal stick in the -plane and the first transmitter (receiver). The center frequency of Ricker waves is uniformly taken as and the amplitude is set to be . The terminal time is taken aswhere denotes the travel distance of the radar. In our simulations we take .
The scattered data were generated by GprMax based on the FDTD method (see [20, 38]). The unbounded exterior domain is truncated by an absorbing boundary condition. The mesh of the forward solver is successively refined until the relative error of the successive measured scattered data is below .
For the -th transmitter, the scattered data are recorded at 1063 time points given byLinear interpolation is adopted to obtain the value of for an arbitrary . In Figure 4, we show the signals of the incident electric field and the total fields irradiated by the first transmitter. The incident field obviously has a launch waveform, as shown in Figure 4(a). The significant changes before are due to the transmitter itself rather than the reflected wave field from the metal stick. In Figure 4(b), one can observe the changes of the signal at about caused by the presence of the object. Subtracting the incident field from the total field, one can obtain the scattered signal, which will be utilized to locate the position of the stick.
(a) Incident field
(b) Total field
For with , , we first plot the functionwhich forms the B-scan image of the object. Note that only the information of the signal of is involved in B-scan for , whereas all signals are required in our indicator equation (15). Figure 5(a) shows the B-scan image of the small object located at (0.2, 0.167). The -coordinate of the location corresponds to the trace number 100. In Figure 5(b), one can see an approximate location of the target using our indicator (15). Numerical experiments show that the Kirchhoff migration method is not only efficient but also robust with respect to polluted data at high noise levels (Figure 6). Note that in this paper the scattered near-field data is polluted bywhere is the noise ratio and is the standard normal distribution.
(a) B-scan image formed by the indicator equation (15)
(b) The Kirchhoff migration image
(a) B-scan image formed by the indicator equation (15)
(b) The Kirchhoff migration image
3.2. Imaging with Experimental Data
In our experiment settings, a metal stick is detected in a sand box by a GSSI Mini HR hand-held ground-penetrating radar; see Figure 7 for the experiment site. The size of the sand box is 0.45 m×0.3m×0.2m, the thickness of the wooden shell is 0.015m, the dielectric properties of the sand are 4, and the dielectric properties of the wooden box are considered approximately the same as the sand. For electric dipoles with TE polarization, the three-dimensional Maxwell equation can be reduced to a two-dimensional model for scalar wave equations. We set up a two-dimensional coordinate system and take the initial position of the source and receiver of the radar at (0.1, 0.1). The metal stick is centered at (0.2, 0.167) with the radius 0.009. The radar rolls along the straight line = 0.1 on the wooden shell and moves forward in the positive -direction at an even pace. In this process there are totally 268 electric dipoles excited by the radar before it reaches the ending point (0.25, 0.1). The temporal function of the incident wave is a Ricker wavelet excited at the frequency 2.6GHz with the amplitude 1.
The B-scan image formed by our experimental data is plotted with Matlab R2013b, as shown in Figure 8(a). Figure 8(b) shows the image using the indicator (12), where the position of the target is precisely located. In comparison with the counterpart from synthetic data (see Figure 5(b)), we find that the inversion of the depth of the stick contains relatively large errors.
(a) B-scan image
(b) Imaging by using the Kirchhoff migration
Next, we want to find two metal sticks with different positions. The centers of the cross section in the -plane are given by (0.15, 0.15) and (0.25, 0.17), respectively. The other parameters in our experiment are the same as before. See Figure 9(b) for the experiment site and Figure 9(a) for the location of the two sticks.
(a) Locations of the first transmitter/receiver and the two small targets
(b) Experiment site
The B-scan images with synthetic and experimental data are shown in Figures 10(a) and 10(b), respectively. The final inversion images based on the Kirchhoff migration are illustrated in Figures 11(a) and 11(b), respectively.
(a) Synthetic data at the noise level 20%
(b) Experimental data
(a) Synthetic data at the noise level 20%
(b) Experimental data
As shown in Figure 10, the B-scan results in both simulation and experiment can indicate a local peak in the position of the metal sticks, verifying the correctness of the inversion scheme. In Figure 11, the local maximum values of the indicator function are obtained approximately at (0.15, 0.15) and (0.25, 0.17), which agree well with the positions of the targets in our model. From the experiment image, it can be concluded that the contour image and the position of the left metal stick agree well with the true scatterers, but the size of the right one is not well reconstructed. Further, there is an untrustworthy imaging area between the two metal sticks, maybe due to the offset of the two maximizes nearby.
4. Imaging Extended Scatterers
The main purpose of this section is to test the Kirchhoff migration for imaging extended scatterers from synthetic data of several incident electric dipoles. As in the previous section, the forward data are obtained by the software GprMax. We consider the second case of our inverse electromagnetic scattering problem formulated at the end of Section 2, where the source positions are uniformly distributed on the circle for some , where . Denote by the number of incident dipoles. The Cartesian coordinates of the -th dipole can be formulated asIn all examples we use receivers equally spaced on . For notational convenience, we use to represent the scattered signal excited by the -th transmitter recorded at the -th receiver. The flight time of an impulsive signal from to a sampling point and then back to is then given byIn this case the imaging function based on Kirchhoff migration scheme can be formulated asIn all of our numerical examples to be reported below, we set . The sampling region is set to be the rectangular domain . The homogeneous background medium is set to be air, with the dielectric constant and the electric conductivity . Unless otherwise stated, the excitation frequency of a Ricker wave is taken as . The amplitude is set to be and the time window is .
In the first example, we excite only one electric dipole located at and measure the scattered field at for . Three types of scatterers will be recovered: two small circular metals, a square-like metal, and an L-shaped metal; see Figures 12(a), 12(c), and 12(e). The last two targets will be referred to as extended scatterers, since their size is much bigger than the incident wavelength. The circular metals are centered at and with the radius 0.01 m. They can be treated as point-like scatterers, because their radiuses are both much smaller than the wavelength corresponding to . The coordinates of the corners of the square-like metal are given by , , , and , and those for the L-shaped metal are , , , , , and . Figures 12(b), 12(d), and 12(f) show the reconstructions from the data of one dipole only. It can be seen that a single source can be used to capture the position of point scatterers and the illuminated boundary of the square-like and L-shaped scatterers. However, the entire shape of extended objects cannot be well reconstructed.
(a) Two small objects
(b) Location of the two small objects
(c) A square-like metal
(d) Imaging of the square-like metal
(e) An L-shaped metal
(f) Imaging of the L-shaped metal
In the second example, we increase the number of incident dipoles for reconstructing the entire shape; see Figure 13. The data of four incident waves without noise can generally produce a full-range image of an extended scatterer. However, there are bigger background disturbances and the scatterer cannot be accurately recovered. Using eight and sixteen incoming sources can image the outline of the scatterer. In Figures 13(e) and 13(f), the maximum values of the indicator are almost distributed on the contour of the object and the imaging is better than the case of eight dipoles. This means that sixteen sources are capable of imaging the scatterer, and therefore it is not necessary to use additional sources. Hence, the data of multiple point sources in different directions need to be measured in order to better capture the shape of the target.
(a) = 4
(b) = 4
(c) = 8
(d) = 8
(e) = 16
(f) = 16
In the third example, we discuss the sensitivity of the imaging quality to the center frequency of incident pulse functions. We use M=16 incident dipoles and N=64 receivers. As shown in Figure 14, the shape of an extended object can be imaged at various frequencies from 1GHz to 10GHz. However, at lower frequencies the target cannot be accurately recovered. Increasing the excitation frequency of the source (which means a shorter wavelength) may lead to an image with higher resolutions.
(a) =1GHz
(b) =2GHz
(c) =3GHz
(d) =4GHz
(e) =5GHz
(f) =10GHz
5. Conclusions
We apply the Kirchhoff migration approach to solve inverse scattering problems for time-dependent electromagnetic waves. The inversion scheme involves only integral calculations and is robust to polluted data at high noisy levels. We use both synthetic and experimental data to examine the performance of such a direct sampling scheme for locating point-like scatterers. The multiple-input-multiple-out scheme is used for imaging extended scatterers from the data generated by the software GprMax. Our experiments show that the data of one electric dipole can correctly reconstruct the location of two small scatterers, and that sixteen incoming sources can perfectly reconstruct the shape of an extended target. In addition, the source irradiated at 10GHz can be used to more clearly reconstruct the geometrical shape than the data of lower frequencies ranging from 1GHz to 5GHz. Our future work consists of imaging extended obstacles from experimental data buried in a random background medium, which is more challenging than the problem in a stationary homogeneous medium.
Data Availability
The data used to support the findings of this study are included within the supplementary information files.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The work of Hongwei Zhou and Ling Ma was supported by the National Key Research and Development Program (Grant No. 2017YFD0600101), the Fundamental Research Funds for the Central Universities (Grant No. 2572019BF08), the China Postdoctoral Science Foundation (Grant No. 2018M640288), and the Heilongjiang Postdoctoral Fund (Grant No. LBH-Z18004). The work of Guanghui Hu is supported by NSFC (Grant No. 11671028) and NSAF (Grant No. U1530401).
Supplementary Materials
We used the GSSI Mini HR hand-held ground-penetrating radar to scan the sand box in Figure 7. The ground-penetrating radar had a transmitting antenna and a receiving antenna, the transmitting antenna transmitted a Ricker wavelet excited at the frequency 2.6GHz with the amplitude 1, and the receiving antenna received data as in the supplementary information files of “data of fig7.CSV” and “data of fig9b”. The difference between “data of fig7.CSV” and “data of fig9b” is that the file “data of fig7.CSV” is for one metal stick, and the file “data of fig9b” is for two metal sticks. (Supplementary Materials)