基于改進(jìn)Krylov子空間算法的井中激電反演

李長(zhǎng)偉，熊 彬，王有學(xué)，羅潤(rùn)林

桂林理工大學(xué)地球科學(xué)學(xué)院，桂林 541004

1 引 言

井中激發(fā)極化法，簡(jiǎn)稱井中激電，作為在鉆孔中進(jìn)行激發(fā)極化測(cè)量的一種非常有效的地球物理勘探方法，在金屬礦勘探中得到了越來(lái)越多的應(yīng)用.井中激電反演包括了電阻率和極化率的反演，極化率數(shù)據(jù)的反演可以在電阻率反演的基礎(chǔ)上得到[1-2］.有限元法由于其對(duì)復(fù)雜模型的處理能力，在地球物理電磁法模擬中得到廣泛的應(yīng)用[3-7］.通過(guò)將電位分解成正常電位和異常電位，Coggon[3］，Lowry等[8］給出了在正演模擬中能消去源點(diǎn)奇異性影響的異常電位法，Zhao等[9］進(jìn)一步指出正常場(chǎng)的電導(dǎo)率應(yīng)為點(diǎn)源處的電導(dǎo)率取值.借鑒異常電位法的思想，Pidlisecky等[10］在有限體積法的正演模擬中給出了消去邊界效應(yīng)和源點(diǎn)奇異性的右端項(xiàng)校正技術(shù).對(duì)于有限元方程，由于其系數(shù)矩陣的對(duì)稱正定陣性，常采用預(yù)處理共軛梯度法（PCG）求解[11-12］.Gauss-Newton法和共軛梯度法是求解反演優(yōu)化問(wèn)題的常用的方法，結(jié)合Krylov子空間算法，可以不進(jìn)行Jacobian矩陣的存儲(chǔ)，直接計(jì)算Jacobian矩陣向量積來(lái)實(shí)現(xiàn)迭代求解[13-14］.有限元法不要求規(guī)則的網(wǎng)格剖分，允許采用非結(jié)構(gòu)化網(wǎng)格來(lái)離散模擬區(qū)域，Rücker等[15-16］采用非結(jié)構(gòu)化網(wǎng)格實(shí)現(xiàn)三維電阻率有限元正反演，提高了對(duì)復(fù)雜模型的處理能力.

本文在正演模擬時(shí)，考慮到井眼影響，采用結(jié)構(gòu)化和非結(jié)構(gòu)化網(wǎng)格相結(jié)合的方式對(duì)模擬區(qū)域剖分.將右端項(xiàng)校正技術(shù)應(yīng)用到有限元法的正演模擬中以減少邊界效應(yīng)和源點(diǎn)奇異性引起的誤差.Krylov子空間法是求解大型稀疏線性方程組的常用方法.實(shí)際勘探中一般會(huì)采用多個(gè)點(diǎn)源電極布置，對(duì)應(yīng)每個(gè)源電極位置都會(huì)形成一個(gè)有限元方程，從而形成了多個(gè)線性方程組.我們針對(duì)這多個(gè)線性方程組系數(shù)矩陣之間差異不大的特點(diǎn)，采用循環(huán)Krylov子空間算法提高多線性方程組的求解效率.反演采用粗網(wǎng)格剖分，用以降低反演問(wèn)題的自由度和不適定性，正演網(wǎng)格在反演網(wǎng)格的基礎(chǔ)上加密得到，以保證正演的計(jì)算精度.在反演迭代中，用不精確PCG算法近似求解模型修正量方程減少了計(jì)算量.根據(jù)剛度矩陣與電導(dǎo)率向量參數(shù)的線性關(guān)系，進(jìn)一步簡(jiǎn)化了Jacobian矩陣向量積的計(jì)算公式，避免了存儲(chǔ)剛度矩陣對(duì)電導(dǎo)率的顯式求導(dǎo)過(guò)程.

2 正演問(wèn)題

2.1 基本原理

2.1.1 點(diǎn)源場(chǎng)基本原理

地下穩(wěn)恒電流產(chǎn)生的電場(chǎng)滿足下面的微分方程[17］：

我們采用有限元方法進(jìn)行正演模擬，對(duì)無(wú)限區(qū)域做截?cái)嗵幚?，施加人工邊界條件

轉(zhuǎn)化為變分問(wèn)題[17］

（3）中第一項(xiàng)是體積分項(xiàng)，第二項(xiàng)是邊界積分項(xiàng).其中r是點(diǎn)源到邊界的距離，n是邊界外法線方向.

用有限元方法對(duì)問(wèn)題離散后得到線性方程組

（4）被稱為有限元方程，其中K是離散正演算子，被稱為剛度矩陣，u是節(jié)點(diǎn)電位未知量，eisrc是一個(gè)點(diǎn)源節(jié)點(diǎn)處分量為1，其它分量為0的向量，稱為有限元方程的右端項(xiàng).

2.1.2 視極化率的計(jì)算

用ρ表示電阻率，η表示極化率，ρs表示視電阻率，ηs表示視極化率，Δu，Δu1，Δu2分別表示極化場(chǎng)、一次場(chǎng)、二次場(chǎng)電位.根據(jù)Seigel理論[1］，等效電阻率為

視極化率為

在正演時(shí)，當(dāng)求得一次場(chǎng)電位后，用等效電阻率代替相應(yīng)的電阻率，再進(jìn)行一次求解，便得到極化場(chǎng)電位，利用（6）式即可求得視極化率.

2.2 網(wǎng)格剖分

結(jié)構(gòu)化網(wǎng)格的節(jié)點(diǎn)排列有序，每個(gè)節(jié)點(diǎn)和鄰點(diǎn)的關(guān)系固定不變，網(wǎng)格生成速度快，數(shù)據(jù)結(jié)構(gòu)簡(jiǎn)單，但是不能方便地進(jìn)行局部網(wǎng)格加密.非結(jié)構(gòu)化網(wǎng)格的節(jié)點(diǎn)編號(hào)命名并無(wú)一定規(guī)則，可以靈活地進(jìn)行網(wǎng)格局部加密，能夠有效地減少網(wǎng)格剖分單元.一般來(lái)說(shuō)，井眼（直徑為數(shù)十厘米）相對(duì)于整個(gè)剖分區(qū)域很小，對(duì)于這種窄狹區(qū)域，非結(jié)構(gòu)化網(wǎng)格為了保證網(wǎng)格質(zhì)量會(huì)導(dǎo)致大量的對(duì)提高精度意義不大的節(jié)點(diǎn)密集于這些區(qū)域，造成計(jì)算量的浪費(fèi)和生成網(wǎng)格的困難[18］.因此在考慮井眼影響時(shí)，我們?cè)诰蹍^(qū)域以井為中心采用放射狀結(jié)構(gòu)化網(wǎng)格，在外圍區(qū)域采用非結(jié)構(gòu)化網(wǎng)格剖分；若不考慮井眼影響時(shí)則對(duì)整個(gè)區(qū)域進(jìn)行非結(jié)構(gòu)化網(wǎng)格剖分.網(wǎng)格剖分見(jiàn)圖1.

圖1 網(wǎng)格剖分：井眼區(qū)域（a），外圍區(qū)域（b）Fig.1 Finite element mesh：borehole（a），outside（b）

2.3 有限元方程的右端項(xiàng)校正

在人工截?cái)噙吔缟喜捎茫?）式中的混合邊界條件可在合理的計(jì)算范圍內(nèi)得到較高的正演模擬精度，但仍需要將邊界選取在足夠遠(yuǎn)的地方以減少邊界效應(yīng)的影響.此外，源點(diǎn)奇異性的存在造成源附近的模擬誤差較大，異常電位法可以有效地降低源點(diǎn)奇異性引起的誤差[8］.異常電位法得到的有限元方程為

式中us為異常電位，up為正常電位，u=up+us，σp為點(diǎn)源處的電導(dǎo)率值.設(shè)源節(jié)點(diǎn)的編號(hào)為s，若源點(diǎn)附近網(wǎng)格剖分足夠精細(xì)，在所有包含點(diǎn)源節(jié)點(diǎn)的網(wǎng)格單元上有σp-σ=0，K（σp-σ）的第s行和第s列的元素均為0.

（7）式等價(jià)于

與原有限元方程相比，（8）式采用了新的右端項(xiàng)，能夠補(bǔ)充正演算子由于邊界取得不夠遠(yuǎn)所丟失的信息，減少邊界效應(yīng)和源點(diǎn)奇異性引起的誤差[10］.均勻半空間情況下有限元?jiǎng)偠染仃嘖（σ）為pσp的線性函數(shù)，K（σp）up與σp的取值無(wú)關(guān).

2.4 多線性方程組求解技術(shù)

實(shí)際的勘探工作往往涉及到多個(gè)點(diǎn)源電極布置，每一個(gè)點(diǎn)源位置對(duì)應(yīng)不同的有限元方程，則需要求解多個(gè)線性方程組.考察變分問(wèn)題（3），設(shè)點(diǎn)源個(gè)數(shù)為nsrc，將有限元方程（4）進(jìn)一步寫(xiě)為

其中K0對(duì)應(yīng)體積分項(xiàng)，ΔKj對(duì)應(yīng)第j個(gè)點(diǎn)源的邊界積分項(xiàng).ΔKj是比K稀疏得多的低秩矩陣，將兩項(xiàng)分開(kāi)并分別進(jìn)行存儲(chǔ)，只需要存儲(chǔ)K0和nsrc個(gè)ΔKj，可減少存儲(chǔ)需求.

多線性方程組（9）的系數(shù)矩陣間只存在邊界積分項(xiàng)的不同，相互間差異不大，可以先選擇其中的一個(gè)線性方程組作為種子系統(tǒng)預(yù)先求解，存儲(chǔ)在求解過(guò)程中生成的子空間信息，在求解其余線性方程組時(shí)，將初始?xì)埩亢退阉鞣较蛟诖俗涌臻g上進(jìn)行投影以加快算法的收斂速度.當(dāng)多線性方程組的右端項(xiàng)靠近，即右端項(xiàng)向量間的夾角小時(shí)，算法有更好的收斂效果.由于（9）的右端項(xiàng)向量相互正交.設(shè)第k個(gè)線性方程組被選為種子系統(tǒng)，為改善右端項(xiàng)的靠近程度，構(gòu)造一個(gè)新的多線性方程組

其中b0= （1，1，…，1）T是一個(gè)分量全為1的向量，bi= （1，1，…，1，0，1，…，1）T是一個(gè)對(duì)應(yīng)點(diǎn)電源位置的分量為0，其余分量全為1的向量，方程（10）和（9）有如下關(guān)系：

（10）式的右端項(xiàng)之間滿足

其中n是系數(shù)矩陣的維數(shù)，ei，ui分別為原多線性方程組（9）的右端項(xiàng)和未知項(xiàng)向量.由于（ΔAi-ΔAk）的范數(shù)很小，式（12）表明越大的n意味著新的多線性方程組的右端項(xiàng)越靠近，因此也應(yīng)該有更好的收斂性.算法的具體實(shí)現(xiàn)見(jiàn)參考文獻(xiàn)[19］.

3 反演問(wèn)題

3.1 目標(biāo)函數(shù)

井中激電的主要反演參數(shù)為電阻率和極化率，本文利用正則化原理和Gauss-Newton法反演.設(shè)m= （m1，m2，…，mM）T為模型參數(shù)向量，d= （d1，d2，…，dN）T為數(shù)據(jù)參數(shù)向量，數(shù)據(jù)誤差為ε1，ε2，…，εN.取m=logσ，d=logρs，這里的σ（=1／ρ），ρs分別代表電導(dǎo)率和視電阻率，ρ為電阻率.給定目標(biāo)函數(shù)為

其中λ為正則化因子，d（m）表示給定模型參數(shù)由正演計(jì)算得到的數(shù)據(jù)，dobs表示實(shí)際測(cè)量數(shù)據(jù)，D=diag（1／εi）為數(shù)據(jù)加權(quán)矩陣，W 為模型光滑性約束矩陣，mref是根據(jù)先驗(yàn)信息給出的參考模型.

利用Gauss-Newton法求目標(biāo)函數(shù)的極小值，每步迭代更新如下：

其中α為線性搜索步長(zhǎng)，Δmk為模型修正量.通過(guò)求解下面的線性方程組得到

其中J=?d／?m，是Jacobian矩陣.方程（15）中光滑性約束W控制了模型單元在空間中的平穩(wěn)變化程度，一般取為模型的離散差分算子.若對(duì)模型做一階正則化約束，即最平模型估計(jì)，W 取為梯度算子.在三維情況下，需要在三個(gè)坐標(biāo)方向上考慮，取

其中αx，αy，αz為三個(gè)坐標(biāo)方向上的加權(quán)，用于控制在不同方向上的光滑程度，需要根據(jù)先驗(yàn)信息進(jìn)行選擇，如令αx取較大值時(shí)意味著解在x方向變化更平緩.

3.2 網(wǎng)格剖分

非結(jié)構(gòu)化網(wǎng)格的生成過(guò)程比較復(fù)雜，TetGen是一款免費(fèi)開(kāi)源的四面體網(wǎng)格生成工具，可以方便地定義分區(qū)網(wǎng)格剖分和設(shè)置網(wǎng)格控制節(jié)點(diǎn)，本文中利用TetGen來(lái)生成非結(jié)構(gòu)化網(wǎng)格.有限元正演的精度與網(wǎng)格剖分精細(xì)程度密切相關(guān).為保證正演結(jié)果的精度，正演計(jì)算的網(wǎng)格需要剖分的足夠精細(xì)；反演網(wǎng)格定義了模型參數(shù)的數(shù)目和區(qū)域的電性劃分，若剖分過(guò)細(xì)，會(huì)造成測(cè)量數(shù)據(jù)個(gè)數(shù)和模型參數(shù)個(gè)數(shù)的比例更小，增大了反演過(guò)程的不適定性，并且增大了計(jì)算量.我們對(duì)反演網(wǎng)格進(jìn)行較粗的剖分，采用結(jié)構(gòu)化三棱柱單元或非結(jié)構(gòu)化四面體單元.在每一個(gè)反演網(wǎng)格單元內(nèi)進(jìn)行加細(xì)剖分得到若干個(gè)正演網(wǎng)格單元，由同一個(gè)反演單元剖分得到的正演網(wǎng)格單元具有相同的模型參數(shù)屬性.正反演采用不同的網(wǎng)格系統(tǒng)，需要處理的一個(gè)重要問(wèn)題是正演和反演過(guò)程中的網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)傳遞.由反演網(wǎng)格加密得到正演網(wǎng)格過(guò)程中，原反演網(wǎng)格的節(jié)點(diǎn)進(jìn)行保留并保持編號(hào)不變，新增加的節(jié)點(diǎn)編號(hào)總是排在后面，這樣的編號(hào)方式使正反演網(wǎng)格之間的數(shù)據(jù)傳遞不需要復(fù)雜的處理進(jìn)行對(duì)應(yīng).

3.3 Jacobian矩陣計(jì)算

利用Krylov子空間法求解方程組（15），只需要計(jì)算Jacobian矩陣向量積，不必要直接計(jì)算和存儲(chǔ)Jacobian矩陣，稱為Jacobian-free Krylov方法[20-22］.利用求導(dǎo)的鏈?zhǔn)椒▌t，對(duì)有限元方程（4）或（8）兩邊求導(dǎo)

于是

記

有

對(duì)任一向量

設(shè)ei表示第i個(gè)分量為1，其余分量為0的單位向量，由于剛度矩陣是電導(dǎo)率向量參數(shù)的線性算子，有

于是

對(duì)任意向量

其中z=K-1v.對(duì)于Jacobian矩陣轉(zhuǎn)置與向量乘積，有

其中x=QTPy.實(shí)際計(jì)算GTv時(shí)利用一次單元?jiǎng)偠染仃嚿蛇^(guò)程便可得到它的所有分量，只需要在每個(gè)有限單元上進(jìn)行積分時(shí)設(shè)置電導(dǎo)率為1，在計(jì)算的同時(shí)乘以相應(yīng)的ziuj，并只對(duì)屬于同一反演單元的正演單元?jiǎng)偠染仃囘M(jìn)行累加.

（26），（28）和（29）式利用剛度矩陣表示，給出了Jacobian矩陣及其轉(zhuǎn)置與向量的乘積的更為方便的計(jì)算方法，避免了顯式地求取和存儲(chǔ)剛度矩陣對(duì)電導(dǎo)率的導(dǎo)數(shù).由于采用了精確的單元積分公式，而不采用數(shù)值積分，在計(jì)算Jv，JTv時(shí)做一次單元積分所花費(fèi)的計(jì)算時(shí)間非常少.程序在Intel?CoreTM2 Duo CPU，2.66GHz，3.00GB內(nèi)存，Windows XP下運(yùn)行，對(duì)506262個(gè)四面體單元做一次單元積分生成總剛度矩陣所用CPU時(shí)間為0.75s.

在多個(gè)點(diǎn)源的情況下，設(shè)nsrc為點(diǎn)源個(gè)數(shù)，可以對(duì)每個(gè)源分別計(jì)算然后進(jìn)行集成或累加，得到Jx，JTx.設(shè)Ji為第i個(gè)點(diǎn)源對(duì)應(yīng)的Jacobian矩陣，則可將向量x進(jìn)行相應(yīng)的分解x= （x1，x2，…，xnsrc）T，有

3.4 反演算法

記

模型修正量方程（15）可寫(xiě)為

Dembo等指出，對(duì)方程組（33）只需要近似求解[23］.借鑒不精確Newton法的思想，我們采用不精確的PCG算法求解方程組（33）：設(shè)置較低的PCG迭代收斂準(zhǔn)則和迭代次數(shù)限制；在計(jì)算矩陣向量積Hv時(shí)，設(shè)置較低的正演計(jì)算精度.

式（14）中線性搜索步長(zhǎng)α應(yīng)滿足 Wolfe準(zhǔn)則[21］，對(duì)不等式

進(jìn)行估計(jì)，若不等式不滿足，則令

再進(jìn)行估計(jì)判斷.其中c為某一常數(shù)，一般取一個(gè)較小的數(shù)，如取為10-4，ρ∈（0，1）為某一常數(shù)，在本文的數(shù)值算例中我們?nèi)ˇ?0.5.

正則化因子λ提供了最佳數(shù)據(jù)擬合和最合理解的折衷，對(duì)問(wèn)題解的性態(tài)起著關(guān)鍵的作用.在測(cè)量數(shù)據(jù)的噪聲水平參數(shù)δ可獲取或能近似估計(jì)的情況下，偏差原理是一種廣為采用的正則化參數(shù)選取策略.其基本思想是好的正則解應(yīng)能使殘量范數(shù)同數(shù)據(jù)噪聲水平相符.記測(cè)量數(shù)據(jù)為dobs，設(shè) ‖dtrue-dobs‖ ≤δ，則合適的解應(yīng)滿足

整個(gè)反演步驟如下：

（1）讀入最大迭代次數(shù)限制、迭代停止準(zhǔn)則等反演參數(shù)，讀入初始模型，選擇較大的值作為初始正則化因子；

（2）利用不精確PCG算法求解模型修正量Δmk；

（3）計(jì)算搜索步長(zhǎng)α；

（4）更新模型參數(shù)mk+1=mk+αΔmk；

（5）估計(jì)殘量，根據(jù)最大迭代次數(shù)限制和迭代收斂準(zhǔn)則判斷，若滿足繼續(xù)第6步，否則轉(zhuǎn)到步驟2；

（6）根據(jù)偏差原理，判斷不等式‖A（m）-dobs‖≤δ是否成立，若不成立，則減小正則化因子λ，轉(zhuǎn)到步驟2；

（7）反演求解結(jié)束，輸出結(jié)果.

3.5 極化率反演

極化率的反演我們采用算子線性近似方法[2，24］

在電阻率反演的基礎(chǔ)上，加入光滑約束和先驗(yàn)信息，求解如下的目標(biāo)函數(shù)極小問(wèn)題：

利用PCG對(duì)（39）式求解，便得到相應(yīng)的極化率參數(shù).其中數(shù)據(jù)加權(quán)D及光滑性限制W和視電阻率反演的選擇方式相同，λ的選擇應(yīng)滿足偏差原理.

4 數(shù)值實(shí)驗(yàn)

4.1 算例1

見(jiàn)圖2，背景電阻率為ρ=200Ωm，其中有一個(gè)低阻異常體，電阻率為ρ1=50Ωm，大小為30m×30m×20m，位于X：[-40m，-10m］，Y：[-40m，-10m］，Z：[50m，70m］處；一個(gè)高阻異常體，電阻率為ρ2=500Ωm，大小為30m×30m×20m，位于X：[10m，40m］，Y：[10m，40m］，Z：[30m，50m］處.

設(shè)分別在井中不同深度A0（z=0m），A1（z=20m），A2（z=50m），A3（z=70m）處固定A 極供電，在地面11條測(cè)線上逐點(diǎn)移動(dòng)M極，在121個(gè)數(shù)據(jù)點(diǎn)上進(jìn)行測(cè)量，得到484個(gè)視電阻率數(shù).計(jì)算區(qū)域選為60m×60m×60m，將計(jì)算區(qū)域離散成12×12×12×2共3456個(gè)三棱柱體反演單元，2197個(gè)節(jié)點(diǎn)，在反演剖分的基礎(chǔ)上進(jìn)行正演網(wǎng)格剖分，得到234204個(gè)正演單元，40181個(gè)正演節(jié)點(diǎn).

給定初始模型m0和先驗(yàn)?zāi)Ｐ蚼ref為ρ=200Ωm的均勻模型，數(shù)據(jù)相對(duì)擬合差定義為

設(shè)置停止準(zhǔn)則為R＜1×10-4，迭代8次達(dá)到收斂要求，所花CPU時(shí)間為187.79s，迭代收斂曲線見(jiàn)圖3.

圖4是Y=-53m，Y=-33m，Y=-16m，Y=16m，Y=33m，Y=53m處的電阻率反演切片圖.從圖上可以明顯看出異常體的存在，異常位置反映基本準(zhǔn)確.由于反演過(guò)程中對(duì)模型加入光滑性約束，反演出的異常體電阻率值與真實(shí)電阻率值有所偏差.

4.2 算例2

計(jì)算區(qū)域和網(wǎng)格剖分及模型大小設(shè)置同算例1.低阻體的電阻率為ρ1=50Ωm，極化率η1=5%，高阻體的電阻率為ρ2=500Ωm，極化率η2=10%，背景電阻率為ρ=200Ωm，極化率為η=1%.

進(jìn)行井-井觀測(cè).在中間鉆孔（X=0，Y=0）中不同深度A0（z=0m），A1（z=20m），A2（z=50m），A3（z=70m）處固定A極供電，在其它8個(gè)鉆孔中共88個(gè)數(shù)據(jù)點(diǎn)上進(jìn)行測(cè)量，8個(gè)測(cè)量鉆孔的井口位置分別為：（X=-50m，Y=-50m）；（X=-50m，Y=50m）；（X=50m，Y=-50m）；（X=50m，Y=50m）；（X=0m，Y=30m）；（X=0m，Y=-30m）；（X=-30m，Y=0m）；（X=30m，Y=0m）.得到352個(gè)視電阻率和視極化率數(shù)據(jù).

給定初始模型m0和先驗(yàn)?zāi)Ｐ蚼ref為ρ=200Ωm的均勻模型，設(shè)置停止準(zhǔn)則為R＜1×10-5，共迭代14次達(dá)到收斂要求.在電阻率反演基礎(chǔ)上進(jìn)行極化率反演，得到XOZ面的電阻率和極化率反演切片見(jiàn)圖5，圖6.從圖上可以看到采用井-井方式，電阻率和極化率的反演效果都比較理想，異常體的形態(tài)和位置反映準(zhǔn)確.

5 結(jié) 論

隨著礦產(chǎn)資源的供需矛盾加劇，人們急需尋找第二找礦空間.井中激發(fā)極化法是深部金屬礦勘探中的一種常用方法.本文系統(tǒng)研究了井中激發(fā)極化法的正反演技術(shù)，在電阻率法正反演的基礎(chǔ)上實(shí)現(xiàn)激發(fā)極化數(shù)據(jù)的正反演.正演采用有限元方法，在正演基礎(chǔ)上實(shí)現(xiàn)了井中激電的正則化反演.文中對(duì)網(wǎng)格剖分方式、高效的正反演算法、極化率的反演和方程組求解技術(shù)等方面進(jìn)行了討論.正演計(jì)算中采用了適合測(cè)井模型的剖分方式，可以處理地形、斜井和考慮井眼影響的模擬問(wèn)題.對(duì)有限元方程進(jìn)行右端項(xiàng)校正使算法在保證計(jì)算精度的前提下有效地減少了模擬區(qū)域的網(wǎng)格剖分?jǐn)?shù).改進(jìn)的多線性方程組求解技術(shù)提高了正演的計(jì)算效率.反演采用Jacobianfree Krylov子空間技術(shù)，不需要計(jì)算和存儲(chǔ)Jacobian矩陣，進(jìn)一步推導(dǎo)了更為簡(jiǎn)便的Jacobian矩陣向量積的計(jì)算方法，避免了剛度矩陣對(duì)電導(dǎo)率求導(dǎo)的顯式計(jì)算和存儲(chǔ).反演的迭代求解過(guò)程中，采用近似求解的思想減少了迭代次數(shù)和計(jì)算時(shí)間.在上述分析的基礎(chǔ)上開(kāi)發(fā)了正反演算法程序，數(shù)值算例驗(yàn)證了算法程序的效率和可靠性.

（References）

[1]Seigel H O.Mathematical formulation and type curves for induced polarization.Geophysics，1959，24（3）：547-565.

[2]Oldenburg D W，Li Y.Inversion of induced polarization data.Geophysics，1994，59（9）：1327-1341.

[3]Coggon J H.Electromagnetic and electrical modeling by the finite element method.Geophysics，1971，36（1）：132-155.

[4]Li Y G，Spitzer K.Three-dimensional DC resistivity forward modelling using finite elements in comparison with finitedifference solutions.Geophys.J.Int.，2002，151（3）：924-934.

[5]Sasaki Y.3-D resistivity inversion using the finite-element method.Geophysics，1994，59（12）：1839-1848.

[6]Pain C C，Herwanger J V，Worthington M H，et al.Effective multidimensional resistivity inversion using finiteelement techniques.Geophys.J.Int.，2002，151（3）：710-728.

[7]阮百堯，熊彬.電導(dǎo)率連續(xù)變化的三維電阻率測(cè)深有限元模擬.地球物理學(xué)報(bào)，2002，45（1）：131-138.Ruan B Y，Xiong B.A finite element modeling of threedimensional resistivity sounding with continuous conductivity.Chinese J.Geophys.（in Chinese），2002，45（1）：131-138.

[8]Lowry T，Allen M B，Shive P N.Singularity removal：a refinement of resistivity modeling techniques.Geophysics，1989，54（6）：766-774.

[9]Zhao S K，Yedlin M J.Some refinements on the finitedifference method for 3-D DC resistivity modeling.Geophysics，1996，61（5）：1301-1307.

[10]Pidlisecky A，Haber E，Knight R.RESINVM3D：A 3D resistivity inversion package.Geophysics，2007，72（2）：H1-H10.

[11]Spitzer K.A 3-D finite-difference algorithm for DC resistivity modelling using conjugate gradient methods.Geophys.J.Int.，1995，123（3）：903-914.

[12]吳小平，汪彤彤.利用共軛梯度算法的電阻率三維有限元正演.地球物理學(xué)報(bào)，2003，46（3）：428-432.Wu X P，Wang T T.A 3-D finite-element resistivity forward modeling using conjugate gradient algorithm.Chinese J.Geophys.（in Chinese），2003，46（3）：428-432.

[13]Loke M H，Chambers J E，Ogilvy R D.Inversion of 2D spectral induced polarization imaging data.Geophysical Prospecting，2006，54（3）：287-301.

[14]吳小平，徐果明.利用共軛梯度法的電阻率三維反演研究.地球物理學(xué)報(bào)，2000，43（3）：420-427.Wu X P，Xu G M.Study on 3-D resistivity reversion using conjugate gradient method.Chinese J.Geophys.（in Chinese），2000，43（3）：420-427.

[15]Rücker C， Günther T， Spitzer K.Three-dimensional modelling and inversion of DC resistivity data incorporating topography-I.modelling.Geophys.J.Int.，2006，166（2）：495-505.

[16]Günther T， Rücker C， Spitzer K.Three-dimensional modeling and inversion of DC resistivity data incorporating topography-II.inversion.Geophys.J.Int.，2006，166：506-517.

[17]徐世浙.地球物理中的有限單元法.北京：科學(xué)出版社，1994：178-188.Xu S Z.The Finite Element Method in Geophysics（in Chinese）.Beijing：Science Press，1994：178-188.

[18]任政勇，湯井田.基于局部加密非結(jié)構(gòu)化網(wǎng)格的三維電阻率法有限元數(shù)值模擬.地球物理學(xué)報(bào)，2009，52（10）：2627-2634.Ren Z Y，Tang J T.Finite element modeling of 3-D DC resistivity using locally refined unstructured meshes.Chinese J.Geophys.（in Chinese），2009，52（10）：2627-2634.

[19]Li C W，Xiong B，Qiang J K，et al.Multiple linear system techniques for 3Dfinite element method modeling of direct current resistivity.Journal of Central South University，2012，19（2）：424-432.

[20]Knoll D A， Keyes D E.Jacobian-free Newton-Krylov methods：a survey of approaches and applications.Journal of Computational Physics，2004，193（2）：357-397.

[21]Nocedal J，Wright S J.Numerical Optimization.New York：Springer，2006.

[22]吳小平，徐果明.電阻率三維反演中偏導(dǎo)數(shù)矩陣的求取和分析.石油地球物理勘探，1999，34（4）：363-372.Wu X P，Xu G M.Derivation and analysis of partial derivative matrix in resistivity 3-D inversion.Oil Geophysical Prospecting （in Chinese），1999，34（4）：363-372.

[23]Dembo R S，Eisenstat S C，Steihaug T.Inexact Newton methods.SIAM J.Numer.Anal.，1982，19（2）：400-408.

[24]阮百堯，村上裕，徐世浙.激發(fā)極化數(shù)據(jù)的最小二乘二維反演方法.地球科學(xué)，1999，24（6）：619-624.Ruan B Y，Yutaka M，Xu S Z.Least square 2-D inversion method for induced polarization data.Earth Science（in Chinese），1999，24（6）：619-624.