Spherical domain magnetotelluric non-structural vector finite element inversion method capable of automatically modeling

Abstract

The invention discloses a spherical domain magnetotelluric non-structural vector finite element inversion method capable of automatically modeling, which comprises the following steps of: automatically establishing a spherical domain initial model through global topographic data, and performing spatial discretization by using a non-structural tetrahedral mesh to realize three-dimensional magnetotelluric non-structural finite element solution of two polarization modes in a spherical domain; constructing a spherical domain magnetotelluric three-dimensional regularization inversion objective function according to the actual observation data and the prediction data; adopting adjoint forward modeling to solve a Jacobian matrix and a product of transpose and a vector of the Jacobian matrix, and utilizing a nonlinear conjugate gradient method to realize iterative inversion of spherical domain magnetotelluric, so as to obtain a final inversion model and complete spherical domain magnetotelluric non-structural vector finite element inversion. According to the method, spherical domain magnetotelluric non-structural vector finite element inversion is realized, errors caused by coordinate projection can be avoided, and a large-scale inversion result is more accurate and more conforms to a real earth model.

Description

Technical Field

[0001] This invention relates to the field of geophysical technology, and more specifically to an automatic modeling method for unstructured vector finite element inversion of spherical magnetotelluric data. Background Technology

[0002] Magnetotelluric (MT) is a geophysical electromagnetic detection technique that uses natural alternating electromagnetic fields as a source to observe the electric and magnetic responses at the Earth's surface, thereby obtaining the distribution of electrical charge distribution near the surface and even the transition zone between the crust and mantle. It has advantages such as large detection depth, no shielding by high-resistivity layers, strong resolution of low-resistivity layers, and strong lateral resolution. Therefore, it is widely used in fields such as the study of deep Earth structure and dynamics, exploration of natural resources such as minerals and geothermal energy, and environmental and geological disaster monitoring.

[0003] With the support of deep structure exploration projects implemented by various countries, a large number of three-dimensional magnetotelluric (MT) field datasets urgently need to be analyzed and interpreted. However, most magnetotelluric forward and inverse studies to date are still based on the assumption of plane waves in the natural incident field and the assumption of a flat surface of the Earth's medium. However, the actual landmass structure is an undulating, spherical shape. Existing studies have shown that for large-scale, long-period magnetotelluric data, the errors caused by coordinate projection and the curvature of the Earth cannot be ignored. Therefore, the traditional Cartesian plane rectangular coordinate system and the three-dimensional magnetotelluric inversion system without topography can no longer meet the requirements for the inversion of three-dimensional magnetotelluric data at the land scale.

[0004] The finite difference method is easy to implement on computers and highly practical, enabling rapid computational simulations. However, it typically divides the computational domain into regular hexahedral meshes, making it poorly adaptable to complex geological environments, arbitrarily undulating terrain, and the curvature of the Earth, thus hindering accurate simulation. In contrast, the unstructured vector finite element method can flexibly handle the modeling of complex structures, providing more accurate simulations of undulating terrain and realistic Earth models. Furthermore, it automatically satisfies the dispersion condition of the electromagnetic field, effectively avoiding the spurious solution problem of the nodal finite element method. However, both methods are currently rarely used for 3D magnetotelluric inversion in the spherical domain. Summary of the Invention

[0005] To address the above technical problems, this invention provides an automatic modeling method for unstructured vector finite element inversion of spherical magnetotelluric data. The specific technical solution is as follows:

[0006] An automated modeling method for unstructured vector finite element inversion of spherical magnetotelluric data includes the following steps:

[0007] Step 1: Automated modeling to obtain the inversion model and initialize the inversion parameters;

[0008] Step 2: Obtain the observation dataset for inversion;

[0009] Step 3, the initial model in step 1 is subjected to spherical domain non-structural vector finite element magnetotelluric forward to obtain predicted data and combined with the observation data in step 2 to build the objective function of the spherical domain magnetotelluric three-dimensional regularization inversion;

[0010] Step 4, the derivation of both sides of the forward equation, wherein the partial derivative of the predicted data to the model parameter is the Jacobian matrix; the product of the Jacobian matrix and its transpose and the vector is obtained by solving the equation using the adjoint forward;

[0011] Step 5, the model update amount is calculated by using the nonlinear conjugate gradient method, and a new inversion model is obtained;

[0012] Step 6, the model forward response and the root mean square error are calculated .

[0013] Step 7, judge whether the inversion reaches the termination condition.

[0014] The present application has the following beneficial effects:

[0015] In consideration of the curvature of the earth, the spherical domain magnetotelluric non-structural vector finite element inversion is realized, which can avoid the error caused by coordinate projection in the traditional magnetotelluric inversion, so that the large-scale inversion result is more accurate and more consistent with the real earth model. BRIEF DESCRIPTION OF DRAWINGS

[0016] Figure 1 The method flowchart of the present application. DETAILED DESCRIPTION

[0017] In order to make the purpose, technical scheme and advantages of the present application clearer and more apparent, the present application will be further described in detail below in combination with the drawings and examples. It should be understood that the specific examples described herein are only used to explain the present application and do not limit the protection scope of the present application.

[0018] In one embodiment, as shown in Figure 1 , an automatic modeling spherical domain magnetotelluric non-structural vector finite element inversion method is provided, comprising the following steps:

[0019] Step 1, automatic modeling to obtain an inversion model (or iterative model) and initialize the inversion parameters; specifically including the following steps:

[0020] Firstly, the terrain corresponding to the region is discretized by triangle to obtain the triangular grid nodes, and then the elevation information of the triangular grid nodes is obtained by linear interpolation using DEM elevation data, and then the 3D spherical domain terrain interface is obtained. Finally, based on the obtained 3D spherical domain terrain interface grid information, the inversion region information and the measuring point information, Gmsh is called again to model the inversion model. At the same time, the conductivity information of the inversion model, the upper and lower limits of the conductivity of the inversion model and the maximum number of iterations are initialized.

[0021] Step 2, obtaining an observation data set for inversion; specifically including the following steps:

[0022] In the embodiment, the observation data obtained after corresponding processing is inverted by using the real part and the imaginary part of the impedance tensor.

[0023] Step 3, performing spherical non-structural vector finite element magnetotelluric forward for the inversion model (or iterative model) in step 1 to obtain predicted data, and combining the observation data in step 2 to construct a spherical magnetotelluric three-dimensional regularization inversion objective function; specifically including the following steps:

[0024] In the embodiment, a regularization inversion objective function based on L2 norm is constructed, which is used for optimization of the inversion model in the regularization inversion process. According to the observation data, the predicted data obtained by forward, and the inversion model, the regularization inversion objective function constructed is:

[0025] ;

[0026] Wherein, represents the discrete model parameter vector of the spherical calculation domain to be solved, that is, the conductivity vector, is the parameter vector (conductivity vector) of the inversion reference model, is the observation data vector (such as impedance tensor component , etc.), is the predicted data vector ( , is the forward operator, is the model weight matrix, is the data weight matrix, is the regularization factor.

[0027] The initial model (or iterative model) in step 1 is subjected to spherical non-structural vector finite element magnetotelluric forward to obtain predicted data, which includes: when the inversion model is subjected to non-structural finite element forward calculation, the electromagnetic field satisfies the following double spin equation:

[0028] ;

[0029] Wherein, respectively refer to the imaginary unit, the angular frequency and the magnetic permeability in vacuum is the electric field intensity with unit V / m, is the magnetic field intensity with unit A / m, is the dielectric conductivity, is the curl operator. According to the Galerkin method, the vector interpolation basis function is applied to the electric field double curl equation in the element space as a weighting function is any differentiable vector, the finite element control equation can be obtained:

[0030]

[0031] where, is the mass matrix of element e, is the stiffness matrix of element e, . Here the superscript indicates that the current physical quantity is the physical quantity inside any one element, the subscript indicates the local number of the first edge, the subscript indicates the local number of the first edge, indicates the local number of the fourth edge, indicates the volume of the tetrahedral element.

[0032] To solve the finite element control equation, the tetrahedral elements discretized by Gmsh need to be assembled into the total stiffness matrix and the boundary conditions need to be applied. At the same time, the full impedance tensor response of the spherical domain magnetotelluric field needs to be solved by solving two types of linear equations under orthogonal polarization field sources respectively: . Here is the total system matrix, is the total mass matrix, is the total stiffness matrix, is the electric field intensity to be solved, is the right end item of the equation after applying the boundary conditions. The two types of orthogonal polarization field sources selected in the spherical domain magnetotelluric field are , where are the latitude and longitude respectively. By solving the equation set by the Pardiso direct solver, the electric field intensity on all edges under the two types of polarization sources is obtained, and the edge magnetic field intensity in all tetrahedral elements is obtained by . Finally, the full impedance tensor response of the magnetotelluric field at the observation point is obtained by interpolation.

[0033] ​​​ ;

[0034] where superscripts 1 and 2 represent two different spherical orthogonal polarization sources, and superscript -1 represents matrix inversion.

[0035] Step 4: Derivation is performed on both sides of the forward equation, where the partial derivative of the predicted data with respect to the model parameters is the Jacobian matrix; the product of the Jacobian matrix and its transpose and a vector is obtained by solving the equation using the adjoint forward; specifically including the following steps:

[0036] In the embodiment, when the Jacobian matrix of the predicted data is calculated, the forward equation is simultaneously derived with respect to the model parameters (i.e., the electrical conductivity ), and the following is obtained:

[0037] ;

[0038] where is the total system matrix, and superscript -1 represents inversion operation. In spherical magnetotelluric forward calculation, the predicted data is obtained by interpolation of the electric field intensity , that is, , where is an interpolation function for interpolating the electric field to the corresponding receiving point position, and thus the Jacobian matrix of the predicted data can be written as:

[0039] ;

[0040] where represents the kth column of the Jacobian matrix, represents the number of discrete tetrahedrons in the inversion model, and since total field-based forward calculation is adopted, there is , and the matrix is defined as:

[0041] ;

[0042] Therefore, the Jacobian matrix and its transpose can be expressed as

[0043] ;

[0044] The product of the Jacobian matrix and its transpose and an arbitrary vector is converted into forward calculation, avoiding the display storage of the Jacobian matrix, which can effectively save the memory. An arbitrary vector is defined as , and two adjoint source vectors

[0045] ;

[0046] Then we get two adjoint equations as follows:

[0047] ;

[0048] where, are the solution vectors of the two adjoint equations respectively.

[0049] Solving the above adjoint equations, we can get the product of the Jacobian matrix and its transpose and any vector :

[0050] ;

[0051] The product of the Jacobian matrix and its transpose and any vector can be obtained by calculating the solutions of the adjoint equations for two polarizations simultaneously.

[0052] ;

[0053] ;

[0054] where, is the interpolation function of the electric field solution of two different polarizations interpolated to the corresponding receiving point position, is the electric field intensity obtained by the spherical MT 3D forward for two different polarizations, and T represents the matrix transpose. This method makes the calculation time of the gradient vector close to single forward, which can greatly improve the efficiency of 3D MT inversion.

[0055] Step 5, using the nonlinear conjugate gradient method to calculate the model update, obtaining a new iteration model; specifically including the following steps:

[0056] In the embodiment, according to the objective function constructed in step 3, the nonlinear conjugate gradient iteration algorithm is constructed, first, define as the data fitting difference of the nth iteration step:

[0057] wherein n is the current iteration step number, and then according to the inversion objective function, the steepest descent direction is obtained,

[0058] ;

[0059] wherein, is the regularization factor of the nth iteration step; denotes the Jacobian matrix of the nth iteration step

[0060] Using the conjugate gradient method, the conjugate gradient descent direction is obtained,

[0061] ;

[0062] The model descent distance at the nth iteration step is thus constructed ,

[0063] ;

[0064] The inversion model is further updated,

[0065] ;

[0066] wherein, is the parameter vector of the inversion model at the nth iteration step, is the parameter vector of the inversion model at the (n+1)th iteration step;

[0067] In the nonlinear conjugate gradient iteration, the Jacobian matrix of the predicted data and its transpose, the model roughness and the regularization factor need to be calculated. The product of the Jacobian matrix and its transpose with an arbitrary vector has been calculated in step 4. In the present embodiment, the model roughness is defined as:

[0068] ;

[0069] wherein,

[0070] ;

[0071] ;

[0072] ;

[0073] represents the number of discrete tetrahedrons in the inversion model, represents the volume of the discrete tetrahedron (denoted as inversion unit ) in the inversion model, represents the number of all inversion units sharing a common vertex with the inversion unit , represents the number of inversion units sharing a common vertex with the inversion unit , represents the difference between the model parameter values of the inversion unit and the inversion unit sharing a common vertex with it (denoted as inversion unit ), represents the distance between the centers of the inversion unit and the inversion unit , represents the distance between the centers of the inversion unit and the inversion unit The ratio of the volume of the body to the volume of the inversion unit body of its common vertex. Wherein the model weight matrix Only need to calculate once after the grid partition is determined, avoiding repeated construction in the inversion iteration process, high computational efficiency.

[0074] In the nonlinear conjugate gradient iteration, the regularization factor plays a role in preventing model overfitting, in this embodiment, in the 0th step iteration, the regularization factor is set to ; and the regularization factor in the 1st step iteration is set to Its expression is:

[0075] ;

[0076] In the formula, is the parameter vector of the inversion model of the 1st step iteration, is the data weight matrix;

[0077] ;

[0078] Wherein, is the ith observation data, is the relative error of the ith observation data, is a small amount to prevent overfitting of relatively small observation data.

[0079] In this embodiment, the cooling method is used to update the regularization factor iteratively:

[0080] ;

[0081] In this way, after giving the product of the Jacobian matrix and its transpose and the vector, the model roughness and the regularization factor, nonlinear conjugate gradient iteration can be performed to update the inversion model.

[0082] Step 6, calculate the model forward response and the root mean square error ; Specifically, the following steps are included:

[0083] In this embodiment, the updated model is reprocessed by the spherical domain magnetotelluric three-dimensional non-structural vector finite element forward to obtain the predicted data at the observation point, that is, the impedance tensor. In order to describe the fitting degree between the predicted data and the observation data, define , that is, the root mean square error of the data, its expression is:

[0084] ;

[0085] Wherein, is the data fitting error, is the data amount of the observation data, is the ith predicted data in the predicted data, is the i-th observation data in the observation data, is the error value corresponding to the i-th observation data in the observation data.

[0086] Step 7, judging whether the inversion reaches a termination condition; specifically including the following steps:

[0087] When the termination condition is reached, it is considered that the inversion has obtained a solution conforming to the true electrical structure of the underground, that is, the current inversion model and response can be output, and the inversion is terminated; in the embodiment, the inversion termination condition is set as:

[0088] the inversion iteration number reaches a preset maximum iteration number;

[0089] the data root mean square error, i.e. is less than a threshold value between two consecutive iterations, expressed as:

[0090] ;

[0091] wherein, subscript k represents the iteration number, represents the threshold value;

[0092] If the inversion termination condition is not reached, steps 3-7 are continued to be iterated until the inversion termination condition is reached.

Claims

1. A self-modeling spherical domain magnetotelluric non-structural vector finite element inversion method, characterized in that: The method comprises the following steps: Step 1, automatically modeling to obtain an inversion model and initialize inversion parameters; Step 2, obtaining an observation data set for inversion; Step 3, performing spherical non-structured vector finite element magnetotelluric forward calculation on the inversion model in step 1 to obtain predicted data, and combining the observation data in step 2 to construct a spherical magnetotelluric three-dimensional regularization inversion objective function; Step 4, deriving both sides of the forward equation, wherein the partial derivative of the predicted data with respect to the model parameters is a Jacobian matrix; the product of the Jacobian matrix and the transpose of the vector is obtained by solving the adjoint forward equation; Step 5, calculating a model update amount by using a nonlinear conjugate gradient method to obtain a new inversion model; Step 6, Calculate model forward response and root mean square error ; Step 7, judging whether the inversion reaches a termination condition.

2. The automatically modelable spherical domain magnetotelluric non-structural vector finite element inversion method according to claim 1, characterized in that: Step 1 specifically comprises the following steps: First, the region corresponding to the terrain is discretized into triangular meshes to obtain triangular mesh nodes, then the elevation data is used to obtain the elevation information of the triangular mesh nodes by linear interpolation, and then the 3D spherical terrain interface is obtained; finally, based on the obtained 3D spherical terrain interface mesh information, inversion region information and measurement point information, Gmsh is called again to model to obtain an inversion model, and the conductivity information of the inversion model, the upper and lower limits of the inversion model conductivity and the maximum number of iterations are initialized.

3. The automatically modelable spherical domain magnetotelluric non-structural vector finite element inversion method according to claim 1, characterized in that: Step 2 specifically comprises the following steps: The processed observation data is obtained, and the real part and the imaginary part of the impedance tensor are used for inversion.

4. The automatically modelable spherical domain magnetotelluric non-structural vector finite element inversion method according to claim 1, characterized in that: Step 3 specifically comprises the following steps: A regularization inversion objective function based on L2 norm is constructed, and the regularization inversion objective function is: ； wherein denotes the discrete model parameter vector of the sphere computational domain to be solved, i.e. the conductivity vector, is the parameter vector of the inversion reference model, is the observation data vector, is the prediction data vector, , is the forward operator, is the model weight matrix, is the data weight matrix, is the regularization factor.

5. The automatically modelable spherical domain magnetotelluric non-structural vector finite element inversion method according to claim 4, characterized in that: In step 3, the initial model in step 1 is subjected to spherical non-structured vector finite element magnetotelluric forward calculation to obtain predicted data, which comprises: When the inversion model is subjected to non-structured finite element forward calculation, the electromagnetic field satisfies the following double curl equation: ； where respectively refer to the imaginary unit, the angular frequency and the magnetic permeability in vacuum; is the electric field intensity, is the magnetic field intensity, is the dielectric permittivity, is the curl operator, the vector interpolation basis functions as a weighting function, the vector identity is applied to the electric field biharmonic equation within the element space is an arbitrary differentiable vector, the finite element governing equation can be obtained: ； wherein, is the mass matrix of element e, is the stiffness matrix of element e, , denotes the physical quantity in the interior of an arbitrary element, the index denotes the local number of the th edge, the index denotes the local number of the th edge, denotes the volume of the tetrahedral element.

6. The automatically modelable spherical domain magnetotelluric non-structural vector finite element inversion method according to claim 5, characterized in that: In step 3, the tetrahedral elements obtained by Gmsh are assembled into the total system matrix and boundary conditions are applied to solve the finite element control equation, and the full impedance tensor response of the spherical domain magnetotelluric field is obtained , respectively, the linear equations under two types of orthogonal polarization field sources are solved: , is the total system matrix, is the total mass matrix, is the total stiffness matrix, is the electric field intensity to be solved, is the right end item of the equation after applying the boundary conditions, and the two types of orthogonal polarization field sources selected in the spherical domain magnetotelluric field are , where are the latitude and longitude, respectively, and the electric field intensity on all edges under two types of polarization sources is obtained by solving the equation set , and the edge magnetic field intensity in all tetrahedral elements is obtained by , and finally the magnetotelluric full impedance tensor response at the observation point is obtained by interpolation ; ； Wherein, the superscripts 1 and 2 represent two different spherical orthogonal polarization sources, and the superscript -1 represents matrix inversion.

7. The automatically modelable spherical domain magnetotelluric non-structural vector finite element inversion method according to claim 6, characterized in that: Step 4 specifically comprises the following steps: When computing the Jacobian matrix of the predicted data, the forward equation simultaneously from both ends derivation, we have ； where, is the total system matrix, in the spherical domain MT forward problem, the predicted data is obtained by interpolating the electric field intensity , i.e. , where is the interpolation function that interpolates the electric field solution to the corresponding receiver location, then the Jacobian matrix of the predicted data is written as: ； wherein, denotes the k-th column of the Jacobian matrix, denotes the number of discrete tetrahedrons in the inverse model, since forward calculations are performed based on the total field, thus , the matrix is defined as: ； Therefore, the Jacobian matrix and its transpose are represented as: ； The product of the Jacobian matrix and its transpose with an arbitrary vector is transformed into a forward calculation, defining an arbitrary vector as Two adjoint source vectors are redefined as : ； Then the following two adjoint forward equations are obtained: ； wherein, respectively denote the solution vectors of the corresponding two adjoint forward equations; Solving the above adjoint forward equation gives the product of the Jacobian matrix and its transpose with any vector : ； The Jacobian matrix and its transpose of the spherical domain MT under the MT are obtained by calculating the product of the results of the two polarization modes and any vector simultaneously with the forward solution. ； ； wherein, is an interpolation function that interpolates the electric field at two different polarizations to the corresponding receiver location, are the electric field strengths obtained from the 3D forward modeling of the spherical domain MT for two different polarizations, and T denotes the matrix transpose.

8. The automatically modelable spherical domain magnetotelluric non-structural vector finite element inversion method according to claim 7, characterized in that: Step 5 specifically comprises the following steps: Based on the objective function constructed in Step 3, a nonlinear conjugate gradient iterative algorithm is constructed. First, define Data fitting difference for the nth iteration step: ； where n is the current iteration step number, and the steepest descent direction is obtained according to the inversion objective function , ； wherein, regularization factor for the n-th iteration; Jn represents the Jacobian matrix for the n-th iteration Using the conjugate gradient method, the conjugate gradient descent direction is obtained : ； Thus the model drop distance at the nth iteration step is constructed : ； Then the inversion model is updated: ； wherein, is the parameter vector of the inversion model for the n-th iteration, is the parameter vector of the inversion model for the n+1-th iteration; When performing nonlinear conjugate gradient iteration, the Jacobian matrix and its transpose of the predicted data, the model roughness and the regularization factor are used, and the model roughness is defined as: ； In the formula, ； ； ； This indicates the number of discrete tetrahedra in the inversion model. In the inversion model, the first The volume of a discrete tetrahedron Representation and Inversion Units The number of all inversion units sharing a vertex. Representation and Inversion Units The inversion unit numbering of the shared vertex Indicates the inversion unit Inversion unit sharing a vertex with it (referred to as inversion unit) The difference in model parameter values ​​between ) Indicates the inversion unit and inversion unit The distance between centers Indicates the inversion unit The ratio of the volume of a given element to the sum of the volumes of its shared inversion units, where the model weight matrix is... Calculate once after the mesh is determined; In the 0th iteration, the regularization factor is set to ; while in the 1st iteration, the regularization factor is set to , which is expressed as: ； wherein is the parameter vector of the inverse model of the 1st iteration, is the data weight matrix; ； wherein, is the ith observation, is the relative error of the ith observation, is a small quantity to prevent overfitting of relatively small observations; The cooling method is used to update the regularization factor iteratively: ； After the product of the Jacobian matrix and its transpose and the vector, the model roughness and the regularization factor are given, the nonlinear conjugate gradient iteration is performed to update the inversion model.

9. The automatically modelable spherical domain magnetotelluric non-structural vector finite element inversion method according to claim 1, characterized in that: Step 6 specifically comprises the following steps: The updated model is used to perform the spherical domain magnetotelluric 3D non-structural vector finite element forward, to obtain the predicted data at the observation point, i.e. impedance tensor, defined as The root mean square error of data, which describes the fitting degree between the predicted data and the observed data, is defined as: ； wherein, is the data fitting difference, is the data quantity of the observation data, is the i-th prediction data in the prediction data, is the i-th observation data in the observation data, is the error value corresponding to the i-th observation data in the observation data.

10. The automatically modelable spherical domain magnetotelluric non-structural vector finite element inversion method according to claim 1, characterized in that: Step 7 specifically comprises the following steps: The inversion termination condition is set as: the number of inversion iterations reaches a preset maximum number of iterations; The data root mean square error is given by The expression is less than a threshold between two consecutive iterations. ； wherein the subscript k denotes the iteration number, is expressed as a threshold value.

Images