Fine imaging method, medium and equipment for ground penetrating radar fusion inversion

Abstract

The present invention relates to the field of ground penetrating radar data processing technology, and specifically to a fine imaging method, medium and equipment for ground penetrating radar fusion inversion, including: using ground penetrating radar to scan a certain detection area to obtain the original GPR data profile; analyzing the wavefront continuity of the GPR data profile, introducing a correction factor to non-destructively separate the GPR data profile to obtain a reflection data profile and a diffraction data profile; constructing an inversion initial model based on the hyperbolic characteristics of the diffraction data profile; constructing an inversion objective function based on the reflection data profile and the diffraction data profile; using the inversion initial model as the inversion initial solution, and using the L-BFGS algorithm to solve the optimization problem of the inversion objective function to obtain a diffraction wave-reflection wave fusion imaging result. The method can effectively improve the imaging effect and accuracy of the ground penetrating radar processing results for tiny diffraction anomalies, optimize the inversion imaging performance, and promote the practical application of ground penetrating radar fine imaging and full waveform inversion imaging.

Description

Technical Field

[0001] The present invention relates to the technical field of ground penetrating radar data processing, and in particular to a fine imaging method, medium and equipment for ground penetrating radar fusion inversion. Background Art

[0002] In the application scenarios of ground penetrating radar, microstructures and layered structures are common observation targets, such as structural damage under the surface of paved roads, steel bar layout in concrete walls, etc. In traditional ground penetrating radar data processing and imaging, there is a lack of a dedicated separation and imaging process for key targets, resulting in weak ground penetrating radar's ability to identify abnormal body targets with crosstalk under complex structural conditions.

[0003] Ground Penetrating Radar (GPR) is an important non-destructive testing and detection method, which is widely used in shallow geophysical detection fields such as engineering detection, environmental exploration and archaeological investigation. It uses the reflection, diffraction and other wave field propagation characteristics of high-frequency electromagnetic waves in underground media to study the parameter distribution of underground media. GPR directly obtains the two-way travel time, amplitude, frequency, phase and other information of electromagnetic waves in the medium. The imaging processing of the full-width ground penetrating radar data is the prerequisite for a comprehensive and objective understanding of the geometric size and physical parameters of the underground medium.

[0004] At present, the imaging processing methods of GPR data include relatively mature inversion methods including travel time tomography inversion, wavefront tomography inversion (WTI) with attenuation coefficient and amplitude as parameters or full amplitude information of wave field as parameters, and full wavefield inversion (FWI). However, these inversion imaging methods are only for GPR reflection field, especially in On-Ground coupling mode. Due to the relative lack of GPR data sets from underground target bodies (compared with Cross-hole mode), the pathological and nonlinear characteristics of GPR reflection field inversion problems are particularly prominent, and the local imaging resolution of underground target bodies needs to be improved urgently. In addition, the GPR inversion for traditional full wavefield information lacks the ability to reconstruct the characteristics of diffracted bodies. Summary of the invention

[0005] The present invention aims to provide a fine imaging method of ground penetrating radar fusion inversion that effectively improves the imaging effect and accuracy of ground penetrating radar processing results for tiny diffraction anomalies. The specific technical scheme is as follows:

[0006] The present invention provides a ground penetrating radar fusion inversion fine imaging method, comprising the following steps:

[0007] S1: Use ground penetrating radar to scan a certain detection area to obtain the original GPR data profile, where: the GPR data profile is composed of a combination of multiple single-channel data obtained by continuous acquisition;

[0008] S2: Analyze the wavefront continuity of the GPR data profile based on the multi-channel coherence properties of the GPR data profile, and introduce a correction factor to use the coherence subtraction method to losslessly separate the GPR data profile to obtain the reflection data profile and the diffraction data profile;

[0009] S3: Constructing the initial inversion model based on the hyperbolic characteristics of the diffraction data profile;

[0010] S4: GPR diffraction wave fusion inversion imaging, specifically:

[0011] S4.1. constructing an inversion objective function based on the reflection data profile and the diffraction data profile;

[0012] S4.2. The inversion initial model is used as the inversion initial solution, and the L-BFGS algorithm is used to solve the optimization problem of the inversion objective function. When the number of inversion iterations in the L-BFGS algorithm is greater than or equal to the first set threshold or the inversion objective function value is less than the second set threshold, the diffraction wave-reflection wave fusion imaging result that can completely invert the underground fine electrical structure of the detection area is obtained.

[0013] Optionally, the S2 includes:

[0014] S2.1. Using the coherence between the signal strength of each channel of the GPR data profile and other channels within a given range as an indicator for analyzing the continuity of the wavefront between channels, the coherence value of each data point on the GPR data profile within the given range is calculated, wherein: the given range includes the sampling time point range and the recording channel range;

[0015] S2.2, introduce the coherent stacking time correction factor to stack the GPR data profile, and obtain the signal strength of the data profile at each data point after stacking, which is the reflection data profile, specifically:

[0016] S2.2.1, construct the coherent superposition time correction factor;

[0017] S2.2.2. Introduce the coherent stacking time correction factor to perform data stacking on the GPR data profile to obtain a complete reflection data profile;

[0018] S2.3, introduce the coherent subtraction amplitude correction factor and the coherent subtraction time shift correction factor to subtract the GPR data profile from the reflection data profile to obtain a completely separated reflection data profile and diffraction data profile, specifically:

[0019] S2.3.1, constructing a coherent subtraction amplitude correction factor and a coherent subtraction time shift correction factor;

[0020] S2.3.2. Based on the coherent subtraction amplitude correction factor and the coherent subtraction time shift correction factor, the coherent subtraction method is used to losslessly separate the GPR data profile to obtain completely separated reflection data profile and diffraction data profile.

[0021] Optionally, in S2.1, the specific formula of the coherence value is as follows:

[0022] ;

[0023] in: is the coordinate of the data point on the GPR data profile, is the track coordinate, is the sampling time point coordinate; For GPR data profile Signal strength at For data points coherence values ​​within a given range; is the sampling time point range; The range of the recording channel.

[0024] Optionally, the S2.2.1 specifically includes:

[0025] ①, Assume that the coherent superposition time correction factor is , after introducing random initialization Perform data stacking on the GPR data profile to obtain the stacked data profile , the specific formula is as follows:

[0026] ;

[0027] in: It is a data point The corresponding correction value is an unknown parameter; The data section after stacking Data points on

[0028] ② Calculate the data profile after stacking The coherence value of each data point is calculated, and the optimization problem of maximizing the coherence is constructed. The specific formula is as follows:

[0029] ;

[0030] ③. Use the particle swarm algorithm to solve the optimization problem of maximizing coherence and obtain the optimal coherent stacking time correction factor corresponding to each point on the stacked data profile. , where: The convergence condition of the particle swarm algorithm is: is greater than the third set threshold, or the number of iterations is greater than the fourth set threshold.

[0031] Optionally, in S2.3.1, constructing the coherent subtraction amplitude correction factor and the coherent subtraction time shift correction factor specifically includes:

[0032] ①, Assume the amplitude correction factor of coherent subtraction is The time shift correction factor of the coherent subtraction method is , after introducing random initialization and , subtract the GPR data profile from the reflection data profile to obtain the diffraction wave field , the specific formula is as follows:

[0033] ;

[0034] in: and They are points The corresponding correction value is an unknown parameter;

[0035] ② Calculate the diffraction wave field The coherence value of each data point in , and the optimization problem of maximizing the coherence ratio is constructed. The specific formula is as follows:

[0036] ;

[0037] in: is the quantity to be solved; The coherence of the diffraction wave data profile is obtained by coherence subtraction. It means when =1, = 0, the coherence of the diffraction wave data profile obtained after random initialization and The value is equivalent to not making a correction to the subtraction;

[0038] ③. The particle swarm algorithm is used to solve this optimization problem, and finally the coherent subtraction amplitude correction factor corresponding to each data point is obtained as follows: The time shift correction factor of the coherent subtraction method is , where: The convergence condition of the particle swarm algorithm is: is greater than the third set threshold, or the number of iterations is greater than the fourth set threshold.

[0039] Optionally, the first threshold is set to 1000-2000, and the second threshold is set to ; The third threshold value is set to 0.8~0.95, and the fourth threshold value is set to 200~400.

[0040] Optionally, the S3 includes:

[0041] S3.1. obtaining geometric information of the hyperbolic feature by circling the hyperbolic feature on the diffraction data profile;

[0042] S3.2, obtaining velocity distribution based on geometric information of the hyperbolic feature, and performing two-dimensional interpolation based on velocity points corresponding to a plurality of hyperbolic features to obtain velocity distribution on a diffraction data profile;

[0043] S3.3. Convert the velocity distribution into electrical structure distribution, and perform two-dimensional interpolation based on multiple hyperbolic characteristic vertices to obtain the electrical structure distribution on the diffraction data profile, that is, invert the initial model.

[0044] Optionally, the specific expression of the inversion objective function in S4.1 is as follows:

[0045] ;

[0046] in: is the diffraction field inversion weight, is the reflection field inversion weight, , is the model parameter vector to be inverted, ,Pick ; is the diffraction field inversion objective function; is the objective function of reflection field inversion, and:

[0047] ;

[0048] in: The observed diffraction wave data vector constructed from the diffraction data profile obtained in S2; The observed reflection wave data vector constructed from the reflection data profile obtained in S2; is the diffraction wave structure model vector to be solved, is the reflection wave structure model vector to be solved; represent In the model response relationship The forward response under represent In the model response relationship The forward response under .

[0049] The present invention also provides a readable storage medium having computer program instructions stored thereon, and when the computer program instructions are executed by a processor, the fine imaging method of ground penetrating radar fusion inversion as described above is implemented.

[0050] The present invention also provides an electronic device, comprising: at least one processor, at least one memory and computer program instructions stored in the memory, when the computer program instructions are executed by the processor, the fine imaging method of ground penetrating radar fusion inversion as described above is implemented.

[0051] Based on the multi-channel coherent attribute properties of ground penetrating radar data, the present invention performs coherent continuity analysis on the collected data, and acquires the diffraction wave field in the collected data in a non-destructive and purely data-driven manner; performs diffraction wave hyperbola feature recognition based on the extracted diffraction wave field, and preliminarily infers the underground velocity structure information by using the corresponding relationship between the geometric attributes of the response wavefront of the anomaly body in the ground penetrating radar data and the dynamic attributes of underground electromagnetic wave propagation, and constructs an initial velocity model; uses the constructed velocity model as the initial model to perform fusion constraint inversion on the separated diffraction wave field and the original collected full wave field. This method can effectively improve the imaging effect and accuracy of the ground penetrating radar processing results for tiny diffraction anomalies, optimize the inversion imaging performance to a certain extent, and promote the practical application of ground penetrating radar fine imaging and full waveform inversion imaging.

[0052] In addition to the above-described purposes, features and advantages, the present invention has other purposes, features and advantages. The present invention will be further described in detail with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] The drawings constituting a part of this application are used to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the drawings:

[0054] Figure 1 It is a technical roadmap of the fine imaging method of ground penetrating radar fusion inversion in an embodiment of the present invention;

[0055] Figure 2 is a schematic diagram of a B-scan collected by a ground penetrating radar in an embodiment of the present invention;

[0056] Figure 3 is a schematic diagram of the superposition of uncorrected GPR data profiles in an embodiment of the present invention;

[0057] Figure 4 is a schematic diagram of superimposed GPR data profiles after correction in an embodiment of the present invention;

[0058] Figure 5 is the technical roadmap of S2 in the embodiment of the present invention;

[0059] Figure 6 is a schematic diagram of the morphology of electromagnetic wave echo scattering of a ground penetrating radar in an embodiment of the present invention;

[0060] Figure 7is a schematic diagram of the geometric shape of the hyperbolic feature on the diffraction data profile in S3.1 in an embodiment of the present invention;

[0061] Figure 8 is the technical roadmap of S3 in the embodiment of the present invention, Figure 8 (a) is a schematic diagram of the diffraction data profile. Figure 8 (b) is a schematic diagram of the effect of the extracted hyperbolic features. Figure 8 (c) is a schematic diagram of the obtained velocity distribution effect. Figure 8 (d) is a schematic diagram of the effect of electrical structure distribution;

[0062] Fig. 9 (a) is the effect diagram of GPR data profile. Fig. 9 (b) is the separation effect diagram of the existing algorithm 1 (diffraction wave separation algorithm based on plane wave deconstruction filtering). Fig. 9 (c) is the separation effect diagram of the existing algorithm 2 (diffraction wave separation algorithm based on common reflection surface element superposition wavefront attribute analysis). Fig. 9 (d) is a separation effect diagram of S2 in an embodiment of the present invention;

[0063] Fig.10 is a comparison diagram of the diffraction wave-reflection wave fusion imaging results of the embodiment of the present invention and the traditional method, Fig.10 (a) is a schematic diagram of a uniform initial model constructed by the traditional method; Fig.10 (b) is a schematic diagram of the initial inversion model in this embodiment; Fig.10 (c) is a schematic diagram of the full wave field inversion effect of the traditional method; Fig.10 (d) is a schematic diagram of the fusion inversion effect in this embodiment; Fig.10 (e) is a schematic diagram of the effect of the traditional method of full wave field on the separate diffraction field imaging; Fig.10 (f) is a schematic diagram of the effect of fusion inversion on separate imaging of the diffraction fields in this embodiment;

[0064] Fig.11 It is a schematic diagram of the process of S4.2 in an embodiment of the present invention.

[0065] Description of Figure Numbers:

[0066] 1 Antenna direction, 2 Electromagnetic wave path, 3 Target object, 4 Air layer, 5 Dielectric layer, 6 Scattering curve. DETAILED DESCRIPTION

[0067] In order to make the purpose, features and advantages of the present invention more obvious and easy to understand, the specific implementation methods of the present invention are described in detail below in conjunction with the accompanying drawings. Several embodiments of the present invention are shown in the accompanying drawings. However, the present invention can be implemented in many different forms and is not limited to the embodiments described herein.

[0068] The embodiments of the present invention are described in detail below with reference to the accompanying drawings, but the present invention can be implemented in many different ways as defined and covered by the claims.

[0069] In one embodiment, if Figure 1 As shown, a fine imaging method of ground penetrating radar fusion inversion includes the following steps:

[0070] The GPR data profile is composed of multiple single-channel data acquired continuously. The data points on the profile can be indexed and located by channel number-time (or sampling point). Indicates. Since electromagnetic waves scatter and reflect in underground abnormal structures when GPR collects data, the recorded signals between each channel on the profile have continuous characteristics, indicating the response of the underground structure, such as Figure 2 When there are abnormal scatterers underground, the corresponding diffraction response signal on the GPR data section will present a characteristic hyperbolic waveform. Compared with the continuous reflection wavefront presented by the horizontal strata and other reflective structures on the GPR data section, this type of hyperbolic waveform has the discontinuous characteristics of multi-channel data. Based on this, the hyperbolic characteristic signal corresponding to the scatterer can be identified by analyzing the continuity of the multi-channel data on the GPR data section, thereby separating the diffraction field signal corresponding to the scatterer.

[0071] When there are abnormal scatterers underground, the corresponding diffraction response signal on the GPR data profile will present a characteristic hyperbolic waveform. Compared with the continuous reflection wavefront presented by reflection structures such as horizontal strata on the GPR data profile, this type of hyperbolic waveform has the discontinuous characteristics of multi-channel data. Based on this, the hyperbolic characteristic signal corresponding to the scatterer can be identified by analyzing the continuity of multi-channel data on the GPR data profile, thereby separating the diffraction field signal corresponding to the scatterer.

[0072] 1. Continuity analysis of GPR data coherence attributes and separation of diffraction fields.

[0073] S1: Use ground penetrating radar to scan a certain detection area to obtain the original GPR data profile, where: the GPR data profile is composed of a combination of multiple single-channel data obtained by continuous acquisition;

[0074] S2: Analyze the wavefront continuity of the GPR data profile based on the multi-channel coherence properties of the GPR data profile, and introduce a correction factor to use the coherence subtraction method to losslessly separate the GPR data profile to obtain the reflection data profile and the diffraction data profile;

[0075] See also Figures 2 to 5 , S2 includes:

[0076] In order to accurately analyze the continuity between each trace on the GPR data profile, this embodiment proposes to quantify it using a coherence index. For simplicity, the GPR data profile is regarded as a two-dimensional array, which is expressed as: , is the point on the GPR data profile The signal strength at the location (i.e., the field strength amplitude recorded by GPR at that location). The similarity between the signal strength of each channel and other channels within a given range is used as an indicator to evaluate the continuity of signals between multiple channels.

[0077] S2.1. Use the coherence between the data of each channel of the GPR data profile and the signal strength of other channels within a given range as an indicator for analyzing the continuity of the wavefront between channels, and calculate the coherence value of each data point on the GPR data profile within a given range, where: the given range includes the sampling time point range and the recording channel range.

[0078] In S2.1, the specific formula for the coherence value is as follows:

[0079] ;

[0080] in: is the coordinate of the data point on the GPR data profile, is the track coordinate, is the sampling time point coordinate; For GPR data profile Signal strength at For data points coherence values ​​within a given range; is the sampling time point range; The range of the recording channel.

[0081] According to the dynamic theory of high-frequency electromagnetic wave propagation in the medium, the signal at each point on the GPR data can be regarded as the linear superposition of the reflection signal and the diffraction signal. Therefore, on the premise of obtaining the reflection signal on the profile, the diffraction signal can be obtained by subtracting the GPR data profile from the reflection profile. One way to obtain the reflection signal is to superimpose the GPR data profile, such as Figure 3 As shown, for each point on the cross section :

[0082] ;

[0083] in: The section after superposition exist The signal strength, Same meaning as above, it is the size of the recording channel range, Figure 2 Have consistent meaning.

[0084] The stacking method defined by this formula can only ensure that the completely horizontal reflected wavefront event axis is reconstructed. However, in actual production scenarios, the reflected signal on the GPR data profile often has a certain inclination. In this case, stacking using this formula cannot completely reconstruct the reflected wave field. Therefore, it is necessary to modify the stacking algorithm, such as Figure 4 shown.

[0085] S2.2, introduce the coherent stacking time correction factor to stack the GPR data profile, and obtain the signal strength of the data profile at each data point after stacking, which is the reflection data profile, specifically:

[0086] S2.2.1, construct the coherent superposition time correction factor;

[0087] S2.2.1 specifically includes:

[0088] As can be seen from S2.1, the coherence between the wavefront events of continuous reflected waves is relatively large, so the coherence can be used as an indicator to determine whether the superposition is completely included in the non-zero tilt wavefront event.

[0089] ①, Assume that the coherent superposition time correction factor is , after introducing random initialization Perform data stacking on the GPR data profile to obtain the stacked data profile , the specific formula is as follows:

[0090] ;

[0091] in: It is a data point The corresponding correction value is an unknown parameter; The data section after stacking Data points on

[0092] ② Calculate the data profile after stacking The coherence value of each data point is calculated, and the optimization problem of maximizing the coherence is constructed. The specific formula is as follows:

[0093] ;

[0094] ③. Use the particle swarm algorithm to solve the optimization problem of maximizing coherence and obtain the optimal coherent stacking time correction factor corresponding to each point on the stacked data profile. , where: The convergence condition of the particle swarm algorithm is: is greater than the third set threshold, or the number of iterations is greater than the fourth set threshold.

[0095] S2.2.2, introduce the coherent stacking time correction factor to perform data stacking on the GPR data profile to obtain a complete reflection data profile; Figure 4 As shown in Figure 1, the stacked event axis is no longer limited to horizontal signals, and has a good reflection of the non-zero dip wavefront on the profile. At this time, the stacking result can be regarded as an accurate reflection profile.

[0096] S2.3, introduce the coherent subtraction amplitude correction factor and the coherent subtraction time shift correction factor to subtract the GPR data profile from the reflection data profile to obtain a completely separated reflection data profile and diffraction data profile, specifically:

[0097] As described in S2.2, the diffraction wave distribution can be approximately obtained by subtracting the GPR data profile from the reflection profile. Without correction, it can be expressed as:

[0098] ;

[0099] in: is the diffraction wave signal, and The meaning is consistent with the previous text. Due to the difference in the principles of diffraction and reflection of high-frequency electromagnetic waves, the intensity of the reflected signal obtained by S2.2 will be significantly greater than the diffraction signal. If the direct subtraction method shown in this formula is used, the obtained diffraction wave signal intensity distribution will be inaccurate; in addition, the extension difference between the diffraction wave and the reflected wave front signal in a single-channel record causes the direct subtraction method shown in this formula to cause deviations in the sampling point direction of the extraction result. Therefore, this embodiment constructs a coherent subtraction amplitude correction factor and coherent subtraction time shift correction factor .

[0100] S2.3.1, constructing a coherent subtraction amplitude correction factor and a coherent subtraction time shift correction factor;

[0101] In S2.3.1, constructing the coherent subtraction amplitude correction factor and the coherent subtraction time shift correction factor specifically includes:

[0102] ①, Assume the amplitude correction factor of coherent subtraction is The time shift correction factor of the coherent subtraction method is , after introducing random initialization and , subtract the GPR data profile from the reflection data profile to obtain the diffraction wave field , the specific formula is as follows:

[0103] ;

[0104] in: and They are points The corresponding correction value is an unknown parameter;

[0105] ② Calculate the diffraction wave field The coherence value of each data point in , and the optimization problem of maximizing the coherence ratio is constructed. The specific formula is as follows:

[0106] ;

[0107] in: is the quantity to be solved; The coherence of the diffraction wave data profile is obtained by coherence subtraction. It means when =1, = 0, the coherence of the diffraction wave data profile obtained after random initialization and The value is equivalent to not making a correction to the subtraction;

[0108] This embodiment uses a particle swarm algorithm to solve this optimization problem and finally obtains the correction factor corresponding to each point on the profile. . Then, by performing coherent subtraction, a completely separated diffraction profile can be obtained.

[0109] ③. The particle swarm algorithm is used to solve this optimization problem, and finally the coherent subtraction amplitude correction factor corresponding to each data point is obtained as follows: The time shift correction factor of the coherent subtraction method is , where: The convergence condition of the particle swarm algorithm is: is greater than the third set threshold, or the number of iterations is greater than the fourth set threshold.

[0110] In this embodiment, the third set threshold is 0.9, and the fourth set threshold is 300.

[0111] S2.3.2. Based on the coherent subtraction amplitude correction factor and the coherent subtraction time shift correction factor, the coherent subtraction method is used to losslessly separate the GPR data profile to obtain completely separated reflection data profile and diffraction data profile.

[0112] S3: Constructing the initial inversion model based on the hyperbolic characteristics of the diffraction data profile;

[0113] Inversion of GPR data can make more accurate judgments on the electrical characteristics and morphology of the detection area. Accurate inversion requires accurate initial model definition. In GPR inversion, the initial model is controlled by the electrical parameters of the medium (mainly conductivity). , dielectric constant and magnetic permeability ). According to Maxwell's equations, different dielectric electrical parameters and different electromagnetic wave frequencies lead to different propagation speeds of electromagnetic waves in the medium. Therefore, this embodiment proposes a method for constructing an initial velocity model using the geometric and dynamic characteristics of diffraction wave signals. The scatterer response in GPR recorded data is mostly in the form of a hyperbola,

[0114] like Figure 6 As shown, is the center projection position of target body 3, is the GPR antenna position (i.e., observation point position); GPR antenna direction 1 is determined by the antenna position Point to the center projection position of target body 3 The line between the GPR antenna position and the target 3 is the electromagnetic wave path 2; is the propagation speed of electromagnetic waves in the dielectric layer 5, satisfying ,in: is the propagation speed of electromagnetic waves in vacuum, and are the relative permittivity and relative magnetic permeability of the dielectric layer 5 respectively; and are the dielectric constant and magnetic permeability of the air layer 5, respectively; is the time required for the electromagnetic wave to propagate from the target 3 to the projection position directly above it, It is the time required for the electromagnetic wave to propagate from the target to the observation point.

[0115] S3 includes:

[0116] Define the target position as , is the depth position in the actual physical scene, ; So combined Figure 7 The above geometric relationship can be used to obtain the standard form of the hyperbola equation:

[0117] ;

[0118] Therefore, by identifying the geometric features of the diffraction wave hyperbola shape, the corresponding dielectric electrical parameters can be determined. This embodiment uses curvature fitting to identify it. For any hyperbola, its vertex position can be known from its geometric properties. The curvature is ,like Figure 7 shown.

[0119] In this embodiment, the diffraction wave profile is extracted The characteristic hyperbola on the surface is circled to obtain its geometric information and further obtain the velocity distribution. On this basis, the velocity points corresponding to multiple characteristic hyperbolas are Perform two-dimensional interpolation to obtain the velocity distribution on the profile; then use the formula (where relative magnetic permeability In most scenarios, it can be regarded as, that is, no magnetic anomaly distribution) converting the velocity distribution into electrical structure distribution , similarly, two-dimensional interpolation is performed on several characteristic hyperbola vertices to obtain the electrical structure distribution on the profile as an accurate inversion initial model, such as Figure 8 shown.

[0120] S3.1. obtaining geometric information of the hyperbolic feature by circling the hyperbolic feature on the diffraction data profile;

[0121] S3.2, obtaining velocity distribution based on geometric information of the hyperbolic feature, and performing two-dimensional interpolation based on velocity points corresponding to a plurality of hyperbolic features to obtain velocity distribution on a diffraction data profile;

[0122] S3.3. Convert the velocity distribution into electrical structure distribution, and perform two-dimensional interpolation based on multiple hyperbolic characteristic vertices to obtain the electrical structure distribution on the diffraction data profile, that is, invert the initial model.

[0123] S4: GPR diffraction wave fusion inversion imaging, specifically:

[0124] Inversion imaging of GPR data requires the construction of an inversion objective function. The traditional inversion problem is expressed as:

[0125] ;

[0126] in: is the observed data vector; is the model parameter vector to be solved by inversion; and The solution dimension of the inversion problem is defined, which is related to the number of parameters to be solved and the amount of observation data; Represents the model parameter vector in the model response relationship Forward response under ; is the random noise in the observed data and can be ignored.

[0127] The inversion objective function is constructed based on the traditional inversion problem:

[0128] ;

[0129] in: is the optimization objective function; It can be roughly regarded as the GPR data profile mentioned above.

[0130] Obviously, the optimization problem given by the inversion objective function only considers the GPR data profile record where the reflected wave and the diffraction field are not separated, and it is difficult to perform fine imaging of the separated diffraction wave information given in S2.3 above. Therefore, this embodiment proposes a fusion inversion imaging objective function, which integrates the reflection field data obtained by coherent superposition and the diffraction field data obtained by separation to perform fusion inversion.

[0131] S4.1. constructing an inversion objective function based on the reflection data profile and the diffraction data profile;

[0132] The specific expression of the inversion objective function in S4.1 is as follows:

[0133] ;

[0134] in: is the diffraction field inversion weight, is the reflection field inversion weight, , is the model parameter vector to be inverted, ,Pick ; is the diffraction field inversion objective function; is the objective function of reflection field inversion, and:

[0135] in: The observed diffraction wave data vector constructed from the diffraction data profile obtained in S2; The observed reflection wave data vector constructed from the reflection data profile obtained in S2; is the diffraction wave structure model vector to be solved, is the reflection wave structure model vector to be solved; represent In the model response relationship The forward response under represent In the model response relationship The forward response under . The model vector to be solved is is the electrical structure dual parameter distribution, that is, the relative dielectric constant and conductivity .

[0136] S4.2. The initial inversion model is used as the initial inversion solution, and the L-BFGS algorithm is used to solve the optimization problem of the inversion objective function. When the number of inversion iterations in the L-BFGS algorithm is greater than or equal to the first set threshold or the inversion objective function value is less than the second set threshold, the diffraction wave-reflection wave fusion imaging result that can completely invert the underground fine electrical structure of the detection area is obtained, that is, the current model finally obtained ,like Fig.10 and Fig.11 As shown, Fig.10 The dotted line in the figure marks the real structure shape and position. It can be seen that compared with the traditional method, the imaging effect and accuracy in this embodiment are better. In this embodiment, the first threshold value is set to 1000~2000, and the second threshold value is set to .

[0137] like Fig. 9 As shown, Fig. 9 The method in this embodiment represented by (d) is compared with Fig. 9 The diffraction wave separation algorithm based on plane wave deconstruction filtering given in (b) of FIG. 11A can better preserve the response information of the diffraction wave vertex position and achieve a nearly lossless separation effect. Fig. 9 The wavefront property analysis diffraction wave separation algorithm based on the superposition of common reflection surface elements given in (c) achieves a relatively close separation effect.

[0138] Based on the multi-channel coherent attribute properties of the ground penetrating radar data, this embodiment performs coherent continuity analysis on the collected data, and obtains the diffraction wave field in the collected data in a non-destructive and purely data-driven manner; based on the extracted diffraction wave field, the hyperbolic feature of the diffraction wave is identified, and the corresponding relationship between the geometric attributes of the response wavefront of the anomaly in the ground penetrating radar data and the dynamic attributes of the underground electromagnetic wave propagation is used to preliminarily infer the underground velocity structure information and construct an initial velocity model; the constructed velocity model is used as the initial model to perform fusion constraint inversion on the separated diffraction wave field and the original collected full wave field. This method can effectively improve the imaging effect and accuracy of the ground penetrating radar processing results for tiny diffraction anomalies, optimize the inversion imaging performance to a certain extent, and promote the practical application of ground penetrating radar fine imaging and full waveform inversion imaging.

[0139] This embodiment also provides a readable storage medium on which computer program instructions are stored. When the computer program instructions are executed by a processor, the fine imaging method of ground penetrating radar fusion inversion as described above is implemented.

[0140] It should be noted that the device embodiments described above are merely schematic, wherein the units described as separate components may or may not be physically separated, and the components displayed as units may or may not be physical units, that is, they may be located in one place, or they may be distributed on multiple network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the scheme of this embodiment. In addition, in the accompanying drawings of the device embodiments provided by the present invention, the connection relationship between the modules indicates that there is a communication connection between them, which may be specifically implemented as one or more communication buses or signal lines. A person of ordinary skill in the art may understand and implement it without paying any creative effort.

[0141] This embodiment also includes an electronic device, including: at least one processor, at least one memory, and computer program instructions stored in the memory. When the computer program instructions are executed by the processor, the fine imaging method of ground penetrating radar fusion inversion as described above is implemented.

[0142] Exemplarily, the computer program may be divided into one or more modules / units, which are stored in the memory and executed by the processor to implement the present invention. The one or more modules / units may be a series of computer program instruction segments capable of implementing specific functions, which are used to describe the execution process of the computer program in the electronic device.

[0143] The electronic device may be a computing device such as a mobile phone, a desktop computer, a notebook, a PDA, a cloud server, etc. The electronic device may include, but is not limited to, a processor and a memory. For example, the electronic device may also include an input / output device, a network access device, a bus, etc.

[0144] The processor may be a central processing unit (CPU), other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or any conventional processor, etc. The processor is the control center of the electronic device, and various interfaces and lines are used to connect various parts of the entire electronic device.

[0145] The memory can be used to store the computer program and / or module, and the processor implements the computer program by running or executing the computer program and / or module stored in the memory, and calling the data stored in the memory. The memory can mainly include a program storage area and a data storage area, wherein the program storage area can store an operating system, an application required for at least one function (such as a sound playback function, an image playback function, etc.), etc.; the data storage area can store data created according to the use of the mobile phone (such as audio data, a phone book, etc.), etc. In addition, the memory can include a high-speed random access memory, and can also include a non-volatile memory, such as a hard disk, a memory, a plug-in hard disk, a smart memory card (Smart Media Card, SMC), a secure digital (SecureDigital, SD) card, a flash card (Flash Card), at least one disk storage device, a flash memory device, or other volatile solid-state storage devices.

[0146] Wherein, if the module / unit integrated in the electronic device is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on such an understanding, the present invention implements all or part of the processes in the above-mentioned embodiment method, and can also be completed by instructing the relevant hardware through a computer program. The computer program can be stored in a readable storage medium, and the computer program can implement the steps of the above-mentioned various method embodiments when executed by the processor. Wherein, the computer program includes computer program code, and the computer program code can be in source code form, object code form, executable file or some intermediate form, etc. The computer-readable medium may include: any entity or device capable of carrying the computer program code, recording medium, U disk, mobile hard disk, disk, optical disk, computer memory, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), electric carrier signal, telecommunication signal and software distribution medium, etc.

[0147] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A fine imaging method of ground penetrating radar fusion inversion, characterized in that: The following steps are involved: S1: Use ground penetrating radar to scan a certain detection area to obtain the original GPR data profile, where: the GPR data profile is composed of a combination of multiple single-channel data obtained by continuous acquisition; S2: Analyze the wavefront continuity of the GPR data profile based on the multi-channel coherence properties of the GPR data profile, and introduce a correction factor to use the coherence subtraction method to losslessly separate the GPR data profile to obtain the reflection data profile and the diffraction data profile; S3: Constructing the initial inversion model based on the hyperbolic characteristics of the diffraction data profile; S4: GPR diffraction wave fusion inversion imaging, specifically: S4.

1. constructing an inversion objective function based on the reflection data profile and the diffraction data profile; S4.

2. The inversion initial model is used as the inversion initial solution, and the L-BFGS algorithm is used to solve the optimization problem of the inversion objective function. When the number of inversion iterations in the L-BFGS algorithm is greater than or equal to the first set threshold or the inversion objective function value is less than the second set threshold, the diffraction wave-reflection wave fusion imaging result that can completely invert the underground fine electrical structure of the detection area is obtained.

2. The fine imaging method of ground penetrating radar fusion inversion according to claim 1 is characterized in that: The S2 includes: S2.

1. Using the coherence between the signal strength of each channel of the GPR data profile and other channels within a given range as an indicator for analyzing the continuity of the wavefront between channels, the coherence value of each data point on the GPR data profile within the given range is calculated, wherein: the given range includes the sampling time point range and the recording channel range; S2.2, introduce the coherent stacking time correction factor to stack the GPR data profile, and obtain the signal strength of the data profile at each data point after stacking, which is the reflection data profile, specifically: S2.2.1, construct the coherent superposition time correction factor; S2.2.

2. Introduce the coherent stacking time correction factor to perform data stacking on the GPR data profile to obtain a complete reflection data profile; S2.3, introduce the coherent subtraction amplitude correction factor and the coherent subtraction time shift correction factor to subtract the GPR data profile from the reflection data profile to obtain a completely separated reflection data profile and diffraction data profile, specifically: S2.3.1, constructing a coherent subtraction amplitude correction factor and a coherent subtraction time shift correction factor; S2.3.

2. Based on the coherent subtraction amplitude correction factor and the coherent subtraction time shift correction factor, the coherent subtraction method is used to losslessly separate the GPR data profile to obtain completely separated reflection data profile and diffraction data profile.

3. The fine imaging method of ground penetrating radar fusion inversion according to claim 2 is characterized in that: In S2.1, the specific formula of the coherence value is as follows: ； in: is the coordinate of the data point on the GPR data profile, is the track coordinate, is the sampling time point coordinate; For GPR data profile Signal strength at For data points coherence values ​​within a given range; is the sampling time point range; The range of the recording channel.

4. The fine imaging method of ground penetrating radar fusion inversion according to claim 3 is characterized in that: S2.2.1 specifically includes: ①, Assume that the coherent superposition time correction factor is , after introducing random initialization Perform data stacking on the GPR data profile to obtain the stacked data profile , the specific formula is as follows: ； in: It is a data point The corresponding correction value is an unknown parameter; The data section after stacking Data points on ② Calculate the data profile after stacking The coherence value of each data point is calculated, and the optimization problem of maximizing the coherence is constructed. The specific formula is as follows: ； ③. Use the particle swarm algorithm to solve the optimization problem of maximizing coherence and obtain the optimal coherent stacking time correction factor corresponding to each point on the stacked data profile. , where: The convergence condition of the particle swarm algorithm is: is greater than a third set threshold, or the number of iterations is greater than a fourth set threshold.

5. The fine imaging method of ground penetrating radar fusion inversion according to claim 4 is characterized in that: In S2.3.1, constructing the coherent subtraction amplitude correction factor and the coherent subtraction time shift correction factor specifically includes: ①, Assume the amplitude correction factor of coherent subtraction is The time shift correction factor of the coherent subtraction method is , after introducing random initialization and , subtract the GPR data profile from the reflection data profile to obtain the diffraction wave field , the specific formula is as follows: ； in: and They are points The corresponding correction value is an unknown parameter; ② Calculate the diffraction wave field The coherence value of each data point in , and the optimization problem of maximizing the coherence ratio is constructed. The specific formula is as follows: ； in: is the quantity to be solved; The coherence of the diffraction wave data profile is obtained by coherence subtraction. It means when =1, = 0, the coherence of the diffraction wave data profile obtained after random initialization and The value is equivalent to not making a correction to the subtraction; ③. The particle swarm algorithm is used to solve this optimization problem, and finally the coherent subtraction amplitude correction factor corresponding to each data point is obtained as follows: The time shift correction factor of the coherent subtraction method is , where: The convergence condition of the particle swarm algorithm is: is greater than the third set threshold, or the number of iterations is greater than the fourth set threshold.

6. The fine imaging method of ground penetrating radar fusion inversion according to claim 5 is characterized in that: The first threshold is set to 1000~2000, and the second threshold is set to ; The third threshold value is set to 0.8~0.95, and the fourth threshold value is set to 200~400.

7. The fine imaging method of ground penetrating radar fusion inversion according to any one of claims 1 to 6, characterized in that: The S3 includes: S3.

1. obtaining geometric information of the hyperbolic feature by circling the hyperbolic feature on the diffraction data profile; S3.2, obtaining velocity distribution based on geometric information of the hyperbolic feature, and performing two-dimensional interpolation based on velocity points corresponding to a plurality of hyperbolic features to obtain velocity distribution on a diffraction data profile; S3.

3. Convert the velocity distribution into electrical structure distribution, and perform two-dimensional interpolation based on multiple hyperbolic characteristic vertices to obtain the electrical structure distribution on the diffraction data profile, that is, invert the initial model.

8. The method for fine imaging of ground penetrating radar fusion inversion according to any one of claims 1 to 6, characterized in that: The specific expression of the inversion objective function in S4.1 is as follows: ； in: is the diffraction field inversion weight, is the reflection field inversion weight, , is the model parameter vector to be inverted, ,Pick ; is the diffraction field inversion objective function; is the objective function of reflection field inversion, and: ； in: The observed diffraction wave data vector constructed from the diffraction data profile obtained in S2; The observed reflection wave data vector constructed from the reflection data profile obtained in S2; is the diffraction wave structure model vector to be solved, is the reflection wave structure model vector to be solved; represent In the model response relationship The forward response under represent In the model response relationship The forward response under .

9. A readable storage medium, characterized in that: Computer program instructions are stored thereon, and when the computer program instructions are executed by a processor, the fine imaging method of ground penetrating radar fusion inversion as claimed in any one of claims 1 to 8 is implemented.

10. An electronic device, characterized in that: include: At least one processor, at least one memory and computer program instructions stored in the memory, when the computer program instructions are executed by the processor, the fine imaging method of ground penetrating radar fusion inversion according to any one of claims 1 to 8.

Images