A Numerical Solution Method and System for Ground Penetrating Radar Based on Divergence-Preserving ADI-FDTD

Abstract

The present invention discloses a numerical solution method and system for ground penetrating radar based on divergence-preserving ADI-FDTD, belonging to the field of efficient numerical simulation of electromagnetic waves in ground penetrating radar. The method includes: listing the differential time domain of the Maxwell's equations of the ground penetrating radar; discretizing the differential time domain into a matrix form; for a multi-pole Debye dispersive medium with a loss term, constructing a constitutive relation between the electric flux density and the electric field strength; transforming the constitutive relation from the frequency domain to the z domain, and then through the transformation between the z domain and the time domain, converting it to the time domain, and substituting it into the matrix form discretized from the differential time domain to obtain a discrete matrix form based on the time domain, and discretizing it into two sub-time step iterations to obtain a basic numerical discretization framework; calculating the numerical discretization equations of each component of the electric field and the magnetic field according to the basic numerical discretization framework, and calculating the electric field strength and the magnetic field strength, and deriving the electric field strength and the magnetic field strength that satisfy divergence preservation.

Description

Technical Field

[0001] The present invention relates to the field of efficient numerical simulation of electromagnetic waves in ground penetrating radar, and specifically relates to a numerical solution method and system for ground penetrating radar based on divergence-preserving ADI-FDTD. Background Art

[0002] Ground Penetrating Radar (GPR) is a widely used geophysical technique that uses antennas to transmit and receive electromagnetic (EM) waves and can detect the material properties and distribution patterns in underground media. When the EM wave passes through the ground, its path, field strength, and waveform are affected by the electrical and geometric properties of the medium. By analyzing the propagation time, amplitude, waveform, and other parameters of the received signal (usually displayed on an oscilloscope) and comparing them with a pre-calculated model, the spatial position, structure, and distribution of underground interfaces or geological structures can be inferred. Generating a database of reflected wave signals corresponding to underground conditions through EM simulation facilitates forward modeling. The Maxwell equations govern the propagation of the EM field through these media, and solving these equations is the key to accurately simulating ground penetrating radar.

[0003] The FDTD method is a widely adopted method for discretizing the Maxwell equations in the time and space domains. By doing so, it enables the EM field to be computed along the time axis in the computational space, effectively capturing transient EM field information. This method incorporates the electromagnetic properties of the target medium into the finite-difference time-domain method using a spatial grid. By assigning appropriate electromagnetic parameters in different regions, various complex structures including inhomogeneity, anisotropy, dispersion characteristics, and nonlinearity can be simulated.

[0004] Despite its practicality, the traditional explicit FDTD method has a key limitation: the time step is restricted by the Courant-Friedrichs-Lewy (CFL) stability condition. In GPR simulation, the dispersion characteristics of the soil and the target, combined with the irregular structural features in certain regions, require the use of a fine, uniform Yee grid to obtain numerical accuracy. However, this leads to two major challenges:

[0005] 1) The time step must be small enough to satisfy the CFL stability condition, resulting in a long computation time;

[0006] 2) Uniform fine grid discretization increases the number of unknowns and significantly increases the memory usage.

[0007] To address these challenges, several unconditionally stable FDTD methods, such as the ADI-FDTD and LOD-FDTD methods, and the DP-ADI-FDTD method, which is a divergence-preserving and unconditionally stable method that loses its divergence-preserving property when solving Maxwell's equations, leading to charge accumulation in the source-free region, have been introduced to solve the EM Particle in Cell (EMPIC) problem. This method was later extended to simulate multi-pole Debye dispersive media. Despite its advantages, the numerical implementation based on the ADE method remains complex and cumbersome.

[0008] In addition, sub-grid techniques have been proposed to improve the computational efficiency of GPR and other EM simulations by eliminating the need for a uniform Yee grid throughout the computational domain. These methods allow for finer discretization in regions where high resolution is required and coarser grids in less critical regions, thus preventing unnecessary oversampling. Previous sub-grid methods can be classified into two types:

[0009] 1) One method uses a time step that satisfies the stability condition of a fine grid over the entire domain, where spatial interpolation is used to exchange EM field information between the coarse grid region and the fine grid region;

[0010] 2) Another method uses different time steps for the coarse grid region and the fine grid region, each time step satisfying its respective CFL stability condition, and employs spatio-temporal interpolation to exchange EM information between the grids.

[0011] Recently, unconditionally stable FDTD methods have been applied to the fine grid region, enabling the entire domain to be solved using a time step that satisfies the CFL stability condition of the coarse grid.

[0012] These sub-grid techniques allow for larger time steps, which significantly improve the computational efficiency compared to traditional hybrid sub-grid methods. Moreover, using a unified time step eliminates the need for time interpolation and simplifies the numerical implementation. However, existing methods still face some challenges, including:

[0013] 1) Mainly focused on two-dimensional models, with limited applicability to complex dispersive media;

[0014] 2) Lack of a unified electromagnetic exchange formula for arbitrary ratios at the coarse-fine grid interface;

[0015] 3) Lack of the divergence-preserving property of the hybrid sub-grid method used in unconditionally stable FDTD methods. Summary of the Invention

[0016] In view of the deficiencies of the prior art, the present invention proposes a numerical solution method and system for ground penetrating radar based on divergence-preserving ADI-FDTD.

[0017] The object of the present invention can be achieved by the following technical solutions:

[0018] In the first aspect of the present invention, a numerical solution method for ground penetrating radar based on divergence-preserving ADI-FDTD is involved, including the following steps:

[0019] List the differential time domain of the Maxwell equations of the ground penetrating radar;

[0020] Discretize the differential time domain of the Maxwell equations into matrix form;

[0021] For a multi-pole Debye dispersive medium with loss terms, construct a constitutive relation between the electric flux density and the electric field strength;

[0022] According to the z-transform theory, transform the constitutive relation from the frequency domain to the z-domain, and then through the transformation between the z-domain and the time domain, transform it to the time domain, and substitute it into the matrix form discretized from the differential time domain of the Maxwell equations to obtain a discrete matrix form based on the time domain;

[0023] Discretize the discrete matrix form based on the time domain into two sub-time step iterations, the first sub-time period , the second sub-time period , derive the numerical iteration formula for the electromagnetic field components, where n is the number of iterations; calculate the electric and magnetic field strength expressions for the first sub-time period and the second sub-time period respectively, and obtain the basic numerical discretization framework of the ZT-DP-ADI-FDTD method for the lossy multi-pole Debye model;

[0024] According to the basic numerical discretization framework, calculate the numerical discretization equations for each component of the electric and magnetic fields, and calculate the electric field strength and the magnetic field strength, and derive the electric field strength and the magnetic field strength that satisfy divergence preservation.

[0025] Optionally, the differential time domain expression of the Maxwell equations is as follows:

[0026] (1)

[0027] (2)

[0028] Where D is the electric flux density, H is the magnetic field strength, E is the electric field strength, μ 0 is the magnetic permeability of free space, and the constitutive relation between the electric flux density D and the electric field E is expressed as D = ε 0 ε r E , ε 0 is the permittivity of free space, ε r is the relative permittivity.

[0029] Optionally, the matrix form is:

[0030] (3)

[0031] where V DP = E DP_x , E DP_y , E DP_z , H DP_x, H DP_y , H DP_z represents the divergence-free electric and magnetic field strength components, V DP The superscripts n and n + 1 of E DP_x , E DP_y , E DP_z respectively represent the electric field strengths in the x, y, and z directions. H DP_x, H DP_y , H DP_z respectively represent the magnetic field strengths in the x, y, and z directions. The matrices A and B are spatial partial derivative matrices, specifically:

[0032] , , where , , , and respectively represent the partial derivatives of the field components in the x, y, and z directions.

[0033] Optionally, the constitutive relation between the electric flux density and the electric field strength is constructed as:

[0034] (4)

[0035] where is the conductivity, is the difference between the relative permittivities at zero frequency and infinite frequency, represents the relative permittivity at infinite frequency; is the pole relaxation time, is the number of poles in the sensitivity response; represents the imaginary part; represents the angular frequency; represents the electric flux density in the frequency domain; represents the electric field strength in the frequency domain.

[0036] Optionally, the transformation of the constitutive relation from the frequency domain to the z-domain and then to the time domain through the transformation between the z-domain and the time domain includes the following steps:

[0037] According to the z-transform theory, transform (4) from the frequency domain to the z-domain to obtain:

[0038] (5)

[0039] where ψ p and Q are auxiliary variables, defined as follows:

[0040] (6)

[0041] (7)

[0042] Through the transformation between the z-domain and the time domain, transform formulas (5), (6), and (7) into the time domain and substitute them into (3) to obtain:

[0043] (8)

[0044] where φ p = ψ p / ε 0 ε ∞ , = Q / ε 0 ε ∞ , and in addition, the ε r in A and B in the matrix becomes ε ∞， is the difference between the relative dielectric constants at zero frequency and infinite frequency; : ψ in the z-domain p , : Q in the z-domain.

[0045] Optionally, the calculation process of the basic numerical discretization framework of the lossy multi-pole Debye model ZT-DP-ADI-FDTD method includes the following steps:

[0046] Discretize formula (8) into two sub-time step iterations. For the first sub-time period :

[0047] (9)

[0048] (10)

[0049] For the second sub - time period :

[0050] (11)

[0051] (12)

[0052] where V = E x , E y , E z , H x , H y , H z , E x 、 E y 、 E z respectively represent the components of the electric field that does not obey the divergence property in the x, y, and z directions, H x 、 H y 、 H z respectively represent the components of the magnetic field that does not obey the divergence property in the x, y, and z directions; 、 φ p 、 V DP 、 The superscripts all represent the number of iteration steps;

[0053] Substituting (10) into (11) gives:

[0054] (13)

[0055] Define the auxiliary variable U = U E U H = U ex , U ey , U ez , Uhx , U hy , U hz , U E 、 U H respectively represent the auxiliary variables of the electric field and the magnetic field, U ex 、 U ey 、 U ez respectively represent the components of the auxiliary variable of the said electric field in the x, y, and z directions, U hx 、 U hy 、 U hz respectively represent the components of the auxiliary variable of the said magnetic field in the x, y, and z directions; where ; then (13) is reformulated as:

[0056] (14)

[0057] Substitute the previous time step of (9) into formula (12) to obtain:

[0058] (15)

[0059] Let U n+1 = V n+1 + V n+1 / 2 , then formula (15) is reformulated as:

[0060] (16).

[0061] Optionally, according to formula (14), obtain U E and U H the numerical discretization equations of the field components:

[0062] (17)

[0063] (18)

[0064] ψ p 、 U E 、 U HThe superscripts all represent the iteration steps, where:

[0065] (19)

[0066] (20)

[0067] Numerical iteration formulas for the electric and magnetic fields are obtained:

[0068] (24a)

[0069] In addition, and The iteration formulas are:

[0070] (24b)

[0071] (24c)

[0072] U ex ,、 U ey 、 U ez 、 E x 、 E y 、 E z The superscripts all represent the iteration steps, 、 、 represent the coordinate indices in the x, y, and z directions respectively;

[0073] Substitute the condition into formula (18) to get:

[0074] (25)

[0075] Obtain The calculation formula is:

[0076] (26)

[0077] Substitute formulas (24a), (24b), (24c), and (26) into formula (9) to obtain the calculation formulas for the E DP and H DP fields that satisfy divergence conservation:

[0078] (27)

[0079] H represents the magnetic field, E represents the electric field, and the superscripts of E and H represent the iteration steps.

[0080] The second aspect of the present invention relates to a numerical solution system for ground penetrating radar based on divergence-preserving ADI-FDTD, including:

[0081] A Maxwell's equations calculation module, configured to list the differential time domain of Maxwell's equations for ground penetrating radar; discretize the differential time domain of the Maxwell's equations into a matrix form;

[0082] A constitutive relation formula construction module for electric flux density and electric field strength, which constructs a constitutive relation formula for electric flux density and electric field strength for a multi-pole Debye dispersion medium with a loss term;

[0083] A z-domain transformation module: configured to, according to the z-transform theory, transform the constitutive relation formula from the frequency domain to the z-domain, and then through the transformation between the z-domain and the time domain, transform it into the time domain, and substitute it into the matrix form discretized from the differential time domain of the Maxwell's equations to obtain a discrete matrix form based on the time domain;

[0084] A basic numerical discretization framework construction module: discretize the discrete matrix form based on the time domain into two sub-time step iterations, the first sub-time period , the second sub-time period , derive the numerical iteration formula for the electromagnetic field components, where n is the number of iterations; calculate the electric field and magnetic field strength expressions for the first sub-time period and the second sub-time period respectively, and obtain the basic numerical discretization framework of the ZT-DP-ADI-FDTD method for the lossy multi-pole Debye model;

[0085] And, a divergence-preserving electric field strength and magnetic field strength solution module: calculate the numerical discretization equations for each component of the electric field and magnetic field according to the basic numerical discretization framework, and calculate the electric field strength and magnetic field strength, and derive the divergence-preserving electric field strength and magnetic field strength.

[0086] The third aspect of the present invention relates to a numerical solution method for ground penetrating radar based on a three-dimensional hybrid sub-grid. The three-dimensional hybrid sub-grid includes a number of coarse grids and fine grids. The excitation source is placed in the coarse grid area. When the electromagnetic wave propagates to the interface area between the coarse grid and the fine grid, interpolate the electromagnetic field values in the coarse grid at the interface into the corresponding fine grid components;

[0087] The method includes the following steps:

[0088] Update the electric field values and magnetic field values of the entire coarse grid using the finite-difference time-domain method;

[0089] Calculate the electric field values and magnetic field values of the fine grid in the fine grid area through the above-mentioned numerical solution method for ground penetrating radar based on divergence-preserving ADI-FDTD;

[0090] Interpolate the electric field and magnetic field values in the coarse grid to the interface of the fine grid at the same interface by using spatial linear interpolation; perform numerical iteration by using the above-mentioned ground penetrating radar numerical solution method based on divergence-preserving ADI-FDTD, and transfer the electric field component and magnetic field component at the interface of the fine grid to the fine grid area;

[0091] After the iteration is completed in the fine grid area, interpolate the values of the fine grid electromagnetic field at the junction and adjacent junctions of the coarse grid and the fine grid back to the electromagnetic field values of the corresponding coarse grid at the junction.

[0092] The fourth aspect of the present invention relates to a ground penetrating radar, including a storage medium storing instructions capable of running the above-mentioned ground penetrating radar numerical solution method based on divergence-preserving ADI-FDTD or the above-mentioned ground penetrating radar numerical solution method based on three-dimensional hybrid sub-grids, or the above-mentioned ground penetrating radar numerical solution system based on divergence-preserving ADI-FDTD.

[0093] Advantages of the present invention:

[0094] 1) The DP-ADI-FDTD method (ZT-DP-ADI-FDTD) based on Z-transform of the present invention can effectively and accurately simulate the electromagnetic wave propagation in lossy multi-level Debye dispersive media due to its unconditional stability and divergence-preserving characteristics. Compared with the ADE method, the Z-transform technology converts the multi-level Debye model from the frequency domain to the Z domain, eliminating the need to store the fields from multiple previous time steps while retaining the core iteration equation. This results in a more flexible and efficient numerical implementation.

[0095] 2) The present invention combines the explicit FDTD method with the implicit ZT-DP-ADI-FDTD method to develop a hybrid sub-grid technology. This method uses the ZT-DP-ADI-FDTD method to calculate the electromagnetic field components in the fine grid area and the FDTD method in the coarse grid area. This hybrid sub-grid technology can adapt to any odd-even grid ratio, further improving the flexibility of the modeling process. Brief description of the drawings

[0096] The present invention will be further described below with reference to the accompanying drawings.

[0097] Figure 1 For Embodiment 2 of the present application, the divergence of the electric flux density calculated by different methods (unit: dB): (a) FDTD method, (b) DP-ADI-FDTD method (CFLN = 1), (c) DP-ADI-FDTD method (CFLN = 3), (d) DP-ADI-FDTD method (CFLN = 5);

[0098] Figure 2For the ZT-DP-ADI-FDTD method in Embodiment 3 of the present application, when the time step is taken as 1, 3, and 5 times the time step of the traditional FDTD at the stability limit (i.e., CFLN = 1, 3, 5), the modulus of the eigenvalue of the amplification matrix M is obtained;

[0099] Figure 3 For the three-dimensional hybrid sub-grid schematic diagram in Embodiment 4 of the present application, where (a) is the three-dimensional interface between the coarse grid and the fine grid, and (b) is the interpolation schematic diagram of E x on the interface.

[0100] Figure 4 For the three-dimensional ground penetrating radar system for multi-target detection in Embodiment 5 of the present application, the transmitting end is Tx and the receiving end is Rx;

[0101] Figure 5 For the time-domain waveform of the receiving point at different Tx-Rx antenna distances in Embodiment 5 of the present application;

[0102] Figure 6 For the time snapshot of calculating the Ez amplitude using the hybrid sub-grid DP-ADI-FDTD method in Embodiment 5 of the present application;

[0103] Figure 7 For the time-domain result of Ez at the receiving antenna in Embodiment 5 of the present application;

[0104] Figure 8 For the relative calculation error of the hybrid sub-grid ADI-FDTD method at different ratios in Embodiment 5 of the present application. Detailed implementation manners

[0105] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without making creative efforts belong to the scope of protection of the present invention.

[0106] Embodiment 1. DP-ADI-FDTD method based on the z-transformed lossy multi-pole Debye model

[0107] The differential time-domain expressions of Maxwell's equations are shown as follows:

[0108] (1)

[0109] (2)

[0110] where D is the electric flux density, H is the magnetic field strength, and E is the electric field strength, is the divergence, μ 0 is the magnetic permeability of free space, and the constitutive relation between the electric flux density D and the electric field E is expressed as D = ε 0 ε r E , ε 0 is the permittivity of free space, ε r is the relative permittivity.

[0111] Discretize (1) and (2) into the following matrix form:

[0112] (3)

[0113] where V DP = E DP_x , E DP_y , E DP_z , H DP_x , H DP_y , H DP_z represents the divergence-preserving electric and magnetic field strength components, V DP The superscripts n and n + 1 of represent the iteration steps; I is the identity matrix; E DP_x , E DP_y , E DP_z are the intensities of the electric field strength in the x, y, and z directions respectively, H DP_x , H DP_y , H DP_z are the intensities of the magnetic field strength in the x, y, and z directions respectively. The matrices A and B are the spatial partial derivative matrices, specifically:

[0114] ,

[0115] where , , , and represent the partial derivatives of the field components in the x, y, and z directions respectively.

[0116] For a multi-pole Debye dispersive medium with loss terms, the constitutive relation between D and E is expressed as:

[0117] (4)

[0118] where is the conductivity, is the difference between the relative permittivities at zero frequency and infinite frequency, represents the relative permittivity at infinite frequency; is the pole relaxation time, is the number of poles in the sensitivity response; represents the imaginary part; represents the angular frequency; represents the electric flux density in the frequency domain; represents the electric field strength in the frequency domain.

[0119] According to the z-transform theory, transforming (4) from the frequency domain to the z-domain gives:

[0120] (5)

[0121] where ψ p and Q are auxiliary variables, defined as follows:

[0122] (6)

[0123] (7)

[0124] By transforming between the z-domain and the time domain, converting equations (5), (6), and (7) to the time domain and substituting into (3), we get:

[0125] (8)

[0126] where φ p = ψ p / ε 0 ε ∞ , = Q / ε 0 ε ∞ , and in addition, ε r in A and B of the matrix becomes ε ∞， is the difference between the relative permittivities at zero frequency and infinite frequency; : ψ in the z-domainp , : Q in the z domain.

[0127] To derive the numerical iteration formula for the electromagnetic field components, (8) is discretized into two sub-time-step iterations

[0128] For the first sub-time period ( ):

[0129] (9)

[0130] (10)

[0131] For the second sub-time period ( ):

[0132] (11)

[0133] (12)

[0134] where V = E x , E y , E z , H x , H y , H z , E x 、 E y 、 E z represent the components of the electric field that does not obey the divergence property in the x, y, and z directions respectively, H x 、 H y 、 H z represent the components of the magnetic field that does not obey the divergence property in the x, y, and z directions respectively; 、 φ p 、 V DP 、 The superscripts of all represent the iteration steps.

[0135] Substituting (10) into (11) gives:

[0136] (13)

[0137] Define the auxiliary variableU = U E U H = U ex , U ey , U ez , U hx , U hy , U hz , U E 、 U H respectively represent the auxiliary variables of the electric field and the magnetic field, U ex 、 U ey 、 U ez respectively represent the components of the auxiliary variable of the electric field in the x, y, and z directions, U hx 、 U hy 、 U hz respectively represent the components of the auxiliary variable of the magnetic field in the x, y, and z directions; where . Then (13) can be reformulated as:

[0138] (14)

[0139] Substituting the previous time step of (9) into (12) gives:

[0140] (15)

[0141] Let U n+1 = V n+1 + V n+1 / 2 , then (15) can be reformulated as:

[0142] (16)

[0143] Finally, formula (14) and formula (16) constitute the basic numerical discretization framework of the lossy multi-pole Debye model ZT-DP-ADI-FDTD method.

[0144] It should be noted that, compared with the traditional ADI-FDTD method, this method introduces a matrix-free right-hand side, significantly improving the computational efficiency. According to Equation (14), we obtain U E and U H The numerical discretization equations for the field components are as follows:

[0145] (17)

[0146] (18)

[0147] ψ p 、 U E 、 U H The superscripts of

[0148] (19)

[0149] (20)

[0150] Taking the x-direction as an example, substituting Equations (18), (19), and (20) into Equation (17), we obtain the numerical iteration equation at the n+1 / 2 time step:

[0151] (21)

[0152] where , , , ; represents the grid size adopted in the y-direction, : represents the magnetic field component in the z-direction, n is the iteration step, : auxiliary variable The x-direction component of 、 、 represent the coordinate indices in the x, y, and z directions respectively. It should be noted that is interpreted as the imaginary part in Equation (4), and as the coordinate index in the y-direction in the remaining equations.

[0153] The field components of

[0154] (22)

[0155] Similarly, and The numerical iteration formula is discretized as follows:

[0156] (23)

[0157] H y represents the magnetic field component in the y - direction; U ex The superscript represents the iteration step number; represents the grid size in the z - direction.

[0158] Where , .

[0159] (24a)

[0160] In addition, and The iteration formula is:

[0161] (24b)

[0162] 4c)

[0163] U ex ,、 U ey 、 U ez 、 E x 、 E y 、 E z The superscripts all represent the iteration step number.

[0164] Next, substituting the condition into Equation (18) gives:

[0165] (25)

[0166] Using a similar method, the calculation formula for can be obtained as:

[0167] (26)

[0168] The numerical iteration formulas for the electric field and magnetic field are given by Equation (24a), Equation (24b), Equation (24c), and Equation (26) respectively. Finally, based on Equation (9), the calculation formulas for the divergence - free electric field E DP and magnetic field H DP are derived:

[0169] (27)

[0170] Among them, the superscripts of E and H represent the iteration steps.

[0171] Example 2 Analysis of the divergence-preserving property of the ZT-DP-ADI-FDTD method in Example 1

[0172] The divergence condition that satisfies Gauss's law is a necessary condition for solving Maxwell's equations. Although the DP-ADI-FDTD method has the property of preserving divergence when dealing with non-dispersive media, which has been numerically verified [D. N. Smithe, J. R. Cary, and J. A. Carlsson, “Divergence preservation in the ADI algorithms for electromagnetics,” J. Comput. Phys., vol. 228, no. 19, pp. 7289-7299, 2009.], whether this property is still maintained when dealing with lossy multi-pole Debye dispersive media needs to be further demonstrated.

[0173] To analyze the divergence-preserving property of the proposed ZT-DP-ADI-FDTD method, (8) is reformulated as:

[0174] (28)

[0175] Where is the current source. Here, similar to the traditional FDTD method, the current source is added at half time steps.

[0176] The right side of formula (28) is further derived as:

[0177] (29)

[0178] The left side of formula (28) is further derived as:

[0179] (30)

[0180] Substituting formula (29) and formula (30) into formula (31), we get:

[0181] (31)

[0182] It can be seen from (3) that (A + B) is a numerical curl operator ( ). Taking the numerical divergence ( ) of (A + B) for any array of field components, for any field component , the following condition exists:

[0183] (32)

[0184] Represents the divergence central difference operator.

[0185] Therefore, the numerical divergence on both sides of Equation (31) gives the following equation:

[0186] (33)

[0187] Substitute the constitutive relation into Equation (33) to get:

[0188] (34)

[0189] It can be seen from Equation (34) that the numerical divergence of D satisfies Gauss's law, indicating that the ZT-DP-ADI-FDTD method proposed in this paper maintains the divergence property when simulating lossy multi-pole Debye dispersive media. In addition, Equation (34) shows that when there is no current source in the computational domain, the numerical divergence of D is zero. This is a key property of EMPIC simulation because it can ensure the prevention of charge accumulation during the simulation.

[0190] To further verify this property, a numerical simulation of homogeneous soil was carried out in a 40×40×40 grid simulation domain. The spatial step used in the simulation is 1 mm, and the dispersion characteristics of the soil are modeled using the second-order Debye model with the following parameters ε ∞ = 3.2, ∆ ε 1 = 0.75, ∆ ε 2 = 0.3, τ 1 = 2.71 ns, τ 2 = 0.108 ns, σ = 0.397´10 -3 S / m. The differential Gaussian source is located at the center of the computational domain, and its time-domain expression is:

[0191] (35)

[0192] where τ = 0.5309 ns, t 0 = 4 τ .

[0193] The formula for the normalized divergence of the electric flux density D is:

[0194] (36)

[0195] Represents the maximum value of the absolute value of the electric flux density D.

[0196] Figure 1 Gives the divergence of D obtained by the traditional FDTD method (CFLN = 1) and the ZT-DP-ADI-FDTD method proposed in this paper (CFLN = 1, 2, 4), where . is the maximum time step allowed by the traditional finite-difference time-domain method. The numerical divergence of the electric flux D obtained by the ZT-DP-ADI-FDTD method is close to that of the traditional FDTD method. Specifically, at the source excitation position, the divergence values generated by both methods are around -30 dB, while in the passive region, the divergence of D drops below -200 dB. This indicates that, similar to the traditional FDTD method, the ZT-DP-ADI-FDTD method proposed in this embodiment satisfies the divergence condition of Gauss's law when numerically solving Maxwell's equations.

[0197] Example 3 Analysis of the unconditional stability characteristics of the ZT-DP-ADI-FDTD method in Example 1

[0198] In this example, the Fourier stability analysis method is used to study the stability of the ZT-DP-ADI-FDTD method. According to formulas (14) and (16), the matrix form of the time stepping of the electromagnetic field from time n to n + 1 can be obtained:

[0199] (37)

[0200] where .

[0201] , ,

[0202] , .

[0203] where respectively represent the x, y, and z components of the auxiliary variables and . The subscripts x 1, x 2, y 1, y 2, z 1, z 2 in 1 and 2 respectively represent the first and second terms of the Debye dispersion; , , , , , ,

[0204] , , , . and both represent a 6×3 matrix.

[0205] The discrete form of the three-dimensional plane wave of the field component in the spatial domain is defined as:

[0206] (38)

[0207] k x , k y , k z represent the wave numbers in the x, y, and z directions respectively; represents the grid size in the x direction.

[0208] Using the central difference method, the spatial partial derivative of the field component is approximated as:

[0209] (39)

[0210] where , ( represent the x, y, and z directions of the coordinate axes respectively), is the wave number in the

[0211] According to the Fourier stability analysis method, as long as the modulus of the eigenvalue of the amplification matrix M is not greater than 1, the stability of the difference equation can be guaranteed. We use a numerical method to solve the modulus of the eigenvalue of matrix M, , and the parameters of the Debye medium are set as ε ∞ = 3.2, ∆ ε 1 = 0.75, ∆ ε 2 = 0.3, τ 1 = 2.71 ns, τ 2 =0.108 ns, σ = 0.397´10 -3 S / m. In addition, the wave number in the matrix is unknown. According to equation (39), is a periodic sine function, and the matrix eigenvalues obtained in [0, π] and [π, 2π] are the same. Therefore, in this embodiment, through Calculate the modulus of the eigenvalues of the M matrix by uniformly sampling in the interval [0, π]. Figure 2 It shows the modulus of the eigenvalues of the amplification matrix M when the time step of the ZT-DP-ADI-FDTD method takes 1, 3, and 5 times the time step of the traditional FDTD method at the stability limit (i.e., CFLN = 1, 3, 5). It can be seen that the modulus of all eigenvalues is within the unit circle, indicating that the proposed method is unconditionally stable.

[0212] Example 4. Three-dimensional hybrid sub-grid technology based on dispersive DP-ADI-FDTD and FDTD methods

[0213] In this example, the ZT-DP-ADI-FDTD method is combined with the traditional FDTD method to propose a hybrid sub-grid modeling and simulation method. In the coarse grid region, the traditional FDTD method is used to calculate the electromagnetic field components, and in the fine grid region, the ZT-DP-ADI-FDTD method is used to calculate the electromagnetic field components. In the three-dimensional sub-grid technology, there are six coarse-fine grid interfaces. For simplicity, this example describes the interpolation process of the E x ( e x ) component and applies a similar interpolation scheme to other electromagnetic field components and interfaces.

[0214] Generally, the excitation source is placed in the coarse grid region. When the electromagnetic wave propagates to the interface region between the coarse grid and the fine grid, the electromagnetic field values in the coarse grid at the interface are interpolated into the corresponding fine grid components. This process can ensure that the fine grid region receives the "excitation source" at the interface. As shown in Figure 3 (a), the initial value of the fine grid electric field component e x at the interface is unknown and cannot be directly calculated using the ZT-DP-ADI-FDTD method. To determine the e x value at the interface, the electric field component E x in the coarse grid at the same interface is interpolated to e x using spatial linear interpolation. Then, numerical iteration is performed using the ZT-DP-ADI-FDTD method to transfer the e x value at the interface to the fine grid region. Figure 3 (b) in Figure 3 shows j F = j A (j f = 1 ) at the interface between the coarse grid and the fine grid of E x distribution, where the subscript F represents the coarse grid coordinates and f represents the fine grid coordinates. In the fine grid region, as shown in the figure, from point to the value can be obtained through the following spatial interpolation expression. It should be noted that in this embodiment, an odd ratio of the coarse grid to the fine grid is used for clarity. However, a similar method can also be adopted for an even ratio.

[0215] Case 1:

[0216] (40)

[0217] Case 2:

[0218] (41)

[0219] where m 1 represents the number of grids in the x - direction, , and M represents the ratio of the coarse grid to the fine grid. and ( , ) are the grid coordinates along the x - and z - directions on the coarse grid interface respectively. Similarly, and ( , ) are the grid coordinates along the x - and z - directions on the fine grid interface. It should be noted that , , and represent the front and back interfaces in the x, y, and z directions of the coarse grid region respectively.

[0220] In Figure 3 (b), the electric field value in the right - hand frame corresponds to Case 1 and is interpolated using formula (40). The value in the left - hand frame corresponds to Case 2 and is interpolated using formula (41). Therefore, the following results can be obtained: , , , , , . Other field values can be obtained by a similar calculation method.

[0221] Subsequently, the iterative expression of needs to be corrected. Taking equation (21) at the n + 1 / 2 time step as an example: ​

[0222] (42)

[0223] Here, a, b, and c have the same definitions as a, b, and c in formula (21). Denotes the index in the z - direction of the fine grid, m 2 Denotes the number of fine grids in the k - direction; Denotes the x - component of the auxiliary variable u of the electric field in the fine grid; Denotes the x - component of the auxiliary variable in the fine - grid region, and the superscript represents the iteration step number.

[0224] For :

[0225] (43)

[0226] Where . The correction equations for the electric field in other directions also have a similar form. Denotes the x - component of the magnetic field in the fine - grid region, and the superscript represents the iteration step number.

[0227] After one iteration of the electromagnetic field components in the fine - grid region, the values of the fine - grid electromagnetic field at the fine - coarse grid junctions and adjacent junctions must be interpolated back to the values of the coarse - grid electromagnetic field at the corresponding junctions. Taking the x - direction as an example, the interpolation formulas for other directions are obtained in a similar way.

[0228] If M is odd,

[0229] (44)

[0230] (45)

[0231] If M is even,

[0232] (46)

[0233] (47)

[0234] Denotes the x - component of the magnetic field in the fine - grid region that obeys the divergence property, and the superscript represents the iteration step number; Denotes the x - component of the electric field in the fine - grid region that obeys the divergence property, and the superscript represents the iteration step number.

[0235] In summary, the numerical implementation of the hybrid sub - grid technique that combines the traditional FDTD method with the ZT - DP - ADI - FDTD method proposed in this embodiment can be summarized as follows:

[0236] 1) Update the values of the electric field (E) and magnetic field (H) for the entire coarse grid using the traditional finite-difference time-domain method;

[0237] 2) Formulas (40), (41), (42), and (43) are used to calculate the values of the fine-grid electric field (e) at the interface;

[0238] 3) Formulas (21), (22), (23), (24a), (24b), and (24c) are used to calculate the values of the fine-grid electric field (e) in the inner fine-grid region;

[0239] 4) Formulas (26) and (27) are used to calculate the values of the magnetic field (h) for the entire fine-grid region;

[0240] 5) Formulas (44) and (45) or (46) and (47) are used to correct the values of the coarse-grid electric field (E) and magnetic field (H), especially at and near the interface between the coarse-grid and fine-grid regions.

[0241] Example 5. Numerical Example

[0242] To verify the accuracy and efficiency of the hybrid sub-grid method proposed in Example 4 for simulating multi-target ground-penetrating radar scenarios, dielectric cylinders and dielectric spheres were considered in the same soil layer. As Figure 4 shown, the sub-grid size was set to 0.6 m × 0.6 m × 0.6 m, and the ratio of the coarse grid to the fine grid was set to 3 (M = 3). Ten layers of CPML were used, and the transmitter and receiver were located 0.18 meters from the soil surface, with a separation of 0.9 meters between them. The centers of the sphere and the cylinder were located 0.9 meters below the soil surface. The relative dielectric constant of the dielectric sphere was 30, and its radius was 0.24 meters, while the relative dielectric constant of the dielectric cylinder was 50, its radius was 0.18 meters, and its length was 0.48 meters. In addition, to demonstrate the flexibility of the proposed method, simulations were performed for different coarse-to-fine grid ratios (M = 3, 4, 5) to verify the method's ability to handle odd and even grid ratios.

[0243] Figure 5 is the time-domain waveform of the normalized electric field Ez at the receiving antenna when the distance between the transmitting antenna and the receiving antenna varies from 0 to 0.78 m. It can be seen that as the Tx-Rx distance increases, the electromagnetic radiation detected by the receiving antenna also increases proportionally. Figure 6 The snapshots of the electric field shown further demonstrate that there is no reflection when the electric field propagates from the coarse-grid region to the fine-grid region, thus verifying the accuracy and effectiveness of the proposed hybrid sub-grid method.

[0244] The calculation results are as Figure 7As shown, it indicates that the results obtained by the hybrid sub-grid discretization method with different coarse-to-fine grid ratios are consistent with the uniform fine-grid discretization method. Figure 8 The relative computational errors of the hybrid sub-grid method at different coarse-to-fine grid ratios are given. As the coarse-to-fine grid ratio increases, the relative computational error decreases, indicating that finer discretization in the target area improves the accuracy of the results, which is in line with the expected behavior.

[0245] Finally, this embodiment compares the computational time and memory usage of the uniform fine-grid method and the hybrid sub-grid method with three different coarse-to-fine grid ratios. As shown in Table 1, compared with the uniform fine-grid method, the hybrid sub-grid method significantly reduces the computational time and memory usage at all grid ratios. These results further verify the accuracy and efficiency of the hybrid sub-grid method. A more in-depth analysis of Table 1 shows that as the coarse-to-fine grid ratio increases, the computational time and memory usage also increase. This is because, with a fixed fine-grid area, a smaller grid size requires solving more unknowns, resulting in higher computational time and memory requirements. Specifically, in the ZT-DP-ADI-FDTD method, numerical operations such as matrix inversion amplify the impact of the increasing number of unknowns on the computational time and memory usage. Combining the data in Figure 8 and Table 1, it can be concluded that although increasing the coarse-to-fine grid ratio can improve computational accuracy, it also leads to a decrease in computational efficiency and an increase in memory consumption. Therefore, when applying the hybrid sub-grid method to the numerical simulation of ground penetrating radar, carefully choosing the ratio of the coarse grid to the fine grid is crucial for achieving the best balance among computational efficiency, memory usage, and accuracy.

[0246] Table 1

[0247]

[0248] In each embodiment of the present application, the superscripts n, n + 1, and n + 1 / 2 of each symbol represent the corresponding iteration steps of the symbol.

[0249] In the description of this specification, the descriptions referring to terms such as "one embodiment", "example", "specific example", etc. mean that the specific features, structures, materials, or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials, or characteristics described can be combined in a suitable manner in any one or more embodiments or examples.

[0250] The basic principles, main features and advantages of the present invention have been shown and described above. Those skilled in the art should understand that the present invention is not limited by the above embodiments, and what is described in the above embodiments and the specification only illustrates the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention will have various changes and improvements, and all these changes and improvements fall within the scope of the present invention claimed.

Claims

1. A numerical solution method for ground penetrating radar based on divergence-preserving ADI-FDTD, characterized in that: The following steps are involved: The differential time domain of Maxwell's equations for ground penetrating radar is listed; the differential time domain of the Maxwell's equations is discretized into a matrix form; for a multipolar Debye dispersive medium with loss terms, a constitutive relationship between electric flux density and electric field intensity is constructed; according to the z-transformation theory, the constitutive relationship is transformed from the frequency domain to the z-domain, and then converted to the time domain through the transformation between the z-domain and the time domain, and substituted into the matrix form discretized from the differential time domain of the Maxwell's equations to obtain a discrete matrix form based on the time domain; the discrete matrix form based on the time domain is discretized into two sub-time step iterations, the first sub-time period , the second sub-period , derive the numerical iteration formula of the electric and magnetic field components, where n is the number of iterations; respectively calculate the electric field and magnetic field intensity expressions of the first sub-time period and the second sub-time period, and obtain the basic numerical discretization framework of the lossy multipole Debye model ZT-DP-ADI-FDTD method; calculate the numerical discretization equations of each component of the electric field and the magnetic field according to the basic numerical discretization framework, and calculate the electric field intensity and the magnetic field intensity, and derive the electric field intensity and the magnetic field intensity that satisfy the divergence preservation; the calculation process of the basic numerical discretization framework of the lossy multipole Debye model ZT-DP-ADI-FDTD method includes the following steps: discretize the discrete matrix form based on the time domain into two sub-time step iterations, for the first sub-time period : For the second sub-period : in V = [ E x , E y , E z , H x , H y , H z ], E x , E y , E z They represent the components of the electric field in the x, y, and z directions that do not obey the divergence characteristics, H x , H y , H z They represent the components of the magnetic field in the x, y, and z directions that do not obey the divergence characteristics; , φ p , V DP , The superscripts of represent the number of iteration steps; Defining auxiliary variables U = [ U E U H ] = [ U ex , U ey , U ez , U hx , U hy , U hz ], U E , U H denote the auxiliary variables of the electric and magnetic fields, respectively, U ex , U ey , U ez denote the components of the auxiliary variable of the electric field in the x, y, and z directions, respectively, U hx , U hy , U hz Respectively represent the components of the auxiliary variables of the magnetic field in the x, y, and z directions; ; then restate it as: make U n+1 = V n+1 + V n+1 / 2 , then it can be restated as: .

2. The numerical solution method of ground penetrating radar based on divergence-preserving ADI-FDTD according to claim 1 is characterized in that: The differential time domain expression of the Maxwell equations is as follows: Where D is the electric flux density, H is the magnetic field intensity, and E is the electric field intensity. μ 0 is the free space permeability, and the constitutive relationship between D and E is expressed as D = ε 0 ε r E , ε 0 is the dielectric constant of free space, ε r is the relative dielectric constant.

3. The numerical solution method of ground penetrating radar based on divergence-preserving ADI-FDTD according to claim 1 is characterized in that: The matrix form is: in V DP =[ E DP_x , E DP_y , E DP_z , H DP_x, H DP_y , H DP_z ] represents the divergence-preserving electric and magnetic field intensity components, V DP The superscript n and n+1 indicate the number of iteration steps; I is the identity matrix; E DP_x , E DP_y , E DP_z are the electric field strengths in the x, y, and z directions, respectively. H DP_x ,H DP_y , H DP_z are the strengths of the magnetic field in the x, y, and z directions respectively. Matrices A and B are spatial partial derivative matrices, specifically: , ,in , , , and denote the partial derivatives of the field components in the x, y, and z directions, respectively.

4. The numerical solution method of ground penetrating radar based on divergence-preserving ADI-FDTD according to claim 1 is characterized in that: The constitutive relationship between electric flux density and electric field intensity is constructed as follows: in, is the conductivity, is the difference between the relative permittivity at zero frequency and infinite frequency, Relative permittivity at infinite frequency; is the extreme relaxation time, is the number of poles in the sensitivity response; represents the imaginary part; represents the angular frequency; Represents the electric flux density in the frequency domain; Represents the electric field strength in the frequency domain.

5. The numerical solution method of ground penetrating radar based on divergence-preserving ADI-FDTD according to claim 4 is characterized in that: The said transforming the constitutive relation from the frequency domain to the z domain, and then converting it to the time domain through the transformation between the z domain and the time domain, comprises the following steps: according to the z-transformation theory, transforming the constitutive relation of the electric flux density and the electric field intensity from the frequency domain to the z domain, and obtaining: where ψ p and Q are auxiliary variables, defined as follows: Through the transformation between z domain and time domain, we can get: in φ p = ψ p / ε 0 ε ∞ , = Q / ε 0 ε ∞ In addition, the matrix A and B ε r becomes ε ∞， is the difference between the relative permittivity at zero frequency and infinite frequency; : In the z domain ψ p , : Q in the z domain.

6. The numerical solution method of ground penetrating radar based on divergence-preserving ADI-FDTD according to claim 1 is characterized in that: get U E and U H Numerical discretization equations for the field components: ψ p , U E , U H The superscripts of represent the number of iteration steps; among them: Numerical iterative formulas for finding electric and magnetic fields: also, and The iteration formula is: U ex ,、 U ey , U ez , E x , E y , E z The superscripts of represent the number of iteration steps. , , Respectively represent the coordinate indexes in the x, y, and z directions; Substitution U H The numerical discretization equations of the field components are obtained: get The calculation formula is: Get the divergence-preserving E DP and H DP The calculation formula of the field is: H represents the magnetic field intensity, E represents the electric field intensity, and the superscripts of E and H represent the number of iteration steps.

7. A ground penetrating radar numerical solution system based on divergence-preserving ADI-FDTD, characterized in that: include: Maxwell equations calculation module, used to list the differential time domain of Maxwell equations of ground penetrating radar; discretize the differential time domain of Maxwell equations into matrix form; constitutive relation construction module of electric flux density and electric field intensity, for multi-pole Debye dispersive media with loss terms, construct the constitutive relation of electric flux density and electric field intensity; z-domain transformation module: used to transform the constitutive relation from frequency domain to z-domain according to z-transformation theory, and then convert it to time domain through the transformation between z-domain and time domain, and substitute it into the matrix form discretized from the differential time domain of Maxwell equations to obtain a discrete matrix form based on time domain; basic numerical discretization framework construction module: discretize the discrete matrix form based on time domain into two sub-time step iterations, the first sub-time period , the second sub-period , derive the numerical iterative formula of the electric and magnetic field components, where n is the number of iterations; calculate the electric field and magnetic field intensity expressions of the first sub-time period and the second sub-time period respectively, and obtain the basic numerical discretization framework of the lossy multipole Debye model ZT-DP-ADI-FDTD method; and, a divergence-preserving electric field intensity and magnetic field intensity solution module: calculate the numerical discretization equations of each component of the electric field and the magnetic field according to the basic numerical discretization framework, and calculate the electric field intensity and the magnetic field intensity, and derive the electric field intensity and the magnetic field intensity that satisfy the divergence preservation; wherein the calculation process of the basic numerical discretization framework of the lossy multipole Debye model ZT-DP-ADI-FDTD method includes the following steps: discretize the discrete matrix form based on the time domain into two sub-time step iterations, for the first sub-time period : For the second sub-period : in V = [ E x , E y , E z , H x , H y , H z ], E x , E y , E z They represent the components of the electric field in the x, y, and z directions that do not obey the divergence characteristics, H x , H y , H z They represent the components of the magnetic field in the x, y, and z directions that do not obey the divergence characteristics; , φ p , V DP , The superscripts of represent the number of iteration steps; Defining auxiliary variables U = [ U E U H ] = [ U ex , U ey , U ez , U hx , U hy , U hz ], U E , U H denote the auxiliary variables of the electric and magnetic fields, respectively, U ex , U ey , U ez denote the components of the auxiliary variable of the electric field in the x, y, and z directions, respectively, U hx , U hy , U hz Respectively represent the components of the auxiliary variables of the magnetic field in the x, y, and z directions; ; then restate it as: make U n+1 = V n+1 + V n +1 / 2 , then it can be restated as: .

8. A ground penetrating radar numerical solution method based on three-dimensional mixed subgrids, characterized in that: The three-dimensional hybrid subgrid includes a plurality of coarse grids and fine grids, and the excitation source is placed in the coarse grid area. When the electromagnetic wave propagates to the interface area between the coarse grid and the fine grid, the electromagnetic field value in the coarse grid at the interface is interpolated into the corresponding fine grid component; the method includes the following steps: using the finite difference time domain method to update the electric field value and magnetic field value of the entire coarse grid; calculating the fine grid electric field value and magnetic field value in the fine grid area by using the numerical solution method for ground penetrating radar based on divergence-preserving ADI-FDTD as described in any one of claims 1 to 6; interpolating the electric field value and magnetic field value in the coarse grid to the interface of the fine grid at the same interface by means of spatial linear interpolation; using the numerical solution method for ground penetrating radar based on divergence-preserving ADI-FDTD as described in any one of claims 1 to 6 to perform numerical iteration, and transferring the electric field component and magnetic field component at the fine grid interface to the fine grid area; after completing the iteration in the fine grid area, interpolating the values ​​of the fine grid electromagnetic field at the junction between the coarse grid and the fine grid and the adjacent junction back to the coarse grid electromagnetic field value at the corresponding junction.

9. A ground penetrating radar, characterized in that: It comprises a storage medium storing instructions capable of running the ground penetrating radar numerical solution method based on divergence preserving ADI-FDTD as described in any one of claims 1 to 6 or the ground penetrating radar numerical solution method based on three-dimensional mixed subgrids as described in claim 8, or the ground penetrating radar numerical solution system based on divergence preserving ADI-FDTD as described in claim 7.

Images