CN113255230A - Gravity model forward modeling method and system based on MQ radial basis function - Google Patents
Gravity model forward modeling method and system based on MQ radial basis function Download PDFInfo
- Publication number
- CN113255230A CN113255230A CN202110665099.7A CN202110665099A CN113255230A CN 113255230 A CN113255230 A CN 113255230A CN 202110665099 A CN202110665099 A CN 202110665099A CN 113255230 A CN113255230 A CN 113255230A
- Authority
- CN
- China
- Prior art keywords
- radial basis
- basis function
- function
- forward modeling
- gravity model
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 101
- 230000005484 gravity Effects 0.000 title claims abstract description 59
- 238000009826 distribution Methods 0.000 claims abstract description 9
- 239000011159 matrix material Substances 0.000 claims description 10
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 abstract description 24
- OLTSGVZGKOFTHZ-UHFFFAOYSA-N P.P.P.P.P.P.P.P.P Chemical compound P.P.P.P.P.P.P.P.P OLTSGVZGKOFTHZ-UHFFFAOYSA-N 0.000 description 29
- 230000008569 process Effects 0.000 description 8
- 230000008901 benefit Effects 0.000 description 7
- 238000004422 calculation algorithm Methods 0.000 description 7
- 238000002474 experimental method Methods 0.000 description 7
- 238000011161 development Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 5
- 238000011160 research Methods 0.000 description 5
- 238000004088 simulation Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- XEEYBQQBJWHFJM-UHFFFAOYSA-N Iron Chemical compound [Fe] XEEYBQQBJWHFJM-UHFFFAOYSA-N 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 2
- 238000013528 artificial neural network Methods 0.000 description 2
- 238000013527 convolutional neural network Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 230000000149 penetrating effect Effects 0.000 description 2
- 238000005293 physical law Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 241000238097 Callinectes sapidus Species 0.000 description 1
- 229910000831 Steel Inorganic materials 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000037396 body weight Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 239000003245 coal Substances 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 1
- 239000010931 gold Substances 0.000 description 1
- 229910052737 gold Inorganic materials 0.000 description 1
- 229910052742 iron Inorganic materials 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000009828 non-uniform distribution Methods 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
- 239000010959 steel Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/04—Architecture, e.g. interconnection topology
- G06N3/048—Activation functions
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- Software Systems (AREA)
- Artificial Intelligence (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Health & Medical Sciences (AREA)
- Mathematical Physics (AREA)
- Computational Linguistics (AREA)
- Molecular Biology (AREA)
- Computing Systems (AREA)
- Biophysics (AREA)
- Biomedical Technology (AREA)
- Data Mining & Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Medical Informatics (AREA)
- Computer Hardware Design (AREA)
- Geometry (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention provides a gravity model forward modeling method and system based on MQ radial basis functions. The method comprises the following steps: gridding the underground space of the working area; adjusting parameters of the radial basis function according to the problems to be solved; the collected actual data (including well data) or the calculated model data are brought into a function equation taking a radial basis function as a shape function, and a weight coefficient vector is solved; and solving the information of the global grid point by using the weight coefficient vector and the points with known information to finally obtain the global gravity field distribution to finish forward modeling. The gravity model forward method and the gravity model forward system based on the MQ radial basis functions provided by the invention adopt a step-by-step calculation method, shorten forward calculation time, better meet the precision requirement, and reduce the number of forward parameters by using the radial basis functions.
Description
Technical Field
The invention relates to the technical field of forward modeling of gravity models, in particular to a forward modeling method and a forward modeling system of a gravity model based on MQ radial basis functions.
Background
There are two main types of methods for establishing gravity model numerical values: a grid-based method and a non-grid method. The mesh-based method is a conventional method in the past, including a difference method, a finite element method, and the like, and has been widely applied to many disciplines and projects. The denser the grid of the difference method is, the more nodes are, the higher the solution precision is. Because the orthogonal grids with uniform sizes are adopted, when the shape of the solved domain boundary is complex, the solving precision of the difference method is limited, and even the solving cannot be carried out. When the finite element method is used for treating the problems of continuous grid reconstruction or large deformation and the like, the finite element grid can be seriously distorted, not only the grid reconstruction is required, but also the understanding precision is seriously influenced [ stretch-breaking, Liu rock, Su. [ theory and application of meshless method [ J ]. mechanical progress, 2009,39(01):1-36 ]. The gridless method is a node-based numerical method, avoids complex mesh subdivision, provides a new method for solving a problem by numerical values [ Lijunjie, rigorous bin. development of gridless method and application in geophysics [ J ]. geophysical development, 2014,29(01): 452-.
The radial basis function method belongs to a meshless method, has infinite smoothness and a derivative of a function about a certain coordinate direction, can be expressed as a weighted linear combination of function values at all discrete points along the direction [ Wanfangzuan, Liao soldier, Xiemang, differential integration method characteristics and improvement [ J ]. computational mechanics reports, 2015,32(06): 765) 771 ], and has a good effect on solving various partial differential equations (such as elliptic, parabolic and hyperbolic). Common radial basis function methods are: local radial basis function method [ C.K.Lee, X.Liu, S.C.Fan.local multiple adaptation for solving boundary value schemes [ J ]. computerized Mechancs, 2003,30(5-6) ], differential product-based radial basis function method [ C.Shu, H.Ding, H.Q.Chen, T.G.Wang.An upward with local RBF-DQ method for mapping of indirect computing flows [ J ]. Computer Methods in Applied Mechanics and Engineering,2004,194(18) ], hermitian radial basis point interpolation method [ Xin Liu, G.R.Liu, Katai, K.Y.Lang.Lam.polar (36) of compressive analysis [ P.P.E.L.P.L.P.P.L.P.P.L.P.P.P.L.P.L.P.P.E.P.P.L.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.L.P.P.P.P.L.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.No.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.No.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.No.P.P.P.No.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.P.
Parameter selection of radial basis functions is a difficulty which always exists along with radial basis development, and Pankaj k.mishra et al discuss ill-conditioned research of solving Poisson's problem by coefficient matrixes composed of several infinite smooth shape functions (Gaussian functions, MQ functions and inverse MQ functions) in detail, disclose the relation between condition number, root mean square error and shape function shape parameters, and provide a basis for selection of radial basis function shape parameters [ Pankaj k.mishra, Sankar k.Nath.on the conversion of Iterative Methods and pseudo-inverse applications in Global processing theory [ J ]. International Journal of Applied and Computational principles, 2017,3(4) ]. Bear positive exceeds points out numerical instability of radial basis functions, and the selection of optimal shape parameters and the variance of kernel functions are closely related [ bear positive exceeds, a plurality of problems in radial basis function approximation research [ D ]. university of redun, 2007 ]. Ahmad Golbai details the difference between the normal shape parameter and the deformed shape parameter, and constructively proposes two methods of deformed shape parameters, namely Binary Shape Parameter (BSP) and mixed shape parameter (HSP), which have better array calculation performance [ Ahmad Golbai, Ehsan Mohbanfar, Hamed Rabei. on the new variable shape parameter sequences for radial bases functions [ J ]. Computational and Applied Mathematics,2015,34(2) ].
In the aspect of selecting the support domain, the selection of the shape function support domain has great influence on the numerical precision of the solution. The common method for selecting the global support domain usually increases the operation cost when solving a large-scale problem, the local circular domain support also faces the troubles of difficult coefficient selection and reduced precision, among them, the better improvement is a node integral unitless Galerkin method [ J ] of windward deviation of integral nodes, 2012,29(02):183-190 ] and an optimization method for a non-latticed function influence domain based on a convolutional neural network [ Liuyuxiang, Wang east, fan present, Chenjian, Housonyang ] optimization research of the non-latticed function influence domain based on the convolutional neural network [ J/OL ]. report of solid mechanics: 1-19[2021-04-27]. https: orf/10.19636/j.cnki.cjsm42-1250/o3.2021.003 ].
The non-uniform distribution point format is combined with the MQ radial basis function to be used for numerically simulating the boundary value problem of the partial differential equation, a satisfactory result is obtained [ Zhu Zi, Zeng Ching, Lu Guang Yin, Jie Geng, the gravity tensor finite element forward modeling of the second-degree body [ J ]. geophysical prospecting and chemical prospecting, 2010,34(05):668-671+685 ] ], particularly, the solution of the elliptical partial differential equation provides theoretical support for solving the geophysical problem by using the radial basis function, and the potential field equation is the elliptical partial differential equation.
Radial basis functions or even meshless methods have relatively few applications in gravity magnetic prospecting and all geophysical methods, however, the corresponding finite element method has a great deal of development [ gold steel harmonize, huxiangyun, supreme, south minggen, healthy men, jinjing qi, liuhui, and complex shape body weight magnetic anomaly isoparametric finite element integration algorithm research [ J ]. petroleum geophysical exploration, 2009,44(02):231 + 239+129+255.] [ Zhu is from the force, Zeng is red, luguanyin, Jie Geng. the gravity tensor finite element forward modeling of the second degree body [ J ]. geophysical prospecting and chemical prospecting, 2010,34(05):668 one 671+685 ], [ Zhaodan, Dujian, Guo Zhiyong, Xuwei, Liu Zhongzi ] magnetic anomaly simulation and detection identification of buried parallel iron pipelines [ J ]. underground space and engineering newspapers, 2020,16(03): 891-896+932.]. The conventional finite element method often has difficulty in obtaining good results when dealing with many practical problems, such as: in the inversion numerical solution of the high-dimensional geophysical property, a discrete point model can be hardly constructed, so that a meshless method has development opportunities.
A unit-free Galerkin method is used for the Queen-King to carry out seismic wave forward modeling [ the King-Yueying-seismic wave forward modeling in the unit-free Galerkin method initially probes the geophysical progress [ J ]. 2007(05):1539-1544 ]. Von German mountain and others make forward modeling of ground penetrating radar based on the unit-less Galerkin method [ Von German mountain, King Honghua, Daidaiwei ] forward modeling of ground penetrating radar based on the unit-less Galerkin method [ J ] geophysical science, 2013,56(01): 298-. The two-dimensional potential field is extended upwards by a Galerkin method by a Naphus, a Lipeng, a Lijing, a grid-free method for the two-dimensional potential field extension is researched [ J ] geophysical prospecting calculation technology, 2017,39(02): 155-. Li Jun et al performed a two-dimensional forward modeling of gravity anomaly [ Li Jun, strictly home ] a meshless method [ J ] in the two-dimensional forward modeling of gravity anomaly [ coal field geology and exploration, 2018,46(06):181 plus 186.] and a meshless solution of magnetotelluric two-dimensional variation problem [ Li Jun, strictly home ] a finite element-radial basis point interpolation method [ J ] in the two-dimensional forward modeling of magnetotelluric [ nonferrous academy of China, 2015,25(05):1314 plus 1324 ].
In summary, the radial basis function method has been well studied in numerical theory, but has been applied to geophysical fields only in a few cases. In view of its excellent numerical simulation capability, it is believed that it has great development potential in the field of geophysical.
The most common and classical method for gravity forward is to use the formula of gravity calculation:
and (3) performing two-dimensional and three-dimensional geological subdivision on the underground space where the geologic body is located, calculating the gravity values of each unit by using a formula (1) and accumulating the gravity values together, thereby finishing the forward modeling of the gravity values of all points. And then, calculating point by point in the observation area to obtain the gravity values of all position points in the area, thus finishing the forward calculation of the gravity. Although the point-by-point calculation method can obtain accurate gravity values of a working area, the time consumption is high, the time requirement of repeated forward inversion in large-scale calculation is not facilitated, and the precision control is complicated.
The method for solving the gravity values of all observation points in a measurement area by adopting a point interpolation radial basis function comprises the following implementation processes: selecting a polynomial function as the linear basis function, e.g. PT(x,y)＝{1,x,y,x2,xy,y2,x3,x2y,y3Calculating the gravity value of a small number of random position points, and then according to the formula c, changing the gravity value to P-1U calculates a coefficient vector c, where U is a series of calculated gravity values corresponding to P (x, y) location points. And then calculated according to the following formula,
U(X,Y)＝PT(X,Y)P-1U (2)
wherein U (X, Y) is the weight of an unknown pointForce value, PT(X, Y) is the corresponding basis function matrix.
Compared with the traditional grid method, the method of point interpolation can obtain the required distribution of the gravity values with less calculation amount. The step of adding point interpolation in the forward and backward evolution process not only can not increase the operation amount, but also can save the operation time, and the positions and the number of discrete points can be artificially selected to carry out multi-scale analysis and local analysis, thereby reducing the number of inversion parameters and the ill-conditioned degree of gravity inversion.
The direct forward modeling with the formula (1) consumes a lot of time, and the coefficient matrix formed in the corresponding inversion process is also extremely large. The inversion process is time-consuming, sometimes even difficult to invert, and meanwhile, the solution difficulty is increased due to too many inversion parameters.
Although the number of solving parameters can be reduced and the solving time is shortened by using the point interpolation function, the accuracy is found to be lower after experiments. Particularly, since the calculation method of the basis function is a power function, the value at a position of the grid far from the origin position may be large, resulting in poor stability. The actual underground space is not different due to the position of the selected coordinate origin, and the value of the actual underground space is only related to the position distance between the geologic body and the observation point. In general, pure point interpolation methods have limited accuracy.
Disclosure of Invention
The technical problem to be solved by the invention is to provide a gravity model forward modeling method and system based on MQ radial basis functions, a step-by-step calculation method is adopted, forward calculation time is shortened, the precision requirement can be well met, and the number of forward parameters is reduced by using the radial basis functions.
In order to solve the technical problem, the invention provides a gravity model forward modeling method based on MQ radial basis functions, which comprises the following steps: gridding the underground space of the working area; adjusting parameters of the radial basis function according to the problems to be solved; the method comprises the steps of substituting collected actual data or calculated model data into a function equation with a radial basis function as a shape function, and solving a weight coefficient vector, wherein the actual data comprises logging data; and solving the information of the global grid point by using the weight coefficient vector and the points with known information to finally obtain the global gravity field distribution to finish forward modeling.
In some embodiments, the radial basis functions include: MQ radial basis functions.
In some embodiments, the MQ radial basis functions have the form:
wherein r isIIn order to calculate the norm of a point and a node, c and theta need to be flexibly selected according to actual conditions.
In some embodiments, gridding the subsurface space of the work area includes: a series of Pi values are defined from surface observation data and well data.
In some embodiments, the background mesh nodes are calculated by affine transformation as follows:
in some embodiments, the radial basis function is expressed as:
wherein R isT(X) is a radial basis function, and b is a coefficient vector.
In some embodiments, the matrix form of r (x) is as follows:
in addition, the invention also provides a gravity model forward modeling system based on the MQ radial basis functions, which comprises: one or more processors; storage means for storing one or more programs; when executed by the one or more processors, cause the one or more processors to implement an MQ radial basis function-based gravity model forward method in accordance with the foregoing.
After adopting such design, the invention has at least the following advantages:
compared with the traditional forward modeling of the gravity model, the forward modeling method has the advantages that the forward modeling calculation time is shortened by adopting a step-by-step calculation method, the precision requirement can be well met, and the number of inversion parameters and the ill-condition of gravity inversion are reduced by utilizing the radial basis function. The method can save the operation cost on a neural network algorithm or other intelligent algorithms which need to be subjected to forward and backward operations for many times;
compared with other interpolation algorithms, the MQ-based radial basis function has the advantage of high similarity with the components of the target function, namely the MQ radial basis function is highly similar to the gravity field generated by a unit gravity source.
Drawings
The foregoing is only an overview of the technical solutions of the present invention, and in order to make the technical solutions of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and the detailed description.
FIG. 1a is a graph of the effect of different values of c on the shape of the MQ function;
FIG. 1b is a graph of the effect of different values of θ on the shape of the MQ function;
FIG. 2 is a theoretical work area layout diagram;
FIG. 3a is a schematic diagram of a uniform node placement method;
FIG. 3b is a schematic diagram of a Chebyshev node placement method;
FIG. 4 is a process for performing the method;
FIG. 5 is a schematic view of a model;
FIG. 6 is a forward model gravity field diagram;
FIG. 7 is a graph of the results of point interpolation;
FIG. 8 is a schematic diagram of a point interpolation error matrix;
FIG. 9 is a graph of radial basis gridless interpolation results;
FIG. 10 is a schematic diagram of a radial basis meshless error matrix;
FIG. 11 is a comparison of results of surface gravity anomalies;
FIG. 12 is a comparison graph of different solution results of trigonometric functions;
FIG. 13 is a comparison graph of the results of a six-point trigonometric function collocation and a radial basis function collocation solution.
Detailed Description
The preferred embodiments of the present invention will be described in conjunction with the accompanying drawings, and it will be understood that they are described herein for the purpose of illustration and explanation and not limitation.
Compared with the traditional forward modeling of the gravity model, the forward modeling calculation time is shortened by adopting a step-by-step calculation method, the accuracy requirement can be better met, and the number of forward modeling parameters is reduced by utilizing the radial basis functions. The method can save the operation cost on a neural network algorithm or other intelligent algorithms which need to be subjected to forward and backward operations for many times; compared with other interpolation algorithms, the MQ-based radial basis function has the advantage of high similarity to the components of the target function, namely the MQ radial basis function is high similar to the gravity field generated by a unit gravity source, and a foundation is provided for later radial basis function inversion.
The most important parts of the proposed method are the selection of the radial basis functions and the selection of the parameters of the radial basis functions. In order to select the comparison well, an observation experiment and a verification experiment are carried out, and then a forward process is carried out by using the selected proper parameters.
Parameter selection of radial basis function
The radial basis function used at this time is a Multi-Quadric function provided by Hardy, and the specific application form is as follows:
in the formula rIFor calculating distance norm of point and node. The values of the shape parameters c and theta need to be discussed, and proper values can better help the accuracy of the radial basis interpolation result. At present, a selection scheme and a comparison method for unifying the optimal values of the values do not exist, and the values of c and theta are required to be flexibly selected according to actual conditions. And the selection of parameters is often determined based on the accuracy of the simulation results, which makes the parameter selection very complicated and difficult to popularize.
For this problem, we have made a corresponding study on equation (3), and the results are shown in fig. 1. As can be seen from fig. 1(a), the influence of c on the MQ function is embodied in weight assignment of the influence domain, and when c is smaller, the MQ function is sharper, so that the influence domain of the node is smaller, resulting in insufficient utilization of node information; when c is too large, although the influence domain of the node can be enlarged, in practical situations, the node often does not have a large influence range due to physical laws, and the calculation amount is increased. This is also the reason why the previous research often supports domain partitioning for nodes. On the premise that the support domain is divided, taking a larger c value can cause a serious long lattice phenomenon, and an ideal result cannot be obtained. And the magnitude of the c value is similar to the upward continuation method of gravity data processing. The study on θ is shown in fig. 1(b), and it can be seen that the selection of θ does not have a significant effect on the shape of the function. In addition, the indexes in the theoretical derivation formula of gravity have strict physical theory, and the index items in the shape function are not suitable to be adjusted greatly. If necessary, the adjustment should be finely adjusted near the index of the corresponding theoretical formula so as to meet the requirement of the overall error.
Second, practical application model of method
The real application environment of the method is described, as shown in fig. 2. The study area is ABB ' A ', AA ' is the surface observation surface in the study area, and a series of P obtained from the ground observation data and the well dataiValues of P fromiThe values are then spread over the study area using grid points, and then forward and inverse calculations are performed.
The arrangement of the grid points provides two methods, as shown in fig. 3. Wherein, the Chebyshev node defined in the solving interval can be calculated by affine transformation according to the following formula for any interval [ a, b ],
the advantage of using the Chebyshev node is that the boundary conditions can be better fitted.
Third, the actual operation process
The actual operation process is shown in fig. 4.
The forward model was established as a horizontal cylinder centered at point (50, -50), with a radius of 10 meters and an infinite y-axis extension, and a residual density of 2g/cm3。
According to the formula of the residual linear density, lambda is sigma, S is sigma, pi R2And further deducing a gravity anomaly expression of each point of the calculation area:
in the formula, delta g is a calculated gravity value, lambda is the linear density of the cylinder, D is the buried depth of the central axis, and x is the horizontal direction horizontal coordinate. The model is shown in fig. 5, the center is a field source (a cylinder extending horizontally and infinitely), and the background node is divided into parts shown by small blue points in the graph.
(2) According to the split nodes and the gravity forward calculation formula, a forward result of the model can be calculated, as shown in fig. 6. As can be seen from the above figure, the forward result is a cluster of concentric circles centered on the center of the anomaly. Although the gravity value of the center is assumed to be 0g.u (similarly at the core) according to the theoretical formula of the potential field, it is reasonably interpolated to 1.9g.u under the influence of the surrounding medium in order to match the actual situation. The subsequent graphs are all expressed by node coordinates, so that the position relation can be analyzed conveniently.
(3) And performing interpolation simulation after the model and forward modeling are completed. An appropriate number of points are artificially selected as interpolation source points, wherein the selected coordinates are ten in total and are expressed by node coordinates, and the node coordinates are as follows: [ 1815915148; 2244566999]THere, theThe reason why the selected nodes have more x coordinates consistent is that: among geological and geophysical problems to be solved, underground information is mostly obtained by drilling. This results in a striped distribution of data, which is chosen to simulate the form in which data is actually acquired.
(4) In order to better show the relevant characteristics of the radial basis function interpolation, numerical simulation of point interpolation is also carried out for comparison analysis, so that the common principle and advantages and disadvantages of the interpolation method are found out, and the method is improved. The expansion of the linear basis function p is:
(5) and calculating a point interpolation result and an error. The experimental interpolation results are shown in fig. 7, and the error is calculated by calculating the absolute value of the difference between the calculated value and the theoretical value. The error amount results are shown in fig. 8. It can be seen from fig. 8 that the result of the point interpolation still has a better depicting ability on the shape of the object, but there are two more serious defects. Firstly, the whole target body can deviate from the real target body position in two directions with large grid numbers (downwards and rightwards deviate after a plurality of experiments). The reason for the deviation is that it is inevitable that the error is large when the calculation is performed near a large number. The root cause is that the condition number of the linear basis function matrix p is too large, and is calculated to be 9179.3. The second is to interpolate a value violating the real physical law, such as a value of a (1, 9) position in the grid. And these problems tend to become particularly pronounced at the border locations of the grid.
(6) The shape function selected for the experiment radial basis interpolation is MQ function, and c is selected according to the empirical parameters and the distribution rule of the background grid
I.e. the product of the grid spacing 12.5 and the empirical parameter 1.789.
The expression for the radial basis function is:
in the formula: rT(X) is a radial basis function, the expression is as in equation 8, and b is a coefficient vector.
The matrix form of R is:
(7) and calculating a radial basis interpolation result and an error. The experimental interpolation results are shown in fig. 9, and the error is calculated by calculating the absolute value of the difference between the calculated value and the theoretical value. The error amount results are shown in fig. 10.
(8) As can be seen from the interpolation results of fig. 8: compared with a point difference method, the radial basis interpolation method has no phenomenon of target body deviation, can better interpolate a target body, and has better fitting performance on the shape. The radial basis interpolation method can also meet the constraint of physical conditions, namely the non-negativity of numerical values, and is a better numerical calculation and geophysical inversion method. The high-precision characteristic in the graph is just a key element required by the geophysical inversion. As can be seen from the radial basis gridless error matrix of fig. 10: the overall error can be controlled within a relatively small range. Due to the characteristics of the adopted shape function, the problem of serious boundary abrupt change does not exist, but the defect is also shown in the figure, and the local extreme value with larger error is in the area near the abnormal body of the model.
(9) In order to better show the results and prepare for the next inversion weight update, the model data of the earth surface and the results obtained by the two methods are shown and analyzed as shown in fig. 11. It can be seen from fig. 11 that the point interpolation has better performance in controlling the overall shape, but has a local deviation, which has a great influence on the accuracy of the geophysical inversion. The radial basis meshless method has excellent performance in interpolation result accuracy at a dense node and a middle node, but an error still exists in the boundary problem.
The invention has the following advantages. Firstly, the radial basis functions are very similar to the gravity value distribution generated by the geologic body unit, and the fitting effect is better according to the principle that the more similar the basis functions and the target functions are. The true distribution position and the abnormal value of the geologic body can be well restored by applying the radial basis function.
Compared with a polynomial-based point interpolation method, the point interpolation method has a larger condition number due to the selection of a coordinate system, so that the stability and the accuracy of the solution are poor. However, the radial basis method is only related to the norm of the position distance between the observation point and the source point, because it is not related to the selection of coordinates. Therefore, the condition number is not increased, and the stability and accuracy of the solution are affected.
Compared with a theoretical subdivision grid forward modeling method, the calculation speed of the radial basis function is greatly improved. Experiments in the previous section can prove that the operation steps are only about one fifth of that of a theoretical subdivision grid forward modeling method.
In addition, since no connection relation exists between the non-grid method point and the point, the region which is considered to have the abnormality can be artificially encrypted, so that the calculation accuracy of the region is higher.
Fitting a trigonometric function and a radial basis function to a differential equation: y ═ 3+ sinx, where x is 0 and y is 0; and solving the x ═ pi and y ═ 0, and summarizing the rule.
The specific implementation scheme is as follows:
(1) and solving to obtain an accurate solution. As a comparative study of the accuracy of the two methods. The exact solution is:
(2) and setting differential intervals for solving, wherein the intervals are designed to be equidistant and the intervals are pi/256.
(3) And selecting appropriate test functions, namely a trigonometric function (2) and a radial basis function (3).
Where n is the number of interpolation points, which is 2 and 6, respectively. I is the number of radial basic coordination points. r is the norm between the coordinates. c and theta are shape parameters. Here c is set to 5 and theta is set to 1.
(4) And selecting the position of the dispensing point. When n is 2 in the experiment, the distribution point is 5 × π/12+ (i-1) × pi/12. Where i takes 0 and 1. When n is 6, the position of the coordination point is selected to be pi/20 + pi/12 + (i-1) × pi/6. In selecting the position, care should be taken to choose symmetry within the domain of definition, which increases the condition number of the matrix. And a stable solution cannot be obtained.
(5) The function is brought into and solved for the differential residual equation (4), the well is solved for
(6) A comparison graph (figure 1) of two-point matching and six-point matching results of trigonometric function interpolation and a comparison graph (figure 2) of six-point trigonometric function matching and radial basis matching results are made, and the root mean square error of the two-point matching and the six-point trigonometric function matching and the radial basis matching is calculated. The formula is as follows:
where RMSD represents the root mean square error, xobsRepresenting an observed value, xrealRepresents a theoretical value.
The root mean square error of the two-point matching points of the trigonometric function is calculated to be 0.6127, the root mean square error of the six-point matching points of the trigonometric function is 0.0913, the root mean square error of the radial basic matching points is 0.0190, and the root mean square error of the radial basic matching points is the minimum.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the present invention in any way, and it will be apparent to those skilled in the art that the above description of the present invention can be applied to various modifications, equivalent variations or modifications without departing from the spirit and scope of the present invention.
Claims (8)
1. A gravity model forward modeling method based on MQ radial basis functions is characterized by comprising the following steps:
gridding the underground space of the working area;
adjusting parameters of the radial basis function according to the problems to be solved;
the method comprises the steps of substituting collected actual data or calculated model data into a function equation with a radial basis function as a shape function, and solving a weight coefficient vector, wherein the actual data comprises logging data;
and solving the information of the global grid point by using the weight coefficient vector and the points with known information to finally obtain the global gravity field distribution to finish forward modeling.
2. The MQ radial basis function-based gravity model forward modeling method according to claim 1, wherein the radial basis functions comprise: MQ radial basis functions.
3. The MQ radial basis function-based gravity model forward method according to claim 2, wherein the MQ radial basis functions have the form:
wherein r isIIn order to calculate the norm of a point and a node, c and theta need to be flexibly selected according to actual conditions.
4. The MQ radial basis function-based gravity model forward modeling method according to claim 1, wherein gridding the subsurface space of the working area comprises:
a series of Pi values are defined from surface observation data and well data.
5. The MQ radial basis function-based gravity model forward modeling method according to claim 4, wherein the background mesh nodes are calculated by affine transformation according to the following formula:
6. the MQ radial basis function-based gravity model forward modeling method according to claim 1, wherein the radial basis function is expressed as:
wherein R isT(X) is a radial basis function, and b is a coefficient vector.
7. The MQ radial basis function-based gravity model forward modeling method according to claim 1, wherein the matrix form of r (x) is as follows:
8. a MQ radial basis function-based gravity model forward modeling system, comprising:
one or more processors;
storage means for storing one or more programs;
when executed by the one or more processors, cause the one or more processors to implement the MQ radial basis function-based gravity model forward method according to any one of claims 1 to 7.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110665099.7A CN113255230B (en) | 2021-06-16 | 2021-06-16 | Gravity model forward modeling method and system based on MQ radial basis function |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110665099.7A CN113255230B (en) | 2021-06-16 | 2021-06-16 | Gravity model forward modeling method and system based on MQ radial basis function |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113255230A true CN113255230A (en) | 2021-08-13 |
| CN113255230B CN113255230B (en) | 2024-02-20 |
Family
ID=77188165
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110665099.7A Active CN113255230B (en) | 2021-06-16 | 2021-06-16 | Gravity model forward modeling method and system based on MQ radial basis function |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113255230B (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| NL2031632B1 (en) * | 2022-04-20 | 2023-11-07 | Chinese Acad Of Geological Sciences | Gravity inversion method and system based on meshfree method |
| US11815647B1 (en) * | 2022-04-20 | 2023-11-14 | Chinese Academy Of Geological Sciences | Gravity inversion method and system based on meshfree method |
| CN117237558A (en) * | 2023-11-10 | 2023-12-15 | 中南大学 | Fracture surface reconstruction method based on variational model and related equipment |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20030060981A1 (en) * | 1999-04-02 | 2003-03-27 | Conoco Inc. | Nonlinear constrained inversion method to determine base of salt interface from gravity and gravity tensor data |
| CN106767780A (en) * | 2016-11-28 | 2017-05-31 | 郑州轻工业学院 | Based on the extension ellipsoid set-membership filtering method that Chebyshev polynomial interopolations are approached |
| CN107391871A (en) * | 2017-08-03 | 2017-11-24 | 中国空气动力研究与发展中心计算空气动力研究所 | A kind of space lattice deformation method based on parallelization RBF |
| CN108205610A (en) * | 2018-01-10 | 2018-06-26 | 河海大学 | Concrete block design of Cooling System method based on quick exact numerical reconfiguration technique |
| CN108490496A (en) * | 2018-03-26 | 2018-09-04 | 中国石油化工股份有限公司 | Gravitational field inversion of Density method based on pseudo-radial basis function neural network |
| CN110443432A (en) * | 2019-08-14 | 2019-11-12 | 中国科学院武汉岩土力学研究所 | A kind of optimization algorithm solving Free Surface of Seepage based on radial basic point interpolation method |
-
2021
- 2021-06-16 CN CN202110665099.7A patent/CN113255230B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20030060981A1 (en) * | 1999-04-02 | 2003-03-27 | Conoco Inc. | Nonlinear constrained inversion method to determine base of salt interface from gravity and gravity tensor data |
| CN106767780A (en) * | 2016-11-28 | 2017-05-31 | 郑州轻工业学院 | Based on the extension ellipsoid set-membership filtering method that Chebyshev polynomial interopolations are approached |
| CN107391871A (en) * | 2017-08-03 | 2017-11-24 | 中国空气动力研究与发展中心计算空气动力研究所 | A kind of space lattice deformation method based on parallelization RBF |
| CN108205610A (en) * | 2018-01-10 | 2018-06-26 | 河海大学 | Concrete block design of Cooling System method based on quick exact numerical reconfiguration technique |
| CN108490496A (en) * | 2018-03-26 | 2018-09-04 | 中国石油化工股份有限公司 | Gravitational field inversion of Density method based on pseudo-radial basis function neural network |
| CN110443432A (en) * | 2019-08-14 | 2019-11-12 | 中国科学院武汉岩土力学研究所 | A kind of optimization algorithm solving Free Surface of Seepage based on radial basic point interpolation method |
Non-Patent Citations (2)
| Title |
|---|
| 李俊杰 等: "重力异常二维正演中的无网格方法", 煤田地质与勘探, vol. 46, no. 6, pages 181 - 186 * |
| 耿美霞;杨庆节;: "应用RBF神经网络反演二维重力密度分布", 石油地球物理勘探, no. 04, pages 651 - 657 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| NL2031632B1 (en) * | 2022-04-20 | 2023-11-07 | Chinese Acad Of Geological Sciences | Gravity inversion method and system based on meshfree method |
| US11815647B1 (en) * | 2022-04-20 | 2023-11-14 | Chinese Academy Of Geological Sciences | Gravity inversion method and system based on meshfree method |
| CN117237558A (en) * | 2023-11-10 | 2023-12-15 | 中南大学 | Fracture surface reconstruction method based on variational model and related equipment |
| CN117237558B (en) * | 2023-11-10 | 2024-02-13 | 中南大学 | Fracture surface reconstruction method based on variational model and related equipment |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113255230B (en) | 2024-02-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113255230A (en) | Gravity model forward modeling method and system based on MQ radial basis function | |
| CN105549106B (en) | A kind of gravity multiple solutions inversion method | |
| CN106199742B (en) | A kind of dimension of Frequency-domain AEM 2.5 band landform inversion method | |
| CN113553748B (en) | Three-dimensional magnetotelluric forward modeling numerical simulation method | |
| CN105701274A (en) | Generation method of three-dimensional local average random field samples of geotechnical parameters | |
| CN111967169B (en) | Two-degree body weight abnormal product decomposition numerical simulation method and device | |
| CN112687001A (en) | Three-dimensional geological structure model random generation and uncertainty analysis method | |
| CN113420487A (en) | Three-dimensional gravity gradient tensor numerical simulation method, device, equipment and medium | |
| CN115292973A (en) | Arbitrarily sampled space wave number domain three-dimensional magnetic field numerical simulation method and system | |
| CN108230452A (en) | A kind of model filling-up hole method based on textures synthesis | |
| CN115201927A (en) | Method for heavy magnetic three-dimensional joint inversion based on cross gradient constraint | |
| CN111597753A (en) | Data depth change characteristic self-adaptive two-dimensional resistivity inversion method and system | |
| CN114692523A (en) | Flow velocity prediction method of self-adaptive high-dimensional hydrodynamic equation based on graph convolution | |
| Liu et al. | Recovery of high frequency wave fields from phase space–based measurements | |
| CN111339688B (en) | Method for solving rocket simulation model time domain equation based on big data parallel algorithm | |
| CN116595827B (en) | Infinite dimension strip shot peening process planning method and system | |
| CN114970289B (en) | Three-dimensional magnetotelluric anisotropy forward modeling numerical simulation method, equipment and medium | |
| CN113779818B (en) | Three-dimensional geologic body electromagnetic field numerical simulation method, device, equipment and medium thereof | |
| CN113673163B (en) | Three-dimensional magnetic abnormal constant rapid forward modeling method, device and computer equipment | |
| CN115755199A (en) | Practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method | |
| CN114036806A (en) | Three-dimensional geothermal field numerical simulation method based on thermal conductivity anisotropic medium | |
| Ashby et al. | Modeling groundwater flow on MPPs | |
| Khan | Inverse problem in ground water: model development | |
| Wang et al. | Hermite radial basis collocation method for unsaturated soil water movement equation | |
| e Sousa et al. | Error estimates and adaptive procedures in geotechnical problems |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |