Nesting finite element magnetotelluric three-dimensional forward modeling method

Abstract

The invention discloses a nested finite element magnetotelluric three-dimensional forward modeling method. The method comprises the following steps: constructing a coarse mesh and fine mesh unstructured tetrahedron three-dimensional model; carrying out attribute assignment on the grid material; performing forward modeling calculation on the global coarse mesh model; the coarse grids are de-mapped to the boundaries of the sub-regions; performing forward modeling calculation on the sub-region fine mesh model; and calculating the magnetotelluric response. According to the method, the magnetotelluric response is calculated in two steps through the global coarse mesh and the core region sub-model fine mesh, so that the calculation time and memory consumption are remarkably reduced. The method lays a foundation for improving the magnetotelluric three-dimensional forward and reverse calculation performance, so that the magnetotelluric method can be better applied to the fields of geoscience, mineral exploration and the like.

Description

Technical Field

[0001] The present invention relates to the field of electromagnetic exploration and numerical simulation, and particularly to a nested finite element three-dimensional magnetotelluric forward simulation method. Background Art

[0002] The magnetotelluric method has become an important geophysical exploration technique for studying the deep structure of the earth, and is widely used in crust-mantle structure detection, mineral resource exploration, and engineering exploration. With the popularization of array-based magnetotelluric observation methods, it has become possible to obtain a large number of data sets, which is crucial for reconstructing the electrical structure of the earth. However, performing three-dimensional inversion on large-scale magnetotelluric data sets on a sufficiently fine grid faces significant computational challenges. Forward simulation, as the basis of three-dimensional inversion, is also the most time-consuming part. Therefore, an accurate and efficient three-dimensional forward algorithm is crucial for the successful inversion of large-scale magnetotelluric data sets.

[0003] When solving the magnetotelluric forward problem, boundary conditions are crucial for ensuring the uniqueness of the solution. In three-dimensional magnetotelluric forward modeling, the top boundary is usually used as the excitation source, and the non-homogeneous Dirichlet condition is often applied. For the bottom and side boundaries, the forward modeling results of one-dimensional / two-dimensional models of a half-space or layered model are usually used as boundary conditions. However, the one-dimensional / two-dimensional models do not fully consider the scattering effects caused by three-dimensional anomalies, resulting in modeling errors near the boundaries, i.e., boundary effects. In three-dimensional magnetotelluric forward modeling, to weaken the boundary effects, the model region is usually extended outward to multiple skin depths. The extended region increases the number of grid cells, thus greatly increasing the computational cost of the forward modeling. Especially when using a direct solver to solve the linear equations, as the number of grid cells increases, the size of the coefficient matrix, the calculation time, and the memory requirement will all increase non-linearly. In the existing three-dimensional forward technologies, the extended boundaries lead to an increase in the forward calculation time and memory requirement, limiting the application of three-dimensional magnetotelluric forward and inverse technologies.

[0004] Chinese Patent Application CN 109977585 A discloses a high-precision magnetotelluric forward simulation method, including: establishing a geophysical model, using Chebyshev point dissection for the geoelectric model; establishing a Chebyshev derivative matrix; using the Chebyshev derivative matrix to convert the partial differential equation of magnetotelluric forward calculation into an algebraic equation; boundary condition processing; solving a large-scale non-symmetric dense complex coefficient linear equation system; and calculating the magnetotelluric response. This patent can only calculate two-dimensional models with structured grids. However, the real geological structure is three-dimensional, and two-dimensional algorithms cannot accurately calculate three-dimensional models with complex structures.

[0005] Chinese Patent Application CN 110068873 A discloses a three-dimensional magnetotelluric forward modeling method based on spherical coordinates, belonging to the technical field of geophysical exploration, including the following steps: a. Establish a magnetotelluric control equation; b. In the spherical coordinate system, divide it into several small inverted quadrangular prism grid cells; c. Set the spherical coordinate model parameters and construct a spherical coordinate staggered grid cell with discrete magnetic fields; d. Perform a loop for a single frequency, and numerically discretize the magnetotelluric control equation in the spherical coordinate staggered grid cell; e. Solve the linear equations; f. Calculate the impedance from the spherical coordinate tensor impedance formula, and then substitute it into the Cagniard apparent resistivity calculation formula to obtain the apparent resistivity and phase at the measuring point. This patent uses a spherical coordinate system to reduce the influence of the earth's curvature on the accuracy of magnetotelluric forward modeling. However, this patent uses quadrangular prism grid cells and cannot simulate any complex geological structures. Summary of the Invention

[0006] The purpose of the present invention is to overcome the deficiencies of the above-mentioned prior art and provide a nested finite element three-dimensional magnetotelluric forward modeling simulation method. This method calculates the magnetotelluric response in two steps through a global coarse grid and a fine grid of the sub-model in the core area, significantly reducing the calculation time and memory consumption of three-dimensional forward modeling.

[0007] To achieve the above purpose, the present invention adopts the following technical solutions: A nested finite element three-dimensional magnetotelluric forward modeling simulation method, including the following steps: S1: Construct a non-structured tetrahedron three-dimensional model of a coarse grid and a fine grid. The first set of grids is a coarse grid for discretizing the entire calculation domain, denoted as . The grid cells in the coarse grid are usually larger and have the largest cell size in the extended area. Grid refinement is not performed around the observation points in the core area. Therefore, the scale of the linear equations generated by the discretization of the coarse grid is small and can be quickly solved by a direct solver to obtain a preliminary approximation of the electric field.

[0008] The second set of grids is used to discretize a small sub-model in the core area, called . This grid is refined in the area containing the observation points. The linear equations of the sub-model are very small, enabling the electric field in the core area to be calculated quickly and accurately.

[0009] S2: Assign grid material properties: Set the conductivity, permittivity, and permeability in the coarse grid and the fine grid.

[0010] S3: Forward modeling calculation of the global coarse grid model: Solve the linear equations generated by the discretization of the global model coarse grid to calculate the electric field on the coarse grid. The model boundary is set several skin depths away from the core area. The linear equations and their boundary conditions involved in this step can be expressed as: Equation (1): , Equation (2): , where is the global coarse - grid finite - element linear system, is the electric field on the coarse grid, is the right - hand side including boundary conditions, is the mapping operator, is the electric field obtained from 1D forward modeling. The superscript "GCM" represents the coarse grid of the discretized global model, and "1D" represents the 1D layered - earth model.

[0011] S4: Map the coarse - grid solution to the sub - region boundary: Project the solution vector on the coarse grid onto the edge - centers of the sub - model boundary to obtain the boundary value , which is expressed by the formula: Equation (3): , where is the projection operator used to project the electric field on the coarse grid onto the fine grid, is the unit tangent vector of the edge, is the sparse interpolation matrix used to interpolate the solution on the coarse grid to the vector electric field along the boundary of the sub - model grid. The Delaunay algorithm is used to quickly locate the tetrahedral elements enclosing the edge - center. When the point is inside the tetrahedral element , the interpolation matrix is defined as: Equation (4): , where and represent the vector basis functions of the e - th tetrahedral element and the j - th edge within the element, respectively.

[0012] S5: Forward calculation of the sub - region fine - grid model: Calculate the electric field on the fine grid by solving the linear equations generated from the finite - element discretization of the sub - model, with its boundary condition being . Here, the boundary is close to the core region. The linear equations involved in this step can be expressed as: Equation (5): , Equation (6): , where is the sub - model fine - grid finite - element linear system, is the electric field on the sub - model fine grid, is a mapping operator, and the superscripts "SFM" and "GCM" represent the sub-model fine grid and the global coarse grid, respectively.

[0013] S6: Calculate the magnetotelluric response. The electric and magnetic fields at the observation points can be calculated through vector basis functions and Faraday's law: Equation (7): , Equation (8): , where: is the amplitude of the electric field of the j-th edge in the tetrahedral element e containing the observation point, is , is the angular frequency, is the magnetic permeability of vacuum, is the vector basis function of the j-th edge in the tetrahedral element e.

[0014] By solving the following equation, the impedance tensor at the observation point can be calculated: Equation (9): , where, , , , are the components of the impedance tensor, , , , , , , , are the horizontal components of the electric field and the magnetic field calculated at the measurement point. The subscripts "x" and "y" represent two horizontal directions, and the subscripts "1" and "2" represent two polarization modes.

[0015] The apparent resistivity and the phase can be calculated through the impedance tensor, that is: Equation (10): , Equation (11): , where, Z is the impedance transfer function, and i and j represent the x and y directions, respectively.

[0016] Preferably, in steps S3 and S5, a direct solver is used to calculate the linear equations. Since the finite element coefficient matrices A corresponding to the two polarization modes at each frequency are the same, only the right-hand side terms are different. Therefore, when using the direct solver to solve the linear system, only the matrix A needs to be decomposed once, thus effectively reducing the computational amount. In addition, by adopting the OpenMP parallel technology, the computational efficiency is further improved during the process of solving the linear system.

[0017] Compared with the patent applications in the background art, the present invention adopts unstructured tetrahedral meshes, which can accurately simulate any arbitrarily complex three-dimensional geological structure. In addition, the present invention proposes a nested finite element method for two sets of unstructured tetrahedral meshes to improve the calculation efficiency and reduce the memory.

[0018] The present invention proposes a nested finite element method for three-dimensional magnetotelluric forward modeling. The present invention solves the limitation of the existing forward methods due to the boundary effect, which leads to a large extended area. Different from the traditional finite element method that solves large-scale linear equations through a globally refined mesh, the nested finite element method only needs to solve two smaller linear systems, thus significantly reducing the calculation time and memory consumption. The present invention lays a foundation for improving the three-dimensional forward and inverse calculation performance of magnetotellurics, so that the magnetotelluric method can be better applied to the fields of geoscience and mineral exploration. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 Flowchart of the nested finite element method for three-dimensional magnetotelluric forward modeling; Figure 2 Schematic diagram of the trapezoidal mountain test model and its electrical parameters; Figure 3 Schematic diagram of extracting the core region sub-model from the globally refined mesh of the trapezoidal mountain test model; Figure 4 Comparison diagram of magnetotelluric responses of the three-dimensional trapezoidal hill model under different forward modeling methods; the first to fourth rows respectively represent the apparent resistivity in the xy mode, the apparent resistivity in the yx mode, the phase in the xy mode, and the phase in the yx mode. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0020] The present invention will be further described below in conjunction with the drawings and embodiments.

[0021] The structures, ratios, sizes, etc. shown in the drawings of this specification are only used to cooperate with the content disclosed in the specification for those skilled in this technology to understand and read, and are not used to limit the limited conditions under which the present invention can be implemented. Therefore, they do not have a substantial technical meaning. Any modification of the structure, change of the proportional relationship, or adjustment of the size, without affecting the effects that the present invention can produce and the purposes that can be achieved, should still fall within the scope covered by the technical content disclosed in the present invention. At the same time, the terms such as "upper", "lower", "left", "right", "middle", and "one" cited in this specification are only for the convenience of clear narration and are not used to limit the scope under which the present invention can be implemented. The change or adjustment of their relative relationships, without substantial change in the technical content, should also be regarded as the scope within which the present invention can be implemented.

[0022] For the flow of the nested finite element method for three-dimensional magnetotelluric forward modeling, see Figure 1。To better illustrate the advantages and purposes of this embodiment and verify the accuracy and superiority of the computational performance of the nested finite element magnetotelluric three-dimensional forward algorithm of the present invention, the embodiment discloses an internationally commonly used test model. In the trapezoidal mountain test model, the algorithm of the present invention is compared with the secondary field finite element method proposed by Xie et al. (2022) and the traditional finite element algorithm, and the accuracy of the nested finite element magnetotelluric forward algorithm is verified.

[0023] For example: The embodiment uses an internationally commonly used three-dimensional trapezoidal mountain model to verify the feasibility and accuracy of the developed algorithm. The background resistivity of this model is 100 Ω·m, and the resistivity of air is 10 8 Ω·m. The relative permittivity and relative permeability of the entire model are both 1. The height of the trapezoidal mountain is 0.45 km, the top is a flat square with a side length of 0.45 km. The bottom is also a square with a side length of 2 km, as Figure 2 shown. The range of the global model is -28 km to 28 km horizontally and -7 km to 14 km vertically. The measurement line is arranged along y = 0, with a total of 27 observation points and a frequency of 2 Hz. In the global model, we generated two sets of unstructured tetrahedral meshes with different element sizes, which are called the global fine mesh model and the global coarse mesh model respectively. The global fine mesh model has 474162 edges, and the global coarse mesh model has 65941 edges. The core area of the simulation is meshed with a fine mesh, as Figure 3 shown. A fine mesh sub-model is extracted with a truncation distance of 1.4 km extended outward from the observation point, as Figure 3 shown. The fine mesh sub-model has 119110 edges.

[0024] The magnetotelluric responses of the trapezoidal mountain model are calculated by two methods: nested finite element and traditional finite element. The calculation time of the traditional finite element algorithm is 105.6 seconds, while the completion time of the nested finite element algorithm is 9.8 seconds, with an acceleration of about 10.8 times. In addition, the peak memory usage of the nested finite element algorithm is 0.82 GB, only 11.3% of the 4.6 GB memory consumption of the traditional finite element algorithm. The experiment shows that the nested finite element magnetotelluric three-dimensional forward algorithm proposed by the present invention significantly improves the calculation efficiency and reduces the memory requirement.

[0025] The apparent resistivity and phase results obtained by different methods are as Figure 4As shown. Among them, reference, FE-GFM, NFE-GCM, and NFE-SFM represent the calculation of the secondary field finite element method proposed by Xie et al. (2022), the traditional finite element calculation in the global fine grid, the global coarse grid calculation in the nested finite element, and the submodel fine grid calculation in the nested finite element, respectively. The results of FE-GFM and NFE-SFM are very close to the reference results. The apparent resistivity and phase of the observation stations near the trapezoidal mountain show significant changes due to the undulating terrain. At the observation stations far from the trapezoidal mountain, the apparent resistivity and phase are close to the values of the homogeneous half-space model (apparent resistivity 100 Ω·m, phase 45°). However, the results of FE-GCM deviate from the reference solution. The experiment shows that using the nested finite element algorithm can achieve a calculation accuracy comparable to that of the traditional finite element method, verifying the effectiveness of the method proposed in the present invention.

[0026] The embodiments demonstrate the effectiveness, superiority, flexibility, and generality of the present invention, and can be adjusted and applied according to different requirements and application scenarios. Therefore, the present invention can be implemented in a variety of different ways, providing a solid foundation and technical support for the application of magnetotelluric methods in various fields.

[0027] Although the specific implementation manners of the present invention are described above in conjunction with the accompanying drawings, it is not a limitation on the protection scope of the present invention. Those skilled in the art should understand that, based on the technical solutions of the present invention, various modifications or deformations that can be made without creative efforts by those skilled in the art are still within the protection scope of the present invention.

Claims

1. A nested finite element magnetotelluric three-dimensional forward modeling method, characterized in that: The following steps are involved: S1: Construct a coarse grid and a fine grid unstructured tetrahedron 3D model; the first set of grids is a coarse grid, which is used to discretize the entire computational domain and is denoted as ; The second set of grids is the sub-model fine grid, which is used to discretize small sub-models in the core area, called ; S2: Mesh material property assignment: setting the conductivity, permittivity and permeability in the coarse and fine meshes; S3: Global coarse grid model forward calculation: solve the linear equations generated by the discretization of the global model coarse grid and calculate the electric field on the coarse grid , model boundary Set a few skin depths away from the core region, the linear equations involved in this step and their boundary conditions are expressed as: Formula (1): , Formula (2): , in, is a global coarse-grid finite element linear system, is the electric field on the coarse grid, is the right-hand side term including the boundary conditions, is a mapping operator, is the electric field obtained by one-dimensional forward modeling, the superscript "GCM" indicates the coarse grid of the discretized global model, and "1D" indicates the one-dimensional layered earth model; S4: Mapping the coarse grid solution to the sub-region boundary: Mapping the coarse grid solution vector to the sub-region boundary Projecting from the coarse mesh to the submodel boundary The edge center of , the formula is: Formula (3): , in, is the projection operator, which is used to project the electric field on the coarse grid onto the fine grid. is the unit tangent vector of the edge, is a sparse interpolation matrix used to interpolate the solution on the coarse grid to the vector electric field along the boundary of the submodel grid; the Delaunay algorithm is used to quickly locate the tetrahedral element surrounding the edge center; when the point is inside the tetrahedral element, the interpolation matrix Defined as: Formula (4): ; in, and Respectively represent the vector basis functions of the e-th tetrahedral unit and the j-th edge in the unit; in, is the projection operator, which is used to project the electric field on the coarse grid onto the fine grid. is the unit tangent vector of the edge, is a sparse interpolation matrix used to interpolate the solution on the coarse grid to the vector electric field along the boundary of the submodel grid; the Delaunay algorithm is used to quickly locate the tetrahedral element surrounding the edge center; when the point is inside the tetrahedral element, the interpolation matrix Defined as: Formula (4): ; in, and Respectively represent the vector basis functions of the e-th tetrahedral unit and the j-th edge in the unit; S5: Sub-region fine mesh model forward calculation: The fine mesh is calculated by solving the linear equations generated by the finite element discretization of the sub-model. The electric field on , and its boundary conditions are , here, the border In the area close to the core, the linear equations involved in this step are expressed as: Formula (5): , Formula (6): , Among them, among them, is the sub-model fine mesh finite element linear system, is the electric field on the submodel fine grid, is the mapping operator, and the superscripts "SFM" and "GCM" denote the submodel fine grid and the global coarse grid, respectively; S6: Calculate the magnetotelluric response. The electric field E and magnetic field H at the observation point are calculated using vector basis functions and Faraday's law: Formula (7): , Formula (8): , in: is the electric field amplitude of the jth edge in the tetrahedral element e containing the observation point, yes , is the angular frequency, is the vacuum permeability, is the vector basis function of the jth edge in tetrahedral element e; The impedance tensor at the observation point is calculated by solving the following equation : Formula (9): , in, , , , are the components of the impedance tensor, , , , , , , , The horizontal components of the electric field and magnetic field calculated at the measuring point. The subscripts "x" and "y" represent two horizontal directions, and the subscripts "1" and "2" represent two polarization modes. Calculating apparent resistivity from the impedance tensor and Phase ,Right now: Formula (10): , Formula (11): , Where Z is the impedance transfer function, i and j represent the x and y directions, Im and Re represent the imaginary and real parts, respectively.

2. The nested finite element magnetotelluric three-dimensional forward modeling method according to claim 1 is characterized in that: In the step S1, the grid cells in the coarse grid have the maximum cell size in the extended area, and no grid refinement is performed around the observation points in the core area.

3. The nested finite element magnetotelluric three-dimensional forward simulation method according to claim 1, characterized in that: In the step S1, the sub-model fine grid is refined in the area containing the observation point.

4. The nested finite element magnetotelluric three-dimensional forward modeling method according to claim 1, characterized in that: In step S1, in steps S3 and S5, a direct solver is used to calculate the linear equations; since the finite element coefficient matrices corresponding to the two polarization modes at each frequency The same, only the right-hand side is different, so when using the direct solver to solve the linear system, only the matrix A decomposition is performed once, which effectively reduces the amount of calculation; in addition, by adopting OpenMP parallel technology, the calculation efficiency is further improved in the process of solving linear systems.

Images