Gpr diffraction and reflection wave separation method using two-dimensional variational mode decomposition

Abstract

The application relates to a GPR diffraction wave and reflection wave separation method using two-dimensional variational modal decomposition, which comprises the following steps: obtaining original underground space ground penetrating radar data for pretreatment, and obtaining pretreated ground penetrating radar data; performing two-dimensional variational modal decomposition on the pretreated data to obtain each order IMF component of the two-dimensional ground penetrating radar data; calculating the frequency-wavenumber spectrum corresponding to each order IMF component; combining the IMF component and the frequency-wavenumber spectrum to analyze the echo type corresponding to each modal component; respectively reconstructing the modal components corresponding to different echo types to obtain ground penetrating radar data after diffraction wave and reflection wave separation; and if the data contains noise, the high-frequency and high-wavenumber modal components corresponding to random noise and random disturbance can be directly removed. According to the distribution characteristic difference of the echo characteristics of different scale target bodies of two-dimensional ground penetrating radar data in the frequency-beam domain, the diffraction wave-reflection wave is separated, and the target geological body is accurately identified and detected.

Description

Technical Field

[0001] This invention belongs to the field of ground-penetrating radar signal processing technology, specifically relating to a data analysis method that uses two-dimensional variational mode decomposition to perform mode decomposition on ground-penetrating radar data and combines frequency-wavenumber spectrum analysis to analyze the wave field type corresponding to each mode. In particular, it relates to a ground-penetrating radar diffraction-reflection wave field separation method based on 2D-VMD. Background Technology

[0002] Ground-penetrating radar (GPR) technology is an important tool in engineering geophysical exploration, widely used for identifying underground targets such as pipelines, cavities, and cracks. In actual exploration, when using broadband recording antennas, two types of effective information are typically collected: reflection information from smooth layers and diffraction information from small-scale, non-uniform, discontinuous geological bodies. However, due to the shielding effect of strong geological reflections, diffraction information is often difficult to extract, thus affecting subsequent interpretation. Therefore, how to separate diffraction and reflection waves from the actual acquired GPR data to achieve the judgment and identification of targets with different characteristics is of significant research importance.

[0003] Based on the acquisition characteristics of ground-penetrating radar (GPR) data, target identification in radar profiles primarily relies on the corresponding diffracted wave echo characteristics. Traditional GPR data processing methods, based on coherence and texture attribute analysis, can improve the resolution of diffracted waves to some extent; however, the shielding effect of strong reflections still causes diffracted waves to become blurred, making identification difficult. Therefore, if diffracted waves can be separated from GPR data to enhance their relevance, small-scale geological anomalies can be effectively identified. Such research is abundant in seismic exploration, typically based on the kinematic and dynamic differences between reflected and diffracted waves to achieve diffracted wave separation. Commonly used diffracted wave separation methods in seismic data processing include: the plane wave destruction filter (PWD) method, which uses the local dip angle difference between adjacent traces of reflected and diffracted waves to suppress reflected waves; the shot gather diffracted wave separation method, combining PWD with spatially variable median filtering algorithms; the Radon transform method based on the kinematic differences between reflected and diffracted waves in common offset traces; and the FK filtering method based on plane wave decomposition.

[0004] Although the electromagnetic waves of ground-penetrating radar (GPR) have similar kinematic propagation laws to the elastic waves of seismic exploration, conventional single-coverage GPR data is often affected by noise. After suppressing strong reflections, the remaining weak diffraction and noise data have a low signal-to-noise ratio. To extract diffraction waves from GPR data to detect small-scale anomalies, Zhaoet.al. proposed a GPR diffraction wave extraction method in 2019 that combines a local plane wave deconstruction filter (PWD) and an online dictionary learning algorithm. This method uses the local plane wave elimination method to estimate and eliminate strong reflections, and a sparse coding algorithm to extract weaker diffraction waves. Chinese patent CN113945987B discloses a method, device, and electronic equipment for detecting geological defects. It uses Seislet transform to shape and regularize GPR wavefield data to obtain separated diffraction wave data, and then performs migration imaging on the separated diffraction wave data. Based on the imaging results, it can accurately detect urban underground geological defects such as goafs and cracks.

[0005] In engineering applications, 2D ground-penetrating radar (GPR) data contains both diffraction wave information from non-uniform and discontinuous geological bodies and reflected wave information from strata. Different targets exhibit different echo characteristics under varying background environments, resulting in varying energy levels of reflected and diffracted waves, often making effective identification of the target geological body difficult. While existing technologies have achieved some degree of separation between reflected and diffracted waves, they primarily focus on suppressing strong reflected waves and highlighting diffracted waves, often resulting in the loss of valuable information and employing complex algorithms with cumbersome steps and limited intuitiveness. If it were possible to effectively separate reflected and diffracted wave pairs without losing any valuable (reflection) information, it would be possible to simultaneously analyze the geological structural information corresponding to both reflected and diffracted waves in GPR data. Summary of the Invention

[0006] The purpose of this invention is to provide a method for separating diffracted and reflected waves from ground-penetrating radar (GPR) data based on two-dimensional variational mode decomposition (2D-VMD). 2D-VMD decomposes the original data into different modal components. The resulting IMF components contain information ranging from low-frequency, low-wavenumber components to high-frequency, high-wavenumber components. Based on the kinematic differences between reflected and diffracted waves, and combined with the frequency-wavenumber spectrum of each modal component, the wavefield type corresponding to each modal component is analyzed. Then, the corresponding modes are reconstructed for different wavefield types, thus obtaining the data after separating the diffracted and reflected waves. This solves the problem of effectively separating diffracted and reflected waves from GPR data, enabling accurate detection of underground anomalies.

[0007] The objective of this invention is achieved through the following technical solution:

[0008] A method for separating GPR diffracted and reflected waves using two-dimensional variational mode decomposition includes the following steps:

[0009] a. Obtain raw underground space ground-penetrating radar data and preprocess it by performing static correction, DC drift removal, and automatic gain control steps to obtain preprocessed ground-penetrating radar data.

[0010] b. Initialize the number of decomposition modes K and the penalty factor α of the 2D-VMD algorithm, and perform two-dimensional variational mode decomposition on the preprocessed data to obtain the intrinsic mode (IMF) components of each order of the two-dimensional ground penetrating radar data.

[0011] c. Calculate the 2D frequency-wavenumber spectrum corresponding to each order of IMF component;

[0012] d. Combine the IMF components in the time-space domain with their frequency-wavenumber spectra to analyze the echo types corresponding to each modal component;

[0013] e. Reconstruct the modal components corresponding to different echo types to obtain ground-penetrating radar data after separation of diffraction and reflection waves; if the data contains noise, the high-frequency and high-wavenumber modal components corresponding to random noise and random disturbances can be directly removed.

[0014] Further, step b involves performing two-dimensional variational mode decomposition on the preprocessed ground-penetrating radar data, specifically including the following steps:

[0015] b1. Define a two-dimensional analytic signal in the frequency domain;

[0016] b2. Establish the two-dimensional variational mode decomposition constraint equations;

[0017] b3. Introduce the quadratic penalty term α and the Lagrange multiplication operator λ to solve the constrained variational problem, transforming the constrained variational problem into an unconstrained problem, and obtain the augmented Lagrange expression;

[0018] b4. The Lagrange multiplier λ is incorporated into the quadratic penalty term, and the alternating multiplier method is used to optimize the calculation of the modal components;

[0019] b5. Perform variational mode decomposition on the original signal to obtain K intrinsic mode components.

[0020] Furthermore, step b1 defines a two-dimensional analytic signal in the frequency domain, specifically as follows:

[0021]

[0022] in, Let ω be the frequency-wavenumber domain representation of the k-th modal component, where ω = (ω1, ω2) represents the frequency-wavenumber domain, {ω k}={ω1,ω2,…,ωk Let} be the center frequency of each modal component, and k (k = 1, 2, ..., K) be the order of the modal component. To obtain the unidirectionality of the analytic signal, one half-plane in the frequency-wavenumber domain is set to zero. The two-dimensional analytic signal with the above Fourier property is:

[0023]

[0024] Where x = (x1, x2) represents the spatiotemporal domain, uk(x) is the spatiotemporal domain representation of the k-th modal component, and u AS k(x) represents the expression by u k The complex analytic signal consisting of (x) and its Hilbert transform, δ( <x,ω k >) represents the Dirac impulse function, * indicates convolution, j(j 2 =-1) is the imaginary unit, and the analytic signal follows ω k The direction calculation is linearly independent.

[0025] Furthermore, in step b2, the two-dimensional variational mode decomposition constraint equations are established, specifically as follows:

[0026]

[0027] For complex analytic signals, multiply by The exponential term is mixed and modulated to the estimated center frequency, and the demodulated gradient L2 norm squared is calculated. To estimate each modal component u k The bandwidth is such that the sum of the estimated bandwidths for each mode is minimized, where αk is the Tikhonov regularization parameter, u k ={u1,u2,…,u k Let f(x) represent the decomposed modal functions, and let f(x) represent the two-dimensional ground-penetrating radar data to be decomposed. The constrained variational problem is solved by introducing a quadratic penalty term α and a Lagrange operator λ.

[0028] Furthermore, in step b3, the resulting augmented Lagrange expression is specifically as follows:

[0029]

[0030] Solving the unconstrained problem becomes:

[0031]

[0032] Furthermore, in step b4, the Lagrange multiplication operator λ is incorporated into the quadratic penalty term, resulting in:

[0033]

[0034] The alternating multiplier method is used to analyze the modal components u. k After optimization calculation, the (n+1)th frequency domain ADMM iteration format of the k-th order IMF component is: (5)

[0035]

[0036]

[0037] in, They represent u respectively k 2D Fourier transforms of f(x) and f(x). This represents the k-th residual of the signal;

[0038] Furthermore, step b5 specifically involves:

[0039] According to formulas (7) and (8), variational mode decomposition is performed on the original signal to obtain K intrinsic mode components. The specific update process of the 2D-VMD algorithm is as follows:

[0040] b51, Initialization

[0041] b52, n = n + 1, are updated alternately in the frequency domain.

[0042] b53. When k = k + 1, repeat until k = K, then update λ.

[0043]

[0044] Where τ is the noise tolerance parameter. The remainder terms are summed in the form of Lagrange multipliers.

[0045] b54. Set the convergence error ε, and repeat formulas (2) and (3) until the following iteration stopping condition is met.

[0046]

[0047] Compared with the prior art, the beneficial effects of the present invention are:

[0048] This invention proposes a ground-penetrating radar (GPR) diffraction-reflection wave separation method based on two-dimensional VMD (Vibration Dynamics). It separates diffraction and reflection waves based on the differences in the frequency-beam domain distribution characteristics of echoes from targets at different scales in two-dimensional GPR data, enabling the identification and detection of target geological bodies. Specifically, it has the following advantages:

[0049] 1. This invention performs 2D-VMD decomposition on the original GPR data and analyzes each modal component and its frequency-wavenumber spectrum by utilizing the difference in distribution characteristics between reflected and diffracted waves in the frequency-wavenumber domain. This allows for more direct determination of the modal components corresponding to diffracted and reflected waves without suppressing effective reflected wave information.

[0050] 2. The modal components after 2D-VMD decomposition contain component information from low frequency and low wavenumber to high frequency and high wavenumber. While separating diffraction waves and reflected waves, it can remove high frequency and high wavenumber modal components corresponding to high frequency random noise and random disturbances in the data, thus achieving the effect of noise reduction. Attached Figure Description

[0051] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0052] Figure 1 The overall flowchart of the GPR diffracted and reflected wave separation method of this invention using two-dimensional variational mode decomposition is shown below;

[0053] Figure 2 Original measured data of underground steel reinforcement pipelines;

[0054] Figure 3 Preprocessed data;

[0055] Figure 4a IMF1 component;

[0056] Figure 4b Frequency-wavenumber spectrum of IMF1;

[0057] Figure 5a IMF2 component;

[0058] Figure 5b Frequency-wavenumber spectrum of IMF2;

[0059] Figure 6a IMF3 component;

[0060] Figure 6b Frequency-wavenumber spectrum of IMF3;

[0061] Figure 7a IMF4 components;

[0062] Figure 7b Frequency-wavenumber spectrum of IMF4;

[0063] Figure 8aIMF5 component;

[0064] Figure 8b Frequency-wavenumber spectrum of IMF5;

[0065] Figure 9 Reconstructing IMF2 and IMF3 yields reconstructed diffraction wave data;

[0066] Figure 10 The reflected wave data are obtained by reconstructing the IMF4 components. Detailed Implementation

[0067] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and not intended to limit it. Furthermore, it should be noted that, for ease of description, the accompanying drawings show only the parts relevant to the present invention, and not all of the structures.

[0068] It should be noted that similar reference numerals and letters in the following figures indicate similar items; therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures. Furthermore, in the description of this invention, terms such as "first," "second," etc., are used only to distinguish descriptions and should not be construed as indicating or implying relative importance.

[0069] This invention relates to a ground-penetrating radar diffraction-reflection wave separation method based on two-dimensional variational mode decomposition, which specifically includes the following steps:

[0070] a. Obtain raw underground space ground-penetrating radar data and preprocess it by performing steps such as static correction and DC drift removal to obtain preprocessed ground-penetrating radar data;

[0071] b. Perform two-dimensional variational mode decomposition (2D-VMD) on the preprocessed data. Initialize the number of decomposition modes K and the penalty factor α for the 2D-VMD algorithm. Define the two-dimensional analytic signal in the frequency domain:

[0072]

[0073] in, Let ω be the frequency-wavenumber domain representation of the k-th modal component, where ω = (ω1, ω2) represents the frequency-wavenumber domain, {ω k}={ω1,ω2,…,ω k Let} represent the center frequency of each modal component, and k (k = 1, 2, ..., K) represent the order of the modal component. To obtain the unidirectionality of the analytic signal, one half-plane in the frequency-wavenumber domain is set to zero. The two-dimensional analytic signal with the above Fourier property is:

[0074]

[0075] Where x = (x1, x2) represents the spacetime domain, u k (x) is the spatiotemporal representation of the k-th modal component, u AS,k (x) represents u k The complex analytic signal consisting of (x) and its Hilbert transform, δ( <x,ω k >) represents the Dirac impulse function, and * denotes convolution. j(j 2 =-1) is the imaginary unit. The analytic signal travels along ω k The direction calculation is linearly independent.

[0076] Establish the two-dimensional variational mode decomposition constraint equations:

[0077]

[0078] For complex analytic signals, multiply by The exponential term is mixed and modulated to the estimated center frequency, and the demodulated gradient L2 norm squared is calculated. To estimate each modal component u k The bandwidth is such that the sum of the estimated bandwidths for each mode is minimized. Where α k u is the Tikhonov regularization parameter. k ={u1,u2,…,u k Let} represent the decomposed modal functions, and f(x) represent the two-dimensional ground-penetrating radar data to be decomposed. The constrained variational problem is solved by introducing a quadratic penalty term α and the Lagrange operator λ.

[0079] By introducing a quadratic penalty term α and the Lagrange operator λ to solve the constrained variational problem, and then transforming the constrained problem into an unconstrained problem, we obtain its augmented Lagrange expression:

[0080]

[0081] Solving the unconstrained problem becomes:

[0082]

[0083] Incorporating the Lagrange multiplication operator λ into the quadratic penalty term, we obtain:

[0084]

[0085] The Alternating Direction Method of Multipliers (ADMM) is used to analyze the modal components u. k For optimization calculations, the (n+1)th frequency domain ADMM iteration format for the k-th order IMF component is:

[0086]

[0087]

[0088] in, They represent u respectively k 2D Fourier transforms of f(x) and f(x). This represents the k-th residual of the signal. According to formulas (7) and (8), variational mode decomposition is performed on the original signal to obtain K intrinsic mode components. The specific update process of the 2D-VMD algorithm is as follows:

[0089] 1. Initialization

[0090] 2. n = n + 1, updated alternately in the frequency domain.

[0091] 3. When k = k + 1, repeat until k = K, then update the value of λ.

[0092]

[0093] Where τ is the noise tolerance parameter. The remainder terms are summed in the form of Lagrange multipliers.

[0094] 4. Set the convergence error ε, and repeat formulas (2) and (3) until the following iteration stopping condition is met.

[0095]

[0096] c. Calculate the 2D frequency-wavenumber spectrum corresponding to each order of IMF component;

[0097] d. Combine the IMF components in the time-space domain with their frequency-wavenumber spectra to analyze the echo types corresponding to each modal component;

[0098] e. Reconstruct the modal components corresponding to different echo types to obtain ground-penetrating radar data after separation of diffraction and reflection waves; if the data contains noise, the high-frequency and high-wavenumber modal components corresponding to random noise and random disturbances can be removed during data reconstruction to improve the signal-to-noise ratio.

[0099] Example 1

[0100] A method for separating GPR diffracted and reflected waves using two-dimensional variational mode decomposition is illustrated in the following flowchart: Figure 1 As shown. To better illustrate the implementation effect of the algorithm, a specific example is given below:

[0101] a. Obtaining Figure 2The original measured data of the underground steel mesh are shown, and the measured data has been preprocessed. The preprocessed data is as follows: Figure 3 As shown;

[0102] b. Perform 2D-VMD decomposition on the preprocessed ground-penetrating radar data to obtain IMF components of each order, such as... Figures 4a-8a As shown;

[0103] c. Calculate the frequency-wavenumber spectrum corresponding to each order of IMF component, such as Figures 4b to 8b ;

[0104] d. Analyze the IMF components of each order and their frequency-wavenumber spectra (Figures 4-8). In Figure 4, IMF1 is the DC component. In Figures 5 and 6, IMF2 and IMF3 correspond to the diffracted wave components generated by the reinforcing bars. In their respective frequency-wavenumber spectra, when the frequency is positive, the wave values ​​of IMF2 and IMF3 are opposite, indicating that the exit angles of the diffracted waves are opposite. In Figure 7, IMF4 corresponds to the reflected wave, and its wave value in the frequency-wavenumber spectrum is concentrated near 0 wavenumber. In Figure 8, the IMF5 component corresponds to high-frequency, high-wavenumber random noise, which can be directly discarded during data reconstruction to improve the signal-to-noise ratio.

[0105] e. Reconstruct IMF2 and IMF3 to obtain the reconstructed diffraction wave data, such as Figure 9 As shown, the reflected wave data is obtained by reconstructing the IMF4 components, as follows: Figure 10 As shown, effective separation of diffracted and reflected waves is achieved without loss of useful information. From the separated diffracted wave data profile, this invention allows for more intuitive analysis of the diffracted wave information generated by the steel mesh, or further, offset imaging to achieve convergence and localization of the diffracted waves, unaffected by interference or shielding from reflected wave energy.

[0106] Note that the above description is merely a preferred embodiment 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 various obvious changes, readjustments, and substitutions can be made 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 which is determined by the scope of the appended claims.

Claims

1. A method for separating GPR diffracted and reflected waves using two-dimensional variational mode decomposition, characterized in that, Includes the following steps: a. Obtain raw underground space ground-penetrating radar data and preprocess it by performing static correction, DC drift removal, and automatic gain control steps to obtain preprocessed ground-penetrating radar data. b. Initialize the number of decomposition modes K and penalty factor of the 2D variational mode decomposition (2D-VMD) algorithm. Two-dimensional variational mode decomposition is performed on the preprocessed data to obtain the intrinsic mode (IMF) components of each order of the two-dimensional ground-penetrating radar data. c. Calculate the 2D frequency-wavenumber spectrum corresponding to each order of IMF component; d. Combine the IMF components in the time-space domain with their frequency-wavenumber spectra to analyze the echo types corresponding to each modal component; e. Reconstruct the modal components corresponding to different echo types to obtain ground-penetrating radar data after separation of diffraction and reflection waves; if the data contains noise, directly remove the high-frequency and high-wavenumber modal components corresponding to random noise and random disturbances.

2. The method for separating GPR diffracted and reflected waves using two-dimensional variational mode decomposition according to claim 1, characterized in that, Step b involves performing two-dimensional variational mode decomposition on the preprocessed ground-penetrating radar data, specifically including the following steps: b1. Define a two-dimensional analytic signal in the frequency domain; b2. Establish the two-dimensional variational mode decomposition constraint equations; b3. Introduce a secondary penalty term. and Lagrange multiplication operators Solve the constrained variational problem, transform the constrained variational problem into an unconstrained problem, and obtain the augmented Lagrange expression; b4. The Lagrange multiplication operator The modal components are optimized by incorporating a secondary penalty term and using the alternating multiplier method. b5. Perform variational mode decomposition on the original signal to obtain K intrinsic mode components.

3. The method for separating GPR diffracted and reflected waves using two-dimensional variational mode decomposition according to claim 2, characterized in that, Step b1, define the two-dimensional analytic signal in the frequency domain, specifically as follows: （1） in, For the first Frequency-wavenumber domain representation of the first modal components, Represents the frequency-wavenumber domain. The center frequencies of each modal component are... The modal component order is... To obtain the unidirectionality of the analytic signal, a half-plane in the frequency-wavenumber domain is set to zero, resulting in a two-dimensional analytic signal with Fourier properties: （2） in, Represents the spacetime domain, Let k be the spatiotemporal representation of the modal component. Indicates by The complex analytic signal formed by its Hilbert transform, For Dirac impulse functions, This represents convolution, where j is the imaginary unit, and the analytic signal follows... The direction calculation is linearly independent.

4. The method for separating GPR diffracted and reflected waves using two-dimensional variational mode decomposition according to claim 3, characterized in that, Step b2, establish the two-dimensional variational mode decomposition constraint equations, specifically: （3） For complex analytic signals, multiply by The exponential term is mixed and modulated to the estimated center frequency, and the demodulated gradient L2 norm squared is calculated. To estimate each modal component The bandwidth is such that the sum of the estimated bandwidths for each mode is minimized, where, For Tikhonov regularization parameters, This represents the two-dimensional ground-penetrating radar data to be decomposed, by introducing a quadratic penalty term. and Lagrange operators Solve the constrained variational problem.

5. The method for separating GPR diffracted and reflected waves using two-dimensional variational mode decomposition according to claim 4, characterized in that, Step b3 yields the following augmented Lagrange expression: （4） Solving the unconstrained problem becomes: （5）。 6. The method for separating GPR diffracted and reflected waves using two-dimensional variational mode decomposition according to claim 5, characterized in that, Step b4, apply the Lagrange multiplication operator Including the secondary penalty term, we get: （6） The alternating multiplier method is used for modal components. Perform optimization calculations, the first The (n+1)th frequency domain ADMM iteration format for the order IMF component is: （7） （8） in, , They represent and 2D Fourier transform; The signal represents the first k One residual.

7. The method for separating GPR diffracted and reflected waves using two-dimensional variational mode decomposition according to claim 6, characterized in that, Step b5, specifically: According to formulas (7) and (8), variational mode decomposition is performed on the original signal to obtain K intrinsic mode components. The specific update process of the 2D-VMD algorithm is as follows: b51, Initialization ; b52 Alternate updates in the frequency domain ; b53, when Repeat until Update （9） in, This is a noise tolerance parameter. The remainder terms are summed in the form of Lagrange multipliers; b54. Set convergence error Repeat steps b52 and b53 until the following iteration stopping condition is met. （10）。

Images