彭 淼,熊 杰

(中國地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院,北京 100083)

基于Matlab的直流電阻率三維數(shù)值模擬及可視化

彭 淼,熊 杰

(中國地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院,北京 100083)

針對直流電阻率法的三維正演問題,采用積分方程法的數(shù)值解法,特別地,對板狀體三維模型做了具體的網(wǎng)格剖分,在Matlab下編程實(shí)現(xiàn)了求解其電位分布的算法.根據(jù)電位分布分別采用中間梯度法,三極剖面法和聯(lián)合剖面法計(jì)算了地面的視電阻率分布.結(jié)果表明,視電阻率異常圖充分顯示了積分方程求解的正確性.對于該結(jié)果,在Matlab中采用G UI編程生成了更為方便的人機(jī)交互式正演界面,該界面具有直觀、高效的特點(diǎn).關(guān)鍵詞：Matlab;積分方程;直流電阻率法;三維正演;G UI

0 引 言

直流電阻率法是地球物理勘探中一種常用的方法.它是以巖礦石導(dǎo)電性差異為基礎(chǔ),通過觀測、研究人工建立的地下穩(wěn)定電流場的分布規(guī)律以達(dá)到找礦和解決其他地質(zhì)問題為目的的一種電法勘探方法[1].研究直流電阻率三維正演問題是電法勘探中的一項(xiàng)重要任務(wù),近30年來國內(nèi)外已有不少學(xué)者對該問題采用有限單元法、有限差分法、積分方程法等數(shù)值模擬方法進(jìn)行了算法研究,成果顯著[2-4].

本文采用積分方程法進(jìn)行直流電阻率的三維正演.對于該正演問題,周熙襄等[1]詳細(xì)闡述了積分方程的近視求解方法,并對具體模型作了正演計(jì)算.他們實(shí)現(xiàn)該算法是用FORTRAN語言編寫的,而本文采用可視化功能更為強(qiáng)大的Matlab來實(shí)現(xiàn).筆者不僅最終成功實(shí)現(xiàn)了直流電阻率的積分方程法三維正演算法的高效計(jì)算,而且還采用G UI編程開發(fā)了正演用戶界面,這使得方便、快捷的人機(jī)交互式正演成為可能.

1 積分方程及其求解方法

假設(shè)在無限大半空間的均勻介質(zhì)中,存在一個三維地質(zhì)體,它的電性與周圍介質(zhì)的電性不同.位于地面的一個電流源供給穩(wěn)定電流 I,均勻介質(zhì)的電阻率為ρ1,三維地質(zhì)體的電阻率為ρ2,表面電荷密度為 q.根據(jù)基本的電場和電磁場定律及邊界條件,可以導(dǎo)出如下積分方程：

式(1)、(2)為電阻率突變界面積累電荷密度 q所滿足的積分方程式.當(dāng)外電源和電阻率分布給定時,求解式(1)、(2)即可得到 q的分布.通過已知的q分布求解式(3)便可以計(jì)算出空間任意點(diǎn)的電位.

對于上面的積分方程,要給出其解析解是很困難的,本文采用數(shù)值計(jì)算的近似解法,其具體步驟如下：

首先,將積分域離散化為有限個小面積的組合.然后,將每一個小面積用一個簡單的平面近似代替.設(shè)小平面內(nèi)的q為常數(shù),這樣,積分方程就變成了一組線性方程組：

式中,qi(或qj)為第i個(或第j個)小平面內(nèi)的積累電荷密度,它在該小平面內(nèi)為一常數(shù),Δsj(或Δs—j)為小面元的面積.

式(4)中,當(dāng)i=1,2,…,n時,可以得到n個線性方程,組成一個 n階線性方程組,其矩陣方程為,

其中,系數(shù)矩陣A的對角線元素為,

非對角線元素為,

右端項(xiàng)為,

然后,對系數(shù)矩陣各元素中的面積分項(xiàng)進(jìn)行計(jì)算,有以下近似表達(dá)式：

將式(3)～(9)聯(lián)立起來求解,就可以得到q的分布,此即為積分方程(1)的數(shù)值近似解.

2 板狀體三維模型

上述積分方程的數(shù)值近似求解方法適用于任意形體的三維模型計(jì)算.本文采用的是板狀體三維模型(見圖1).圖1中上部的長方體代表虛源.

圖1 板狀體三維模型示意圖及表面網(wǎng)格剖分

對于一個長、寬、高分別為A、B、C的板狀體,本文采取以下方式劃分表面的小面元：把長方體的每條棱劃分成N段,其中 a=A/N,b=B/N,c= C/N,那么共出現(xiàn)3種不同的小面元,a×b,b×c和c×a,一共有6×N2個小面元.這樣劃分的優(yōu)點(diǎn)是避免了奇異面元的討論,從而能很方便地實(shí)現(xiàn)算法.

由于每一個劃分得足夠小的面元內(nèi)的積累電荷密度 qi都可以視為常數(shù),故對于式(4)、(5)、(6)中的積分項(xiàng),可以采用下列近似公式求解：

通過對板狀體上頂面、底面、前側(cè)面、后側(cè)面、左側(cè)面以及右側(cè)面的分類討論,計(jì)算系數(shù)矩陣和右端項(xiàng),并求解線性方程組可得到q的分布,然后利用下式計(jì)算地表任意點(diǎn)的電位：

最后,本文在Matlab的編程環(huán)境中實(shí)現(xiàn)了計(jì)算地表電位分布的算法,并將程序模塊保存于相應(yīng)文件中.

3 視電阻率的計(jì)算及算例

為了檢驗(yàn)程序的有效性和正確性,本文分別采用中間梯度法、三極剖面法以及聯(lián)合剖面法計(jì)算得到了地表的視電阻率分布.計(jì)算時均采用以下公式：

對于三極剖面法和聯(lián)合剖面法,其裝置系數(shù)為,

而對于中間梯度法,除計(jì)算主測線上的視電阻率值外,還可計(jì)算在離開主測線上一定距離且與之平行的旁測線上的視電阻率值,其裝置系數(shù)為：

中間梯度法板狀體三維模型正演主測線視電阻率對比如圖2所示.

實(shí)線參數(shù)如下：頂部埋深為5 m的板狀體：長B=20 m,寬A =10 m,高C=10 m;直流電源位于地面主測線,相距AB= 420 m;電流I=20 A;圍巖電阻率ρ1=1 000Ωm;板狀體電阻率ρ2=20Ωm;測量電極MN=2 m.其他曲線參數(shù)除長和圍巖電阻率按圖例改變外,其余不變.圖2 中間梯度法板狀體三維模型正演主測線視電阻率對比

從圖2的對比可以看出,其視電阻率曲線完全符合理論解釋,即：都在模型中點(diǎn)出現(xiàn)低阻異常極小值;都在遠(yuǎn)離模型的兩側(cè)接近圍巖真實(shí)電阻率;靠近模型兩側(cè)都出現(xiàn)小的極大值.這說明正演結(jié)果是穩(wěn)定可靠的.

由于三維正演可以得到地面任意位置的電位,根據(jù)該點(diǎn)電位可以計(jì)算地表視電阻率的分布.圖3、圖4和圖5是分別采用中間梯度法,三極剖面法以及聯(lián)合剖面法計(jì)算得到的地面視電阻率的分布圖.這3張圖都很直觀地顯示了板狀體的大致位置.其中：圖3中心的低阻異常與板狀體中心位置吻合;圖4和圖5的對比也進(jìn)一步說明了正演結(jié)果的正確性.

參數(shù)設(shè)置：頂部埋深為5 m的板狀體,長B=20 m,寬A= 10 m,高C=10 m;電流I=20 A;圍巖電阻率ρ1=1 000 Ωm;板狀體電阻率ρ2=20Ωm;測量電極MN=2 m.圖3 中間梯度法板狀體三維模型正演地表視電阻率平面圖

參數(shù)設(shè)置：頂部埋深為10 m的板狀體,長B=10 m,寬A= 10 m,高C=6 m;AM=8 m;MN=2 m;電流I=20 A;圍巖電阻率ρ1=500Ωm;板狀體電阻率ρ2=50Ωm.圖4 三極剖面法板狀體三維模型正演地表視電阻率平面圖

參數(shù)設(shè)置：頂部埋深為10 m的板狀體,長B=10 m,寬A= 10 m,高C=6 m;AM=8 m;MN=2 m;電流I=20 A;圍巖電阻率ρ1=500Ωm;板狀體電阻率ρ2=50Ωm.圖5 聯(lián)合剖面法板狀體三維模型正演地表視電阻率平面圖

需要說明的是,以上圖示結(jié)果的正演計(jì)算均采用6×10×10的網(wǎng)格剖分,與更為精細(xì)的6×15×15的網(wǎng)格剖分僅存在0.38%的平均誤差,而以上圖示的數(shù)值計(jì)算和成圖的時間只需4～8 s(計(jì)算機(jī)的配置為：英特爾酷睿2雙核E4500CPU,2 G內(nèi)存).因此,該程序在Matlab平臺上的運(yùn)行是高效的.

4 GUI界面設(shè)計(jì)

對于上述正演模擬,我們通過修改模型參數(shù)可以獲得若干種正演結(jié)果.但是,每次修改參數(shù)都需要在程序代碼中進(jìn)行,這樣既費(fèi)時又費(fèi)事,給正演研究過程帶來諸多不便.為解決這一問題,我們在Matlab平臺上采用G UI編程設(shè)計(jì)了人機(jī)交互式正演界面.

人機(jī)交互式正演界面的右側(cè)是參數(shù)輸入面板和操作控制臺,其中集成了中間梯度法、三極剖面法以及聯(lián)合剖面法3種方法求解地表視電阻率.用戶只需依次輸入想要改變的參數(shù)值,然后左鍵單擊正演模擬按鈕便能完成從數(shù)值模擬到成圖的所有過程.人機(jī)交互界面的左側(cè)是成圖區(qū),包含了模型的三維立體圖以及正演結(jié)果圖.其中模型的三維立體圖可以分二維和三維變換顯示.采用G UI編程用戶從正演開始到成圖結(jié)束的等待時間只需4～8 s,且正演過程能夠重復(fù)進(jìn)行.圖6顯示了采用中間梯度法正演模擬完成后的用戶界面.

圖6 中間梯度法板狀體三維模型人機(jī)交互式正演界面圖

5 結(jié) 語

本文利用經(jīng)典的解三維問題的積分方程法進(jìn)行數(shù)值模擬,目的不在于對該算法進(jìn)行突破創(chuàng)新,而在于提出一種能在Matlab平臺上解決電法正演問題的思路.同時,我們開發(fā)的人機(jī)交互式正演界面只是一個初步嘗試,卻已具備簡潔、直觀、快速等諸多優(yōu)點(diǎn),相信通過引入更多的正演程序模塊以及G UI編程的進(jìn)一步拓展,其能夠發(fā)展成為功能更為強(qiáng)大的正演模擬系統(tǒng).

[1]周熙襄,鐘本善.電法勘探數(shù)值模擬技術(shù)[M].成都：四川科學(xué)技術(shù)出版社,1986.

[2]李金銘.地電場與電法勘探[M].北京：地質(zhì)出版社,2004. [3]傅良魁.應(yīng)用地球物理教程[M].北京：地質(zhì)出版社,1991. [4]江玉樂,雷 宛.地球物理數(shù)據(jù)處理教程[M].北京：地質(zhì)出版社,2006.

[5]顧觀文.電阻率/激電測深二維人機(jī)交互正演模擬[J].物探化探計(jì)算技術(shù),2007,29(S1)：89-105.

[6]李曉昌.在Matlab平臺上實(shí)現(xiàn)可控源音頻大地電磁反演數(shù)據(jù)三維可視化顯示[J].物探化探計(jì)算技術(shù),2007,29 (S1)：68-104

3-D Numerical Simulation and Visualization of DC Resistivity Based on Matlab

PENG Miao,XIONG Jie

(School of Geophysics and Information Technology,China University of Geosciences,Beijing 100083,China)

Integral equation method was used in this paper to solve the 3-Dforward problemof DC resistivity.Especially,specific gridding dissection was made on the 3-D tabular model and an algorithm for solving its potential distribution was realized by programming on Matlab.Central gradient array,tri-electrode profile and combined profile survey were utilized to calculate its apparent resistivity distribution on ground according to the potential distribution.The results indicate that this solution of the integral equation is correct according to apparent resistivity anomalyfigure.Besides,a human interface for more convenient forward stimulation was created by G UI programming on Matlab.This interface is intuitive and efficient.

Matlab;integral equation;DC resistivity;three-dimensional forward simulation;G UI

P631.3+25

：A

2010-07-25.

彭 淼(1963—),男,博士研究生,從事電法勘探技術(shù)研究.

1004-5422(2010)03-0213-04