Magnetotelluric three-dimensional inversion method with conductivity continuously changing in blocks

Abstract

The invention relates to the field of geophysics, and discloses a conductivity block continuous change magnetotelluric three-dimensional inversion method, and the method comprises the steps: adding a model constraint term based on the conductivity of a non-structural tetrahedral mesh node into a target function, and obtaining an inversion target function expression; enabling the norm of the conductivity variation of the unit and the surrounding unit to be minimum, and solving a model constraint term to obtain the model constraint term based on the tetrahedral mesh node conductivity; carrying out gradient solution on the inversion objective function expression; specific expressions of elements in stiffness matrix gradients under four different conditions are solved, gradient components of stiffness matrixes corresponding to all nodes are solved, and then the objective function gradient is solved in magnetotelluric inversion; according to the method, updating of the node conductivity model during three-dimensional underground medium inversion is achieved, unknown numbers generated in the inversion process are greatly reduced, the uniqueness of a solution is improved, meanwhile, the time needed by the inversion process is greatly shortened, and a better inversion effect is obtained.

Description

Technical Field

[0001] The invention relates to the field of geophysics, and in particular to a three-dimensional magnetotelluric inversion method for block-by-block continuous variation of 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 magnetotelluric three-dimensional inversion, conductivity parameters are usually established on tetrahedral units, that is, a grid unit has a fixed conductivity value. Since the inversion unknowns are the number of model units, and the corresponding number of grid units in magnetotelluric three-dimensional inversion is huge, there are many inversion unknowns and the solution is highly multi-solutionable.

[0003] The present invention adopts a parameterized assignment method of establishing conductivity on tetrahedral unit nodes in magnetotelluric inversion. Since 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 magnetotelluric three-dimensional inversion algorithm based on 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 a model term 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 a better inversion effect.

[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 new stiffness matrix gradients, and is characterized in that the method includes:

[0008] S1: 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, and 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;

[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 items for the node resistivity parameters, and obtain the model constraint items 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, so as to realize 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 difference term is used to fit the observed data; is the model constraint, 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 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 The random number is Gaussian noise; the data weighting matrix can be written as W d =diag{1 / ε1,1 / ε2,…,1 / ε N}, when the standard deviation ε iThe larger the value, the smaller its weight in the data weighting matrix. i The smaller it is, the greater its weight in the data weighting matrix is, and the full tensor impedance response is selected as the observed 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 items for the node resistivity parameters, and obtaining the model constraint items based on the conductivity of the tetrahedral mesh nodes comprises:

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

[0023]

[0024] The solution of 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 mesh 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 ith bottom tetrahedron, and are the conductivity parameter of the jth 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 expression of the objective function of magnetotelluric inversion is:

[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 solving and converting into matrix form, we get the forward linear equations;

[0037] The partial derivatives of m are calculated for both sides of the forward linear equations, 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] Preferably, according to the upper and lower limits of conductivity Where a and b are constants, simplifying them yields:

[0040]

[0041] so It can be further deduced as:

[0042]

[0043] Among them, p and q represent the row index and column index corresponding to the element, 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, N(j) represents the number of inversion units that share the same vertex with inversion unit node j, k represents the sequence number of the inversion unit that shares the same vertex with inversion node j, Ω k represents 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.

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

[0045]

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

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

[0048]

[0049] 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 tetrahedral unit are 1, 2, 3, and 4 respectively, the gradient integral term of the stiffness matrix is ​​L j The expressions are L1, L2, L3, and L4 respectively.

[0050] 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, it also includes:

[0051] 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.

[0052] 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, comprises: 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 the same 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 achieve the solution of the gradient of the objective function in magnetotelluric inversion.

[0053] The beneficial effects of the present invention are embodied in the following aspects: the specific content of the present invention is to solve the problem that in the traditional three-dimensional magnetotelluric 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 three-dimensional magnetotelluric 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 of establishing the conductivity on the tetrahedral unit nodes in the magnetotelluric inversion. Because the number of tetrahedral unit nodes in the 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 of continuous change of the conductivity of the underground medium to obtain a better inversion effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] In order to more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the following is a brief introduction to the drawings required for the specific embodiments or the description of the prior art. In all the drawings, similar elements or parts are generally identified by similar reference numerals. In the drawings, the elements or parts are not necessarily drawn according to the actual scale.

[0055] Figure 1 It is a flow chart of a method for three-dimensional magnetotelluric inversion with continuous block variation of conductivity provided by an embodiment of the present invention.

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

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

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

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

[0060] Figure 6 It 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.

[0061] Figure 7 It is a comparison diagram of the RMS convergence curve during the inversion iteration process provided in Example 1 of the present invention.

[0062] 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.

[0063] Fig. 9 It is a schematic diagram of the high and low resistance anomaly model under the valley terrain provided by Example 2 of the present invention.

[0064] Fig.10 This is the local encryption of the inversion grid near the measuring point provided by Example 2 of the present invention.

[0065] Fig.11 It 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.

[0066] Fig.12 It is a comparison diagram of the RMS convergence curve during the inversion iteration process provided in Example 2 of the present invention.

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

[0068] In order to make the technical personnel in the technical field better understand the scheme of the present invention, the technical scheme in the embodiment of the present invention will be clearly and completely described below in combination with the drawings in the embodiment of the present invention. Obviously, the described embodiment is only a part of the embodiment of the present invention, not all the embodiments. Based on the embodiment of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0069] The terms "first", "second", etc. in the specification and claims of the present invention and the above-mentioned drawings are used to distinguish different objects rather than to describe a specific order. In addition, the terms "including" and "having" and any variations thereof are intended to cover non-exclusive inclusions. For example, a process, method, device, product or end including a series of steps or units is not limited to the listed steps or units, but may optionally include steps or units that are not listed, or may optionally include other steps or units inherent to these processes, methods, products or ends.

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

[0071] The present invention adopts a parameterized assignment method of establishing conductivity on tetrahedral unit nodes in magnetotelluric inversion. Since 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 magnetotelluric three-dimensional inversion algorithm based on continuous change of conductivity block is proposed.

[0072] The present invention aims to provide a three-dimensional magnetotelluric inversion method with block-by-block continuous variation of conductivity. By constructing a model term 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 a better inversion effect.

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

[0074] S1: 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, including the data fitting difference term and the three-dimensional model constraint term.

[0075] 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, and the inversion objective function expression based on the node model parameters is obtained as follows:

[0076]

[0077] Among them, in the above formula, The difference term is used to fit the observed data; is the model constraint, 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 mesh nodes; M is the number of model parameters, λ is the regularization factor, which is 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:

[0078]

[0079] where d=[d1,d2,…,dN ] 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 A random number, i.e. Gaussian noise;

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

[0081]

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

[0083] 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.

[0084] In the embodiment of the present invention, in the continuous space, the model constraint term based on the L2 norm can be written as:

[0085]

[0086] The solution of the discretized model constraints is:

[0087]

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

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

[0090] In an embodiment of the present invention, the model parameters at a node are obtained, and all adjacent tetrahedral units containing the node are found; the model constraint items of each tetrahedral unit are 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 items of all adjacent inversion units at the node, the model constraint items of the node are obtained by taking the average value. The expression of the model constraint term based on the unit node conductivity parameter is:

[0091]

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

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

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

[0095]

[0096] Then the expression of the objective function of magnetotelluric inversion is:

[0097]

[0098] 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 derivatives of m on both sides of the forward linear equations, and combine with the expression of the sparse stiffness matrix in to further deduce The specific calculation formula of the component is:

[0099]

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

[0101]

[0102] so It can be further deduced as:

[0103]

[0104] Among them, p and q represent the row index and column index corresponding to the element, 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, N(j) represents the number of inversion units that share the same vertex with inversion unit node j, k represents the sequence number of the inversion unit that shares the same vertex with inversion node j, Ω k represents 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:

[0105]

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

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

[0108]

[0109] 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 tetrahedral unit are 1, 2, 3, and 4 respectively, the gradient integral term of the stiffness matrix is ​​L j The expressions are L1, L2, L3,

[0110] L4.

[0111] 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, so as to realize the solution of the objective function gradient in magnetotelluric inversion.

[0112] 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. ijThe specific expression of each element in is obtained; the gradients of the stiffness matrices of the N(j) inversion units that share the same 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 achieve the solution of the gradient of the objective function in magnetotelluric inversion.

[0113] Example 1

[0114] The present application embodiment takes a three-dimensional flat terrain model as an example. Figure 4 As shown in the figure, the size of the two high and low resistance cubes in the model are 10×10×5km, the top surface is buried at a depth of 2km, 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 Within the horizontal range, the point distance in the horizontal direction is 2.5km, the total number of measurement points is 169, and the observation frequency is logarithmic (common logarithmic) and evenly distributed at 7 frequency points of 0.01-10Hz. For the real and imaginary parts of each component of the full tensor impedance, a standard deviation of 5% is added Gaussian random noise with a standard deviation of 0.005 is added to the tilt sub-data as synthetic data. Set the forward calculation area 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 in 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 continuous cooling number is 20 times, and the minimum allowable value is 10 -14 ; The 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 perform inversion meshing on the model. The meshes under the two inversion algorithms are consistent, with 451048 grid cells and 72415 nodes. The local encryption of the inversion mesh near the measuring point is shown in Figure 2. Figure 5 shown.

[0115] 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 the regularization factor dropped to 4.12739×10 -5, the number of inversion unknowns is 291432, and the entire iteration process takes 5.37 hours. When the conductivity block continuous change inversion algorithm is used for inversion, after 27 inversion iterations, the number of continuous cooling times of the regularization factor reaches the preset maximum value, and the inversion iteration process ends. The initial value of the regularization factor is set to 10. At the end of the iteration, the regularization factor drops to 0.37561, and the RMS value is 1.00647. The number of inversion unknowns is 55696, and the entire iteration process takes 1.49 hours.

[0116] Table 1 Parameters of the two algorithms at the end of inversion

[0117]

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

[0119] Figure 8 This is the inversion result diagram of the high and low resistivity anomaly model under flat terrain. The inversion results under the conductivity block uniform algorithm and the conductivity block continuous change algorithm are compared and analyzed. From the sections of the inversion results at y=0km and z=4.5km, it can be seen that the inversion effects of the two anomalies under the two algorithms are relatively good, and the size and position are inverted relatively accurately. Both can well invert the low-resistance anomaly with a conductivity close to 10Ω@m. For the inversion of high-resistance anomalies, the conductivity block uniform algorithm can invert the high-resistance anomaly with a conductivity of about 300Ω@m, while the conductivity block continuous change algorithm can invert the high-resistance anomaly with a conductivity of more than 500Ω@m, and the inversion recovery effect is significantly better.

[0120] Example 2

[0121] This embodiment takes the terrain model as an example. Fig. 9 It can be seen that there is a 10Ω·m low-resistance anomaly directly below the valley terrain, and the vertical distance between its upper interface and the lowest point of the valley is 3km. There is a 2000Ω@m high-resistance anomaly directly below the peak terrain, and the vertical distance between its upper interface and the peak is 6km. The background conductivity is 100Ω@m. The size of the two anomalies is 10×10×5km. The total number of measurement points is 40, evenly distributed on the surface with a horizontal range of X:-12km~12km, Y:-21km~21km, and the point spacing is 6km. There are 16 observation frequencies, at 0.001 The 10Hz interval is evenly distributed in logarithmic intervals, namely 0.001Hz, 0.00185Hz, 0.00342Hz, 0.00631Hz, 0.0117Hz, 0.0216Hz, 0.0398Hz, 0.0736Hz, 0.1360Hz, 0.2512Hz, 0.4642Hz, 0.8577Hz, 1.5849Hz, 2.9287Hz, 5.4117Hz and 10Hz. The standard deviation of 3% is added to the real and imaginary data of each component in the full tensor impedance in the observed data. random noise to be used as synthetic data.

[0122] The inversion calculation of the model is carried out by using the traditional block uniform inversion algorithm and the conductivity block continuous change inversion algorithm in this paper. The initialization model and parameters are consistent under the inversion of 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 model is meshed using the open source Tetgen software. 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 encryption of the inversion grid near the measuring point is shown in Figure 2. Fig.10 shown.

[0123] 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 was 10 at the beginning of the iteration. When the iteration ended, the regularization factor dropped to 7.32488×10 -10 , the number of inversion unknowns is 398073, and the entire iteration process takes 17.75 hours. After 99 inversion iterations of the conductivity block continuous change inversion algorithm, the RMS value dropped to 1.00489, and the inversion iteration process ended. The initial value of the regularization factor is 10, and the regularization factor drops to The number of inversion unknowns is 65890, and the entire iterative process takes 8.77 hours.

[0124] Table 2 Parameters of the two algorithms at the end of inversion

[0125]

[0126]

[0127] Fig.11 is a comparison of the cooling curves of the regularization factor during the inversion iteration process. Fig.12 It is a comparison chart of the RMS convergence curve during the inversion iteration process. Fig.13 This is the inversion result diagram of the high and low resistivity anomaly model under terrain. The inversion results of the conductivity block uniform algorithm and the conductivity block continuous change algorithm are compared and analyzed, and the sections of the inversion results at x = 0km and z = 7km are analyzed. It can be seen that the inversion effect of the conductivity block continuous inversion algorithm is much better than that of the conductivity block uniform algorithm. The conductivity of the low resistivity anomaly can be restored to 10Ω·m, and the high resistivity anomaly can also invert the conductivity area above 200Ω@m. In general, when considering the influence of terrain, the conductivity block continuous inversion algorithm can better restore the high resistivity anomaly and low resistivity anomaly for the underground than the conductivity block uniform inversion algorithm. The inversion iterative process generates fewer unknowns, which greatly shortens the time required for the inversion iterative process.

[0128] The beneficial effects of the present invention are embodied in the following aspects: the specific content of the present invention is to solve the problem that in the traditional three-dimensional magnetotelluric 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 three-dimensional magnetotelluric 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 of establishing the conductivity on the tetrahedral unit nodes in the magnetotelluric inversion. Because the number of tetrahedral unit nodes in the 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 of continuous change of the conductivity of the underground medium to obtain a better inversion effect.

[0129] The above is only an embodiment of the present invention, and the common knowledge such as the known specific technical schemes or characteristics in the scheme is not described in detail here; it should be pointed out that for those skilled in the art, several modifications and improvements can be made without departing from the scheme of the present invention, which should also be regarded as the protection scope of the present invention, and these will not affect the effect of the present invention and the practicality of the patent. The protection scope required by this application shall be based on the content of its claims, and the specific implementation methods and other records in the specification can be used to interpret the content of the claims.

Claims

1. A three-dimensional magnetotelluric inversion method for continuous change of conductivity in blocks, characterized in that: The method comprises: S1: 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, and 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; 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 items for the node resistivity parameters, and obtain the model constraint items 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, so as to realize the solution of the objective function gradient in magnetotelluric inversion.

2. The method for three-dimensional magnetotelluric inversion of conductivity block-by-block continuous variation according to claim 1, characterized in that: 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 difference term is used to fit the observed data; is the model constraint, 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.

3. The method for three-dimensional magnetotelluric inversion of conductivity block-by-block continuous variation according to claim 2 is characterized by: 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 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 A random number, i.e. Gaussian noise; The data weighting matrix can be written as W d =diag{1 / ε1,1 / ε2,…,1 / ε N }, when the standard deviation ε i The larger the value, the smaller its weight in the data weighting matrix. i The smaller it is, the greater its weight in the data weighting matrix is, and the full tensor impedance response is selected as the observed data, where the impedance tensor response expression is: In the formula, superscripts 1 and 2 represent two polarization directions.

4. The method for three-dimensional magnetotelluric inversion of conductivity block-by-block continuous variation according to claim 3 is characterized in that: 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 conductivity of the tetrahedral mesh nodes include: In continuous space, the model constraint based on the L2 norm can be written as: The solution of the discretized model constraints is: Among them, C m is the model covariance matrix, which constrains the model smoothness, W m is the model weight matrix.

5. The method for three-dimensional magnetotelluric inversion of conductivity block-by-block continuous variation according to claim 4, characterized in that: The step of 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 includes: Get the model parameters at a node and find all adjacent tetrahedral elements containing this node; The model constraint item of each tetrahedral unit is obtained 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 up 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 The conductivity parameters of the jth node of unit i and the gradient of the linear interpolation basis function, respectively, L2 represents the two-normal number of the oriented vector.

6. A three-dimensional inversion method for conductivity continuous change in earth electromagnetic inversion according to claim 5, 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 expression of the objective function of magnetotelluric inversion is: Where m r is the prior model, is the data covariance matrix; The objective function is partial derivative with respect to m, and Further solving and converting into matrix form, we get the forward linear equations; The partial derivatives of m are calculated for both sides of the forward linear equations, and combined with the expression of the sparse stiffness matrix in, we can further deduce The specific calculation formula of the component is:

7. The method for three-dimensional magnetotelluric inversion of conductivity block-by-block continuous variation according to claim 6, 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 corresponding to the element, 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, N(j) represents the number of inversion units that share the same vertex with inversion unit node j, k represents the sequence number of the inversion unit that shares the same vertex with inversion node j, Ω k represents 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.

8. The method for three-dimensional magnetotelluric inversion of conductivity block-by-block continuous variation according to claim 7, characterized in that: After performing gradient solving for 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 coordinates; 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 tetrahedral unit are 1, 2, 3, and 4 respectively, the gradient integral term of the stiffness matrix is ​​L j The expressions of L1, L2, L3, L 4。 9. The method for three-dimensional magnetotelluric inversion of conductivity block-by-block continuous variation according to claim 8, characterized in that: After solving the stiffness matrix gradient of the inversion unit with the same vertex of 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.

10. The method for three-dimensional magnetotelluric inversion of conductivity block-by-block continuous variation according to claim 9, characterized in that: 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 magnetotelluric inversion, includes: Solve the stiffness matrix gradient D for four different cases ij The specific expression of each element in ; Calculate the gradients of the stiffness matrices of N(j) inversion units that share the same vertex with the inversion unit node j and sum them up 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