An Automatic Framework Using Space-Time Processing and TR-MUSIC for Subsurface and Through-Wall Multitarget Imaging

We present an automatic framework combined space-time signal processing with Time Reversal electromagnetic (EM) inversion for subsurface and through-wall multitarget imaging using electromagnetic waves. This framework is composed of a frequencywavenumber (FK) filter to suppress direct wave and medium bounce, a FK migration algorithm to automatically estimate the number of targets and identify target regions, which can be used to reduce the computational complexity of the following imaging algorithm, and a EM inversion algorithm using Time Reversal Multiple Signal Classification (TR-MUSIC) to reconstruct hidden objects. The feasibility of the framework is demonstrated with simulated data generated by GPRMAX.


Introduction
Electromagnetic waves play a significant role in subsurface and through-wall imaging (TWI) due to their strong penetrability.In addition, systems based on electromagnetic waves have a variety of practical applications including surveillance and reconnaissance, rescue missions, search for suspects and hostages, archaeological applications, and demining [1,2].
In general, the remote sensing problems [3] contain the following five main methodological issues: (1) hidden object detection and localization; (2) medium material recognition; (3) hidden object size estimation; (4) hidden object shape estimation; (5) hidden object material recognition.One approach for the reconstruction of hidden objects is represented by electromagnetic field inversion [4], and an alternative is represented by techniques based on the pattern recognition approach [5].However, most of the research activity has been focused on the hidden object detection and localization issues, and the hidden object imaging is the signification step.So many imaging algorithms have been proposed in the past few years.In this paper, we will focus on the hidden object detection and localization.
Previously, several effective imaging algorithms that take into account the medium reflection, signal distortion, and delay effects have been proposed, such as subspace based method [2], nonlinear inverse scattering algorithms, and linear inverse scattering algorithms [6].The nonlinear inverse scattering algorithms can perform a quantitative reconstruction of the hidden targets, however, which are still limited to implement in practical applications due to its iterative computation thus making it very time-consuming [7].Aiming to this problem, the linear inverse scattering algorithms have been proposed to implement the detection and imaging of the hidden object.Linearized inversion schemes based on diffraction tomography, which require much less computational efforts, are well suited for in practical applications since they can be implemented with fast Fourier transform [8].However, these imaging algorithms have not described the preprocessing step, nor have taken into account the multiple targets scenario.To address the above problems, in this paper, we would present the preprocessing techniques and consider the multitarget scenario and we will adopt the subspace-based method.
The detection and location of the hidden objects are determined by applying inversion schemes to the scattering data.Nevertheless, the quality of imaging largely depends on a number of interference factors among the raw data [9].The major factors are (1) the antenna direct wave and medium bounce; (2) the migration of targets information due to medium roughness; (3) the diffracted wave on the lithology mutations point of the radar profile.Therefore, International Journal of Antennas and Propagation there is a need to adopt some appropriate data processing techniques to suppress the interference signal and enhance the reflected wave from the target.However, this is a quite challenging task because of the complicated EM scattering phenomenon.
The scattering signal will overlap between the target scattering signal and the antenna direct wave when the target object is located in the shallow underground region or the region near the building wall.In such circumstances, the direct wave would have a great influence on the object scattering signal.Therefore, this is a key step to suppress direct wave.In order to solve this problem, many signal processing approaches have been developed, and most of them have a good effect.They include advanced approaches for median filter, average filter [10], wavelet transform [11]; but their performance needs a large number of data or a big computational cost.So we propose to apply a FK filter approach to suppress the influence of the direct wave in this paper.An FK filter is a 2-dimensional process in frequency-Wavenumber domain and utilizes a relationship between the velocity of the incidence along a spatial sampling line and the slope of the spectrum in the frequency-wavenumber domain.It will implement a separation of incident waves that arrive from different directions by using this relationship [12].
Generally, the radar profile is different from geological space because of some factors, such as wave propagation characteristics, the choice of observation system, underground structure, and radar wave recording mode.To remove the influence of factors mentioned above, there is a need to adopt migration methods for the purpose of hidden objects location and imaging.At present, the migration methods in data processing technology can be classified into two main classes.The first one exploits the wave equation to back-propagate the recorded wave field with the aim of achieving migration imaging, such as FK migration [13][14][15], Kirchhoff migration [16], and finite difference migration [17].The second class exploits the ray theory with the aim of making the reflected wave automatically return to its true spatial location.
In this paper, we present an automatic framework combined space-time signal processing with Time Reversal electromagnetic (EM) inversion for subsurface and throughwall multitarget imaging using electromagnetic waves.In particular, at the beginning, a FK filter approach for filtering out the direct wave from the raw data is used.The underlying processing method with the aim of making the reflected wave automatically returns to its true spatial location by FK migration algorithm.Simultaneously automatically estimate the number of targets and identify target regions, which can be used to reduce the computational complexity of the following imaging algorithm.Subsequently, the TR-MUSIC method is applied for an efficient localization and imaging of multiple hidden targets with a high resolution.
This paper is organized as follows.In Section 2, we describe the preprocessing method.Section 3 presents the formulation of the TR-MUSIC algorithm based on the subsurface and through-wall imaging.Then, several examples of multitarget localization and imaging are presented in Section 4, and the conclusions are made in Section 5.

Preprocessing Technique Theory
2.1.FK Filter Algorithm.Frequency-wavenumber filter is used to suppress interference waves in the frequency-wavenumber spectrum.The radar profile is a 2D signal and can be represented as E(t, x), and the frequency spectrum of the signal is where x is the horizontal coordinate, t is the two-way travel time along the depth direction, and k x is the wavenumber vectors in the horizontal component.Then, use an FK filter as where H(ω, k x ) is a transfer function of filter to be presented in the following and G(ω, k x ) is a filtered spectrum.By applying 2-D inverse Fourier transform to G(ω, k x ), we can obtain a time-spatial radar profile suppressing the interference wave component as [12]: In general, the spectrum component of the direct wave and medium bounce in the FK domain is distributed along the frequency axis, spatial frequency k x = 0.In this paper, the main interference waves of simulation data are the direct wave and medium bounce, so the filter to be applied in the FK domain is given by In the following sections, we present the FK migration algorithm and describe the relative formulas of this method.

Frequency-Wavenumber Migration. FK migration algo-
rithm is used to make the object signal back to its true spatial location and estimate the number of targets and identify target regions.The waves propagate at velocity v in the medium.For practical scenario, we assume that the transmitter and the corresponding receiver are collocated.While the waves travel a two-leg path in reality, mathematically we could assume that the waves travel only from the target to the receiver at half velocity v/2.The relative formulas of FK migration algorithm applied to the half velocity v/2 are listed as follows: where v is the velocity of migration and k z is the wavenumber vectors in the vertical component.At the outset, by applying 2-D Fourier transform to E(t, x), we can obtain a frequencywavenumber spatial record data.And then, B(k x , k z ) is yielding through inserting (1) into the relation (6).Finally, we obtain the inverse Fourier transform E(x, z, 0) of B(k x , k z ) by using (7).E(x, z, 0) is the reflection point migrated image:

Mathematical Formulation
3.1.Through-Wall Model.The geometry of the problem is described in Figure 1.Multiple objects are hidden behind a known homogeneous wall with dielectric permittivity ε s , electric conductivity σ s , and thickness d.The first and third layers are freespace with dielectric permittivity ε 0 , and the magnetic permeability is the same everywhere and equal to the one in free space μ 0 .We consider a 2-D geometry where both the wall and the objects are assumed to be infinitely long and invariant along the y direction.The transmitting and receiving antennas are assumed to be ideal dipoles and center at the same position R ant at a distance of z t in front of the homogeneous wall located in the first layer medium, where ant = 1, 2, 3, . . ., N, R ant = (X ant , Z ant ).A set of M targets are located at position r obj in an inaccessible region behind the wall, for obj = 1, 2, 3, . . ., M, r obj = (x obj , z obj ), M < N.
TR-MUSIC algorithm exploits the multistatic response matrix K in order to image the object and is given by [18] where τ obj is the scattering strength of the obj-th target and g obj (ω) is a vector associated with the layered medium Green function of the obj-th target for the homogeneous wall and is given by Due to reciprocity, K ji (ω) = K i j (ω), so that K is a N × N symmetric matrix.The multistatic response matrix K is a function of the geometry and physics of wave propagation between the antenna elements and the objects [19].Where G is the layered medium Green function, the spectrum form of G can be written as [2]: where T is the transmission coefficient of the wall where k 1 , k 2 are the wavenumbers in the freespace and in the wall, respectively.k x is the horizontal component of the wavenumber vector, and k 1z , k 2z are the vertical components of the wavenumber vector, k 1 and k 2 , respectively.In order to efficiently solve (10), in this paper, the saddle point method is employed for the asymptotic evaluation.Let then the stationary phase point is given by The solution of the above equation is the stationary phase point k x0 .Assuming that Φ(k x ) is a high oscillating function and by employing the saddle point method, the layered medium Green function can be derived as [18]: where Φ is the second derivative of Φ, As long as M < N, the N × N multistatic response matrix K is rank deficient.This characteristic of the response matrix forms the basis for applying subspace-based signal processing techniques, so that we exploit the singular value decomposition (SVD) of the response matrix, and is given by where the matrix Σ is diagonal matrix, in which diagonal elements are nonnegative and known as singular values.U, V are N × N orthogonal matrices.Since K is rank deficient, the space of the data is composed of signal subspace and null subspace.In this paper, we assume that there is L nonzero singular value of the response matrix, videlicet; the first L columns of U correspond to the vector g obj (ω), while the remaining N − L columns span the null subspace.Then, the locations of the objects can be determined through the timereversal MUSIC pseudospectrum where • denotes the inner product and u p is the pth column vector of the matrix U. When r matches with the true location of the target and it shows a peak at the target location.
The mathematical formulation of the subsurface model is presented in the following section, and it is noted that the spectrum form of the Green function is changed.

Subsurface Model.
The aim of this section is to describe the mathematical formulation of the multitarget subsurface imaging based on the through-wall imaging theory.The planar interface z = 0 divides the space into two half spaces (see Figure 2).The upper half space is air with dielectric permittivity and magnetic permeability equal to ε 0 and μ 0 , respectively.The lower half space is homogeneous soil with relative dielectric permittivity ε b , magnetic permeability μ 0 , and conductivity σ b .The transceiver antenna is located in the upper half space at a distance of D1 from the ground surface with an interelement spacing D2, while the targets are buried in the lower half space.In this section, we still consider the 2-D structure where the objects are infinitely long and invariant along the y direction.Furthermore, we assume that the location of the targets and antennas based on the subsurface imaging is the same as the through-wall model.
In the mathematical formulation of the through-wall imaging, we consider the Green function form of the threelayered background medium (air-wall-air).In the following, the Green function for the two-layered background medium (air-soil) is applied for the subsurface imaging, so that the spectrum form of the Green function can be rewritten as [20] G R ant , r obj , ω where T(k x ) is the transmission coefficient, where k 1 , k 2 are the wavenumbers in the freespace and in the homogeneous soil, respectively.The same step, the saddle point method is employed for an efficient evaluation of the above integral equation, and the Green function can be derived as where k x0 is the stationary phase point, where Φ is the second derivative of Φ, and the other steps are the same as previous through-wall mode.

Results and Discussion
In this section, some numerical simulation results are presented in order to show the effectiveness of the automatic algorithm.The reference scenario and the measurement configuration are the same as already defined in the preceding section.The exact scattered field is computed in time domain by means of the GPRMAX software, and a line source, which radiates a Ricker wavelet with a central frequency equal to 1 GHz, is used.

Through-Wall Case.
In order to show the effectiveness of the automatic algorithm for through-wall radar imaging (TWRI), a numerical result for the imaging of two targets behind a homogeneous wall is presented in this section.The measurement configuration for the system is shown in Figure 1.
The radar system scans the region of interest along a line parallel to the wall in x direction form −1 m to 1 m with a step 0.02 m at a distance of 0.1 m from the front wall.The dielectric constant, conductivity, and thickness of the wall are ε s = 4, σ s = 0.01 S/m, and d = 0.2 m.

Preprocessing Results of Through-Wall Model. Figure 3(a) depicts the raw radar data in the presence of additive
Gaussian noise with a signal-to-noise ratio (SNR) of 40 dB for two targets in the time domain.The first square object has a relative dielectric permittivity of 9 and a conductivity of 0.05 S/m.And the second has a relative dielectric permittivity of 12 and a conductivity of 0.05 S/m.The two square targets whose sizes are 0.05 m by 0.05 m are centered at (−0.4 m, −1.2 m) and (0.4 m, −0.9 m), respectively.Between about 0 and 1.2 ns, we can observe the multiple reflections occurring in the antenna.These reflections interfere seriously in the focus of the target reflections.Due to the transceiver closed to the homogeneous wall, the direct wave and the front wall echoes almost reach the transceiver simultaneously.At around 4.3 ns, a second-order reflection between the antenna and the homogeneous wall is visible.Then, the hyperbola originating from the obj1 target appears at around 11.9 ns, while the hyperbola of the obj2 target is appearing at about 8.9 ns.
In Figure 3(a), compared with the direct wave, the curves presented by targets are relatively weak.Therefore, in order to suppress the influences of the direct wave and medium bounce, the FK filter method is applied for filtering out the major antenna effects and the homogeneous wall surface reflection.In Figure 3(b), the 2-D radar data after filtering are presented in the time domain.We can see that the antenna direct wave is effectively removed and the targets reflected waves are enhanced.The remaining oscillations visible in the processed data are an artifact of the inverse Fourier transform applied to data collected within a limited frequency range.could find that these objects are clearly identified and we also can automatically estimate the number of targets, which can be used to reduce the computational complexity of the following imaging algorithm.Moreover, it can be seen that the shape of the buried objects along the scanning axis is distinguished.It would help us to extract the information about the position, depth, and shape of the buried objects from the processed data.After the preprocessing step, the imaging results of the multitarget are shown in the next section.

TR-MUSIC Reconstruction Results of Through-Wall
Model.In this section, we are devoted to presenting the reconstruction results of the TR-MUSIC algorithm when applied to the scattered field data achieved by the proposed filtering procedure.In the simulations, images of certain canonical targets are shown revealing their region of localization.For reducing the computational complexity, attention is restricted in a limited region, which is estimated by the migration image.And then, we assume a 2-D geometry where the wall is homogeneous medium and the model ignores the multiple scattering between the targets.In order to show the contribution of FK migration algorithm, the next, we will show that the TR-MUSIC reconstruction results utilize the scattered field data without and with FK migration procedure, respectively.The first sets of reconstructions are concerned with the eigenvalues of the multistatic response matrix of two targets.Figure 4(a) plots the eigenvalues of the multistatic response matrix for the two targets without migration procedure.For the computation of the TR-MUSIC pseudospectrum, the break point is chosen as the fifth point where the other singular values drop to zero. Figure 4(b) depicts the eigenvalues distribution of the multistatic response matrix for the two targets after migration procedure, and the break point is also chosen as the fifth point.By using the eigenvectors corresponding to the eigenvalues after the break point, the TR-MUSIC imaging results of the two targets are shown in Figure 5. Figure 5(a) shows the TR-MUSIC reconstruction image of the throughwall model adopted directly on the scattered field data without FK migration procedure.As it can be seen, the image is not an accurate reconstruction of the two objects.By using the scattered data after migration procedure, the imaging result of the two objects is shown in Figure 5(b).From this figure, we could see that the automatic framework is able to locate and determine the extended targets.

Subsurface Imaging
Case.This section is devoted to present some numerical tests with the aim of confirming the effectiveness and performance of the automatic algorithm as applied to subsurface imaging problems.The measurement configuration is taking into account the interface and the losses in the soil (background medium).The soil has a relative dielectric permittivity of 3 and a conductivity of 0.01 S/m, and we assume that the soil is homogeneous.With respect to the through-wall model, the subsurface imaging model is taking into account a two-layered background medium.The measurement configuration for the system is shown in Figure 2. The radar system scans the region of interest along a line parallel to the air-soil interface in x direction form −1 m to 1 m with a step 0.02 m at a distance of 0.1 m from the air/soil interface.

Preprocessing Results of Subsurface
Model.The 2-D raw ground penetrating radar simulation data in the presence of additive Gaussian noise with a signal-to-noise ratio (SNR) of 45 dB for the buried objects in the time domain is presented in Figure 6(a).It can be seen that the multiple reflections occurring in the antenna are appearing in between about 0 and 1.0 ns.These reflections distort the signal for later times as well.At around 2.0 ns, the soil surface reflection is visible.Then, the hyperbola originating from the obj1 target appears at around 11.6 ns, while the hyperbola of the obj2 target is appearing at about 8.2 ns.The first square object has a relative dielectric permittivity of 8 and a conductivity of 0.005 S/m.It is sized 0.05 m by 0.05 m, the abscissa of its center is −0.5 m, and the depth of its center is −0.8 m, while the second circular cylinder target with diameter 0.05 m is centered at (0.5 m, −0.5 m).Its relative dielectric permittivity and conductivity are 12 and 0.005 S/m, respectively.Moreover, it can be noted that, with respect to the antenna direct wave, the buried objects reflected waves are relatively weak.So there is a need to adopt some appropriate data processing techniques to suppress the antenna direct wave and enhance the reflected wave from the targets.In Figure 6(b), the 2-D radar data after applying the FK filter method is presented in the time domain.It can be noted that the antenna direct wave is effectively removed and the targets reflected waves are enhanced.
The FK migration image result of the two buried objects is depicted in Figure 6(c).With respect to Figure 3(c), there is a good effect of the reconstruction.From this figure, we could find that the upper side and the lower edge of the buried objects are rightly localized.Moreover, it can be seen that the shape of the buried objects along the scanning axis is distinguished.It would help us to extract the information about the position, depth, and shape of the buried objects from the processed data.
In the following section, we present the TR-MUSIC reconstruction results by being applied to the scattered field after the preprocessing procedure.

TR-MUSIC Reconstruction Results of Subsurface Model.
This section is devoted to illustrate the reconstruction results by applying the migration data and without migration data.The reconstructions are performed by the TR-MUSIC algorithm applied to the subsurface model.The first step of the reconstruction aims to choose a break point from the eigenvalues of the multistatic response matrix of the two buried targets.For the two buried objects in this simulation, the singular values are shown in Figure 7. Figure 7(a) plots the eigenvalues of the multistatic response matrix for the two buried targets without migration procedure.For the computation of the TR-MUSIC pseudospectrum, the break point is chosen as the sixth point where the other singular values drop to zero. Figure 7(b) depicts the eigenvalues of the multistatic response matrix for the two buried targets after migration procedure, and the break point is also chosen as the sixth point.Figure 8(a) shows the TR-MUSIC reconstruction image of the subsurface model adopted directly on the scattered field data without FK migration procedure.As it can be seen, the image is not an accurate reconstruction of the two buried objects.By using the scattered data after migration procedure, the imaging result of the two buried objects is shown in Figure 8(b).From this figure, we could see that the automatic framework is able to locate and determine the extended buried targets.

Concluding Remarks
This paper has presented an automatic framework combined space-time signal processing techniques with TR-MUSIC algorithm for subsurface and through-wall multitarget imaging using electromagnetic waves.The space-time signal processing techniques include the FK filter method to suppress the antenna direct wave and medium bounce, and then the FK migration algorithm to automatically estimate the number of targets and identify target regions, which can be used to reduce the computational complexity of the following imaging algorithm.And finally, by employing the null space of the multistatic response matrix, simultaneous detection and localization of multiple targets can be achieved by TR-MUSIC algorithm.The experimental results based on simulated scattered field data show that the presented automatic framework could obtain a good estimation of the positions of the hidden targets, and TR-MUSIC imaging algorithm is not an accurate reconstruction of the hidden objects by applying the scattered field data without FK migration procedure, so that the presented framework is an effective algorithm for locating and determining the hidden targets.However, the scenarios assume that the background medium is homogeneous and it ignores the mutual interactions among these hidden objects.Further research will focus on improving the proposed algorithm to adapt to the inhomogeneous background medium and the mutual interactions scenario.

Figure 3 (Figure 4 :
Figure 4: Eigenvalues of the multistatic response matrix of through-wall imaging targets: (a) on the scattered data without migration; (b) on the scattered data after migration.

Figure 5 :Figure 6 :Figure 7 :
Figure 5: Image obtained with the through-wall radar scattered field data by TR-MUSIC algorithm; colors represent the pseudospectrum value: (a) on the scattered data without migration procedure; (b) on the scattered data after migration procedure.

Figure 8 :
Figure 8: Image obtained with the GPR radar scattered field data by TR-MUSIC algorithm; colors represent the pseudospectrum value: (a) on the scattered data without migration procedure; (b) on the scattered data after migration procedure.