李清坤, 曹 輝, 曹禮剛, 邢雯淋, 張朝暉

(成都理工大學(xué) 地球勘探與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，成都 610059)

0 引言

作為一種高靈活性且實(shí)用性的地球物理勘探方法，直流電阻率法常用于各大工程領(lǐng)域中,因此直流電阻率法的正反演研究工作一直備受關(guān)注[1-3]。自20世紀(jì)70年代以來(lái)，有限元法被應(yīng)用到地球物理正演計(jì)算。Coggon[4]提出的方法為有限單元法在地球物理中的應(yīng)用奠定了正演理論基礎(chǔ)；BIBBY[5]在Coggon工作的基礎(chǔ)上給出了軸對(duì)稱(chēng)體直流電阻率模型計(jì)算方法，這大大地增加了直流電法在復(fù)雜地形的應(yīng)用的可能性；Dey等[6]將有限差分方法引入點(diǎn)源二維任意地電剖面電阻率法正演計(jì)算中；周熙襄等[7]在Dey等提出的方法上加入了混合邊界條件進(jìn)行正演計(jì)算并給出試算結(jié)果，提高了正演計(jì)算的效率和精度，進(jìn)一步證明了有限差分法的實(shí)用價(jià)值；徐世浙[8]在書(shū)中詳細(xì)介紹了對(duì)于直流電阻率二維的正演有限元法數(shù)值計(jì)算過(guò)程；底青云[9]運(yùn)用直流電場(chǎng)有限元方法進(jìn)行了二維線(xiàn)源正演模擬，探討了改進(jìn)復(fù)雜結(jié)構(gòu)下的有限元正演模擬的實(shí)效性；阮百堯等[10]提出了具有高精度的有限元方法數(shù)值解，解決了二維地電斷面電阻率測(cè)深的變分問(wèn)題；Greenhalgh等[11]介紹了一種基于高斯正交網(wǎng)格(GQG)的新型電阻率建模方法的理論，可以在具有任意地形條件下對(duì)2.5D/3D各向異性模型中計(jì)算出電勢(shì)；湯井田等[12]運(yùn)用非結(jié)構(gòu)化三角形網(wǎng)格,對(duì)起伏地形下2.5D直流電阻率進(jìn)行正演模擬,并利用比較法來(lái)消除地形的影響，獲得了良好的效果。

在反演研究中，Loke和Barker[13-14]采用擬牛頓方法對(duì)高斯-牛頓方法進(jìn)行了改進(jìn)，使得反演的精度得到了保證；阮百堯等[15]提出了一種二維直流電阻率測(cè)深的快速反演方法，在實(shí)際的生產(chǎn)工作取得了很好的效果；吳小平等[16]采用共軛梯度迭代法，實(shí)現(xiàn)了直流電阻率測(cè)量數(shù)據(jù)的三維最小構(gòu)造反演；Singh等[17]提出了基于神經(jīng)網(wǎng)絡(luò)的二維電阻率測(cè)深快速優(yōu)化反演算法，并應(yīng)用實(shí)際的電阻率測(cè)深數(shù)據(jù)反演中，使得數(shù)據(jù)解釋更為快速、高效、準(zhǔn)確；馮德山等[18]采用有限單元法對(duì)超高密度電法的全四極裝置進(jìn)行正演模擬，采用廣義最小二乘正則化反演,較好地刻畫(huà)出復(fù)雜地形條件下的異常體的形態(tài)及地下電性的分布；吳小平等[19]基于非結(jié)構(gòu)化四面體網(wǎng)格，采用不完全Gauss-Newton反演法，實(shí)現(xiàn)了對(duì)任意地形條件下偶極-偶極視電阻率數(shù)據(jù)的三維反演；周勇等[20]使用了基于模型敏感度信息的非結(jié)構(gòu)化自適應(yīng)網(wǎng)格算法來(lái)生成滿(mǎn)足反演要求的高質(zhì)量網(wǎng)格，并使用了穩(wěn)定的雙共軛梯度法來(lái)求解Gauss-Newton方程以實(shí)現(xiàn)三維直流電阻率穩(wěn)定反演法。

筆者在前人研究的基礎(chǔ)上，研究了起伏地形條件下直流電阻率2.5D反演算法，為直流電阻率法的數(shù)據(jù)處理解釋提供理論支持。首先，將有限體積法引入二維正演線(xiàn)性方程組中求解，提高了正演數(shù)值模擬的精度和速度；然后，將反演中常用的正則化Tikhonov反演思想引入其中，并在反演算法中采用多線(xiàn)程并行計(jì)算求解矩陣，大大節(jié)省了內(nèi)存和時(shí)間的消耗，最后，對(duì)pole-dipole和dipole-pole和dipole-dipole裝置的反演效果進(jìn)行了對(duì)比分析。

1 正演算法

假設(shè)地下全空間電導(dǎo)率為σ，點(diǎn)電流源位于地表，則滿(mǎn)足直流電阻率法的2.5維偏微分方程(PDEs)及其邊界條件(BC)公式為式(1)[21]。

(1)

式中：φ為電位向量,V；σ為電導(dǎo)率,S/m；I為供電電流強(qiáng)度,A·m-3；rs+為供電電極正極位置；rs-為供電電極負(fù)極位置。在直流電阻率法中，使用接收電極對(duì)電位場(chǎng)中的差異進(jìn)行采樣，以收集觀(guān)測(cè)數(shù)據(jù)。為了模擬該偏微分方程(PDEs)(或一組PDEs，如果存在多個(gè)電流源位置)，必須將該方程離散化到計(jì)算網(wǎng)格上。

此方程在地表空氣邊界滿(mǎn)足Newman邊界條件，在無(wú)窮遠(yuǎn)處滿(mǎn)足混合邊界條件，通過(guò)使用交錯(cuò)網(wǎng)格有限體積法(FVM)，我們將系統(tǒng)離散化為：

(2)

Aφ=-q

(3)

其中

(4)

式中：A為系數(shù)矩陣；q為場(chǎng)源向量，采用不完全的LU分解法解方程[22-23]，即可得到電位U和視電阻率值ρs

(5)

其中：ε為接收電極處電位向量；ΔU為電位差；K為裝置系數(shù)。

2 反演算法

地球物理中反演一直是一項(xiàng)復(fù)雜而艱巨的任務(wù)，由于反演問(wèn)題本身就存在多解性和不確定性，基于此筆者采用Tikhonov正則化反演方法[24-26]，建立目標(biāo)函數(shù)為式(6)。

φ(m)=φd(m)+βφm(m)

(6)

式中：φ為正則化因子;φ(m)為數(shù)據(jù)目標(biāo)函數(shù)項(xiàng)；φm(m)為模型目標(biāo)函數(shù)項(xiàng)。

這里采用基于模型與先驗(yàn)信息mref的最小二范數(shù)

(7)

式中：Wm為模型權(quán)重，是x、y和z方向上一階平滑度的矩陣。每個(gè)Wm矩陣可以用整數(shù)離散表示。正則化Wm是應(yīng)用為每個(gè)單元格具有不同權(quán)重的標(biāo)量或?qū)蔷仃嚰訖?quán)總和α*的矩陣。

將數(shù)據(jù)目標(biāo)函數(shù)寫(xiě)為式(8)[24]。

(8)

式中：Wd為數(shù)據(jù)的權(quán)重矩陣，該矩陣為對(duì)角矩陣，其元素為Wdii=1∈i，其中εi是第i個(gè)標(biāo)準(zhǔn)偏差；F(m)為正演函數(shù)；m為模型參數(shù)；d為觀(guān)測(cè)所得的數(shù)據(jù)。

則其梯度為：

(9)

式中：J[m]為靈敏度矩陣;J[m]ij為第i個(gè)基準(zhǔn)點(diǎn)相對(duì)于第j個(gè)模型參數(shù)的變化方式；在第k次迭代中，從模型mk開(kāi)始搜索擾動(dòng)δm以減少目標(biāo)函數(shù)，可以給出反演模型所需要的公式：

F[mk+δm]≈F[mk]+J[mk]

(10)

求解該擾動(dòng)，將梯度設(shè)為0

(11)

更新的模型由mk+1=mk+γδm給出，其中γ∈(0,1]，該系數(shù)可以通過(guò)一維搜索找到，其默認(rèn)值設(shè)為“1”。直到達(dá)到合適的停止標(biāo)準(zhǔn)方可停止迭代優(yōu)化過(guò)程。反演流程圖如圖1所示。

圖1 反演并行算法流程圖

3 靈敏度

對(duì)于靈敏度的求解，我們參考文獻(xiàn)[24]，直流電阻率法的正演問(wèn)題為對(duì)形如式(12)的偏微分方程進(jìn)行求解：

c(m,u)=A(m)u-q=0

(12)

其中:m為正演模型；u表示需要求解的場(chǎng)。對(duì)于給定的模型m，可以計(jì)算出場(chǎng)u(m)，由于求解得到的數(shù)據(jù)是場(chǎng)的子集，因此可以由線(xiàn)性投影P來(lái)從模型中創(chuàng)建投影數(shù)據(jù)

dpred=Pu(m)

(13)

其中P是字段在數(shù)據(jù)空間上的投影。

對(duì)于反演問(wèn)題，主要是解決如何通過(guò)改變模型來(lái)轉(zhuǎn)換數(shù)據(jù)。將式(13)泰勒展開(kāi)得式(14)。

O(h2‖v‖)

(14)

根據(jù)式(14)，改變模型得到相應(yīng)反演數(shù)據(jù)，也是我們所求的解的方程?；诖穗x散化，我們得出離散化空間中的靈敏度。靈敏度矩陣可以寫(xiě)成式(15)。

J=-PA-1G

(15)

從而可以對(duì)偏微分方程(PDE)進(jìn)行求導(dǎo)

(16)

將式(16)代入式(15)中可以得到靈敏度矩陣為式(17)。

J=-P(uc(m,u))-1mc(m,u)

(17)

4 正確性驗(yàn)證

為了驗(yàn)證該直流電阻率法正演算法的準(zhǔn)確性與可行性，設(shè)置了一個(gè)場(chǎng)源在地面的均勻全空間模型，其電阻率為100 Ω·m。圖2為均勻全空間模型偶極裝置數(shù)值模擬結(jié)果，從圖3可知，在靠近坐標(biāo)軸原點(diǎn)區(qū)域，解析解與數(shù)值計(jì)算的結(jié)果出入較大，其余計(jì)算區(qū)域的正演模擬結(jié)果與解析解吻合較好，總體平均誤差小于0.1%。證明了該正演算法是準(zhǔn)確可靠的，可以為反演的模型做進(jìn)一步地準(zhǔn)備。

圖2 偶極裝置解析解與數(shù)值模擬計(jì)算結(jié)果對(duì)比圖

圖3 解析解與數(shù)值模擬計(jì)算結(jié)果誤差圖

5 算例分析

5.1 算例一

構(gòu)建的模型如圖4所示,長(zhǎng)為200 m，深為60 m，包含高阻和低阻不同形狀異常體的二層帶地形模型，淺層為1 000 Ω·m的圓形高阻體，異常體的中心坐標(biāo)(100 m，-20 m)，其半徑為10 m；另設(shè)兩個(gè)10 Ω·m的長(zhǎng)方形形高阻體，異常體的長(zhǎng)為40 m，寬10 m，圍巖電阻率為100 Ω·m。網(wǎng)格剖分密度為1 m×1 m，測(cè)線(xiàn)沿著X軸方向布設(shè)20個(gè)供電電極，測(cè)點(diǎn)布設(shè)間隔為10 m。AB極距為10 m，MN極距為10 m。

圖4 平地形模型

如圖3建立的數(shù)值計(jì)算模型，進(jìn)行直流電阻率法的四極和三極裝置進(jìn)行正演數(shù)值計(jì)算。首先對(duì)模型分別采用偶極-偶極(dipole-dipole)、單極-偶極(pole-dipole)、偶極-單極(dipole-pole)、三種布極方式，對(duì)構(gòu)建的模型進(jìn)行直流電阻率法正演計(jì)算，得到視電阻率擬剖面(圖5)?？梢钥吹剑琩ipole-dipole方式能較好分辨出高阻體和低阻體的存在，pole-dipole和dipole-pole方式雖能反映出分別存在高低阻體，但對(duì)于埋深和方位的理解存在較大偏差。

圖5 視電阻率擬剖面

圖6為視電阻率直方圖。由圖6可以看到，三種布極方式描述的視電阻率直方圖能夠較好地反映出視電阻率的集中分布情況，對(duì)于高阻體的存在給出了明確的信息，但對(duì)于低阻體的反映則不明顯。

圖6 視電阻率直方圖

采用靈敏度加權(quán)反演后，得到靈敏度如圖7所示。由圖7可以看到，三種布極方式對(duì)于淺層有著較高的靈敏度，隨著地層變深而靈敏度降低，當(dāng)存在高阻體時(shí)，靈敏度存在較明顯地降低，而存在低阻體時(shí)，靈敏度降低相對(duì)周?chē)鷩鷰r存在明顯地變緩。

圖8為恢復(fù)電導(dǎo)率反演結(jié)果，圖8中黑色實(shí)心倒三角為布設(shè)的20個(gè)電極。本次通過(guò)建立模型進(jìn)行數(shù)值模擬，正演過(guò)后的視電阻率帶入反演算法中，通過(guò)對(duì)比反演成像之后的結(jié)果圖與原始模型可以發(fā)現(xiàn)，建立的模型與反演結(jié)果吻合得很好，證明了本文算法的可靠性。可以看到，對(duì)于高阻體，三種布極方式均能正確反演出來(lái)，其中dipole-dipole方式反演結(jié)果較為理想。而對(duì)于pole-dipole和dipole-pole低阻體則在圖8反演圖中可以看出，對(duì)低阻體的形態(tài)有所偏差。在平地形模型下，dipole-dipole裝置的視電阻率異常也是明顯大于三極裝置的視電阻率異常，dipole-dipole裝置的分辨率比pole-dipole和dipole-pole裝置的分辨率高。

5.2 算例二

異常體組合模型的三個(gè)電阻率分別是1 000 Ω·m、10 Ω·m和1 000 Ω·m。背景電阻率為均勻介質(zhì)值為100 Ω·m，圖9為起伏地形電阻率模型。在凹地形的一側(cè)埋有1 000 Ω·m的條形高阻體。長(zhǎng)為60 m，寬為10 m，另一側(cè)凸地形同樣埋深條形高阻體，中間設(shè)有一個(gè)10 Ω·m的圓形低阻體。對(duì)比圖11三種裝置的反演結(jié)果圖，從這三個(gè)裝置的反演結(jié)果，可以看出三種布極方式都可以準(zhǔn)確地反映異常體的埋藏深度和形態(tài)。通過(guò)比較不難發(fā)現(xiàn)dipole-dipole裝置反演得到的異常幅值大小以及邊界范圍與真實(shí)模型更為貼近。對(duì)于低阻體的刻畫(huà)相對(duì)于diploe-pole和pole-dipole更為真實(shí)準(zhǔn)確。

圖9 起伏地形電阻率模型

圖10 靈敏度分析

圖11 反演結(jié)果

圖10表明在相同的條件下，dipole-dipole裝置比pole-dipole和dipole-pole裝置對(duì)異常體更靈敏，異常響應(yīng)明顯較大。在近地表對(duì)異常的響應(yīng)，dipole-dipole裝置相對(duì)比pole-dipole和dipole-pole裝置更為集中，主要表現(xiàn)為富集在地表能量的響應(yīng)更強(qiáng)。結(jié)合圖8與圖11可以看出，地形對(duì)三種布極方式的影響會(huì)使得裝置對(duì)地下異常體的形態(tài)變形，從而產(chǎn)生假的異常給后續(xù)的解釋帶來(lái)困難，dipole-dipole裝置反映異常的能力均強(qiáng)于pole-dipole和dipole-pole裝置，且隨著地表場(chǎng)源的位置靠近異常體，dipole-dipole裝置的異常幅值明顯增大，而pole-dipole和dipole-pole裝置的異常幅值相對(duì)較小，在凹地形異常體埋深較淺，三種布極方式對(duì)異常體的形態(tài)都可以較為明顯刻畫(huà)出來(lái)，在凸地形中由于異常體的埋深較深，dipole-dipole相對(duì)于pole-dipole和dipole-pole表現(xiàn)更為出色，對(duì)起伏地形中低阻異常體也可以看出，dipole-dipole具有更好的分辨率以及對(duì)異常體形態(tài)的刻畫(huà)更為準(zhǔn)確。

6 結(jié)論

基于前人的基本理論，利用有限體積法對(duì)起伏地形的直流電阻率法進(jìn)行正演，計(jì)算出大小相同，位置、電阻率不同的高阻和低阻圓形異常體的視電阻率值。利用加權(quán)正則化Tikhonov反演方法實(shí)現(xiàn)起伏地形2.5D電阻率反演成像算法，模型計(jì)算結(jié)果表明，該反演算法較穩(wěn)定、快速，反演結(jié)果良好；在正演算法中，利用多線(xiàn)程并行計(jì)算大型稠密矩陣，為反演算法中快速回代求解提供了條件；反演算法中構(gòu)建了帶地形和高、低阻異常體的二層模型，實(shí)現(xiàn)了靈敏度加權(quán)的2.5D反演方法，研究了直流電阻率法三種布極方式對(duì)反演結(jié)果的影響，得出結(jié)論不同的布極方式對(duì)地下導(dǎo)體反演的效果不同，三種布極方式反演結(jié)果表明，dipole-dipole裝置對(duì)地下異常體的響應(yīng)更為明顯，探測(cè)的能力更好，分辨率也同樣比dipole-pole、pole-dipole兩種布極方式要好。