CN113255230A - 基于mq径向基函数的重力模型正演方法及系统 - Google Patents
基于mq径向基函数的重力模型正演方法及系统 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
本发明提供了一种基于MQ径向基函数的重力模型正演方法及系统。该方法包括:将工作区域的地下空间进行网格化;根据需要解决的问题,调整好径向基函数的参数;将采集的实际数据(包含井中数据)或计算的模型数据带入由径向基函数作为形函数的函数方程中,并求解出权系数向量;用权系数向量及已知信息的点,求解出全局的网格点的信息,最终得到全局的重力场分布,完成正演。本发明提供的基于MQ径向基函数的重力模型正演方法及系统采用分步计算方法,缩短了正演计算时间,又能较好地满足精度要求,利用径向基函数还能减少正演参数的数量。
Description
技术领域
本发明涉及重力模型正演技术领域,特别是涉及一种基于MQ径向基函数的重力模型正演方法及系统。
背景技术
建立重力模型数值方法有两大类:基于网格的方法和无网格的方法。基于网格的方法是过去传统的方法,包括差分法和有限元法等,已经广泛应用于许多学科和工程。差分法的网格越密,节点越多,求解的精度就越高。由于其采用的是均匀大小的正交网格,在求解域边界形状比较复杂时,差分法求解的精度就会受到限制,甚至不能求解。有限元法在处理需要不断进行网格重构或者大变形等问题的时候,有限元网格可能会产生严重扭曲,不仅需要网格重构,而且严重地影响了解的精度[张雄,刘岩,马上.无网格法的理论及应用[J].力学进展,2009,39(01):1-36.]。无网格法是一种基于节点的数值方法,它避免了繁琐的网格剖分,为数值求解问题提供了新方法[李俊杰,严家斌.无网格法进展及在地球物理学中的应用[J].地球物理学进展,2014,29(01):452-461.],在多个领域都表现出优越性。
径向基函数法属于无网格法,它所拥有的无限光滑加上函数关于某个坐标方向的导数,可以表示为沿该方向所有离散点处函数值的加权线性组合[汪芳宗,廖小兵,谢雄.微分求积法的特性及其改进[J].计算力学学报,2015,32(06):765-771.],它在解决各类偏微分方程(如:椭圆型、抛物型和双曲型)问题具有较好的效果。常用的径向基函数方法有:局部径向基函数法[C.K.Lee,X.Liu,S.C.Fan.Local multiquadric approximation forsolving boundary value problems[J].Computational Mechanics,2003,30(5-6).]、基于微分求积的径向基函数法[C.Shu,H.Ding,H.Q.Chen,T.G.Wang.An upwind local RBF-DQ method for simulation of inviscid compressible flows[J].Computer Methodsin Applied Mechanics and Engineering,2004,194(18).]、埃尔米特径向基点插值法[Xin Liu,G.R.Liu,Kang Tai,K.Y.Lam.Radial point interpolation collocationmethod(RPICM)for the solution of nonlinear poisson problems[J].ComputationalMechanics,2005,36(4).]。
径向基函数的参数选择是一个伴随径向基发展一直存在的难点,PankajK.Mishra等人详细论述了几种无限光滑的形函数(高斯函数、MQ函数及逆MQ函数)构成的系数矩阵求解泊松问题的病态性研究,揭露了条件数、均方根误差与形函数形状参数的关系,为径向基函数形状参数的选择提供了依据[Pankaj K.Mishra,Sankar K.Nath.On theConvergence of Iterative Methods and Pseudoinverse Approaches in GlobalMeshless Collocation[J].International Journal of Applied and ComputationalMathematics,2017,3(4).]。熊正超指出了径向基函数数值计算不稳定性,最优形状参数的选择和核函数的方差有着密切的联系[熊正超.径向基函数逼近中的若干问题研究[D].复旦大学,2007.]。Ahmad Golbabai详述了常形状参数和变形状参数的区别,并建设性地提出了二元形状参数(BSP)和混合形状参数(HSP)两种变形状参数的方法,都有较好的数组计算表现[Ahmad Golbabai,Ehsan Mohebianfar,Hamed Rabiei.On the new variable shapeparameter strategies for radial basis functions[J].Computational and AppliedMathematics,2015,34(2).]。
支持域的选择方面,形函数支持域的选择对求解数值精度有很大的影响。常用的选择全局支持域的方法在求解大型问题时往往会加重其运算成本,局部的圆域支撑又面临着系数难以选择和精度下降的困扰,其中较好的改进是在求解对流占优问题时采用的迎风偏移局部支撑域[王晓东,欧阳洁,王玉龙,蒋涛.积分节点迎风偏移的节点积分无单元Galerkin方法[J].计算物理,2012,29(02):183-190.]和基于卷积神经网络对无网格形函数影响域的优化方法[刘宇翔,王东东,樊礼恒,陈健,侯松阳.基于卷积神经网络的无网格形函数影响域优化研究[J/OL].固体力学学报:1-19[2021-04-27].https://doi.org/10.19636/j.cnki.cjsm42-1250/o3.2021.003.]。
将非均匀配点格式与MQ径向基函数相结合用于数值模拟偏微分方程边值问题,得到了令人满意的结果[朱自强,曾思红,鲁光银,严文婕.二度体的重力张量有限元正演模拟[J].物探与化探,2010,34(05):668-671+685.],特别是对椭圆型偏微分方程的求解,为地球物理问题使用径向基函数求解提供了理论支撑,因为位场方程为椭圆型的偏微分方程。
径向基函数乃至无网格方法在重磁勘探及所有地球物理方法的应用都是相当少的,但与之对应的有限元方法却有着长足的发展[金钢燮,胡祥云,超敬来,南明根,康正男,金敬七,刘慧.复杂形体重磁异常的等参数有限元积分算法研究[J].石油地球物理勘探,2009,44(02):231-239+129+255.][朱自强,曾思红,鲁光银,严文婕.二度体的重力张量有限元正演模拟[J].物探与化探,2010,34(05):668-671+685.][赵丹丹,杜坚,郭智勇,徐伟,刘忠祥.埋地平行铁质管线磁异常模拟与探测识别[J].地下空间与工程学报,2020,16(03):891-896+932.]。传统的有限元法在很多的实际问题处理时往往难以取得良好的结果,比如:在地球物理高维的反演数值求解中,几乎是不可能构造离散点模型,这使得无网格方法有了发展的机会。
王月英用无单元Galerkin法进行了地震波正演模拟[王月英.地震波正演模拟中无单元Galerkin法初探[J].地球物理学进展,2007(05):1539-1544.]。冯德山等人基于无单元Galerkin法做出了探地雷达的正演模拟[冯德山,王洪华,戴前伟.基于无单元Galerkin法探地雷达正演模拟[J].地球物理学报,2013,56(01):298-308.]。孔倩等人使用Galerkin方法对二维位场进行向上延拓[孔倩,李鹏,李晶.二维位场延拓的无网格方法研究[J].物探化探计算技术,2017,39(02):155-160.]。李俊杰等人基于径向基函数进行了重力异常的二维正演[李俊杰,严家斌.重力异常二维正演中的无网格方法[J].煤田地质与勘探,2018,46(06):181-186.]以及大地电磁二维变分问题的无网格解法[李俊杰,严家斌.大地电磁二维正演中的有限元-径向基点插值法[J].中国有色金属学报,2015,25(05):1314-1324.]。
总而言之,径向基函数方法在数值理论上已经有了较为完善的研究,但在地球物理上的应用还较少。鉴于其优良的数值模拟能力,相信它在地球物理领域定有较大的发展潜力。
重力正演最常用和经典的方法就是利用万有引力计算公式:
对地质体所在的地下空间做二维及三维的地质剖分,然后每一个单元都用公式(1)来计算重力值再累加在一起,就完成了各点的重力值正演。再在观测区域逐点计算得到区域所有位置点的重力值,如此完成重力的正演计算。该逐点计算方法虽然可以得到工作区准确的重力值,但是却极其耗时,不利于在大规模计算时反复正反演的时间要求,精度控制也较为繁琐。
采用点插值径向基函数求取测区所有观测点重力值方法,其实施过程为:选取一个多项式函数作为线性基函数,如PT(x,y)={1,x,y,x2,xy,y2,x3,x2y,y3},计算出少量的随机位置点的重力值,然后根据公式c=P-1U算出系数向量c,此处U是一系列对应P(x,y)位置点的已求出的重力值。再根据下面的公式计算,
U(X,Y)=PT(X,Y)P-1U (2)
其中U(X,Y)为未知点的重力值,PT(X,Y)为对应的基函数矩阵。
这种点插值的方法相较于传统的网格方法而言,可以用较少的运算量就得到我们所需要的重力值的分布。在正反演过程中加上点插值这一步不仅不会加大运算量,反而可以节省运算时间,而且可以人为的选择离散点的位置和数量做多尺度的分析及局部分析,从而减少了反演参数的数量和重力反演的病态程度。
用公式(1)直接正演需要消耗大量的时间,对应的反演过程中所构成的系数矩阵也是极其庞大。求逆的过程也是十分的耗时,有时甚至难以求逆,同时由于反演的参数过多也会增加求解难度。
使用了点插值函数后虽然可以减少其求解参数的数量,且求解的时间缩短,但经过实验后发现其精确度较低。特别是由于基函数的计算方法是幂函数,使得在网格距离原点位置较远处的值会很大,导致了稳定性较差。而实际的地下空间是不会因为所选择的坐标原点的位置而产生不同,其值只会与地质体和观测点的位置距离有关。总之,纯粹的点插值方法精度有限。
发明内容
本发明要解决的技术问题是提供一种基于MQ径向基函数的重力模型正演方法及系统,采用分步计算方法,缩短了正演计算时间,又能较好地满足精度要求,利用径向基函数减少了正演参数的数量。
为解决上述技术问题,本发明提供了一种基于MQ径向基函数的重力模型正演方法,所述方法包括:将工作区域的地下空间进行网格化;根据需要解决的问题,调整好径向基函数的参数;将采集的实际数据或计算的模型数据带入由径向基函数作为形函数的函数方程中,并求解出权系数向量,其中,实际数据包括测井数据;用权系数向量及已知信息的点,求解出全局的网格点的信息,最终得到全局的重力场分布,完成正演。
在一些实施方式中,径向基函数包括:MQ径向基函数。
在一些实施方式中,MQ径向基函数具有如下形式:
其中,rI为计算点与节点的范数,c和θ需要根据实际情况灵活地选择参数。
在一些实施方式中,将工作区域的地下空间进行网格化,包括:由地面观测数据及井中数据定义为一系列Pi值。
在一些实施方式中,背景网格节点由如下公式进行仿射变换计算得出:
在一些实施方式中,径向基函数的表达式为:
其中,RT(X)为径向基函数,b为系数向量。
在一些实施方式中,R(X)的矩阵形式如下:
此外,本发明还提供了一种基于MQ径向基函数的重力模型正演系统,所述系统包括:一个或多个处理器;存储装置,用于存储一个或多个程序;当所述一个或多个程序被所述一个或多个处理器执行,使得所述一个或多个处理器实现根据前文所述的基于MQ径向基函数的重力模型正演方法。
采用这样的设计后,本发明至少具有以下优点:
本发明相较于传统的重力模型正演而言,采用分步计算方法,缩短了正演的计算时间,又能较好地满足精度的要求,利用径向基函数减少了反演参数的数量和重力反演的病态。它能在需要多次正反演的神经网络算法或其它的智能算法上节省运行成本;
同时相较于其它的插值算法,基于MQ的径向基函数还拥有与目标函数成分高度相似的优点,即MQ径向基函数与单位重力源产生的重力场高度相似。
附图说明
上述仅是本发明技术方案的概述,为了能够更清楚了解本发明的技术手段,以下结合附图与具体实施方式对本发明作进一步的详细说明。
图1a是c值不同对MQ函数形状的影响;
图1b是θ值不同对MQ函数形状的影响;
图2是理论工区布点图;
图3a是均匀节点布置方法的示意图;
图3b是Chebyshev节点布置方法的示意图;
图4是方法的执行过程;
图5是模型示意图;
图6是模型正演重力场图;
图7是点插值结果图;
图8是点插值法误差矩阵的示意图;
图9是径向基无网格法插值结果图;
图10是径向基无网格法误差矩阵的示意图;
图11是地表重力异常结果对比图;
图12是三角函数不同配点解结果对比图;
图13是六点三角函数配点与径向基函数配点解结果对比图。
具体实施方式
以下结合附图对本发明的优选实施例进行说明,应当理解,此处所描述的优选实施例仅用于说明和解释本发明,并不用于限定本发明。
本发明涉及了一种基于MQ径向基函数的重力数据正演的近似计算方法,本发明相较于传统的重力模型正演而言,采用分步计算方法,缩短了正演计算时间,又能较好地满足精度要求,利用径向基函数减少了正演参数的数量。它能在需要多次正反演的神经网络算法或其它智能算法上节省运行成本;同时相较于其它的插值算法,基于MQ的径向基函数还拥有与目标函数成分高度相似的优点,即MQ径向基函数与单位重力源产生的重力场高度相似,为之后的径向基函数反演提供了基础。
所提出的方法最主要的部分就是径向基函数的选择以及径向基函数参数的选择。为了较好地选择对比,我们做了观察实验和验证实验,然后用选取的合适的参数做正演流程。
一、径向基函数的参数选取
本次用到的径向基函数为Hardy提出的Multi-Quadric函数,应用的具体形式为:
式中rI为计算点与节点的距离范数。需要讨论的是形状参数c及θ的取值,合适的取值对径向基插值结果的准确性有较好的帮助。目前尚没有统一其最优值的选择方案及比对方法,c与θ的取值往往需要根据实际情况灵活地选择参数。而且选择参数往往是基于模拟结果的准确性决定,这导致了参数选择变得十分复杂且难以推广。
针对此问题我们对公式(3)做了相应的研究,结果如图1。从图1(a)中可以看出,c对MQ函数的影响体现在影响域的权值分配,当c较小时MQ函数较为尖锐,使得节点的影响域较小,造成对节点信息利用的不充分;当c过大时,虽然能使节点的影响域变大,但在实际情况中节点往往因为物理规律而并不具备很大的影响范围,又会增大计算量。这也是前人的研究中往往会对节点进行支持域划分的原因。在已经划分了支持域的前提下,取用较大的c值会造成严重的龙格现象,而得不到理想的结果。同时c值的大小类似于重力数据处理的向上延拓方法。对于θ的研究如图1(b)所示,从图中可以看出θ的选择对于函数的形状不会产生明显的影响。再加之对于重力的理论推导公式中的指数都有严格的物理理论,不宜对形函数中的指数项进行较大的调节。如需调节也应在对应理论公式指数的附近进行细微地调整,以满足整体的误差要求。
二、方法的实际应用模型
先介绍所用方法的现实应用环境,如图2所示。研究区域为ABB'A',在研究区域中AA’为地表观测面,现有由地面观测数据及井中数据得到的一系列Pi值,由这些Pi值再用网格点遍布研究区域,然后进行正反演计算。
网格点的布置提供两种方法,如图3所示。其中定义在求解区间的Chebyshev节点对于任意区间[a,b],可以用如下公式进行仿射变换计算得出,
运用Chebyshev节点的优势是能更好地拟合边界条件。
三、实际操作过程
实际操作过程参见图4所示。
建立的正演模型是一个中心位于点(50,-50),半径为10米,y轴无限延伸的水平圆柱体,圆柱体的剩余密度为2g/cm3。
根据剩余线密度公式λ=σ.S=σ.πR2,进而推导出计算区域各点的重力异常表达式:
式中△g为计算重力值,λ为圆柱体线密度,D为中轴线埋深,x为水平方向横坐标。模型如图5所示,中心为场源(水平无限延伸的圆柱体),背景节点剖分如图中的各小蓝点所示。
(2)根据剖分的节点及重力正演计算公式,可以计算出模型的正演结果,如图6所示。从上图中可以看出,正演结果为以异常体中心为圆心的一簇同心圆。虽然根据位场的理论公式,假设中心的重力值为0g.u(类似地核处),但为了和实际情况相匹配,受周围介质的影响将其合理插值为1.9g.u。之后的图都以节点坐标来表示,便于分析位置关系。
(3)模型及正演完成之后进行插值模拟。人为选择出适量的点作为插值的源点,在此选择的坐标共十个,表示方式是以节点坐标表示,它们是:[1815915148;2244566999]T,这里选择的节点有较多的x坐标一致的原因是:在需要解决的地质地球物理问题中,地下信息多采用钻井的方式取得。这就会使数据呈现条带式的分布,这种选择方法就是为了模拟实际采集数据的形式。
(4)为了更好地展示径向基函数插值的相关特征,同时也做了点插值的数值模拟用于对比分析,以找出插值法的共同原理和优缺点以及改进方法。线性基函数p的展开式为:
(5)计算点插值结果及误差。实验插值结果如图7,误差的计算方法是计算值与理论值之差的绝对值。误差量结果如图8所示。从图8中可以看出点插值的结果在目标体形状上还是有着较好的刻画能力,但是存在两个比较严重的缺陷。一是整个目标体会朝网格编号大的两个方向偏离真实的目标体位置(经过多次试验都是向下和右两个方向偏离)。偏离的原因是在大的编号附近进行计算时难免会使误差较大。根本原因是线性基函数矩阵p的条件数过大,经过计算得条件数为9179.3。二是插值出违反现实物理规律的值,比如网格中的(1,9)位置的值。而且这些问题往往会在网格的边界位置变得特别明显。
径向基函数的表达式为:
式中:RT(X)为径向基函数,表达式如公式8,b为系数向量。
R的矩阵形式为:
(7)计算出径向基插值结果及误差。实验插值结果如图9,误差的计算方法是计算值与理论值之差的绝对值。误差量结果如图10所示。
(8)从图8的插值结果可以看出:径向基插值法相较于点差值法而言不存在目标体偏离的现象,能够较好地插值出目标体,且形状也有较好的拟合性。径向基插值法同时也能满足物理条件的约束,即数值的非负性,是一种较好的数值计算和地球物理反演方法。图中的高精度特性正是地球物理反演需要的关键要素。从图10的径向基无网格法误差矩阵中可以看出:在整体的误差都能控制在相对较小的范围内。由于采用的形函数的特性,不存在严重的边界突变问题,但在图中同样也表现出了其不足的一面,从图中可以看出误差较大的地方极值在模型异常体附近地区。
(9)为了更好地展示结果和为下一步工作反演权值更新做准备,我们将地表的模型数据及两种方法得到的结果进行了展示及分析如图11。从图11中可以看出点插值在整体形状的把控上具有较好的表现,但局部基本都会有偏差,这对于地球物理反演的准确性造成较大的影响。而径向基无网格法在节点较密处及节点中部插值结果的准确性都有着优秀的表现,但在边界问题上依旧会存在误差。
本发明主要有以下几点优势。首先径向基函数与地质体单元产生的重力值分布十分相近,根据基函数与目标函数越相似其拟合效果越好的原则。运用径向基函数能较好地还原地质体的真实分布位置及异常值。
相较于基于多项式的点插值方法,点插值的方法由于坐标系的选择会有较大的条件数导致解的稳定性及精确性变差。但径向基的方法由于其与坐标的选择无关,只与观测点与源点的位置距离的范数有关。所以并不会使其条件数变大,进而影响解的稳定性与精确性。
径向基函数相较于理论的剖分网格正演方法,其计算速度具有较大的提升。上节中的实验可以证明,其运算步骤只有理论的剖分网格正演方法的五分之一左右。
此外,由于无网格方法点与点之间没有任何的连接关系,所以我们可以人为地加密我们认为存在异常的区域,使得该区域的计算精度更高。
用三角函数和径向基函数配点法对微分方程:y″=3+sinx在边界条件x=0,y=0;x=π,y=0进行求解,并总结其规律。
具体实施方案:
(1)求解出其精确解。作为两种方法准确性的对比研究。其精确解为:
(2)设置求解的差分间隔,间隔设计为等距,间距为π/256。
(3)选取合适的试函数,分别为三角函数(2)和径向基函数(3)。
其中n为插值点数,在此分别取2、6。I为径向基配点数。r为配点间的范数。c和θ为形状参数。在此设置为c为5和θ设置为1。
(4)选取配点的位置。在实验中选择为n为2时,配点位置为5*π/12+(i-1)*pi/12。其中i取0和1。n为6时,配点位置选为π/20+π/12+(i-1)*π/6。在选择位置时要注意别选择在定义域内对称,这样的选择会增大矩阵的条件数。进而无法得到稳定的解。
(5)将函数带入并求解微分的残值方程(4),井求解
(6)做出三角函数插值的两点配点与六点配点结果的对比图(图1),及六点三角函数配点和径向基配点的对比图(图2),并计算出它们的均方根误差。公式为:
其中RMSD代表均方根误差,xobs代表观测值,xreal代表理论值。
计算得出三角函数两点配点的均方根误差为0.6127,三角函数六点配点的均方根误差为0.0913,径向基配点的均方根误差为0.0190,径向基配点的均方根误差最小。
以上所述,仅是本发明的较佳实施例而已,并非对本发明作任何形式上的限制,本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰,均落在本发明的保护范围内。
Claims (8)
1.一种基于MQ径向基函数的重力模型正演方法,其特征在于,包括:
将工作区域的地下空间进行网格化;
根据需要解决的问题,调整好径向基函数的参数;
将采集的实际数据或计算的模型数据带入由径向基函数作为形函数的函数方程中,并求解出权系数向量,其中,实际数据包括测井数据;
用权系数向量及已知信息的点,求解出全局的网格点的信息,最终得到全局的重力场分布,完成正演。
2.根据权利要求1所述的基于MQ径向基函数的重力模型正演方法,其特征在于,径向基函数包括:MQ径向基函数。
4.根据权利要求1所述的基于MQ径向基函数的重力模型正演方法,其特征在于,将工作区域的地下空间进行网格化,包括:
由地面观测数据及井中数据定义为一系列Pi值。
8.一种基于MQ径向基函数的重力模型正演系统,其特征在于,包括:
一个或多个处理器;
存储装置,用于存储一个或多个程序;
当所述一个或多个程序被所述一个或多个处理器执行,使得所述一个或多个处理器实现根据权利要求1至7任意一项所述的基于MQ径向基函数的重力模型正演方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110665099.7A CN113255230B (zh) | 2021-06-16 | 2021-06-16 | 基于mq径向基函数的重力模型正演方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110665099.7A CN113255230B (zh) | 2021-06-16 | 2021-06-16 | 基于mq径向基函数的重力模型正演方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113255230A true CN113255230A (zh) | 2021-08-13 |
CN113255230B CN113255230B (zh) | 2024-02-20 |
Family
ID=77188165
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110665099.7A Active CN113255230B (zh) | 2021-06-16 | 2021-06-16 | 基于mq径向基函数的重力模型正演方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113255230B (zh) |
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 (zh) * | 2023-11-10 | 2023-12-15 | 中南大学 | 一种基于变分模型的断裂面重建方法及相关设备 |
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 (zh) * | 2016-11-28 | 2017-05-31 | 郑州轻工业学院 | 基于Chebyshev多项式插值逼近的扩展椭球集员滤波方法 |
CN107391871A (zh) * | 2017-08-03 | 2017-11-24 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于并行化径向基函数的空间网格变形方法 |
CN108205610A (zh) * | 2018-01-10 | 2018-06-26 | 河海大学 | 基于快速精确数值重构技术的混凝土块冷却系统设计方法 |
CN108490496A (zh) * | 2018-03-26 | 2018-09-04 | 中国石油化工股份有限公司 | 基于拟径向基函数神经网络的重力场密度反演方法 |
CN110443432A (zh) * | 2019-08-14 | 2019-11-12 | 中国科学院武汉岩土力学研究所 | 一种基于径向基点插值法求解渗流自由面的优化算法 |
-
2021
- 2021-06-16 CN CN202110665099.7A patent/CN113255230B/zh 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 (zh) * | 2016-11-28 | 2017-05-31 | 郑州轻工业学院 | 基于Chebyshev多项式插值逼近的扩展椭球集员滤波方法 |
CN107391871A (zh) * | 2017-08-03 | 2017-11-24 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于并行化径向基函数的空间网格变形方法 |
CN108205610A (zh) * | 2018-01-10 | 2018-06-26 | 河海大学 | 基于快速精确数值重构技术的混凝土块冷却系统设计方法 |
CN108490496A (zh) * | 2018-03-26 | 2018-09-04 | 中国石油化工股份有限公司 | 基于拟径向基函数神经网络的重力场密度反演方法 |
CN110443432A (zh) * | 2019-08-14 | 2019-11-12 | 中国科学院武汉岩土力学研究所 | 一种基于径向基点插值法求解渗流自由面的优化算法 |
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 (zh) * | 2023-11-10 | 2023-12-15 | 中南大学 | 一种基于变分模型的断裂面重建方法及相关设备 |
CN117237558B (zh) * | 2023-11-10 | 2024-02-13 | 中南大学 | 一种基于变分模型的断裂面重建方法及相关设备 |
Also Published As
Publication number | Publication date |
---|---|
CN113255230B (zh) | 2024-02-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113255230A (zh) | 基于mq径向基函数的重力模型正演方法及系统 | |
CN106199742B (zh) | 一种频率域航空电磁法2.5维带地形反演方法 | |
CN113553748B (zh) | 一种三维大地电磁正演数值模拟方法 | |
CN106683185B (zh) | 一种基于大数据的高精度曲面建模方法 | |
CN105701274A (zh) | 一种岩土参数三维局部平均随机场样本的生成方法 | |
CN113420487B (zh) | 一种三维重力梯度张量数值模拟方法、装置、设备和介质 | |
CN112687001B (zh) | 三维地质结构模型随机生成及不确定性分析方法 | |
CN111967169B (zh) | 二度体重力异常积分解数值模拟方法和装置 | |
CN115292973A (zh) | 一种任意采样的空间波数域三维磁场数值模拟方法及系统 | |
CN111597753A (zh) | 数据深度变化特征自适应的二维电阻率反演方法及系统 | |
CN114692523A (zh) | 基于图卷积的自适应高维流体动力方程的流速预测方法 | |
CN111339688B (zh) | 基于大数据并行算法求解火箭仿真模型时域方程的方法 | |
CN116595827B (zh) | 无限维度条带喷丸成形工艺规划方法和系统 | |
CN114970289B (zh) | 三维大地电磁各向异性正演数值模拟方法、设备及介质 | |
CN113779818B (zh) | 三维地质体其电磁场数值模拟方法、装置、设备及介质 | |
CN113673163B (zh) | 一种三维磁异常数快速正演方法、装置和计算机设备 | |
CN115755199A (zh) | 一种实用的非结构网格三维电磁反演平滑正则化方法 | |
CN108491671A (zh) | 一种土性参数六节点三角形随机场单元样本的生成方法 | |
CN114036806A (zh) | 基于热导率各向异性介质的三维地温场数值模拟方法 | |
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 | |
Li et al. | Fine numerical analysis of the crack-tip position for a Mumford–Shah minimizer | |
CN116911146B (zh) | 三维重力场全息数值模拟及cpu-gpu加速方法 |
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 |