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

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

0 引言

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

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

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

1 直流電阻率法三維正演

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

(1)

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

(2)

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

(3)

KU=P

(4)

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

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

2.1 正則化反演

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

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

(5)

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

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

(6)

式中：Wd為數(shù)據(jù)權(quán)重矩陣；A為正演算子。穩(wěn)定因子能對(duì)模型解的空間進(jì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為先驗(yàn)?zāi)Ｐ汀?/p>

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

(8)

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

(9)

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

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

Δm=mk+1-mk

(10)

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

(11)

(12)

(13)

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

(14)

2.3 BICGSTAB算法

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

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

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

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

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

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

(15)

5個(gè)單元形成5個(gè)方程組

(16)

式(16)可寫為矩陣形式

Fq=ω

(17)

式中

q=(FTF)-1FTω

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

(18)

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

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

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

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

(2)計(jì)算J；

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

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

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

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

3 模型驗(yàn)證

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

圖3 理論模型

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

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

圖5a為未進(jìn)行網(wǎng)格優(yōu)化的反演結(jié)果，由于網(wǎng)格過于粗糙，對(duì)異常體邊界刻畫不清晰，導(dǎo)致異常體形態(tài)偏離真實(shí)模型；場(chǎng)源附近過于粗糙的網(wǎng)格導(dǎo)致了嚴(yán)重的奇異性問題，具體表現(xiàn)在場(chǎng)源附近反演出多個(gè)低阻假異常；反演電阻率大于理論模型。圖5b為采用本文所述網(wǎng)格優(yōu)化策略進(jìn)行一次優(yōu)化得到的結(jié)果，可見對(duì)模型靈敏度較大的區(qū)域進(jìn)行網(wǎng)格細(xì)化后，地表附近未見明顯異常，反演的異常體位置及電阻率更接近真實(shí)模型。從表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ǔ)上得到了進(jìn)一步改善，表現(xiàn)為邊界的分辨率更高，異常體形態(tài)也更接近真實(shí)情況。

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

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

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

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

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

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

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

圖6 直流電阻率勘探測(cè)線部署圖

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

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

5 結(jié)束語

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