An adaptive finite element method for 3D electromagnetic potential forward modeling
By transforming Maxwell's equations into governing equations for scalar and vector potentials, and utilizing the normal current density discontinuity error estimation and an octet refinement strategy, the non-uniqueness and convergence problems of finite element numerical solutions in three-dimensional electromagnetic exploration are solved, achieving efficient mesh adaptive optimization and accurate electromagnetic field calculation.
Patent Information
- Application Number
- CN202310155976.5
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2023-02-23
- Publication Date
- 2026-02-13
- Estimated Expiration
- 2043-02-23
AI Technical Summary
Existing technologies suffer from non-uniqueness and convergence issues in finite element numerical solutions in 3D electromagnetic exploration, especially in frequency domain geophysical EM forward modeling, resulting in high computational costs and insufficient accuracy.
The Maxwell equations are transformed into governing equations for scalar and vector potentials using the Coulomb gauge. Gradient terms of Lagrange multipliers are introduced, and the error generated by the discontinuity of normal current density is used as the posterior error estimate. An octet refinement strategy is adopted to refine the mesh, and an adaptive finite element method is constructed.
It improves the stability and accuracy of finite element numerical solutions, reduces the number of iterations, lowers computational costs, and achieves efficient mesh adaptive optimization.
Smart Images
Figure CN116227286B_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The present application relates to the field of forward technology of geophysical electromagnetic exploration method, and particularly relates to a self-adaptive finite element method for three-dimensional electromagnetic potential forward modeling. BACKGROUND
[0002] Developing a numerical forward modeling method that can accurately calculate the induced electromagnetic (EM) field and understand the propagation of the EM field on the ground is necessary for the inversion and interpretation of three-dimensional geophysical EM data. In order to solve the traditional EM Helmholtz equation, the vector finite element (FE) method is applied to solve the electric diffusion equation (Borner et al., 2008; Liu et al., 2008; Farquharson and Miensopust 2011; Ren et al. 2013) and the magnetic field (Franke et al. 2007) diffusion equation (varisuha, 2018). In order to obtain a high-precision numerical solution with as little computational cost as possible, the adaptive finite element forward modeling technology has been rapidly developed. However, when solving the frequency domain geophysical EM forward problem with finite elements, the solving efficiency of the finite element system is still limited by the non-uniqueness of the solution.
[0003] The divergence correction based on the non-dispersive magnetic flux or current density can reduce the non-uniqueness of the solution and accelerate the convergence of the finite element numerical solution. This divergence correction technique is very effective in improving the finite element numerical solution efficiency of passive and active electromagnetic systems (Mackie et al., 1994; Smith, 1996; Druskin et al., 1999 and LaBrecque, 1999). For the electric field diffusion equation, this divergence correction based on the current conservation in the region uses the electric scalar potential as the divergence correction term, and uses the electric scalar potential gradient to correct the numerical solution of the electric field diffusion equation. According to the role of the electric scalar potential term in the process of finite element numerical solution, the divergence correction can be divided into explicit and implicit types.
[0004] The explicit divergence correction scheme is inserted into the iterative solution process to correct the numerical solution (Schwarzbach et al. 2010). The correction interval determines the effect of improving the solution of the problem, which needs to be preset in the specific case of the electromagnetic system and there is no fixed standard (Liu and Yin, 2013; Zhou et al. 2021). The implicit divergence correction scheme is to establish a vector-scalar potential system with Coulomb gauge (Biro and Preis, 1989; Haber, 2000, 2001; Badea, 2001; Everett et al. 2001), which means that no additional divergence correction is needed. The electric scalar potential is still a correction term in the divergence correction process. The applied Coulomb gauge ensures that the electromagnetic potential does not contain divergent spurious modes, and further limits the non-uniqueness of the vector potential (Badea, 2001; Ansari, 2014).
[0005] In order to further enhance the integrity of the A-φ system (control equation) to improve the solving efficiency of the control equation, and solve the problem that the current divergence constraint for calculating the "posterior error" (the posterior error is used to determine the unit that needs to be refined, and finally optimize the grid) is lacking in the current "electromagnetic potential control equation", we develop a goal-oriented adaptive finite element scheme for the Coulomb gauge A-φ system. SUMMARY
[0006] The present application aims to provide an adaptive finite element method for three-dimensional electromagnetic potential forward modeling to solve the problems mentioned in the background art.
[0007] To achieve the above purpose, the present application provides the following technical scheme:
[0008] An adaptive finite element method for three-dimensional electromagnetic potential forward modeling, comprising the following steps:
[0009] S1, converting Maxwell's equations into control equations of scalar potential and vector potential by using Coulomb gauge;
[0010] S2, introducing a gradient term of Lagrange multiplier into the control equation, discretizing the solution region, and constructing a finite element linear equation system;
[0011] S3, solving the finite element linear equation system;
[0012] S4, using the error generated by the discontinuity of the normal current density as the posterior error estimator, dual-weighting the posterior error estimator and determining the target grid that needs to be refined;
[0013] S5, using an octant refinement strategy to encrypt the target grid;
[0014] S6, using the encrypted grid for finite element forward modeling calculation.
[0015] Preferably, the step S1 specifically comprises:
[0016] Take the time harmonic factor as e iωt Under quasi-static conditions, the propagation of electromagnetic fields satisfies Maxwell's equations:
[0017]
[0018]
[0019] Substitute equation (2) into equation (1) to obtain the electric field double spinor equation:
[0020]
[0021] The electromagnetic field is represented by vector-scalar potential:
[0022]
[0023]
[0024] Equation (3) is represented by electromagnetic potential:
[0025]
[0026] The divergence of the current density in the passive region is zero, For the active region The electromagnetic equation satisfies:
[0027]
[0028] The Coulomb gauge is represented as a separate equation as follows:
[0029]
[0030] Where E is the electric field intensity, H is the magnetic field intensity, σ is the conductivity of the underground dielectric, μ0 is the magnetic permeability of free space, ω is the angular frequency, and the current density J s Is the field source of the induced electromagnetic field.
[0031] Preferably, the step S2 specifically comprises:
[0032] The gradient λ term of the Lagrange multiplier introduced in the control equation is represented as:
[0033]
[0034] Using the vector shape function N and the node shape function N, the vector potential and the scalar potential are discretized as follows:
[0035]
[0036]
[0037] where n A and n φ are the number of edges and nodes of the element, respectively; the above equation is discretized as follows:
[0038]
[0039] where, is the approximate solution of A, λ, Ω is the domain of the calculation model, i and j are the edge numbers (i, j = 1, 2, 3,..., N edges ), and l and k are the node numbers (l, k = 1, 2, 3,..., N nodes );
[0040] The current divergence term is discretized as follows:
[0041]
[0042] The condition that the divergence of the vector potential in the Coulomb gauge is zero is discretized as follows:
[0043]
[0044] Combining equations (12), (13) and (14), the finite element linear equation set is expressed as:
[0045]
[0046]
[0047] where the term related to the excitation source in the right end of equation (13) and equation (14) is put into the vector S.
[0048] Preferably, in step S4, the weighted posterior error estimator of the i-th element is expressed as:
[0049] η i = η w,i · η e,i (17)
[0050] η e,i is the posterior error estimator, which is expressed as:
[0051]
[0052] where the k-th element face F k of the i-th element is considered with the normal n error, and ε k+ and εk- respectively, are the normal errors inside and outside the element, and the error size is L2 norm;
[0053] η w,i is the weighted coefficient, and is expressed as:
[0054]
[0055] where w + and w - are the influence functions inside and outside the element, respectively, and are obtained by solving the dual weak form of the electromagnetic potential equation.
[0056] Preferably, in the step S5, the octahedral refinement strategy is: taking the midpoint of each edge of the old tetrahedron as the base point, connecting the adjacent base points to form a cross section, and dividing the old tetrahedron into eight new tetrahedrons.
[0057] The beneficial effects of the present application are:
[0058] (1) The present scheme implicitly integrates the divergence correction into the finite element system equation, and uses the Lagrange compensation term to enhance the continuity of the vector potential and the stability in the solving process;
[0059] (2) The present scheme determines the error generated by the discontinuity of the normal current density as the posteriori error estimator; compared with the other three kinds of errors, when the degrees of freedom are similar, the discontinuity of the normal current density can obtain the highest precision numerical solution at the observation point;
[0060] (3) The present scheme uses the octahedral refinement strategy to encrypt the target grid, which can reduce the number of iterations required by the equal division refinement, and the grid self-adapting process does not show obvious instability, so that the grid can reach the optimal state earlier, that is, the calculation efficiency can be improved by reducing the calculation amount of solving the additional system equation caused by the grid self-adapting. BRIEF DESCRIPTION OF DRAWINGS
[0061] Figure 1 A flow chart of a three-dimensional electromagnetic potential forward modeling self-adapting finite element method provided by the embodiment of the present application;
[0062] Figure 2 A tetrahedral octahedral refinement schematic diagram;
[0063] Figure 3 A three-layer layered model schematic diagram based on MT forward modeling;
[0064] Figure 4 MT self-adapting forward response error reduction processes based on different posteriori error estimators in the three-layer layered model; (a) and (b) respectively represent the error reduction processes of apparent resistivity and phase;
[0065] Figure 5 Fig. 19 is a schematic diagram of a 3D electrical heterogeneous model for MT;
[0066] Figure 6 Fig. 20 is a diagram of MT adaptive forward response error reduction curves of four kinds of posterior error estimates in the 3D electrical heterogeneous model; (a) and (b) represent error reduction processes of apparent resistivity and phase, respectively;
[0067] Figure 7 Fig. 21 is a diagram of mesh adaptation results of different posterior error estimates of the 3D electrical heterogeneous model for MT; (a)-(d) represent meshes corresponding to Jn, Bn, Bt and residual, respectively, and subsequent labels 1, 2 and 3 represent mesh refinement results of the 1st, 10th and 20th iterations, respectively; white boxes represent abnormal blocks;
[0068] Figure 8 Fig. 22 is a diagram of mesh adaptation effects of eight-subdivision and two-subdivision in the 3D electrical heterogeneous model for MT; (a) and (b) represent decreasing processes of response error with mesh iteration when the mesh growth rate is 1%, (c) and (d) represent decreasing processes when the mesh growth rate is 3%, (a) and (c) represent apparent resistivity error, and (b) and (d) represent phase error;
[0069] Figure 9 Fig. 23 is a diagram of a 3D electrical heterogeneous model with a line source;
[0070] Figure 10 Fig. 24 is a diagram of adaptive forward response error of different posterior error estimates in the line source 3D electrical heterogeneous model; (a) and (b) represent error reduction processes of apparent resistivity and phase, respectively;
[0071] Figure 11 Fig. 25 is a diagram of mesh adaptation results of different posterior error estimates of the line source 3D electrical heterogeneous model; (a)-(d) represent meshes corresponding to Jn, Bn, Bt and residual, respectively, and subsequent labels 1, 2 and 3 represent mesh refinement results of the 1st, 10th and 20th iterations, respectively; white boxes represent abnormal blocks. DETAILED DESCRIPTION
[0072] The application will be further described in detail below in combination with the drawings and embodiments:
[0073] As shown in the drawings, an adaptive finite element method for 3D electromagnetic potential forward modeling comprises the following steps: Figure 1 S1, converting Maxwell equations into control equations of scalar potential and vector potential by using Coulomb gauge;
[0074] Specifically, it comprises:
[0075]
[0076] Taking the time harmonic factor as e iωt Under quasi-static conditions, the propagation of electromagnetic fields satisfies Maxwell's equations:
[0077]
[0078]
[0079] Substituting equation (2) into equation (1), the electric field double-rotor equation is obtained:
[0080]
[0081] The electromagnetic field is expressed by vector-scalar potential as:
[0082]
[0083]
[0084] Equation (3) is expressed by electromagnetic potential as:
[0085]
[0086] The divergence of the current density in the passive region is zero, For the active region The electromagnetic equation satisfies:
[0087]
[0088] The Coulomb gauge is expressed as a separate equation as follows:
[0089]
[0090] Where E is the electric field intensity, H is the magnetic field intensity, σ is the conductivity of the dielectric, μ0 is the magnetic permeability of free space, ω is the angular frequency, and the current density J s Is the field source of the induced electromagnetic field.
[0091] S2, the gradient term of the Lagrange multiplier introduced in the control equation, the solving area is discretized, and the finite element linear equation group is constructed;
[0092] Specifically, it includes:
[0093] In order to enhance the continuity of the numerical solution of A in the unit body and enhance its stability in the solving process, the gradient λ term of the Lagrange multiplier is introduced in the control equation, which is expressed as:
[0094]
[0095] Using the vector shape function N and the node shape function N, the vector potential and the scalar potential are discretized as follows:
[0096]
[0097]
[0098] where n A and n φ are the number of edges and nodes of the unit cell, respectively; the above equation is discretized as follows:
[0099]
[0100] where, is the approximate solution of A, λ, Ω is the domain of the calculation model, i and j are the edge numbers (i, j = 1, 2, 3,..., N edges ), and l and k are the node numbers (l, k = 1, 2, 3,..., N nodes );
[0101] The current divergence term is discretized as follows:
[0102]
[0103] The condition that the divergence of the vector potential in the Coulomb gauge is zero is discretized as follows:
[0104]
[0105] In combination with equations (12), (13) and (14), the finite element linear equation set is represented as:
[0106]
[0107]
[0108] where the term related to the excitation source in the right end of equation (13) and equation (14) is put into the vector S.
[0109] S3, solving the finite element linear equation set;
[0110] Specifically, based on the initial three-dimensional tetrahedral mesh partitioning, the finite element linear equation set is solved.
[0111] S4, using the error generated by the discontinuity of the normal current density as the posterior error estimator, dual-weighted posterior error estimator and determining the target mesh that needs to be refined;
[0112] Specifically, the weighted posterior error estimator of the i-th unit is represented as:
[0113] η i = η w,i · η e,i (17)
[0114] η e,i is the posterior error estimator, for a tetrahedron, the error on its four faces is added, when measuring the tangential discontinuity of a physical quantity, denoted as:
[0115]
[0116] where F k is the normal n error of the reference face, ε k+ and ε k- are the normal error inside and outside the element, the error size is the L2 norm;
[0117] η w,i is the weighting coefficient, denoted as:
[0118]
[0119] where w + and w - are the influence functions inside and outside the element, which are obtained by solving the dual weak form of the electromagnetic potential equation, that is:
[0120] Equation (15) can be expressed as where Φ represents the vector basis function, and L is a square integrable function, which can be expressed as:
[0121]
[0122] By solving the dual weak form of the electromagnetic potential equation, the influence function can be obtained, and its dual form is as follows:
[0123] B * (w,Φ)=L(Φ) (21)
[0124] where B*(,) is the dual form of B(,) of the electromagnetic potential equation set, and due to its characteristics, B*(,) is equal to B(,).
[0125] S5, using an octree refinement strategy to encrypt the target mesh;
[0126] Specifically, in order to quickly refine to meet the accuracy requirements of numerical solution, the midpoint of each edge of the old tetrahedron is taken as the base point, the adjacent base points are connected to form a cross section, and the old tetrahedron is divided into eight new tetrahedrons. Moreover, in order to avoid that some tetrahedrons are too small in the adaptive refinement process, a height threshold h t Set the minimum cell volume threshold of the regular tetrahedron Considering that the minimum cell size affects the accurate calculation of electromagnetic wave propagation, it is necessary to set a reference skin depth δ, and the cell height threshold is δ / 8.
[0127] S6. Use the encrypted mesh to perform finite element forward modeling calculations.
[0128] Specifically, the mesh obtained after refinement in step S5 is recorded as the initial mesh, and steps S3-S5 are repeated until a finite element solution that satisfies the given error condition is obtained.
[0129] Numerical Examples Based on MT Forward Modeling
[0130] Numerical Example 1
[0131] Establish a three-layered model ( Figure 3 To test the influence of different posterior error estimates on the grid adaptation process of the magnetotelluric system, the resistivity of the three layers was 100 Ω·m, 10 Ω·m, and 100 Ω·m, respectively. The thickness of the first and second layers was 4 km, and the last layer was a uniform space extending downwards. The size of the finite element discrete model was 60 km × 60 km × 1700 km, and the air layer thickness was set to 20 km. We calculated the magnetotelluric (MT) response at 1 Hz. During the grid adaptive optimization process, based on the overall posterior error, the top 1% of error values were marked, and the corresponding element cells were bisected and refined. Four posterior error estimates were selected: normal current discontinuity Jn, residual Res, normal magnetic flux density Bn, and tangential magnetic flux density Bt. The error reduction process of the MT adaptive forward response for the four posterior error estimates was obtained. Figure 4 With each mesh update, the degrees of freedom (DoF) of the system equations will change accordingly. The faster the error decreases and the smaller the mismatch value, the more effective the adaptive optimization of the mesh.
[0132] In the mesh adaptive optimization process using the magnetic flux density (B) discontinuity as the posterior error estimate, whether in the normal or tangential direction, the MT response error at the observation points did not show a significant decreasing trend. However, the errors corresponding to the residuals and normal current discontinuities decreased with increasing iterations. Among these, the error curve for Jn showed the best decreasing trend. In the results for apparent resistivity and phase, the error curve for Jn decreased faster than the residual curve. When the number of iterations is the same, the normal current density discontinuity can effectively refine more elements, resulting in more degrees of freedom in the corresponding finite element system equations and higher MT response accuracy.
[0133] Numerical Example 2
[0134] Forward modeling of MT based on a more complex three-dimensional electrical inhomogeneity model is used for optimal posterior error estimation in test grid adaptation. In the three-dimensional electrical inhomogeneity model, two 10 Ω·m anomalous blocks are buried in a 100 Ω·m homogeneous half-space medium underground. Figure 5). The top depths of the two anomalous blocks are 4 km and 6 km, respectively, and the horizontal distance between them is 2 km. The size of the two anomalous blocks is 2 km x 2 km x 2 km. The measurement line is placed on the ground surface directly above the blocks.
[0135] The cells with error size in the top 1% are halved, and the grid adaptation iteration response error curves corresponding to the four posterior error estimates are shown in Figure 6 . The normal current discontinuity can be used as an error estimate to achieve the fastest response error reduction rate. The response error reduction trend of the vector potential residual is not as good as that of the tangential magnetic induction intensity discontinuity. In addition, the residual as an error indicator is difficult to accurately indicate the area that has a greater impact on the accuracy of the numerical solution, which leads to a limited number of correctly refined cells after each optimization iteration. As an error estimate, the tangential magnetic induction intensity discontinuity can maintain the same error curve decline trend as the normal current discontinuity in the early stages of the optimization process, but it cannot maintain this trend in the middle and later stages of the optimization iteration.
[0136] The impact of different posterior error estimates on grid adaptation optimization can be confirmed by the grid changes with iterations Figure 7 . The Jn discontinuity is used as a posterior error estimate to correctly indicate the cells to be refined in the near-surface and anomalous block areas. The posterior error estimate of the B discontinuity can only indicate the cells to be refined in the near-surface area, but it is too concentrated. They also incorrectly mark part of the air. In the residual grid adaptation optimization process, the grid refinement degree is not high for both the near-surface area and the anomalous block area. The grid change results of the residual show that it is difficult to accurately identify the factors affecting the accuracy of MT response data as a posterior error estimate.
[0137] The grid adaptation effects of octree refinement and halving refinement are shown by the response error changes after optimization iterations Figure 8 . When the dof is similar, the two refinement strategies can achieve similar MT response accuracy. According to the size of the posterior error, whether the top 1% or the top 3% of the cells are marked for refinement, the error curve of octree refinement and halving refinement has basically the same decline trend.
[0138] Numerical example based on line source forward modeling
[0139] Numerical example one
[0140] A three-dimensional block anomaly model with a line source is used Figure 9The effectiveness of the adaptive finite element scheme in solving the measured electromagnetic potential problem is tested. There is a 120 m x 200 m x 400 m block in a homogeneous subsurface space with a resistivity of 50 Ω·m, and a 100-meter-deep block is buried on the top of the block. A 100-meter-long grounding wire is laid on the ground at a horizontal distance of 900 meters from the block, and a 1A current with a frequency of 3Hz is transmitted.
[0141] For the forward modeling of the target-oriented adaptive finite element A-φ system with a line source, the influence of the above four error estimates in the grid adaptation process is evaluated by comparing the error curves of the response data. Figure 10 The error curves of the four kinds of posterior error estimates all show a downward trend, indicating that they can effectively indicate the factors that have a greater impact on the accuracy of the 3Hz line source forward simulation response data.
[0142] Among the four error curves, the error curve of Bn has no obvious downward trend, and the error curve of Bt only has an obvious downward trend in the early iterations. The residual curve has a significant decreasing trend, and the change trend of the real part of the electric field is more obvious than that of the imaginary part of the electric field. However, this trend is not as good as that of Jn. As the grid optimization process proceeds, Jn gradually becomes the most prominent, and its error reduction rate is significantly faster than that of the other three error indicators. Specifically, Jn has more degrees of freedom in one iteration, or when the degrees of freedom are similar, its error is significantly smaller.
[0143] The grid changes corresponding to the four posterior error estimates in the grid adaptation process are shown in Figure 11 Based on the error estimates of Bt or Bn discontinuity, part of the air can be easily marked as an area that needs to be refined. Their grid refinement is carried out around the line source and the upper and lower areas of the ground surface, and their indication of the abnormal block area is limited. Whether tangential or normal, the error estimates of the magnetic induction intensity B are difficult to clearly indicate those factors that reduce the accuracy of the response data. The residual-based error estimate cannot effectively mark the area near the line source as an area to be refined, and its grid refinement effect in the near-surface area is not ideal. The grid refinement of Jn is reflected on the current line source, the near-surface area and the abnormal block area. Moreover, the grid refinement of Jn is not obviously concentrated in a certain area.
[0144] The octree refinement and the binary refinement can obtain similar numerical solution accuracy, which indicates that the advantage of the octree refinement is that it can refine the elements to the desired state with fewer iteration times. Considering the grid adaptation optimization iteration as an optimization problem that generates elements that minimize the numerical solution error of the measuring point response data, the octree refinement and the binary refinement can be regarded as the difference in the step size of finding the extreme value in the optimization process, and the octree refinement strategy is a larger step optimization scheme.
[0145] The above-mentioned are only embodiments of the present application, and common technical solutions and / or common knowledge of the scheme are not described in detail. It should be pointed out that, for those skilled in the art, without departing from the technical solutions of the present application, a number of modifications and improvements can be made, which should also be considered as the protection scope of the present application, and these will not affect the effect and practicality of the patent. The protection scope claimed in the present application should be subject to the content of its claims, and the specific implementation mode and the like recorded in the specification can be used to explain the content of the claims.
Claims
1. An adaptive finite element method for three-dimensional electromagnetic potential forward modeling, characterized in that, Includes the following steps: S1. Use the Coulomb gauge to transform Maxwell's equations into governing equations for scalar and vector bits; S2. The gradient term of the Lagrange multiplier is introduced into the governing equation to discretize the solution domain and construct a system of linear equations for the finite element method. S3. Solve the aforementioned system of linear finite element equations; S4. Use the error caused by the discontinuity of normal current density as the posterior error estimate, dual-weighted posterior error estimate and determine the target mesh that needs to be refined. In step S4, the weighted posterior error estimate of the i-th grid cell is expressed as: ， η e,i It is the posterior error estimator, expressed as: ， Wherein, the k-th element surface Fk of the i-th element is the reference surface considering the normal error n, ε k+ and ε k- These are the normal errors inside and outside the element, respectively, and the magnitude of the error is the L2 norm. η w,i These are weighting coefficients, expressed as: ， w + and w - These are the influence functions inside and outside the unit, respectively, obtained by solving the dual weak form of the electromagnetic potential equation; S5. The target mesh is refined using an octet refinement strategy. In step S5, the octet refinement strategy is as follows: using the midpoints of each edge of the old tetrahedron as base points, adjacent base points are connected to form a cross section, dividing the old tetrahedron into eight new tetrahedrons. To avoid some tetrahedrons becoming too small during adaptive refinement, a height threshold (ht) is used to set the minimum cell volume threshold for the regular tetrahedron. Considering that the minimum cell size affects the accurate calculation of electromagnetic wave propagation, it is necessary to set it with reference to the skin depth δ, and the cell height threshold is δ / 8; S6. Use the encrypted mesh to perform finite element forward modeling calculations.
2. The adaptive finite element method for three-dimensional electromagnetic potential forward modeling according to claim 1, characterized in that, Step S1 specifically includes: Take the time harmonic variation factor as Under quasi-static conditions, the propagation of the electromagnetic field satisfies Maxwell's equations: ， Substituting equation (2) into equation (1), we obtain the double curl equation of the electric field: ， Electromagnetic fields can be represented by a vector-scalar potential as follows: ， Equation (3) can be expressed in terms of electromagnetic potential as follows: ， The current density divergence is zero in the passive region, ▽·J=0; for the active region, ▽·J=▽·J s Electromagnetic type satisfies: ， The Coulomb gauge is represented as a single equation as follows: ， Where E is the electric field strength, H is the magnetic field strength, and σ is the conductivity of the underground dielectric. ω is the permeability of free space, ω is the angular frequency, and J is the current density. s It is the source of induced electromagnetic fields.
3. The adaptive finite element method for three-dimensional electromagnetic potential forward modeling according to claim 2, characterized in that: Step S2 specifically includes: The gradient ▽λ term of the Lagrange multipliers is introduced into the governing equations, and is expressed as: ， Using the vector shape function N j and node shape function N k The vector potential and scalar potential are discretized as follows: ， in, and Let be the number of edges and nodes of the unit cell, respectively; the above equations are discretized as follows: ， in, yes The approximate solution, where Ω represents all domains of the computational model, and i and j are the number of edges, where i,j = 1, 2, 3, ..., N. edges l and k are the number of nodes, where l,k = 1, 2, 3, ..., N nodes ; The current divergence term is discretized as follows: （13）， The conditional discretization of the vector potential divergence in the Coulomb gauge to be zero is as follows: ， Combining equations (12), (13), and (14), the finite element linear equation system can be expressed as: ， Among them, the terms on the right-hand side of equations (13) and (14) related to the excitation source The terms are placed into vector S.