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

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

0 引言

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

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

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

1 直流電阻率法三維正演

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

(1)

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

(2)

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

(3)

KU=P

(4)

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

2 非結(jié)構(gòu)化自適應(yīng)網(wǎng)格反演

2.1 正則化反演

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

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

(5)

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

目標函數(shù)的數(shù)據(jù)擬合泛函為

(6)

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

(7)

式中：αS、αx、αy、αz為比例系數(shù)；mapr為先驗?zāi)Ｐ汀?/p>

Zhdanov[37]給出穩(wěn)定因子的一般形式為

(8)

式中Wm為模型權(quán)重矩陣。由此可以得到總目標函數(shù)

(9)

2.2 基于高斯—牛頓的最優(yōu)化解

高斯—牛頓法求解最優(yōu)化目標函數(shù)[38-39]的模型改進量

Δm=mk+1-mk

(10)

求解目標函數(shù)對Δm的一階導(dǎo)數(shù)，并令其為0，得到

(11)

(12)

(13)

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

(14)

2.3 BICGSTAB算法

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

2.4 非結(jié)構(gòu)化網(wǎng)格粗糙度計算

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

圖1 非結(jié)構(gòu)化網(wǎng)格單元

ωi的值可以通過一個與x、y、z有關(guān)的線性函數(shù)擬合

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

(15)

5個單元形成5個方程組

(16)

式(16)可寫為矩陣形式

Fq=ω

(17)

式中

q=(FTF)-1FTω

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

(18)

2.5 自適應(yīng)網(wǎng)格的生成

設(shè)N為數(shù)據(jù)點數(shù)，定義模型靈敏度為

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

(1)生成初始網(wǎng)格，設(shè)置優(yōu)化比例、最小單元；

(2)計算J；

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

(4)設(shè)置優(yōu)化單元質(zhì)量標準；

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

圖2 網(wǎng)格優(yōu)化流程

3 模型驗證

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

圖3 理論模型

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

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

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

因此，反演結(jié)果在一次優(yōu)化的基礎(chǔ)上得到了進一步改善，表現(xiàn)為邊界的分辨率更高，異常體形態(tài)也更接近真實情況。

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

表1 不同網(wǎng)格條件下反演次數(shù)及誤差

圖4 網(wǎng)格優(yōu)化前(a)、后(b)對比

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

4 實測數(shù)據(jù)反演

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

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

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

圖7 實測數(shù)據(jù)反演結(jié)果

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

5 結(jié)束語

本文采用基于模型靈敏度的網(wǎng)格優(yōu)化算法，生成了高質(zhì)量的反演網(wǎng)格，進而有效提高了反演的效果。將本文算法運用于理論模型及實際地下溶洞探測，均取得了理想的反演結(jié)果?；谀Ｐ挽`敏度的網(wǎng)格優(yōu)化算法從反演需求出發(fā)進行模型剖分，為高效、精確地進行三維反演提供了一種新思路。