A method for suppressing multiple waves of ground penetrating radar based on one-wave guiding weighted matching

Abstract

This invention discloses a weighted matched multiple suppression method for ground-penetrating radar (GPR) based on primary wave guidance, belonging to the field of GPR data processing technology. The method first performs regularized predictive deconvolution on preprocessed radar data to obtain an initial multiple model. Then, based on the energy distribution of the coarse primary wave estimate, a weighted mask and a diagonal weighted matrix are constructed. These are explicitly introduced into the weighted least squares solution of the matched filter to obtain the matched filter. Finally, the multiple estimation is optimized using this filter and subtracted from the original data to obtain the suppression result. This invention applies soft constraints to the matched filter through a primary wave-guided weighted mask, reducing the risk of excessive subtraction in the primary wave-dominant region and improving the directional suppression capability of the multiple wave-dominant region. It effectively suppresses multiples while maintaining the structural integrity of the primary wave, exhibiting advantages such as robust processing and low parameter sensitivity.

Description

Technical Field

[0001] This invention belongs to the field of ground penetrating radar data processing technology, and specifically relates to a weighted matching ground penetrating radar multiple wave suppression method based on primary wave guidance. Background Technology

[0002] Ground penetrating radar (GPR), as a highly efficient geophysical exploration technology, works by emitting high-frequency electromagnetic waves into the ground and receiving reflected echoes from interfaces between different media, thereby locating and identifying concealed underground targets. Due to its advantages such as high resolution, non-destructive operation, and ease of use, GPR technology is widely used in infrastructure exploration, geological surveys, and engineering quality inspection.

[0003] However, in practical engineering applications, GPR data quality is susceptible to various factors. Among them, multiple interference is one of the long-standing key problems in ground-penetrating radar data processing. Multiples often overlap with primary reflection events in both time and space, thereby reducing the signal-to-noise ratio of the radar profile and affecting the accurate identification of underground targets.

[0004] Several mature methods have been developed for addressing the multiple reflection problem in seismic data processing. Predictive deconvolution and adaptive matched filter subtraction are widely used for multiple reflection suppression and have formed a relatively systematic theoretical framework. With the continuous expansion of the application of ground penetrating radar in near-surface detection, these methods are also gradually being introduced into GPR data processing to reduce the impact of multiple reflections on imaging results.

[0005] However, due to the wide bandwidth and strong inhomogeneity of GPR data, and the high degree of superposition between primary and secondary waves, traditional methods often face problems of insufficient stability and excessive subtraction in practical applications, especially under conditions of high noise or model mismatch. To improve the numerical stability of the secondary wave suppression process, regularization has been introduced into solving inverse problems such as predictive deconvolution and matched filtering. Although regularization methods have improved the stability of secondary wave suppression to some extent, existing improvement methods mostly focus on optimizing the predictive operator or matched filter itself, lacking explicit constraints on preserving the primary wave structure, and still struggling to achieve an ideal balance between secondary wave suppression and primary wave protection.

[0006] For GPR data, the continuity and stability of the primary wave reflection structure are of equal importance to subsequent target identification and interpretation as multiple wave suppression. Therefore, it is necessary to introduce constraint information that can distinguish between the primary and multiple wave dominant regions during multiple wave suppression to achieve a reasonable balance between multiple wave suppression and primary wave protection. Summary of the Invention

[0007] To address the aforementioned technical problems, this invention provides a weighted matching ground-penetrating radar multiple suppression method based on primary wave guidance, aiming to suppress multiple residuals, avoid excessive subtraction of effective signals, and improve the signal-to-noise ratio and interpretability of radar profiles.

[0008] To achieve the above objectives, the technical solution of the present invention is as follows: A method for suppressing multiple waves from a weighted matching ground-penetrating radar based on primary wave guidance includes the following steps: Step S1: Collect ground-penetrating radar data and preprocess the data to obtain preprocessed radar data; Step S2: Perform regularized predictive deconvolution processing on the preprocessed radar data to obtain the initial multiple wave model; Step S3: Construct a coarse first-wave estimate based on the initial multiple wave model, generate a weighted mask according to the energy distribution characteristics of the coarse first-wave estimate, and construct a diagonal weighted matrix based on the weighted mask; Step S4: Explicitly incorporate the diagonal weighting matrix into the weighted least squares solution process of the matched filter to obtain the matched filter; Step S5: Optimize the initial multiple model using the matched filter to obtain an updated multiple estimate, and subtract the updated multiple estimate from the preprocessed radar data to obtain the result after multiple suppression.

[0009] In the above scheme, step S2 specifically includes: Step S21: Let the preprocessed radar data be... The predicted step size is α The prediction filter of length n is The predicted value is then expressed as: ; in, For the first One filter coefficient; This refers to the index number of the filter coefficients; Step S22: Define prediction error And construct the error energy objective function under the minimum mean square error: ; Step S23: Based on the least squares criterion, derive the Wiener-Hough equation by minimizing the error energy. Its matrix form is as follows: ; in, R For the reason autocorrelation function The constructed Toplitz matrix, The cross-correlation vector; Step S24: Introduce the first regularization parameter Construct the regularization equation: ; in, It is the identity matrix; Step S25: Solve the above regularization equation to obtain a stable prediction filter c; Step S26: Use the prediction filter c to process the preprocessed radar data. Perform convolution to obtain the initial multiple wave model. : .

[0010] In the above scheme, step S3 specifically includes: Step S31: Perform global amplitude matching on the initial multiple wave model using a global scalar factor, and subtract the amplitude-matched multiple wave model from the preprocessed radar data to obtain a rough first wave estimate; Step S32: Perform Hilbert transform on the coarse first-wave estimate to obtain its envelope and then normalize it; Step S33: Generate a binary mask based on the normalized envelope through threshold segmentation, and construct a weight function based on the binary mask; Step S34: Perform a moving average smoothing process on the weight function to obtain a smoothed weight function, and then construct a diagonal weighted matrix.

[0011] In a further technical solution, the specific calculation method for step S31 is as follows: ; ; in, As a global scalar factor, For the initial multiple wave model, For preprocessed radar data, This is a rough estimate of the first wave.

[0012] In a further technical solution, the specific calculation method for step S32 is as follows: ; ; in, For envelope, Let Hilbert transform function, For normalized envelope.

[0013] In a further technical solution, the specific calculation method for step S33 is as follows: ; ; in, It is a binary mask. For the preset threshold, For the weight function, This refers to the weighting parameters for the primary wave region.

[0014] In a further technical solution, the specific calculation method for step S34 is as follows: ; ; in, For smoothing weight function, The length of the moving average filter. This indicates the current time point, i.e., the center position of the sliding window; Relative to a point in time The offset; This represents the weight value at the nth offset within a window centered at n; It is a diagonal weighted matrix.

[0015] In the above scheme, step S4 specifically includes: The initial multiple wave model Constructed as a convolution matrix M, based on the diagonal weighting matrix Establish the weighted least squares objective function: ; Introducing a second regularization parameter The matched filter is obtained by solving the following weighted regularization equation: ; in, The convolution matrix constructed for the initial multiple wave model. It is a diagonal weighted matrix. This is the second regularization parameter. It is the identity matrix. For preprocessed radar data, This is a matched filter.

[0016] In the above scheme, step S5 specifically includes: The matched filter obtained using step S4 For the initial multiple wave model Perform convolution to obtain the optimized wave multiple estimation. and from the preprocessed radar data Subtracting this estimate from the middle yields the results of multiple wave suppression. .

[0017] Through the above technical solution, the weighted matching ground-penetrating radar multiple suppression method based on primary wave guidance provided by the present invention has the following beneficial effects: 1. Reconstruct the multiple-wave suppression process to improve the systematic nature of the processing mechanism. The traditional linear processing flow of "prediction-subtraction" is reconstructed into a closed-loop optimization flow of "prediction-estimation-weighted inversion-adaptive subtraction". Structural prior is introduced in the matched filter solution stage to improve multiple wave suppression from the algorithm mechanism level and avoid local optimization that only stays at the parameter or operator level.

[0018] 2. Protect the primary wave in-phase axis and reduce the risk of excessive subtraction of the effective signal. A weighted mask is constructed using the primary wave energy distribution. In the primary wave-dominant region, the fitting weight of the matched filter is reduced to avoid overfitting and over-subtraction of the effective signal, thus maintaining the continuity and stability of the primary wave phase axis.

[0019] 3. Enhance the fitting ability of multiple waves and improve the suppression effect of multiple waves. Maintaining a high weight in the multiple wave-dominant region enhances the matching filter's fitting ability in that region, thereby achieving directional suppression of the multiple wave and reducing residual energy and artifacts.

[0020] 4. Avoid the shortcomings of traditional methods and reduce the probability of effective wave damage. By embedding a weighted mask into the inversion solution process, the filter is able to distinguish between the primary and secondary wave-dominant regions during the solution stage, thereby reducing the risk of excessive subtraction and effective wave damage commonly found in traditional methods.

[0021] 5. Improve algorithm robustness and reduce parameter sensitivity. Through regularization constraints and weighted smoothing, the invention maintains stable processing performance even under conditions of strong noise or local overlap between primary and secondary waves. When the core parameters fluctuate within ±20%, the signal-to-noise ratio improvement fluctuates by less than ±0.21 dB, and the structural similarity index fluctuates by less than ±0.0015, indicating that the invention is insensitive to parameter changes.

[0022] 6. Improved quantitative evaluation indicators, superior to existing methods. The measured data and numerical simulation results show that the present invention outperforms the traditional predictive deconvolution method and L1 norm matching method in terms of quantitative evaluation indicators such as signal-to-noise ratio improvement and structural similarity index, and achieves a more reasonable balance between multiple wave suppression and primary wave structure preservation. Attached Figure Description

[0023] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below.

[0024] Figure 1 This is a schematic flowchart of a weighted matching ground-penetrating radar multiple wave suppression method based on primary wave guidance disclosed in an embodiment of the present invention; Figure 2 This is the original radar data measured in practice; Figure 3 The results are shown by different methods; (a) is the predicted deconvolution, (b) is L1 norm matching, and (c) is the weighted matching method of this invention. Figure 4 The simulation data and processing results are as follows: (a) is the original radar data of the simulation, (b) is the result after prediction deconvolution processing, (c) is the result after L1 norm matching processing, and (d) is the result after weighted matching processing of this invention. Figure 5 A comparative analysis of the weighted masking effect of different methods is presented; where (a) is the predicted deconvolution, (b) is the unweighted matched filtering, and (c) is the weighted matching method of this invention. Detailed Implementation

[0025] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention.

[0026] This invention provides a method for suppressing multiple waves from a weighted matching ground-penetrating radar based on primary wave guidance, such as... Figure 1 As shown, the method includes the following steps: I. Data Acquisition and Preprocessing First, raw radar data of the underground target area is collected using ground-penetrating radar (GPR). This radar data is a time-series signal arranged by channels, denoted as... , where L is the total number of lanes.

[0027] Each radar data point undergoes preprocessing, specifically including: DC removal: Subtract the mean value from each signal to eliminate system bias.

[0028] Background removal: Calculate the average channel for all channels and subtract the average channel from the original data for each channel to reduce horizontal continuous interference.

[0029] Bandpass filtering: Based on the frequency band range of the effective signal from the ground penetrating radar, the data is bandpass filtered to remove high-frequency noise and low-frequency drift.

[0030] After preprocessing, the preprocessed radar data is output and denoted as... .

[0031] II. Regularized Predictive Deconvolution for Obtaining the Initial Multiple Wave Model For preprocessed data The initial multiple wave model is obtained using a regularized prediction deconvolution method. The specific process is as follows: Step S21: Let the prediction step size be... The prediction filter length is The prediction filter is denoted as Then the predicted value Represented as: ; in, For the first One filter coefficient; This refers to the index number of the filter coefficients; Step S22: Define prediction error The objective function for error energy is: ; Step S23: Based on the least squares criterion, minimize the error energy to derive the Wiener-Hough equation: ; From the input signal Find the autocorrelation function The predictive filter can be obtained by solving the matrix equation. The above formula can be further written as: ; in, For the reason autocorrelation function The constructed Toplitz matrix, This is a cross-correlation vector.

[0032] Step S24: To solve the autocorrelation matrix To solve ill-conditioned unstable problems, Tikhonov regularization is introduced by adding a small first regularization parameter. By moving to the main diagonal of the matrix, we improve its condition number and construct a regularization equation: ; in, for Identity matrix.

[0033] Step S25: Solve the equation to obtain a stable prediction filter. .

[0034] Step S26: Using a prediction filter For raw data Perform convolution to obtain the initial multiple wave model. : ; In actual calculations, the convolution results need to be truncated or padded with zeros to maintain consistent length.

[0035] III. Constructing a coarse first-wave estimate and generating a weighted mask The goal of this step is to obtain the primary wave energy distribution characteristics and construct a weighting matrix to guide the matched filtering.

[0036] Step S31: Perform global amplitude matching on the initial multiple wave model using a global scalar factor, and subtract the amplitude-matched multiple wave model from the preprocessed radar data to obtain a rough first wave estimate. Details are as follows: Introducing a global scalar factor The initial multiple wave model was tested using the least squares criterion. Perform overall amplitude matching to quickly obtain the energy distribution characteristics of the primary wave, the expression of which is: ; This factor enables the scaled multiple model to match the energy level of the multiple components in the original data.

[0037] Subtracting the amplitude-matched multiple wave model from the original radar data yields a rough first wave estimate. : ; The estimate does not need to be very precise; it only needs to reflect the approximate energy distribution of the primary wave.

[0038] Step S32: Perform a Hilbert transform on the coarse first-wave estimate to obtain its envelope and then normalize it; specifically as follows: Rough first wave estimation Perform a Hilbert transform and calculate its envelope. : ; Then the envelope is normalized: ; in, For envelope, Let Hilbert transform function, For normalized envelope.

[0039] The normalized envelope value is compressed to between 0 and 1, reflecting the relative strength distribution of primary wave energy on the time axis.

[0040] Step S33: Generate a binary mask based on the normalized envelope through threshold segmentation, and construct a weight function based on the binary mask; specifically as follows: Set preset threshold (In this embodiment, we take) Generate a binary mask : ; The region with a value of 1 indicates a strong primary wave energy region, while the region with a value of 0 indicates a multiple wave or noise-dominated region.

[0041] To avoid excessive subtraction in the primary wave region during matched filtering, an initial weighting function is constructed based on the binary mask, assigning lower weights to regions where primary wave energy is dominant. ; in, The weighting parameter for the first wave region (taken in this embodiment) This function results in a lower weight for the primary wave dominance region. The region dominated by multiple waves maintains a high weight (weight 1).

[0042] It is important to emphasize that the weighted mask is not directly constructed based on the predicted multiple waves themselves, but rather derived from the amplitude envelope of a coarse estimate of the first wave. Its role is merely to impose soft constraints on the matched filtering process. Therefore, even if there is local energy overlap between the first and multiple waves, the mask will not forcibly classify the event as either a first or multiple wave, thus avoiding systematic misclassification problems. This first wave estimate is only used to extract energy distribution features rather than precise separation results; its purpose is to provide spatial constraints for the matched filtering process. Even if there are some misclassifications in local areas, the smoothing processing of the weighted mask and the introduction of regularization terms can effectively suppress the amplified impact of misclassifications on the filtering results, thereby ensuring the overall robustness of the algorithm.

[0043] Step S34: Perform a moving average smoothing process on the weight function to obtain a smoothed weight function, and then construct a diagonal weighted matrix, as follows: To ensure the continuity of the weighting function over time, for Use a length of Smoothing is performed using a moving average filter: ; Smoothing eliminates abrupt changes in weights at the boundaries, ensuring the continuity of the weighting function.

[0044] Finally, a diagonal weighted matrix is ​​constructed based on the smoothed weight function: ; in, For smoothing weight function, The length of the moving average filter. This indicates the current time point, i.e., the center position of the sliding window; Relative to a point in time The offset; This represents the weight value at the nth offset within a window centered at n; It is a diagonal weighted matrix.

[0045] This weighting matrix will be used in the weighted matched filter solution in subsequent steps.

[0046] IV. Solving the Matched Filter Using Weighted Matched Filtering This step is the core of the method, and its goal is to solve for the optimal matched filter under the guidance of a weighted mask.

[0047] First, the initial multiple wave model Constructed as a convolution matrix Then, the weighted least squares criterion is introduced to establish the objective function: ; in, Let be the matched filter to be solved.

[0048] To obtain a robust matched filter, the derivative of the objective function is taken and set to zero, and a second regularization parameter is introduced. The weighted regularization equation is obtained as follows: ; Solving this equation yields the optimal matched filter. .

[0049] In this step, the diagonal weighted matrix The function is as follows: In the region dominated by the primary wave, the matching filter's fit is weakened due to its lower weight, thus avoiding overfitting of the effective signal; in the region dominated by the secondary wave, the matching filter has a stronger fitting ability due to its higher weight, thus achieving directional suppression of the secondary wave.

[0050] V. Multiple Wave Suppression and Result Output Find the matched filter Then, it is applied to the initial multiple wave model. This yields optimized, more accurate multiple wave estimates: ; Finally, the optimized multiple wave estimates were extracted from the original radar record. Subtracting from the middle, we obtain the output signal after removing multiple waves: ; Output the radar data after the multiple wave suppression. This is used for subsequent underground target identification and interpretation.

[0051] In summary, the weighted matching method based on primary wave guidance proposed in this invention achieves a key breakthrough over the traditional prediction and subtraction framework. By introducing a primary wave energy sensing mechanism, a new intelligent processing flow is constructed, encompassing prediction, estimation, guided matching, and adaptive subtraction. This algorithm innovatively utilizes a weighted mask to deeply integrate prior information from primary wave estimation into the regularization solution process of the matched filter, thereby accurately guiding the algorithm to focus on the multiple wave development region and effectively avoiding the risk of effective wave damage.

[0052] Experiment and Results Analysis A single evaluation index is insufficient to fully reflect the effect of multiple wave suppression. Therefore, this invention employs a combination of quantitative indicators and profile structure comparison to comprehensively evaluate the processing performance of different methods. The classic predictive deconvolution method is selected for comparison to highlight the improved effect of the algorithm in this invention.

[0053] 1. Definition of Evaluation Indicators To objectively and quantitatively evaluate the processing effects of different multiple wave suppression methods, this invention evaluates the processing results from two aspects: signal energy improvement and structure preservation. Specifically, signal-to-noise ratio (SNR), SNR improvement (ΔSNR), and structural similarity index (SSIM) are used as evaluation indicators, and qualitative visual analysis of radar profiles is combined to comprehensively evaluate the algorithm performance.

[0054] The signal-to-noise ratio (SNR) is defined as the ratio of the average power in the signal region to that in the noise region, and its calculation formula is as follows: ; in, and These represent the average power of samples within the signal window and noise window, respectively. The signal region is selected within a time-space window containing major primary wave reflection events, while the noise region is selected within an area that does not contain significant reflection events and is dominated by multiple waves and background interference. To ensure fairness in comparing different methods, all methods use the same signal and noise window settings.

[0055] The signal-to-noise ratio (SNR) improvement is used to measure the change in overall signal quality before and after multiple suppression. It is defined as the difference between the processed SNR and the unprocessed SNR, i.e.: ; SSIM is used to compare the relative preservation of primary wave structure by different methods, rather than absolute consistency with the real underground model. This index takes into account the similarity between brightness, contrast and structural information, and its value ranges from [0,1]. The closer the value is to 1, the more similar the two images are in structure.

[0056] It should be noted that the structural similarity index (SSIM) is only applicable when a reference structure exists. Therefore, this invention uses SSIM to evaluate the processing results only in numerical simulation data, with the reference image being the known primary wave field in the simulation model. For measured data, due to the lack of a real primary wave reference structure, this invention does not use SSIM as a quantitative evaluation index, but mainly relies on the qualitative comparison of the signal-to-noise ratio improvement and radar profile structure for analysis.

[0057] 2. Standardized parameter description To ensure fairness in comparing different methods, this invention uses uniform settings for matched filtering parameters and weighted mask parameters in all experiments, as shown in Table 1. The prediction step size in predictive deconvolution... The operator length n is selected based on the main period characteristics of multiple waves in different data, and its value will be explained in the corresponding experimental section.

[0058] Table 1 Experimental parameter settings

[0059] 3. Parameter sensitivity analysis To verify the robustness of the proposed algorithm to the core parameters, this invention conducts a ±20% fluctuation test on key adjustable parameters to analyze the impact of parameter changes on the processing results. The test parameters and ranges are as follows: mask threshold T (0.24 / 0.3 / 0.36), prediction stage regularization parameter... (8e-4 / 1e-3 / 1.2e-3), Regularization parameters for the matching phase (8e-4 / 1e-3 / 1.2e-3), Mask smoothing window length (12 / 15 / 18) First wave weight (0.08 / 0.1 / 0.12).

[0060] Numerical simulation data test results show that when the parameters fluctuate within the above range, the method of the present invention... The maximum fluctuation amplitude was ±0.21dB, and the maximum fluctuation amplitude of SSIM was ±0.0015. This result shows that the algorithm of this invention is less sensitive to the core parameters. Even with small deviations in the parameters, it can still maintain stable multi-wave suppression effect and single-wave fidelity performance, verifying the robustness of the algorithm.

[0061] 4. Analysis of measured data results To verify the applicability of the proposed method in practical engineering scenarios, multiple wave suppression experiments were conducted using ground-penetrating radar data from road inspection. This data contained significant multiple wave interference, which affected the interpretation of underground structures. Figure 2 The original radar profile is presented, and it can be observed that the phase axes of multiple waves are more obvious in the shallow and locally deep regions.

[0062] Based on the main period characteristics of multiple waves in the measured road data, the prediction step size is determined. The parameters are set to 10ns, the prediction operator length n is set to 17ns, and the other parameters related to matched filtering and weighted mask are all set according to the unified settings given in Table 1. Figure 3 The processing results of the actual test data are given.

[0063] For the measured data, Figure 3 The predicted deconvolution and shown in (a) Figure 3 The L1 norm matching shown in (b) weakens the secondary echo energy in the 20–30 ns interval to some extent, but at the same time, it affects the primary wave reflection near 20 ns, resulting in local amplitude weakening and decreased continuity of the primary wave phase axis. In contrast, Figure 3 The method of the present invention shown in (c) effectively suppresses secondary echoes while better preserving the structural characteristics of primary wave reflections. The primary wave phase axis is more continuous and stable, the local reflection pattern is preserved, and no obvious excessive subtraction phenomenon is observed. A comprehensive comparison shows that although the three methods achieve similar improvements in overall energy performance, the method of the present invention achieves a more reasonable balance between primary wave protection and secondary wave suppression, thereby improving the structural interpretability of the radar profile.

[0064] To further conduct quantitative evaluation, Table 2 presents the signal-to-noise ratio (SNR) and the improvement in SNR after processing with different methods. The signal window for the measured data was set to a time range of 7-17 ns and a channel range of 18-28 ns. Fixed windows were used for both methods to ensure fairness in the comparison.

[0065] Table 2 Quantitative Analysis Results

[0066] The results show that the method of the present invention is superior to the traditional predictive deconvolution method and L1 norm matching method in terms of signal-to-noise ratio improvement. This indicates that under the same parameter settings and processing flow, the introduction of a weighted matching mechanism guided by a primary wave helps to reduce the impact of multiple waves and background interference on the effective signal, thereby improving the overall signal quality of the measured radar data.

[0067] It should be noted that, due to limitations in field conditions, the measured data lacks completely independent true information about the underground environment. This invention supplements and verifies this deficiency through numerical simulation and comparison with multiple methods.

[0068] 5. Analysis of Numerical Simulation Results To further verify the effectiveness of the proposed method and its ability to preserve the primary wave structure under controllable conditions, this invention constructs a three-dimensional numerical model based on the open-source electromagnetic simulation software gprMax, and compares and analyzes the suppression effects of different methods on multiple waves. Unlike measured data, the propagation paths and wavefield characteristics of primary and multiple waves can be clearly distinguished in the numerical simulation, thus providing a reliable reference for the mechanism verification of the algorithm's performance.

[0069] The simulation model is a 3D forward model (2.0m × 2.0m × 0.003m), with perfectly matched layer (PML) absorbing boundary conditions. The excitation source is a Hertzian dipole in the z-direction, 0.05m above the ground surface. The receiving and transmitting antennas are co-located, moving along the x-axis in 0.02m steps, collecting 90 data channels. The time window is 60ns, and the sampling rate is 0.1ns. Three typical dielectric materials are used in the model to simulate the actual detection environment: the upper layer is air (relative permittivity...). electrical conductivity S / m, thickness 0.2m); the middle surface layer is a concrete layer ( =8, =0.01S / m, thickness 0.3m); the lower layer is soil medium ( =5, =0.001S / m, thickness 1.5m). To simulate an underground target, an infinitely long ideal conductor (PEC) cylinder with a radius of 0.075m was embedded in the soil layer at a distance of 0.7m from the antenna.

[0070] Based on the known wavelet characteristics and multiple wave propagation paths in the simulation model, the prediction step size is determined. The operator length n was set to 6ns and 12ns respectively, and the other parameters related to matched filtering and weighted mask were all set according to the unified settings given in Table 1.

[0071] Figure 4 (a) shows the original radar data obtained from the simulation, in which multiple wave interference is evident in the shallow and some deep regions. Figure 4 As can be seen in (b), traditional predictive deconvolution has a certain effect on weakening the energy of multiple waves, but local remnants still exist. Figure 4 As shown in (c), compared to the predicted deconvolution, the L1 criterion is slightly unstable in some regions, with significant attenuation of target reflection and a slight degradation in lateral continuity. In contrast, Figure 4The method of the present invention shown in (d) suppresses multiple waves while making the target reflection event clearer, and its waveform continuity and energy concentration are well preserved.

[0072] Since the primary wave field is known in the simulation data, in addition to the improvement in signal-to-noise ratio, the structural similarity index is further used to quantitatively evaluate the results of different methods. The signal window for the numerical simulation data is 0-26 ns in time and 1-90 channels, while the noise window is 26-60 ns in time and 1-90 channels. Fixed windows are used for both methods to ensure fairness in the comparison. The quantitative analysis results are shown in Table 3.

[0073] Table 3 Quantitative Analysis Results

[0074] Numerical simulation results show that both predictive deconvolution and the method of this invention can effectively suppress multiple wave interference and maintain the main structural features of the primary wave, with a generally high SSIM value within the signal region. In comparison, the method of this invention performs slightly better in terms of signal-to-noise ratio improvement and structural similarity, indicating that the weighting strategy guided by the primary wave helps to achieve a more reasonable balance between multiple wave suppression and primary wave structure preservation.

[0075] 6. Comparative Analysis of Weighted Masking Effects To further clarify the role of the weighted mask in the method proposed in this invention, a comparative analysis was conducted on the matched filtering results with and without the introduction of the weighted mask based on numerical simulation data. In this comparative experiment, the predicted deconvolution parameters, matched filter length, and regularization parameters were kept consistent. The only variable was whether or not a primary-guided weighted mask was introduced during the matched filtering solution process, thus ensuring that the comparison results could truly reflect the influence of the weighted mask on the matched filtering process. Figure 5 A comparison is given between the matched filtering results without introducing a weighted mask and the processing results of the method of this invention.

[0076] Figure 5 A comparison of the results obtained by different methods on the same numerical simulation data is presented. It can be seen that... Figure 5 The predictive deconvolution method shown in (a) weakens the multiple wave energy to some extent, but significant residues remain. In contrast, Figure 5 The unweighted matched filter in (b) achieved a more significant effect in further suppressing multiple waves, but due to the lack of constraints on the primary wave-dominant region, its matching process exhibited slight overfitting in local areas, manifested as residual artifacts in the background and longitudinal instability.

[0077] This invention introduces a primary-wave guided weighted mask during matched filtering, applying soft constraints to the primary-wave energy-dominant region, effectively suppressing the excessive subtraction problem in unweighted matched filtering. Since the weighted mask is introduced only as a soft constraint, combined with smoothing and regularization terms, its local misclassification does not systematically drive the improvement of the matching results, thus avoiding the amplified impact of mask bias on the results. Figure 5 As can be seen in (c), the primary wave phase axis maintains better continuity and smoothness in the transverse direction, while the background area is more uniform and no obvious residual artifacts appear.

[0078] The comparative results show that predictive deconvolution has limited effectiveness in weakening multiples, while unweighted matched filtering, although further enhancing multiple suppression, is prone to overfitting in local regions. The method of this invention, by introducing a weighted mask guided by the primary wave, imposes constraints on the matching process, effectively avoiding excessive subtraction, making the primary wave phase axis more continuous and smooth, and the background region more stable, demonstrating its advantage in balancing primary wave protection and multiple wave suppression.

[0079] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined in this invention may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for suppressing multiple waves from a weighted matching ground-penetrating radar based on primary wave guidance, characterized in that, Includes the following steps: Step S1: Collect ground-penetrating radar data and preprocess the data to obtain preprocessed radar data; Step S2: Perform regularized predictive deconvolution processing on the preprocessed radar data to obtain the initial multiple wave model; Step S3: Construct a coarse first-wave estimate based on the initial multiple wave model, generate a weighted mask according to the energy distribution characteristics of the coarse first-wave estimate, and construct a diagonal weighted matrix based on the weighted mask; Step S4: Explicitly incorporate the diagonal weighting matrix into the weighted least squares solution process of the matched filter to obtain the matched filter; Step S5: Optimize the initial multiple model using the matched filter to obtain an updated multiple estimate, and subtract the updated multiple estimate from the preprocessed radar data to obtain the result after multiple suppression.

2. The method according to claim 1, characterized in that, Step S2 specifically includes: Step S21: Let the preprocessed radar data be... The predicted step size is α The prediction filter of length n is The predicted value is then expressed as: ； in, For the first One filter coefficient; This refers to the index number of the filter coefficients; Step S22: Define prediction error And construct the error energy objective function under the minimum mean square error: ； Step S23: Based on the least squares criterion, derive the Wiener-Hough equation by minimizing the error energy. Its matrix form is as follows: ； in, R For the reason autocorrelation function The constructed Toplitz matrix, The cross-correlation vector; Step S24: Introduce the first regularization parameter Construct the regularization equation: ； in, It is the identity matrix; Step S25: Solve the above regularization equation to obtain a stable prediction filter c; Step S26: Use the prediction filter c to process the preprocessed radar data. Perform convolution to obtain the initial multiple wave model. : 。 3. The method according to claim 1, characterized in that, Step S3 specifically includes: Step S31: Perform global amplitude matching on the initial multiple wave model using a global scalar factor, and subtract the amplitude-matched multiple wave model from the preprocessed radar data to obtain a rough first wave estimate; Step S32: Perform Hilbert transform on the coarse first-wave estimate to obtain its envelope and then normalize it; Step S33: Generate a binary mask based on the normalized envelope through threshold segmentation, and construct a weight function based on the binary mask; Step S34: Perform a moving average smoothing process on the weight function to obtain a smoothed weight function, and then construct a diagonal weighted matrix.

4. The method according to claim 1, characterized in that, The specific calculation method for step S31 is as follows: ； ； in, As a global scalar factor, For the initial multiple wave model, For preprocessed radar data, This is a rough estimate of the first wave.

5. The method according to claim 4, characterized in that, The specific calculation method for step S32 is as follows: ； ； in, For envelope, Let Hilbert transform function, For normalized envelope.

6. The method according to claim 5, characterized in that, The specific calculation method for step S33 is as follows: ； ； in, It is a binary mask. For the preset threshold, For the weight function, This refers to the weighting parameters for the primary wave region.

7. The method according to claim 6, characterized in that, The specific calculation method for step S34 is as follows: ； ； in, For smoothing weight function, The length of the moving average filter. This indicates the current time point, i.e., the center position of the sliding window; Relative to a point in time The offset; This represents the weight value at the nth offset within a window centered at n; It is a diagonal weighted matrix.

8. The method according to claim 1, characterized in that, Step S4 specifically includes: The initial multiple wave model Constructed as a convolution matrix M, based on the diagonal weighting matrix Establish the weighted least squares objective function: ； Introducing a second regularization parameter The matched filter is obtained by solving the following weighted regularization equation: ； in, The convolution matrix constructed for the initial multiple wave model. It is a diagonal weighted matrix. This is the second regularization parameter. It is the identity matrix. For preprocessed radar data, This is a matched filter.

9. The method according to claim 1, characterized in that, Step S5 specifically includes: The matched filter obtained using step S4 For the initial multiple wave model Perform convolution to obtain the optimized wave multiple estimation. and from the preprocessed radar data Subtracting this estimate from the middle yields the results of multiple wave suppression. .

Images