## 4 Numerical Experiment  
### Effectiveness Verification of MWT-MSSA Method  
To verify the efficiency of the MWT-MSSA GPR diffraction separation method proposed in this work, we first perform testing and analysis on simulated GPR data that feature small-scale anomalies. The model dimensions were set to 8.0 m \(\times\) 2.5 m, with two undulating reflection interfaces located at approximately depths of 0.9 m and 1.4 m, respectively. The designed relative dielectric permittivity and electric conductivity of the layers are shown in Figure 3a. Three voids with diameters of roughly 0.1 m are specifically designed to be located in the second, third layers and at the interfaces. The simulation was implemented using the finite difference time-domain method (FDTD) with a Ricker wavelet source at a central frequency of 400 MHz and a time window of 50 ns. The common- offset configuration was employed for data acquisition, with the grid size at 0.005 m and the trace interval at 0.05 m. Figure 3b shows two nearly linear reflection events that correspond to the different media layers and three hyperbolic diffraction events that are produced by  
Figure 2: Flowchart of the proposed MWT-MSSA algorithm of GPR diffraction separation. the voids. It is noteworthy that the hyperbolic diffraction exhibits a unique shape when the voids are further away from the interfaces, as shown in the left and middle diffraction. However, when the voids are closer to or coincide with the interfaces, the energy at the vertex and wings of the hyperbolic diffraction is severely obscured by reflection, which maintains a stronger energy due to the presence of interfaces. This can evidently be shown in the right-side diffraction, where the reflection caused by interference could complicate the accurate identification and extraction of the diffraction. The simulated data shown in Figure 3b are first processed using MWT to separate the reflection and diffraction components into distinct wavelet transform coefficients. To overcome the ambiguities of MWT coefficients in describing the spatial and temporal distribution, we employ an inverse wavelet transform to reconstruct GPR profiles that could precisely depict the characteristics of the wavelet transform coefficients across all levels. Figure 4 presents the coefficient profiles corresponding to the first, second, and third levels after the wavelet transformation was performed on the data shown in Figure 3b. The biorthogonal wavelet basis function bior4.4 is employed herein for wavelet transform. A detailed comparison of results using other alternative wavelet basis functions will be discussed in Section 4.2. As illustrated in Figure 4a,e,i, the majority of the reflective components in the synthetic data tend to concentrate toward the low-frequency approximation coefficients \(d_{LLk}\) at each level, while the diffraction data are predominantly distributed in the high-frequency vertical detailed coefficients \(d_{HLk}\) at each level, as shown in Figure 4c,g,k. After the first-level wavelet transform, the faint energy at the wings of the diffraction is decomposed into the first-level high-frequency vertical detailed coefficient \(d_{HL1}\), whereas the remaining energy of both reflection and diffraction are retained in the first level of low-frequency approximation coefficient \(d_{LL1}\). Following the second-level wavelet transform, the majority of the diffractive energy previously contained in the first level of low-frequency approximation coefficient is now redistributed into the second level of high-frequency vertical detailed coefficient, as shown in Figure 4g, leaving only faint energy at the diffraction vertexes within the second level of low-frequency approximation coefficient, as evidenced in Figure 4e,f,h. Similarly, as presented in Figure 4i-l, the third level wavelet transform further decomposes the diffractive vertex energy (i.e., initially present in the second level low-frequency approximation coefficient) into the third level of high-frequency vertical detailed coefficient.