## 2 Methodology
### BP Imaging Algorithm  
Previous studies on GPR imaging algorithms mainly focused on the back-projection (BP) algorithm, the Kirchhoff Migration algorithm, the F-K Migration algorithm, the Diffraction Tomography (DT) algorithm, the Range Migration (RM) algorithm, and the Reverse Time Migration (RTM) algorithm, and others [42; 43]. Among them, the BP algorithm fully considers the effect of media stratification in the imaging geometry scene. The refraction of electromagnetic waves at the dielectric interface can be compensated for more precisely and adapted to the road structure [43], so it has great application value in road detection [44].  
The BP algorithm is derived from Computer Tomography (CT) imaging technology [45], which was initially one of the imaging methods of Synthetic Aperture Radar (SAR). Since the working principle of GPR is the same as that of SAR, the imaging method of GPR can learn from that of SAR [46]. The principle of the BP algorithm is the coherent superposition of the target scattered energy in the time domain of the received signal to obtain high-resolution imaging results [47; 48; 49].  
The principle of BP algorithm imaging is shown in Figure 5. The upper layer is air with a permeability of 1 and a relative permittivity of 1, while the lower layer is a uniform dielectric with a permeability of \(\mu_{1}\) and a relative dielectric constant of \(\varepsilon_{1}\). Radar antennas send and receive signals at \(N\) positions. For the nth antenna, its coordinate is set as \((x_{n},-h)\), and its electromagnetic wave is refracted at \(x_{r}\). A random point \(A(x,z)\) in the medium is an example for analysis. According to [PERSON]'s law and geometric relations, Equations (6)-(8) can be obtained.  
Figure 4: GPR image of the underground cavity model.  
Figure 3: Three-level wavelet decomposition.  
\[\frac{\sin\theta_{0}}{\sin\theta_{1}}=\sqrt{\varepsilon_{1}} \tag{6}\] \[\sin\theta_{0}=\frac{x_{r}-x_{n}}{\sqrt{h^{2}+\left(x_{r}-x_{n} \right)^{2}}}\] (7) \[\sin\theta_{1}=\frac{x-x_{r}}{\sqrt{z^{2}+\left(x-x_{r}\right)^{2}}} \tag{8}\]  
where \(\theta_{0}\) and \(\theta_{1}\) are the incidence angle and refraction angle, respectively. Equation (9) can be obtained by combining Equations (6)-(8):  
\[\left(x_{r}-x_{n}\right)^{2}[z^{2}+\left(x-x_{r}\right)^{2}]=\varepsilon_{1} \left(x-x_{r}\right)^{2}[h^{2}+\left(x_{r}-x_{n}\right)^{2}] \tag{9}\]  
According to this equation, the position of the refraction point \(x_{r}\) can be calculated, and then the round-trip propagation time \(T_{A,\;n}\) of the electromagnetic wave emitted by the antenna at point \(A\) can be obtained, as shown in Equation (10):  
\[T_{A,n}=\frac{2}{c}(\sqrt{h^{2}+\left(x_{r}-x_{n}\right)^{2}}+\sqrt{\varepsilon _{1}}\sqrt{z^{2}+\left(x-x_{r}\right)^{2}}) \tag{10}\]  
where \(c\) is the speed of light.  
The amplitude of the nth antenna at point \(A\) can be obtained by combining \(T_{A,\;n}\) and the GPR detection data, which is shown in Equation (11):  
\[U_{A,n}=\left.I_{n}(t)\right|_{t=T_{A,n}} \tag{11}\]  
where \(I_{n}(t)\) indicates the A-Scan data at the nth antenna position at time \(t\).  
The final response amplitude of imaging point \(A\) can be obtained by summing the amplitude of point \(A\) with the antenna at each position on the measuring line. The above steps are repeated for each point throughout the imaging area to attain the final image result. Taking Figure 4b as an example, the effect of BP imaging is shown in Figure 6.  
Figure 5: Principle of BP algorithm imaging.