An Improved FK Migration Method for Ground Penetrating Radar Signals Based on Phase Information

Abstract

This invention relates to signal migration processing methods, and discloses an improved FK migration method for ground-penetrating radar (GPR) signals based on phase information. The method includes the following steps: preprocessing the acquired complex domain echo data to decompose it into amplitude and phase information; segmenting the amplitude information into blocks according to the time-space dimension and fitting an amplitude attenuation model; calculating the phase change rate along the time and space dimensions to obtain the phase fluctuation rate variance; correcting the wavefield propagation path based on real-time electromagnetic wave velocity; performing a two-dimensional Fourier transform on the preprocessed echo signal to convert it to the frequency-wavenumber domain, and correcting it to obtain a new frequency-wavenumber domain echo signal; selecting different interpolation methods for completion; and performing a two-dimensional inverse Fourier transform to obtain the GPR migration imaging result. The method disclosed in this invention achieves a progressive improvement in noise suppression, signal enhancement, and imaging optimization, significantly enhancing the geological information interpretation capability of GPR data.

Description

Technical Field

[0001] This invention relates to signal offset processing methods, and particularly to an improved FK offset method for ground penetrating radar signals based on phase information. Background Technology

[0002] Underground pipeline networks are the core of infrastructure such as urban water supply, gas transmission, and oilfield pipeline management, and their accurate detection relies on ground-penetrating radar (GPR) technology. GPR obtains detailed information about underground structures by emitting electromagnetic waves and receiving signals reflected back from underground targets. Considering the deviations caused by uneven ground in the detection environment during underground pipeline detection, migration algorithms are needed to make the boundaries and details of underground geological bodies clearer. This can effectively handle wavefield imaging of complex geological structures such as folds, faults, and multi-layered media, improving imaging resolution. Common GPR migration algorithms include frequency-wavenumber domain migration (FK migration), Kirchhoff integral migration, and reverse time migration (RTM). Among them, FK migration, with its low computational complexity, adaptability to large-scale data, and strong anti-interference ability, can meet the engineering applications of underground pipeline detection.

[0003] Traditional FK migration algorithms, a key technology in ground-penetrating radar (GPR) imaging, convert time-space echo data to the frequency-wavenumber domain, perform wavefield extrapolation at a fixed wave velocity, select appropriate interpolation methods for calculation, and finally convert the data back from the frequency-wavenumber domain to the spatial domain to obtain the migrated result. This enables the repositioning of tilted reflective interfaces and is a core method for locating underground targets. However, it has significant drawbacks in practical applications: electromagnetic waves attenuate exponentially in underground media, and the echo amplitude of deep pipe networks is greatly reduced, failing to retain target signal characteristics. This results in blurred outlines of deep pipe networks and low imaging accuracy for deep targets. Underground media are usually layered, and changes in water content lead to significant differences in wave velocity. Using a fixed average velocity for FK migration can produce incorrect repositioning in areas with drastic velocity changes, leading to structural distortion or residual diffraction waves. Summary of the Invention

[0004] To address the aforementioned technical problems, this invention provides an improved FK migration method for ground-penetrating radar signals based on phase information. By processing amplitude and phase information, dynamically adjusting interpolation weights, and combining a single-channel three-dimensional imaging process, high-precision detection of underground pipe networks can be achieved.

[0005] To achieve the above objectives, the technical solution of the present invention is as follows:

[0006] An improved FK migration method for ground-penetrating radar signals based on phase information includes the following steps:

[0007] Step 1: Preprocess the complex domain echo data collected by the single-channel sweep frequency ground penetrating radar along the preset survey line, and then decompose the preprocessed echo signal to obtain amplitude information and phase information.

[0008] Step 2: Divide the amplitude information into blocks according to the time-space dimension, calculate the average amplitude of each amplitude block and the average amplitude difference between the current amplitude block and the adjacent amplitude blocks, and take the maximum amplitude difference to obtain the amplitude attenuation rate; combine the physical law of electromagnetic wave exponential decay to obtain the amplitude attenuation model, and set the attenuation threshold based on the maximum value of the amplitude attenuation model.

[0009] Step 3: Calculate the phase change rate along the time and space dimensions of the phase information and take the mean to obtain the variance of the phase fluctuation rate;

[0010] Step 4: Combining the amplitude decay rate and the variance of the phase fluctuation rate, the basic dielectric constant of the medium and the correction coefficient are introduced to obtain the real-time electromagnetic wave velocity; the wave field propagation path is corrected based on the real-time electromagnetic wave velocity to obtain the corrected time wavenumber.

[0011] Step 5: Perform a two-dimensional Fourier transform on the preprocessed echo signal in the horizontal-time dimension, and convert it to the frequency-wavenumber domain to obtain the echo signal in the frequency-wavenumber domain. Then, based on the corrected time wavenumber, obtain a new echo signal in the frequency-wavenumber domain.

[0012] Step 6: Select different interpolation methods to complete the echo signal in the new frequency-wavenumber domain based on the relationship between the maximum amplitude difference and the attenuation threshold. Finally, perform a two-dimensional inverse Fourier transform on the completed frequency-wavenumber domain echo signal to obtain the ground-penetrating radar migration imaging result.

[0013] In the above scheme, step 1, the preprocessing includes DC removal, background elimination, mean filtering, and time gain compensation.

[0014] In the above scheme, in step 1, the preprocessed echo signal Decomposed into amplitude information With phase information The formula is:

[0015] ;

[0016] ;

[0017] in, For the real part of the complex number, It represents the imaginary part of a complex number.

[0018] In the above scheme, the specific method of step 2 is as follows:

[0019] Amplitude information Divided according to the time-space dimension into Block amplitude block, In the time and space domains, the first... The first time block, the first Amplitude blocks within a spatial block; in the time dimension For time block indexes, each time block contains Each sampling point; in the spatial dimension, For spatial block indexes, each spatial block contains One sampling point;

[0020] For each block Calculate the average amplitude within the block. The formula is , It is used to mark the time direction in small blocks. The first, spatial direction The index of each sampling point;

[0021] Calculate the difference between the average amplitude of the current amplitude block and the average amplitude of its adjacent amplitude blocks to the right and below, and take the maximum amplitude difference. ; Utilizing the maximum amplitude difference The amplitude decay rate was calculated. :

[0022] ;

[0023] in, For time step, For spatial step size, The standard spatial attenuation coefficient, The standard time decay factor; The spatial attenuation coefficient, This is the time decay coefficient;

[0024] Based on the amplitude of the electromagnetic wave at the initial position Propagation time of electromagnetic waves Depth of the target Spatial attenuation coefficient Time decay coefficient The amplitude attenuation model is obtained. ;

[0025] The attenuation threshold is obtained based on the maximum value of the attenuation model. , , This represents the adjustment coefficient.

[0026] In the above scheme, the specific method for step 3 is as follows:

[0027] For phase information, in the time dimension, the time step is calculated as follows: phase change rate :

[0028] ,

[0029] in, Indicates time step Later, in time ,space Phase value at; Indicates the current time ,space Phase value at;

[0030] In the spatial dimension, the computational spatial step size is... phase change rate :

[0031] ,

[0032] in, Represents spatial step size Later, in time ,space Phase value at;

[0033] Taking the average of the two yields the average rate of phase change in both time and space dimensions. This reflects the overall fluctuation of phase over time and space:

[0034] ,

[0035] Then the variance of the phase fluctuation rate is calculated. :

[0036] ;

[0037] in, This represents the total number of time samples. Indicates the index of the time sampling point. Indicates the spatial sampling point index. This represents the total number of spatial samples. for The global mean.

[0038] In the above scheme, in step 4, the real-time electromagnetic wave velocity The calculation formula is as follows:

[0039] ;

[0040] in, At the speed of light, , For correction factor, The fundamental dielectric constant of the dielectric. The amplitude decay rate, The variance of the phase fluctuation rate, is the relative permittivity that varies with time and space.

[0041] In the above scheme, the corrected time wavenumber calculation method in step 4 is as follows:

[0042] ;

[0043] in, This is the corrected time wavenumber; For real-time electromagnetic wave velocity, At the speed of light, For space wavenumber, For reference wavenumber.

[0044] In the above scheme, step 5, the method for performing a two-dimensional Fourier transform of the preprocessed echo signal in the horizontal-time dimension to convert it to the frequency-wavenumber domain is as follows:

[0045] ;

[0046] in, For space wavenumber, For time wavenumber, This is the preprocessed echo signal; Indicates the index of the time sampling point. Indicates the spatial sampling point index. This represents the echo signal in the frequency-wavenumber domain obtained after the two-dimensional Fourier transform;

[0047] Based on the corrected time wavenumber Obtain the echo signal in the new frequency-wavenumber domain .

[0048] In the above scheme, in step 6, based on the maximum amplitude difference With decay threshold The following are methods for completing the echo signal in the new frequency-wavenumber domain by selecting different interpolation methods based on the relationship:

[0049] when When using linear interpolation, the following is obtained: ;

[0050] in, This represents the result of linear interpolation, which is the echo signal in the frequency-wavenumber domain after linear interpolation completion. Represents the echo signal in the frequency-wavenumber domain The Middle Average amplitude data of each data block Represents the echo signal in the frequency-wavenumber domain Zhongyu The adjacent first Average amplitude data of each data block The relation coefficients for the current data block;

[0051] when When using cubic spline interpolation, let the given... Echo signal in the frequency-wavenumber domain interpolation nodes The cubic spline interpolation function is as follows:

[0052] ,

[0053] in, This represents the result of cubic spline interpolation, which is the echo signal in the frequency-wavenumber domain after the cubic spline interpolation is completed. For the segmented interval index of cubic spline interpolation, , Representing the spatial-temporal dimension variables in the frequency-wavenumber domain, This represents the x-coordinate of the node in the cubic spline interpolation. , , , Represents the coefficients of a cubic spline function;

[0054] Will and Unified representation of the echo signal in the completed frequency-wavenumber domain .

[0055] In the above scheme, step 6, the method for performing a two-dimensional inverse Fourier transform on the completed frequency-wavenumber domain echo signal is as follows:

[0056] ;

[0057] in, This indicates the final ground-penetrating radar offset imaging result. This represents the echo signal in the completed frequency-wavenumber domain. For space wavenumber, The corrected time wavenumber, This represents the two-dimensional inverse Fourier transform.

[0058] Through the above technical solution, the improved FK migration method for ground-penetrating radar signals based on phase information provided by the present invention has the following beneficial effects:

[0059] 1. In the preprocessing of this invention, the DC removal operation can accurately eliminate DC component interference introduced by instruments or the environment in the original data, making the signal baseline return to stability; the background elimination algorithm, based on the statistical characteristics of the signal, can identify and remove background clutter unrelated to underground targets, reducing the masking effect of noise on the effective signal, and making the target reflection initially prominent; the mean filtering targets false reflections caused by multipath propagation in ground penetrating radar, and cuts off the propagation path of such interference through a wavefield separation strategy to avoid misleading target interpretation; the time gain compensation, based on the signal amplitude distribution, intelligently strengthens weak reflection signals and moderately suppresses strong interference, balancing the signal performance of targets at different depths and with different materials, so that the reflection information of each layer is presented more evenly. After this series of preprocessing steps, the noise floor of the original data is significantly reduced.

[0060] 2. This invention adjusts the interpolation strategy based on the attenuation characteristics of phase and amplitude information. It fits an amplitude attenuation model by combining the physical law of exponential attenuation of electromagnetic waves and sets an attenuation threshold. The interpolation method is selected based on the relationship between the maximum amplitude difference and the attenuation threshold to better preserve and recover useful information in the data. Simultaneously, the interpolation logic is adjusted based on the attenuation characteristics of phase information, relying more on phase interpolation in regions with large amplitude variations. The amplitude and phase information of each frequency component are obtained. In the FK transform space, an attenuation threshold can be set to delineate the interpolation region based on the amplitude attenuation characteristics of the signal. In regions with small amplitude attenuation, conventional amplitude plus phase interpolation is used. In regions with large amplitude attenuation, the influence of amplitude interpolation is reduced, and phase interpolation is relied upon more heavily to avoid interpolation distortion of deep target signals due to amplitude attenuation.

[0061] 3. The improved FK offset of this invention calculates the real-time electromagnetic wave velocity by combining the amplitude attenuation rate and the variance of the phase fluctuation rate, and then corrects the wavefield propagation path to obtain the corrected time wavenumber. By optimizing the wavenumber domain filtering strategy and phase correction model, it addresses the signal dispersion and boundary blurring problems remaining after basic processing, achieving more accurate wavefield extension and imaging focusing. The precise focusing depth fusion of preprocessing and improved FK offset achieves a progressive improvement in noise suppression, signal enhancement, and imaging optimization, significantly enhancing the geological information interpretation capability of ground-penetrating radar data. Attached Figure Description

[0062] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below.

[0063] Figure 1 This is a schematic diagram of an improved FK offset method for ground-penetrating radar signals based on phase information disclosed in an embodiment of the present invention;

[0064] Figure 2 It consists of a single-channel swept-frequency ground-penetrating radar system;

[0065] Figure 3 Images before and after preprocessing; (a) before preprocessing, (b) after preprocessing;

[0066] Figure 4 The figure shows a comparison between the basic FK offset and the improved method. (a) shows the result of the traditional FK offset, and (b) shows the result of the improved FK offset of this invention.

[0067] In the diagram, 1. Detection vehicle; 2. Frequency-sweeping ground-penetrating radar board; 3. Single-channel antenna; 4. External trigger wheel; 5. Main unit software; 6. Pipeline marked area. Detailed Implementation

[0068] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention.

[0069] This invention provides an improved FK migration method for ground-penetrating radar signals based on phase information, such as... Figure 1 As shown, it includes the following steps:

[0070] Step 1: Preprocess the complex domain echo data collected by the single-channel sweep frequency ground penetrating radar along the preset survey line, and then decompose the preprocessed echo signal to obtain amplitude information and phase information.

[0071] like Figure 2 As shown, the single-channel swept-frequency ground-penetrating radar includes a detection trolley 1. The bottom of the trolley 1 is equipped with a swept-frequency ground-penetrating radar plate 2 and a single-channel antenna 3. The front wheels are external trigger wheels 4, which move along a preset survey line within the marked area 6 of the pipeline. Data acquisition and control are performed by the host software 5. The external trigger wheels 4 record the distance traveled by the detection trolley 1, and the host software 5 simultaneously acquires complex domain echo signals. ( This represents the index of the time sampling point, in nanoseconds (ns). The index of the spatial sampling point (in meters) is used to sequentially perform DC removal, background elimination, mean filtering, and time gain compensation on the original echo signal to eliminate noise and interference.

[0072] 1. DC removal

[0073] Eliminating the DC component in the signal improves the performance and measurement accuracy of the radar system. The formula is as follows: ,in This represents the number of time sampling points; This is the echo signal after DC removal.

[0074] 2. Background Removal

[0075] To remove static noise and background interference from radar data in order to clearly identify underground targets, the mean noise value of target-free areas is extracted. The formula is ; The echo signal after background removal;

[0076] 3. Mean filtering

[0077] for For each sampling point, a neighborhood window is selected within a certain range. The average value of all sampling points within this window is calculated, and this average value replaces the original sampling point's value. This method effectively suppresses high-frequency noise and low-frequency interference in the data, making the target signal more prominent and clear, resulting in a filtered echo signal. .

[0078] 4. Time gain compensation

[0079] Time-gain compensated radar (TGC) is used to adjust the gain of the data based on the propagation characteristics of radar waves in underground media, thereby compensating for energy attenuation and improving the signal-to-noise ratio. The TGC slope is adjusted. (dB / ns), start time and termination time For each sampling point, determine whether it falls within the TGC's effective range based on its corresponding time. If it does, proceed according to the formula... Calculate dB gain Then use the formula Convert it to a linear scaling factor The amplitude of the output signal is calculated by directly multiplying the two values; if the amplitude is outside the range, the amplitude of the signal at that sampling point is kept constant to obtain the amplified echo signal. .

[0080] The preprocessed echo signal Decomposed into amplitude information With phase information The formula is:

[0081] ;

[0082] ;

[0083] in, For the real part of the complex number, It represents the imaginary part of a complex number.

[0084] Step 2: Divide the amplitude information into blocks according to the time-space dimension, calculate the average amplitude of each amplitude block and the average amplitude difference between the current amplitude block and the adjacent amplitude blocks, and take the maximum amplitude difference to obtain the amplitude attenuation rate; combine the physical law of electromagnetic wave exponential decay to obtain the amplitude attenuation model, and set the attenuation threshold based on the maximum value of the amplitude attenuation model.

[0085] The specific method is as follows:

[0086] Amplitude information Divided according to the time-space dimension into Block amplitude block, In the time and space domains, the first... The first time block, the first Amplitude blocks within a spatial block; in the time dimension For time block indexes, each time block contains Each sampling point; in the spatial dimension, For spatial block indexes, each spatial block contains One sampling point;

[0087] For each block Calculate the average amplitude within the block. The formula is , It is used to mark the time direction in small blocks. The first, spatial direction The index of each sampling point;

[0088] Calculate the difference between the average amplitude of the current amplitude block and the average amplitude of its adjacent amplitude blocks to the right and below, and take the maximum amplitude difference. ; Utilizing the maximum amplitude difference The amplitude decay rate was calculated. :

[0089] ;

[0090] in, For time step, For spatial step size, The standard spatial attenuation coefficient, The standard time decay factor; The spatial attenuation coefficient, This is the time decay coefficient;

[0091] Based on the amplitude of the electromagnetic wave at the initial position Propagation time of electromagnetic waves Depth of the target Spatial attenuation coefficient Time decay coefficient The amplitude attenuation model is obtained. ;

[0092] The attenuation threshold is obtained based on the maximum value of the attenuation model. , , This represents the adjustment coefficient.

[0093] Step 3: Calculate the phase change rate along the time and space dimensions of the phase information and take the mean to obtain the variance of the phase fluctuation rate.

[0094] Calculating the phase difference along both the time and space dimensions can reflect the degree of phase fluctuation. The specific method is as follows:

[0095] For phase information, in the time dimension, the time step is calculated as follows: phase change rate :

[0096] ,

[0097] in, Indicates time step Later, in time ,space Phase value at; Indicates the current time ,space Phase value at; time step Spatial step size Adjustments were made based on the distribution of media and pipelines in the measured data.

[0098] In the spatial dimension, the computational spatial step size is... phase change rate :

[0099] ,

[0100] in, Represents spatial step size Later, in time ,space Phase value at;

[0101] Taking the average of the two yields the average rate of phase change in both time and space dimensions. This reflects the overall fluctuation of phase over time and space:

[0102] ,

[0103] Then the variance of the phase fluctuation rate is calculated. :

[0104] ;

[0105] in, This represents the total number of time samples. Indicates the index of the time sampling point. Indicates the spatial sampling point index. This represents the total number of spatial samples. for The global mean.

[0106] Step 4: Combining the amplitude decay rate and the variance of the phase fluctuation rate, the basic dielectric constant of the medium and the correction coefficient are introduced to obtain the real-time electromagnetic wave velocity; the wave field propagation path is corrected based on the real-time electromagnetic wave velocity to obtain the corrected time wavenumber.

[0107] Because the underground medium is non-uniform, wave velocity varies. Wave velocity obtained by combining phase information more closely reflects the actual underground conditions. Real-time electromagnetic wave velocity. The calculation formula is as follows:

[0108] ;

[0109] in, At the speed of light, , For correction factor, The fundamental dielectric constant of the dielectric. The amplitude decay rate, The variance of the phase fluctuation rate, is the relative permittivity that varies with time and space.

[0110] The corrected method for calculating the time wavenumber is as follows:

[0111] ;

[0112] in, The corrected time wavenumber is a parameter that describes the frequency characteristics of the wave field in the time dimension in the frequency-wavenumber domain. After real-time electromagnetic wave velocity correction, it can accurately reflect the actual propagation law of the wave field in the underground non-uniform medium. For real-time electromagnetic wave velocity, At the speed of light, For space wavenumber, For reference wavenumber.

[0113] Step 5: Perform a two-dimensional Fourier transform on the preprocessed echo signal in the horizontal-time dimension to convert it to the frequency-wavenumber domain, and then obtain the echo signal in the frequency-wavenumber domain based on the corrected time-wavenumber.

[0114] Preprocessed echo signal The method for performing a two-dimensional Fourier transform in the horizontal-time dimension to convert to the frequency-wavenumber domain is as follows:

[0115] ;

[0116] in, For space wavenumber, For time wavenumber, This is the preprocessed echo signal; Indicates the index of the time sampling point. Indicates the spatial sampling point index. This represents the echo signal in the frequency-wavenumber domain obtained after the two-dimensional Fourier transform;

[0117] Based on the corrected time wavenumber Obtain the echo signal in the new frequency-wavenumber domain .

[0118] Step 6: Select different interpolation methods to complete the echo signal in the new frequency-wavenumber domain based on the relationship between the maximum amplitude difference and the attenuation threshold. Finally, perform a two-dimensional inverse Fourier transform on the completed frequency-wavenumber domain echo signal to obtain the ground-penetrating radar migration imaging result.

[0119] Based on the maximum amplitude difference With decay threshold The following are methods for completing the echo signal in the new frequency-wavenumber domain by selecting different interpolation methods based on the relationship:

[0120] when At that time, it was determined to be a region with relatively gentle decay, combined with the relation coefficient of the current data block. Using linear interpolation, we obtain ;

[0121] in, This represents the result of linear interpolation, which is the echo signal in the frequency-wavenumber domain after linear interpolation completion. Represents the echo signal in the frequency-wavenumber domain The Middle Average amplitude data of each data block Represents the echo signal in the frequency-wavenumber domain Zhongyu The adjacent first Average amplitude data of each data block The relation coefficients for the current data block;

[0122] when When the signal is identified as a region of strong attenuation, cubic spline interpolation is used. The core of this method is to construct a cubic polynomial so that the interpolation function is not only continuous in terms of function value at the nodes, but also in terms of its first and second derivatives. This allows for a smooth and accurate fit to signals with strong attenuation and abrupt changes.

[0123] Given Echo signal in the frequency-wavenumber domain interpolation nodes The cubic spline interpolation function is as follows:

[0124] ,

[0125] in, This represents the result of cubic spline interpolation, which is the echo signal in the frequency-wavenumber domain after the cubic spline interpolation is completed. For the segmented interval index of cubic spline interpolation, , Representing the spatial-temporal dimension variables in the frequency-wavenumber domain, This represents the x-coordinate of the node in the cubic spline interpolation. , , , The coefficients of the cubic spline function are represented by the coefficients. The interpolation curve is smooth and accurate by solving for continuity of the function at the nodes, continuity of the first derivative, and continuity of the second derivative.

[0126] Will and Unified representation of the echo signal in the completed frequency-wavenumber domain .

[0127] The method for performing a two-dimensional inverse Fourier transform on the completed frequency-wavenumber domain echo signal is as follows:

[0128] ;

[0129] in, This indicates the final ground-penetrating radar offset imaging result. This represents the echo signal in the completed frequency-wavenumber domain. For space wavenumber, The corrected time wavenumber, This represents the two-dimensional inverse Fourier transform.

[0130] This step enables precise imaging of underground targets, effectively eliminating the effects of medium inhomogeneity and wave field attenuation.

[0131] To demonstrate the effectiveness of the method of this invention, a 10-meter-long area of ​​gray granite floor tiles in a school building was selected as the detection scenario. This area contains numerous manhole covers and underground pipes, and a single-channel sweep frequency ground-penetrating radar system was used for detection. Figure 3 Image (a) shows the imaging results of the raw data collected. Figure 3 Image (b) shows the preprocessed imaging result. Then, the preprocessed imaging result is processed using both the traditional FK migration method and the improved FK migration method of this invention. Figure 4As can be seen in (a), the target outline of the underground reflected signal is blurred and locally broken. Because traditional algorithms use a fixed wave velocity, they cannot adapt to the wave velocity variations in non-uniform media, leading to reflected wave deviation. Furthermore, for high dynamic range data acquired by frequency-sweeping ground-penetrating radar, the offset effect is poor, and traditional FK offset is not applicable; Figure 4 In (b), the target outline is sharper and more continuous, the reflection point imaging is more concentrated, there is no obvious trailing, the dynamic wave velocity accurately corrects the propagation path, and the target outline is clearer.

[0132] Figure 4 In (a), the background shows obvious clutter (false signals), which is a direct result of wavefield migration error caused by a fixed wave velocity; while Figure 4 In (b), the background is cleaner; this type of error was eliminated through dynamic wave velocity and targeted interpolation, resulting in a significantly higher imaging signal-to-noise ratio. Figure 4 (a)

[0133] Figure 4 In the middle (a) deep region, the signal is almost completely submerged by noise, making it difficult to distinguish the target; while Figure 4 (b) The deep signal identification is higher because the improved algorithm can dynamically compensate for the energy attenuation of deep electromagnetic waves and repair the signal in the strong attenuation area, thus avoiding the problem of blurred deep imaging in traditional algorithms.

[0134] In summary, the improved algorithm overcomes the shortcomings of traditional fixed wave velocity by adapting to the non-uniformity of the underground medium, resulting in better imaging accuracy and clarity, and facilitating subsequent identification and analysis of underground targets.

[0135] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A phase information based ground penetrating radar signal improved FK migration method, characterized in that, The method comprises the following steps: Step 1, pre-processing the complex domain echo data collected by the single-channel sweep frequency ground penetrating radar along the preset survey line, and then decomposing the pre-processed echo signal to obtain amplitude information and phase information; Step 2: block processing the amplitude information according to time-space dimensions, calculating the average amplitude of each amplitude block and the average amplitude difference between the current amplitude block and the adjacent amplitude block, taking the maximum amplitude difference, and obtaining the amplitude attenuation rate by using the maximum amplitude difference; combining the electromagnetic wave exponential attenuation physical law to obtain the amplitude attenuation model, and setting the attenuation threshold based on the maximum value of the amplitude attenuation model; Step 3: calculating the phase change rate along the time and space dimensions respectively and taking the average to obtain the phase fluctuation rate variance; Step 4, combining the amplitude attenuation rate and the phase fluctuation rate variance, introducing the medium basic dielectric constant and the correction coefficient to obtain the real-time electromagnetic wave velocity; correcting the wave field propagation path based on the real-time electromagnetic wave velocity to obtain the corrected time wave number; Step 5: performing two-dimensional Fourier transform on the pre-processed echo signal in the horizontal plane-time dimension, converting to the frequency-wave number domain to obtain the echo signal in the frequency-wave number domain, and then obtaining the new echo signal in the frequency-wave number domain based on the corrected time wave number; Step 6, according to the relationship between the maximum amplitude difference and the attenuation threshold, different interpolation methods are selected to complete the new echo signal in the frequency-wave number domain, and finally the two-dimensional inverse Fourier transform is performed on the completed echo signal in the frequency-wave number domain to obtain the ground penetrating radar migration imaging result.

2. The method of claim 1, wherein, In step 1, the pre-processing includes removing direct current, background elimination, mean filtering and time gain compensation.

3. The method of claim 1, wherein, In step 1, the pre-processed echo signal is decomposed into amplitude information and phase information , according to the formula: ； ； wherein is the real part of the complex number, is the imaginary part of the complex number.

4. The method of claim 1, wherein, The specific method of step 2 is as follows: amplitude information divided into amplitude blocks, denotes an amplitude block in the time and spatial domain, the th time block, the th spatial block; in the time dimension, is a time block index, each time block contains sample points; in the spatial dimension, is a spatial block index, each spatial block contains sample points; For each block , the average amplitude in the block is calculated , the formula is , is the index used to mark the time direction of the first , the space direction of the first sampling point in the small block. Calculate the difference of the average magnitude of the current magnitude block and the right, lower adjacent magnitude blocks, and take the maximum magnitude difference ; use the maximum magnitude difference to calculate the magnitude decay rate : ； wherein, is a time step, is a space step, is a standard spatial attenuation coefficient, is a standard temporal attenuation coefficient; is a spatial attenuation coefficient, is a temporal attenuation coefficient; According to the amplitude of the electromagnetic wave at the initial position , the propagation time of the electromagnetic wave , the depth of the target , the spatial attenuation coefficient , the time attenuation coefficient , the amplitude attenuation model is obtained ; obtaining an attenuation threshold based on a maximum value of the attenuation model , , denotes an adjustment factor.

5. The method of claim 4, wherein, The specific method of step 3 is as follows: For the phase information, in the time dimension, the phase change rate is calculated as : ， in, Indicates time step Later, in time ,space Phase value at; Indicates the current time ,space Phase value at; In the spatial dimension, the spatial step is calculated as the rate of change of phase : ， wherein, denotes a spatial step later in time , the phase value at a spatial place; Taking the average of both gives the average rate of phase change with time and spatial dimensions , reflecting the overall fluctuation of phase with time and space: ， Further, the phase fluctuation rate variance is calculated : ； wherein, is the total number of time samples, denotes the time sample point index, denotes the spatial sample point index, is the total number of spatial samples, is the global mean of .

6. The method of claim 1, wherein, In step 4, real-time electromagnetic wave velocity The calculation formula is as follows: ； wherein, c is the speed of light, , is a correction factor, is the medium base dielectric constant, is the amplitude decay rate, is the phase fluctuation rate variance, is the relative dielectric constant as a function of time and space.

7. The method of claim 1, wherein, In step 4, the method for calculating the corrected time wave number is as follows: ； wherein, is the corrected time wavenumber; is the real-time electromagnetic wave velocity, is the speed of light, is the spatial wavenumber, is the reference wavenumber.

8. The method of claim 1, wherein, In step 5, the method for converting the pre-processed echo signal to the frequency-wave number domain by two-dimensional Fourier transform in the horizontal plane-time dimension is as follows: ； wherein, is the spatial wave number, is the temporal wave number, is the pre-processed echo signal; denotes the temporal sample point index, denotes the spatial sample point index, denotes the echo signal in the frequency-wave number domain after two-dimensional Fourier transformation; based on the corrected time-wavenumber obtaining a new frequency-wavenumber domain echo signal .

9. The method of claim 1, wherein, In step 6, the different interpolation methods are selected according to the relationship between the maximum amplitude difference and the decay threshold The method of completing the new frequency-wavenumber domain echo signal is as follows: When linear interpolation method, we get ; wherein, represents a linear interpolation result, and is an echo signal in the frequency-wavenumber domain after linear interpolation completion; represents an average amplitude data of the i-th data block in the echo signal in the frequency-wavenumber domain represents an average amplitude data of the i-th data block in the echo signal in the frequency-wavenumber domain represents an average amplitude data of the i-th data block in the echo signal in the frequency-wavenumber domain represents an average amplitude data of the i-th data block in the echo signal in the frequency-wavenumber domain represents an average amplitude data of the i-th data block in the echo signal in the frequency-wavenumber domain represents an average amplitude data of the i-th data block in the echo signal in the frequency-wavenumber domain represents an average amplitude data of the i-th data block in the echo signal in the frequency-wavenumber domain is a correlation coefficient of the current data block; When a cubic spline interpolation method is used, given a frequency-wave number domain echo signal interpolation nodes , the cubic spline interpolation function is as follows: ， wherein, represents the result of cubic spline interpolation, is the echo signal in the frequency-wavenumber domain after the cubic spline interpolation is completed; is the index of the segment interval of the cubic spline interpolation, , represents the space-time dimension variable in the frequency-wavenumber domain, represents the node abscissa of the cubic spline interpolation, , , , represents the coefficient of the cubic spline function; Will be And The echo signal in the frequency-wavenumber domain after completion is uniformly represented as .

10. The method of claim 1, wherein, In step 6, the method for performing two-dimensional inverse Fourier transform on the completed echo signal in the frequency-wave number domain is as follows: ； wherein, represents the final ground penetrating radar migration imaging result, represents the completed frequency-wavenumber domain echo signal, is the spatial wavenumber, is the corrected time wavenumber, represents the two-dimensional inverse Fourier transform.

Images