Data Interpretation Technology of GPR Survey Based on Variational Mode Decomposition

Featured Application

The instantaneous frequency analysis method of the GPR was established based on VMD. The proposed analysis method was used to solve the time-frequency spectrum by taking one simulation GPR signal. The proposed analysis method was applied to detect the geological hazards of the highway surveys.

Abstract

Data interpretation is the crucial scientific component that influences the inspection accuracy of ground penetrating radar (GPR). Developing algorithms for interpreting GPR data is a research focus of increasing interest. The problem of algorithms for interpreting GPR data is unresolved. To this end, this study proposes a sophisticated algorithm for interpreting GPR data with the aim of improving the inspection resolution. The algorithm is formulated by integrating variational mode decomposition (VMD) and Hilbert–Huang transform techniques. With this method, the intrinsic mode function of the GPR data is first produced using the VMD of the data, followed by obtaining the instantaneous frequency by using the Hilbert–Huang transform to analyze the intrinsic mode functions. The instantaneous frequency data can be decomposed into three frequency attributes, including frequency division section, time-frequency section, and space frequency section, which constitute a platform to gain insight into the nature of the GPR data, such that the inspected media components can be examined. The effectiveness of the proposed method on a synthetic signal from a GPR forward model was studied, with the multi-resolution performance being tested. Inspecting the media of a highroad by analyzing the GPR data, with the abnormal characteristics being designated, validated the applicability of the proposed method.

1. Introduction

The spatial distribution of underground media may be characterized by data interpretation of the effective reflection electromagnetic waves of ground penetrating radar (GPR). The extraction of echo signals is one of the most critical technologies in GPR measurement data interpretation, which is directly related to the detection accuracy towards the underground target of interest [1].

The collected GPR signal typically possesses non-stationary and random characteristics due to the complexity of the underground medium and the interference factors during the measurement. In general, Fourier transform the power spectrum estimation are the representative theories in GPR signal processing, whereas they are suitable in nature for application in stationary signal analysis [2]. With these methods, the frequency components of the whole signal can be produced, but the joint time and frequency resolution of the signal cannot be simultaneously obtained. For non-stationary GPR signals, the change in signal characteristics with time is essential in reflecting the variation in the underground medium being investigated. To this end, the instantaneous frequency characteristic of GPR signals plays an important role in the recognition of underground media [3,4]. The instantaneous frequency of a non-stationary signal is commonly examined while using time-frequency analysis methods, typically short-time Fourier transform and wavelet transform. These methods are competent for obtaining the energy intensity of the instantaneous frequency at different frequencies and times [5,6,7]. However, the above methods are based on Fourier analysis theory, so they are limited by the disadvantages of Fourier analysis.

Conventional GPR is a broadband signal, while the instantaneous frequency analysis method is based on the Hilbert transform and it is unable to separate signals at different scales [8,9]. Since Norden E. Huang proposed the Hilbert-Huang Transform (HHT), the instantaneous spectrum of the intrinsic mode function (IMF) of GPR signals can be obtained by using the empirical mode decomposition (EMD) [10,11,12]. Due to some problems in the EMD methods, such as mode aliasing and endpoint effects, the HHT method and its use are limited to a certain extent [13,14,15,16]. A noise-assisted data analysis method, the ensemble empirical mode decomposition (EEMD) is applied in order to overcome this problem, which represents one improvement of EMD. By adding finite white noise to the investigated signal, the EEMD method can automatically eliminate the mode mixing automatically [17]. The reconstructed signal from IMFs includes residual noise and different realizations of signal plus noise may produce different number of modes with EEMD. Torres proposed the complete ensemble empirical mode decomposition (CEEMD) in order to overcome these situations [18]. Furthermore, some scholars have proposed that the EEMD and the CEEMD could both be used in GPR signal analysis [19,20,21,22]. However, these methods also have limitations, for example, the lack of mathematical theory and the degrees of freedom will reduce the algorithm’s robustness, which leave room for theoretical development and the improvement of the robustness of the decomposition [23].

Recently, the variational mode decomposition (VMD), which encompasses multiple adaptive Wiener filter groups, was discovered to have good robustness with the ability to overcome the disadvantages of mode aliasing, small end effects, and pseudo components when compared to other conventional algorithms [23,24]. The VMD algorithm has been applied to some fields, such as mechanical fault diagnosis or structural damage identification, and it has achieved significant results [25,26,27,28]. VMD has also been applied to analyze seismic data [29,30]. However, relevant reports on utilizing VMD to analyze GPR data are lacking. Therefore, investigating the application of VMD theory in GPR data processing is not only academically significant, but it will also improve the detection effect of GPR in engineering.

In this study, the VMD algorithm was used to calculate the IMF, and the Hilbert transform was used to compute the instantaneous frequency of the IMF. The instantaneous frequency of the GPR was established based on VMD. The VMD, EMD, and EEMD algorithms were used to solve the time-frequency spectrum by taking one synthetic signal and one simulation GPR signal for analysis. The effect of the analysis was that it compared these three methods against each other based on the time-frequency spectrum. Finally, the VMD method was applied to detect the geological hazards of the highway surveys and to test the effects of its application.

2. Theory and Methodology

2.1. The Instantaneous Frequency Based on VMD Theory

2.1.1. VMD Theory

VMD theory’s overall framework, which includes the construction and solution of both variational models and problems, is a variational problem. The variational model is introduced to solve the signal decomposition when the VMD method is used to obtain the IMF components. The constraints are used to search for the optimal solution of the variational model. The center frequency and bandwidth of every modal function component are iteratively updated to achieve signal decomposition. Finally, the frequency band of the signal is adaptively decomposed to obtain the narrow-band IMF of the predetermined scale [23].

In the construction of variational models, the mode function $u_k(t)$ was taken as the Hilbert transform and every analytic signal was obtained. Subsequently, the central frequency of prediction was $e^{-j{\omega }_{k}t}$ and the frequency spectrum of every mode was modulated to the corresponding frequency baseband. Finally, the $L^2$ norm of the modulation signal was solved and the signal bandwidth of every mode was evaluated. The variational problem is the solution of the K mode function $u_k\left(t\right)$, and the sum of the estimated bandwidth of every mode is minimal when the constraining condition is the sum of every mode equal to the input signal. The signal decomposed to N IMF classifications. The corresponding variational problem is defined, as follows:

(1)

Every IMF component was taken as the Hilbert transform, so that the analytic signal could be obtained.

$$\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)$$

(1)

where $\delta \left(t\right)$ is the Dirichlet function and $u_k\left(t\right)$ is K IMF.

(2)

The analytic signal of the estimated central frequency can be changed to the baseband by using frequency shift.

$$\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}$$

(2)

where ${\omega }_k$ is the central frequency.

(3)

The $L^2$ norm of the module signal was solved, and the bandwidth of every mode was estimated. The variational problem was constructed, as follows:

$$\{\begin{array}{c}\underset{\left\{u_k\right\},\left\{w_k\right\}}{\min }\left\{{\displaystyle \sum _k{\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}\Vert }^2}\right\}\\ s.t.{\displaystyle \sum _k{u}_k=f\left(t\right)}\end{array}$$

(3)

For solving the variational problem, the original signal was decomposed into K narrow mode components to obtain the optimal solution to the constrained variational problem. The augmented Lagrangian can be expressed, as follows:

$$\begin{array}{r}L\left(\left\{u_k\right\},\left\{{\omega }_k\right\},\lambda \right)=\alpha {{\displaystyle \sum _k\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}\Vert }}_2^2\\ +{\Vert f\left(t\right)-{\displaystyle \sum _k{u}_k\left(t\right)}\Vert }_2^2+\langle \lambda \left(t\right),f\left(t\right)-{\displaystyle \sum _k{u}_k\left(t\right)}\rangle \end{array}$$

(4)

where $\alpha$ is the penalty factor, $\lambda \left(t\right)$ is the Lagrange multiplier, $\Vert \cdot {\Vert }_2^2$ is the norm, and $\langle \cdot \rangle$ is the inner-product operator.

The alternative direction multiplier operator was used. ${\omega }_{{}_k}^{}$, $u_{{}_k}^{}$, and $\lambda$ were updated, so the variational problem could be solved.

2.1.2. The Instantaneous Frequency Spectrum Solution

The solution condition of instantaneous frequency needs the IMF components that can be obtained with VMD. When the input signal is decomposed into several IMFs that are based on VMD, the solution of instantaneous frequency can be implemented. The Hilbert transform is used to construct the analytical signal, and the instantaneous frequency and amplitude of each modal function are obtained. The Hilbert instantaneous frequency spectrum is constructed based on the instantaneous frequency and the amplitude of each modal function. Figure 1 shows the process of instantaneous frequency spectrum implement.

In the solution of instantaneous frequency, the original signal f (t) can be decomposed into n number of IMF while using the VMD method:

$$f(t)={\displaystyle \sum _{k=1}^n{u}_k\left(t\right)}$$

(5)

The Hilbert transform was performed on IMF. The transform equation is expressed as:

$$H\left[u_k\left(t\right)\right]=\frac{1}{\pi }P{\displaystyle {\int }_{-\infty }^{+\infty }\frac{u_k\left(\tau \right)}{t-\tau }}$$

(6)

where $P$ is the Cauchy principal value and $\tau$ is time.

The analytical signal can be constructed using Formula (6) and IMF. The formula is expressed as:

$$A\left[u_k\left(t\right)\right]=u_k\left(t\right)+iH\left[u_k\left(t\right)\right]$$

(7)

In the signal’s Equation (7), its amplitude and phase is:

$$a_k\left(t\right)=\sqrt{u_k^2\left(t\right)+H^2\left[u_k\left(t\right)\right]}$$

(8)

$${\theta }_k\left(t\right)=\arctan \frac{H\left[u_k\left(t\right)\right]}{u_k\left(t\right)}$$

(9)

The phase of the signal was derived. Then, the instantaneous frequency of the signal can be obtained:

$${\omega }_k\left(t\right)=\frac{1}{2}\frac{d{\theta }_k\left(t\right)}{dt}$$

(10)

where ${\omega }_k\left(t\right)$ is the instantaneous frequency.

Then, $f\left(t\right)$ can be expressed as the following:

$$f\left(t\right)=\mathrm{Re}\left({\displaystyle \sum _{k=1}^n{a}_k\left(t\right)e^{j{\varphi }_k\left(t\right)}}\right)$$

(11)

Formula (11) was expanded and the Hilbert instantaneous frequency is expressed as the following:

$$H\left(t,\omega \right)=\mathrm{Re}\left({\displaystyle \sum _{k=1}^n{a}_k\left(t\right)e^{j{\displaystyle \int {\omega }_k\left(t\right)dt}}}\right)$$

(12)

where $H\left(t,\omega \right)$ is the Hilbert instantaneous frequency spectrum and $\omega$ is the frequency.

Based on VMD, the Hilbert instantaneous frequency spectrum can be solved from Formula. $H\left(t,\omega \right)$ describes the amplitude of the signal in the entire frequency domain with the time and frequency of the transformation. In Fourier analysis, the energy representation at a certain frequency is a harmonic wave that exists for the entire time. By contrast, in Hilbert’s instantaneous frequency analysis, the frequency energy at some point may indicate its occurrence in the entire time. A frequency wave appears and it represents the probability of occurrence.

2.2. Time-Frequency Analysis Method of GPR

The radargram is the synthesis of the time domain signal. The instantaneous frequency spectrum method cannot analyze the radargram very well. The instantaneous frequency analysis method of GPR needs extending from single time domain signal. During the computation process of the GPR time-frequency analysis, every trace data was decomposed by VMD into the IMF. The analytic signal was constructed using the Hilbert transformation. Subsequently, the instantaneous frequency and amplitude was obtained. Finally, we obtained the Hilbert spectrum. Recording the data of two-dimensional GPR along a measuring line can be expressed as $f\left(x,t\right)$. The reflected pulse spectrum of GPR by Hilbert transform can be expressed as:

$$H\left(x,t,\omega \right)=\mathrm{Re}\left({\displaystyle \sum _{k=1}^n{a}_k\left(x,t\right)e^{j{\displaystyle \int {\omega }_k\left(x,t\right)dt}}}\right)$$

(13)

The GPR input data of the time-frequency analysis is profile data. According to Formula (13), the GPR data of the time-frequency analysis is three-dimensional (3D) data. Thus, 3D GPR data can be divided into three frequency attributes. Figure 2 shows the analysis process.

When 3D GPR data frequency ω is fixed, $H\left(x,t\right)$ is the division frequency section:

$$H\left(x,t\right)=H\left(x,t,\omega \right)|\omega =const$$

(14)

The division frequency section is a filter cross-sectional view, as different frequency bands for the filtered scan and the reflected signal can be investigated in each frequency band. The division frequency section is the most common method for frequency analysis of GPR.

When the 3D GPR data location x is fixed, $H\left(t,\omega \right)$ is the time frequency section:

$$H\left(t,\omega \right)=H\left(x,t,\omega \right)|x=const$$

(15)

The time frequency analysis section records the frequency characteristics and the reflected wave distributed in the frequency axis.

When the 3D GPR data time t is fixed, $H\left(x,\omega \right)$ is the space frequency section:

$$H\left(x,\omega \right)=H\left(x,t,\omega \right)|t=const$$

(16)

The space frequency section shows the spatial frequency analysis of GPR data. The frequency components correspond to the changes in the reflection signal of certain depth along the survey profile.

In the three frequency attributes, the transverse profile is the division frequency section, the longitudinal profile is the time frequency section, and the horizontal profile is the space frequency section. The time-frequency analysis method of GPR is based on the three frequency attributes.

3. Comparison of VMD, EMD, and EEMD

A synthetic signal s is taken as an example in order to test the performance of VMD decomposition for signal. The synthetic signal s contains the two high frequency components under the background of the low frequency component. The high frequency components need to be separated from background for testing the effect of the VMD decomposition. The synthetic signal s is overlapped by s₁ and s₂, where s₁ is the sinusoidal periodic function and s₂ is the segmented sine function (Figure 3a). The function is as follows:

$$s_1=\sin \left(20\pi t\right)$$

(17)

$$s_2=\{\begin{array}{c}0.2\times \sin \left(300\pi t\right)when0.05\le t\le 0.1\\ 0.4\times \sin \left(100\pi t\right)when0.15\le t\le 0.25\\ 0other\end{array}$$

(18)

$$s=s_1+s_2$$

(19)

The synthetic signal s was decomposed while using VMD (Figure 3a). Figure 4b–d shows the results of the decomposition.

According to the results of Figure 3b–d, the signal s was decomposed to obtain three IMFs. The three IMFs correspond to the three sine functions of signal s. IMF1 corresponds to s₁, IMF2 corresponds to s₂ in 0.15~0.25 s, and IMF3 corresponds to s₂ in 0.05~0.1 s. The results show that the VMD can properly separate IMFs from the signals.

In order to study the performance advantages of VMD, the signal s was decomposed using VMD, EMD, and EEMD, respectively, and the instantaneous spectrum of the signal s was obtained by using these three methods (Figure 4).

The results for the instantaneous frequency of s in Figure 4 show the performances of these three methods. In Figure 4a, the EMD method only obtained a component at a frequency of 10 Hz, but failed to get the components at the 50 Hz and 150 Hz frequencies, most especially the components at the 150 Hz frequency. In Figure 4b, EEMD can separate the different components into the three frequencies. However, the component at 10 Hz frequency shows some fluctuations in the 0.15~0.25 s range, and the other two components have some bending corresponding to the frequency section. We found that EEMD can separate the different components, although the accuracy is not high. According to Figure 4c, the VMD method clearly separates the three different components with accurate frequencies for each time period. Thus, VMD’s effectiveness in the elimination of mode aliasing is significantly superior to that of the other two methods.

During the mode decomposition process, the endpoint effect is usually a serious problem in EMD, since it causes the energy to be changed after decomposition. The signals u(i) decomposition generates and the IMF energy can be expressed, as follows:

$$E=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{u}^2\left(i\right)}}{n}}$$

(20)

where E is the IMF energy of n samples, u(i) is the signal, and n is the sample number of signals.

The energy leakage can be reflected by defining the deviation between the total energy of the IMF and the original signal energy. The evaluation index can be expressed, as follows:

$$\xi =\frac{\left|\sqrt{{\displaystyle \sum _{i=1}^k{E}_i^2}}-E_u\right|}{E_u}$$

(21)

where E_u is the signal energy, E_i is the IMF energy of i number, and k is the total number of IMFs.

According to

Table 1. Evaluation index of energy leakage.

	ξ
EMD	EEMD	VMD
0.065	0.0125	0.0119

, the energy leakage when using VMD is smaller than that when using EMD or EEMD, while the end effect of using VMD is not obvious. The VMD is composed of a plurality of adaptive Wiener filters. It avoids the error accumulation of EMD and EEMD that is caused by the uncertainty of the extreme values at the endpoints during recursive components. Therefore, it is superior to the other two methods regarding the endpoint effect.

4. Time Frequency Analysis of GPR Simulation Signal

4.1. Time Frequency Analysis of VMD

A simulation model was created in order to evaluate the contribution of the time frequency analysis based on Matlab 2016 platform (The MathWorks, Natick, MA, USA). During the simulation signal acquisition of GPR, the size of the assumed structural model is shown in Figure 5a. There was one steel cylinder, where the center coordinates were (2.5, 1.5) and the cylinder radius was 0.2 m. For the simulation parameters, the diectric permittivity and electrical conductivity of the concrete were 10 and 0.001 S/m, while the diectric permittivity and electrical conductivity of the steel were 1 and 106 S/m. The computation area was 5.0 m × 3.0 m. The grid size was 0.002 m and the length of the time window was 40 ns. The Ricker wave was considered as input excitation and the central frequency was 0.5 GHz. The transceiver distance was 0.04 m and the number of the simulation sampling trace was 230. Along the direction of the x axis, the 230 time-series GPR signal was synthesized and Figure 5b shows the scanner result of the simulation with the finite difference time domain (FDTD).

In Figure 5b, the reflected signal of the direct wave is significantly shown and the reinforcement reflex arc is also displayed in the center of the concrete. However, the reflected arc was relatively weak and it was not clearly distinguished in the internal structure. Based on the data of the GPR simulation scan, the three-dimensional body can be constituted after the time frequency analysis. Figure 6 shows the results.

Figure 6 shows the 3D data of GPR $H\left(x,t,\omega \right)$. The data structure consists of two parts, in accordance with the direct wave and reflection wave signal of reinforcement, and the concrete data effectively displays internal structure information. However, the internal details of the structure are not well displayed in the 3D GPR data and the body of the 3D GPR data needs slicing for analysis. From what is mentioned above, it can be concluded that the transverse profile is the division frequency section, the longitudinal profile is the time frequency section, and the horizontal profile is the space frequency section. The sliced 3D data body of GPR was processed to study the internal structure detection effects with the time frequency analysis of GPR.

First, the data volume $H\left(x,t,\omega \right)$ was transversely cut to obtain a frequency division section. Figure 7 shows the frequency division section $H\left(x,t\right)$, with a central frequency of 0.5 GHz. In Figure 7, the main frequency contour displays the reflex arcs of the hyperbola below the direct wave and the reflected signal of the reinforcement as compared with Figure 5b, which has been significantly enhanced. The effect of the abnormality detection significantly improved.

Secondly, the GPR data $H\left(x,t,\omega \right)$ processed longitudinal cutting and the time frequency section $H\left(t,\omega \right)$ could be obtained. Figure 8 displays the time frequency section at the center of the GPR line. In Figure 8, the two extreme points are 4 ns and 20 ns. The two time positions correspond to the travel times of the direct wave and reinforcement reflection wave, with both frequencies in the 0.5 GHz GPR transmitter main frequency.

Finally, the GPR data body $H\left(x,t,\omega \right)$ processed the horizontal slice. The space frequency section $H\left(x,\omega \right)$ could be obtained. In Figure 9, the space frequency section shows the layer at time 20 ns. In Figure 9, the extreme range in the x-axis is about 2.1 m~2.5 m, in relation to the bar diameter. The frequency of the position was 0.5 GHz, which was mainly transmitted at the frequency of GPR.

According to the above analysis, it can be found that the proposed method is a very useful tool for the interpretation of GPR data. The method not only provides information of the detection targets in different views, but it also displays internal information of the detection target in detail. In these three kinds of sections, the frequency division section can be used as a comparison with the original image, which is the most common method of frequency analysis in GPR, while the other two sections have secondary analytical verification and other view functions.

4.2. Time Frequency of VMD Comparison with Other Methods

For further verifying the advantage of the proposed method, VMD comparison with other methods were studied for signal interpretation of GPR and they used a detection of three-layer concrete pavement as an example. See Figure 10a for the GPR forward sketch. The model measured 2.5 m × 0.45 m; the uppermost layer was air, the middle layer was concrete, and the bottom layer was dry sand. There was one incline facet between the concrete and dry sand. The depth of the concrete was between 0.15 m and 0.25 m. The permittivity and conductivity of the concrete were 6 and 0.005 S/m, respectively, while the permittivity and conductivity of the dry sand were 3 and 0.0001 S/m respectively. During the simulation of GPR with FDTD, the computational grid size was 0.0025 m, the length of the time window was 12 ns, and the absorbing boundary was an unsplit field perfectly matched layer (UPML). The central frequency of the transceiver antenna was set to 900 MHz, the transceiver antenna distance was 0.025 m, and the work mode was self-excitation and self-receiving. From the left 0.1 m position to the right 2.3 m at 0.02 m intervals, a total of 115 trace data was obtained. Figure 10b shows the scanner result of the simulation.

In Figure 10b, each IMF can be decomposed using VMD and the frequency division section can be obtained while using the Hilbert transform. Figure 11 shows the frequency division section of the simulation signal of GPR at 1.0 and 1.2 GHz frequencies.

According to Figure 11, the frequency division section at 1.0 GHz clearly shows the ground-coupled direct wave (Figure 11a), and the slice of the time-frequency spectrum at 1.2 GHz clearly shows the reflected wave between the concrete and the dry sand interface (Figure 11b). It can be seen that the time-frequency spectrum separates the different GPR components that are based on VMD.

EMD and EEMD were used to calculate the time-frequency spectrum of the GPR signal simulation in order to further study the effects of VMD on GPR signals, respectively. Figure 12 shows the time-frequency spectrum at 1.0 GHz and 1.2 GHz based on EMD, and Figure 13 shows the time-frequency spectrum at 1.0 GHz and 1.2 GHz based on EEMD.

In Figure 12, the frequency division section at 1.0 GHz has a weak performance on the direct wave (Figure 12a). However, the frequency division section at 1.2 GHz does not show a reflected wave on the interface between the concrete and sand (Figure 12b), but there is still direct wave information present. The time-frequency spectrum that was obtained using EMD is not clear on the separation between the different components of the radar signal.

Figure 13a shows the frequency division sections at 1.0 GHz based on EEMD. The frequency division section at 1.0 GHz shows direct waves. However, the time-frequency spectrum at 1.2 GHz displays reflection at the interface of concrete and sand (Figure 13b), but the effect is not obvious. There are direct waves and other information in Figure 13b. EEMD cannot decompose each of the frequency components of the simulation GPR signal. VMD is better than the other two methods in separating the different components of the GPR signal among the three methods for solving the time-frequency spectrum analysis.

5. The Analysis of Real Exploration GPR Signals

The Yunnan Province of China plans to build a highway through a mining area and the project required a geophysical survey. The detection device used was a GSSI SIR-3000 GPR. The central frequency of the antenna was 100 MHz. The length of the detection section was 15.0 m and acquisition space interval was 0.050 m. The cross section had 300 trace records, with each trace containing 512 sampling points. Figure 14a shows the survey records of the GPR.

Based on the original records of the field (Figure 14a), the mining area was not clear enough. The detection data was processed using VMD Time Frequency in order to test VMD’s effect. During processing GPR data, there are several processing steps that are applied. Firstly, the preprocessing of original GPR data is performed to improve the observation of embedded cavities in the radargrams. Mean denoising is used for filtering random noise, and removing DC is used for decreasing the direct current effect. Secondly, every trace time domain data was decomposed by VMD into the IMFs after the preprocessing. Subsequently, the analytic signal was constructed using Hilbert transformation that is based on IMFs, and the instantaneous frequency and amplitude was obtained at every trace. Finally, the Hilbert spectrum of every trace was integrated, and the time-frequency spectrum of the detection signal was obtained while using the proposed method. The 3D GPR data can be sliced on frequency division plane at 80 kHz. Figure 14b shows the frequency division section of the signal at 80 kHz.

In Figure 14b, it can be seen that there were obvious anomalies that were located in the survey line at 4.5 m, which are obviously different from other regions. The region with anomalies is displayed within elliptic area in Figure 14b. It is a verified fact that there were caverns at the 4.5 m location of the survey line and 7.0 m in depth after drilling. In Figure 14b, the quality of the data section is significantly improved when compared with Figure 14a. The non-uniformity recognition of the medium with the time frequency section is more distinct than the original time domain section. The local anomaly features are highlighted, and the goaf is clearly visible from the time frequency section. In Figure 15, frequency division sections at 80 kHz were also obtained with EMD and EEMD. However the region with anomalies does not display well when comparing with Figure 14b. VMD is also better than the other two methods in application. The evidence proves that the proposed method has significant applicability for processing GPR signals.

6. Conclusions

Based on the VMD and the Hilbert transform theory, the VMD technique was introduced into the GPR data interpretation process, and the instantaneous frequency analysis method of the GPR signal was also established. The analysis results showed that the proposed method has better outcomes than the traditional methods. Finally, the method was applied to a real engineering exploration. The conclusions are as follows:

(1)

The VMD method is better at overcoming the influences of mode aliasing and endpoint effects as compared to the traditional methods EMD and EEMD.

(2)

The time-frequency spectrum of the GPR signal can be obtained using VMD. The slice section of the 3D time-frequency data body that was constructed from the GPR scan line data displays the information of the detection target from multiple angles.

(3)

VMD separates the IMF of the GPR signal and has obvious advantages when compared to the traditional method EMD and EEMD.

(4)

According to the time-frequency spectrum of the GPR signal based on VMD, the processing method can highlight the abnormal areas and it shows good applicability in the discovery of the cavern anomalies.