An improved deconvolution method for suppressing multiple waves of ground penetrating radar

Abstract

The present invention discloses a method for improving deconvolution to suppress multiple waves of ground-penetrating radar, comprising the following steps: 1. Reading ground-penetrating radar data and preprocessing the data, including: processing the radar data by channel difference; normalizing the radar data; and filtering the radar data; and 2. Improving predictive deconvolution, including: locating the specific position of the reflected wave signal using a local peak method; fitting a reflected wave hyperbola based on the located signal position; controlling the intensity factor of the predictive deconvolution using an estimated intensity parameter; suppressing multiple waves in the ground-penetrating radar data using the predictive deconvolution principle, and outputting the ground-penetrating radar data after eliminating the interference of the multiple waves. The improved predictive deconvolution method designed by the present invention automatically fits a hyperbola based on the similarity of the multiple wave shapes and accurately controls the position of the deconvolution to eliminate the multiple waves, eliminating the interference of the multiple waves and retaining the hyperbola corresponding to the weak target as much as possible.

Description

Technical Field

[0001] The invention belongs to the field of intelligent information processing, and in particular is a method for improving deconvolution to suppress multiple waves of ground penetrating radar. Background Art

[0002] Currently, the mainstream methods for suppressing GPR multiple waves can be divided into three categories.

[0003] One type requires advance knowledge of the relevant information of the ground-penetrating radar's transmitted waves. The multiple wave model is inferred by inverting the radar data based on the wave equation or simulating electromagnetic waves. The multiple wave model is then adaptively subtracted from the original radar data to achieve the purpose of suppressing the multiple waves. This is called the predictive subtraction method based on the wave equation (Niu Binhua et al., Progress in Multiple Wave Suppression Technology Based on Wave Equations, 2003; Wu Di et al., Research on Data-Driven Multiple Wave Attenuation Methods, 2008; Xie Songlei, Research on Complex Structure Multiple Wave Suppression Methods Based on Wave Equations, 2013).

[0004] One type of method that reconstructs the full waveforms of primaries and multiples through iterative inversion is called sparse inversion method (VANGROENESTIJN et al., Estimation of primaries by sparse inversion from passiveseismic data, 2010; FHCYPMA et al., Estimating primaries by sparse inversion, a generalized approach, 2013; YANG Xuming et al., Adaptive sparse inversion multiple suppression method, 2020; BAO Peinan et al., Interlayer multiple suppression method based on iterative inversion, 2021).

[0005] Another type of method, called filtering, relies solely on the characteristic differences between primary and multiple waves to filter out multiple interference using filters (Zhang Chuncheng et al., Research on Shallow Ground Penetrating Radar Clutter Suppression and Synthetic Aperture Imaging Based on Hyperbolic Characteristics, 2005; Li Qiang et al., Sea Clutter Suppression Method Based on Radon Transform at High-Swept Capes, 2007; Zhu Hongxiang et al., Crustal Structure Estimation in Sedimentary Basins - Elimination of Receiver Function Multiple Reverberation Using Predictive Deconvolution Methods, 2018; Ma Jitao, Comparative Analysis of Multiple Suppression Algorithms Based on Three-Dimensional Parabolic Radon Transform, 2022).

[0006] The defects of existing methods:

[0007] (1) Since the performance of the deconvolution algorithm depends on the precise setting of algorithm parameters, practical engineering applications have high requirements for parameter selection;

[0008] (2) When the multiple waves overlap with the effective reflection signal waveforms of other target objects, it is very difficult to solve the optimal prediction filter factor. The current method cannot effectively remove the interference of the multiple waves. Summary of the Invention

[0009] This invention addresses the limitations of predictive deconvolution techniques in obtaining optimal predictive filter factors for regular hyperbolic reflection signals. It provides an improved deconvolution method for suppressing multiples in ground-penetrating radar (GPR). This method employs a hyperbola-fitting approach to extract the reflection signal hyperbola and the multiple hyperbola within a specified area. By comparing the fitted hyperbolas, the multiples' positions are precisely located, thereby constraining the deconvolution delay parameter. This method accurately and conveniently obtains the optimal predictive filter factor in the presence of interference, facilitating subsequent suppression of multiples and extracting valuable subsurface information, providing accurate insights for geological interpretation.

[0010] The present invention provides a method for improving deconvolution to suppress ground penetrating radar multiple waves, comprising the following steps:

[0011] Step 1: Read the ground penetrating radar data and preprocess the data. The specific steps include:

[0012] (1.1) Differential processing of radar data;

[0013] (1.2) Normalize radar data;

[0014] (1.3) Filtering and processing radar data;

[0015] Step 2: Improve the prediction deconvolution. The specific steps include:

[0016] (2.1) Locate the specific position of the reflected wave signal using the local peak method;

[0017] (2.2) Fitting the reflected wave hyperbola according to the located signal position;

[0018] (2.3) Use the estimated strength parameter to control the strength factor of the predictive deconvolution;

[0019] (2.4) The multiple waves in the GPR data are suppressed by predictive deconvolution, and the GPR data with the multiple wave interference eliminated is output.

[0020] The beneficial effects of the method of the present invention are:

[0021] (1) Preprocessing the input hyperbola GPR data can significantly eliminate irrelevant information such as possible noise interference, retain and enhance the hyperbola features, and improve the reliability and stability of the fitted hyperbola.

[0022] (2) The improved predictive deconvolution method designed in the present invention automatically fits a hyperbola based on the similarity of the shapes of multiple waves and accurately controls the position of the deconvolution to eliminate the multiple waves, thereby eliminating the interference of multiple waves and retaining the hyperbola corresponding to the weak target as much as possible. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] Figure 1 Flowchart for the improved predictive deconvolution algorithm. DETAILED DESCRIPTION

[0024] The present invention will be further described below with reference to the embodiments and drawings, but the present invention is not limited thereto.

[0025] Example

[0026] A method for suppressing multiple waves of ground penetrating radar by improving deconvolution comprises the following steps:

[0027] Step 1: Read the ground penetrating radar data and preprocess the data. The specific steps include:

[0028] (1.1) Differential processing of radar data;

[0029] Assume that i is the number of GPR data channels, parameter t represents the number of time samples per channel, and d i (t) is the i-th GPR data, and the mathematical expression of the input GPR data D(t) is shown in formula (1):

[0030] D(t)={d1(t),...,d i (t),...,d L (t)},t=1,2,...,T (1)

[0031] D(t) is differentially processed to reduce the irregular fluctuations between data and eliminate possible noise interference as much as possible to make the reflection curve more stable. This embodiment uses the second-order difference idea of ​​image processing to perform differential processing on each data of the ground penetrating radar data to obtain the radar data u after differential processing. i (t), which can better highlight the curve shape of the reflected wave compared to the original data.

[0032] (1.2) Normalize radar data;

[0033] Normalization is a way to simplify calculations, which transforms dimensional expressions into dimensionless ones, making them into scalar quantities with a fixed standard form. This embodiment uses the Min-Max normalization method to transform the ground penetrating radar data u i (t) is scaled to [0,1], and n is obtained. i (t) is the normalized GPR data.

[0034] (1.3) Filtering and processing radar data;

[0035] Data filtering can reduce noise interference and lower the impact of data quality degradation caused by data loss. If a point is significantly different from its surroundings, it is infected by noise and the abrupt data point needs to be removed.

[0036] First, the normalized ground penetrating radar data n i (t) Adopt the neighborhood averaging method to obtain the radar data m after neighborhood averaging processing i (t), t=1,2,...T;

[0037] Next, a filter is created by estimating the local mean and variance around each data point of the GPR data. i (t) Perform filtering to obtain the filtered GPR data of channel i, as shown in formula (2):

[0038]

[0039] W(t)={w1(t),...,w i (t),...w L (t)}, t=1,2,...,T, W(t) is the GPR data after filtering, The Wiener filter is used for filtering, which can adapt to the local variance of the image. When the variance is large, almost no smoothing is performed; when the variance is small, more smoothing is performed.

[0040] Step 2: Improve the prediction deconvolution. The algorithm flow chart is as follows: Figure 1 As shown, the specific steps of the improved method include:

[0041] (2.1) Locate the specific position of the reflected wave signal using the local peak method;

[0042] After preprocessing the data, the specific location of the reflected wave needs to be located in order to constrain the location area of ​​the multiple waves eliminated by deconvolution. The reflected waves are mostly parabolic or hyperbolic, so finding the peak becomes the breakthrough point for locating the object position.

[0043] This embodiment uses the local peak method to locate the exact position of the reflected wave. First, for each data w i (t) Perform first-order difference calculation to obtain g i (t), and then find the possible peak value, see formula (3):

[0044] P i ={g i (t)||g i (t)-β|>3σ 2} (3)

[0045] In formula (3), β and σ 2 g i (t) mean and variance, and finally, find P i The position of the maximum value is used for hyperbola fitting, see formula (4):

[0046] v i ={t|g i (t) = maxP i} (4)

[0047] Through the above algorithm processing, the position l of the reflected wave is obtained i Information, where l i This is key information and requires the use of location l i To fit and optimize the hyperbola expression of the reflected wave, and minimize the error between the fitted hyperbola and the actual reflected wave, let the set of position coordinates l i ={(i,v i )|1≤i≤L,1≤v i ≤T} is the final calculated value.

[0048] (2.2) Fitting the reflected wave hyperbola according to the located signal position;

[0049] Use the local peak method to determine the reflected wave position l i Finally, to further infer the shape of the reflected wave, it is necessary to fit a hyperbola. The fitted hyperbola can infer the position and shape of the multiple reflected waves, so that the interference of the multiple waves can be directly subtracted from the original data in the later stage.

[0050] Let the parameter set of the hyperbola be h = {h1, h2, h3, h4}, where h1 is the fitting angle; h2 is the average value of the range of the offset i of the data to be fitted; h3 is the average value of the offset i of the data to be fitted; h4 is the position w of the peak of the data to be fitted i The average value of

[0051] It is known that the vertical coordinate of the data to be fitted is v i , assuming that the vertical coordinate of the fitting data is y i , according to the least square theorem, the hyperbola fitting can be transformed into a constrained optimization problem, see formula (5):

[0052]

[0053] in, Indicates the optimized parameter value. It can be seen that the fitted hyperbola corresponds to each data position See formula (6):

[0054]

[0055] (2.3) Use the estimated strength parameter to control the strength factor of the predictive deconvolution;

[0056] Assume that the reflection signal area is located from the ζth channel to the εth channel, and the time sampling arrive The multiple wave hyperbola is obtained by scanning D(t). Assuming that the scanning step is λ, that is, the first scanning position is from the ζth track to the εth track, and the time sampling is arrive Each scan can obtain a fitted hyperbola. Due to the periodicity of multiple waves, when multiple waves are successfully obtained, the difference between their shape and the hyperbola fitted by the reflection signal is small. Therefore, by comparing the hyperbola fitted by the reflection signal with the hyperbola fitted by each scan, it can be determined whether the hyperbola obtained by the scan is a multiple wave.

[0057] Assume that the reflected signal fits the hyperbola as: Where i is the number of channels, are the hyperbola parameters for fitting the reflection signal, and the hyperbola for fitting the m-th multiple wave is: Fit hyperbola parameters to the reflection signal;

[0058] First, calculate the fitting degree of the hyperbola corresponding to the reflected signal and the hyperbola corresponding to the multiple wave, as shown in formula (7):

[0059]

[0060] When the fitting degree between the hyperbola corresponding to the reflected signal and the hyperbola corresponding to the multiple wave satisfies the following conditions, see formula (8), it is determined that the multiple wave is successfully located;

[0061] δm<T (8)

[0062] Where T is the judgment threshold. When the fitting degree is less than the judgment threshold, it can be determined that the multiple waves corresponding to the hyperbola are successfully located.

[0063] Then extract the predicted deconvolution prediction step size α for subsequent deconvolution operations, and calculate it as shown in formula (9):

[0064]

[0065] (2.4) Suppressing multiple waves in the GPR data through predictive deconvolution and outputting GPR data after eliminating multiple wave interference;

[0066] The core problem of predictive deconvolution is to design the deconvolution factor s(t), as shown in formula (10):

[0067]

[0068] Where α is the deconvolution prediction step size, ρ is the prediction filter length, and c(t) = [c(0), c(1), ..., c(ρ)] is the prediction filter factor, which can be obtained based on the least squares theorem;

[0069] In order to avoid the instability caused by the amplitude of the wavelet amplitude spectrum being zero or close to zero at a certain frequency, it is necessary to pre-whiten the wavelet amplitude spectrum during the solution process;

[0070] Finally, the deconvolution factor s(t) is combined with the i-th channel data d i (t) Convolution, we can get the i-th channel data q after removing the multiple wave interference i (t), see formula (11):

[0071]

[0072] Through the above steps, we finally achieved the goal of completely suppressing the multiple waves of strong targets and making the primary wave signals of weak targets clearer, which is more convenient for subsequent geological research and analysis, and also provides some research ideas and methods worthy of reference for handling related work.

Claims

1. A method for improving deconvolution to suppress ground penetrating radar multiple waves, characterized in that: The following steps are involved: Step 1: Read the ground penetrating radar data and preprocess the data. The specific steps include: (1.1) Differential processing of radar data; Assume that i is the number of GPR data channels, parameter t represents the number of time samples per channel, and d i (t) is the i-th GPR data, and the mathematical expression of the input GPR data D(t) is shown in formula (1): D(t)={d1(t),...,d i (t),...,d L (t)},t=1,2,...,T (1) Perform differential processing on D(t), and use the second-order differential idea of ​​image processing to perform differential processing on each data of the ground penetrating radar data to obtain the ground penetrating radar data u after differential processing. i (t); (1.2) Normalize radar data; The Min-Max normalization method is used to convert the ground penetrating radar data u i (t) is scaled to [0,1], and n is obtained. i (t) is the normalized ground penetrating radar data; (1.3) Filtering and processing radar data; First, the normalized ground penetrating radar data n i (t) Adopt the neighborhood averaging method to obtain the radar data m after neighborhood averaging processing i (t), t=1,2,...T; Next, a filter is created by estimating the local mean and variance around each data point of the GPR data. i (t) Perform filtering to obtain the filtered GPR data of channel i, as shown in formula (2): W(t)={w1(t),...,w i (t),...w L (t)}, t=1,2,...,T, W(t) is the GPR data after filtering, The Wiener filter is used for filtering, which can adapt to the local variance of the image. When the variance is large, almost no smoothing is performed; when the variance is small, more smoothing is performed. Step 2: Improve the prediction deconvolution. The specific steps include: (2.1) Locate the specific position of the reflected wave signal using the local peak method; First, for each data w i (t) Perform first-order difference calculation to obtain g i (t), and then find the possible peak value, see formula (3): P i ={g i (t)||g i (t)-β|>3σ 2 } (3) In formula (3), β and σ 2 g i The mean and variance of (t); Finally, find P i The position of the maximum value is used for hyperbola fitting, see formula (4): v i ({t|g i (t)6maxP i } (4) Through the above algorithm processing, the position l of the reflected wave is obtained i Information, where l i This is key information and requires the use of location l i To fit and optimize the hyperbola expression of the reflected wave, and minimize the error between the fitted hyperbola and the actual reflected wave, let the set of position coordinates l i ={(i,v i )|1≤i≤L,1≤v i ≤T} is the final calculated value; (2.2) Fitting the reflected wave hyperbola according to the located signal position; Let the parameter set of the hyperbola be h = {h1, h2, h3, h4}, where h1 is the fitting angle; h2 is the average value of the range of the offset i of the data to be fitted; h3 is the average value of the offset i of the data to be fitted; h4 is the position w of the peak of the data to be fitted i The average value of It is known that the vertical coordinate of the data to be fitted is v i , assuming that the vertical coordinate of the fitting data is y i , according to the least square theorem, the hyperbola fitting can be transformed into a constrained optimization problem, see formula (5): in, Indicates the optimized parameter value. It can be seen that the fitted hyperbola corresponds to each data position See formula (6): (2.3) Use the estimated strength parameter to control the strength factor of the predictive deconvolution; Assume that the reflection signal area is located from the ζth channel to the εth channel, and the time sampling arrive The multiple wave hyperbola is obtained by scanning D(t). Assuming that the scanning step is λ, that is, the first scanning position is from the ζth track to the εth track, and the time sampling is arrive Each scan can obtain a fitted hyperbola. Due to the periodicity of multiple waves, when multiple waves are successfully obtained, the difference between their shape and the hyperbola fitted by the reflection signal is small. Therefore, by comparing the hyperbola fitted by the reflection signal with the hyperbola fitted by each scan, it can be determined whether the hyperbola obtained by the scan is a multiple wave. Assume that the reflected signal fits the hyperbola as: Where i is the number of channels, are the hyperbola parameters for fitting the reflection signal, and the hyperbola for fitting the m-th multiple wave is: Fit hyperbola parameters to the reflection signal; First, calculate the fitting degree of the hyperbola corresponding to the reflected signal and the hyperbola corresponding to the multiple wave, as shown in formula (7): When the fitting degree between the hyperbola corresponding to the reflected signal and the hyperbola corresponding to the multiple wave satisfies the following conditions, see formula (8), it is determined that the multiple wave is successfully located; d m <T (8) Where T is the judgment threshold. When the fitting degree is less than the judgment threshold, it can be determined that the multiple waves corresponding to the hyperbola are successfully located. Then, the prediction deconvolution step size α is extracted for the subsequent deconvolution operation, and the calculation is shown in formula (9): (2.4) Suppressing multiple waves in the GPR data through predictive deconvolution and outputting GPR data after eliminating multiple wave interference; The core problem of predictive deconvolution is to design the deconvolution factor s(t), as shown in formula (10): Where α is the deconvolution prediction step size, ρ is the prediction filter length, and c(t) = [c(0), c(1), ..., c(ρ)] is the prediction filter factor, which can be obtained based on the least squares theorem; Finally, the deconvolution factor s(t) is combined with the i-th channel data d i (t) Convolution, we can get the i-th channel data q after removing the multiple wave interference i (t), see formula (11): Ultimately, the goal is achieved that the multiple waves of strong targets can be completely suppressed and the primary wave signals of weak targets are clearer.

Images