Adaptive 3D fast inversion imaging method for magnetotelluric data

Abstract

The application discloses a self-adapting three-dimensional fast inversion imaging method for magnetotelluric data, which adopts a non-structural tetrahedral grid to divide a model; adopts an adaptive non-structural finite element method of a Coulomb standard potential to calculate a forward response of the model; performs adaptive encryption on the forward grid; and adopts a limited memory quasi-Newton method (L-BFGS) to perform inversion, so as to obtain an inversion imaging result. The application adopts a separated forward and inversion adaptive grid system, ensures the forward precision, does not need to increase the parameters of the inversion grid, improves the calculation efficiency and the calculation precision, introduces an MPI strategy to perform parallel operation on the frequency, further improves the calculation efficiency, adopts the non-structural tetrahedral grid division, improves the applicability of the three-dimensional inversion imaging of the magnetotelluric data in a complex undulating area, and has important practical values in deep earth exploration and magnetotelluric data inversion imaging for complex rugged mineral resources exploration.

Description

TECHNICAL FIELD

[0001] The application relates to a self-adaptive three-dimensional fast inversion imaging method for magnetotelluric data, and belongs to the field of magnetotelluric measurement for deep earth exploration and mineral resource exploration in complex mountainous areas. BACKGROUND

[0002] Magnetotelluric method is a classical geophysical method for studying the conductivity distribution of the earth interior. Thanks to the development of computing technology, three-dimensional inversion of magnetotelluric data has made great progress. The quality and computational efficiency of the inversion result are directly determined by the model grid discretization used in three-dimensional inversion, so improving the quality of the forward and inversion grid is a key link.

[0003] In recent years, domestic and foreign scholars have done a lot of research on the optimization design of the grid in forward and inversion, among which the refinement of the grid in forward is mainly studied, mainly including grid self-adaptive refinement based on posterior error estimation and grid refinement based on prior information constraint, and the refinement of the grid in inversion is less studied.

[0004] Non-structured grid can be flexibly divided into complex geometric models and is increasingly used in forward and inversion of geophysical electromagnetic data. A single grid is usually divided according to the user's experience, which limits the flexibility of grid division and reduces the computational efficiency in practical application, so the non-structured grid inversion technology has problems of low computational efficiency and poor adaptability to complex structures and undulating terrain areas in three-dimensional magnetotelluric inversion. SUMMARY

[0005] The application aims to provide a self-adaptive three-dimensional fast inversion imaging method for magnetotelluric data, which can improve the flexibility and computational efficiency of non-structured grid division and improve the applicability to complex structures and undulating terrain areas.

[0006] In order to achieve the above-mentioned purpose, the application provides a self-adaptive three-dimensional fast inversion imaging method for magnetotelluric data, which comprises the following steps:

[0007] S1, dividing the model by using a non-structured tetrahedral grid;

[0008] S2, calculating the forward response of the model by using a non-structured finite element method of Coulomb standard potential;

[0009] S3, adaptively encrypting the forward grid;

[0010] S4, performing inversion by using a limited memory quasi-Newton method (L-BFGS) to obtain an inversion imaging result.

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

[0012] S1.1, first establish the poly file of mesh partition, mainly including the control point coordinates of geometric model, line segment, face and other elements, the numbering of each region and the volume constraint value;

[0013] S1.2, using non-structured mesh partition software TetGen, input poly file, generate node, element, neigh, face four files, node file is the node list file of non-structured grid, that is, it contains the numbering and three-dimensional coordinates (x, y, z) of all nodes of the grid, element file is tetrahedron list file, which contains the node list and attribute number corresponding to each tetrahedron, neigh file is the tetrahedron number adjacent to each four face, face file is the list file of triangular face contained in the grid;

[0014] S1.3, renumber the generated tetrahedral element attribute, the tetrahedral attribute in the inversion region is numbered in order (1, 2, 3, …), that is, each tetrahedron is an inversion parameter, a new tetrahedral file is generated, which is used for subsequent forward and inversion;

[0015] S1.4, assign resistivity initial value to each region of the grid, generate resistivity file, where the inversion region is marked with a number greater than 0, and the non-inversion region is marked as 0.

[0016] Further, the specific process of S2 is:

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

[0018]

[0019] Where i is the imaginary unit, ω is the angular frequency, μ0 is the vacuum dielectric permeability, and σ is the 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 φ, the formula is as follows:

[0021]

[0022] S2.3, using Coulomb gauge condition: Substitute formula (3) and (4) into formula (1) and (2) to obtain Helmholtz equation about A-φ potential:

[0023]

[0024] S2.4, the model area is divided into tetrahedral elements, interpolation is carried out by using tetrahedral element shape function, each element is integrated, and each element integral is synthesized to obtain a total linear equation group:

[0025] Ku=0(7);

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

[0027] S2.5, the boundary conditions are put into the linear equation group (7), a linear equation group solver Pardiso is selected to solve, the finite element solution of A, φ potential is obtained, and the finite element solution of electromagnetic field is obtained according to formula (3) and (4).

[0028] Further, the specific process of S3 is:

[0029] S3.1, the elements near the receiving point are encrypted, the element where the receiving point is located is found, whether the volume meets the requirement of the set receiving point element volume is judged, if not, the encryption is carried out in the subsequent grid refinement; the above steps are repeated until the volume of the element meets the requirement;

[0030] S3.2, each element of the grid is cycled, the minimum value of the distance from the element to all receiving points is calculated, and the maximum volume that the element should meet is calculated according to the distance size and the skin depth size, and the calculation formula is as follows:

[0031]

[0032] Wherein, δ e is the skin depth, λ is the control factor, and r is the minimum distance from the element to the receiving point; if the volume of the element exceeds the maximum volume constraint, the encryption is carried out in the subsequent grid refinement; the above steps are repeated until the volume of all elements meets the requirement;

[0033] S3.3, since an independent forward and inverse double grid system is adopted, each element that needs to be inverted in the inversion grid is marked as an independent element, and after the grid refinement, the parameters of each refined element in the inversion grid are transmitted to each sub-element after the refinement of the element.

[0034] Further, the specific process of S4 is:

[0035] S4.1, the inversion objective function is constructed by using the Tikhonov regularization method, and the objective functional Φ is represented as:

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

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

[0038] S4.2, the gradient of the objective function (9) is obtained:

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

[0040] wherein, J is a Jacobian matrix, obtained by adjoint equation method during forward modeling;

[0041] S4.3, the objective function is solved by iteration method, and the model vector is updated constantly, so that the objective function reaches a minimum, and the model update iteration form is represented as:

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

[0043] wherein, α k is a search step, controlling the correction size of the model, and p k is a search direction, controlling the correction direction of the model.

[0044] Further, the S2 and S3 further comprise parallel computing of the forward response of the model and the Jacobian matrix based on the MPI strategy, and the computing process is as follows: the total number of groups is calculated according to the frequency, one group for one frequency, and numbered in turn; the main process distributes the computing tasks to the sub-processes in groups in turn, and if the number of sub-processes is insufficient, it waits for the remaining sub-processes to run and then continues to distribute, and this process is completed through message passing between the main process and the sub-processes; after receiving the sub-tasks, each sub-process performs adaptive encryption on 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 calculation results of the model response and the Jacobian matrix, releases the state of all processes, and MPI ends.

[0045] This invention employs a separate forward and inverse adaptive grid system, ensuring forward modeling accuracy without increasing the parameters of the inverse grid, thus improving both computational efficiency and accuracy. The introduction of an MPI strategy for parallel frequency operations further enhances computational efficiency. Furthermore, the use of unstructured tetrahedral mesh generation improves the applicability of 3D inversion imaging of magnetotelluric data in complex and undulating regions. This invention improves the flexibility and computational efficiency of unstructured mesh generation, as well as its applicability to complex geological structures and undulating terrain, making it of significant practical value in geomagnetic measurements for deep-earth exploration and complex, rugged mineral resource exploration. Attached Figure Description

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

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

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

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

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

[0051] Figure 6 This is a graph showing the root mean square (RMS) fitting and objective function value curves in an embodiment of the present invention. Detailed Implementation

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

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

[0054] S1. The model is divided using an unstructured tetrahedral mesh;

[0055] S2. The forward response of the model is calculated using the unstructured finite element method with Coulomb gauge potential.

[0056] S3. Adaptively refine the forward mesh;

[0057] S4. The inversion is performed using the finite-memory quasi-Newton method (L-BFGS) to obtain the inversion imaging results.

[0058] Example: To test the feasibility of this imaging method, for example... Figure 2The model of multiple anomalies with undulating terrain is shown to be inverted and tested. The model is a sloping terrain with a low-resistivity anomaly and a high-resistivity anomaly. The resistivity of the low-resistivity anomaly is 10Ωm, the resistivity of the high-resistivity anomaly is 3000Ωm, the resistivity of the surrounding rock is 100Ωm, and the geometric dimensions of the two anomalies are 8km×8km×4km. Figure 2 (a) shows the planar location of the anomaly in the model and the distribution of magnetotelluric measurement points. The distance between magnetotelluric measurement points and the distance between lines are both 2km. There are a total of 7 measurement lines, with 13 measurement points on each line. The distribution of measurement points completely covers the range of the two anomalies, for a total of 91 measurement points. Figure 2 (b) in the figure is a vertical depth slice of the model. The slope terrain rises 1 km in elevation from west to east, distributed between -3 km and 3 km in the east-west direction. The top of the low-resistivity anomaly is buried at a depth of 2 km, and the top of the high-resistivity anomaly is buried at a depth of 3 km. The magnetotelluric observation frequencies range from 4 Hz to 0.01 Hz, with a total of 10 frequencies. 3% Gaussian random noise was added to the model response data to simulate the measured data.

[0059] The model is discretized using an unstructured tetrahedral mesh. The model inversion mesh is as follows: Figure 3 As shown, it contains a total of 42,639 nodes and 259,558 elements, with a total of 145,525 inversion parameters. Due to the use of independent forward and inversion mesh systems, the forward mesh is adaptively refined from the inversion mesh. The forward mesh of this embodiment is as follows: Figure 4 As shown, Figure 4 (a) and (b) in the diagram correspond to forward modeling meshes at frequencies of 4Hz and 0.14Hz, respectively. The 4Hz mesh, after refinement, contains 63,066 nodes and 369,043 core elements, while the 0.14Hz mesh, after refinement, contains 45,577 nodes and 272,708 elements. From... Figure 4 As can be seen, the mesh refinement is inconsistent at different frequencies. At lower frequencies, the mesh is sparser, mainly with interpolation refinement near the measurement points, while at higher frequencies, the mesh is denser, which is consistent with the propagation law of electromagnetic waves.

[0060] Figure 5 This is a comparison between the final inversion result of the undulating terrain model and the actual model, where... Figure 5 (a) in the model is the real model. Figure 5 (b) in the figure shows the final inversion result. As can be seen from the figure, the inversion result recovers the spatial location and morphology of the two anomalies under undulating terrain quite well. In particular, the resistivity and morphology of the low-resistivity anomaly are closer to the real model, indicating that the inversion result is accurate and reliable. This embodiment demonstrates that the present invention is accurate and feasible, capable of inverting and imaging target bodies under undulating terrain. Due to the use of unstructured mesh partitioning, it is beneficial for partitioning complex undulating terrain models, and can save a large number of mesh units compared to regular meshes.

[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 decreased from 2.91 to below 0.92, and the objective function value decreased from 6.92 × 10⁴ to 3.09 × 10³. Both curves converged stably, indicating that the inverted imaging results can stably converge to near the true model. Because the observation frequency was 10, 11 threads were used for parallel computation, and the inversion time was 4.1 hours. This embodiment verifies that the present invention has high computational efficiency.

Claims

1. A self-adaptive 3D fast inversion imaging method for magnetotelluric data, characterized in that, It comprises the following steps: S1, using unstructured tetrahedral mesh to subdivide the model; S2, using the unstructured finite element method of Coulomb norm potential to calculate the forward response of the model; S3, self-adaptive encryption of the forward grid; S4, using the limited memory quasi-Newton method L-BFGS inversion to obtain the inversion imaging result; The specific process of S1 is: S1.1, establish the poly file of grid subdivision, including the control point coordinates of the geometric model, line segment, face element, the number of each region and the volume constraint value; S1.2, using unstructured grid subdivision software TetGen, input the poly file, generate node, element, neigh, face four files, wherein the node file is the node list file of unstructured grid, containing all the node numbers and three-dimensional coordinates (x, y, z) of the grid, the element file is the tetrahedron list file, containing the node list and attribute number corresponding to each tetrahedron, the neigh file is the tetrahedron number connected with each four face, and the face file is the list file of triangular faces contained in the grid; S1.3, renumber the generated tetrahedral element attribute, the tetrahedral attribute in the inversion region is numbered in order, that is, each tetrahedron is an inversion parameter, and a new tetrahedral file is generated for subsequent forward and inversion; S1.4, assign the initial resistivity value to each region of the grid, generate resistivity file, wherein the inversion region is marked with a number greater than 0, and the non-inversion region is marked as 0; The specific process of S3 is: S3.1, encrypt the unit near the receiving point, find out the unit where the receiving point is located, judge whether the volume meets the requirement of the set receiving point unit volume, if not, encrypt it in the subsequent grid refinement; Repeat the above steps until the volume of the unit meets the requirement; S3.2, loop each unit of the grid, find the minimum value of the distance from the unit to all receiving points, and calculate the maximum volume that the unit should meet according to the distance and skin depth, the calculation formula is as follows: （8）； wherein, is a skin depth, is a control factor, is a minimum distance from the unit to the receiving point; if the unit volume exceeds this maximum volume constraint, it is refined in a subsequent mesh refinement; the above steps are repeated until the volume of all units meets the requirements. S3.3, since an independent forward and inversion double grid system is adopted, each unit that needs to be inverted in the inversion grid is marked as an independent unit, and after grid refinement, the parameters of each refined unit in the inversion grid are transmitted to each subunit after refinement.

2. The adaptive 3D fast inversion imaging method for MT data according to claim 1, wherein, The specific process of S2 is: S2.1, the differential equations satisfied by electric field E and magnetic field H are as follows: （1）； （2）； wherein is the imaginary unit, is the angular frequency, is the vacuum permeability, is the electrical conductivity; S2.2 According to the definition of Coulomb potential, the electric field E and the magnetic field H can be expressed as a magnetic vector potential A and an electric scalar potential respectively. The formula is as follows: （3）； （4）； S2.3, Using the Coulomb gauge condition: Substituting equations (3) and (4) into equations (1) and (2), we obtain the Helmholtz equation for the potential: Substituting equations (3) and (4) into equations (1) and (2), we obtain the Helmholtz equation for the potential: （5）； （6）； S2.4, the model region is subdivided into tetrahedral elements, interpolation is carried out by using tetrahedral element shape function, integration is carried out for each element, and the total linear equation set is obtained by synthesizing the integration of each element: （7）； where K is the stiffness matrix, u is the unknown A, potential; S2.5, the boundary conditions are put into the linear equation set (7), a linear equation set solver Pardiso is selected to solve, A, The finite element solution of the potential is obtained, and the finite element solution of the electromagnetic field is obtained according to formulas (3) and (4).

3. The adaptive 3D fast inversion imaging method for MT data according to claim 1, wherein, The specific process of S4 is: S4.1, Constructing the inversion objective function by the Tikhonov regularization method, the objective functional is expressed as: （9）； where m is the N-dimensional model parameter vector, generally the log resistivity values; is the prior or initial model parameter vector; R is the roughness operator matrix; μ is the Lagrange multiplier to balance the model constraint term and the data fitting term; W is the data covariance matrix, a diagonal matrix composed of data error terms; d is the observation data vector; F(m) is the forward response corresponding to the model m; S4.2, the gradient of the objective function (9) is obtained: （10）； wherein, J is the Jacobian matrix, which is obtained by the adjoint equation method during the forward modeling; S4.3, the objective function is solved by iteration method, and the model vector is updated constantly to make the objective function reach the minimum, and the model update iteration form is represented as: （11）； wherein, is a search step, controlling the size of the modification of the model, is a search direction, controlling the direction of the modification of the model.

4. The adaptive 3D fast inversion imaging method for MT data according to claim 1, characterized in that, The S2 and S3 further comprise 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 group for one frequency, and numbered in sequence; the main process distributes the calculation tasks to the sub-processes in sequence, and if the number of sub-processes is insufficient, the remaining sub-processes are waited to run and then the distribution is continued, and the process is completed through message transmission between the main process and the sub-processes; after receiving the sub-tasks, each sub-process performs adaptive encryption on the forward grid 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 calculation results of the model response and Jacobian matrix, releases the state of all processes, and the MPI ends.

Images