## 2 Materials and Methods
### Ensemble Empirical Mode Decomposition  
In 1998, [PERSON] et al. [21] proposed EMD for Hilbert transform with physical meanings to analyze nonlinear and non-stationary signals in multiple resolutions. It offers a self-adaptive approach to decompose any complicated signal series into a finite number of IMFs, which focus on the local time scale and involve only one oscillation mode in each waveform. The instantaneous frequency defined in an IMF corresponds well with the local physics of the studied system, which would have an expected performance in interpreting short-distance features from GPR radargrams. According to the definition, an arbitrary signal \(X(Z)\) and its Hilbert transform \(Y(Z)\) form the complex conjugate pair to obtain the analytical signal \(P(Z)\). \[P(Z)=X(Z)+iY(Z)=a(Z)e^{i\theta(Z)} \tag{4}\]  
where, \(i\) represents the imaginary unit, \(a(Z)=\sqrt{X^{2}(Z)+Y^{2}(Z)}\) is the signal amplitude, and \(\theta(Z)=arctan(\frac{Y(Z)}{X(Z)})\) is the signal phase. Therefore, the instantaneous frequency is defined as:  
\[\omega=\frac{d\theta(Z)}{dZ} \tag{5}\]  
In order to handle the issue that the instantaneous frequency has no precise calculation approach, signal data are limited by narrow bands to provide physical meanings for the instantaneous frequency, and these narrow bands can be obtained by EMD procedure. Processing a single signal series, the EMD method applies the following procedures:  
1. Firstly record the initial data, \(r_{1}=A(Z)\), and initialize from \(i=1\), \(j=1\);
2. Find the upper and lower envelopes (\(U_{j},L_{j}\)) of the data by cubic spline interpolation and then calculate their means, \(m_{j}=(U_{j}+L_{j})/2\);
3. Setting \(h_{i,1}=r_{i}\), obtain the rest part by removing the means \(m_{j}\), \(h_{i,j+1}=h_{i,j}-m_{j}\), \(j=j+1\);
Repeat the procedure ii. and iii. until \(h_{i,j+1}\) satisfies the decomposition threshold, and we can obtain the \(ith\) IMF as \(C_{i}=h_{i,j+1}\). In [PERSON]'s method, the threshold parameter \(SD\) should be between 0.2 and 0.3;  
\[SD=\sum_{k=1}^{j}\frac{m_{k}^{2}}{h_{i,j}^{2}} \tag{6}\]  
\(i=i+1\), calculate the rest part of signals \(r_{i}=r_{i-1}-C_{i-1}\). Repeat the above steps ii. to iv. to obtain diverse IMFs and the final residue data \(r_{n}\). \[A(Z)=\sum_{i=1}^{n-1}C_{i}+r_{n} \tag{7}\]  
After the EMD algorithm, the analytical signal for Hilbert transform with physical meanings is expressed as:  
\[P(Z)=\sum_{i=1}^{n-1}C_{i}e^{i\int\omega_{i}(Z)dZ} \tag{8}\]  
However, the significant drawback of EMD, mode mixing appearing frequently, has depressed its performance and restricted the applications. This effect either mingles signals of widely disparate scales in an IMF, or residents similar-scale signals in different IMFs. A noise assisted method, EEMD, was then developed by [PERSON] and [PERSON] [24] to overcome the problem of mode mixing. In the EEMD approach, the true IMFs are defined as the ensemble mean of multiple decomposition trials, each with initialized data of the original signals plus a random white noise with narrow amplitude, and this thereby improves mode separating as other researches demonstrate [27; 28]. The EEMD algorithm can be summarized as follows:  
Initializing from \(i=1\), generate a random white noise and add it to pre-processed signals, \(A_{i}(Z)=A(Z)+noise_{i}\). The standard deviation of noise amplitudes is suggested to be \(1/5\) that of the signals;  
Apply EMD to decompose \(A_{i}(Z)\) into different IMFs \(C_{i,1},C_{i,2},...,C_{i,n-1}\);  
\(i=i+1\), repeat the above two steps for \(N\) times, with the minimum value of \(N\) should satisfy Equation (9), where \(\varepsilon_{n}\) is the final standard deviation of error arising from added noises, and \(\varepsilon\) is the standard deviation of noise amplitude;  
\[ln\varepsilon_{n}+\frac{\varepsilon}{2}lnN=0 \tag{9}\]  
Obtain the true decomposition results by calculating the ensemble means of all corresponding IMFs.