Ground penetrating radar forward method based on finite field

Abstract

The application discloses a ground penetrating radar forward method based on a finite field, which comprises the following steps: constructing a ground penetrating radar finite field forward model; setting the positions of transmitting antennas and receiving antennas, the center frequency of the transmitting antennas, and the attenuation proportion of the finite field boundary; calculating the finite field range corresponding to the transmitting antennas; calculating the ground penetrating radar finite field forward parameters by adding a perfect matched layer; loading a pulse source for the set transmitting antenna position, updating the simulated area electromagnetic field value at the same time; updating the electromagnetic field auxiliary field; increasing the time step, repeating the above steps until the numerical simulation of the whole time is completed; obtaining the wave field data of the receiving antennas by using the set receiving antenna positions; repeating the above steps until all the transmitting antennas complete the excitation processing, thereby completing the ground penetrating radar finite field forward calculation and obtaining the radar profile; and the calculation cost of the method is reduced and the calculation efficiency is improved.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration technology, specifically relating to a ground-penetrating radar forward modeling method based on a finite domain. Background Technology

[0002] Ground-penetrating radar (GPR) is characterized by its non-destructive, efficient, and high-precision nature, making it an important shallow geophysical exploration method. It is widely used in engineering exploration, hydrological environment, disaster relief, pipeline detection, road inspection, and other fields. The transmitting and receiving antennas of GPR are packaged together and move along the ground surface, enabling large-scale data acquisition in a short time. The high efficiency and high sampling rate of the acquisition method result in a huge dataset, requiring significant manpower and computing costs for interpretation. Therefore, higher demands are placed on the efficiency of data processing.

[0003] GPR forward modeling, as a tool, simulates the radar wave propagation process and obtains theoretical observation data through numerical methods, which helps to understand the response characteristics of radar to different media. In addition, as the basis of reverse time migration imaging and full waveform inversion, forward modeling directly affects the accuracy and speed of the interpretation of measured data. Improving the speed of forward modeling is of great significance for promoting the timeliness of data interpretation to meet the engineering requirements.

[0004] Currently, commonly used forward modeling methods for GPR mainly include the Finite Difference Method (FDTD), the Finite Element Method (FETD), and the Finite Volume Method (FVTD). Among them, FETD and FVTD can use unstructured meshes, which makes spatial discretization more flexible and suitable for solving irregular models. However, FETD requires solving the inverse of a large stiffness matrix at each time step, resulting in high computational complexity and making parallel computing difficult. FVTD requires a large number of meshes to ensure numerical accuracy, which leads to a large amount of computation. FDTD transforms the Maxwell equations in differential form into a set of difference equations, which has the advantages of easy implementation, fast operation speed, and high-efficiency parallel computing. It is the most widely used method in GPR forward modeling. GPR forward modeling is a typical multi-source problem. After each point source excitation, the wave field of the simulation region needs to be calculated. When the forward model is large or the mesh is dense, the computational memory and computation time increase sharply, making it difficult to perform forward modeling on microcomputers with limited running memory.

[0005] Considering that GPR has the characteristics of high antenna frequency, the detection objects are mostly lossy media, the radar wave attenuation is fast, the detection depth and the range of influence are very limited, and most of the effective signal is concentrated near the antenna. When performing forward modeling for areas far from the antenna, the wave field intensity can be ignored and has little significance for the detection results. Instead, it causes a serious waste of computing resources. When performing forward modeling for larger models, the computing memory is too large and it is difficult to implement on a microcomputer. At the same time, the problem of long forward modeling time also makes it difficult for data processing steps such as migration and inversion to meet the timeliness required by engineering.

[0006] In summary, current GPR forward modeling methods often suffer from high computational costs and relatively low computational efficiency when dealing with large-scale models. Summary of the Invention

[0007] The purpose of this invention is to provide a finite-domain-based ground-penetrating radar forward modeling method that reduces computational costs and improves computational efficiency while ensuring accuracy.

[0008] The ground-penetrating radar forward modeling method based on a finite domain provided by this invention includes the following steps:

[0009] S1. Construct a finite-domain forward model for ground-penetrating radar;

[0010] S2. Set the positions of the ground-penetrating radar transmitting and receiving antennas, the center frequency of the transmitting antenna, and the finite domain boundary attenuation ratio;

[0011] S3. Using the antenna position and center frequency set in step S2, and the physical property parameters of the forward model constructed in step S1, calculate the finite domain range corresponding to the transmitting antenna;

[0012] S4. Using the finite domain range obtained in step S3, calculate the finite domain forward modeling parameters of ground penetrating radar by adding a perfect matching layer;

[0013] S5. Apply a loading pulse source to the transmitting antenna position set in step S2, and update the electromagnetic field value of the simulated area at the same time;

[0014] S6. Update the electromagnetic auxiliary field in the perfectly matched layer region;

[0015] S7. Add time steps and repeat steps S5-S6 above until the numerical simulation of the entire time is completed;

[0016] S8. Using the position of the receiving antenna set in step S2, obtain the wave field data of the receiving antenna;

[0017] S9. Repeat steps S3-S8 above until all transmitting antennas have completed the excitation process, and then complete the finite-domain forward modeling calculation of the ground penetrating radar to obtain the radar profile.

[0018] Step S1, which involves constructing a finite-domain forward model for ground-penetrating radar, specifically includes:

[0019] Establish a finite domain (FD) forward model for ground penetrating radar (GPR) based on application requirements;

[0020] The physical properties of the forward model include relative permittivity and conductivity;

[0021] Step S2, which involves setting the positions of the ground-penetrating radar transmitting and receiving antennas, the center frequency of the transmitting antenna, and setting the finite-domain boundary attenuation ratio, specifically includes:

[0022] Setting the location of the GPR transmitting antenna specifically includes:

[0023] Set the number, starting position, spacing, and ending position of the transmitting antennas;

[0024] Setting the GPR receiving antenna position specifically includes:

[0025] Set the number, start position, spacing, and end position of the receiving antennas, and the transmit / receive distance between the transmitting and receiving antennas;

[0026] Setting the center frequency of the transmitting antenna, specifically including:

[0027] The higher the center frequency of the transmitting antenna, the faster the electromagnetic waves attenuate and the shallower the detection depth.

[0028] Different transmit antenna frequencies need to be selected based on the model depth. The frequency is determined using the following formula:

[0029]

[0030] Where f represents the center frequency of the radar transmitting antenna; x represents the spatial resolution; ε r Represents the relative permittivity of the surrounding rock;

[0031] The finite field boundary decay ratio R0 affects the acceleration effect. The larger R0 is, the better the acceleration effect, but the error also increases accordingly.

[0032] The optimal range of R0 values ​​for different antenna frequencies was determined experimentally, specifically including:

[0033] When the antenna frequency is 100MHz, the optimal value range for R0 is 0.01 to 0.08.

[0034] When the antenna frequency is 200MHz, the optimal value range for R0 is 0.025 to 0.1.

[0035] When the antenna frequency is 400 MHz, the optimal value range of R0 is 0.06 - 0.15;

[0036] When the antenna frequency is 900 MHz, the optimal value range of R0 is 0.09 - 0.19;

[0037] Using the antenna position and center frequency set in step S2 and the physical property parameters of the forward model constructed in step S1, calculate the finite domain range corresponding to the transmitting antenna in step S3, specifically including:

[0038] (3 - 1) Calculation of the transmission coefficient:

[0039] Use the following formula to calculate the transmission coefficient T1 when the radar wave vertically enters below the transmitting antenna:

[0040]

[0041] where ε1 represents the relative permittivity of air; ε2 represents the relative permittivity of the formation below the transmitting antenna;

[0042] (3 - 2) Calculation of the medium absorption attenuation ratio:

[0043] Starting from the transmitting antenna, use the following formula to calculate the attenuation ratio A caused by medium absorption when propagating different distances in the set direction N :

[0044]

[0045] where A N represents the attenuation ratio caused by medium absorption when the radar wave propagates different distances in the set direction; Δx represents the spatial step of the grid; if a certain point on the ground is N grids away from the transmitting antenna, the propagation distance of the electromagnetic wave is equal to NΔx; α i represents the minimum attenuation coefficient in the vertical direction at the position of the i-th grid from the transmitting antenna in the model, and is calculated using the following formula:

[0046]

[0047] where ω represents the angular frequency, and the calculation formula is as follows:

[0048] ω = 2πf

[0049] where f represents the center frequency of the transmitting antenna; ε represents the permittivity of the medium, with the unit of F / m; μ represents the magnetic permeability of the medium, with the unit of H / m; σ represents the conductivity of the medium, with the unit of S / m;

[0050] (3 - 3) Calculation of spherical spreading attenuation:

[0051] Calculate the attenuation B caused by spherical diffusion propagating different distances in the set direction N :

[0052] (1) Two-dimensional case:

[0053]

[0054] (2) Three-dimensional case:

[0055]

[0056] where d tr represents the distance between the receiving antenna and the transmitting antenna; d N represents the distance that the electromagnetic wave propagates N grid cells from the transmitting antenna, and the calculation formula is as follows:

[0057] d N = NΔx

[0058] (3-4) Calculation of the total attenuation ratio:

[0059] Use the following formula to calculate the total attenuation ratio R at different distances propagated from the transmitting antenna N :

[0060] R N = T1·A N ·B N

[0061] (3-5) Attenuation ratio screening condition:

[0062] Use the following formula to represent the selection condition:

[0063] R N ≤ R0

[0064] where R0 represents the set attenuation ratio of the finite domain boundary;

[0065] (3-6) Calculation of the finite domain boundary:

[0066] Select the points that meet the attenuation ratio screening condition, and select the point closest to the transmitting antenna from the selected points as the finite domain boundary in the set direction;

[0067] (3-7) Calculation of the maximum finite domain boundary:

[0068] Use the following formula to calculate the maximum finite domain boundary corresponding to the recording time length constraint condition in the set direction:

[0069]

[0070] where t NIt represents the time for electromagnetic waves to propagate through N grids; i represents the cumulative number of times in the calculation process, starting from 1 and gradually increasing to N; v i It represents the electromagnetic wave propagation speed in the vertical direction at a point where the model is i grids away from the transmitting antenna; Δx represents the spatial step size;

[0071] The selection condition is expressed by the following formula:

[0072] t N ≥T / 2

[0073] Select the points that meet the above conditions, and select the point closest to the transmitting antenna from the selected points as the maximum finite domain boundary;

[0074] If the finite domain boundary is greater than the maximum finite domain boundary, change the corresponding finite domain boundary to the position of the maximum finite domain boundary;

[0075] (3-8) Determination of the final finite domain range:

[0076] Repeat the above steps (3-2)-(3-7) until the finite domain boundaries in each direction are calculated, and then determine the final finite domain range;

[0077] For the finite domain range calculated by step S3 in step S4, by adding a perfectly matched layer, calculate the forward modeling parameters of the ground penetrating radar finite domain, specifically including:

[0078] Before exciting the transmitting antenna and calculating the field strength components, according to the known data, calculate the parameters required for forward modeling of the finite domain obtained in step S3, including the convolutional perfectly matched layer (CPML) absorbing boundary and the finite difference time domain (FDTD) in the forward modeling method;

[0079] (4-1) Calculation of the auxiliary absorption coefficient in the CPML region:

[0080] Before performing forward modeling, it is necessary to set CPML outside the simulation region and set the parameters inside CPML, specifically including:

[0081] Calculate the coordinate stretching factor s using the following formula i :

[0082]

[0083] where κ i represents the auxiliary absorption coefficient, which is used to improve the absorption effect of PML on surface waves; σ i represents the conductivity parameter in the i direction inside the CPML layer; αi It represents the auxiliary absorption coefficient, which is used to improve the absorption effect of PML on low-frequency components; j represents the imaginary unit; ω represents the angular frequency; ε0 represents the permittivity of vacuum;

[0084] The thickness of CPML is finite, and σ, κ, and α can vary monotonically in the layer;

[0085] The following formula is used to describe the setting of the three parameters in the z direction:

[0086]

[0087] where z0 represents the interface between CPML and the simulation region; m represents a parameter; d represents the thickness of CPML; σ max represents the maximum conductivity within CPML, and the calculation formula is as follows:

[0088]

[0089] where ε r represents the relative permittivity of the CPML layer; δ represents the cell size;

[0090]

[0091] where κ max represents the maximum value of κ within CPML;

[0092]

[0093] where α max represents the maximum value of α within CPML;

[0094] (4-2) Calculation of forward parameters of the field strength component update equation:

[0095] According to the two-dimensional Maxwell equation and the constitutive equation, the decoupling equation of the two-dimensional time-domain TM mode is described by the following formula:

[0096]

[0097] where x represents the x direction in the rectangular coordinate system; y represents the y direction in the rectangular coordinate system; z represents the z direction in the rectangular coordinate system; H x represents the component of the magnetic field strength in the x direction; H y represents the component of the magnetic field strength in the y direction; E z represents the component of the electric field strength in the z direction; ε represents the permittivity of the medium; t represents time; σ represents the conductivity of the medium; J z represents the applied current source; μ represents the permeability of the medium;

[0098] Expanding the above formula in the Yee grid gives the update equations for the field strength components as follows:

[0099]

[0100] where the superscript defining the field strength represents time, the subscript represents the component direction, and the content inside the parentheses represents the position in the grid; represents the x - direction magnetic field component at the position (i,j) at the time n + 1 / 2; represents the x - direction magnetic field component at the position (i,j) at the time n - 1 / 2; represents the z - direction electric field component at the position (i,j) at the time n; represents the z - direction electric field component at the position (i,j + 1) at the time n; Δy represents the spatial interval of the grid in the y - direction; CQ(m) represents the medium parameter at the spatial position where the magnetic field or electric field on the left - hand side of the above formula is located, and the calculation formula is as follows:

[0101]

[0102] where, Δt represents the time interval; μ(m) represents the magnetic permeability of the medium at the position m;

[0103]

[0104] where, represents the y - direction magnetic field component at the position (i,j) at the time n + 1 / 2; represents the y - direction magnetic field component at the position (i,j) at the time n - 1 / 2; represents the z - direction electric field component at the position (i + 1,j) at the time n; Δx represents the spatial interval of the grid in the x - direction;

[0105]

[0106] where, represents the z - direction electric field component at the position (i,j) at the time n + 1; represents the y - direction magnetic field component at the position (i + 1,j) at the time n + 1 / 2; represents the y - direction magnetic field component at the position (i,j) at the time n + 1 / 2; represents the x - direction magnetic field component at the position (i,j + 1) at the time n + 1 / 2; represents the x - direction magnetic field component at the position (i,j) at the time n + 1 / 2; represents the applied current density in the z - direction at the position (i,j) at the time n + 1 / 2; CA(m) represents the update equation coefficient, and the calculation formula is as follows:

[0107]

[0108] Among them, σ(m) represents the dielectric conductivity at position m; ε(m) represents the dielectric permittivity at position m;

[0109] CB(m) represents the coefficient of the update equation, and the calculation formula is as follows:

[0110]

[0111] Perform the processing of loading the pulse source for the position of the transmitting antenna set in step S2 in step S5, and at the same time update the electromagnetic field values in the simulation area, which specifically includes:

[0112] For the position of the transmitting antenna set in step S2, the loading of the pulse source is represented by the following formula:

[0113]

[0114] Among them, i s , j s represent the position of the transmitting antenna in the grid; represents the z-direction electric field intensity component at the position of the transmitting antenna at the (n + 1)th moment; represents the z-direction current intensity of the current source at the (n + 1 / 2)th moment;

[0115] Update the magnetic field strength component by using the following formula;

[0116]

[0117]

[0118] Update the electric field component by using the following formula:

[0119]

[0120] Update the electromagnetic field auxiliary field in the perfectly matched layer region in step S6, which specifically includes: perform the update processing of the electromagnetic field auxiliary field for CPML, which specifically includes:

[0121] (6 - 1) Calculate the auxiliary expression of the magnetic field update equation in CPML:

[0122] Calculate the auxiliary expression of the magnetic field update equation in CPML by using the following formula:

[0123]

[0124]

[0125] Among them, b i and a i The calculation formulas of are as follows:

[0126]

[0127]

[0128] (6 - 2) Calculate the magnetic field in the CPML region:

[0129] Use the following formula to calculate the magnetic field in the CPML region:

[0130]

[0131]

[0132] (6 - 3) Calculate the auxiliary expression of the electric field update equation in the CPML:

[0133] Use the following formula to calculate the auxiliary expression of the electric field update equation in the CPML:

[0134]

[0135]

[0136] (6 - 4) Calculate the electric field in the CPML region:

[0137] Use the following formula to calculate the electric field in the CPML region:

[0138]

[0139] The step of obtaining the wave field data of the receiving antenna by using the position of the receiving antenna set in step S2 described in step S8 specifically includes:

[0140] Output the E data at the position of the receiving antenna from the full wave field information according to the position of the receiving antenna, z as the radar record received by the receiving antenna after the current excitation of the transmitting antenna;

[0141] The forward - sounding method of ground - penetrating radar based on a finite domain provided by the present invention combines the attenuation characteristics of GPR high - frequency electromagnetic waves, limits the GPR forward - sounding region near the radar antenna, and reduces the calculation scale. In addition, the method of the present invention proposes a constraint condition according to the receiving time length to prevent the situation of an overly large finite domain, further ensuring the speed of finite - domain forward - sounding. Compared with the traditional GPR forward - sounding, the method of the present invention can effectively reduce the calculation scale, reduce the operation memory, allow more channels to be parallel, and significantly improve the forward - sounding efficiency on the premise of ensuring the accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0142] Figure 1 It is a schematic diagram of the method flow of the method of the present invention.

[0143] Figure 2 Schematic diagram for comparison between the simulation area of the method of the present invention and that of the conventional method

[0144] Figure 3 Schematic diagram for calculating the finite domain boundary of the method of the present invention

[0145] Figure 4 Schematic diagram of the relative permittivity and conductivity of Model 1 in the method of the present invention

[0146] Figure 5 Schematic diagram of the forward modeling time-consuming and relative error curves of Model 1 in the method of the present invention under different attenuation ratios (R0)

[0147] Figure 6 Schematic diagram of the relative permittivity and conductivity of Model 2 in the method of the present invention: Figure 6 (a) Schematic diagram of the relative permittivity distribution of the model Figure 6 (b) Schematic diagram of the conductivity distribution of the model

[0148] Figure 7 Schematic diagram of the forward modeling result and the finite domain forward modeling error of Model 2 in the method of the present invention: Figure 7 (a) Schematic diagram of the forward modeling result Figure 7 (b) Schematic diagram of the error between the forward modeling result and the conventional forward modeling result

[0149] Figure 8 Schematic diagram for comparison between the 94th trace conventional forward modeling and the finite domain forward modeling wave fields of Model 2 in the method of the present invention

[0150] Figure 9 Schematic diagram of Model 3 in the method of the present invention: Figure 9 (a) Schematic diagram of the relative permittivity distribution of the model Figure 9 (b) Schematic diagram of the conductivity distribution of the model

[0151] Figure 10 Schematic diagram of the relative permittivity distribution profiles of three survey lines of Model 3 in the method of the present invention and the corresponding finite domain forward modeling results: Figure 10 (a) Schematic diagram of the relative permittivity profile of the survey line at x = 4m Figure 10 (b) Schematic diagram of the relative permittivity profile of the survey line at x = 4m Figure 10 (c) Schematic diagram of the relative permittivity profile of the survey line at x = 8m Figure 10 (d) Schematic diagram of the finite domain forward modeling result of the survey line at x = 8m Figure 10 (e) Schematic diagram of the finite domain forward modeling result of the survey line at x = 12m Figure 10 (f) Schematic diagram of the finite domain forward modeling result of the survey line at x = 12m Detailed implementation mode

[0152] As Figure 1 shown in the following is the schematic diagram of the method flow of the method of the present invention: The forward ground penetrating radar method based on finite field provided by the present invention includes the following steps:

[0153] S1. Construct a finite domain forward model of ground penetrating radar; specifically including:

[0154] Establish a finite domain (Finite Domain, FD) forward model of ground penetrating radar (GPR) according to application requirements;

[0155] The physical property parameters of the forward model include relative permittivity and conductivity;

[0156] S2. Set the positions of the transmitting antenna and the receiving antenna of the ground penetrating radar, and the center frequency of the transmitting antenna, and at the same time set the finite domain boundary attenuation ratio; specifically including:

[0157] Set the position of the GPR transmitting antenna, specifically including:

[0158] Set the number, starting position, spacing, and ending position of the transmitting antenna;

[0159] Set the position of the GPR receiving antenna, specifically including:

[0160] Set the number, starting position, spacing, and ending position of the receiving antenna, and the transceiver distance between the transmitting antenna and the receiving antenna;

[0161] Set the center frequency of the transmitting antenna, specifically including:

[0162] The higher the center frequency of the transmitting antenna, the faster the electromagnetic wave attenuates and the shallower the detection depth;

[0163] It is necessary to select different transmitting antenna frequencies according to the model depth, and use the following formula to determine the frequency:

[0164]

[0165] where, f represents the center frequency of the radar transmitting antenna; x represents the spatial resolution; ε r represents the relative permittivity of the surrounding rock;

[0166] The finite domain boundary attenuation ratio R0 affects the acceleration effect. The larger R0 is, the better the corresponding acceleration effect is, and the misunderstanding error also increases accordingly;

[0167] Determine the optimal value range of R0 at different antenna frequencies through experiments, specifically including:

[0168] When the antenna frequency is 100 MHz, the optimal value range of R0 is 0.01 to 0.08;

[0169] When the antenna frequency is 200 MHz, the optimal value range of R0 is 0.025 to 0.1;

[0170] When the antenna frequency is 400 MHz, the optimal value range of R0 is 0.06 to 0.15;

[0171] When the antenna frequency is 900 MHz, the optimal value range of R0 is 0.09 to 0.19;

[0172] S3. Calculate the finite domain range corresponding to the transmitting antenna by using the antenna position and center frequency set in step S2, and the physical property parameters of the forward model constructed in step S1; specifically including:

[0173] (3-1) Calculation of the transmission coefficient:

[0174] Use the following formula to calculate the transmission coefficient T1 when the radar wave vertically enters below the transmitting antenna:

[0175]

[0176] where, ε1 represents the relative dielectric constant of air; ε2 represents the relative dielectric constant of the formation below the transmitting antenna;

[0177] (3-2) Calculation of the attenuation ratio:

[0178] As Figure 2 shown is the comparison schematic diagram of the simulation area of the method of the present invention and the simulation area of the conventional method: Starting from the transmitting antenna, use the following formula to calculate the attenuation ratio A of the propagation medium absorption in the set direction N :

[0179]

[0180] where, A N represents the attenuation ratio caused by the absorption of the radar wave in the propagation medium at different distances in the set direction; Δx represents the spatial step of the grid; if a certain point on the ground surface is N grids away from the transmitting antenna, the propagation distance of the electromagnetic wave is equal to NΔx; α i represents the minimum attenuation coefficient in the vertical direction at the position of the i-th grid of the model away from the transmitting antenna, and is calculated by the following formula:

[0181]

[0182] where, ω represents the angular frequency, and the calculation formula is as follows:

[0183] ω = 2πf

[0184] Among them, f represents the center frequency of the transmitting antenna;

[0185] ε represents the dielectric constant of the medium, with the unit of F / m; μ represents the magnetic permeability of the medium, with the unit of H / m; σ represents the conductivity of the medium, with the unit of S / m;

[0186] (3-3) Calculation of attenuation:

[0187] Calculate the attenuation B caused by spherical diffusion of propagation at different distances in the set direction N :

[0188] (1) Two-dimensional case:

[0189]

[0190] (2) Three-dimensional case:

[0191]

[0192] Among them, d tr represents the distance between the receiving antenna and the transmitting antenna; d N represents the distance that the electromagnetic wave travels N grid distances from the transmitting antenna, and the calculation formula is as follows:

[0193] d N = NΔx

[0194] (3-4) Calculation of the total attenuation ratio:

[0195] Use the following formula to calculate the total attenuation ratio R at different distances from the transmitting antenna N :

[0196] R N = T1·A N ·B N

[0197] (3-5) Attenuation ratio screening condition:

[0198] Use the following formula to represent the selection condition:

[0199] R N ≤ R0

[0200] Among them, R0 represents the set attenuation ratio of the finite domain boundary;

[0201] (3-6) Calculation of the finite domain boundary:

[0202] As Figure 3 shown is the schematic diagram of the finite domain boundary calculation of the method of the present invention:

[0203] Select points that meet the attenuation ratio screening criteria, and select the point closest to the transmitting antenna from the selected points as the boundary of the finite domain in the set direction;

[0204] (3-7) Calculation of the maximum finite domain boundary:

[0205] Use the following formula to calculate the maximum finite domain boundary corresponding to the recording time length constraint condition in the set direction:

[0206]

[0207] where, t N represents the time for electromagnetic waves to propagate through N grids; i represents the cumulative number of times in the calculation process, starting from 1 and gradually increasing to N; v[[ID=1s]] i represents the electromagnetic wave propagation speed in the vertical direction at a point where the model is i grids away from the transmitting antenna; Δx represents the spatial step size;

[0208] Use the following formula to represent the selection condition:

[0209] t N ≥T / 2

[0210] Select points that meet the above conditions, and select the point closest to the transmitting antenna from the selected points as the maximum finite domain boundary;

[0211] If the finite domain boundary is greater than the maximum finite domain boundary, change the corresponding finite domain boundary to the position of the maximum finite domain boundary;

[0212] (3-8) Determination of the final finite domain range:

[0213] Repeat the above steps (3-2)-(3-7) until the finite domain boundaries in each direction are calculated, and then determine the final finite domain range;

[0214] S4. Using the finite domain range calculated in step S3, calculate the forward modeling parameters of the ground penetrating radar by adding a perfectly matched layer; specifically including:

[0215] Before exciting the transmitting antenna and calculating the field strength components, calculate the parameters required for forward modeling of the finite domain obtained in step S3 according to the known data, including the convolutional perfectly matched layer (CPML) absorption boundary and the finite difference time domain (FDTD) in the forward modeling method;

[0216] (4-1) Calculation of the auxiliary absorption coefficient in the CPML region:

[0217] Before performing the forward simulation, it is necessary to set CPML outside the simulation area and set the parameters inside CPML, specifically including:

[0218] The coordinate stretching factor s is calculated using the following formula i :

[0219]

[0220] where κ i represents the auxiliary absorption coefficient, which is used to improve the absorption effect of PML on surface waves; σ i represents the conductivity parameter in the i direction within the CPML layer; α i represents the auxiliary absorption coefficient, which is used to improve the absorption effect of PML on low-frequency components; j represents the imaginary unit; ω represents the angular frequency; ε0 represents the permittivity of vacuum;

[0221] The thickness of CPML is finite, and σ, κ, and α can vary monotonically in the layer;

[0222] The settings of the three parameters in the z direction are described using the following formula:

[0223]

[0224] where z0 represents the interface between CPML and the simulation area; m represents a parameter, which takes the value of 4 in the method of the present invention; d represents the thickness of CPML; σ max represents the maximum conductivity within CPML, and the calculation formula is as follows:

[0225]

[0226] where ε r represents the relative permittivity of the CPML layer; δ represents the cell size;

[0227]

[0228] where κ max represents the maximum value of κ within CPML, which takes the value of 5 in the method of the present invention;

[0229]

[0230] where α max represents the maximum value of α within CPML, which takes the value of 0.008 in the method of the present invention;

[0231] (4-2) Calculation of forward simulation parameters for the field strength component update equation:

[0232] According to the two-dimensional Maxwell equation and the constitutive equation, the decoupling equation for the two-dimensional time-domain TM mode is described using the following formula:

[0233]

[0234] Where x represents the x-direction in the rectangular coordinate system; y represents the y-direction in the rectangular coordinate system; z represents the z-direction in the rectangular coordinate system; H x H represents the component of the magnetic field strength in the x-direction; y E represents the component of the magnetic field strength in the y-direction. z ε represents the component of the electric field intensity in the z-direction; ε represents the dielectric constant of the medium; t represents time; σ represents the conductivity of the medium; J z The applied current source is represented by μ; the magnetic permeability of the medium is represented by μ.

[0235] Expanding the above formula in the Yee grid yields the update equations for the field strength components, as shown below:

[0236]

[0237] In this definition, the superscript of the field strength indicates time, the subscript indicates the direction of the component, and the part inside the parentheses indicates the position in the grid. This represents the x-direction magnetic field component at position (i,j) at time n+1 / 2; This represents the x-direction magnetic field component at position (i,j) at time n-1 / 2; Let z represent the electric field component in the z-direction at position (i,j) at time n; Let represent the z-direction electric field component at time n (i, j+1); Δy represents the spatial spacing of the grid in the y-direction; CQ(m) represents the medium parameter at the spatial location of the magnetic or electric field on the left side of the above formula, and the calculation formula is shown below:

[0238]

[0239] Where Δt represents the time interval; μ(m) represents the permeability of the medium at position m;

[0240]

[0241] in, This represents the y-direction magnetic field component at position (i,j) at time n+1 / 2; This represents the y-direction magnetic field component at position (i,j) at time n-1 / 2; Δx represents the z-direction electric field component at position (i+1,j) at time n; Δx represents the spatial spacing of the grid in the x-direction.

[0242]

[0243] in, Let z represent the electric field component in the z-direction at position (i,j) at time n+1; This represents the y-direction magnetic field component at position (i+1,j) at time n+1 / 2; This represents the y-direction magnetic field component at position (i,j) at time n+1 / 2; This represents the x-direction magnetic field component at position (i, j+1) at time n+1 / 2; This represents the x-direction magnetic field component at position (i,j) at time n+1 / 2; Let represent the applied current density in the z-direction at position (i,j) at time n+1 / 2; CA(m) represents the coefficients of the update equation, calculated as follows:

[0244]

[0245] Where σ(m) represents the dielectric conductivity at position m; ε(m) represents the dielectric constant at position m;

[0246] CB(m) represents the coefficients of the update equation, and the calculation formula is shown below:

[0247]

[0248] S5. Apply a loading pulse source to the transmitting antenna position set in step S2, and update the electromagnetic field values ​​in the simulated region; specifically including:

[0249] For the antenna position set in step S2, the loading pulse source is represented by the following formula:

[0250]

[0251] Among them, i s ,j s Indicates the position of the transmitting antenna in the grid; This represents the z-direction electric field intensity component at the position of the transmitting antenna at time n+1; This represents the current intensity in the z-direction of the current source at time n+1 / 2;

[0252] The magnetic field strength components are updated using the following formula;

[0253]

[0254]

[0255] The electric field components are updated using the following formula:

[0256]

[0257] S6. Update the electromagnetic auxiliary field in the perfectly matched layer region; specifically, this includes: updating the electromagnetic auxiliary field for CPML, specifically including: (6-1) calculating the auxiliary expression of the magnetic field update equation in CPML:

[0258] The auxiliary expression for the CPML internal magnetic field update equation is calculated using the following formula:

[0259]

[0260]

[0261] Among them, b i and a i The calculation formulas are as follows:

[0262]

[0263]

[0264] (6-2) Calculate the magnetic field within the CPML region:

[0265] The magnetic field within the CPML region is calculated using the following formula:

[0266]

[0267]

[0268] (6-3) Calculate the auxiliary expression for the electric field update equation within CPML:

[0269] The auxiliary expression for the internal electric field update equation of CPML is calculated using the following formula:

[0270]

[0271]

[0272] (6-4) Calculate the electric field within the CPML region:

[0273] The electric field within the CPML region is calculated using the following formula:

[0274]

[0275] S7. Add time steps and repeat steps S5-S6 above until the numerical simulation of the entire time is completed;

[0276] S8. Using the position of the receiving antenna set in step S2, obtain the wave field data of the receiving antenna; specifically including:

[0277] Based on the position of the receiving antenna, output the E value of the receiving antenna position from the full-wave field information. z The data serves as the radar record received by the receiving antenna after the current excitation of the transmitting antenna.

[0278] S9. Repeat steps S3-S8 above until all transmitting antennas have completed the excitation process, and then complete the finite-domain forward modeling calculation of the ground penetrating radar to obtain the radar profile.

[0279] The method of this invention is described in detail through the following sets of numerical experiments, specifically including:

[0280] 1. Numerical experiment to explore the optimal range of values ​​for the boundary attenuation ratio R0 of a finite field:

[0281] To explore the optimal value of the boundary attenuation ratio R0 in a finite field to balance accuracy and speed requirements, a model is established... Figure 4 Numerical experiments were conducted on the 1m×5m undulating interface model shown to test R0. Figure 4 The diagram shown is a schematic representation of the relative permittivity and conductivity of Model 1 in the method of this invention.

[0282] The simulation area has a grid size of 200×1000, a spatial sampling interval of 0.005m, and a 10-layer CPML absorbing boundary. The transmitting antenna is 0.05m above the ground, with 197 observation points evenly distributed and a transmit / receive distance of 8cm. A 900MHz Ricker wavelet is used as the pulse source, the recording time T is 16ns, and the time step is 0.01ns. The traditional FDTD forward modeling takes 1270s.

[0283] Finite-domain forward modeling experiments were conducted using different values ​​of R0, increasing the value of R0 from 0.04 to 0.25, for a total of 85 experiments. The time taken for each finite-domain forward modeling was recorded, and the results were processed to remove direct waves before calculating the relative error compared to traditional FDTD forward modeling records. The results are as follows: Figure 5 As shown, Figure 5 The figure shown is a schematic diagram of the forward modeling time and relative error curves of Model 1 under different attenuation ratios (R0) in the method of the present invention;

[0284] exist Figure 5In the figure, the solid line represents the curve of relative error in finite field forward modeling as a function of R0, and the dashed line represents the curve of forward modeling runtime as a function of R0. As R0 increases, the finite field range decreases, the algorithm runtime decreases in a near-exponential manner, and the relative error increases in a near-exponential manner. To balance the accuracy and speed requirements of the algorithm, the interval shown by the dashed line in the figure is truncated. On the left side of the interval, the relative error of the algorithm is very small, but the runtime is too long, and the algorithm acceleration effect is very poor. On the right side of the interval, the runtime of the algorithm decreases slowly, but the error increases sharply, which cannot meet the accuracy requirements of forward modeling. Considering both accuracy and speed, R0 (0.09~0.19) in the middle of the dashed line in the figure is selected as a suitable attenuation ratio.

[0285] Different application scenarios require antennas of different frequencies for detection. Using the same method, we conducted multiple experiments on the range of R0 values ​​at several frequencies, and obtained the results shown in the table below:

[0286] Table 1 Optimal range of R0 values ​​at several antenna frequencies

[0287] 100MHz 0.01～0.08 200MHz 0.025～0.1 400MHz 0.06～0.15 900MHz 0.09～0.19

[0288] 2. Numerical experiment with time constraints:

[0289] To verify the effect of the recording time length constraint on the finite field size constraint, a system is established as follows: Figure 6 The 1m×5m road defect model shown is as follows: Figure 6 The diagram shows the relative permittivity and conductivity of Model 2 in the method of this invention: air layer 0.05m, underground medium depth 0.95m. Figure 6 (a) is a schematic diagram of the relative permittivity distribution of the model. Figure 6 (b) is a schematic diagram of the electrical conductivity (S / m) distribution of the model, from left to right representing three types of road defects: cracks (0.75m), cavities (2.15m), and subsidence (3.65m to 3.93m);

[0290] A 900MHz Ricker wavelet was used as the pulse source. The transmitting antenna was 0.05m above the ground, and 200 observation points were evenly distributed. The recording time T was 15ns, the time step was 0.01ns, the grid size of the simulation area was 200×1000, and the spatial step was 0.005m. Conventional FDTD forward modeling and FD-FDTD forward modeling constrained by the recording time were used, with the attenuation ratio R0 set to 0. Figure 7 The figure shows the forward modeling results of Model 2 in the method of this invention and a schematic diagram of the finite field forward modeling error. Figure 7 (a) is a schematic diagram of the forward modeling results; Figure 7 (b) is a schematic diagram showing the error between the forward modeling results and the conventional forward modeling results;

[0291] 1) Error manifestation:

[0292] The radar records obtained by FD-FDTD and FDTD are very similar. Figure 7 In (b), the FD-FDTD with recording time constraints only has some error in the larger recording time range, with a maximum error of 0.0656. The overall error distribution is small and the amplitude is small, which can be ignored relative to the amplitude of the effective wave (35.69). To investigate the error of the radar response to anomalies in the data, the relative error was calculated after removing the direct wave for both. The relative error was calculated using the following formula:

[0293]

[0294] Where ||·|| denotes norm calculation; R fd R represents radar records used in finite-field forward modeling. tra This represents the radar record from the conventional forward modeling; based on the above formula, the calculated result is 0.00934%.

[0295] 2) Cost analysis:

[0296] Two-dimensional forward modeling only calculates a finite domain truncation simulation region in the horizontal direction, using the number of grids in the horizontal direction (x-direction) of the simulation region to simply replace the memory representation; FDTD requires wavefield calculation for the entire model in each forward modeling process, with 1000 grids in the horizontal direction; FD-FDTD, which records the time length constraint, has an average of 463 horizontal grids in the simulated region of the antenna each time, effectively reducing the scale of forward modeling; in terms of running time, FD-FDTD takes only 297s, far less than FDTD forward modeling takes 1206s, achieving a speedup of 4.06;

[0297] The experimental results above show that when R0 is too small, the recording time constraint can prevent the finite field from being too large. The constraint calculates the finite field with the maximum range, which can better balance the accuracy and speed requirements and ensure the acceleration performance of the finite field method.

[0298] 3. Error Source Analysis Numerical Experiment:

[0299] To investigate the sources of error, the 94th error, which has a relatively large error in the forward modeling results of Model 2, was selected for error source analysis.

[0300] like Figure 8 The diagram shown is a comparison of the conventional forward modeling and the finite domain forward modeling wavefield of the 94th trace of Model 2 in the method of this invention: Figure 8The image shows wavefield snapshots of the traditional FDTD at t = 5 ns, t = 7.5 ns, and t = 9 ns, as well as wavefield snapshots of the finite-domain FDTD. At t = 7.5 ns, the wave from the traditional FDTD just reaches the collapse anomaly on the right side of the model and generates a diffracted wave. However, this anomaly does not exist in the forward modeling region of the finite-domain FDTD, and the diffracted wave is not recorded in the receiving antenna. This is an important reason for the finite-domain forward modeling error of the 94th channel.

[0301] The above experiments demonstrate that the main source of error in finite-domain forward modeling is the phenomenon that the actual propagation distance of radar waves exceeds the finite-domain range. When the transmitting antenna is a certain distance from the ground surface, the spherical diffusion characteristic of electromagnetic waves will cause some electromagnetic waves to propagate a certain distance in the air and then penetrate into the underground medium. However, the energy of this portion of electromagnetic waves is very small, and the energy of the diffracted waves generated by anomalies in distant regions is almost completely attenuated when they reach the receiving antenna. Moreover, in actual detection, radar antennas are mostly shielded antennas, and the transmitted electromagnetic pulses are transmitted almost vertically downwards to the ground surface. The received signals are all reflected echoes from underground. Therefore, when using the finite-domain method for forward modeling, the error caused by the small portion of echoes that propagate a certain distance in the air and then enter the underground medium outside the finite domain can be ignored.

[0302] Furthermore, the visible errors at the upper part of the simulated region in the wavefield snapshots at 7.5 ns and 9 ns are due to the PML absorption error caused by the process of loading PML after truncating the finite domain. Combining the error pattern with the difference pattern of the wavefield snapshots, it can be inferred that the error of FD-FDTD compared to FDTD mainly lies in the difference in the absorption process of PML between the two. FD-FDTD absorbs the wave earlier than FDTD, but this error is only calculated relative to the forward modeling result of FDTD. Since there is also a certain error in the absorption of PML during the forward modeling process of FDTD, and this error value is very small, does not belong to the effective detection information, and has almost no interference to the effective wave, it can be ignored.

[0303] 4. Finite-field forward modeling experiment for large-scale models:

[0304] A publicly available high-resolution three-dimensional hydrological facies dataset was used, sourced from a well-maintained gravel quarry near Herten (southwest Germany); such as Figure 9 The diagram shown is of model 3 in the method of this invention: the model size is 16m × 10m × 7m, and the physical property parameters of the slices at positions X = 0, 4, 8, 12, and 16m are as follows. Figure 9 As shown; Figure 9 This is a schematic diagram of model 3 in the method of the present invention; Figure 9 (a) is a schematic diagram of the relative permittivity distribution of the model; Figure 9(b) is a schematic diagram of the conductivity distribution of the model; the grid size of the model is 639×401×288, and the spatial step is 0.025m; a 25cm thick air layer is established above the model; three measurement lines are established at x=4, 8, and 12m, with 96 measurement points evenly distributed on each measurement line, 0.1m above the ground surface, using a 100MHz Ricker wavelet as the excitation source, with a transmit / receive distance of 0.3m; the time step is 0.025ns, the recording time is set to 200ns, and the total time is 8000 steps;

[0305] As numerical simulations demand higher precision, finer meshes, and larger computational domains, the computational efficiency of CPUs, which primarily rely on serial computation, is gradually becoming insufficient to meet practical needs. GPUs, with their greater number of stream processors, are well-suited for computing large numbers of relatively simple tasks and are currently the mainstream acceleration method. However, GPU memory is significantly limited compared to CPU memory. For example, a personal computer with a 12th Gen Intel® Core™ i7-12700 CPU and an NVIDIA GeForce GTX 1650 SUPER GPU has only 4GB of video memory. If this large-scale model is performed on this computer, with the addition of 10 CPML layers, the memory required to store the wave field and model parameters alone reaches 5.73GB, far exceeding the computer's capacity.

[0306] The finite field algorithm can solve the above problems well; with R0 = 0.04, the forward modeling calculation of a large three-dimensional model of three survey lines was completed on a microcomputer in 244 hours; Figure 10 This is a schematic diagram showing the relative permittivity distribution profiles of the three measurement lines of Model 3 in the method of this invention and the corresponding finite domain forward modeling results; Figure 10 (a) is a schematic diagram of the relative permittivity profile of the survey line at x = 4m; Figure 10 (b) is a schematic diagram of the relative permittivity profile of the survey line at x = 8m; Figure 10 (c) is a schematic diagram of the relative permittivity profile of the survey line at x = 12m; Figure 10 (d) is a schematic diagram of the finite-domain forward modeling results of the survey line at x = 4m; Figure 10 (e) is a schematic diagram of the finite-domain forward modeling results of the survey line at x = 8m; Figure 10 (f) is a schematic diagram of the finite domain forward modeling results of the survey line at x = 12m; in the ground penetrating radar forward modeling profile, the interface information between the layers of the model is perfectly reflected, and the deep reflection wave information is consistent with the characteristics of the deep part of the model containing a large amount of sand and gravel;

[0307] To compare the speed difference between the FD-FDTD and FDTD models, FD-FDTD forward modeling was performed on a workstation using an Intel(R) Xeon(R) Platinum 8168 CPU and an NVIDIA Quadro P6000 GPU. The computational memory and time statistics for both are shown in the table below. Conventional forward modeling calculates the wavefield for the entire model, requiring approximately 15.63 GB of computational memory per channel, while FD-FDTD requires an average of only 3.36 GB per channel, reducing the average computational memory requirement to 0.21 GB, allowing for more parallel threads. This significant reduction in computational scale enabled FD-FDTD to complete the forward modeling of three survey lines in just 79 hours; whereas FDTD takes 7 hours to calculate one observation point, and calculating the 288 observation points across three survey lines is estimated to require 2016 hours.

[0308]

[0309] The above experiments verified the powerful advantages of the finite field algorithm in forward modeling of large models. Performing wave field calculations within a finite field reduced the average memory usage for a single excitation to about 1 / 5, allowing for more threads to perform calculations, and ultimately increasing the running speed by 25.5 times. The significant reduction in memory usage enabled the successful completion of forward modeling of large 3D models on microcomputers that could not handle the memory requirements of conventional FDTD calculations.

Claims

1. A ground-penetrating radar forward modeling method based on a finite domain, comprising the following steps: S1. Construct a finite-domain forward model for ground-penetrating radar; S2. Set the positions of the ground-penetrating radar transmitting and receiving antennas, the center frequency of the transmitting antenna, and the finite domain boundary attenuation ratio; S3. Using the antenna position and center frequency set in step S2, and the physical property parameters of the forward model constructed in step S1, calculate the finite domain range corresponding to the transmitting antenna; specifically including: (3-1) Calculation of transmission coefficient: The transmission coefficient of radar waves entering perpendicularly below the transmitting antenna is calculated using the following formula. : in, Indicates the relative permittivity of air; This represents the relative permittivity of the stratum beneath the transmitting antenna; (3-2) Calculation of the dielectric absorption attenuation ratio: Starting from the transmitting antenna, the attenuation ratio caused by medium absorption at different distances propagating in a set direction is calculated using the following formula. : in, This indicates the attenuation ratio caused by medium absorption at different distances as the radar wave propagates in a set direction; This represents the spatial step size of the grid; if a point on the Earth's surface is a distance from the transmitting antenna... If there are 1 grid, then the propagation distance of the electromagnetic wave is equal to 1 / 2 grid. ; Indicates the distance between the model and the transmitting antenna. The minimum attenuation coefficient of each grid position along the vertical direction is calculated using the following formula: in, The angular frequency is expressed by the following formula: in, Indicates the center frequency of the transmitting antenna; The dielectric constant of a medium, expressed in units of 1000 kJ / m². ; The magnetic permeability of a medium is expressed in units of 1000 ppm. ; The conductivity of a medium is expressed in units of 1000 kJ / m². ; (3-3) Calculation of spherical diffusion attenuation: Calculate the attenuation caused by spherical diffusion at different distances propagating in a given direction. : (1) Two-dimensional case: (2) Three-dimensional situation: in, This indicates the distance between the receiving antenna and the transmitting antenna; This indicates that electromagnetic waves propagate from the transmitting antenna. The distance between grid cells is calculated using the following formula: (3-4) Calculation of the total attenuation ratio: The total attenuation ratio at different propagation distances from the transmitting antenna is calculated using the following formula. : (3-5) Attenuation ratio screening criteria: The selection criteria are expressed by the following formula: in, This represents the set attenuation ratio at the boundary of the finite field; (3-6) Calculation of the boundary of a finite field: Select points that meet the attenuation ratio screening conditions, and choose the point closest to the transmitting antenna from the selected points as the boundary of the finite domain in the set direction; (3-7) Calculation of the boundary of the largest finite field: The maximum finite field boundary corresponding to the time length constraint in the set direction is calculated using the following formula: in, Indicates electromagnetic wave propagation The time for each grid; This indicates the number of accumulations during the calculation process, from... Start by taking values, gradually increasing to ; Indicates the distance between the model and the transmitting antenna. The electromagnetic wave propagation speed at each point in the grid in the vertical direction; Indicates the spatial step size; The selection criteria are expressed by the following formula: Select points that satisfy the above conditions, and choose the point closest to the transmitting antenna from the selected points as the boundary of the maximum finite field; If the boundary of a finite field is greater than the boundary of the largest finite field, the corresponding boundary of the finite field is changed to the position of the boundary of the largest finite field. (3-8) Determining the final range of the finite field: Repeat steps (3-2)-(3-7) above until the finite field boundary in each direction is calculated, and then determine the final finite field range; S4. Using the finite domain range calculated in step S3, calculate the finite domain forward modeling parameters of ground penetrating radar by adding a perfect matching layer; S5. Apply a loading pulse source to the transmitting antenna position set in step S2, and update the electromagnetic field value of the simulated area at the same time; S6. Update the electromagnetic auxiliary field in the perfectly matched layer region; S7. Add time steps and repeat steps S5-S6 above until the numerical simulation of the entire time is completed; S8. Using the position of the receiving antenna set in step S2, obtain the wave field data of the receiving antenna; S9. Repeat steps S3-S8 above until all transmitting antennas have completed the excitation process, thereby completing the finite-domain forward modeling calculation of the ground penetrating radar and obtaining the radar profile.

2. The ground-penetrating radar forward modeling method based on a finite domain according to claim 1, characterized in that... Step S1, which involves constructing a finite-domain forward model for ground-penetrating radar, specifically includes: Establish a finite-domain forward model for ground-penetrating radar based on application requirements; The physical properties of the forward model include relative permittivity and conductivity.

3. The ground-penetrating radar forward modeling method based on a finite domain according to claim 2, characterized in that... Step S2, which involves setting the positions of the ground-penetrating radar transmitting and receiving antennas, the center frequency of the transmitting antenna, and setting the finite-domain boundary attenuation ratio, specifically includes: Setting the location of the GPR transmitting antenna specifically includes: Set the number, starting position, spacing, and ending position of the transmitting antennas; Setting the GPR receiving antenna position specifically includes: Set the number, start position, spacing, and end position of the receiving antennas, as well as the transmit / receive distance between the transmitting and receiving antennas; Setting the center frequency of the transmitting antenna, specifically including: The higher the center frequency of the transmitting antenna, the faster the electromagnetic waves attenuate and the shallower the detection depth. Different transmit antenna frequencies need to be selected based on the model depth. The frequency is determined using the following formula: in, Indicates the center frequency of the radar transmitting antenna; Indicates spatial resolution; Represents the relative permittivity of the surrounding rock; Finite domain boundary attenuation ratio Affects the acceleration effect. The larger the value, the better the acceleration effect, but the misunderstanding error also increases accordingly. Experiments were conducted to determine the antenna frequencies at different frequencies. The optimal range of values ​​for includes: When the antenna frequency is 100MHz The optimal value range is 0.01 to 0.08; When the antenna frequency is 200MHz The optimal value range is 0.025~0.1; When the antenna frequency is 400MHz The optimal value range is 0.06 to 0.15; When the antenna frequency is 900MHz The optimal value range is 0.09 to 0.

19.

4. The ground-penetrating radar forward modeling method based on a finite domain according to claim 3, characterized in that... Step S4, which involves using the finite domain range calculated in step S3 and adding a perfect matching layer to calculate the finite domain forward modeling parameters of ground-penetrating radar, specifically includes: Before exciting the transmitting antenna and calculating the field strength components, the parameters required for forward modeling of the finite field obtained in step S3 are calculated based on the known data, including the absorption boundary of the fully matched convolutional layer CPML and the finite difference method FDTD in the forward modeling method. Specifically, this includes the calculation of the auxiliary absorption coefficient within the CPML region and the calculation of the forward modeling parameters of the field strength component update equation.

5. The ground-penetrating radar forward modeling method based on a finite domain according to claim 4, characterized in that... The calculation of the auxiliary absorption coefficient within the CPML region specifically includes: Before performing forward modeling, it is necessary to configure CPML outside the simulation region and set the parameters inside CPML, specifically including: The coordinate scaling factor is calculated using the following formula. : in, This represents the auxiliary absorption coefficient, used to improve the absorption effect of PML on surface waves; Indicates the CPML layer directional conductivity parameters; This represents the auxiliary absorption coefficient, used to improve the absorption effect of PML on low-frequency components; Represents the imaginary unit; Indicates angular frequency; The dielectric constant of vacuum; CPML has a limited thickness. , , It can change monotonically within a layer; The following formula describes the situation in In terms of direction, three parameters need to be set: in, This represents the interface between CPML and the simulation region; Indicates parameters; Indicates the thickness of CPML; The maximum conductivity within the CPML is expressed by the following formula: in, Indicates the relative permittivity of the CPML layer; Indicates cell size; in, Indicates the content within CPML The maximum value; in, Indicates the content within CPML The maximum value.

6. The ground-penetrating radar forward modeling method based on a finite domain according to claim 5, characterized in that... The calculation of the forward modeling parameters of the field strength component update equation includes: Based on the two-dimensional Maxwell equations and constitutive equations, the decoupling equations for the two-dimensional time-domain TM mode are described by the following formulas: Among them, the department direction; Indicates the magnetic field strength at Components in direction; Indicates the magnetic field strength at Components in direction; Indicates the electric field strength at Components in direction; Indicates the dielectric constant of the medium; Indicates time; Indicates the electrical conductivity of the medium; Indicates the applied current source; Indicates the magnetic permeability of the medium; Expanding the above formula in the Yee grid yields the update equations for the field strength components, as shown below: In this definition, the superscript of the field strength indicates time, the subscript indicates the direction of the component, and the part inside the parentheses indicates the position in the grid. express time Location Directional magnetic field components; express time Location Directional magnetic field components; express time Location Directional electric field components; express time Location directional electric field components; Indicates the grid in Spatial interval in direction; The medium parameter representing the spatial location of the magnetic or electric field on the left side of the above formula is calculated as follows: in, Indicates a time interval; express The permeability of the medium at a given location; in, express time Location Directional magnetic field components; express time Location Directional magnetic field components; express time Location directional electric field components; Indicates the grid in Spatial interval in direction; in, express time Location directional electric field components; express time Location Directional magnetic field components; express time Location Directional magnetic field components; express time Location Directional magnetic field components; express time Location Directional magnetic field components; express time Location The direction of the applied current density; The coefficients of the updated equation are calculated using the following formula: in, express The dielectric conductivity at the location; express The dielectric constant of the medium at a given location; The coefficients of the update equation are represented by the following formula: 。 7. The ground-penetrating radar forward modeling method based on a finite domain according to claim 6, characterized in that... Step S5, which involves loading a pulse source for the antenna position set in step S2 and updating the electromagnetic field values ​​in the simulated region, specifically includes: For the antenna position set in step S2, the loading pulse source is represented by the following formula: in, Indicates the position of the transmitting antenna in the grid; express The position of the transmitting antenna at all times directional electric field intensity components; express The current source at all times Directional current intensity; The magnetic field strength components are updated using the following formula; The electric field components are updated using the following formula: 。 8. The ground-penetrating radar forward modeling method based on a finite domain according to claim 7, characterized in that... Step S6, which involves updating the electromagnetic auxiliary field in the perfectly matched layer region, specifically includes: The electromagnetic field auxiliary field update process for CPML includes: (6-1) Calculate the auxiliary expression for the CPML internal magnetic field update equation: The auxiliary expression for the CPML internal magnetic field update equation is calculated using the following formula: in, and The calculation formulas are as follows: (6-2) Calculate the magnetic field within the CPML region: The magnetic field within the CPML region is calculated using the following formula: (6-3) Calculate the auxiliary expression for the internal electric field update equation of CPML: The auxiliary expression for the internal electric field update equation of CPML is calculated using the following formula: (6-4) Calculate the electric field within the CPML region: The electric field within the CPML region is calculated using the following formula: 。 9. The ground-penetrating radar forward modeling method based on a finite domain according to claim 8, characterized in that... Step S8, which involves using the position of the receiving antenna set in step S2 to obtain the wave field data of the receiving antenna, specifically includes: Based on the position of the receiving antenna, output the location of the receiving antenna from the full-wave field information. The data serves as the radar record received by the receiving antenna after the current excitation of the transmitting antenna.

Images