Self-adaptive three-dimensional fast inversion imaging method for magnetotelluric data

Abstract

The invention discloses a self-adaptive three-dimensional fast inversion imaging method for magnetotelluric data. The method comprises the following steps of: subdividing a model by adopting a non-structural tetrahedral mesh; calculating forward modeling response of the model by adopting a self-adaptive non-structural finite element method of coulomb specification potential; carrying out adaptive encryption on the forward modeling grid; and carrying out inversion by adopting a finite memory quasi-Newton method (L-BFGS) to obtain an inversion imaging result. According to the method, the separated forward and reverse adaptive grid system is adopted, the forward precision is ensured, meanwhile, parameters of an inversion grid do not need to be increased, the calculation efficiency and the calculation precision are improved, the MPI strategy is introduced to conduct parallel operation on the frequency, and the calculation efficiency is further improved; by adopting the non-structural tetrahedral mesh generation, the applicability of the three-dimensional inversion imaging of the magnetotelluric data in a complex fluctuating region is improved, and the method has important practical value in the aspects of geomagnetic data inversion imaging of deep ground exploration and complex rugged mineral resource exploration and the like.

Description

Technical Field

[0001] The present invention relates to an adaptive three-dimensional rapid inversion imaging method of magnetotelluric data, and belongs to the field of geomagnetic measurement for deep earth detection and mineral resource exploration in complex mountainous areas. Background Art

[0002] Magnetotellurics (MT) is a classic geophysical method for studying the distribution of electrical conductivity within the Earth's interior. Thanks to advances in computing technology, significant progress has been made in the three-dimensional inversion of MT data. The discretization of the model grid used in 3D inversion directly determines the quality of the inversion results and computational efficiency. Therefore, improving the quality of the forward and inversion grids is a key component.

[0003] In recent years, domestic and foreign scholars have done a lot of research on the optimization design of grids in forward and inversion. Among them, the most research is on the grid refinement problem in forward modeling, mainly grid adaptive refinement based on posterior error estimation and grid refinement based on prior information constraints, while there are fewer studies on grid refinement in inversion.

[0004] Unstructured grids can flexibly divide complex geometric models and are increasingly used in the forward and inversion of geophysical electromagnetic data. Usually, a single grid is used based on the user's experience. In practical applications, this limits the flexibility of the grid division and reduces its computational efficiency. Therefore, unstructured grid inversion technology has the problems of low computational efficiency in three-dimensional magnetotelluric inversion and poor adaptability to complex structures and undulating terrain areas. Summary of the Invention

[0005] The object of the present invention is to provide an adaptive three-dimensional rapid inversion imaging method for magnetotelluric data, which can improve the flexibility and computational efficiency of unstructured grid subdivision and improve the applicability to complex structures and undulating terrain areas.

[0006] To achieve the above object, the present invention provides an adaptive three-dimensional rapid inversion imaging method for magnetotelluric data, comprising the following steps:

[0007] S1. Use unstructured tetrahedral mesh to divide the model;

[0008] S2. Calculate the forward response of the model using the unstructured finite element method using Coulomb gauge potential;

[0009] S3, adaptively encrypt the forward modeling grid;

[0010] S4. Use the limited-memory quasi-Newton method (L-BFGS) to invert and obtain the inversion imaging result.

[0011] Furthermore, the specific process of S1 is as follows:

[0012] S1.1. First, create a poly file for meshing, which mainly includes the coordinates of the control points, line segments, surfaces and other elements of the geometric model, the number of each area and the volume constraint value;

[0013] S1.2. Use the unstructured mesh generation software TetGen to input the poly file and generate four files: node, element, neigh, and face. The node file is the node list file of the unstructured mesh, which contains the numbers and three-dimensional coordinates (x, y, z) of all the nodes in the mesh. The element file is the tetrahedron list file, which contains the node list and attribute numbers corresponding to each tetrahedron. The neigh file contains the numbers of the tetrahedrons adjacent to each tetrahedron. The face file is the list file of the triangular faces contained in the mesh.

[0014] S1.3. Renumber the generated tetrahedron unit attributes. The tetrahedron attributes in the inversion area are numbered in sequence (1, 2, 3, ...). That is, each tetrahedron is an inversion parameter. Generate a new tetrahedron file for subsequent forward modeling and inversion.

[0015] S1.4. Assign an initial resistivity value to each area of ​​the grid and generate a resistivity file, in which the inversion area is marked with a number greater than 0, and the non-inversion area is marked with 0.

[0016] Furthermore, the specific process of S2 is as follows:

[0017] S2.1. The differential equations satisfied by the electric field (E) and the magnetic field (H) are as follows:

[0018]

[0019] Where i is the imaginary unit, ω is the angular frequency, μ0 is the magnetic permeability of the vacuum medium, and σ is the electrical conductivity;

[0020] S2.2. According to the definition of Coulomb potential, the electric field (E) and magnetic field (H) are expressed as a magnetic vector potential A and an electric scalar potential φ, as follows:

[0021]

[0022] S2.3. Use Coulomb norm conditions: Substituting equations (3) and (4) into equations (1) and (2), we obtain the Helmholtz equation for the A-φ potential:

[0023]

[0024] S2.4. Divide the model area into tetrahedral elements, interpolate using the tetrahedral element shape function, integrate each element, and synthesize the integrals of each element to obtain the overall linear equation system:

[0025] Ku=0(7);

[0026] Where K is the stiffness matrix, u is the A, φ potential to be solved;

[0027] S2.5. Place the boundary conditions into the linear equation system (7), select the linear equation system solver Pardiso to solve, obtain the finite element solution of A and φ potential, and then obtain the finite element solution of the electromagnetic field according to formulas (3) and (4).

[0028] Furthermore, the specific process of S3 is as follows:

[0029] S3.1. Densify the cells near the receiving point, find the cell where the receiving point is located, and determine whether its volume is less than the set receiving point cell volume requirement. If not, densify it in the subsequent mesh refinement. Repeat the above steps until the cell volume meets the requirement.

[0030] S3.2. Loop through each cell of the grid and find the minimum distance from the cell to all receiving points. Based on the distance and skin depth, calculate the maximum volume that the cell should satisfy. The calculation formula is as follows:

[0031]

[0032] Among them, δ e is the skin depth, λ is the control factor, and r is the minimum distance from the unit to the receiving point. If the unit volume exceeds the maximum volume constraint, it will be encrypted in the subsequent mesh refinement. Repeat the above steps until the volume of all units meets the requirements.

[0033] S3.3. Since an independent forward and inversion dual-grid system is used, each unit that needs to be inverted in the inversion grid is marked as an independent unit. After the grid is refined, the parameters of each refined unit in the inversion grid are transferred to the refined sub-units of the unit.

[0034] Furthermore, the specific process of S4 is as follows:

[0035] S4.1. The inversion objective function is constructed by the Tikhonov regularization method. The objective functional Φ is expressed as:

[0036] Φ(m)=μ||R(m-m0)|| 2 +||W(dF(m))|| 2 (9);

[0037] Where m is the N-dimensional model parameter vector, generally the logarithmic resistivity value; m0 is the prior or initial model parameter vector; R is the roughness operator matrix; μ is the Lagrange multiplier used to balance the model constraint terms and data fitting terms; W is the data covariance matrix, a diagonal matrix composed of data error terms; d is the observed data vector; F(m) is the forward response corresponding to model m;

[0038] S4.2. Obtain the gradient of the objective function (9) and obtain:

[0039] g(m)=-2J T W T W(dF(m))+2μR T R(m-m0) (10);

[0040] Where J is the Jacobian matrix, which is obtained by the adjoint equation method during forward modeling;

[0041] S4.3. Solve the objective function through iteration and continuously update the model vector to minimize the objective function. The model update iteration form is expressed as:

[0042] m k+1 =m k +α k p k (11);

[0043] Among them, α k is the search step size, which controls the correction size of the model, p k It is the search direction and controls the correction direction of the model.

[0044] Furthermore, S2 and S3 also include parallel calculation of the forward response and Jacobian matrix of the model based on the MPI strategy. The calculation process is as follows: the total number of groups is calculated according to the frequency, one group for each frequency, and numbered in sequence; the main process distributes the calculation task groups to the sub-processes in sequence. If the number of sub-processes is insufficient, it waits for the remaining sub-processes to complete the operation before continuing the distribution. This process is completed through message passing between the main process and the sub-processes; after receiving the sub-task, each sub-process adaptively encrypts the forward grid of the task group and calls the PARDISO solver to calculate the forward response of the model, and then sends the calculation results to the main process; until all sub-tasks are completed, the main process obtains all the calculation results, outputs the model response and Jacobian matrix calculation results, releases all process states, and ends MPI.

[0045] This invention utilizes a separate forward and inverse adaptive grid system to ensure forward modeling accuracy without increasing inversion grid parameters, thereby improving computational efficiency and accuracy. The introduction of an MPI strategy for frequency parallelization further enhances computational efficiency. The use of unstructured tetrahedral meshing improves the applicability of three-dimensional inversion imaging of magnetotelluric data in complex, rugged areas. This invention enhances the flexibility and computational efficiency of unstructured meshing, as well as its applicability to complex structures and rugged terrain. It has significant practical value in deep earth exploration and geomagnetic surveying for complex and rugged mineral resource exploration. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 It is a schematic diagram of the workflow of the method of the present invention;

[0047] Figure 2 is a schematic diagram of a terrain model in an embodiment of the present invention;

[0048] Figure 3 is a schematic diagram of a terrain model inversion grid in an embodiment of the present invention;

[0049] Figure 4 is a schematic diagram of a forward grid of a terrain model in an embodiment of the present invention;

[0050] Figure 5 is a comparison diagram of the terrain model inversion result and the real model in an embodiment of the present invention;

[0051] Figure 6 3 is a graph showing the root mean square fitting RMS and the objective function value in an embodiment of the present invention. DETAILED DESCRIPTION

[0052] The present invention will be further described below with reference to the accompanying drawings.

[0053] like Figure 1 As shown, an adaptive three-dimensional rapid inversion imaging method for magnetotelluric data includes the following steps:

[0054] S1. Use unstructured tetrahedral mesh to divide the model;

[0055] S2. Calculate the forward response of the model using the unstructured finite element method using Coulomb gauge potential;

[0056] S3, adaptively encrypt the forward modeling grid;

[0057] S4. Use the limited-memory quasi-Newton method (L-BFGS) to invert and obtain the inversion imaging result.

[0058] Example: In order to test the feasibility of the imaging method, Figure 2The inversion calculation is performed on the model of multiple anomalies with undulating terrain shown in the figure. The model is a low-resistance anomaly and a high-resistance anomaly under an inclined terrain. The resistivity of the low-resistance anomaly is 10Ωm, the resistivity of the high-resistance anomaly is 3000Ωm, and the resistivity of the surrounding rock is 100Ωm. The geometric dimensions of the two anomalies are 8km×8km×4km. Figure 2 (a) shows the plane position of the anomaly in the model and the distribution of magnetotelluric measurement points. The magnetotelluric measurement point and line distances are both 2 km, with a total of 7 measurement lines and 13 measurement points per line. The measurement point distribution completely covers the range of the two anomalies, totaling 91 measurement points. Figure 2 (b) shows a vertical depth slice of the model. The slope terrain rises 1 km from west to east, and is distributed between -3 km and 3 km east-west. The top of the low-resistance anomaly is buried at a depth of 2 km, and the top of the high-resistance anomaly is buried at a depth of 3 km. Magnetotelluric observations were conducted at 10 frequencies, ranging from 4 Hz to 0.01 Hz. 3% Gaussian random noise was added to the model response data to simulate measured data.

[0059] The model is discretized using unstructured tetrahedral mesh. The model inversion mesh is as follows Figure 3 As shown, it contains 42639 nodes and 259558 units, of which the inversion parameters are 145525. Since an independent forward and inverse grid system is used, the forward grid is obtained by adaptively encrypting the inversion grid. The forward grid of this embodiment is as follows Figure 4 As shown, Figure 4 (a) and (b) in the figure correspond to the forward modeling grids of 4 Hz and 0.14 Hz respectively. The 4 Hz grid contains 63,066 nodes and 369,043 core elements after refinement, while the 0.14 Hz grid contains 45,577 nodes and 272,708 core elements after refinement. Figure 4 It can be seen that the degree of mesh encryption for forward modeling at different frequencies is inconsistent. When the frequency is low, the mesh is sparse, and interpolation encryption is mainly performed near the measurement point. When the frequency is high, the mesh is dense, which is consistent with the propagation law of electromagnetic waves.

[0060] Figure 5 is the comparison between the final inversion result of the undulating terrain model and the real model, where Figure 5 (a) in the figure is the real model. Figure 5 (b) in the figure is the final inversion result. As can be seen from the figure, the inversion results have restored the spatial position and morphology of the two anomalies under the undulating terrain well, especially the resistivity and morphology of the low-resistance anomaly are closer to the real model, indicating that the inversion result is accurate and reliable. This embodiment shows that the present invention is accurate and feasible, and can perform inversion imaging of the target body under the undulating terrain. Due to the use of unstructured grid subdivision, it is conducive to the subdivision of complex undulating terrain models, and can save a large number of grid cells relative to regular grids;

[0061] like Figure 6 The figure shows the root mean square (RMS) fitting and objective function of the data during the model inversion process. As can be seen from the figure, after 47 iterations, the RMS decreases from 2.91 to within 0.92, and the objective function value decreases from 6.92×104 to 3.09×103. Both curves converge stably, indicating that the inversion imaging results can stably converge to the true model. Because the number of observation frequencies is 10 and 11 threads are used for parallel computing, the inversion time required is 4.1 hours. This example verifies the high computational efficiency of the present invention.

Claims

1. An adaptive three-dimensional rapid inversion imaging method for magnetotelluric data, characterized in that: The steps include: S1. Use unstructured tetrahedral mesh to divide the model; S2. Calculate the forward response of the model using the unstructured finite element method using Coulomb gauge potential; S3, adaptively encrypt the forward modeling grid; S4. Use the limited memory quasi-Newton method L-BFGS inversion to obtain the inversion imaging result.

2. The adaptive three-dimensional rapid inversion imaging method of magnetotelluric data according to claim 1, characterized in that: The specific process of S1 is: S1.

1. Create a poly file for meshing, including the coordinates of the control points, line segments, surface elements of the geometric model, the numbers of each area, and the volume constraint values; S1.

2. Use the unstructured mesh generation software TetGen to input the poly file and generate four files: node, element, neigh, and face. The node file is the node list file of the unstructured mesh, which contains the numbers and three-dimensional coordinates (x, y, z) of all the nodes in the mesh. The element file is the tetrahedron list file, which contains the node list and attribute numbers corresponding to each tetrahedron. The neigh file is the number of the tetrahedrons adjacent to each tetrahedron. The face file is the list file of the triangular faces contained in the mesh. S1.

3. Renumber the generated tetrahedron unit attributes. The tetrahedron attributes in the inversion area are numbered sequentially, that is, each tetrahedron is an inversion parameter. A new tetrahedron file is generated for subsequent forward modeling and inversion. S1.

4. Assign an initial resistivity value to each area of ​​the grid and generate a resistivity file, in which the inversion area is marked with a number greater than 0, and the non-inversion area is marked with 0.

3. The adaptive three-dimensional rapid inversion imaging method of magnetotelluric data according to claim 1, characterized in that: The specific process of S2 is: S2.

1. The differential equations satisfied by the electric field E and the magnetic field H are as follows: Where i is the imaginary unit, ω is the angular frequency, μ0 is the magnetic permeability of the vacuum medium, and σ is the electrical conductivity; S2.

2. According to the definition of Coulomb potential, the electric field E and magnetic field H are expressed as a magnetic vector potential A and an electric scalar potential φ, as follows: S2.

3. Use Coulomb norm conditions: Substituting equations (3) and (4) into equations (1) and (2), we obtain the Helmholtz equation for the A-φ potential: S2.

4. Divide the model area into tetrahedral elements, interpolate using the tetrahedral element shape function, integrate each element, and synthesize the integrals of each element to obtain the overall linear equation system: Ku=0(7); Where K is the stiffness matrix, u is the A, φ potential to be solved; S2.

5. Place the boundary conditions into the linear equation system (7), select the linear equation system solver Pardiso to solve, obtain the finite element solution of A and φ potential, and then obtain the finite element solution of the electromagnetic field according to formulas (3) and (4).

4. The adaptive three-dimensional rapid inversion imaging method of magnetotelluric data according to claim 1, characterized in that: The specific process of S3 is as follows: S3.

1. Densify the cells near the receiving point, find the cell where the receiving point is located, and determine whether its volume is less than the set receiving point cell volume requirement. If not, densify it in the subsequent mesh refinement. Repeat the above steps until the cell volume meets the requirement. S3.

2. Loop through each cell of the grid and find the minimum distance from the cell to all receiving points. Based on the distance and skin depth, calculate the maximum volume that the cell should satisfy. The calculation formula is as follows: Among them, δ e is the skin depth, λ is the control factor, and r is the minimum distance from the unit to the receiving point. If the unit volume exceeds the maximum volume constraint, it will be encrypted in the subsequent mesh refinement. Repeat the above steps until the volume of all units meets the requirements. S3.

3. Since an independent forward and inversion dual-grid system is used, each unit that needs to be inverted in the inversion grid is marked as an independent unit. After the grid is refined, the parameters of each refined unit in the inversion grid are transferred to the refined sub-units of the unit.

5. The adaptive three-dimensional rapid inversion imaging method of magnetotelluric data according to claim 1, characterized in that: The specific process of S4 is: S4.

1. The inversion objective function is constructed by the Tikhonov regularization method. The objective functional Φ is expressed as: Φ(m)=μ||R(m-m0)|| 2 +||W(d-F(m))|| 2 (9); Where m is the N-dimensional model parameter vector, generally the logarithmic resistivity value; m0 is the prior or initial model parameter vector; R is the roughness operator matrix; μ is the Lagrange multiplier used to balance the model constraint terms and data fitting terms; W is the data covariance matrix, a diagonal matrix composed of data error terms; d is the observed data vector; F(m) is the forward response corresponding to model m; S4.

2. Obtain the gradient of the objective function (9) and obtain: g(m)=-2J T W T W(d-F(m))+2μR T R(m-m0) (10); Where J is the Jacobian matrix, which is obtained by the adjoint equation method during forward modeling; S4.

3. Solve the objective function through iteration and continuously update the model vector to minimize the objective function. The model update iteration form is expressed as: m k+1 =m k +α k p k (11); Among them, α k is the search step size, which controls the correction size of the model, p k It is the search direction and controls the correction direction of the model.

6. The adaptive three-dimensional rapid inversion imaging method of magnetotelluric data according to claim 1, characterized in that: Said S2 and S3 also include parallel calculation of the forward response and Jacobian matrix of the model based on the MPI strategy, and the calculation process is as follows: the total number of groups is calculated according to the frequency, one frequency is one group, and they are numbered in sequence; The main process groups the computational tasks and assigns them to the subprocesses in sequence. If the number of subprocesses is insufficient, it waits for the remaining subprocesses to complete before continuing to assign tasks. This process is accomplished through message passing between the main process and the subprocesses. After receiving a subtask, each subprocess adaptively encrypts the forward modeling grid of the task group and calls the PARDISO solver to calculate the forward response of the model. The calculation results are then sent to the main process. Until all subtasks are completed, the main process obtains all calculation results, outputs the model response and Jacobian matrix calculation results, releases all process states, and MPI ends.

Images