Efficient computing method and system for ground penetrating imaging based on CM24 high-order time domain finite difference and high-order PML

Abstract

The application provides a ground penetrating imaging efficient calculation method and system based on CM24 high-order time domain finite difference and high-order PML, relates to the technical field of geophysical exploration and electromagnetic wave imaging, and the method comprises the following steps: performing forward wave field propagation calculation on a wide-band corrected electromagnetic field update equation to obtain forward simulation data of each receiving point; according to the forward simulation data, time-reversing recorded echo data and then loading the time-reversed echo data to corresponding receiving point positions as excitation sources to perform reverse wave field propagation calculation and obtain reverse propagation wave fields; and according to the forward propagation wave fields and the reverse propagation wave fields, multiplying the forward wave field values and the reverse wave field values at each time step and accumulating the product results of all time steps to generate ground penetrating radar reverse time migration imaging results reflecting underground medium discontinuous interfaces and target positions. The application solves the problems of performance decline of a traditional absorbing boundary under a coarse grid condition and introduction of an additional reflection by the absorbing boundary.

Description

Technical Field

[0001] This invention relates to the field of structural health monitoring and simulation analysis technology for hydraulic engineering, and in particular to an efficient calculation method and system for ground-penetrating imaging based on CM24 high-order time-domain finite difference and high-order PML. Background Technology

[0002] In a ground-penetrating radar (GPR) detection scenario, a transmitting antenna, a receiving antenna, and a multi-layered underground medium are involved. This multi-layered medium includes an air layer, dielectric layers with different dielectric parameters, and embedded targets. The transmitting antenna emits broadband pulsed electromagnetic waves into the underground medium. These electromagnetic waves generate reflected echoes at the interfaces between the dielectric layers and at the embedded targets, which are then received by the receiving antenna. This detection process involves calculations of the forward and reverse propagation of electromagnetic waves within a continuous space.

[0003] Traditional finite-difference time-domain methods for simulating electromagnetic wave propagation are subject to strict limitations on mesh size due to numerical dispersion problems. To ensure phase accuracy, the mesh size is typically required to be no larger than 1 / 20 of the wavelength. This limitation leads to a significant increase in memory consumption and low computational efficiency in simulations of large-area detection scenarios.

[0004] Existing high-order temporal finite difference methods show a significant decrease in absorption performance at the absorption boundaries of conventional perfectly matched convolutional layers under coarse-grid resolution. Using a fine-grid in the absorption boundary region introduces additional reflection interference at the interface between the coarse and fine grids, further affecting the final reverse-time migration imaging quality.

[0005] Furthermore, existing high-order finite-difference time-domain methods typically employ coefficient selection strategies aimed at minimizing numerical dispersion error at a single frequency point. Since ground-penetrating radar employs broadband pulse excitation, coefficient selection at a single frequency point is insufficient to minimize the overall error across the entire operating frequency band, resulting in inadequate phase accuracy in broadband electromagnetic wave propagation simulations. Summary of the Invention

[0006] This invention provides an efficient calculation method and system for ground-penetrating imaging based on CM24 high-order finite-difference time-domain and high-order PML, which solves the problems of performance degradation of traditional absorbing boundaries and additional reflection introduced by the interface under coarse grid conditions.

[0007] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows: Firstly, an efficient calculation method for ground-penetrating imaging based on CM24 high-order finite-difference time-domain and high-order PML is provided, the method comprising: Based on the detection scenario, the calculation area, spatial grid size, dielectric parameters of the medium, excitation source type and center frequency are set, and the positions of the transmitting and receiving antennas are determined to establish a calculation model of the ground penetrating radar detection area. Based on the computational model, an electromagnetic field update equation containing fourth-order spatial difference is constructed to obtain a computational framework for low-numerical dispersive wave field propagation under coarse-grid resolution conditions. Based on the aforementioned computational framework and the spectral distribution characteristics of the ground-penetrating radar, coefficients that minimize the overall numerical dispersion error within the operating frequency band are selected to generate the broadband corrected electromagnetic field update equation. Based on the broadband corrected electromagnetic field update equation, by introducing a higher-order scaling transformation factor, a higher-order convolution fully matched layer absorption boundary condition is constructed to generate a truncated boundary absorption model under coarse grid conditions. Based on the truncated boundary absorption model, the forward wave field propagation calculation is performed on the broadband corrected electromagnetic field update equation to obtain forward modeling data for each receiving point. Based on the forward modeling data, the recorded echo data is reversed in time and reloaded as an excitation source to the corresponding receiving point position to perform reverse wave field propagation calculation and obtain the reverse propagation wave field. Based on the forward propagation wavefield and the reverse propagation wavefield, at each time step, the forward wavefield value is multiplied by the reverse wavefield value after time inversion, and the product results of all time steps are accumulated to generate a ground-penetrating radar reverse time migration imaging result that reflects the discontinuous interface of the underground medium and the target position.

[0008] Secondly, a high-efficiency computing system for ground-penetrating imaging based on CM24 high-order finite-difference time-domain and high-order PML includes: The model building module is used to set the calculation area, spatial grid size, medium electromagnetic parameters and center frequency according to the detection scenario, and to determine the position of the transmitting and receiving antennas to build a ground penetrating radar calculation model. The equation construction module is used to construct an electromagnetic field update equation containing fourth-order spatial difference based on the computational model, and obtain a computational framework for low numerical dispersive wave field propagation. The coefficient optimization module is used to select the coefficients that minimize the overall numerical dispersion error within the working frequency band based on the computational framework and the characteristics of the working spectrum distribution, and to generate the broadband corrected electromagnetic field update equation. The boundary construction module is used to construct the high-order convolution fully matched layer absorption boundary conditions based on the broadband modified electromagnetic field update equation and through a high-order scaling transformation factor, thereby generating a truncated boundary absorption model. The wavefield calculation and imaging module is used to calculate the forward wavefield propagation based on the truncated boundary absorption model to obtain forward modeling data. The data is then reversed in time and used as an excitation source to reload the data to the receiving point to calculate the reverse wavefield propagation to obtain the reverse propagation wavefield. Finally, the forward and reverse propagation wavefields are multiplied at each time step and the product of all time steps is accumulated to generate the ground penetrating radar reverse time migration imaging result.

[0009] The above-described solution of the present invention has at least the following beneficial effects: Based on the spectral distribution characteristics of ground-penetrating radar, coefficients that minimize the overall numerical dispersion error within the operating frequency band are selected to generate a broadband-corrected electromagnetic field update equation. This operation enables high-precision electromagnetic wave propagation simulation under broadband pulse excitation, improves broadband phase accuracy, and solves the problem of insufficient broadband phase accuracy caused by selecting a single frequency point.

[0010] By employing a high-order scaling transform factor, a high-order convolutional perfectly matched layer absorbing boundary condition is constructed. Under coarse-grid resolution, this boundary condition effectively enhances the absorption of evanescent and low-frequency wave components, significantly suppresses reflection interference at truncated boundaries, and solves the problems of performance degradation of traditional absorbing boundaries and the introduction of additional reflections at interfaces under coarse-grid conditions, ensuring the accuracy and stability of forward modeling and reverse-time migration imaging.

[0011] By combining wideband coefficient selection with higher-order absorption boundary conditions, high-quality ground-penetrating radar forward modeling and reverse-time migration imaging were achieved under coarse-grid conditions with a 1 / 4 wavelength grid resolution. This method significantly improves computational efficiency and reduces memory consumption while maintaining good imaging quality, overcoming the problems of excessively high grid resolution requirements and large computational resource consumption associated with traditional methods. Attached Figure Description

[0012] Figure 1 This is a flowchart of an efficient calculation method for ground-penetrating imaging based on CM24 high-order temporal finite difference and high-order PML provided by an embodiment of the present invention.

[0013] Figure 2 The electric field component in the CM24 method A schematic diagram of the field components involved in the update.

[0014] Figure 3 The embodiments of the present invention provide the transmitted waveform spectral distribution and the corresponding frequency numerical dispersion error.

[0015] Figure 4 Wideband dispersion error diagrams corresponding to different coefficients at a spatial resolution of 4.

[0016] Figure 5 The reflection error map of the CM24-HO-CPML proposed in this invention at the observation point location.

[0017] Figure 6 Imaging model diagram.

[0018] Figure 7 The RTM imaging results of this invention.

[0019] Figure 8A comparison of the RTM imaging results and CM24-CPML imaging results of this invention. Detailed Implementation

[0020] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0021] like Figure 1 As shown, this embodiment provides an efficient calculation method for ground-penetrating imaging based on CM24 high-order finite-difference time-domain and high-order PML, including the following steps: Step S1: Based on the detection scenario, set the calculation area, spatial grid size, dielectric parameters of the medium, excitation source type and center frequency, determine the position of the transmitting and receiving antennas, and establish a calculation model of the ground penetrating radar detection area.

[0022] Step S2: Based on the computational model, construct the electromagnetic field update equation containing fourth-order spatial difference to obtain the computational framework for low-numerical dispersive wave field propagation under coarse grid resolution conditions.

[0023] Step S3: Based on the calculation framework and the operating spectrum distribution characteristics of the ground penetrating radar, select the coefficients that minimize the overall numerical dispersion error within the operating frequency band, and generate the broadband corrected electromagnetic field update equation.

[0024] Step S4: Based on the broadband corrected electromagnetic field update equation, by introducing a higher-order scaling transformation factor, construct a higher-order convolution fully matched layer absorption boundary condition, and generate a truncated boundary absorption model under coarse mesh conditions.

[0025] Step S5: Based on the truncated boundary absorption model, perform forward wave field propagation calculations on the broadband corrected electromagnetic field update equations to obtain forward simulation data for each receiving point.

[0026] Step S6: Based on the forward simulation data, the recorded echo data is reversed in time and reloaded to the corresponding receiving point position as the excitation source to perform reverse wave field propagation calculation and obtain the reverse propagation wave field.

[0027] Step S7: Based on the forward propagation wave field and the reverse propagation wave field, multiply the forward wave field value with the reverse wave field value after time inversion at each time step, and accumulate the product results of all time steps to generate a ground-penetrating radar reverse time migration imaging result that reflects the discontinuous interface of the underground medium and the target position.

[0028] The beneficial effects of this embodiment are as follows: Step S2 constructs an electromagnetic field update equation containing fourth-order spatial difference. Combined with step S3, which selects coefficients based on the spectral distribution characteristics of the ground-penetrating radar (GPR) to minimize the overall numerical dispersion error within the working frequency band, a broadband-corrected electromagnetic field update equation is generated. This operation achieves high-precision electromagnetic wave propagation simulation under broadband pulse excitation, effectively improving phase accuracy within the broadband band and solving the problem of insufficient broadband phase accuracy caused by the single-frequency-point coefficient selection strategy in existing technologies. Step S4 constructs a higher-order scaling transform factor and a higher-order convolution fully matched layer absorption boundary condition, effectively enhancing the absorption capability of evanescent waves and low-frequency wave components under coarse-grid resolution conditions. This significantly suppresses reflection interference at truncated boundaries, solving the problems of decreased performance of traditional CPML absorption boundaries under coarse-grid conditions and the introduction of additional reflections at the coarse-fine grid interface. Through the combination of the above broadband correction and higher-order absorption boundary conditions, high-quality GPR forward modeling and reverse-time migration imaging are achieved under coarse-grid conditions with a 1 / 4 wavelength grid resolution. While maintaining good imaging quality, computational efficiency is significantly improved, and memory consumption is significantly reduced. In practical applications, the specific implementation process is as follows: Step S1: Establish a calculation model of the ground-penetrating radar detection area: A computational model of the ground-penetrating radar detection area was established to provide a basic computational framework for subsequent electromagnetic field update equation construction and wave field propagation calculations. The specific implementation process is as follows: Step S11: Based on the detection scenario, set the path along... direction, direction and The three-dimensional computational area size along the direction. In a typical application scenario of health monitoring of hydraulic engineering structures, the ground-penetrating radar detection area is along... The length of the direction (horizontal survey line direction) is set according to the detection requirements, for example, 2.0 m; along The width of the direction (horizontal direction perpendicular to the survey line) is determined based on the antenna's movement range and the target's distribution range, for example, set to 0.5 m; along The height of the direction (depth direction) needs to cover the entire area from the surface air layer to below the depth of the target, for example, set to 1.5 m. Set it at the periphery of the calculation area. Layers (e.g., 8 to 12 layers) of high-order convolutions perfectly match the layer absorption boundaries, which surround the entire computational region to achieve absorption truncation of outgoing traveling waves.

[0029] Step S12: Based on the 3D computational region dimensions set in step S11, set the electromagnetic parameters of the air layer and different dielectric layers, and determine the position and shape of the embedded target. The air layer is located at the top of the computational region, with a thickness of, for example, 0.2 m, and its relative permittivity... Take 1.0, conductivity Take 0.0 S / m. The underground medium layer is set according to the actual geological conditions: the relative permittivity of the first medium layer (such as sandy soil layer). The conductivity is between 6.0 and 9.0. Take values ​​from 0.001 to 0.01 S / m; the relative permittivity of the second dielectric layer (such as a clay layer). Take values ​​from 12.0 to 20.0, conductivity The relative permittivity is set to 0.01 to 0.05 S / m. The embedded target is located in the underground medium layer and can be rectangular, circular, or irregular in shape. and conductivity Based on the target material (e.g., cavity extraction) Metals are used for high conductivity, concrete is used for... (Up to 8.0) confirmed.

[0030] Step S13: Based on the above-mentioned electromagnetic parameters of the medium and the embedded target characteristics, set the spatial grid size. This embodiment uses a 1 / 4 wavelength grid resolution, i.e., the spatial grid step size. , and All are set to one-quarter of the minimum wavelength in the medium. The center frequency is... (For example, 500MHz), the corresponding wavelength in vacuum is (in The speed of light in a vacuum (m / s), minimum wavelength in the medium , where is the wavelength of the maximum frequency with a magnitude of -40 dB in the normalized spectrum of the radiation field within the calculation region. Spatial grid step size. Pick ,in R In this embodiment, the number of grids per minimum wavelength is set as follows. (i.e., 1 / 4 wavelength resolution). Simultaneously, the time step is determined to satisfy the Courant-Friedrichs-Lewy (CFL) stability condition.

[0031] Step S14: Determine the positions of the transmitting and receiving antennas and the measurement point spacing in bistatic radar observation mode. The transmitting antenna Tx and the receiving antenna Rx are arranged along the measurement line with a fixed offset, and the measurement points are along... direction The antennas are evenly distributed. Both transmitting and receiving antennas are placed in the air layer, at a height of, for example, 0.05 m above the ground surface. After all antenna parameters are set, a complete calculation model of the ground-penetrating radar detection area is established.

[0032] Through steps S11 to S14 above, this step establishes a complete ground-penetrating radar detection calculation model that includes the three-dimensional calculation region, medium distribution, embedded target, grid resolution, and transceiver antenna configuration, providing an accurate spatial discretization basis for the subsequent construction of electromagnetic field update equations.

[0033] Step S2, construct the CM24 electromagnetic field update equations containing fourth-order spatial differences: like Figure 2 As shown, in the differential template structure of the CM24 method, the electric field... E z Taking the component as an example, a CM24 electromagnetic field update equation containing fourth-order spatial difference is constructed to obtain a wave field propagation calculation framework with low numerical dispersion characteristics under coarse grid resolution. The key feature of the CM24 method is that it uses three integrator switching circuits and a cross-shaped integrator loop to update the electric field, and minimizes numerical dispersion by optimizing the coefficients in the combined loop.

[0034] Step S21: Based on the spatial grid size in the computational model, construct the electric field components containing fourth-order spatial differences. The time-step expression. The time-domain Maxwell curl equation in Cartesian coordinates is: ; in, Indicates along Electric field intensity components in the direction (unit: V / m). Indicates along Magnetic field strength components in the direction (unit: A / m). Indicates along Magnetic field strength components in the direction (unit: A / m). The dielectric constant of the medium (unit: F / m).

[0035] The time derivative is discretized using the central difference method, with the time step denoted as . The spatial grid step size is denoted as follows: and electric field components Defined in integer time steps ( ,in (Time step index), magnetic field components and Defined in half-integer time steps ( ).

[0036] For electric field components The spatial partial derivatives are approximated by the CM24 method using a weighted combination of central difference operators with different spans, which are compact difference operators spanning one grid step (span). p=1), and the extended difference operator spanning three grid steps (span) p =3). Along direction and y The spatial center difference of the direction is approximately: ; ; Thus, the electric field components are obtained. The CM24 time step expression: ; in, and These represent the grid points ( i , j )superior, time, right x direction and right y Directional difference operator. The difference coefficients (dimensionless) spanning one grid step. The difference coefficients (dimensionless) span three grid steps. The difference coefficients (dimensionless) span three grid steps. The difference coefficients (dimensionless) are for cross-shaped integral switching.

[0037] ; ; Step S22: Using the aforementioned spatial center difference operator, calculate the update equations for the electric and magnetic fields respectively. During the calculation, updating the electric field requires 24 magnetic field components from the surrounding 16 grids, and updating the magnetic field requires 24 electric field components from the surrounding 9 grids. Therefore, special boundary treatment is required at the boundary of the computational domain (near the inner side of the PML absorbing boundary).

[0038] Step S23: For the CM24 fourth-order spatial difference scheme, its numerical wavenumber Physical wavenumber in continuous space The relationship between them is determined by the following dispersion relation: ; in, Indicates along Numerical wavenumber in direction (unit: rad / m). Indicates along Actual physical wavenumber in direction (unit: rad / m). Indicates along Spatial grid step size in the direction (unit: m). , , Let be the coefficient, satisfying .along Numerical wavenumber of direction The expressions are similar. By selecting appropriate... , The value of makes the numerical wavenumber Approximating the physical wavenumber over a relatively large wavenumber range This reduces numerical dispersion error.

[0039] Step S24: Obtained in steps S21 to S23 Based on the component update expression and dispersion relation, a computational framework for low-numerical-dispersion wavefield propagation under coarse-grid resolution is obtained. This framework includes... The positional relationships of all field points involved in component updates.

[0040] Step S25: Based on the aforementioned low numerical dispersive wave field propagation calculation framework, and following the same derivation method, sequentially obtain the remaining five electromagnetic field components (in the three-dimensional case, they are...). , , , , The time-step formula is used to complete the construction of the complete CM24 electromagnetic field update equation, which includes all six electromagnetic field components.

[0041] In three dimensions, spatial coordinates are expanded to... The corresponding six field component update equations all adopt the same fourth-order spatial difference scheme to ensure consistent numerical accuracy throughout the entire three-dimensional computational space.

[0042] Through the above operations in this step, the CM24 electromagnetic field update equation with low numerical dispersion characteristics under coarse grid (1 / 4 wavelength resolution) conditions was successfully constructed, laying the core algorithm foundation for realizing efficient ground-penetrating radar wave field propagation simulation.

[0043] Step S3, Update the electromagnetic field equations after broadband correction: like Figure 3 Based on the spectral distribution characteristics of GPR and the numerical dispersion error curve of the CM24 algorithm, and based on the dispersion relationship of the CM24 computational framework and the spectral distribution characteristics of the ground penetrating radar operating frequency band, a wideband overall optimization strategy is used to select the coefficient combination that minimizes the numerical dispersion error throughout the entire operating frequency band. , This generates a broadband-corrected electromagnetic field update equation. For example... Figure 4 As shown, at a fixed grid resolution RWhen the coefficient is 4, wide-band dispersion error (WBDE) is used to measure the overall numerical dispersion error level corresponding to different coefficient combinations.

[0044] Step S31: Based on the CM24 calculation framework and the spectral distribution characteristics of ground-penetrating radar, obtain the broadband dispersion error measurement index at a fixed spatial resolution. (The spatial resolution is then set to a fixed value.) Under the condition of (i.e., 4 grids per minimum wavelength), define the direction of propagation ( (The direction of electromagnetic wave propagation is relative to) Angle between axes The angle relative to the z-axis direction The average numerical dispersion error on ) is:

[0045] For 3D mode, propagation direction Coverage 0 to Scope Coverage 0 to The range. Represents frequency Average numerical dispersion error (dimensionless) over the entire propagation space, integral sign Indicates the direction angle of propagation From 0 to Definite integral operations.

[0046] Step S32: Calculate the single frequency based on the above broadband dispersion error measurement index. The corresponding average numerical dispersion error across the entire propagation space In practical calculations, the direction angle is discretized. (For example, with) step size (Sampling is performed), and the above definite integral is replaced by numerical summation: ; in, This represents the total number of sampling points for the propagation direction angle. This indicates a summation operation.

[0047] Step S33: Based on the average numerical dispersion error Normalized amplitude with electric field spectrum Construct the objective function for the overall numerical dispersion error across the wideband: ; in, Indicating the combination of coefficients The overall numerical dispersion error index (dimensionless) for the wideband. Indicates the lower limit frequency of the operating frequency band (unit: Hz). Indicates the upper limit frequency of the operating frequency band (unit: Hz). The normalized amplitude (dimensionless) represents the frequency domain definite integral operation.

[0048] The excitation source is usually selected with a center frequency of The Ricker wavelet, or first-order Gaussian pulse derivative, has the following time-domain expression: ; in, express Excitation current density at time (unit: A / m) 2 ), Indicates the center frequency (unit: Hz). Indicates time (unit: seconds). Indicates the peak pulse time (unit: seconds, usually taken as...). ), Represented by natural constant An exponential function with base 0.05. The spectral amplitude of the Ricker wavelet. Obtained through Fourier transform, and then normalized. .

[0049] By searching in the above coefficient parameter space, Reaching the minimum value By combining these, the optimal differential coefficients after broadband correction can be obtained. In an exemplary preferred embodiment, for Spatial resolution, optimal coefficient combination can make the center frequency The overall numerical dispersion error across the entire operating frequency band on both sides is minimized.

[0050] The coefficients obtained from the above optimization Substituting these values ​​into the electromagnetic field update equation constructed in step S2 generates the broadband-corrected electromagnetic field update equation. This broadband coefficient correction solves the problem of insufficient broadband phase accuracy caused by traditional methods optimizing only a single frequency point, ensuring that the wave field propagation simulation under broadband pulse excitation by ground-penetrating radar maintains high phase accuracy throughout the entire operating frequency band.

[0051] Step S4, High-order convolutional perfectly matched layer (HO-CPML) absorbs boundary conditions: like Figure 5As shown, the absorption effect of CM24-HO-CPML is verified (this verification uses a simplified test scenario independent of the main simulation). Based on the broadband modified electromagnetic field update equation, by introducing a high-order scaling transformation factor, a high-order convolutional perfectly matched layer (HO-CPML) absorption boundary condition compatible with the fourth-order spatial difference scheme of CM24 is constructed, generating a truncated boundary absorption model with excellent absorption performance under coarse mesh conditions.

[0052] Step S41: Based on the broadband corrected electromagnetic field update equation, adjust the scaling transformation factor of the standard convolution perfectly matched layer. Perform factorization. In the frequency domain, the scaling coordinate transformation factor of CPML is expressed as: ; in, Indicates along direction( The scaling factor (dimensionless) of ). express The real stretching factor in the direction (dimensionless) ), express Conductivity distribution in the direction (unit: S / m). express Frequency shift factor in direction (unit: S / m) Represents the imaginary unit ( ), Angular frequency (unit: rad / s) Represents the vacuum permittivity ( F / m).

[0053] This embodiment uses the reciprocal of the standard CPML scaling factor. Transform it into the product of two terms:

[0054] Step S42: Calculate the poles and corresponding residues of the complex frequency function based on the product form of the two terms. It has two first-order poles in the complex frequency plane, located at: ; in, and They represent Two single poles in the complex frequency plane. The corresponding residues are respectively ( k =1,2): ; in, Indicates the first pole The corresponding residue, Indicates the second pole The corresponding residue.

[0055] By using the inverse Fourier transform, the frequency domain Transforming to the time domain, its time-domain impulse response is: ; in, Indicates along Temporal convolution active function in the direction (unit: s) -1 ), Dirac function, Represents the unit step function. Indicates time (unit: seconds).

[0056] Step S43: Based on the two first-order poles and their corresponding residues, construct the recursive relationship for the auxiliary terms. In the numerical implementation of the finite-difference time domain, the electric field components... Convolution operation in free space Updated using the following recursive relation:

[0057] ; in, , Indicates the first time step edge , Auxiliary convolution variables in the direction, This represents the recursive decay coefficient (dimensionless). Represents the recursive coefficient (the unit depends on the specific implementation). , Indicates the first The spatial partial derivative of the magnetic field at each time step.

[0058] For the higher-order CPML in this embodiment, the recursion coefficients and The expression is as follows: ; in, Indicates the time step (unit: s). Indicates along Spatial grid step size in the direction (unit: m).

[0059] By absorbing the boundary layer , and Extending the polynomial order to , and A high-order polynomial distribution (e.g., second- or third-order polynomial distribution) is adopted along the depth direction of the absorbing boundary, forming a high-order CPML (HO-CPML) absorbing boundary condition: ; in, Represents the depth coordinates (in meters) measured from the boundary of the PML region. This indicates the total thickness of the PML layer (unit: m). express The maximum value at the outer boundary of PML. express The maximum value at the outer boundary of PML. , They represent , The order of the polynomial distribution.

[0060] Integrating the aforementioned HO-CPML convolution terms into the broadband-corrected CM24 electromagnetic field update equation generated in step S3 completes the construction of the truncated boundary absorption model under coarse-grid conditions. Through the high-order CPML processing in this step, the absorption boundary effectively enhances its absorption capacity for evanescent and low-frequency wave components under coarse-grid resolution, significantly suppresses reflection interference at the truncated boundary, and avoids additional reflection problems at the interface between coarse and fine grids.

[0061] Step S5, forward wave field propagation calculation: Using the truncated boundary absorption model constructed in step S4 and the broadband corrected electromagnetic field update equation generated in step S3, forward wave field propagation calculations are performed to obtain forward simulation data for each receiving point location.

[0062] Step S51: Based on the truncated boundary absorbing model, set multiple layers of absorbing boundaries around the computational region. In this embodiment, the computational region along... direction, direction and Set up around the direction Layer HO-CPML absorbing boundary (e.g.) (8 to 12 layers), the thickness of the absorption boundary is (in (Spatial grid step size). Medium parameters within the absorbing boundary. , Set the higher-order polynomial distribution according to step S43.

[0063] Step S52: Based on the updated electromagnetic field equations after the multilayer absorption boundary and broadband correction, set up a current source with a specific polarization direction and center frequency. As the excitation source. The excitation source uses a center frequency of The Ricker subwavelength (e.g., 500MHz) is expressed in time domain as described in step S33. The spatial location of the excitation source is applied at the grid point where the transmitting antenna Tx is located, with the polarization direction along... Axial direction. The excitation source is loaded in the form of electric field components at each time step. Add a current source term to the updated equation.

[0064] Step S53: Perform time-domain finite-difference time-step calculation based on the excitation source. (The calculation is performed in time steps.) Increasing from 1 to (Total number of time steps), at each time step, update all six electromagnetic field components sequentially. , , , , , This completes the calculation of the forward wave field propagation from the start of excitation to the propagation of the electromagnetic wave to all spatial points within the computational region. Total simulation time... The value of needs to ensure that the electromagnetic wave can propagate from the transmitting antenna to the target at the deepest point in the computational region and return to the receiving antenna. During the wave field propagation calculation, the auxiliary convolution variable within the PML absorbing layer... Synchronized updates.

[0065] Step S54: Record the forward simulation data of the entire computational domain and the electric field components at each receiving point Rx location. Time series data were obtained to acquire forward modeling data at coarse grid resolution. For multiple receiving antennas deployed on the ground (e.g., uniformly distributed along the survey line)... (1 receiving point), extract the location of each receiving point at all time steps from the simulation results. Field values ​​form the received data matrix. ,in Indicates the receiving time index. Indicates the receiver point index.

[0066] This step yields complete and effectively suppressed boundary reflection interference forward modeling data at coarse grid resolution, providing high-quality input data for subsequent reverse time migration imaging.

[0067] Step S6, Calculation of reverse wave field propagation: The recorded echo data in the forward simulation data obtained in step S5 is time-reversed and used as an excitation source to be reloaded to the corresponding receiving point position. The reverse wave field propagation calculation is then performed to obtain the reverse propagation wave field.

[0068] Step S61: Based on the forward modeling data, extract the echo data recorded at each receiving point (i.e., the electric field component containing information on reflection and scattering from the subsurface medium). Since the HO-CPML absorption boundary in step S4 has significantly suppressed boundary reflection, the recorded time series of each receiving point can be directly used as echo data.

[0069] Step S62: Reverse the echo data according to time. The specific operation is as follows: If the receiving point... The time series recorded during forward propagation is ,in Then its time-reversed sequence is The time-reversed echo data is used as a new excitation source, maintaining the same physical parameter settings as the forward propagation, including the same medium electromagnetic parameter distribution, the same spatial grid size, and the same time step.

[0070] Step S63: Reload the excitation source obtained after time reversal to the corresponding receiving point position. Specifically, reload the excitation source obtained after time reversal to the corresponding receiving point position. Time-reversed data at each receiving point location As the excitation source at this location, the loading method is the same as that of the excitation source of the transmitting antenna in forward propagation (in the electric field component). (A current source term is added to the update equation), but at this time the spatial location of the excitation source is distributed at all receiving points rather than at a single transmitting antenna.

[0071] Step S64: Using the same broadband modified electromagnetic field update equation and HO-CPML truncation boundary absorption model as in Step S5, and with the time-reversed echo data as the excitation source, perform reverse wave field propagation calculations from the last time step to the initial time step (i.e., time index from...). Decrease to 1) to obtain the complete backpropagation wavefield. The backpropagation wavefield is recorded as the electric field components at each time step of each spatial grid point. .

[0072] Through this step, the echo information recorded at the receiving point is successfully re-injected into the computational space through time reversal and backpropagation. This allows the backpropagation wavefield to focus naturally when it encounters discontinuous interfaces in the underground medium and embedded targets, providing backpropagation wavefield data for cross-correlation calculations of reverse time migration imaging.

[0073] Step S7, Reverse Time Migration (RTM) Imaging Results: The forward propagation wavefield obtained in step S5 and the reverse propagation wavefield obtained in step S6 are subjected to zero-delay cross-correlation imaging processing to generate ground-penetrating radar reverse time migration (RTM) imaging results that reflect the discontinuous interface of the underground medium and the target location.

[0074] Step S71: Traverse every spatial grid point within the computational region. For the three-dimensional case, traverse all... Grid index.

[0075] Step S72: For each spatial grid point Perform forward propagation of wavefield values ​​at all time steps. With the backpropagation wave field value Time-step product operation. Forward propagating wave field. The electric field component values ​​at each time step of each grid point are taken from the forward modeling process in step S5, and the backpropagation wave field is... The corresponding data is taken from the backpropagation simulation process in step S6.

[0076] Step S73: Sum the product results of all time steps to obtain the spatial grid point. The zero-delay cross-correlation imaging value at that location is given by the imaging condition formula: ; in, Represents spatial grid points The reverse time-shifted imaging value at the location (dimensionless or with energy dimension). Indicates the index of time steps From 0 to Summation operation, Indicates the first A time step at a spatial grid point The forward propagation electric field component at the location (unit: V / m). Indicates the first A time step at a spatial grid point The value of the back-propagating electric field component at the location (unit: V / m). The above zero-delay cross-correlation imaging condition utilizes the coherence of the forward and reverse wave fields at the same time. At the interface of the underground medium discontinuity and the location of the target, the forward-propagating incident wave field and the reverse-propagating reflected wave field are synchronized in time, and the cross-correlation result produces a high-amplitude response.

[0077] Step S74: Calculate the zero-delay cross-correlation results of all spatial grid points. This generates the final ground-penetrating radar reverse-time migration imaging profile. The imaging profile is presented in spatial coordinates. The horizontal and vertical axes (for a two-dimensional section), or with In a three-dimensional coordinate system (for a three-dimensional data volume), the imaging values The size is rendered in grayscale or color to reflect the spatial location and morphology of discontinuous interfaces (such as layer interfaces and fracture surfaces) and abnormal targets (such as cavities, pipes, cracks, etc.) in the underground medium.

[0078] In a preferred embodiment, the imaging results Subsequent processing is performed, including but not limited to: normalization (mapping the imaging values ​​to the [0,1] interval), Laplacian filtering (enhancing interface continuity), gain compensation (compensating for deep energy attenuation), etc., to improve the interpretability of the imaging profile.

[0079] like Figure 7 As shown, Figure 7 (a) The VV polarization RTM imaging results clearly depict the layered interfaces between region B and region C, and between region C and region D. Meanwhile, within region B, three and two rock targets are clearly visible at the two red circle markers, respectively. Figure 7 (b) shows that the VH polarization RTM imaging results exhibit strong scattered energy ripples within the yellow circle area, which can provide a valid reference for the identification of the rock piles in this area using VV polarization imaging results. This imaging result can clearly reflect the discontinuities of the underground medium, the location and morphology of buried targets and anomalies.

[0080] like Figure 8 As shown, Figure 8 (a) shows the imaging results of the method (CM24-HO-CPML) in this embodiment. Figure 8 (b) Image results using CM24 combined with traditional CPML. As can be seen in the comparison, when using traditional CPML, obvious boundary reflection fringes appear in the image, imaging artifacts exist above the pile of stones in the region, and the interface between region C and region D is relatively blurred. This is mainly due to the lower grid resolution. R At a resolution of 4, traditional CPML absorption at the boundary exhibits significant reflection errors, thus affecting the quality of RTM imaging. The method in this embodiment effectively suppresses boundary reflections using CM24-HO-CPML, resulting in clearer and more accurate imaging results.

[0081] Through zero-delay cross-correlation imaging processing in this step, the information from the forward and reverse wave fields is successfully fused into a high-resolution imaging result that reflects the subsurface structure, accurately characterizing the discontinuous interfaces and target volume distribution of the subsurface medium without relying on a priori velocity models.

[0082] Embodiments of the present invention also provide a high-efficiency computational system for ground-penetrating imaging based on CM24 high-order finite-difference time-domain and high-order PML, including: The model building module is used to set the calculation area, spatial grid size, medium electromagnetic parameters and center frequency according to the detection scenario, and to determine the position of the transmitting and receiving antennas to build a ground penetrating radar calculation model. The equation construction module is used to construct an electromagnetic field update equation containing fourth-order spatial difference based on the computational model, and obtain a computational framework for low numerical dispersive wave field propagation. The coefficient optimization module is used to select the coefficients that minimize the overall numerical dispersion error within the working frequency band based on the computational framework and the characteristics of the working spectrum distribution, and to generate the broadband corrected electromagnetic field update equation. The boundary construction module is used to construct the high-order convolution fully matched layer absorption boundary conditions based on the broadband modified electromagnetic field update equation and through a high-order scaling transformation factor, thereby generating a truncated boundary absorption model. The wavefield calculation and imaging module is used to calculate the forward wavefield propagation based on the truncated boundary absorption model to obtain forward modeling data. The data is then reversed in time and used as an excitation source to reload the data to the receiving point to calculate the reverse wavefield propagation to obtain the reverse propagation wavefield. Finally, the forward and reverse propagation wavefields are multiplied at each time step and the product of all time steps is accumulated to generate the ground penetrating radar reverse time migration imaging result.

[0083] It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

[0084] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A high-efficiency calculation method for ground-penetrating imaging based on CM24 high-order finite-difference time-domain and high-order PML, characterized in that, The method includes: Based on the detection scenario, the calculation area, spatial grid size, dielectric parameters of the medium, excitation source type and center frequency are set, and the positions of the transmitting and receiving antennas are determined to establish a calculation model of the ground penetrating radar detection area. Based on the computational model, an electromagnetic field update equation containing fourth-order spatial difference is constructed to obtain a computational framework for low-numerical dispersive wave field propagation under coarse-grid resolution conditions. Based on the aforementioned computational framework and the spectral distribution characteristics of the ground-penetrating radar, coefficients that minimize the overall numerical dispersion error within the operating frequency band are selected to generate the broadband corrected electromagnetic field update equation. Based on the broadband corrected electromagnetic field update equation, by introducing a higher-order scaling transformation factor, a higher-order convolution fully matched layer absorption boundary condition is constructed to generate a truncated boundary absorption model under coarse grid conditions. Based on the truncated boundary absorption model, the forward wave field propagation calculation is performed on the broadband corrected electromagnetic field update equation to obtain forward modeling data for each receiving point. Based on the forward modeling data, the recorded echo data is reversed in time and reloaded as an excitation source to the corresponding receiving point position to perform reverse wave field propagation calculation and obtain the reverse propagation wave field. Based on the forward propagation wavefield and the reverse propagation wavefield, at each time step, the forward wavefield value is multiplied by the reverse wavefield value after time inversion, and the product results of all time steps are accumulated to generate a ground-penetrating radar reverse time migration imaging result that reflects the discontinuous interface of the underground medium and the target position.

2. The method according to claim 1, characterized in that, Based on the aforementioned computational model, an electromagnetic field update equation incorporating fourth-order spatial difference is constructed to obtain a computational framework for low-numerical-dispersion wavefield propagation under coarse-grid resolution conditions, specifically including: Based on the spatial grid size in the computational model, a time-step expression for the electric and magnetic field components containing fourth-order spatial differences is constructed. Based on the time step expression of the electric field component, a three-dimensional non-dissipative wavenumber is achieved by adding a cross-shaped integral loop on a two-dimensional template using the spatial center difference operator. Based on the three-dimensional non-dissipative wavenumber, the phase accuracy is improved while maintaining the fourth-order spatial difference, and a low-numerical dispersive wave field propagation calculation framework under coarse grid resolution is obtained. Based on the low numerical dispersive wave field propagation calculation framework, the time step formulas for all electromagnetic field components are derived and obtained, thus completing the construction of the complete electromagnetic field update equations.

3. The method according to claim 2, characterized in that, Based on the computational framework and the spectral distribution characteristics of the ground-penetrating radar, coefficients that minimize the overall numerical dispersion error within the operating frequency band are selected to generate the broadband corrected electromagnetic field update equation, specifically including: Based on the computational framework and the spectral distribution characteristics of the ground penetrating radar, a broadband dispersion error measurement index under a fixed spatial resolution is obtained. Based on the broadband dispersion error measurement index, calculate the average numerical dispersion error of a single frequency throughout the entire propagation space; Based on the magnitude of the average numerical dispersion error and the normalized amplitude of the electric field spectrum, coefficients that minimize the overall numerical dispersion error within the working frequency band are selected to generate the broadband corrected electromagnetic field update equation.

4. The method according to claim 3, characterized in that, Based on the broadband corrected electromagnetic field update equation, by introducing a higher-order scaling transformation factor, a higher-order convolutional perfectly matched layer absorbing boundary condition is constructed to generate a truncated boundary absorbing model under coarse-grid conditions, specifically including: Based on the broadband corrected electromagnetic field update equation, the scaling transformation factor of the standard convolution fully matched layer is expanded into a product of two terms; Based on the product of the two terms, obtain the two first-order poles and the corresponding residues of the poles for the corresponding scaling transformation factor; Based on the two first-order poles and their corresponding residues, an auxiliary term is constructed to realize the recursive convolution calculation of the absorbing boundary, making it compatible with the fourth-order spatial difference scheme, and generating a truncated boundary absorbing model under coarse grid conditions.

5. The method according to claim 4, characterized in that, Based on the truncated boundary absorption model, the forward wave field propagation calculation is performed on the broadband corrected electromagnetic field update equation to obtain forward modeling data for each receiving point, specifically including: According to the truncated boundary absorption model, multiple layers of absorption boundaries are set around the calculation area to ensure that electromagnetic waves are effectively absorbed when they propagate to the absorption boundary during the simulation time. Based on the multilayer absorption boundary and the broadband modified electromagnetic field update equation, a current source with polarization direction and center frequency is set as the excitation source. Based on the excitation source, perform forward wave field propagation calculations to obtain forward modeling data of the entire calculation area and recorded data of each receiving point; Based on the recorded data, the electric field components at the observation point locations are extracted to obtain forward modeling data at a coarse grid resolution.

6. The method according to claim 5, characterized in that, Based on the forward modeling data, the recorded echo data is reversed in time and reloaded as the excitation source to the corresponding receiving point position. Reverse wavefield propagation calculations are then performed to obtain the reverse propagation wavefield, specifically including: Based on the forward modeling data, extract the echo data recorded at each receiving point; Based on the echo data, the echo data is reversed in time and used as a new excitation source, with the medium type set to vacuum type; Based on the new excitation source, it is reloaded to the corresponding receiving point position, and the reverse wave field propagation calculation is performed using the same broadband modified electromagnetic field update equation and truncated boundary absorption model to obtain the reverse propagation wave field.

7. The method according to claim 6, characterized in that, Based on the forward propagation wavefield and the reverse propagation wavefield, at each time step, the forward wavefield value is multiplied by the time-inverted reverse wavefield value, and the product results of all time steps are accumulated to generate a ground-penetrating radar reverse-time migration imaging result reflecting the discontinuity interface of the underground medium and the target location, specifically including: Based on the forward propagation wave field and the reverse propagation wave field, traverse every spatial grid point within the calculation area; Based on each spatial grid point, at each time step, the forward wave field value is multiplied by the reverse wave field value after time inversion; Based on the multiplication result, the product results of all time steps are summed to obtain the zero-delay cross-correlation calculation result of the spatial grid point; Based on the zero-delay cross-correlation calculation results, ground-penetrating radar reverse time migration imaging results reflecting the discontinuous interface of the underground medium and the target location are generated.

8. The method according to claim 7, characterized in that, Based on the detection scenario, the computational region, spatial grid size, dielectric parameters, excitation source type, and center frequency are set, and the positions of the transmitting and receiving antennas are determined. A computational model of the ground-penetrating radar detection area is then established, specifically including: Based on the detection scenario, set along x direction, y direction and z The dimensions of the three-dimensional computational region in the direction are determined, and multiple absorbing boundaries are set around the computational region; Based on the dimensions of the three-dimensional calculation region, the dielectric parameters of different dielectric layers are set, and the position and shape of the embedded target are determined. Based on the electromagnetic parameters and the maximum frequency of the transmitted electric field, the spatial grid size corresponding to 1 / 4 wavelength grid resolution is set, and the positions of the transmitting and receiving antennas and the measurement point intervals in the bistatic radar observation mode are determined to establish a calculation model of the ground penetrating radar detection area.

9. The method according to claim 8, characterized in that, The process involves constructing an electromagnetic field update equation with fourth-order spatial difference based on the computational model, obtaining a computational framework for low-numerical dispersive wave field propagation under coarse-grid resolution conditions, and further includes specific update operations for the electromagnetic field components, specifically including: Based on the aforementioned computational framework, the update expressions for the electromagnetic field components at specific moments on the grid points are obtained; Based on the update expression and the coefficients to be optimized, the time-step update of the electromagnetic field components is completed, and a calculation framework for low numerical dispersive wave field propagation under coarse grid resolution is obtained.

10. A high-efficiency computing system for ground-penetrating imaging based on CM24 high-order finite-difference time-domain and high-order PML, characterized in that, The system is used to perform the method as described in any one of claims 1 to 9, comprising: The model building module is used to set the calculation area, spatial grid size, dielectric parameters and center frequency according to the detection scenario, and to determine the position of the transmitting and receiving antennas to build a ground penetrating radar calculation model. The equation construction module is used to construct an electromagnetic field update equation containing fourth-order spatial difference based on the computational model, and obtain a computational framework for low numerical dispersive wave field propagation. The coefficient optimization module is used to select the coefficients that minimize the overall numerical dispersion error within the working frequency band based on the computational framework and the characteristics of the working spectrum distribution, and to generate the broadband corrected electromagnetic field update equation. The boundary construction module is used to construct the high-order convolution fully matched layer absorption boundary conditions based on the broadband optimized electromagnetic field update equation and through a high-order scaling transformation factor, thereby generating a truncated boundary absorption model. The wavefield calculation and imaging module is used to calculate the forward wavefield propagation based on the truncated boundary absorption model to obtain forward modeling data. The data is then reversed in time and used as an excitation source to reload the data to the receiving point to calculate the reverse wavefield propagation to obtain the reverse propagation wavefield. Finally, the forward and reverse propagation wavefields are multiplied at each time step and the product of all time steps is accumulated to generate the ground penetrating radar reverse time migration imaging result.

Images