Electromagnetic detection inversion is a typical ill-posed problem, and it is easy to cause multiple solutions in the inversion results. Ill-posedness is an inherent feature of the inversion itself. It is difficult to overcome the essential difficulty without additional information for the solution. The solution to this problem is An effective method is to study constrained inversion. In this paper, the current mainstream Gauss-Newton-conjugate gradient method (GN-CG) is used to directly impose constraints on the inversion objective function, and the value range of the dielectric resistivity is used as the prior information and point penalty function method outside the constraints. Compared with the conventional three-dimensional resistivity inversion objective function, the objective function of the inequality constraints is added, which can suppress the ambiguity of the inversion theoretically. The test results of various theoretical models show that the three-dimensional well-ground resistivity inversion algorithm based on inequality constraints can effectively improve the accuracy of the inversion results, and the method of imposing inequality constraints by the penalty function method is realistic and effective.
Abstract
The inversion of electromagnetic detection data is a typical ill-posed problem and is prone to cause a multiplicity of solutions of the inversion results. The ill-posedness is an inherent characteristic of inversion and is difficult to overcome without additional information. An effective way to solve this problem is constrained inversion. In this study, the Gauss-Newton - conjugate gradient (GN-CG) method was used to directly impose constraints on the inversion objective function. Specifically, the dielectric resistivity range was introduced into the inversion objective function as the prior information and constraints using the exterior penalty function method. Compared with the conventional three-dimensional resistivity inversion objective function, the objective function with inequality constraints can suppress the multiplicity of solutions in theory. As revealed by the testing results of various theoretical models, the three-dimensional borehole-to-surface resistivity inversion algorithm based on inequality constraints effectively improves the precision of inversion results, and the way of imposing inequality constraints using the penalty function method is feasible and effective.
Preface
Electromagnetic prospecting is based on the electrical differences between different rocks in the earth's crust (such as differences in electrical conductivity, magnetic permeability, dielectric properties, and electrochemical properties), and studies the temporal and spatial distribution of artificial or natural electromagnetic fields through observation. In order to achieve the purpose of searching for underground useful mineral resources, ascertaining underground geological structures and solving geological problems [ 1 ] . Electromagnetic prospecting methods are mainly divided into direct current method based on geometric sounding principle, magnetotelluric method (MT/AMT/CSAMT) based on frequency domain sounding principle and transient electromagnetic method (TEM) based on time domain sounding principle [ 2 ] etc. The DC resistivity method is one of the most classic methods of electrical exploration, and it has been widely and effectively applied in the exploration of metal minerals, non-metallic minerals, coal fields, oil and gas, and groundwater, and has been extended to hydrology, engineering, environment, archaeology, etc. and other detection fields that are closely related to national economic construction and human social life [ 3 ] . The DC resistivity method has a variety of flexible observation methods, such as the electrical profile method of single-pole, dipole and multi-pole devices, and the electrical sounding method. The electrical method, the power supply device is always placed in the deep part of the borehole to make it close to the object to be detected, so as to increase the current intensity or the abnormal response received [ 4-5 ] , the well-ground resistivity method is mainly used in metal mines Secondary resource exploration and prediction of reservoir boundaries [ 6 ] , etc., usually the acquisition device collects data on the grid along parallel lines with different electrode arrays, and uses 3D inversion algorithm for interpretation. With the development of computer and numerical calculation technology, the three-dimensional forward and inversion algorithm of resistivity is used in grid subdivision (structured [ 7 ]⇓ ⇓ ⇓ ⇓ - 12 ] , unstructured [ 13 ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ - 21 ] ), numerical methods (finite difference method [ 22 ⇓ ⇓ ⇓ ⇓ - 27 ] , finite element method [ 28 - 29 ] ) , the solution of the objective function (Gauss-Newton method (GN) [ 3 , 19 , 30 ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ - 37 ] , quasi-Newton method (QN) [ 38 ⇓ ⇓⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ - 54 ] , nonlinear conjugate gradient method (NLCG) [ 45 , 51 , 55 ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ ⇓ - 63 ] , etc. ) .
Since the grid of structured units cannot adapt to the 3D inversion of resistivity in complex terrain, the unstructured finite element method has achieved great success in 3D numerical simulation of resistivity in complex terrain. The advantages of allowing local encryption and being able to simulate complex geometric models have greatly improved the efficiency of 3D unstructured finite element solutions. In the case of the same calculation accuracy, compared with the structured regular grid, the calculation time and storage capacity of the unstructured grid can be reduced by about an order of magnitude [ 3 ] . In the three-dimensional inversion of resistivity, due to the large number of inversion parameters and the large amount of data, the Jacobian matrix (partial derivative matrix) is difficult to overcome the huge amount of calculation and storage requirements. Many scholars have studied the inversion method that can avoid calculating the Jacobian matrix. algorithm, Zhang et al. [ 26 ] and Wu Xiaoping [ 36 ] introduced the conjugate gradient method, which solved the two problems of resistivity 3D forward modeling efficiency and 3D inversion Jacobian matrix calculation and storage, and realized relatively fast and effective resistivity 3D inversion. The optimization methods used in 3D data inversion mainly include nonlinear conjugate gradient method (NLCG), Gauss-Newton method (GN) and quasi-Newton method (QN). NLCG and QN only need the gradient information of the objective function, and do not need to explicitly calculate the sensitivity matrix. Since only the first-order derivative (gradient) information of the objective function is used, there is only linear convergence, and the model needs to be updated many times to converge to the ideal value. , while the GN method uses the information of the first-order derivative (Gradient) and the second-order derivative (Hessian) of the objective function at the same time, showing approximately quadratic convergence, and the number of model updates required is far less than that of the NLCG and QN methods. The normal equation of the sub-Gauss-Newton iteration is solved by the preconditioned conjugate gradient method (PCG), which not only avoids displaying the solution and storing the sensitivity matrix, but also reduces the number of PCG iterations to reduce the calculation time of the inversion. Due to the nonlinear characteristics of the inversion, this inexact solution not only reduces the amount of calculation for the linearized iterative inversion, but also avoids the over-solution of the normal equations in each GN iteration, so there is an advantage in the inversion search process. Chance to jump out of local minima[33]。
An effective way to suppress the multi-solutions of 3D inversion of resistivity is to develop constrained inversion and add the existing prior information into the inversion process by mathematical means to improve the inversion accuracy. Sasaki Y [ 21 ] proposed a smooth constraint, so that the resistivity difference between adjacent grids is extremely small and the transition is smooth. This constraint has natural rationality and has become a commonly used constraint form in resistivity inversion; Kim et al [ 64 ] defined a new parameter vector to represent the model parameters according to the imposition form of smooth constraints, and applied the inequality constraints to the inversion equation, and achieved good results; Li et al[65] used constrained optimization in The interior point penalty function method adds the inequality constraints into the objective function, which improves the accuracy of inversion imaging; Huang Junge et al. [66] and Wan Xinlin et al. [ 67 ] respectively improved the roughness matrix to improve the effect of 3D inversion ; Liu Bin et al . [ 68 ] introduced adaptive weighted smoothness constraints in the inversion objective function to control the deep resolution and discussed the addition of inequalities and prior structural constraints in their literature [ 69-70 ] , eliminating the inversion The false anomalies and redundant structures in the imaging obviously suppress the multi-solution problem of resistivity detection and inversion.
Electromagnetic detection inversion is a typical ill-posed problem, and it is easy to cause multiple solutions in the inversion results. Illness is an inherent feature of the inversion itself, and it is difficult to overcome the essential difficulty without additional information for the solution [ 71 ] . Using related technology to introduce prior information to achieve high-precision imaging of underground structures is one of the important development directions of electromagnetic methods in the future [ 72 ] . The variation range of resistivity is a very important prior information, how to apply inequality constraints to the inversion equation is a key problem. The method adopted by Kim et al. [ 64 ] is an empirical imposition method. In the inversion process, the constraint item is cumbersome and not easy to understand; Li et al. [ 65 ] adopted the interior point penalty function method, adding logarithmic form to the objective function The barrier function, this application form is simple, but the interior point method requires that the selection of the initial value must be within the feasible region, and the initial point should be a feasible point far from the constraint boundary, too close to a certain constraint boundary will cause difficulty in solving . If the constraints are complex, it is difficult to select the initial point, and the interior point method can only solve the situation where the constraint is an inequality, and the exterior point method can solve the situation where the constraint is a mixture of equality and inequality, and the exterior point method has no requirements for the initial point. . Combining the previous research and the problems existing in the above inversion, this paper introduces the inequality constraint representing the variation range of the model parameters into the objective function as a point penalty function outside the prior information [73], and selects the well-ground observation method , unstructured The combination of tetrahedral finite element method and Gauss Newton method (GN) realized the three-dimensional constrained inversion algorithm of well-ground resistivity method, and verified the stability and effectiveness of the developed algorithm by synthesizing observation data with typical three-dimensional theoretical models .
1. 3D forward modeling theory
1.1 Governing equations
The partial differential equation and its boundary conditions satisfied by the total potential of the point current source field in the three-dimensional DC resistivity method are [ 9 , 16 , 74 ] :
∇ ⋅ ( p(x,y,z)∇u(x,y,z))=−4 p.mohAId(rA),inΩ∂u∂n=0,onCs∂u∂n+cosθru=0,onC∞
(1)
where: σ is the electrical conductivity of the rock; u is the electric potential; I is the current intensity; δ is the Dirac function ; r A is the position of the source point;ohA= 2 p
; If the source is underground,ohA= 4 p
, n is the boundaryC∞
The outer normal vector on the boundary; r is the distance from the point power source to the boundary point; Γ s is the natural boundary condition;C∞
is the infinite boundary; θ is the angle between the outer normal vector n of the point on the infinite boundary and the distance r from the field source to the boundary point . The integral equation of the variational problem corresponding to formula (1) can be derived by using the weighted residual method:
F(u)=∫Oh[12p(∇u)2−4 p.mohAd(A)u]dOh+∫C∞12pu2cos(r,n)rdCdF(u)=0
(2)
The calculation area adopts tetrahedron subdivision and linear difference, and finally forms a large sparse symmetric linear equation system. The matrix expression is as follows:
Ku=P
(3)
where: K is an×n
symmetric matrix of ; u is an n × 1 column vector, which represents the potential vector on the three-dimensional grid node; P is a column vector containing field source information; P = (0,…, I ,…,0), in order to save memory, use Row compression storage (compressed sparse row or compressed row storage) processes the coefficient matrix K , and uses the symmetric hyper-relaxation preconditioned conjugate gradient algorithm [ 11 , 20 ] (SSORPCG) to solve the equations (3).
1.2 Algorithm Verification
In order to verify the correctness of the forward modeling program in this paper, a model of a buried sphere in a uniform half space is designed in this paper, as 1a, the model is a spherical anomalous body embedded in a uniform half space [ 16 ] , and the radius of the sphere is R =2.25 m, the coordinates of the center of the sphere are (0,0,-4.5), the resistivity value of the spherer0=1Ω⋅m
, the resistivity value of the formation surrounding rock isr1=10Ω⋅m
, the coordinates of the point current source A are at (-5,0,0), the magnitude of the current is 1 A, the measuring electrodes M are along the x direction, and the spacing is 0.25 m. The size of the entire calculation area is 500 m×500 m×300 m, the size of the target area is 100 m×100 m×100 m, the total number of grid nodes is 273 582, and the total number of grid units is 1 625 838. The grid division effect is shown in Figure 1b ~1d, the encryption processing is carried out at the source point, measuring point and sphere respectively. Using the pole-pole device, the potential value at the measuring point is calculated by the finite element forward modeling program in this paper, and compared with the analytical solution, as solution. The maximum error is less than 0.34% ( see Figure 3 ).
2. 3D inversion theory
2.1 Objective function
According to Tikhonov regularization theory, the inversion objective function in the sense of least squares is [ 3 , 19 , 31 ] :
Φ ( m ) =Phid( m ) + lPhim(m)=(Wd(dobs−F(m)))22+ l(Wm(m−mref))22
(4)
Where: F ( m ) is the forward response function; m is the model parameter(mi,i=1,2,…,M)
, d obs(dj, j = 1.2 , … , N)
is the observation data, W d isN×N
order observation data error weighting matrix,Wd=diag(1/p1,1/p2,…,1/pn)
,pi
is the standard deviation of the i -th observed data; W m isM×M
The first-order model weighting matrix is usually defined by the discrete difference operator of the model unit, and the first-order regularization constraint is generally taken. λ is a regularization factor, which is used to balance the weight of data fitting and model smoothness; m ref is a reference model containing prior information of model parameters. The inversion observation parameter d obs is the potential value of the pole-pole, and the model parameter is the conductivity value of the unitp
, due to the large variation range, in order to stabilize the inversion, the logarithm is usually used to calibrate the observation data and model parameters, that is,d=lnPhiobs,m=lnσ
. For the application of smooth constraints, please refer to reference [ 62 ], for the principle of applying inequality constraints, please refer to references [ 62-63 ], and to construct an inversion objective function with inequality constraints in the conventional resistivity inversion objective function, see formula (5):
Φ ( m ) =Phid( m ) + lPhim(m)+μP(m)=∥Wd(dobs−F(m))∥22+λ∥Wm(m−mref)∥22+μ[(min(m−mmin),0)2+μ(min(mmax−m),0)2]
(5)
In the formula:mmin,mmax
is the upper and lower limits of grid resistivity for each unit. Equation (5) is the objective function of the three-dimensional resistivity inversion with inequality constraints in this paper. Compared with the conventional three-dimensional resistivity inversion objective function, this method adds inequality constraints. The expression is simple, easy to understand, and easy to understand. accomplish. At present, the optimization methods commonly used in 3D resistivity inversion include quasi - Newton method (QN), Gauss - Newton method (GN), nonlinear conjugate gradient method (NLCG) and the variants of these methods [ 33,75 ] . The second-order Taylor expansion is performed on the objective function (5), and the normal equation of the Gauss-Newton method (Gauss-Newton, GN) is obtained by ignoring the higher-order terms:
H(m)△(m)=−g(m)
(6)
Among them: g ( m ) is the first-order derivative of the objective function (gradient matrix, Gradient), H ( m ) is the second-order approximate derivative of the objective function (Hessian matrix, Hessian), respectively expressed as follows:
g(m)=∇Φ(m)=−2JT(m)WTdWd[dobs−F(m)]+2λWTmWm(m−mref)+ 2μ[min(m−mmin,0)−min(mmax−m,0)]H(m)=JT(m)WTdWdJ( m ) + lWTmWm+ mWc
(7)
in:J(m)
is the sensitivity matrix ( the first derivative of the forward response function F ( m )),Wc
is an inequality constraint matrix: its elementsWc(i,j)
The values are as follows:
Wc⎛⎝⎜i,j⎞⎠⎟=⎧⎩⎨⎪⎪0, mmin≤m≤mmax1, m<mmin1, m>mmax
(8)
In the GN method, in order to obtain the modification amount of the model, the method equation (6) needs to be solved for each iteration. It can be seen from the equation (7) that the key to the calculation of the gradient matrix g and the Hessian matrix H lies in the calculation of the sensitivity matrix J , directly Solving equation (6) requires explicit calculation and storage of the sensitivity matrix and Hessian matrix, and involves a large amount of forward calculation, which consumes a large amount of computing memory and computing time, and the existing computing resources are difficult to meet the requirements. In order to avoid the Jacobian For the direct calculation of the matrix, the conjugate gradient method is usually used to iteratively solve the formula (6), forming the GN-CG method. In the CG algorithm, it is only necessary to calculate and store the Jacobian matrix or the product J x andJTy
, so there is no need to consider the storage problem of the Jacobian matrix, and J x andJTy
The calculation of is equivalent to solving a forward modeling problem, which greatly speeds up the calculation speed of the inversion. Zhang et al. [ 26 ] in their paper calculated J x andJTy
A detailed analysis was carried out.
2.2 Model update
Using the CG method to obtain the solution Δ of the normal equation (6)m
Finally, a linear search is required to obtain the model update amount of each GN iteration, and the update formula is:
mk+1=mk+ a D m
(9)
Among them: α is the update step size of the model. The accurate line search of the step size α often requires a lot of calculations. Generally, the non-exact one-dimensional linear search method is used to obtain the step size. The common strategy is to use the Wofe-Powell criterion to perform non-exact linear The update step α is obtained by searching , and the Wofe-Powell criterion includes sufficient descent conditions and curvature conditions, and the formula is:
F (mk+akpk)≤Φ(mk)+c1ak∇Φ(mk)Tpk∇Φ(mk+akpk)Tpk≥c2∇Φ(mk)Tpk
(10)
in:Phi
is the forward operator,mk
is the model parameter of the k -th inversion iteration;ak
is the iteration step;pk
is the search direction; c 1 and c 2 are constants, satisfying0<c1<c2<1
, please refer to [ 73 ] for the specific search method . In general, the smaller c2 is , the more accurate the linear search is,c2=0.1
yields a fairly accurate linear search, takingc2=0.9
A rather weak linear search is obtained, the smaller c 2 is, the longer the search time is, according to the previous research [ 45 ] , usually the GN method and the L_BFGS method generally takec1=10−4,c2=0.9
。
Social Plugin