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

藍澤鸞，張志勇，鄧居智，周 峰，李 曼

（東華理工大學放射性地質與勘探技術國防重點學科實驗室，南昌 330013）

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

藍澤鸞，張志勇，鄧居智，周 峰，李 曼

（東華理工大學放射性地質與勘探技術國防重點學科實驗室，南昌 330013）

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

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

0 引言

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

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

1 直流電阻率正演原理

2．5D直流電阻率所滿足的數學模型歸結為點源二維問題，其經過傅里葉變換后，將變成波數域二維問題［21］：

式（1）中：區域Ω為二維研究區域；Γs、Γ∞是二維區域的邊界，Γs為區域Ω的地表邊界，Γ∞為區域Ω的地下邊界；σ為電導率；r為點電源到邊界的距離；n為無窮遠邊界的外法向方向；k為波數；K0、K1分別為零階、一階第二類貝塞爾函數。

將公式（1）轉化為變分問題：

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

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

式中：r為點的位置；km是波數；gm是加權系數［21］（表1）。

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

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

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

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

反演是資料解釋的重要環節，為了能夠很好地對異常體邊界進行識別，研究了基于TV井地直流電阻率正則化反演技術。然而目標函數的最優化求解與正則化因子的求取是正則化反演的關鍵，反演采用L曲線法，實現正則化因子的求取；利用CG方法求解反演目標函數，提高反演過程的穩定性與計算效率。

2．1 基于TV正則化反演

Tikhonov［23］正則化為：

其中：P（α，m）為目標函數，α是正則化因子；φd（m）也稱數據目標函數，通常采用實測數據與預測數據的距離表示：

其中：mapr為先驗模型；＾Wm為模型權重矩陣。如果以式（7）為標準形式，則通過引入矩陣＾We可以得到不同穩定因子的統一表示形式［33］。

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

其中β為一個小的數。分析知，式（9）作為穩定因子將突出模型梯度變化，能夠突出邊界識別。

為了利用式（9）統一表示穩定因子，由式（9）構造模型權重函數為式（10）。

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

2．2 共軛梯度法（CG）

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

通過線性搜索使目標函數最小化的方式來求解下降步長：

然后通過對目標函數進行一系列的化簡，得到下降步長式（17）。

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

對于公式（5），正則化因子α是數據目標函數與模型目標的權重系數，它決定了反演結果的好壞。α的過大與過小都會影響反演效果，如何正確地選擇正則化因子是反演的關鍵。

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

圖2 L曲線Fig．2 L－curve

曲線曲率的表達式為：

其中：ρ＝φd（m）；η＝φm（m）；η′為η對α的一階導數。

3 算例分析

3．1 均勻半空間模型

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

圖3 數值解與解析解電位對比曲線圖Fig．3 Potential contrast curve of numerical solution and analytical solution

3．2 模型I

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

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

3．3 組合模型

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

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

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

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

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

4 結論

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

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

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

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

圖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近似快速三維反演井地電法數據［J］．球物理學進展，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］王志剛，何展翔，魏文博，等．井地電法的準解析近似三維反演研究［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］．物探化探計算技術，2007，29（5）：425－436．

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

［12］呂玉增，阮百堯，黃俊革．直流電井間三維直接成像［J］．物探化探計算技術，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］沈平，強建科，李永軍，等．井間視電阻率的幾何成像方法［J］．中南大學學報：自然科學版，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］岳建華，劉志新．井－地三維電阻率成像技術［J］．地球物理學進展，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］．煤田地質與勘探，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］．北京大學學報：自然科學版，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］．地球物理學進展，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，FATEMI 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］．北京：科學出版社，1994．

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

［22］吳小平，徐果明，李時燦，等．利用不完全Cholesky共軛梯度法求解點源三維地電場［J］．地球物理學報，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

國家自然科學基金（41304055，41304056）；東華理工大學研究生創新基金（DYCA13026）；江西省研究生創新基金（YC2013－S177）

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