## 3 Methodology  
### Principle of Multi-Level Wavelet Transform  
Two-dimensional (2D) time-domain GPR data are typically represented as \(d(x,\,t)\), \(x=1,\,2,\ldots,\,N_{x}\), \(t=1,\,2,\,\ldots,\,N_{t}\), where \(N_{x}\) represents the total number of traces, and \(N_{t}\) represents the number of samples (e.g., when collecting 10 m of data with an offset of 0.02 m and sampling points of 1024, the final acquired data will have \(N_{x}\) traces range from 1 to 500 and \(N_{t}\) samples range from 1 to 1024). According to wavelet transform theory, the first-level 2D discrete wavelet transform can be expressed as [34; 35]:  
\[\left\{\begin{array}{l}LL_{1}(x,\,t)=\sum\limits_{m}\sum\limits_{n}h[m]h[n]d( 2x-m,2t-n)\\ LH_{1}(x,\,t)=\sum\limits_{m}\sum\limits_{n}h[m]g[n]d(2x-m,2t-n)\\ HL_{1}(x,\,t)=\sum\limits_{m}\sum\limits_{n}g[m]h[n]d(2x-m,2t-n)\\ HH_{1}(x,\,t)=\sum\limits_{m}\sum\limits_{n}g[m]g[n]d(2x-m,2t-n)\\ \end{array}\right. \tag{1}\]  
where \(LL_{1}\), \(LH_{1}\), \(HL_{1}\), and \(HH_{1}\) represent the low-frequency approximation coefficient, high-frequency horizontal, vertical and diagonal detail coefficient of the first-level wavelet transform of \(d(x,\,t)\), respectively. \(h[m]\) and \(g[n]\) represent the filter coefficients used to extract the low-frequency components and the high-frequency components, respectively. \(d(2x-m,2t-n)\) denotes the local sampling value. \(m\) and \(n\) are the filter indices in horizontal and vertical directions, respectively.  
The MWT builds upon the results of each previous wavelet transformation level. The low-frequency approximation coefficients \(LL_{i-1}\) derived from the previous level are further decomposed using wavelet transform. This iterative process is repeated until the desired number of decomposition levels is achieved [36].  
MWT generally involves applying wavelet decomposition iteratively to the low-frequency approximation coefficients \(LL_{i}\) obtained from the previous transformation level,  
Figure 1: Simulated GPR profile (**a**) and its singular spectrum at the dominant frequency (**b**), with the diffraction profiles separated by the mean value method (**c**) and the cumulative energy approach (**d**), respectively.  
until the specified number of decomposition levels are achieved. Based on Equation (1), the \(i\)-th level wavelet transform of \(d(x,t)\) can be expressed as [37; 38]:  
\[\left\{\begin{array}{l}\begin{array}{l}\text{{LL}}_{i}(x,t)=\sum\limits_{m} \sum\limits_{n}h[m]h[n]LL_{i-1}(2x-m,2t-n)\\ \text{{LH}}_{i}(x,t)=\sum\limits_{m}\sum\limits_{n}h[m]g[n]LL_{i-1}(2x-m,2t-n) \\ \text{{HL}}_{i}(x,t)=\sum\limits_{m}\sum\limits_{n}g[m]h[n]LL_{i-1}(2x-m,2t-n) \\ \text{{HH}}_{i}(x,t)=\sum\limits_{m}\sum\limits_{n}g[m]g[n]LL_{i-1}(2x-m,2t-n) \end{array}\right.\end{array}\right. \tag{2}\]  
Upon completing \(k\)-level wavelet transforms, \(k\) sets low-frequency approximation coefficients, high-frequency horizontal, vertical and diagonal detailed coefficients can be obtained. To visually illustrate the distribution characteristics of the GPR wavefield at each level of the wavelet transform, the inverse transform of the wavelet coefficients is applied, as described in Equation (2). This inverse transform generates time-domain GPR coefficient profiles \(d_{LLk}(x,t)\), \(d_{{LH}}(x,t)\), \(d_{{HL}}(x,t)\), and \(d_{{HH}}(x,t)\), which provide a detailed representation of the wavelet transform coefficients across each decomposition level.  
\[\left\{\begin{array}{l}d_{{LL}_{k}}(x,t)=\sum\limits_{m}\sum\limits_{n}h[m]h[ n]LL_{k}(2x-m,2t-n)\\ d_{{LH}_{k}}(x,t)=\sum\limits_{m}\sum\limits_{n}h[m]g[n]LH_{k}(2x-m,2t-n)\\ d_{{HL}_{k}}(x,t)=\sum\limits_{m}\sum\limits_{n}g[m]h[n]HL_{k}(2x-m,2t-n)\\ d_{{HH}_{k}}(x,t)=\sum\limits_{m}\sum\limits_{n}g[m]g[n]HH_{k}(2x-m,2t-n)\end{array}\right. \tag{3}\]  
Without loss of generality, the subsequent discussion refers to \(d_{{\ast}{\ast}_{k}}(x,\ t)\) as the time-domain GPR coefficient profiles that are derived from the inverse transform of the wavelet coefficients at each decomposition level.