Regional decomposition magnetotelluric three-dimensional forward modeling method and system

Abstract

The present invention belongs to the technical field of electromagnetic exploration and numerical simulation. A regional decomposition magnetotelluric three-dimensional forward modeling method and system are provided, which introduces overlapping DDM to decompose the large-scale MT problem into multiple smaller subdomains, thereby reducing the computational complexity; a direct solver is used to solve each small linear system in the subdomain, which not only avoids the convergence problem of the iterative solver, but also reuses the LU decomposition results in the Schwarz iteration process to speed up the calculation speed; a multi-grid scheme is integrated to accelerate convergence, and good initial conditions are provided through coarse grid solutions, thereby reducing the number of iterations and further improving the computational efficiency; the use of unstructured grid technology enhances the adaptability of complex geological structures, thereby improving the simulation accuracy; the present invention develops an airborne magnetotelluric three-dimensional forward modeling technology with fast computing speed and low memory requirements, which provides a basis for theoretical research on airborne magnetotelluric and three-dimensional data inversion.

Description

Technical Field

[0001] The invention relates to the technical field of electromagnetic exploration and numerical simulation, and in particular to a regional decomposed magnetotelluric three-dimensional forward simulation method and system. Background Art

[0002] The statements in this section merely provide background art related to the present invention and do not necessarily constitute prior art.

[0003] As the discovery of shallow and exposed ore bodies is gradually exhausted, exploration work is gradually extending to greater depths, and the difficulty is also significantly increased. In order to meet global resource needs, it is urgent to develop exploration technologies suitable for deep targets. The airborne magnetotelluric method is a geophysical technology that uses natural electromagnetic sources for detection. It measures magnetic field changes through airborne receivers to obtain electrical information about underground structures. Airborne magnetotelluric not only has a greater detection depth, but is also more efficient than ground-based magnetotelluric methods and is suitable for areas with significant terrain undulations. After decades of development, thanks to advances in electromagnetic theory, data processing and electronic technology, airborne magnetotelluric technology has made important breakthroughs.

[0004] In recent years, a variety of natural field airborne magnetotelluric measurement systems have been introduced and have shown great potential in deep mineral exploration. As an important tool for airborne magnetotelluric data interpretation, three-dimensional forward modeling and inversion technology obtain underground resistivity distribution by processing observation data. The computational efficiency of three-dimensional forward modeling determines the performance of three-dimensional inversion. With the continuous advancement of high-performance computing technology, electromagnetic forward modeling has made significant progress.

[0005] At present, the mainstream three-dimensional magnetotelluric forward modeling methods include integral equation method, finite difference method, finite volume method, finite element method and their hybrid methods. The vector finite element method supports the use of unstructured tetrahedral grids and has high flexibility to simulate complex geological structures. Therefore, it is widely used in three-dimensional electromagnetic forward modeling and inversion. The scale of airborne magnetotelluric forward modeling problems is usually much larger than that of ground-based magnetotelluric. The main reasons include: (1) the coverage of aerial measurements is larger; (2) the frequency of airborne magnetotelluric measurements is higher (30Hz~20000Hz), which is more susceptible to complex terrain and shallow geological structures. Accurate simulation of the electromagnetic response of mountainous terrain and complex geology requires extremely fine grid modeling, which forms a large-scale linear equation system with millions of unknowns, consuming huge computing resources.

[0006] When solving large linear systems of equations in MT forward modeling, either iterative solvers or direct solvers can be used. The advantage of iterative solvers is their low memory consumption, but their convergence speed depends largely on the condition number of the coefficient matrix, the choice of preconditioner, and the initial guess. As the frequency decreases and the matrix size increases, the linear system of Helmholtz equations for the electric field obtained by finite element discretization becomes more ill-conditioned, resulting in an increase in the condition number, making it difficult for iterative solvers to converge. Although direct solvers provide accurate and stable solutions, they usually consume much more memory than iterative solvers. As the problem size increases, the memory and computing time requirements increase nonlinearly. Therefore, whether iterative or direct solvers are used, the solution of large-scale airborne MT forward problems faces significant limitations. Summary of the invention

[0007] In order to address the deficiencies of the prior art, the present invention provides a method and system for three-dimensional magnetotelluric forward modeling by regional decomposition, introduces overlapping domain decomposition methods (Domain Decomposition Methods, DDM), decomposes a large-scale magnetotelluric (MT) forward modeling problem into multiple smaller subdomains, and reduces the computational complexity; uses a direct solver to solve each small linear system in the subdomain, avoids the convergence problem of the iterative solver, and reuses the LU decomposition results in the Schwarz iteration process, thereby speeding up the calculation speed; integrates a multi-grid scheme to accelerate convergence, and provides good initial conditions through a coarse grid solution, thereby reducing the number of iterations and further improving the computational efficiency; uses unstructured grid technology to enhance the adaptability of complex geological structures, thereby improving the simulation accuracy.

[0008] In order to achieve the above object, the present invention adopts the following technical solution:

[0009] In a first aspect, the present invention provides a regional decomposed magnetotelluric three-dimensional forward modeling method.

[0010] A regional decomposition magnetotelluric three-dimensional forward modeling method includes the following processes:

[0011] Generate overlapping sub-regions in the unstructured grid of the computational domain for forward simulation according to each tetrahedral unit in the tetrahedral grid of the computational domain;

[0012] Performing vector finite element discretization on the overlapping sub-regions to obtain a finite element linear equation group for each overlapping sub-region;

[0013] According to the vector finite element discretization results, the boundary conditions of the overlapping sub-regions are set, and the finite element linear equations of all overlapping sub-regions are solved;

[0014] Update the global solution by combining the solutions in all overlapping subregions;

[0015] It is determined whether the convergence condition is reached according to the global solution. If the convergence condition is reached and the measuring point is in a tetrahedral unit, the electric field component and the magnetic field component at the measuring point are calculated to obtain the impedance tensor at the observation point and the apparent resistivity of the airborne magnetotelluric system. If the convergence condition is not reached, the boundary condition of the overlapping sub-region is reset, and the iteration is continued until the convergence condition is reached.

[0016] In a second aspect, the present invention provides a regional decomposed magnetotelluric three-dimensional forward modeling system.

[0017] A regional decomposition magnetotelluric three-dimensional forward modeling system, comprising:

[0018] The overlapping sub-region generating unit is configured to: generate overlapping sub-regions in the unstructured grid of the computational domain of the forward simulation;

[0019] The vector finite element discretization unit is configured to: perform vector finite element discretization on the overlapping sub-regions to obtain a finite element linear equation group for each overlapping sub-region;

[0020] The finite element linear equation solving unit is configured to: set the boundary conditions of the overlapping sub-regions according to the vector finite element discretization results, and solve the finite element linear equations of all overlapping sub-regions;

[0021] The global solution updating unit is configured to: update the global solution by combining the solutions in all overlapping sub-regions;

[0022] The forward calculation unit is configured to: determine whether the convergence condition is reached according to the global solution, if the convergence condition is reached, and when the measuring point is in a tetrahedral unit, calculate the electric field component and the magnetic field component at the measuring point to obtain the impedance tensor at the observation point and the apparent resistivity of the airborne magnetotelluric system; if the convergence condition is not reached, reset the overlapping sub-region boundary condition, and continue to iterate until the convergence condition is reached.

[0023] Compared with the prior art, the present invention has the following beneficial effects:

[0024] 1. The present invention proposes a regional decomposition aeromagnetotelluric three-dimensional forward modeling method. For large-scale three-dimensional aeromagnetotelluric forward modeling problems, the additive Schwarz overlapping domain decomposition method is used to divide the large-scale problem into multiple small-scale sub-domains, effectively reducing the computational complexity.

[0025] 2. The present invention uses a direct solver in the subdomain, avoiding the convergence problem of the traditional iterative solver, and reuses the LU decomposition results in the Schwarz iteration, greatly reducing the consumption of computing resources; for areas where parameters remain unchanged during the inversion process (such as air layers or filling areas), the LU decomposition results can be reused, further reducing the burden of inversion calculations.

[0026] 3. The multi-grid method is introduced in the Schwarz iteration, and the coarse grid solution is used as the initial condition to accelerate the convergence speed and reduce the number of iterations; the unstructured tetrahedral grid is used to discretize the complex geological structure, which enhances the algorithm's adaptability to complex terrain and simulation accuracy; compared with the traditional finite element method, the regional decomposition magnetotelluric three-dimensional forward modeling method greatly reduces memory usage and calculation time.

[0027] 4. The present invention is suitable for large-scale aeromagnetotelluric forward and inversion calculations, which effectively breaks through the limitations of existing technologies in quickly solving large-scale aeromagnetotelluric forward modeling. This method can maintain high accuracy under complex terrain and large-scale computing conditions, significantly improve computing efficiency and reduce memory consumption, providing a solid technical foundation for aeromagnetotelluric three-dimensional forward and inversion, and contributing to the further promotion and application of this technology in multiple application fields.

[0028] Advantages of additional aspects of the present invention will be given in part in the following description, and in part will become obvious from the following description, or will be learned through practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] The accompanying drawings in the specification, which constitute a part of the present invention, are used to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute improper limitations on the present invention.

[0030] Figure 1 A flow chart of a regional decomposed magnetotelluric three-dimensional forward modeling method provided in Example 1 of the present invention;

[0031] Figure 2 The conductivity distribution and fine grid diagram of the COMMEMI-3D3 model provided in Example 1 of the present invention;

[0032] Figure 3 A grid partition diagram of the COMMEMI-3D3 model provided in Example 1 of the present invention;

[0033] Figure 4 Convergence diagram of Schwarz iteration under different overlapping layers and initial boundary conditions provided in Example 1 of the present invention;

[0034] Figure 5The calculation results of the regional decomposition algorithm under different overlapping layers and initial boundary conditions provided in Example 1 of the present invention are as follows: Figure 5 The upper sub-graphs in the figure are the apparent resistivity. Figure 5 The sub-graphs in the middle and lower layers are the relative errors of the apparent resistivity calculated by the regional decomposition algorithm and the traditional finite element algorithm;

[0035] Figure 6 A schematic diagram of a regional decomposed magnetotelluric three-dimensional forward modeling system provided in Example 2 of the present invention. DETAILED DESCRIPTION

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

[0037] It should be noted that the following detailed descriptions are exemplary and are intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meanings as those commonly understood by those skilled in the art to which the present invention belongs.

[0038] In the absence of conflict, the embodiments of the present invention and the features of the embodiments may be combined with each other.

[0039] Embodiment 1:

[0040] In order to accurately simulate the electromagnetic field characteristics in complex geological structures and solve the problem of low computational efficiency of large-scale three-dimensional airborne magnetotelluric forward modeling, the present invention has developed a new finite element domain decomposition method, which combines the additive overlapping Schwarz method, the multigrid scheme and the direct solver, and uses partitioned unstructured tetrahedral grids. This method solves the computational bottleneck of large linear equations by decomposing large forward problems into small problems in subdomains; the present invention uses unstructured tetrahedral grids to discretize complex geological structures and divides the entire computational domain into overlapping subdomains. The linear equations obtained in each subdomain are directly solved. The solver is connected to the solution, so as to achieve fast calculation, and the matrix decomposition results can be reused; finally, the solution of the global problem is gradually approached by additive Schwarz iteration. During the additive Schwarz iteration, there is no need to exchange data between subdomains. The global solution is obtained by using the solution of each subdomain, and then used as the subdomain boundary condition of the next iteration; in order to accelerate the convergence speed of the additive Schwarz iteration, the present invention combines the multi-grid scheme and uses the coarse grid solution instead of the analytical solution of the layered earth model as the initial subdomain boundary condition; the accuracy and superiority of the algorithm are verified by the layered model and complex geological model embodiments. Numerical results show that compared with the traditional finite element method, the regional decomposition magnetotelluric three-dimensional forward modeling method of the present invention significantly reduces the computational memory requirements and time while maintaining accuracy.

[0041] Specifically, Figure 1 As shown, the following process is included:

[0042] S1: Construct the governing equations for forward modeling of aeromagnetotelluric.

[0043] The governing equation of the forward modeling of aeromagnetotelluric can be expressed as:

[0044] (1);

[0045] in, and represent the electric field and magnetic field respectively, is the angular frequency, is the frequency, parameter , and represent the electrical conductivity, magnetic permeability and dielectric constant respectively, represents the curl operator.

[0046] In the computational domain The border Imposing nonhomogeneous Dirichlet boundary conditions on :

[0047] (2);

[0048] in, Indicates at the border For the forward electromagnetic problem in the air, two polarization modes must be solved for each frequency, namely and The source excitation is applied on the upper boundary, and the electric fields on other boundaries are obtained by forward modeling of the one-dimensional background model.

[0049] S2: Constructing overlapping sub-region models in unstructured grids.

[0050] The grid partitioning method affects the convergence of Schwarz iteration. Regular grid partitioning helps to accelerate the iterative convergence of the domain decomposition method. However, in earth electromagnetic simulation, local refinement is usually used to improve accuracy, resulting in non-uniform unstructured grids in different subdomains. In order to balance the number of elements in each subdomain, ideal regular partitioning is often not achievable. The commonly used partitioning tool METIS does not consider the geometric shape of the subdomain, and the grid partitioning result is irregular in shape. The present invention proposes an algorithm to achieve more regular grid partitioning. The centroid of the tetrahedron is divided into multiple subdomains by recursively calling K-means clustering to ensure that the number of points contained in each subdomain is less than the specified maximum value. The K-means clustering algorithm considers the spatial distance between elements, thereby achieving more compact and regular grid partitioning.

[0051] More specifically, the process includes:

[0052] (1) Generate computational domain Tetrahedral mesh of ;

[0053] (2) The computational domain The tetrahedron in is divided into non-overlapping subdomains, , ensuring that their union covers the entire computational domain , that is, satisfying: ;

[0054] (3) By extending each subdomain The boundaries of the Each expansion operation adds a layer of tetrahedrons, and the number of expansion operations is used to represent the final overlap width;

[0055] (4) Define each overlapping subdomain through Boolean operations The border , where the subdomain The inner boundary of For subdomain The intersection with the rest, outer boundary The boundary between the subdomain and the overall domain The intersection of;

[0056] (5) Renumber the edges of each subdomain and establish an edge mapping operator between the subdomain and the global domain.

[0057] S3: Vector finite element discretization of overlapping sub-regions.

[0058] The first-order vector finite element method is used to discretize the airborne magnetotelluric forward problem. The Galerkin method is used and the weak integral equation of the overlapping subdomain can be expressed as:

[0059] (3);

[0060] in, represents the first-order vector basis function.

[0061] By using the first-order vector Green's theorem, equation (3) can be written in matrix form:

[0062] (4);

[0063] in, , and They are the stiffness matrix associated with the dual curl operator, the mass matrix associated with the conductivity, and the mass matrix associated with the dielectric constant. The internal definitions are:

[0064] (5);

[0065] (6);

[0066] (7);

[0067] in, and Representative , No. The vector basis functions of the edges, and is The inner edge number.

[0068] Assembling the elements in all overlapping subdomains together, we obtain the following linear equation system:

[0069] (8);

[0070] in, is a size of × The sparse complex matrix of and The size is The vector representing The electric field unknowns and right-hand side terms of the subdomains. Indicates subdomain The number of edges in .

[0071] S4: Set the overlapping sub-region boundary conditions.

[0072] The simplest initial subdomain boundary conditions are constructed by the analytical solution of the one-dimensional model, which matches the boundary conditions of the global problem. The present invention proposes a more efficient method to use a coarse grid solution to solve the global problem, reducing the number of iterations by using better initial conditions. If the number of cells in the coarse grid is still large, the multigrid technique is combined with regional decomposition. First, the regional decomposition algorithm is applied to calculate the solution on the coarse grid, and then the coarse grid solution is used to obtain the initial boundary conditions of the fine grid. Finally, the regional decomposition algorithm is used to calculate the solution on the fine grid.

[0073] If the number of iterations is 1, the boundary conditions are initialized for each subdomain , which can be expressed as:

[0074] (9);

[0075] in, It is The initial amplitudes of the boundaries of the overlapping sub-regions are obtained through one-dimensional forward analytical solution or global coarse grid solution.

[0076] If the number of iterations is greater than or equal to 2, the boundary conditions of the subdomain are updated using the solution of the previous iteration, which can be expressed as:

[0077] (10);

[0078] in, is the boundary of the ith overlapping subregion The solution of the iteration, It is The boundaries of the overlapping sub-regions The solution of the iteration.

[0079] S5: Solve the finite element linear equations for all subdomains.

[0080] Since the direct solver has advantages in stability, accuracy and handling of multiple right-hand side terms, the present invention uses the direct solver to solve the finite element linear equations of the subdomain. In the forward modeling of aeromagnetotelluric, two polarization modes need to be solved, the matrix In the Schwarz iteration process, only the right-hand side remains unchanged. Therefore, a linear system matrix decomposition is performed in the first iteration, and only back-substitution is required to obtain the solution in subsequent iterations.

[0081] Using the direct solver, only one matrix decomposition is required in the Schwarz iteration. The linear system solution can be obtained in just one back iteration, which significantly reduces the computational cost. In contrast, the iterative solver needs to solve Sublinear system, where The invention adopts OpenMP thread parallel computing to solve the linear system with direct solver, thus improving the computing efficiency.

[0082] S6: Update the global solution.

[0083] Since the solution error in the overlapping area is usually large, especially near the internal boundary of the subdomain, this patent uses the solution of the non-overlapping domain to update the global solution, and takes the arithmetic mean of the solution on the overlapping boundary. This strategy reduces the global solution update error caused by the solution error in the overlapping area, thereby accelerating the convergence of the iterative domain decomposition algorithm.

[0084] In this implementation, by combining all non-overlapping subdomains Update the global solution :

[0085] (11);

[0086] in, Indicates The global solution in the iteration, Indicates In the overlapping sub-areas Iterative solution.

[0087] S7: Convergence judgment of Schwarz iteration.

[0088] Calculate the current iteration Compared with the previous iteration The L2 norm of the standardized error between is expressed as:

[0089] (12);

[0090] in, Indicates The global solution in iterations.

[0091] If the error If the value is less than the preset convergence threshold, the iteration is considered to have converged and the solution process is terminated. Otherwise, it returns to S4 and continues to iterate until convergence or the maximum number of iterations is reached; through the Schwarz iteration from S4 to S7, the aeromagnetotelluric forward problem can be quickly solved and an accurate global solution can be gradually obtained in each subdomain.

[0092] S8: Calculate the aeromagnetotelluric response data.

[0093] The electric and magnetic field components at the measuring point are calculated using vector basis functions and Faraday's law. When the measurement point is in a tetrahedral unit, the calculation formulas for the electric field component and the magnetic field component are:

[0094] (13);

[0095] (14);

[0096] in, represents the discrete electric field defined at the edge center, represents the tetrahedron containing the observation point, Local numbers for the edges of the tetrahedron.

[0097] Taking the mode of measuring the magnetic field in the air and the electric field at the ground base station as an example, once the electric field and magnetic field are known, the impedance tensor can be calculated at the observation point by the following equation:

[0098] (15);

[0099] in, , , , , , , , are the horizontal components of the electric field and magnetic field calculated at the measuring point. The subscripts “1” and “2” represent two polarization modes. "and" ” represent ground base stations and aerial receivers respectively, , , and are the components of the impedance tensor respectively.

[0100] Apparent resistivity of airborne magnetotelluric system The calculation is as follows:

[0101] (16);

[0102] in, represents the absolute value of the determinant, is the impedance tensor, is the angular frequency, is the magnetic permeability.

[0103] In order to better illustrate the advantages and purposes of the present invention and verify the correctness and accuracy of the regional decomposition aeromagnetotelluric three-dimensional forward algorithm, this implementation method is verified using the international reference model COMMEMI-3D3 model as an example, and the regional decomposition aeromagnetotelluric three-dimensional forward algorithm of the present invention is compared with the traditional finite element magnetotelluric three-dimensional forward algorithm to verify the accuracy of the algorithm.

[0104] The challenge of the COMMEMI-3D3 model simulation is that the model has a high conductivity contrast. The model contains a set of shallow conductors with resistivities of 30 Ωm and 300 Ωm, and two highly conductive anomalies with resistivities of 0.1 Ωm and 0.3 Ωm. These anomalies are buried in a one-dimensional layered background consisting of three layers with resistivities of 103 Ωm, 104 Ωm and 10 Ωm, respectively, with layer thicknesses of 1 km, 2.5 km and 22.5 km, respectively, and an air layer thickness of 10 km. The survey area is from -2.5 km to 2.5 km in the x and y directions, with a station spacing of 100 m. The flight altitude is 100 m, and the ground base station is located at (-20,0,0) (in km). The model is discretized into two sets of grids: the coarse grid consists of 86575 tetrahedra and 102461 edges, and the fine grid consists of 719349 tetrahedra and 846197 edges. The conductivity distribution of the model and the fine grid are shown in Figure 2. Figure 2 shown.

[0105] The maximum number of tetrahedrons in each subdomain of the fine grid is set to 20,000. The fine grid model is divided into 59 subdomains, while the coarse grid model is divided into 6 subdomains. Figure 3 shown.

[0106] In order to evaluate the effects of different overlap widths (number of layers) and initial boundary conditions on the accuracy and efficiency of the algorithm, six parameter settings were considered in the COMMEMI-3D3 model, with overlap widths of 1, 2, and 3 layers, and initial boundary conditions derived from the analytical solution of the one-dimensional layered model. and the 3D global coarse grid solution .

[0107] Figure 4 The convergence of the Schwarz iteration of the fine grid model under different overlap widths and initial boundary conditions is shown. Table 1 lists the calculation statistics of different methods. Taking the traditional finite element method as a reference, the normalized L2 norm error of the apparent resistivity at all observation points is defined as ,in and They represent the results obtained by the regional decomposition algorithm and directly solving the global finite element forward equations, respectively, and are used to measure the simulation accuracy of the regional decomposition algorithm.

[0108] Table 1: Statistics of the COMMEMI-3D3 model

[0109]

[0110] From Table 1 and Figure 4 The following results can be observed: (1) The peak memory usage of the domain decomposition algorithm is significantly lower than that of the traditional finite element method. The computational time of the domain decomposition algorithm is also shorter than that of the traditional finite element method; (2) Increasing the overlap width (number of layers) helps the Schwarz iteration converge, but it also increases the degrees of freedom and computational scale of the subdomain, resulting in a rapid increase in computational time and memory usage. Excessive overlap width may cause the domain decomposition algorithm to lose its computational performance advantage; (3) The introduction of the multigrid scheme significantly improves the convergence speed of the Schwarz iteration. The coarse grid solution provides an effective initial boundary condition for the fine grid simulation, thereby reducing the number of required iterations. This method reduces the computational time, and the modification of the boundary conditions does not increase memory consumption; (4) After reaching the convergence threshold, the accuracy of the domain decomposition algorithm is comparable to that of the traditional finite element method, and the normalized L2 norm error is less than 0.56%.

[0111] The present invention plots the apparent resistivity results calculated by FE-DDM under different numbers of overlapping layers and initial boundary conditions and compares them with traditional FE simulations. Figure 5 As shown, is the analytical solution of the one-dimensional layered model, The three-dimensional global coarse grid solution is given with 1 layer, 2 layers and 3 layers respectively. The results accurately identify the conductors with a maximum relative error of less than 4%. This shows that the FE-DDM algorithm can effectively simulate the three-dimensional airborne magnetotelluric response with high conductivity contrast under different initial boundary conditions and overlap widths. In summary, considering the computational time and memory consumption, the Schwarz regional decomposition algorithm can be used to use a 1-layer overlapping tetrahedral grid combined with a global coarse grid solution as the initial boundary for simulation.

[0112] Embodiment 2:

[0113] like Figure 6 As shown, this implementation provides a regional decomposition magnetotelluric three-dimensional forward modeling system, including:

[0114] The overlapping sub-region generating unit is configured to: generate overlapping sub-regions in the unstructured grid of the computational domain of the forward simulation;

[0115] The vector finite element discretization unit is configured to: perform vector finite element discretization on the overlapping sub-regions to obtain a finite element linear equation group for each overlapping sub-region;

[0116] The finite element linear equation solving unit is configured to: set the boundary conditions of the overlapping sub-regions according to the vector finite element discretization results, and solve the finite element linear equations of all overlapping sub-regions;

[0117] The global solution updating unit is configured to: update the global solution by combining the solutions in all overlapping sub-regions;

[0118] The forward calculation unit is configured to: determine whether the convergence condition is reached according to the global solution, if the convergence condition is reached and when the measuring point is at the e When the tetrahedral unit is in the grid, the electric field component and the magnetic field component at the measuring point are calculated to obtain the impedance tensor at the observation point and the apparent resistivity of the airborne magnetotelluric system; if the convergence condition is not met, the boundary condition of the overlapping sub-region is reset and the iteration is continued until the convergence condition is met.

[0119] It is understandable that the above-mentioned units can be separately or completely combined into one or several other units to constitute, or one (some) of the units can be further divided into multiple smaller units in function to constitute, which can achieve the same operation without affecting the realization of the technical effects of the embodiments of the present application. The above-mentioned units are divided based on logical functions. In practical applications, the functions of one unit can also be implemented by multiple units, or the functions of multiple units can be implemented by one unit. In other embodiments of the present application, the system may also include other units. In practical applications, these functions can also be implemented with the assistance of other units, and can be implemented by the collaboration of multiple units.

[0120] According to another embodiment of the present application, the system described in this embodiment can be constructed, and the method of Example 1 of the present application can be implemented by running a computer program (including program code) capable of executing the steps involved in the corresponding method described in Example 1 on a general-purpose computing device such as a computer, which includes processing elements and storage elements such as a central processing unit (CPU), random access memory (RAM), and read-only memory (ROM). The computer program can be recorded on, for example, a computer-readable recording medium, and loaded into the above-mentioned computing device through the computer-readable recording medium and run therein.

[0121] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A regional decomposition magnetotelluric three-dimensional forward modeling method, characterized in that: The process includes: Generate overlapping sub-regions in the unstructured grid of the computational domain for forward simulation according to each tetrahedral unit in the tetrahedral grid of the computational domain; Performing vector finite element discretization on the overlapping sub-regions to obtain a finite element linear equation group for each overlapping sub-region; According to the vector finite element discretization results, the boundary conditions of the overlapping sub-regions are set, and the finite element linear equations of all overlapping sub-regions are solved; Update the global solution by combining the solutions in all overlapping subregions; Judging whether the convergence condition is reached according to the global solution, if the convergence condition is reached and when the measuring point is in a tetrahedral unit, calculating the electric field component and the magnetic field component at the measuring point, and obtaining the impedance tensor at the observation point and the apparent resistivity of the airborne magnetotelluric system; If the convergence condition is not met, reset the boundary conditions of the overlapping sub-regions and continue iterating until the convergence condition is met; Set overlapping sub-region boundary conditions, including: If the number of iterations is 1, the boundary conditions are initialized for each subdomain , expressed as: ； in, It is The initial amplitude of the boundaries of the overlapping sub-regions is obtained by one-dimensional forward analytical solution or global coarse grid solution; If the number of iterations is greater than or equal to 2, the boundary conditions of the subdomain are updated using the solution of the previous iteration, expressed as: ； in, It is The boundaries of the overlapping sub-regions The solution of the iteration, It is The boundaries of the overlapping sub-regions The solution of the iteration.

2. The regional decomposition magnetotelluric three-dimensional forward modeling method according to claim 1, characterized in that: Based on the tetrahedral elements in the tetrahedral mesh of the computational domain, the overlapping sub-regions in the unstructured mesh of the computational domain for forward simulation are generated, including: Generate Computed Domain Tetrahedral mesh of Compute the domain The tetrahedron in is divided into non-overlapping subdomains , the union of all non-overlapping subdomains covers the computational domain ,satisfy: ; By extending each subdomain The boundaries of the , each expansion operation will add a layer of tetrahedron, and the number of expansion operations is used to represent the final overlap width; Each overlapping sub-region is defined by a Boolean operation The border , where the subdomain The inner boundary of For subdomain The intersection with the rest, outer boundary The boundary between the subdomain and the overall domain The intersection of; The edges of each subdomain are renumbered, and an edge mapping operator is established between the subdomain and the global domain.

3. The regional decomposition magnetotelluric three-dimensional forward modeling method according to claim 1, characterized in that: The overlapping sub-regions are discretized by vector finite element to obtain a finite element linear equation group of each overlapping sub-region, including: The first-order vector finite element method is used to discretize the airborne magnetotelluric forward problem. The Galerkin method is used and the weak integral equation of the overlapping sub-region is: ； in, represents the first-order vector basis function, is the angular frequency, is the frequency, parameter , and represent the electrical conductivity, magnetic permeability and dielectric constant respectively, is the electric field amplitude on the edge, represents the curl operator; By using the first-order vector Green's theorem, written in matrix form: ； in, , and They are the stiffness matrix associated with the dual curl operator, the mass matrix associated with the conductivity, and the mass matrix associated with the dielectric constant; Assembling the elements in all overlapping subregions together, we obtain the following system of linear equations: ； in, is a size of The sparse complex matrix of and The sizes are The vector representing The unknown electric field and the right-hand side term of the subdomain, the variable Indicates subdomain The number of edges in .

4. The regional decomposition magnetotelluric three-dimensional forward modeling method according to claim 3, characterized in that: Stiffness matrix associated with the double curl operator , the mass matrix related to conductivity and the mass matrix of the dielectric constant term The elements in the given tetrahedron The internal definitions are , and : ； ； ； in, and Representative , No. The vector basis functions of the edges, and In the tetrahedron The inner edge number.

5. The regional decomposition magnetotelluric three-dimensional forward modeling method according to claim 1, characterized in that: Solve the finite element linear equations for all overlapping subdomains, including: Using a direct solver to solve the finite element linear equations for the subdomains.

6. The regional decomposition magnetotelluric three-dimensional forward modeling method according to claim 1, characterized in that: Judging whether the convergence condition is reached according to the global solution includes: Calculate the current iteration Compared with the previous iteration The L2 norm of the standardized error between , if the error If the value is less than the preset convergence threshold, the iteration is considered to have converged and reached the convergence condition; otherwise, the iteration continues until the convergence condition is reached or the maximum number of iterations is reached.

7. The regional decomposition magnetotelluric three-dimensional forward modeling method according to claim 1, characterized in that: When the measuring point is When there are four tetrahedral units, the electric field component and magnetic field component of the measuring point are: ； ； in, represents the discrete electric field defined at the edge center, represents the tetrahedron containing the observation point, is the local number of the tetrahedron edges, represents the curl operator.

8. The regional decomposition magnetotelluric three-dimensional forward modeling method according to claim 1, characterized in that: Impedance tensor at the observation point ,include: ； in, and represent ground base stations and aerial receivers respectively, , , , , , , , are the horizontal components of the electric field and magnetic field calculated at the measuring point, , , and are the components of the impedance tensor respectively; Apparent resistivity of airborne magnetotelluric system ,include: ； represents the absolute value of the determinant, is the angular frequency, is the magnetic permeability.

9. A regional decomposition magnetotelluric three-dimensional forward modeling system, characterized in that: include: The overlapping sub-region generating unit is configured to: generate overlapping sub-regions in the unstructured grid of the computational domain of the forward simulation; The vector finite element discretization unit is configured to: perform vector finite element discretization on the overlapping sub-regions to obtain a finite element linear equation group for each overlapping sub-region; The finite element linear equation solving unit is configured to: set the boundary conditions of the overlapping sub-regions according to the vector finite element discretization results, and solve the finite element linear equations of all overlapping sub-regions; The global solution updating unit is configured to: update the global solution by combining the solutions in all overlapping sub-regions; The forward calculation unit is configured to: determine whether the convergence condition is reached according to the global solution, if the convergence condition is reached and when the measuring point is in a tetrahedral unit, calculate the electric field component and the magnetic field component at the measuring point to obtain the impedance tensor at the observation point and the apparent resistivity of the airborne magnetotelluric system; If the convergence condition is not met, reset the boundary conditions of the overlapping sub-regions and continue iterating until the convergence condition is met; Set overlapping sub-region boundary conditions, including: If the number of iterations is 1, the boundary conditions are initialized for each subdomain , expressed as: ； in, It is The initial amplitude of the boundaries of the overlapping sub-regions is obtained by one-dimensional forward analytical solution or global coarse grid solution; If the number of iterations is greater than or equal to 2, the boundary conditions of the subdomain are updated using the solution of the previous iteration, expressed as: ； in, It is The boundaries of the overlapping sub-regions The solution of the iteration, It is The boundaries of the overlapping sub-regions The solution of the iteration.

Images