Method for 3d anisotropy resistivity imaging based on non-structural finite elements

Abstract

The application discloses a kind of resistivity three-dimensional anisotropy imaging methods based on non-structure finite element, belong to the direct current method imaging technical field in geological exploration, comprising: constructing detection model and its anisotropy control equation in direct current field for anisotropic medium, wherein anisotropy control equation includes resistivity tensor, and is converted into finite element integral equation with boundary condition, using non-structure tetrahedron unit to carry out grid division discretization to detection model and calculate each unit node node potential, node potential is brought into finite element integral equation and is forward simulated;Regularization objective function is constructed, and in inversion regularization objective function, model term is converted into constraint domain related to gradient filter operator and upper and lower limit constraint, minimize inversion regularization objective function to update model parameter, obtain resistivity three-dimensional anisotropy imaging result, so as to ensure inversion accuracy, greatly reduce data volume and data acquisition cost.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of direct current imaging in geological exploration, and particularly relates to a resistivity three-dimensional anisotropy imaging method based on non-structural finite elements. BACKGROUND

[0002] At present, the underground medium is generally regarded as isotropic in the field of geophysical inversion, and this assumption has good applicability through long-term experiments. However, more and more studies show that the underground medium is not simply isotropic, and the results of inversion using isotropic theory have many false anomalies, which makes people get a wrong geological understanding. Anisotropy characteristics exist widely in nature, especially in places with developed bedding, fractures and underground water-rich cavities. These places often have resources or potential hazards, and are the areas that need to be focused on. The study of anisotropy has been started for decades, and has developed rapidly in the electromagnetic field, and has been widely used in resource exploration, potential geological hazards and other aspects.

[0003] Direct current resistivity method is a commonly used geophysical exploration method in the electromagnetic field. Due to the inherent characteristics of direct current method, the detection range is mostly the near-surface area. The introduction of anisotropy concept into inversion imaging will undoubtedly increase the imaging accuracy and has more practicality in resource exploration and engineering application. The study of direct current anisotropy has a long history, but due to the increase of inversion parameters, the non-uniqueness of inversion results is more obvious, which brings great difficulties and challenges to inversion work. In the early stage, people used simple models such as half-space or one-dimensional anisotropic layered medium to study the anisotropy law, and gave the analytical solution of the related model, which provided a reliable basis for the later study of direct current anisotropy. With the increase of medium layers, the model also tends to be complex. Although one-dimensional anisotropy inversion work is carried out, due to the inherent non-uniqueness of inversion and the increase of parameters, the results are no longer reliable.

[0004] With the increase of dimensions and the complexity of geological structure, people cannot directly get the corresponding solution, and have to use numerical simulation method for approximate solution. In the study of two-dimensional anisotropy, finite difference is used for forward modeling, and regularization factor is added in the inversion process to reduce the artifacts of inversion results and greatly improve the calculation speed.

[0005] Because the electrical anisotropy has directionality, the conventional two-dimensional data acquisition method has great limitations and cannot obtain effective anisotropy information, and therefore it is necessary to consider a three-dimensional observation device to bring higher quality data results and to carry out three-dimensional anisotropy research of direct current resistivity. In order to obtain better data information, well data is particularly important, and well data can better invert the current flow direction and reconstruct the geological structure. The use of unstructured grids can finely subdivide complex terrains, and a large number of researchers have carried out a large amount of theoretical research through well, well-ground and well-well data, analyzed the influence of different data combinations on anisotropy inversion results, and obtained better application in actual engineering, and the results are better than isotropy.

[0006] At present, non-destructive detection has become an important research hotspot, and when well-to-well measurement cannot be carried out in some specific areas, the non-destructive detection method becomes more important. This method does not need to destroy the geological structure and can save a lot of cost. However, it is difficult to invert the multiple parameters of anisotropy only with surface data. SUMMARY

[0007] In view of the above, the purpose of the present application is to provide a resistivity three-dimensional anisotropy imaging method based on a non-structural finite element, to solve the technical problem of poor resistivity three-dimensional anisotropy imaging effect in the current direct current resistivity anisotropy non-destructive detection.

[0008] To achieve the above-mentioned purpose of the application, the resistivity three-dimensional anisotropy imaging method based on the non-structural finite element provided by the embodiment comprises the following steps:

[0009] A detection model and its corresponding anisotropy control equation in a direct current electric field are constructed for anisotropic media in a detection area, wherein the anisotropy control equation contains a resistivity tensor, the anisotropy control equation is decomposed by using Galerkin, a finite element integral equation with a boundary condition is obtained by simultaneously introducing a Dirichlet boundary condition, a non-structural tetrahedral element is used to discretize and subdivide the detection model, and the node potential on each element node is calculated, and the node potential is brought into the finite element integral equation for forward modeling;

[0010] An inversion regularization objective function is constructed, the model term is converted into a constraint domain related to a gradient filter operator and upper and lower limit constraints in the inversion regularization objective function, the gradient calculation of the regularization objective function is realized through an adjoint forward modeling calculation and a gradient filtering calculation, the model parameters are updated by minimizing the inversion regularization objective function, and then the resistivity three-dimensional anisotropy imaging result is obtained.

[0011] Preferably, the anisotropy control equation is expressed as:

[0012] ;

[0013] wherein, denotes the electric potential, denotes the electric current, denotes the receiving point position, source point denotes the source point position, denotes the distance between the receiving point position and the source point position, denotes the resistivity tensor, symbol denotes the gradient, denotes the Dirac operator.

[0014] Preferably, the Galerkin method is used to decompose the anisotropic control equation, and the Dirichlet boundary condition is introduced to obtain the finite element integral equation with boundary conditions, which is expressed as:

[0015] ;

[0016] wherein, denotes the calculated detection area, superscript T denotes the transpose, denotes the electric potential U of the finite element integral equation function.

[0017] Preferably, the non-structural tetrahedral element is used to mesh the detection model for discretization and to calculate the node potential at each element node, including:

[0018] Any point p inside the tetrahedral element has a potential value of U p The line connecting point p and the four nodes divides the tetrahedron into four small tetrahedral elements. According to the definition of three-dimensional natural coordinates, the ratio of the volume of the small tetrahedral element composed of point p and the other three nodes excluding the k th node to the volume of the large tetrahedral element is the volume function L k When calculating U p , an interpolation basis function N k is constructed. N k = L k The following relationship is obtained:

[0019] ;

[0020] wherein, V e denotes the volume of the tetrahedral element, V k denotes the volume of the small tetrahedral element composed of point p and any three nodes in the tetrahedron excluding the k th node;

[0021] Based on this, the node potential on each unit node is calculated. , is represented as:

[0022] ;

[0023] in, U p Represents any point within a tetrahedral element. p The potential value at that location, N k Indicates the first k Interpolation basis functions for each node, U k Represents the tetrahedral element of the first k The node potential of each node, k The values ​​are 1, 2, 3, and 4.

[0024] Preferably, an inversion regularization objective function is constructed, and the model terms in the inversion regularization objective function are transformed into a constraint domain related to the gradient filtering operator and upper and lower bound constraints, including:

[0025] ;

[0026] in, The data covariance matrix is ​​a diagonal matrix related to the data error. This represents the observation data obtained at the receiving point. This represents the predicted data obtained through forward modeling. Represents the regularization factor used in the data fitting term during balanced inversion. and model items Smoothness, symbol Denotes the square of the 2-norm. The constraint domain related to the gradient filtering operator and upper and lower bound constraints is represented as:

[0027] ;

[0028] Among them, model roughness Related to gradient filtering operators, Represents model terms, , Indicates the first j Model parameters for each tetrahedral element. and They represent The lower and upper limits, This represents the initial model term.

[0029] Preferably, the formula used in gradient filtering calculation is:

[0030] ;

[0031] in, Represents the gradient and sensitivity matrix Derivative of the constrained domain model Model roughness metric , This represents the normalized residual, and .

[0032] Preferably, the model roughness metric (C) w -1 ) T Constructed using the inverse distance interpolation operator, it is represented as:

[0033] ;

[0034] ;

[0035] in, Indicates the first i The roughness of a tetrahedral element, where M represents the total number of tetrahedral elements. Indicates the relationship with the first i The number of adjacent tetrahedral elements of a tetrahedral element. p This is a parameter characterizing the smoothness of the model, with a value ranging from 0 to 1. Indicates the relationship with the first i The adjacent tetrahedral unit of the first j The distance between tetrahedral units of p The power of 0 indicates that it is related to the power of 1. i The roughness value contributed by non-adjacent tetrahedral elements to a given tetrahedral element is 0. Indicates the adjacent first j The tetrahedral unit for the first i The contribution of roughness to each tetrahedral element.

[0036] Preferably, a finite-memory quasi-Newton method is used to minimize the inversion regularization objective function, and then update the model parameters.

[0037] Compared with the prior art, the beneficial effects of the present invention include at least the following:

[0038] The method of the present application adopts finite element non-structural tetrahedron to accurately mesh the detection model of complex terrain, locally meshes the key calculation area, reduces the number of meshes while improving the calculation efficiency, adopts the total field method, and applies Dirichlet boundary condition at the medium boundary to realize high-precision forward modeling. In the inversion process, the model regularization term in the conventional objective function is changed, a new anisotropy tomography technology based on gradient filtering is used, which can reduce false anomalies in the inversion process, especially can constrain the multi-solution problem caused by multiple parameters in anisotropy inversion, and the research on underground anisotropy abnormal body shows that the method can accurately restore the anisotropy characteristics of the underground abnormal body. Compared with existing anisotropy technology, the method of the present application does not need well data, greatly reduces the data amount and data acquisition cost. BRIEF DESCRIPTION OF DRAWINGS

[0039] In order to more clearly illustrate the technical solutions of the embodiments of the present application or the prior art, the drawings needed to be used in the embodiments or the prior art description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0040] Figure 1 is a flow chart of the three-dimensional anisotropy imaging method of resistivity provided by the embodiment based on non-structural finite element;

[0041] Figure 2 is a structural schematic diagram of the non-structural tetrahedral element provided by the embodiment;

[0042] Figure 3 is a schematic diagram of anisotropy detection observation model provided by the embodiment, wherein the black dots represent receiving points, the red dots represent transmitting sources, and the square area below the transmitting source represents an abnormal body;

[0043] Figure 4 is an arbitrary anisotropy resistivity inversion result provided by the embodiment, wherein the white square box represents the true angle area, and the blue color represents the inversion result;

[0044] Figure 5 is an arbitrary anisotropy angle inversion result provided by the embodiment, wherein the white square box represents the true angle area, and the red color represents the inversion result. DETAILED DESCRIPTION

[0045] In order to make the purpose, technical solutions and advantages of the present application more clear, the present application will be further described in detail below with reference to the drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application, and do not limit the protection scope of the present application.

[0046] The inventive concept of the present application is that, in order to solve the technical problems that the current direct current resistivity anisotropy nondestructive detection research is less, and the direct current resistivity anisotropy detection effect is poor in the nondestructive detection based on only surface data, the present application provides a resistivity three-dimensional anisotropy imaging method based on unstructured finite elements, which uses an unstructured tetrahedral grid to finely discretize structures of different scales, realizes high-precision forward and adjoint forward calculation based on a node-based finite element method, converts model terms to a new constraint domain in the objective function, constructs a gradient filtering term when calculating the gradient to prevent false anomalies caused by dramatic changes in parameters, so as to achieve the purpose of constraining parameters, and updates the inversion model by using the (limited memory quasi-Newton) L-BFGS algorithm to realize three-dimensional inversion imaging of direct current anisotropy. Theoretical examples show that the method in the present application can efficiently image the three-dimensional anisotropy structure of the underground under the condition of only surface data, effectively restores the anisotropy characteristics of the geological body, and reduces false anomalies compared with the traditional method, thereby providing a reliable technical solution for three-dimensional direct current method work on complex media.

[0047] Based on the above inventive concept, as shown in Figure 1 , the embodiment provides a resistivity three-dimensional anisotropy imaging method based on unstructured finite elements, which comprises the following steps:

[0048] Step 1, a detection model for anisotropic media in a detection area and its corresponding anisotropy control equation in a direct current electric field are constructed, wherein the anisotropy control equation contains a resistivity tensor.

[0049] In the embodiment, a detection model for anisotropic media in a detection area is constructed, wherein the detection model contains a receiving point, a transmitting source, and an anomaly body. An anisotropy control equation for the anisotropic media in the detection model in a direct current electric field is also constructed, which is expressed as:

[0050] ;(1)

[0051] wherein, represents an electric potential, represents an electric current, represents a receiving point position, a source point represents a source point position, represents a distance between the receiving point position and the source point position, represents a resistivity tensor, and a symbol represents a gradient, represents a Dirac operator.

[0052] The conductivity tensor contains 6 independent components σ xy , σ yx , σ xz , σ zx, σ yz , σ zy , and is a symmetric positive definite component, i.e.

[0053] ; (2)

[0054] wherein x , y, z represent three directions of an orthogonal coordinate system, and since the conductivity tensor is positive definite and symmetric, σ xy = σ yx , σ xz = σ zx , and σ yz = σ zy . Therefore, the general electric anisotropy conductivity tensor can be calculated by six independent components.

[0055] The conductivity tensor can also be obtained by three times of Euler rotation of the principal axis conductivity tensor σ x , σ y , and σ z , i.e.

[0056] ; (3)

[0057] wherein R is a rotation matrix, R = R x R y R z , and the superscript T represents the transpose of a matrix, wherein R x , R y , and R z are three components of the rotation matrix R, and the rotation matrix R is expressed as:

[0058] ; (4)

[0059] wherein α, β, and γ respectively represent rotation angles around x, y, and z.

[0060] Step 2, the Galerkin method is used to decompose the anisotropic control equation, and a Dirichlet boundary condition is introduced to obtain a finite element integral equation with a boundary condition.

[0061] In the embodiment, for the anisotropic control equation, the Galerkin method is used to decompose the anisotropic control equation, and a Dirichlet boundary condition is introduced to obtain an equivalent finite element integral equation, which is expressed as:

[0062] ; (5)

[0063] wherein represents a detection area to be calculated, the superscript T represents the transpose, represents a potentialU The finite element integral equation function.

[0064] Step 3: The detection model is meshed and discretized using unstructured tetrahedral elements, and the node potential on each element node is calculated. The node potential is then substituted into the finite element integral equation for forward modeling.

[0065] In this embodiment, unstructured tetrahedral elements are used to mesh and discretize the detection model. The field to be determined is placed on the element nodes, and the nodal potential on each element node is calculated. , is represented as:

[0066] (6)

[0067] in, U p Represents any point within a tetrahedral element. p The potential value at that point, N k Indicates the first k The interpolation basis function of the node. U k Represents the tetrahedral element of the first generation. k The node potential of each node, k The values ​​are 1, 2, 3, and 4.

[0068] like Figure 2 The tetrahedral element shown has a potential value of p taken inside it. U p The lines connecting point p to the four nodes further divide the tetrahedron into four smaller tetrahedral elements. According to the definition of three-dimensional natural coordinates, the line connecting point p to the four nodes (excluding the first node) is defined as... k The ratio of the volume of the smaller tetrahedron formed by the other three nodes to the volume of the larger tetrahedron is a function of the volume. L k According to formula (6) in calculation U p Interpolation basis functions need to be constructed. N k When using linear interpolation N k = L k Then we get the following relationship:

[0069] (7)

[0070] in, V e This represents the volume of a tetrahedral unit cell. V k This indicates that point p is in relation to the tetrahedron except for the first... kAny three nodes outside the node form a small tetrahedron unit volume;

[0071] The formula (6) is brought into the formula (5), and after integration for each discrete unit, matrix assembly is carried out, and the following formula is obtained:

[0072] ; (8)

[0073] Wherein, K is a coefficient matrix, U is the potential to be calculated, and P is a source term. Since the coefficient matrix K is a symmetric sparse matrix, CRS (Compressed Row Storage) is used for compressed storage to reduce memory consumption. By solving the above formula (8), the potential value at any place in the calculation area is calculated, and high-precision forward simulation is realized.

[0074] S4, the inversion regularization objective function is constructed, and the model term in the inversion regularization objective function is converted into a constraint domain related to the gradient filter operator and the upper and lower limit constraint.

[0075] In the embodiment, for the inversion problem, in order to make the inversion robust, and considering that there is inherent sensitivity difference between adjacent tetrahedral elements of different volumes, when constructing the inversion regularization objective function, a gradient filtering regularization scheme is adopted, that is, the model term in the inversion regularization objective function is converted into a constraint domain related to the gradient filter operator and the upper and lower limit constraint, and the specific inversion regularization objective function is represented as:

[0076] ; (9)

[0077] Wherein, The data covariance matrix is a diagonal matrix related to the data error, The observed data obtained at the receiving point is represented as The predicted data calculated by forward simulation is represented as The regularization factor is used to balance the smoothness of the data fitting term and the model term in the inversion, the symbol represents the square of 2-norm, The constraint domain related to the gradient filter operator and the upper and lower limit constraint is represented as

[0078] ; (10)

[0079] Wherein, the model roughness is related to the gradient filter operator, The model term is represented as The initial model term is represented as , The first jmodel parameters of a tetrahedron unit, and respectively represent the lower limit and the upper limit, the model parameter m in anisotropy is written as , , represent the components of m in the resistivity of the three principal axes x, y, z and the rotation angles α, β, γ.

[0080] S5, the gradient calculation of the regularization objective function is realized by the accompanying forward simulation calculation and the gradient filter calculation.

[0081] In the embodiment, in order to calculate the minimum objective function (9), the gradient of is needed, according to the sensitivity definition and the application of the chain rule, the Jacobian matrix is obtained as follows:

[0082] ; (11)

[0083] Suppose the sensitivity matrix , the constraint domain model derivative , and the model roughness , the formula (11) can be simplified as:

[0084] ; (12)

[0085] The gradient of the formula (9) is:

[0086] ; (13)

[0087] wherein, represents the gradient, represents the normalized residual, and , . From the formula (9), it can be obtained that:

[0088] ; (14)

[0089] wherein, represents the model parameter of the i-th tetrahedron unit, represents the model term of the i-th tetrahedron unit, wherein, .

[0090] The most important thing in the gradient filter method is to find (C w -1 ) T . From the formula (13), this term directly acts on the transpose product of the sensitivity and the data residual, is used for smoothing the sensitivity and realizing the gradient filter, and is directly constructed by using the inverse distance interpolation operator, that is:

[0091] ; (15)

[0092] ; (16)

[0093] wherein, represents the roughness of the i th tetrahedron unit, M represents the total number of tetrahedron units, represents the number of tetrahedron units adjacent to the i th tetrahedron unit, p is a parameter representing the smoothness of the model, and the value range is 0-1, represents the distance i between the j th tetrahedron unit adjacent to the th tetrahedron unit, p power of the distance i , 0 represents that the roughness value contributed by other tetrahedron units not adjacent to the th tetrahedron unit is 0, j represents the contribution value of the adjacent i th tetrahedron unit to the roughness of the th tetrahedron unit.

[0094] S6, while the forward simulation calculation and the gradient filtering calculation are performed, the gradient calculation of the regularization objective function is performed. The model parameters are updated by minimizing the inversion regularization objective function, and then the resistivity three-dimensional anisotropy imaging result is obtained.

[0095] As can be seen from formula (13), in order to obtain the gradient g of the objective function, the transpose product of the sensitivity matrix J and the residual v needs to be calculated, and the sensitivity matrix J can be derived as:

[0096] ; (17)

[0097] wherein, L is an interpolation operator, and by taking the partial derivative of both sides of formula (8), the following formula is obtained:

[0098] ; (18)

[0099] wherein, , , .

[0100] Substitute formula (18) into formula (17) to obtain the transpose of the sensitivity matrix, and the product of the transpose of the sensitivity matrix and the vector J T v can be obtained by the following formula:

[0101] ; (19)

[0102] A new parameter u is introduced d :

[0103] ; (20)

[0104] After a simple transformation, we have

[0105] ; (21)

[0106] According to the definition of L T v, we only need to solve the linear equations (21) once to get u d . The limited memory quasi-Newton (L-BFGS) method is used to minimize the objective function, and the iterative optimization of the objective function is performed until the root mean square error (RMS) of the data is less than 1 or the maximum number of iterations is reached. u d is output, and the model parameters are updated according to u d , and the resistivity inversion result is obtained for imaging.

[0107] The embodiment also provides experimental examples to illustrate the effect of the above-mentioned method of the application. The model observed in the method is a two-level model, as shown in Figure 3 , the sources are uniformly distributed in the survey area, and all survey points receive the signal when each source emits the signal. An anisotropic anomaly body exists directly below the survey points, with an upper top buried at a depth of 10 m. The size of the anomaly body is 20 m*20 m*10 m, and the principal axis resistivity is ρ x / ρ y / ρ z =10 / 50 / 100Ωm, and the angles α / β / γ=10 / 20 / 30°.

[0108] Figure 4 The principal axis resistivity inversion result is shown in Figure 4 , where a, b, and c are the xoz sections of the x axis, y axis, and z axis at y=0 m. As can be seen from the figure, the three principal axis resistivities are well recovered and close to the true resistivity values, Figure 4 , where d, e, and f are the xoy sections of the x axis, y axis, and z axis at z=15 m. The resistivity values in the range of the anomaly body are well recovered. As can be seen from the above figures, the method can well recover the resistivity characteristics of the anomaly body for any anisotropy inversion.

[0109] Figure 5 The arbitrary anisotropy angle inversion result is shown in Figure 5 , where a, b, and c are the xoz sections of α, β, and γ at y=0 m, Figure 5 , where d, e, and f are the xoy sections of α, β, and γ at z=15 m. As can be known from the true model, the three angles gradually increase. The inversion results from the two direction sections can show that the method has good constraint on the angle change.

[0110] From the inversion results above, it can be concluded that the method proposed in this paper has good effect in any anisotropy inversion, not only can restore the resistivity value, but also can restore the angle information to the greatest extent. Compared with isotropic inversion, although the inversion parameters increase by 6 times, the method proposed in this paper still only needs to use the ground data to obtain the anisotropy information of the underground abnormal body, and has strong applicability.

[0111] The above specific embodiments have described the technical solutions and beneficial effects of the present application in detail. It should be understood that the above description is only the most preferred embodiment of the present application and is not intended to limit the present application. Any modification, supplement and equivalent replacement made within the principle range of the present application shall be included in the protection scope of the present application.

Claims

1. A direct current resistivity 3D anisotropy imaging method based on non-structural finite elements, characterized in that, For direct current resistivity anisotropy nondestructive detection, only under the condition of surface data, the underground three-dimensional anisotropic structure can be efficiently imaged, including the following steps: A detection model and its corresponding anisotropy control equation in the direct current field are constructed for anisotropic medium in the detection area, wherein the anisotropy control equation contains a resistivity tensor, which is expressed as: ； wherein, denotes the electric potential, denotes the electric current, denotes the receiver position, source position r e denotes the source position, denotes the distance between the receiver position and the source position, denotes the resistivity tensor, symbol denotes the gradient, denotes the Dirac operator; The anisotropy control equation is decomposed by using the Galerkin method, and a finite element integral equation with boundary conditions is obtained by introducing Dirichlet boundary conditions, which is expressed as: ； wherein represents the calculated detection region, the superscript T represents the transpose, F U represents the finite element integral equation function of the potential U ​​ A non-structural tetrahedral element is used to discretize the detection model and calculate the node potential at each element node, including: Take any point p inside the tetrahedral element; its potential value is U p The lines connecting point p to the four nodes divide the tetrahedron into four smaller tetrahedral elements. According to the definition of three-dimensional natural coordinates, the line connecting point p to the four nodes (excluding the first node) divides the tetrahedron into four smaller tetrahedral elements. k The ratio of the volume of the smaller tetrahedron formed by the other three nodes to the volume of the larger tetrahedron is a function of the volume. L k In calculation U p Interpolation basis functions need to be constructed. N k When using linear interpolation N k = L k Then we get the following relationship: ； wherein, V e represents the volume of a tetrahedron unit, V k represents the volume of a tetrahedron unit formed by the point p and any three nodes of the tetrahedron except the k first node. Based on this, the node potential on each unit node is calculated is expressed as: ； wherein, U p denotes the potential value at an arbitrary point within the tetrahedral element, p N k denotes the interpolation basis function of the k th node, U k denotes the node potential of the k th node in the tetrahedral element, k taking values 1, 2, 3, 4;​ The node potential is brought into the finite element integral equation for forward modeling; An inversion regularization objective function is constructed, and in the inversion regularization objective function, the model term is converted into a constraint domain related to the gradient filter operator and the upper and lower limit constraints, the gradient calculation of the regularization objective function is realized through the adjoint forward simulation calculation and the gradient filter calculation, the model parameters are updated by minimizing the inversion regularization objective function, and then the resistivity three-dimensional anisotropy imaging result is obtained; Wherein, the inversion regularization objective function is constructed, and in the inversion regularization objective function, the model term is converted into a constraint domain related to the gradient filter operator and the upper and lower limit constraints, including: ； where, denotes the data covariance matrix, which is a diagonal matrix related to the data error, denotes the observed data obtained at the receiving points, denotes the predicted data calculated by forward modeling, denotes the regularization factor, which is used to balance the smoothness of the data fitting term and the model term in the inversion, and the symbol denotes the square of the 2-norm, denotes the constraint domain related to the gradient filter operator and the upper and lower bound constraints, and is expressed as: ； where the model roughness is related to the gradient filter operator, is the model term, , is the model parameter of the j th tetrahedral element, and are the lower and upper bounds of , respectively, is the initial model term, and the model parameter m in anisotropy is written as , is the component of m in the x, y, z principal axis direction resistivity and the rotation angle a, b, g direction.

2. The nonstructural finite element based 3D DC resistivity anisotropy imaging method of claim 1, wherein, When the gradient filter calculation is performed, the formula used is: ； wherein, denotes the gradient, the sensitivity matrix , the constraint region model derivative , the model roughness measure , denotes the normalized residual, and .

3. The nonstructural finite element based 3D DC resistivity anisotropy imaging method of claim 1, wherein, Model roughness metric (C w -1 ) T is constructed using inverse distance interpolation operator, denoted as: ； ； in, Indicates the first i The roughness of a tetrahedral element, where M represents the total number of tetrahedral elements. Indicates the relationship with the first i The number of adjacent tetrahedral elements of a tetrahedral element. p This is a parameter characterizing the smoothness of the model, with a value ranging from 0 to 1. Indicates the relationship with the first i The adjacent tetrahedral unit of the first j The distance between tetrahedral units of p The power of 0 indicates that it is related to the power of 1. i The roughness value contributed by non-adjacent tetrahedral elements to a given tetrahedral element is 0. Indicates the adjacent first j The tetrahedral unit for the first i The contribution of roughness to each tetrahedral element.

4. The nonstructural finite element based 3D DC resistivity anisotropy imaging method of claim 1, wherein, The limited memory quasi-Newton method is used to minimize the inversion regularization objective function, and then the model parameters are updated.

Images