Orthogonal panel mesh based complex ground surface ground penetrating radar imaging method

Abstract

The application discloses a complex ground surface ground penetrating radar amplitude-preserving imaging method based on an orthogonal skin mesh, and comprises the following steps: 1) constructing a ground surface model according to ground fluctuation data; constructing an electromagnetic wave underground propagation speed model according to geological exploration data; 2) generating an orthogonal skin mesh by using a Ryskin-Leal method; 3) discretizing the orthogonal skin mesh by using a finite difference time domain method; 4) performing reverse extrapolation on the received ground penetrating radar signal, and reversely propagating the electromagnetic wave field to an initial time; in the reverse extrapolation process, the spatial position of a reflection interface is calculated through a cross-correlation imaging condition, and an underground imaging profile is obtained; 5) calculating a target function gradient, and solving an inversion problem by using a conjugate gradient method to iteratively update the electromagnetic wave underground propagation speed model; and 6) generating an imaging profile of an underground structure according to the iteratively converged electromagnetic wave underground propagation speed model. The application significantly improves the precision, stability and anti-noise performance of underground imaging.

Description

Technical Field

[0001] This invention relates to the field of radar imaging, specifically to a method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted grids. Background Technology

[0002] As geological exploration demands extend to complex environments, Ground Penetrating Radar (GPR) technology, with its high resolution and excellent penetration capabilities, is increasingly widely used in the field of underground structure detection, particularly demonstrating unique advantages in shallow geological surveys, archaeological excavations, and engineering inspections. To improve detection accuracy, numerical simulation methods are commonly used for forward modeling and imaging of the GPR wavefield. Among these methods, body-fitted meshes, as a mesh generation method that can effectively fit arbitrary undulating surface morphologies, significantly reduce terrain discretization errors and improve simulation accuracy under complex surface conditions compared to traditional stepped meshes.

[0003] Existing technologies suffer from at least the following drawbacks: First, traditional migration imaging methods are often based on the assumption of a horizontal surface, making them poorly adaptable to highly undulating terrain and unable to accurately reconstruct the spatial distribution of underground structures, leading to a decline in image quality. Second, in real, complex surface environments, radar wave propagation paths are significantly affected by irregular terrain, resulting in multiple refractions and scattering phenomena. This non-target energy mixes into the received signal, making it difficult to distinguish between effective and interference signals, thus reducing the efficiency of target reflection information extraction. Furthermore, due to the lack of accurate joint modeling methods for real surface terrain and underground medium parameters, existing methods struggle to achieve high-precision inversion of underground electromagnetic properties. In addressing these issues, researchers have attempted to improve image quality by introducing complex medium models, multi-scale algorithms, or pre-stack denoising techniques, but all suffer from problems such as high computational load, poor stability, or limited applicability. Especially when dealing with the combined effects of non-flat terrain and non-uniform media, traditional methods often struggle to balance imaging accuracy and computational efficiency. Summary of the Invention

[0004] The purpose of this invention is to provide a method for amplitude-preserving imaging of complex surfaces using ground-penetrating radar based on orthogonal body-fitted meshes, comprising the following steps:

[0005] 1) Construct a surface model based on ground undulation data; construct an electromagnetic wave underground propagation velocity model based on geological exploration data;

[0006] 2) The Ryskin–Leal method is used to generate orthogonal body meshes, thereby converting the coordinates of the physical domain to the coordinates of the computational domain, and mapping the velocity model and medium parameters of the physical domain to the computational domain;

[0007] 3) The orthogonal body mesh is discretized using the finite difference time-domain method to simulate the propagation of electromagnetic waves in the underground medium;

[0008] 4) In the computational domain, the received ground-penetrating radar signal is back extrapolated to propagate the electromagnetic wave field back to the initial time. During the back extrapolation process, the spatial position of the reflection interface is calculated based on the propagation time of the electromagnetic wave and the underground propagation velocity model, and the underground imaging profile is obtained by cross-correlation imaging conditions.

[0009] 5) Based on the orthogonal body mesh, the gradient of the objective function is calculated by iterative inversion method, and the conjugate gradient method is used to solve the optimization inversion problem. The electromagnetic wave underground propagation velocity model is iteratively updated until the imaging results are stable and the reflection event and interface location no longer show significant improvement.

[0010] 6) Based on the electromagnetic wave underground propagation velocity model after iterative convergence, extract the imaging profiles of all measuring points and superimpose them to generate the imaging profile of the underground structure.

[0011] Furthermore, the surface model is constructed based on ground undulation data;

[0012] The electromagnetic wave underground propagation velocity model was constructed based on geological exploration data.

[0013] Furthermore, in step 2), the steps of generating orthogonal body meshes using the Ryskin–Leal method include:

[0014] 2.1) Introduce a uniform equidistant grid into the computational domain and combine it with measured data to obtain the undulation height of the physical domain surface to reflect the actual surface conditions; the measured data are obtained by GPS / RTK surveying or topographic mapping.

[0015] 2.2) The nodes in the computational domain are mapped to the corresponding undulation heights in the physical domain (x,z) by piecewise linear interpolation, ensuring that the grid nodes are attached to the surface and medium interface;

[0016] 2.3) Iterative optimization of nodes within the mesh is performed using a distortion function based on the scale factor;

[0017] 2.4) Determine whether the node update magnitude meets the convergence criterion. If yes, output an orthogonal body-fitted mesh; otherwise, return to step 2.3. , For the node update magnitude; This is a preset threshold.

[0018] Furthermore, iterative optimization of the nodes inside the mesh refers to performing a weighted average of the nodes inside the mesh based on the distortion function, so that the boundary and interface nodes remain fixed.

[0019] Furthermore, the distortion function is as follows:

[0020] (1)

[0021] (2)

[0022] (3)

[0023] In the formula, , for and Orthogonal scale factors of direction; n and g represent respectively , The relative smoothness index. , It is the average scale.

[0024] Furthermore, in step 3), the step of discretizing the orthogonal body mesh using the finite difference time-domain method includes:

[0025] 3.1) Construct Maxwell's equations in a rectangular coordinate system, that is:

[0026] (4)

[0027] In the formula, H x and H z E represents the magnetic field components in the x and z directions. y Let ε be the electric field intensity in the y direction, σ be the permittivity, μ be the conductivity, and t be the time.

[0028] 3.2) Transform Maxwell's equations in rectangular coordinates to curvilinear coordinates. , The form below;

[0029] The curvilinear coordinate system is established by taking the origin at the starting point of the computational domain surface. The direction extends along the tangent to the surface undulation curve. The direction is taken as the normal direction perpendicular to the ground surface, thus constructing a curvilinear coordinate system that adapts to the actual terrain;

[0030] Curvilinear coordinate system , Maxwell's equations under ( ) are as follows:

[0031] (5)

[0032] In the formula, E z Let be the electric field intensity in the z-direction; J and det represent the Jacobian matrix and Jacobian determinant, respectively.

[0033] 3.3) Solve equation (5) using a fully staggered mesh to achieve discretization of the orthogonal body mesh; a fully staggered mesh means that different components of the electric field and magnetic field are defined on the same mesh, so that each variable is defined at a different position on the mesh.

[0034] Furthermore, the objective function As shown below:

[0035] (6)

[0036] in, Represents a simulated record; It represents the difference between the actual observation record and the background field record; M is the number of excitations by the transmitting antenna; N is the number of receivers; T is the recording duration; Let j be the spatial location of the receiving point. By introducing an orthogonal body-fitted mesh to calculate the objective function, the numerical calculation error caused by undulating terrain can be effectively reduced, thus improving the calculation accuracy.

[0037] Furthermore, the gradient g of the objective function is shown below:

[0038] (7)

[0039] Among them, H x and H z E represents the magnetic field components in the x and z directions. y Let ε be the electric field intensity in the y-direction, σ be the permittivity, μ be the conductivity, and t be the time. The orthogonal body-fitted mesh, through precise mesh generation, ensures high accuracy and reliability in gradient calculations. By performing detailed calculations on this mesh, the final imaging results can accurately reflect the electromagnetic properties of the underground structure.

[0040] Furthermore, in step 5), equation (6) is solved by the conjugate gradient method, with minimizing the objective function as the optimization objective, so as to finally obtain the perturbation value of the underground relative permittivity; this perturbation value is used to reflect the reflection characteristics of different materials in the underground structure, including the position, shape and electromagnetic properties of the object.

[0041] Furthermore, during the iterative update of the electromagnetic wave underground propagation velocity model, the influence of model parameter perturbation on the electromagnetic field response is characterized by the perturbation equation, and the relationship between the observation data residuals and model correction is established, providing a basis for the calculation of the objective function gradient and the iterative update of the model.

[0042] The perturbation equation is as follows:

[0043] (8)

[0044] In the formula, Ey Let be the electric field intensity in the y-direction. σ is the perturbation of the dielectric constant, μ is the conductivity, and t is the permeability.

[0045] The technical effects of this invention are undeniable. This invention combines body-fitting mesh technology with an amplitude-preserving imaging method to improve the accuracy and stability of ground-penetrating radar in underground detection and imaging processes.

[0046] By employing a body-fitted mesh, this invention can effectively solve the stepped discretization error caused by traditional regular meshes when processing undulating surfaces, ensuring that the mesh can accurately fit the interface between the surface and the subsurface medium, thereby avoiding interference from false scattering.

[0047] Meanwhile, this invention incorporates a amplitude-preserving imaging method, which, by inverting the collected radar data, can extract the true information of the underground medium and generate more accurate underground imaging results.

[0048] Specifically, this invention performs precise coordinate transformation between the computational and physical domains, ensuring that grid nodes can match surface undulations and medium interfaces, thereby improving the stability and computational accuracy of numerical simulations. The amplitude-preserving imaging method further enhances imaging quality through forward and backward propagation of wavefields and iterative optimization, making the inversion results closer to the actual underground structure.

[0049] The combined approach of this invention can generate high-resolution subsurface images in complex geological environments, especially in areas close to the surface, where the imaging effect is clearer and more consistent with the actual reflectivity model. Furthermore, the method of this invention has a strong noise suppression effect, and can still obtain relatively accurate imaging results even under low signal-to-noise ratio conditions.

[0050] In summary, this invention, by combining body-fitted meshes with a amplitude-preserving imaging method, significantly improves the accuracy, stability, and noise resistance of underground imaging. It is particularly suitable for underground exploration and imaging in complex geological environments and has broad application prospects. Attached Figure Description

[0051] Figure 1 This is a schematic diagram illustrating the conversion between the physical domain and the computational domain of this invention;

[0052] Figure 2 This is a flowchart of the imaging algorithm of the present invention;

[0053] Figure 3 This invention establishes a small rock block model and a perturbation model obtained from smoothing calculations;

[0054] Figure 4 This is a schematic diagram of orthogonal body mesh generation and simulation data obtained in this invention.

[0055] Figure 5 These are the imaging results corresponding to the small rock block model of this invention.

[0056] Figure 6 It is the perturbation model obtained by multi-layer geological model and smoothed model calculation in this invention.

[0057] Figure 7 These are radar data with different noise levels and the imaging results of the corresponding amplitude-preserving imaging method according to the present invention. Detailed Implementation

[0058] The present invention will be further described below with reference to embodiments, but it should not be construed that the scope of the present invention is limited to the following embodiments. Various substitutions and modifications made based on ordinary technical knowledge and common practices in the art without departing from the above-described technical concept of the present invention should be included within the scope of protection of the present invention.

[0059] Example 1:

[0060] See Figures 1 to 7 A method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted meshes includes the following steps:

[0061] 1) Construct a surface model based on ground undulation data; construct an electromagnetic wave underground propagation velocity model based on geological exploration data;

[0062] 2) The Ryskin–Leal method is used to generate orthogonal body meshes, thereby converting the coordinates of the physical domain to the coordinates of the computational domain, and mapping the velocity model and medium parameters of the physical domain to the computational domain;

[0063] 3) The orthogonal body mesh is discretized using the finite difference time-domain method to simulate the propagation of electromagnetic waves in the underground medium;

[0064] 4) In the computational domain, the received ground-penetrating radar signal is back extrapolated to propagate the electromagnetic wave field back to the initial time. During the back extrapolation process, the spatial position of the reflection interface is calculated based on the propagation time of the electromagnetic wave and the underground propagation velocity model, and the underground imaging profile is obtained by cross-correlation imaging conditions.

[0065] 5) Based on the orthogonal body mesh, the gradient of the objective function is calculated by iterative inversion method, and the conjugate gradient method is used to solve the optimization inversion problem. The electromagnetic wave underground propagation velocity model is iteratively updated until the imaging results are stable and the reflection event and interface location no longer show significant improvement.

[0066] 6) Based on the electromagnetic wave underground propagation velocity model after iterative convergence, extract the imaging profiles of all measuring points and superimpose them to generate the imaging profile of the underground structure.

[0067] Example 2:

[0068] The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body mesh is the same as that in Example 1. In addition, the surface model is constructed based on ground undulation data.

[0069] The electromagnetic wave underground propagation velocity model was constructed based on geological exploration data.

[0070] Example 3:

[0071] A method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body mesh, with technical content the same as any one of embodiments 1-2, further comprising the following steps in step 2), which include generating the orthogonal body mesh using the Ryskin–Leal method:

[0072] 2.1) Introduce a uniform equidistant grid into the computational domain and combine it with measured data to obtain the undulation height of the physical domain surface to reflect the actual surface conditions; the measured data are obtained by GPS / RTK surveying or topographic mapping.

[0073] 2.2) The nodes in the computational domain are mapped to the corresponding undulation heights in the physical domain (x,z) by piecewise linear interpolation, ensuring that the grid nodes are attached to the surface and medium interface;

[0074] 2.3) Iterative optimization of nodes within the mesh is performed using a distortion function based on the scale factor;

[0075] 2.4) Determine whether the node update magnitude meets the convergence criterion. If yes, output an orthogonal body mesh; otherwise, return to step 2.3. , For the node update magnitude; This is a preset threshold.

[0076] Example 4:

[0077] The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body mesh is the same as any one of Examples 1-3. Further, iterative optimization of the nodes inside the mesh refers to performing a weighted average of the nodes inside the mesh according to the distortion function, so that the boundary and interface nodes remain fixed.

[0078] Example 5:

[0079] A method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted meshes, with the same technical content as any one of Examples 1-4, further wherein the distortion function is as follows:

[0080] (1)

[0081] (2)

[0082] (3)

[0083] In the formula, , for and Orthogonal scale factors of direction; n and g represent respectively , The relative smoothness index. , It is the average scale.

[0084] Example 6:

[0085] A method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body mesh, with technical content identical to any one of embodiments 1-5, further comprising the following step 3), which involves discretizing the orthogonal body mesh using the finite difference time-domain method:

[0086] 3.1) Construct Maxwell's equations in a rectangular coordinate system, that is:

[0087] (4)

[0088] In the formula, H x and H z E represents the magnetic field components in the x and z directions. y Let ε be the electric field intensity in the y direction, σ be the permittivity, μ be the conductivity, and t be the time.

[0089] 3.2) Transform Maxwell's equations in rectangular coordinates to curvilinear coordinates. , The form below;

[0090] The curvilinear coordinate system is established by taking the origin at the starting point of the computational domain surface. The direction extends along the tangent to the surface undulation curve. The direction is taken as the normal direction perpendicular to the ground surface, thus constructing a curvilinear coordinate system that adapts to the actual terrain;

[0091] Curvilinear coordinate system , Maxwell's equations under ( ) are as follows:

[0092] (5)

[0093] In the formula, E z Let be the electric field intensity in the z-direction; J and det represent the Jacobian matrix and Jacobian determinant, respectively.

[0094] 3.3) Solve equation (5) using a fully staggered mesh to achieve discretization of the orthogonal body mesh; a fully staggered mesh means that different components of the electric field and magnetic field are defined on the same mesh, so that each variable is defined at a different position on the mesh.

[0095] Example 7:

[0096] A method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted meshes, with technical content identical to any one of Examples 1-6, further comprising the objective function... As shown below:

[0097] (6)

[0098] in, Represents a simulated record; It represents the difference between the actual observation record and the background field record; M is the number of excitations by the transmitting antenna; N is the number of receivers; T is the recording duration; It represents the spatial location of the j-th receiving point. By introducing an orthogonal body-fitted mesh to calculate the objective function, the numerical calculation error caused by undulating terrain can be effectively reduced, thus improving the calculation accuracy.

[0099] Example 8:

[0100] A method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted meshes is provided, with the technical content being the same as any one of Examples 1-7. Furthermore, the gradient g of the objective function is as follows:

[0101] (7)

[0102] Among them, H x and H z E represents the magnetic field components in the x and z directions. y Let ε be the electric field intensity in the y-direction, σ be the permittivity, μ be the conductivity, and t be the time. The orthogonal body-fitted mesh, through precise mesh generation, ensures high accuracy and reliability in gradient calculations. By performing detailed calculations on this mesh, the final imaging results can accurately reflect the electromagnetic properties of the underground structure.

[0103] Example 9:

[0104] The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body mesh is the same as any one of embodiments 1-8. Further, in step 5), equation (6) is calculated and solved by the conjugate gradient method, with minimizing the objective function as the optimization objective, so as to finally obtain the perturbation value of the underground relative permittivity; the perturbation value is used to reflect the reflection characteristics of different materials in the underground structure, including the position, shape and electromagnetic properties of the object.

[0105] Example 10:

[0106] A method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body mesh is provided, with the same technical content as any one of Examples 1-9. Furthermore, during the iterative update of the electromagnetic wave underground propagation velocity model, the influence of model parameter perturbation on electromagnetic field response is characterized by perturbation equations, and the relationship between observation data residuals and model correction is established, providing a basis for objective function gradient calculation and model iterative update.

[0107] The perturbation equation is as follows:

[0108] (8)

[0109] In the formula, E y Let be the electric field intensity in the y-direction. σ is the perturbation of the dielectric constant, μ is the conductivity, and t is the permeability.

[0110] Example 11:

[0111] A method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted meshes is described below:

[0112] Orthogonal body mesh generation algorithm

[0113] This invention provides a body-fitting mesh generation method. When undulating terrain exists, the discrete mesh used in numerical simulation needs to closely fit the undulating boundaries to avoid spurious scattering caused by stepped discretization. Such a mesh is a body-fitting mesh. The smoothness of the mesh interior and the orthogonality of the boundaries ensure the stability and accuracy of the numerical simulation. The generation of orthogonal body-fitting meshes is actually a mathematical transformation process, that is, transforming an arbitrary-shaped physical domain into a rectangular computational domain (such as...). Figure 1 (As shown). In this transformation, there is a one-to-one correspondence between points, which can be mathematically represented as:

[0114] (1)

[0115] To solve partial differential equations in the physical domain within the computational domain, we need to convert the partial derivatives in the physical domain to their partial derivatives in the computational domain. In a volume-fitted mesh, the partial derivatives in physical space can be converted to their partial derivatives in computational space using the inverse of the Jacobian matrix:

[0116] The mapped Jacobian matrix J and its determinant detJ can be used to convert physical domain partial derivatives into computational domain partial derivatives.

[0117] (2)

[0118] This paper employs a two-dimensional orthogonal body-fitting mesh generation strategy based on the Ryskin–Leal method. This method first introduces a uniform, equidistant mesh into the computational domain and calculates the undulation height of the physical domain's surface. Subsequently, through piecewise linear interpolation, nodes in the computational domain are mapped to the corresponding heights in the physical domain (x,z), ensuring that the mesh nodes are attached to the surface and medium interface.

[0119] Building upon this, to improve mesh smoothness and orthogonality, a distortion function based on the scale factor is introduced to iteratively optimize the internal nodes. Definition and The orthogonal scale factor for the direction is:

[0120] (3)

[0121] The average scales in the corresponding directions are denoted as follows: , Then, the relative smoothness index is defined as follows:

[0122] (4)

[0123] Then, the distortion function is constructed:

[0124] (5)

[0125] Using Gauss-Seidel iteration, the coordinates of internal nodes are updated by a weighted average based on a distortion function, keeping the boundary and interface nodes fixed. If the node update magnitude satisfies the convergence criterion... If the mesh iteration converges, then the mesh is considered to have converged. This method can effectively generate two-dimensional body-fitted meshes with good orthogonality and smoothness.

[0126] 2. Finite-Difference Time-Domain Method for Ground Penetrating Radar under Body-Fit Mesh

[0127] When solving electromagnetic field problems involving GPR, the constitutive relations of the medium are usually considered, and the calculation is simplified using a two-dimensional transverse magnetic (TM) model. Maxwell's equations in Cartesian coordinates are as follows:

[0128] (6)

[0129] Where Hx and Hz are the magnetic field components in the x and z directions, Ey is the electric field intensity in the y direction, ε is the permittivity, σ is the conductivity, μ is the permeability, and t is time.

[0130] In a body-fitted mesh, the above partial derivatives need to be converted to curvilinear coordinates. , The original equation can be written in the form of () using formula 2:

[0131] (7)

[0132] When solving the wave equation in a curvilinear coordinate system, the Standard Staggered Grid (SSG) is no longer applicable because each variable needs to have its spatial derivatives calculated in two directions at the same grid point. Therefore, this invention uses a Fully Staggered Grid (FSG) to solve equation (3). The core idea of ​​the FSG grid is that different components of the electric and magnetic fields are defined alternately on the same grid, so that each variable is defined at a different position on the grid, thereby avoiding interpolation operations and improving computational efficiency and accuracy.

[0133] 3. Ground Penetrating Radar Amplitude Preservation Imaging Algorithm with Body-Fit Mesh

[0134] Ground-penetrating radar's amplitude-preserving imaging algorithm combines backpropagation and the least squares method to achieve high-precision underground imaging in complex geological environments. By substituting the first two terms in equation (6), we can obtain the following formula:

[0135] (8)

[0136] When the medium is disturbed by ∆ε, the change in the wave field ∆Ey can be expressed as:

[0137] (9)

[0138] By subtracting equations (6) and (7) and ignoring higher-order minor quantities, we obtain:

[0139] (10)

[0140] Let x = ∆ε represent the change in dielectric constant, and we obtain the perturbation equation:

[0141] (9)

[0142] By iteratively updating x, its objective function is defined as:

[0143] (11)

[0144] Where ∆dcal represents the simulated record, which must satisfy equation (5); ∆dobs is the difference between the actual observation record and the background field record; M is the number of excitations of the transmitting antenna; N is the number of receivers; T is the recording duration; and rj is the spatial location of the j-th receiving point. Using the Lagrange multiplier method, the gradient formula of the objective function is:

[0145] (12)

[0146] Where Hx and Hz are the magnetic field components in the x and z directions, respectively, Ey is the electric field intensity in the y direction, ε is the permittivity, σ is the conductivity, μ is the permeability, and t is time. The orthogonal body-fitted mesh, through precise mesh generation, ensures high accuracy and reliability in gradient calculations. By performing detailed calculations on this mesh, the final imaging results can accurately reflect the electromagnetic properties of the underground structure.

[0147] Finally, the perturbation value of the underground relative permittivity is obtained by calculating using the conjugate gradient method, with minimizing the objective function as the optimization objective. This perturbation value can effectively reflect the reflection characteristics of different materials in the underground structure, including the location, shape, and electromagnetic properties of the object.

[0148] Figure 2 A flowchart of the imaging algorithm of the present invention is shown, which includes the following steps:

[0149] Step S1, Construct a body-fitted mesh model: Based on ground undulation data, construct an accurate surface model; based on geological exploration data, construct a subsurface velocity model as the basis for wave field propagation;

[0150] Step S2, conversion between computational and physical domains: The body-fitted mesh generated by the Ryskin–Leal method is used to convert the coordinates of the physical domain to the coordinates of the computational domain; the velocity model, medium parameters, etc. in the physical domain are mapped to the computational domain to adapt to the numerical calculation of the finite difference method.

[0151] Step S3, Forward Wave Field Calculation: The computational domain is efficiently discretized, and the propagation of electromagnetic waves in the underground medium is simulated using the finite difference time-domain method.

[0152] Step S4, Reverse Wavefield Calculation: The collected ground-penetrating radar signals are back extrapolated until the electromagnetic wavefield returns to its initial state. During this process, the location of the reflecting interface is calculated based on the propagation time and velocity model of electromagnetic waves to obtain the underground imaging profile. Finally, the imaging profiles of all measuring points are superimposed and integrated to obtain the imaging result.

[0153] Step S5, Iterative Update: Based on the calculated gradient, update the velocity model by solving the inversion problem using the conjugate gradient method.

[0154] Step S6, Imaging result output: Generate an imaging profile of the underground structure based on the final updated velocity model.

[0155] Example 12:

[0156] The verification of the amplitude-preserving imaging method for complex surface ground-penetrating radar based on orthogonal body mesh is as follows:

[0157] To verify the advantages of body-fitted meshes over regular meshes in data interpretation, a multi-block model was established as follows, with the actual model shown in Figure [image missing]. Figure 3 (a) shows a model with dimensions of 1.0m × 2.0m, a time interval dt = 0.016ns, a time step of 900, and a surface with a certain degree of undulation. A smoothed real model is used as the initial model for the migration process. To obtain the initial model, a Gaussian filter with a template size of 40 × 40 and a standard deviation of 4 is used for smoothing; these two parameters determine the degree of smoothness relative to the real model. The perturbation model calculated from the real model and the smoothed model is shown below. Figure 3 As shown in b.

[0158] The schematic diagram of orthogonal body mesh generation and the simulation data of this invention are as follows: Figure 4 As shown, the radar acquisition center frequency is 1.5 GHz, and 100 channels of ground-penetrating radar data in single-transmit, single-receive mode are uniformly acquired at the ground surface. This study uses an amplitude-preserving imaging algorithm as the data interpretation algorithm. A smoothing model is used as the initial velocity model, and calculations are performed using a conventional reverse time migration (RTM) imaging algorithm, an amplitude-preserving imaging algorithm under a regular grid, and an amplitude-preserving imaging algorithm under a body-fitted grid.

[0159] In this experiment, the maximum number of iterations in the calculation was 50. Figure 5 'a' represents the RTM result; Figure 5 b and Figure 5 c represents the amplitude-preserving imaging results in the computational and physical domains obtained using the body-fitted mesh. Figure 5 d represents the amplitude-preserving imaging result of a traditional regular grid.

[0160] The results show that RTM cannot capture the dielectric constant variations of small rock fragments, and the location and shape of the fragments are not clearly displayed. In contrast, amplitude-preserving imaging can capture the dielectric constant variations, and the imaging results are consistent with the dielectric constant variations of the near-real model. Compared with the results from rectangular meshes, the near-surface portion of the amplitude-preserving imaging results from body-fitted meshes is clearer, and the morphology and properties of subsurface anomalies are closer to the real reflectivity model.

[0161] To further verify the noise resistance of the amplitude-preserving imaging method under body-fitted mesh, this paper constructs a multi-layer geological model with undulating surface, the real model of which is as follows: Figure 6 (a) shows a model with dimensions of 1.0m × 3.0m, a time interval dt = 0.01ns, a time step of 2000, and a surface with a certain degree of undulation. A Gaussian filter of the template size is used for smoothing; the perturbation model calculated from the real model and the smoothed model is as follows. Figure 6 As shown in b.

[0162] The radar's center frequency was 1.5 GHz, and 100 channels of ground-penetrating radar data were uniformly collected at the ground surface in single-transmit, single-receive mode. In the numerical simulation tests of this study, Gaussian white noise with different noise levels relative to the effective signal was added to the collected data before imaging calculations were performed. Figure 7 (a) and Figure 7 (c) Received data after adding different levels of noise. The effective signal is relatively clear when the signal-to-noise ratio is -5dB. Figure 7 (a)); The signal-to-noise ratio is −10dB, and the effective signal is barely visible. Figure 7 (c)). In this experiment, a smoothing model is used as the initial velocity model, and amplitude-preserving imaging under body-fitted mesh is used for calculation.

[0163] Figure 7 (b) and Figure 7 (d) Results of amplitude-preserving imaging at the corresponding noise level are given. At a signal-to-noise ratio of -5 dB, noise has virtually no impact on the imaging results. Figure 7 (b) At a signal-to-noise ratio of −10 dB, noise has a significant impact on the imaging results, with many cluttered artifacts appearing in the background, but anomalous objects can still be identified relatively well. Figure 7 (d)). The numerical test results show that the algorithm has a good noise suppression effect.

Claims

1. A method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted meshes, characterized in that, Includes the following steps: S1) Construct a surface model based on ground undulation data; Based on geological exploration data, a model of the underground propagation velocity of electromagnetic waves was constructed. S2) The Ryskin–Leal method is used to generate an orthogonal body mesh, thereby converting the coordinates of the physical domain to the coordinates of the computational domain, and mapping the velocity model and medium parameters of the physical domain to the computational domain; S3) The orthogonal body mesh is discretized using the finite difference time-domain method to simulate the propagation of electromagnetic waves in the underground medium; S4) In the computational domain, the received ground-penetrating radar signal is extrapolated backward to propagate the electromagnetic wave field back to the initial time. In the reverse extrapolation process, based on the propagation time of electromagnetic waves and the underground propagation velocity model, the spatial location of the reflection interface is calculated through cross-correlation imaging conditions to obtain the underground imaging profile. S5) Based on the orthogonal body mesh, the gradient of the objective function is calculated by the iterative inversion method, and the conjugate gradient method is used to solve the optimization inversion problem. The electromagnetic wave underground propagation velocity model is iteratively updated until the imaging results are stable. S6) Based on the electromagnetic wave underground propagation velocity model after iterative convergence, extract the imaging profiles of all measuring points and superimpose them to generate the imaging profile of the underground structure.

2. The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted mesh according to claim 1, characterized in that, The surface model is constructed based on ground undulation data; The electromagnetic wave underground propagation velocity model was constructed based on geological exploration data.

3. The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted mesh according to claim 1, characterized in that, Step S2), the steps of generating an orthogonal body mesh using the Ryskin–Leal method include: S2.1) Introduce a uniform equidistant grid into the computational domain and combine it with measured data to obtain the undulation height of the physical domain surface to reflect the actual surface conditions; the measured data is obtained by GPS / RTK surveying or topographic mapping. S2.2) The nodes in the computational domain are mapped to the corresponding undulation heights in the physical domain (x,z) by piecewise linear interpolation, ensuring that the grid nodes are attached to the surface and medium interface; S2.3) Iterative optimization of nodes within the mesh is performed using a distortion function based on the scale factor; S2.4) Determine whether the node update magnitude satisfies the convergence criterion. If yes, output an orthogonal body-fitted mesh; otherwise, return to step S2.

3. , For the node update magnitude; This is a preset threshold.

4. The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted mesh according to claim 3, characterized in that, Iterative optimization of nodes within a mesh refers to applying a weighted average to the nodes within the mesh based on a distortion function, thereby keeping the boundary and interface nodes fixed.

5. The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted mesh according to claim 3, characterized in that, The distortion function is shown below: (1) (2) (3) In the formula, , for and Orthogonal scale factors of direction; n and g represent respectively , The relative smoothness index; f is the distortion function; , It is the average scale.

6. The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body mesh according to claim 1, characterized in that, Step S3), the discretization of the orthogonal body mesh using the finite difference time-domain method includes: S3.1) Construct Maxwell's equations in a rectangular coordinate system, that is: (4) In the formula, H x and H z E represents the magnetic field components in the x and z directions. y Let ε be the electric field intensity in the y direction, σ be the permittivity, μ be the conductivity, and t be the time. S3.2) Transform Maxwell's equations in rectangular coordinates to curvilinear coordinates. , The form below; The curvilinear coordinate system is established by taking the origin at the starting point of the computational domain surface. The direction extends along the tangent to the surface undulation curve. The direction is taken as the normal direction perpendicular to the ground surface, thus constructing a curvilinear coordinate system that adapts to the actual terrain; Curvilinear coordinate system , Maxwell's equations under ( ) are as follows: (5) In the formula, E z Let be the electric field intensity in the z-direction; J and det represent the Jacobian matrix and Jacobian determinant, respectively. S3.3) The equation (5) is solved by using a fully staggered mesh to achieve the discretization of the orthogonal body mesh; the fully staggered mesh means that the different components of the electric field and the magnetic field are defined on the same mesh, so that each variable is defined at a different position on the mesh.

7. The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body-fitted mesh according to claim 1, characterized in that, objective function As shown below: (6) in, Represents a simulated record; It represents the difference between the actual observation record and the background field record; M is the number of excitations by the transmitting antenna; N is the number of receivers; T is the recording duration; It is the spatial location of the j-th receiving point.

8. The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body mesh according to claim 1, characterized in that, The gradient g of the objective function is shown below: (7) Among them, H x and H z E represents the magnetic field components in the x and z directions. y Let t be the electric field intensity in the y direction and t be time.

9. The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body mesh as described in claim 1, characterized in that, In step S5), equation (6) is solved by the conjugate gradient method, with minimizing the objective function as the optimization objective, so as to finally obtain the perturbation value of the underground relative permittivity; this perturbation value is used to reflect the reflection characteristics of different materials in the underground structure.

10. The method for amplitude-preserving imaging of complex surface ground-penetrating radar based on orthogonal body mesh according to claim 1, characterized in that, During the iterative update of the electromagnetic wave underground propagation velocity model, the influence of model parameter perturbation on electromagnetic field response is characterized by perturbation equations, and the relationship between observation data residuals and model correction is established. The perturbation equation is as follows: (8) In the formula, E y Let be the electric field intensity in the y-direction. σ is the perturbation of the dielectric constant, μ is the conductivity, and t is the permeability.

Images