Ground Penetrating Radar Signal Denoising Method Combining Two-Dimensional VMD and DT-CWT

Abstract

The present invention relates to a ground-penetrating radar signal denoising method that combines two-dimensional VMD with DT-CWT. The method first performs one-dimensional VMD on the original two-dimensional ground-penetrating radar noisy data in the time dimension, removes the intrinsic modal components corresponding to the noise according to corresponding judgment criteria, and then reconstructs the data. The reconstructed data is then subjected to one-dimensional VMD in the spatial dimension, and the intrinsic modal components corresponding to the noise are removed and then reconstructed. The reconstructed data after the two-dimensional VMD denoising process is subjected to a discrete dual-tree complex wavelet transform, and a corresponding threshold formula is selected for threshold denoising. The data is then reconstructed by inverse dual-tree complex wavelet transform to obtain the denoised data. This method takes into account the noise distribution of the data in the spatial dimension, performs two-dimensional VMD denoising on the original noisy data, and effectively improves the data signal-to-noise ratio. At the same time, to further improve the data signal-to-noise ratio, a dual-tree complex wavelet transform that can achieve precise reconstruction is used to perform secondary denoising on the data.

Description

Technical Field

[0001] The present invention belongs to the technical field of ground penetrating radar signal processing, and specifically relates to a time-frequency analysis of ground penetrating radar data using variational mode decomposition (VMD) and dual-tree complex wavelet transform (DT-CWT), and more particularly to a ground penetrating radar signal denoising technology using a combination of two-dimensional VMD and DT-CWT. Background Art

[0002] Ground-penetrating radar (GPR) is a geophysical exploration method that uses high-frequency pulsed electromagnetic waves to detect and locate subsurface geological bodies. It boasts strong anti-interference capabilities, flexible and convenient operation, high-resolution imaging, intuitive results, and high efficiency and non-destructiveness. It is widely used in geological surveys and engineering investigations. In actual exploration, wideband recording is used to obtain more reflected wave information. This means that while the effective signal is recorded, various interfering noise signals are also recorded, resulting in a decrease in the signal-to-noise ratio (SNR) of the GPR data. Therefore, noise reduction processing must be performed on the collected GPR data to improve the SNR and data quality, providing clear and reliable GPR profiles for subsequent geological interpretation.

[0003] Because ground-penetrating radar signals exhibit typical time-varying and non-stationary characteristics, conventional time-domain or frequency-domain analysis is inappropriate and cannot achieve ideal denoising results. Time-frequency analysis is a signal processing method that describes the time-varying frequency of a signal. It maps a one-dimensional time signal onto a two-dimensional time-frequency plane and analyzes the signal in the time-frequency domain. This method comprehensively reflects the signal's time-frequency characteristics and is particularly suitable for processing non-stationary signals. Conventional time-frequency analysis methods include the short-time Fourier transform (STFT), wavelet transform (WT), empirical mode decomposition (EMD), its improved EMD-like methods, and variational mode decomposition (VMD). While the short-time Fourier transform (STFT) can extract frequency characteristics of local time segments by selecting a window function, this method cannot achieve both high time-frequency resolution and high resolution. The wavelet transform (WT) can analyze signals at different time and frequency scales by shifting and scaling the scaling function and the wavelet function, exhibiting superior time-frequency characteristics. The wavelet's automatic zooming function can effectively separate sudden changes and high-frequency noise from the signal, achieving signal denoising. However, wavelet bases are empirically derived. When using wavelet transforms for signal time-frequency analysis and noise reduction, the selection of wavelet basis functions, decomposition levels, and threshold functions directly affects the results, making them non-adaptive. Empirical mode decomposition (EMD) decomposes signals based on the time-scale characteristics of the data, eliminating the need for pre-defined basis functions and enabling adaptive signal analysis. However, the modal components (IMFs) derived from EMD decomposition lack a strict mathematical definition, and this method is susceptible to noise, resulting in modal aliasing and end-point effects. While EMD decomposition has been continuously researched and improved, it still has some shortcomings. Variational mode decomposition (VMD), based on EMD, defines intrinsic mode components (IMFs) with more stringent constraints and uses an iterative approach to search for the optimal solution to the variational model, thereby determining the center frequency and bandwidth of each IMF component. This method has a comprehensive mathematical theory and is an adaptive and fully non-recursive method for modal variation and signal time-frequency analysis.

[0004] During actual ground-penetrating radar (GPR) detection, the collected GPR data is relatively complex due to the complexity of the detection environment and electromagnetic wave propagation processes; the diverse, non-uniform, and highly random distribution of underground media; the small geometric scale of the detection target; and factors such as natural interference and cluttered echoes from man-made structures. In actual data processing, a single time-frequency analysis method often fails to achieve effective denoising results, necessitating a combination of different time-frequency analysis methods to analyze and reduce noise in GPR signals. Chinese patent CN113238190A proposes a GPR echo signal denoising method based on EMD combined with wavelet thresholding. This method improves the signal-to-noise ratio (SNR) of GPR data by combining EMD decomposition with wavelet thresholding techniques. The combination of different time-frequency analysis denoising methods has also been applied to signal denoising in other fields. Qiao Yun et al. published a seismic signal denoising method based on VMD and improved wavelet thresholding. They proposed combining VMD with the improved wavelet thresholding method and applying it to seismic signal denoising. By constructing a new threshold function, they improved the shortcomings of traditional soft and hard thresholding functions and, combined with the good adaptability of the VMD method, further improved the signal-to-noise ratio of seismic signals. Liu Chong et al. published a method for denoising partial discharge signals by combining VMD with improved wavelet thresholding. Combining the respective advantages of VMD and wavelet thresholding, they first filtered out periodic narrowband interference signals through VMD decomposition and then used the improved wavelet thresholding method to filter out Gaussian white noise.

[0005] According to the characteristics of ground-penetrating radar (GPR) data, effective signals at common depth points between adjacent channels exhibit strong correlation, while noise exhibits no correlation. For two-dimensional GPR data, the noise contained in the data, such as random noise and noise caused by medium randomness, often exists in both the temporal and spatial dimensions. When combining various denoising methods to denoise measured data, previous studies have focused solely on the temporal dimension, suppressing noise in one-dimensional time signals while ignoring noise in the spatial dimension. Furthermore, while wavelet analysis has excellent localization in the time-frequency domain and can effectively analyze the time-frequency characteristics of non-stationary signals, conventional discrete wavelet transforms (DWTs) typically use real wavelet basis functions, which can only extract signal information from the perspective of amplitude and lack phase information. Complex wavelets, on the other hand, can extract information from both amplitude and phase perspectives. However, conventional complex wavelets are generally continuous wavelets with unique phase-frequency characteristics, and cannot accurately reconstruct signals. Summary of the Invention

[0006] The purpose of the present invention is to provide a ground penetrating radar signal denoising method combining two-dimensional VMD and DT-CWT for noisy two-dimensional ground penetrating radar data, so as to solve the problem of performing two-dimensional VMD and DT-CWT denoising on the original data in the time dimension and spatial dimension respectively according to the distribution characteristics of the two-dimensional ground penetrating radar noise, which can effectively improve the signal-to-noise ratio of the data.

[0007] The purpose of the present invention is achieved through the following technical solutions:

[0008] A ground penetrating radar signal denoising method combining two-dimensional VMD and DT-CWT includes the following steps:

[0009] a. Obtain the original GPR data and represent the original two-dimensional data with M time sampling points and N as s(x, t), where x represents the horizontal direction and t represents the time direction. The time sampling sequence of the i-th (i=1,…,M) channel can be recorded as s(x i ,t), the spatial sampling sequence at the jth (j=1,…,N) time sample point is recorded as s(x,t j );

[0010] b. Perform the first dimension VMD on the time dimension, initialize the decomposition mode number K and penalty factor α of the VMD algorithm, and loop for each s(x i ,t) Perform one-dimensional VMD decomposition in the time dimension to obtain the intrinsic mode (IMF) components of each signal. After completing the one-dimensional VMD decomposition of all channels, the IMF components of all channels are obtained, which are recorded as imf1;

[0011] c. Calculate the correlation coefficients between each IMF component of imf1 and the original data, and denote the correlation coefficient of the kth (k = 1, 2, ..., K) IMF component as R k , where the largest correlation coefficient is R max , define R k <R max The IMF component of / 10 is noise, R k >R max The IMF component of / 10 is a valid signal. The IMF component corresponding to the valid signal is screened and the signal is reconstructed to obtain the two-dimensional data s1(x, t);

[0012] d. Perform the second dimension VMD on the spatial dimension. Similar to step b, initialize the decomposition mode number K' and penalty factor α' of the VMD algorithm, cycle the time sampling points, and for each time sampling point s(x,t j) Perform one-dimensional VMD decomposition in the spatial dimension to obtain the IMF components of each time-space sampling signal, complete the one-dimensional VMD decomposition of all time sampling points, and obtain the IMF components of all time sampling points, which are recorded as imf2;

[0013] e. Calculate the correlation coefficients of each order IMF component of imf2 and the original data respectively, and denote the mutual correlation coefficient of the mth (m=1,2,…,K') IMF component as R m , where the largest correlation coefficient is R' max Similar to step c, define R m <R' max The IMF component of / 10 is noise, R m >R' max The IMF component of / 10 is a valid signal. The IMF component corresponding to the valid signal is screened and the signal is reconstructed to obtain the two-dimensional data s2(x,t);

[0014] f. Use dual-tree complex wavelet transform (DT-CWT) to perform threshold denoising on the reconstructed data s2(x, t) to obtain the denoised data s3(x, t).

[0015] Furthermore, in step b, one-dimensional VMD is performed on each channel of the original data s(x, t) in the time dimension, where the i-th channel s(x i ,t) includes the following steps:

[0016] b1. Establish the variational mode decomposition constraint equation,

[0017]

[0018] Where u k is the kth IMF component, ω k for u k The center frequency, K is the mode number, j(j 2 =-1) is the imaginary unit, ∑u k =s(x i ,t) is a constraint condition, which ensures that the sum of all modal components after VMD decomposition is the original signal;

[0019]

[0020] Indicated by u k The complex analytic signal composed of its Hill transform, δ(t) is the Dirac impulse function; * represents convolution, multiplying the complex analytic signal by The term adjusts the spectrum of each modal component to its corresponding baseband and calculates the L2 norm square of the demodulated gradient To estimate each modal component u kbandwidth;

[0021] b2. To find the optimal solution to the constrained variational problem, the Lagrangian multiplication operator λ(t) and the second-order penalty factor α are introduced to transform the constrained variational problem (1) into an unconstrained problem. Its Lagrangian expansion expression is:

[0022]

[0023] The Lagrange multiplier λ(t) can ensure the strictness of the constraints, and the second-order penalty factor α can ensure the accuracy of signal reconstruction in Gaussian noise environment;

[0024] b3. Use the alternating direction multiplier method (ADMM) to continuously update each modal component and its center frequency, and finally obtain the saddle point of the unconstrained model, which is the optimal solution to the original problem. Combined with the Parseval Fourier isometric transform principle, the n+1th frequency domain ADMM iteration format of the kth order IMF component is:

[0025]

[0026]

[0027]

[0028] Where ω represents the frequency. s(x i ,t), Fourier transform of λ(t); τ is the noise tolerance parameter, set the convergence error ε, initialize Iterative update according to the iterative format (3) to (5) ω k , When the following iterative convergence conditions are met, the iteration is stopped and the corresponding k-th order IMF component is output ω k ;

[0029]

[0030] Furthermore, in step f, dual dual-tree complex wavelet transform is used to further denoise the reconstructed data after two-dimensional VMD denoising. The specific steps are as follows:

[0031] f1. Construct a one-dimensional dual-tree complex wavelet function. Assume that the dual-tree complex wavelet consists of two parallel wavelet trees A and B, which respectively constitute the real part and the imaginary part. The scaling function corresponding to tree A is and the wavelet function ψ h (t), the scaling function corresponding to tree B and the wavelet function ψ g(t), construct complex wavelet

[0032] ψ(t)=ψ h (t)+jψ g (t) (7)

[0033] In order to make its spectrum unilateral, the wavelet function ψ is required h (t) and ψ g (t) are Hilbert transform pairs, that is, they satisfy

[0034] ψ g (t)=H{ψ h (t)} (8)

[0035] Here, H{·} represents the Hilbert transform.

[0036] f2. Construct a dual double-tree complex wavelet filter. Let h0(n) and h1(n) be the low-pass and high-pass filters corresponding to the real part tree A, respectively. Let g0(n) and g1(n) be the low-pass and high-pass filters corresponding to the imaginary part tree B, respectively. There needs to be a half-sample delay between the low-pass filters of tree A and tree B, that is,

[0037] g0(n)=h0(n-0.5) (9)

[0038] Correspondingly, in the frequency domain,

[0039]

[0040] Where G0(ω) and H0(ω) are the Fourier transforms of g0(n) and h0(n), respectively, which are two conjugate orthogonal low-pass filter banks. Formula (10) is called the half-frame shift condition, which is equivalent to increasing the sampling of the low-pass signal to twice the original sampling at each scale;

[0041] f3. Set the number of decomposition layers and perform a two-dimensional dual-tree complex wavelet transform on the reconstructed data s2(x, t) after VMD denoising. The two-dimensional dual-tree complex wavelet transform is obtained by the tensor product of the one-dimensional dual-tree complex wavelet. The decomposition is equivalent to performing a one-dimensional dual-tree complex wavelet transform on the rows and columns of the two-dimensional data respectively, that is, performing a one-dimensional dual-tree complex wavelet transform on the time dimension and the space dimension of the cascade of s2(x, t);

[0042] f4, wavelet coefficient threshold processing, set the threshold T, select the threshold function, and perform wavelet soft threshold denoising on the wavelet coefficients of each layer. Subtract T from the wavelet coefficients whose absolute value is greater than the threshold T, and set the wavelet coefficients whose absolute value is less than the threshold T to zero. The soft threshold function is expressed as:

[0043]

[0044] The threshold T is

[0045]

[0046] d u,v represents the vth wavelet coefficient of the uth layer, and V is the number of wavelet coefficients in the uth layer; is the processed wavelet coefficient, and the noise standard deviation σ can be estimated using the empirical formula, that is,

[0047]

[0048] f5. Data reconstruction: perform two-dimensional dual-tree complex wavelet inverse transform on the processed low-frequency wavelet coefficients and high-frequency wavelet coefficients to obtain the denoised data s3(x,t).

[0049] Compared with the prior art, the present invention has the following beneficial effects:

[0050] The ground penetrating radar signal denoising method combining two-dimensional VMD and DT-CWT in the present invention performs two-dimensional VMD and DT-CWT denoising on the original data in the time dimension and space dimension respectively according to the distribution characteristics of two-dimensional ground penetrating radar noise, which can effectively improve the signal-to-noise ratio of the data.

[0051] It has the following advantages:

[0052] 1. By performing two-dimensional VMD on the original data in both the time and space dimensions, the present invention can better analyze the characteristics of the two-dimensional ground penetrating radar signal in both dimensions. At the same time, the time sampling sequence and the space sampling sequence are denoised in both the time and space dimensions, thus avoiding the insufficient noise suppression caused by conventional one-dimensional VMD denoising only in the time dimension.

[0053] 2. The present invention constructs a dual tree discrete complex wavelet and a corresponding conjugate orthogonal filter pair that are mutually Hilbert transform pairs, so that the phase-frequency characteristics of the wavelet are optimally matched with those of the processed signal. The characteristics of the signal can be analyzed from the perspectives of amplitude and phase, and the dual tree discrete complex wavelet transform can achieve accurate reconstruction of the signal.

[0054] 3. The present invention combines the respective advantages of VMD and dual-tree complex wavelet threshold denoising to perform denoising on GPR data. It can not only analyze the time-frequency characteristics of two-dimensional signals from different angles, but also further improve the signal-to-noise ratio of GPR data. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the embodiments. It should be understood that the following drawings only illustrate certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other relevant drawings can be obtained based on these drawings without paying any creative work.

[0056] Figure 1 Overall flow chart of the ground penetrating radar signal denoising method combining two-dimensional VMD and DT-CWT in the present invention;

[0057] Figure 2 Homogeneous layered medium model;

[0058] Figure 3 Forward simulation data of uniform layered media;

[0059] Figure 4 Layered random equivalent medium model;

[0060] Figure 5 Layered random equivalent medium forward simulation data

[0061] Figure 6 1D-VMD reconstructed data, SNR = 9.36;

[0062] Figure 7 2D-VMD reconstructed data, SNR = 12.47;

[0063] Figure 8 2D-VMD+DT-CWT reconstructed data, SNR=14.85. DETAILED DESCRIPTION

[0064] The present invention will be further described below in conjunction with embodiment:

[0065] The present invention will be further described in detail below with reference to the accompanying drawings and examples. It will be understood that the specific embodiments described herein are intended only to illustrate the present invention and are not intended to limit the present invention. It should also be noted that, for ease of description, the accompanying drawings only illustrate portions relevant to the present invention, not all structures.

[0066] It should be noted that similar reference numerals and letters represent similar items in the following drawings. Therefore, once an item is defined in one drawing, it does not need to be further defined or explained in subsequent drawings. At the same time, in the description of the present invention, the terms "first", "second", etc. are used only to distinguish the description and should not be understood as indicating or implying relative importance.

[0067] The ground penetrating radar signal denoising method combining two-dimensional VMD and DT-CWT of the present invention specifically comprises the following steps:

[0068] a. Obtain the original GPR data and represent the original two-dimensional data with M time sampling points and N as s(x,t), where x represents the horizontal direction and t represents the time direction. The time sampling sequence of the i-th (i=1,…,M) channel can be recorded as s(x i ,t), the spatial sampling sequence at the jth (j=1,…,N) time sample point is recorded as s(x,t j );

[0069] b. Perform VMD on the time dimension. Initialize the decomposition mode number K and penalty factor α of the VMD algorithm. Repeat the loop. For each channel s(x i ,t) performs one-dimensional VMD decomposition in the time dimension to obtain the intrinsic mode (IMF) components of each signal. i ,t)’s i-th VMD decomposition:

[0070] Establish the variational mode decomposition constraint equation,

[0071]

[0072] Where u k is the kth IMF component, ω k for u k The center frequency, K is the mode number, j(j 2 =-1) is the imaginary unit. ∑u k =s(x i ,t) is a constraint condition, which ensures that the sum of all modal components after VMD decomposition is the original signal.

[0073]

[0074] Indicated by u k The complex analytic signal composed of its Hill transform, δ(t) is the Dirac impulse function; * represents convolution. Multiply the complex analytic signal by The term adjusts the spectrum of each modal component to its corresponding baseband and calculates the L2 norm square of the demodulated gradient To estimate each modal component u k bandwidth;

[0075] In order to find the optimal solution of the constrained variational problem, the Lagrangian multiplication operator λ(t) and the second-order penalty factor α are introduced to transform the constrained variational problem (1) into an unconstrained problem. Its Lagrangian expansion expression is:

[0076]

[0077] The Lagrange multiplier λ(t) can ensure the strictness of the constraints, and the second-order penalty factor α can ensure the accuracy of signal reconstruction in Gaussian noise environment;

[0078] The alternating direction multiplier method (ADMM) is used to continuously update the modal components and their center frequencies, and finally the saddle point of the unconstrained model is obtained, which is the optimal solution to the original problem. Combined with the Parseval Fourier isometric transform principle, the n+1th frequency domain ADMM iteration format of the kth order IMF component is:

[0079]

[0080]

[0081]

[0082] Where ω represents the frequency. s(x i ,t), Fourier transform of λ(t). τ is the noise tolerance parameter. Set the convergence error ε and initialize Iterative update according to the iterative format (3) to (5) ω k , When the following iterative convergence conditions are met, the iteration is stopped and the corresponding k-th order IMF component is output ω k ;

[0083]

[0084] Complete the one-dimensional VMD decomposition of all channels and obtain the IMF components of all orders of all channels, which are recorded as imf1;

[0085] c. Calculate the correlation coefficients between each IMF component of imf1 and the original data, and denote the correlation coefficient of the kth (k = 1, 2, ..., K) IMF component as R k , where the largest correlation coefficient is R max . Define R k <R max The IMF component of / 10 is noise, R k >R max The IMF component of / 10 is a valid signal. The IMF component corresponding to the valid signal is screened and the signal is reconstructed to obtain the two-dimensional data s1(x,t);

[0086] d. Perform VMD on the second dimension of the spatial dimension. Similar to step b, initialize the decomposition mode number K' and penalty factor α' of the VMD algorithm. Cycle time sampling points, for each time sampling point s(x,tj ) Perform one-dimensional VMD decomposition in the spatial dimension to obtain the IMF components of each time-space sampling signal, complete the one-dimensional VMD decomposition of all time sampling points, and obtain the IMF components of all time sampling points, which are recorded as imf2;

[0087] e. Calculate the correlation coefficients of each order IMF component of imf2 and the original data respectively, and denote the mutual correlation coefficient of the mth (m=1,2,…,K') IMF component as R m , where the largest correlation coefficient is R' max Similar to step c, define R m <R' max The IMF component of / 10 is noise, R m >R' max The IMF component of / 10 is a valid signal. The IMF component corresponding to the valid signal is screened and the signal is reconstructed to obtain the two-dimensional data s2(x,t);

[0088] f. Use dual-tree complex wavelet transform to perform threshold denoising on the reconstructed data s2(x,t).

[0089] Construct a one-dimensional dual-tree complex wavelet function. Assume that the dual-tree complex wavelet consists of two parallel wavelet trees A and B, which respectively form its real part and imaginary part. The scaling function corresponding to tree A is and the wavelet function ψ h (t), the scaling function corresponding to tree B and the wavelet function ψ g (t). Construct complex wavelet

[0090] ψ(t)=ψ h (t)+jψ g (t) (7)

[0091] In order to make its spectrum unilateral, the wavelet function ψ is required h (t) and ψ g (t) are Hilbert transform pairs, that is, they satisfy

[0092] ψ g (t)=H{ψ h (t)} (8)

[0093] Here, H{·} represents the Hilbert transform.

[0094] Construct a dual-tree complex wavelet filter. Let h0(n) and h1(n) be the low-pass and high-pass filters corresponding to the real part tree A, and g0(n) and g1(n) be the low-pass and high-pass filters corresponding to the imaginary part tree B. There needs to be a half-sample delay between the low-pass filters of tree A and tree B, that is,

[0095] g0(n)=h0(n-0.5) (9)

[0096] Correspondingly, in the frequency domain,

[0097]

[0098] Where G0(ω) and H0(ω) are the Fourier transforms of g0(n) and h0(n), respectively, which are two conjugate orthogonal low-pass filter banks. Equation (10) is the half-frame shift condition, which is equivalent to increasing the sampling of the low-pass signal by a factor of two at each scale.

[0099] Set the number of decomposition levels and perform a two-dimensional dual-tree complex wavelet transform on the reconstructed data s2(x,t) after VMD denoising. The two-dimensional dual-tree complex wavelet transform is obtained by tensor product of one-dimensional dual-tree complex wavelets. The decomposition is equivalent to performing a one-dimensional dual-tree complex wavelet transform on the rows and columns of the two-dimensional data, that is, performing a cascade of one-dimensional dual-tree complex wavelet transforms in the time and space dimensions on s2(x,t).

[0100] Wavelet coefficient threshold processing. Set the threshold T, select the threshold function, and perform wavelet soft threshold denoising on the wavelet coefficients of each layer. Subtract T from the wavelet coefficients whose absolute value is greater than the threshold T, and set the wavelet coefficients whose absolute value is less than the threshold T to zero. The soft threshold function is expressed as

[0101]

[0102] The threshold T is

[0103]

[0104] d u,v represents the vth wavelet coefficient of the uth layer, and V is the number of wavelet coefficients in the uth layer; is the processed wavelet coefficient. The noise standard deviation σ can be estimated using the empirical formula, namely

[0105]

[0106] Data reconstruction: Perform a two-dimensional dual-tree complex wavelet inverse transform on the processed low-frequency wavelet coefficients and high-frequency wavelet coefficients to obtain the denoised data s3(x,t).

[0107] Example 1

[0108] The overall flow chart of the ground penetrating radar signal denoising method combining two-dimensional VMD and DT-CWT is as follows: Figure 1 To better illustrate the implementation effect of the algorithm process, a set of specific examples are given below:

[0109] a. Establishing a layered medium with uniform medium ( Figure 2), its finite difference forward simulation data ( Figure 3 ) as the noise-free data s0(x,t). Establish a layered random equivalent medium model ( Figure 4 ), and its finite difference forward simulation data is used as the original noisy data s(x,t)( Figure 5 );

[0110] b. Perform VMD on the first dimension of time dimension for s(x,t). Initialize the decomposition mode number K and penalty factor α of VMD algorithm. Loop, for each s(x i ,t) Perform one-dimensional VMD decomposition in the time dimension to obtain the intrinsic mode (IMF) components of each signal. After completing the one-dimensional VMD decomposition of all channels, the IMF components of all channels are obtained, which are recorded as imf1;

[0111] c. Calculate the correlation coefficients between each IMF component of imf1 and the original data, and denote the correlation coefficient of the kth (k = 1, 2, ..., K) IMF component as R k , where the largest correlation coefficient is R max . Define R k <R max The IMF component of / 10 is noise, R k >R max / 10 IMF components are valid signals. Filter the IMF components corresponding to the valid signals and reconstruct the signals to obtain the two-dimensional data s1(x,t)( Figure 6 ). According to the signal-to-noise ratio calculation formula The signal-to-noise ratio of the reconstructed data s1(x,t) after the first-dimensional VMD is SNR=9.37;

[0112] d. Perform the second dimension VMD on the spatial dimension of s1(x,t). Similar to step b, initialize the decomposition mode number K' and penalty factor α' of the VMD algorithm. Cycle time sampling point, for each time sampling point s1(x,t j ) Perform one-dimensional VMD decomposition in the spatial dimension to obtain the IMF components of each time-space sampling signal, complete the one-dimensional VMD decomposition of all time sampling points, and obtain the IMF components of all time sampling points, which are recorded as imf2;

[0113] e. Calculate the correlation coefficients of each order IMF component of imf2 and s1(x, t) respectively, and denote the mutual correlation coefficient of the mth (m=1,2,…,K') IMF component as R m , where the largest correlation coefficient is R' max Similar to step c, define R m <R' max The IMF component of / 10 is noise, Rm >R' max / 10 IMF components are valid signals. Filter the IMF components corresponding to the valid signals and reconstruct the signals to obtain the two-dimensional data s2(x,t)( Figure 7 ), its signal-to-noise ratio SNR=12.47;

[0114] f. Use dual-tree complex wavelet transform to perform threshold denoising on the reconstructed data s2(x,t) to obtain the denoised data s3(x,t)( Figure 8 ), its signal-to-noise ratio SNR=14.85.

[0115] In this example, for noisy data ( Figure 5 ) performs conventional one-dimensional VMD decomposition only in the time dimension, and the reconstructed data after denoising ( Figure 6 )The signal-to-noise ratio is 9.37, and then the two-dimensional VMD decomposition in time and space is performed ( Figure 7 ), the signal-to-noise ratio of the reconstructed data after denoising is 12.47. Then, the data after 2D VMD denoising is subjected to 2D dual-tree complex wavelet transform and threshold denoising to reconstruct the data ( Figure 8 ) signal-to-noise ratio is 14.85, further improving the data signal-to-noise ratio. This shows that compared to conventional one-dimensional VMD or single denoising methods, the GPR signal denoising method combining two-dimensional VMD and DT-CWT proposed in this paper has better denoising effects and can produce GPR data with a better signal-to-noise ratio.

[0116] Note that the above are only preferred embodiments of the present invention and the technical principles employed. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described herein, and that various obvious changes, readjustments, and substitutions can be made by those skilled in the art without departing from the scope of protection of the present invention. Therefore, although the present invention has been described in detail through the above embodiments, the present invention is not limited to the above embodiments and may include many other equivalent embodiments without departing from the concept of the present invention. The scope of the present invention is determined by the scope of the appended claims.

Claims

1. A ground penetrating radar signal denoising method combining two-dimensional VMD and DT-CWT, characterized in that: The following steps are involved: a. Obtain the original GPR data and represent the original two-dimensional data with M time sampling points and N as s(x, t), where x represents the horizontal direction and t represents the time direction. The time sampling sequence of the i-th (i=1,…,M) channel can be recorded as s(x i ,t), the spatial sampling sequence at the jth (j=1,…,N) time sample point is recorded as s(x,t j ); b. Perform the first dimension VMD on the time dimension, initialize the decomposition mode number K and penalty factor α of the VMD algorithm, and loop for each s(x i ,t) Perform one-dimensional VMD decomposition in the time dimension to obtain the intrinsic mode (IMF) components of each signal. After completing the one-dimensional VMD decomposition of all channels, the IMF components of all channels are obtained, which are recorded as imf1; c. Calculate the correlation coefficients between each IMF component of imf1 and the original data, and denote the correlation coefficient of the kth (k = 1, 2, ..., K) IMF component as R k , where the largest correlation coefficient is R max , define R k <R max The IMF component of / 10 is noise, R k >R max The IMF component of / 10 is a valid signal. The IMF component corresponding to the valid signal is screened and the signal is reconstructed to obtain the two-dimensional data s1(x, t); d. Perform the second dimension VMD on the spatial dimension. Similar to step b, initialize the decomposition mode number K' and penalty factor α' of the VMD algorithm, cycle the time sampling points, and for each time sampling point s(x,t j ) Perform one-dimensional VMD decomposition in the spatial dimension to obtain the IMF components of each time-space sampling signal, complete the one-dimensional VMD decomposition of all time sampling points, and obtain the IMF components of all time sampling points, which are recorded as imf2; e. Calculate the correlation coefficients of each order IMF component of imf2 and the original data respectively, and denote the mutual correlation coefficient of the mth (m=1,2,…,K') IMF component as R m , where the largest correlation coefficient is R' max , similar to step c, define R m <R' max The IMF component of / 10 is noise, R m >R' max The IMF component of / 10 is a valid signal. The IMF component corresponding to the valid signal is screened and the signal is reconstructed to obtain the two-dimensional data s2(x,t); f. Use dual-tree complex wavelet transform (DT-CWT) to perform threshold denoising on the reconstructed data s2(x, t) to obtain the denoised data s3(x, t).

2. The method for denoising ground penetrating radar signals by combining two-dimensional VMD and DT-CWT according to claim 1, characterized in that: Step b: perform one-dimensional VMD on each channel of the original data s(x,t) in the time dimension, where the i-th channel s(x i ,t) includes the following steps: b1. Establish variational mode decomposition constraint equation Where u k is the kth IMF component, ω k for u k The center frequency, K is the mode number, j(j 2 =-1) is the imaginary unit, ∑u k =s(x i ,t) is a constraint condition, which ensures that the sum of all modal components after VMD decomposition is the original signal; Indicated by u k The complex analytic signal composed of its Hill transform, δ(t) is the Dirac impulse function; * represents convolution, multiplying the complex analytic signal by The term adjusts the spectrum of each modal component to its corresponding baseband and calculates the L2 norm square of the demodulated gradient To estimate each modal component u k bandwidth; b2. To find the optimal solution to the constrained variational problem, the Lagrangian multiplication operator λ(t) and the second-order penalty factor α are introduced to transform the constrained variational problem (1) into an unconstrained problem. Its Lagrangian expansion expression is: The Lagrange multiplier λ(t) can ensure the strictness of the constraints, and the second-order penalty factor α can ensure the accuracy of signal reconstruction in Gaussian noise environment; b3. Use the alternating direction multiplier method (ADMM) to continuously update each modal component and its center frequency, and finally obtain the saddle point of the unconstrained model, which is the optimal solution to the original problem. Combined with the Parseval Fourier isometric transform principle, the n+1th frequency domain ADMM iteration format of the kth order IMF component is: Where ω represents the frequency, s(x i ,t), Fourier transform of λ(t); τ is the noise tolerance parameter, set the convergence error ε, initialize Iterative update according to the iterative format (3) to (5) ω k , When the following iterative convergence conditions are met, the iteration is stopped and the corresponding k-th order IMF component is output ω k ; 3. The method for denoising ground penetrating radar signals by combining two-dimensional VMD and DT-CWT according to claim 1, characterized in that: Step f, further denoising the reconstructed data after 2D VMD denoising using dual dual-tree complex wavelet transform, the specific steps are as follows: f1. Construct a one-dimensional dual-tree complex wavelet function. Assume that the dual-tree complex wavelet consists of two parallel wavelet trees A and B, which respectively constitute the real part and the imaginary part. The scaling function corresponding to tree A is and the wavelet function ψ h (t), the scaling function corresponding to tree B and the wavelet function ψ g (t), construct complex wavelet ψ(t)=ψ h (t)+jψ g (t) (7) In order to make its spectrum unilateral, the wavelet function ψ is required h (t) and ψ g (t) are Hilbert transform pairs, that is, they satisfy ψ g (t)=H{ψ h (t)} (8) Where H{·} represents the Hilbert transform, f2. Construct a dual double-tree complex wavelet filter. Let h0(n) and h1(n) be the low-pass and high-pass filters corresponding to the real part tree A, respectively. Let g0(n) and g1(n) be the low-pass and high-pass filters corresponding to the imaginary part tree B, respectively. There needs to be a half-sample delay between the low-pass filters of tree A and tree B, that is, g0(n)=h0(n-0.5) (9) Correspondingly, in the frequency domain, Where G0(ω) and H0(ω) are the Fourier transforms of g0(n) and h0(n), respectively, which are two conjugate orthogonal low-pass filter banks. Formula (10) is called the half-frame shift condition, which is equivalent to increasing the sampling of the low-pass signal to twice the original sampling at each scale; f3. Set the number of decomposition layers and perform a two-dimensional dual-tree complex wavelet transform on the reconstructed data s2(x, t) after VMD denoising. The two-dimensional dual-tree complex wavelet transform is obtained by the tensor product of the one-dimensional dual-tree complex wavelet. The decomposition is equivalent to performing a one-dimensional dual-tree complex wavelet transform on the rows and columns of the two-dimensional data respectively, that is, performing a one-dimensional dual-tree complex wavelet transform on the time dimension and the space dimension of the cascade of s2(x, t); f4, wavelet coefficient threshold processing, set the threshold T, select the threshold function, and perform wavelet soft threshold denoising on the wavelet coefficients of each layer. Subtract T from the wavelet coefficients whose absolute value is greater than the threshold T, and set the wavelet coefficients whose absolute value is less than the threshold T to zero. The soft threshold function is expressed as: The threshold T is d u,v represents the vth wavelet coefficient of the uth layer, and V is the number of wavelet coefficients in the uth layer; is the processed wavelet coefficient, and the noise standard deviation σ can be estimated using the empirical formula, that is, f5. Data reconstruction: perform two-dimensional dual-tree complex wavelet inverse transform on the processed low-frequency wavelet coefficients and high-frequency wavelet coefficients to obtain the denoised data s3(x,t).

Images