A three-dimensional magnetotelluric inversion method based on continuous conductivity block variation

Abstract

The invention relates to the field of geophysics and discloses a three-dimensional magnetotelluric inversion method for block-by-block continuous variation of conductivity. The method includes adding a model constraint term based on the conductivity of an unstructured tetrahedral grid node to an objective function to obtain an inversion objective function expression. The model constraint term is solved by minimizing the norm of the conductivity variation between a unit and surrounding units to obtain a model constraint term based on the conductivity of the tetrahedral grid node. The inversion objective function expression is gradient-solved. Specific expressions of elements in the stiffness matrix gradient under four different conditions are solved and the gradient components of the stiffness matrix corresponding to all nodes are solved, thereby achieving the solution of the objective function gradient in magnetotelluric inversion. The method realizes the updating of the node conductivity model during three-dimensional underground medium inversion, greatly reduces the unknowns generated in the inversion process, improves the uniqueness of the solution, and greatly shortens the time required for the inversion process to obtain a better inversion effect.

Description

Technical Field

[0001] The present invention relates to the field of geophysics, and in particular to a three-dimensional magnetotelluric inversion method for block-by-block continuous variation of electrical conductivity. Background Art

[0002] Magnetotelluric inversion is a quantitative interpretation method that infers the distribution of underground conductivity based on electromagnetic field data observed on the surface. In three-dimensional magnetotelluric inversion, conductivity parameters are usually established on tetrahedral units, that is, a fixed conductivity value is assigned to each grid unit. Since the inversion unknowns are the number of model units, and the corresponding number of grid units in three-dimensional magnetotelluric inversion is huge, there are many inversion unknowns and the solution is highly multi-solutionable.

[0003] The present invention adopts a parameterized assignment method that establishes conductivity on tetrahedral unit nodes in magnetotelluric inversion. Because the number of tetrahedral unit nodes in unstructured grids is much smaller than the number of tetrahedral units, this conductivity assignment method can greatly reduce the inversion unknowns and improve the uniqueness of the solution. At the same time, it adapts to the situation of continuous change of underground medium conductivity and obtains better inversion effect. Therefore, a three-dimensional magnetotelluric inversion algorithm based on the continuous change of conductivity block is proposed. Summary of the Invention

[0004] The present invention aims to provide a three-dimensional magnetotelluric inversion method with block-by-block continuous variation of conductivity. By constructing model terms based on unit nodes and solving a new stiffness matrix gradient, the node conductivity model is updated during the three-dimensional underground medium inversion, which greatly reduces the unknowns generated in the inversion process, improves the uniqueness of the solution, and greatly shortens the time required for the inversion process, thereby obtaining better inversion results.

[0005] In order to achieve the above object, the present invention provides the following method:

[0006] The present invention provides a three-dimensional magnetotelluric inversion method for continuous block variation of conductivity:

[0007] The method includes solving model constraints based on node conductivity parameters and solving a new stiffness matrix gradient, and is characterized in that the method includes:

[0008] S1: The idea of ​​Tikhonov regularized inversion is introduced, and a model constraint term based on the conductivity of the unstructured tetrahedral grid nodes is added to the objective function. The inversion objective function expression based on the node model parameters is obtained, including the data fitting error term and the three-dimensional model constraint term;

[0009] S2: Minimize the norm of the conductivity change between the unit and the surrounding units, constrain the smoothness of the model, solve the model constraint terms for the node resistivity parameters, and obtain the model constraint terms based on the conductivity of the tetrahedral mesh nodes;

[0010] S3: directly constructing a three-dimensional network model through the model constraint items based on the conductivity of the tetrahedral mesh nodes, and calculating the three-dimensional model constraint items at the nodes of the three-dimensional network model;

[0011] S4: performing a gradient solution on the inversion objective function expression according to the data fitting difference term and the three-dimensional model constraint term;

[0012] S5: Solve the specific expression of each element in the stiffness matrix gradient under four different conditions, and solve the gradient components of the stiffness matrix corresponding to all nodes, thereby realizing the solution of the objective function gradient in magnetotelluric inversion.

[0013] Preferably, the idea of ​​Tikhonov regularized inversion is introduced, and a model constraint term based on the conductivity of the unstructured tetrahedral grid nodes is added to the objective function, and the inversion objective function expression based on the node model parameters is obtained as follows:

[0014]

[0015] Among them, in the above formula, The data fitting difference term is used to fit the observed data; It is a model constraint item, providing the model’s prior information so that the model maintains a certain degree of smoothness; m=[m1,m2,…,m M ] T is the inversion model parameter vector, that is, the conductivity model parameter vector established on the tetrahedral grid nodes; M is the number of model parameters, and λ is the regularization factor used to balance the data fitting term and the model constraint term.

[0016] Preferably, the expression of the data fitting difference term is as follows:

[0017]

[0018] where d=[d1,d2,…,d N ] T is the observation data vector, N is the size of the observation data, d is the full tensor impedance data or full tensor impedance data and dipole vector impedance, F(m) is the forward response, W d is an N×N order data weighting matrix containing data errors, which represents the weight of each data in the inversion process. The added noise is expected to be 0 and the standard deviation is ε i The random number is Gaussian noise; the data weight matrix is ​​written as W d =diag{1 / ε1,1 / ε2,…,1 / ε N}, when the standard deviation ε iThe larger it is, the smaller its weight in the data weighting matrix is. i The smaller it is, the greater its weight in the data weighting matrix is. The full tensor impedance response is selected as the observation data, where the impedance tensor response expression is:

[0019]

[0020] In the formula, superscripts 1 and 2 represent two polarization directions.

[0021] Preferably, the step of minimizing the norm of the conductivity variation between the unit and the surrounding units, constraining the smoothness of the model, solving the model constraint terms for the node resistivity parameters, and obtaining the model constraint terms based on the tetrahedral mesh node conductivity includes:

[0022] In continuous space, the model constraint based on the L2 norm is written as:

[0023]

[0024] The solution to the discretized model constraints is:

[0025]

[0026] Among them, C m is the model covariance matrix, constraining the model smoothness, W m is the model weight matrix.

[0027] Preferably, the step of directly constructing a three-dimensional network model through the model constraint items based on the conductivity of the tetrahedral grid nodes and calculating the three-dimensional model constraint items at the nodes of the three-dimensional network model includes: obtaining the model parameters at a node and finding all adjacent tetrahedral units containing this node; summing the product of the model parameters of all tetrahedral units containing the node and the gradient of the linear interpolation basis function to obtain the model constraint items of each tetrahedral unit; summing the model constraint items of all adjacent inversion units at the node and taking the average value to obtain the model constraint item of the node. The expression of the model constraint term based on the unit node conductivity parameter is:

[0028]

[0029] Where M is the number of tetrahedral unit nodes in the underground area, V i is the volume of the i-th bottom tetrahedron, and are the conductivity parameter of the j-th node of unit i and the gradient of the linear interpolation basis function, respectively, and L2 represents the second norm of the orientation quantity.

[0030] Preferably, the step of performing a gradient solution on the inversion objective function expression according to the data fitting difference term and the three-dimensional model constraint term comprises:

[0031] Provide a priori model and add it to the constraints of the three-dimensional model:

[0032]

[0033] Then the objective function of magnetotelluric inversion is expressed as:

[0034]

[0035] where m r is the prior model, is the data covariance matrix;

[0036] The objective function is partial derivative with respect to m, and Further solve and convert into matrix form to obtain the forward linear equations;

[0037] The partial derivatives of m on both sides of the forward linear equations are obtained, and combined with the expression of the sparse stiffness matrix in, we can further deduce The specific calculation formula of the component is:

[0038]

[0039] Among them, p and q represent the row index and column index of the element respectively, j is a node in the inversion area, σ j represents the conductivity of the inversion unit node j, Ω j is the area covered by all inversion units that share the same vertex with inversion unit node j, k represents the serial number of the inversion unit that shares the same vertex with inversion node j, Ω k Indicates the area where the kth inversion unit is located, L j It is the interpolation basis function of the nodes at the common vertices of the inversion unit.

[0040] Preferably, according to the upper and lower limits of conductivity Where a and b are constants, simplifying them yields:

[0041]

[0042] so It can be further deduced as:

[0043]

[0044] Among them, p and q represent the row index and column index of the element respectively, j is a node in the inversion area, σj represents the conductivity of the inversion unit node j, Ω j is the area covered by all inversion units that share the same vertex with the inversion unit node j, N(j) represents the number of inversion units that share the same vertex with the inversion unit node j, k represents the sequence number of the inversion unit that shares the same vertex with the inversion node j, Ω k Indicates the area where the kth inversion unit is located, L j It is the interpolation basis function of the nodes at the common vertices of the inversion unit.

[0045] Preferably, after performing gradient solving for the inversion objective function expression, the method further comprises: solving the gradient of the stiffness matrix of the inversion unit having the same vertex as the inversion unit node according to the unit integral formula of the volume coordinate; solving the stiffness matrix gradient of the kth inversion unit Ω having the same vertex as the inversion unit node j; k In the equation, the gradient of the stiffness matrix can be expressed as:

[0046]

[0047] In the formula, f ij =a i a j +b i b j +c i c j ;

[0048] According to the unit integration formula of volume coordinates:

[0049]

[0050] The kth inversion unit contains four nodes, and the expression of the stiffness matrix gradient will be different after the inversion unit node j falls on any node; when the node numbers of the inversion unit node j in the tetrahedron unit are 1, 2, 3, and 4 respectively, the L in the gradient integral term of the stiffness matrix is j The expressions are L1, L2, L3, and L4 respectively.

[0051] Preferably, after solving the stiffness matrix gradient of the inversion unit of the common vertex of the inversion unit nodes according to the unit integration formula of the volume coordinate, the method further includes:

[0052] The complex sensitivity matrix is ​​converted into the form of a real sensitivity matrix and then the formula of the complex sensitivity matrix is ​​substituted into the calculation formula of the objective function fitting difference gradient to solve the objective function gradient.

[0053] Preferably, the step of solving the specific expression of each element in the stiffness matrix gradient under four different conditions and solving the gradient components of the stiffness matrix corresponding to all nodes, thereby realizing the step of solving the gradient of the objective function in the magnetotelluric inversion, includes: solving the stiffness matrix gradient D under four different conditions ij The specific expression of each element in is obtained; the gradients of the stiffness matrices of the N(j) inversion units that share a vertex with the inversion unit node j are summed to obtain the gradient of the stiffness matrix at the inversion unit node j; the gradients of the stiffness matrices at all inversion unit nodes are solved to realize the solution of the gradient of the objective function in magnetotelluric inversion.

[0054] The beneficial effects of the present invention are embodied in the following aspects: the specific content of the present invention is to address the problem that in traditional magnetotelluric three-dimensional inversion, the conductivity parameters are established on tetrahedral units. Since the inversion unknowns are the number of model units, and the corresponding number of grid units in the magnetotelluric three-dimensional inversion is huge, which leads to a large number of inversion unknowns and a strong multi-solution problem. Therefore, it is proposed to adopt a parameterized assignment method in magnetotelluric inversion that establishes the conductivity on the tetrahedral unit nodes. Because the number of tetrahedral unit nodes in an unstructured grid is much smaller than the number of tetrahedral units, this conductivity assignment method can greatly reduce the inversion unknowns, improve the uniqueness of the solution, and adapt to the situation where the conductivity of the underground medium changes continuously, thereby obtaining a better inversion effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the following briefly describes the drawings required for the specific embodiments or the description of the prior art. Similar elements or parts are generally identified by similar reference numerals throughout the drawings. Elements or parts in the drawings are not necessarily drawn to scale.

[0056] Figure 1 The present invention provides a flow chart of a three-dimensional magnetotelluric inversion method for block-by-block continuous variation of conductivity.

[0057] Figure 2 This is a schematic diagram of an unstructured triangular network model provided by an embodiment of the present invention.

[0058] Figure 3 It is a schematic diagram of the inversion process provided by an embodiment of the present invention.

[0059] Figure 4 This is a schematic diagram of a high and low resistance anomaly model under flat terrain provided by Example 1 of the present invention.

[0060] Figure 5 This is a diagram of the local densification of the inversion grid near the measuring point provided in Example 1 of the present invention.

[0061] Figure 6 This is a comparison diagram of the cooling curves of the regularization factor during the inversion iteration process provided in Example 1 of the present invention.

[0062] Figure 7 This is a comparison diagram of the RMS convergence curve during the inversion iteration process provided by Example 1 of the present invention.

[0063] Figure 8 This is an inversion result diagram of the high and low resistivity anomaly model under flat terrain provided by Example 1 of the present invention.

[0064] Figure 9 This is a schematic diagram of a high and low resistivity anomaly model in a valley terrain provided by Example 2 of the present invention.

[0065] Figure 10 This is the local densification of the inversion grid near the measuring point provided by Example 2 of the present invention.

[0066] Figure 11 This is a comparison diagram of the cooling curves of the regularization factor during the inversion iteration process provided in Example 2 of the present invention.

[0067] Figure 12 This is a comparison diagram of the RMS convergence curve during the inversion iteration process provided by Example 2 of the present invention.

[0068] Figure 13 This is an inversion result diagram of the high and low resistivity anomaly model under valley terrain provided by Example 2 of the present invention. DETAILED DESCRIPTION

[0069] In order to help those skilled in the art better understand the present invention, the following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0070] The terms "first," "second," and so on, in the description and claims of the present invention and the accompanying drawings are used to distinguish between different items, not to describe a specific order. Furthermore, the terms "including," "having," and any variations thereof, are intended to cover non-exclusive inclusions. For example, a process, method, apparatus, product, or end comprising a series of steps or elements is not limited to the listed steps or elements but may optionally include steps or elements not listed therein, or may optionally include other steps or elements inherent to such process, method, product, or end.

[0071] References herein to "embodiments" mean that a particular feature, structure, or characteristic described in connection with the embodiments may be included in at least one embodiment of the present invention. The appearance of this phrase in various places in the specification does not necessarily refer to the same embodiment, nor does it constitute a separate or alternative embodiment that is mutually exclusive of other embodiments. It is understood, both explicitly and implicitly, by those skilled in the art that the embodiments described herein may be combined with other embodiments.

[0072] The present invention adopts a parameterized assignment method that establishes conductivity on tetrahedral unit nodes in magnetotelluric inversion. Because the number of tetrahedral unit nodes in unstructured grids is much smaller than the number of tetrahedral units, this conductivity assignment method can greatly reduce the inversion unknowns and improve the uniqueness of the solution. At the same time, it adapts to the situation of continuous change of underground medium conductivity and obtains better inversion effect. Therefore, a three-dimensional magnetotelluric inversion algorithm based on the continuous change of conductivity block is proposed.

[0073] The present invention aims to provide a three-dimensional magnetotelluric inversion method with block-by-block continuous variation of conductivity. By constructing model terms based on unit nodes and solving a new stiffness matrix gradient, the node conductivity model is updated during the three-dimensional underground medium inversion, which greatly reduces the unknowns generated in the inversion process, improves the uniqueness of the solution, and greatly shortens the time required for the inversion process, thereby obtaining better inversion results.

[0074] The specific embodiment of the present invention provides a three-dimensional magnetotelluric inversion method for continuous change of conductivity block. Figure 1 、 Figure 2 and Figure 3 As shown, the following steps are included:

[0075] S1: The idea of ​​Tikhonov regularized inversion is introduced, and a model constraint term based on the conductivity of the unstructured tetrahedral grid nodes is added to the objective function. The inversion objective function expression based on the node model parameters is obtained, including the data fitting difference term and the three-dimensional model constraint term.

[0076] In the embodiment of the present invention, the idea of ​​Tikhonov regularized inversion is introduced, and a model constraint term based on the conductivity of the unstructured tetrahedral grid nodes is added to the objective function. The inversion objective function expression based on the node model parameters is obtained as follows:

[0077]

[0078] Among them, in the above formula, The data fitting difference term is used to fit the observed data; It is a model constraint item, providing the model’s prior information so that the model maintains a certain degree of smoothness; m=[m1,m2,…,m M ] Tis the inversion model parameter vector, that is, the conductivity model parameter vector established on the tetrahedral mesh nodes; M is the number of model parameters, and λ is the regularization factor used to balance the data fitting term and the model constraint term to ensure that these two terms are not overfitted; the expression of the data fitting difference term is as follows:

[0079]

[0080] where d=[d1,d2,…,d N ] T is the observation data vector, N is the size of the observation data, d can be the full tensor impedance data or the full tensor impedance data and the dipole vector impedance, and F(m) is the forward response. d is an N×N order data weighting matrix containing data errors, which represents the weight of each data in the inversion process. The added noise is expected to be 0 and the standard deviation is ε i Random numbers, that is, Gaussian noise;

[0081] The data weighting matrix is ​​written as W d =diag{1 / ε1,1 / ε2,…,1 / ε N}, when the standard deviation ε i The larger it is, the smaller its weight in the data weighting matrix is. i The smaller it is, the greater its weight in the data weighting matrix is. The full tensor impedance response is selected as the observation data, where the impedance tensor response expression is:

[0082]

[0083] In the formula, superscripts 1 and 2 represent two polarization directions.

[0084] S2: Minimize the norm of the conductivity change between the unit and the surrounding units, constrain the smoothness of the model, solve the model constraints for the node resistivity parameters, and obtain the model constraints based on the conductivity of the tetrahedral mesh nodes.

[0085] In an embodiment of the present invention, in a continuous space, the model constraint term based on the L2 norm is written as:

[0086]

[0087] The solution to the discretized model constraints is:

[0088]

[0089] Among them, C m is the model covariance matrix, constraining the model smoothness, W m is the model weight matrix.

[0090] S3: A three-dimensional network model is directly constructed by using model constraints based on the conductivity of tetrahedral mesh nodes, and three-dimensional model constraints at the nodes of the three-dimensional network model are calculated.

[0091] In an embodiment of the present invention, the model parameters at a node are obtained, and all adjacent tetrahedral units containing this node are found; the model constraint term of each tetrahedral unit is obtained by summing the product of the model parameters of all tetrahedral units containing the node and the gradient of the linear interpolation basis function; after summing the model constraint terms of all adjacent inversion units at the node, the model constraint term of the node is obtained by taking the average value. The expression of the model constraint term based on the unit node conductivity parameter is:

[0092]

[0093] Where M is the number of tetrahedral unit nodes in the underground area, V i is the volume of the i-th bottom tetrahedron, and are the conductivity parameter of the j-th node of unit i and the gradient of the linear interpolation basis function, respectively, and L2 represents the second norm of the orientation quantity.

[0094] S4: Gradient solve the inversion objective function expression based on the data fitting difference term and the 3D model constraint term.

[0095] In an embodiment of the present invention, a priori model is provided and added to the constraints of the three-dimensional model:

[0096]

[0097] Then the objective function of magnetotelluric inversion is expressed as:

[0098]

[0099] where m r is the prior model, is the data covariance matrix; the objective function is partial derivative with respect to m, and Further solve and transform into matrix form to obtain the forward linear equations; calculate the partial derivative of m on both sides of the forward linear equations, and combine the expression of the sparse stiffness matrix in to further deduce The specific calculation formula of the component is:

[0100]

[0101] Among them, p and q represent the row index and column index of the element respectively, j is a node in the inversion area, σ jrepresents the conductivity of the inversion unit node j, Ω j is the area covered by all inversion units that share the same vertex with inversion unit node j, k represents the serial number of the inversion unit that shares the same vertex with inversion node j, Ω k Indicates the area where the kth inversion unit is located, L j It is the interpolation basis function of the nodes at the common vertices of the inversion unit.

[0102] According to the upper and lower limits of conductivity Where a and b are constants, simplifying them yields:

[0103]

[0104] so It can be further deduced as:

[0105]

[0106] Among them, p and q represent the row index and column index of the element respectively, j is a node in the inversion area, σ j represents the conductivity of the inversion unit node j, Ω j is the area covered by all inversion units that share the same vertex with the inversion unit node j, N(j) represents the number of inversion units that share the same vertex with the inversion unit node j, k represents the sequence number of the inversion unit that shares the same vertex with the inversion node j, Ω k Indicates the area where the kth inversion unit is located, L j is the interpolation basis function of the node at the common vertex of the inversion unit; in the kth inversion unit Ω that shares the same vertex with the inversion unit node j k In the equation, the gradient of the stiffness matrix can be expressed as:

[0107]

[0108] In the formula, f ij =a i a j +b i b j +c i c j ;

[0109] According to the unit integration formula of volume coordinates:

[0110]

[0111] The kth inversion unit contains four nodes, and the expression of the stiffness matrix gradient will be different after the inversion unit node j falls on any node; when the node numbers of the inversion unit node j in the tetrahedron unit are 1, 2, 3, and 4 respectively, the L in the gradient integral term of the stiffness matrix is jThe expressions are L1, L2, L3, and L4 respectively.

[0112] S5: Solve the specific expression of each element in the stiffness matrix gradient under four different conditions, and solve the gradient components of the stiffness matrix corresponding to all nodes, thereby realizing the solution of the objective function gradient in magnetotelluric inversion.

[0113] In the embodiment of the present invention, the complex sensitivity matrix is ​​converted into the form of a real sensitivity matrix and then the formula of the complex sensitivity matrix is ​​substituted into the calculation formula of the objective function fitting difference gradient to solve the objective function gradient; the stiffness matrix gradient D under four different conditions is solved. ij The specific expression of each element in is obtained; the gradients of the stiffness matrices of the N(j) inversion units that share a vertex with the inversion unit node j are summed to obtain the gradient of the stiffness matrix at the inversion unit node j; the gradients of the stiffness matrices at all inversion unit nodes are solved to realize the solution of the gradient of the objective function in magnetotelluric inversion.

[0114] Example 1

[0115] The present application embodiment takes a three-dimensional flat terrain model as an example. Figure 4 As shown in the figure, the dimensions of the two high-resistance and low-resistance cubes in the model are both 10×10×5 km, the top surface is buried at a depth of 2 km, the conductivity is 1000Ω·m and 10Ω·m respectively, and the background conductivity is 100Ω·m. The measuring points are evenly distributed on the surface [-15km, 15km] 2 In the horizontal range, the point distance in the horizontal direction is 2.5 km, the total number of measurement points is 169, and the observation frequency is logarithmic (common logarithmic) and evenly spaced at 7 frequency points from 0.01 to 10 Hz. For the real and imaginary parts of each component of the full tensor impedance, the standard deviation is added as Gaussian random noise with a standard deviation of 0.005 is added to the inclination data as synthetic data. The forward calculation area is set to Ω = [-130km, 130km] 3 , the inversion area is the entire underground forward modeling area. At the same time, the traditional block uniform inversion algorithm and the conductivity block continuous change inversion algorithm of this paper are used to carry out the inversion calculation of the model. The initialization model and parameters are consistent under the inversion of the two algorithms. The uniform half-space model is set as the initial model of the inversion, the conductivity is 100Ωm, the initial regularization factor is 10, the cooling factor is 1.2, the number of continuous cooling times is 20 times, and the minimum allowable value is 10 -14The RMS target value is set to 0.995; the maximum number of inversion iterations is set to 160. The open source Tetgen software is used to generate the inversion mesh for the model. The meshes under the two inversion algorithms are consistent, with 451048 grid cells and 72415 nodes. The local refinement of the inversion mesh near the measurement point is shown in the figure. Figure 5 shown.

[0116] At the end of the inversion, the parameters of the two algorithms are shown in Table 1. After 84 inversion iterations, the RMS value of the conductivity block uniform inversion algorithm dropped to 0.99433, and the inversion iteration process ended. The initial value of the regularization factor was set to 10, and at the end of the iteration, the regularization factor dropped to 4.12739×10 -5 The number of inversion unknowns was 291,432, and the entire iterative process took 5.37 hours. When using the block-by-block continuous variation conductivity inversion algorithm, after 27 inversion iterations, the number of consecutive cooling cycles of the regularization factor reached the preset maximum value, terminating the inversion iteration. The regularization factor was initially set to 10, and by the end of the iteration, it had dropped to 0.37561, with an RMS value of 1.00647. The number of inversion unknowns was 55,696, and the entire iterative process took 1.49 hours.

[0117] Table 1 Parameters corresponding to the two algorithms at the end of inversion

[0118]

[0119] Figure 6 is a comparison of the cooling curves of the regularization factor during the inversion iteration process. Figure 7 It is a comparison of the RMS convergence curve during the inversion iteration process.

[0120] Figure 8 The inversion results for the high- and low-resistance anomaly model on flat terrain are shown in Figure 2. A comparative analysis of the inversion results using the conductivity block uniform algorithm and the conductivity block continuous variation algorithm reveals that both algorithms perform well in inverting the two anomalies, accurately retrieving their size and location. Low-resistance anomalies with conductivity close to 10 Ω·m can be accurately inverted using both algorithms. For high-resistance anomalies, the conductivity block uniform algorithm can invert high-resistance anomalies with conductivity around 300 Ω·m, while the conductivity block continuous variation algorithm can invert high-resistance anomalies with conductivity above 500 Ω·m, achieving significantly better inversion recovery.

[0121] Example 2

[0122] This embodiment takes the terrain model as an example. Figure 9A low-resistance anomaly with a resistance of 10 Ω·m is located directly beneath the valley, with its top interface 3 km vertically from the valley's lowest point. A high-resistance anomaly with a resistance of 2000 Ω·m is located directly beneath the peak, with its top interface 6 km vertically from the peak's summit. The background conductivity is 100 Ω·m. Both anomalies measure 10 × 10 × 5 km. There are 40 measurement points, evenly distributed across the surface in a horizontal range of -12 km to 12 km in the X direction and -21 km to 21 km in the Y direction, with a spacing of 6 km. There are 16 observation frequencies, which are equally spaced in the range of 0.001 to 10 Hz, namely 0.001 Hz, 0.00185 Hz, 0.00342 Hz, 0.00631 Hz, 0.0117 Hz, 0.0216 Hz, 0.0398 Hz, 0.0736 Hz, 0.1360 Hz, 0.2512 Hz, 0.4642 Hz, 0.8577 Hz, 1.5849 Hz, 2.9287 Hz, 5.4117 Hz and 10 Hz. The standard deviation of the real and imaginary data of each component of the full tensor impedance is added to the observation data. random noise to be used as synthetic data.

[0123] The inversion calculation of the model is carried out using both the traditional block uniform inversion algorithm and the conductivity block continuous variation inversion algorithm proposed in this paper. The initialization model and parameters are consistent under the two algorithms. The uniform half-space model is set as the initial inversion model, the conductivity is 100Ω·m, the initial regularization factor is 10, the cooling factor is 1.2, and the minimum allowable value is 10 -14 The RMS target value is set to 1.005; the maximum number of inversion iterations is set to 160. The open source Tetgen software is used to perform inversion meshing on the model. The meshes under the two inversion algorithms are consistent, with 463356 grid cells and 72894 nodes. The inversion area during the inversion process is the underground area of ​​the entire model. The local densification of the inversion mesh near the measuring point is shown in the figure below. Figure 10 shown.

[0124] At the end of the inversion, the parameters of the two algorithms are shown in Table 2. After 137 inversion iterations, the RMS value of the conductivity block uniform inversion algorithm dropped to 1.00488, and the inversion iteration ended. The initial value of the regularization factor at the beginning of the iteration was 10. When the iteration ended, the regularization factor dropped to 7.32488×10 -10, the number of inversion unknowns is 398073, and the entire iterative process takes 17.75 hours. After 99 inversion iterations, the RMS value of the conductivity block continuous change inversion algorithm dropped to 1.00489, and the inversion iteration process ended. The initial value of the regularization factor is 10, and the regularization factor dropped to 5.19189×10 at the end of the iteration. -7 The number of inversion unknowns is 65890, and the entire iterative process takes 8.77 hours.

[0125] Table 2 Parameters corresponding to the two algorithms at the end of inversion

[0126]

[0127] Figure 11 is a comparison of the cooling curves of the regularization factor during the inversion iteration process. Figure 12 It is a comparison of the RMS convergence curve during the inversion iteration process. Figure 13 Figure 1 shows the inversion results for a model of high and low resistivity anomalies under terrain considerations. The inversion results using the conductivity block uniform algorithm and the conductivity block continuous variation algorithm are compared and analyzed, with slices of the inversion results at x = 0 km and z = 7 km analyzed. It can be seen that the conductivity block continuous inversion algorithm achieves significantly better inversion results than the conductivity block uniform algorithm, recovering conductivity values ​​down to 10 Ω·m for low-resistance anomalies and even inverting conductivity values ​​above 200 Ω·m for high-resistance anomalies. Overall, when considering the effects of terrain, the conductivity block continuous inversion algorithm better recovers high and low resistivity anomalies in the subsurface than the conductivity block uniform inversion algorithm. The inversion iteration process generates fewer unknowns, significantly shortening the inversion iteration time.

[0128] The beneficial effects of the present invention are embodied in the following aspects: the specific content of the present invention is to address the problem that in traditional magnetotelluric three-dimensional inversion, the conductivity parameters are established on tetrahedral units. Since the inversion unknowns are the number of model units, and the corresponding number of grid units in the magnetotelluric three-dimensional inversion is huge, which leads to a large number of inversion unknowns and a strong multi-solution problem. Therefore, it is proposed to adopt a parameterized assignment method in magnetotelluric inversion that establishes the conductivity on the tetrahedral unit nodes. Because the number of tetrahedral unit nodes in an unstructured grid is much smaller than the number of tetrahedral units, this conductivity assignment method can greatly reduce the inversion unknowns, improve the uniqueness of the solution, and adapt to the situation where the conductivity of the underground medium changes continuously, thereby obtaining a better inversion effect.

[0129] The above description is merely an embodiment of the present invention. Common knowledge such as the specific technical solutions or features of the solutions is not described in detail here. It should be noted that those skilled in the art may make several modifications and improvements without departing from the solution of the present invention, and these modifications and improvements should also be considered as the scope of protection of the present invention. These modifications and improvements will not affect the effects of the present invention and the practicality of the patent. The scope of protection claimed in this application shall be based on the content of the claims, and the specific embodiments and other descriptions in the specification may be used to interpret the content of the claims.

Claims

1. A three-dimensional magnetotelluric inversion method for continuous block variation of conductivity, characterized by: The method comprises: S1: The idea of ​​Tikhonov regularized inversion is introduced, and a model constraint term based on the conductivity of the unstructured tetrahedral grid nodes is added to the objective function. The inversion objective function expression based on the node model parameters is obtained, including the data fitting error term and the three-dimensional model constraint term; S2: Minimize the norm of the conductivity change between the unit and the surrounding units, constrain the smoothness of the model, solve the model constraint terms for the node resistivity parameters, and obtain the model constraint terms based on the conductivity of the tetrahedral mesh nodes; S3: directly constructing a three-dimensional network model through the model constraint items based on the conductivity of the tetrahedral mesh nodes, and calculating the three-dimensional model constraint items at the nodes of the three-dimensional network model; S4: performing a gradient solution on the inversion objective function expression according to the data fitting difference term and the three-dimensional model constraint term; S5: Solve the specific expression of each element in the stiffness matrix gradient under four different conditions, and solve the gradient components of the stiffness matrix corresponding to all nodes, thereby solving the gradient of the objective function in the magnetotelluric inversion; The idea of ​​Tikhonov regularized inversion is introduced, and the model constraint term based on the conductivity of the unstructured tetrahedral grid nodes is added to the objective function. The inversion objective function expression based on the node model parameters is obtained as follows: Among them, in the above formula, The data fitting difference term is used to fit the observed data; It is a model constraint item, providing the model’s prior information so that the model maintains a certain degree of smoothness; m=[m1,m2,…,m M ] T is the inversion model parameter vector, that is, the conductivity model parameter vector established on the tetrahedral grid nodes; M is the number of model parameters, and λ is the regularization factor used to balance the data fitting term and the model constraint term; The expression of the data fitting difference term is as follows: where d=[d1,d2,…,d N ] T is the observation data vector, N is the size of the observation data, d is the full tensor impedance data or full tensor impedance data and dipole vector impedance, F(m) is the forward response, W d is an N×N order data weighting matrix containing data errors, which represents the weight of each data in the inversion process. The added noise is expected to be 0 and the standard deviation is ε i Random numbers, that is, Gaussian noise; The data weighting matrix is ​​written as W d =diag{1 / ε1,1 / ε2,…,1 / ε N }, when the standard deviation ε i The larger it is, the smaller its weight in the data weighting matrix is. i The smaller it is, the greater its weight in the data weighting matrix is. The full tensor impedance response is selected as the observation data, where the impedance tensor response expression is: In the formula, superscripts 1 and 2 represent two polarization directions; The steps of minimizing the norm of the conductivity variation between the unit and the surrounding units, constraining the smoothness of the model, solving the model constraint items for the node resistivity parameters, and obtaining the model constraint items based on the tetrahedral mesh node conductivity include: In continuous space, the model constraint based on the L2 norm is written as: The solution to the discretized model constraints is: Among them, C m is the model covariance matrix, constraining the model smoothness, W m is the model weight matrix; The step of directly constructing a three-dimensional network model by using the model constraint items based on the conductivity of the tetrahedral mesh nodes and calculating the three-dimensional model constraint items at the nodes of the three-dimensional network model includes: Get the model parameters at a node and find all adjacent tetrahedral elements containing this node; Obtain a model constraint item of each tetrahedral unit by summing the product of the model parameters of all tetrahedral units including the node and the gradient of the linear interpolation basis function; After summing the model constraints of all adjacent inversion units at the node, the average value is taken to obtain the model constraint term of the node. The expression of the model constraint term based on the unit node conductivity parameter is: Where M is the number of tetrahedral unit nodes in the underground area, V i is the volume of the i-th bottom tetrahedron, and are the conductivity parameter of the j-th node of unit i and the gradient of the linear interpolation basis function, respectively, and L2 represents the second norm of the orientation quantity.

2. The method for three-dimensional magnetotelluric inversion of conductivity block-by-block continuous variation according to claim 1, characterized in that: The step of performing a gradient solution on the inversion objective function expression according to the data fitting difference term and the three-dimensional model constraint term comprises: Provide a priori model and add it to the constraints of the three-dimensional model: Then the objective function of magnetotelluric inversion is expressed as: where m r is the prior model, is the data covariance matrix; The objective function is partial derivative with respect to m, and Further solve and convert into matrix form to obtain the forward linear equations; The partial derivatives of m on both sides of the forward linear equations are obtained, and combined with the expression of the sparse stiffness matrix in, we can further deduce The specific calculation formula of the component is: Among them, p and q represent the row index and column index of the element respectively, j is a node in the inversion area, σ j represents the conductivity of the inversion unit node j, Ω j is the area covered by all inversion units that share the same vertex with inversion unit node j, k represents the serial number of the inversion unit that shares the same vertex with inversion node j, Ω k Indicates the area where the kth inversion unit is located, L j It is the interpolation basis function of the nodes at the common vertices of the inversion unit.

3. The method for three-dimensional magnetotelluric inversion with continuous block variation of conductivity according to claim 2, characterized in that: According to the upper and lower limits of conductivity Where a and b are constants, simplifying them yields: so It can be further deduced as: Among them, p and q represent the row index and column index of the element respectively, j is a node in the inversion area, σ j represents the conductivity of the inversion unit node j, Ω j is the area covered by all inversion units that share the same vertex with the inversion unit node j, N(j) represents the number of inversion units that share the same vertex with the inversion unit node j, k represents the sequence number of the inversion unit that shares the same vertex with the inversion node j, Ω k Indicates the area where the kth inversion unit is located, L j It is the interpolation basis function of the nodes at the common vertices of the inversion unit.

4. The method for three-dimensional magnetotelluric inversion with block-by-block continuous conductivity variation according to claim 3, characterized in that: After performing gradient solving on the inversion objective function expression, the method further includes: solving the gradient of the stiffness matrix of the inversion unit of the common vertex of the inversion unit nodes according to the unit integration formula of the volume coordinate; In the kth inversion unit Ω that shares the same vertex with the inversion unit node j k In the equation, the gradient of the stiffness matrix can be expressed as: In the formula, f ij =a i a j +b i b j +c i c j ; According to the unit integration formula of volume coordinates: The kth inversion unit contains four nodes, and the expression of the stiffness matrix gradient will be different after the inversion unit node j falls on any node; when the node numbers of the inversion unit node j in the tetrahedron unit are 1, 2, 3, and 4 respectively, the L in the gradient integral term of the stiffness matrix is j The expressions are L1, L2, L3, and L4 respectively.

5. The method for three-dimensional magnetotelluric inversion of conductivity block-by-block continuous variation according to claim 4, characterized in that: After solving the stiffness matrix gradient of the inversion unit with the same vertex as the inversion unit node according to the unit integration formula of the volume coordinate, it also includes: The complex sensitivity matrix is ​​converted into the form of a real sensitivity matrix and then the formula of the complex sensitivity matrix is ​​substituted into the calculation formula of the objective function fitting difference gradient to solve the objective function gradient.

6. The method for three-dimensional magnetotelluric inversion of conductivity block-by-block continuous variation according to claim 5, characterized in that: The steps of solving the specific expression of each element in the stiffness matrix gradient under four different conditions and solving the gradient components of the stiffness matrix corresponding to all nodes, thereby solving the gradient of the objective function in the magnetotelluric inversion, include: Solve the stiffness matrix gradient D under four different conditions ij The specific expression of each element in; Calculate the gradient of the stiffness matrix of N(j) inversion units that share the same vertex with the inversion unit node j and sum them to obtain the gradient of the stiffness matrix at the inversion unit node j; The gradient of the stiffness matrix at all inversion unit nodes is solved to achieve the solution of the gradient of the objective function in magnetotelluric inversion.

Images