李靖煒，謝興兵，高文龍

(長江大學 地球物理與石油資源學院，油氣資源與勘探技術教育部重點實驗室，湖北 武漢 430100)

0 引言

直流電阻率法是電法勘探中使用較廣，較成熟的一種物探方法[1-2]。其通過人工建立的穩(wěn)定電場探測巖礦石反應的電學性質差異，來獲取地質結構信息[3-13]。隨著儀器發(fā)展和勘探需求的進步，三維正反演是適應復雜地質條件下勘探的必要技術手段，且通過三維正演算法，對地質模型進行三維可視化處理，可以很好的幫助研究人員理解地質模型電阻率的空間分布。近半個世紀以來，國內外已有很多學者使用有限元法(FE，F(xiàn)inite Element)、有限差分法(FD，F(xiàn)inite Difference)、積分方程法(IE，Integral Equation)等數(shù)值模擬方法，包括一次場法和異常電位法[7，10]，結構化網格和非結構化網格等對該問題進行了研究，并取得豐富的成果[3-19]。

本文基于有限差分法實現(xiàn)了直流電阻率的三維正演模擬及可視化，有限差分數(shù)值模擬方法具有較快的計算速度、編程實現(xiàn)簡單、占用的物理存儲空間相對較小等特點，同時也是應用較廣泛、研究較多的一類數(shù)值模擬方法。對于使用有限差分法進行直流電法三維正演問題，Dey等[3]對于該模擬問題進行了有限差分法三維數(shù)值模擬，但其求解差分格式的線性方程時效率很低。周熙攘等[5]對于該模擬問題詳細論述了有限差分法數(shù)值模擬方法及該方法的優(yōu)缺點，并使用Fortran語言做了正演計算。Zhao等[10]和吳小平等[12,20]對該三維模擬問題優(yōu)化了有限差分大型稀疏方程組的求解算法。鄧正棟等[21]對該問題進行了點源三維有限差分及異常電位法[7，10，22]的詳細推導，并做了一維驗證。王志剛，何展翔，魏文博等[23]對于該模擬問題進行了算法及求解方法等若干討論。朱占升等[17]通過MATLAB平臺對電法勘探數(shù)據(jù)進行了三維可視化處理，但其使用的面繪制很難對空間內部等進行多視角觀察。鑒于此，本文從異常電位所滿足的微分方程出發(fā)，采用有限差分算法進行直流電阻率的三維數(shù)值計算，在MATLAB下編程使得數(shù)值計算及三維可視化同時進行。在此基礎上，期望借助于MATLAB GUI編程實現(xiàn)具有三維視電阻率多角度顯示功能的直流電法三維數(shù)值模擬系統(tǒng)，幫助研究人員能夠更加方便、直觀和快捷地理解不同地球物理模型的電阻率的空間分布及其變化規(guī)律。

1 控制方程

對于均勻半無限空間中點源分布電流場的計算問題，一般需要把計算區(qū)域限定在有限的空間，即在一個有限的求解區(qū)內求解對應的方程，結合物理和數(shù)學的要求，對相應的邊界上規(guī)定邊值，確定方程的唯一解。

求解的空間模型見圖1。圖1中的Ⅰ和Ⅱ區(qū)域分別代表不同電導率分布的異常體空間和背景空間，使用異常電位法，將電位分為背景電位U0和異常電位Ua，即U=U0+Ua，σa表示實際電導率σ減去背景電導率σ0，即σa=σ-σ0，T0，T1…T5依次分別代表空間模型的上邊界、左邊界、前邊界、右邊界、下邊界和后邊界。

圖1 三維直流電法求解空間模型

在三維笛卡爾坐標系中，若地下均勻半空間中場源是一個坐標位于地面(x0，y0，z0)點，電流強度為I的點電流源，通過推導可得到異常電位Ua滿足的微分方程(1)：

?·(σ?Ua)=-?·(σa?U0)

(1)

邊界條件滿足方程(2)：

(2)

其中：σ=σ(x，y，z)，U=U(x，y，z)。在正演計算中U0，σa和σ和邊界條件是給定的，只有Ua是未知數(shù)，所以異常電位Ua的解是唯一的。

2 有限差分法數(shù)值模擬

點電流源的電位在三維空間滿足的方程及定解條件在數(shù)學上稱為求解偏微分方程的邊值問題。本文采用有限差分法解決此邊值問題。有限差分法是用差商代替微分，把理論上的電場微分方程變換為適于數(shù)值計算的差分形式，形成能夠數(shù)值計算的差分方程。

在電法勘探三維正演模擬中根據(jù)有限差分中剖分網格形式的不同，可以分為交錯網格[24]、常規(guī)網格、無網格差分[25]等。綜合分析算法的計算成本及對硬件的要求，本文選擇使用常規(guī)六面體網格進行有限差分。

2.1 剖分網格

為便于計算，需要對求解區(qū)域進行網格剖分，將電位U賦值于六面體網格節(jié)點上。數(shù)值計算中使用符號i、j和k對每個節(jié)點編號，區(qū)分表現(xiàn)為網格節(jié)點沿X、Y和Z軸方向的編號，令hi=Δxi、hj=Δyj和hk=Δzk示網格節(jié)點U(xi，yj，zk)沿X、Y和Z軸正方向的步長，網格節(jié)點電位Ui，j，z=U(xi，yj，zk)；單元電導率σi，j，z=σ(xi，yj，zk)。網格剖分樣式見圖2，對于小體積異常體，對異常體部分進行細小的網格剖分，對異常體之外的部分使用常規(guī)剖分。

圖2 網格示意圖

2.2 差分格式

通過對異常電位微分方程(1)的差分變換，得到線性方程(3)：

[A]·{Ua}=-[Aa]·{U0}

(3)

在數(shù)值計算中，需要注意邊界內的節(jié)點電位與邊界上的節(jié)點電位的計算。當節(jié)點(i，j，k)位于邊界內時，線性方程(3)中的矩陣 [A]、[Aa]和{Ua}為

[A]=[ai，j，kai-1，j，kai，j-1，kai，j，k-1ai+1，j，kai,j+1，kai，j，k+1]

(4)

式中：

hi表示沿i方向的步長，

ai，j，k=-[ai-1，j，k+ai，j-1，k+ai，j，k-1+ai+1，j，k+ai,j+1，k+ai，j，k+1]，

其余(i，j，k)相鄰節(jié)點系數(shù)與此類似。

{Ua}=[Uai，j，kUai-1，j，kUai，j-1，kUai，j，k-1Uai+1，j，kUai,j+1，kUai，j，k+1]T

(5)

[Aa]=[aai，j，kaai-1，j，kaai，j-1，kaai，j，k-1aai+1，j，kaai,j+1，kaai，j，k+1]

(6)

式中：

aai，j，k=-(aai-1，j，k+aai，j-1，k+aai，j，k-1+aai+1，j，k+aai,j+1，k+aai，j，k+1)，

其余(i，j，k)點相鄰系數(shù)形式與其相似。

當節(jié)點(i，j，k)位于地面T0邊界上時，因為地面節(jié)點是半個單元，空間中不存在(i，j，k-1)節(jié)點，同時場源位于地面，其滿足鏡像法原理[5]，所以在地面邊界上，其相應的節(jié)點點電流源的電流強度放大一倍，線性方程(3)中的矩陣 [A]、[Aa] 和 {Ua} 為

[A]=[ai，j，kai-1，j，kai，j-1，kai+1，j，kai,j+1，kai，j，k+1]

(7)

式中：

ai，j，k=-(ai-1，j，k+ai，j-1，k+ai+1，j，k+ai,j+1，k+ai，j，k+1)，

其余(i，j，k)相鄰系數(shù)形式與其相似。

同理，{Ua} 中沒有Uai，j，k-1節(jié)點，其中aai，j，k+1為

{Ua}=[Uai，j，kUai-1，j，kUai，j-1，kUai+1，j，kUai,j+1，kUai，j，k+1]T

(8)

[Aa]=[aai，j，kaai-1，j，kaai，j-1，kaai+1，j，kaai，j+1，kaai，j，k+1]

(9)

式中：

aai，j，k=-(aai-1，j，k+aai，j-1，k+aai+1，j，k+aai,j+1，k+aai，j，k+1)，

其余(i，j，k)相鄰系數(shù)形式與其相似。

當節(jié)點(i，j，k)位于T1，T2，…，T5邊界上時，此類邊界面上不存在場源，進行數(shù)值計算時所設異常模型距離邊界面相對較遠，所以令異常電位Uai，j，k=0，令邊界面上的電位等于U0。

在電導率不同的分界面上，因為差分過程中是對節(jié)點及其相鄰節(jié)點進行差商代替，所以在差分方程中電位在此分界上是自然連續(xù)的。

結合參數(shù)矩陣(4)～(9)求解線性方程(3)即可確定異常場電位，最后將有限差分法求得的異常場電位Ua與解析解求得的正常場電位U0相加，可以近似得到實際電場電位U。

2.3 線性方程組的求解

本文選取逐次超松弛(Successive Over Relaxati-on， SOR)迭代法[4]求解線性方程(3)。點電流源三維數(shù)值模擬離散網格的節(jié)點總數(shù)一般可達百萬級別，形成的矩陣為巨型稀疏陣，線性方程組的階數(shù)很高，若直接求取，雖然速度較快，但占用內存量較大，累計誤差也較嚴重。若選用常規(guī)迭代法求解，占用內存少，編制程序較簡單，也可以克服誤差的累積，但收斂速度較慢。在保持迭代的優(yōu)點基礎上，為加快收斂速度，選取SOR迭代法。超松弛迭代法首先對每一個節(jié)點上的電位賦予初值(本文異常電位Ua初值賦予0)，而后依照節(jié)點編號依次由差分方程求解各節(jié)點的電位，從而每一個節(jié)點通過超松弛迭代計算該節(jié)點的電位值：

(10)

3 算例

為了驗證算法的正確性，通過使用二級裝置，得到三維視電阻率數(shù)據(jù)，然后，三維數(shù)值模擬結果與解析解一維視電阻率曲線對比并檢驗算法是否正確。

模型Ⅰ：空間范圍100 m×100 m×100 m，點源位于地面T0中心(空間模型見圖1)，背景電阻率300 Ω·m，低阻球形異常體電阻率30 Ω·m，半徑10 m，球心埋深50 m(點源正下方)。

通過數(shù)值計算求得三維視電阻率數(shù)據(jù)，過點源取一垂直切片，由圖3可見，紅線框標識了異常體大致位置(點源位于異常體正上方)，從圖中可以明顯看出低阻異常體位置以及視電阻率的連續(xù)性變化，這可以說明算法是適合的；通過數(shù)值解與解析解對比驗證，得到一維視電阻率對比曲線，由圖4可見，從視電阻率ρs曲線上來看，數(shù)值解具有較好的對稱性(圖4a)，從相對誤差來看，其相對誤差在0.5%以內(圖4b)，說明三維數(shù)值計算的精度相對較高，用于研究異常體的變化規(guī)律和可視化是可行的。

圖3 差分數(shù)值解三維視電阻率切片

圖4 差分數(shù)值解與解析解對比

模型Ⅱ：空間范圍100 m×100 m×100 m，點源位于地面T0中心，背景電阻率300 Ω·m，使用兩個立方體異常模型，低阻30 Ω·m，高阻300 Ω·m，異常體邊長均為10 m，體中心埋深均為50 m。

為進一步驗證算法，使用模型Ⅱ的組合異常體模型。異常體位于點源正下方的兩側，通過二級裝置計算得到三維視電阻率，平行X軸取過點源與兩個異常體的切片見圖5，紅線框標識了異常體大致位置。從圖5中可以明顯分辨出低阻體和高阻體，及其空間位置和異常體視電阻率的形態(tài)，這也說明當空間存在組合異常體時，該算法依然是可行的。

圖5 組合異常體三維視電阻率切片

4 可視化及MATLAB GUI界面

可視化即對數(shù)值計算的三維數(shù)據(jù)體進行三維可視化，圖6是算例1中低阻球形異常體視電阻率3D等值線可視化圖，其主要使用MATLAB的函數(shù)“contourslice”對視電阻率三維數(shù)據(jù)進行可視化處理。結合MATLAB GUI開發(fā)應用，即使用GUI圖形對象組成用戶界面，通過這種用戶界面，三維數(shù)值模擬的命令和對程序的控制是通過“選擇”各種圖形對象來實現(xiàn)的。

圖6 低阻異常體視電阻率三維等值線可視化

4.1 GUI開發(fā)環(huán)境

本文使用MATLAB平臺，其GUIDE(GUI Development Environment)提供了一套可視化的創(chuàng)建圖形窗口(GUI)的工具，包括布局編輯器(Layout Edtor)、幾何排列工具(Alignment Tool)、屬性編輯器(Property Inspector)、對象瀏覽器(Object Browser)及菜單編輯器(Menu Editor)，使用GUIDE可方便創(chuàng)建GUI應用程序。在MATLAB中，GUI的設計是以M文件的編程模式完成的，它能夠按照用戶設計的GUI布局，將GUI的布局代碼存儲在FIG文件中，同時產生一個M文件用于存儲調用函數(shù)，在M文件中不再蘊含GUI的布局代碼。利用這一框架編制自己的應用程序，在開發(fā)應用程序時代碼量大大縮小。

4.2 開發(fā)流程

GUI圖形界面的功能，還是要經過一定的設計思路和計算方法，由特定的程序來完成。開發(fā)GUI應用流程圖見圖7(Z軸表示深度，二極裝置等效于AO的距離，為更貼合地下三維空間，筆者使用三維笛卡爾坐標表示)，為了實現(xiàn)程序的功能，需要通過GUIDE完成操作界面的布局及控件的標識，還需要在運行程序前編寫有限差分三維模擬算法代碼，完成程序中變量的賦值、輸入輸出、計算及繪圖等工作，最終實現(xiàn)數(shù)據(jù)可視化。

圖7 GUI應用開發(fā)流程圖

4.3 應用

運行程序(平臺MATLAB)，通過命令窗口直接輸入文件名或用openfig，open或hgload命令即可運轉GUI程序。

在模擬實驗中，通過點擊不同的菜單欄選項，可以得到對應的視電阻率三維可視化圖。不同的菜單對應不同的模擬實驗，部分三維直流視電阻率可視化程序操作界面見圖8。同時，使用者可在對應的文本框中設置電阻率參數(shù)等，這樣使得操作者不必對算法有深入的認識，可以通過參數(shù)設置等對其快速的模擬，這也是GUI界面的優(yōu)點。相反，該應用程序也有一些小缺點是可移植性差，對于MATLAB開發(fā)的GUI界面程序，其依賴于MATLAB Mathworks工作平臺，即使包裝成“.exe”Windows程序，其仍需要安裝MATLAB相關輔助程序。

圖8 三維直流視電阻率可視化程序界面

5 結論

本文使用經典的有限差分法、異常電位法和SOR迭代求解三維數(shù)值模擬問題，論述了數(shù)值計算過程及GUI開發(fā)過程，目的不是追求在算法上進行創(chuàng)新和突破，而是在于通過具有豐富的可視化包的MATLAB平臺進行視電阻率三維數(shù)值模擬及可視化，來避免使用第三方軟件建?；虺蓤D帶來的不便。同時，由三維正演模型，可以得到很多有意義的三維直流電法響應的相關結論。

1)通過算例分析表明，基于MATLAB使用有限差分對視電阻率進行三維數(shù)值模擬及可視化，不僅是可行的，而且算法簡單有效，同時三維數(shù)值計算的結果也具有很高的精度。

2)通過對離散異常電位微分方程得到的系數(shù)矩陣研究，使用SOR迭代法，其差分格式通過MATLAB編程不僅不用考慮對大型稀疏系數(shù)矩陣進行壓縮儲存，而且計算速度較快。

此外，筆者編寫的三維可視化GUI應用程序只是一個初級階段，但已具備方便快速的計算、多角度直觀的三維可視化、簡約的操作界面等諸多優(yōu)點，有信心通過加入更多的三維數(shù)值計算模塊以及對MATLAB可視化編程和GUI編程的深層次開拓，其能夠形成功能更為全面的三維數(shù)值模擬可視化系統(tǒng)。