A 3D inversion method based on direct current method constrained by seismic data

Abstract

The present invention discloses a three-dimensional inversion method for direct current electrical methods constrained by seismic data, which relates to the field of geophysical exploration and includes processing seismic data and inverting resistivity data. Using seismic wave layer velocity information to constrain high-density electrical inversion can play a complementary role to a certain extent, facilitating the use of three-dimensional seismic image information to supplement the resolution of resistivity inversion, reflecting the boundaries of underground electrical anomalies, and further characterizing the boundary characteristics of the anomalies through the integration of geological and geophysical information. Seismic data constraints are applied to three-dimensional resistivity inversion, and electrical data under undulating terrain are numerically simulated and inverted using the finite volume method and the Gauss-Newton inversion method. Seismic anisotropic velocity tensor constraints are also introduced for inversion, thereby reducing the multi-solution problem of single-method inversion, obtaining a more accurate three-dimensional underground electrical structure, and improving the resolution of the inversion for vertical boundaries.

Description

Technical Field

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

[0002] Electrical prospecting is a classic geophysical exploration method. Usually, the results of high-density electrical inversion are generally output as two-dimensional apparent resistivity profiles. However, for complex geological conditions, two-dimensional inversion cannot meet exploration requirements. 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 prospecting itself need to be solved urgently: the resistivity inversion algorithm needs to be improved, and the requirements for algorithm and hardware conditions for three-dimensional inversion calculation are increased; the vertical resolution of the underground resistivity distribution results of high-density electrical data inversion is not high, making it difficult to distinguish the boundaries of anomalies, which is insufficient for grasping the boundaries of anomalies.

[0003] Joint inversion is currently one of the primary methods for improving inversion quality. Joint inversion involves the combined application of multiple geophysical observations in geophysical inversion, leveraging the interrelationships between the petrophysical and geometric parameters of geological bodies to jointly invert the same subsurface geological and geophysical model. However, while joint inversion can improve the interpretation accuracy of single electromagnetic data by leveraging other geophysical exploration data, it also requires a significant increase in computational effort, placing higher demands on computer computing power and memory storage, making joint inversion relatively less practical.

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

[0005] To address the aforementioned challenges, the present invention provides a three-dimensional inversion method for direct current (DC) electrical resistivity constrained by seismic data. Using seismic wave layer velocity to constrain high-density electrical inversion can complement high-density electrical inversion to a certain extent, facilitating the use of image information to supplement the vertical resolution of resistivity images, invert the boundaries of underground electrical anomalies, and further characterize the boundary characteristics of the anomalies through the integration of geological and geophysical information. This invention studies three-dimensional inversion of direct current (DC) apparent resistivity constrained by seismic wave layer velocity. The present invention applies seismic data constraints to three-dimensional resistivity inversion. Finite volume numerical simulation and Gauss-Newton inversion methods are used to numerically simulate and invert electrical data over undulating terrain. Seismic wave layer velocity tensor structure constraints are also introduced for inversion. This reduces the multi-solution problem of single-method inversion, resulting in a more accurate three-dimensional underground electrical structure and improved vertical boundary resolution of the inversion.

[0006] Therefore, the present invention adopts the above-mentioned seismic data-constrained DC electrical method 3D inversion method, which has the following beneficial effects:

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

[0008] (2) The present invention introduces structural constraints driven by seismic data, breaking the equivalence dilemma of electrical inversion, reducing the model space solution set by 40%-60%, and significantly 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 memory usage by 35%.

[0010] (4) In the present invention, GPU accelerates structural tensor calculation and eigendecomposition, and the operation speed of key modules is increased by 8-12 times compared with 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 with 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 This is the overall technical roadmap for 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 the forward modeling of the three-dimensional finite volume method using the unit center method in an embodiment of the present invention;

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

[0018] Figure 4 3D pseudo cross-sectional view of apparent resistivity of the 3D resistivity model in an embodiment of the present invention;

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

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

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

[0022] Figure 8 The velocity model and the maximum principal eigenvalue at x=0 in the embodiment of the present invention are Cross-section comparison diagram;

[0023] Figure 9 The velocity model and the maximum principal eigenvalue at y=0 in the embodiment of the present invention are Cross-section comparison diagram;

[0024] Figure 10 The velocity model and the maximum principal eigenvalue at z = -20m in the embodiment of the present invention are shown in FIG. Cross-section comparison diagram;

[0025] Figure 11 This 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] Figure 12 This 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] Figure 13 This 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] Figure 14 This 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] Figure 15 Graph showing the variation of the root mean square error with the number of inversion iterations in an embodiment of the present invention;

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

[0031] Figure 17 This is a three-dimensional resistivity inversion cross-section diagram at y=0m based on 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 with reference to the accompanying drawings and embodiments.

[0033] Unless otherwise defined, technical or scientific terms used in the present invention shall have the same meaning as commonly understood by one of ordinary skill in the art to which the present invention belongs.

[0034] The words “include” or “comprising” and similar words 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 accompanying drawings. It 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 stipulated and limited, the terms such as “attachment” should be understood in a broad sense. For example, it can be a fixed connection, a detachable connection, or an integral whole; 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 the 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 associated with the electrical exploration data, remove the air layer, and calculate the structural tensor to obtain structural tensor data. This data is then accelerated through GPU and CPU hybrid computing to calculate the structural tensor and perform eigendecomposition.

[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. The three-dimensional discrete gradient operator is defined as:

[0041] ;

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

[0043] ;

[0044] Gaussian kernel 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] Where, 、 and are the gradients in each direction, is the structure tensor matrix corresponding to node (i, j, k), is the convolution operator, is a two-dimensional Gaussian difference operator, To control the Gaussian function;

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

[0050] ;

[0051] ;

[0052] Structure tensor matrix For one The real pair matrix of is decomposed into three eigenvalues ( 、 、 ) and three pairwise orthogonal eigenvectors ( 、 、 ), where the maximum eigenvalue corresponds to The eigenvector of represents the direction pointing to the strongest local gradient energy, perpendicular to the structure boundary;

[0053] In a continuous medium, perpendicular to the stratum interface; in sharp edges, Perpendicular to the fault plane; second eigenvector In and The orthogonal plane represents the secondary change direction of the local structure, and the third eigenvector points in the direction of minimum local variation, usually parallel to the main extension direction of the structure, , indicating strong anisotropy; , indicating isotropy.

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

[0055] GPU-CPU heterogeneous parallel computing: Computationally intensive tasks such as structural tensor extraction and feature decomposition are accelerated by the GPU, while the core inversion logic is scheduled by CPU multi-threading; multi-scale Gaussian filtering replaces the traditional diffusion equation, reducing the complexity of three-dimensional similarity calculations 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 to 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 algorithm stability.

[0058] Memory and efficiency balance:

[0059] The SimPEG stiffness matrix is reconstructed through sparse matrix compression storage (CSR format), reducing 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).

[0060] The calculation process is as follows:

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

[0062] 1.2: The gradient calculation stage uses a parallel implementation of the 3D central difference method. CUDA thread blocks are responsible for differential operations in a 32×32×32 voxel area. Shared memory is used to cache local data to reduce global memory access latency.

[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 (SMs) 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 a cascade of three one-dimensional convolution operations, significantly reducing the computational complexity.

[0064] During the construction of the structural tensor, element-level operations on the six components are performed using the GPU's Tensor Cores, performing matrix multiplication and addition fusion calculations. FP16 mixed-precision mode is used to increase computational throughput to 83 TFLOPS. The Gaussian filtering stage leverages the separable convolution feature to decompose the three-dimensional filter into three axial one-dimensional convolutions. Combined with the 72MB on-chip L2 cache, this enables filtering at a speed of 120 million voxels per second. The eigendecomposition stage calls a batch function to reconstruct the matrix operations of the entire three-dimensional mesh into a four-dimensional tensor. Leveraging the 4070's 24MB L2 cache to achieve a 98% cache hit rate, the eigendecomposition time for a 512^3 mesh is kept within 9.3 seconds, 17 times faster than the EPYC 7k62's OpenBLAS multi-threaded implementation.

[0065] The weight generation stage uses JIT-compiled CUDA kernels to parallelize conditional logic. A 2.0x anisotropic enhancement factor is used in coherent regions, while a 0.5x suppression factor is used in edge regions. Masking operations are accelerated using the 4070's 192 texture mapping units. Measured data shows that a complete 512^3 geological model is processed in just 28 seconds, with 92% GPU computing and 8% CPU-assisted processing. Total system power consumption remains stable at 320W, achieving a computational energy efficiency of 8.75 GFLOPS / W. This hybrid architecture leverages the advantages of parallel computing and the multi-core data feed capabilities of the EPYC processor, providing a production-grade solution for large-scale 3D geological modeling.

[0066] 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 geological model boundaries and invalid area interference, providing input data for anisotropic inversion; define the anisotropy coefficient , through the 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 earthquake characteristics.

[0067] In step 2, by defining Gaussian filters for horizontal and vertical structures and , calculate similarity , where the diffusion tensor Determined by the eigenvector and weight:

[0068] ;

[0069] Similarity Defined as:

[0070] ;

[0071] The calculation is divided into the following three steps, using a horizontal priority filtering strategy to enhance the response of horizontally continuous structures:

[0072] (1) Perform transverse Gaussian filtering along the xy plane. The anisotropic Gaussian kernel is shown as follows:

[0073] ;

[0074] Here , indicating smoothing only in the xy plane;

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

[0076] ;

[0077] Here , indicating filtering along the z direction;

[0078] (3) Final smoothing:

[0079]

[0080] Where, Represents the multiplication of elements at corresponding positions of the matrix; is three-dimensional image data with dimensions (x, y, z); for( , )'s Gaussian kernel parameters, Therefore A three-dimensional Gaussian filter with parameters , which filters only in the x and y directions and keeps the original value in the z direction; for( , ) three-dimensional isotropic Gaussian filter for smoothing in all directions;

[0081] The calculation of is divided into the following three steps, which enhance the response of vertical discontinuities by longitudinal-first filtering:

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

[0083] ;

[0084] ②Here , indicating smoothing only on the yz plane;

[0085] ③ Perform horizontal Gaussian filtering along the x direction:

[0086] ;

[0087] Here , indicating filtering along the x direction;

[0088] ④Final smoothing:

[0089] ;

[0090] Where, Represents the multiplication of elements at corresponding positions in the matrix.

[0091] In the mask processing process 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:

[0092] ;

[0093] In the formula Represented as a binary matrix related to spatial positions), it 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 operations:

[0094] ;

[0095] Where, is the original data matrix is the original data matrix, Expressed as The binary matrix generated by the function is the same as the original data matrix The same dimensions, Represented by the mask matrix The result of selectively retaining the original data D is: Represents the multiplication of the elements at corresponding positions in the matrix to ensure that the values in the invalid area are reset to zero to avoid affecting subsequent calculations.

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

[0097] ;

[0098] A dynamic mask is introduced in the structural similarity calculation to keep the boundary clear by applying the mask multiple times. The calculation is expressed as:

[0099] Initial data mask:

[0100] ;

[0101] Gaussian filter mask:

[0102] ;

[0103] Secondary mask:

[0104] ;

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

[0106] Step 3: Define the anisotropy coefficient by referring to the structural characteristics of the correlation between the eigenvalue ratios corresponding to the main feature direction and the normal feature direction, and store and visualize the anisotropy weights in a modular form;

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

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

[0109] ;

[0110] ;

[0111] Where k is a positive real number used to adjust the eigenvalue ratio 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 similarity is determined When the value of edge_threshold is less than the boundary judgment condition, it is an edge area, and the smoothing intensity is reduced in the edge area to retain the sharp edge;

[0114] When similarity is determined The value greater than the coherence judgment condition is coherence_threshold, which is a structural coherence area. In the coherence area, 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 pre-processed by the mask, and the anisotropic weight after mask processing is expressed as:

[0118] ;

[0119] in is the weight calculation method, are the eigenvalue and eigenvector of the structure tensor respectively, and a and b are the assigned weight coefficients.

[0120] The structure tensor is extracted based on the seismic wave velocity, and the anisotropic regularization weight matrix is constructed; the direction of the main eigenvector ( ) constrains the smoothing direction of the resistivity model to achieve spatial consistency control of the electrical boundary and seismic interface.

[0121] Based on the SimPEG inversion framework, a new seismic structure regularization module (Regularization.SeismicStructure) is added. It extracts the directional characteristics of the stratum interface through the seismic wave velocity structure tensor and constructs anisotropic smoothing constraints. This improves the vertical resolution by 35% compared to 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, achieving efficient coupling of seismic characteristics and electrical models, reducing 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, and obtain a set of discrete equations by discretizing the integral equations. The required apparent resistivity is obtained by solving the discrete equations.

[0124] The theoretical basis and basic equations of the point source field in the three-dimensional high-density electrical numerical simulation are based on a steady current field. Assume that there is a point 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. The relationship between these physical quantities is:

[0125] (1)

[0126] (2)

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

[0128] (3)

[0129] Solving equation (3) yields:

[0130] (4)

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

[0132] (5)

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

[0134] (6)

[0135] In the three-dimensional numerical simulation, it is assumed that the position of the point charge e is , the charge density q is:

[0136] (7)

[0137] in is the Dirac function, and the current intensity is

[0138] (8)

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

[0140] (9)

[0141] For the dipole power supply, the point power supply and point power supply The currents are and , the basic equation satisfied by the three-dimensional electric field potential is:

[0142] (10)

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

[0144] (11)

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

[0146] (12)

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

[0148] (13)

[0149] In step S1, in the grid division 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 recorded as , the source term of the control equation discretized in the control volume is:

[0150] ;

[0151] Where u represents the potential, represents the conductivity of the medium, The source term represents the injection or reception of current per unit volume;

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

[0153] ;

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

[0155] ;

[0156] Within the control volume Indicates volume, is the surface area, and the area differentiation of the right equation is simplified to:

[0157] ;

[0158] because , = , , further simplified to:

[0159] ;

[0160] , , are the lengths of the control volume in the x, y, and z directions, respectively, where Expressed as:

[0161] ;

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

[0163] ;

[0164] In the above formula:

[0165] ;

[0166] therefore:

[0167] ;

[0168] ;

[0169] ;

[0170] The source term is:

[0171] .

[0172] Step S2: Obtain resistivity inversion under seismic data constraints using the Gauss-Newton inversion framework. The linear system is solved using the conjugate gradient method (PCG). Inversion is performed using an adaptive regularization parameter under the L-curve criterion. The root mean square error (RMS) is used to evaluate the inversion effect.

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

[0174] In SimPEG's Inversion.BaseInvProblem class, the L-curve criterion is extended to drive the regularization parameter adaptation mechanism, dynamically adjusting the smoothing strength and data fitting weight, avoiding manual trial and error parameter adjustment, and improving the stability of the inversion results by 25%.

[0175] Designing a beta cooling strategy , 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%.

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

[0177] ;

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

[0179] ;

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

[0181] ;

[0182] in:

[0183] ;

[0184] ;

[0185] ;

[0186] is the model vector, is the data vector, is the forward response of the model vector m, is the data item weight matrix related to the measured data, For the The standard deviation of the number;

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

[0188] ;

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

[0190] ;

[0191] 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:

[0192] ;

[0193] Sensitivity Matrix Calculate as element 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:

[0194]

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

[0196] ;

[0197] Jacobi preconditioning is used to improve the matrix condition number. The matrix condition number (Condition Number) is defined as:

[0198] ;

[0199] in, Represents the matrix norm (usually the spectral norm):

[0200] The larger the condition number, the closer the matrix is to being singular (irreversible), the more sensitive it is to errors in numerical solutions, and the slower the iterative method converges.

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

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

[0203]

[0204] The Hessian matrix is approximately

[0205]

[0206] 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).

[0207] 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:

[0208]

[0209] Ideally, the new matrix condition number , thereby accelerating the convergence of the iterative method.

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

[0211]

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

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

[0214] ;

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

[0216] The constructed preconditioner M is:

[0217] ;

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

[0219] ;

[0220] The preconditioner M can be expressed as:

[0221] ;

[0222] The gradient is calculated as follows:

[0223] ;

[0224] 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 , ;

[0225] in, is the Jacobian matrix, and its elements are expressed as:

[0226] ;

[0227] Where, Expressed as A forward response, Indicates the model parameters, which is expressed as The observation data for Partial derivatives of model parameters

[0228] The updated approximate objective function is:

[0229] ;

[0230] In the formula is the data item weight matrix related to the measured data;

[0231] right about Taking the derivative and setting the reciprocal to zero gives us the linear system:

[0232] ;

[0233] The model update formula is:

[0234] ;

[0235] in: ;

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

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

[0238] ;

[0239] Where, represents the gradient vector, represents the Hessian matrix;

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

[0241] ;

[0242] Where, represents the residual after preconditioning, represents the initial residual, Indicates the initial search direction,

[0243] Iterative update step size Compute the quadratic function that is minimized for J :

[0244] ;

[0245] Along direction The minimum point satisfies ,get:

[0246] ;

[0247] The numerator is the inner product of the preconditioned residual, which measures the current residual energy; represents the residual of the k-th iteration, 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 The direction of maximum drop, Indicates the search direction of the kth iteration;

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

[0249] ;

[0250] ;

[0251] Step along the search direction , 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:

[0252] ;

[0253] renew Weighing the energy ratio of the new and old residuals and adjusting the search direction yields:

[0254] ;

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

[0256] ;

[0257] Where, represents the adjustment coefficient of the conjugate direction;

[0258] 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:

[0259] ;

[0260] Where, represents the initial residual The norm of , as the convergence criterion.

[0261] 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:

[0262] ;

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

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

[0265] ;

[0266] along with Changes, the curve shows the trade-off under different regularization strengths:

[0267] For small Region: The data fitting residual is small, the solution norm is large, the curve is approximately horizontal, and the residual changes slowly;

[0268] For big Region: The solution norm is small, the residual increases rapidly, the curve is approximately vertical, and the solution norm changes slowly;

[0269] For the inflection point region, the curvature at the inflection point is the largest, which corresponds to the optimal ;

[0270] The inflection point corresponds to the point with the maximum curvature. The calculation steps are:

[0271] Parametric curves: treating the L-curve as a function ,in

[0272] ;

[0273] Calculate curvature:

[0274] ;

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

[0276] The L-curve criterion is used to realize the regularization parameter adaptation, and the data fitting difference is obtained. and model complexity The trade-off curve selects the optimal , for the generated anisotropic weights, traverse different value, run the inversion to calculate the residual and solution norm, and finally identify the inflection point of the L-curve to ensure that the inversion model fits the data and maintains the structural characteristics. Estimated to be:

[0277]

[0278] Where, Indicates poor data fit, and the denominator represents the anisotropic characteristics of model complexity

[0279] ;

[0280] in is the cooling factor, represents the regularization parameter of the kth iteration, The above setting is to use strong regularization (large ) to ensure stability, and then gradually reduce it to improve resolution.

[0281] The inversion is defined as the RMS value between the model response and the measured data:

[0282] ;

[0283] Where, Represents the measured data, Represents the forward calculation data, For the The standard deviation of the data, The total number of data.

[0284] 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;

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

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

[0287] ;

[0288] The anisotropy strength is:

[0289] ;

[0290] The final x-direction weight is:

[0291] ;

[0292] Where, 、 and Represents the standard basis vector (x, y, z axis unit vector)

[0293] For the direction vector , the regularization term in one direction is decomposed into:

[0294] ;

[0295] Where, is the absolute value sum of the direction vector, is the model norm regularization parameter, , Represents the regularization parameter of the gradient term in the x, y, and z directions:

[0296] ;

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

[0298] ;

[0299] Where, is the regularization term, is the model vector.

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

[0301] Model instance calculation

[0302] This proposal proposes a seismic data-constrained 3D inversion algorithm based on the DC electrical method, which is implemented based on the maximum principal eigenvector of the structural tensor. To verify the effectiveness of this guidance method and the agreement between 3D Gauss-Newton inversion and 3D data-constrained inversion, a simple 3D subsurface model containing high-resistivity and low-resistivity anomalies was constructed and inverted.

[0303] Establishment of resistivity model: First, the present invention established the following test Figure 3 The three-dimensional resistivity model of a uniform half-space underground body with high and low resistivity anomalies is shown to test whether constrained inversion can accurately recover the target volume boundary. The surface terrain of the model is a basin 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-resistivity area with a resistivity of 10Ω·m (yellow). The resistivity of the underground background part is 100Ω·m, and the air resistivity is At the same time, the model takes into account the influence of terrain on the forward data and removes the air layer above the terrain.

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

[0305] First, a main survey line 1 with a length of 176 m from x = -88 m to 88 m was arranged at y = 0 to survey the abnormal display of the two anomalies in the xoz plane.

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

[0307] 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-resistivity layered anomalies.

[0308] For the three survey lines above, a dipole-dipole setup was selected, using a 1A current source, an electrode spacing of 12m, and six receivers. To eliminate the effects of topography, each electrode was redefined to be positioned above the terrain. Three-dimensional finite volume numerical simulations were performed using a modified version of the open-source software package SimPEG to calculate the received signal data for the dipole-dipole setup under a 1A current source. These data were then added with random Gaussian noise with a mean of zero and a standard deviation of 10% and used as synthetic data for inversion.

[0309] 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 )、 ( Figure 5) 2D cross-sections. Both the 3D pseudo-cross-sections and the 2D forward cross-sections clearly reflect high and low resistivity anomalies, while also significantly impacting the identification of multiple targets.

[0310] Build velocity model

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

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

[0313] 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 region refinement in grid setting to establish the three-dimensional inversion grid.

[0314] Define the basic unit size as 4 as the benchmark 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 units; the middle refined area: the unit length is halved to 2, containing 30 units, and the total length is 60 units. The y-direction maintains the same structure as the x-direction to form a symmetrical partition with the same total length and number of units. The z-direction is divided into two areas: the lower coarse grid area: the grid length is 4, containing 10 units, and the total length is 40; the upper refined area: the length is 2, containing 20 units, and the total length is 40 units. After the above refined and coarse parts of the grid are generated, a total of Taking the influence of the terrain into consideration, the remaining valid calculation grids are 100416 after removing the part above the terrain.

[0315] Calculating constraint weights for seismic data

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

[0317] The structure tensor field of the three-dimensional image reveals the geological structure characteristics of different regions through the shape of the characteristic ellipsoid and the eigenvalues of the structure tensor in the three-dimensional model. When the model is in a homogeneous medium area, the three eigenvalues 、 、 Both approach zero, at which point the characteristic ellipsoid degenerates into a tiny point-like structure near the origin, reflecting the lack of directional change 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 It is significantly larger than the other two eigenvalues, forming a cigar-shaped ellipsoid extending along the direction of the main eigenvector, with its major axis perpendicular to the geological interface. The surface structure area is characterized by the first two eigenvalues. 、 Close to and much larger than The ellipsoid is compressed into a disk shape with its normal perpendicular to the bedding plane. At this point, the weight generation module strengthens the constraint strength in both directions within the plane. The quantitative analysis of the characteristic ellipsoid begins with the calculation of the gradient field. Spatial differential information is obtained through Gaussian derivative filtering, and then a structural tensor matrix containing six independent components is constructed. The eigenvalues and eigenvectors extracted after eigendecomposition fully describe the local geometric characteristics.

[0318] Figure 8 、 Figure 9 and Figure 10 The model and maximum principal eigenvalue of the guide image at x=0, y=0 and depth z=-20m are In the cross-sectional comparison diagram, 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 closer to the surface.

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

[0320] As can be seen in the figure, the anisotropic weights increase from the center to the edge of the target, ultimately reaching 1 at the background, effectively reflecting both the center and edges of the target. Regarding the anisotropic weights, the weight coefficients defined above are assigned as a=2 and b=0.5. To distinguish the boundaries between the two targets, the anisotropic weights at the center of the anomaly are assigned a smaller weight, while the anisotropic weights at the target boundary effectively reflect the target boundary recognition.

[0321] 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 as shown in the figure below. Figure 15As shown in the figure, it can be seen that the root mean square errors of both inversion methods fluctuate and drop 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.

[0322] The comparison between the Gauss-Newton inversion and image-guided inversion results and the correct model on the y=0 interface is shown in Figure 2. Figure 16 、 Figure 17 As shown in the figure, compared to the single-method inversion results generated by traditional Gauss-Newton inversion, the 3D image-guided inversion more clearly depicts the two target volumes, demonstrating reliability in identifying multiple targets. Furthermore, judging by the internal high-resistance / low-resistance morphology and extension of the two target volumes, the center and boundaries of the high-resistance volume are accurately depicted, improving vertical resolution. The center, upper and lower limits of the low-resistance volume are also more accurately depicted than with traditional inversion, and are generally consistent with the true model.

[0323] Example 1

[0324] The methods 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:

[0325] Theoretical application:

[0326] 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 the collaborative inversion of electrical and seismic data is established, promoting the development of multi-method data fusion theory.

[0327] 2. 3D inversion algorithm optimization: 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 electromagnetics and gravity, improving the computational efficiency and stability of complex models.

[0328] 3. Structural Constraint Inversion Theory: Through the directional coupling of seismic image feature extraction and electrical models, it provides a new approach to solving the multi-solution problem of single-method inversion. It is particularly suitable for the refined characterization of structural features such as layered media and fault boundaries.

[0329] Engineering Applications: Compatible with SimPEG standard input formats (Survey, Mesh), it supports seamless integration with existing electrical and seismic processing workflows (such as ResIPy and SeisPy), reducing user migration costs. It is also compatible with the terrain correction module, enabling joint modeling of seismic and electrical data on undulating surfaces through the mesh deformation interface (Mesh.TensorMesh), addressing the limitations of the default flat terrain assumption.

[0330] 1. Mineral Resource Exploration: This technology is suitable for precise 3D electrical structure imaging of metal ore bodies and oil and gas reservoirs. It can be combined with seismic velocity model constraints to enhance ore body boundary identification and reduce drilling verification costs. For example, in skarn deposits, combined analysis of resistivity anomalies and seismic velocity mutations can effectively delineate mineralized alteration zones.

[0331] 2. Hydrogeological survey: Used for the detailed characterization of the spatial distribution of aquifers and the boundaries of groundwater pollution plumes. Combined with seismic reflection interface constraints, it can distinguish the contact relationship between low-resistance aquifers and high-resistance bedrock, thereby improving the accuracy of water resource assessment.

[0332] 3. Engineering geological survey:

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

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

[0335] 4. Environment and Hazard Geology:

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

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

[0338] 5. Urban underground space development: 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.

[0339] Therefore, the present invention adopts the above-mentioned three-dimensional inversion method of DC electrical 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 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.

[0340] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not 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 solutions of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A three-dimensional inversion method constrained by direct current (DC) method based on seismic data, characterized by: Including the processing of seismic data and the inversion processing of resistivity data; The processing of seismic data includes the following steps: Step 1: Process the 3D seismic velocity data or 3D seismic images associated with the electrical exploration data, remove the air layer, and calculate the structural tensor to obtain structural tensor data. This data is then accelerated through GPU and CPU hybrid computing to calculate the structural tensor and perform eigendecomposition. Step 2: Using the multi-scale Gaussian filtering method, weighted data is constructed by judging the correlation, and the mask processing method is used to solve the problem of blurred terrain boundaries and invalid area interference in the geological model, providing input data for anisotropic inversion; Step 3: Define the anisotropy coefficient by referring to the structural characteristics of the correlation between the eigenvalue ratios corresponding to the main feature direction and the normal feature direction, and store and visualize the anisotropy weights 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, and obtain a set of discrete equations by discretizing the integral equations. The required apparent resistivity is obtained by solving the discrete equations. Step S2: Obtain resistivity inversion under seismic data constraints using the Gauss-Newton inversion framework. The linear system is solved using the conjugate gradient method (PCG). Inversion is performed using an adaptive regularization parameter under the L-curve criterion. The root mean square error (RMS) 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. The method for three-dimensional inversion based on 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. The three-dimensional discrete gradient operator is defined as: ； Remove high-frequency noise using Gaussian smoothing: ； Gaussian kernel 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: ； Where, 、 and are the gradients in each direction, is the structure tensor matrix corresponding to node (i, j, k), is the convolution operator, is a two-dimensional Gaussian difference operator, To control the Gaussian function; The expression of the Gaussian difference operator is: ； ； Structure tensor matrix For one The real pair matrix of is decomposed into three eigenvalues ( 、 、 ) and three pairwise orthogonal eigenvectors ( 、 、 ), where the maximum eigenvalue corresponds to The eigenvector of represents the direction pointing to the strongest local gradient energy, perpendicular to the structure boundary; In a continuous medium, perpendicular to the stratum interface; in sharp edges, Perpendicular to the fault plane; second eigenvector In and The orthogonal plane represents the secondary change direction of the local structure, and the third eigenvector points in the direction of minimum local variation, usually parallel to the main extension direction of the structure, , indicating strong anisotropy; , indicating isotropy.

3. The three-dimensional inversion method of DC electrical method constrained by seismic data according to claim 2, characterized in that: In step 1, GPU and CPU hybrid computing accelerates structural tensor calculation and eigendecomposition based on the graphics card hardware environment. The calculation process is as follows: 1.1: Migrate velocity model data from host memory to graphics card memory through the asynchronous transmission channel of the CuPy library, and use the unified memory architecture of CUDA 8.0 to achieve zero-copy data exchange between the CPU and GPU; 1.2: The gradient calculation stage uses a parallel implementation of the 3D central difference method. CUDA thread blocks are responsible for differential operations in a 32×32×32 voxel area. Shared memory is used to cache local data to reduce global memory access latency. 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 (SMs) for parallel execution. The element-level operations corresponding to the components are vectorized through the GPU broadcast mechanism. Combined with the separable convolution feature of the Gaussian filter, the three-dimensional convolution is decomposed into a cascade of three one-dimensional convolution operations.

4. The method of three-dimensional inversion constrained by direct current method according to claim 3, characterized in that: In step 2, by defining Gaussian filters for horizontal and vertical structures and , calculate similarity , where the diffusion tensor Determined by the eigenvector and weight: ； Similarity Defined as: ； The calculation is divided into the following three steps, using a horizontal priority filtering strategy to enhance the response of horizontally continuous structures: (1) Perform transverse Gaussian filtering along the xy plane. The anisotropic Gaussian kernel is shown as follows: ； Here , indicating smoothing only in the xy plane; (2) Perform longitudinal Gaussian filtering along the z direction: ； Here , indicating filtering along the z direction; (3) Final smoothing: Where, Represents the multiplication of elements at corresponding positions of the matrix; is three-dimensional image data with dimensions (x, y, z); for( , ), Therefore A three-dimensional Gaussian filter with parameters , which filters only in the x and y directions and keeps the original value in the z direction; for( , ) three-dimensional isotropic Gaussian filter for smoothing in all directions; The calculation of is divided into the following three steps, which enhance the response of vertical discontinuities by longitudinal-first filtering: ① Perform longitudinal Gaussian filtering along the yz plane: ； Here , indicating smoothing only on the yz plane; ② Perform horizontal Gaussian filtering along the x direction: ； Here , indicating filtering along the x direction; ③Final smoothing: Where, Represents the multiplication of elements at corresponding positions in the matrix.

5. The three-dimensional inversion method of DC method constrained by seismic data according to claim 4, characterized in that: In the mask processing process of step 2, the essence of the mask is a binary matrix, which is used to mark the valid area of the data space. The mathematical form is expressed as: ； In the formula Represented as a binary matrix with spatial position correlation, data is selectively retained or suppressed through element-wise multiplication: ； Where, is the original data matrix, Expressed as The binary matrix generated by the function is the same as the original data matrix The same dimensions, Represents the mask matrix The result of selectively retaining the original data D; The calculation of the mask introduced in the gradient calculation is: ； A dynamic mask is introduced in the structural similarity calculation to keep the boundary clear by applying the mask multiple times. The calculation is expressed as: Initial data mask: ； Gaussian filter mask: ； Secondary mask: 。 6. The three-dimensional inversion method of DC electrical 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 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 similarity is determined When the value is less than the boundary judgment condition, it is an edge area. In the edge area, the smoothing intensity is reduced to retain the sharp edge. The value of the boundary judgment condition edge_threshold is set to a constant. When similarity is determined The value greater than the coherence judgment condition is the structural coherence area, in which the smoothing along the main direction is enhanced and the noise is suppressed; the value of the coherence judgment condition coherence_threshold is set to a constant; 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: ； in is the weight calculation method, are the eigenvalue and eigenvector of the structure tensor respectively, and a and b are the assigned weight coefficients.

7. The 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 recorded as , the source term of the control equation discretized in the control volume is: ； Where u represents the potential, σ represents the dielectric conductivity, The source term represents the current injection or reception per unit volume; Find the volume integral of the above equation within the control volume: ； Applying Gauss's theorem on the left side yields: ； Within the control volume Indicates volume, is the surface area, and the area differentiation of the right equation is simplified to: ； because , = , , further simplified to: ； , , are the lengths of the control volume in the x, y, and z directions, respectively, where Expressed as: ； is the harmonic mean of the conductivity of the two control volumes in the x direction: ； In the above formula: ； therefore: ； ； ； The source term is: 。 8. The 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 , ； in, is the Jacobian matrix, and its elements are expressed as: ； Where, Expressed as A forward response, Indicates the model parameters, For the The observation data for Partial derivatives of model parameters; The updated approximate objective function is: ； In the formula is the data item weight matrix with the measured data; right about Taking the derivative and setting the reciprocal to zero gives us the linear system: ； The model update formula is: ； in: ; The obtained linear system is solved by the conjugate gradient method (CG) combined with the preconditioning method - 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: ； Where, represents the gradient vector, 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: ； Where, represents the residual after preconditioning, represents the initial residual, Indicates the initial search direction, Iterative update step size Compute the quadratic function that is minimized for J : ； Along direction The minimum point satisfies , get the step length : ； The numerator is the inner product of the preconditioned residual, which measures the current residual energy; represents the residual of the kth iteration; 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. Indicates the search direction of the kth iteration; The subsequent model update and residual update are: ； ； Step along the search direction , 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: ； renew Weighing the energy ratio of the new and old residuals and adjusting the search direction yields: ； Combine the current preconditioned residual with the previous step direction to generate the conjugate direction: ； Where, 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, represents the initial residual The norm of , as the convergence criterion.

9. The 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 observation data; L is the regularization matrix, λ is the regularization parameter; The L-curve is defined as a parameterized curve: ； along with Changes, the curve shows the trade-off under different regularization strengths: For small Region: The data fitting residual is small, the solution norm is large, the curve is approximately horizontal, and the residual changes slowly; For big 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 is the largest, which corresponds to the optimal ; The inflection point corresponds to the point with the maximum curvature. The calculation steps are: Parametric curves: treating the L-curve as a function ,in ； Calculate 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 is obtained. and model complexity The trade-off curve selects the optimal β, and for the generated anisotropic weights, traverses different value, run the inversion to calculate the residual and solution norm, and finally identify the inflection point of the L-curve to ensure that the inversion model fits the data while maintaining the structural characteristics. The initial β estimate is: Where, Indicates poor data fit, and the denominator represents the anisotropic characteristics of model complexity ； in is the cooling factor, represents the regularization parameter of the kth iteration, β For cooling strategy; The inversion is defined as the RMS value between the model response and the measured data: ； Where, Represents the measured data, Represents the forward calculation data, For the The standard deviation of the data, The total number of data.

10. The method for three-dimensional inversion constrained by direct current method according to claim 9, wherein in step S3, generating a regularization term containing the weight of the seismic data constraint 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: ； Where, 、 and Represents the standard basis vectors, x, y, z axis unit vectors; For the direction vector , the regularization term in one direction is decomposed into: ； Where, is the absolute value sum of the direction vector, is the model norm regularization parameter, , Represents the regularization parameter of the gradient term in the x, y, and z directions; ； Anisotropic regularization is the sum of three orthogonal regularization terms: ； Where, is the regularization term, is the model vector.

Images