Search

US-20260127338-A1 - 3D INVERSION METHOD OF DC METHOD USING SEISMIC DATA CONSTRAINTS

US20260127338A1US 20260127338 A1US20260127338 A1US 20260127338A1US-20260127338-A1

Abstract

A 3D inversion method of DC method using seismic data constraints is disclosed, comprising processing of seismic data and inversion processing of resistivity data; by utilizing seismic wave interval velocity information to constrain the high-density electrical inversion, which can play a complementary role to a certain extent, it is beneficial to supplement the resolution of resistivity inversion with 3D seismic image information, reflect a boundary of subsurface electrical anomalous body, and further characterizes the boundary characteristics of anomalous body through comprehensive geological-geophysical analysis. Seismic data constraints are applied to 3D resistivity inversion, performing numerical simulation and inversion of electrical data under undulating terrain using the finite volume method and Gauss-Newton inversion method. Furthermore, the seismic anisotropic velocity tensor constraint is introduced for inversion to mitigate the non-uniqueness inherent in single-method inversion, yielding a more accurate 3D subsurface electrical structure and improving the resolution of longitudinal boundaries.

Inventors

  • Zhenwei GUO
  • Yanyi WANG
  • Zhuo Liu
  • Yanjun Chen
  • Jianxin Liu
  • Bochen Wang
  • Chengxi WANG

Assignees

  • CENTRAL SOUTH UNIVERSITY

Dates

Publication Date
20260507
Application Date
20260106
Priority Date
20250319

Claims (10)

  1. 1 . A three-dimensional (3D) inversion method of a direct current (DC) method using seismic data constraints, comprising a processing of seismic data and an inversion of resistivity data; wherein the processing of the seismic data comprises the following steps: step 1: processing 3D seismic velocity data or a 3D seismic image related to electric exploration data, removing an air layer part, performing a calculation of a structure tensor to obtain structure tensor data, and accelerating a calculation of the structure tensor and eigendecomposition through a hybrid calculation of a graphics processing unit (GPU) and a central processing unit (CPU); step 2: adopting a multi-scale Gaussian filtering method, constructing weight data through correlation determinations, solving problems of fuzzy terrain boundary and invalid region interference of a geological model by a mask processing method, thereby providing input data for anisotropic inversion; and step 3: calculating and defining an anisotropy coefficient with reference to a structural feature of the correlation by an eigenvalue ratio corresponding to a main feature direction and a normal feature direction, and simultaneously outputting a storage and a visualization of anisotropy weight in a modular form; wherein the processing of the inversion of resistivity data comprises the following steps: S1: performing a calculation by adopting a finite volume method, discretizing a calculation region into non-overlapping control volumes or control voxels, applying conservation equations to integrate the control volumes, obtaining a set of discrete equations by discretizing the integral equations, and obtaining required apparent resistivity by solving the discrete equations; S2: obtaining resistivity inversion under seismic data constraint by updating a framework of Gauss-Newton inversion, solving a linear system by a preconditioned conjugate gradient (PCG) method, performing an inversion with an adaptive regularization parameter guided by a L-curve criterion, and evaluating an inversion effect by using a root mean square error; S3: receiving the anisotropic weight stored from the seismic data calculation in step 3 through a modular interface, and forming a new regularization term to constrain the inversion; and S4: performing a model generation and calculation, and visualizing results of resistivity inversion under the constraint of seismic data.
  2. 2 . The 3D inversion method of the DC method using seismic data constraints according to claim 1 , wherein in step 1, a calculation method of structure tensor data is as follows: calculating a gradient field by using a central difference combined with a Gaussian derivative kernel, and a 3D discrete gradient operator is defined as: { I x [ i , j , k ] = I [ i + 1 , j , k ] - I [ i - 1 , j , k ] 2 ⁢ Δ ⁢ x I y [ i , j , k ] = I [ i , j + 1 , k ] - I [ i , j - 1 , k ] 2 ⁢ Δ ⁢ y I z [ i , j , k ] = I [ i , j , k + 1 ] - I [ i , j , k - 1 ] 2 ⁢ Δ ⁢ z ; eliminating a high-frequency noise by Gaussian smoothing: = G σ * I x , = G z * I z , = G σ * I z ; a discrete form of Gaussian kernel G, is: G σ [ n ] = 1 ( 2 ⁢ π ) 1 / 2 ⁢ σ ⁢ e n 2 2 ⁢ σ 2 ; wherein the structure tensor is a second-order tensor matrix derived from a gradient, and a mathematical form of the 3D structure tensor is: S σ = G σ * [ I x 2 I x ⁢ I y I x ⁢ I z I x ⁢ I y I y 2 I y ⁢ I z I x ⁢ I z I y ⁢ I z I z 2 ] = [ S xx S xy S xz S xy S yy S yz S xz S yz S zz ] ; where I x , I y and I z are gradients in each direction, S σ is a structure tensor matrix corresponding to a node (i,j,k), *is a convolution operator, G σ is a two-dimensional Gaussian difference operator, and σ controls a Gaussian function; wherein an expression of the Gaussian difference operator is: G σ ( x , y , z ) = 1 ( 2 ⁢ π ⁢ σ ) 3 / 2 ⁢ G x ⁢ G y ⁢ G z ⁢ e - ( x 2 2 ⁢ σ x 2 + y 2 2 ⁢ σ y 2 + z 2 2 ⁢ σ z 2 ) ; S σ = λ 1 ⁢ ν 1 ⁢ ν 1 T + λ 2 ⁢ ν 2 ⁢ ν 2 T + λ 3 ⁢ ν 3 ⁢ ν 3 T ⁢ ( λ 1 > λ 2 > λ 3 ) ; the structure tensor matrix S σ is a 3×3 real matrix, which is decomposed into three eigenvalues (λ 1 ,λ 2 ,λ 3 ) and three pairwise orthogonal eigenvectors (v 1 ,v 2 ,v 3 ), wherein the eigenvector v 1 corresponding to a largest eigenvalue denotes a direction pointing to a maximum local gradient energy, which is perpendicular to a structural boundary; wherein in a continuous medium, v 1 is perpendicular to a formation interface; in a sharp edge, v 1 is perpendicular to a fault surface; the second eigenvector v 2 is in a plane orthogonal to v 1 , which denotes a secondary change direction of a local structure, the third eigenvector v 3 points to a direction with a minimum local change, which is typically parallel to a main extension direction of the structure, λ 1 >>λ 3 , denoteing strong anisotropy; λ 1 ≈λ 3 , denoteing isotropy.
  3. 3 . The 3D inversion method of the DC method using seismic data constraints according to claim 2 , wherein in step 1, the accelerating the calculation of the structure tensor and eigendecomposition through the hybrid calculation of the GPU and CPU is performed on a hardware environment with a graphics card, and a calculation process is as follows: 1.1: migrating velocity model data from host memory to video memory of the graphics card through an asynchronous transmission channel of a CuPy library, and achieving a zero-copy data interaction between CPU and GPU by using a unified memory architecture of CUDA8.0; 1.2: utilizing a parallelized 3D central difference method for a gradient computation, with compute unified device architecture (CUDA) thread blocks handling differential operations on 32×32×32 voxel regions, and employing shared memory caching of local data to reduce global memory access latency; and 1.3: during a construction process of the structure tensor, mapping calculations of six independent components (Txx, Txy, Txz, Tyy, Tyz, and Tzz) to different streaming multiprocessors (SM) for a parallel execution, vectorizing element-wise operations corresponding to the components through a broadcast mechanism of GPU, and utilizing separable convolution characteristics of Gaussian filter to decompose the 3D convolution into a cascade of three one-dimensional (1D) convolutions.
  4. 4 . The 3D inversion method of the DC method using seismic data constraints according to claim 3 , wherein in step 2, a similarity S is calculated by defining the Gaussian filters S 1 and S 2 of transverse structures and longitudinal structures, wherein a diffusion tensor D is determined by a feature vector and weight: ∂ m ∂ t = ∇ · ( D ⁢ ∇ m ) ; wherein the similarity S is defined as: S = S 2 S 1 + ϵ ; wherein a calculation of S 1 is divided into the following three steps, wherein a response of a horizontal continuous structure is enhanced by a transverse-priority filtering strategy: (1) performing a transverse Gaussian filtering along an x-y plane, wherein the anisotropic Gaussian kernel is shown as follows: I temp ⁢ 1 = G σ 1 * I σ 1 = ( 2 , 2 , 0 ) ; here σ x =σ y =2, σ z =0, which denotes only smoothing in the x-y plane; (2) performing a longitudinal Gaussian filtering along a z direction: I temp ⁢ 2 = G 2 * I temp ⁢ 1 σ 2 = ( 0 , 0 , 8 ) ; here σ x =σ y =0, σ z =8, which denotes filtering along the z direction; (3) final smoothing: S 1 = G σ 3 * ( I temp ⁢ 1 ⊙ I temp ⁢ 2 ) σ 3 = 2 ; where ⊙ denotes a multiplication of corresponding position elements of the matrix; wherein I is 3D image data, with a dimension of (x, y, z); σ 1 is a Gaussian kernel parameter of (σ x =2, σ y =2,), G σ 1 is a 3D Gaussian filter with σ 1 as the parameter, only filtering in the x and y directions, and maintaining an original value in the z direction; and G σ 3 is a 3D isotropic Gaussian filter for (σ x =2, σ y =2,), smoothing in all directions; wherein a calculation of S 2 is divided into the following three steps, wherein a response of the vertical continuous structure is enhanced by a longitudinal-priority filtering strategy: {circle around (1)} performing a longitudinal Gaussian filtering along a y=z plane: I temp ⁢ 3 = G σ 4 * I σ 4 = ( 0 , 8 , 8 ) ; here σ y =σ z =8, σ x =0, which denotes only smoothing in the y=z plane; (2) performing a transverse Gaussian filtering along an x direction: I temp ⁢ 4 = G 5 * I temp ⁢ 3 σ 5 = ( 2 , 0 , 0 ) ; here σ y =σ z =0, σ x =2, which denotes filtering along the x direction; (3) final smoothing: S 2 = G σ 6 * ( I temp ⁢ 3 ⊙ I temp ⁢ 4 ) σ 6 = 8 ; where ⊙ denotes a multiplication of corresponding position elements of the matrix.
  5. 5 . The 3D inversion method of the DC method using seismic data constraints according to claim 4 , wherein in the mask processing process of step 2, an essence of the mask is a binary matrix, configured to mark valid regions of data space, mathematically expressed as: M ⁡ ( x ) = { 1 , x ∈ valid ⁢ region 0 , x ∈ invalid ⁢ region ; where M(x) denotes a spatial position-dependent binary matrix, configured to selectively retain or suppress data by performing an element-wise multiplication: D masked = D ⊙ M ; where D is an original data matrix, M is denoted as a binary matrix generated by the function M(x), which is the same as a dimension of the original data matrix D, and D masked denotes a result of the selective retention of the original data D through the mask matrix M; wherein a calculation of introducing a mask in the gradient calculation is: ∇ I masked = ∇ I ⊙ M ; wherein a dynamic mask is introduced in a structural similarity calculation, and the boundary is kept clear through multiple applications of the mask, and a calculation is: an initial data mask: D 0 = D ⊙ M ; a gaussian filter mask: D 1 = G σ ( D 0 ) ; a secondary mask: D 2 masked = D 1 ⊙ M .
  6. 6 . The 3D inversion method of the DC method using seismic data constraints according to claim 5 , wherein in step 3, a process of anisotropy and weight calculation is as follows: calculating and defining an anisotropy coefficient by an eigenvalue ratio corresponding to a main feature direction and a normal feature direction: λ r = λ 1 λ 3 + ϵ ; α = 1 - e - k ⁢ λ r ; where k is a positive real number, configured to adjust an influence of an eigenvalue ratio λ r on the anisotropy coefficient, k→0: wherein the anisotropy effect almost disappears, and the model degenerates into isotropic smoothing, wherein a larger k results in more sensitive regional response; a smaller k weakens anisotropic differences, and a value of k is [0.5,2.0]. obtaining weight values of the following three constraints based on the determination of the calculated similarity; when the similarity S is less than a value of a boundary determination condition, it is determined as an edge region, and a smoothing intensity is reduced in the edge region, and a sharp boundary is retained; wherein a value edge_threshold of the boundary determination condition is set to a constant; when the similarity S is greater than a value of a coherence determination condition, it is determined as a structurally coherent region, and smoothing along the main direction is enhanced in a coherent region to suppress noise; wherein a value coherence_threshold of the coherence determination condition is set to a constant; other parts are set to 1: ω = { a ⁡ ( 1 + α ) S > coherence_threshold b ⁡ ( 1 - α ) S < edge_threshold 1 otherwise ; wherein input data of weight calculation has been preprocessed by the mask, and anisotropic weight after mask processing is expressed as: W ⁡ ( x ) = f ⁡ ( ∇ I masked , λ i , v i ) · M ⁡ ( x ) ; where ƒ is a weight calculation method, λ i and v i are an eigenvalue and an eigenvector of the structure tensor, respectively, and a and b are given weight coefficients.
  7. 7 . The 3D inversion method of DC method using seismic data constraints according to claim 6 , wherein in S1, according to a theory of finite volume method, a control volume of a unit where a node (i, j, k) is located is recorded as V i,j,k , and a discrete source term of a control equation in the control volume is: ∇ · ( σ ⁢ ∇ u ) = - s ; where u denotes a potential, σ denotes a conductivity of a medium, and −s is defined as a source term denoting current injection or reception per unit volume; performing a volume integral of the above equation in the control volume: 1 ❘ "\[LeftBracketingBar]" V i , j , k ❘ "\[RightBracketingBar]" ⁢ ∫ V i , j , k ∇ · ( σ ⁢ ∇ u ) ⁢ dV = - 1 ❘ "\[LeftBracketingBar]" V i , j , k ❘ "\[RightBracketingBar]" ⁢ ∫ V i , j , k s ⁢ dV ; applying the Gauss theorem to the left side expression to obtain: 1 ❘ "\[LeftBracketingBar]" V i , j , k ❘ "\[RightBracketingBar]" ⁢ ∫ V i , j , k ∇ · ( σ ⁢ ∇ u ) ⁢ dV = 1 ❘ "\[LeftBracketingBar]" V i , j , k ❘ "\[RightBracketingBar]" ⁢ ∫ ∂ V i , j , k ( σ ⁢ ∇ u ) · n ⁢ dS ; wherein in the control volume, V i,j,k denotes a volume, ∂V i,j,k is a surface area, and performing surface integration and simplification on the right side expression to obtain: ∫ ∂ V i , j , k ( σ ⁢ ∇ u ) · n ⁢ dS = ∫ S i + 1 2 σ ⁢ ∇ u ⁢ dS - ∫ S i - 1 2 σ ⁢ ∇ u ⁢ dS + ∫ S j + 1 2 σ ⁢ ∇ u ⁢ dS - ∫ S j - 1 2 σ ⁢ ∇ u ⁢ dS + ∫ S k + 1 2 σ ⁢ ∇ u ⁢ dS - ∫ S k - 1 2 σ ⁢ ∇ u ⁢ dS ; since ⁢ S i + 1 2 = S i - 1 2 , S j + 1 2 = S j - 1 2 , and S k + 1 2 = S k - 1 2 , further ⁢ simplifying ⁢ to ⁢ obtain : 1 ❘ "\[LeftBracketingBar]" V i , j , k ❘ "\[RightBracketingBar]" ⁢ ∫ ∂ V i , j , k ( σ ⁢ ∇ u ) · n ⁢ dS = ( σ ⁢ ∂ u ∂ x ) i + 1 2 , j , k - ( σ ⁢ ∂ u ∂ x ) i - 1 2 , j , k h i x + ( σ ⁢ ∂ u ∂ y ) i , j + 1 2 , k - ( σ ⁢ ∂ u ∂ y ) i , j - 1 2 , k h i y + ( σ ⁢ ∂ u ∂ z ) i , j , k + 1 2 - ( σ ⁢ ∂ u ∂ z ) i , j , k - 1 2 h i z ; h i x , h i y , and ⁢ h i z are lengths of the control volume in the x, y, and z directions, respectively, wherein ( σ ⁢ ∂ u ∂ x ) i + 1 2 , j , k is expressed as: ( σ ⁢ ∂ u ∂ x ) i + 1 2 , j , k = σ i + 1 2 , j , k ( ∂ u ∂ x ) i + 1 2 , j , k ; σ i + 1 2 , j , k is a harmonic mean of the two control volumes of the conductivity in the x direction: σ i + 1 2 , j , k = h i + 1 2 x ( ∫ x i x i + 1 σ - 1 ( x , y , z ) ⁢ dx ) - 1 ; in the above equation: h i + 1 2 x = x i + 1 - x i = ( h i + 1 x + h i x ) / 2 ; thus: ( σ ⁢ ∂ u ∂ x ) i + 1 2 , j , k = h i + 1 2 x ( h i x 2 ⁢ σ i , j , k + h i + 1 x 2 ⁢ σ i + 1 , j , k ) - 1 ⁢ ( u i + 1 , j , k - u i , j , k h i + 1 2 x ) ; ( σ ⁢ ∂ u ∂ x ) i - 1 2 , j , k = h i + 1 2 x ( h i x 2 ⁢ σ i , j , k + h i - 1 x 2 ⁢ σ i - 1 , j , k ) - 1 ⁢ ( u i , j , k - u i - 1 , j , k h i - 1 2 x ) ; 1 ❘ "\[LeftBracketingBar]" V i , j , k ❘ "\[RightBracketingBar]" ⁢ ∫ ∂ V i , j , k ( σ ⁢ ∇ u ) · n ⁢ dS = σ i + 1 2 , j , k ( u i + 1 , j , k - u i , j , k h i + 1 2 x ) - σ i - 1 2 , j , k ( u i , j , k - u i - 1 , j , k h i - 1 2 x ) h i x + σ i , j + 1 2 , k ⁢ ( u i , j + 1 , k - u i , j , k h j + 1 2 y ) - σ i , j - 1 2 , k ( u i , j , k - u i , j - 1 , k h j - 1 2 y ) h i y + σ i , j , k + 1 2 ( u i , j , k + 1 - u i , j , k h k + 1 2 z ) - σ i , j , k - 1 2 ( u i , j , k - u i , j , k - 1 h k - 1 2 z ) h i z = - 1 ❘ "\[LeftBracketingBar]" V i , j , k ❘ "\[RightBracketingBar]" ⁢ ∫ V i , j , k s ⁢ dV ; wherein the source term is: 1 ❘ "\[LeftBracketingBar]" V i , j , k ❘ "\[RightBracketingBar]" ⁢ ∫ V i , j , k s ⁢ dV = { I ( x a , y a , z a ) ∈ V i , j , k 0 ( x a , y a , z a ) ∉ V i , j , k .
  8. 8 . The 3D inversion method of DC method using seismic data constraints according to claim 7 , wherein in S 2 , a core of the Gauss-Newton inversion method is to perform a first-order Taylor expansion to a nonlinear forward response F(m), wherein a current model is m k , F ⁡ ( m k + 1 ) = F ⁡ ( m k ) + J k ( m k + 1 - m k ) ; where J k is a Jacobian matrix, the element is denoted by: J i ⁢ j = ∂ F i d ⁢ m j ; where F i is denoted as an i th forward response, m; denotes a j th model parameter, and J ij is a partial derivative of an i th observation data to the j th model parameter; wherein an updated approximate objective function is: Φ ⁡ ( m k + 1 ) = 1 2 ⁢  W d ( d - F ⁡ ( m k ) - W d ⁢ J k ( m k + 1 - m k ) )  2 + β 2 ⁢  W m ( m k + 1 - m r ⁢ e ⁢ f )  2 ; where W d is a data-weighting matrix for the measured data; computing a derivative of Φ(m k+1 ) with respect to m k+1 and setting the derivative to zero, to obtain a linear system: ( J k T ⁢ W d T ⁢ W d ⁢ J k + β ⁢ W m T ⁢ W m ) ⁢ ( m k + 1 - m k ) = J k T ⁢ W d T ⁢ W d ( d - F ⁡ ( m k ) ) - β ⁢ W m T ⁢ W m ( m k - m k + 1 ) ; wherein a model update equation is: m k + 1 = m k + αδ ⁢ m ; wherein : δ ⁢ m = ( J k T ⁢ W d T ⁢ W d ⁢ J k + β ⁢ W m T ⁢ W m ) - 1 ; the resulting linear system is solved by a preconditioned conjugate gradient method (PCG), which is a conjugate gradient (CG) method combined with a precondition: initializing parameters by setting an initial model update to zero, and an initial residual is equal to a gradient vector: δ ⁢ m 0 = 0 r 0 = g - H ⁢ δ ⁢ m 0 ; where g denotes a gradient vector, and H denotes a Hessian matrix; accelerating convergence by applying a preconditioned matrix M to residual; and setting an initial search direction as the preconditioned residual: z 0 = M - 1 ⁢ r 0 p 0 = z 0 ; where z 0 denotes a preconditioned residual, r 0 denotes an initial residual, and p 0 denotes an initial search direction, an iteratively updated step size α k is calculated by minimizing a quadratic function ƒ(δm): f ⁡ ( δ ⁢ m ) = 1 2 ⁢ ( δ ⁢ m T ⁢ H ⁢ δ ⁢ m - g T ⁢ δ ⁢ m ) ; wherein ∂ f ∂ α = 0 is satisfied at a minimum point along direction p k , and a step size α k is obtained: α k = r k T ⁢ z k p k T ⁢ Hp k ; wherein a molecule is an inner product of the preconditioned residual, configured to measure current residual energy; r k denotes a residual of a k th iteration; z k denotes a preconditioned residual, to adjust a distribution of the residual through the preconditioned matrix M, thus improving a condition number of the linear system, and accelerating the convergence; a denominator is an energy of the search direction under the Hessian metric, and p k denotes a search direction of the k th iteration; wherein subsequent model updates and residual updates are: δ ⁢ m k + 1 = δ ⁢ m k + α k ⁢ P k ; r k + 1 = r k - α k ⁢ P k ; stepping α k along the search direction, updating the model correction; applying preconditions and calculating new directions, processing a new residual by the preconditioner, and eliminating an ill-conditionality of the Hessian matrix: z k + 1 = M - 1 ⁢ r k + 1 ; updating β k to measure an old and new residual energy ratio, and adjusting the search direction to obtain: β k = r k + 1 T ⁢ z k + 1 r k T ⁢ z k ; generating a conjugate direction by combining the current preconditioned residual with the previous direction: p k + 1 = z k + 1 + β k ⁢ p k ; where β k denotes an adjustment coefficient of the conjugate direction; setting a termination condition to terminate the iteration when a residual norm is reduced to 1% of the initial residual to balance a computational efficiency and an accuracy, thus avoiding excessive iterations:  r k   g  < t ⁢ o ⁢ lCG ; where ∥g∥ denotes a norm of the initial residual r 0 =g, configured as a convergence benchmark.
  9. 9 . The 3D inversion method of the DC method using seismic data constraints according to claim 8 , wherein in S2, in the process of the adaptive regularization, a L-curve criterion is used to determine an optimal value of regularization parameter, and an objective function of regularization problem can be defined as: min x {  Ax - b  2 2 + λ 2 ⁢  Lx  2 2 } ; wherein A is a coefficient matrix, b is an observation data; L is a regularization matrix, and λ is a regularization parameter; wherein L-curve is defined as a parametric curve: ( log ⁢  Ax λ - b  2 , log ⁢  Lx λ  2 ) ; wherein as λ changes, the curve shows a trade-off under different regularization strengths: for a small λ region: a data fitting residual is small, a solution norm is large, the curve is approximately horizontal, and the residual change is slow; for a large λ region: the solution norm is small, the residual increases rapidly, the curve is approximately vertical, and the solution norm changes slowly; for an inflection point region, a curvature is largest at an inflection point that balances data fitting with solution smoothness, corresponding to an optimal λ; the inflection point corresponds to a point with the largest curvature, wherein calculation steps are as follows: parameterizing curve: defining the L-curve as a function η(ρ), wherein ρ = log ⁢  Ax λ - b  ; η = log ⁢  Lx λ  calculating a curvature: κ ⁡ ( λ ) = ρ′η″ - ρ″η′ ( ( ρ′ ) 2 + ( η′ ) 2 ) 3 / 2 ; selecting a maximum curvature point: the corresponding λ is an optimal value; wherein a L-curve criterion is used to achieve the adaptive regularization parameter, an optimal β is selected by a trade-off curve between a data fitting difference Φ d and a model complexity Φ m , wherein for the generated anisotropic weights, the different λ values are traversed, and the inversion calculation residual and solution norm are executed, finally, the L-curve inflection point is identified to ensure that the inversion model not only fits the data but also maintains the structural characteristics, wherein the initial β is estimated to be: β 0 = Φ d ( ( ρ′ ) 2 + ( η′ ) 2 ) 3 / 2 ; where Φ d is a difference of data fitting, and a denominator denotes an anisotropic feature of the model complexity; β k + 1 = β k γ ; wherein γ is a cooling factor, β k denotes a regularization parameter of the k th iteration, and β is a cooling strategy; the inversion is defined by the root mean square (RMS) between the model response and the measured data: RMS = ∑ i = 1 N [ ( d i obs - d i fwd ) σ i ] 2 N ; where d i obs denotes measured i th data, d i fwd denotes i th data of a forward calculation, σ i is a measurement standard deviation of the i th data, and N is a total number of data.
  10. 10 . The 3D inversion method of the DC method using seismic data constraints according to claim 9 , wherein in step S3, the regularization term generation with seismic data constraint weight comprises the following steps: a projection of a main feature vector in each direction is denoted as: { d x = ❘ "\[LeftBracketingBar]" V 1 · e x ❘ "\[RightBracketingBar]" d y = ❘ "\[LeftBracketingBar]" V 1 · e y ❘ "\[RightBracketingBar]" d z = ❘ "\[LeftBracketingBar]" V 1 · e z ❘ "\[RightBracketingBar]" ; an anisotropic strength is: A = 1 - e - k ⁢ λ 1 λ 3 + ϵ ; a final x-direction weight is: w x = { d x = 2 ⁢ ( 1 + d x ) S > τ coh d x = 0.5 ( 1 - d x ) S < τ edge 1 else ; where e x , e y , and e z denote standard basis vectors, which are unit vectors in the x, y, and z axes; for a direction vector v=(a, b, c), a one-direction regularization term decomposition is denoted as: φ v ( m ) = α s ⁢  m  2 + a ⁢ α x ∑ ❘ "\[LeftBracketingBar]" v i ❘ "\[RightBracketingBar]" ⁢  W x ⁢ ∂ x m  2 + b ⁢ α y ∑ ❘ "\[LeftBracketingBar]" v i ❘ "\[RightBracketingBar]" ⁢  W y ⁢ ∂ y m  2 + c ⁢ α x ∑ ❘ "\[LeftBracketingBar]" v i ❘ "\[RightBracketingBar]" ⁢  W z ⁢ ∂ z m  2 ; where Σ|v i | is a sum of an absolute value of the direction vector, α s is a regularization parameter of a model norm, α x , α y , and α z denote regularization parameters of gradient terms in the x, y, z directions; ∑ ❘ "\[LeftBracketingBar]" v i ❘ "\[RightBracketingBar]" = ❘ "\[LeftBracketingBar]" a ❘ "\[RightBracketingBar]" + ❘ "\[LeftBracketingBar]" b ❘ "\[RightBracketingBar]" + ❘ "\[LeftBracketingBar]" c ❘ "\[RightBracketingBar]" ; anisotropic regularization is formed by summing three orthogonal direction regularization terms: φ m ( m ) = φ x ( m ) + φ y ( m ) + φ z ( m ) ; where φ m (m) is a regularization term, and m is a model vector.

Description

TECHNICAL FIELD The present disclosure relates to the field of geophysical exploration, in particular to a three-dimensional (3D) inversion method of a direct current (DC) method using seismic data constraints. BACKGROUND Electrical exploration is a classic geophysical exploration method. While high-density electrical method results typically produce two-dimensional apparent resistivity profiles. Nevertheless, for complex geological conditions, two-dimensional inversion is difficult to meet the exploration requirements, and 3D inversion is required to accurately grasp the geometric position and physical characteristics of the target. There are two key challenges with the high-density electrical method: firstly, the resistivity inversion algorithm needs to be improved, and the calculation of 3D inversion requires improved algorithm and hardware conditions; secondly, the longitudinal resolution of subsurface resistivity distributions derived from high-density electrical data inversion is insufficient, making it difficult to delineate anomaly boundaries, which makes it challenging to accurately characterize the anomaly boundaries. Joint inversion is one of the main methods to improve the quality of inversion at present. The joint inversion refers to the joint application of multiple geophysical observation data in geophysical inversion, and the same subsurface geological and geophysical model is jointly inverted through the relationship between rock physical properties and geometric parameters of geological bodies. 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 forward higher requirements for computer computing power and memory storage, thereby reducing the practical applicability of joint inversion. Therefore, it is necessary to provide a 3D inversion method of the DC method using seismic data constraints to solve the above problems. SUMMARY In order to solve the above problems, the present disclosure provides a 3D inversion method of a DC method using seismic data constraints. By utilizing seismic velocity information to constrain the high-density electrical inversion, which can play a complementary role to a certain extent, it is beneficial to supplement the longitudinal resolution of resistivity image information, invert a boundary of an subsurface electrical anomalous body, and further characterize the boundary characteristics of the anomalous body through comprehensive geological-geophysical information. The present disclosure studies 3D inversion of DC apparent resistivity constrained by seismic wave interval velocity, which applies seismic data constraints to 3D resistivity inversion, performs numerical simulation and inversion of electrical data under undulating terrain using the finite volume method and Gauss-Newton inversion method. Furthermore, the seismic anisotropic velocity tensor constraint is introduced for inversion to mitigate the non-uniqueness inherent in single-method inversion, yielding a more accurate 3D subsurface electrical structure and improving the resolution of longitudinal boundaries. Therefore, the present disclosure adopts the 3D inversion method of the DC method using seismic data constraints, and has the following beneficial effects: (1) In the present disclosure, the anisotropic characteristics of the formation interface are extracted through the seismic interval velocity structure tensor, the direction adaptive regularization constraint is constructed, thereby improving the vertical resolution of the resistivity inversion by more than 30% compared with the conventional method, and the anomalous body boundary positioning error≤5%.(2) In the present disclosure, structural constraint terms driven by seismic data are introduced, which breaks the equivalence dilemma of electrical inversion, reduces the model space solution set by 40%-60%, while yielding notably effective false anomaly suppression.(3) In the present disclosure, a finite volume method-Gauss-Newton mixed framework is employed, which enhances the stiffness matrix sparsity by 20% and reduces memory usage by 35%.(4) In the present disclosure, the GPU accelerates structure tensor calculation and eigendecomposition, increasing the operational velocity of the key module by 8-12 times compared to CPU execution.(5) In the present disclosure, the preconditioned conjugate gradient method (PCG) is combined with the adaptive regularization parameter strategy, which reduces the number of iterations by 50%-70%.(6) In the present disclosure, the visualization and storage of both seismic data structure gradient calculations and weight calculations are implemented in a modular manner, thereby ensuring the readability and stability of the technical methodology.(7) The present disclosure supports 3D model inversion at a million-unit level (grid scale ≥1,