張博漢，張懷榜

?

基于MatlabGUI的多球體重力正演

張博漢1，張懷榜2，3

( 1.長(zhǎng)江大學(xué) 地球物理與石油資源學(xué)院，湖北 武漢 430100；2.成都理工大學(xué) 地球物理學(xué)院，四川 成都 610059；3.中石化石油工程地球物理有限公司 勝利分公司，山東 東營(yíng) 257086)

為了更好地解決重力勘探中的正演、反演問(wèn)題，提高重力勘探資料的解釋精度，在MatlabGUI程序開發(fā)環(huán)境下對(duì)地球物理多個(gè)球體重力異常正演問(wèn)題進(jìn)行了研究與程序編寫，重點(diǎn)研究了一個(gè)均勻球體、兩個(gè)及兩個(gè)以上均勻球體模型的正演模擬，使一些較復(fù)雜的重力問(wèn)題能夠通過(guò)多球體模型近似模擬得以解決。該正演模擬程序界面簡(jiǎn)潔，操作方便，可交互輸入球體模型深度、密度、半徑等參數(shù)，模擬結(jié)果可進(jìn)行三維、二維等多方式顯示，三維曲面可任意調(diào)整觀測(cè)視角，二維曲線可任意選擇測(cè)線位置進(jìn)行顯示。本程序的模擬結(jié)果可以為復(fù)雜球體重力勘探資料解釋與反演提供參考，也是球體重力問(wèn)題正演研究的有效手段。

多球體模型；重力異常正演；MatlabGUI程序模擬

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

1 引 言

正演是重力勘探的重要手段，要對(duì)重力資料進(jìn)行正確的地質(zhì)解釋，首先必須搞清楚已知地質(zhì)體的異常特征、數(shù)值大小及其分布和變化規(guī)律，也就是要先解決重力勘探的正演問(wèn)題，一些復(fù)雜的地質(zhì)現(xiàn)象可以視為多個(gè)簡(jiǎn)單地質(zhì)體的綜合效應(yīng)，將復(fù)雜問(wèn)題簡(jiǎn)單化，分解復(fù)雜問(wèn)題為多個(gè)簡(jiǎn)單問(wèn)題是解決重力問(wèn)題的重要思路與方法。在實(shí)際工作中，一些等軸形狀的地質(zhì)體，如礦巢、礦囊、巖株和其他的等軸的塊狀體都可以近似當(dāng)作球體來(lái)研究，特別當(dāng)球體尺度與埋深的比值較小時(shí)效果更好，而基于Matlab的正演模擬特別適合球體重力問(wèn)題的研究[1-4]，因?yàn)樗幸粋€(gè)GUI(Graphical User Interface)可視化編程平臺(tái)，界面友好直觀、交互方便，而且可即時(shí)顯示，非常便于模擬結(jié)果的檢驗(yàn)和修改。因此，本文在研究重力異常模擬方法的基礎(chǔ)上編寫了MatlabGUI球體正演程序，從而為重力問(wèn)題的研究提供服務(wù)。

2 球體重力正演

球體正演模型可分為單個(gè)球體、兩個(gè)球體、多個(gè)球體模型等，球體越多模擬越復(fù)雜，結(jié)果也越接近實(shí)際情況，越能解決實(shí)際問(wèn)題[5-7]。下面分別進(jìn)行討論。

2.1 單個(gè)均勻球體計(jì)算

為使計(jì)算簡(jiǎn)化，將坐標(biāo)原點(diǎn)O選在球心在地面的投影點(diǎn)上，對(duì)于均勻球體，設(shè)球心埋藏深度為H(m)，半徑為R(m)，剩余密度為Δρ(g/cm3)，在其外部空間任一點(diǎn)所引起的重力異常等價(jià)于球心質(zhì)點(diǎn)的重力場(chǎng)的作用，這可使計(jì)算簡(jiǎn)化；對(duì)于非均勻球體，計(jì)算則較為復(fù)雜，采用多次迭代數(shù)值計(jì)算方法才能逐步逼近實(shí)際模型。

(1)

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

(2)

式(2)中：Δg(x,y)為布格重力異常(mGal)，G為萬(wàn)有引力常數(shù)，G=6.67×10-11m3/kg·s2。

根據(jù)坐標(biāo)點(diǎn)繪制的Δg重力異常是一個(gè)三維立體曲面。

當(dāng)沿一條線觀測(cè)時(shí)，式(2)得到簡(jiǎn)化，得到一條二維重力異常曲線。

分析式(2)可以知道單個(gè)球體重力異常基本特征：

1)重力異常最大值位于球心正上方(此時(shí)假設(shè)球心投影為原點(diǎn))，重力異常最大值為

(3)

在解反問(wèn)題時(shí)，此特征用來(lái)確定球心在地面的投影位置，式(3)還說(shuō)明重力異常與球體剩余質(zhì)量成正比，與埋深的平方成反比。

2)式(2)中因含有x,y，故異常值相對(duì)球心在地面投影點(diǎn)對(duì)稱分布，隨著x,y的無(wú)限增大，趨于0，由球心投影點(diǎn)向外其異常值的導(dǎo)數(shù)由大逐漸變小，異常值等值線為一系列由密到疏的同心圓。

(4)

(5)

在解正問(wèn)題時(shí)，式(5)是繪制Δg曲線的簡(jiǎn)便方法；解反問(wèn)題時(shí)，該式則可用來(lái)推算球心深度。

用式(2)還可以進(jìn)一步計(jì)算重力位的二次偏微商Vxz、Vxy、Vzz，Vxz、Vxy、Vzz是重力異常Δg沿兩個(gè)水平方向和垂直方向的變化率，可以提高重力異常的分辨率：

(6)

(7)

(8)

2.2 兩個(gè)及兩個(gè)以上均勻球體的計(jì)算

兩個(gè)均勻球體模型，設(shè)兩球球心埋藏深度分別為H1、H2，半徑為R1、R2,剩余密度為Δρ1、Δρ2，球體中心連線中點(diǎn)在地面的投影為坐標(biāo)原點(diǎn)O，則在地面任一點(diǎn)P(x,y,0)的重力異常值可視為兩球各自在地面重力作用的疊加，計(jì)算公式為

Δg(x,y)=

(9)

式中：(x10,y10)、(x20,y20)分別為兩球的球心在地面投影點(diǎn)的坐標(biāo)。

二次偏微商的計(jì)算公式如下：

Vxz=

(10)

Vxy=

(11)

Vzz=

(12)

兩個(gè)以上球體的模型，可通過(guò)設(shè)置合適的坐標(biāo)原點(diǎn)O進(jìn)行類似計(jì)算。

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

Matlab是美國(guó)MathWorks公司于1984年推出的計(jì)算繪圖功能強(qiáng)大的高級(jí)程序編寫平臺(tái)，是Matrix Laboratory的縮寫，是進(jìn)行矩陣運(yùn)算、符號(hào)運(yùn)算、數(shù)字信號(hào)處理、數(shù)學(xué)物理分析模擬的有利工具，目前應(yīng)用領(lǐng)域比較廣泛。GUI是Matlab圖形交互的一個(gè)程序編寫界面，在該界面下編寫程序，界面直觀，操作簡(jiǎn)便，方便檢查和修改，程序編制運(yùn)算效率較高，下面具體說(shuō)明程序編制主要過(guò)程。

3.1 重力異常正演程序界面布設(shè)

以Matlab2014a為例，在新建菜單“圖形用戶界面”中選擇新建GUI，打開可視化界面編制面板(圖1)，選用Edit Text、Static Text、PushButton、ToggleButton等控件設(shè)計(jì)重力異常正演程序界面，此程序中主要有15個(gè)控件，設(shè)置每個(gè)控件屬性(圖2)，設(shè)置完成保存程序?yàn)锽allGravityModelling.fig。

圖1 GUI程序設(shè)計(jì)面板Fig.1 GUI programming interface

圖2 控件屬性設(shè)置界面Fig.2 Control property inspector interface

3.2 正演程序編寫

在面板界面上方點(diǎn)擊編輯器按鈕(或右鍵點(diǎn)擊按鈕“查看回調(diào)”下的“CallBack”)進(jìn)入回調(diào)函數(shù)編輯界面，回調(diào)主函數(shù)為

function varargout = oneballgravity(varargin)

主函數(shù)由Matlab自動(dòng)生成，不需要進(jìn)行編輯，用戶只編輯控件按鈕對(duì)應(yīng)的子函數(shù)即可，在每個(gè)控件子函數(shù)后續(xù)行中編寫運(yùn)行程序，圖3是計(jì)算二次偏微商Vzz的子函數(shù)程序編寫示例。主函數(shù)中varargin是子函數(shù)調(diào)用參數(shù)，用戶點(diǎn)擊不同的操作按鈕，varargin會(huì)產(chǎn)生不同的參數(shù)值，主函數(shù)根據(jù)varargin參數(shù)調(diào)用不同控件子函數(shù)正演計(jì)算。

3.3 程序運(yùn)行及結(jié)果顯示

MatlabGUI界面、操作程序運(yùn)行及函數(shù)調(diào)用過(guò)程：

1)在GUI面板上方點(diǎn)擊運(yùn)行按鈕(或用快捷鍵Ctrl+T)，運(yùn)行程序，出現(xiàn)圖4操作界面；

圖3 二次偏微商Vzz函數(shù)子程序Fig.3 The sub-rountine of second partial derivative Vzz

2)填入模型參數(shù)，如球體半徑為50 m，埋深230 m，密度為2 g/cm3等，然后點(diǎn)擊要計(jì)算的按鈕，如Vxy按鈕；

3)主函數(shù)oneballgravity (varargin)根據(jù)varargin返回值調(diào)用與所點(diǎn)擊控件按鈕相對(duì)應(yīng)的子函數(shù)計(jì)算重力異常；

4)根據(jù)用戶需要和編寫程序輸出計(jì)算結(jié)果。圖5是一個(gè)均勻球體正演二次偏微商Vxy的結(jié)果，圖6是兩個(gè)均勻球體正演二次偏微商Vxz的結(jié)果。Vxy表示重力異常水平分量gx在y方向的變化率，單位為mGal/m，Vxz表示重力異常水平分量gx在z方向的變化率，單位為mGal/m。

計(jì)算結(jié)果可以看出，本程序不但可以計(jì)算1個(gè)球體、2個(gè)球體、多個(gè)球體及不同半徑球體模型的正演結(jié)果，而且可以將計(jì)算結(jié)果在3D空間顯示，也可以繪制一條二維重力異常曲線(圖7)，顯示角度可以任意調(diào)整。本重力異常正演程序?yàn)榻鉀Q一些復(fù)雜重力問(wèn)題提供了有力的正演分析工具，可有效解決由于球體模型數(shù)量少，正演結(jié)果與實(shí)際誤差較大的問(wèn)題。

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

圖5 一個(gè)球體Vxy計(jì)算結(jié)果Fig.5 One-ball Vxy modelling result

圖6 兩個(gè)球體Vxz計(jì)算結(jié)果Fig.6 Two-ball Vxz modelling result

圖7 二維重力異常Δg曲線Fig.7 2D gravity anomaly Δg curve

3.4 模擬結(jié)果的正確性驗(yàn)證

計(jì)算結(jié)果的正確性驗(yàn)證是正演模擬的重要環(huán)節(jié)，本程序設(shè)計(jì)了兩種檢驗(yàn)方式確保模擬結(jié)果的正確性和可靠性：

1)基于重力公式的正向檢驗(yàn)方式，首先根據(jù)公式計(jì)算出模型中某個(gè)點(diǎn)的理論數(shù)值，然后將計(jì)算結(jié)果與模擬結(jié)果進(jìn)行比較，算出模擬誤差，誤差超過(guò)標(biāo)準(zhǔn)限視為模擬失敗。對(duì)本程序1個(gè)球體Δg的驗(yàn)證：當(dāng)x=500 m時(shí)，理論值是0.048 908 mGal，實(shí)際模擬結(jié)果是0.048 908 mGal，誤差為0，模擬結(jié)果是正確的。

2)基于反演的驗(yàn)證方式，根據(jù)模擬數(shù)據(jù)反演球體的中心點(diǎn)位置、球心埋深、球體剩余質(zhì)量等參數(shù)，將反演結(jié)果與實(shí)際模型進(jìn)行對(duì)比，誤差限超標(biāo)視為模擬失敗。對(duì)本程序1個(gè)球體Δg反演①模擬數(shù)據(jù)Δg最大值位于x=0點(diǎn)，與理論值符合。②模擬數(shù)據(jù)Δgmax=0.954 736。

根據(jù)式(5)反算出球體的埋深為200.35 m，與模型埋深200 m的誤差為0.35 m，此誤差不是模擬引起的，是由于觀測(cè)數(shù)據(jù)點(diǎn)不夠密所致。

通過(guò)球體剩余質(zhì)量等參數(shù)驗(yàn)證也得到相似的結(jié)果，綜合認(rèn)為：模擬是正確的。

4 結(jié)論與建議

MatlabGUI可視化界面平臺(tái)上編制了一個(gè)均勻球體、兩個(gè)及兩個(gè)以上均勻球體的重力異常正演模擬程序，通過(guò)模型測(cè)試與實(shí)際應(yīng)用，得到以下認(rèn)識(shí)：

1)基于MatlabGUI的多球體重力正演模擬操作簡(jiǎn)便、界面友好，可進(jìn)行三維、二維不同的觀測(cè)方式進(jìn)行模擬計(jì)算，結(jié)果可以選擇不同觀測(cè)面、不同觀測(cè)線進(jìn)行多維度顯示，為多角度重力異常問(wèn)題研究提供了有效手段。

2)MatlabGUI重力異常程序正演與Visual C++、Fortran、QT語(yǔ)言模擬相比，具有結(jié)果顯示更直觀、視覺效果更好、運(yùn)行修改更簡(jiǎn)潔方便的優(yōu)勢(shì)，更適合重力問(wèn)題的研究。

下一步，筆者所在的研究團(tuán)隊(duì)將開展均勻水平圓柱體、鉛錘臺(tái)階及水平物質(zhì)半平面、傾斜脈等規(guī)則模型的重力異常正演問(wèn)題研究及程序編寫，為重力正演、反演問(wèn)題提供更強(qiáng)的技術(shù)支撐。

[1]張劍，師學(xué)明，劉夢(mèng)花.基于MATLAB開發(fā)環(huán)境的球體重力正演[J].工程地球物理學(xué)報(bào)，2007，4(5)：460-464.

[2]陳義群，陳華.基于MATLAB的工程物探軟件快速開發(fā)[J].地球物理學(xué)進(jìn)展，2004，19(4)：802-806.

[3]徐佳，朱魯，翟培合.基于Voxler平臺(tái)的電法數(shù)據(jù)三維可視化[J].工程地球物理學(xué)報(bào)，2014，11(6)，772-775.

[4]夏媛媛，趙民，藏歌,等.正演模擬技術(shù)在解釋反演中的應(yīng)用[J].工程地球物理學(xué)報(bào)，2014，11(6)，842-846.

[5]長(zhǎng)春地質(zhì)學(xué)院重力教研室.重力勘探[M].北京：地質(zhì)出版社，1980.

[6]陳善.重力勘探[M].北京：地質(zhì)出版社，1988.

[7]張勝業(yè)，潘玉玲.應(yīng)用地球物理原理[M].武漢：中國(guó)地質(zhì)大學(xué)出版社，2004.

[8]陳垚光，毛濤濤，王正林.精通MATLABGUI設(shè)計(jì)(第三版)[M].北京：電子工業(yè)出版社，2013.

[9]羅華飛.MATLABGUI程序設(shè)計(jì)學(xué)習(xí)手記[M].北京：北京航空航天大學(xué)，2011.

[10]趙麗萍，單波，丁曉英.基于MATLAB的標(biāo)準(zhǔn)靜力觸探數(shù)據(jù)的輸入和輸出[J].工程地球物理學(xué)報(bào)，2014，11(1)，101-105.

[11]宋葉志，賈東永.MATLAB數(shù)值分析與應(yīng)用[M].北京：機(jī)械工業(yè)出版社，2009.

[12]張志涌，劉瑞楨，楊程櫻.掌握和精通MATLAB[M].北京：北京航空航天大學(xué)出版社，1997.

[13]劉衛(wèi)國(guó).MATLAB程序設(shè)計(jì)教程[M].北京：中國(guó)水利水電出版社，2005.

[14]童孝忠，柳建新.MATLAB程序設(shè)計(jì)及在地球物理中的應(yīng)用[M].長(zhǎng)沙：中南大學(xué)出版社，2013.

[15]飛思科技產(chǎn)品研發(fā)中心.MATLAB基礎(chǔ)與提高[M].北京：電子工業(yè)出版社，2006.

[16]胡飛，石瑞平，陳建國(guó).MAPGIS在圖元的物理重排問(wèn)題[J].工程地球物理學(xué)報(bào)，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

國(guó)家重大專項(xiàng)大型油氣田及煤成氣開發(fā)(編號(hào)：2016ZX05005005)；長(zhǎng)江大學(xué)大學(xué)生創(chuàng)新創(chuàng)業(yè)基金

張博漢(1995-)，男，本科學(xué)生，主要學(xué)習(xí)和研究方向是計(jì)算地球物理學(xué)。E-mail：18086457768@163.com

張懷榜(1966-)，男，高級(jí)工程師，博士研究生，主要從事石油地球物理勘探方法研究。E-mail：zbhzhb@163.com

P631.1

A

2016-05-31