Resistivity three-dimensional anisotropic imaging method based on non-structural finite element

Abstract

The invention discloses a resistivity three-dimensional anisotropic imaging method based on a non-structural finite element, which belongs to the technical field of direct current method imaging in geological exploration, and comprises the following steps: constructing a detection model and an anisotropic control equation of the detection model in a direct current electric field for an anisotropic medium, converting the detection model into a finite element integral equation with boundary conditions, performing mesh generation discretization on the detection model by adopting a non-structural tetrahedron element, calculating a node potential on each element node, and substituting the node potentials into the finite element integral equation for forward modeling; according to the method, an inversion regularization objective function is constructed, a model item is converted into a constraint domain related to a gradient filtering operator and upper and lower limit constraints in the inversion regularization objective function, the inversion regularization objective function is minimized to update model parameters, and a resistivity three-dimensional anisotropic imaging result is obtained, so that the inversion accuracy is ensured, and the resistivity three-dimensional anisotropic imaging efficiency is improved. And the data volume and the data acquisition cost are greatly reduced.

Description

Technical Field

[0001] The invention belongs to the technical field of direct current electrical imaging in geological exploration, and in particular relates to a resistivity three-dimensional anisotropic imaging method based on non-structural finite elements. Background Art

[0002] At present, the field of geophysical inversion generally regards underground media as isotropic. This assumption has been proven to be applicable through long-term experiments. However, more and more studies have shown that underground media are not simply isotropic. The results of inversion using isotropic theory have many false anomalies, which make people get a wrong geological understanding. Anisotropic characteristics are widely present in nature, especially in areas with developed stratification, cracks, and cavities rich in groundwater. These places are often rich in resources or have potential hazards, and are areas that need to be detected. The study of anisotropy began decades ago and has developed rapidly in the electromagnetic field. It has been widely used in resource exploration, potential geological hazards, etc.

[0003] The DC resistivity method is a commonly used geophysical exploration method in the electromagnetic field. Due to the inherent characteristics of the DC resistivity method, the detection range is mostly near the surface area. Incorporating the concept of anisotropy into the inversion imaging will undoubtedly increase the accuracy of imaging and be more practical in resource exploration and engineering applications. The study of DC anisotropy has a long history, but due to the increase in inversion parameters, the non-uniqueness of the inversion results has become more obvious, which has brought great difficulties and challenges to the inversion work. In the early days, people used simple models such as half-space or one-dimensional anisotropic layered media to study the anisotropy law, and gave analytical solutions to related models, which provided a reliable basis for subsequent DC anisotropy research. As the number of medium layers increases, the model also tends to be complex. Although one-dimensional anisotropic inversion work is performed, the results are no longer reliable due to the inherent non-uniqueness of the inversion and the increase in parameters.

[0004] As the dimension increases and the geological structure becomes more complex, people cannot directly obtain the corresponding solution and have to use numerical simulation methods for approximate solutions. In the study of two-dimensional anisotropy, finite differences are used for forward modeling, and regularization factors are added during the inversion process to reduce the artifacts of the inversion results and greatly improve the calculation speed.

[0005] Since electrical anisotropy is directional, conventional two-dimensional data acquisition methods have great limitations and cannot obtain effective anisotropy information. Therefore, it is necessary to consider the use of three-dimensional observation devices to bring higher quality data results and conduct three-dimensional anisotropy research on DC resistivity. In order to obtain better data information, well data is particularly important. Well data can better invert the current flow direction and reconstruct the geological structure. The use of unstructured grids can finely divide complex terrain. A large number of researchers have conducted a lot of theoretical research on well, well-ground and well-well data, analyzing the impact of different data combinations on anisotropic inversion results, and have been well applied in actual engineering, and the results obtained are better than isotropy.

[0006] Currently, non-destructive detection has become an important research hotspot. When inter-well measurements cannot be carried out in certain specific areas, non-destructive detection methods become more important. This method does not require the destruction of geological structures and can save a lot of costs. However, it is extremely difficult to invert anisotropic multi-parameters using only surface data. Summary of the invention

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

[0008] To achieve the above-mentioned object of the invention, an embodiment provides a resistivity three-dimensional anisotropy imaging method based on non-structural finite elements, comprising the following steps: A detection model and its corresponding anisotropic control equation in a DC electric field are constructed for the anisotropic medium in the detection area, wherein the anisotropic control equation contains the resistivity tensor. The anisotropic control equation is decomposed by Galerkin, and the Dirichlet boundary condition is introduced to obtain the finite element integral equation with boundary conditions. The detection model is meshed and discretized using unstructured tetrahedral units, and the node potential on each unit node is calculated. The node potential is brought into the finite element integral equation for forward simulation. An inversion regularization objective function is constructed, and the model terms are converted into constraint domains related to the 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 by accompanying forward simulation calculation and gradient filtering calculation. The inversion regularization objective function is minimized to update the model parameters, and then the three-dimensional anisotropic resistivity imaging results are obtained.

[0009] Preferably, the anisotropy governing equation is expressed as: ; in, represents the potential, represents the current, Indicates the location of the receiving point, the source point represents the source location, Represents the distance between the receiving point and the source point. Represents the resistivity tensor, symbol represents the gradient, represents the Dirac operator.

[0010] Preferably, the Galerkin decomposition of the anisotropic control equation is adopted, and the Dirichlet boundary condition is introduced to obtain the finite element integral equation with boundary conditions, which is expressed as: ; in, Indicates the calculated detection area, the superscript T represents transpose, Indicates potential U Finite element integral equation function.

[0011] Preferably, the detection model is meshed and discretized using unstructured tetrahedral units and the node potential on each unit node is calculated, including: Take any point p inside the tetrahedral unit, and its potential value is U p The lines connecting point p and the four nodes divide the tetrahedron into four small tetrahedral units. According to the three-dimensional natural coordinate definition, the lines connecting point p and the four nodes without k The ratio of the volume of the small tetrahedron formed by the other three nodes of the node to the volume of the large tetrahedron is the volume function L k , in the calculation U p Interpolation basis functions need to be constructed N k , when linear interpolation is used N k = L k , we get the following relationship: ; in, V e represents the volume of tetrahedral unit, V k represents the distance between point p and the tetrahedron except the k The small tetrahedral unit volume composed of any three nodes outside the nodes; Based on this, the node potential on each unit node is calculated , expressed as: ; in, U p Represents any point in the tetrahedral unitp The potential value at N k Indicates k The interpolation basis functions of nodes, U k Indicates the tetrahedral unit k The node potential of each node, k The value is 1, 2, 3, 4.

[0012] Preferably, an inversion regularization objective function is constructed, and the model term is converted into a constraint domain related to the gradient filter operator and the upper and lower limit constraints in the inversion regularization objective function, including: ; in, represents the data covariance matrix, which is a diagonal matrix associated with the data error, represents the observation data obtained at the receiving point, represents the predicted data obtained through forward simulation calculation, Represents the regularization factor, which is used to balance the data fitting term in the inversion and model terms The smoothness of represents the square of the 2-norm, Represents the constraint domain related to the gradient filter operator and the upper and lower limit constraints, expressed as: ; Among them, the model roughness Related to the gradient filter operator, represents the model term, , Indicates j The model parameters of the tetrahedral elements are: and Respectively The lower and upper limits of represents the initial model terms.

[0013] Preferably, the formula used in gradient filtering calculation is: ; in, represents the gradient, sensitivity matrix , the constraint domain model derivative , model roughness metric , represents the normalized residual, and .

[0014] Preferably, the model roughness metric (C w -1 ) TIt is constructed using the inverse distance interpolation operator, expressed as: ; ; in, Indicates i The roughness of tetrahedral units, M represents the total number of tetrahedral units, Indicates i The number of tetrahedral units adjacent to a tetrahedral unit, p It is a parameter that characterizes the smoothness of the model, and its value range is 0-1. Indicates i The adjacent tetrahedral unit j The distance between tetrahedral elements of p Power, 0 means the same as the i The roughness value contributed by other tetrahedral elements that are not adjacent to the tetrahedral element is 0. Indicates the adjacent j The tetrahedral unit is i The roughness contribution of the tetrahedral elements.

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

[0016] Compared with the prior art, the present invention has the following beneficial effects: The method of the present invention uses finite element unstructured tetrahedron to accurately mesh the detection model of complex terrain, and locally encrypts the grid in the key calculation area, thereby reducing the number of grids and improving the calculation efficiency. The method of the present invention adopts the total field method and applies the Dirichlet boundary condition at the medium boundary to achieve high-precision forward simulation. The inversion process changes the model regularization term in the conventional objective function and uses a new direct current anisotropic tomography technology based on gradient filtering. This method can reduce false anomalies in the inversion process, especially can constrain the multi-solution problem caused by multiple parameters in anisotropic inversion. The study of underground anisotropic anomalies shows that this method can accurately restore the anisotropic characteristics of underground anomalies. Compared with the existing anisotropic technology, the method of the present invention does not require well data, which greatly reduces the amount of data and data acquisition costs. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0018] Figure 1 It is a flowchart of a resistivity three-dimensional anisotropy imaging method based on non-structured finite elements provided in an embodiment; Figure 2 is a schematic structural diagram of a non-structural tetrahedron unit provided in an embodiment; Figure 3 is a schematic diagram of an anisotropic detection observation model provided in an embodiment, wherein a black dot represents a receiving point, a red dot represents a transmitting source, and a square area below the transmitting source represents an abnormal body; Figure 4 is an arbitrary anisotropic resistivity inversion result provided by the embodiment, wherein the white box represents the true angle area, and the blue represents the inversion result; Figure 5 : is an arbitrary anisotropic angle inversion result provided by the embodiment, wherein the white box represents the true angle area, and the red box represents the inversion result. DETAILED DESCRIPTION

[0019] To make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific implementation methods described herein are only used to explain the present invention and do not limit the scope of protection of the present invention.

[0020] The inventive concept of the present invention is: in view of the technical problems that there are few studies on nondestructive detection of DC resistivity anisotropy and the detection effect of DC resistivity anisotropy is not good in nondestructive detection based on surface data only, the present invention provides a resistivity three-dimensional anisotropic imaging method based on non-structured finite elements, which uses unstructured tetrahedral grids to finely discretize structures of different scales, realizes high-precision forward modeling and adjoint forward modeling calculations based on node-based finite element method, transforms model terms into new constraint domains in the objective function, constructs gradient filter terms when gradients are obtained, prevents drastic changes in parameters from causing false anomalies, thereby achieving the purpose of constraining parameters, and uses (limited memory quasi-Newton) L-BFGS algorithm to update the inversion model to realize DC anisotropic three-dimensional inversion imaging. Theoretical examples show that the method can efficiently image underground three-dimensional anisotropic structures with only surface data, effectively restores the anisotropic characteristics of geological bodies, reduces false anomalies compared with traditional methods, and provides a reliable technical solution for carrying out three-dimensional DC electrical method work on complex media.

[0021] Based on the above invention concept, Figure 1 As shown, an embodiment provides a resistivity three-dimensional anisotropy imaging method based on a non-structured finite element, comprising the following steps: Step 1: construct a detection model for the anisotropic medium in the detection area and its corresponding anisotropic control equation in a DC electric field, wherein the anisotropic control equation contains a resistivity tensor.

[0022] In the embodiment, a detection model is constructed for the anisotropic medium in the detection area, wherein the detection model includes a receiving point, a transmitting source, and an abnormal body, etc. An anisotropic control equation in a DC electric field is also constructed for the anisotropic medium in the detection model, which is expressed as: ; (1) in, represents the potential, represents the current, Indicates the location of the receiving point, the source point represents the source location, Represents the distance between the receiving point and the source point. Represents the resistivity tensor, symbol represents the gradient, represents the Dirac operator.

[0023] Conductivity tensor Contains 6 independent components σ xy , σ yx , σ xz , σ zx , σ yz , σ zy ,and is a symmetric positive definite component, that is: ; (2) in x , y,z Represent the three directions of the rectangular coordinate system. Since the conductivity tensor is positive definite symmetric, σ xy =σ yx ,σ xz =σ zx , σ yz =σ zy Therefore, the conductivity tensor for general electrical anisotropy can be calculated from six independent components.

[0024] The conductivity tensor can also be expressed by the principal axis conductivity tensor σ x , σ y , σ z After 3 Euler rotations, we get: ; (3) Where R is the rotation matrix, R=R x R y R z , the superscript T represents the transpose of the matrix, where Rx , R y , R z are the three components of the rotation matrix R, which is expressed as: ; (4) Among them, α, β, and γ represent the rotation angles around x, y, and z respectively.

[0025] Step 2: Use the Galerkin method to decompose the anisotropic control equation and introduce the Dirichlet boundary condition to obtain the finite element integral equation with boundary conditions.

[0026] In the embodiment, the Galerkin method is used to decompose the anisotropic control equation, and the Dirichlet boundary condition is introduced to obtain an equivalent finite element integral equation, which is expressed as: ; (5) in, Indicates the calculated detection area, the superscript T represents transpose, Indicates potential U The finite element integral equation function.

[0027] Step 3: Use unstructured tetrahedral units to mesh and discretize the detection model and calculate the node potential on each unit node, and bring the node potential into the finite element integral equation for forward simulation.

[0028] In the embodiment, the detection model is meshed and discretized using unstructured tetrahedral units, the field to be calculated is placed on the unit nodes, and the node potential on each unit node is calculated. , expressed as: ; (6) in, U p Represents any point in the tetrahedral unit p The potential value at N k Indicates k The interpolation basis functions of the nodes, U k Indicates the tetrahedral unit k The node potential of each node, k The value is 1, 2, 3, 4.

[0029] like Figure 2 The tetrahedral unit shown in the figure has a potential value of U p The lines connecting point p and the four nodes divide the tetrahedron into four small tetrahedral units. According to the three-dimensional natural coordinate definition, the lines connecting point p and the four nodes withoutk The ratio of the volume of the small tetrahedron formed by the other three nodes of the node to the volume of the large tetrahedron is the volume function L k , according to formula (6) in the calculation U p Interpolation basis functions need to be constructed N k , when linear interpolation is used N k = L k , we get the following relationship: ; (7) in, V e represents the volume of tetrahedral unit, V k represents the distance between point p and the tetrahedron except the k The small tetrahedral unit volume composed of any three nodes outside the nodes; Substituting formula (6) into formula (5) and performing matrix assembly after integrating each discrete unit, we can obtain: ; (8) Among them, K is the coefficient matrix, U is the potential to be calculated, and P is the source term. Since the coefficient matrix K is a symmetric sparse matrix, CRS (Compressed Row Storage) is used to compress the storage to reduce memory consumption. By solving the above equation (8), the potential value at any point in the calculation area is solved to achieve high-precision forward simulation.

[0030] S4, construct an inversion regularization objective function, and transform the model terms in the inversion regularization objective function into constraint domains related to the gradient filter operator and upper and lower limit constraints.

[0031] In the embodiment, for the inversion problem, in order to make the inversion robust and considering the inherent sensitivity difference between adjacent tetrahedral elements of different volumes, a gradient filtering regularization scheme is adopted when constructing the inversion regularization objective function, that is, the model term in the inversion regularization objective function is converted into a constraint domain related to the gradient filtering operator and the upper and lower limit constraints. The specific inversion regularization objective function is expressed as: ; (9) in, represents the data covariance matrix, which is a diagonal matrix associated with the data error, represents the observation data obtained at the receiving point, represents the predicted data obtained through forward simulation calculation, Represents the regularization factor, which is used to balance the data fitting term in the inversion and model terms The smoothness of represents the square of the 2-norm, Represents the constraint domain related to the gradient filter operator and the upper and lower limit constraints, expressed as: ; (10) Among them, the model roughness Related to the gradient filter operator, represents the model term, represents the initial model term, , Indicates j The model parameters of the tetrahedral elements are: and Respectively The lower and upper limits of the anisotropic model parameter m are written as , It represents the resistivity of m in the three principal axis directions of x, y, and z, as well as the components of the rotation angles α, β, and γ.

[0032] S5, regularized objective function gradient calculation is realized by adjoint forward simulation calculation and gradient filtering calculation.

[0033] In the embodiment, in order to calculate the minimization objective function (9), it is necessary to calculate The gradient of , according to the sensitivity definition and applying the chain rule, can be obtained as for: ; (11) Assuming the sensitivity matrix , the constraint domain model derivative , model roughness metric , formula (11) can be simplified to obtain: ; (12) The gradient of formula (9) is: ; (13) in, represents the gradient, represents the normalized residual, and , From formula (9), we can get: ; (14) in, represents the model parameters of the i-th tetrahedral unit, represents the model term of the i-th tetrahedral unit, where .

[0034] The most important thing in the gradient filtering method is to find (C w -1 ) T From formula (13), we can see that this term directly acts on the transposed product of sensitivity and data residual, which is used to smooth the sensitivity and realize gradient filtering. It is directly constructed using the inverse distance interpolation operator, that is: ; (15) ; (16) in, Indicates i The roughness of tetrahedral units, M represents the total number of tetrahedral units, Indicates i The number of tetrahedral units adjacent to a tetrahedral unit, p It is a parameter that characterizes the smoothness of the model, and its value range is 0-1. Indicates i The adjacent tetrahedral unit j The distance between tetrahedral elements of p Power, 0 means the same as the i The roughness value contributed by other tetrahedral elements that are not adjacent to the tetrahedral element is 0. Indicates the adjacent j The tetrahedral unit is i The roughness contribution of the tetrahedral elements.

[0035] S6, while performing forward simulation calculation and gradient filtering calculation, the gradient of the regularized objective function is calculated. The inversion regularized objective function is minimized to update the model parameters, and then the resistivity three-dimensional anisotropic imaging results are obtained.

[0036] It can be seen from formula (13) that in order to obtain the gradient g of the objective function, it is necessary to calculate the transposed product of the sensitivity matrix J and the residual v. The sensitivity matrix J can be derived as: ; (17) Where L is the interpolation operator. By taking partial derivatives on both sides of formula (8), we get: ; (18) in, , , .

[0037] Substituting formula (18) into formula (17) yields the transpose of the sensitivity matrix. The product of the transpose of the sensitivity matrix and the vector J is T v can be obtained by the following formula: ; (19) Introduce a new parameter u d : ; (20) After a simple transformation, we can get: ;(twenty one) According to L T The definition of v is that we only need to solve the linear equation system (21) once to get u d The objective function is minimized using the limited memory quasi-Newton (L-BFGS) method, and the objective function is iteratively optimized until the root mean square error (RMS) of the data is less than 1 or the maximum number of iterations is reached. Output u d , according to u d The model parameters are updated to obtain resistivity inversion results for imaging.

[0038] The embodiment also provides an experimental example to illustrate the effect of the above method of the present invention. The observation model of the method in this paper adopts a secondary model, such as Figure 3 As shown, the sources are evenly distributed in the measurement area. When each source transmits a signal, all measurement points receive it. An anisotropic anomaly exists directly below the measurement point. The top surface is buried 10m deep. The size of the anomaly is 20m*20m*10m. The main axis resistivity is ρ x / ρ y / ρ z =10 / 50 / 100Ωm, angle α / β / γ=10 / 20 / 30°.

[0039] Figure 4 is the main axis resistivity inversion result, Figure 4 a, b and c are the xoz sections of the x-axis, y-axis and z-axis at y=0m. It can be seen from the figure that the resistivity of the three main axes is well restored and close to the true resistivity value. Figure 4 In the figure, d, e and f are the xoy sections of the x-axis, y-axis and z-axis at z=15m respectively. The resistivity of the abnormal body is well restored. From the above figures, it can be seen that the method in this paper can well restore the resistivity characteristics of the abnormal body for any anisotropic inversion.

[0040] Figure 5 is the inversion result of arbitrary anisotropy angle, Figure 5 a, b and c are the xoz sections of α, β and γ at y=0m, respectively. Figure 5 In the figure, d, e and f are the xoy sections of α, β and γ at z=15m respectively. From the real model, we can know that the three angles gradually increase. From the inversion results of the slices in two directions, we can see that the method in this paper has good constraints on the angle changes.

[0041] From the above inversion results, it can be concluded that the method proposed in this paper has a good effect in any anisotropic inversion. It can not only restore the resistivity value, but also restore the angle information to the greatest extent. Compared with the isotropic inversion, although the inversion parameters are increased by 6 times, the method proposed in this paper still only needs to use ground data to obtain the anisotropic information of underground anomalies, which has strong applicability.

[0042] The specific implementation methods described above provide a detailed description of the technical solutions and beneficial effects of the present invention. It should be understood that the above is only the most preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, supplements and equivalent substitutions made within the scope of the principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A resistivity three-dimensional anisotropy imaging method based on non-structured finite elements, characterized in that: The following steps are involved: A detection model and its corresponding anisotropic control equation in a DC electric field are constructed for the anisotropic medium in the detection area, wherein the anisotropic control equation contains the resistivity tensor. The anisotropic control equation is decomposed by the Galerkin method, and the Dirichlet boundary condition is introduced to obtain the finite element integral equation with boundary conditions. The detection model is meshed and discretized using unstructured tetrahedral units, and the node potential on each unit node is calculated. The node potential is brought into the finite element integral equation for forward simulation. An inversion regularization objective function is constructed, and the model terms are converted into constraint domains related to the 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 by accompanying forward simulation calculation and gradient filtering calculation. The inversion regularization objective function is minimized to update the model parameters, and then the three-dimensional anisotropic resistivity imaging results are obtained.

2. The resistivity three-dimensional anisotropy imaging method based on non-structured finite elements according to claim 1, characterized in that: The anisotropic governing equation is expressed as: ； in, represents the potential, represents the current, Indicates the location of the receiving point, the source point represents the source location, Represents the distance between the receiving point and the source point. Represents the resistivity tensor, symbol represents the gradient, represents the Dirac operator.

3. The resistivity three-dimensional anisotropy imaging method based on non-structured finite elements according to claim 2, characterized in that: The Galerkin decomposition of the anisotropic control equation and the introduction of the Dirichlet boundary condition are used to obtain the finite element integral equation with boundary conditions, which can be expressed as: ； in, Indicates the calculated detection area, the superscript T represents transpose, Indicates potential U The finite element integral equation function.

4. The resistivity three-dimensional anisotropy imaging method based on non-structured finite elements according to claim 2, characterized in that: The detection model is meshed and discretized using unstructured tetrahedral elements and the node potential on each element node is calculated, including: Take any point p inside the tetrahedral unit, and its potential value is U p The lines connecting point p and the four nodes divide the tetrahedron into four small tetrahedral units. According to the three-dimensional natural coordinate definition, the lines connecting point p and the four nodes without k The ratio of the volume of the small tetrahedron formed by the other three nodes of the node to the volume of the large tetrahedron is the volume function L k , in the calculation U p Interpolation basis functions need to be constructed N k , when linear interpolation is used N k = L k , we get the following relationship: ； in, V e represents the volume of tetrahedral unit, V k represents the distance between point p and the tetrahedron except the k The small tetrahedral unit volume composed of any three nodes outside the nodes; Based on this, the node potential on each unit node is calculated , expressed as: ； in, U p Represents any point in the tetrahedral unit p The potential value at N k Indicates k The interpolation basis functions of nodes, U k Indicates the tetrahedral unit k The node potential of each node, k The value is 1, 2, 3, 4.

5. The resistivity three-dimensional anisotropy imaging method based on non-structured finite elements according to claim 1, characterized in that: Construct an inversion regularization objective function and transform the model terms in the inversion regularization objective function into constraint domains related to the gradient filter operator and upper and lower limit constraints, including: ； in, represents the data covariance matrix, which is a diagonal matrix associated with the data error, represents the observation data obtained at the receiving point, represents the predicted data obtained through forward simulation calculation, Represents the regularization factor, which is used to balance the data fitting term in the inversion and model terms The smoothness of represents the square of the 2-norm, Represents the constraint domain related to the gradient filter operator and the upper and lower limit constraints, expressed as: ； Among them, the model roughness Related to the gradient filter operator, represents the model term, , Indicates j The model parameters of the tetrahedral elements are: and Respectively The lower and upper limits of represents the initial model terms.

6. The resistivity three-dimensional anisotropy imaging method based on non-structural finite elements according to claim 5, characterized in that: When calculating the gradient filter, the formula used is: ； in, represents the gradient, sensitivity matrix , the constraint domain model derivative , model roughness metric , represents the normalized residual, and .

7. The resistivity three-dimensional anisotropy imaging method based on non-structured finite elements according to claim 6, characterized in that: Model roughness metric (C w -1 ) T It is constructed using the inverse distance interpolation operator, expressed as: ； ； in, Indicates i The roughness of tetrahedral units, M represents the total number of tetrahedral units, Indicates i The number of tetrahedral units adjacent to a tetrahedral unit, p It is a parameter that characterizes the smoothness of the model, and its value range is 0-1. Indicates i The adjacent tetrahedral unit j The distance between tetrahedral elements of p Power, 0 means the same as the i The roughness value contributed by other tetrahedral elements that are not adjacent to the tetrahedral element is 0. Indicates the adjacent j The tetrahedral unit is i The roughness contribution of the tetrahedral elements.

8. The resistivity three-dimensional anisotropy imaging method based on non-structured finite elements according to claim 1, characterized in that: A limited memory quasi-Newton method is used to minimize the inverse regularization objective function and then update the model parameters.

Images