NL2031632B1

NL2031632B1 - Gravity inversion method and system based on meshfree method

Info

Publication number: NL2031632B1
Application number: NL2031632A
Authority: NL
Inventors: Hu Jinmin; Zhang Yongqian; Huang Yao; Liu Yan
Original assignee: Chinese Acad Of Geological Sciences
Priority date: 2022-04-20
Filing date: 2022-04-20
Publication date: 2023-11-07

Abstract

A gravity inversion method and system based a meshfree method are provided. An appropriate approximate function construction method is selected; appropriate parameters of the hybrid, radial basis function are selected, with an appropriate judgment method; an equation construction mode is selected; weighting is carried out on the distance of the hybrid radial basis function according to ore body occurrence and form in prior information; known underground density information is loaded, and a system of equations is constructed with a hybrid radial basis function meshfree point collocation method; the system. of equations is solved, and global estimated density distribution is solved with a coefficient matrix generated by matching the obtained coefficient vector through a global background mesh; observation data is loaded, and inversion is carried out by using a pre—optimal conjugate gradient method and taking the estimated density distribution as a constraint; underground density distribution is obtained to complete inversion.

Description

P1286 /NLpd

GRAVITY INVERSION METHOD AND SYSTEM BASED ON MESHEFREE METHOD

TECHNICAL FIELD

The disclosure relates to the technical field of gravity mod- el inversion, and in particular to a gravity model inversion meth- od and system based on a meshfree method.

BACKGROUND ART

Gravity prospecting is a geophysical method which is very im- portant in various fields such as mineral exploration, oil and gas exploration, and structural interpretation. Gravity inversion is one of the most important means in gravity prospecting data pro- cessing. In the current gravity inversion, the superior methods include a simulated annealing method based on random search, a neural network algorithm, Bayesian inversion based on statistics, a preconditioned conjugate gradient method based on nonlinear so- lution and the like. Relatively speaking, the study on geophysical inversion using meshfree methods is still in the beginning stage.

The meshfree methods have become a new class of computational methods with considerable success, but their application in the field of geophysics is still in the beginning of development.

From application examples, they have been applied to various branches of geophysics, indicating that they can be well applied to the problem of geophysics. For the forward and inversion of gravity, two-dimensional forward of gravity can already be solved with the meshfree methods, but there are not good solutions for the inversion. The disclosure uses meshfree local Petrov-Galerkin method in combination with hybrid radial basis functions to solve the problem of gravity inversion.

SUMMARY

The technical problem to be solved by the disclosure is to provide a gravity model inversion method and system based on a hy- brid basis function meshfree method. By adopting the hybrid radial basis function meshfree method, density field estimation is car- ried out on known prior information through the hybrid basis meshfree method, and then the estimated density field serves as initial density for inversion. With the method, the prior infor- mation can be efficiently utilized, so that the accuracy of gravi- ty inversion is improved.

In order to solve the above-mentioned technical problem, the disclosure provides a gravity model inversion method and system based on a meshfree method. The steps specifically include: a hy- brid radial basis function formed by a multiguadric (MQ) radial basis function and a shape parameter independent cubic kernel function is used; appropriate parameters of the hybrid radial ba- sis function are selected by using an appropriate judgment method; in the method, root mean square deviation (RMSD) and degree of freedom are selected for judgment; norm weighting is carried out on the distance norm of the hybrid radial basis function according to orebody occurrence and form in the prior geological infor- mation; known underground density information is loaded, and a system of equations is constructed by using a hybrid radial point interpolation method; inversion is carried out by using a precon- ditioned conjugate gradient method and taking estimated density distribution as a constraint; and underground density distribution is obtained to complete inversion.

The hybrid basis function may be used for interpolation of discrete data and can get rid of the constraints of the mesh. The cubic kernel function is introduced, so that the situation that the stability of the equation is affected by too large condition number is well avoided, and the accuracy is improved.

The disclosure has the advantages that the gravity inversion method and system based on the meshfree method are provided. With the hybrid radial point collocation method, the prior geclogical information can be efficiently utilized, and the known orebody form and occurrence can be brought into the constraint of inver- sion. With the radial distance as the only variable for calculat- ing the estimated density distribution, the calculation complexity can be simplified, meanwhile, the physical rule constraint can al- so be well met, and thus the efficiency and accuracy of inversion are improved.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing is merely an overview of the technical solution of the disclosure, in order to clearly understand the technical means of the disclosure, the disclosure will be further described in detail below in combination with the drawings and specific im- plementation modes.

FIG. 1 illustrates the influence of different values of pa- rameter c on the shape of an MQ function.

FIG. 2 illustrates the influence of different values of pa- rameter 8 on the shape of an MQ function.

FIG. 3 illustrates the influence of different values of pa- rameter £# on the shape of an MQ function.

FIG. 4 illustrates a shape diagram of a cubic function.

FIG. 5 illustrates a shape diagram of a hybrid basis func- tion.

FIG. 6 illustrates a diagram of the variations in the root mean square error with shape parameters.

FIG. 7 illustrates a diagram of the variations in the condi- tion numbers with shape parameters.

FIG. 8 illustrates a diagram of the variations in the condi- tional numbers with degree of freedom.

FIG. 9 illustrates a diagram of the variations in the root mean square error with degree of freedom.

FIG. 10 illustrates a weighting schematic diagram of a kernel function.

FIG. 11 illustrates a schematic diagram of a model.

FIG. 12 illustrates an estimation model diagram adopting a meshfree method.

FIG. 13 illustrates a forward result diagram of a model.

FIG. 14 illustrates a traditional inversion result diagram.

FIG. 15 illustrates an inversion result diagram with prior information estimation adopting a meshfree method.

DETAILED DESCRIPTION OF THE EMBODIMENTS

The preferred embodiments of the disclosure will be described below with reference to the drawings, and it is to be understood that the preferred embodiments described herein are for the pur- pose of illustrating and explaining the disclosure only and are not intended to limit the disclosure.

The first step: parameter selection of a radial basis func- tion

The radial basis function used herein is an improved near function of a MQ function proposed by Hardy, and the specific form of the MQ function is: 2 2 = g (x)= (er) +c) (7)

In the formula, r: is the distance norm of the calculation point and the node. So far, no uniform and effective parameter se- lection methods are available. To this end, it is necessary to se- lect and evaluate parameters. Aiming at the above-mentioned prob- lem, corresponding studies have been conducted on formula (1), with the result shown in FIGs. 1 to 3. As can be seen from FIG. 1, the influence of c on the MQ function is embodied in two aspects.

On the one hand, it can change the sharpness of the radial basis function, i.e., change the weight distribution inside the influ- ence domain. On the other hand, its magnitude determines the mag- nitude of an extreme value of the radial basis function, and c is exponentially related to the maximum value.

As can be seen from FIG. 2, the influence of © on the MQ function is reflected in overall rise and fall, the smaller 6 is, the larger all values of the entire function are. As can be seen from FIG. 3, the influence of ese on the MQ function is reflected in a change to its smoothness. The smaller its value is, the smoother the obtained radial basis function is, the influence of e is dif- ferent from that of c, it can independently change the smoothness without influencing the magnitude of an extreme value.

The MQ function of the general form results in solving an ill-conditioned system of equations, and small errors in the data may result in significant large errors in interpolation solution.

In order to reduce the limitation, a traditional MQ function and a cubic kernel function independent of shape parameters are combined to obtain a hybrid basis function. The specific form is: br) =all+ (er) + pr, (8}

The shape of the cubic kernel function is shown as FIG. 4.

The shape of the hybrid kernel function is shown as FIG. 5. As can be seen from FIG. 5, a determines the component ratio of the MQ radial basis function, and also determines the extreme value of the center. B determines the components of the cubic kernel func- 5 tion. As with the previous property, the larger the value of eg, the sharper the obtained radial basis function is.

The second step: evaluation method for parameter selection

Both the MQ radial basis function and the hybrid basis func- tion have the problems of excessive parameters and difficulty in selection. For this, some parameters are introduced to evaluate whether the arranged node spacing and the selection of the shape parameters of the basis function are excellent or not. 1. Influence of shape parameters on root mean square error: shape parameters have a great influence on the root mean square error of the solution value, the variations in the root mean square error with shape parameters shown in FIG. 6. The overall error of the hybrid basis function used herein is obviously lower than that of the traditional method. 2. Influence of shape parameters on condition number: the variations in the condition numbers with shape parameters is shown in FIG. 7. As can be seen from the figure, the value of & can be selected, so that the root mean square error and equation singu- larity are relatively small at the same time. 3. Influence of degree of freedom on condition number: the variations in the conditional numbers with degree of freedom is shown in FIG. 8. As can be seen from the figure, in the tradition- al method, when the degree of freedom is relatively small, the condition number obviously increases along with the increase of the degree of freedom, and then tends to be stable. The hybrid ba- sis function is obviously more excellent than the traditional method in the aspect of influence of degree of freedom. 4. Influence of degree of freedom on root mean sguare error: the variations in the root mean square error with degree of free- dom is shown in FIG. 9. As can be seen from the figure, when the degree of freedom is relatively small, the root mean square errors can be reduced along with increase of the degree of freedom in all the three methods. However, the root means square error of the method based on gaussian function and MQ function will increase sharply when the degree of freedom reaches a certain value. This condition does not occur with the improved hybrid basis function.

Therefore, the hybrid basis function can better adapt to the change of the degree of freedom.

The third step: selection of an equation construction mode

The weighted residual method is a mathematical method capable of directly solving an approximate solution from a differential equation, and is a powerful method for constructing approximate solutions of ordinary differential equations and partial differen- tial equations and discrete system equations. The governing equa- tion for the studied problem can be expressed as follows:

To = f(x) xeQ

B(u(x))=g(x) xel (2)

Where L is the differential operator within the solution do- main ©, B is the differential operator on the boundary T, and u (x) is the field variable within the solution domain. The margin form of the equation is obtained according to formula (2):

IE =I(u(x))- f(x)=0 xe

R. =B(u(x))-g(x)=0 xel (2)

The margin equation is constructed in the integral form: [wr do + [vr dr =0

Q ï (4) w and v are weight functions. When the formula meets the spe- cific conditions that the basis function is continuous at a cer- tain order, the weight functions and the basis function are line- arly independent, etc., an approximate solution converging to an exact solution can be obtained.

The used hybrid basis function generated based on the MQ ra- dial basis function already contains the response to distance. Af- ter the proportioning problem of the weight functions is reduced, the equation is higher in flexibility, and combination with inver- sion becomes more feasible and efficient.

The fourth step: re-weighted radial basis function

In actual underground geclegic bodies, especially in the field of metal mineral exploration, ore bodies often have obvious structural characteristics. In order to solve the above-mentioned problem, a distance re-weighting method is used to enable the method to be more suitable for estimation of underground geologic ore bodies.

First, in the aspect of occurrence, when the occurrence of the geologic bodies is consistent with that of the basis function, the optimal calculation result can be obtained. The specific form of the re-weighted distance is:

I. d =|Ax, A*S*G*R*G'| Ax, Av | (5) a ind —sin@ cosd (6) eo op (7) where Ax and Ay represent the horizontal and vertical dis- tances, respectively, in the two-dimensional case. ò is a scaling parameter, G represents a rotation matrix, and © is a rotation pa- rameter. R is a stretching matrix, and p is a stretching parame- ter. A specific influence mode can be seen in FIG. 10. In summary, the scaling parameter is to change the smoothness degree of the basis function, the stretching matrix is to match the approximate form of the estimated ore bodies, and the rotation matrix is to estimate the approximate tendency occurrence of the ore bodies.

Actual operation process

A workplace with the depth of 500 m and the width of 1000 m is established, the background mesh spacing is 20 m, and a hypo- thetical geologic ore body is established in the work area. A rec- tangular two-dimensional ore body with the length of 140 m and the width of 80 m is used herein, and its residual density is 5 g/cm’.

See FIG. 11 for specific representations.

The gravity inversion method based on the meshfree method is characterized in that the prior information is reasonably utilized and is introduced into the inversion process. Therefore, it is as- sumed that a certain amount of scatter density information has been obtained by methods such as well logging. Then the hybrid ba-

sis function with the appropriate parameters is chosen to be sub- stituted into the equation:

U) ~u*(X)=Y R(X)b, = R7(X)b i=l (8)

In the formula, R'(X) is a weighting hybrid basis function, with the expression shown in formula (8) and the weighed distance shown in formula (5). b is a coefficient vector. The obtained co- efficient vector is substituted into global background to solve estimate ore body distribution, as shown in FIG. 12. Compared with prior constraint in the traditional method, the method pays more attention to overall grasp.

The forward value of the model is calculated according to the forward calculation formula, as shown in FIG. 13.

Ag = & = aj ZED

WIE CT ©)

As can be seen from the forward result diagram of the model, the forward value is the maximum at the position above the ore body, and is gradually reduced along with the distance away from the ore body, so that the basic distribution rule is met. Mean- while, as the buried depth of the orebody is large and the trans- verse extension is insufficient, the inversion difficulty increas- es.

In inversion, a preconditioned conjugate gradient method is mainly used as a solving method. For gravity imaging inversion, the condition number of the coefficient matrix is quite large, so that the convergence of a conjugate gradient method is influenced.

In order to solve the problem, a preconditioned matrix is added to make it more stable. The specific form is:

PG GAm = PG’ Ad (10) where P is the preconditioned matrix and is taken as the re- ciprocal of the square of the depth. G is a system matrix, Am is a model correction, and Ad is observation data. The characteristic values of PG'G may be centrally distributed on the diagonal line, and the condition number of the equation is optimized.

The traditional inversion objective function basically satis-

fies formula (11).

Hm) = ig — g(m)) +4, Wim -H| + A, T(m) (17)

LL 2 r(m)=>Y > [(m —m) +&]? =H (12)

Where g°®* is an observed value, g (m) is an estimated value,

Wis a model weighting matrix, H is prior underground information,

T is a stable functional, and A is a regularization parameter. In the meshfree gravity inversion, the prior information and stable functional need to be unified in the same matrix formed by a radi- al basis function.

A traditional inversion result without constraint of prior information, as shown in FIG. 14, and an inversion result with prior model estimation adopting a meshfree method, as shown in

FIG. 15, are subjected to contrastive analysis. As can be seen from the figures, the inversion result without prior information shows that its center basically conforms to the position of the model, but it is in a concentric oval shape and is insufficient in focusing capacity, so that the shape of the whole orebody is much larger than the real shape. While the inversion result with esti- mation adopting the meshfree method can be well matched with the real shape of the model, and the defect of excessive smoothness is avoided.

Claims

CONCLUSIONS

1. Gravity inversion method and system based on a mesh-free method, in which a hybrid basis function mesh-free method is applied, density field estimation is performed on known prior information via the hybrid basis mesh-free method, and then the estimated density field serves as initial tial density for inversion, in particular as follows: (1) a suitable approximate function construction method is selected, and a hybrid radial basis function formed by an MQ radial basis function and a shape parameter-independent cubic kernel function is used; (2) Appropriate parameters of the hybrid radial basis function are selected by using an appropriate assessment method, and the selected parameters of the hybrid radial basis function are the weight of each component and the flatness of the components of the MQ radial basis function; {3) a comparison construction mode is selected, a strong construction is selected in the method, because the radial basis function has a form similar to that of a weak weight function, the operation can be simplified; standard weighting is carried out on the distance standard of the hybrid radial basis function according to the occurrence and shape of ore in the previous information; (4) known subsurface density information is loaded and a system of equations is constructed using a hybrid radial base point collocation method; (5) the system of equations is solved, and global estimated density distribution is solved using a coefficient matrix generated by matching the obtained coefficient vector with a global background mesh; (6) observation data is loaded, and inversion is performed by using a preconditioned conjugate gradient method and taking the estimated density distribution as a constraint; (7) subsurface density distribution is obtained to complete inversion.

The method according to claim 1, wherein the hybrid radial base point collocation method is applied to achieve the purpose that the prior geological information is efficiently used, and the prior shape and appearance of the ore body can be brought within the constraint of inversion .

The method of claim 1, wherein with the radial distance as the only variable for calculating the estimated density distribution, the calculation complexity can be simplified while also satisfying the physical control constraint, and thus the inversion efficiency and accuracy can therefore be effectively improved.