Reflectivity-Consistent Sparse Blind Deconvolution for Denoising and Calibration of Multichannel GPR Volume Images

Vehicle-mounted multichannel ground penetrating radar (MC-GPR) is a revolutionary technology that facilitates the acquisition of volume images by arranging multiple antennas; however, its images are highly affected by noise due to different antenna characteristics. This study proposes reflectivity-consistent sparse blind deconvolution (RC-SBD) for appropriate denoising of ground penetrating radar (GPR) volume images. RC-SBD interprets the observed waveform as the convolution of the emitted wavelets and reflectivity, plus stationary clutter such as reflections from the vehicle itself. The method obtains denoised reflectivity by estimating the wavelets and clutter. The key feature of RC-SBD is that it extends the existing SBD method to 3-D, and introduces an assumption of reflectivity smoothness in the horizontal direction, expressed by the total variation (TV) regularization term. The estimation is formulated as a minimization problem involving $\ell _{2}$ and $\ell _{1}$ norms and is optimized using the Split–Bregman algorithm. Trade-off hyperparameters of the objective function are optimized via Bayesian optimization, maximizing the kurtosis of the calibrated volume image. Validation with synthetic data demonstrates accurate wavelet estimation and significant denoising of the volume image. Real-world data application further reveals considerable improvements in the channel-depth cross section, providing a clear visualization of structures like rebar and steel plates. Notably, the calibrated image remains stable across diverse datasets, including earthwork and bridge sections, showcasing the versatility and reliability of the proposed methodology.


Reflectivity-Consistent Sparse Blind Deconvolution for Denoising and Calibration of Multichannel GPR Volume Images
Takanori Imai and Tsukasa Mizutani , Member, IEEE Abstract-Vehicle-mounted multichannel ground penetrating radar (MC-GPR) is a revolutionary technology that facilitates the acquisition of volume images by arranging multiple antennas; however, its images are highly affected by noise due to different antenna characteristics.This study proposes reflectivity-consistent sparse blind deconvolution (RC-SBD) for appropriate denoising of ground penetrating radar (GPR) volume images.RC-SBD interprets the observed waveform as the convolution of the emitted wavelets and reflectivity, plus stationary clutter such as reflections from the vehicle itself.The method obtains denoised reflectivity by estimating the wavelets and clutter.The key feature of RC-SBD is that it extends the existing SBD method to 3-D, and introduces an assumption of reflectivity smoothness in the horizontal direction, expressed by the total variation (TV) regularization term.The estimation is formulated as a minimization problem involving ℓ 2 and ℓ 1 norms and is optimized using the Split-Bregman algorithm.Tradeoff hyperparameters of the objective function are optimized via Bayesian optimization, maximizing the kurtosis of the calibrated volume image.Validation with synthetic data demonstrates accurate wavelet estimation and significant denoising of the volume image.Real-world data application further reveals considerable improvements in the channel-depth cross section, providing a clear visualization of structures like rebar and steel plates.Notably, the calibrated image remains stable across diverse datasets, including earthwork and bridge sections, showcasing the versatility and reliability of the proposed methodology.

I. INTRODUCTION
G ROUND penetrating radar (GPR) constitutes a techno- logical approach that enables the efficient acquisition of subsurface information by transmitting electromagnetic waves into the ground and measuring their reflections.With a long history of application in civil engineering, GPR has been instrumental in detecting buried pipes [1], [2], voids [3], [4], and damage [5], [6].The increasing demand for inspecting substantial volumes of infrastructure has fostered the need for high-speed and comprehensive surveys.Consequently, research into a system known as vehicle-mounted multichannel Takanori Imai is with the Department of Civil Engineering, Graduate School of Engineering, University of Tokyo, Tokyo 113-8656, Japan (e-mail: imai-takanori@g.ecc.u-tokyo.ac.jp).
Digital Object Identifier 10.1109/TGRS.2023.3317846GPR (MC-GPR) has seen considerable progression in recent years [7], [8], [9], [10], [11], [12], [13], [14].This system is characterized by the arrangement of multiple antennas orthogonal to the scan (or driving) direction, enabling efficient generation of 3-D images of the subsurface.Despite its innovative nature, the analysis of volume images, which are often described as the spatial amalgamation of multiple 2-D (B-scan) images [15], presents certain complexities.In particular, noise that appears streaky or random is frequently evident in the images generated by MC-GPR applications.This noise is largely attributable to the utilization of distinct antennas across individual channels.Past research applications have primarily centered around objects that are relatively conspicuous in the image, such as medium-sized cavities and buried pipes.
As such, the fundamental solution to noise has often been moved to the margins of these studies.Rigorous channelto-channel calibration of MC-GPR will be indispensable for exploring novel applications, such as the detection of small voids, rebar, and microdamage.The theory behind the noise generation can be succinctly explained through straightforward mathematical equations.Let us denote f and w as the vectors of observed waves and transmitted wavelets, respectively, and r as the vector representing subsurface reflectivity.For vehicle-mounted types, although not universally introduced in other studies, the stationary clutter g, which does not change with observation locations, should be considered.The clutter g represents direct waves and reflections from the vehicle body itself, and is a constant reflection independent of the exploration position.These variables are related by the following equation: where x is the observation point in the scan direction, c is the channel number, t is the reflection time depth and N is the Gaussian noise with a variation of σ 2 .The above relation is also shown as a conceptual diagram in Fig. 1.Given conditions of ideal propagation devoid of frequency-dependent attenuation, observations are made through the convolution of the wavelets and subsurface reflectivity.Due to minor disparities in wavelets between channels and the existence of the clutter, a discontinuous image is obtained in the channel direction, even if the reflectivity is smooth.Deconvolution by wavelets serves as an effective approach for inter-channel calibration and can be executed as follows:  where ˆdenotes the discrete Fourier transform (DFT) of the vectors.When the transmitted wavelet w is known, methods such as a Wiener filter [16], [17], [18] or deterministic deconvolution [19] prove efficient.However, the exact wavelet is often not known in various circumstances.Such an issue, termed blind deconvolution, generally represents an ill-posed problem that seeks to find both unknown reflectivity r and wavelets w.The goal of deconvolution varies, e.g., for calibration (as in this study) or input electric field estimation for full-waveform inversion (FWI) [20], [21], [22], or one may be interested in the obtained high-resolution reflectivity itself.In the realm of seismic data processing, methods based on the presumption of the wavelet being of minimum phase (spiking deconvolution) have been proposed [23], [24].Nevertheless, the wavelet is frequently nonminimum phase in GPR, leading to the proposition of deconvolution methods premised on the mixed-phase assumption [25], [26], [27].
Sparse modeling has gained momentum in the realms of image processing and machine learning as a technique underpinning the presumption of sparsity in information and facilitating the extraction of its essence [28], [29].This concept has found its way into the GPR field, where it is commonly recognized as sparse blind deconvolution (SBD).The SBD method seeks to identify higher-resolution reflectivity and enhanced wavelets by presupposing that subsurface reflectivity is sparse, an assumption grounded in practical reasoning.Past studies have demonstrated robust recovery of nonminimum phase wavelets and high-resolution reflectivity [30], [31], [32], [33].Li [31] proposed a technique that alternately updates the wavelet and reflectivity for superior results.They introduced a regularization term that integrates high order statistics [34] and the ℓ p norm, the effectiveness of which was corroborated using simulated data.Jazayeri et al. [33] further put forward an alternating update technique to effectuate more robust deconvolution.By solving the ℓ 2 − ℓ 1 minimization problem with unconstrained basis pursuit denoising (UBPDN) [35], their method realized a higher resolution of experimental GPR waveforms.
Despite their efficacy, the aforementioned methods have been designed for B-scan images acquired with a pair of antennas, while methods suitable for volume images measured with multichannel antenna pairs have yet to be proposed.A plausible approach would be to apply the SBD method independently for each channel and obtain the wavelets, but  OBSERVATION MODEL this presents challenges such as phase inconsistency in the depth direction.In other words, phenomena such as the unevenness or irregularity in the depth of the road surface can occur.Consequently, an appropriate 3-D extension of the method is needed.
This study puts forth an extended alternating update SBD method, adapted to three dimensions, for the extraction of wavelets from MC-GPR for volume image calibration.The optimization equation incorporates the assumption of consistent subsurface reflectivity.Specifically, the horizontal total variation (TV) of the reflectivity is integrated into the objective function.This study terms the proposed method reflectivity consistent SBD (RC-SBD).Additionally, the study acknowledges and addresses the stationary clutter engendered by vehicle-mounted radar.Section II elucidates the formulation and algorithm devised to solve the problem.Section III validates the results of the proposed method using synthetic data, while Section IV presents its application to real-world data.Conclusions are given in the last section.

II. METHODOLOGY A. Proposed Model
The variables utilized for constructing the model are summarized in Table I.These variables are either vectors or multidimensional arrays, simplistically referred to as tensors.The elements n d , n s , and n c represent the length of a volume image in the depth direction, scan direction, and number of channels, respectively.Data acquisition by onvehicle MC-GPR can be conceptualized according to the following equation: When W is constituted by a Toeplitz form, the first term can typically represent a linear convolution.Specifically, we adopt the stepped frequency continuous wave (SF-CW) system, which observes data in the frequency domain.Therefore, for each channel c, the tensor W is assumed to be a circulant matrix constituted as follows: (4) In this study, tensor products such as W R are defined to be the matrix product for each channel c.It follows that the calculation below is performed: ( The second term denotes a stationary clutter tensor that includes direct waves and reflections from mounting devices.We assume that this clutter is stationary and does not undergo substantial variations in the scan direction. The SBD of GPR, proposed in a previous study [33], identified sparse reflectivity and wavelet by solving the following equation: min R(:,:,c) Here, ∥x∥ p represents the ℓ p norm, which is defined by . λ r and λ w stand for the regularization hyperparameters.As simultaneous resolution of the minimization problem concerning R and W proved to be challenging, they split (6) into two equations, and resolved the problem for R and W alternately.In other words, reflectivity is initially obtained given the initial wavelet, followed by updating the wavelet utilizing the updated reflectivity.Reiteration of the aforementioned procedure led to minimizers, and the method was shown to perform well for sparse reflectivity estimation and transmitted wavelet estimation in the single-channel case.
We extend the method to 3-D by introducing the assumption that the reflectivity is smooth in the scan and channel directions.This is accomplished by solving the optimization problem min The above equation is also split into two following equations and solved alternately: The newly introduced regularization hyperparameters are denoted by λ tv and λ g .The first term common to ( 8) and ( 9) represents a fidelity term grounded on the observation model (3), serving to weaken the constraint.λ tv symbolizes TV, renowned as one of the most efficacious image denoising techniques [36], [37].The TV minimization method finds application within the GPR field for the purpose of image restoration or noise reduction [38], [39], and ensuring the stability of the FWI [40].
The TV term functions to minimize the magnitude of the derivative, essentially smoothing the image.Nevertheless, since it is analogous to the ℓ 1 regularization term, it has the capacity to preserve meaningful derivatives.Therefore, it is possible to maintain the effects of the ℓ 1 regularization term R, which induces sparsity in the reflectivity, while simultaneously ensuring smoothness of R. Two variants of the TV indices exist: isotropic and anisotropic.The latter, which is less prone to induce image distortion [41], is employed in our study and can be represented by the following equation: ∇ s and ∇ c are the derivative operator in the scan and channel direction.Here, the forward difference operator with Neumann boundary conditions is adopted, i.e.,

B. Split Bregman Algorithm
In the proposed optimization problem (8), the involvement of nondifferentiable functions such as the TV and the ℓ 1 -norm makes the problem challenging to solve.
Assuming that → R N is a convex and differentiable function, we consider the minimization problem given by ( 14) for a specific u ∈ R N .The constraints in ( 14) are relaxed to transform the problem into ( 15) The constraints are weakened by By introducing a auxiliary variable d ∈ R N , the problem can be reformulated as expressed in (16), and further transformed into The constraints are also relaxed in this step.It is known that by applying the Bregman iteration and then using the "adding the noise back" method, (17) can be further rewritten into an iterative formula as in Here, k denotes the iteration number and d g ∈ R N is a variable associated with the subgradient.Equation (18) can be "split" into two subproblems, as given in In (20), u is decoupled from the L1 norm, and if E is differentiable, the problem becomes simple to solve.Equation ( 21) is known to have a closed-form solution, given by Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

C. Iterative Optimization
Now, we consider a specific solution to the main problem using.In Section II-B, it was necessary to introduce two variables for the splitting optimization: an auxiliary variable and a variable related to the subgradient (they correspond to d k and d k g , respectively).To apply the split Bregman algorithm, new variables are introduced for the derivative in each direction of the image and for the ℓ 1 norm (Table II).Using these variables, main minimization problem (8) can also be expressed in the form of iteration as in ( 18) and ( 19) The solutions to the above subproblems must be followed sequentially.First, ( 24) is a simultaneous minimization problem for R and G.This can be solved using the gradient method, since L is fortunately convex and has only one solution.The partial derivative of L with respect to G is and ∂L/∂ G = 0 has a closed-form solution of The partial derivative with respect to R is and the gradient can be calculated using G obtained in (28).The next step is minimization with respect to alternative variables (25).This has the following closed-form solution with shrinkage: where, T = (( . The last step is the minimization with respect to the subgradient-related variable (26).Analogous to (19), they can be updated by the following equations: Using the sparse and smooth reflectivity R obtained from the above iterative process, the transmitted wavelet W can be estimated based on (9).This is done in a similar way to the derivation of the Wiener filter.The updating equation for each channel c is as follows: where * denotes the conjugate.The above procedure can solve the main problem (7) and is summarized in Algorithm 1. i rg and µ represent the iteration number and the learning rate of the gradient method, respectively.The gradient descent method is shown as an example in Algorithm 1, but methods such as momentum and Adam [42], which are expected to have fast or stable convergence, can of course be adopted.For the initial value of W , the experimentally observed direct waves can be used.The variables introduced for the split Bregman algorithm are initialized as tensors with all zero elements.

III. VALIDATION BY SYNTHETIC DATA A. Data Generation
The time signal of wavelets and clutter significantly varies depending on their vehicular installation, making these shapes exceedingly challenging to accurately determine by experiment.Consequently, we employed synthetic data with known true values for methodological validation.
The creation of synthetic data is based on the SF-CW approach, in accordance with the real-world data acquisition methods discussed later in Section IV.This method involves transmitting a sinusoidal wave with a frequency that changes at a fixed step, while measuring the amplitude and phase Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.(36) until convergence of W of the reflected waves.The discretized frequency used for transmission and observation is from 50 to 3030 MHz, with a step size of 20 MHz.The conversion to the time domain is conducted using the inverse discrete Fourier transform (IDFT).

Algorithm 1 Transmitted Wave Estimation
The wavelets were generated across 29 channels.Three distinct types of volume images were synthesized, each of which includes manually generated layered and hyperbolic responses.These were generated by the convolution of transmitted waves with 3-D synthetic sparse data in the frequency domain, followed by the addition of stationary clutter.20 dB Gaussian noise was also added to the volume image.For the algorithm's initial wavelets, random errors were imposed on the amplitude and phase spectra of the true transmitted waves.Within the algorithm, zero-padding increased the time waveform to n d = 2 10 points.Furthermore, the following normalization was applied to the wavelets at the onset of each iteration to avoid the explicit solution of W = 0: Two further modifications were implemented to stabilize the results.The first concerns a constraint on the updating of (36).The calibration process involves multiplying by the inverse of the wavelet, a procedure that can become unstable when F − G is small.To mitigate this, only the top half of the F − G amplitudes were used for updating at each frequency.The second modification was the application of smoothing in the frequency domain, which was accomplished using a Gaussian window.After convergence, the observed data F was calibrated in the frequency domain based on (2).We want to note that the calibration data obtained from ( 2) is determined in the frequency domain, and thus, it does not necessarily yield a sparse waveform.

B. Selection of Hyperparameters
The primary optimization problem (7) includes hyperparameters such as λ tv , λ r , λ w , λ , along with the hidden Split Bregman parameter α.These hyperparameters determine the trade-off between the terms of the objective function.Some of the hyperparameters affect the estimation results, hence careful consideration was needed.Typically fixed, the Split Bregman parameter is assigned a value of α = 5 in this case.Furthermore, λ tv is set to 1, as we discerned that excessively small λ tv values lead to nonsmooth calibrated images.The remaining hyperparameters exert a substantial influence on the estimation outcomes, thus searched in a certain range.Kurtosis is frequently employed as a metric for desirable waveforms [27], [43], [44].This study follows suit, utilizing the maximization of kurtosis of the calibrated image as a condition for the optimal hyperparameters, i.e., the maximization of the subsequent objective function: In this equation, µ and σ represent the mean and standard deviation of the random variable x, where the calibrated volume image substitutes for x.The validity of using kurtosis as an indicator is briefly validated here.Fig. 2 presents scatter plots of kurtosis against root mean squared error (RMSE) of the calibrated volume image and the wavelets for each set of hyperparameters.Each black point signifies a combination of λ r , λ g = 0.01, 0.1, 1, 10 and λ w = 10 −6 , 10 −5 , 10 −4 , 10 −3 , 10 −2 , translating to a total of 80 grid search outcomes.Within this range, the RMSE of the volume image exhibits a distinctly negative correlation.Similarly, the correlation coefficient between kurtosis and wavelets error also demonstrated a negative correlation at −0.77.Consequently, the maximization of kurtosis emerges as an apt indicator for determining the hyperparameters that yield accurate volume image with the wavelets also being well estimated at this point.
For practical considerations and to ascertain optimal parameters with fewer trials compared to grid search, a search utilizing Bayesian optimization was attempted [45].λ r and λ g were sought in the range from 0.01-10, and λ w was searched in the range from 10 −6 -10 −2 .The optimized results of a total of 30 trials are shown in Table III.Also, kurtosis and RMSE for this parameter set are indicated by the red dots in Fig. 2. Convergence to the vicinity of the optimal range by grid search was substantiated.

C. Calibration Results
Fig. 3 shows the results of the wavelet estimation using parameters derived from Bayesian optimization.In the time domain, the predicted waveform nearly aligns with the true waveform, signifying precise estimation [Fig.3(a)].Fig. 3(b) demonstrates the RMSE between the true wavelet and the estimated wavelet for each channel.The decreasing errors across all channels demonstrate the effectiveness of multichannel wavelet estimation using the proposed approach.
estimation result for the stationary clutter is exhibited in Fig. 4, where Fig. 4(a) represents the true and the estimated clutter, and Fig. 4(b) shows the residual between the two.While simple ℓ 2 norm adjustment does not yield a perfect estimate, the method manages to capture a general trend.Fig. 4(b) retains horizontal streaks around the travel time of 2, 5, and 10 ns.This can be attributed to the difficulty in distinguishing a strong and horizontally uniform layered response from stationary clutter in the scan direction.One way to mitigate this estimation error is to use a variety of data, including complex reflections.Fig. 5 presents the slices of the volume image after calibration.Fig. 5(a) is calibrated by the true wavelets and stationary clutter and Fig. 5(b) comprises the image calibrated according to (2), solely using the initial wavelets.Due to errors in the wavelets and deconvolution performed without eliminating the stationary clutter, Fig. 5(b) exhibits significant noise.This noise, resulting from the channel characteristics, is especially noticeable in the scan-channel (horizontal) and channel-depth cross sections, obscuring the subsurface responses.Although the scan-depth cross section appears cleaner than the two above as it utilizes a single antenna, it is still impacted by the stripe noise.Fig. 5(c) illustrates the volume image calibrated by the estimated wavelets and stationary clutter.All three sections reveal that noise has been effectively eliminated.This denoising effect is particularly conspicuous in the channel-depth cross section, where the two hyperbolic responses are clearly discernible.
Also for reference, the peak signal-to-noise ratio (PSNR) expressed by the following equation is shown above the volume data: where MSE is the Mean Squared Error of the noiseless and the target data, and R is the peak value of the data.Compared to the initial wavelet calibration, PSNR was improved by about 7 dB, and the noise improvement was confirmed in terms of the index.It can be concluded that the proposed method can precisely estimate the wavelets and clutter from the observed data based on the model in (3).Furthermore, calibration using these estimates yields smooth and highly discernible data.

IV. APPLICATION TO REAL-WORLD DATA A. Measurement System
The proposed method was applied to real-world data, collected from two vehicle-mounted MC-GPRs, as illustrated in Fig. 6.The vehicle in Fig. 6(a) was designated as Vehicle-1, while the one depicted in Fig. 6(b) was referred to as Vehicle-2.
These vehicles were equipped with radar systems that used the (SF-CW) method for the transmission and reception of electromagnetic waves.The operating frequency band spanned from 50 to 3030 MHz, with 20 MHz intervals.The data acquisition interval was set to 7 cm in the scan direction and 7.5 cm in the channel direction.These MC-GPRs were equipped with 29 and 25 channels respectively, allowing for data acquisition across widths of 2.1 and 1.8 m.The radar was of an air-coupled type, with the antenna positioned 15 cm above the ground.Direct waves were recorded with the antenna oriented toward the air in a preliminary experiment.These waves were then subtracted from the observed data before being used as inputs to the algorithm.The direct waves also provided the initial value for the transmitted wavelets.
Vehicle-1 was used to collect data from experimental environments that simulated the conditions of an earthwork road and a bridge slab, as shown in Fig. 7.The bridge slab site was embedded with steel bars at regular intervals of either 15 or 30 cm and selectively integrated steel plates to expose the characteristics of the floor slab's substructure.Conversely, Vehicle-2 was used on national roads located within Saitama Prefecture, which incorporated diverse sections including both earthwork and bridges.
Hyperparameters were determined by Bayesian optimization based on kurtosis maximization, similar to the previous chapter.Three training and three validation datasets were prepared for each vehicle.Each dataset contained 50 points in the scan direction, corresponding to a size of 3.5 m.The wavelets and clutter were estimated using the training data, and  these estimates were then used to calibrate the validation data according to (2).Hyperparameters maximizing the kurtosis of the calibrated validation data were sought.The Split Bregman and TV parameters were fixed to α = 10 and λ tv = 1.λ r and λ g were explored in the range from 0.01-10, and λ w was explored in the range from 10 −6 -10 −2 .The hyperparameter sets chosen through 30 iterations are presented in Table IV.Fig. 8 provides a comparative depiction of the waveforms of both the estimated and initial wavelets for representative channels, corresponding to the optimal parameter sets.
The computations were programed and executed in MAT-LAB on a computer equipped with an Intel Core i9-9980HK CPU and an NVIDIA GeForce RTX 2080 GPU.The average computation time for each pair of hyperparameters was  approximately 200 s, amounting to a total computation time of 2 h.While the computation may seem time-consuming, it is important to note that the estimated wavelets and clutter for frequency domain calibration can be utilized for any observed data, provided that the installation conditions on the vehicle remain unchanged.Therefore, computation time does not present a significant obstacle for practical applications.

B. Calibration Results
Fig. 9 presents a comparison between the volume image calibrated with direct waves and those calibrated with estimated wavelets by the RC-SBD for Vehicle-1.The top row of the figure displays results calibrated by subtracting direct waves and dividing by the same, while the bottom row presents results from the proposed method, calibrated by subtracting both direct waves and the estimated clutter, and dividing by the estimated wavelets.Fig. 9(a) and (d) provide an overall view of the volume image, revealing a significant improvement in the image quality, with the stripe noise observed in the synthetic data now eliminated.
The efficacy of the RC-SBD is particularly noticeable in the channel depth slice [Fig.9 hyperbolic response from densely placed rebars.Steel plates integrated into the underside of the bridge slab become visible around 8 ns in the depth direction.
The scan-channel cross section, or horizontal cross section, also sees significant improvement.Fig. 9(c) and (f) show the depth at which the rebar reflection was detected.Striped noise at 0.6 and 1.4 m in the channel direction overlap with reflections from the main bars in Fig. 9(c), obscuring them.However, these stripes are mitigated in Fig. 9(f), revealing clear, periodic rebar responses.The removal of antenna characteristics for each channel significantly improved the noise across all cross sections, resulting in more consistent volume image.
The application for even longer data is presented in Fig. 10.These figures show the results of calibrating in-service road image of around 1400 scan points (100 m) in the scan direction.For this vehicle, calibration with direct waves introduces unwanted low-frequency noise into the waveform.In contrast, calibration using the estimated wavelets and clutter effectively eliminates the presence of unwanted waves and enhances visibility.As mentioned previously, this volume image includes both the bridge and earthwork sections.Fig. 10(b) and (e) represent channel-depth slices of the bridge section.The hyperbolic curve caused by the rebar is clearly visible through 5-6 ns depth.Also, at around the depth of 7 ns, a response presumed to be a girder under the slab can be seen.Fig. 10(c) and (f) the slices of the earthwork section, where the boundary response of the road surface and the base is easily discernable following clutter removal.Notably, the quality of the image improves continuously over the entire 100 m span.Even with training data as short as 3.5 × 3 = 10.5 m, the estimates prove versatile enough to be applied to both earthwork and bridge sections.

V. CONCLUSION
This study introduced the RC-SBD methodology for the calibration and denoising of MC-GPR volume image.This approach represents an evolution of the existing SBD algorithm, previously used for B-scan images.An integral part of this evolution is the incorporation of a TV regularization term, which articulates the assumption of smooth subsurface reflectivity.In the context of synthetic data, the RC-SBD method proved highly efficient.The wavelets estimated by the method deviated minimally from the ground truth, which translated into a substantial enhancement in the quality of the resulting volume image.A key factor in this result was the automatic hyperparameter estimation via Bayesian optimization, which was aimed at maximization of kurtosis.This process successfully yielded optimal wavelets.Furthermore, the methodology demonstrated significant robustness when applied to real-world data.It shed light on the responses of rebar, steel plate, and strata in channel-depth cross sections, which had previously been unobserved.The process of subtraction and division using the estimated clutter and transmitted wavelets resulted in marked improvements in the volume image with a range of 100 m.This included both earthwork and bridge sections.A significant advantage of this approach is that the transmitted wave and clutter estimated for a specific vehicle can be generalized to any observed data, as long as the radar installation conditions remain constant.This insight suggests that the RC-SBD method offers an efficient and flexible tool for the calibration and denoising of MC-GPR volume image.
While this study has made significant progress in enhancing the calibration and denoising of MC-GPR volume image through the RC-SBD method, it also leaves several important aspects unaddressed.
1) The study approximated the amplitude of stationary clutter using the ℓ 2 norm.However, this approximation is a coarse estimate and is highly reliant on the training data used.This approach, while effective as a first approximation, may be insufficient for more accurate or nuanced applications.The introduction of more rigorous estimation techniques, such as statistical or stochastic methods, may provide a more robust and accurate measure of stationary clutter amplitude.
2) The current RC-SBD method does not accommodate the time-varying attenuation of the transmitted wave.This limitation is a departure from the physical reality of wave propagation, where attenuation is often frequency-dependent and varies with time.
The existing methodology essentially attempts to find an average wavelet that results in an optimal calibrated volume image.Incorporating a term representing frequency-dependent attenuation could potentially improve the accuracy of deconvolution or calibration procedures in future applications of the method.

Manuscript received 27
June 2023; revised 6 August 2023 and 12 September 2023; accepted 17 September 2023.Date of publication 22 September 2023; date of current version 9 October 2023.The work of Tsukasa Mizutani was supported by JST FOREST Program under Grant JPMJFR215R.(Corresponding author: Tsukasa Mizutani.)

2 .
Scatter plots of kurtosis of the calibrated image against.(a) RMSE of the calibrated volume image.(b) RMSE of the wavelets.

Fig. 3 .
Fig. 3. RC-SBD results of synthetic data regarding wavelet estimation (channel no.20).(a) Comparison of the truth, estimated, and initial wavelet.(b) RMSE convergence of the estimated wavelets by 15 × iteration.

Fig. 4 .
Fig. 4. Comparison of the stationary clutter in the time domain.(a) Truth and estimated for the representative channels.(b) Residual between the truth and the estimated.

Fig. 5 .
Fig. 5. Comparison of calibrated volume image.(a) Calibrated by the true transmitted wave and stationary clutter.(b) Calibrated by the initial transmitted wave.(c) Calibrated by the estimated transmitted wave and stationary clutter.

Fig. 7 .
Fig. 7. Pictures of the bridge slab field.(a) Panomaric view.(b) Rebars and the steel plate installed inside the slab.

Fig. 9 .
Fig. 9. Comparison of the calibrated volume image (vehicle-1).Top row (a)-(c) calibrated by the direct waves.Bottom row (d)-(f) calibrated by the estimated wavelets and the clutter.(a) and (d) Overall view of volume image.(b) and (e) Channel-depth slice images.(c) and (f) Scan-channel slice images.

Fig. 10 .
Fig. 10.Comparison of the calibrated volume image (vehicle-2).Top row (a)-(c) calibrated by the direct waves.Bottom row (d)-(f) calibrated by the estimated wavelets and the clutter.(a) and (d) Overall view of volume image.(b), (c), (e) and (f) Channel-depth slice images.
(a) and (d)], which utilizes multiple antennas for observation.Around 5 ns in the depth direction, rebar reflections can be seen, displaying the overlap of the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I DEFINITIONS
AND DESCRIPTIONS FOR THE

TABLE II DEFINITIONS
AND DESCRIPTIONS FOR THE SPLIT BREGMAN ALGORITHM prox is a shrinkage function: prox(X, λ) = sgn(X ) ⊙ max (|X | − λ, 0), where sgn and ⊙ refers to the signum function and the Hadamard product, respectively.

TABLE III OPTIMIZED
HYPERPARAMETERS FOR THE SYNTHETIC DATA

TABLE IV OPTIMIZED
HYPERPARAMETERS FOR THE REAL-WORLD DATA