基于非結(jié)構(gòu)網(wǎng)格的電阻率三維帶地形反演

吳小平， 劉洋， 王威

1 中國(guó)科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院 地震與地球內(nèi)部物理實(shí)驗(yàn)室， 合肥 230026 2 蒙城地球物理國(guó)家野外科學(xué)觀測(cè)研究站， 蒙城 233500 3 中海石油(中國(guó))有限公司深圳分公司研究院， 廣州 510240

?

基于非結(jié)構(gòu)網(wǎng)格的電阻率三維帶地形反演

吳小平1，2， 劉洋1，2， 王威3

1 中國(guó)科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院 地震與地球內(nèi)部物理實(shí)驗(yàn)室， 合肥 230026 2 蒙城地球物理國(guó)家野外科學(xué)觀測(cè)研究站， 蒙城 233500 3 中海石油(中國(guó))有限公司深圳分公司研究院， 廣州 510240

地表起伏地形在野外礦產(chǎn)資源勘察中不可避免，其對(duì)直流電阻率法勘探影響巨大.近年來(lái)，電阻率三維正演取得諸多進(jìn)展，特別是應(yīng)用非結(jié)構(gòu)網(wǎng)格我們能夠進(jìn)行任意復(fù)雜地形和幾何模型的電阻率三維數(shù)值模擬，但面向?qū)嶋H應(yīng)用的起伏地形下電阻率三維反演依然困難.本文基于非結(jié)構(gòu)化四面體網(wǎng)格，并考慮到應(yīng)用GPS/GNSS時(shí)，區(qū)域地球物理調(diào)查中可非規(guī)則布設(shè)測(cè)網(wǎng)的實(shí)際特點(diǎn)，實(shí)現(xiàn)了任意地形(平坦或起伏)條件下、任意布設(shè)的偶極-偶極視電阻率數(shù)據(jù)的不完全Gauss-Newton三維反演.合成數(shù)據(jù)的反演結(jié)果表明了方法的有效性，可應(yīng)用于復(fù)雜野外環(huán)境下的三維電法勘探.

電阻率； 三維反演； 非結(jié)構(gòu)網(wǎng)格； 地形

1 引言

直流電阻率法是電法勘探的最經(jīng)典方法之一，在金屬礦產(chǎn)、非金屬礦產(chǎn)、煤田、油氣及地下水等資源、能源勘探中獲得廣泛而有效的應(yīng)用，并已拓展到水文、工程、環(huán)境、考古等與國(guó)民經(jīng)濟(jì)建設(shè)、人類社會(huì)生活密切相關(guān)的探測(cè)領(lǐng)域(傅良魁，1983；吳小平和汪彤彤，2002).電法勘探具靈活性及勘探費(fèi)用低的顯著特點(diǎn)，演變出諸如單極、偶極及多極裝置的電剖面法、電測(cè)深方法等多種觀測(cè)方式.特別是20世紀(jì)末發(fā)明的高密度電法儀器，可以高效地獲得大量地電觀測(cè)數(shù)據(jù)，使得地下精細(xì)結(jié)構(gòu)的電阻率三維反演成像成為可能，并迅速成為該領(lǐng)域國(guó)內(nèi)外學(xué)者關(guān)注的前沿課題(Park and Van Gregory，1991；Li and Oldenburg，1994；Dabas et al.，1994；Sasaki，1994；Zhang et al.，1995；Loke and Barker，1996；何繼善，1997；底青云等，1997；吳小平和汪彤彤，2002；黃俊革等，2004).

隨著計(jì)算機(jī)及數(shù)值計(jì)算技術(shù)的發(fā)展，電阻率三維正、反演取得很大進(jìn)展，預(yù)條件共軛梯度、多重網(wǎng)格等算法的引入大大提高了電阻率三維有限差分法、有限單元法數(shù)值模擬的計(jì)算速度和精度(Spitzer，1995；Zhou and Greenhalgh，2001；Li and Spitzer，2002；Wu，2003；Wu et al.，2003；Lu et al.，2010；Ren and Tang，2010；Pan and Tang，2014；Ren and Tang，2014)，為電阻率三維快速反演奠定基礎(chǔ).在電阻率三維反演中，由于反演參數(shù)多和數(shù)據(jù)量大，雅可比矩陣(偏導(dǎo)數(shù)矩陣)巨大的計(jì)算量和存儲(chǔ)需求難于克服.許多學(xué)者研究可避開(kāi)計(jì)算雅可比矩陣的反演算法，如基于Born近似的三維反演方法(Li and Oldenburg，1994)，避開(kāi)了雅可比矩陣計(jì)算，簡(jiǎn)單快速地得到地下介質(zhì)導(dǎo)電率的三維分布.Born近似在地下介質(zhì)導(dǎo)電率變化不大時(shí)比較精確，地下介質(zhì)越復(fù)雜，則精確程度越低，反演分辨率受到很大限制.Zhang等(1995)、吳小平和徐果明(2000)引入共軛梯度方法，解決了電阻率三維有限差分正演效率及三維反演中雅可比矩陣求取、存儲(chǔ)兩方面問(wèn)題，實(shí)現(xiàn)了比較快速、有效的電阻率真三維反演.宛新林等(2005)通過(guò)對(duì)光滑度矩陣元素進(jìn)行適當(dāng)改進(jìn)，使之適用于各種情況下光滑度矩陣的求取，并用最小二乘正交分解法求解反演方程，有效地提高了計(jì)算速度.李長(zhǎng)偉等(2012)基于改進(jìn)Krylov子空間算法實(shí)現(xiàn)了井中激電快速反演；戴前偉等(2012)、劉斌等(2012a)研究了不同約束條件下的電阻率三維反演.Pidlisecky等(2007)采用預(yù)條件雙共軛梯度法求解正演中的大型線性方程組，開(kāi)發(fā)了電阻率三維不完全Gauss-Newton反演程序.電阻率三維反演也逐漸在水文、工程地球物理(Park，1998；Chambers et al.，2011；劉斌等，2012b)、環(huán)境地質(zhì)(Weller et al.，1999；Chambers et al.，2002；Bentley and Gharibi，2004；)、考古(Papadopoulos et al.，2006，2007)以及永久凍土探測(cè)(R?dder and Kneisel，2012；Kneisel et al.，2014)等領(lǐng)域獲得應(yīng)用.

以上所述電阻率三維反演及其應(yīng)用均是基于平坦地形.而實(shí)際勘探中，地形影響不可避免，并對(duì)電阻率反演結(jié)果可能造成難于預(yù)料的偏差.Oppliger(1984)、Holcombe和Jiracek(1984)分別利用積分方程、有限元法進(jìn)行了電阻率三維地形改正，Tsourlos等(1999)研究了二維情況下各種電極排列的地形影響及其改正.他們的結(jié)果表明：做地形改正只能是非常近似的，地下結(jié)構(gòu)稍復(fù)雜，這種改正誤差就很大.只有將地形同時(shí)帶入反演算法中，才能精確消除地形影響及其對(duì)反演結(jié)果的偏差.Yi等(2001)、吳小平(2005)分別利用有限元和有限差分?jǐn)?shù)值模擬實(shí)現(xiàn)了帶地形電阻率三維反演，由于采用結(jié)構(gòu)化的網(wǎng)格單元剖分，不能適應(yīng)于復(fù)雜地形下的電阻率三維反演.近年來(lái)，非結(jié)構(gòu)有限元方法在復(fù)雜地形電阻率三維數(shù)值模擬中取得極大成功(Rücker et al.，2006；Ren and Tang，2010，2014；Wang et al.，2013)，非結(jié)構(gòu)化網(wǎng)格具有單元質(zhì)量可控、允許局部加密、能夠模擬復(fù)雜幾何模型等優(yōu)點(diǎn)，使得三維非結(jié)構(gòu)有限單元求解效率大幅提高，同等計(jì)算精度情況下，非結(jié)構(gòu)化網(wǎng)格相對(duì)于結(jié)構(gòu)化規(guī)則網(wǎng)格，計(jì)算時(shí)間和存儲(chǔ)量均可下降約一個(gè)數(shù)量級(jí).在此基礎(chǔ)上，Günther和Rücker(2006)實(shí)現(xiàn)了電阻率帶地形三維反演，是該領(lǐng)域的一個(gè)重要進(jìn)展.

電阻率三維觀測(cè)方式也是值得研究的重要方面，它不僅關(guān)系到野外觀測(cè)的效率和成本，也關(guān)系到數(shù)據(jù)的分辨力和反演效果(吳小平和汪彤彤，2002).傳統(tǒng)的地面物探往往都需要先進(jìn)行人工物探測(cè)量網(wǎng)布設(shè)，這是一項(xiàng)繁重的野外工作，其效率低、成本高.在山區(qū)、森林覆蓋地區(qū)，通視條件差，國(guó)家控制點(diǎn)少且遠(yuǎn)離工區(qū)，使得用傳統(tǒng)的經(jīng)緯儀來(lái)測(cè)定電法觀測(cè)點(diǎn)的坐標(biāo)幾乎不可能(林君和石磊，1996；陳清禮和胡家華，1998).GPS/GNSS(Global Navigation Satellite System)導(dǎo)航與定位技術(shù)快速發(fā)展(陳俊勇等，2007)，中國(guó)的北斗衛(wèi)星導(dǎo)航系統(tǒng)也日趨成熟(楊元喜，2010)，手持型設(shè)備的定位精度也在不斷提高(劉少杰等，2006；李國(guó)防和閆新亮，2010；王克曉等，2012)，為電法勘探測(cè)網(wǎng)布設(shè)提供了一種全新的工作方法(劉會(huì)敏和趙灝，2013)，在進(jìn)行區(qū)域地球物理調(diào)查時(shí)，可形成適應(yīng)于工區(qū)地形、環(huán)境條件的任意靈活三維測(cè)量電極布設(shè)方式，而不必嚴(yán)格按設(shè)定的規(guī)則測(cè)網(wǎng)順序進(jìn)行.然而，由于任意布極三維測(cè)量的觀測(cè)數(shù)據(jù)結(jié)構(gòu)復(fù)雜，其三維反演面臨諸多挑戰(zhàn)，Günther和Rücker(2006)的電阻率帶地形三維反演沒(méi)有這方面的解決方案，也未見(jiàn)有針對(duì)性的其他工作發(fā)表.

本文將任意布設(shè)電極的三維觀測(cè)方式與電阻率三維非結(jié)構(gòu)有限單元數(shù)值模擬相結(jié)合，實(shí)現(xiàn)任意地形(平坦或起伏)條件下、任意布設(shè)的偶極-偶極視電阻率觀測(cè)數(shù)據(jù)的不完全Gauss-Newton三維反演.這是野外數(shù)據(jù)三維反演解釋的關(guān)鍵，具有顯著的實(shí)際意義.

2 電阻率三維非結(jié)構(gòu)有限單元正演

雙電極供電時(shí)的直流電阻率法控制方程可以表示為：

φ

其中，φ為電勢(shì)，供電電流為I，正電極位置為rs+，負(fù)電極位置為rs-，σ為電導(dǎo)率.用有限元離散求解區(qū)域，可將此式寫成矩陣表達(dá)形式(Wangetal.,2013)：

A(σ)u=q,

(1)

u是表示三維網(wǎng)格節(jié)點(diǎn)上電勢(shì)分布的向量,A表示為拉普拉斯正演算子，實(shí)際是一個(gè)與網(wǎng)格剖分及電導(dǎo)率分布有關(guān)的對(duì)稱、稀疏矩陣，q是一個(gè)描述源分布的向量.

正演問(wèn)題實(shí)際就是在已知區(qū)域中電導(dǎo)率和源分布的情況下，求解電勢(shì)的分布：

u=A-1(σ)q.

(2)

我們實(shí)現(xiàn)了電阻率三維非結(jié)構(gòu)有限單元數(shù)值模擬(Wangetal.，2013)，本文利用三維非結(jié)構(gòu)有限單元可對(duì)任意電阻率模型進(jìn)行正演計(jì)算.圖1a顯示均勻半空間中有一半徑為35m的低阻球體模型及其非結(jié)構(gòu)網(wǎng)格剖分，圖1b為供電極及測(cè)量電極分布的示意圖.點(diǎn)A(30，500)為供電點(diǎn)，21個(gè)測(cè)量電極位于井中測(cè)線z=900～1000 m(x=-70 m)上，間隔10 m.背景電阻率為500 Ωm，球體電阻率為5 Ωm.該球體模型的電位分布有解析解(Large，1971；Zhdanovand Keller，1994)，可用于檢驗(yàn)電阻率三維正演數(shù)值解的可靠性.圖2是三維非結(jié)構(gòu)有限單元數(shù)值模擬結(jié)果與解析解的對(duì)比，平均誤差為0.9%，表明我們的電阻率三維非結(jié)構(gòu)有限單元數(shù)值模擬結(jié)果是可靠的，獲得較高的計(jì)算精度.

圖3是平坦地形下的高阻體模型示意圖，地下高阻異常體電阻率2000 Ωm，異常尺寸20 m×20 m×20 m，頂部埋深10 m，背景電阻率200 Ωm，地表測(cè)區(qū)80 m×80 m，相鄰電極間距離均為4 m，即有21×21個(gè)測(cè)量電極，供電極距AB=100 m.圖4是非結(jié)構(gòu)有限元三維正演計(jì)算獲得的中間梯度視電阻率平面圖，圖中明顯可見(jiàn)高阻異常響應(yīng).

圖1 球體模型(a) 非結(jié)構(gòu)網(wǎng)格剖分; (b) 供電極及測(cè)量電極分布示意圖.Fig.1 Sphere model(a) Unstructured grids; (b) Schematic diagram of electrodes as source and receivers.

圖2 球體模型的數(shù)值解與解析解對(duì)比Fig.2 Comparison between numerical and analytical solutions of the sphere model

實(shí)際測(cè)量中，地形影響不可避免.非結(jié)構(gòu)有限元方法在三維地形的數(shù)值模擬中具有結(jié)構(gòu)化網(wǎng)格無(wú)可比擬的優(yōu)勢(shì).圖5是一個(gè)用凸起半球模擬山峰地形下的高阻異常體及其非結(jié)構(gòu)網(wǎng)格剖分，異常體參數(shù)與圖3中的模型相同.圖6為山峰地形下高阻體模型的中間梯度視電阻率平面圖(AB=100 m)，結(jié)果表明：地形影響非常大，反映的主體響應(yīng)為低阻異常，基本掩蓋了地下為高阻異常體的真實(shí)結(jié)構(gòu).

圖3 高阻體模型示意圖Fig.3 Schematic diagram of the high resistivity model

圖4 高阻體模型的視電阻率平面圖(AB=100 m)Fig.4 Apparent resistivity of the high resistivity model (AB=100 m)

3 電阻率三維不完全Gauss-Newton反演

直流電阻率不完全Gauss-Newton反演中，觀測(cè)數(shù)據(jù)d一般是模型三維剖分網(wǎng)格節(jié)點(diǎn)上電位u的子集，即僅在若干個(gè)有限的位置上進(jìn)行觀測(cè)，有：

d=Qu=QA-1(σ)q,

(3)

其中，Q是u到d的映射矩陣.反演的目的就是從這些有限位置觀測(cè)到的電位d得到地下的電導(dǎo)率σ(或電阻率)分布.通常情況下，為了防止反演不穩(wěn)定而出現(xiàn)偏差比較大的電導(dǎo)率(如不符合物理意義的負(fù)值)，一般取反演參數(shù)m為ln(σ).經(jīng)典的最小二乘法反演是對(duì)以下泛函求極值：

Φd=‖D(d(m)-dobs)‖2,

(4)

一般采用加上模型約束的正則化方法，變成無(wú)約束的最優(yōu)化問(wèn)題，其目標(biāo)泛函為：

(5)

其中，d(m)為計(jì)算數(shù)據(jù)，dobs為觀測(cè)值，β為正則化參數(shù)，mref為參考模型.W為模型加權(quán)矩陣，一般取一階正則化約束，用梯度算子保證模型的光滑度(TikhonovandArsenin，1977)，擴(kuò)展到三維：

圖5 山峰地形下高阻體模型的非結(jié)構(gòu)網(wǎng)格剖分圖Fig.5 Unstructured grids of the high resistivity model under a hill

W= (αsI+αxGx+αyGy+αzGz)T

×(αsI+αxGx+αyGy+αzGz),

其中，Gx，Gy，Gz為三個(gè)方向的梯度算子；αx，αy，αz為三個(gè)方向的光滑度加權(quán)，αs為最小權(quán)重.一般來(lái)說(shuō)W矩陣可根據(jù)先驗(yàn)信息進(jìn)行調(diào)整，比如Pidlisecky等(2007)在W矩陣右端乘上隨深度z遞減的1/z2矩陣，以模擬地表測(cè)量隨深度的衰減.本文不涉及對(duì)W矩陣的再調(diào)整，并取αx=αy=αz=1，αs=0.01.Gauss-Newton方法的法方程形式為：

(JTDTDJ+βWTW)Δm=

-[JTDTD(QA-1q-dobs)+βWTW(m-mref)],

(6)

令

H=JTDTDJ+βWTW，

g=JTDTD(QA-1q-dobs)+βWTW(m-mref),

(6)式可簡(jiǎn)寫為：

HΔm=-g,

(7)

Haber等(2000)將J表示為：

J=-QA-1B,

(8)

其中，B=?(A(m)u)/?m，B是與模型參數(shù)和網(wǎng)格劃分有關(guān)的一個(gè)矩陣.Dembo等(1982)指出，對(duì)方程組(7)只需要近似求解，借鑒不完全Gauss-Newton法，采用不精確的預(yù)條件共軛梯度算法(PCG)求解方程組(7)，即設(shè)置較低的PCG迭代收斂準(zhǔn)則和迭代次數(shù)限制，預(yù)條件矩陣為WTW(Haber，2005).PCG迭代過(guò)程中需要計(jì)算矩陣向量積Hv，只需要處理J與向量v的乘積Jv或JTv即可，可按如下步驟：

(1) 先計(jì)算一個(gè)中間向量w=Bv；

(2) 然后求解線性系統(tǒng)y=A-1w；

(3) 之后y向量左乘投影矩陣Q即可.

由于反演過(guò)程中需要反復(fù)求解三維正演的線性系統(tǒng)y=A-1w來(lái)完成雅可比矩陣與向量的乘積，因而求解Gauss-Newton反演的法方程非常耗時(shí).同樣基于不完全Gauss-Newton思想，在計(jì)算矩陣向量積Hv時(shí)，設(shè)置較低的正演計(jì)算精度.JTv的計(jì)算類似，詳見(jiàn)相關(guān)文獻(xiàn)(吳小平和徐果明，1999；李長(zhǎng)偉等，2012)，大大提高了電阻率三維不完全Gauss-Newton反演算法的效率.

引入上述方法之后，簡(jiǎn)要的不完全Gauss-Newton反演流程如下：

(1) 給定初始模型和參考模型，均取為均勻半空間.

(2) 利用PCG算法不精確求解方程(7)，獲得模型修正量Δm.

(3) 確定線性搜索步長(zhǎng)α：

(4) 更新模型：

mi+1=mi+αΔm.

(5) 判斷是否繼續(xù)迭代：

計(jì)算目標(biāo)函數(shù)并估計(jì)數(shù)據(jù)殘差，若滿足精度要求則反演結(jié)束并輸出結(jié)果；若不滿足，重復(fù)步驟(2)—(4)進(jìn)行迭代.對(duì)于節(jié)點(diǎn)數(shù)為250000左右的非結(jié)構(gòu)網(wǎng)格，約5次迭代可得到較好的結(jié)果，總迭代次數(shù)一般不超過(guò)10次.反演程序由IntelFortran編譯運(yùn)行，在PC上計(jì)算耗時(shí)30min左右(CPU為i5-3.2GHz，內(nèi)存為8G).

圖6 山峰地形下高阻體模型的視電阻率平面圖(AB=100 m)Fig.6 Apparent resistivity of the high resistivity model under a hill (AB=100 m)

4 任意偶極-偶極數(shù)據(jù)的電阻率三維反演4.1 任意偶極-偶極觀測(cè)方式

實(shí)際上，由于非結(jié)構(gòu)網(wǎng)格單元的建模優(yōu)勢(shì)，可以針對(duì)不同地形特點(diǎn)，設(shè)計(jì)該地形條件下最合理的觀測(cè)方式，即對(duì)任意形式的偶極-偶極測(cè)量數(shù)據(jù)，無(wú)論觀測(cè)方式簡(jiǎn)單或復(fù)雜，均可用于本文任意地形條件下的三維電阻率反演.

4.2 平坦地形下高阻體模型的反演

模型設(shè)置如前文圖3中模型，即平坦地形下有一立方體高阻異常，背景電阻率為200 Ωm，異常體電阻率為2000 Ωm.反演結(jié)果如圖8所示，該結(jié)果表明地下存在高阻異常，與圖3所示的模型示意圖相當(dāng)一致，定位準(zhǔn)確，反演獲得很好的效果.

4.3 平坦地形下多異常體的反演

為了驗(yàn)證反演算法的正確性與分辨率，反演更為復(fù)雜一些的模型很有必要.圖9為平坦地形下并排兩個(gè)立方體異常模型，其中一個(gè)為低阻體，電阻率20 Ωm；另一個(gè)為高阻體，電阻率2000 Ωm，頂部埋深均為10 m，異常間隔20 m，背景電阻率為200 Ωm.

反演結(jié)果如圖10所示，清楚地反映了地下低阻與高阻異常體并存，橫向分辨率很好，并且沒(méi)有其他的多余結(jié)構(gòu)會(huì)影響解釋結(jié)果.與圖9所示的真實(shí)模型示意圖相比，定位準(zhǔn)確，吻合好，反演結(jié)果很好地分辨了地下真實(shí)的電性結(jié)構(gòu).

4.4 山峰地形下高阻體模型的反演

模型設(shè)置如前文圖5中模型，山峰地形下有一立方體高阻異常，背景電阻率為200 Ωm，異常體電阻率為2000 Ωm.倘若不考慮山峰地形的影響，直接使用平坦地形下的電阻率三維反演，則結(jié)果如圖11所示，顯示地下為一低阻異常體，與圖5的真實(shí)模型完全不一致.可見(jiàn)，地形造成的影響非常大，甚至獲得相反的異常結(jié)果.

吳小平(2005)的結(jié)果表明，無(wú)論是在數(shù)據(jù)空間還是在模型空間進(jìn)行地形校正，起伏地形條件下的電阻率三維反演均不能完全消除地形影響，必須在反演中直接引入地形信息，進(jìn)行帶地形的電阻率三維反演才有可能.圖12為該模型的帶地形電阻率三維反演結(jié)果，清晰地反映了山峰地形下的高阻異常體，與圖5的真實(shí)模型相當(dāng)一致，有效地消除了地形影響.

4.5 復(fù)雜起伏地形條件下的電阻率三維反演

實(shí)際野外勘察工作中，山峰、山谷起伏不定，地形條件十分復(fù)雜.隨著GPS/GNSS技術(shù)的日益成熟，復(fù)雜環(huán)境中觀測(cè)點(diǎn)的準(zhǔn)確定位變得簡(jiǎn)單易行，為靈活、高效的野外電阻率三維測(cè)量創(chuàng)造了有利條件.

圖7 隨機(jī)觀測(cè)裝置示意圖Fig.7 Schematic of the random acquisition

圖8 高阻體模型的反演結(jié)果Fig.8 Inversion of the high resistivity model

圖9 低阻和高阻共存的組合模型示意圖Fig.9 Schematic diagram of the model with low and high resistivity anomalous bodies

圖10 低阻和高阻異常體的組合模型反演結(jié)果Fig.10 Inversion of the combination model with a low resistivity and a high resistivity anomalous bodies

圖11 山峰地形下高阻體模型的反演結(jié)果，反演中沒(méi)有考慮地形Fig.11 Inversion of the high resistivity model under a hill，ignoring the influence of the topography

結(jié)合本文實(shí)現(xiàn)的任意偶極-偶極視電阻率觀測(cè)數(shù)據(jù)的帶地形三維反演，將在野外實(shí)際數(shù)據(jù)的三維反演解釋中發(fā)揮重要作用.

圖13是模擬野外實(shí)際情況構(gòu)建的一個(gè)相對(duì)復(fù)雜的起伏地形及其非結(jié)構(gòu)網(wǎng)格剖分.我們分別反演該復(fù)雜地形下的低阻異常體和高阻異常體三維模型，低阻異常體電阻率為20 Ωm，高阻異常體電阻率為2000 Ωm，背景電阻率均為200 Ωm，其他模型參數(shù)見(jiàn)圖3.當(dāng)起伏地形下分別存在低阻異常體和高阻異常體時(shí)，帶地形電阻率三維反演結(jié)果見(jiàn)圖14、圖15，均清晰地反映了地下的真實(shí)電阻率結(jié)構(gòu)，獲得很好的反演效果，表明本文基于非結(jié)構(gòu)網(wǎng)格的電阻率三維反演方法在野外復(fù)雜地形條件下亦可獲得可靠結(jié)果.

4.6 數(shù)據(jù)含10%高斯噪聲的三維反演

模型設(shè)置與4.2節(jié)相同，即平坦地形下有一立方體高阻異常，背景電阻率為200 Ωm，異常體電阻率為2000 Ωm.在生成合成數(shù)據(jù)時(shí)，將默認(rèn)引入的5%高斯噪聲增加為10%高斯噪聲.反演結(jié)果如圖16所示，同樣取得較好的反演結(jié)果，電阻率結(jié)構(gòu)清晰并且無(wú)多余的虛假信息.與圖8的反演結(jié)果相比，只是在高阻異常中心處的電阻率值略有減小，表明本文的電阻率三維反演結(jié)果可靠，反演方法對(duì)觀測(cè)數(shù)據(jù)的噪聲依賴性不大.

圖12 山峰地形下高阻體模型的帶地形反演結(jié)果Fig.12 Inversion of the high resistivity model under a hill，considering the influence of the topography

圖14 復(fù)雜起伏地形下存在低阻異常體的帶地形反演結(jié)果Fig.14 Inversion of the low resistivity model under a complex topography，considering the influence of the topography

圖15 復(fù)雜起伏地形下存在高阻異常體的帶地形反演結(jié)果Fig.15 Inversion of the high resistivity model under a complex topography，considering the influence of the topography

圖16 引入10%測(cè)量誤差時(shí)高阻體模型的反演結(jié)果Fig.16 Inversion of the high resistivity model with a measurement error of 10%

圖13 復(fù)雜起伏地形的非結(jié)構(gòu)網(wǎng)格剖分示意圖Fig.13 Unstructured grids of a complex topography

5 結(jié)論

隨著GPS/GNSS技術(shù)的日益成熟，野外區(qū)域地球物理調(diào)查中觀測(cè)點(diǎn)的準(zhǔn)確定位變得簡(jiǎn)單易行，為靈活、高效的野外電阻率三維測(cè)量創(chuàng)造了有利條件.結(jié)合電阻率三維非結(jié)構(gòu)有限單元數(shù)值模擬，本文實(shí)現(xiàn)了任意偶極-偶極視電阻率觀測(cè)數(shù)據(jù)的帶地形三維反演，理論模型反演以及模擬野外實(shí)際地形的復(fù)雜模型反演均取得很好的效果，模擬實(shí)際數(shù)據(jù)帶高斯噪聲的三維反演也取得可靠的結(jié)果，將對(duì)野外實(shí)際數(shù)據(jù)的三維反演解釋向前推進(jìn)關(guān)鍵一步.

Armijo L. 1966. Minimization of functions having Lipschitz continuous first partial derivatives.Pac.J.Math., 16(1): 1-3. Bentley L R, Gharibi M. 2004. Two- and three-dimensional electrical resistivity imaging at aheterogeneous remediation site.Geophysics, 69(3): 674-680.

Chambers J E, Ogilvy R D, Kuras O, et al. 2002. 3D electrical imaging of known targets at a controlled environmental test site.Environ.Geol.,41(6): 690-704.

Chambers J E, Wilkinson P B, Kuras O, et al. 2011. Three-dimensional geophysical anatomy of an active landslide in Lias Groupmudrocks, Cleveland Basin, UK.Geomorphology, 125(4): 472-484.

Chen J Y, Dang Y M, Cheng P F. 2007. Development and progress in GNSS.JournalofGeodesyandGeodynamics(in Chinese), 27(5): 1-4.

Chen Q L, Hu J H. 1998. The use of portable GPS in electrical prospecting.GeophysicalProspectingforPetroleum(in Chinese), 37(4): 128-133. Dabas M, Tabbagh A, Tabbagh J. 1994. 3-D inversion in subsurface electrical surveying-I. Theory.Geophys.J.Int.,119(3): 975-990.

Dai Q W, Chai X C, Chen D P. 2012. 3D DC resistivity inversion based on damped Gauss-Newton method.ChineseJournalofEngineeringGeophysics(in Chinese), 9(4): 375-379.

Dembo R S, Eisenstat S C, Steihaug T. 1982. Inexact newton methods.SIAMJ.Numer.Anal., 19(2): 400-408.

Di Q Y, Wang M Y, Yan S M, et al. 1997. The application of the high density resistivity method for the seawave-proof dam in Zhuhai-harbour.ProgressinGeophysics(in Chinese), 12(2): 79-88.

Fu L K. 1983. Electrical Prospecting Tutorial (in Chinese). Beijing:Geological Publishing House. Günther T, Rücker C, Spitzer K. 2006. Three-dimensional modelling and inversion of DC resistivity data incorporating topography—II. Inversion.Geophys.J.Int.,166(2): 506-517. Haber E, Ascher U M, Oldenburg D. 2000. On optimization techniques for solving nonlinear inverse problems.InverseProblems, 16(5): 1263-1280. Haber E. 2005. Quasi-newton methods for large-scale electromagnetic inverse problems.InverseProblems, 21(1): 305-323.

He J S. 1997. Development and prospect of electrical prospecting method.ChineseJ.Geophys.(ActaGeophysicaSinica)(in Chinese), 40(S1): 308-316. Holcombe H T, Jiracek G R. 1984. Three-dimensional terrain corrections in resistivity surveys.Geophysics, 49(4): 439-452.

Huang J G, Ruan B Y, Bao G S. 2004. Resistivity inversion on 3-D section based on FEM.JournalofCentralSouthUniversity(ScienceandTechnology)(in Chinese), 35(2): 295-299.

Kneisel C, Emmert A, K?stl J. 2014. Application of 3D electrical resistivity imaging for mapping frozen ground conditions exemplified by three case studies.Geomorphology, 210: 71-82.

Large D B. 1971. Electric potential near a spherical body in a conducting half-space.Geophysics, 36(4): 763-767.

Li C W, Xiong B, Wang Y X, et al. 2012. Inversion of IP logging based on improved Krylov subspace methods.ChineseJ.Geophys. (in Chinese), 55(11): 3862-3869, doi: 10.6038/j.issn.0001-5733.2012.11.034.

Li G F, Yan X L. 2010. Positioning accuracy test of portable GPS receiver and its application in geological exploration of mineral.GeotechnicalEngineeringWorld(in Chinese), 1(4): 380-384.

Li Y G, Oldenburg D W. 1994. Inversion of 3-D DC resistivity data using an approximate inverse mapping.Geophys.J.Int.,116(3): 527-537.

Li Y G, Spitzer K. 2002. Three-dimensional DC resistivity forward modelling using finite elements in comparison with finite-difference solutions.Geophys.J.Int.,151(3): 924-934.

Lin J, Shi L. 1996. Prospect for application of GPS survey in geophysical exploring.EquipmentforEarthScience(in Chinese), 4: 1-5.Liu B, Li S C, Li S C, et al. 2012a. 3D electrical resistivity inversion with least-squares method based on inequality constraint and its computation efficiency optimization.ChineseJ.Geophys. (in Chinese), 55(1): 260-268, doi: 10.6038/j.issn.0001-5733.2012.01.025.

Liu B, Nie L C, Li S C, et al. 2012b. 3D electrical resistivity inversion tomography with spatial structural constraint.ChineseJournalofRockMechanicsandEngineering(in Chinese), 31(11): 2258-2268.

Liu H M, Zhao H. 2013. The use of portable GPS in electrical and magnetic prospecting for network layout.West-ChinaExplorationEngineering(in Chinese), 25(8): 86-88.Liu S J, Song Z C, Liu G. 2006. Error analysis of handheld GPS positioning in field geological survey.SurveyingandMappingofSichuan(in Chinese), 29(1): 11-14.

Loke M H, Barker R D. 1996. Practical techniques for 3D resistivity surveys and data inversion.GeophysicalProspecting, 44(3): 499-523.Lu J J, Wu X P, Spitzer K. 2010. Algebraic multigrid method for 3D DC resistivity modelling.ChineseJ.Geophys.,53(3): 700-707.

Oppliger G L. 1984. Three-dimensional terrain corrections for mise-à-la-masse and magnetometric resistivity surveys.Geophysics, 49(10): 1718-1729.

Pan K J, Tang J T. 2014. 2. 5-D and 3-D DC resistivity modelling using an extrapolation cascadic multigrid method.Geophys.J.Int., 197(3): 1459-1470.

Papadopoulos N G, Tsourlos P, Tsokas G N, et al. 2006. Two-dimensional and three-dimensional resistivity imaging in archaeological site investigation.Archaeol.Prospect.,13(3): 163-181.

Papadopoulos N G, Tsourlos P, Tsokas G N, et al. 2007. Efficient ERT measuring and inversion strategies for 3D imaging of buried antiquities.NearSurf.Geophys.,5(6): 349-361.

Park S K, Van Gregory P. 1991. Inversion of pole-pole data for 3-D resistivity structure beneath arrays of electrodes.Geophysics, 56(7): 951-960.

Park S. 1998. Fluid migration in the vadose zone from 3-D inversion of resistivity monitoring data.Geophysics, 63(1): 41-51.

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

Ren Z Y, Tang J T. 2010. 3D direct current resistivity modeling with unstructured mesh by adaptive finite-element method.Geophysics, 75(1): H7-H17.

Ren Z Y, Tang J T. 2014. A goal-oriented adaptive finite-element approach for multi-electrode resistivity system.Geophys.J.Int.,199(1): 136-145.

R?dder T, Kneisel C. 2012. Permafrost mapping using quasi-3D resistivity imaging,Murtèl,Swiss Alps.NearSurfaceGeophysics, 10(2): 117-127.Rücker C, Günther T, Spitzer K. 2006. Three-dimensional modelling and inversion of dc resistivity data incorporating topography—I. Modelling.Geophys.J.Int.,166(2): 495-505.

Sasaki Y. 1994. 3-D resistivity inversion using the finite-element method.Geophysics, 59(12): 1839-1848. Spitzer K. 1995. A 3-D finite-difference algorithm for DC resistivity modelling using conjugate gradient methods.Geophys.J.Int.,123(3): 903-914.Tikhonov A N, Arsenin V A. 1977. Solutions of Ill-posed Problems. Washington:W. H. Winston and Sons.

Tsourlos P I, Szymanski J E, Tsokas G N. 1999. The effect of terrain topography on commonly used resistivity arrays.Geophysics, 64(5): 1357-1363.

Wan X L, Xi D Y, Gao E G, et al. 2005. 3-D resistivity inversion by the least-squares QR factorization method under improved smoothness constraint condition.ChineseJ.Geophys. (in Chinese), 48(2): 439-444, doi: 10.3321/j.issn:0001-5733.2005.02.030.

Wang K X, Li F Y, Liu H L, et al. 2012. Handheld GPS positioning accuracy and error.GNSSWorldofChina(in Chinese), 36(6): 83-86.

Wang W, Wu X P, Spitzer K. 2013. Three-dimensional DC anisotropic resistivity modelling using finite elements on unstructured grids.Geophys.J.Int.,193(2): 734-746. Weller A, Frangos W, Seichter M. 1999. Three-dimensional inversion of induced polarization data from simulated waste.J.Appl.Geophys., 41(1): 31-47.

Wu X P, Xu G M. 1999. Derivation and analysis of partial derivative matrix in resistivity 3-D inversion.OilGeophysicalProspecting(in Chinese), 34(4): 363-372.

Wu X P, Xu G M. 2000. Study on 3-D resistivity inversion using conjugate gradient method.ChineseJ.Geophys. (in Chinese), 43(3): 420-427, doi: 10.3321/j.issn:0001-5733.2000.03.016.

Wu X P, Wang T T. 2002. Progress of study on 3D resistivity inversion methods.ProgressinGeophysics(in Chinese), 17(1): 156-162, doi: 10.3969/j.issn.1004-2903.2002.01.023.

Wu X P. 2003. A 3-D finite-element algorithm for DC resistivity modelling using the shifted incomplete Cholesky conjugate gradient method.Geophys.J.Int.,154(3): 947-956.

Wu X P, Xiao Y F, Qi C, et al. 2003. Computations of secondary potential for 3D DC resistivity modelling using an incomplete Choleski conjugate-gradient method.GeophysicalProspecting, 51(6): 567-577.

Wu X P. 2005. 3-D resistivity inversion under the condition of uneven terrain.ChineseJ.Geophys. (in Chinese), 48(4): 932-936, doi: 10.3321/j.issn:0001-5733.2005.04.028.

Yang Y X. 2010. Progress, contribution and challenges of compass/Beidou satellite navigation system.ActaGeodaeticaetCartographicaSinica(in Chinese), 39(1): 1-6. Yi M J, Kim J H, Song Y, et al. 2001. Three-dimensional imaging of subsurface structures using resistivity data.GeophysicalProspecting, 49(4): 483-497.

Zhang J, Mackie R L, Madden T R. 1995. 3-Dresistivity forward modeling and inversion using conjugate gradients.Geophysics, 60(5): 1313-1325.

Zhdanov M S, Keller G V. 1994. The Geo-Electrical Methods in Geophysical Exploration (Vol. 31). Netherlands: Elsevier Science Limited. Zhou B, Greenhalgh S A. 2001. Finite element three-dimensional direct current resistivity modelling: accuracy and efficiency considerations.Geophys.J.Int.,145(3): 679-688.

附中文參考文獻(xiàn)

陳俊勇, 黨亞明, 程鵬飛. 2007. 全球?qū)Ш叫l(wèi)星系統(tǒng)的進(jìn)展. 大地測(cè)量與地球動(dòng)力學(xué), 27(5): 1-4.

陳清禮, 胡家華. 1998. 手持式GPS在電法勘探中的應(yīng)用. 石油物探, 37(4): 128-133.

戴前偉, 柴新朝, 陳德鵬. 2012. 基于阻尼型高斯牛頓法的三維直流電阻率反演. 工程地球物理學(xué)報(bào), 9(4): 375-379.

底青云, 王妙月, 嚴(yán)壽民等. 1997. 高密度電阻率法在珠海某防波堤工程中的應(yīng)用. 地球物理學(xué)進(jìn)展, 12(2): 79-88.

傅良魁. 1983. 電法勘探教程. 北京: 地質(zhì)出版社.

何繼善. 1997. 電法勘探的發(fā)展和展望. 地球物理學(xué)報(bào), 40(增刊): 308-316.

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

李長(zhǎng)偉, 熊彬, 王有學(xué)等. 2012. 基于改進(jìn)Krylov子空間算法的井中激電反演. 地球物理學(xué)報(bào), 55(11): 3862-3869, doi: 10.6038/j.issn.0001-5733.2012.11.034.

李國(guó)防, 閆新亮. 2010. 手持GPS定位精度測(cè)試及其在礦產(chǎn)勘查中的應(yīng)用. 礦產(chǎn)勘查, 1(4): 380-384.

林君, 石磊. 1996. GPS在地球物理勘探中的應(yīng)用展望. 地學(xué)儀器, 4: 1-5.

劉斌, 李術(shù)才, 李樹(shù)忱等. 2012a. 基于不等式約束的最小二乘法三維電阻率反演及其算法優(yōu)化. 地球物理學(xué)報(bào), 55(1): 260-268, doi: 10.6038/j.issn.0001-5733.2012.01.025.

劉斌, 聶利超, 李術(shù)才等. 2012b. 三維電阻率空間結(jié)構(gòu)約束反演成像方法. 巖石力學(xué)與工程學(xué)報(bào), 31(11): 2258-2268.

劉會(huì)敏, 趙灝. 2013. 手持GPS在電、磁法勘探測(cè)網(wǎng)布設(shè)工作中的應(yīng)用. 西部探礦工程, 25(8): 86-88.

劉少杰, 宋在超, 劉剛. 2006. 野外地質(zhì)測(cè)量中手持GPS定位的誤差分析. 四川測(cè)繪, 29(1): 11-14.

宛新林, 席道瑛, 高爾根等. 2005. 用改進(jìn)的光滑約束最小二乘正交分解法實(shí)現(xiàn)電阻率三維反演. 地球物理學(xué)報(bào), 48(2): 439-444, doi: 10.3321/j.issn:0001-5733.2005.02.030.

王克曉, 李鳳友, 劉煥玲等. 2012. 手持GPS定位精度與誤差的研究. 全球定位系統(tǒng), 36(6): 83-86.

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

吳小平, 徐果明. 2000. 利用共軛梯度法的電阻率三維反演研究. 地球物理學(xué)報(bào), 43(3): 420-427, doi: 10.3321/j.issn:0001-5733.2000.03.016.

吳小平, 汪彤彤. 2002. 電阻率三維反演方法研究進(jìn)展. 地球物理學(xué)進(jìn)展, 17(1): 156-162, doi: 10.3969/j.issn.1004-2903.2002.01.023.

吳小平. 2005. 非平坦地形條件下電阻率三維反演. 地球物理學(xué)報(bào), 48(4): 932-936, doi: 10.3321/j.issn:0001-5733.2005.04.028.

楊元喜. 2010. 北斗衛(wèi)星導(dǎo)航系統(tǒng)的進(jìn)展、貢獻(xiàn)與挑戰(zhàn). 測(cè)繪學(xué)報(bào), 39(1): 1-6.

(本文編輯 何燕)

3D resistivity inversion incorporating topography based on unstructured meshes

WU Xiao-Ping1，2， LIU Yang1，2， WANG Wei3

1LaboratoryofSeismologyandPhysicsofEarth’sInterior,SchoolofEarthandSpaceSciences,UniversityofScienceandTechnologyofChina,Hefei230026,China2NationalGeophysicalObservatoryatMengcheng,Mengcheng233500,China3ResearchInstitute,CNOOCLtd-Shenzhen,Guangzhou510240,China

Surface topographies have a great influence for the direct current (DC) resistivity methods, which cannot be avoided in actual mineral explorations. 3D DC resistivity forward modeling is available in recent years, especially for arbitrary topography and complicated subsurface structures using unstructured grids. However, surface topography is still a challenge for 3D interpretation in realistic applications, which may cause significant error in the 3D resistivity inversion without topography. Additionally it is a hard work to lay measurement points on regular observation network in complex terrains and the corresponding data cannot be simulated on ordinary structured grid. Therefore, 3D resistivity inversion incorporating topography based on unstructured meshes is necessary.We use unstructured finite element method for 3D resistivity forward modeling in order to simulate arbitrary topography and complicated subsurface structures. Our modeling result for a sphere model shows high accuracy in comparison with analytical solution. On the basis, we implement an inexact Gauss-Newton inversion for dipole-dipole configuration on arbitrary surface topography. With the development of GPS/GNSS technique, it is not necessary to lay measurement points on regular observation network exactly in the field survey. The inversion method developed in this paper can inverse the resistivity data from arbitrary dipole-dipole measurements, which is more convenient for 3D interpretation in realistic applications.A random acquisition system for arbitrary dipole-dipole measurements is designed, including 16 dipoles as transmitted electrodes and 100 random diploes as receiver electrodes, i.e. 1600 random dipole-dipole apparent resistivities. Firstly, flat terrain models are used to verify our 3D resistivity inversion for the random dipole-dipole apparent resistivities data, obtaining the inverted model in good agreement with subsurface geoelectrical structures. Then a high resistivity model under a mountain ridge is simulated to show the significant influence from surface topography. The 3D resistivity inversion obtains a low resistivity structure if the topography is ignored, showing a wrong subsurface structure. Our 3D resistivity inversion incorporating topography based on unstructured meshes, in which the topography is directly incorporated into the inversion algorithm, obtains the true high resistivity structure under a mountain ridge. The 3D inversions for models with complicated topography also turn out to be very successful. All dipole-dipole apparent resistivities data for synthetic examples above are generated with 5% Gaussian errors. The 3D resistivity inversion for synthetic data with 10% Gaussian errors is presented finally. Good result shows the 3D resistivity inversion algorithm in this study is very robust.It becomes simple and practicable for the location of measurement points in field geophysical survey using modern GPS/GNSS technique, providing favorable conditions for flexible and efficient 3D resistivity field measurements. In combination to 3D resistivity modeling using unstructured finite element method, we implement an inexact Gauss-Newton 3D resistivity inversion for the random dipole-dipole measurements on arbitrary surface topographies. Synthetic examples show our 3D inversion routines obtain good results for theoretical model and simulated realistic model with complicated topography. The 3D resistivity inversion for synthetic data with 10% Gaussian errors converges stably and the result is also reliable. The 3D resistivity inversion incorporating topography based on unstructured meshes in this paper promotes a key step towards the 3D field interpretation in realistic applications.

Resistivity； Three-dimensional inversion； Unstructured mesh； Topography

國(guó)家自然科學(xué)基金(41374076，41130420)、國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)(2014AA06A610、2012AA09A201)、天然氣水合物資源勘查與試采工程國(guó)家專項(xiàng)(GZH201400305)和國(guó)家重大科學(xué)儀器設(shè)備開(kāi)發(fā)專項(xiàng)項(xiàng)目任務(wù)(2011YQ05006008)聯(lián)合資助.

吳小平，男，1967年生，教授，博士生導(dǎo)師，研究方向?yàn)殡姶艤y(cè)深和地球內(nèi)部物理. E-mail:wxp@ustc.edu.cn

10.6038/cjg20150808.

10.6038/cjg20150808

P631

2014-12-22，2015-06-30收修定稿

吳小平， 劉洋， 王威. 2015. 基于非結(jié)構(gòu)網(wǎng)格的電阻率三維帶地形反演. 地球物理學(xué)報(bào)，58(8):2706-2717,

Wu X P, Liu Y, Wang W. 2015. 3D resistivity inversion incorporating topography based on unstructured meshes.ChineseJ.Geophys. (in Chinese),58(8):2706-2717,doi:10.6038/cjg20150808.