Ground penetrating radar numerical solution method and system based on divergence-preserving ADI-FDTD

Abstract

The invention discloses a ground penetrating radar numerical solution method and system based on divergence-keeping ADI-FDTD, and belongs to the field of efficient numerical simulation of electromagnetic waves in ground penetrating radars. The method comprises the following steps: listing a differential time domain of a Maxwell equation set of the ground penetrating radar; discretizing the differential time domain into a matrix form; constructing a constitutive relation between the electric flux density and the electric field intensity for the multipolar Debye dispersive medium with the loss item; transforming the constitutive relational expression from a frequency domain to a z domain, transforming the constitutive relational expression into a time domain through transformation between the z domain and the time domain, substituting the constitutive relational expression into a matrix form discretized by the differential time domain to obtain a discrete matrix form based on the time domain, and discretizing the discrete matrix form into two sub-time step iterations to obtain a basic numerical discretization framework; and calculating a numerical discretization equation of each component of the electric field and the magnetic field according to the basic numerical discretization framework, calculating the electric field intensity and the magnetic field intensity, and deducing the electric field intensity and the magnetic field intensity meeting the divergence keeping degree.

Description

Technical Field

[0001] The present invention relates to the field of efficient numerical simulation of electromagnetic waves in ground penetrating radar, and in particular 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 technology that uses antennas to transmit and receive electromagnetic (EM) waves, which can detect material properties and distribution patterns in underground media. When EM waves pass through the ground, their paths, field strengths, and waveforms 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 pre-calculated models, the spatial position, structure, and distribution of underground interfaces or geological structures can be inferred. EM simulation generates a database of reflected wave signals corresponding to underground conditions to facilitate forward modeling. Maxwell's equations control the propagation of EM fields through these media, and solving these equations is the key to accurately simulate ground penetrating radar.

[0003] The FDTD method is a widely adopted method for discretizing Maxwell's equations in both time and space domains. By doing so, it enables the computation of spatial EM fields along the time axis, 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, and 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 simulations, the dispersive properties of soils and targets, combined with the irregular structural features in certain areas, require the use of fine, uniform Yee meshes to achieve numerical accuracy. However, this leads to two major challenges: 1) The time step must be small enough to satisfy the CFL stability condition, resulting in longer calculation time; 2) Uniformly fine grid discretization increases the number of unknowns and significantly increases memory usage.

[0005] To address these challenges, several unconditionally stable FDTD methods, such as ADI-FDTD and LOD-FDTD methods, DP-ADI-FDTD method is a divergence-preserving, unconditionally stable method, which loses the divergence-preserving property when solving Maxwell equations, resulting in charge accumulation in the passive region, was introduced to solve the EM particle in cell (EMPIC) problem. This method was later generalized to simulate multipolar Debye dispersive media. Despite its advantages, the numerical implementation based on the ADE method is still complex and cumbersome.

[0006] In addition, subgridding 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 areas where high resolution is required, while using coarser grids in less critical areas, thus preventing unnecessary oversampling. Previous subgridding methods can be divided into two types: 1) One method uses a time step that satisfies the stability condition of the fine grid over the entire domain, where spatial interpolation is used to exchange EM field information between coarse-grid and fine-grid regions; 2) Another method uses different time steps for coarse-grid and fine-grid regions, each time step satisfies its own CFL stability condition, and uses spatiotemporal interpolation to exchange EM information between grids.

[0007] Recently, unconditionally stable FDTD methods have been applied to the fine-grid region, allowing the entire domain to be solved using time steps that satisfy the coarse-grid CFL stability condition.

[0008] These subgrid techniques allow for larger time steps, which greatly improves computational efficiency compared to traditional hybrid subgrid methods. They also use a uniform time step, eliminating the need for time interpolation and simplifying numerical implementation. However, existing methods still face some challenges, including: 1) It mainly focuses on two-dimensional models, which have limited applicability to complex dispersive media; 2) Lack of a unified, arbitrarily scaled electromagnetic exchange formula at the coarse-fine grid boundary; 3) Lack of the divergence preserving property of the hybrid subgrid method used in the unconditionally stable FDTD method. Summary of the invention

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

[0010] The purpose of the present invention can be achieved through the following technical solutions: The first aspect of the present invention relates to a numerical solution method for ground penetrating radar based on divergence-preserving ADI-FDTD, comprising the following steps: List the time-domain derivatives of Maxwell's equations for ground penetrating radar; Discretize the differential time domain of the Maxwell equations into a matrix form; For multipolar Debye dispersive media with loss terms, the constitutive relation between electric flux density and electric field intensity is constructed; According to the z-transformation theory, the constitutive relation 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 obtained by the differential time-domain discretization of the Maxwell 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, n is the number of iterations; calculate the electric 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 ZT-DP-ADI-FDTD method for the lossy multipole Debye model; The numerical discretization equations of the electric field and magnetic field components are calculated according to the basic numerical discretization framework, and the electric field intensity and magnetic field intensity are calculated to derive the electric field intensity and magnetic field intensity that satisfy the divergence preservation.

[0011] Optionally, the differential time-domain expression of the Maxwell equations is as follows: (1) (2) Where D is the electric flux density, H is the magnetic field intensity, E is the electric field intensity, μ0 is the free space magnetic permeability, and the constitutive relationship between the electric flux density D and the electric field E is expressed as D = ε0ε r E, ε0 is the dielectric constant of free space, ε r is the relative dielectric constant.

[0012] Optionally, the matrix form is: (3) 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 intensity components, V DP The superscripts n and n+1 indicate the number of iterations; I is the identity matrix; EDP_x , E DP_y , E DP_z are the electric field strengths in the x, y, and z directions, respectively, and 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.

[0013] Optionally, the constitutive relationship between electric flux density and electric field intensity is constructed as follows: (4) 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.

[0014] Optionally, the transforming of the constitutive relation from the frequency domain to the z domain, and then converting to the time domain through the transformation between the z domain and the time domain, comprises the following steps: According to the z-transform theory, (4) is transformed from the frequency domain to the z-domain to obtain: (5) where ψ p and Q are auxiliary variables, defined as follows: (6) (7) By transforming between the z domain and the time domain, formula (5), formula (6) and formula (7) are converted to the time domain and substituted into (3), we get: (8) where φ p =ψ p / ε0ε ∞ , = Q / ε0ε∞ , in addition, the ε in matrices A and B r becomes ε ∞， is the difference between the relative permittivity at zero frequency and infinite frequency; : ψ in z domain p , : Q in the z domain.

[0015] 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: Formula (8) is discretized into two sub-time steps. For the first sub-time period : (9) (10) For the second sub-period :

[0016] Where 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, respectively. 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; Substituting (10) into (11), we obtain: (13) Define auxiliary variable U = [U E U H ] = [U ex , U ey , U ez , U hx , U hy , U hz ],U E , U H Denote auxiliary variables of electric and magnetic fields, U ex , Uey , U ez Respectively represent the auxiliary variables of the electric field in the x, y, and z directions, 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 (13) can be restated as: (14) Substitute the previous time step of (9) into formula (12) to obtain: (15) Let U n+1 = V n+1 +V n+1 / 2 , then formula (15) can be restated as: (16).

[0017] Optionally, according to formula (14), U E and U H Numerical discretization equations for the field components:

[0018] ψ p , U E , U H The superscripts of represent the number of iteration steps, where:

[0019] Numerical iterative formulas for finding electric and magnetic fields: (24a) also, and The iteration formula is:

[0020] U ex ,、U ey , U ez 、E x 、E y 、E z The superscripts of represent the number of iteration steps. , , Represents the coordinate index in the x, y, and z directions respectively; The conditions Substituting into formula (18) we get: (25) get The calculation formula is: (26) Substituting formula (24a), formula (24b), formula (24c) and formula (26) into formula (9), we can obtain E that satisfies the divergence preservation. DP and H DP The calculation formula of the field is: (27) H represents magnetic field, E represents electric field, and the superscripts of E and H represent the number of iteration steps.

[0021] The second aspect of the present invention relates to a ground penetrating radar numerical solution system based on divergence-preserving ADI-FDTD, comprising: A Maxwell equations calculation module is used to list the differential time domain of the Maxwell equations of the ground penetrating radar; discretize the differential time domain of the Maxwell equations into a matrix form; The module for constructing the constitutive relation between electric flux density and electric field intensity is used to construct the constitutive relation between electric flux density and electric field intensity for multipolar Debye dispersive media with loss terms. Z-domain transformation module: used for transforming the constitutive relation from the frequency domain to the z-domain according to the z-transformation theory, and then converting it to the time domain through the transformation between the z-domain and the time domain, and substituting it into the matrix form obtained by the differential time-domain discretization of the Maxwell equations to obtain a discrete matrix form based on the time domain; Basic numerical discretization framework building module: discretize the time-domain based discrete matrix form 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, n is the number of iterations; calculate the electric 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 ZT-DP-ADI-FDTD method for the lossy multipole Debye model; And, a divergence-preserving electric field strength and magnetic field strength solution module: calculates the numerical discretization equations of each component of the electric field and the magnetic field according to the basic numerical discretization framework, and calculates the electric field strength and the magnetic field strength, and derives the electric field strength and the magnetic field strength that satisfy the divergence preservation.

[0022] The third aspect of the present invention relates to a numerical solution method for ground penetrating radar based on a three-dimensional hybrid subgrid, wherein the three-dimensional hybrid subgrid includes a plurality of coarse grids and fine grids, and an 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 comprises the following steps: The finite-difference time-domain method is used to update the electric and magnetic field values ​​of the entire coarse grid; The fine grid electric field value and magnetic field value in the fine grid area are calculated by the above-mentioned divergence-preserving ADI-FDTD-based ground penetrating radar numerical solution method; The electric field value and the magnetic field value in the coarse grid are interpolated to the interface in the fine grid on the same interface by means of spatial linear interpolation; the above-mentioned divergence-preserving ADI-FDTD-based ground penetrating radar numerical solution method is used for numerical iteration to transfer the electric field component and the magnetic field component at the fine grid interface to the fine grid area; After the iteration is completed in the fine grid area, the values ​​of the fine grid electromagnetic field at the boundary between the coarse grid and the fine grid and at the adjacent boundary are interpolated back to the coarse grid electromagnetic field values ​​at the corresponding boundary.

[0023] 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 divergence-preserving ADI-FDTD-based ground penetrating radar numerical solution method or the above-mentioned three-dimensional mixed sub-grid-based ground penetrating radar numerical solution method, or the above-mentioned divergence-preserving ADI-FDTD-based ground penetrating radar numerical solution system.

[0024] Beneficial effects of the present invention: 1) The Z-transform DP-ADI-FDTD method (ZT-DP-ADI-FDTD) method of the present invention can effectively and accurately simulate electromagnetic wave propagation in lossy multi-level Debye dispersive media due to its unconditional stability and dispersion-preserving properties. Compared with the ADE method, the Z-transform technique converts the multi-level Debye model from the frequency domain to the Z domain, eliminating the need to store fields from multiple previous time steps while retaining the core iterative equations. This leads to a more flexible and efficient numerical implementation.

[0025] 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 uses 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

[0026] The present invention will be further described below in conjunction with the accompanying drawings.

[0027] Figure 1The electric flux density divergence (in dB) calculated by different methods in Example 2 of the present application: (a) FDTD method, (b) DP-ADI-FDTD method (CFLN = 1), (c) DP-ADI-FDTD method (CFLN = 3), (d) DP-ADI-FDTD method (CFLN = 5); Figure 2 The modulus of the eigenvalues ​​of the matrix M is amplified when the time step of the ZT-DP-ADI-FDTD method in Example 3 of the present application is 1, 3, and 5 times the time step of the traditional FDTD under the stability limit (ie, CFLN = 1, 3, and 5); Figure 3 Schematic diagram of a three-dimensional mixed subgrid in Example 4 of the present application, wherein (a) a three-dimensional interface between a coarse grid and a fine grid, (b) an E x Interpolation diagram of .

[0028] Figure 4 The three-dimensional ground penetrating radar system for multi-target detection in Example 5 of the present application, the transmitting end is Tx, and the receiving end is Rx; Figure 5 The time domain waveforms of the receiving points at different Tx-Rx antenna distances in Embodiment 5 of the present application; Figure 6 It is a time snapshot of calculating the Ez amplitude using the hybrid subgrid DP-ADI-FDTD method in Example 5 of the present application; Figure 7 is the time domain result of Ez at the receiving antenna in Example 5 of the present application; Figure 8 It is the relative calculation error of the mixed sub-grid ADI-FDTD method at different scales in Example 5 of the present application. DETAILED DESCRIPTION

[0029] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0030] Example 1. DP-ADI-FDTD method based on z-transform lossy multipole Debye model The differential time-domain expression of Maxwell's equations is as follows: (1) (2) Where D is the electric flux density, H is the magnetic field intensity, and E is the electric field intensity. is the divergence, μ0 is the free space permeability, and the constitutive relationship between the electric flux density D and the electric field E is expressed as D = ε0ε r E, ε0 is the dielectric constant of free space, ε r is the relative dielectric constant.

[0031] Discretize (1) and (2) into the following matrix form: (3) 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 intensity components, V DP The superscripts n and n+1 indicate the number of iterations; 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, and 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: ,

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

[0033] For a multipolar Debye dispersive medium with loss terms, the constitutive relationship between D and E is expressed as: (4) 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.

[0034] According to the z-transform theory, (4) is transformed from the frequency domain to the z-domain to obtain: (5) where ψ p and Q are auxiliary variables, defined as follows:

[0035] By transforming between the z domain and the time domain, formula (5), formula (6) and formula (7) are converted to the time domain and substituted into (3), we get: (8) where φ p =ψ p / ε0ε ∞ , = Q / ε0ε ∞ , in addition, the ε in matrices A and B r becomes ε ∞， is the difference between the relative permittivity at zero frequency and infinite frequency; : ψ in z domain p , : Q in the z domain.

[0036] In order to derive the numerical iterative formula for the electric and magnetic field components, (8) is discretized into two sub-time steps: For the first sub-period ( ): (9) (10) For the second sub-period ( ): (11) (12) Where 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, respectively. 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 represent the number of iteration steps.

[0037] Substituting (10) into (11), we obtain: (13) Define auxiliary variable U = [U E U H ] = [U ex , U ey , U ez , U hx , U hy , U hz ],U E , U H Denote auxiliary variables of electric and magnetic fields, U ex , U ey , U ez Respectively represent the auxiliary variables of the electric field in the x, y, and z directions, 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 (13) can be restated as: (14) Substituting the previous time step of (9) into (12) yields: (15) Let U n+1 = V n+1+Vn+1 / 2 , then (15) can be restated as: (16) Finally, equations (14) and (16) constitute the basic numerical discretization framework of the ZT-DP-ADI-FDTD method for the lossy multi-pole Debye model.

[0038] It is worth noting that compared with the traditional ADI-FDTD method, this method introduces a matrix-free right-hand side, which significantly improves the computational efficiency. According to formula (14), we get U E and U H The numerical discretization equations of the field components are as follows:

[0039] ψ p , U E , U H The superscripts of represent the number of iteration steps; among them: (19) (20) Taking the x direction as an example, substitute formula (18), formula (19) and formula (20) into formula (17) to obtain the numerical iteration equation at the n+1 / 2 time step: (twenty one) in , ; Indicates the grid size used in the y direction, : represents the magnetic field component in the z direction, n is the number of iteration steps, : Auxiliary variables The x-direction component of , n is the number of iteration steps; , , Represents the coordinate index in the x, y and z directions respectively. It should be noted that It is interpreted as the imaginary part in formula (4) and as the coordinate index in the y direction in the other formulas.

[0040] The field components of are calculated as follows: (twenty two) akin, and The numerical iteration formula is discretized as: (twenty three) H y represents the magnetic field component in the y direction; U ex The superscript of indicates the number of iteration steps; Indicates the grid size in the z direction.

[0041] in , .

[0042] (24a) also, and The iteration formula is: (24b) (24c) U ex ,、U ey , U ez 、E x 、E y 、E z The superscripts represent the number of iteration steps.

[0043] Next, the condition Substituting into formula (18) we get: (25) Using a similar method, we can get The calculation formula is: (26) The numerical iteration formulas of the electric field and magnetic field are given by formula (24a), formula (24b), formula (24c) and formula (26), respectively. Finally, based on formula (9), the electric field E that satisfies the divergence preservation is derived. DP and magnetic field H DP The calculation formula is: (27) The superscripts E and H represent the number of iterations.

[0044] Example 2 Analysis of the divergence preservation characteristics of the ZT-DP-ADI-FDTD method in Example 1 Satisfying the divergence condition of 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, it has been numerically verified [DN Smithe, JRCary, and JA Carlsson, "Divergence preservation in the ADI algorithms forelectromagnetics," J. Comput. Phys., vol. 228, no. 19, pp. 7289-7299, 2009.], whether this property is still maintained when dealing with lossy multipole Debye dispersive media needs further proof.

[0045] In order to analyze the dispersion-preserving property of the proposed ZT-DP-ADI-FDTD method, (8) is reformulated as: (28) in is a current source. As with the traditional FDTD method, a current source is added at half time.

[0046] The right side of formula (28) is further derived as: (29) The left side of formula (28) is further deduced as: (30) Substituting formula (29) and formula (30) into formula (31), we obtain: (31) From (3), we can see that (A + B) is a numerical curl operator ( ), taking the numerical divergence of (A + B) for an arbitrary array of field components ( ), for any field component , the following conditions exist: (32) represents the divergence central difference operator.

[0047] Therefore, the numerical divergence on both sides of equation (31) can be obtained as follows:

[0048] From formula (34), it can be seen 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 characteristics when simulating lossy multipole Debye dispersive media. In addition, formula (34) shows that the numerical divergence of D is zero when there is no current source in the computational domain. This is a key feature of EMPIC simulation because it ensures that charge accumulation is prevented during the simulation.

[0049] To further verify this characteristic, a numerical simulation of uniform soil was performed in a 40×40×40 grid simulation domain. The spatial step size used in the simulation was 1 mm, and the dispersion characteristics of the soil were 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: (35) where τ = 0.5309 ns, t0 = 4τ.

[0050] The formula for the normalized divergence of the electric flux density D is: (36) Indicates the maximum absolute value of the electric flux density D.

[0051] Figure 1 The divergence of D obtained by the traditional FDTD method (CFLN = 1) and the proposed ZT-DP-ADI-FDTD method (CFLN = 1, 2, 4) is given, where . tmax 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 ​​produced by both methods are around -30 dB, while in the passive area, the divergence of D drops below -200 dB. This shows 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 the Maxwell equations.

[0052] Example 3 Analysis of the unconditional stability characteristics of the ZT-DP-ADI-FDTD method in Example 1 This embodiment uses the Fourier stability analysis method to study the stability of the ZT-DP-ADI-FDTD method. According to formula (14) and formula (16), the matrix form of the time step of the electromagnetic field in the time period from n to n+1 can be obtained: (37) in .

[0053]

[0054]

[0055] in Respectively represent: auxiliary variables and The x, y, and z components of the chromatic aberration are shown in Table 1. The 1 and 2 in the subscripts x1, x2, y1, y2, z1, and z2 represent the first and second terms of the Debye dispersion, respectively.

[0056] The three-dimensional plane wave of the field component is discretely defined in the spatial domain as: (38) k x , k y , k z represent the wave numbers in the x, y, and z directions respectively; Indicates the grid size in the x direction.

[0057] Using the central difference method, the spatial partial derivatives of the field components are approximated as: (39) in , ( Represent the x, y, and z directions of the coordinate axis respectively), yes Direction of the wave number.

[0058] According to the Fourier stability analysis method, as long as the modulus of the eigenvalue of the amplified matrix M is not greater than 1, the stability of the differential equation can be guaranteed. We use a numerical method to solve the modulus of the eigenvalue of the matrix M. , the parameter of the Debye medium is set to ε ∞ = 3.2, Δε1= 0.75, Δε2= 0.3, τ1= 2.71 ns, τ2= 0.108ns, σ = 0.397×10 -3 S / m. In addition, the wave number in the matrix is unknown. According to equation (39), is a periodic sinusoidal function, and The matrix eigenvalues ​​obtained at [0, π] and [π, 2π] are the same, so this embodiment uses Compute the modulus of the eigenvalues ​​of the M matrix uniformly on [0, π]. Figure 2 The figure shows the modulus of the eigenvalues ​​of the matrix M when the time step of the ZT-DP-ADI-FDTD method is 1, 3, and 5 times the time step of the traditional FDTD in the stability limit (i.e., CFLN = 1, 3, and 5). It can be seen that the modulus of all eigenvalues ​​is within the unit circle, which shows that the proposed method is unconditionally stable.

[0059] Example 4. Three-dimensional hybrid subgrid technology based on dispersive DP-ADI-FDTD and FDTD methods In this embodiment, the ZT-DP-ADI-FDTD method is combined with the traditional FDTD method to propose a hybrid sub-grid modeling and simulation method. The traditional FDTD method is used to calculate the electromagnetic field components in the coarse grid area, and the ZT-DP-ADI-FDTD method is used to calculate the electromagnetic field components in the fine grid area. In the three-dimensional sub-grid technology, there are six coarse-fine grid interfaces. For simplicity, this embodiment describes one of the interfaces. x (e x ) components and applies a similar interpolation scheme to other electromagnetic field components and interfaces.

[0060] Usually, the excitation source is placed in the coarse mesh area. When the electromagnetic wave propagates to the interface area between the coarse mesh and the fine mesh, the electromagnetic field value in the coarse mesh at the interface is interpolated to the corresponding fine mesh component. This process ensures that the fine mesh area receives the "excitation source" at the interface. Figure 3 As shown in (a), the fine grid electric field component e at the interface x The initial value of is unknown and cannot be directly calculated using the ZT-DP-ADI-FDTD method. xThe electric field component E in the coarse grid at the same interface is calculated by spatial linear interpolation. x Interpolate to e x Then, the ZT-DP-ADI-FDTD method was used to perform numerical iterations and the e x The value is transferred to the fine grid area. Figure 3 (b) shows Figure 3 (a) in j F = j A (j f =1) at the interface between the coarse and fine grids x Distribution, where the subscript F represents the coarse grid coordinates and f represents the fine grid coordinates. In the fine grid area, as shown in the figure, from point arrive The value of can be obtained by the following spatial interpolation expression. It should be noted that an odd ratio of coarse grid to fine grid is used in this embodiment for clarity. However, a similar method can be used for even ratios.

[0061] Scenario 1:

[0062]

[0063] Scenario 2:

[0064] (41).

[0065] Where m1 represents the number of grids in the x direction, , M represents the ratio of coarse grid to fine grid. and ( , ) are the grid coordinates along the x and z directions on the coarse grid interface. Similarly, and ( , ) are the grid coordinates along the x and z directions on the fine grid interface. It is worth noting that , ,and Represent the front and rear interfaces in the x, y and z directions of the coarse mesh area, respectively.

[0066] exist Figure 3 In (b), the electric field values ​​in the right box correspond to case 1 and are interpolated using formula (40). The values ​​in the left box correspond to case 2 and are interpolated using formula (41). Therefore, the following results can be obtained: , , , , , Other field values ​​can be obtained by similar calculation methods.

[0067] Then, it is necessary to correct Iterative expression of . Take equation (21) at the n+1 / 2 time step as an example: (42) The definitions of a, b, and c here are the same as those in formula (21). represents the index of the fine grid in the z direction, and m2 represents the number of fine grids in the k direction; The x-direction component of the auxiliary variable u representing the electric field in the fine grid; Represents auxiliary variables in the fine grid region The x-component of , and the superscript indicates the number of iterations.

[0068] for :

[0069] in The correction equations for the electric field in other directions have similar forms. represents the x-component of the magnetic field in the fine mesh area, and the superscript represents the number of iterations.

[0070] After the electromagnetic field components in the fine grid area complete one iteration, the values ​​of the fine grid electromagnetic field at the boundary between the coarse and fine grids and the adjacent boundary must be interpolated back to the coarse grid electromagnetic field values ​​at the corresponding boundary. Taking the x direction as an example, the interpolation formulas in other directions use similar methods.

[0071] If M is an odd number,

[0072] If M is even,

[0073] It represents the x-component of the magnetic field in the fine mesh area that obeys the divergence characteristic, and the superscript represents the number of iteration steps; It represents the x-component of the electric field in the fine mesh region that obeys the divergence characteristic. The superscript represents the number of iterations.

[0074] In summary, the numerical implementation of the hybrid subgrid technology that combines the traditional FDTD method with the ZT-DP-ADI-FDTD method proposed in this embodiment can be summarized as follows: 1) Use the traditional finite-difference time-domain method to update the electric field (E) and magnetic field (H) values ​​of the entire coarse grid; 2) Formulas (40), (41), (42), and (43) are used to calculate the fine grid electric field (e) value at the interface; 3) Formulas (21), (22), (23), (24a), (24b), and (24c) are used to calculate the fine grid electric field (e) value within the inner fine grid region; 4) Formulas (26) and (27) are used to calculate the magnetic field (h) value of the entire fine grid area; 5) Equations (44) and (45) or (46) and (47) are used to correct the coarse-grid electric (E) and magnetic (H) field values, especially at and near the interface between the coarse-grid and fine-grid regions.

[0075] Example 5. Numerical example In order to verify the accuracy and efficiency of the hybrid subgrid method proposed in Example 4 in simulating multi-target GPR scenarios, dielectric cylinders and dielectric spheres were considered in the same soil layer. Figure 4 As shown in the figure, the subgrid size is set to 0.6 m × 0.6 m × 0.6 m, and the ratio of coarse grid to fine grid is set to 3 (M = 3). A 10-layer CPML is used, and the transmitter and receiver are located 0.18 m above the soil surface with a spacing of 0.9 m between them. The center of the sphere and cylinder is 0.9 m below the soil surface. The relative permittivity of the dielectric sphere is 30 and the radius is 0.24 m, while the relative permittivity of the dielectric cylinder is 50, the radius is 0.18 m, and the length is 0.48 m. In addition, to demonstrate the flexibility of the proposed method, simulations with different coarse and fine grid ratios (M = 3, 4, 5) are performed to verify the ability of the method to handle odd and even grid ratios.

[0076] Figure 5 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 in the range of 0 ~ 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 electric field snapshots shown further demonstrate that there is no reflection when the electric field propagates from the coarse mesh region to the fine mesh region, thus verifying the accuracy and effectiveness of the proposed hybrid sub-grid method.

[0077] The calculation results are as follows Figure 7 As shown in Figure 3, it is shown that the results obtained by the mixed subgrid discretization method with different coarse-fine grid ratios are consistent with the uniform fine grid discretization method. Figure 8 The relative computational errors of the mixed subgrid method for different coarse-fine grid ratios are given. As the coarse-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 consistent with the expected behavior.

[0078] Finally, this example compares the computation time and memory usage of the uniform fine grid method and the hybrid subgrid method with three different coarse-fine grid ratios. As shown in Table 1, the hybrid subgrid method significantly reduces the computation time and memory usage at all grid ratios compared to the uniform fine grid method. These results further verify the accuracy and efficiency of the hybrid subgrid method. A deeper analysis of Table 1 shows that as the coarse-fine grid ratio increases, the computation 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 computation time and memory requirements. Specifically, in the ZT-DP-ADI-FDTD method, numerical operations such as matrix inversion amplify the impact of the increase in the number of unknowns on computation time and memory usage. Combining the data in Figure 8 and Table 1, it can be concluded that although increasing the coarse-fine grid ratio can improve the computational accuracy, it also leads to reduced computational efficiency and increased memory consumption. Therefore, when applying the hybrid subgrid method to ground penetrating radar numerical simulations, carefully selecting the ratio of coarse grid to fine grid is crucial to achieve the best balance between computational efficiency, memory usage, and accuracy.

[0079] Table 1

[0080] In each embodiment of the present application, the superscript n, n+1, n+1 / 2 of each symbol indicates the number of iteration steps corresponding to the symbol.

[0081] In the description of this specification, the description with reference to the terms "one embodiment", "example", "specific example", etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representation of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any one or more embodiments or examples in a suitable manner.

[0082] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments, and the above embodiments and descriptions are only for explaining the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, and these changes and improvements all fall within the scope of the present invention to be protected.

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 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 ZT-DP-ADI-FDTD method of the lossy multipole Debye model; calculate the numerical discretization equations of each component of the electric field and 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.

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, E is the electric field intensity, μ0 is the free space magnetic 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: 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 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, and H DP_x, HDP_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: where φ p =ψ p / ε0ε ∞ , = Q / ε0ε ∞ , in addition, the ε in matrices A and B r becomes ε ∞， is the difference between the relative permittivity at zero frequency and infinite frequency; : ψ in 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 5 is characterized in that: 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: discretizing 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 : Where V = [Ex, Ey, Ez, Hx, Hy, Hz], Ex, Ey, Ez represent the components of the electric field in the x, y, z directions that do not obey the divergence characteristics, Hx, Hy, Hz represent the components of the magnetic field in the x, y, z directions that do not obey the divergence characteristics, respectively; ,φp ,VDP, The superscripts of represent the number of iteration steps; Define auxiliary variables U = [UE UH] = [Uex, Uey, Uez, Uhx, Uhy, Uhz], UE and UH represent auxiliary variables of electric field and magnetic field respectively, Uex, Uey and Uez represent components of auxiliary variables of electric field in x, y and z directions respectively, and Uhx, Uhy and Uhz represent components of auxiliary variables of magnetic field in x, y and z directions respectively; ; then restate it as: 。 7. The numerical solution method of ground penetrating radar based on divergence-preserving ADI-FDTD according to claim 6 is characterized in that: Get U E And the numerical discretization equations of the UH field components: The superscripts of ψp, UE, and UH all represent the number of iteration steps; among them: Numerical iterative formulas for finding electric and magnetic fields: 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; Substitute U H The numerical discretization equations of the field components are obtained: get The calculation formula is: Get E that satisfies divergence preservation 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.

8. 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 ZT-DP-ADI-FDTD method of the lossy multipole Debye model; 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 meet the divergence preservation.

9. 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 7; 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 7 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.

10. 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 7 or the ground penetrating radar numerical solution method based on three-dimensional mixed subgrids as described in claim 9, or the ground penetrating radar numerical solution system based on divergence preserving ADI-FDTD as described in claim 8.

Images