A method for constructing an initial model for ground-penetrating radar imaging

Abstract

This invention relates to the field of ground-penetrating radar (GPR) detection technology, and particularly to a method for constructing an initial model for GPR imaging. The method includes: acquiring GPR data and preprocessing it; performing frequency-wavenumber domain migration imaging under multi-velocity scanning to obtain migration imaging results at different velocities; constructing a two-dimensional focusing evaluation index system to evaluate all migration imaging results and obtain the optimal velocity and optimal focusing imaging result; extracting subsurface medium structural features based on the optimal focusing imaging result, and constructing an initial GPR imaging model by combining the physical mapping relationship between electromagnetic wave velocity and relative permittivity. This invention solves the technical problems of traditional initial model construction relying on human experience, insufficient prior information, and low matching degree between the model and the actual subsurface structure. The final output initial model can be directly used as input for subsequent GPR imaging or inversion algorithms, effectively improving the focusing of subsequent imaging and the stability and reliability of inversion results.

Description

Technical Field

[0001] This invention relates to the field of ground-penetrating radar (GPR) detection technology, and in particular to a method for constructing an initial model for GPR imaging. Background Technology

[0002] Ground penetrating radar (GPR) is a non-destructive testing technique that utilizes the propagation characteristics of electromagnetic waves in underground media. It is widely used in fields such as underground structure detection, engineering inspection, resource exploration, and archaeology. In GPR imaging and wave equation-based inversion processes, the initial model of the electromagnetic parameters of the underground medium (such as propagation velocity and dielectric constant) has a significant impact on the imaging results and inversion stability.

[0003] In practical applications, due to the complex structure and significant variations in physical properties of underground media, coupled with a lack of accurate prior information, constructing a reasonable initial model becomes a crucial issue restricting the imaging and inversion effects of ground-penetrating radar. If the initial model differs significantly from the actual media, it can easily lead to poor focusing of the imaging results, structural distortion, and even cause the inversion process to fail to converge or get trapped in local optima.

[0004] The existing methods for building initial models are as follows: Model building methods based on experience or simplified assumptions: In current engineering applications, common initial model building methods mainly rely on human experience or simplified medium assumptions. For example, they assign values ​​based on the empirical dielectric constant of typical materials, or assume that the underground medium has a horizontally layered structure and use uniform background parameters. These methods are simple to implement, but they are highly dependent on the operator's experience, make it difficult to reflect the complex structural characteristics of the actual underground medium, and result in poor model adaptability and stability.

[0005] Model building methods based on external information or prior data: Some methods construct an initial model of the subsurface medium by introducing borehole data, geological profiles, or other geophysical exploration results. Although such methods can improve the rationality of the model under certain conditions, their applicability is limited by the availability and consistency of external data, and they are difficult to promote and apply when reliable prior information is lacking or when field conditions are limited.

[0006] Parameter estimation methods based on single imaging or inversion results: Existing technologies also include methods that analyze ground-penetrating radar (GPR) imaging results to estimate the propagation velocity or dielectric constant of the medium, and then construct an initial model accordingly. For example, propagation parameters are selected based on a single migration imaging result or a specific evaluation index. However, such methods typically rely on a single parameter or a single imaging result, making it difficult to comprehensively reflect the changing characteristics of imaging focusing under different parameter conditions. They are also susceptible to noise, anomalies, or imaging assumptions, resulting in limited model stability and reliability.

[0007] Model building methods based on inversion algorithms: Some studies have attempted to directly obtain electromagnetic parameter models of subsurface media using inversion algorithms. However, these methods are usually highly dependent on the initial model and have high computational complexity. They are also prone to unstable results if the initial model is unreasonable. Therefore, they are inherently inseparable from the problem of initial model construction.

[0008] Based on the above, it is necessary to provide a method for constructing an initial model for ground-penetrating radar imaging in order to improve the rationality and applicability of the initial model. Summary of the Invention

[0009] The main objective of this invention is to provide a method for constructing an initial model for ground-penetrating radar imaging. The specific technical solution is as follows: A method for constructing an initial model for ground-penetrating radar imaging includes the following steps: Step S1: Acquire ground-penetrating radar data and perform preprocessing; Step S2: Perform frequency-wavenumber domain migration imaging under multi-velocity scanning to obtain migration imaging results at different velocities; Step S3: Construct a two-dimensional focusing evaluation index system, evaluate all offset imaging results, and obtain the optimal speed and optimal focusing imaging results; Step S4: Extract the structural features of the underground medium based on the optimal focusing imaging results, combine the physical mapping relationship between electromagnetic wave velocity and relative permittivity, complete the assignment of the relative permittivity of the underground medium and the region filling, and finally construct the initial model of ground penetrating radar imaging.

[0010] Furthermore, step S1 includes: Step S1.1: Acquire two-dimensional ground-penetrating radar profile data ,in: The spatial coordinates of the survey line direction represent the horizontal position of the underground medium; This is the two-way propagation time of electromagnetic waves in the underground medium; Step S1.2: Process the two-dimensional profile data from the ground-penetrating radar. Remove background: ; in: The data is a two-dimensional ground-penetrating radar profile after background removal. This represents the total number of sampling points along the survey line direction; For the direction of the survey line Spatial coordinates of each sampling point ; Step S1.3: Perform power function gain processing on the background-removed ground-penetrating radar 2D profile data: ; in: This is the two-dimensional ground-penetrating radar profile data after power function gain processing; This is the gain coefficient; Step S1.4: Perform Fourier transform on the ground-penetrating radar two-dimensional profile data after power function gain processing: ; in: Ground-penetrating radar data in the frequency-wavenumber domain; The spatial wavenumber along the survey line represents the wavefield variation characteristics in the horizontal direction. ω is the angular frequency of the electromagnetic wave; The imaginary unit, ; It is a natural constant.

[0011] Furthermore, step S2 includes: Step S2.1: Preset a discrete velocity set that covers the range of electromagnetic wave propagation speeds in actual underground media. ,in, The total number of discrete velocities. For the first Each scan speed, ; Step S2.2: For any scan rate Establish the frequency-wavenumber domain offset dispersion relation: ; in: The vertical wavenumber characterizes the wavefield variation characteristics in the vertical depth direction. Step S2.3, for Make a judgment: when When ≤0, it indicates the selected scanning speed. If the requirements are not met, return to step S2.2, reselect the scanning speed, and repeat the process. Calculate; when If the value is greater than 0, proceed to step S2.4; Step S2.4: Calculate the scanning speed Corresponding phase shift operator : ; in: The vertical depth coordinates of the underground medium; It is a natural constant; Step S2.5: Transfer the ground-penetrating radar data in the frequency-wavenumber domain. Corresponding scanning speed Phase shift operator Multiply the data, then transform it back to the time-space domain using a two-dimensional inverse Fourier transform to obtain the frequency-wavenumber domain migration imaging result at this scanning speed. The calculation formula is as follows: ; in: For scanning speed Ground-penetrating radar offset imaging results; Step S2.6: For the velocity set For each scan velocity, steps S2.2 to S2.5 are executed sequentially to complete frequency-wavenumber domain migration imaging across the entire velocity range, ultimately yielding a migration imaging result set that corresponds one-to-one with the velocity set. .

[0012] Furthermore, step S3 includes: Step S3.1: Process the offset imaging results set Each imaging result Preprocessing is required; Step S3.2: Calculate the radial energy concentration (REC) and the in-phase axis sharpness index (HSI) based on the preprocessed imaging results, and fuse the radial energy concentration (REC) and the in-phase axis sharpness index (HSI) to obtain the focus quality factor (FQF). Step S3.3: Select the optimal speed and optimal focusing imaging result based on the focusing quality factor FQF.

[0013] Furthermore, step S3.1 includes: ① Set of offset imaging results Each imaging result An analytical signal is constructed using Hilbert transform, and the amplitude envelope of the imaging result is extracted. The calculation formula is as follows: ; ; in: For imaging results The corresponding analytical signal; For Hilbert transform operators; The amplitude envelope matrix of the imaging result; Modular operation for complex numbers; ②The amplitude envelope matrix Normalization is performed: ; in: This is the normalized amplitude envelope matrix; The amplitude envelope matrix The minimum value in; The amplitude envelope matrix The maximum value in.

[0014] Furthermore, in step S3.2, the radial energy concentration REC is calculated as follows: ; in: For scanning speed Radial energy concentration of the imaging results; For the imaging results in the horizontal direction The number of sampling points on; For the imaging results in the vertical depth direction Number of sampling points on; The in-phase axis sharpness index (HSI) is calculated as follows: ① Calculate the lateral gradient: ; in: This represents the gradient value of the normalized amplitude envelope matrix in the horizontal direction. Horizontal direction The sampling interval on; Indicates spatial location The normalized amplitude envelope matrix at the point; Indicates spatial location The normalized amplitude envelope matrix at the point; ② Calculate the sharpness index of the in-phase axis: ; in: For scanning speed The in-phase axis sharpness index of the imaging results; for The gradient mean; for The gradient standard deviation; The focus quality factor (FQF) is calculated as follows: ① Normalize the sharpness index of the in-phase axis: ; in: For scanning speed The normalized in-phase axis sharpness index; For velocity set The minimum value of the sharpness index of all in-phase axes; For velocity set The maximum value of the sharpness index of all in-phase axes; ② Calculate the focusing quality factor FQF: ; in: For scanning speed The focus quality factor of the imaging results.

[0015] Furthermore, step S3.3 specifically involves: using the Focused Quality Factor (FQF) as the core evaluation criterion, and evaluating the velocity set. The FQF values ​​corresponding to all scanning speeds are analyzed, and the optimal speed is automatically selected. The imaging result corresponding to the optimal speed is the one with the best focusing effect across the entire speed range, i.e., the optimal focusing imaging result. The rules for selecting the optimal speed are as follows: ①If the FQF value is in the velocity set If there exists a single global maximum value, then the scanning speed corresponding to that maximum value is the optimal speed, denoted as . The formula is: ; in: To achieve the optimal speed; This represents the operator that finds the independent variable that maximizes the function; ②If the FQF value is in the velocity set If there are multiple local maxima, then the scanning speed corresponding to each local maximum is taken as the optimal speed, denoted as the optimal speed set. , The number of local maxima; ; in: This represents the optimal speed set.

[0016] Furthermore, in step S4, based on the optimal focused imaging results, the structural features of the underground medium are extracted using the Canny edge detection algorithm: the normalized amplitude envelope matrix corresponding to the optimal velocity is... As input, Gaussian smoothing, gradient calculation, non-maximum suppression, and double thresholding are performed sequentially to extract the structural edge set of the subsurface medium; as detailed below: ① The normalized amplitude envelope matrix of the input Gaussian smoothing is performed using the following formula: ; in: This is the normalized amplitude envelope matrix after Gaussian smoothing; This is a two-dimensional convolution operation; It is a two-dimensional Gaussian kernel. The standard deviation of the Gaussian kernel; ② Calculate the horizontal and vertical gradients of the normalized amplitude envelope matrix after Gaussian smoothing, and then fuse the horizontal and vertical gradients to obtain the gradient magnitude and gradient direction. The calculation formula is as follows: ; ; ; ; in: This represents the horizontal gradient value of the normalized amplitude envelope matrix after Gaussian smoothing; This represents the vertical gradient value of the normalized amplitude envelope matrix after Gaussian smoothing. Vertical depth direction z The sampling interval is a fixed constant. The gradient magnitude represents the intensity of the edge. The gradient direction represents the orientation of the edge. Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point; Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point; Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point; Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point; ③ Gradient magnitude along the gradient direction By performing point-by-point judgment, only local maxima along the gradient direction are retained, and gradient values ​​at non-edge points are suppressed to obtain the gradient magnitude after non-maximum suppression. ④ Set a high threshold and low threshold , > The gradient magnitude after non-maximum suppression is compared with the high threshold. and low threshold Comparison: Gradient magnitude greater than The points are strong edge points and are directly identified as the structural edge of the underground medium; Gradient magnitude is between and The points between them are weak edge points. Only when a point is connected to a strong edge point is it determined to be a structural edge of the underground medium. Gradient magnitude less than Points that are not edge points are directly removed; Through the above operations, the final set of structural edges of the subsurface medium is obtained, denoted as... ,in This indicates that the location is a structural boundary. This indicates that the location is a non-structural boundary.

[0017] Furthermore, in step S4, the physical mapping relationship between electromagnetic wave velocity and relative permittivity is as follows: ; The formula for calculating the relative permittivity, obtained by transforming the above equation, is as follows: ; in: This represents the actual propagation speed of electromagnetic waves in underground media; The speed at which electromagnetic waves propagate in a vacuum; The relative permittivity of the underground medium; Optimal speed Substituting into the above equation, we obtain the relative permittivity value corresponding to the optimal velocity, denoted as . : ; If an optimal velocity set exists For each optimal speed , Calculate the relative permittivity corresponding to the optimal velocity. .

[0018] Furthermore, in step S4, based on the structural edge set The underground space is divided into multiple independent medium regions, and then the relative permittivity is considered. The relative permittivity of each region was assigned and the background was filled, thus constructing an initial model of the relative permittivity of the subsurface medium. That is, the initial model for ground-penetrating radar imaging, specifically: From the set of structural edges The enclosed region is denoted as Each closed region represents an independent underground medium. For each closed region Assign it the relative permittivity obtained by the corresponding optimal velocity conversion. If multiple optimal velocities exist, assign corresponding velocities to different regions based on the target of the offset focusing. ; For the background region outside the closed region, the relative permittivity value corresponding to the scan rate of the relatively stable portion of the focus quality factor (FQF) is selected as the background fill value. ; The final initial model of the relative permittivity of the subsurface medium for: .

[0019] The beneficial effects achieved by this solution are: This invention uses ground-penetrating radar (GPR) measured data as the core input, without relying on precise prior medium parameters or external geological data. Through three core steps—frequency-wavenumber domain (FK) migration imaging under multi-velocity scanning, calculation of two-dimensional focusing evaluation index and selection of optimal velocity, and Canny edge detection and construction of an initial dielectric constant model—it automatically completes the construction of an initial model of the relative dielectric constant of the subsurface medium under data-driven conditions. This solves the technical problems of traditional initial model construction relying on human experience, insufficient prior information, and low matching degree between the model and the actual subsurface structure. The final output initial model can be directly used as input for subsequent GPR imaging or inversion algorithms, effectively improving the focusing of subsequent imaging and the stability and reliability of inversion results. Attached Figure Description

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

[0021] Figure 1 This is a flowchart illustrating the method for constructing the initial model for ground-penetrating radar imaging in an embodiment of the present invention. Figure 2 This is a flowchart of discrete velocity migration imaging in an embodiment of the present invention; Figure 3 This is a flowchart of the focused evaluation analysis and initial model construction in an embodiment of the present invention; Figure 4 This is a model of undulating strata and isolated boulders designed in an embodiment of the present invention; Figure 5 The data in this embodiment of the invention is the original data after gaining; Figure 6 This is a diagram showing the optimal velocity offset result; Figure 7 This is a gradient magnitude plot; Figure 8 Image showing the edge detection results; Figure 9 The diagram shows a comparison between the empirical model and the model of this invention, where: a) is the uniform model, b) is the gradient model, and c) is the model constructed in this invention.

[0022] The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0023] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0024] Example See Figure 1 This invention proposes a method for constructing an initial model for ground-penetrating radar imaging, comprising the following steps: Step S1: Acquire ground-penetrating radar data and perform preprocessing; including: Step S1.1: Acquire two-dimensional ground-penetrating radar profile data ,in: The spatial coordinates of the survey line direction represent the horizontal position of the underground medium; It is the two-way propagation time of electromagnetic waves in the underground medium, that is, the total propagation time of electromagnetic waves from the transmitting antenna, after being reflected by the underground medium, and back to the receiving antenna.

[0025] Step S1.2: Process the two-dimensional profile data from the ground-penetrating radar. Remove background: ; in: The data is a two-dimensional ground-penetrating radar profile after background removal. This represents the total number of sampling points along the survey line direction; For the direction of the survey line Spatial coordinates of each sampling point ; The core operation of background removal is to perform background removal on each two-way propagation time. The mean of all sampling points in the corresponding horizontal direction is calculated, and then the mean is subtracted from the original data to eliminate background noise in the horizontal direction.

[0026] Step S1.3: Considering the geometric diffusion attenuation and medium absorption attenuation during the propagation of electromagnetic waves underground, the reflected signal energy of deep media is weak and easily masked. Therefore, power function gain processing is applied to the background-removed ground-penetrating radar two-dimensional profile data to enhance the weak deep reflection signal and balance the energy distribution of the entire profile. ; in: This is the two-dimensional ground-penetrating radar profile data after power function gain processing; is the gain coefficient, a constant set according to the actual detection scenario, which can be adjusted to match the attenuation characteristics of different underground media.

[0027] Step S1.4: To convert the radar data from the time-space domain to the frequency-wavenumber domain and meet the mathematical requirements of frequency-wavenumber (FK) migration imaging, a Fourier transform is performed on the ground-penetrating radar two-dimensional profile data after power function gain processing. The calculation formula is as follows: ; in: Ground-penetrating radar data in the frequency-wavenumber domain; The spatial wavenumber along the survey line represents the wavefield variation characteristics in the horizontal direction. ω is the angular frequency of the electromagnetic wave; The imaginary unit, , It is a natural constant.

[0028] Step S2: Perform frequency-wavenumber domain migration imaging under multi-velocity scanning to obtain migration imaging results at different velocities; See Figure 2 Specifically, it includes the following steps: Step S2.1: Preset a discrete velocity set that covers the range of electromagnetic wave propagation speeds in actual underground media. ,in, The total number of discrete velocities. For the first Each scan speed, .

[0029] Step S2.2: For any scan rate Establish the frequency-wavenumber domain offset dispersion relation: ; in: The vertical wavenumber characterizes the wavefield variation characteristics in the vertical depth direction. The scanning velocity parameter is any discrete velocity from the preset velocity set, representing the assumed underground propagation velocity of electromagnetic waves in this migration imaging.

[0030] Step S2.3, for Make a judgment: when When ≤0, it indicates the selected scanning speed. If the requirements are not met, return to step S2.2, reselect the scanning speed, and repeat the process. Calculate; when If the value is greater than 0, proceed to step S2.4.

[0031] Step S2.4: The phase shift operator is the core operator for FK migration imaging. Its function is to perform phase correction on the data in the frequency-wavenumber domain to restore the actual reflection characteristics of the subsurface medium. For each scanning velocity in the preset velocity set... Construct a corresponding phase shift operator separately. The calculation formula is: ; in: The vertical depth coordinates of the underground medium; It is a natural constant.

[0032] Step S2.5: Transfer the ground-penetrating radar data in the frequency-wavenumber domain. Corresponding scanning speed Phase shift operator Multiply the data, then transform it back to the time-space domain using a two-dimensional inverse Fourier transform to obtain the frequency-wavenumber domain migration imaging result at this scanning speed. The calculation formula is as follows: ; in: For scanning speed The ground-penetrating radar offset imaging results are presented below, forming a two-dimensional matrix representing different horizontal positions. x Different vertical depths z The intensity of the reflected signal at that location.

[0033] Step S2.6: For the velocity set For each scan velocity, steps S2.2 to S2.5 are executed sequentially to complete frequency-wavenumber domain migration imaging across the entire velocity range, ultimately yielding a migration imaging result set that corresponds one-to-one with the velocity set. .

[0034] Step S3: Construct a two-dimensional focusing evaluation index system, evaluate all offset imaging results, and obtain the optimal speed and optimal focusing imaging results; The core objective of this step is to construct a two-dimensional focusing evaluation index system specifically for ground-penetrating radar FK migration imaging, to quantitatively evaluate all imaging results obtained in step S2, and to objectively and automatically select the imaging result with the best focusing characteristics and the corresponding optimal speed, thus solving the problem of traditional methods relying on manual judgment of focusing effect.

[0035] The focusing effect of offset imaging directly reflects the degree of matching between the imaging results and the actual underground structure: when the scanning speed... When the imaging speed approaches the actual electromagnetic wave propagation speed of the underground medium, the imaging results show good focusing, characterized by convergence of the reflection phase axis, concentrated energy, and high edge sharpness; when the scanning speed... When the speed deviates from the true speed, the imaging results are poorly focused, manifested as widening of the reflection phase axis, energy dispersion, and low edge sharpness.

[0036] Based on the physical characteristics of the focusing effect described above, this step first preprocesses the imaging results, then constructs two basic evaluation indicators: Radial Energy Concentration (REC) and In-phase Axis Sharpness Index (HSI). Finally, these are fused to obtain the Focus Quality Factor (FQF) as the core evaluation indicator. The optimal speed is automatically selected through the quantization value of FQF. (See [link to relevant documentation]). Figure 3 Specifically, it includes the following steps: Step S3.1: Process the offset imaging results set Each imaging result Preprocessing is performed to eliminate phase interference and normalize the data to a uniform range to ensure the comparability of subsequent evaluation indicators. This includes: ① Set of offset imaging results Each imaging result An analytical signal is constructed using Hilbert transform, and the amplitude envelope of the imaging result is extracted to eliminate the interference of phase on the evaluation of focusing effect. The calculation formula is as follows: ; ; in: For imaging results The corresponding analytical signal; For Hilbert transform operators; The amplitude envelope matrix of the imaging result retains only amplitude information and eliminates phase interference; Modular operation for complex numbers; ②The amplitude envelope matrix Normalization is performed to map the data to the [0,1] interval, eliminating amplitude differences in imaging results at different scanning speeds and ensuring the comparability of evaluation metrics across the entire speed range. ; in: This is the normalized amplitude envelope matrix, which serves as the basis for subsequent calculations of focused evaluation indicators; The amplitude envelope matrix The minimum value in; The amplitude envelope matrix The maximum value in.

[0037] Step S3.2: Calculate the radial energy concentration (REC) and the in-phase axis sharpness index (HSI) based on the preprocessed imaging results, and fuse the REC and HSI to obtain the focus quality factor (FQF); specifically: Radial energy concentration (REC) characterizes the focusing effect from the energy distribution dimension, quantifying the degree of concentration of reflected signal energy radially. The better the focusing effect, the more concentrated the energy, and the closer the REC value is to 1; the worse the focusing effect, the more dispersed the energy, and the closer the REC value is to 0. The radial energy concentration REC is calculated as follows: ; in: For scanning speed The radial energy concentration of the imaging result, with a value range of [0,1]; For the imaging results in the horizontal direction Number of sampling points on; For the imaging results in the vertical depth direction The number of sampling points.

[0038] The in-phase axis sharpness index (HSI) characterizes the focusing effect from the edge gradient dimension, quantifying the edge sharpness of the reflecting in-phase axis. The better the focusing effect, the more concentrated the gradient at the edge of the in-phase axis, and the larger the HSI value; the worse the focusing effect, the more diffuse the gradient, and the smaller the HSI value. After subsequent HSI normalization, the better the focusing effect and the more concentrated the energy, the closer the normalized value is to 0. The in-phase axis sharpness index (HSI) is calculated as follows: ① Calculate the lateral gradient: ; in: The gradient value of the normalized amplitude envelope matrix in the horizontal direction represents the rate of change of the edge of the in-phase axis in the horizontal direction. Horizontal direction The sampling interval is a fixed constant. Indicates spatial location The normalized amplitude envelope matrix at the point; Indicates spatial location The normalized amplitude envelope matrix at the point; ② Calculate the sharpness index of the in-phase axis: ; in: For scanning speed The in-phase axis sharpness index of the imaging result is a positive real number; the larger the value, the higher the edge sharpness. for The gradient mean, ; for The gradient standard deviation, .

[0039] The Focus Quality Factor (FQF) is a fusion evaluation index that normalizes and fuses the Radial Energy Concentration (REC) and In-Axis Sharpness Index (HSI), while simultaneously constraining energy concentration characteristics and edge stability characteristics. It outputs a single scalar value ranging from [0,1]. A larger value indicates better focusing of the imaging result and is the core basis for subsequent optimal velocity selection. The Focus Quality Factor (FQF) is calculated as follows: ① Normalize the sharpness index of the in-phase axis: ; in: For scanning speed The normalized in-phase axis sharpness index is defined with a value range of [0,1]. For velocity set The minimum value of the sharpness index of all in-phase axes; For velocity set The maximum value of the sharpness index of all in-phase axes.

[0040] ② Calculate the focusing quality factor FQF: ; in: For scanning speed The focusing quality factor of the imaging result has a value range of [0,1]. The larger the value, the better the focusing effect.

[0041] Step S3.3: Select the optimal speed and optimal focusing imaging result based on the focusing quality factor FQF.

[0042] Using the Focused Quality Factor (FQF) as the core evaluation criterion, the velocity set is evaluated. The FQF values ​​corresponding to all scanning speeds are analyzed, and the optimal speed is automatically selected. The imaging result corresponding to the optimal speed is the one with the best focusing effect across the entire speed range, i.e., the optimal focusing imaging result. The rules for selecting the optimal speed are as follows: ①If the FQF value is in the velocity set If there exists a single global maximum value, then the scanning speed corresponding to that maximum value is the optimal speed, denoted as . The formula is: ; in: The optimal speed is the scanning speed that best represents the actual electromagnetic wave propagation speed in the underground medium. This represents the operator that finds the independent variable that maximizes the function; ②If the FQF value is in the velocity set If there are multiple local maxima, then the scanning speed corresponding to each local maximum is taken as the optimal speed, denoted as the optimal speed set. , The number of local maxima; ; in: This represents the optimal speed set.

[0043] After selecting the optimal velocity (or the optimal velocity set), extract the corresponding imaging result as the base image for the subsequent step S4.

[0044] Step S4: Extract the structural features of the underground medium based on the optimal focusing imaging results, combine the physical mapping relationship between electromagnetic wave velocity and relative permittivity, complete the assignment of the relative permittivity of the underground medium and the region filling, and finally construct the initial model of ground penetrating radar imaging.

[0045] The core objective of this step is to extract the true structural boundaries of the subsurface medium based on the optimal focused imaging results selected in step S3, using Canny edge detection. Then, combining this with the physical mapping relationship between electromagnetic wave velocity and relative permittivity, the relative permittivity of the subsurface medium is assigned and the region is filled, ultimately constructing an initial ground-penetrating radar imaging model with structural constraints and property indications. The specific operations are divided into three stages: Canny edge detection, determination of the velocity-permeability mapping relationship, and initial model construction and filling. See also... Figure 3 The details are as follows: Canny edge detection: Canny edge detection is a high-precision edge extraction algorithm that effectively suppresses spurious scattering signals and improves the continuity and integrity of structural boundaries. Based on optimal focusing imaging results, the Canny edge detection algorithm is used to extract subsurface medium structural features: the normalized amplitude envelope matrix corresponding to the optimal velocity is used... As input, Gaussian smoothing, gradient calculation, non-maximum suppression, and double thresholding are performed sequentially to extract the structural edge set of the subsurface medium; as detailed below: ① The normalized amplitude envelope matrix of the input Gaussian smoothing is applied to eliminate noise interference in the image and prevent noise from being falsely detected as edges. The calculation formula is as follows: ; in: This is the normalized amplitude envelope matrix after Gaussian smoothing; This is a two-dimensional convolution operation; It is a two-dimensional Gaussian kernel, and it is a smoothing operator. ; is the standard deviation of the Gaussian kernel, and is a constant set according to the noise level of the imaging results.

[0046] ② Calculate the horizontal and vertical gradients of the normalized amplitude envelope matrix after Gaussian smoothing, and then fuse the horizontal and vertical gradients to obtain the gradient magnitude and gradient direction. The calculation formula is as follows: ; ; ; ; in: This represents the horizontal gradient value of the normalized amplitude envelope matrix after Gaussian smoothing; This represents the vertical gradient value of the normalized amplitude envelope matrix after Gaussian smoothing. Vertical depth direction z The sampling interval is a fixed constant. The gradient magnitude represents the intensity of the edge. The gradient direction represents the orientation of the edge. Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point; Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point; Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point; Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point.

[0047] ③ Gradient magnitude along the gradient direction By performing point-by-point judgment, only the local maximum value in the gradient direction is retained, and the gradient value of non-edge points is suppressed to refine the edge, eliminate the edge blurring problem, and obtain the gradient magnitude after non-maximum suppression.

[0048] ④ Set a high threshold and low threshold , > The gradient magnitude after non-maximum suppression is compared with the high threshold. and low threshold Comparison: Gradient magnitude greater than The points are strong edge points and are directly identified as the structural edge of the underground medium; Gradient magnitude is between and The points between them are weak edge points. Only when a point is connected to a strong edge point is it determined to be a structural edge of the underground medium. Gradient magnitude less than Points that are not edge points are directly removed; Through the above operations, the final set of structural edges of the subsurface medium is obtained, denoted as... ,in This indicates that the location is a structural boundary. This indicates that the location is a non-structural boundary.

[0049] The velocity-dielectric constant mapping relationship is determined as follows: The propagation speed of electromagnetic waves in underground media has a fixed physical mapping relationship with the relative permittivity of the medium. This relationship provides a theoretical basis for converting the optimal velocity into a medium property parameter (relative permittivity). The physical mapping relationship between electromagnetic wave velocity and relative permittivity is as follows: ; The formula for calculating the relative permittivity, obtained by transforming the above equation, is as follows: ; in: This represents the actual propagation speed of electromagnetic waves in underground media; Let be the speed of electromagnetic waves in a vacuum, and be a physical constant. ≈3×10⁸ m / s; The relative permittivity of the underground medium is a core physical property parameter for ground-penetrating radar imaging and inversion. Optimal speed Substituting into the above equation, we obtain the relative permittivity value corresponding to the optimal velocity, denoted as . : ; in: The relative permittivity corresponding to the optimal velocity characterizes the core physical properties of the underground medium.

[0050] If an optimal velocity set exists For each optimal speed , Calculate the relative permittivity corresponding to the optimal velocity. .

[0051] Initial model building and population: Based on structural edge sets The underground space is divided into multiple independent medium regions, and then the relative permittivity is considered. The relative permittivity of each region was assigned and the background was filled, thus constructing an initial model of the relative permittivity of the subsurface medium. That is, the initial model for ground-penetrating radar imaging, specifically: From the set of structural edges The enclosed region is denoted as Each closed region represents an independent underground medium. For each closed region Assign it the relative permittivity obtained by the corresponding optimal velocity conversion. If multiple optimal velocities exist, assign corresponding velocities to different regions based on the target of the offset focusing. ; For the background region outside the closed region, the relative permittivity value corresponding to the scan rate of the relatively stable portion of the focus quality factor (FQF) is selected as the background fill value. The background fill value ensures the continuity and rationality of the model; The final initial model for the relative permittivity of the subsurface medium is as follows: ; in: The initial model for the relative permittivity of the underground medium is a two-dimensional matrix representing different horizontal positions. Different vertical depths The relative permittivity value at that location.

[0052] To verify the effectiveness of the method of the present invention, the following simulation experiment was conducted: A set of undulating strata and isolated rock models were designed, such as Figure 4 As shown: its horizontal length is 20m and its depth is 10m. The model consists of three layers, with the upper layer being a layered medium with a thickness of 3m and a relative permittivity of [missing value]. The intermediate layer is approximately 1.5m thick, and its contact surface with the upper and lower strata is undulating. The relative permittivity of this layer is... The thickness of the bottom stratum is approximately 5.5m, and its relative permittivity is... The underlying medium contains two irregular clay anomalies with a relative permittivity of [value missing]. The depth of the anomaly is approximately 7-9m. This model considers the shielding effect of the irregular anomaly in the undulating strata with high ohmic loss, in order to simulate the gravel exploration under strong attenuation in multiple strata that is prone to occur in actual exploration.

[0053] The gprMax algorithm, implemented using the finite-difference time-domain algorithm, was used to perform forward modeling of the undulating strata boulder model. An antenna center frequency of 100 MHz and a time window of 200 ns were used. A total of 970 A-scans were acquired with a channel spacing of 0.02 m. The amplified data obtained from the raw data is shown below. Figure 5 As shown; The wavefield obtained from the forward modeling was migrated and edge detected. The optimal velocity migration result is shown in [link to relevant documentation]. Figure 6 See gradient magnitude results. Figure 7 See edge detection results. Figure 8 .

[0054] Figure 9 The core model comparison results of this numerical simulation experiment directly compare the effects of the initial model constructed by the traditional model construction method and the method of this invention, clearly demonstrating the significant advantages of this invention in the detailed characterization of underground media structure and physical properties.

[0055] Traditional model building methods are all based on simplified assumptions made by human experience and do not incorporate the true wavefield characteristics of ground-penetrating radar data. As a result, the models they construct deviate significantly from the actual structure and physical property distribution of the subsurface medium. The specific defects are as follows: Uniform models (e.g.) Figure 9 (a) completely ignores the stratification, property differentiation, and local anomalies of the underground medium, simplifying the complex underground space into a single homogeneous medium. This approach is only applicable to extremely simple detection scenarios. In actual engineering, the underground medium inevitably exhibits changes in structure and properties. Such models lead to severe distortion of subsequent imaging results and difficulty in converging the inversion process, which is one of the core reasons for the low inversion accuracy of traditional ground-penetrating radar.

[0056] Gradient models (such as) Figure 9 b) of the model considers the variation trend of physical properties with depth, but it can only depict continuous and smooth gradient changes. It cannot reflect the discontinuous structural boundaries of the subsurface medium, nor can it depict local abrupt changes in physical properties. For exploration scenarios with complex subsurface structures and significant changes in physical properties, this model still suffers from serious structural information loss, resulting in extremely poor detail depiction capabilities in subsequent imaging and inversion results.

[0057] This invention uses ground-penetrating radar measured / simulated data as the sole driving force. It obtains the optimal velocity matching the real medium through multi-velocity scanning and focusing evaluation, and combines this with Canny edge detection to extract precise structural boundaries. The final initial model is entirely based on the real wavefield characteristics reflected by the data, without any artificial simplification assumptions. It highly matches the real structure and physical property distribution of the underground medium. Compared with traditional models, its core advantages lie in two dimensions: structural characterization and physical property distribution, as detailed below: 1. Accurately characterize the discontinuous structural boundaries of underground media. Figure 9 As can be clearly seen in (c), the model constructed by this invention can accurately reproduce the stratification interfaces, structural abrupt changes in the depth direction, and the contour boundaries of local anomalies in the subsurface medium. The boundary positions are clear and the continuity is good, consistent with the preset real forward model ( Figure 4The structural features are highly consistent. This advantage stems from the structural extraction path of the present invention, which is "optimal focused imaging + Canny edge detection": the optimal focused imaging result eliminates the problems of in-phase axis broadening and energy dispersion caused by velocity deviation, providing a clear basic image for edge detection; Canny edge detection, combined with Gaussian smoothing, non-maximum suppression, and dual threshold connection, effectively suppresses false scattering signals, ensuring the accuracy and continuity of structural boundary extraction.

[0058] The traditional model ( Figure 9 a) and b) in the above have no ability to delineate structural boundaries and cannot reflect the actual structural characteristics of the underground medium at all.

[0059] 2. Accurately reproduce the physical property distribution characteristics of underground media. The distribution of relative permittivity directly reflects the differences in physical properties of underground media. The model constructed in this invention ( Figure 9 c) can perform property matching based on wave field data, assigning relative permittivity values ​​that match the actual properties to different structural regions, which not only characterizes the differences in properties in different regions, but also ensures the rationality of properties in the same region.

[0060] In traditional models, the uniform model ( Figure 9 a)) The physical properties are consistent throughout the entire space, gradient model ( Figure 9 b) in the above can only reflect the overall gradient change and cannot restore the actual physical property distribution of the underground medium, which deviates greatly from the real situation.

[0061] 3. Adaptable to complex underground media detection scenarios Figure 9 The results in (c) demonstrate that the method of this invention does not rely on prior geological information; it can construct an initial model of complex structures and diverse physical properties of underground media using only ground-penetrating radar data. This solves the technical challenge of traditional models being difficult to reasonably define in scenarios with complex underground media structures and significant changes in physical properties. Whether it is a layered rock and soil mass or an underground space with local anomalies, this invention can accurately characterize its structural and physical property features, providing high-quality initial input for subsequent fine imaging and inversion.

[0062] The above description is only a preferred embodiment of the present invention and does not limit the scope of protection of the present invention. All equivalent structural transformations made under the inventive concept of the present invention using the contents of the present invention specification and drawings, or direct / indirect applications in other related technical fields, are included within the scope of protection of the present invention.

Claims

1. A method for constructing an initial model for ground-penetrating radar imaging, characterized in that, Includes the following steps: Step S1: Acquire ground-penetrating radar data and perform preprocessing; Step S2: Perform frequency-wavenumber domain migration imaging under multi-velocity scanning to obtain migration imaging results at different velocities; Step S3: Construct a two-dimensional focusing evaluation index system, evaluate all offset imaging results, and obtain the optimal speed and optimal focusing imaging results; Step S4: Extract the structural features of the underground medium based on the optimal focusing imaging results, combine the physical mapping relationship between electromagnetic wave velocity and relative permittivity, complete the assignment of the relative permittivity of the underground medium and the region filling, and finally construct the initial model of ground penetrating radar imaging.

2. The method for constructing an initial model for ground-penetrating radar imaging according to claim 1, characterized in that, Step S1 includes: Step S1.1: Acquire two-dimensional ground-penetrating radar profile data ,in: The spatial coordinates of the survey line direction represent the horizontal position of the underground medium; This is the two-way propagation time of electromagnetic waves in the underground medium; Step S1.2: Process the two-dimensional profile data from the ground-penetrating radar. Remove the background: ； in: The data is a two-dimensional ground-penetrating radar profile after background removal. This represents the total number of sampling points along the survey line direction; For the direction of the survey line Spatial coordinates of each sampling point ; Step S1.3: Perform power function gain processing on the background-removed ground-penetrating radar 2D profile data: ； in: This is the two-dimensional ground-penetrating radar profile data after power function gain processing; This is the gain coefficient; Step S1.4: Perform Fourier transform on the ground-penetrating radar two-dimensional profile data after power function gain processing: ； in: Ground-penetrating radar data in the frequency-wavenumber domain; The spatial wavenumber along the survey line represents the wavefield variation characteristics in the horizontal direction. ω is the angular frequency of the electromagnetic wave; The imaginary unit, ; It is a natural constant.

3. The method for constructing an initial model for ground-penetrating radar imaging according to claim 2, characterized in that, Step S2 includes: Step S2.1: Preset a discrete velocity set that covers the range of electromagnetic wave propagation speeds in actual underground media. ,in, The total number of discrete velocities. For the first Each scan speed, ; Step S2.2: For any scan rate Establish the frequency-wavenumber domain offset dispersion relation: ； in: The vertical wavenumber characterizes the wavefield variation characteristics in the vertical depth direction. Step S2.3, for Make a judgment: when When ≤0, it indicates the selected scanning speed. If the requirements are not met, return to step S2.2, reselect the scanning speed, and repeat the process. Calculate; when If the value is greater than 0, proceed to step S2.4; Step S2.4: Calculate the scanning speed Corresponding phase shift operator : ； in: The vertical depth coordinates of the underground medium; It is a natural constant; Step S2.5: Transfer the ground-penetrating radar data in the frequency-wavenumber domain. Corresponding scanning speed Phase shift operator Multiply the data, then transform it back to the time-space domain using a two-dimensional inverse Fourier transform to obtain the frequency-wavenumber domain migration imaging result at this scanning speed. The calculation formula is as follows: ； in: For scanning speed Ground-penetrating radar offset imaging results; Step S2.6: For the velocity set For each scan velocity, steps S2.2 to S2.5 are executed sequentially to complete frequency-wavenumber domain migration imaging across the entire velocity range, ultimately yielding a migration imaging result set that corresponds one-to-one with the velocity set. .

4. The method for constructing an initial model for ground-penetrating radar imaging according to claim 3, characterized in that, Step S3 includes: Step S3.1: Process the offset imaging results set Each imaging result Preprocessing is required; Step S3.2: Calculate the radial energy concentration (REC) and the in-phase axis sharpness index (HSI) based on the preprocessed imaging results, and fuse the radial energy concentration (REC) and the in-phase axis sharpness index (HSI) to obtain the focus quality factor (FQF). Step S3.3: Select the optimal speed and optimal focusing imaging result based on the focusing quality factor FQF.

5. The method for constructing an initial model for ground-penetrating radar imaging according to claim 4, characterized in that, Step S3.1 includes: ① Set of offset imaging results Each imaging result An analytical signal is constructed using Hilbert transform, and the amplitude envelope of the imaging result is extracted. The calculation formula is as follows: ； ； in: For imaging results The corresponding analytical signal; For Hilbert transform operators; The amplitude envelope matrix of the imaging result; Modular operation for complex numbers; ②The amplitude envelope matrix Normalization is performed: ； in: This is the normalized amplitude envelope matrix; The amplitude envelope matrix The minimum value in; The amplitude envelope matrix The maximum value in.

6. The method for constructing an initial model for ground-penetrating radar imaging according to claim 5, characterized in that, In step S3.2, the radial energy concentration REC is calculated as follows: ； in: For scanning speed Radial energy concentration of the imaging results; For the imaging results in the horizontal direction The number of sampling points on; For the imaging results in the vertical depth direction The number of sampling points on; The in-phase axis sharpness index (HSI) is calculated as follows: ① Calculate the lateral gradient: ； in: This represents the gradient value of the normalized amplitude envelope matrix in the horizontal direction. Horizontal direction The sampling interval on; Indicates spatial location The normalized amplitude envelope matrix at the point; Indicates spatial location The normalized amplitude envelope matrix at the point; ② Calculate the sharpness index of the in-phase axis: ； in: For scanning speed The in-phase axis sharpness index of the imaging results; for The gradient mean; for The gradient standard deviation; The focus quality factor (FQF) is calculated as follows: ① Normalize the sharpness index of the in-phase axis: ； in: For scanning speed The normalized in-phase axis sharpness index; For velocity set The minimum value of the sharpness index of all in-phase axes; For velocity set The maximum value of the sharpness index of all in-phase axes; ② Calculate the focusing quality factor FQF: ； in: For scanning speed The focus quality factor of the imaging results.

7. The method for constructing an initial model for ground-penetrating radar imaging according to claim 6, characterized in that, Step S3.3 specifically involves: using the Focused Quality Factor (FQF) as the core evaluation criterion, and evaluating the velocity set. The FQF values ​​corresponding to all scanning speeds are analyzed, and the optimal speed is automatically selected. The imaging result corresponding to the optimal speed is the one with the best focusing effect across the entire speed range, i.e., the optimal focusing imaging result. The rules for selecting the optimal speed are as follows: ①If the FQF value is in the velocity set If there exists a single global maximum value, then the scanning speed corresponding to that maximum value is the optimal speed, denoted as . The formula is: ； in: To achieve the optimal speed; This represents the operator that finds the independent variable that maximizes the function; ②If the FQF value is in the velocity set If there are multiple local maxima, then the scanning speed corresponding to each local maximum is taken as the optimal speed, denoted as the optimal speed set. , The number of local maxima; ； in: This represents the optimal speed set.

8. The method for constructing an initial model for ground-penetrating radar imaging according to claim 7, characterized in that, In step S4, based on the optimal focused imaging results, the structural features of the subsurface medium are extracted using the Canny edge detection algorithm: the normalized amplitude envelope matrix corresponding to the optimal velocity is... As input, Gaussian smoothing, gradient calculation, non-maximum suppression, and double thresholding are performed sequentially to extract the structural edge set of the subsurface medium; as detailed below: ① The normalized amplitude envelope matrix of the input Gaussian smoothing is performed using the following formula: ； in: This is the normalized amplitude envelope matrix after Gaussian smoothing; This is a two-dimensional convolution operation; It is a two-dimensional Gaussian kernel. The standard deviation of the Gaussian kernel; ② Calculate the horizontal and vertical gradients of the normalized amplitude envelope matrix after Gaussian smoothing, and then fuse the horizontal and vertical gradients to obtain the gradient magnitude and gradient direction. The calculation formula is as follows: ； ； ； ； in: This represents the horizontal gradient value of the normalized amplitude envelope matrix after Gaussian smoothing; This represents the vertical gradient value of the normalized amplitude envelope matrix after Gaussian smoothing. Vertical depth direction z The sampling interval is a fixed constant. The gradient magnitude represents the intensity of the edge. The gradient direction represents the orientation of the edge. Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point; Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point; Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point; Indicates spatial location The Gaussian smoothed normalized amplitude envelope matrix at the point; ③ Gradient magnitude along the gradient direction By performing point-by-point judgment, only local maxima along the gradient direction are retained, and gradient values ​​at non-edge points are suppressed to obtain the gradient magnitude after non-maximum suppression. ④ Set a high threshold and low threshold , > The gradient magnitude after non-maximum suppression is compared with the high threshold. and low threshold Comparison: Gradient magnitude greater than The points are strong edge points and are directly identified as the structural edge of the underground medium; Gradient magnitude is between and The points between them are weak edge points. Only when a point is connected to a strong edge point is it determined to be a structural edge of the underground medium. Gradient magnitude less than Points that are not edge points are directly removed; Through the above operations, the final set of structural edges of the subsurface medium is obtained, denoted as... ,in This indicates that the location is a structural boundary. This indicates that the location is a non-structural boundary.

9. The method for constructing an initial model for ground-penetrating radar imaging according to claim 8, characterized in that, In step S4, the physical mapping relationship between electromagnetic wave velocity and relative permittivity is as follows: ； The formula for calculating the relative permittivity, obtained by transforming the above equation, is as follows: ； in: This represents the actual propagation speed of electromagnetic waves in underground media; The speed at which electromagnetic waves propagate in a vacuum; The relative permittivity of the underground medium; Optimal speed Substituting into the above equation, we obtain the relative permittivity value corresponding to the optimal velocity, denoted as . : ； If an optimal velocity set exists For each optimal speed , Calculate the relative permittivity corresponding to the optimal velocity. .

10. The method for constructing an initial model for ground-penetrating radar imaging according to claim 9, characterized in that, In step S4, based on the structural edge set The underground space is divided into multiple independent medium regions, and then the relative permittivity is considered. The relative permittivity of each region was assigned and the background was filled, thus constructing an initial model of the relative permittivity of the subsurface medium. That is, the initial model for ground-penetrating radar imaging, specifically: From the set of structural edges The enclosed region is denoted as Each closed region represents an independent underground medium. For each closed region Assign it the relative permittivity obtained by the corresponding optimal velocity conversion. If multiple optimal velocities exist, assign corresponding velocities to different regions based on the target of the offset focusing. ; For the background region outside the closed region, the relative permittivity value corresponding to the scan rate of the relatively stable portion of the focus quality factor (FQF) is selected as the background fill value. ; The final initial model of the relative permittivity of the subsurface medium for: 。

Images