Matlab軟件包在地球物理數(shù)值計算與輔助教學中的應用

童孝忠，王 濤，柳建新

（1.中南大學 地球科學與信息物理學院，長沙 410083；2.有色資源與地質災害探查湖南省重點實驗室，長沙 410083）

0 前言

Matlab軟件包是美國MathWorks公司研發(fā)的一款軟件產(chǎn)品，是目前在國際上被廣泛接受和使用的數(shù)值計算軟件。對于數(shù)值計算問題，用其它程序設計語言編程求解非常麻煩，并且需要具備專門的數(shù)學知識及一定的程序設計技能，而用Matlab編程，往往只要少數(shù)幾個語句即可完成求解任務，具有編程效率高、使用方便等特點。Matlab軟件包這種強大的數(shù)值計算能力，使其成為地球物理科學計算方面的首選解題工具［1－5］。

Matlab軟件包為用戶開發(fā)圖形界面提供了一個方便、高效的集成環(huán)境—圖形用戶開發(fā)環(huán)境（Graphical User Interface Development Environment，GUIDE）。它將所有的GUI支持的用戶控件都集成在一起，并提供界面外觀、屬性和行為響應方式的設置方法。GUIDE將用戶設計好的GUI界面保存在一個FIG（圖形）資源文件中，同時還能夠生成包含GUI初始化和組件界面布局控制代碼的M文件，用戶可以采用這一框架編制自己的應用程序。將Matlab的圖形用戶界面設計（GUI）融入到課件制作中，能夠輕松實現(xiàn)地球物理課程的輔助教學［6－7］。

作者在本文針對地球物理數(shù)值計算和教學過程中的難點和重點，借助Matlab軟件包，給出具體的實例。

1 Matlab在數(shù)值計算中的應用

1.1 重力異常正演

重力正演問題就是給定地下某種地質體的形狀、產(chǎn)狀和剩余密度等，通過理論計算求取它在地面上或空間范圍內引起的異常大小、特征和變換規(guī)律等。這里只分析球體模型重力異常的正演計算。

設球心的埋藏深度為D，球的半徑為R，剩余密度為σ，則它的剩余質量為M＝。將坐標原點選在球心在地面的投影處，則球體的重力異常表達式為［8］：

式中 G ＝6.67×10－11m3／kg·s2＝6.67×10－2（g.u.）m2／t；D 的 單 位 為 m；M 的 單 位 為 t（噸）。

假設R＝50m，D＝100m，σ＝1t／m3，編制的Matlab程序代碼如下：

程序執(zhí)行后的結果如圖1所示。

1.2 磁異常正演

球體磁場的正演公式為［9］式（2）與式（3）。

式中 R為球體中心埋深；m為球體磁矩，且m＝MV（M為磁化強度，V 為球體體積）；I為磁化傾角；A′為觀測剖面與磁化強度水平投影夾角。

假設球體中心埋深R＝15m，半徑r＝10m，k＝0.015（SI），當?shù)卮艌鯞＝50 000nT。建立球體磁異常正演計算的Matlab函數(shù)文件MAG＿sphere＿FWD.m，程序如下：

function［Hax，Hay，Za，Delta＿T］＝MAG＿sphere＿FWD（A，I）

%磁性球體磁場正演，輸入?yún)?shù)A為觀測剖面與磁化強度水平投影夾角，I為磁化傾角

當A′＝45°，I＝0°時，在 MATLAB命令窗口調用函數(shù)文件

? ［Hax，Hay，Za，Delta＿T］＝MAG＿sphere＿FWD（pi／4，0）；

圖示正演計算的磁異常的結果如圖2所示。

1.3 電磁法正演計算

這里討論層狀介質的大地電磁測深正演計算。假定地電剖面是水平分層均勻的，如果地電剖面共有N 層，則共有2 N－1個參數(shù)，hi（i＝1，2，…，N－1），ρi（i＝1，2，…，N）。hi、ρi分別代表第i層的厚度和電阻率。對于上述的一維層狀介質模型，計算視電阻率ρa和相位φa的公式為式（4）［10］。

式中 Z1表示第一層地表的波阻抗；μ為導磁率；ω＝為角頻率；Z可用公式（5）的遞推公式計算。1

建立層狀介質大地電磁正演計算的Matlab函數(shù)文件MT1D＿FWD.m，程序如下：

當模型參數(shù)為ρ1＝100Ω·m ，ρ2＝600Ω·m，h1＝1 800m時，在 MATLAB命令窗口調用函數(shù)文件：

?［pc，ph］＝ MT1D＿FWD（［100，600］，1800）；

圖示正演計算結果如圖3所示。

2 Matlab在教學中的應用

將Matlab的圖形用戶界面設計融入到課件制作中，通過Matlab強大的數(shù)值計算功能和GUI的可視化界面，能夠輕松地實現(xiàn)地球物理課程的輔助教學。

下面給出磁異常向上延拓GUI制作的操作步驟，設計的GUI用戶界面如圖4所示。

（1）打開GUI設計窗口，添加有關控件對象。在MATLAB命令窗口輸入命令guide，打開GUI設計窗口。單擊GUI設計窗口控件工具欄中的Axes按鈕，并在圖形窗口中拖出一個矩形框，調整好大小和位置；再添加3個按鈕、2個靜態(tài)文本框和1個可編輯文本框，并調整好大小和位置。

（2）利用屬性編輯器，設置圖形對象的屬性。打開屬性編輯器，修改控件對象的屬性及屬性值。將3個按鈕的String屬性設置為Za、Ha和Delta＿T；將2個靜態(tài)文本框的String屬性設置為“球體磁異常延拓用戶界面”和“上延高度”。

（3）添加“模型”菜單項。雙擊圖上空白處，直接在屬性查看器中將menu和toolbar的屬性修改為figure，然后打開菜單編輯器，新建一個菜單項，它的Label屬性設為“模型”，在剛建的菜單項下建立兩個子菜單項，其Label屬性分別設為“單個球體”和“組合球體”。

（4）保存圖形用戶界面。選擇File菜單中的Save命令，將設計的圖形用戶界面保存為 Mag＿Extent.fig，同時生成 Mag＿Extent.m文件。

（5）編寫代碼，實現(xiàn)控件及菜單功能。在函數(shù)文件Mag＿Extent.m添加相應代碼，省去自帶的注釋部分，程序代碼如下：

圖3 一維模型的大地電磁測深響應Fig.3 Magnetotelluric responses of one－dimensional model

圖4 磁異常延拓結果圖Fig.4 The results of magnetic anomaly extension

磁異常計算函數(shù)ProlongateForward.m的程序代碼如下：

運行圖形用戶界面，若選擇組合球體模型，上延高度取為10m，即可得如圖4所示的計算結果。

3 結論

（1）Matlab軟件包在解決實際的數(shù)值計算問題中，具有使用更為簡便、語句功能更強的特點，適合用于計算地球物理的研究與應用。

（2）利用Matlab軟件包，通過具體應用實例，給出了直觀形象的圖形界面，豐富了教學內容，達到了一定的教學效果，為課程教學方法和手段的改革探索了新的思路。

［1］ALAN W.Geophysica：MATLAB－based softwarefor the simulation，display and processing of nearsurface geophysical data［J］.Computers ＆ Geosciences，2002，28（6）：751－762.

［2］馮治權.MATLAB及其在重力基點網(wǎng)平差中的應用［J］.物探化探計算技術，2003，25（4）：336－339.

［3］李曉昌.在MATLAB平臺上實現(xiàn)可控源音頻大地電磁反演數(shù)據(jù)三維可視化顯示［J］.物探化探計算技術，2007，29（增刊）：68－71.

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

［5］柳建新，籍煒，劉穎，等.基于 Matlab與Fortran混合編程的一維fCSEM正演可視化［J］.物探化探計算技術，2012，34（2）：224－228.

［6］郭良輝，孟小紅，石磊，等.MATLAB在“應用重力學”課程教學中的應用［J］.科技信息，2010（3）：12－13.

［7］邵小桃，郭勇，李一玫.“電磁場與電磁波”課程的Matlab輔助教學［J］.電氣電子教學學報，2010，32（5）：111－113.

［8］曾華霖.重力場與重力勘探［M］.北京：地質出版社，2005.

［9］管志寧.地磁場與磁力勘探［M］.北京：地質出版社，2005.

［10］柳建新，童孝忠，郭榮文，等.大地電磁測深法勘探—資料處理、反演與解釋［M］.北京：科學出版社，2012.