基于MatlabGUI的多球体重力正演

张博汉，张怀榜



基于MatlabGUI的多球体重力正演

张博汉1，张怀榜2，3

( 1.长江大学 地球物理与石油资源学院，湖北 武汉 430100；2.成都理工大学 地球物理学院，四川 成都 610059；3.中石化石油工程地球物理有限公司 胜利分公司，山东 东营 257086)

为了更好地解决重力勘探中的正演、反演问题，提高重力勘探资料的解释精度，在MatlabGUI程序开发环境下对地球物理多个球体重力异常正演问题进行了研究与程序编写，重点研究了一个均匀球体、两个及两个以上均匀球体模型的正演模拟，使一些较复杂的重力问题能够通过多球体模型近似模拟得以解决。该正演模拟程序界面简洁，操作方便，可交互输入球体模型深度、密度、半径等参数，模拟结果可进行三维、二维等多方式显示，三维曲面可任意调整观测视角，二维曲线可任意选择测线位置进行显示。本程序的模拟结果可以为复杂球体重力勘探资料解释与反演提供参考，也是球体重力问题正演研究的有效手段。

多球体模型；重力异常正演；MatlabGUI程序模拟

( 1.GeophysicsandOilResourcesInstitute,YangtzeUniversity,WuhanHubei430100,China;2.GeophysicsInstitute,ChengduUniversityofTechnology,ChengduSichuang610059,China;3.ShengliBranch,SINOPECGeophysicalCorporation,DongyingShandong257086,China)

1 引 言

正演是重力勘探的重要手段，要对重力资料进行正确的地质解释，首先必须搞清楚已知地质体的异常特征、数值大小及其分布和变化规律，也就是要先解决重力勘探的正演问题，一些复杂的地质现象可以视为多个简单地质体的综合效应，将复杂问题简单化，分解复杂问题为多个简单问题是解决重力问题的重要思路与方法。在实际工作中，一些等轴形状的地质体，如矿巢、矿囊、岩株和其他的等轴的块状体都可以近似当作球体来研究，特别当球体尺度与埋深的比值较小时效果更好，而基于Matlab的正演模拟特别适合球体重力问题的研究[1-4]，因为它有一个GUI(Graphical User Interface)可视化编程平台，界面友好直观、交互方便，而且可即时显示，非常便于模拟结果的检验和修改。因此，本文在研究重力异常模拟方法的基础上编写了MatlabGUI球体正演程序，从而为重力问题的研究提供服务。

2 球体重力正演

球体正演模型可分为单个球体、两个球体、多个球体模型等，球体越多模拟越复杂，结果也越接近实际情况，越能解决实际问题[5-7]。下面分别进行讨论。

2.1 单个均匀球体计算

为使计算简化，将坐标原点O选在球心在地面的投影点上，对于均匀球体，设球心埋藏深度为H(m)，半径为R(m)，剩余密度为Δρ(g/cm3)，在其外部空间任一点所引起的重力异常等价于球心质点的重力场的作用，这可使计算简化；对于非均匀球体，计算则较为复杂，采用多次迭代数值计算方法才能逐步逼近实际模型。

(1)

地表任意位置P(x,y,0)的重力异常为：

(2)

式(2)中：Δg(x,y)为布格重力异常(mGal)，G为万有引力常数，G=6.67×10-11m3/kg·s2。

根据坐标点绘制的Δg重力异常是一个三维立体曲面。

当沿一条线观测时，式(2)得到简化，得到一条二维重力异常曲线。

分析式(2)可以知道单个球体重力异常基本特征：

1)重力异常最大值位于球心正上方(此时假设球心投影为原点)，重力异常最大值为

(3)

在解反问题时，此特征用来确定球心在地面的投影位置，式(3)还说明重力异常与球体剩余质量成正比，与埋深的平方成反比。

2)式(2)中因含有x,y，故异常值相对球心在地面投影点对称分布，随着x,y的无限增大，趋于0，由球心投影点向外其异常值的导数由大逐渐变小，异常值等值线为一系列由密到疏的同心圆。

(4)

(5)

在解正问题时，式(5)是绘制Δg曲线的简便方法；解反问题时，该式则可用来推算球心深度。

用式(2)还可以进一步计算重力位的二次偏微商Vxz、Vxy、Vzz，Vxz、Vxy、Vzz是重力异常Δg沿两个水平方向和垂直方向的变化率，可以提高重力异常的分辨率：

(6)

(7)

(8)

2.2 两个及两个以上均匀球体的计算

两个均匀球体模型，设两球球心埋藏深度分别为H1、H2，半径为R1、R2,剩余密度为Δρ1、Δρ2，球体中心连线中点在地面的投影为坐标原点O，则在地面任一点P(x,y,0)的重力异常值可视为两球各自在地面重力作用的叠加，计算公式为

Δg(x,y)=

(9)

式中：(x10,y10)、(x20,y20)分别为两球的球心在地面投影点的坐标。

二次偏微商的计算公式如下：

Vxz=

(10)

Vxy=

(11)

Vzz=

(12)

两个以上球体的模型，可通过设置合适的坐标原点O进行类似计算。

3 MatlabGUI可视化正演程序[8-16]

Matlab是美国MathWorks公司于1984年推出的计算绘图功能强大的高级程序编写平台，是Matrix Laboratory的缩写，是进行矩阵运算、符号运算、数字信号处理、数学物理分析模拟的有利工具，目前应用领域比较广泛。GUI是Matlab图形交互的一个程序编写界面，在该界面下编写程序，界面直观，操作简便，方便检查和修改，程序编制运算效率较高，下面具体说明程序编制主要过程。

3.1 重力异常正演程序界面布设

以Matlab2014a为例，在新建菜单“图形用户界面”中选择新建GUI，打开可视化界面编制面板(图1)，选用Edit Text、Static Text、PushButton、ToggleButton等控件设计重力异常正演程序界面，此程序中主要有15个控件，设置每个控件属性(图2)，设置完成保存程序为BallGravityModelling.fig。

图1 GUI程序设计面板Fig.1 GUI programming interface

图2 控件属性设置界面Fig.2 Control property inspector interface

3.2 正演程序编写

在面板界面上方点击编辑器按钮(或右键点击按钮“查看回调”下的“CallBack”)进入回调函数编辑界面，回调主函数为

function varargout = oneballgravity(varargin)

主函数由Matlab自动生成，不需要进行编辑，用户只编辑控件按钮对应的子函数即可，在每个控件子函数后续行中编写运行程序，图3是计算二次偏微商Vzz的子函数程序编写示例。主函数中varargin是子函数调用参数，用户点击不同的操作按钮，varargin会产生不同的参数值，主函数根据varargin参数调用不同控件子函数正演计算。

3.3 程序运行及结果显示

MatlabGUI界面、操作程序运行及函数调用过程：

1)在GUI面板上方点击运行按钮(或用快捷键Ctrl+T)，运行程序，出现图4操作界面；

图3 二次偏微商Vzz函数子程序Fig.3 The sub-rountine of second partial derivative Vzz

2)填入模型参数，如球体半径为50 m，埋深230 m，密度为2 g/cm3等，然后点击要计算的按钮，如Vxy按钮；

3)主函数oneballgravity (varargin)根据varargin返回值调用与所点击控件按钮相对应的子函数计算重力异常；

4)根据用户需要和编写程序输出计算结果。图5是一个均匀球体正演二次偏微商Vxy的结果，图6是两个均匀球体正演二次偏微商Vxz的结果。Vxy表示重力异常水平分量gx在y方向的变化率，单位为mGal/m，Vxz表示重力异常水平分量gx在z方向的变化率，单位为mGal/m。

计算结果可以看出，本程序不但可以计算1个球体、2个球体、多个球体及不同半径球体模型的正演结果，而且可以将计算结果在3D空间显示，也可以绘制一条二维重力异常曲线(图7)，显示角度可以任意调整。本重力异常正演程序为解决一些复杂重力问题提供了有力的正演分析工具，可有效解决由于球体模型数量少，正演结果与实际误差较大的问题。

图4 程序操作界面Fig.4 Programming interface

图5 一个球体Vxy计算结果Fig.5 One-ball Vxy modelling result

图6 两个球体Vxz计算结果Fig.6 Two-ball Vxz modelling result

图7 二维重力异常Δg曲线Fig.7 2D gravity anomaly Δg curve

3.4 模拟结果的正确性验证

计算结果的正确性验证是正演模拟的重要环节，本程序设计了两种检验方式确保模拟结果的正确性和可靠性：

1)基于重力公式的正向检验方式，首先根据公式计算出模型中某个点的理论数值，然后将计算结果与模拟结果进行比较，算出模拟误差，误差超过标准限视为模拟失败。对本程序1个球体Δg的验证：当x=500 m时，理论值是0.048 908 mGal，实际模拟结果是0.048 908 mGal，误差为0，模拟结果是正确的。

2)基于反演的验证方式，根据模拟数据反演球体的中心点位置、球心埋深、球体剩余质量等参数，将反演结果与实际模型进行对比，误差限超标视为模拟失败。对本程序1个球体Δg反演①模拟数据Δg最大值位于x=0点，与理论值符合。②模拟数据Δgmax=0.954 736。

根据式(5)反算出球体的埋深为200.35 m，与模型埋深200 m的误差为0.35 m，此误差不是模拟引起的，是由于观测数据点不够密所致。

通过球体剩余质量等参数验证也得到相似的结果，综合认为：模拟是正确的。

4 结论与建议

MatlabGUI可视化界面平台上编制了一个均匀球体、两个及两个以上均匀球体的重力异常正演模拟程序，通过模型测试与实际应用，得到以下认识：

1)基于MatlabGUI的多球体重力正演模拟操作简便、界面友好，可进行三维、二维不同的观测方式进行模拟计算，结果可以选择不同观测面、不同观测线进行多维度显示，为多角度重力异常问题研究提供了有效手段。

2)MatlabGUI重力异常程序正演与Visual C++、Fortran、QT语言模拟相比，具有结果显示更直观、视觉效果更好、运行修改更简洁方便的优势，更适合重力问题的研究。

下一步，笔者所在的研究团队将开展均匀水平圆柱体、铅锤台阶及水平物质半平面、倾斜脉等规则模型的重力异常正演问题研究及程序编写，为重力正演、反演问题提供更强的技术支撑。

[1]张剑，师学明，刘梦花.基于MATLAB开发环境的球体重力正演[J].工程地球物理学报，2007，4(5)：460-464.

[2]陈义群，陈华.基于MATLAB的工程物探软件快速开发[J].地球物理学进展，2004，19(4)：802-806.

[3]徐佳，朱鲁，翟培合.基于Voxler平台的电法数据三维可视化[J].工程地球物理学报，2014，11(6)，772-775.

[4]夏媛媛，赵民，藏歌,等.正演模拟技术在解释反演中的应用[J].工程地球物理学报，2014，11(6)，842-846.

[5]长春地质学院重力教研室.重力勘探[M].北京：地质出版社，1980.

[6]陈善.重力勘探[M].北京：地质出版社，1988.

[7]张胜业，潘玉玲.应用地球物理原理[M].武汉：中国地质大学出版社，2004.

[8]陈垚光，毛涛涛，王正林.精通MATLABGUI设计(第三版)[M].北京：电子工业出版社，2013.

[9]罗华飞.MATLABGUI程序设计学习手记[M].北京：北京航空航天大学，2011.

[10]赵丽萍，单波，丁晓英.基于MATLAB的标准静力触探数据的输入和输出[J].工程地球物理学报，2014，11(1)，101-105.

[11]宋叶志，贾东永.MATLAB数值分析与应用[M].北京：机械工业出版社，2009.

[12]张志涌，刘瑞桢，杨程樱.掌握和精通MATLAB[M].北京：北京航空航天大学出版社，1997.

[13]刘卫国.MATLAB程序设计教程[M].北京：中国水利水电出版社，2005.

[14]童孝忠，柳建新.MATLAB程序设计及在地球物理中的应用[M].长沙：中南大学出版社，2013.

[15]飞思科技产品研发中心.MATLAB基础与提高[M].北京：电子工业出版社，2006.

[16]胡飞，石瑞平，陈建国.MAPGIS在图元的物理重排问题[J].工程地球物理学报，2005，2(3)：235-238.

On Multi-sphere Gravity Forward Modeling Based on MatlabGUI

Zhang Bohan1, Zhang Huaibang2,3

In order to solve the problems of forward and inverse gravity method, and for improving the precision of gravity data interpretation, this paper based on the MatlabGUI development and environment discussed gravity anomaly forward modeling of multiple spheres and how to write the processing MatlabGUI codes. The main researching results focused on the forward modeling problems of one homogeneous sphere, two homogeneous spheres and several homogeneous spheres. It makes some of the more complex gravity problems be solved by multi-sphere gravity simulating approximately. The forward simulation program has a simple interface, and its operation is convenient. The sphere model depth, media density, sphere radius and other parameters can be input interactively. The simulation results can be displayed in the three-dimensional space or in a curve. The 3D figures can be observed in any direction by adjusting the showing angle. A line gravity curve can be displayed conveniently by inputting the line parameters. The procedure simulation results can give us some guides or arouse us ideas for interpreting complex sphere gravity data and doing gravity inversion, The MatlabGUI gravity modellings are also wonderful tools for researching sphere gravity forward problems.

multi-sphere model; gravity anomaly forward modeling; MatlabGUI program simulation

1672—7940(2016)05—0580—06

10.3969/j.issn.1672-7940.2016.05.004

国家重大专项大型油气田及煤成气开发(编号：2016ZX05005005)；长江大学大学生创新创业基金

张博汉(1995-)，男，本科学生，主要学习和研究方向是计算地球物理学。E-mail：18086457768@163.com

张怀榜(1966-)，男，高级工程师，博士研究生，主要从事石油地球物理勘探方法研究。E-mail：zbhzhb@163.com

P631.1

A

2016-05-31