## 1 Introduction  
Ground Penetrating Radar (GPR) is a fast, high-resolution, nondestructive detection technology which has been widely used in various fields such as building structure investigation [1; 2], environmental research [3], engineering detection [4], etc. Due to the complex and varied distribution of the underground medium, multiple waves, including diffraction waves, are usually produced when electromagnetic waves encounter a contrast of electrical properties. In most cases, the GPR measurement area is filled with various electromagnetic interferences; the noise brought on by the surrounding environment will contaminate the collected GPR observed signal, such that the GPR profile is unable to accurately reflect the underground structure, thus requiring methods to improve the imaging solution [5]. Migration typically provides an efficient method for imaging complicated structures. It can properly position the reflection waves and focus the diffraction waves backward to their proper spatial position [6; 7], and the distribution of underground structures can be obtained effectively [8] using finite-difference migration [9], Kirchhoff integration migration [10],F-K migration [11], reverse time migration (RTM) [12; 13; 14], etc. The RTM method, based on two-way wave equations, properly utilizes full-wavefield information and possesses the unique advantages of high precision, accurate phase, and propagation in any specific direction [15; 16], which is suitable for imaging complex structures [17]. However, the RTM method remains only the conjugate transpose of the forward operator, not its inverse, and the result of the RTM method does not represent the reflection coefficient. For example, RTM images often contain various undesirable artifacts when sensor bandwidth and coverage are limited, which can affect the imaging resolution of underground structures including possible anomalies. Furthermore, conventional RTM results usually contain only information on the location of underground structures without providing sufficiently precise amplitude information, which makes it impossible to obtain an accurate velocity distribution of the subsurface medium [18]. Least-squares reverse-time migration (LSRTM) solves RTM imaging as a least-squares inversion problem and achieves an imaging result that best matches GPR observation data by implementing an iterative algorithm [19]. Based on the observational data space, the method is generally employed on ray class and wave equation continuation operators to perform high-resolution subsurface imaging. However, in cases of severely degraded data quality, LSRTM imaging results contain artifacts that usually cause inversion instability and slow convergence of the inversion [20]. Since LSRTM involves solving the linear inverse problem, several different approaches have been developed to improve LSRTM imaging. For instance, quadratic and non-quadratic regularization (Cauchy norm) have been incorporated into LSRTM to improve resolution while suppressing spatial artifacts [21]. Qp-LSRTM has been employed to rectify distortion caused by strong subsurface attenuation. The merit of this approach compared with conventional RTM and LSRTM is that Qp-LSRTM compensates for the amplitude loss due to attenuation and produces images with better balanced amplitudes and more resolution below highly attenuative layers [22]. Furthermore, Tikhonov regularization is employed alongside LSRTM for ensuring the stability of the inversion. However, Tikhonov regularization usually tends to over-smooth the boundary between the target region and the background region [23]. Instead, approximate Total Variation (TV) regularization based on image reconstruction helps to recover the discontinuous image efficiently and has many other advantages. For example: limiting the smoothness of the image, eliminating the stair effect, avoiding edge blurring, etc. [24].