地-井時域激電擬牛頓法反演問題研究

潘和平 孟慶鑫

(1.中國地質(zhì)大學(xué)地球物理與空間信息學(xué)院，湖北 武漢 430074；2.石家莊經(jīng)濟學(xué)院勘查技術(shù)與工程學(xué)院，河北 石家莊 050031)

引 言

井中激電地表-井中工作方式(地-井(Induced Polarization method,IP)IP)可實現(xiàn)多方位供電觀測，接收裝備置于井下接近礦體，獲得異常較大，對鉆孔地況要求低，工作效率高，觀測信息較之井中直流電阻率法(井中Direct Current Resistivity method,DC)更豐富，解釋判斷方法原理直觀簡便，是金屬礦普查或危機礦山勘探工作中重要的井中物探方法，取得了較好的應(yīng)用效果[1-3].

井中DC/IP反演研究主要針對地面及井-地測量方式，Li等[4]針對地面與井中觀測激電數(shù)據(jù)進行三維反演研究，Weller等[5]研究了激電數(shù)據(jù)的三維反演問題，Zhou等[6]完成了鉆孔電阻率快速成像，Wilkinson等[7]研究了礦井電阻率成像問題；國內(nèi)相關(guān)反演研究也取得了很大成果，阮百堯等[8]、吳小平等[9]分別研究了地面觀測方式下的激電與電阻率法的反演技術(shù)，呂玉增等[10]、沈平等[11]研究了井間電阻率成像問題，岳建華等[12]、安然等[13]進行了井-地電阻率法成像研究.對于地-井激電的反演，由于觀測數(shù)據(jù)量少，很難進行成像反演，主要采用切線法、正演擬合法等，來確定埋深、距井距離等參數(shù)以判斷礦體方位等[14].

本文應(yīng)用有限差分法和解析式法模擬三維地電模型中地-井激電方位測量的異常響應(yīng)，基于正演模擬，采用理論解析式建立反演目標(biāo)函數(shù)，利用最優(yōu)化擬牛頓法[15-16]反演目標(biāo)體方位等參數(shù).

1 地-井時域激電三維正演

1.1 解析法與數(shù)值模擬法

正演是反演的基礎(chǔ)，時域激電正演模擬主要由解析法和數(shù)值模擬法實現(xiàn).解析法是針對球體、橢球體等規(guī)則形體，通過數(shù)學(xué)解法推導(dǎo)出電場分布情況，結(jié)果精確嚴(yán)格，一般作為理論解和參考比對依據(jù)[1].數(shù)值模擬利用數(shù)值方法對描述地球物理現(xiàn)象的偏微分方程及相應(yīng)的邊界條件進行近似求解模擬，適用性強效率高.黃智輝等[17]應(yīng)用解析法對極化橢球體的地-井測量情況進行了闡述，文獻[18]和[19]分別采用有限差分法和有限元法對地-井IP進行正演模擬并闡述分析了其異常特征；聶在平等[20-21]應(yīng)用有限元法進行了井中電磁響應(yīng)數(shù)值模擬研究.

本文正演主要基于有限差分法，以等效電阻率原理為出發(fā)點，采用解析解與數(shù)值解相結(jié)合的方式對三維地電模型的地-井IP進行正演模擬，正演所涉及的稀疏矩陣方程計算選用穩(wěn)定雙共軛梯度法[22]，減少迭代次數(shù)提高了效率.離散差分方程的推導(dǎo)過程及截斷邊界條件處理方法見文獻[23]，本文僅列出穩(wěn)流電場在求解三維區(qū)域內(nèi)的七點差分格式

(1)

式中：

ao-1，p，q=(σo-1，p，q+σo，p，q)/[(Δxo+Δxo-1)Δxo-1]；

ao，p-1，q=(σo，p-1，q+σo，p，q)/[(Δyp+Δyp-1)Δyp-1]；

ao，p，q-1=(σo，p，q-1+σo，p，q)/[(Δzq+Δzq-1)Δzq-1]；

ao+1，p，q=(σo+1，p，q+σo，p，q)/[(Δxo+Δxo-1)Δxo]；

ao，p+1，q=(σo，p+1，q+σo，p，q)/[(Δyp+Δyp-1)Δyp]；

ao，p，q+1=(σo，p，q+1+σo，p，q)/[(Δzq+Δzq-1)Δzq]；

ao，p，q=-[ao-1，p，q+ao+1，p，q+ao，p-1，q+ao，p+1，q+

ao，p，q-1+ao，p，q+1]；

此外，本文還采用了穩(wěn)定點電流場下井旁球狀極化體的解析式作為有限差分模擬的參考依據(jù)和反演工作中的正演計算，其表達式可見文獻[1].

1.2 實際算例與地-井時域IP異常的一般特征規(guī)律

如圖1(a)所示，縮小比例模型取均勻三維網(wǎng)格剖分，由離散網(wǎng)格節(jié)點上的電性值代替整個連續(xù)地電模型，中心垂直虛線為垂直鉆孔，A0、A1、A2、A3、A4分別為場源布設(shè)于井口方位、主方位、反方位和兩個輔助方位，虛線立方體表示井旁異常體.

算例1地電模型如圖1(b)所示，無極化均質(zhì)大地電阻率100 Ω·m，尺寸100 m×100 m×100 m

(a)

(b)

(c)圖1 三維地電模型示意圖

(三維網(wǎng)格剖分?jǐn)?shù)101×101×101)；井旁直徑8 m極化球體，球心坐標(biāo)(x=9 m，y=0 m，z=-39 m)，電阻率10 Ω·m，極化率50%；場源位于A1方位坐標(biāo)(x=11 m，y=0 m，z=0 m)，電流強度20 A；井Zk平面坐標(biāo)(x=0 m，y=0 m).

由圖2正演模擬結(jié)果對比可知：有限差分解所得曲線形態(tài)、數(shù)據(jù)結(jié)果都與解析解極為近似，顯示出有限差分法在地-井IP模擬方面的適用性和準(zhǔn)確性；但模擬計算中仍存在誤差，這是因為有限差分只是利用有限的離散節(jié)點代替連續(xù)空間模型，誤差會隨測點接近異常體而逐步變大，但考慮到電位差的數(shù)值大小，這種誤差是允許的.

(a) 一次場電位差/V (b) 二次場電位差/V圖2 有限差分解與解析解對比圖

算例2地電模型如圖1(c)所示，鉆孔位于原點(0 m，0 m，0 m)，無極化均質(zhì)大地電阻率100 Ω·m；井旁有一極化立方低阻體，邊長8 m，中心坐標(biāo)(8 m，3 m，-28 m)；五方位測量方式各方位長度L=20 m；接收電極MN距均為1 m.

(a) (b) (c) (d)圖3 井旁低阻極化立方體有限差分正演模擬結(jié)果

正演結(jié)果如圖3所示：各方位的視電阻率極小值和視極化率極大值大致都在靠近異常體中心深度的觀測位置上；并因激發(fā)方位差異，異常曲線出現(xiàn)不同特征(各方位曲線特征與文獻[19]利用有限元法模擬結(jié)果相近).分析二次場電位差可知各方位測量曲線形態(tài)與視極化率曲線相近，但對于該地電模型，井口方位所測二次場幅值最大，顯示出不同方位對于異常體不同的極化效果[1，17].總之，不同方位的激電測量結(jié)果差別明顯，可用來推斷目標(biāo)體位置等信息；視電阻率和一次場電位差主要反應(yīng)場源對異常體的激發(fā)情況，視極化率和二次場電位差主要反應(yīng)場源對異常體的極化情況.

2 擬牛頓最優(yōu)化反演方法

本節(jié)基于前文正演模擬中對地-井IP測量特征的認(rèn)識，采用穩(wěn)流電場中規(guī)則形體解析式進行反演目標(biāo)函數(shù)和偏導(dǎo)數(shù)矩陣的構(gòu)建，應(yīng)用擬牛頓法建立最優(yōu)化擬合反演方法.

2.1 反演基本思路

首先基于最小二乘原則建立反演目標(biāo)函數(shù)[16]

(2)

本文以極化球體為例，選取其理論解析式作為目標(biāo)函數(shù)，均勻無極化介質(zhì)中井旁球體極化場電位解析式見文獻[1].經(jīng)試算可知：雖然極化球體解析式涉及不同系數(shù)的高次勒讓德多項式，但可用少數(shù)幾個多項式近似代替；目標(biāo)函數(shù)可視為一個多變量高次函數(shù)，基于此的反演工作也就相當(dāng)于在限定條件下利用最優(yōu)化法求解一個多變量的高次函數(shù).

2.2 共軛梯度法與變尺度法

牛頓法優(yōu)點是收斂速度快，但需要計算Hessian陣(H)的逆矩陣H-1，為此研究者采用共軛梯度法和變尺度法對其進行了改進.兩種方法都是選擇一個擬矩陣，使其盡可能近似于H-1，避免直接求取逆矩陣[24].

共軛梯度(Conjugate Gradient，CG)法的基本原理是將共軛性和梯度法相結(jié)合，利用所求模型參數(shù)向量在已知迭代點上的梯度方向形成一系列線性無關(guān)彼此成H共軛的向量ρM，并沿該方向進行搜索，形成每次迭代遞推格式，求出目標(biāo)函數(shù)的極小點，獲取所求模型參數(shù)的最優(yōu)化解.其迭代式可寫為

v(M)=v(1)+β1ρ1+β2ρ2+…+βM-1ρM-1.

(3)

CG法仍要計算Hessian陣，研究者提出許多變種算法，其中較簡便的如Fletcher-Reeves(FR)法，根據(jù)式(4)將基本共軛法進一步簡化為FRCG法[25].

(4)

(5)

(6)

(7)

2.3 具體反演方法

本文的反演試算采用FRCG法(Fletcher-Reeves Conjugate-Gradient)和BFGS法(Boyden-Fletcher-Goldfarb-Shanno)，觀測數(shù)據(jù)為地-井IP方位(選取主、反、兩個輔助方位)測量的一次場電位差和極化場電位差，分別對井旁異常體的一次場和極化場進行反演，待反演模型參數(shù)為：異常體電阻率(極化電阻率ρ*或一次場電阻率ρ2)，異常體中心坐標(biāo)(xp，yp，zp)，球體半徑a；即五個待求參數(shù).最優(yōu)化反演方法存在多解性，地-井IP所測數(shù)據(jù)量又較少，基于對地-井方位測量特征的認(rèn)識，適當(dāng)加入一定的限制條件對所要反演的參數(shù)向量進行約束，可使反演結(jié)果更精確.圖4為FRCG和BFGS最優(yōu)化反演流程圖.

(a) FRCG最優(yōu)化法反演程序框圖 (b) BFGS最優(yōu)化法反演程序框圖圖4 最優(yōu)化法反演程序框圖

3 反演算例與分析

本節(jié)應(yīng)用擬牛頓最優(yōu)化反演方法，對三維地電模型下井旁球體的方位、埋深等參數(shù)進行正演擬合反演，并對該方法的可行性和適用性進行簡要分析.

3.1 試算算例

算例1如圖1(c)所示，地電模型為非極化均質(zhì)大地(電阻率100 Ω·m)；垂直井位于原點，井旁有一半徑4 m的極化球體，中心坐標(biāo)(xp=8 m，yp=7 m，zp=-15 m)，球體電阻率ρ2=10 Ω·m，極化率50%；五方位測量場源為A0、A1、A2、A3、A4，電流強度20 A，方位長度L=20 m.

地-井IP觀測數(shù)據(jù)分別為一次場電位差和極化場電位差，由極化球體的解析式求解得出，并加入1%的隨機噪音.應(yīng)用最優(yōu)化法進行反演，利用極化球體解析式進行正演計算，分別對一次場和極化場電位差數(shù)據(jù)進行反演.

(a) 主剖面 (b) 輔助剖面 (c) 主剖面 (d) 輔助剖面圖5 最優(yōu)化反演擬合結(jié)果

反演初始模型參數(shù)為：球心坐標(biāo)(xp=4 m，yp=4 m，zp=-5 m)，ρ2=20 Ω·m，ρ*=30 Ω·m，a=2 m.反演擬合結(jié)果如圖5所示.

應(yīng)用BFGS法進行反演(迭代次數(shù)20)，一次場各參數(shù)反演結(jié)果：ρ2=10.325 Ω·m，a=4.329 m，中心坐標(biāo)為(xp=7.992 m，yp=7.091 m，zp=-14.838 m)； 極化總場各參數(shù)反演結(jié)果：ρ*=20.807 Ω·m，a=4.281 m，中心坐標(biāo)為(xp=8.101 m，yp=7.130 m，zp=-14.937 m).

應(yīng)用FRCG法進行反演(迭代次數(shù)20)，一次場各參數(shù)反演結(jié)果：ρ2=10.441 Ω·m，a=4.294 m，中心坐標(biāo)為(xp=8.209 m，yp=7.056 m，zp=-14.931 m)；極化總場各參數(shù)反演結(jié)果：ρ*=20.750 Ω·m，a=4.321 m，中心坐標(biāo)為(xp=8.223 m，yp=7.103 m，zp=-14.815 m).

算例2仍采用圖1(c)所示地電模型，地-井IP五方位觀測數(shù)據(jù)由有限差分法進行正演計算(地電模型為40 m×40 m×40 m，網(wǎng)格節(jié)點為81×81×81)得出，分別為一次場電位差和極化場電位差，并同樣基于有限差分正演計算進行反演，模型參數(shù)與前同.由于是均勻網(wǎng)格節(jié)點組成，每次反演所得方位數(shù)據(jù)均轉(zhuǎn)化在與之最近的網(wǎng)格節(jié)點上.

反演初始模型參數(shù)為：球心坐標(biāo)(xp=4 m，yp=4 m，zp=-5 m)，ρ2=15 Ω·m，ρ*=30 Ω·m，a=3 m.反演擬合結(jié)果如圖6所示.

(a) 主剖面 (b) 輔助剖面 (c) 主剖面 (d) 輔助剖面圖6 最優(yōu)化反演擬合結(jié)果

應(yīng)用BFGS法進行反演(迭代次數(shù)6)，一次場各參數(shù)反演結(jié)果：a=4.5 m，ρ2=10.95 Ω·m，中心坐標(biāo)為(xp=8 m，yp=7.5 m，zp=-13.5 m)；極化總場反演結(jié)果：ρ*=21.62 Ω·m，a=4.5 m，中心坐標(biāo)為(xp=8 m，yp=7.5 m，zp=-14 m).

應(yīng)用FRCG法進行反演(迭代次數(shù)6)，一次場各參數(shù)反演結(jié)果：a=4.5 m，ρ2=11.41 Ω·m，中心坐標(biāo)為(xp=8 m，yp=7.5 m，zp=-13.5 m)；極化總場參數(shù)反演結(jié)果：ρ*=20.96 Ω·m，a=4.5 m，中心坐標(biāo)為(xp=8.5 m，yp=7 m，zp=-14 m).需要說明，反演結(jié)果都是小數(shù)，而求解區(qū)域為均勻網(wǎng)格剖分，坐標(biāo)點和球體半徑參數(shù)都變?yōu)榕c之最鄰近的節(jié)點數(shù)值，在數(shù)值上都根據(jù)剖分情況有所取舍，于是造成兩種擬牛頓方法反演結(jié)果近似，反演擬合數(shù)據(jù)曲線幾乎重合，實際上各個參數(shù)計算結(jié)果都略有差別.

3.2 算例分析

從上述算例反演擬合情況來看，所得曲線趨勢和觀測值保持一致，擬合效果較為理想，同時對待求參數(shù)的反演結(jié)果也與模型設(shè)定較為接近，說明本文選取反演方法可達到較好的擬合效果，具有一定的適用性；同時具體反演參數(shù)結(jié)果和所設(shè)模型仍有一定差別，因為最優(yōu)化法的多解性，在所選限定條件內(nèi)，還有能達到相近結(jié)果的參數(shù)組合；經(jīng)過試算，適當(dāng)加入限定條件，反演參數(shù)結(jié)果可進一步優(yōu)化.

4 結(jié) 論

本文利用有限差分法，以極化球體解析式為比對，對三維地電模型中地-井時域激電異常特征進行了正演模擬；應(yīng)用共軛法、變尺度最優(yōu)化法進行了基于解析法正演與有限差分正演的反演研究，獲得了如下結(jié)論和認(rèn)識：

應(yīng)用擬牛頓法(BFGS法、FRCG法)，在適當(dāng)?shù)南薅l件下，利用地-井激電方位測量數(shù)據(jù)，對井旁規(guī)則形體或類似規(guī)則形體進行多參數(shù)的最優(yōu)化反演是可行的，獲得的反演結(jié)果較好，擬合精度較高.

擬牛頓最優(yōu)化反演法擁有較好的收斂性，具有一定的適用性，但對于反演初值設(shè)定要求較高，在進行反演工作時按情況選取較為接近結(jié)果的數(shù)據(jù)可以在保證擬合精度的前提下減少反演迭代次數(shù).

實際地質(zhì)條件下的井旁目標(biāo)體往往不是規(guī)則形體，電性和大小等參數(shù)也難以加入合適的限定條件，故需研究適用性更廣的解釋模型；此外，多參數(shù)反演較為復(fù)雜，各個待求參數(shù)對整個觀測數(shù)據(jù)的影響有較大差異，對此進行研究可進一步提高改進反演方法技術(shù)；考慮到實際工作中地質(zhì)情況較復(fù)雜，應(yīng)適當(dāng)采取人機互動的方式進行擬合反演.

[1] 蔡柏林，黃智輝，谷守民.井中激發(fā)極化法[M].北京:地質(zhì)出版社，1983．

[2] STEPHEN T，MUDGE.Radial resistivity/IP surveys using a down-hole current electrode[J].Exploration Geophysics，2004，35:188-193．

[3] 高長榮.井中激電在西霞礦區(qū)的應(yīng)用[J].物探與化探，2007，31(S1):98-101．

GAO Changrong.The application of the borehole IP method in the XiXia ore district[J].Geophysical & Geochemical Exploration.2007，31(S1):98-101.(in Chinese)

[4] LI Y G，OLDENBURG D W.3-D inversion of induced polarization data[J].Geophysics，2000，65(6):1931-1945．

[5] WELLER A，F(xiàn)RANGOS W，SEICHTER M.3-D inversion of IP data from simulated waste[J].Journal of Applied Geophysics，2000，44:67-83．

[6] ZHOU B，GREENHALGH S A.Rapid 2D/3D crosshole resistivity imaging using the analytic sensitivity function[J].Geophysics，2002，67(2):755-765．

[7] WILKINSON P B，CHAMBERS J E，et al.Extreme sensitivity of crosshole electrial resistivity tomography measurements to geometric errors[J].Geophy J Int，2008，173:49-62．

[8] 阮百堯，村上裕，徐世浙.激發(fā)極化數(shù)據(jù)的最小二乘二維反演方法[J].地球科學(xué)，1999，24(6):619-624．

RUAN Baiyao，YUTAKA M，XU Shizhe.Least square 2-D inversion method for induced polarization data[J].Geophysical & Geochemical Exploration,1999，24(6):619-624.(in Chinese)

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

WU Xiaoping，XU Guoming.Study on 3-D resistivity inversion using conjugate gradient method[J].Chinese J Geophys，2000，43(3):420-470.(in Chinese)

[10] 呂玉增，阮百堯，黃俊革.直流電井間三維直接成像[J].物探化探計算技術(shù)，2003，25(1):60-64．

LV Yuzeng，RUAN Baiyao，HUANG Junge.The 3-D immediate corsshole tomography with direct current[J].Computing Techniques for Geophysical & Geochemical Exploration.2003，25(1):60-64.(in Chinese)

[11] 沈 平，強建科，李永軍，等.井間視電阻率的幾何成像方法[J].中南大學(xué)學(xué)報:自然科學(xué)版，2010，41(3):1079-1084．

SHEN Ping，QIANG Jianke，LI Yongjun，et al. Geometry image method of crosshole apparent resistivity[J].Journal of Central South University:Science and Technology，2010，41(3):1079-1084.(in Chinese)

[12] 岳建華，劉志新.井-地三維電阻率成像技術(shù)[J].地球物理學(xué)進展，2005，20(2):407-411．

YUE Jianhua，LIU Zhixin.Three dimension resistivity tomography of mine ground[J].Progress in Geophysics，2005，20(2):407-411.(in Chinese)

[13] 安 然，李桐林，徐凱軍.井-地三維電阻率反演研究[J].地球物理學(xué)進展，2007，22(1):247-249．

AN Ran，LI Tonglin，XU Kaijun.Well-surface 3-D resistivity inversion[J].Progress in Geophysics，2007，22(1):247-249.(in Chinese)

[14] 周 峰，潘和平，吳國平，等.井中激電地-井方式井旁球體正反演.物探與化探，2008，32(3):321-325．

ZHOU Feng，PAN Heping，WU Guoping，et al.Forward and inversion of sphere beside well in then method of ground-well induced polarization[J].Geophysical & Geochemical Exploration，2008，32(2):321-325.(in Chinese)

[15] 潘和平，黃智輝.巖性和煤質(zhì)最優(yōu)化變尺度法分析[J].物探與化探，1991，15(3):168-173．

PAN Heping，HUANG Zhihui.Analysis of lithology and coal quality with the optimized scale transformation method[J].Geophysical & Geochemical Exploration，1991，15(3):168-173.(in Chinese)

[16] 潘和平，馬火林，蔡柏林，等.地球物理測井與井中物探[M].北京:科學(xué)出版社，2009．

[17]黃智輝.井中激電地-井方式方位測量資料解釋問題的探討[J].物探與化探，1979，3(3):22-30．

HUANG Zhihui.Discussion on data interpretation of azimuth measurement in the surface-hole mode of IP borehole[J].Geophysical & Geochemical Exploration.1979，3(3):22-30.(in Chinese)

[18] 周 峰，潘和平，文國軍，等.基于有限差分的井中激發(fā)極化法正演[J].電波科學(xué)學(xué)報，2010，25(4):785-791．

ZHOU Feng，PAN Heping，WEN Guojun，et al.Three-dimensional forward modeling of induced of polartization borehole using finite difference method[J].Chinese Journal of Radio Science，2010，25(4):785-791.(in Chinese)

[19] 呂玉增，阮百堯，彭蘇萍.地-井方位激電觀測異常特征研究[J].地球物理學(xué)進展，2012，27(1):201-216．

Lü Yuzeng，RUAN Baiyao，PENG Suping.A study on anomaly of surface-borehole direction induced polarization survey[J].Progress in Geophysics，2012，27(1):201-216.(in Chinese)

[20] 孫向陽，聶在平，李愛勇，等.用于電磁感應(yīng)建模的一種快速有效計算方法[J].電波科學(xué)學(xué)報，2008，23(5):932-936．

SUN Xiangyang，NIE Zaiping，LI Aiyong，et al.A fast and effective computational method of modeling the electromagnetic induction[J].Chinese Journal of Radio Science，2008，23(5):932-936.(in Chinese)

[21] 孫向陽，聶在平，李愛勇，等.用高階疊層有限元法計算隨鉆測井的三維電磁響應(yīng)[J].電波科學(xué)學(xué)報，2009，24(2):273-279．

SUN Xiangyang，NIE Zaiping，LI Aiyong，et al.The modeling of logging while drilling tool’s 3-D electromagnetic response using the high order hierarchical vector finite element method[J].Chinese Journal of Radio Science，2009，24(2):273-279.(in Chinese)

[22] 吳建平，王正華，李曉梅.稀疏線性方程組的高效求解與并行計算[M].長沙：湖南科技出版社，2004．

[23] 羅延鐘，張桂青.電子計算機在電法勘探中的應(yīng)用[M].武漢:武漢地質(zhì)學(xué)院出版社，1987．

[24] 王家映.地球物理反演理論[M].北京:高等教育出版社，2002．

[25] 袁亞湘，孫文瑜.最優(yōu)化理論與方法[M].2版.北京：科學(xué)出版社，1997．