基于MatlabGUI的多球體重力正演

張博漢，張懷榜

?

基于MatlabGUI的多球體重力正演

張博漢1，張懷榜2，3

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

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

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

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

1 引 言

正演是重力勘探的重要手段，要對重力資料進行正確的地質解釋，首先必須搞清楚已知地質體的異常特征、數值大小及其分布和變化規(guī)律，也就是要先解決重力勘探的正演問題，一些復雜的地質現象可以視為多個簡單地質體的綜合效應，將復雜問題簡單化，分解復雜問題為多個簡單問題是解決重力問題的重要思路與方法。在實際工作中，一些等軸形狀的地質體，如礦巢、礦囊、巖株和其他的等軸的塊狀體都可以近似當作球體來研究，特別當球體尺度與埋深的比值較小時效果更好，而基于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)可以知道單個球體重力異?；咎卣鳎?/p>

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自動生成，不需要進行編輯，用戶只編輯控件按鈕對應的子函數即可，在每個控件子函數后續(xù)行中編寫運行程序，圖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 模擬結果的正確性驗證

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

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語言模擬相比，具有結果顯示更直觀、視覺效果更好、運行修改更簡潔方便的優(yōu)勢，更適合重力問題的研究。

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

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

[2]陳義群，陳華.基于MATLAB的工程物探軟件快速開發(fā)[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]張勝業(yè)，潘玉玲.應用地球物理原理[M].武漢：中國地質大學出版社，2004.

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

[9]羅華飛.MATLABGUI程序設計學習手記[M].北京：北京航空航天大學，2011.

[10]趙麗萍，單波，丁曉英.基于MATLAB的標準靜力觸探數據的輸入和輸出[J].工程地球物理學報，2014，11(1)，101-105.

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

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

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

[14]童孝忠，柳建新.MATLAB程序設計及在地球物理中的應用[M].長沙：中南大學出版社，2013.

[15]飛思科技產品研發(fā)中心.MATLAB基礎與提高[M].北京：電子工業(yè)出版社，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

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

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

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

P631.1

A

2016-05-31