基于二維數(shù)據(jù)的加權(quán)非線性共軛梯度三維反演

韓 雪，強建科，魯 凱，李俊營，滿開峰，毛先成

?

基于二維數(shù)據(jù)的加權(quán)非線性共軛梯度三維反演

韓 雪1，強建科2,3，魯 凱2，李俊營2，滿開峰2，毛先成2,3

( 1.廣東省地質(zhì)物探工程勘察院，廣東 廣州 510800；2.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083；3.中南大學(xué) 有色金屬成礦預(yù)測教育部重點實驗室，湖南 長沙 410083)

針對當(dāng)前電阻率數(shù)據(jù)的解釋仍以二維反演為主，而地下為三維結(jié)構(gòu)的情況，研究起伏地形條件下三維源二維采集方式獲取的多條測線測深數(shù)據(jù)的三維反演問題。正演算法采用三棱柱剖分的有限單元法，可以模擬起伏地形下的復(fù)雜模型；在求解大型線性方程組時采用了不完全喬利斯基分解回代技術(shù)，大大提高了三維多點電源的正演計算速度；在計算偏導(dǎo)數(shù)矩陣過程中，依據(jù)互換原理，只需有效組合正演時每個點電源所對應(yīng)的節(jié)點電位便可得到偏導(dǎo)數(shù)矩陣，節(jié)約大量計算時間。其次三維反演算法采用加權(quán)正則化共軛梯度法，抑制了嚴(yán)重病態(tài)的非線性問題。模型計算結(jié)果表明：加權(quán)非線性共軛梯度反演算法既保證了反演穩(wěn)定收斂，又能適合起伏地形反演，達到了滿意的效果。

電阻率三維反演；起伏地形；共軛梯度反演；二維數(shù)據(jù)三維反演

1 引 言

直流電阻率法在探測巖溶、采空區(qū)等方面應(yīng)用廣泛[1,2]，然而在實際勘查中經(jīng)常會遇到地形起伏情況，實測數(shù)據(jù)會發(fā)生嚴(yán)重的畸變，這對后期數(shù)據(jù)的反演解釋工作帶來了非常大的困難。對于三維反演，偏導(dǎo)數(shù)矩陣的計算時間和計算精度，直接影響反演計算的時間和精度，因此，偏導(dǎo)數(shù)矩陣求取是反演問題的關(guān)鍵之一。在電阻率一維反演中，偏導(dǎo)數(shù)矩陣常常采用差分法近似獲得，由于模型參數(shù)少，即使多次調(diào)用正演程序，總體速度仍然很快，但對于二、三維反演來說，模型參數(shù)很多，差分法顯然是不合適的。為了快速方便地求出精度較高的偏導(dǎo)數(shù)矩陣，國內(nèi)外的地球物理學(xué)家進行了不懈的努力。Tripp等(1984)[3]提出了一種利用電位函數(shù)與模型參數(shù)間的簡單關(guān)系來計算偏導(dǎo)數(shù)的近似方法，當(dāng)時是用在二維地電斷面的反演問題中，取得了很好的效果。以后，很多地球物理學(xué)家的偏導(dǎo)數(shù)計算大都采用這個辦法。阮百堯等(1999，2001)[4,5]對二維偏導(dǎo)數(shù)矩陣的求取做了詳細論述；黃俊革(2003)[6]利用互換定律在三維方面的應(yīng)用作了研究，提高了計算速度；Papadopoulos等(2011)[7]改進了雅可比矩陣的計算量，去除對反演沒有重要貢獻的參數(shù)，在保證精度的基礎(chǔ)上，減少內(nèi)存，提高反演速度。本文采用互換定律來計算偏導(dǎo)數(shù)矩陣。

Park等(1991)[8]發(fā)表了有限差分的三維反演算法，開始了比較系統(tǒng)的反演方法研究。Sasaki(1994)[9]闡述了基于有限單元法的三維電阻率反演，但對深部單元的分辨能力不強，反演得到的電阻率值與實際模型有較大差別。Zhang等(1995)[10]使用共軛梯度算法對三維模型進行了正反演研究。吳小平等(1999，2000，2005)[11-13]、劉海飛等(2011)[14]研究了E-SCAN測量方式在非平坦地形上的三維正反演問題，取得了一些好的效果，但在實際工程勘探中，受成本和周期的影響，獲取三維勘探原始數(shù)據(jù)存在較大困難，而且E-SCAN測量方式是基于二極電位測量裝置的，分辨率低于梯度裝置，目前野外勘探仍然以多剖面的二維測量為主。黃俊革等(2004)[15]對三維電阻率反演作了改進，提出“體積因子”、將全局反演改為局部反演等措施，提高了反演的速度和精度。強建科等(2007，2013)[16,17]提出三棱柱單元剖分的有限元方法可以模擬大傾角的復(fù)雜地形，同時改進了偏導(dǎo)數(shù)矩陣，解決了反演深度偏淺的弊病。Gunther等(2006)[18]實現(xiàn)了任意地形的三維電阻率反演，正演采用的是非結(jié)構(gòu)化網(wǎng)格技術(shù)的有限單元法。Oldenborger等(2009)[19]利用點擴散函數(shù)實現(xiàn)了三維電阻率成像反演。

本文研究的是目前野外勘探常采用的三維源二維采集方式獲取的三維數(shù)據(jù)，依此數(shù)據(jù)進行三維起伏地形反演算法。

2 三維反演思路

研究反演的第一步是正演問題。本文采用三維起伏地形有限元正演算法[16]，地下空間采用三棱柱(五面體)單元剖分，這種單元可以模擬起伏地形，插值函數(shù)選為三線性插值函數(shù)，外邊界采用混合邊界[20]。

在實現(xiàn)三維反演中，采用固定網(wǎng)格大小，只考慮視電阻率對地下各個網(wǎng)格電導(dǎo)率的響應(yīng)關(guān)系，可以減少未知參數(shù)數(shù)量，便于求解。

在通常情況下，直流電阻率法的正演問題可以表示成如下形式，即

A(m)=d

(1)

式(1)中：A是非線性正演算子，這里用有限元法實現(xiàn)；m是模型參數(shù)電導(dǎo)率，單位為S/m；d是觀測數(shù)據(jù)，即視電阻率值，單位為Ω·m。

而反演過程就是根據(jù)非線性正演算子A和觀測數(shù)據(jù)矢量d求出模型參數(shù)m。

m=A-1d

(2)

由于式(2)是病態(tài)的，即解是非惟一和不穩(wěn)定的，因此要根據(jù)正則化原理求解病態(tài)問題，很多學(xué)者將正則化反演方法應(yīng)用于直流電的三維反演[21-27]，目標(biāo)函數(shù)方程如下：

Pα(m,d)=φ(m)+αS(m)=min

(3)

rn=A(mn)-d

(4)

(5)

(6)

(7)

(8)

3 偏導(dǎo)數(shù)矩陣

在共軛梯度求解中，需要知道偏導(dǎo)數(shù)矩陣，即視電阻率對模型單元電導(dǎo)率的靈敏度矩陣

F=[aij]

其中i=1,2,…,Ns,

j=1,2,…,Mm

(9)

式(9)中，ρsi是實測得到的視電阻率數(shù)據(jù)，共有Ns個；σj是第j個網(wǎng)格電導(dǎo)率，共有Mm個模型參數(shù)。

以上偏導(dǎo)數(shù)矩陣計算，可歸結(jié)到網(wǎng)格節(jié)點電位對模型中各個網(wǎng)格電導(dǎo)率的偏導(dǎo)數(shù)矩陣計算問題。

例如，對于偶極-偶極裝置來說，視電阻率計算公式為

(10)

對上式兩端求電導(dǎo)率偏導(dǎo)數(shù)，有

(11)

由此看出，任何排列視電阻率對模型電導(dǎo)率的偏導(dǎo)數(shù)都可用點電流源的電位對電導(dǎo)率的偏導(dǎo)數(shù)來表達。

在正演計算時，有限元的單元集成方程為

KU=P

(12)

式(12)中，K為系數(shù)矩陣，U為單元節(jié)點電位，P為電源項和邊界條件。

在方程(12)兩邊，對地下節(jié)點K標(biāo)識的網(wǎng)格單元電導(dǎo)率σK求偏導(dǎo)數(shù)，由于場源項P與地下電導(dǎo)率無關(guān)，故求導(dǎo)結(jié)果為

(13)

比較方程(12)和(13)，可將方程(13)理解為，地面節(jié)點電位U對任意網(wǎng)格單元電導(dǎo)率σK的偏導(dǎo)數(shù)，可以由方程(13)右端項(新場源項)求得，而該新場源項僅與該單元(三棱柱)矩陣系數(shù)和6個節(jié)點電位有關(guān)，因此，把新場源項可以看作位于這6個節(jié)點上的6個微小點電流源：

(14)

式(14)中Kij(i=1，2，…，6；j=1，2，…，6)為三棱柱單元的剛度矩陣之諸元素。當(dāng)對求解區(qū)的網(wǎng)格進行剖分時，所用步長是相等的，即求解區(qū)采用均勻剖分，這時，式(14)的系數(shù)矩陣對于所有單元都是相同的。這給計算帶來了方便。

(15)

將(14)式代入(15)式，有：

(16)

有了地面節(jié)點的電位對模型電導(dǎo)率的一階偏導(dǎo)數(shù)，按照(11)式可以組合出偶極裝置的視電阻率對電導(dǎo)率的偏導(dǎo)數(shù)矩陣F，即雅可比矩陣。獲得雅可比矩陣F后，就可進行反演成像計算了。

從上述計算步驟可以看出，在形成雅可比矩陣的過程中，其計算量是很小的，因為上述步驟中用到的不同供電點對應(yīng)的所有網(wǎng)格節(jié)點電位，都來自于正演計算中保留下來的；單元剛度矩陣數(shù)據(jù)也是正演時保存的。計算雅可比矩陣只是利用了正演的計算結(jié)果，通過換算得到的。當(dāng)然，如果只需地面節(jié)點電位的正演結(jié)果，不需要保存整個斷面空間的電位函數(shù)值。因此，為了計算雅可比矩陣需要更多的計算機內(nèi)存空間以保存每一個地面電極供電時的斷面空間的電位值，而這對于現(xiàn)在的計算機條件來說，是可以做到的。也就是說，現(xiàn)在只需作一次正演就可以得到雅可比矩陣，使得反演所需的時間大大地減少了，這種求雅可比矩陣的方法使得高密度電法的三維反演成像成為可能，為保障高密度電法的地質(zhì)效果提供了條件。

4 三維反演算例

下面的算例中均采用偶極—偶極裝置。

模型1：簡單模型如圖1(a)所示，均勻半空間中有一個4 m×4 m×4 m的低阻立方體模型，電阻率為10 Ω·m， 頂埋深為3 m， 中心埋深5 m，圍巖電阻率為100 Ω·m，測量點距1 m，線距1 m，每條剖面上25個供電點，共設(shè)計5條剖面。通過正演獲得“實測數(shù)據(jù)”。在反演中，初始模型電阻率取“實測視電阻率數(shù)據(jù)”的平均值作為背景電阻率。

圖1(b)是偶極-偶極裝置三維正演垂直切片圖，可以看出5條視電阻率等值線異常呈標(biāo)準(zhǔn)的“八字形”低阻特點，且前后左右對稱。在低阻體正上方三條測線(Y=-1 m、Y=0 m、Y=1 m)視電阻率差異較大，兩邊測線(Y=-2 m、Y=2 m)視電阻率差異變小。

圖1(c)是三維反演結(jié)果，經(jīng)過約1小時8次迭代，擬合誤差小于5 %，這5條斷面圖很好地還原了低阻體的大小、形狀和埋深，說明該算法的收斂性很好。

模型2：復(fù)雜模型如圖2(a)所示，均勻半空間中有一個4 m×4 m×4 m的高阻立方體模型和一個長寬高為6 m×4 m×6 m低阻長方體模型，高阻立方體電阻率為500 Ω·m，頂埋深為2 m，低阻立方體電阻率為10 Ω·m，頂埋深為4 m，兩個異常相隔4 m，圍巖電阻率為100 Ω·m，測量點距1 m，線距1 m，每條剖面上31個供電點，共設(shè)計5條剖面。模擬“實測數(shù)據(jù)”由三維正演得出。初始模型電阻率取“實測視電阻率數(shù)據(jù)”的平均值作為背景電阻率。

圖2 復(fù)雜模型三維反演結(jié)果Fig.2 3D inversion result of complex model

圖3 起伏地形低阻體三維正反演結(jié)果Fig.3 The result of forward calculation and inversion under undulating terrain 3D low resistance body model

對上述復(fù)雜模型模擬的視電阻率數(shù)據(jù)進行三維反演，評估該算法對高低阻組合模型的分辨能力。圖2(b)是三維反演結(jié)果水平切片，沿縱深方向取4個切片，從圖中可看出，在Z=3 m切片上，左側(cè)有一個高電阻率異常，其異常中心與真實模型一致；在Z=5 m、7 m切片中，右側(cè)有一個低電阻率異常，其異常中心與真實模型一致。由此得出，該反演算法能夠從復(fù)雜的組合視電阻率異常中反演出不同深度的高、低阻體，成像效果明顯。

圖2(c)是三維反演結(jié)果縱向切片，從5條斷面圖上看出，在X方向上，該算法能夠清楚地識別間隔4 m的兩個高低電阻率異常，且異常中心與真實模型一致，成像效果很好。圖2(d)對Y=0 m斷面和Z=5 m切片上實測數(shù)據(jù)與預(yù)測模型數(shù)據(jù)進行對比，二者視電阻率等值線形態(tài)一致，說明預(yù)測模型與真實模型非常接近。

模型3：在模型1的基礎(chǔ)上增加一個三維山脊地形，如圖3(a)所示，坡度約30°，斜坡長8 m，山脊高4 m，其他參數(shù)同模型1。

圖3(b)為5條測線上的視電阻率等值線斷面圖，從圖上看出，三維山脊對視電阻率異常產(chǎn)生了非常嚴(yán)重的干擾，致使山脊下視電阻率值變得很大，其視電阻率為圍巖電阻率的3倍多，兩側(cè)為低阻異常。與圖1(b)比較看出，山脊產(chǎn)生的虛假高視電阻率異常完全破壞了低阻體產(chǎn)生的“八字形”視電阻率異常特征。

圖3(c)為起伏地形三維反演結(jié)果，從圖中看出，反演算法正確反演出了低阻體的位置、形狀和大小，反演結(jié)果與地下真實模型相符，沒有得出多余的構(gòu)造。由此得出，直接帶入地形進行三維反演才能夠有效地解決由地形產(chǎn)生的虛假異常，從而得到地下真實的電性分布特征，同時證明該算法反演成像穩(wěn)定可靠。

5 結(jié) 論

偏導(dǎo)數(shù)矩陣求解是電阻率三維反演的一個關(guān)鍵問題，根據(jù)穩(wěn)定電流場的互換原理，推導(dǎo)出含有起伏地形信息的雅可比矩陣公式，通過組合各節(jié)點的電位，能夠快速得到雅可比矩陣，節(jié)約大量計算時間。在前人研究的基礎(chǔ)上，提出加權(quán)正則化非線性共扼梯度方法，優(yōu)化了權(quán)重系數(shù)和正則化因子，從而實現(xiàn)起伏地形電阻率三維反演算法，算例表明該方法成像效果良好。在求解大型線性方程組時采用了不完全喬利斯基分解回代技術(shù)，大大提高了三維多點電源的正演計算速度問題。

[1]肖敏,陳昌彥,白朝旭,等.北京地區(qū)淺層采空區(qū)高密度電法探測應(yīng)用分析[J].工程地球物理學(xué)報,2014,11(1):29-35.

[2]袁忠明,韋乙杰.高密度電阻率法偶極裝置在某道路工程巖溶探測中的應(yīng)用[J].工程地球物理學(xué)報,2015,12(1):72-76.

[3] Tripp A C, Hohmann G W, Swift C M. Two-dimensional resistivity inversion[J].Geophysics,1984,49 (10):1 708-1 717.

[4]阮百堯,村上裕,徐世浙.電阻率激發(fā)極化法數(shù)據(jù)的二維反演程序[J].物探化探計算技術(shù),1999,21(2): 116-125.

[5]阮百堯.視電阻率對模型電阻率的偏導(dǎo)數(shù)矩陣計算方法[J].地質(zhì)與勘探,2001,37(6):39-41.

[6]黃俊革.三維電阻率/激化率有限單元正演模擬與反演成像[D].湖南:中南大學(xué),2003.

[7]Papadopoulos N G, Tsourlos P, Papazachos C, et al.An algorithm for fast 3D inversion of sur-face electrical resistivity tomography data: application on imaging buried antiquities[J]. Geophysical Prospecting,2011,59(3):557-575.

[8]Park S K, Van G P.Inversions of pole-pole data for 3-D resistivity structure beneath arrays of electrodes[J].Geophysics,1991,56(7):951-960.

[9]Sasaki Y.3-D resistivity inversion using the finite-element method[J]. Geophysics,1994,59(12): 1 839-1 848.

[10]Zhang J, Mackie R L, Madden T R. 3-D resistivity forward modeling and inversion using conjugate gradients[J]. Geophysics,1995,60(5):1 313-1 325.

[11]吳小平,徐果明.電阻率三維反演中偏導(dǎo)數(shù)矩陣的求取與分析[J].石油地球物理勘探,1999,34(4): 363-373.

[12]吳小平,徐果明.利用共軛梯度法的電阻率三維反演研究[J].地球物理學(xué)報,2000,43(3):420-427.

[13]吳小平.非平坦地形條件下電阻率三維反演[J].地球物理學(xué)報,2005,48(4):932-936.

[14]劉海飛,柳建新,郭榮文,等.起伏地形三維激電連續(xù)介質(zhì)模型快速反演[J].吉林大學(xué)學(xué)報(地球科學(xué)版),2011,41(4):1 212-1 218.

[15]黃俊革,阮百堯,鮑光淑.基于有限單元法的三維地電斷面電阻率反演[J].中南大學(xué)學(xué)報(自然科學(xué)版), 2004,35(2):295-299.

[16]強建科,羅延鐘.三維地形直流電阻率有限元法模擬[J].地球物理學(xué)報,2007,50(5):1 606-1 613.

[17]Qiang Jianke, Han Xue, Dai Shikun.3D DC Resistivity Inversion with Topography Based on Regularized Conjugate Gradient Method[J]. International Journal of Geophysics,2013.Article ID 931876, 9 pages, 2013. doi:10.1155/2013/931876.

[18]Günther T, Rücker C, Spitzer K. Three-dimensional modelling and inversion of dc resistivity data incorporating topography-II. Inversion[J]. Geophys.J.Int.,2006,166(2):506-517.

[19]Oldenborger G A, Routh P S. The point-spread function measure of resolution for the 3-D electrical resistivity experiment[J]. Geophys.J.Int.,2009,176(2):405-414.

[20]Dey A, Morrison H F. Resistivity modeling for arbitrarily shaped three-dimensional structures [J]. Geophysics,1979,44(4):753-780.

[21]Ellis R G, Oldenburg D W.T he pole-pole 3-D DC resistivity inverse problem: a conjugate gradient approach[J]. Geophys.J.Int.,1994,119(1):187-194.

[22]LaBrecque D J, Miletto M, Daily W, et al. The effects of noise on Occam's inversion of resistivity tomography data[J]. Geophysics,1996,61(2):538-548.

[23]Yi M J, Kim J H, Song Y, et al. Three-dimensional imaging of subsurface structures using resistivity data[J]. Geophysical Prospecting,2001,49(4):483-497.

[24]Pain C C, Herwanger J V, Worthington M H, et al. Effective multidimensional resistivity inversion using finite-element techniques[J]. Geophys.J.Int.,2002,151(3):710-728.

[25]Pidlisecky A, Haber E, Knight R. RESINVM3D: A 3D resistivity inversion package[J]. Geophysics, 2007,72(2):H1-H10.

[26]Marescot L, Palma Lopes S, Green A G. Nonlinear inversion of geoelectric data acquired across 3D objects using a finite-element approach[J]. Geophysics,2008,73(3):F121-F133.

[27]Papadopoulos N G, Yi M J, Kim J H, et al. Geophysical investigation of tumuli by means of surface 3D electrical resistivity tomography[J]. Journal of Applied Geophysics,2010,70(3): 192-205.

[28]Zhdanov M S. Geophysical inverse theory and regularization problems[M].New York:Elsevier, 2002.

The 3D Inversion with Non-linear Conjugate Gradient Based on 2D DC Resistivity Data

Han Xue1，Qiang Jianke2,3，Lu Kai2，Li Junying2，Man Kaifeng2，Mao Xiancheng2,3

( 1.Geologic&GeophysicalEngineeringExplorationInstituteofGuangdong,GuangzhouGuangdong510800,China;2.SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,ChangshaHunan410083,China;3.KeyLaboratoryofMetallogenicPredictionofNonferrousMetals,MinistryofEducation,CentralSouthUniversity,ChangshaHunan410083,China)

In view of the current resistivity data interpretation, the 2D inversion is still taken as the main point, but actually the 3D structure is used more underground. This paper presents a research on the 3D inversion for the multiple line sounding data over rugged terrain area. Triangular-prism mesh finite element method is employed in forward algorithm, which can be used to simulate the complex models with the topography. In solving large linear equations, incomplete Cholesky decomposition and back substitution technique are used, which can significantly improve the speed of 3D multi-source modeling. According to the interchangeable principle, the partial derivative matrix could be acquired by combining the potential of forward simulation, which saves large amount of time. The weighted regularized conjugate gradient is applied to 3D inversion, which solves the severely ill-posed non-linear problem. The model calculation results show that this method could both guarantee the stable inversion convergence and be suitable for rugged topography inversion to reach satisfactory results.

3D resistivity inversion; rugged terrain; conjugate gradient inversion; 3D inversion with 2D data

1672—7940(2016)05—0561—09

10.3969/j.issn.1672-7940.2016.05.001

國家自然科學(xué)基金(編號：41472301；41174104))

韓 雪(1988-)，女，助理工程師，碩士，主要從事電法勘探研究。E-mail:kahanxue@126.com)

強建科(1967-)，男，副教授，博士，主要從事電磁法正反演研究。E-mail:qiangjianke@163.com

P631.3 文獻標(biāo)碼： A

2016-05-28