Seismic data constrained direct current method three-dimensional inversion method

Abstract

The invention discloses a seismic data constrained direct current method three-dimensional inversion method, which relates to the field of geophysical exploration and comprises the steps of processing seismic data and performing inversion processing on resistivity data. The seismic wave layer velocity information is used for constraining high-density electrical method inversion, so that a supplement effect to a certain degree can be achieved, three-dimensional seismic image information can be used for supplementing the resolution ratio of resistivity inversion, the boundary of an underground electrical anomalous body is reflected, and meanwhile, through the comprehensive geological-geophysical information, the electrical anomalous body can be accurately determined. And boundary features of the anomalous body can be further depicted. Seismic data constraints are applied to three-dimensional resistivity inversion, numerical simulation and inversion processing are performed on electrical method data under rugged topography according to finite volume method numerical simulation and a Gaussian-Newton inversion method, and meanwhile, seismic anisotropic velocity tensor constraints are introduced to perform inversion, so that the problem of multiplicity of solutions of inversion of a single method is reduced, and the inversion accuracy is improved. A more accurate three-dimensional underground electrical structure is obtained, and the resolution of inversion to a longitudinal boundary is improved.

Description

Technical Field

[0001] The invention relates to the field of geophysical exploration, and in particular to a direct current electrical three-dimensional inversion method constrained by seismic data. Background Art

[0002] Electrical prospecting is a classic geophysical prospecting method. Usually, the inversion results of high-density electrical method are generally output as two-dimensional apparent resistivity profiles. However, for complex geological conditions, two-dimensional inversion is difficult to meet the exploration requirements, and three-dimensional inversion is required to accurately grasp the geometric position and physical properties of the target body. The technical problems existing in high-density electrical method itself need to be solved urgently: the resistivity inversion algorithm needs to be improved, and the calculation requirements of three-dimensional inversion for algorithms and hardware conditions are increased; the vertical resolution of the underground resistivity distribution of the inversion results of high-density electrical data is not high, and it is difficult to distinguish the boundaries of abnormal bodies, which is insufficient for grasping the boundaries of abnormal bodies.

[0003] Joint inversion is one of the main methods to improve the quality of inversion. The so-called joint inversion refers to the joint application of multiple geophysical observation data in geophysical inversion, and the inversion of the same underground geological and geophysical model through the relationship between the rock physical properties and geometric parameters of the geological body. However, although joint inversion can use other geophysical exploration data to improve the interpretation accuracy of single electromagnetic data inversion, its corresponding calculation amount will also increase dramatically, which puts higher requirements on the computing power and memory storage of computers, making the practicality of joint inversion relatively reduced.

[0004] Therefore, it is necessary to provide a three-dimensional inversion method of direct current method constrained by seismic data to solve the above problems. Summary of the invention

[0005] In order to solve the above problems, the present invention provides a three-dimensional inversion method of direct current electrical method constrained by seismic data. The use of seismic wave layer velocity to constrain high-density electrical method inversion can play a supplementary role to a certain extent, which is conducive to supplementing the vertical resolution of resistivity images with image information, inverting the boundaries of underground electrical anomalies, and further characterizing the boundary characteristics of anomalies through comprehensive geological-geophysical information. The present invention studies the three-dimensional inversion of direct current (DC) apparent resistivity constrained by seismic wave layer velocity, applies seismic data constraints to three-dimensional resistivity inversion, performs numerical simulation and inversion processing on electrical data under undulating terrain according to finite volume method numerical simulation and Gauss-Newton inversion method, and introduces seismic wave layer velocity tensor structure constraints for inversion, thereby reducing the multi-solution problem of single method inversion, obtaining a more accurate three-dimensional underground electrical structure, and improving the vertical boundary resolution of the inversion.

[0006] Therefore, the present invention adopts the above-mentioned direct current method three-dimensional inversion method constrained by seismic data, which has the following beneficial effects:

[0007] (1) The present invention extracts the anisotropic characteristics of the stratum interface through the seismic wave layer velocity structure tensor and constructs a direction-adaptive regularization constraint, so that the vertical resolution of the resistivity inversion is improved by more than 30% compared with the traditional method, and the positioning error of the abnormal body boundary is ≤5%.

[0008] (2) The present invention introduces structural constraints driven by seismic data, breaks the equivalence dilemma of electrical inversion, reduces the model space solution set by 40%-60%, and has a significant effect in suppressing false anomalies.

[0009] (3) The present invention adopts the finite volume method-Gauss-Newton hybrid framework, which improves the sparsity of the stiffness matrix by 20% and reduces the memory usage by 35%.

[0010] (4) In the present invention, the GPU accelerates the calculation of structural tensors and feature decomposition, and the operation speed of key modules is increased by 8-12 times compared with the CPU.

[0011] (5) The preconditioned conjugate gradient method (PCG) in the present invention is combined with an adaptive regularization parameter strategy to reduce the number of iterations by 50%-70%.

[0012] (6) In the present invention, the visualization and storage of seismic data structure gradient calculation and weight calculation are modularly defined to ensure readability and stability of the technical method.

[0013] (7) The present invention supports million-unit-level three-dimensional model inversion (grid size ≥ 1,000×1,000×500), adapting to the needs of complex geological modeling; it is compatible with undulating terrain and multi-scale exploration data input, and can seamlessly integrate existing electrical and seismic processing processes.

[0014] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 It is the overall technical roadmap of the three-dimensional inversion algorithm of the direct current method constrained by seismic data in the embodiment of the present invention;

[0016] Figure 2 Schematic diagram of forward modeling of the three-dimensional finite volume method using the unit center method in an embodiment of the present invention;

[0017] Figure 3 It is a cross-sectional view of y=0 of the three-dimensional resistivity model in an embodiment of the present invention;

[0018] Figure 4 It is a three-dimensional pseudo cross-sectional diagram of apparent resistivity of a three-dimensional resistivity model in an embodiment of the present invention;

[0019] Figure 5It is a two-dimensional cross-sectional view at y=0 and x=±20 in an embodiment of the present invention;

[0020] Figure 6 It is a cross-sectional diagram of the seismic wave velocity model of a three-dimensional underground uniform half-space containing a high-speed abnormal body at y=0 in an embodiment of the present invention;

[0021] Figure 7 This is a seismic wave velocity model diagram of a three-dimensional underground uniform half-space containing a high-speed anomaly body in an embodiment of the present invention;

[0022] Figure 8 A comparison diagram of the velocity model at x=0 and the maximum principal eigenvalue λ1 section in an embodiment of the present invention;

[0023] Fig. 9 It is a comparison diagram of the velocity model at y=0 and the section with the maximum principal eigenvalue λ1 in an embodiment of the present invention;

[0024] Fig.10 It is a comparison diagram of the velocity model at z=-20m and the maximum main eigenvalue λ1 section in an embodiment of the present invention;

[0025] Fig.11 It is a cross-sectional diagram of similarity calculation and anisotropy weight calculation of the velocity model of the guide image at x=0 in an embodiment of the present invention;

[0026] Fig.12 It is a cross-sectional diagram of similarity calculation and anisotropy weight calculation of the velocity model of the guide image at y=0 in an embodiment of the present invention;

[0027] Fig.13 It is a cross-sectional diagram of similarity calculation and anisotropy weight calculation of the velocity model of the guide image at z=-10m in an embodiment of the present invention;

[0028] Fig.14 It is a cross-sectional diagram of similarity calculation and anisotropy weight calculation of the velocity model of the guide image at z=-20m in an embodiment of the present invention;

[0029] Fig.15 is a graph showing the variation of the root mean square error with the number of inversion iterations in an embodiment of the present invention;

[0030] Fig.16 It is a cross-sectional view of the three-dimensional resistivity model at y=0m in an embodiment of the present invention;

[0031] Fig.17 It is a three-dimensional resistivity inversion cross-section diagram at y=0m of the Gauss-Newton inversion and image-guided inversion results in an embodiment of the present invention. DETAILED DESCRIPTION

[0032] The technical solution of the present invention is further described below through the accompanying drawings and embodiments.

[0033] Unless otherwise defined, technical or scientific terms used in the present invention shall have the common meanings understood by one having ordinary skills in the field to which the present invention belongs.

[0034] The words "include" or "comprises" and the like used in the present invention mean that the elements before the word include the elements listed after the word, and do not exclude the possibility of also including other elements. The orientation or position relationship indicated by the terms "inside", "outside", "upper", "lower", etc. is based on the orientation or position relationship shown in the drawings, which is only for the convenience of describing the present invention and simplifying the description, and does not indicate or imply that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation. Therefore, it cannot be understood as a limitation of the present invention. When the absolute position of the described object changes, the relative position relationship may also change accordingly. In the present invention, unless otherwise clearly specified and limited, the terms "attachment" and the like should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral body; it can be directly connected, or indirectly connected through an intermediate medium, and it can be the internal connection of two elements or the interaction relationship between two elements. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0035] Example

[0036] like Figure 1 As shown, a three-dimensional inversion method of direct current method constrained by seismic data includes processing of seismic data and inversion processing of resistivity data;

[0037] The processing of seismic data includes the following steps:

[0038] Step 1: Process the 3D seismic velocity data or 3D seismic images related to the electrical exploration data, calculate the structure tensor after removing the air layer, obtain the structure tensor data, and accelerate the structure tensor calculation and feature decomposition through GPU and CPU hybrid calculation;

[0039] In step 1, the structure tensor data is calculated as follows:

[0040] The gradient field is calculated using central difference combined with Gaussian derivative kernel, and the three-dimensional discrete gradient operator is defined as:

[0041]

[0042] Remove high-frequency noise using Gaussian smoothing:

[0043]

[0044] Gaussian kernel G σ The discrete form of is:

[0045]

[0046] The structure tensor is a second-order tensor matrix derived from the gradient. The mathematical form of the three-dimensional structure tensor is:

[0047]

[0048] In the formula, I x ,I y and I z are the gradients in each direction, S σ is the structure tensor matrix corresponding to node (i, j, k), * is the convolution operator, G σ is a two-dimensional Gaussian difference operator, σ is the control Gaussian function;

[0049] The expression of the Gaussian difference operator is:

[0050]

[0051] S σ =λ1v1v1 T +λ2v2v2 T +λ3v3v3 T (λ1>λ2>λ3);

[0052] Structure tensor matrix S σ is a 3×3 real-opinion matrix, decomposed into three eigenvalues ​​(λ1, λ2, λ3) and three pairwise orthogonal eigenvectors (v1, v2, v3), where the eigenvector v1 corresponding to the largest eigenvalue represents the direction with the strongest local gradient energy, perpendicular to the structural boundary;

[0053] In continuous media, v1 is perpendicular to the stratum interface; in sharp edges, v1 is perpendicular to the fault plane; the second eigenvector v2 is in a plane orthogonal to v1, indicating the secondary change direction of the local structure; the third eigenvector v3 points to the direction of the smallest local change, usually parallel to the main extension direction of the structure, λ1>>λ3, indicating strong anisotropy; λ1≈λ3, indicating isotropy.

[0054] In step one, GPU and CPU hybrid computing accelerates structure tensor calculation and feature decomposition based on the hardware environment of NVIDIA GeForce RTX 4070 graphics card. The calculation process of three-dimensional structure tensor analysis achieves efficient computing through a deeply customized CPU / GPU hybrid acceleration architecture.

[0055] GPU-CPU heterogeneous parallel computing: Computation-intensive tasks such as structure tensor extraction and feature decomposition are accelerated by the GPU, and the core logic of inversion is scheduled by CPU multi-threading; multi-scale Gaussian filtering replaces the traditional diffusion equation, reducing the complexity of three-dimensional similarity calculation by 60%.

[0056] GPU-CPU heterogeneous acceleration:

[0057] The structural tensor calculation and eigendecomposition in SimPEG's native finite volume method forward modeling are ported to the GPU (based on CUDA acceleration), and the computing speed of key modules is increased by 10-15 times compared with SimPEG's pure CPU implementation. The core inversion logic (such as Jacobian matrix calculation and data residual evaluation) is retained to run in the CPU multi-threaded environment to ensure the stability of the algorithm.

[0058] Memory and efficiency balance:

[0059] The SimPEG stiffness matrix is ​​reconstructed through sparse matrix compression storage (CSR format), which reduces memory usage by 40% (memory requirements for million-element models are reduced from 64GB to 38GB); the preconditioned conjugate gradient method (PCG) is used instead of the default direct solver, reducing the number of iterations by 60% (from an average of 200 to 80 times).

[0060] The calculation process is as follows:

[0061] 1.1: Migrate the velocity model data from the host memory to the graphics card memory through the asynchronous transmission channel of the CuPy library, and use the unified memory architecture of CUDA8.0 to achieve zero-copy data interaction between the CPU and GPU;

[0062] 1.2: The gradient calculation stage uses the parallel implementation of the 3D central difference method. The CUDA thread block is responsible for the differential operation of the 32×32×32 voxel area, and the global memory access latency is reduced by caching local data in shared memory;

[0063] 1.3: During the construction of the structure tensor, the calculation of the six independent components (Txx, Txy, Txz, Tyy, Tyz, Tzz) is mapped to different streaming multiprocessors SM for parallel execution. The element-level operations corresponding to the components are vectorized through the GPU broadcast mechanism. Combined with the separable convolution characteristics of the Gaussian filter, the three-dimensional convolution is decomposed into three one-dimensional convolution cascade operations, significantly reducing the computational complexity.

[0064] During the construction of the structural tensor, the element-level operations of the six components are performed through the GPU's Tensor Core to perform matrix multiplication and addition fusion calculations, and the FP16 mixed precision mode is used to increase the computing throughput to 83TFLOPS. The Gaussian filtering stage uses the separable convolution feature to decompose the three-dimensional filtering into three axial one-dimensional convolutions. With the 72MB on-chip storage of the L2 cache, the filtering speed reaches 120 million voxels per second. The feature decomposition link calls the batch function to reconstruct the matrix operation of the entire three-dimensional grid into a four-dimensional tensor. With the help of the 4070's 24MB L2 cache to achieve a 98% cache hit rate, the feature decomposition time of the 512^3 grid is controlled within 9.3 seconds, which is 17 times faster than the OpenBLAS multi-threaded implementation of EPYC 7k62.

[0065] The weight generation stage uses JIT-compiled CUDA kernels to parallelize conditional logic, using a 2.0x anisotropic enhancement factor in coherent areas and a 0.5x suppression factor in edge areas, and using 4070's 192 texture mapping units to accelerate mask operations. Actual measured data shows that it takes only 28 seconds to fully process a 512^3 geological model, of which 92% is GPU computing and 8% is CPU-assisted processing. The power consumption of the entire machine is stable at 320W, and the computing energy efficiency ratio reaches 8.75GFLOPS / W. This hybrid architecture gives full play to the advantages of parallel computing and EPYC's multi-core data supply capabilities, providing a production-level solution for large-scale 3D geological modeling.

[0066] Step 2: Use the multi-scale Gaussian filtering method to construct weight data through correlation judgment, solve the problem of blurred boundaries and invalid area interference of geological models through mask processing method, and provide input data for anisotropic inversion; define the anisotropic coefficient By eigenvalue ratio Quantify structural strength; construct a three-state weight function (edge ​​suppression / structure enhancement / default smoothing) to achieve dynamic regularized strength control based on seismic characteristics.

[0067] In step 2, the similarity S is calculated by defining Gaussian filters S1 and S2 for the horizontal and vertical structures, where the diffusion tensor D is determined by the feature vector and the weight:

[0068]

[0069] The similarity S is defined as:

[0070]

[0071] The calculation of S1 is divided into the following three steps, which enhances the response of horizontally continuous structures through a lateral-first filtering strategy:

[0072] (1) Perform lateral Gaussian filtering along the xy plane, and the anisotropic Gaussian kernel is shown as follows:

[0073]

[0074] Here σ x =σ y =2,σ z =0, indicating smoothing only in the xy plane;

[0075] (2) Perform longitudinal Gaussian filtering along the z direction:

[0076]

[0077] Here σ x =σ y =0,σ z =8, indicating filtering along the z direction;

[0078] (3) Final smoothing:

[0079]

[0080] σ3=2

[0081] Where ⊙ represents the multiplication of the corresponding elements of the matrix; I is the three-dimensional image data with dimensions (x, y, z); σ1 is (σ x =2,σ y =2,σ z =0), It is a three-dimensional Gaussian filter with σ1 as parameter, filtering is performed only in the x and y directions, and the z direction keeps the original value; is(σ x =2,σ y =2,σ z =2) three-dimensional isotropic Gaussian filter for smoothing in all directions;

[0082] The calculation of S2 is divided into the following three steps, which enhance the response of vertical discontinuities through longitudinal priority filtering:

[0083] ① Perform longitudinal Gaussian filtering along the yz plane:

[0084]

[0085] ② Here σ y =σ z =8,σ x =0, indicating smoothing only in the yz plane;

[0086] ③ Perform lateral Gaussian filtering along the x direction:

[0087]

[0088] Here σ y =σ z =0,σ x =2, indicating filtering along the x direction;

[0089] ④Final smoothing:

[0090]

[0091] Where ⊙ represents the multiplication of elements at corresponding positions in the matrix.

[0092] In the mask processing of step 2, the essence of the mask is a binary matrix (0 / 1 matrix), which is used to mark the valid area of ​​the data space and is highly consistent with terrain processing. The mathematical form is:

[0093]

[0094] Where M(x) is represented as a binary matrix related to spatial position, which is used to mark which positions (x) in the data space belong to the valid area (retain data) and which belong to the invalid area (suppress data), and selectively retain or suppress data through element multiplication operation:

[0095] D masked =D⊙M;

[0096] Where D is the original data matrix, D is the original data matrix, M is the binary matrix generated by the M(x) function, and has the same dimension as the original data matrix D. masked It represents the result of selectively retaining the original data D through the mask matrix M. ⊙ represents the multiplication of the elements at the corresponding positions of the matrix to ensure that the values ​​of the invalid areas are reset to zero to avoid affecting the subsequent calculations.

[0097] In order to suppress the gradient noise in the invalid area (air layer), the calculation of the mask is introduced in the gradient calculation as follows:

[0098]

[0099] A dynamic mask is introduced in the calculation of structural similarity. The boundary is kept clear by applying the mask multiple times. The calculation is expressed as:

[0100] Initial data mask:

[0101] D0=D⊙M;

[0102] Gaussian filter mask:

[0103] D1=G σ (D0);

[0104] Secondary mask:

[0105]

[0106] Through a multi-level and dynamic mask application strategy, problems such as blurred geological model boundaries and interference from invalid areas are effectively solved, providing high-quality input data for subsequent anisotropic inversion.

[0107] Step 3: Define the anisotropy coefficient by calculating the structural characteristics of the eigenvalue ratio corresponding to the main feature direction and the normal feature direction with reference to the correlation, and output the storage and visualization of the anisotropy weight in a modular form;

[0108] In step 3, the anisotropy and weight calculation process is as follows:

[0109] The anisotropy coefficient is defined by calculating the ratio of the eigenvalues ​​corresponding to the main feature direction and the normal feature direction:

[0110]

[0111] Where k is a positive real number used to adjust the eigenvalue ratio λ r Impact on the anisotropy coefficient, k→0: the anisotropy effect almost disappears, and the model degenerates into isotropic smoothness. The larger the k, the more sensitive the response of the region. The smaller the k, the weaker the anisotropy difference. The value of k is [0.5,2.0];

[0112] The weight values ​​of the following three constraints are obtained by calculating the similarity;

[0113] When the similarity judgment S is less than the value of the boundary judgment condition, edge_threshold is the edge area, and the smoothing strength is reduced in the edge area to retain the sharp edge;

[0114] When the similarity judgment S is greater than the coherence judgment condition value coherence_threshold, it is a structural coherent area. In the coherent area, the smoothing along the main direction is enhanced to suppress noise;

[0115] The other parts are set to 1:

[0116]

[0117] The inheritance of implicit masks is also taken into account in the generation of anisotropic weights: the input data for weight calculation has been preprocessed by the mask, and the anisotropic weight after mask processing is expressed as:

[0118]

[0119] Where f is the weight calculation method, λ i ,v iare the eigenvalue and eigenvector of the structure tensor respectively, and a and b are the assigned weight coefficients.

[0120] The structural tensor is extracted based on the seismic wave interlayer velocity, and the anisotropic regularized weight matrix is ​​constructed. The smoothing direction of the resistivity model is constrained by the direction of the main eigenvector (v1) to achieve spatial consistency control of the electrical boundary and the seismic interface.

[0121] Based on the SimPEG inversion framework, a new seismic structure regularization module (Regularization.SeismicStructure) is added to extract the directional characteristics of the formation interface through the seismic wave velocity structure tensor and construct anisotropic smoothing constraints, which improves the vertical resolution by 35% compared with the SimPEG default inversion (the measured fault boundary positioning error is ≤3%).

[0122] A multi-scale Gaussian filter-structure tensor fusion algorithm is proposed to replace the original isotropic diffusion filter of SimPEG, realize the efficient coupling of seismic characteristics and electrical models, and reduce false anomalies by 50%. The inversion processing of resistivity data includes the following steps:

[0123] Step S1: Calculate using the finite volume method, discretize the calculation area into non-overlapping control volumes or control elements, apply the conservation equation to integrate on the control volume, obtain a set of discrete equations by discretizing the integral equations, and obtain the required apparent resistivity by solving the discrete equations;

[0124] The theoretical basis and basic equations of the point power field in the three-dimensional high-density electrical numerical simulation are in a stable current field. Assume that there is a point power source with a current intensity of I at point A on the ground (the default setting in SimPEG is 1A): j is the current density, E is the electric field intensity, U is the total potential, and σ is the dielectric conductivity. Then the relationship between these physical quantities is:

[0125]

[0126] j=σE (2)

[0127] Substituting (2) into (1) we get:

[0128]

[0129] By solving equation (3), we can get:

[0130]

[0131] According to the principle of charge conservation, the relationship between current density j and charge density q is:

[0132]

[0133] The potential u and charge density equation are obtained as follows:

[0134]

[0135] In the three-dimensional numerical simulation, it is assumed that the position of the point charge e is (x0,y 0, z0), the charge density q is:

[0136] q=eδ(x-x0)δ(y-y0)δ(z-z0) (7)

[0137] Where δ is the Dirac function and the current intensity is

[0138]

[0139] The basic equation of the three-dimensional electric field potential is obtained from the point power source:

[0140]

[0141] For a dipole power source, the point power source A(x A ,y A ,z A ) and point power supply B(x B ,y B ,z B ) are I and -I respectively, and the basic equation satisfied by the three-dimensional electric field potential is:

[0142]

[0143] For the boundary value, on the surface G s The normal component of the current density is zero:

[0144] ∈Γ S (11)

[0145] Assuming that the abnormal body has no effect on the potential on the infinite boundary, at the infinite boundary G s The potential is obtained:

[0146] ∈Γ ∞ (12)

[0147] R is the distance from the point source to the boundary. The derivative of (12) is:

[0148] ∈Γ ∞ (13)

[0149] In step S1, in the grid generation and processing, the structured grid of the unit center method is used to discretize the three-dimensional calculation area, such as Figure 2As shown, according to the control volume method theory, the control volume of the unit where the node (i, j, k) is located is denoted by V i,j,k , the source term of the control equation discretized in the control volume is:

[0150]

[0151] In the formula, u represents the potential, σ represents the conductivity of the medium, and -s represents the source term, which represents the current injection or reception per unit volume;

[0152] Solve the volume integral of the above equation within the control volume:

[0153]

[0154] Applying Gauss's theorem on the left side yields:

[0155]

[0156] Control volume V i,j,k Indicates volume, is the surface area, and the area can be differentiated by the right formula to obtain:

[0157]

[0158] because Further simplifying:

[0159]

[0160] are the lengths of the control volume in the x, y, and z directions, respectively. It is expressed as:

[0161]

[0162] is the harmonic mean of the conductivity in the x direction of the two control volumes:

[0163]

[0164] In the above formula:

[0165]

[0166] therefore:

[0167]

[0168] The source term is:

[0169]

[0170] Step S2: The resistivity inversion under the constraint of seismic data is obtained by updating the framework of Gauss-Newton inversion. The linear system is solved by the conjugate gradient method PCG. At the same time, the adaptive regularization parameter under the reference of L-curve criterion is used for inversion. The root mean square error is used to evaluate the inversion effect.

[0171] The hybrid solution framework of the finite volume method (FVM) and Gauss-Newton method realizes the rapid generation and sparse storage of the forward stiffness matrix; the conjugate gradient method (PCG) optimization based on Jacobi preconditioning is combined with the dynamic adjustment strategy of the adaptive regularization parameters based on the L-curve criterion.

[0172] In the Inversion.BaseInvProblem class of SimPEG, the L-curve criterion is extended to drive the regularization parameter adaptive mechanism, dynamically adjust the smoothing intensity and data fitting weight, avoid manual trial and error parameter adjustment, and improve the stability of the inversion results by 25%.

[0173] Design β cooling strategy (β = β0 / γ n ), strong regularization is used to suppress noise in the early stage, and weak regularization is used to restore details in the later stage. Compared with the SimPEG fixed β value strategy, the root mean square error (RMS) is reduced by 18%.

[0174] In step S2, during model discretization and mesh generation, the governing equation is the current continuity equation:

[0175]

[0176] The finite volume method (FVM) is used for spatial discretization to generate the stiffness matrix A and obtain the linear system:

[0177] A(σ)φ=q;

[0178] In geophysical inversion, the simplest objective function is

[0179]

[0180] in:

[0181] m=[m1,m2,...,m M ] T ;

[0182] d=[d1,d2,...,d N ] T ;

[0183]

[0184] m is the model vector, d is the data vector, F(m) is the forward response of the model vector m, W dis the data item weight matrix with measured data, σ i is the standard deviation of the ith number;

[0185] The inversion objective function after adding Tikhonov regularization constraint is:

[0186]

[0187] The overall objective function used in the present invention is obtained by referring to Gauss-Newton inversion:

[0188]

[0189] The data fitting term Part of the calculation uses the L2 norm of the weighted residual, the regularization term Set as a combination of the minimum model constraint and the weighted directional smoothing constraint:

[0190]

[0191] The sensitivity matrix J[m] is calculated as element J ij It reflects the sensitivity of model parameters to data. Applying weights can balance the update amount of different parameters. The sensitivity matrix is ​​expressed as:

[0192]

[0193] The over-update of deep units is suppressed by dynamically adjusting the sensitivity weights, which are set as:

[0194]

[0195] Jacobi preconditioning is used to improve the matrix condition number. The matrix condition number is defined as:

[0196] κ(A)=||A||·||A -1 ||;

[0197] Here, ||·|| represents the matrix norm (usually the spectral norm):

[0198] The larger the condition number, the closer the matrix is ​​to being singular (irreversible), the more sensitive it is to errors during numerical solution, and the slower the iterative method converges;

[0199] The smaller the condition number, the closer the matrix is ​​to being "well-conditioned", the faster the iterative method converges and the more stable the solution.

[0200] In the inverse problem, each step of the Gauss-Newton method requires solving a linear system:

[0201]

[0202] The Hessian matrix is ​​approximated as

[0203]

[0204] Among them, the coefficient matrix (Hessian matrix) is usually ill-conditioned (high condition number), which leads to slow convergence of the conjugate gradient method (CG).

[0205] The goal of the preconditioner is to transform the original linear system into an equivalent but easier to solve form by constructing the preconditioner matrix M:

[0206] M -1 Ax=M -1 b;

[0207] Ideally, the new matrix M -1 The condition number κ(M -1 A)<<κ(A), thus accelerating the convergence of the iterative method.

[0208] The Jacobi preconditioner is the simplest preconditioner and is constructed as follows:

[0209] M=diag(A) -1 ;

[0210] That is, take the inverse of the diagonal elements of matrix A to form a diagonal matrix.

[0211] The constructed Jacobi preconditioner first calculates the diagonal elements h of the Hessian matrix ii , and then construct a diagonal matrix M as the preconditioner, whose elements are the reciprocals of the diagonal elements:

[0212]

[0213] Where the first term is the weighted column square sum of the sensitivity matrix J, and the second term is the regularization matrix W m The column sum of squares of is multiplied by β;

[0214] The constructed preconditioner M is:

[0215]

[0216] The above Hessian matrix can be further expressed as:

[0217]

[0218] The preconditioner M can be expressed as:

[0219]

[0220] The gradient is calculated as follows:

[0221]

[0222] The core of the Gauss-Newton method is to perform a first-order Taylor expansion on the nonlinear forward response F(m). The current model is m k ,F(m k+1 )=F(m k )+J k (m k+1 -m k );

[0223] Among them, J k is the Jacobian matrix, and its elements are expressed as:

[0224]

[0225] In the formula, F i Denotes the i-th forward response, m j Represents the jth model parameter, which is expressed as the partial derivative of the i-th observation data with respect to the j-th model parameter

[0226] The updated approximate objective function is:

[0227]

[0228] Where W d is the data item weight matrix with measured data;

[0229] For Φ(m k+1 )About m k+1 Taking the derivative and setting the reciprocal to zero gives us the linear system:

[0230]

[0231] The model update formula is:

[0232] m k+1 =m k +αδm;

[0233] in:

[0234] The obtained linear system is solved by the conjugate gradient method (CG) combined with the preconditioner-conjugate gradient method (PCG):

[0235] The initialization parameters are that the initial model update amount is 0, and the initial residual is equal to the gradient vector:

[0236]

[0237] In the formula, g represents the gradient vector, H represents the Hessian matrix;

[0238] Apply the preconditioning matrix M to the residual to accelerate convergence; set the initial search direction to the residual after preconditioning:

[0239]

[0240] In the formula, z0 represents the residual after preconditioning, r0 represents the initial residual, and p0 represents the initial search direction.

[0241] The step size of iterative update is α k Calculate the J-minimized quadratic function f(δm):

[0242]

[0243] Along direction p k The minimum point satisfies get:

[0244]

[0245] The numerator is the inner product of the preconditioned residual, which measures the current residual energy; r k represents the residual of the kth iteration, z k represents the residual after preconditioning. The distribution of the residual is adjusted by the preconditioning matrix M to improve the condition number of the linear system and accelerate convergence. The denominator is the energy of the search direction under the Hessian metric, so that the objective function is along p k The direction of maximum drop, p k Indicates the search direction of the kth iteration;

[0246] The subsequent model update and residual update are:

[0247] δm k+1 =δm k +α k p k ;

[0248] r k+1 =r k -α k p k ;

[0249] Step α along the search direction k , update the model correction; apply preconditioning and calculate the new direction, and process the new residuals through preconditioning to eliminate the ill-conditioning of the Hessian matrix:

[0250] z k+1 =M -1 r k+1 ;

[0251] Update Beta k Weighing the energy ratio of the new and old residuals, adjusting the search direction yields:

[0252]

[0253] Combine the current preconditioned residual with the direction of the previous step to generate the conjugate direction:

[0254] p k+1 =z k+1 +β k p k ;

[0255] In the formula, β k represents the adjustment coefficient of the conjugate direction;

[0256] The termination condition is set to terminate the iteration when the residual norm drops to 1% of the initial residual to balance computational efficiency and accuracy and avoid excessive iterations:

[0257]

[0258] In the formula, ||g|| represents the norm of the initial residual r0=g, which serves as the convergence criterion.

[0259] In the regularization parameter adaptation process of step S2, the L-curve criterion is used to determine the optimal value of the regularization parameter. The objective function of the regularization problem can be defined as:

[0260]

[0261] Where: A is the coefficient matrix, b is the observed data; L is the regularization matrix, λ is the regularization parameter;

[0262] The L-curve is defined as a parameterized curve:

[0263] (log||Ax λ -b||2,log||Lx λ ||2);

[0264] As λ changes, the curve shows the trade-off under different regularization strengths:

[0265] For the small λ region: the data fitting residual is small, the solution norm is large, the curve is approximately horizontal, and the residual changes slowly;

[0266] For large λ region: the solution norm is small, the residual increases rapidly, the curve is approximately vertical, and the solution norm changes slowly;

[0267] For the inflection point region, the curvature at the inflection point of the balance data fit and the smoothness of the solution is the largest, corresponding to the optimal λ;

[0268] The inflection point corresponds to the point with the largest curvature. The calculation steps are:

[0269] Parameterized curve: Consider the L-curve as a function η(ρ), where

[0270]

[0271] Calculate the curvature:

[0272]

[0273] Select the point of maximum curvature: the corresponding λ is the optimal value;

[0274] The L-curve criterion is used to realize the regularization parameter adaptation, and the data fitting difference Φ d With model complexity Φ m The trade-off curve selects the optimal β. For the generated anisotropic weights, traverse different λ values, run the inversion to calculate the residual and solution norm, and finally identify the L-curve inflection point to ensure that the inversion model fits the data while maintaining the structural characteristics. The initial β estimate is:

[0275]

[0276] In the formula, Φ d represents the poor fit of the data, and the denominator represents the anisotropic characteristics of the model complexity

[0277]

[0278] Where γ is the cooling factor, β k represents the regularization parameter of the kth iteration, and β is the cooling strategy. The purpose of the above settings is to use strong regularization (large β) to ensure stability in the early stage, and gradually reduce it in the later stage to improve the resolution.

[0279] The inversion is defined as the root mean square value (RMS) between the model response and the measured data:

[0280]

[0281] In the formula, represents the measured i-th data, represents the i-th data of the forward calculation, σ i is the measurement standard deviation of the i-th data, and N is the total number of data.

[0282] Step S3: the modular interface receives the anisotropic weights calculated and stored from the seismic data in step 3 to form a new regularization term to constrain inversion;

[0283] In step S3, the generation of the regularization term containing the seismic data constraint weights includes the following steps:

[0284] The projection of the main eigenvector in each direction is expressed as:

[0285]

[0286] The anisotropy strength is:

[0287]

[0288] The final x-direction weight is:

[0289]

[0290] In the formula, e x 、e y and e z Represents the standard basis vectors (x, y, z axis unit vectors)

[0291] For the direction vector v = (a, b, c), the regularization term in a single direction is decomposed into:

[0292]

[0293] In the formula, ∑|v i | is the absolute value and direction vector, a s is the model norm regularization parameter, α x , α y , α z Represents the regularization parameter of the gradient term in the x, y, and z directions:

[0294] ∑|v i |=|a|+|b|+|c|;

[0295] Anisotropic regularization is the sum of three orthogonal regularization terms:

[0296]

[0297] In the formula, is the regularization term, and m is the model vector.

[0298] Step S4: Generate and calculate the model, and visualize the results of resistivity inversion under the constraints of seismic data.

[0299] Model instance calculation

[0300] The proposed 3D inversion algorithm of DC method constrained by seismic data is based on the maximum principal eigenvector of the structure tensor. In order to verify the effectiveness of the guided method and the degree of fit between 3D Gauss-Newton inversion and 3D data constrained inversion, a simple 3D underground model with high-resistance and low-resistance anomalies was established in this section for inversion experiments.

[0301] Establishment of resistivity model: First, the present invention tests and establishes Figure 3The three-dimensional underground uniform half-space resistivity model with high and low resistivity anomalies is shown to test whether the constrained inversion can accurately restore the boundary of the target body. The surface terrain of the model is a basin terrain with high surroundings and low center. The model consists of a high-resistance sphere with a resistivity of 1000Ω·m (purple part) and a low-resistance area with a resistivity of 10Ω·m (yellow). The resistivity of the underground background part is 100Ω·m, and the air resistivity is 10 8 Ω·m, and the model takes into account the influence of terrain on the forward data and removes the air layer above the terrain.

[0302] The test network used in this solution test includes three test lines:

[0303] First, a main survey line 1 with a length of 176 m from x = -88 m to 88 m was arranged at y = 0 in the east-west direction to survey the abnormal display of the two abnormal bodies in the xoz plane;

[0304] Secondly, an auxiliary north-south survey line 2 with a length of 176 m from y = -88 m to 88 m was arranged at x = -20 m to survey the display of high-resistance spherical anomalies;

[0305] Finally, an auxiliary north-south survey line 3 with a length of 176 m from y = -88 m to 88 m was arranged at x = 20 m to survey the display of low-resistance layered anomalies.

[0306] For the above three survey lines, the device selected is a dipole-dipole device, using a 1A current source, the electrode spacing is set to 12m, and the number of receivers is set to 6. In order to eliminate the influence of terrain, each electrode is redefined to fall on the terrain. The modified open source software package SimPEG is used for three-dimensional finite volume numerical simulation to calculate the received signal data under the 1A current source under the dipole-dipole device. These data are added with random Gaussian noise with a mean of 0 and a standard deviation of 10% and will be inverted as synthetic data.

[0307] At the same time, the forward modeling line data is converted into a two-dimensional cross-sectional diagram to verify the reliability of the forward modeling data, and a three-dimensional pseudo cross-sectional diagram ( Figure 4 ), y = 0 and x = ± 20 ( Figure 5 ) The visualization part of the forward modeling data including the 2D cross-section diagram. It can be seen from both the 3D pseudo cross-section diagram and the 2D forward cross-section diagram that both high and low resistivity anomalies are reflected more obviously, and at the same time, the identification of multiple targets will not be greatly affected.

[0308] Build velocity model

[0309] First, the present invention tests and establishes Figure 6The three-dimensional underground uniform half-space seismic wave velocity model containing high-speed anomalies is shown to test whether the constrained inversion can accurately restore the boundary of the target body. The velocity model consists of a high-speed sphere (purple part) with a layer velocity of 200m / s and a high-speed layered cube with a layer velocity of 200m / s (yellow). The underground background layer velocity is 50m / s. At the same time, the model takes into account the influence of the terrain on the data and removes the air layer above the terrain. In order to verify whether the model is established correctly, the following is generated Figure 7 The 3D visualization image shown is used as a guide image for constraint.

[0310] Three-dimensional Gauss-Newton inversion of resistivity

[0311] In the numerical simulation experiment, since the approximate position of the target body is known, in order to reduce the number of grids and improve grid efficiency, this scheme adopts the method of intermediate area refinement in grid setting to establish the three-dimensional inversion grid.

[0312] Define the basic unit size as 4 as the reference length of the grid partition, and divide the x-direction into three areas: the coarse grid area on both sides: the cell length is 4, each containing 15 cells; the middle refined area: the cell length is halved to 2, containing 30 cells, and the total length is 60 cells. The y-direction maintains the same structure as the x-direction to form a symmetrical partition with the same total length and number of cells. The z-direction is divided into two areas: the lower coarse grid area: the grid length is 4, containing 10 cells, and the total length is 40; the upper refined area: the length is 2, containing 20 cells, and the total length is 40 cells. After the above refined and coarse part of the grid generation, a total of 60*60*30=108000 grids are generated. Considering the influence of the terrain, the remaining effective calculation grids are 100416 after removing the part above the terrain.

[0313] Calculating Constraint Weights for Seismic Data

[0314] This scheme regards the velocity model 3D image shown as prior underground structure information to guide the inversion of high-density electrical data. Its 3D image visualization is as follows: Figure 6 shown.

[0315] The structure tensor field of the three-dimensional image reveals the geological structural characteristics of different regions in the three-dimensional model through the shape of the characteristic ellipsoid and the eigenvalues ​​of the structure tensor. When the model is in a uniform medium area, the three eigenvalues ​​λ1, λ2, and λ3 are all close to zero. At this time, the characteristic ellipsoid degenerates into a tiny point structure near the origin, reflecting the lack of directional changes in the area. The corresponding code excludes non-geological interference through the active area mask. In the edge or linear structure area, the maximum eigenvalue λ1 is significantly larger than the other two eigenvalues, forming a cigar-shaped ellipsoid extending along the direction of the main eigenvector, and its major axis is perpendicular to the geological interface. The planar structure area is characterized by the first two eigenvalues ​​λ1 and λ2 being similar and much larger than λ3. The ellipsoid is compressed into a disk shape with the normal direction perpendicular to the bedding plane. At this time, the weight generation module will enhance the constraint strength in two directions in the plane. The quantitative analysis of the characteristic ellipsoid begins with the calculation of the gradient field. The spatial differential information is obtained through Gaussian derivative filtering, and then a structure tensor matrix containing six independent components is constructed. The eigenvalues ​​and eigenvectors extracted after eigendecomposition fully describe the local geometric characteristics.

[0316] Figure 8 , Fig. 9 and Fig.10 They are the cross-sectional comparison diagrams of the model and the maximum principal eigenvalue λ1 of the guide image at x=0, y=0 and depth z=-20m. The eigenvalues ​​at the edge of the target body are larger, while the eigenvalues ​​corresponding to the flat area are close to 0. In the surface terrain area, the eigenvalues ​​increase as they approach the surface.

[0317] Fig.11 , Fig.12 , Fig.13 and Fig.14 They are the similarity and anisotropic weight cross-sections of the guide image at x=0, y=0 and depth z=-10, -20m. The similarity and anisotropic weight are calculated using explicit and implicit masks to avoid the influence of terrain on the calculation.

[0318] As can be seen from the figure, the anisotropic weights show an increasing trend from the center to the edge of the target body, and finally extend to the background equal to 1, which corresponds to the good reflection of the center and edge of the target body; in terms of anisotropic weight performance, the weight coefficients defined by the above weights are assigned as a=2, b=0.5. To distinguish the boundaries of the two target bodies, the anisotropic weights at the center of the abnormal body are assigned a smaller weight, and the anisotropic weights on the boundary of the target body reflect the boundary recognition of the target body well.

[0319] In order to compare the inversion results, the constrained inversion and the traditional Gauss-Newton inversion will use the same inversion grid. The root mean square error (DRMS) is compared with the number of inversion iterations. Fig.15As shown, it can be seen that the root mean square errors of both inversion methods fluctuate down to within one percent to reach the convergence condition. The number of inversion iterations of the constrained inversion is more than that of the traditional inversion, but the time taken by the constrained inversion is slightly faster than that of the traditional inversion.

[0320] The comparison between the Gauss-Newton inversion and image-guided inversion results and the correct model on the y=0 interface is shown in Fig.16 , Fig.17 As shown in the figure, compared with the single method inversion results generated by traditional Gauss-Newton inversion, the three-dimensional image-guided inversion more clearly depicts the two target bodies, and is reliable for the identification of multiple target bodies. At the same time, from the perspective of the internal high-resistance / low-resistance morphology and extension of the two target bodies, the center and boundary of the high-resistance body are accurately depicted, which improves the vertical resolution; the center and upper and lower limits of the low-resistance body are more accurately depicted than traditional inversion, and are basically consistent with the real model.

[0321] Embodiment 1

[0322] The method and software of this solution can be applied to the three-dimensional precise inversion of geophysical exploration. The application part can be divided into two parts: theoretical application field and engineering application field as follows:

[0323] Theoretical application:

[0324] 1. Multi-source geophysical data fusion theory: This algorithm provides a new theoretical framework for geophysical joint inversion. By introducing the seismic wave velocity structure tensor to construct anisotropic regularization constraints, a mathematical and physical model for electrical-seismic data collaborative inversion is established, which promotes the development of multi-method data fusion theory.

[0325] 2. Optimization of three-dimensional inversion algorithm: The proposed finite volume method-Gauss-Newton hybrid solution framework, GPU-accelerated structure tensor calculation and adaptive regularization parameter strategy can be extended to other geophysical inversion fields such as electromagnetic method and gravity, improving the computational efficiency and stability of complex models.

[0326] 3. Structural constrained inversion theory: Through the directional coupling of seismic image feature extraction and electrical model, it provides a new idea for solving the multi-solution problem of single-method inversion, which is especially suitable for the refined characterization of structural features such as layered media and fault boundaries.

[0327] Engineering application: Compatible with SimPEG standard input formats (Survey, Mesh), support seamless connection with existing electrical and seismic processing processes (such as ResIPy, SeisPy), and reduce user migration costs; compatible with terrain correction modules, through the mesh deformation interface (Mesh.TensorMesh), realize the joint modeling of seismic and electrical data under undulating surface conditions, and solve the limitations of the default flat terrain assumption.

[0328] 1. Mineral resource exploration: It is suitable for accurate imaging of the three-dimensional electrical structure of metal ore bodies and oil and gas reservoirs. It can be combined with seismic velocity model constraints to improve the ability to identify ore body boundaries and reduce drilling verification costs. For example, in skarn deposits, the mineralized alteration zone can be effectively delineated through the joint analysis of resistivity anomalies and seismic wave velocity mutations.

[0329] 2. Hydrogeological survey: It is used for the detailed characterization of the spatial distribution of aquifers and the boundaries of groundwater pollution plumes. Combined with the constraints of seismic reflection interfaces, it can distinguish the contact relationship between low-resistance aquifers and high-resistance bedrock, and improve the accuracy of water resource assessment.

[0330] 3. Engineering geological survey:

[0331] Tunnels and underground engineering: Identify geological hazards such as karst pipes and fault zones, enhance the vertical resolution of resistivity inversion through seismic layer constraints, and assist in engineering geological stratification.

[0332] Dam foundation and slope stability: Detect weak interlayers and seepage channels, and evaluate rock mass integrity by combining the anisotropic characteristics of seismic wave velocity.

[0333] 4. Environment and Hazard Geology:

[0334] Contaminated site detection: define the spread of pollutants, use seismic interface constraints to distinguish artificial fill from native strata, and improve the vertical positioning accuracy of contaminated boundaries.

[0335] Goaf and ground fissures: Accurately identify underground cavities and rock dislocation zones through spatial coupling of resistivity anomalies and seismic wave velocity drop zones.

[0336] 5. Urban underground space development: It is suitable for shallow electrical structure detection in subway tunnels, pipeline corridors, etc. Combined with seismic profile data constraints, it can effectively distinguish the interface between artificial structures and natural strata and reduce the impact of urban noise interference.

[0337] Therefore, the present invention adopts the above-mentioned direct current electrical three-dimensional inversion method constrained by seismic data, applies the seismic data constraint to the three-dimensional resistivity inversion, performs numerical simulation and inversion processing on the electrical data under the undulating terrain according to the finite volume method numerical simulation and the Gauss-Newton inversion method, and introduces the seismic wave layer velocity tensor structure constraint for inversion, thereby reducing the multi-solution problem of inversion by a single method, obtaining a more accurate three-dimensional underground electrical structure, and improving the vertical boundary resolution of the inversion.

[0338] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solution of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solution to deviate from the spirit and scope of the technical solution of the present invention.

Claims

1. A three-dimensional inversion method of direct current method constrained by seismic data, characterized by: Including processing of seismic data and inversion of resistivity data; The processing of seismic data includes the following steps: Step 1: Process the 3D seismic velocity data or 3D seismic images related to the electrical exploration data, calculate the structure tensor after removing the air layer, obtain the structure tensor data, and accelerate the structure tensor calculation and feature decomposition through GPU and CPU hybrid calculation; Step 2: Use the multi-scale Gaussian filtering method to construct weight data through correlation judgment, and use the mask processing method to solve the problem of blurred terrain boundaries and invalid area interference in the geological model to provide input data for anisotropic inversion; Step 3: Define the anisotropy coefficient by calculating the structural characteristics of the eigenvalue ratio corresponding to the main feature direction and the normal feature direction with reference to the correlation, and output the storage and visualization of the anisotropy weight in a modular form; The inversion processing of resistivity data includes the following steps: Step S1: Calculate using the finite volume method, discretize the calculation area into non-overlapping control volumes or control elements, apply the conservation equation to integrate on the control volume, obtain a set of discrete equations by discretizing the integral equations, and obtain the required apparent resistivity by solving the discrete equations; Step S2: The resistivity inversion under the constraint of seismic data is obtained by updating the framework of Gauss-Newton inversion. The linear system is solved by the conjugate gradient method PCG. At the same time, the adaptive regularization parameter under the reference of L-curve criterion is used for inversion. The root mean square error is used to evaluate the inversion effect. Step S3: the modular interface receives the anisotropic weights calculated and stored from the seismic data in step 3 to form a new regularization term to constrain inversion; Step S4: Generate and calculate the model and visualize the results of resistivity inversion under the constraints of seismic data.

2. A three-dimensional inversion method of direct current method constrained by seismic data according to claim 1, characterized in that: In step 1, the structure tensor data is calculated as follows: The gradient field is calculated using central difference combined with Gaussian derivative kernel, and the three-dimensional discrete gradient operator is defined as: Remove high-frequency noise using Gaussian smoothing: Gaussian kernel G σ The discrete form of is: The structure tensor is a second-order tensor matrix derived from the gradient. The mathematical form of the three-dimensional structure tensor is: In the formula, I x ,I y and I z are the gradients in each direction, S σ is the structure tensor matrix corresponding to the node (i, j, k), * is the convolution operator, G σ is a two-dimensional Gaussian difference operator, σ is the control Gaussian function; The expression of the Gaussian difference operator is: S σ =λ1v1v1 T +λ2v2v2 T +λ3v3v3 T (λ1>λ2>λ3); Structure tensor matrix S σ is a 3×3 real-opinion matrix, decomposed into three eigenvalues ​​(λ1, λ2, λ3) and three pairwise orthogonal eigenvectors (v1, v2, v3), where the eigenvector v1 corresponding to the largest eigenvalue represents the direction with the strongest local gradient energy, perpendicular to the structural boundary; In continuous media, v1 is perpendicular to the stratum interface; in sharp edges, v1 is perpendicular to the fault plane; the second eigenvector v2 is in a plane orthogonal to v1, indicating the secondary change direction of the local structure; the third eigenvector v3 points to the direction of the smallest local change, usually parallel to the main extension direction of the structure, λ1>>λ3, indicating strong anisotropy; λ1≈λ3, indicating isotropy.

3. A three-dimensional inversion method of direct current method constrained by seismic data according to claim 2, characterized in that: In step 1, GPU and CPU hybrid computing accelerates the structure tensor calculation and feature decomposition based on the hardware environment of NVIDIA GeForce RTX4070 graphics card. The calculation process is as follows: 1.1: Migrate the velocity model data from the host memory to the graphics card memory through the asynchronous transmission channel of the CuPy library, and use the unified memory architecture of CUDA8.0 to achieve zero-copy data interaction between the CPU and GPU; 1.2: The gradient calculation stage uses the parallel implementation of the 3D central difference method. The CUDA thread block is responsible for the differential operation of the 32×32×32 voxel area, and the global memory access latency is reduced by caching local data in shared memory; 1.3: During the construction of the structure tensor, the calculation of the six independent components (Txx, Txy, Txz, Tyy, Tyz, Tzz) is mapped to different streaming multiprocessors SM for parallel execution. The element-level operations corresponding to the components are vectorized through the GPU broadcast mechanism. Combined with the separable convolution characteristics of the Gaussian filter, the three-dimensional convolution is decomposed into three cascade operations of one-dimensional convolution.

4. A three-dimensional inversion method of direct current method constrained by seismic data according to claim 3, characterized in that: In step 2, the similarity S is calculated by defining Gaussian filters S1 and S2 for the horizontal and vertical structures, where the diffusion tensor D is determined by the feature vector and the weight: The similarity S is defined as: The calculation of S1 is divided into the following three steps, which enhances the response of horizontally continuous structures through a lateral-first filtering strategy: (1) Perform lateral Gaussian filtering along the xy plane. The anisotropic Gaussian kernel is shown as follows: Here σ x =σ y =2,σ z =0, indicating smoothing only in the xy plane; (2) Perform longitudinal Gaussian filtering along the z direction: Here σ x =σ y =0,σ z =8, indicating filtering along the z direction; (3) Final smoothing: Where ⊙ represents the multiplication of the corresponding elements of the matrix; I is the three-dimensional image data with dimensions (x, y, z); σ1 is (σ x =2,σ y =2,σ z =0), It is a three-dimensional Gaussian filter with σ1 as parameter, filtering is performed only in the x and y directions, and the z direction keeps the original value; is(σ x =2,σ y =2,σ z =2) three-dimensional isotropic Gaussian filter for smoothing in all directions; The calculation of S2 is divided into the following three steps, which enhance the response of vertical discontinuities through longitudinal priority filtering: ① Perform longitudinal Gaussian filtering along the yz plane: Here σ y =σ z =8,σ x =0, indicating smoothing only in the yz plane; ② Perform lateral Gaussian filtering along the x direction: Here σ y =σ z =0,σ x =2, indicating filtering along the x direction; ③Final smoothing: σ6=8 Where ⊙ represents the multiplication of elements at corresponding positions in the matrix.

5. A three-dimensional inversion method of direct current method constrained by seismic data according to claim 4, characterized in that: In the mask processing of step 2, the essence of the mask is a binary matrix (0 / 1 matrix), which is used to mark the valid area of ​​the data space. The mathematical form is: Where M(x) is represented as a binary matrix associated with spatial position, and data is selectively retained or suppressed through element-wise multiplication operation: D masked =D⊙M; Where D is the original data matrix, M represents the binary matrix generated by the M(x) function, which has the same dimension as the original data matrix D. masked It represents the result of selectively retaining the original data D through the mask matrix M; The calculation of the mask introduced in the gradient calculation is: A dynamic mask is introduced in the calculation of structural similarity. The boundary is kept clear by applying the mask multiple times. The calculation is expressed as: Initial data mask: D0=D⊙M; Gaussian filter mask: D1=G σ (D0); Secondary mask: <h2 style=";text-align:left;direction:ltr">D2<h2 style=";text-align:left;direction:ltr"> masked <h2 style=";text-align:left;direction:ltr"> (D1⊙M) 6. A three-dimensional inversion method of direct current method constrained by seismic data according to claim 5, characterized in that: In step 3, the anisotropy and weight calculation process is as follows: The anisotropy coefficient is defined by calculating the ratio of the eigenvalues ​​corresponding to the main feature direction and the normal feature direction: Where k is a positive real number used to adjust the eigenvalue ratio λ r Impact on the anisotropy coefficient, k→0: the anisotropy effect almost disappears, and the model degenerates into isotropic smoothness. The larger the k, the more sensitive the response of the region. The smaller the k, the weaker the anisotropy difference. The value of k is [0.5,2.0]; The weight values ​​of the following three constraints are obtained by calculating the similarity; When the similarity judgment S is less than the value of the boundary judgment condition, edge_threshold (which can be set as a constant) is the edge area, and the smoothing strength is reduced in the edge area to retain the sharp edge; When the similarity judgment S is greater than the coherence judgment condition value coherence_threshold (a constant that can be set), it is a structural coherent area, and the smoothing along the main direction is enhanced in the coherent area to suppress noise; The other parts are set to 1: The input data for weight calculation has been preprocessed by mask, and the anisotropic weight after mask processing is expressed as: Where f is the weight calculation method, λ i ,v i are the eigenvalue and eigenvector of the structure tensor respectively, and a and b are the assigned weight coefficients.

7. A three-dimensional inversion method of direct current method constrained by seismic data according to claim 6, characterized in that: In step S1, according to the control volume method theory, the control volume of the unit where the node (i, j, k) is located is denoted by V i,j,k , the source term of the control equation discretized in the control volume is: In the formula, u represents the potential, σ represents the dielectric conductivity, and -s is the source term representing the current injection or reception per unit volume; Solve the volume integral of the above equation within the control volume: Applying Gauss's theorem on the left side yields: Control volume V i,j,k Indicates volume, is the surface area, and the area can be differentiated by the right formula to obtain: because Further simplifying: are the lengths of the control volume in the x, y, and z directions, respectively. It is expressed as: is the harmonic mean of the conductivity in the x direction of the two control volumes: In the above formula: therefore: The source term is:

8. A three-dimensional inversion method of direct current method constrained by seismic data according to claim 7, characterized in that: In step S2, the core of the Gauss-Newton method is to perform a first-order Taylor expansion on the nonlinear forward response F(m). The current model is m k , F(m k+1 )=F(m k )+J k (m k+1 -m k ); Among them, J k is the Jacobian matrix, and its elements are expressed as: In the formula, F i Denotes the i-th forward response, m j represents the jth model parameter, J ij is the partial derivative of the i-th observation data with respect to the j-th model parameter; The updated approximate objective function is: Where W d is the data item weight matrix with measured data; For Φ(m k+1 )About m k+1 Taking the derivative and setting the reciprocal to zero gives us the linear system: The model update formula is: m k+1 =m k +αδm; in: The obtained linear system is solved by the conjugate gradient method (CG) combined with the preconditioner-conjugate gradient method (PCG): The initialization parameters are that the initial model update amount is 0, and the initial residual is equal to the gradient vector: In the formula, g represents the gradient vector, H represents the Hessian matrix; Apply the preconditioning matrix M to the residual to accelerate convergence; set the initial search direction to the residual after preconditioning: In the formula, z0 represents the residual after preconditioning, r0 represents the initial residual, and p0 represents the initial search direction. The step size of iterative update is α k Calculate J to minimize the quadratic function f(δm: Along direction p k The minimum point satisfies Get the step size α k : The numerator is the inner product of the preconditioned residual, which measures the current residual energy; r k represents the residual of the kth iteration; z k represents the residual after preconditioning. The distribution of the residual is adjusted by the preconditioning matrix M to improve the condition number of the linear system and accelerate convergence. The denominator is the energy of the search direction under the Hessian metric, p k Indicates the search direction of the kth iteration; The subsequent model update and residual update are: δm k+1 =δm k +a k p k ; r k+1 ＝r k -α k p k ； Step α along the search direction k , update the model correction; apply preconditioning and calculate the new direction, and process the new residuals through preconditioning to eliminate the ill-conditioning of the Hessian matrix: z k+1 =M -1 r k+1 ; Update Beta k Weighing the energy ratio of the new and old residuals, adjusting the search direction yields: Combine the current preconditioned residual with the previous step direction to generate the conjugate direction: p k+1 =from k+1 +β k p k ; In the formula, β k represents the adjustment coefficient of the conjugate direction; The termination condition is set to terminate the iteration when the residual norm drops to 1% of the initial residual to balance computational efficiency and accuracy and avoid excessive iterations: Where g represents the norm of the initial residual r0=g, which serves as the convergence criterion.

9. A three-dimensional inversion method of direct current method constrained by seismic data according to claim 8, characterized in that: In the regularization parameter adaptation process of step S2, the L-curve criterion is used to determine the optimal value of the regularization parameter. The objective function of the regularization problem can be defined as: Where: A is the coefficient matrix, b is the observed data; L is the regularization matrix, λ is the regularization parameter; The L-curve is defined as a parameterized curve: (log||Ax λ -b||2,log||Lx λ ||2); As λ changes, the curve shows the trade-off under different regularization strengths: For the small λ region: the data fitting residual is small, the solution norm is large, the curve is approximately horizontal, and the residual changes slowly; For large λ region: the solution norm is small, the residual increases rapidly, the curve is approximately vertical, and the solution norm changes slowly; For the inflection point region, the curvature at the inflection point of the balance data fit and the smoothness of the solution is the largest, corresponding to the optimal λ; The inflection point corresponds to the point with the largest curvature. The calculation steps are: Parameterized curve: Consider the L-curve as a function η(ρ), where Calculate the curvature: Select the point of maximum curvature: the corresponding λ is the optimal value; The L-curve criterion is used to realize the regularization parameter adaptation, and the data fitting difference Φ d With model complexity Φ m The trade-off curve selects the optimal β. For the generated anisotropic weights, traverse different λ values, run the inversion to calculate the residual and solution norm, and finally identify the L-curve inflection point to ensure that the inversion model fits the data while maintaining the structural characteristics. The initial β estimate is: In the formula, Φ d represents the poor fit of the data, and the denominator represents the anisotropic characteristics of the model complexity Where γ is the cooling factor, β k represents the regularization parameter of the kth iteration, β is the cooling strategy; The inversion is defined as the root mean square value (RMS) between the model response and the measured data: In the formula, represents the measured i-th data, represents the i-th data of the forward calculation, σ i is the measurement standard deviation of the i-th data, and N is the total number of data.

10. A three-dimensional inversion method of direct current method constrained by seismic data according to claim 9, characterized in that: in step S3, generating a regularization term containing seismic data constraint weights comprises the following steps: The projection of the main eigenvector in each direction is expressed as: The anisotropy strength is: The final x-direction weight is: In the formula, e x 、e y and e z Represents the standard basis vector (x, y, z axis unit vector) For the direction vector v = a, b, c, the regularization term in a single direction is decomposed as: In the formula, ∑|v i | is the absolute value and direction vector, α s is the model norm regularization parameter, α x , α y , α z Represents the regularization parameter of the gradient term in the x, y, and z directions; ∑|v i |=|a|+|b|+|c|; Anisotropic regularization is the sum of three orthogonal regularization terms: In the formula, is the regularization term, and m is the model vector.

Images