藍(lán)澤鸞，張志勇，鄧居智，周 峰，李 曼

（東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南昌 330013）

基于TV井地2．5D直流電阻率正則化反演

藍(lán)澤鸞，張志勇，鄧居智，周 峰，李 曼

（東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南昌 330013）

針對(duì)目前地面直流電阻率法所面臨的勘探深度有限的問題，開展了基于TV井地2．5D直流電阻率正則化反演研究。首先從井地直流電滿足的邊值問題出發(fā)，采用三角單元二次插值有限元法推導(dǎo)了2．5D井地電阻率正演公式。為了提高計(jì)算效率，對(duì)邊界采用近似處理、以及利用圖論理論的矩陣重排與填入元分析方法，實(shí)現(xiàn)了稀疏矩陣的直接解法，大大提高了計(jì)算效率。為了能夠提高反演結(jié)果的穩(wěn)定以及對(duì)異常邊界分辨能力，采用L曲線法求取正則化因子、TV的穩(wěn)定因子以及共軛梯度方法進(jìn)行反演目標(biāo)函數(shù)的求解，研究表明，反演結(jié)果高效、準(zhǔn)確。最后對(duì)不同井地觀測(cè)方式的反演結(jié)果進(jìn)行對(duì)比分析，得出井地直流電阻率異常特征和分布規(guī)律。

井地直流電法；共軛梯度法；有限單元法；TV穩(wěn)定因子；L曲線

0 引言

直流電阻率法已被廣泛用于礦產(chǎn)資源勘探、環(huán)境與工程物探中，已成為物探常用方法。傳統(tǒng)電阻率測(cè)量主要采用的是固定裝置（pole-pole、dipole-dipole、Wenner以及Schlumberger）進(jìn)行地表觀測(cè)，從而得到電阻率曲線或擬斷面圖，其大致能顯示地下電阻率在垂直方向或水平方向的變化。但是在某些特殊情況下，由于測(cè)區(qū)的限制或電阻率在水平方向差異很大，導(dǎo)致傳統(tǒng)的電法勘探受限而無法勘探到深部的目標(biāo)體，因此開展井地電阻率成像技術(shù)或反演已變成目前的研究重點(diǎn)［1－6］。

井地直流電阻率反演研究主要是針對(duì)井－地－井、井－地、井－井等觀測(cè)方式，Zhou［7］等利用解析靈敏度函數(shù)實(shí)現(xiàn)了快速二維／三維井間電阻率成像；Wikinson［8］等重點(diǎn)討論了礦井直流電阻率成像技術(shù)問題；國(guó)內(nèi)有關(guān)井地直流電阻率反演研究同樣取得了相當(dāng)大的成果，王志剛等［9－11］采用了Born近似實(shí)現(xiàn)了井地電法數(shù)據(jù)的三維反演、井地電法的準(zhǔn)解析近似三維反演研究，以及利用積分方程法實(shí)現(xiàn)了井地電法的三維模擬的研究；呂玉增等［12］實(shí)現(xiàn)了直流電井間三維直接成像；沈平等［13］研究了井間視電阻率的幾何成像方法；岳建華等［14］開展了井－地三維電阻率成像技術(shù)的研究；徐凱軍等［15］實(shí)現(xiàn)了基于共軛梯度法的垂直有限線源的三維電阻率反演；柯敢攀等［16］開展了井地電法的三維正反演研究；安然等［17］進(jìn)行了井地電阻率三維反演研究。對(duì)目前來講，地下情況復(fù)雜，現(xiàn)有的電阻率成像技術(shù)在一定程度難以滿足現(xiàn)有的需求，這里采用了全變差正則化方法，該方法具有一定的邊界識(shí)別能力，是Rudin等［18］在1992提出，主要應(yīng)用于圖像邊界的識(shí)別問題；后來Wang等［19］利用混合正則化方法研究了地震波的反演；Acar等［20］對(duì)該方法進(jìn)行了改進(jìn)，分析該方法的特性。因此借鑒這些理論，采用改進(jìn)的全變差穩(wěn)定因子應(yīng)用到井地直流電阻率反演中，為了提高反演結(jié)果的穩(wěn)定性，運(yùn)用基于L曲線方法求取正則化因子方法，數(shù)值計(jì)算表明，反演結(jié)果準(zhǔn)確、可靠。

1 直流電阻率正演原理

2．5D直流電阻率所滿足的數(shù)學(xué)模型歸結(jié)為點(diǎn)源二維問題，其經(jīng)過傅里葉變換后，將變成波數(shù)域二維問題［21］：

式（1）中：區(qū)域Ω為二維研究區(qū)域；Γs、?！奘嵌S區(qū)域的邊界，Γs為區(qū)域Ω的地表邊界，?！逓閰^(qū)域Ω的地下邊界；σ為電導(dǎo)率；r為點(diǎn)電源到邊界的距離；n為無窮遠(yuǎn)邊界的外法向方向；k為波數(shù)；K0、K1分別為零階、一階第二類貝塞爾函數(shù)。

將公式（1）轉(zhuǎn)化為變分問題：

對(duì)式（2）求變分并令其為零，得到線性方程組（3）。

通過求解線性方程組（3），得到波數(shù)域的電位U，然后通過傅氏反變換得到電位。

式中：r為點(diǎn)的位置；km是波數(shù)；gm是加權(quán)系數(shù)［21］（表1）。

表1 傅氏變換系數(shù)km和gmTab．1 Fourier transform coefficients，kmand gm

鑒于各種地表的地球物理觀測(cè)方法，觀測(cè)數(shù)據(jù)對(duì)淺部分辨率高、向深部逐漸降低，基于這樣的認(rèn)識(shí)，設(shè)計(jì)了基于二叉樹的收縮網(wǎng)格（圖1）剖分算法。該算法采用三角單元剖分，具有較好地?cái)M合地表地形與地下地質(zhì)體的能力。為了進(jìn)一步提高正演計(jì)算的效率和精度，作者采用了基于圖論理論的矩陣重排與填入元分析算法，實(shí)現(xiàn)了基于TV井地2．5D直流電阻率正則化反演算法。該算法采用了改進(jìn)的Cholesky［22］分解算法，其只需要對(duì)總體系數(shù)矩陣進(jìn)行一次分解，然后對(duì)于不同的右端向量進(jìn)行回代即可。正演計(jì)算邊界對(duì)測(cè)點(diǎn)的影響大，通常測(cè)點(diǎn)距離邊界越近，邊界的影響就越大，為了減少邊界對(duì)測(cè)量的影響，采用了邊界近似，主要是將區(qū)域分成目標(biāo)區(qū)和邊界區(qū)，目標(biāo)區(qū)為異常體存在的區(qū)域，其網(wǎng)格步長(zhǎng)采用等距剖分；邊界區(qū)的網(wǎng)格步長(zhǎng)采用指數(shù)遞增的形式來模擬無窮邊界。

圖1 三角單元剖分示意圖Fig．1 Sketch map of triangular element subdivision

2 基于TV井地直流電阻率正則化反演

反演是資料解釋的重要環(huán)節(jié)，為了能夠很好地對(duì)異常體邊界進(jìn)行識(shí)別，研究了基于TV井地直流電阻率正則化反演技術(shù)。然而目標(biāo)函數(shù)的最優(yōu)化求解與正則化因子的求取是正則化反演的關(guān)鍵，反演采用L曲線法，實(shí)現(xiàn)正則化因子的求?。焕肅G方法求解反演目標(biāo)函數(shù)，提高反演過程的穩(wěn)定性與計(jì)算效率。

2．1 基于TV正則化反演

Tikhonov［23］正則化為：

其中：P（α，m）為目標(biāo)函數(shù)，α是正則化因子；φd（m）也稱數(shù)據(jù)目標(biāo)函數(shù)，通常采用實(shí)測(cè)數(shù)據(jù)與預(yù)測(cè)數(shù)據(jù)的距離表示：

其中：mapr為先驗(yàn)?zāi)Ｐ停唬轜m為模型權(quán)重矩陣。如果以式（7）為標(biāo)準(zhǔn)形式，則通過引入矩陣＾We可以得到不同穩(wěn)定因子的統(tǒng)一表示形式［33］。

為了能夠有效地識(shí)別邊界，Rudin提出了Total varivation（TV）因子

其中β為一個(gè)小的數(shù)。分析知，式（9）作為穩(wěn)定因子將突出模型梯度變化，能夠突出邊界識(shí)別。

為了利用式（9）統(tǒng)一表示穩(wěn)定因子，由式（9）構(gòu)造模型權(quán)重函數(shù)為式（10）。

其中：m（r）為模型參數(shù)分布函數(shù)；▽m（r）為模型參數(shù)梯度；e為與計(jì)算機(jī)數(shù)值精度有關(guān)的一個(gè)很小正數(shù)。然后將公式（6）、式（7）、式（10）帶入到目標(biāo)函數(shù)式（5）中，即可得到基于TV穩(wěn)定因了的正則化反演目標(biāo)函數(shù)。

2．2 共軛梯度法（CG）

采用共軛梯度法［15，22］求解公式（12）最優(yōu)化問題，其基本思想跟最速下降法一樣：

通過線性搜索使目標(biāo)函數(shù)最小化的方式來求解下降步長(zhǎng)：

然后通過對(duì)目標(biāo)函數(shù)進(jìn)行一系列的化簡(jiǎn)，得到下降步長(zhǎng)式（17）。

2．3 L曲線求取正則化因子

對(duì)于公式（5），正則化因子α是數(shù)據(jù)目標(biāo)函數(shù)與模型目標(biāo)的權(quán)重系數(shù)，它決定了反演結(jié)果的好壞。α的過大與過小都會(huì)影響反演效果，如何正確地選擇正則化因子是反演的關(guān)鍵。

作者采用由Hansen等［24］提出的L曲線法求取正則化因子，該方法是一種基于數(shù)據(jù)誤差水平未知的啟發(fā)式選取正則化因子的方法，它是通過對(duì)φd（m）與φm（m）在雙對(duì)數(shù)下進(jìn)行刻畫，其曲線形狀為“L”，其正則化因子為曲線曲率最大值，一般是“隅角”對(duì)應(yīng)的位置，如圖2所示。

圖2 L曲線Fig．2 L－curve

曲線曲率的表達(dá)式為：

其中：ρ＝φd（m）；η＝φm（m）；η′為η對(duì)α的一階導(dǎo)數(shù)。

3 算例分析

3．1 均勻半空間模型

為了驗(yàn)證算法的正確性，設(shè)計(jì)均勻半空間電阻率為100Ω·m，計(jì)算了理論電位曲線并與數(shù)值模擬結(jié)果進(jìn)行對(duì)比，計(jì)算結(jié)果如圖3所示，相對(duì)誤差在2%以內(nèi)。

圖3 數(shù)值解與解析解電位對(duì)比曲線圖Fig．3 Potential contrast curve of numerical solution and analytical solution

3．2 模型I

設(shè)置如圖4所示的地電模型，即均勻半空間存在一個(gè)傾斜低阻異常體，異常體的電阻率大小為50Ω·m，背景電阻率的大小為100Ω·m。設(shè)計(jì)了井－地－井的觀測(cè)方式，井所在的位置如圖5虛線所示，并且反演數(shù)據(jù)由正演數(shù)據(jù)加入1%、10%、20%以及30%隨機(jī)噪聲得到，對(duì)比不同隨機(jī)噪聲的反演結(jié)果，如圖5所示。

圖5為不同噪聲水平的反演結(jié)果，隨著噪聲逐漸增強(qiáng)，反演結(jié)果不精細(xì)，反演得到的電阻率斷面逐漸變得模糊，很難分辨出異常體的傾斜方向以及大致范圍。由圖5（a）可以看出，低阻異常體清晰可見，其傾斜的方向也能夠反應(yīng)出異常體傾向，異常范圍集中，且背景電阻率與真實(shí)電阻率相近。隨著噪聲水平的增加，異常范圍相應(yīng)地反映出異常體所在的位置，但是異常體的傾向難以判別，且反演斷面圖存在一些假異常，背景電阻率與真實(shí)電阻率值相差較大，井所在的位置反演結(jié)果隨著噪聲水平的增加，影響較大。

3．3 組合模型

在均勻半空間存在兩個(gè)異常體（圖6），M1為高阻異常體，電阻率值為200Ω·m；M2為低阻異常體，電阻率值為50Ω·m，背景電阻率值為100 Ω·m。設(shè)計(jì)兩種井地的觀測(cè)方式，分別為井－地－井與地－井－地，井所在的位置如圖7虛線所示，并且反演數(shù)據(jù)由正演數(shù)據(jù)加入1%、10%、20%以及30%隨機(jī)噪聲得到，通過對(duì)比兩種觀測(cè)方式不同隨機(jī)噪聲的反演結(jié)果，其反演結(jié)果如圖7所示。

圖7為不同噪聲水平的反演結(jié)果，其反演結(jié)果都能大致反映出異常所在的位置，但是隨著噪聲增加，反演結(jié)果不精細(xì)，且井所在的位置的反演結(jié)果誤差較大。從圖7可以看出，當(dāng)隨機(jī)噪聲為10%、20%和30%時(shí)，反演斷面圖上存在小區(qū)域的假異常，圍巖電阻率與真實(shí)電阻率誤差較大；隨機(jī)噪聲為1%時(shí)，背景電阻率與真實(shí)電阻率相近，異常體的反演相對(duì)集中，其值與真實(shí)值相近。

圖8為不同噪聲水平的反演結(jié)果，反演結(jié)果都能大致反應(yīng)出異常所在的位置，但是隨著噪聲增加，反演結(jié)果不精細(xì)，且井所在的位置的反演結(jié)果誤差較大。從反演斷面圖可以看出，當(dāng)隨機(jī)噪聲為10%、20%和30%時(shí)，反演斷面圖上存在小區(qū)域的假異常，圍巖電阻率與真實(shí)電阻率誤差較大；隨機(jī)噪聲為1%時(shí)，背景電阻率與真實(shí)電阻率相近，異常體的反演相對(duì)集中，其值與真實(shí)值相近。

對(duì)比圖7與圖8反演結(jié)果發(fā)現(xiàn)，井－地－井測(cè)量裝置相對(duì)單井反演效果更明顯，異常范圍更集中，能夠很好地反映出異常的邊界，易于圈定異常體的范圍。隨著噪聲水平的增加，井－地－井測(cè)量方式相對(duì)單井反演結(jié)果影響較小，易于圈定異常的范圍。

圖4 地電模型斷面圖Fig．4 Sketch of the model

4 結(jié)論

作者在前人的基礎(chǔ)上，開發(fā)了基于TV井地2．5D直流電阻率正則化反演算法，并進(jìn)行模型試算，試算表明，程序計(jì)算正確，算法穩(wěn)定高效。研究工作得到以下成果：

1）井－地－井測(cè)量方式相對(duì)單井反演效果更明顯，異常范圍更集中，能夠很好地反映出異常的邊界，易于圈定異常體的范圍。

2）基于L曲線求取正則化因子的方法以及TV穩(wěn)定因子的選擇，提高了反演結(jié)果的穩(wěn)定性。

3）近似的邊界處理以及基于圖論理論的矩陣重排，提高了正演計(jì)算效率及精度。

圖5 不同噪聲水平的井－地－井電阻率反演斷面圖Fig．5 Different noise level of the well－to－well resistivity inversion section

圖6 地電模型斷面圖Fig．6 Sketch of the model

圖7 不同噪聲水平的地－井－地反演斷面圖Fig．7 Different noise level of the Surface－Borehole－Surface resistivity inversion section

圖8 不同噪聲水平的井－地－井反演斷面圖Fig．8 Different noise level of the well－to－well resistivity inversion section

［1］SMTH，N．C，VOZOFF，K．Two-dimensional dc resistivity inversing for dipole－dipole data［J］．IEEETrans．Geosci．Remote Sensing，1984，GE－22 （02）：21－28．

［2］ZHOU B，GREENHALGH，S A．A synthetic study on cross-hole resistivity imaging with different electrode arrsys［J］．Geophys，1997，28（2）：1－5．

［3］ZHANG，J，MACKIE，R，MADDEN，T．3－D resistivity forward modeling and inversion using conjugate gra-dients［J］．Geophysics，1995，60（6）：1313－1325．

［4］OLDENBURG，D W，MCGILLICRAY，P R，ELLIS，R G．Generalised subspace methods for largescale inverse problem［J］．Geophys．1993，114（4）：12 －20．

［5］LI Y G，OLDENBURG，D W．Approximate inverse mapping in dc resistivity problems［J］．Geophys．1992，109（5）：343－362．

［6］MAURIELLOP，PATELLA，D．Resistivity anomaly imaging by probability tomography［J］．Geophys．1999，47（2）：411－429．

［7］BING ZHOU STEWART A，Greenhalgh．Rapid 2-D／3-D crosshole resistivity imaging using the analytic sensitivity function［J］．GEOPHYSICS，2002，67（3）：755－765

［8］WILKINSON P B，CHAMBERS J E，et al．Extreme sensitivity of cross-hole electrical resistivity tomography measurements to geometric errors［J］．Geophys．2008，173（3）：49－62．

［9］王志剛，何展翔，魏文博．Born近似快速三維反演井地電法數(shù)據(jù)［J］．球物理學(xué)進(jìn)展，2007，22（02）：508－513．

WANG ZH G，HE ZH X，WEI W B．Fast 3Dinversion borehole ground electrical method data based on born on approximation［J］．progress in geophsycis，2007，22（02）：508－513．（In Chinese）

［10］王志剛，何展翔，魏文博，等．井地電法的準(zhǔn)解析近似三維反演研究［J］．石油地球物理勘探，2007，42（2）：220－225．

WANG ZH G，HE ZH X，WEI W B，et al．3－D quasi－analytic approximate inversion of borehole－to－surface electric data［J］．OGP，2007，42（2）：220－225．（In Chinese）

［11］王志剛，何展翔，魏文博．積分方程法三維模擬井地電法并行算法研究［J］．物探化探計(jì)算技術(shù)，2007，29（5）：425－436．

WANG ZH G，HE ZH X，WEI W B．Parallel algorithm research of integral equation method to 3Dbore －h(huán)ole surface electromagnetic modeling［J］．Computing techniques for geophysical and geochemical exploration，2007，29（5）：425－436．（In Chinese）

［12］呂玉增，阮百堯，黃俊革．直流電井間三維直接成像［J］．物探化探計(jì)算技術(shù)，2003，25（01）：60－64．

LU Y Z，RUAN B Y，HUANG J G．The 3－D immediate crosshole tomography with direct current［J］．2003，25（01）：60－64．（In Chinese）

［13］沈平，強(qiáng)建科，李永軍，等．井間視電阻率的幾何成像方法［J］．中南大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，41（3）：1079 －1084．

SHEN P，QIANG J K，LI Y J，et al．Geometry image method of crosshole apparent resistivity［J］．Journal of Central South University（Sicence and Technology），2010，41（3）：1079－1084．（In Chinese）

［14］岳建華，劉志新．井－地三維電阻率成像技術(shù)［J］．地球物理學(xué)進(jìn)展，2005，20（2）：407－411．

YUE J H，LIU ZH X．Three dimension resistivity tomography of mine ground［J］．progress in geophusics，2005，20（2）：407－411．（In Chinese）

［15］徐凱軍，李桐林，張輝，等．基于共軛梯度法的垂直有限源三維電阻率反演［J］．煤田地質(zhì)與勘探，2006，34（3）：69－71．

XU K J，LI T L，ZHANG H，et al．3Dresistivity inversion of vertical finite line source using conjugate gradients［J］．Coal geology ＆exploration，2006，34（3）：69 －71．（In Chinese）

［16］柯敢攀，黃清華．井地電法的三維正反演研究［J］．北京大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，45（03）：264－272．

KE G P，HUANG Q H．3DForward and inversion problems of borehole－to－surface electrical method ［J］．Acta scientiarum Naturalium Universitatis Pekinensis2009，45（03）：264－272．（In Chinese）

［17］安然，李桐林，徐凱軍．井地三維電阻率反演研究［J］．地球物理學(xué)進(jìn)展，2007，22（1）：247－249．

AN R，LI T L，XU K J．Wellsurface 3Dresistivity inversion［J］．Progress in geophsycis，2007，22（1）：247 －249．（In Chinese）

［18］RUDIN L I，OSHER S，F(xiàn)ATEMI E．Nonliner total variation based noise removal algorithms［J］．Phsica D：Nonliner Phenomena，1992，60（1－4）：259－268．

［19］WANG Y F，CUI Y，YANG C C．Hybrid regularization methods for seismic reflectiveity inversion［J］．Int，J．Geomath，2011，2（1）：87－112．

［20］ACAR R，VOGEL C R．Analysis of total variation penalty methods［J］．Inverse Problems，1994，（10）：1217－1229．

［21］徐世浙．地球物理中的有限單元法［M］．北京：科學(xué)出版社，1994．

XU SH ZH．The finite element method in geophysics ［M］．Beijing：Science Press，1994．（In Chinese）

［22］吳小平，徐果明，李時(shí)燦，等．利用不完全Cholesky共軛梯度法求解點(diǎn)源三維地電場(chǎng)［J］．地球物理學(xué)報(bào)，1998，41（06）：848－855．

WU X P，XU G M，LI SH C，et al．The calculation of three－dimensional geoelectric field of point source by incomplete cholesky conjugate gradient method［J］．Chinese Journal Geophysics．1998，41（06）：848－855．（In Chinese）

［23］TIKHONOV A N，ARSENIN V Y．Solution of ill－Posed Problems［M］．Washington．DC：V H Winston and Sons，1977．

［24］HANSEN，C．Rank－deficient and discrete ill－posed problems［M］．Numerical aspects of liner inversion：Department of mathematical modeling，Technical University of Denmark Lyngby，1998．

2．5D inversion of borehole-surface DC resistivity based on the total variation function

LAN Ze-luan，ZHANG Zhi-yong，DENG Ju-zhi，ZHOU Feng，LI Man
（Key Laboratory of Radioactive Geology and Exploration Technology Fundamental Science for National Defense，East China Institute of Technology，Nanchang 330013，China）

The 2．5Dregularization inversion for borehole－to－ground DC resistivity has been implemented in this paper，which is based on total variations（TV）to facing the problem that the ground DC resistivity exploration depth was limited．To start from the well satisfied boundary value problem，using the triangular unit finite element quadratic interpolation derived to 2．5Dforward modeling formulation．While the approximate boundary treatment，as well as the use of graph theory matrix rearrangement and fill element analysis method to achieve a direct solution of sparse matrices，greatly improving the computational efficiency．In order to improve the stability of the inversion results and the ability to distinguish the exception of boundaries，using L－curve to achieve regularization factor，stability factor of the TV and the conjugate gradient method for solving the objective function，which have shown that efficiency and accuracy．Finally，the inversion results in different borehole－toground observation way were on the basis of contrastive analysis to reach the anomaly characteristics and distribution rule of borehole－to－ground dc resistivity．

borehole-surface electrical method；conjugate gradient method；finite element method；total variation；L-curve

P 631．3

A

10．3969／j．issn．1001-1749．2015．03．02

1001-1749（2015）03-0273-07

2014-10-19 改回日期：2014-12-28

國(guó)家自然科學(xué)基金（41304055，41304056）；東華理工大學(xué)研究生創(chuàng)新基金（DYCA13026）；江西省研究生創(chuàng)新基金（YC2013－S177）

藍(lán)澤鸞（1988－），女，碩士，主要從事電法數(shù)值模擬研究，E-mail：ZFECIT＠163．com。