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

(成都理工大學 地球勘探與信息技術教育部重點實驗室，成都 610059)

0 引言

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

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

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

1 正演算法

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

(1)

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

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

(2)

Aφ=-q

(3)

其中

(4)

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

(5)

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

2 反演算法

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

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

(6)

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

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

(7)

式中：Wm為模型權重，是x、y和z方向上一階平滑度的矩陣。每個Wm矩陣可以用整數(shù)離散表示。正則化Wm是應用為每個單元格具有不同權重的標量或對角矩陣加權總和α*的矩陣。

將數(shù)據(jù)目標函數(shù)寫為式(8)[24]。

(8)

式中：Wd為數(shù)據(jù)的權重矩陣，該矩陣為對角矩陣，其元素為Wdii=1∈i，其中εi是第i個標準偏差；F(m)為正演函數(shù)；m為模型參數(shù)；d為觀測所得的數(shù)據(jù)。

則其梯度為：

(9)

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

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

(10)

求解該擾動，將梯度設為0

(11)

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

圖1 反演并行算法流程圖

3 靈敏度

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

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

(12)

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

dpred=Pu(m)

(13)

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

對于反演問題，主要是解決如何通過改變模型來轉換數(shù)據(jù)。將式(13)泰勒展開得式(14)。

O(h2‖v‖)

(14)

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

J=-PA-1G

(15)

從而可以對偏微分方程(PDE)進行求導

(16)

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

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

(17)

4 正確性驗證

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

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

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

5 算例分析

5.1 算例一

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

圖4 平地形模型

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

圖5 視電阻率擬剖面

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

圖6 視電阻率直方圖

采用靈敏度加權反演后，得到靈敏度如圖7所示。由圖7可以看到，三種布極方式對于淺層有著較高的靈敏度，隨著地層變深而靈敏度降低，當存在高阻體時，靈敏度存在較明顯地降低，而存在低阻體時，靈敏度降低相對周圍圍巖存在明顯地變緩。

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

5.2 算例二

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

圖9 起伏地形電阻率模型

圖10 靈敏度分析

圖11 反演結果

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

6 結論

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