A hybrid grid-based high-order finite element forward modeling method for magnetotelluric

Abstract

The present invention patent discloses a method for high-order finite element forward modeling of magnetotellurics based on hybrid grids, which specifically relates to the field of geophysical electromagnetic three-dimensional forward modeling. The Delaunay triangulation algorithm is used to divide the earth's surface into triangular grids, and triangular grids are extended on the triangular grids to form triangular prism units; the areas adjacent to the triangular prism units are discretized into tetrahedral units, and the underground electrical model is discretized using the constructed hybrid grid. Finite element unit analysis is performed through high-order shape functions to discretize the magnetotelluric wave equations and form a magnetotelluric high-order finite element system; the hybrid grid finite element system is solved to obtain the magnetotelluric response signal. The technical solution of the present invention solves the problem of increased computational cost caused by the use of traditional unstructured networks in MT forward modeling under high-frequency conditions. The hybrid grid can adapt to complex geoelectric models with strong terrain fluctuations, which requires less computational cost than the use of traditional unstructured units.

Description

Technical Field

[0001] The present invention relates to the field of geophysical electromagnetic three-dimensional forward modeling, and in particular to a magnetotelluric high-order finite element forward modeling method based on a hybrid grid. Background Art

[0002] Unstructured tetrahedral meshes are capable of accurately depicting complex subsurface anomalies and topographic relief because they can theoretically fit the shape of any geological volume. This type of mesh has become a useful tool for discretizing geoelectrical models using the finite element method in numerical simulations of three-dimensional geophysical electromagnetic fields and has been widely used in magnetotelluric forward modeling. However, unstructured meshes often need to be refined into extremely small elements to avoid pathological conditions. This refinement often results in excessive redundant elements and reduces the computational efficiency of finite element method (FEM) forward modeling. Summary of the Invention

[0003] The present invention aims to provide a hybrid grid-based high-order finite element forward modeling method for magnetotelluric (MT), which solves the problem of increased computational cost caused by using traditional unstructured networks in MT forward modeling at high frequencies.

[0004] In order to achieve the above-mentioned purpose, the technical solution of the present invention is as follows: a hybrid grid-based high-order finite element forward simulation method for magnetotelluric, the simulation method is as follows: the Delaunay triangulation algorithm is used to freely divide the earth's surface into triangular grids, and the triangular grids are extended in the vertical direction of the triangular grids to form triangular prism units; the area adjacent to the triangular prism units is discretized into tetrahedral units, so that the upper and lower interfaces of the triangular prism units are coupled together through triangular planes, and the size of the tetrahedral units is based on the size of the coupled triangles and gradually increases toward the side away from the triangle; the constructed hybrid grid is used to discretize the underground electrical model, and finite element unit analysis is performed through high-order shape functions to discretize the magnetotelluric wave equations to form a magnetotelluric high-order finite element system; the hybrid grid finite element system is solved to obtain the magnetotelluric response signal.

[0005] Furthermore, the distance that the triangular mesh extends is 1 to 2 times the attachment depth, and the vertical dimension of each subdivision layer should be less than 1 / 3 of the attachment depth.

[0006] Furthermore, a second-order shape function is used to discretize the triangular prism element, wherein the shape function of the triangular prism element is composed of a vector shape function and a scalar node shape function, and the second-order triangular prism element has 15 node shape functions;

[0007] The node shape functions at the six vertices of the triangular prism element are:

[0008]

[0009]

[0010] The shape functions of the nodes at the midpoints of the triangle sides are:

[0011]

[0012]

[0013] The shape functions of the nodes at the midpoints of the rectangle sides are:

[0014]

[0015] Vector shape functions are derived from node shape functions; their relationship is:

[0016]

[0017] where edge ij is associated with node j.

[0018] Compared with the existing technology, this solution has the following beneficial effects:

[0019] 1. This scheme uses triangular prism elements to discretize the near-surface region, which can significantly reduce the number of elements (NE) required. This is due to the mesh quality constraints of tetrahedrons. It can also improve the computational efficiency of high-frequency data.

[0020] 2. This solution uses different geoelectrical models to analyze the accuracy and efficiency of hybrid grids in calculating magnetotelluric responses. The three-dimensional layered model and the DTM1 (Dublin Test Model 1) model demonstrate their accuracy and applicability to general geoelectrical models. Hybrid grids can construct finite element system matrices with fewer degrees of freedom and achieve higher solution accuracy. Their superiority is also demonstrated in a terrain model example.

[0021] 3. Triangular prism elements offer the flexibility of lateral refinement controlled by triangles. Since rectangular elements are not affected by aspect ratio, they also offer stability in element quality. Therefore, the application of this hybrid grid is not limited to magnetotelluric forward modeling but can also be used for other geophysical electromagnetic computational problems. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 is a schematic diagram of the generation of a triangular prism unit in an embodiment;

[0023] Figure 2 Schematic diagram of a triangular prism unit in an embodiment;

[0024] Figure 32 is a schematic diagram of a pyramid unit in an embodiment;

[0025] Figure 4 Schematic diagram of a second-order triangular prism unit in an embodiment;

[0026] Figure 5 is a schematic diagram of a three-dimensional layered model in an embodiment;

[0027] Figure 6 Schematic diagram of different grid encryption levels at measuring points of a three-dimensional layered model in an embodiment;

[0028] Figure 7 3D layered model at different grid encryption points in the embodiment;

[0029] Figure 8 3. FIG. 3 is a graph showing the change in apparent resistivity and phase with frequency obtained by using hybrid grids of different orders in a three-dimensional layered model in an embodiment;

[0030] Figure 9 1 is a graph showing the variation of apparent resistivity and phase with frequency obtained by using different grids in a three-dimensional layered model in an embodiment;

[0031] Figure 10 is a sparsity distribution diagram of the system matrix of different grids in the embodiment;

[0032] Figure 11 This is a schematic diagram of the DTM1 model in the embodiment;

[0033] Figure 12 1 is a graph showing the calculation results of apparent resistivity and phase obtained under the DTM1 model using different grids at frequencies of 0.1 Hz, 1 Hz, and 10 Hz, respectively;

[0034] Figure 13 is a slope terrain model diagram in the embodiment;

[0035] Figure 14 This is a result diagram of the xy components of apparent resistivity and phase calculated using different grids when the slope terrain model has a different inclination angle in the embodiment;

[0036] Figure 15 1 is a result diagram of the yx components of the apparent resistivity and phase calculated using different grids when the slope terrain model and the inclination angle change in the embodiment;

[0037] Figure 16 is a three-dimensional geoelectric model diagram of complex undulating terrain in the embodiment;

[0038] Figure 171 is a result diagram of the xy and yx components of the MT response at survey line 1 under a three-dimensional geoelectric model of complex undulating terrain in the embodiment;

[0039] Figure 18 3 is a result diagram of the xy and yx components of the MT response at the survey line 2 under the three-dimensional geoelectric model of complex undulating terrain in the embodiment. DETAILED DESCRIPTION

[0040] The present invention will be further described in detail below through specific embodiments:

[0041] Theoretical introduction:

[0042] Ignore the displacement current in the quasi-static field and assume that the time harmonic variation factor is e iωt , the governing equation of the induced electromagnetic field is written as:

[0043]

[0044]

[0045] Where E and H are the induced electric and magnetic fields, respectively, ω is the angular frequency, and σ and μ are the electrical conductivity and magnetic permeability of the electromagnetic medium, respectively.

[0046] After taking the curl of formula (1), substitute formula (2) into the formula (1) after taking the curl, and the Helmholtz equation of the electric field E is obtained as follows:

[0047]

[0048] Governing equations:

[0049] The electric field can be decomposed into the magnetic vector potential A and the electric scalar potential This can reduce the ill-conditioning of the curl-curl electric field equation, speed up the iterative solver to solve the finite element system, and obtain higher numerical accuracy. The relationship between the potential and the field can be expressed as:

[0050]

[0051]

[0052] The magnetic vector potential equation can be written as follows:

[0053]

[0054] Combining J = σE, substituting formula (4) into the Gaussian divergence equation, the current divergence is zero, which can be expressed as:

[0055]

[0056] When both the curl and divergence of A are given, A is uniquely defined as a constant. Therefore, in order to obtain a unique solution, the Coulomb gauge condition is added:

[0057]

[0058] The Coulomb gauge is applied to enforce the divergence limit and enhance the uniqueness of the finite element (FE) equation. Based on formula (6) derived from Ampere's law, we add Lagrange multiplier terms to eliminate the external current density divergence. Formulas (6) and (7) can be rewritten as follows:

[0059]

[0060]

[0061] In equation (10), -iωσ is continuous and thus specifies the divergence of the magnetic vector potential.

[0062] For geophysical electromagnetic field propagation problems, we use homogeneous Dirichlet boundary conditions at the truncated model boundaries.

[0063] Solve a system of linear equations

[0064] We use the Galerkin method to discretize formula (9), and the discretized expression is as follows:

[0065]

[0066] Where N represents the vector shape function. and Respectively and The approximate solution of Ω represents the entire computational domain. Γ and γ represent the outer and inner boundaries. Similarly, the discrete forms of formulas (7) and (8) are as follows:

[0067]

[0068]

[0069] Where N is a scalar shape function. We use a natural electromagnetic field source to form the right side of the system. When the polarization direction of the magnetic field source is in the y direction, the top interface of the air is loaded with a magnetic field source H0 = (0, 1, 0), and the boundary condition of the outer boundary parallel to the x-axis is The boundary conditions for the outer boundary being perpendicular to the x-axis and the bottom boundary being perpendicular to the x-axis are When the polarization direction of the magnetic field source is in the x direction, the air top interface is loaded with a magnetic field source H0 = (1, 0, 0), and the boundary condition of the outer boundary perpendicular to the x axis is The boundary conditions for the outer and bottom boundaries parallel to the x-axis are The boundary conditions are The relevant boundary conditions are added through the second term in Equation (11). The relationship between the magnetic field H and A is given as follows:

[0070]

[0071] After the outer boundary is discretized, it is assumed that one surface of the i-th unit is at the top interface of the air, denoted as Γ i Combining formula (15) and adding the second term in formula (12) to the MT field source, the expression is as follows:

[0072]

[0073] The computational domain is subdivided into smaller subdomain cells. The vector and scalar potentials for each cell in the grid are represented by piecewise polynomial basis functions, which are expressed as follows:

[0074]

[0075]

[0076]

[0077] Among them, N j represents the vector shape function on the j-th edge in the unit, N j represents the scalar shape function of the jth node in the element, nedge and nnode represent the number of edges and nodes in the element, respectively. Combining Equations (12) to (14) and (17) to (19), the matrix form of the system equation after finite element discretization is as follows:

[0078]

[0079]

[0080] where Γ air is the air top interface, i,j=1,…,n edge and l,k=1,…,n node We use the MUMPS direct solver to solve System equations.

[0081] Example

[0082] A hybrid-grid-based high-order finite element forward modeling method for magnetotelluric (MT) is constructed as follows: the Earth's surface is freely divided into a triangular mesh using the Delaunay triangulation algorithm. The triangular mesh is extended perpendicularly to form the top layer of triangular prism elements. The triangular mesh is extended for a distance of 1 to 2 times the attachment depth, and the vertical dimension of each subdivision layer should be less than 1 / 3 of the attachment depth. The region adjacent to the triangular prism element (i.e., the deep region below it and the air layer above it) is discretized into tetrahedral elements. The upper and lower interfaces of the triangular prism element are coupled via the upper and lower triangular planes of the triangular prism element. The size of the tetrahedral element is based on the size of the coupled triangle and gradually increases away from the triangle. The constructed hybrid-grid discretized underground electrical model is used, and finite element analysis is performed using high-order shape functions to discretize the MT wave equation and form a high-order MT finite element system. This hybrid-grid finite element system is solved to obtain the MT response signal.

[0083] The shape function of the triangular prism element (such as Figure 3 The triangular prism element consists of two triangles at the top and bottom and three rectangles connected at the sides.

[0084] The Whitney-1 form of the basis function is selected as the vector shape function of the triangle and is expressed as W. The corresponding edge is determined by the node number represented by the subscript of W. The two nodes at both ends determine the edge represented by W. Then W can be expressed as:

[0085]

[0086] Among them, ξ i (i=1,2,3) is the plane triangle node shape function established by Lagrange interpolation, and the Whitney-1 shape function represents the area coordinate of the triangular prism unit. Represents the gradient vector of the triangular prism element. The calculation of these shape functions can directly refer to the unit analysis of the triangle in the two-dimensional case, but the three-dimensional shape function of the triangular prism cannot refer to the two-dimensional case. It is necessary to calculate the local coordinate ζ in the height direction (such as Figure 2 As shown) is added to the shape function, we get:

[0087] The bottom triangle vector shape function—that is, the bottom vector function—is expressed as:

[0088] N1=(1-ζ)W 23

[0089] N2=(1-ζ)W 31

[0090] N3=(1-ζ)W 12 (twenty two)

[0091] The top triangle vector shape function—the top vector function—is expressed as:

[0092] N4=ζW 23

[0093] N5=ζW 31

[0094] N6=ζW 12 (twenty three)

[0095] The vector shape function of the rectangle side, that is, the volume vector function, is expressed as:

[0096]

[0097] The scalar node shape function of the triangular prism element is:

[0098]

[0099] Shape function of a pyramidal prism element:

[0100] The expression of the vector shape function of the quadrilateral at the bottom of the pyramid prism unit is:

[0101]

[0102]

[0103]

[0104]

[0105] In the formula, N p Represents the vector shape function of the pyramid unit, N p represents the node shape function of the pyramid unit. The expression of the vector shape function of the hypotenuse of the pyramid unit is:

[0106]

[0107]

[0108]

[0109]

[0110] The expression of the node shape function of the pyramid unit is as follows:

[0111]

[0112]

[0113]

[0114] The principle of this technical solution: The geoelectric model can be divided into the near-surface region and the deep-surface region. The near-surface region is discretized by triangular prism units and can be regarded as the boundary layer near the surface (such as Figure 1 The deep region is discretized by tetrahedrons located below the boundary layer. The triangular prisms and tetrahedrons can share a triangular face, which perfectly combines the triangular prism mesh and the tetrahedron mesh (as shown in Figure 1 Other software that can generate tetrahedral, prism, hexahedral, etc. meshes can also be used to generate hybrid meshes, such as Gmsh.

[0115] Unit Analysis of Triangular Prisms and Pyramids:

[0116] Triangular prism unit (such as Figure 2 As shown, a triangular prism element has six nodes and nine edges. The shape functions (defined in the coordinate system (ξ1, ξ2, ξ)) consist of vector shape functions and node shape functions. A triangular prism element consists of two triangles at the top and bottom, and three surrounding rectangles. See the shape functions of the triangular prism element and the shape functions of the pyramidal prism element for vector and node shape functions.

[0117] The shape function of the pyramid unit (such as Figure 3 As shown, a pyramid element has five nodes and eight edges. Node and edge numbers are defined in the coordinate system (ξ1,ξ2,ξ). It consists of vector shape functions and node shape functions. A pyramid prism element consists of a rectangle at the base and four surrounding triangles. For more information about the shape functions of a pyramid element, see the shape functions of a triangular prism element and the shape functions of a pyramid prism element.

[0118] Shape function of the second-order triangular prism element:

[0119] Secondary triangular prism unit such as Figure 4 As shown. The quadratic triangular prism element has 15 nodes, 18 edge variables, and 10 surface variables. i (i=1,2,...,28) is the variable in the unit. Figure 4 As shown, there are 15 node shape functions that define the node and edge numbers in the coordinate system (ξ1,ξ2,ξ).

[0120] The node shape functions at the six vertices of the triangular prism element are:

[0121]

[0122] The shape functions of the nodes at the midpoints of the triangle sides are:

[0123]

[0124] The shape functions of the nodes at the midpoints of the rectangle sides are:

[0125]

[0126] The vector shape functions are derived from those node shape functions. Their relationship is:

[0127]

[0128] Where edge ij is associated with node j. Vector shape functions are on edges, while node shape functions are on nodes.

[0129] by Figure 4 To illustrate, the vector shape function N 1,14 It is a vector function acting on edge 1-14 (edge ​​1-14 is the edge between node 1 and node 14) and has directionality.

[0130] Because the vector shape function has directionality, when the relationship between the vector shape function and the node shape function is written as When defining edge ij as being associated with node j, for example, edge 1-14 is associated with node 14; when the relationship between the vector shape function and the node shape function is written as When defining edge ji to be associated with node i, for example, edge 14-1 is associated with node 1.

[0131] According to equations (21) to (24), the second-order vector shape function of the triangular prism can be easily derived. Taking nodes 1, 7, and 11 as an example:

[0132]

[0133]

[0134]

[0135] The derivation methods for other second-order vector shape functions are the same as above. We can classify vector shape functions into the following categories.

[0136] The vector shape functions for the bottom and top triangle edges can be expressed as:

[0137]

[0138] The vector shape function of the rectangle edge (volume vector function) can be expressed as:

[0139]

[0140] The vector shape functions of the bottom and top triangle faces can be expressed as:

[0141]

[0142] The vector shape function of a rectangular face can be expressed as:

[0143]

[0144]

[0145] Advantages of hybrid grids:

[0146] The advantages of triangular prism elements in finite element magnetotelluric forward modeling are demonstrated by numerical examples of layered models. Apparent resistivity and phase are responses, which helps us analyze the accuracy of the calculation.

[0147] To simplify the analysis, the highest frequency is used to calculate the minimum adhesion depth, which is used to stretch the triangular prism elements in the near-surface region. Although this size is redundant for calculating electromagnetic fields at lower frequencies, it does not affect the superiority of the triangular prism compared to the tetrahedron.

[0148] Meshing of 3D layered models:

[0149] The layered model is divided into three layers (such as Figure 5 (As shown): The resistivities of the first, second, and third layers are 200Ωm, 1000Ωm, and 200Ωm, respectively. The thickness of the first and second layers is 0.5km. Measurement points are located at Y = 0m and Z = 0m, with a spacing of 300m along the X direction from -1200m to 1200m. The entire computational model is 20km × 20km × 70km.

[0150] The size of the grid cells can be divided into horizontal and vertical dimensions, and is controlled by the maximum element size (MES) and the growth rate. These dimensions are primarily controlled by the discretization of triangular prisms in the near-surface region. The size of the tetrahedrons in the deep region is determined by the base size of the triangular prisms.

[0151] For the near-surface region, the horizontal size is controlled by the triangles used to mesh the ground. We use the Triangle Maximum Unit Size (TriMES) and the Triangle Growth Rate (TriGR) to generate these triangles. We use the Prism Growth Rate (PriGR) to control the multiplier by which the vertical size of the prisms grows and to control the total thickness of the prism elements.

[0152] For the deep region, the size of the top tetrahedron coplanar with the coupling plane is controlled by the size of the coupled triangles. The sizes of the other tetrahedra are controlled by the tetrahedron growth rate (TetGR).

[0153] We use the four types of grids listed in Table 1 to discretize the three-dimensional layered model.

[0154] Table 1:

[0155]

[0156] In Table 1, MES represents the maximum element size. The horizontal size is controlled by the surface's triangular MES (Tri MES) and the triangle growth rate (TriGR). TriGR represents the maximum element growth rate for triangular element size. The vertical size is controlled by the prism layer thickness and the tetrahedron growth rate (TetGR). TetGR represents the maximum element growth rate for tetrahedron element size. The prism growth rate (PriGR) represents the prism element growth rate perpendicular to the prism layer thickness.

[0157] As shown in Table 1 above, Grid A uses the hybrid grid of the embodiment, with the upper portion of the hybrid grid being a triangular prism grid and the lower portion being a tetrahedron grid. The horizontal size is determined by the size of the top triangular unit cell. The maximum unit size of the triangle is 1500 meters, and the triangle growth rate is 1.3. The vertical size of the hybrid grid is determined by both the vertical sizes of the triangular prisms and tetrahedrons. The vertical size of the first prism layer is 8 meters, and the prism growth rate is equal to the tetrahedron growth rate. The total number of prism layers is 10. Below the prism layer, the vertical size of the tetrahedron elements is controlled by the maximum unit size of the triangle and the tetrahedron growth rate. The tetrahedron growth rate in Grid A is 1.5.

[0158] Meshes B1, B2, and B3 all use tetrahedral meshes, and the only difference between these meshes is the thickness of the meshes. Mesh B1 is the coarsest mesh. The horizontal size of the mesh B1 unit is the same as the triangular prism unit, and the maximum unit size of the triangle is 1500m. Under the premise of effective mesh division, the vertical size of the first layer of mesh B1 unit is about 500m to obtain good mesh quality. The thickness of meshes B2 and B3 is smaller than that of mesh B1. The maximum unit size of the triangle and the first layer thickness of mesh B2 are 250m and about 200m respectively. The maximum unit size of the triangle and the first layer thickness of mesh B3 are 100m and about 65m respectively. The mesh refinement of the four meshes at the measuring points is compared (such as Figure 6 As shown, (a) is the mesh refinement at the measuring point of grid A, (b) is the mesh refinement at the measuring point of grid B1. (c) is the mesh refinement at the measuring point of grid B2. (d) is the mesh refinement at the measuring point of grid B3) and the mesh diagram of the cross section (as shown Figure 7 As shown, (a) is a cross-section of mesh A. The dark portion represents a triangular prism mesh, and the remaining portion represents a tetrahedral mesh. (b) is a cross-section of mesh B1. (c) is a cross-section of mesh B2. (d) is a cross-section of mesh B3.

[0159] Higher-order shape functions:

[0160] In the frequency range of 1E5 Hz to 1E-1 Hz, the MT response is calculated using a hybrid grid of different orders. A total of 26 frequency points are calculated, of which 10 frequency points are equally spaced between 1E5 Hz and 1E4 Hz, and the remaining frequency points are equally spaced between 1E4 Hz and 1E-1 Hz. The apparent resistivity and phase change with frequency are shown in Figure 2. Figure 8 (The apparent resistivity and phase are shown in (a) and (b), respectively. The errors in apparent resistivity and phase are shown in (c) and (d), respectively. The absolute errors are obtained by normalizing the semi-analytical solution.) When the shape function is first-order, the maximum absolute errors in apparent resistivity and phase are 27.8510% and 6.8170%, respectively. When the shape function is second-order, the maximum absolute errors in apparent resistivity and phase are 2.0968% and 0.4063%, respectively. All calculations were performed in MATLAB 2019, using an Intel Core i5-10400 CPU @ 2.90 GHz and 32 GB of RAM.

[0161] Mesh quality analysis:

[0162] The MT responses were calculated in the frequency range of 1E4 Hz to 1E-1 Hz. Figure 9 The frequency variations of apparent resistivity and phase are shown in (a) and (b), respectively. The errors in apparent resistivity and phase are shown in (c) and (d), respectively. The absolute errors are normalized using the semi-analytical solution. Within this frequency range, the accuracy of apparent resistivity and phase calculated using different grids varies significantly. The root mean square error (RMSE) is used to measure the overall error in apparent resistivity and phase, with reference values ​​calculated using the one-dimensional layered model formula (semi-analytical solution).

[0163] For meshes consisting only of tetrahedrons, the smaller the mesh size, the smaller the calculated response error. The vertical size of the tetrahedral element is a key factor affecting the calculation accuracy. Considering the quality of the tetrahedral element, the horizontal size of the tetrahedral element needs to be consistent with its vertical size. However, this will lead to a sharp increase in the number of tetrahedral elements (NE) and degrees of freedom, and the calculation time will increase with the increase in degrees of freedom.

[0164] Hybrid meshes with triangular prism elements can overcome this limitation. High accuracy can be achieved even with triangular prism elements with an aspect ratio greater than 100. When accuracy is not limited by aspect ratio, the number of required elements is greatly reduced, as is the number of DoFs required. The resulting NE, DoF, and mean absolute error (MAE) are shown in Table 2 below.

[0165] Table 2:

[0166]

[0167] Note: The number of elements (NE) represents the total number of elements in a hybrid or tetrahedral mesh. The number of DoFs and the number of nonzeros (NNZ) reflect the basic characteristics of the FE system matrix. Based on the semi-analytical solution, the mean absolute error is used to measure the error between different mesh generation results. This example provides the memory usage for four types of mesh calculations.

[0168] Table 2 compares the number of finite element system elements, degrees of freedom, calculation time, and numerical solution accuracy of different meshes. A higher order helps to improve the accuracy of the numerical solution of the prismatic element. Comparing the calculation results of mesh A and mesh B1, it can be seen that mesh A has more degrees of freedom and a more accurate numerical solution. Comparison of the calculation results of meshes A and B2 shows that even if the number of DoFs of the tetrahedral mesh is 8 times that of the hybrid mesh, the accuracy of the hybrid mesh is much higher than that of the tetrahedral mesh. Comparison of the results of meshes A and B3 shows that when the number of DoFs of the tetrahedral mesh is 56 times that of the hybrid mesh, the accuracy of the tetrahedral mesh can be at the same level as the hybrid mesh. This shows that the mesh type not only affects the size and number of degrees of freedom of the system matrix, but may also affect the solution characteristics of the system matrix.

[0169] Sparsity distribution of system matrices of different grids (such as Figure 10 As shown in the figure (a) through (d) represent the distribution of the system matrices generated by meshes A, B1, B2, and B3, respectively), this demonstrates that the hybrid mesh's system matrix exhibits superior solution properties. For tetrahedral meshes, the clustering of the system matrix gradually decreases as the number of degrees of freedom increases. While the hybrid mesh's degrees of freedom lie between those of meshes B1 and B2, its system matrix exhibits significantly better clustering than the other two tetrahedral meshes.

[0170] DTM1 model analysis:

[0171] The hybrid grid mainly uses triangular prism elements to process the near-surface area, which is often related to the calculation accuracy of high-frequency data. In order to verify the role of hybrid grid in improving the calculation accuracy of high-frequency data solutions, we increased the complexity of the model based on the three-dimensional layered model and selected the DTM1 model to demonstrate the advantages of triangular prism elements. A survey line was selected to record the response of apparent resistivity and phase (such as Figure 11 shown), where Figure 11 The middle (a) view is a vertical section of the DTM1 model. Figure 7 The middle (b) view is a horizontal section of the DTM1 model, and the survey line is shown in the horizontal section.

[0172] We set up four different meshes in total, namely the hybrid mesh of the embodiment and three tetrahedral meshes. The parameters of the different meshes are shown in Table 3 below.

[0173] Table 3:

[0174]

[0175] As can be seen from Table 3 above, in the hybrid mesh, TriMES, which controls the horizontal element size, is set to 2000m, and the thickness of the first layer of triangular prism elements is 318m. In the tetrahedral mesh, TriMES gradually changes from 2000m to 500m. The size of the first layer of tetrahedral elements is approximately equal to TriMES on the surface. In all four meshes, the growth rate of triangular and tetrahedral elements TriGR and TetGR is set to 1.5. The system equations established for the tetrahedral mesh with TriMES of 1000m and the hybrid mesh with TriMES of 2000m have the same degrees of freedom. The tetrahedral mesh with TriMES of 500m can be regarded as a dense tetrahedral mesh.

[0176] The calculation results of apparent resistivity and phase at frequencies of 0.1Hz, 1Hz and 10Hz are as follows Figure 12 As shown, Figure 12 (a) and (d) represent the apparent resistivity and phase results at a frequency of 10 Hz, (b) and (e) at a frequency of 1 Hz, and (c) and (f) at a frequency of 0.1 Hz. The Nam code response was obtained using the FEM method. The Mackie code response was obtained using the finite difference method. The wsinv3dmt and mt3dinv responses were obtained using the FD code. At a frequency of 10 Hz, the denser the tetrahedral mesh, the smaller the oscillations in the apparent resistivity and phase curves, but neither is as smooth as the response curves obtained with the hybrid mesh. At a frequency of 1 Hz, the oscillations in the tetrahedral mesh response curve decrease, but the smoothness still lags behind that of the hybrid mesh response curve. At a frequency of 0.1 Hz, the smoothness of the tetrahedral mesh response curve is comparable to that of the hybrid mesh. When calculating high-frequency electromagnetic fields from magnetotellurics, the discretization equations using hybrid meshes offer greater stability and higher accuracy.

[0177] When both the hybrid and tetrahedral meshes have horizontal dimensions of approximately 2000 m, the hybrid mesh exhibits more degrees of freedom and greater response stability, as the triangular prism elements in the hybrid mesh can reduce their vertical dimensions without degrading mesh quality. When the tetrahedral mesh is refined, its horizontal dimension is reduced by approximately 1000 m and its vertical dimension is set to approximately 1000 m, ensuring that the tetrahedral mesh has a similar number of degrees of freedom (DoF) as the hybrid mesh. The apparent resistivity and phase curves calculated using this dense tetrahedral mesh still exhibit significant oscillations. The tetrahedral mesh is further refined, with both horizontal and vertical dimensions set to approximately 500 m. At this point, the number of DoFs is nearly double that of the hybrid mesh, but the results remain unstable. This demonstrates that the hybrid mesh is indeed more efficient than a pure tetrahedral mesh for computing high-frequency data.

[0178] Slope terrain model analysis:

[0179] like Figure 13 As shown, the slope angles in this embodiment are 10, 20, and 30 degrees. Due to the volume effect of electromagnetic methods, calculating high-frequency electromagnetic responses is more difficult when the terrain is highly undulating. We analyzed the improvement in hybrid grid computational efficiency under different terrain conditions by adjusting the terrain undulation angle.

[0180] We set three mesh types: hybrid mesh, general tetrahedral mesh, and refined tetrahedral mesh. The parameters of different meshes are shown in Table 4.

[0181] Table 4:

[0182]

[0183] As shown in Table 4 above, the hybrid mesh and the standard tetrahedral mesh have the same horizontal element size, with a TriMES value of 120 m for both. The reinforced tetrahedral mesh has a smaller horizontal element size, with a TriMES value of 50 m. The thickness of the first layer of triangular prism elements in the hybrid mesh is 15 m. The vertical dimensions of the first tetrahedron in the other two tetrahedral meshes are approximately equal to their horizontal dimensions. In these meshes, the mesh growth rate for both triangles and tetrahedrons is 1.5. The system matrices generated by the hybrid mesh and tetrahedral meshes with the same horizontal element size have a similar number of degrees of freedom, while the system matrix generated by the tetrahedral mesh with the smallest horizontal element size has almost twice the number of degrees of freedom.

[0184] When the tilt angle changes, we calculate the xy components under different types of grids ( Figure 14 ) and the yx component ( Figure 15 ) response. For the xy component ( Figure 14 ), the apparent resistivity resolution is comparable, with slight oscillations observed in the results of general tetrahedral meshes. The differences in the results from different meshes are evident in the phase results, but the nature of the oscillations remains largely unchanged with tilt angle. However, hybrid meshes offer very stable results for both apparent resistivity and phase, across various tilt models.

[0185] The results of the yx component agree with those of the xy component in three key respects. First, the results oscillate, but the hybrid grid results have the smallest oscillation. Second, the oscillation differences between the results of the different grids barely change with the inclination of the slope terrain model. Third, the oscillation differences are even more pronounced in the phase results.

[0186] The difference between the yx component and the xy component is that the response of the yx component has a more violent oscillation (by comparing Figure 14 and Figure 15),in Figure 14 and Figure 15 (a1), (b1), and (c1) represent the apparent resistivity calculation results, while (a2), (b2), and (c2) represent the phase calculation results. (a), (b), and (c) represent the calculation results for tilt angles of 10, 20, and 30 degrees, respectively. The slope varies along the X direction, causing greater error in the yx component. Therefore, larger errors may occur in directions with greater terrain undulation.

[0187] Terrain relief model analysis:

[0188] In the high frequency range, the hybrid grid can obtain high-precision electromagnetic responses regardless of whether it is flat terrain or undulating terrain with slopes. In order to further verify the adaptability of the hybrid grid to undulating terrain, we designed a three-dimensional geoelectric model of complex undulating terrain ( Figure 16 (a) shows the undulating terrain and two survey lines. (b) shows the entire 3D model. The accuracy of the magnetotelluric responses calculated using a hybrid mesh and a tetrahedral mesh is compared with reference to the actual terrain. The calculation frequency is 1E4 Hz.

[0189] We set three mesh types: hybrid mesh, general tetrahedral mesh, and refined tetrahedral mesh. The parameters of different meshes are shown in Table 5.

[0190] Table 5:

[0191]

[0192] Table 5 shows that the hybrid mesh and the standard tetrahedral mesh have the same horizontal element size, with a TriMES of 100 m. The refined tetrahedral mesh has a smaller horizontal element size, with a TriMES of 20 m. The TriGR for all three meshes is 1.5. The first-layer mesh thicknesses for the three meshes are 15 m, 80 m, and 20 m, respectively. To ensure that the system matrices generated by the hybrid mesh and the standard tetrahedral mesh have similar DoFs, the growth rates for the hybrid mesh and the standard tetrahedral mesh are set to 1.5 and 1.34, respectively. The TetGR for the refined tetrahedral mesh is 1.5.

[0193] The MT response results of the xy and yx components on the two survey lines are as follows: Figure 17 、 18 As shown, (a1) and (b1) represent apparent resistivity, (a2) and (b2) represent phase results. Figure 17), the results from the general tetrahedral mesh are the least accurate. The hybrid mesh results are closer to those from the refined tetrahedral mesh. The refined tetrahedral mesh results exhibit significant oscillations in the phase response. This demonstrates that the hybrid mesh can achieve stable results while constructing a system matrix with a lower number of degrees of freedom.

[0194] The response results of measurement line 2 ( Figure 18 ) also shows that when the number of DoFs is similar, hybrid meshes provide higher accuracy than standard tetrahedral meshes. Refined tetrahedral meshes yield even more accurate results, but the system matrix has nearly 2.6 times the DoF of the hybrid mesh. Furthermore, the phase results for the refined tetrahedral mesh exhibit stronger oscillations than those for the hybrid mesh.

[0195] The above are only embodiments of the present invention, and common knowledge such as the specific structure and / or characteristics of the scheme are not described in detail here. It should be pointed out that for those skilled in the art, without departing from the structure of the present invention, several variations and improvements can be made, which should also be regarded as the scope of protection of the present invention, and these will not affect the effect of the implementation of the present invention and the practicality of the patent. The scope of protection required by this application shall be based on the content of its claims, and the specific implementation methods and other records in the specification can be used to interpret the content of the claims.

Claims

1. A hybrid grid-based high-order finite element forward modeling method for magnetotelluric, characterized by: The simulation method is as follows: The Delaunay triangulation algorithm is used to freely divide the earth's surface into triangular grids, and the triangular grids are extended in the vertical direction of the triangular grid to form triangular prism units; The area adjacent to the triangular prism unit is discretized into tetrahedral units, thereby coupling the upper and lower interfaces of the triangular prism unit through the triangular plane. The size of the tetrahedral unit is based on the size of the coupled triangle and gradually increases away from the triangle. Using the constructed hybrid grid discretized underground electrical model, finite element unit analysis is performed through high-order shape functions to discretize the magnetotelluric wave equation and form a magnetotelluric high-order finite element system. Solve the hybrid grid finite element system to obtain the magnetotelluric response signal; A second-order shape function is used to discretize the triangular prism element, wherein the shape function of the triangular prism element is composed of a vector shape function and a scalar node shape function, and the second-order triangular prism element has 15 node shape functions; Among them, the node shape functions at the six vertices of the triangular prism unit are: The shape functions of the nodes at the midpoints of the triangle sides are: The shape functions of the nodes at the midpoints of the rectangle sides are: in, i Represents a node, N i Indicates the unit i The node shape function of each node; Vector shape functions are derived from node shape functions; their relationship is: in, For the edge - j The vector shape function of For the unit j The node shape function of the nodes, the edge - j With node j associated.

2. The hybrid grid-based magnetotelluric high-order finite element forward modeling method according to claim 1, characterized in that: The distance that the triangular mesh extends is 1 to 2 times the attachment depth, and the vertical dimension of each subdivision layer should be less than 1 / 3 of the attachment depth.

Images