GPR echo signal denoising method based on improved WTD and VMD

Abstract

The invention relates to the technical field of geophysical exploration signal processing, and discloses a GPR echo signal denoising method based on improved WTD and VMD. Comprising the following steps: adaptively optimizing a decomposition modal number and a penalty factor of the VMD by utilizing an HHO algorithm; performing 2D-VMD decomposition on the original GPR data by using the optimal parameter to obtain a plurality of IMFs; dividing the IMF into a signal IMF and a noise IMF through a correlation coefficient; a wavelet threshold parameter is optimized and improved by using a PSO algorithm, and the noise IMF is processed; and after reconstruction, direct waves are removed by adopting an averaging method. Through a dual parameter optimization mechanism, an improved threshold function with continuity and no deviation and a differential threshold processing strategy, adaptive high-precision denoising of a GPR signal is realized, the signal-to-noise ratio is effectively improved, a weak target signal is reserved, direct wave interference is remarkably suppressed, and the method has the advantages of high robustness and high robustness. The method has the advantages of high parameter adaptive capacity, high denoising precision, high engineering practicability and the like.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration signal processing technology, and in particular to a method for denoising GPR echo signals based on improved WTD and VMD. Background Technology

[0002] Ground penetrating radar (GPR) is a geophysical exploration device that uses high-frequency electromagnetic waves to detect underground structures. Due to its advantages such as short detection cycle, high resolution, and non-destructive testing, it is widely used in fields such as road defect detection, tunnel lining assessment, and urban underground space exploration. However, in actual detection, GPR echo signals are often interfered with by various types of noise, mainly including: direct waves (composed of antenna-coupled waves and surface direct waves), system electronic noise (such as amplifier noise and sampling quantization noise), and background clutter (caused by medium inhomogeneity). These noises reduce the signal-to-noise ratio of the GPR signal, obscure target reflection characteristics, and severely affect subsequent data interpretation.

[0003] In existing research, wavelet threshold denoising (WTD) and variational mode decomposition (VMD) are two commonly used GPR signal denoising methods.

[0004] WTD utilizes wavelet multi-resolution analysis and time-frequency localization characteristics to suppress noise by thresholding the high-frequency coefficients of the signal, which can balance noise suppression and signal detail preservation to a certain extent. However, the traditional WTD method has the following shortcomings: (1) Limitations of the threshold function: The hard threshold function is discontinuous and easily causes reconstruction oscillations, while the soft threshold function is continuous but can lead to excessive smoothing of the signal; (2) Inaccurate threshold setting: High-frequency wavelet coefficients often contain effective signals, but traditional methods treat them all as noise for threshold calculation, resulting in thresholds that are too high or too low, affecting the denoising effect; (3) Single threshold processing strategy: Existing methods use the same threshold for high-frequency coefficients with different decomposition levels and directions, without considering the characteristics of noise decay with the number of decomposition levels, which can easily lead to over-suppression.

[0005] Virtual Mode Decomposition (VMD) achieves efficient frequency component separation by decomposing a signal into several adaptive mode functions, thus effectively solving the mode aliasing problem in Empirical Mode Decomposition (EMD) and gradually becoming another important means of denoising GPR signals. However, the decomposition quality of VMD is highly dependent on the selection of two key parameters: the number of decomposed modes and the penalty factor. Traditional methods often rely on manual experience or trial and error to determine these parameters, making it difficult to guarantee the decomposition effect and easily leading to problems such as under-decomposition of modes, redundant modes, or mode aliasing, resulting in a difficulty in simultaneously achieving noise suppression and signal detail preservation.

[0006] In summary, existing GPR signal denoising methods have not yet achieved an effective balance between parameter adaptability, noise suppression accuracy, and signal detail preservation. There is an urgent need for a joint denoising method that can achieve parameter adaptive optimization, balance noise suppression and signal fidelity, and effectively remove direct waves. Summary of the Invention

[0007] To address the aforementioned technical problems, this invention provides a GPR echo signal denoising method based on improved WTD and VMD, which aims to adaptively optimize denoising parameters, improve signal-to-noise separation accuracy, preserve details of weak target signals, and effectively suppress direct wave interference.

[0008] To achieve the above objectives, the technical solution of the present invention is as follows: A method for denoising GPR echo signals based on improved WTD and VMD includes the following steps: Step 1: Obtain two-dimensional noisy GPR data, preprocess it to obtain one-dimensional noisy data for parameter optimization; Step 2: Using the HHO algorithm, with the goal of minimizing the envelope entropy, optimize the one-dimensional noisy data to obtain the optimal VMD parameter combination, including the optimal number of decomposition modes K and the optimal penalty factor α. Step 3: Using the optimal number of decomposition modes K and the optimal penalty factor α, perform two-dimensional variational mode decomposition on the original two-dimensional noisy GPR data to obtain K IMFs; Step 4: Calculate the correlation coefficient between each IMF and the original two-dimensional noisy GPR data, and classify the K IMFs into signal IMFs and noise IMFs according to the preset threshold; Step 5: For each noisy IMF, use the PSO algorithm with the generalized cross-validation (GCV) criterion as the fitness function to optimize and improve the parameters of WTD, and use the optimized parameters to perform improved wavelet threshold denoising on the noisy IMF to obtain the processed noisy IMF. Step 6: Superimpose the processed noise IMF and the signal IMF to reconstruct the signal and obtain the reconstructed signal; Step 7: Perform direct-to-direct-wave removal processing on the reconstructed signal to obtain the final denoised signal.

[0009] In the above scheme, step 2 includes: Step 2.1: Initialize the HHO algorithm parameters, including population size, maximum number of iterations, search space boundary, and randomly generate the initial population, with each individual representing a candidate [K, α] combination; Step 2.2: For each individual in the population, use the corresponding [K, α] to perform VMD decomposition on the one-dimensional noisy data to obtain K IMFs, calculate the envelope entropy of each IMF, and take the minimum envelope entropy value as the fitness value of the individual; Step 2.3: Based on the exploration and development strategy of the HHO algorithm, iteratively update the population position until the maximum number of iterations or the fitness value converges, and output the globally optimal individual position, that is, the optimal number of decomposition modes K and the optimal penalty factor α.

[0010] In the above scheme, the correlation coefficient is calculated in step 4 as follows: ; in, For the first The correlation coefficients between each IMF component and the original two-dimensional noisy GPR data. This is the original two-dimensional noisy GPR data; This represents the average value of the original two-dimensional noisy GPR data; For the first One IMF component; For the first The average value of each IMF component; This represents the number of sampling points.

[0011] In the above scheme, in step 4, the preset threshold is 0.5. IMFs with a correlation coefficient greater than the threshold are determined as signal IMFs, and those with a correlation coefficient less than the threshold are determined as noise IMFs.

[0012] In a further technical solution, step 2.2, the formula for calculating the envelope entropy is as follows: ; in, For the normalized envelope signal, For the first The envelope entropy of an IMF, This is the length of the envelope sequence, i.e., the signal length.

[0013] In a further technical solution, in step 2.3, the escape energy... The strategy used to control the switching of the HHO algorithm has a value that decays with iteration, calculated as follows: ; in, Use random initial values. This represents the current iteration number. This represents the maximum number of iterations. when When it enters the exploration phase, when It will then enter the utilization stage.

[0014] In the above scheme, step 5 specifically includes the following methods: Step 5.1: For each noise IMF, perform multi-scale discrete wavelet decomposition to obtain the high-frequency coefficients of each layer; Step 5.2: Initialize the PSO algorithm parameters and set the search dimension, which corresponds to the thresholds in the horizontal, vertical, and diagonal directions of each layer and the threshold function adjustment factor of each layer; use minimizing the GCV function as the fitness function, iteratively update the particle position, and find the optimal thresholds in the horizontal, vertical, and diagonal directions of each layer and the optimal adjustment factor of each layer. Step 5.3: Improve threshold processing The optimal thresholds for each layer in the horizontal, vertical, and diagonal directions, as well as the optimal adjustment factors for each layer, obtained through PSO optimization, are substituted into the improved wavelet threshold function and used for thresholding of the wavelet high-frequency coefficients of the corresponding layer. Step 5.4: Wavelet Reconstruction Finally, the high-frequency coefficients and low-frequency coefficients after thresholding are reconstructed by wavelet to obtain the processed noise IMF.

[0015] In a further technical solution, step 5.2, the fitness function is as follows: ; in, The number of high-frequency wavelet coefficients. This represents the number of wavelet coefficients that are set to 0 after thresholding. These are wavelet high-frequency coefficients. For the threshold, To estimate the wavelet coefficients.

[0016] In a further technical solution, the improved wavelet threshold function in step 5.3 is expressed as follows: ; in, These are wavelet high-frequency coefficients. For the threshold, To estimate wavelet coefficients, As a regulating factor, and .

[0017] In the above scheme, the method of step 7 is as follows: ; in, To reconstruct the signal in Pixel value at that location, for The pixel value after noise reduction processing The number of rows, i.e., the number of time sampling points. This refers to the number of columns, or the number of scan channels.

[0018] Through the above technical solution, the GPR echo signal denoising method based on improved WTD and VMD provided by the present invention has the following beneficial effects: (1) Strong parameter self-adaptation capability, freeing it from dependence on human experience. This invention constructs a dual-parameter optimization mechanism by introducing the Harris Eagle Optimization (HHO) algorithm to adaptively optimize the number of decomposition modes K and the penalty factor α in Variational Mode Decomposition (VMD), and simultaneously employing the Particle Swarm Optimization (PSO) algorithm to adaptively select the threshold and threshold function adjustment factor for wavelet threshold denoising. This mechanism overcomes the limitations of traditional methods that rely on manual experience or trial-and-error to determine parameters, enabling the algorithm to automatically adjust parameters based on the statistical characteristics of the GPR signal. This significantly improves the algorithm's generalization ability and robustness under different detection scenarios and noise environments.

[0019] (2) The denoising accuracy is significantly improved and the signal-to-noise separation effect is excellent. The improved wavelet threshold function proposed in this invention combines continuity, asymptoticity, and unbiasedness, solving the dual defects of traditional hard threshold functions (which are prone to reconstruction oscillations due to discontinuity) and soft threshold functions (which are prone to excessive signal smoothing due to constant bias). Simultaneously, the introduction of an adjustment factor makes the threshold function flexible, allowing dynamic adjustment between hard and soft thresholds based on signal characteristics. Furthermore, a differentiated thresholding strategy is adopted, setting thresholds and processing high-frequency coefficients independently for different decomposition levels and different directions (horizontal, vertical, diagonal), fully considering the directional anisotropy of noise distribution. Simulation experiments show that, compared with traditional soft and hard thresholding denoising algorithms, the HVPI algorithm improves the signal-to-noise ratio by [percentage missing]. and The mean square error was reduced respectively and ; (3) The signal details are well preserved, and weak targets are protected. This invention distinguishes between signal IMFs and noise IMFs obtained from VMD decomposition using the correlation coefficient method. Instead of simply removing the noise IMFs, it employs an improved wavelet threshold denoising process. This innovative strategy effectively preserves the weak target reflection signals submerged in the noise IMFs, avoiding the signal detail loss problem caused by directly removing noise components in traditional methods, thus achieving an effective balance between noise suppression and signal fidelity.

[0020] (4) Direct wave suppression is effective, and shallow targets are clearly distinguishable. This invention employs the mean method as the final step in the denoising process. Utilizing the spatial stability of direct waves, direct wave suppression is achieved by subtracting the mean value from each row. This method is computationally simple and highly effective, effectively eliminating the masking effect of strong-energy direct waves on shallow target signals, making shallow reflection characteristics clearly discernible in the denoised signal, and improving the accuracy of GPR data interpretation.

[0021] (5) The algorithm has a complete process and strong engineering applicability. This invention organically combines HHO parameter optimization, 2D-VMD decomposition, correlation coefficient classification, PSO-optimized wavelet threshold denoising, improved threshold function processing, and mean-based direct wave removal to form a complete adaptive denoising process for GPR signals. Each module has a clearly defined function and is tightly integrated. It can be used as the core algorithm of data processing software in a laboratory environment, or it can be ported to an embedded platform for real-time processing of GPR data, demonstrating strong engineering practical value. Attached Figure Description

[0022] 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.

[0023] Figure 1 This is a schematic diagram of a GPR echo signal denoising method based on improved WTD and VMD disclosed in an embodiment of the present invention. Figure 2 A flowchart for global optimization of the decomposition mode number and penalty factor of VMD using HHO; Figure 3 A flowchart illustrating the parameter optimization and improvement of WTD using the PSO algorithm. Detailed Implementation

[0024] 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.

[0025] This invention provides a GPR echo signal denoising method based on improved WTD and VMD (i.e., HVPI algorithm), the process of which is as follows: Figure 1As shown. The specific implementation steps are as follows: Step 1: Data Preprocessing Input raw two-dimensional GPR data Where I represents the number of time sampling points (rows) and J represents the number of spatial scanning channels (columns). To reduce the computational complexity of subsequent HHO parameter optimization, the two-dimensional data is averaged column-wise to obtain a one-dimensional average data vector, i.e., one-dimensional noisy data. This one-dimensional data retains the main frequency characteristics of the original signal and is used for fitness evaluation of the HHO algorithm.

[0026] Step 2: Adaptive optimization of VMD parameters based on HHO This step utilizes the Harris Eagle Optimization (HHO) algorithm to globally optimize the number of decomposition modes K and the penalty factor α in VMD. The flowchart is as follows: Figure 2 As shown. The specific implementation is as follows: Step 2.1: Initialize HHO parameters Set the maximum number of iterations. Search dimensions (correspond and Upper bound of search space Set as The lower realm for Randomly generate the initial population Each individual This represents a potential solution in the search space, namely , Population size.

[0027] Step 2.2: Calculate the fitness value: For each potential solution VMD decomposition is performed on the input one-dimensional noisy data to obtain Each intrinsic mode function (IMF) is then calculated. Fitness value and determine prey location This represents the current optimal solution. The fitness function is defined as the minimum value of the envelope entropy of each IMF, as shown below: .

[0028] The formula for calculating envelope entropy is as follows: ; in, For the normalized envelope signal, For the first The envelope entropy of each modal component, The length of the envelope sequence. A smaller envelope entropy indicates stronger regularity of modal components and less noise. Record the current globally optimal individual, i.e., the individual with the lowest fitness.

[0029] Step 2.3: Iterative Optimization For each iteration, update the escape energy. The value used to control the phase switching of the HHO algorithm decays with each iteration, as shown below: ; in, Use random initial values. This represents the current iteration number.

[0030] when When it enters the exploration phase, when It will then enter the utilization stage.

[0031] (1) Exploration phase (global search): when At that time, based on random numbers Choose two eagle flock position update strategies.

[0032] when At that time, the position is updated according to strategy 1, as shown below: ; in, For randomly selected individuals, This is a random number used to control the movement step size.

[0033] when At that time, the position is updated according to strategy 2, as shown below: ; in, The average position of the population. It is a random number. and These are the upper and lower bounds of the search space, respectively.

[0034] (2) Utilization stage (partial development): when At that time, based on the escape energy and random number Four attack strategies were selected, and the introduction of these four strategies significantly improved the flexibility of local search.

[0035] when and When the prey still has energy but the probability of escape is low, a soft encirclement strategy is executed, as shown below: ; in, This represents the random jump intensity, used to simulate the random jumping behavior of a flock of eagles when prey escapes. It is a random number.

[0036] when and When the prey is exhausted and unable to escape, a hard encirclement strategy is employed, directly tightening the encirclement to attack, as shown below: ; when and When the prey is likely to escape, a soft encirclement and approach dive are executed. The flock of eagles uses Levy Flight to simulate the prey's irregular movements and optimize the attack path, as shown below: ; in, For search dimensions; For random step size vectors ( (V) is used to adjust the amplitude of Levi's flight; The random step size for the Levy flight is used to enhance the diversity of the local search; For the fitness function, select a better solution.

[0037] when and When the prey's energy is low but escape is still possible, a hard-hit approach with a dive is executed. The eagle flock searches intensively near the prey, combined with Levi's flight optimization, as shown below: .

[0038] Step 2.4: Iteration Termination When the maximum number of iterations is reached or the fitness function value converges, the algorithm terminates its iteration and outputs the global optimal solution. That is, the optimal decomposition mode number obtained by optimization. and optimal penalty factor Otherwise, return to step 2.2 and continue iterating.

[0039] Step 3: 2D-VMD Decomposition The optimal decomposition mode number obtained in step 2 and optimal penalty factor Two-dimensional variational mode decomposition (2D-VMD) is performed on the original two-dimensional noisy GPR data. 2D-VMD decomposes the two-dimensional signal into... A narrowband two-dimensional intrinsic mode function (2D-IMF) is obtained, where each IMF represents a component at a different frequency scale. After decomposition, the following is obtained: Two-dimensional IMF components.

[0040] Step 4: Differentiate IMFs using the correlation coefficient method Calculate the correlation coefficients between each IMF and the original two-dimensional noisy data: ; in, For the first The correlation coefficients between each IMF component and the original two-dimensional noisy GPR data. This is the original two-dimensional noisy GPR data; This represents the average value of the original two-dimensional noisy GPR data; For the first One IMF component; For the first The average value of each IMF component; This represents the number of sampling points.

[0041] Depending on the actual signal-to-noise ratio, the threshold can be adjusted between 0.3 and 0.7. In this embodiment, the threshold is set to 0.5, and IMFs with correlation coefficients greater than the threshold are determined to be signal IMFs, while those with correlation coefficients less than the threshold are determined to be noise IMFs.

[0042] Step 5: Improved wavelet thresholding for noise IMF denoising Since noisy IMFs still contain some effective signal components, improved wavelet thresholding denoising is performed on them instead of direct removal, in order to further preserve the details of the target signal. Specifically: multi-scale discrete wavelet decomposition is performed on each noisy IMF. The PSO algorithm is used with the Generalized Cross Validation (GCV) criterion as the fitness function to optimize the thresholds and threshold function adjustment factors in each horizontal, vertical, and diagonal direction. These are then substituted into the improved wavelet threshold function for thresholding, followed by wavelet reconstruction to obtain the processed noisy IMF. The specific process is as follows: Step 5.1: Multi-scale discrete wavelet decomposition The db4 wavelet basis was selected, and the number of decomposition levels was set to 3. Two-dimensional wavelet decomposition was performed on each noise IMF component to obtain three levels of high-frequency coefficients in the horizontal, vertical, and diagonal directions, and one level of low-frequency coefficients.

[0043] Step 5.2: Adaptive optimization of PSO parameters The Particle Swarm Optimization (PSO) algorithm is used to optimize the threshold and threshold function adjustment factor for high-frequency coefficients in each layer and direction. The PSO algorithm flowchart is shown below. Figure 3 As shown, the process includes the following: First, initialize the PSO algorithm parameters: population size is set to 50; maximum number of iterations T is set to 300; inertia weight... The learning factor was gradually reduced from 0.9 to 0.2 during the search process, balancing exploration and development; Since the optimization parameters include thresholds for the horizontal, vertical, and diagonal directions of each layer. , and and threshold function adjustment factors for each layer Therefore, the search dimension is equal to four times the number of wavelet decomposition layers, i.e. .

[0044] The fitness value of each particle is then evaluated, with the GCV function minimized as the fitness function, as shown below: ; in, The number of high-frequency wavelet coefficients. This represents the number of wavelet coefficients that are set to 0 after thresholding. These are wavelet high-frequency coefficients. For the threshold, To estimate the wavelet coefficients.

[0045] By minimizing the GCV value, the threshold λ can be asymptotically optimally selected, and its denoising effect approaches the theoretical optimal solution as the amount of data increases.

[0046] After obtaining the fitness value of each particle, the optimal position of that individual particle can be obtained. and global optimal position Next, the velocity and position of the particles are updated using the following formula and constrained within the boundary range.

[0047] ; ; in, The velocity of an individual particle; The position of an individual particle; Indicates the particle number.

[0048] Repeat the above steps until the maximum number of iterations is reached or the fitness function value converges. This will give you the thresholds for the horizontal, vertical, and diagonal directions of each layer, as well as the threshold function adjustment factor.

[0049] Step 5.3: Improve threshold processing The parameters obtained from PSO optimization , and With regulatory factors Substitute it into the improved threshold function and use it for thresholding the high-frequency coefficients of the corresponding layer wavelet.

[0050] Based on the problems of discontinuity and fixed deviation in traditional threshold functions, it is known that: the threshold function should be continuous to prevent oscillations in the reconstructed signal; the wavelet coefficients before and after thresholding should not deviate to prevent excessive smoothing of the signal or loss of the target signal; the actual situation that coefficients above the threshold still contain some noise should be considered, and they should not be directly retained without processing; an adjustment factor should be introduced into the threshold function to make it flexible and adaptive, able to adapt to the characteristics of different signals and improve the denoising effect. Based on this, this invention proposes an improved wavelet threshold function, the expression of which is: ; in, These are wavelet high-frequency coefficients. For the threshold, To estimate wavelet coefficients, As a regulating factor, and .

[0051] This function has the following properties: (1) Continuity: Improved threshold function in It is continuous at the point of intersection, thus solving the discontinuity problem of the hard threshold function; (2) Asymptoticity: The asymptote of the improved threshold function is It has asymptoticity; (3) Unbiasedness: Improved threshold function and There is no deviation between them, thus eliminating the constant deviation of the soft threshold function; (4) Flexibility: By changing the regulatory factors The value can be selected from the soft threshold ( ) and hard threshold ( The threshold function can be dynamically adjusted between different ranges to adapt to the characteristics of different signals.

[0052] Step 5.4: Wavelet Reconstruction Finally, the high-frequency coefficients and low-frequency coefficients after thresholding are reconstructed by wavelet to obtain the processed noise IMF.

[0053] Perform the above operations on all noisy IMFs.

[0054] Step 6: Signal Reconstruction The processed noise IMF obtained in step 5 is superimposed with the signal IMF determined in step 4 to obtain the reconstructed signal.

[0055] Step 7: Using the mean method to obtain the direct wave The reconstructed signal is subjected to direct-to-transmission processing to obtain the final denoised signal. The calculation formula is as follows: ; in, To reconstruct the signal in Pixel value at that location, for The pixel value after noise reduction processing The number of rows, i.e., the number of time sampling points. This refers to the number of columns, or the number of scan channels.

[0056] Simulation experiments verify: To verify the effectiveness of this invention, a comparative experiment was conducted using simulated GPR data and traditional denoising methods. The simulated data was generated based on the Finite-Domain Difference (FDTD) method, simulating a scenario with underground targets containing cavities, and Gaussian white noise and direct wave interference were added to achieve an input signal-to-noise ratio (SNR_in) of -5 dB. Soft-threshold wavelet denoising, hard-threshold wavelet denoising, VMD denoising (with manually selected parameters), and the HVPI algorithm of this invention were used for processing. The denoising effect was evaluated using the signal-to-noise ratio (SNR) and mean square error (MSE).

[0057] Experimental parameters: HHO population size 30, 50 iterations; PSO population size 50, 300 iterations; db4 wavelet, 3-level decomposition; correlation coefficient threshold 0.5. The processing results are shown in Table 1.

[0058] Table 1. Comparison of denoising performance of different methods

[0059] As shown in Table 1, the HVPI algorithm of this invention has the highest output signal-to-noise ratio and the smallest mean square error, significantly outperforming traditional methods. Experimental results demonstrate that this invention can effectively suppress noise, preserve signal details, and successfully remove direct waves, verifying the algorithm's superiority.

[0060] 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 herein 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 denoising GPR echo signals based on improved WTD and VMD, characterized in that, Includes the following steps: Step 1: Obtain two-dimensional noisy GPR data, preprocess it to obtain one-dimensional noisy data for parameter optimization; Step 2: Using the HHO algorithm, with the goal of minimizing the envelope entropy, optimize the one-dimensional noisy data to obtain the optimal VMD parameter combination, including the optimal number of decomposition modes K and the optimal penalty factor α. Step 3: Using the optimal number of decomposition modes K and the optimal penalty factor α, perform two-dimensional variational mode decomposition on the original two-dimensional noisy GPR data to obtain K IMFs; Step 4: Calculate the correlation coefficient between each IMF and the original two-dimensional noisy GPR data, and classify the K IMFs into signal IMFs and noise IMFs according to the preset threshold; Step 5: For each noisy IMF, use the PSO algorithm with the generalized cross-validation (GCV) criterion as the fitness function to optimize and improve the parameters of WTD, and use the optimized parameters to perform improved wavelet threshold denoising on the noisy IMF to obtain the processed noisy IMF. Step 6: Superimpose the processed noise IMF and the signal IMF to reconstruct the signal and obtain the reconstructed signal; Step 7: Perform direct-to-direct-wave removal processing on the reconstructed signal to obtain the final denoised signal.

2. The GPR echo signal denoising method based on improved WTD and VMD according to claim 1, characterized in that, Step 2 includes: Step 2.1: Initialize the HHO algorithm parameters, including population size, maximum number of iterations, search space boundary, and randomly generate the initial population, with each individual representing a candidate [K, α] combination; Step 2.2: For each individual in the population, use the corresponding [K, α] to perform VMD decomposition on the one-dimensional noisy data to obtain K IMFs, calculate the envelope entropy of each IMF, and take the minimum envelope entropy value as the fitness value of the individual; Step 2.3: Based on the exploration and development strategy of the HHO algorithm, iteratively update the population position until the maximum number of iterations or the fitness value converges, and output the globally optimal individual position, that is, the optimal number of decomposition modes K and the optimal penalty factor α.

3. The GPR echo signal denoising method based on improved WTD and VMD according to claim 1, characterized in that, In step 4, the correlation coefficient is calculated as follows: ； in, For the first The correlation coefficients between each IMF component and the original two-dimensional noisy GPR data. This is the original two-dimensional noisy GPR data; This represents the average value of the original two-dimensional noisy GPR data; For the first One IMF component; For the first The average value of each IMF component; This represents the number of sampling points.

4. The GPR echo signal denoising method based on improved WTD and VMD according to claim 1, characterized in that, In step 4, the preset threshold is 0.

5. IMFs with a correlation coefficient greater than the threshold are identified as signal IMFs, and those with a correlation coefficient less than the threshold are identified as noise IMFs.

5. A method for denoising GPR echo signals based on improved WTD and VMD according to claim 2, characterized in that, In step 2.2, the formula for calculating the envelope entropy is as follows: ； in, For the normalized envelope signal, For the first The envelope entropy of an IMF, This is the length of the envelope sequence, i.e., the signal length.

6. The GPR echo signal denoising method based on improved WTD and VMD according to claim 2, characterized in that, In step 2.3, escape energy The strategy used to control the switching of the HHO algorithm has a value that decays with iteration, calculated as follows: ； in, Use random initial values. This represents the current iteration number. This represents the maximum number of iterations. when When it enters the exploration phase, when It will then enter the utilization stage.

7. The GPR echo signal denoising method based on improved WTD and VMD according to claim 1, characterized in that, The method in step 5 specifically includes: Step 5.1: For each noise IMF, perform multi-scale discrete wavelet decomposition to obtain the high-frequency coefficients of each layer; Step 5.2: Initialize the PSO algorithm parameters and set the search dimension, which corresponds to the thresholds in the horizontal, vertical, and diagonal directions of each layer and the threshold function adjustment factor of each layer; use minimizing the GCV function as the fitness function, iteratively update the particle position, and find the optimal thresholds in the horizontal, vertical, and diagonal directions of each layer and the optimal adjustment factor of each layer. Step 5.3: Improve threshold processing The optimal thresholds for each layer in the horizontal, vertical, and diagonal directions, as well as the optimal adjustment factors for each layer, obtained through PSO optimization, are substituted into the improved wavelet threshold function and used for thresholding of the wavelet high-frequency coefficients of the corresponding layer. Step 5.4: Wavelet Reconstruction Finally, the high-frequency coefficients and low-frequency coefficients after thresholding are reconstructed by wavelet to obtain the processed noise IMF.

8. A method for denoising GPR echo signals based on improved WTD and VMD according to claim 7, characterized in that, In step 5.2, the fitness function is as follows: ； in, The number of high-frequency wavelet coefficients. This represents the number of wavelet coefficients that are set to 0 after thresholding. These are wavelet high-frequency coefficients. For the threshold, To estimate the wavelet coefficients.

9. A method for denoising GPR echo signals based on improved WTD and VMD according to claim 7, characterized in that, The improved wavelet threshold function in step 5.3 is expressed as follows: ； in, These are wavelet high-frequency coefficients. For the threshold, To estimate wavelet coefficients, As a regulating factor, and .

10. A method for denoising GPR echo signals based on improved WTD and VMD according to claim 1, characterized in that, The method for step 7 is as follows: ； in, To reconstruct the signal in Pixel value at that location, for The pixel value after noise reduction processing The number of rows, i.e., the number of time sampling points. This refers to the number of columns, or the number of scan channels.

Images