自適應網格直流電阻率三維反演

周 勇 張志勇* 張大洲 謝尚平 馬 鑫 余鵬洲

(①東華理工大學地球物理與測控技術學院，江西南昌 330013； ②中南大學地球科學與信息物理學院，湖南長沙 410083)

0 引言

直流電阻率法是資源和工程勘探中廣泛應用的一種地球物理方法[1-5]。有限單元法能夠精細刻畫地下真實的復雜地質體，且計算精度高，因而受到廣泛關注[6-11]。區域剖分是有限單元法計算的重要環節，直接影響計算結果。傳統的結構化網格難以解決構造復雜模型邊界及解決源場的奇異性問題，往往造成較大的視電阻率誤差[12-13]。為獲得更精確的結果，必須以增加網格為代價，這無疑會成倍地增加計算量及內存需求。非結構化網格通過網格加密技術在關注區域(如場源、觀測點、異常體、地形)進行網格細化，可以有效地處理異常體邊界及場源點的奇異性問題，進而得到更為精確的解[14]。

采用有限階精度的插值函數進行有限單元正演時，網格越精細得到的結果往往越精確，然而這樣無疑增加了反演的不適定性。通過網格數量與質量的控制可以有效增加反演的穩定性[15]。在地震彈性波層析成像中，基于單元所通過的射線數量進行網格優化的技術得到了應用[16]； 通過組合具有相似物性的鄰近單元以減少單元數量的技術，提高了醫學中電阻率圖像重建的質量[17]； 非線性反演中，基于靈敏度信息的多尺度自適應優化方法很好地避免了總目標函數陷入局部極小值的風險[18]。上述“粗化”類網格優化方法旨在通過減少單元數以減少反演不適定性，按照一定標準“細化”并得到高質量網格是提高反演穩定性的另一種思路。物性梯度可以反映異常體邊界，以物性梯度為依據逐步細化網格，能更好地反演出異常體形狀[19]； 將精度逐漸增加的三重網格分別用于二次場和總場模擬以及反演，在保證正演精度的情況下有效提高了反演效率[20]。此外，小波多尺度變換技術也受到了關注[21]。

自適應網格技術在地球物理反演中逐漸成為熱點，但網格優化缺少評價標準，因此本文研究了依據模型靈敏度的自適應網格生成技術，并開展了三維直流電阻率反演。模型靈敏度反映模型改變量對于所有觀測數據集的影響，若較小的模型改變量造成了較大觀測數據的改變時，則應對該反演區域進行網格優化，同時通過設置最小單元體積以避免產生增加反演不確定性的過小網格。該方法在反演中構建了非結構化網格條件下的最小結構穩定因子；使用高斯—牛頓法最優化求解目標函數； 通過穩定雙共軛梯度法(BICGSTAB)迭代求解高斯—牛頓系統。理論與實測數據反演證明了本文方法的有效性。

1 直流電阻率法三維正演

直流電法的邊值問題可以表示為[8]

(1)

式中：u為電位；r表示位置矢量；r0為電流源位置；δ為狄拉克函數；σ為電導率；I為電流強度；Γ0表示地表界面；ΓR表示無窮遠處邊界；n為邊界外法線方向矢量； 〈·，·〉表示兩個矢量的夾角。邊值問題等價于下式泛函極值問題[5]

(2)

采用有限單元法[22]求解式(2)的極值問題

(3)

KU=P

(4)

求解式(4)可得節點電位，進而計算視電阻率。

2 非結構化自適應網格反演

2.1 正則化反演

正則化是求解不適定性地球物理反問題的重要手段，Tikhonov正則化公式[23]為

Φ=φ(d)+αφ(m)

(5)

式中：Φ為總目標函數，反演就是求解Φ極小值的過程；φ(d)為數據誤差函數，d為數據向量;φ(m)為穩定因子，m為模型參數向量； α為正則化因子，是權衡數據擬合誤差與穩定因子的一個折中系數。α的選取方式，一般包括經驗定值法、無偏風險估計預測方法(UPRE)[24]、廣義交叉驗證方法(GCV)[25-26]、 “L”型曲線法[27]、自適應求取方法[28]等。本文采用一種修正的“L”型曲線方法[29]。

目標函數的數據擬合泛函為

(6)

式中：Wd為數據權重矩陣；A為正演算子。穩定因子能對模型解的空間進行限制，以減少多解性，常用的穩定因子有平滑模型約束穩定因子[30](一階導數、二階導數)、最小結構穩定因子[31-32]、最小支持穩定因子[33-34]和最小梯度支持穩定因子[35-36]等。本文采用最小結構穩定因子

(7)

式中：αS、αx、αy、αz為比例系數；mapr為先驗模型。

Zhdanov[37]給出穩定因子的一般形式為

(8)

式中Wm為模型權重矩陣。由此可以得到總目標函數

(9)

2.2 基于高斯—牛頓的最優化解

高斯—牛頓法求解最優化目標函數[38-39]的模型改進量

Δm=mk+1-mk

(10)

求解目標函數對Δm的一階導數，并令其為0，得到

(11)

(12)

(13)

則式(11)可寫成高斯—牛頓方程

(14)

2.3 BICGSTAB算法

BICGSTAB算法是一種基于雙邊Lanczos算法和殘差正交子空間的迭代方法，在Krylov子空間中選取序列，并選擇合適的搜索方向及步長，使誤差最小。該算法收斂規則、平滑且只需要較小的迭代次數，在求解大型線性方程組時優勢明顯[42]。具體過程參見文獻[43]。

2.4 非結構化網格粗糙度計算

Sethian等[44]提出，單元的空間梯度可以通過方向梯度的線性組合得到。在二維條件下，空間梯度通過計算單元中心點與其相鄰單元中心點之間的方向導數線性組合得到。拓展到三維，在非結構化網格中任選四個相鄰單元，假設其中心點(xi,yi,zi)的某一物性參數為ωi，則中心單元的空間梯度等于ω0對相鄰四個單元中心方向上的導數之和(圖1)。

圖1 非結構化網格單元

ωi的值可以通過一個與x、y、z有關的線性函數擬合

ω(x,y,z)=ax+by+cz+d

(15)

5個單元形成5個方程組

(16)

式(16)可寫為矩陣形式

Fq=ω

(17)

式中

q=(FTF)-1FTω

使用最小二乘方法得到超定方程的解，中心單元中任意一點p0處的梯度為

(18)

2.5 自適應網格的生成

設N為數據點數，定義模型靈敏度為

Sj是模型參數mj改變對整體數據影響的綜合響應，反之也反映了整個數據集對mj改變的識別能力。不考慮數據方差與正則化，Sj實質上是矩陣H的主對角元素的平方根。一般而言，靈敏度矩陣的奇異值與Hessian的特征值是網格質量評判的依據。顯然，如果Hessian是主對角占優的矩陣，則有利于反問題的優化求解[15]。基于Sj進行網格優化實質上是對Hessian矩陣的主對角元素進行優化，所以通過Sj對網格進行優化可以生成適用于反演的高質量網格，計算流程如圖2，運算過程如下：

(1)生成初始網格，設置優化比例、最小單元；

(2)計算J；

(3)分析向量S的平均值方差，并按比例選擇待優化單元；

(4)設置優化單元質量標準；

(5)判斷網格質量，不滿足要求則返回步驟(2)，滿足要求則退出。

圖2 網格優化流程

3 模型驗證

設地下15m處有一體積為20m×20m×20m的低阻異常體(圖3)，電阻率為10Ω·m，地層背景電阻率為100Ω·m，地面布設15×15個接收電極，極距為5m，電極按E-SCAN[4]觀測方法布設。

圖3 理論模型

比較圖4a與圖4b可見，在場源及異常體周邊進行了網格細化，這些區域內很小的模型改變量往往會造成模型靈敏度大幅度的改變； 在距離異常體較遠的區域并未進行網格加密，因為這些區域內的模型靈敏度值往往較低，無需優化。隨著距離的增加，網格逐漸粗糙的過程也更符合真實的地電模型。

圖5所示為不同網格條件下的反演結果及反演迭代過程中相關因素變化曲線(右)對比。

圖5a為未進行網格優化的反演結果，由于網格過于粗糙，對異常體邊界刻畫不清晰，導致異常體形態偏離真實模型；場源附近過于粗糙的網格導致了嚴重的奇異性問題，具體表現在場源附近反演出多個低阻假異常；反演電阻率大于理論模型。圖5b為采用本文所述網格優化策略進行一次優化得到的結果，可見對模型靈敏度較大的區域進行網格細化后，地表附近未見明顯異常，反演的異常體位置及電阻率更接近真實模型。從表1可見，網格優化一次后，均方根誤差明顯下降，從3.60降至2.38。兩次優化后的反演結果(圖5c)顯示，第二次優化網格過程中，網格增加數量較第一次少，說明網格優化趨于穩定。

因此，反演結果在一次優化的基礎上得到了進一步改善，表現為邊界的分辨率更高，異常體形態也更接近真實情況。

使用粗網格進行反演時(圖5a)，迭代前期數據誤差迅速下降，迭代中、后期出現穩定因子與數據誤差同時上升的情況。分析可能的原因是網格過于粗糙，較低的正演精度直接影響反演結果，導致迭代中后期數據誤差隨正則化因子的減小而增大。該方案由于網格數量較少，反演所需時間較短，但反演效果不理想。運用本文所述網格優化策略進行兩次網格細化(圖5c)則有效改善了這一問題。隨著迭代進行，均方根誤差隨穩定因子的增加而下降，在迭代后期逐步趨于穩定。網格的細化意味著計算量的增加，通過設置網格細化閾值避免了產生過多過細網格，以兼顧計算效率。

表1 不同網格條件下反演次數及誤差

圖4 網格優化前(a)、后(b)對比

圖5 不同網格條件下的反演結果(左)及反演迭代過程中相關因素變化曲線(右)對比

4 實測數據反演

將本文所述算法運用于某區巖溶探測數據的反演。布設南北、東西方向測線各10條(圖6)。采用溫納和溫納—施倫貝爾[4]兩種工作裝置共采集有效測點12775個。

根據數據反演結果，得到如下的結論：反演區域發育3個低阻異常體(圖7)。①號異常體推測為被低阻物質填充的完整溶洞； ②號異常體部分出露地表，底部深度為-40m，推測為被低阻物質填充的、出露地表的漏斗狀溶洞； ③號異常體同樣出露地表，底部深度為-23m，推測為被低阻物質填充的漏斗狀溶洞。鉆探記錄顯示這3個異常體是被低阻物質填充的溶洞反映，異常體周圍的高阻體是致密灰巖。

圖6 直流電阻率勘探測線部署圖

圖7 實測數據反演結果

綜上分析，利用本文所述算法結合地質資料，大致確定地下3個被低阻物質填充的溶洞，其規模和形態與實際鉆井信息吻合，證實了本文所述算法的正確性。

5 結束語

本文采用基于模型靈敏度的網格優化算法，生成了高質量的反演網格，進而有效提高了反演的效果。將本文算法運用于理論模型及實際地下溶洞探測，均取得了理想的反演結果。基于模型靈敏度的網格優化算法從反演需求出發進行模型剖分，為高效、精確地進行三維反演提供了一種新思路。