潘和平 孟慶鑫

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

引 言

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

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

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

1 地-井時(shí)域激電三維正演

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

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

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

(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)定點(diǎn)電流場(chǎng)下井旁球狀極化體的解析式作為有限差分模擬的參考依據(jù)和反演工作中的正演計(jì)算，其表達(dá)式可見(jiàn)文獻(xiàn)[1].

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

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

算例1地電模型如圖1(b)所示，無(wú)極化均質(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%；場(chǎng)源位于A1方位坐標(biāo)(x=11 m，y=0 m，z=0 m)，電流強(qiáng)度20 A；井Zk平面坐標(biāo)(x=0 m，y=0 m).

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

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

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

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

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

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

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

2.1 反演基本思路

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

(2)

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

2.2 共軛梯度法與變尺度法

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

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

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

(3)

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

(4)

(5)

(6)

(7)

2.3 具體反演方法

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

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

3 反演算例與分析

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

3.1 試算算例

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

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

(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法進(jìn)行反演(迭代次數(shù)20)，一次場(chǎng)各參數(shù)反演結(jié)果：ρ2=10.325 Ω·m，a=4.329 m，中心坐標(biāo)為(xp=7.992 m，yp=7.091 m，zp=-14.838 m)； 極化總場(chǎng)各參數(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法進(jìn)行反演(迭代次數(shù)20)，一次場(chǎng)各參數(shù)反演結(jié)果：ρ2=10.441 Ω·m，a=4.294 m，中心坐標(biāo)為(xp=8.209 m，yp=7.056 m，zp=-14.931 m)；極化總場(chǎng)各參數(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五方位觀測(cè)數(shù)據(jù)由有限差分法進(jìn)行正演計(jì)算(地電模型為40 m×40 m×40 m，網(wǎng)格節(jié)點(diǎn)為81×81×81)得出，分別為一次場(chǎng)電位差和極化場(chǎng)電位差，并同樣基于有限差分正演計(jì)算進(jìn)行反演，模型參數(shù)與前同.由于是均勻網(wǎng)格節(jié)點(diǎn)組成，每次反演所得方位數(shù)據(jù)均轉(zhuǎn)化在與之最近的網(wǎng)格節(jié)點(diǎn)上.

反演初始模型參數(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法進(jìn)行反演(迭代次數(shù)6)，一次場(chǎng)各參數(shù)反演結(jié)果：a=4.5 m，ρ2=10.95 Ω·m，中心坐標(biāo)為(xp=8 m，yp=7.5 m，zp=-13.5 m)；極化總場(chǎng)反演結(jié)果：ρ*=21.62 Ω·m，a=4.5 m，中心坐標(biāo)為(xp=8 m，yp=7.5 m，zp=-14 m).

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

3.2 算例分析

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

4 結(jié) 論

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

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

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

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

[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] 高長(zhǎng)榮.井中激電在西霞礦區(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é)報(bào)，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].物探化探計(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] 沈 平，強(qiáng)建科，李永軍，等.井間視電阻率的幾何成像方法[J].中南大學(xué)學(xué)報(bào):自然科學(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é)進(jìn)展，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é)進(jìn)展，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] 周 峰，潘和平，吳國(guó)平，等.井中激電地-井方式井旁球體正反演.物探與化探，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] 潘和平，馬火林，蔡柏林，等.地球物理測(cè)井與井中物探[M].北京:科學(xué)出版社，2009．

[17]黃智輝.井中激電地-井方式方位測(cè)量資料解釋問(wèn)題的探討[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] 周 峰，潘和平，文國(guó)軍，等.基于有限差分的井中激發(fā)極化法正演[J].電波科學(xué)學(xué)報(bào)，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] 呂玉增，阮百堯，彭蘇萍.地-井方位激電觀測(cè)異常特征研究[J].地球物理學(xué)進(jìn)展，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)，聶在平，李愛(ài)勇，等.用于電磁感應(yīng)建模的一種快速有效計(jì)算方法[J].電波科學(xué)學(xué)報(bào)，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)，聶在平，李愛(ài)勇，等.用高階疊層有限元法計(jì)算隨鉆測(cè)井的三維電磁響應(yīng)[J].電波科學(xué)學(xué)報(bào)，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] 吳建平，王正華，李曉梅.稀疏線性方程組的高效求解與并行計(jì)算[M].長(zhǎng)沙：湖南科技出版社，2004．

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

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

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