Mine ground penetrating radar clutter suppression method based on improved SSA-VMD-NSST

Abstract

The application discloses a mine ground penetrating radar clutter suppression method based on improved SSA-VMD-NSST, and belongs to the field of signal processing, and comprises the following steps: S1: inputting a pretreated ground penetrating radar spectrum containing clutter, obtaining the average trajectory of the spectrum, optimizing the decomposition number and the penalty factor of VMD decomposition of the average trajectory by using the sparrow search algorithm, obtaining the optimal parameter combination of VMD, and the fitness function is the minimum envelope entropy; S2: decomposing the spectrum by using the optimal parameters, selecting the VMD component with high similarity according to the Pearson correlation coefficient of the modal component and the original spectrum, and recombining the VMD component into a new spectrum; and performing improved NSST operation on the synthesized spectrum according to the threshold value based on the adaptive change of the scale direction and the flexible threshold function, so as to obtain the final denoised spectrum. The application can not only effectively remove the clutter with a specific center frequency, but also suppress the random clutter in the ground penetrating radar spectrum, and improve the denoising performance of the ground penetrating radar spectrum.

Description

Technical Field

[0001] This invention belongs to the field of signal processing, and in particular relates to a method for suppressing clutter from underground ground-penetrating radar in mines based on an improved SSA-VMD-NSST. Background Technology

[0002] Ground-penetrating radar (GPR) detects underground strata by receiving reflected echoes from subsurface targets. Due to its high detection speed and real-time imaging capabilities, it is widely used in archaeology, underground energy exploration, and road and bridge quality inspection. However, in actual detection processes, the underground medium is complex and variable, electromagnetic waves attenuate rapidly, and there is interference from human-made facilities and clutter from the instrument itself or the system. Consequently, the signals received by GPR are inevitably mixed with various types of clutter. Therefore, noise reduction of GPR data is an essential step. Its purpose is to remove clutter, highlight effective feature information, and obtain data with a high signal-to-noise ratio for subsequent data interpretation.

[0003] Denoising GPR data has been extensively studied, with methods primarily involving converting the data into different domains. However, these methods also have their own inherent limitations. In recent years, Variational Mode Decomposition (VMD) has been widely applied in signal processing. In 2014, Dragomiretskiy et al. proposed a multi-component adaptive decomposition method—VMD. Since VMD itself does not exhibit endpoint effects and can adaptively separate two harmonic signals with similar frequencies in a composite signal, this method has been widely used in GPR data denoising. However, the decomposition effect of VMD is affected by the number of modal components and the penalty factor; improper selection of these two factors can lead to a low degree of decomposition.

[0004] Wavelet transform is the mainstream algorithm for GPR image denoising; however, it is rarely used due to its inability to provide anisotropic singularities and optimal approximation. To overcome these shortcomings, researchers have successively proposed ultrawavelet algorithms, such as curvelet, shearlet, and non-subsampled shearlet transform (NSST). These ultrawavelet algorithms have been widely applied to ground-penetrating radar (GPR) image denoising. Currently, many ultrawavelet threshold function studies mainly involve deleting ultrawavelet coefficients smaller than the threshold and retaining those larger than the threshold. Based on this, inverse ultrawavelet transform is used to obtain the denoised GPR image, but the selection method is too simplistic. Therefore, the choice of threshold function is a decisive factor affecting the ultrawavelet denoising effect. In addition, traditional ultrawavelet thresholding mainly utilizes the sparsity of GPR data in the ultrawavelet domain to suppress clutter, while ignoring the distribution characteristics of effective signals and random clutter at different scales and directions in the ultrawavelet domain. This leads to the loss of a large amount of effective signal while removing clutter. Summary of the Invention

[0005] Objective: To address the shortcomings of existing technologies, this invention provides an improved SSA-VMD-NSST method for suppressing clutter in underground ground-penetrating radar (GPR) systems. This invention primarily involves two stages: First, optimizing VMD parameters using a sparrow search algorithm to achieve an optimal combination of mode decomposition numbers and penalty factors. Based on this, Pearson correlation coefficients are used for screening and recombining to form a new spectrum, thereby removing clutter at specific center frequencies from the spectrum. Second, using a scale-direction adaptive threshold based on L2 norm squared constraints and a flexible threshold function, the spectrum obtained in the first stage is subjected to NSST operation to remove random clutter.

[0006] The combined denoising of the two stages effectively addresses the shortcomings of existing ground-penetrating radar (GPR) data denoising methods, improves the signal-to-noise ratio and target detection performance of GPR data, and enhances the practicality and reliability of the algorithm.

[0007] To achieve the above objectives, the present invention employs the following techniques:

[0008] A method for suppressing clutter from underground ground-penetrating radar in mines based on an improved SSA-VMD-NSST method includes the following steps:

[0009] (1) Input the preprocessed ground-penetrating radar image P1 containing clutter, perform grayscale processing, and obtain the average trajectory t of the image;

[0010] (2) Use the sparrow optimization algorithm to optimize two important parameters of VMD decomposition of t: penalty factor and number of decompositions;

[0011] (3) Perform VMD decomposition on P1 based on the obtained optimal parameters to obtain several signal components containing noisy information;

[0012] (4) Based on the Pearson correlation coefficient, several components are compared with the original signal, some components with high similarity are selected, and they are combined into a new spectrum P2.

[0013] (5) Apply NSST operation to P2 using a threshold based on adaptive change in scale direction and a flexible threshold function to complete further denoising work and obtain the final denoised map P3.

[0014] Compared with the prior art, the present invention has the following advantages:

[0015] 1. This invention adaptively determines the VMD decomposition parameters using a sparrow search algorithm to obtain the optimal parameter combination for the VMD algorithm. The VMD signal decomposition method based on sparrow search optimization is applied to the field of ground penetrating radar signal de-interference. This method improves the efficiency of signal decomposition, reduces the error of decomposition and reconstruction, and has certain advantages in signal decomposition.

[0016] 2. To overcome the discontinuity problem of hard thresholds and reduce the bias of soft thresholds in wavelet thresholding, a novel flexible threshold function is proposed. It is not only continuous over its domain R, but also converges quickly to the actual value when the coefficients are greater than the threshold; when the coefficients are less than the threshold, non-zero operations effectively preserve the useful signal, and the included adjustable parameters enhance its adaptability.

[0017] 3. Traditional wavelet threshold definitions do not change with direction and scale, resulting in the loss of much effective signal while denoising. NSST, as a novel multi-scale, multi-directional time-frequency analysis method, possesses optimal sparse representation capabilities, localization features, and direction sensitivity. This invention decomposes multiple sets of noisy ground-penetrating radar signals using NSST, calculates the square of the L2 norm of the NSST domain coefficients at each scale and direction, and rearranges the data. It is found that the square of the normalized L2 norm continuously decreases with changing direction, exhibiting the most significant change. Therefore, a scale- and direction-adaptive threshold calculation method based on L2 norm square constraints is proposed, solving the problem that traditional wavelet thresholds cannot change with scale and direction.

[0018] 4. This invention combines optimized VMD with improved NSST, which can not only effectively remove GPR clutter with a specific center frequency, but also remove random clutter. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in this invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0020] Figure 1 This is the overall flowchart of the present invention;

[0021] Figure 2 This is a detailed flowchart of the SSA-VMD of the present invention;

[0022] Figure 3 This is a flowchart of the improved NSST operation of the present invention;

[0023] Figure 4 This is a cross-sectional view of Embodiment 1 of the present invention;

[0024] Figure 5 This is a diagram showing the effect of NSST three-level decomposition of Embodiment 1 of the present invention;

[0025] Figure 6This is a graph showing the trend of L1 norm, L2 norm, and the square of L2 norm after normalization and rearrangement of NSST coefficients in different directions at different scales in Embodiment 1 of the present invention.

[0026] Figure 7 This is a graph showing different threshold selection functions for this invention;

[0027] Figure 8 The images show the denoising effects of different methods in Embodiment 2 of the present invention, where (a) is wavelet transform (hard threshold), (b) is wavelet transform (soft threshold), (c) is EMD, (d) is CEEMD, (e) is VMD, (f) is NSST (global threshold), (g) is the proposed algorithm, and (h) is the original spectrum. Figure 9 These are the TNR values ​​under different methods in Embodiment 2 of the present invention. Detailed Implementation

[0028] The technical solution of the present invention will be further described below with reference to the accompanying drawings, but it is not limited thereto. Any modifications or equivalent substitutions to the technical solution of the present invention that do not depart from the spirit and scope of the technical solution of the present invention should be covered within the protection scope of the present invention.

[0029] This invention addresses the problems in existing ground-penetrating radar (GPR) denoising techniques, such as the difficulty in selecting VMD parameters, the inability of wavelet transform thresholds to vary with scale and direction, and the overly simplistic threshold selection function. It proposes a two-stage GPR data denoising method. First, a sparrow search algorithm is used to adaptively determine the VMD parameters. Second, a scale- and direction-adaptive threshold and a flexible threshold function are selected to improve the performance of wavelet transform denoising. Finally, an improved SSA-VMD-NSST clutter suppression method for underground ground-penetrating radar in mines is proposed. Figure 1 As shown in the flowchart, the following technical solution is adopted:

[0030] Step 1: Input the ground-penetrating radar image P1 containing clutter after preprocessing (one-dimensional bandpass filtering, AGC, background denoising, etc.), perform grayscale processing so that the range of each pixel in the image is between 0 and 255, and calculate the average trajectory t of the image.

[0031] Step 2 utilizes the sparrow optimization algorithm to optimize two important parameters of VMD decomposition of t: the penalty factor α and the number of modes k; where the number of modes k is the number of IMF sequences of modal components, and the penalty factor α is the data feature corresponding to each IMF sequence of modal components.

[0032] The specific process of VMD decomposition of t in Step 2 includes the following steps:

[0033] The constrained variational model involved in VMD is:

[0034]

[0035] In the formula {u k}={u1,u2,…,u K} and {w k}={w1,w2,…,w K} represents the k-th modal component and its corresponding center frequency, for a total of K modes; δ(t) is the unit impulse function; j represents the imaginary unit; * represents convolution operation; This is a partial derivative operation; f is the target signal.

[0036] To facilitate the solution, a penalty factor α and a Lagrange multiplier λ are introduced to solve the variational constraint problem

[25] . The resulting augmented Lagrange expression is as follows:

[0037]

[0038] The Alternate Direction Method of Multiplers (ADMM) is used to iteratively solve the saddle point of equation (2), and u is iteratively updated in the frequency domain. k w k and λ, as shown below:

[0039] VMD decomposes the signal into K modal components, as follows:

[0040] (1) Initialization λ 1 And n is 0;

[0041] (2)u k and w k Iterative updates are performed using equations (3) and (4) respectively.

[0042]

[0043]

[0044] ∧ represents the Fourier transform; n is the number of iterations;

[0045] (3) Update λ using equation (5).

[0046]

[0047] In the formula, τ is the fidelity coefficient;

[0048] (4) Repeat steps 2 and 3 until the iteration termination condition is met. The termination condition is given by equation (6).

[0049]

[0050] In the formula, ε is the discrimination precision, and ε>0; the superscript n is the number of iteration steps, and the subscript k represents the current mode number;

[0051] (5) Output K modal components.

[0052] The specific process of optimizing VMD parameters based on sparrow search in Step 2 is as follows: Figure 2 As shown, it includes the following steps:

[0053] (1) Determine the number of modes and the numerical range of the penalty factor, and initialize the parameters; whereby a population of n sparrows can be represented by the following formula:

[0054] X = [x1, x2, ..., x] n ] T ,x i =[x i,1 ,x i,2 ...x i,d (7)

[0055] In the formula: d is the dimension of the variable.

[0056] (2) Calculate the fitness of individual sparrows using the envelope entropy function defined by the following formula, and rank the fitness. The magnitude of the envelope entropy is closely related to the sparsity of the IMF component signal. The more periodic impact components related to clutter contained in the IMF component, the greater the sparsity of the signal, and the smaller the envelope entropy value.

[0057]

[0058] Where N is the number of sampling points of the IMF component signal, a(j) is the envelope amplitude signal of the k modal components decomposed by VMD after Hilbert demodulation, and P j It is the probability distribution sequence obtained by normalizing a(j), and the entropy value of the probability distribution sequence is the envelope entropy; E p The smaller the value, the sparser the signal, the less clutter it contains, and the more obvious the periodic clutter distribution in the original signal; conversely, the larger the value, the more sparse the signal. p The larger the value, the weaker the signal sparsity and the more clutter it contains.

[0059] The fitness value f0 of a sparrow is expressed as follows:

[0060]

[0061] (3) Update the positions of the discoverer, follower, and scout using the following formula:

[0062] The formula for updating the discoverer's location is:

[0063]

[0064] In the formula, This represents the position of the i-th sparrow in the j-th dimension at the current iteration number. α is a random number between (0,1), G is the maximum number of iterations, Q is a random number following a normal distribution, L is a d-dimensional matrix of all 1s, and R2 and ST represent the alarm value and the safety threshold, respectively.

[0065] When R2 < ST, it means there are no predators nearby; otherwise, predators are found nearby.

[0066] Followers monitor the discoverer's movement and compete for position. When a follower fails to gain a position, it flies to another location, and its update formula is:

[0067]

[0068] In the formula, x w x represents the worst position in the current population. p This is the best position currently occupied by the discoverer; A is a d-dimensional matrix, where each element is randomly assigned a value of 1 or -1, and A + =A T (AA T ) -1 When i > n / 2, it means that the follower has not obtained food and needs to go to a more distant place to search for food.

[0069] The scout's location update formula is:

[0070]

[0071] In the formula, x b β represents the optimal position of the current population; β represents the step size control parameter; k is a random number between [-1, 1]; ε is a constant used to avoid a denominator of 0; f i f w f g These represent the current fitness value, the worst fitness value globally, and the best fitness value, respectively.

[0072] Where: when f i ≠f g When the signal is 0, it indicates that the sparrow is on the edge of the population, aware of the danger, and needs to move closer to the center of the population to reduce the risk of being preyed upon; conversely, when the signal is 0, it indicates that the sparrow in the center of the population is aware of the danger and needs to escape from its current position.

[0073] (4) Calculate the fitness value and update the sparrow's position.

[0074] (5) Determine whether the maximum number of iterations has been reached. If it is, execute (6); otherwise, repeat steps (2)-(4).

[0075] (6) Output the parameters corresponding to the optimal sparrow position and end the iteration.

[0076] Step 3: Perform VMD decomposition on P1 based on the obtained optimal parameters to obtain several signal components containing noisy information.

[0077] Step 4: Based on the Pearson correlation coefficient, several components are compared with the original signal, and some components with similarity greater than the threshold are selected and combined into a new spectrum P2.

[0078] The Pearson correlation coefficient involved in step 4 is defined as follows:

[0079]

[0080] In the formula, x and y represent the IMF component and the original spectrum, respectively; E(x) represents the expectation of x.

[0081] Pearson correlation coefficient threshold ρ new Defined as follows, if the Pearson correlation coefficient is greater than ρ new If a component is found to be a valid IMF component, it is retained and reconstructed to obtain a new synthetic spectrum; otherwise, it is discarded.

[0082]

[0083] In the formula, ρ max ρ min These are the maximum and minimum values ​​in the Pearson correlation coefficient sequence between these IMF components and the original signal, respectively.

[0084] Step 5 applies an adaptive threshold based on scale direction and a flexible threshold function to P2 to perform NSST operation, completing further denoising work and obtaining the final denoised map P3.

[0085] In Step 5, the process of performing the NSST operation on P2 based on the scale-oriented adaptive threshold and the flexible threshold function is as follows: Figure 3 As shown, the process for NSST operations is as follows:

[0086] (1) At each decomposition level, a non-subsample pyramid filter (NSPF) is used to obtain low-frequency and high-frequency components. Then, NSPF decomposition is performed again, and the low-frequency components that capture singularities in the spectrum are decomposed iteratively. NSPF yields j+1 sub-spectrums of the same size as the source spectrum and j high-frequency spectra, where j is the decomposition level. Then, l-level multi-directional decomposition is performed.

[0087] (2) Each high-frequency component is decomposed into 21 coefficient components of the same size using a shearlet filter box (SFB). These filters are formed by avoiding subsampling to satisfy the shift invariance property.

[0088] (3) After the above two steps, the low-frequency coefficient and the high-frequency coefficient are obtained.

[0089] (4) Retain low-frequency coefficients and perform thresholding based on scale-direction adaptive changes and flexible thresholding functions on high-frequency coefficients. The defined thresholds and threshold selection functions are shown in the following equations:

[0090] The threshold expression is:

[0091]

[0092] In the formula, d p,q This represents the coefficient matrix of the noisy spectrum in the q-th direction at the p-th scale of the NSST, where N is the number of sampling points. It is a constant greater than 0. For about d p,q The square of the 2-norm.

[0093] To characterize the distribution characteristics of effective signals and random clutter at various scales and directions after NSST decomposition of GPR data, we statistically analyzed the NSST coefficients of multiple sets of noisy data after NSST transformation. Example 1 is provided below. Figure 4 The decomposition scale is 3, and the direction matrix is ​​[3 3 4]. The decomposition effect is as follows: Figure 5 There are 8 directions at the first and second scales, and 16 directions at the third scale, for a total of 32 directions.

[0094] We calculated the variation trends of the L1 norm, L2 norm, and the square of the L2 norm after normalized reordering of the NSST coefficients in different scales and directions in Example 1, such as... Figure 6 After normalizing and reordering at each scale, it was found that the square of the L2 norm of the NSST coefficients at each scale continuously decreases with direction, and the trend of the square of the L2 norm is more obvious than that of the L1 and L2 norms. This pattern has been verified in multiple examples, indicating that the square of the L2 norm of the NSST coefficients can more sensitively characterize the distribution variation of the effective signal in different directions. Therefore, this invention chooses the square of the L2 norm to construct the threshold.

[0095] For any fixed scale, the larger the square of the L2 norm of the NSST coefficients in a certain direction, the more effective signal there is in that direction. The smaller the scale, the smaller the threshold value can be used to protect the effective signal in this direction; conversely, a larger threshold value is used to remove random clutter. Simultaneously, as the scale increases, the effective signal decreases, and ln(p+1) gradually increases, requiring a larger threshold value to suppress clutter signals. In summary, the proposed thresholding method effectively protects the effective signal while suppressing random clutter as much as possible.

[0096] The threshold selection function expression is:

[0097]

[0098] In the formula, n is a positive integer; a, β, and γ represent constants with values ​​greater than 0; the value of β is between [0, 0.5]; and the value of γ is between [0, 1].

[0099] Figure 7 A function is selected for different thresholds, where β and γ in the flexible threshold function take values ​​of 0.1 and 0.1, respectively.

[0100] (5) The NSST inverse transform is applied to the low-frequency and high-frequency coefficients after thresholding to obtain the fused spectrum.

[0101] Figure 8 The image shows a comparison of denoising effects under different methods in Example 2. In the preprocessing of Example 2, the bandpass filter bandwidth was 30–300 MHz, the number of AGC control points was 18, and the gain extremum was 50. In SSA initialization, the population size was 10, the maximum number of iterations was 50, the discoverer ratio was 0.7, the warning ratio was 0.2, and the warning value was 0.6. The optimal parameter combination for the penalty factor and the number of decompositions was (1266, 12).

[0102] To utilize objective metrics to measure the performance of different clutter suppression methods, a referenceless image evaluation metric—Target to Noise Ratio (TNR)—is introduced as an evaluation standard, defined as follows:

[0103]

[0104] Where T and N represent the target region and clutter region, respectively; Q T and Q N X(p,q) represents the number of pixels in the target region and the clutter region, respectively; X(p,q) represents the pixel value. Therefore, TNR reflects the contrast between the target and clutter regions. The greater the contrast, the more obvious the denoising effect of the method. Figure 9 The TNR values ​​for Example 2 are shown under different methods.

[0105] In summary, only this invention effectively suppresses deep clutter while simultaneously enhancing the energy of the target reflected wave group, restoring the specific characteristics of the deep waveform signal, and suppressing the appearance of triangular textures. Therefore, this also indirectly explains why the proposed algorithm achieves the highest TNR (Time-to-Noise Ratio) index.

[0106] The above-described embodiments further illustrate the purpose, technical solution, and advantages of the present invention. It should be understood that the above-described embodiments are merely preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made to the present invention within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for suppressing clutter from underground ground-penetrating radar in mines based on an improved SSA-VMD-NSST, characterized in that: The ground-penetrating radar clutter suppression method is performed according to the following steps: (1) Input the preprocessed ground-penetrating radar image P1 containing clutter, perform grayscale processing, and obtain the average trajectory t of the image; (2) The sparrow search algorithm is used to optimize two important parameters of VMD decomposition of t: the penalty factor and the number of decompositions. The specific steps of the algorithm include: Step 1: Determine the number of decompositions and the range of penalty factors, and initialize parameters, including the maximum number of iterations, sparrow population size, discoverer ratio, scout ratio, and warning value; Step 2: Calculate the fitness of individual sparrows in the population using the envelope entropy function defined below, and then rank the fitness scores: (1) Where N is the number of sampling points, It is the envelope signal obtained by demodulating the k modal components of VMD using Hilbert. Through calculation The probability distribution sequence obtained by normalization is called the envelope entropy by calculating the entropy value of the probability distribution sequence. Step 3: Update the positions of the discoverer, follower, and scout using the following formula: The formula for updating the discoverer's location is: (2) In the formula, This represents the information about the j-th dimension position of the i-th sparrow at the current iteration number. It is a random number between (0,1), G is the maximum number of iterations, Q is a random number that follows a normal distribution, L is a d-dimensional matrix of all 1s, and R2 and ST represent the alarm value and the safety threshold, respectively. The follower position update formula is: (3) In the formula, This indicates the worst position of the current population. This is the best position currently occupied by the discoverer; It is a d-dimensional matrix, where each element is randomly assigned the value 1 or -1, and ;when When the time is right, it means that the followers have not yet obtained food and need to go to a more distant place to search for food; The scout's location update formula is: (4) In the formula, This indicates the optimal position of the current population; Indicates the step size control parameter; It is a random number in the range [-1, 1]. It is a constant used to avoid a denominator of 0; These represent the current fitness value, the worst-case fitness value, and the best-case fitness value, respectively. Step 4: Calculate the fitness value and update the sparrow's position; Step 5: Determine if the maximum number of iterations has been reached. If so, proceed to Step 6; otherwise, repeat steps 2-4. Step 6: Output the parameters corresponding to the optimal sparrow position and end the iteration; (3) Perform VMD decomposition on P1 based on the obtained optimal parameters to obtain several signal components containing noisy information; (4) Based on the Pearson correlation coefficient, several components are compared with the original signal, some components with high similarity are selected, and they are combined into a new spectrum P2; (5) Apply NSST operation to P2 using a threshold based on adaptive change in scale direction and a flexible threshold function to complete further denoising work and obtain the final denoised map P3.

2. The method for suppressing clutter from underground ground-penetrating radar in mines based on the improved SSA-VMD-NSST according to claim 1, characterized in that: The Pearson correlation coefficient defined in step 4 is: (5) In the formula, These represent the VMD two-dimensional map components and the original map, respectively. express The expectation.

3. The method for suppressing clutter from underground ground-penetrating radar in mines based on the improved SSA-VMD-NSST according to claim 1, characterized in that: The improved scale orientation adaptive threshold in step 5 is: (6) In the formula, This represents the coefficient matrix of the noisy spectrum in the q-th direction at the p-th scale of the NSST, where N is the number of sampling points. It is a constant greater than 0. For about The square of the 2-norm.

4. The method for suppressing clutter from underground ground-penetrating radar in mines based on the improved SSA-VMD-NSST according to claim 3, characterized in that: The improved flexible threshold function in step 5 is as follows: (7) In the formula, n is a positive integer; a, β, and γ represent constants with values ​​greater than 0.

Images