Three-dimensional ground penetrating radar dispersive medium discontinuous finite element forward modeling method and system and storage medium

Abstract

The invention discloses a three-dimensional ground penetrating radar dispersive medium discontinuous finite element forward modeling method and system and a storage medium, and the method comprises the following steps: S1, building a forward modeling model, and setting parameters at a boundary; s2, inputting parameters, and setting the positions of an excitation point and a receiving point; s3, performing mesh generation on the three-dimensional dispersive medium model; s4, solving numerical flux of six components of the three-dimensional electromagnetic field; s5, updating an auxiliary field variable of the three-dimensional area; s6, updating electromagnetic field components of the whole region; s7, time steps are added, and the steps are repeated until three-dimensional forward modeling simulation of the current time step is completed; s8, repeating the steps until all excitation is completed; and S9, outputting the data. According to the method, the three-dimensional full-region electromagnetic field evolution calculation logic is optimized, and the multi-scale modeling requirement of a large-scale complex underground structure is met. The finally output three-dimensional radar profile can present the depth, form and spatial distribution of the underground target in a three-dimensional manner, and detection blind areas are reduced.

Description

Technical Field

[0001] This invention relates to the field of physical detection technology, and in particular, to a three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method. Background Technology

[0002] Ground penetrating radar (GPR), as a high-precision non-destructive exploration technology, is widely used in fields such as oil exploration, geological disaster early warning, and cultural heritage protection. It obtains information about the structure of underground materials by emitting high-frequency electromagnetic waves and receiving reflected signals. Because the propagation of electromagnetic waves underground involves solving complex wave equations and multi-physics coupling, and because actual exploration faces challenges such as formation absorption attenuation and random noise interference, GPR forward modeling research is of great significance for constructing the mapping relationship between geoelectric models and radar response, optimizing time-shift data processing algorithms, and breaking through inversion imaging technology.

[0003] Compared to the limitations of two-dimensional GPR, which can only provide a single survey profile and is prone to decision-making bias, three-dimensional GPR can display the characteristics of underground targets in a three-dimensional way through dense grid acquisition and volume data reconstruction, reducing blind spots in detection. However, existing methods have problems such as low accuracy and incomplete consideration of the impact on boundaries. Summary of the Invention

[0004] This invention aims to at least solve one of the technical problems existing in the prior art. To this end, this invention proposes a three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0006] A three-dimensional ground-penetrating radar discontinuous medium dispersive element forward modeling method includes the following steps: S1, establishing a forward model of the three-dimensional dispersive medium and setting the CFS-PML parameters at the boundary; S2, inputting the physical property parameters of the dispersive medium and setting the excitation and receiving point positions; S3, meshing the three-dimensional dispersive medium model; S4, solving for the numerical flux of the six components of the three-dimensional electromagnetic field; S5, updating the auxiliary field variables of the three-dimensional CFS-PML region; S6, updating the electromagnetic field components of the entire region; S7, adding a time step and repeating steps S4-S6 until the three-dimensional forward modeling simulation of the current time step is completed; S8, repeating steps S4-S7 until all excitations are completed; S9, outputting data to form a radar profile of the three-dimensional dispersive medium.

[0007] Furthermore, the physical property parameter of the dispersive medium is the relative permittivity ε. r ,

[0008] Where ε ∞τ is the relative permittivity at infinite frequency, Δε is the difference between the relative permittivity at static zero frequency and at infinite frequency, τ is the relaxation time, j is the imaginary unit, and ω is the angular frequency.

[0009] Furthermore, the numerical flux in step S4 is calculated using the following formula:

[0010]

[0011] fluxH x fluxH y fluxH z These represent the numerical fluxes of the magnetic field in the x, y, and z directions, respectively.

[0012] fluxE x fluxE y fluxE z These are the numerical fluxes of the electric field in the x, y, and z directions, respectively.

[0013] κ e v h , κ h v e For numerical flux coefficients;

[0014] n x n y and n z It is the projection of the external normal vector of the tetrahedron onto the x, y, z directions;

[0015] H x H y H z These are the component values ​​of the magnetic field of the current unit in the x, y, and z directions, respectively;

[0016] E x E y E z These are the components of the electric field of the current cell in the x, y, and z directions, respectively.

[0017] These are the components of the magnetic field in the x, y, and z directions of the adjacent units, respectively.

[0018] These are the components of the electric field in the x, y, and z directions of the adjacent cells, respectively.

[0019] HH and EE represent the jump values ​​of the normal components of the magnetic field and electric field, respectively.

[0020] Furthermore, the jump values ​​of the normal components of the magnetic and electric fields can be calculated using the following formula:

[0021]

[0022] Furthermore, the auxiliary field variables include an electric field auxiliary variable P and a magnetic field auxiliary variable Q, wherein the electric field auxiliary variable P includes P0. x The magnetic field auxiliary variable Q includes Q0 x The P x and Q x Through such

[0023] The following formula has been updated:

[0024] P x =ε ∞ [(σ y -σ x -α z )P x,1 -(α y +σ x )P x,2 -(α x +σ x )P x,3 ]+P x,4 +P x,5 +P x,6 ;

[0025] Q x =μ[(σ y -σ x -α z )Q x,1 -(α y +σ x )Q x,2 -(α x +σ x )Q x,3 ] / ε0;

[0026] Among them, P x and Q x These are the electric field components E x and magnetic field component H x Auxiliary field variables; ε ∞ σ is the relative permittivity at infinite frequency; x and σ y These represent the conductivity in the x and y directions within the CFS-PML layer, respectively; α x α y and α z ε0 represents the complex frequency shift parameters in the x, y, and z directions within the CFS-PML layer, respectively; μ is the permeability; ε0 is the permittivity in vacuum; P0 is the dielectric constant in vacuum. x,1 P x,2 P x,3 P x,4 P x,5 P x,6These are the electric field components E x The six auxiliary field variable components, Q x,1 Q x,2 Q x,3 These are the magnetic field components H. x The three auxiliary field variable components.

[0027] Furthermore, the electromagnetic field components of the entire region are updated according to the following formula:

[0028]

[0029] in, Let D be the partial derivative, t be time, and D be the partial derivative. k For the current tetrahedron, M is the mass matrix of the element; S x S y and S z , respectively, are the unit rigidity matrices in the x, y, and z directions; F is the boundary element matrix; J is the integral term of the current source; z The loading stimulus source; dv is the vector composed of the values ​​of the basis functions at this element node; the superscript T indicates the transpose of the vector; dv represents the volume element; Δε is the difference in relative permittivity between the static zero frequency and the infinite frequency; σ is the conductivity.

[0030] Furthermore, the CFS-PML parameters at the boundary include the maximum value α of the complex frequency shift parameter at the outer boundary of the CFS-PML. max The maximum conductivity σ at the outer boundary of CFS-PML max The PML thickness d0 can be obtained using the following formula: σ i and α i :

[0031]

[0032] σ i α represents the conductivity along the i-direction within the CFS-PML layer. i α is the complex frequency shift parameter in the i-direction within the CFS-PML layer; max The maximum value of the complex frequency shift parameter at the outer boundary of CFS-PML; d i It is the distance to the innermost layer i of the PML in the direction; d0 is the thickness of the PML, σ max is the maximum conductivity at the outer boundary of CFS-PML, where m is the exponential order.

[0033] Furthermore, the maximum conductivity σ at the outer boundary of the CFS-PML max It is obtained by calculation using the following formula:

[0034] σmax =-(m+1)ln(R0) / (2ηd0ε r )

[0035] Where m is the exponential order; R0 is the expected reflection error; η is the wave impedance of the PML layer; ε r is the relative permittivity of the medium.

[0036] This invention also provides a three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling system, comprising: a modeling and parameter setting module for establishing a forward model of the three-dimensional dispersive medium and setting CFS-PML parameters at the boundaries; an input setting module for inputting dispersive medium physical property parameters and setting the excitation and receiving point positions; a parameter calculation module for meshing the three-dimensional dispersive medium model; a numerical flux solving module for solving the numerical flux of the six components of the three-dimensional electromagnetic field; an auxiliary field variable update module for updating the auxiliary field variables of the three-dimensional CFS-PML region; an electromagnetic field component update module for updating the electromagnetic field components of the entire region; a cyclic simulation module for increasing the time step and completing the three-dimensional forward simulation of the current time step; a cyclic excitation module for completing all excitations; and an output module for outputting data to form a radar profile of the three-dimensional dispersive medium.

[0037] The present invention also provides a storage medium, wherein the computer-readable storage medium stores the computer program, and the computer program, when executed by the processor, implements the three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method.

[0038] The present invention has the following beneficial effects:

[0039] This invention effectively solves the core problems of imperfect coupling of absorbing boundaries and limited multi-scale modeling accuracy in 3D ground-penetrating radar (GPR) forward modeling of dispersive media by constructing a 3D dispersive medium forward model and integrating CFS-PML boundary processing. By precisely setting CFS-PML parameters and updating the corresponding auxiliary field variables, non-physical reflection interference is significantly reduced, ensuring the efficiency and stability of boundary processing. The calculation logic for the evolution of the electromagnetic field across the entire 3D region is optimized through step-by-step iterative updates of electromagnetic field components and numerical flux, improving simulation accuracy and adapting to the multi-scale modeling needs of large-scale complex underground structures. The final output 3D radar profile can present the depth, morphology, and spatial distribution of underground targets in a three-dimensional manner, reducing detection blind spots and providing reliable support for geoelectric model construction, optimization of time-shift data processing algorithms, and breakthroughs in inversion imaging technology. This significantly enhances the application value and decision-making accuracy of GPR in practical exploration scenarios.

[0040] In addition to the objectives, features, and advantages described above, the present invention has other objectives, features, and advantages. The invention will now be described in further detail with reference to the figures. Attached Figure Description

[0041] The accompanying drawings, which form part of this application, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings:

[0042] Figure 1 This is a schematic diagram of the overall process of the present invention;

[0043] Figure 2 This is a three-dimensional geological model of a goaf area in one embodiment;

[0044] Figure 3 This is a forward-modeling slice of a three-dimensional geological model of a goaf area in one of the embodiments. Detailed Implementation

[0045] It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

[0046] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0047] It should be noted that all directional indications (such as up, down, left, right, front, back, etc.) in the embodiments of the present invention are only used to explain the relative positional relationship and movement of each component in a certain specific posture (as shown in the figure). If the specific posture changes, the directional indication will also change accordingly.

[0048] Furthermore, the use of terms such as "first" and "second" in this invention is for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined with "first" or "second" may explicitly or implicitly include at least one of that feature. Additionally, the technical solutions of the various embodiments can be combined with each other, but only on the basis of being achievable by those skilled in the art. When the combination of technical solutions is contradictory or impossible to implement, such a combination of technical solutions should be considered non-existent and not within the scope of protection claimed by this invention.

[0049] Please refer to Figure 1 The present invention provides a preferred embodiment of a three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method, comprising steps S1, S2, S3, S4, S5, S6, S7, S8 and S9.

[0050] S1. Establish a forward model of the three-dimensional dispersive medium and set the CFS-PML parameters at the boundary. The forward model of the three-dimensional dispersive medium is a virtual model used to simulate the wave / field propagation response of a medium with dispersive characteristics in three-dimensional space, and is used for forward calculations.

[0051] S2: Input the physical properties of the dispersive medium and set the excitation and receiving points. The excitation and receiving points are the radar's transmission and reception points, and need to be preset.

[0052] S3 meshes the 3D dispersive medium model. Meshing divides the 3D model into multiple tetrahedrons. For example... Figure 2 As shown, a three-dimensional dispersive medium model of one embodiment is meshed.

[0053] S4, solve for the numerical flux of the six components of the three-dimensional electromagnetic field.

[0054] S5, update the auxiliary field variables of the three-dimensional CFS-PML region.

[0055] S6 updates the electromagnetic field components across the entire region.

[0056] S7, add a time step, and repeat steps S4-S6 until the 3D forward simulation of the current time step is completed. That is, each time step must undergo a round of steps S4-S6 until all time steps have undergone a round of steps S4-S6.

[0057] S8. Repeat steps S4-S7 until all excitations are completed, meaning all excitation sources at all excitation points have been activated. This means each excitation point needs to undergo one round of steps S4-S7, until all excitation points have undergone one round of steps S4-S7.

[0058] S9 outputs data to form a radar profile of a three-dimensional dispersive medium. For example... Figure 3 The image shown is a radar cross-section obtained in one embodiment.

[0059] This invention effectively solves the core problems of imperfect coupling of absorbing boundaries and limited multi-scale modeling accuracy in 3D ground-penetrating radar (GPR) forward modeling of dispersive media by constructing a 3D dispersive medium forward model and integrating CFS-PML boundary processing. By precisely setting CFS-PML parameters and updating the corresponding auxiliary field variables, non-physical reflection interference is significantly reduced, ensuring the efficiency and stability of boundary processing. The calculation logic for the evolution of the electromagnetic field across the entire 3D region is optimized through step-by-step iterative updates of electromagnetic field components and numerical flux, improving simulation accuracy and adapting to the multi-scale modeling needs of large-scale complex underground structures. The final output 3D radar profile can present the depth, morphology, and spatial distribution of underground targets in a three-dimensional manner, reducing detection blind spots and providing reliable support for geoelectric model construction, optimization of time-shift data processing algorithms, and breakthroughs in inversion imaging technology. This significantly enhances the application value and decision-making accuracy of GPR in practical exploration scenarios. This invention addresses the core issues of imperfect coupling of absorbing boundaries and limited accuracy of multi-scale modeling in forward modeling of dispersive media using 3D ground-penetrating radar. By clarifying the constituent parameters and calculation methods of the relative permittivity of the dispersive medium, it accurately quantifies the frequency dependence characteristics of the medium. Combined with the systematic solution of the numerical flux of the six components of the 3D electromagnetic field and the accurate calculation of the jump value of the normal component, it provides reliable data support for the evolution of the electromagnetic field and effectively avoids waveform distortion and target positioning deviation.

[0060] In some embodiments of the present invention, the dispersive medium physical property parameter is the relative permittivity ε. r ,

[0061] Where ε ∞ τ is the relative permittivity at infinite frequency, Δε is the difference between the relative permittivity at static zero frequency and at infinite frequency, τ is the relaxation time, j is the imaginary unit, and ω is the angular frequency.

[0062] The physical property parameter of the dispersive medium is clearly defined as the relative permittivity, which is composed of the relative permittivity at infinite frequency, the difference between the relative permittivity at static zero frequency and infinite frequency, and the relaxation time. This accurately quantifies the frequency dependence characteristics of the dispersive medium, providing a physical property basis that fits the actual medium for subsequent electromagnetic field evolution calculations. It ensures the accuracy of electromagnetic wave propagation simulation in the dispersive medium and avoids simulation deviations caused by ambiguity in the definition of physical property parameters.

[0063] In some embodiments of the present invention, the numerical flux in step S4 is calculated using the following formula:

[0064]

[0065] fluxH x fluxH y fluxH zThese represent the numerical fluxes of the magnetic field in the x, y, and z directions, respectively.

[0066] fluxE x fluxE y fluxE z These are the numerical fluxes of the electric field in the x, y, and z directions, respectively.

[0067] κ e v h , κ h v e For different numerical flux coefficients;

[0068] n x n y and n z It is the projection of the external normal vector of the tetrahedron onto the x, y, z directions;

[0069] H x H y H z These are the component values ​​of the magnetic field of the current unit in the x, y, and z directions, respectively;

[0070] E x E y E z These are the components of the electric field of the current cell in the x, y, and z directions, respectively.

[0071] These are the components of the magnetic field in the x, y, and z directions of the adjacent units, respectively.

[0072] These are the components of the electric field in the x, y, and z directions of the adjacent cells, respectively.

[0073] HH and EE represent the jump values ​​of the normal components of the magnetic field and electric field, respectively. The numerical flux of the six components of the three-dimensional electromagnetic field is solved using specific formulas. By combining key parameters such as the numerical flux coefficient, the external normal vector projection of the tetrahedron, and the electromagnetic field component values ​​of adjacent elements, and by introducing the jump values ​​of the normal components of the magnetic and electric fields, accurate calculation of the electromagnetic field component flux is achieved. This provides a reliable numerical basis for the evolution and updating of the electromagnetic field, improves the accuracy of three-dimensional electromagnetic field propagation simulation, and meets the characterization requirements of electromagnetic wave propagation characteristics in complex dispersive media.

[0074] In some embodiments of the present invention, the jump values ​​of the normal components of the magnetic field and electric field can be calculated using the following formula:

[0075]

[0076] The specific calculation formulas for the jump values ​​of the normal components of the magnetic field and electric field are given, which provide key parameter support for solving the numerical flux in claim 3, ensure the integrity and accuracy of the numerical flux calculation, avoid electromagnetic field propagation simulation errors caused by deviations in the calculation of the jump values ​​of the normal components, and further improve the numerical calculation reliability of the entire forward modeling method.

[0077] In some embodiments of the present invention, the auxiliary field variables include an electric field auxiliary variable P and a magnetic field auxiliary variable Q, wherein the electric field auxiliary variable P includes P x P y P z The magnetic field auxiliary variable Q includes Q0 x Q y Q z ;

[0078] P x =ε ∞ [(σ y -σ x -α z )P x,1 -(α y +σ x )P x,2 -(α x +σ x )P x,3 ]+P x,4 +P x,5 +P x,6 ;

[0079] Q x =μ[(σ y -σ x -α z )Q x,1 -(α y +σ x )Q x,2 -(α x +σ x )Q x,3 ] / ε0;

[0080] P y =ε ∞ [(σ z -σ y -α x )P y,1 -(α z +σ y )P y,2 -(α y +σ y )P y,3 ]+P y,4 +P y,5 +P y,6 ;

[0081] Q y =μ[(σ z -σ y -α x )Q y,1 -(α z +σ y )Q y,2 -(α y +σ y )Q y,3 ] / ε0;

[0082] P z =ε ∞ [(σ x -σ z -α y )P z,1 -(α x +σ z )P z,2 -(α z +σ z )P z,3 ]+P z,4 +P z,5 +P z,6 ;

[0083] Q z =μ[(σ x -σ z -α y )Q z,1 -(α x +σ z )Q z,2 -(α z +σ z )Q z,3 ] / ε0.

[0084] P x and Q x These are the electric field components E x and magnetic field component H x Auxiliary field variables;

[0085] P y and Q y These are the electric field components E y and magnetic field component H y Auxiliary field variables;

[0086] P z and Q z These are the electric field components E z and magnetic field component H z Auxiliary field variables.

[0087] ε ∞σ is the relative permittivity at infinite frequency; x σ y σ z These represent the conductivity in the x, y, and z directions within the CFS-PML layer, respectively; α x α y and α z ε0 represents the complex frequency shift parameters in the x, y, and z directions within the CFS-PML layer, respectively; μ is the permeability; ε0 is the permittivity in vacuum; P0 is the dielectric constant in vacuum. x,1 P x,2 P x,3 P x,4 P x,5 P x,6 These are the electric field components E x The six auxiliary field variable components, Q x,1 Q x,2 Q x,3 These are the magnetic field components H. x The three auxiliary field variable components.

[0088] P y,1 P y,2 P y,3 P y,4 P y,5 P y,6 These are the electric field components E y The six auxiliary field variable components, Q y,1 Q y,2 Q y,3 These are the magnetic field components H. y The three auxiliary field variable components.

[0089] P z,1 P z,2 P z,3 P z,4 P z,5 P z,6 These are the electric field components E z The six auxiliary field variable components, Q z,1 Q z,2 Q z,3 These are the magnetic field components H. z The three auxiliary field variable components. σ is the conductivity.

[0090] Because it is a dispersive medium, the relative permittivity ε r It is frequency dependent, so there are 6 auxiliary variables for each electric field component and 3 auxiliary variables for each magnetic field component.

[0091] The auxiliary variables can be updated according to the following formula:

[0092]

[0093]

[0094] Integrating the above formula and substituting the current value into the parameter on the right-hand side yields the updated value. (See formula...) In the middle, P on the right side of the equal sign x,1 and E x All values ​​are currently known. After substituting them, the expression is integrated to obtain the updated P. x,1 .

[0095] To address the frequency dependence of the relative permittivity of dispersive media, the settings and related parameter associations of the electric field auxiliary variables P (6 for each electric field component) and magnetic field auxiliary variables Q (3 for each magnetic field component) in the CFS-PML region were clarified. This improved the coupling mechanism between CFS-PML and dispersive media, ensured the stability of boundary treatment in dispersive media environments, and reduced the interference of non-physical reflections on forward modeling results.

[0096] In some embodiments of the present invention, step S6 updates the electromagnetic field components of the entire region, which mainly refers to updating the electromagnetic field components of the simulation region and the CFS-PML boundary. The specific update equation is as follows:

[0097]

[0098]

[0099] Integrating the above formula and substituting the current value into the parameter on the right-hand side yields the updated value. (See formula...) In the middle, H on the right side of the equal sign x Substitute the current known value into the equation, integrate the expression, and obtain the updated P. x,1 .

[0100] Let D be the partial derivative, t be time, and D be the partial derivative. k For the current tetrahedron, M is the mass matrix of the element; S x S y and S z , respectively, are the unit rigidity matrices in the x, y, and z directions; F is the boundary element matrix; J is the integral term of the current source. z The loading stimulus source; Let be the vector composed of the values ​​of the basis functions at this element node; the superscript T indicates the transpose of the vector; dv represents the volume element. The update equations for the electromagnetic field components across the entire region (including the simulation region and the CFS-PML boundary) are clearly defined. By introducing parameters such as the element mass matrix, unit stiffness matrix, boundary element matrix, and current source integral term, a systematic and accurate update of the electromagnetic field components is achieved. This ensures the continuity and accuracy of the electromagnetic field evolution process throughout the entire computational domain, taking into account both the propagation characteristics of the simulation region and the absorption characteristics of the boundary region, thus improving the overall reliability of the forward modeling simulation.

[0101] In some embodiments of the present invention, the CFS-PML parameters at the boundary include the maximum value α of the complex frequency shift parameter at the outer boundary of the CFS-PML. max The maximum conductivity σ at the outer boundary of CFS-PML max The PML thickness d0 can be obtained using the following formula: σ i and α i :

[0102]

[0103] σ i The conductivity in the i-direction within the CFS-PML layer is σ. x σ y σ z ;α i Let α be the complex frequency shift parameter in the i-direction within the CFS-PML layer. x α y and α z ;α max The maximum value of the complex frequency shift parameter at the outer boundary of CFS-PML; d i It is the distance to the innermost layer i of the PML in the direction; d0 is the thickness of the PML.

[0104] σ max is the maximum conductivity at the outer boundary of CFS-PML, where m is the exponential order.

[0105] The specific formulas for establishing the three-dimensional dispersive geological forward model and setting CFS-PML parameters are refined. By accurately defining the calculation methods of key parameters such as conductivity and complex frequency shift parameters, especially the method for determining the maximum value of the complex frequency shift parameter, the absorption performance of CFS-PML for evanescent waves and low-frequency waves is significantly optimized, further reducing non-physical reflections at the boundary, ensuring the accuracy of broadband electromagnetic wave propagation simulation, and adapting to the boundary processing requirements of three-dimensional dispersive medium forward modeling.

[0106] In some embodiments of the present invention, the maximum conductivity σ at the outer boundary of the CFS-PML is... max It is obtained by calculation using the following formula:

[0107] σ max =-(m+1)ln(R0) / (2ηd0ε r ).

[0108] Where m is the exponential order; R0 is the expected reflection error; η is the wave impedance of the PML layer; ε r Let be the relative permittivity of the medium. A specific formula for calculating the maximum conductivity at the outer boundary of CFS-PML is provided. Combined with parameters such as exponential order, expected reflection error, and wave impedance, a quantitative design for the maximum conductivity is achieved. This ensures that the absorption effect at the CFS-PML boundary meets the expected standard, avoiding insufficient absorption performance due to unreasonable conductivity parameter settings. This further improves the scientific rigor and effectiveness of boundary treatment, providing boundary parameter assurance for high-precision forward modeling.

[0109] This invention also provides a three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling system, comprising: a modeling and parameter setting module for establishing a forward model of the three-dimensional dispersive medium and setting CFS-PML parameters at the boundaries; an input setting module for inputting dispersive medium physical property parameters and setting the excitation and receiving point positions; a parameter calculation module for meshing the three-dimensional dispersive medium model; a numerical flux solving module for solving the numerical flux of the six components of the three-dimensional electromagnetic field; an auxiliary field variable update module for updating the auxiliary field variables of the three-dimensional CFS-PML region; an electromagnetic field component update module for updating the electromagnetic field components of the entire region; a cyclic simulation module for increasing the time step and completing the three-dimensional forward simulation of the current time step; a cyclic excitation module for completing all excitations; and an output module for outputting data to form a radar profile of the three-dimensional dispersive medium. The forward modeling method is broken down into multiple functional modules such as modeling and parameter setting, input setting, and parameter calculation. This realizes modular division of labor and collaborative operation in the forward modeling process of 3D ground-penetrating radar dispersive media, improves the maintainability and scalability of the forward modeling system, and ensures the orderliness and efficiency of the forward modeling process by each module performing specific functions. This facilitates subsequent system optimization and function upgrades and adapts to the forward modeling simulation needs of different scenarios.

[0110] This invention also provides a storage medium, a computer-readable storage medium storing the computer program, which, when executed by the processor, implements the three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method. Storing the three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method in the form of a computer program on a readable storage medium achieves the reusability and portability of the forward modeling method, facilitating its deployment and execution on different computing devices, lowering the threshold for practical application, and providing convenient conditions for rapidly conducting forward simulations in actual exploration scenarios, thus promoting and applying this technology in engineering practice.

[0111] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method, characterized in that, Includes the following steps: S1, Establish a forward model of the three-dimensional dispersive medium and set the CFS-PML parameters at the boundary; S2, input the physical properties of the dispersive medium, and set the positions of the excitation point and the receiver point; S3, meshing the three-dimensional dispersive medium model; S4, solve for the numerical flux of the six components of the three-dimensional electromagnetic field; S5, update the auxiliary field variables of the three-dimensional CFS-PML region; S6 updates the electromagnetic field components across the entire region; S7, add a time step, repeat steps S4-S6 until the three-dimensional forward simulation of the current time step is completed; S8. Repeat steps S4-S7 until all excitations are complete. S9 outputs data to form a radar profile of a three-dimensional dispersive medium.

2. The three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method according to claim 1, characterized in that, The physical property parameter of the dispersive medium is the relative permittivity ε. r , Where ε ∞ τ is the relative permittivity at infinite frequency, Δε is the difference between the relative permittivity at static zero frequency and at infinite frequency, τ is the relaxation time, j is the imaginary unit, and ω is the angular frequency.

3. The three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method according to claim 1, characterized in that, The numerical flux in step S4 is calculated using the following formula: fluxH x fluxH y fluxH z These represent the numerical fluxes of the magnetic field in the x, y, and z directions, respectively. fluxE x fluxE y fluxE z These are the numerical fluxes of the electric field in the x, y, and z directions, respectively. κ e v h , κ h v e For numerical flux coefficients; n x n y and n z It is the projection of the external normal vector of the tetrahedron onto the x, y, z directions; H x H y H z These are the component values ​​of the magnetic field of the current unit in the x, y, and z directions, respectively; E x E y E z These are the components of the electric field of the current cell in the x, y, and z directions, respectively. These are the components of the magnetic field in the x, y, and z directions of the adjacent units, respectively. These are the components of the electric field in the x, y, and z directions of the adjacent cells, respectively. HH and EE represent the jump values ​​of the normal components of the magnetic field and electric field, respectively.

4. The three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method according to claim 3, characterized in that, The jump values ​​of the normal components of the magnetic and electric fields can be calculated using the following formula:

5. The three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method according to claim 1, characterized in that, The auxiliary field variables include an electric field auxiliary variable P and a magnetic field auxiliary variable Q. The electric field auxiliary variable P includes P0. x The magnetic field auxiliary variable Q includes Q0 x The P x and Q x Update using the following formula: P x =e ∞ [(s y -s x -a z )P x,1 -(a y +s x )P x,2 -(a x +s x )P x,3 ]+P x,4 +P x,5 +P x,6 ; Q x =μ[(σ y -s x -a z )Q x,1 -(a y +s x )Q x,2 -(a x +s x )Q x,3 ] / ε0; Among them, P x and Q x These are the electric field components E x and magnetic field component H x Auxiliary field variables; ε ∞ σ is the relative permittivity at infinite frequency; x and σ y These represent the conductivity in the x and y directions within the CFS-PML layer, respectively; α x α y and α z ε0 represents the complex frequency shift parameters in the x, y, and z directions within the CFS-PML layer, respectively; μ is the permeability; ε0 is the permittivity in vacuum; P0 is the dielectric constant in vacuum. x,1 P x,2 P x,3 P x,4 P x,5 P x,6 These are the electric field components E x The six auxiliary field variable components, Q x,1 Q x,2 Q x,3 These are the magnetic field components H. x The three auxiliary field variable components.

6. The three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method according to claim 1, characterized in that, The electromagnetic field components of the entire region are updated according to the following formula: in, Let D be the partial derivative, t be time, and D be the partial derivative. k For the current tetrahedron, M is the mass matrix of the element; S x S y and S z , respectively, are the unit rigidity matrices in the x, y, and z directions; F is the boundary element matrix; J is the integral term of the current source; z The loading stimulus source; dv is the vector composed of the values ​​of the basis functions at this element node; the superscript T indicates the transpose of the vector; dv represents the volume element; Δε is the difference in relative permittivity between the static zero frequency and the infinite frequency; σ is the conductivity.

7. The three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method according to claim 1, characterized in that, The CFS-PML parameters at the boundary include the maximum value α of the complex frequency shift parameter at the outer boundary of the CFS-PML. max The maximum conductivity σ at the outer boundary of CFS-PML max The PML thickness d0 can be obtained using the following formula: σ i and α i : σ i α represents the conductivity along the i-direction within the CFS-PML layer. i α is the complex frequency shift parameter in the i-direction within the CFS-PML layer; max The maximum value of the complex frequency shift parameter at the outer boundary of CFS-PML; d i It is the distance to the innermost layer i of the PML in the direction; d0 is the thickness of the PML, σ max is the maximum conductivity at the outer boundary of CFS-PML, where m is the exponential order.

8. The three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method according to claim 7, characterized in that, The maximum conductivity σ at the outer boundary of the CFS-PML max σ is obtained by calculating using the following formula: max =-(m+1)ln(R0) / (2ηd0ε r ) Where m is the exponential order; R0 is the expected reflection error; η is the wave impedance of the PML layer; ε r is the relative permittivity of the medium.

9. A three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling system, characterized in that, include: The modeling and parameter setting module is used to establish a forward model of a three-dimensional dispersive medium and set the CFS-PML parameters at the boundary. The input settings module is used to input the physical properties of the dispersive medium and set the positions of the excitation and receiving points; The parameter calculation module is used to mesh the three-dimensional dispersive medium model; The numerical flux solver module is used to solve for the numerical flux of the six components of the three-dimensional electromagnetic field. The auxiliary field variable update module is used to update the auxiliary field variables of the three-dimensional CFS-PML region; Electromagnetic field component update module, used to update the electromagnetic field components of the entire area; The loop simulation module is used to increase the time step and complete the three-dimensional forward simulation of the current time step; The cyclic excitation module is used to complete all excitations; The output module is used to output data and form a radar profile of a three-dimensional dispersive medium.

10. A storage medium, wherein the computer-readable storage medium stores the computer program, characterized in that, When the computer program is executed by the processor, it implements the three-dimensional ground-penetrating radar dispersive medium discontinuous finite element forward modeling method as described in any one of claims 1 to 8.

Images