Synthetic aperture radar (SAR) processing techniques are commonly used in radar, sonar, seismic and ultrasound imaging^[1,2]. To reduce the false alarm rate of target identification and improve the detection efficiency, it is desirable to improve the imaging resolution and reduce the computation cost.

Among SAR processing techniques, the migration is a family of effective algorithms originated from seismic imaging^[3]. In these algorithms, the reflection signatures from recorded data are refocused back to the true target location and physical shape. For ground penetrating radar (GPR) employing, conventional SAR imaging techniques, in particular the back-projection (BP) method, can be used to achieve the image refocusing by the coherent summation of the back-scattered data after a proper phase compensation process^[4]. However, it is computationally expensive to apply the “delay-and-sum” approach to compensate the phase delay by solving a fourth order equation^[5].

The use of migration algorithms for GPR imaging has become popular in recent years due to the high computational efficiency of the fast Fourier transform (FFT) which enables low execution time at the cost of increased memory usage^[6,7]. However, for multi-layer medium scenarios like standoff GPR imaging, traditional Omega-k algorithm cannot be directly utilized because it requires the velocity of the electromagnetic (EM) wave to be constant. Gazdag introduced the phase shift migration (PSM) algorithm to GPR, which also operates in the wave-number domain^[7]. The PSM considers the variation of the wave propagation speed in the different mediums. However, this approach is slower than Omega-k because it calculates the wave propagation at every incremental depth repetitively and the FFT is not involved in the depth direction.

In this paper, we extend the original Omega-k algorithm to the multi-layer medium scenario, at the same time, reserving its high computational efficiency. For multi-layer geometry in which layers of homogeneous media are stacked vertically, PSM is first applied recursively to each layer boundary, and then, Omega-k is used within the layer. The algorithm exploits fast Omega-k interpolation in the spatial wave-number domain instead of iterative calculating and can be implemented for multi-layer medium imaging with known each layer depth and permittivity. For multi-layer imaging, the wave velocity estimation is necessary and need to be implemented first for lateral image reconstructions. For the standoff GPRs, methods like Hough Transform or template matching^[8,9] are no longer capable due to strong reflection and other clutter from the interface. Here, we propose a target-region-extraction method to extract the target region first and to estimate the velocity within the region with higher accuracy.

The remainder of this paper is organized as follows: The extended Omega-k (EWK) algorithm for multi-layer medium is introduced in Sect. 1. Section 2 introduces the data validation platform and presents the velocity estimation result. In Sect. 3, various experimental imaging results for buried metallic targets in sand is presented. The conclusion is given in Sect. 4.

Considering the imaging problem illustrated in Fig. 1, in which an antenna array is placed in medium 1 (e.g., air) and the target resides in medium 2 (e.g., sand). At each antenna element, a transmitter illuminates the target region and a coherent receiver measures the reflected signal. Refraction and reflection happen on layer interface due to permittivity differences between layers. An averaging method is applied to reduce sand reflection before image generation.

Migration algorithms are simplified to be an inverse problem with the assumption that the targets are the sources of the received waves. Here, we assume that all scattering targets are located in down half space medium 2, and the back reflected signal is recorded along each sampling position at the height , where and are the horizontal position of the transceiver. The response at the transceiver will simply be the superposition of each point on the target multiplied by the round trip phase to that point. Then, we apply FFT over and is,

Skx,ky,k1,Z1=∬srx',y',k1,Z1e-ikxx'+kyy'dx'dy'

where indicates the array position; is wave-number of first layer, and represents the spatial wave-number of and respectively. We assume that the medium is isotropic and homogeneous. Therefore, we can use a single velocity for the wave propagation within each layer. For a homogeneous and lossless medium, it is obvious that , where is speed of light in free-space and is the relative electric permittivity of the medium. By utilizing the principle of stationary phase^[10], Eq. 1 can be solved as,

Skx,ky,k1,Z1=ρx,y,z∙e-ikxx-ikyy-ikz1z+Z1

where is the target position, is the reflectivity function which is simply the ratio of reflected field to incident filed, and is the projection of the wave number in medium 1 along the z-axis

kz1=4k12-kx2-ky2

Suppose we shift down the transceiver along z-axis by , the FFT at new height is

Skx,ky,k1,Z1+∆z=ρ∙e-ikxx-ikyy-ikz1z+Z1+∆z

By substituting Eq. 2 into Eq. 4, we can summarize the principle of PSM. In the wave-number domain, the wave of any depth of can be extrapolated from depth of by multiplication with a complex exponential:

Skx,ky,k1,Z1+∆z=Skx,ky,k1,Z1×eikz1∆z

Finally, an inverse Fourier transform (IFFT) can be applied to obtain the extrapolated wave field at in space-time coordinates.

Ix,y,Z1+∆z=∭Skx,ky,k1,Z1×eikxx+kyy+kz1∆zdkxdkydk1

By proper selecting and using Eq. 5, we can create the image increment-by-increment until the entire depth is imaged. To precisely generate an image, the increment is generally chosen to be smaller than the range resolution. When for a two-layer problem, the field cannot be extrapolated directly with different , referring to the medium relative permittivity, but the extrapolation method within each layer is similar, which is,

Ix,y,Z2+∆z=∭Skx,ky,k2,Z2×eikxx+kyy+kz2∆zdkxdkydk2

where, and denote the field at depth and , respectively. The transmission factors of different media are generally a complex function, dependent on both incident angle and the impedance of the media^[8]. However, for most common GPRs, the direction of the transducer limits the emitted and received wave fields to a relatively small angle interval. Therefore, we assume that the transmission factors are approximately independent of incident angle, so that the wave fields directly above and below an interface are proportional. Thus, the wave field at the interface can be calculated from the wave field measured at as follows:

Skx,ky,k2,Z2∝Skx,ky,k1,Z1⋅eikz1d1

where, is the depth of first layer. By Eqs. 2-6, PSM can be chosen for the imaging of a two-layer model and be implemented for multi-layer medium imaging.

It is a long execution time to create the image through PSM method. Note that for each layer, the wave velocity of the propagating medium is constant under the homogeneous assumptions. Therefore, Omega-k can be of particular interest in each layer since FFT and IFFT approach is significantly faster. According to Eq. 3, we have,

2kl=kx2+ky2+kZl2

thus, we get,

dkldkZl=αkx,ky,kZl=kZlkx2+ky2+kZl2

By applying Eq. 8, we shift down the wave field at the top of the second layer, which equals to that the antenna is sampling at layer interface for the targets within second layer. Therefore, by substituting Eq. 8 into Eq. 7, we can use Omega-k algorithm to generate the image within each layer as follows,

Ilx,y,z=∭S∙αkx,ky,kZl∙eikxx+kyy+kZlzdkxdkydkZl

where S is the wave field at each layer interface. As we know that in frequency wave-number space, if and are defined as uniform, the regular FFT algorithm cannot be directly evaluated on as it is non-uniform. Therefore, a re-interpolation from the k to domain is essential. The imaging procedure for a two-layer medium scenario is summarized in Fig.2a. With known depth and relative permittivity of each layer, the method can be extended to multi-layer medium imaging.

The size of the input data is given by the number of frequency samples (normally also the number of depth pixels), P, the number of measurement positions, M and N. For traditional back-projection method, the complexity is ), which is far more than EWK and PSM. For PSM and EWK, the asymptotic complexities of FFT operations are the same. However, for PSM, the complexity to iterative generate the image pixel-by-pixel is , whereas for EWK, the execution complexity by FFT operation is), where L is the medium layer number. When L is limited, it is obvious that EWK tends to have a lower complexity than PSM with regard to the increased P, as . Here, we give an example of computational time of PSM and EWK, respectively, shown in Fig. 2b, where P varies from 50 to 300, with and . It is obviously that EWK has a much faster execution speed than the PSM. The execution time is clearly consistent with the aforementioned complexity analysis.

The relative permittivity of the medium is an essential value that determines the wave transmission velocity in the respective medium. Unlike down-looking radar, for a two-layer medium, strong boundary reflection, transmission loss and wave scattering make it difficult to obtain the target curve directly for wave velocity estimation. In this paper, we propose a velocity estimation method based on the vertex region extraction; the reflected signal at the vertex region is highly sensitive to the change in wave velocity due to the multiple mediums. The procedure can be summarized as three steps:

i. Establish image data matrix referring to velocity change .

Radar transmission speed in the second layer is assumed to be and is velocity interval. Typical underground medium may involve different relative permittivity. Choose a set of B-scan data and substitute into (7) in sequence. Then, we generate the image data matrix towards.

ii. Extract Target region.

In the data matrix, we define the relative difference function between the normalized average value and the maximum value of as below,

pvxi,zj=Pmaxxi,zj-maxAvexi,zj,Gx,z,v1,Gx,z,vN

The target image can be expressed as the point spread function in Sinc function. The difference between the focused image and the defocused image is the peak value. A wrong velocity will blur the image, along with the decrease of the peak amplitude value of the image. For a cluttered area after applying Eq. 12, approximates to zero, because the clutter energy value has little relations with within the medium layer. However, for the target region, the energy is the strongest with the correct wave velocity. Therefore, this valve will be significantly strengthened after Eq. 12. We define the normalized as,

Pvxi,zj=pvxi,zjmaxabspvxi,zj

Then, extract the target region as:

Cxi,zj=1Pvxi,zj≥T10Pvxi,zj<T1

The region where corresponds to the curve at vertex region in the radar image, where the main component is the target echo hyperbola with very little clutter influence.

iii. Velocity Curve.

For the cluttered region, the amplitude along the axis decreases or rises monotonically, thus the vertex point typically arrives ator for monotonic properties; while for the target region, the amplitude appears to be a parabolic curve, the vertex point normally arrives at the central curve. We find the maximum point, and plot the energy curve refers to the. The maximum value position refers to the correct wave velocity, where the amplitude of the target is best focused after image processing.

To do a realistic validation that EWK can produce high quality image with less time, an experiment with a two-layer structure was conducted as shown in Fig. 3. The ability to detect and image buried objects will depend on soil transmission and scattering. Soils with grain sizes larger than the wave length have low transmittance due to scattering. Dry quartz sand has high transmission in the millimeter wave range due to its fairly uniform grain size, which is smaller than the wavelength of mm Waves.

The diagram shown in Fig. 3 is a simplified view of the main elements and the architecture of the prototype system. The prototype consists of six components: (1) Amplified-multiplier (AMC-i) module; (2) mixer-amplified-multiplier module; (3) Vector Network Analyzer (VNA); (4) software controlling interface; (5) W band antenna; and (6) 2D scanning platform. All components are connected with phase stable cables.

The VDI module extends the frequency range of the 26.5 GHz PNA up to 75~110 GHz via AMC-i. A high direction coupler directs the transmitting signal and reference and measurement receive signals to form the S11 measurements via MixAMC-i. A single 23-dB gain horn antenna with 0.22 rad beam-width is used to transmit and receive the signal. Finally, the transmission coefficient was recorded by the VNA as the scattered data from the targets. Targets of different shapes are tested in the experiments. The ground medium is white sand. The calibrated data are stored and made available in the frequency domain in MATLAB data file format.

Before the experiment, velocity is firstly estimated. The X shape is buried in sand for 15 cm, and the matrix is generated by using EWK method. The threshold we chose is . Note that the parameter is the threshold to pick the target region out from the completely scanning scenario. In our paper, we use Gray Edge Cutting Method to find a proper to extract the strongest target region for velocity estimation^[11]. Our speed ranges from . The signal energy measured about the target area and clutter area versus is shown in Figs. 4a and 4b. The trend line for the energy measured about the clutter region increases monotonically, whereas the curve for the energy about the target is hyperbolic. The curve crests at . Therefore, the corresponding relative permittivity is .

We have conducted a couple of measurements to validate the proposed algorithm. In the first experiment, a resolution plate was buried in sand for 8 cm. The image results are shown in Fig. 5. It is evident from Fig.5a and 5b that the images produced by PSM and EWK are visually similar, and both results are clearly discernible in the images. The calculation time for PSM and EWK are 642 s and 215 s, respectively. The shape of the plate is clearly visible even for the 4 mm slot. The 2D images defocuses a little because the resolution plate is not perfectly placed horizontally. To demonstrate the alignment between target and sand surface, a depth trace through the central region of the targets is drawn in Fig.5d and the result gives indications of the target depth. The first peak is the system coupling delay; the second peak is the sand reflection, which also represents the distance from transceiver towards sand surface is 40.2 cm; and the third peak is the target position arrives at 51.76 cm, the depth of the measured target is 7.7 cm, which agrees to the buried depth. The reflections from sand surface can also be seen with significant intensity values of -7.7dB in the image.

In the second experiment, X shape is buried in sand with depth of 10 cm, and O shape is put on the surface of the sand. We covered the sand surface with fresh-cut dry grass. The imaging results are shown in Fig. 6 and both targets and their specific features are clearly visible. As expected, deeper targets have stronger transmission loss, and the loss at 10 cm is around 8 dB in comparison with the O shape. The grass has strong effect on the O shape image because of the attenuation and reflection.

Unlike SAR image in free-space, the transmission loss in the standoff GPR is proportional to the buried depth. We pick out the maximum pixel of O shape and X shape, respectively, and plot the depth traces in Fig. 6d. The focused position (maximum energy) of O shape is at 52 cm. The focused position of X shape is around 62 cm with maximum energy value of -8.9dB. However, the side-lobe of O shape at 62 cm is around -15.56 dB. This energy causes a strong ghost in Fig. 6b under the image dynamic range of -30 dB.

Afterwards, we re-buried the X shape target in sand for 15 cm and removed the grass. As shown in Fig. 7, the image quality is greatly improved without grass, the maximum energy of X shape is -11 dB due to higher transmission loss. Large bandwidth not only provides superior range resolution, also it helps to reduce the interference from the environment at a given range and to provide better azimuth resolution.

We compare our work with current state art of work^[12]. Our system is demonstrated without focusing lens, which is more easily for practical deployment. Furthermore, the depth we can detect is deeper thanthat in Ref. [12]. In this paper, we have shown imaging of targets under laboratory conditions with dry sand. Thus, a system with higher dynamic range or lower frequency could be of great interest in future investigations of GPR imaging in practical scenarios with higher transmission losses and air/soil reflection losses.

In summary, by integrating the Omega-k method into the PSM, this paper presents a new EWK method for multi-layer medium imaging along with a velocity estimation method for two-layer medium. The algorithm can achieve imagery using Omega-k method with transmission velocity changes. We have shown theoretically that the algorithm has lowered computational complexity compared to the PSM algorithm with increased number of range intervals. The algorithm is formulated in 3D and been implemented for both 2D and 3D applications on synthetic data from a W band GPR prototype. Both 2D and 3D images of the data were presented. Images formed using the algorithm provide well-focusing features of the targets, even very close scattering mechanisms can be easily distinguished. The new method offers a promising faster solution with high fidelity, which provides the potential for various applications.