劉永亮， 李桐林*， 胡英才， 高文， 朱成， 柴倫煒

1 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 1300262 核工業(yè)北京地質(zhì)研究院， 北京 100029

?

快速擬線性近似方法及三維頻譜激電反演研究

劉永亮1， 李桐林1*， 胡英才2， 高文1， 朱成1， 柴倫煒1

1 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 1300262 核工業(yè)北京地質(zhì)研究院， 北京 100029

Zhdanov提出的擬線性（QL）近似方法有效解決了積分方程計(jì)算速度慢、占用內(nèi)存多的缺陷，但因?yàn)槿S頻譜激電（SIP）換源頻繁，每次換源都需要重新計(jì)算剖分單元的并矢格林函數(shù)和一次場(chǎng)，所以用該方法模擬三維SIP的速度仍然緩慢.本文根據(jù)一次場(chǎng)及格林函數(shù)的空間對(duì)稱性，提出了一種適用于多源電磁法的快速Q(mào)L近似正演方法.在此基礎(chǔ)上，實(shí)現(xiàn)了Cole-Cole參數(shù)范圍約束的三維SIP共軛梯度反演方法.理論試算結(jié)果表明：快速Q(mào)L近似方法適用于三維SIP正演模擬，計(jì)算速度較QL近似方法有了很大提高；Cole-Cole參數(shù)范圍約束的三維SIP共軛梯度反演方法對(duì)零頻電阻率、極化率等參數(shù)的反演效果良好，將異常區(qū)域剖分成2000塊時(shí)迭代一次只需約0.2 s的時(shí)間，占用內(nèi)存大約為140 MB.

快速; 擬線性; Cole-Cole模型; 共軛梯度; 范圍約束函數(shù)

1 引言

隨著找礦難度的加大以及電磁法儀器和計(jì)算機(jī)技術(shù)的快速發(fā)展，人們對(duì)電磁法解決地質(zhì)問(wèn)題的期待已越來(lái)越高，傳統(tǒng)的一維、二維解釋模型已不能滿足對(duì)復(fù)雜地質(zhì)模型的刻畫(huà)，三維電磁法正、反演技術(shù)已提到日程，對(duì)此，國(guó)內(nèi)外研究學(xué)者作了大量研究工作（如Eaton,1989; Madden and Mackie,1989; Smith and Booker, 1991;Lee and Xie,1993; Oristaglio et al.,1993 ;Pellerin et al.,1993; Nekut,1994; Torres-Verdin and Habashy,1995; Xie and Lee,1995; Newman and Alumbaugh,1996; Oldenburg et al.,2013; 張輝,2006; 徐凱軍,2007; 李建平,2008; 翁愛(ài)華等,2012等）.但是，三維電磁法內(nèi)存消耗大、計(jì)算速度慢，目前還很難投入到實(shí)際生產(chǎn)應(yīng)用當(dāng)中.Zhdanov和Fang（1996a,1996b）提出擬線性（QL）近似法，通過(guò)對(duì)比積分方程法、QL近似方法和Born近似方法對(duì)不同模型的計(jì)算結(jié)果， 驗(yàn)證了QL近似方法的計(jì)算精度，有效解決了積分方程計(jì)算速度慢、占用內(nèi)存多的缺陷.Zhdanov等（2000）實(shí)現(xiàn)了三維擬線性反演，對(duì)日本Minamikayabe地區(qū)7個(gè)不同頻率（96、48、24、12、6、3和1.5）161個(gè)觀測(cè)點(diǎn)的MT數(shù)據(jù)進(jìn)行反演，其結(jié)果與Takasugi等（1992）的二維反演結(jié)果相近，但是該三維反演結(jié)果提供了更為詳細(xì)的地質(zhì)結(jié)構(gòu)，為解釋該地區(qū)高導(dǎo)有利含水區(qū)提供了依據(jù).

Zhdanov和Tartaras （2002）提出的局部擬線性（LQL）近似方法，能夠在不失計(jì)算精度的前提下，為多源電磁正、反演模擬及應(yīng)用節(jié)省大量時(shí)間.在此基礎(chǔ)上Zhdanov實(shí)現(xiàn)了三維電磁反演，該反演能夠準(zhǔn)確地反演出呈階梯狀逐漸變深的異常體.利用該反演對(duì)加拿大某海灣含鎳-鈷-鉻硫礦區(qū)的HEM數(shù)據(jù)進(jìn)行了反演，從反演結(jié)果來(lái)看，該區(qū)存在兩種不同的巖性，研究區(qū)中心存在低阻異常體并向南東、北西向延伸，該結(jié)果與Oldenburg等（1998）在該區(qū)的三維重力反演結(jié)果有很好的對(duì)應(yīng)關(guān)系.

Yoshioka和Zhdanov（2005）在LQL近似方法的基礎(chǔ)上實(shí)現(xiàn)了基于Cole-Cole模型的三維頻譜激電（SIP）非線性共軛梯度反演，對(duì)理論模型進(jìn)行了反演試算，但只在零頻電阻率上取得了良好的反演效果，其他參數(shù)反演效果不佳.

與其他電磁法不同，SIP方法的發(fā)射和觀測(cè)方式比較特殊，需要不斷更換發(fā)射與接收的位置，通過(guò)變換收發(fā)距離來(lái)達(dá)到測(cè)深的目的.由于換源頻繁，QL近似方法用于三維SIP正演需要不斷重復(fù)計(jì)算反射系數(shù)，并且每次換源都需要重新計(jì)算單元一次場(chǎng)，正演的速度受到了嚴(yán)重的制約.

本文首先利用并矢格林函數(shù)和一次場(chǎng)的空間對(duì)稱性，提出了一種適用于多源電磁法正演的快速Q(mào)L近似方法，該方法只需在發(fā)射源第一次發(fā)射時(shí)計(jì)算并存儲(chǔ)單元之間的并矢格林函數(shù)和一次場(chǎng)，在之后的每次換源正演計(jì)算中不必重新計(jì)算.另外，算法將并矢格林函數(shù)系數(shù)存儲(chǔ)成Toeplitz矩陣，利用Toeplitz矩陣的性質(zhì)減少不必要的計(jì)算.然后，基于文中提出的快速Q(mào)L正演方法，實(shí)現(xiàn)了Cole-Cole模型參數(shù)范圍約束下的三維SIP反演；最后，通過(guò)模型試算，驗(yàn)證了正演算法的正確性和反演算法的有效性.

2 正演

2.1 快速Q(mào)L近似理論

在非均勻介質(zhì)內(nèi)部電磁場(chǎng)總場(chǎng)E等于一次場(chǎng)與二次場(chǎng)之和:

E=En+Ea,

（1）

En為發(fā)射電流在背景區(qū)域內(nèi)部產(chǎn)生的一次場(chǎng)，Ea是由不均勻電導(dǎo)率的分布導(dǎo)致電流產(chǎn)生的二次場(chǎng).一次場(chǎng)為接地導(dǎo)線發(fā)射源在均勻半空間產(chǎn)生的一次場(chǎng).二次場(chǎng)由Zhdanov和Fang（1996a,1996b）提出的擬線性近似方法求得：

×En（r）dv,

（2）

（3）

（4）

在異常區(qū)域內(nèi)反射系數(shù)是緩慢變化的（Zhdanov和Fang,1996a）,因此近似地，可以把（4）式寫(xiě)為

（5）

然后可以得到

（6）

（7）

這里，J為剖分的塊數(shù)，ΔEα β γ為

（8）

*代表復(fù)數(shù)的共軛，ΔEα μ ν為ΔEα β γ的共軛.

求解（7）式，可得反射系數(shù)求解公式：

（9）

其中

（10）

（11）

其中，EB為Born近似解，En為一次場(chǎng)，Gn為格林系數(shù).

任意一單元的反射系數(shù)，可以由所有單元的異常電導(dǎo)率、一次場(chǎng)以及所有單元對(duì)當(dāng)前單元的格林系數(shù)進(jìn)行簡(jiǎn)單的線性運(yùn)算求得，其中格林系數(shù)的計(jì)算量龐大，耗費(fèi)時(shí)間較多.當(dāng)對(duì)異常體在x、y、z三個(gè)方向上均勻地進(jìn)行剖分，分別為nx,ny,nz,則對(duì)x方向，其格林系數(shù)矩陣可以寫(xiě)為

（12）

其中Γ（i,j）所代表的物理意義是第j個(gè)單元作為源項(xiàng)在第i個(gè)單元中心處產(chǎn)生的場(chǎng)，即第j個(gè)單元對(duì)第i個(gè)單元的格林系數(shù).如果第i個(gè)單元和第j個(gè)單元的相對(duì)位置不變，那么Γ（i,j）的值是相同的，例如Γ（1,2）=Γ（2,3），因此該矩陣滿足Toeplitz矩陣的結(jié)構(gòu).對(duì)于該矩陣，只需要存儲(chǔ)第一行和第一列的元素就可以表示出整個(gè)矩陣.

并矢格林函數(shù)代表的是源在接收點(diǎn)處產(chǎn)生的場(chǎng)，并且滿足互換定理（deLugao，1996），例如

（13）

因?yàn)閳?chǎng)具有方向性，因此上式兩端加上絕對(duì)值才成立.

y方向和x方向的格林系數(shù)矩陣具有相同的結(jié)構(gòu)，但z方向不滿足這樣的結(jié)構(gòu)，因此格林系數(shù)矩陣整體是一個(gè)二重Toeplitz矩陣.

（14）

對(duì)y方向，各個(gè)分量的對(duì)稱性質(zhì)為

（15）因此快速Q(mào)L近似方法只計(jì)算、存儲(chǔ)（12）式矩陣中第一行的元素，盡管比QL近似方法多占用一些內(nèi)存，但其計(jì)算量卻是QL近似方法的1/（nxny）.

雖然換源需要重新計(jì)算反射系數(shù)，但每次換源計(jì)算反射系數(shù)所需要的格林系數(shù)矩陣是相同的，因此只需要一次計(jì)算格林系數(shù)矩陣并將其存儲(chǔ)就可以完成所有換源反射系數(shù)的計(jì)算.

單元一次場(chǎng)關(guān)于發(fā)射源中心點(diǎn)空間對(duì)稱，并且只要發(fā)射源與接收點(diǎn)的相對(duì)位置相同，其場(chǎng)值是相等的，因此利用快速Q(mào)L近似方法模擬多源電磁法時(shí)，單元一次場(chǎng)的計(jì)算及存儲(chǔ)方式與格林函數(shù)相同，這里不再重述.

根據(jù)Pelton等（1978）的結(jié)論，巖礦石的復(fù)電阻率可由Cole-Cole模型描述：

（16）

式中，ρ（iω）為復(fù)電阻率，ρ0、η、c、τ分別為零頻電阻率、極化率、頻率相關(guān)系數(shù)、時(shí)間常數(shù)，ω=2πf，f為發(fā)射頻率.

將（16）式的復(fù)電阻率形式代入（2）式中，就可以實(shí)現(xiàn)考慮激電效應(yīng)的三維復(fù)電阻率正演.

2.2 正演模擬

為驗(yàn)證快速Q(mào)L近似方法的準(zhǔn)確性，給出三維正演模型如圖1所示，通過(guò)偶極-偶極觀測(cè)方式進(jìn)行觀測(cè)，虛線y=0為測(cè)線，發(fā)射電流為10A，通過(guò)兩個(gè)頻率0.1 Hz和64 Hz進(jìn)行計(jì)算.圖2給出積分方程法（李建平，2008）、快速Q(mào)L近似方法計(jì)算的二次場(chǎng)曲線對(duì)比結(jié)果.從圖中可以看出快速Q(mào)L近似方法計(jì)算精度較高.

表1給出積分方程法、QL近似方法及快速Q(mào)L近似方法計(jì)算時(shí)間及占用內(nèi)存對(duì)比結(jié)果，可以看到快速Q(mào)L近似方法的計(jì)算速度較快，而且隨著剖分塊數(shù)增加、接收點(diǎn)數(shù)和換源次數(shù)增多，快速Q(mào)L近似方法在計(jì)算速度上的優(yōu)勢(shì)越加明顯.

通過(guò)以上驗(yàn)證及對(duì)比，快速Q(mào)L近似方法計(jì)算準(zhǔn)確，其精度和占用內(nèi)存情況完全適用于反演，并且在模擬多源電磁正演時(shí)，計(jì)算速度較QL近似有了很大提高.

圖1 三維模型

圖2 積分方程法和快速擬線性近似法的三維SIP正演二次場(chǎng)結(jié)果對(duì)比，（a、b） 發(fā)射頻率為0.1 Hz; （c、d） 發(fā)射頻率為64 Hz

表1 積分方程法、QL近似法和快速Q(mào)L近似法計(jì)算時(shí)間和占用內(nèi)存對(duì)比表

3 參數(shù)范圍約束下的三維頻譜激電反演

3.1 反演理論

該反演分為兩步，第一步反演物質(zhì)性質(zhì)系數(shù)m（Zhdanov and Tartaras，2002），第二步反演Cole-Cole模型的四個(gè)參數(shù).

3.1.1 反演物質(zhì)性質(zhì)系數(shù)m

假設(shè)我們通過(guò)擬線性近似方法已經(jīng)得到了偶極-偶極裝置下的二次場(chǎng)，重寫(xiě)（2）式：

（17）

其中m為反映地下物質(zhì)性質(zhì)的參數(shù)，有

（18）

由于三維反演的多解性尤為嚴(yán)重，在這里我們采用光滑約束反演方法.建立的目標(biāo)函數(shù)如下：φ=（Wdd-WdGm）*（Wdd-WdGm）+λLm*RTRm,

（19）這里*表示復(fù)數(shù)的共軛，Wd為數(shù)據(jù)方差向量，λL為拉格朗日因子，R為光滑度矩陣,G=GE（En（r））.上式對(duì)m求導(dǎo)取極小得到反演方程：

（20）

其中Δd為觀測(cè)數(shù)據(jù)與正演理論數(shù)據(jù)的殘差向量，Δm為模型修正量.

擬合精度定義為

（21）

RMS為擬合精度，k為數(shù)據(jù)個(gè)數(shù).

求解上述反演問(wèn)題，本文采用了共軛梯度法，反演流程如下：

（1） 給定初始模型m0及拉格朗日因子λL；

（4）hpi=p·tdi；

（6） 更新模型，mi=mi+tpi·tdi，計(jì)算擬合差，達(dá)到給定精度則跳出循環(huán)，否則繼續(xù)；

（7）ri+1=ri+tpi·hpi；

（9）tdi+1=ri+1+βi·tdi；

（10） 令i=i+1,然后跳到（4）.

在反演過(guò)程中拉格朗日因子λL可以采取自適應(yīng)方式進(jìn)行選擇（劉云鶴和殷長(zhǎng)春，2013），具體如下：開(kāi)始給定一個(gè)相對(duì)較大的λL值進(jìn)行反演迭代，當(dāng)反演步長(zhǎng)小到一定程度時(shí)，將λL值變?yōu)樵瓉?lái)的十分之一，當(dāng)λL值或擬合差小于閥值時(shí)終止迭代.

3.1.2Cole-Cole模型參數(shù)的確定

在反演Cole-Cole模型參數(shù)時(shí)，本文進(jìn)行逐塊反演.根據(jù)擬線性近似理論：

λ（rj）En（rj）=G（m（r））,

（22）

由此可以近似計(jì)算出反射系數(shù)λ，然后通過(guò)（18）式計(jì)算每個(gè)單元的異常復(fù)電導(dǎo)率Δσ（ω）.

（23）

這里σ（ω）=σn+Δσ（ω），σn為背景電導(dǎo)率.

對(duì)于（23）式，通過(guò)四個(gè)頻率列出非線性方程組，利用共軛梯度迭代求解出Cole-Cole模型參數(shù).考慮到Cole-Cole模型參數(shù)均為正實(shí)數(shù)，且在一定范圍內(nèi)變化，在這里引入Kim等（1999）、Commer和Newman（2008）給出的約束函數(shù)：

（24）

其中x代表Cole-Cole模型參數(shù)，Δx為反演參數(shù)修正量，x0為初始模型.這樣就將模型x限定在（xmin,xmax）范圍內(nèi)，能夠有效改善反演的多解性.

3.2 理論模型反演

反演模型如圖1所示，反演區(qū)域如圖3所示，地面設(shè)置260個(gè)觀測(cè)點(diǎn)，將反演區(qū)域均勻剖分成1920塊，采用偶極-偶極觀測(cè)方式如圖3所示，在13條線上進(jìn)行觀測(cè)，每次發(fā)射的接收點(diǎn)數(shù)為12個(gè)，每條線換源22次.對(duì)該模型通過(guò)0.1 Hz、1 Hz、4 Hz及16 Hz四個(gè)頻率進(jìn)行正演模擬，得到的地面觀測(cè)總場(chǎng)作為觀測(cè)數(shù)據(jù)，加入1%的高斯噪聲進(jìn)行反演.通過(guò)反演試算，測(cè)試該反演的速度、占用內(nèi)存大小、穩(wěn)定性及對(duì)高阻體和低阻體的反演效果.

圖3 三維反演區(qū)域（a）及數(shù)據(jù)觀測(cè)方式（b）

反演物質(zhì)性質(zhì)系數(shù)m時(shí)，給定的初始模型為均勻半空間模型，m（r）=0，r=1,2,…,n；反演Cole-Cole模型參數(shù)時(shí)，給定初始模型為ρ0=100 Ωm、η=10-5、c=0.1、τ=10 s.由于Cole-Cole模型各參數(shù)變化范圍大致為η=0～0.98、τ=10-3～5×103s、c=0.1～0.6、ρ0=10-4～105Ωm，因此給定各參數(shù)范圍約束的最大（max）值分別為ρ0max=105Ωm、ηmax=0.98、cmax=0.6、τmax=5×103s；最?。╩in）值分別為ρ0min=10-4Ωm、ηmin=0、cmin=0.1、τmin=10-3s.

給定初始拉格朗日因子值為λL=10-2，隨著反演的進(jìn)行，λL隨著迭代步長(zhǎng)變小而逐漸減小，當(dāng)達(dá)到給定精度或步長(zhǎng)低于閥值時(shí)終止迭代.設(shè)定的精度為10-5，步長(zhǎng)閥值設(shè)定為10-4.

從圖4反演結(jié)果來(lái)看，該反演能夠在零頻電阻率、極化率及頻率相關(guān)系數(shù)上有較高的分辨率，時(shí)間常數(shù)的分辨率差；低阻異常體的反演效果優(yōu)于高阻異常體的反演效果；所有Cole-Cole模型參數(shù)的反演結(jié)果都在自身合理的變化范圍內(nèi)，因?yàn)榉秶s束函數(shù)強(qiáng)行剔除了不合理的解，是改善反演多解性的重要手段.

從靈敏性角度分析:零頻電阻率的靈敏度高于其他三個(gè)參數(shù)；極化率、頻率相關(guān)系數(shù)和時(shí)間常數(shù)的靈敏度受頻率影響，頻率越高，極化率的靈敏度越高，而頻率相關(guān)系數(shù)和時(shí)間常數(shù)的靈敏度越低；Cole-Cole參數(shù)反演靈敏性的失衡造成了時(shí)間常數(shù)的反演結(jié)果較差，無(wú)法準(zhǔn)確反演出異常體的時(shí)間常數(shù).

從圖5的收斂曲線中可以看出反演穩(wěn)定，只是共軛梯度迭代次數(shù)較多.表2給出反演m的時(shí)間及占用內(nèi)存情況；在反演Cole-Cole模型參數(shù)時(shí)，反演迅速，一般需要幾十秒至幾分鐘的時(shí)間.雖然共軛梯度法的迭代次數(shù)較多，但總體反演時(shí)間短，在三維反演領(lǐng)域內(nèi)，該反演在速度上具有較大優(yōu)勢(shì)，而且占用的內(nèi)存小，有利于實(shí)際生產(chǎn)應(yīng)用.

表2 反演m的時(shí)間及占用的內(nèi)存

本文所有正反演計(jì)算均在PC機(jī)上進(jìn)行，計(jì)算機(jī)內(nèi)核為AMD Athlon（tm） II 631 Quad-Core Processor 2.60 GHz, 2.99 GB 內(nèi)存.

圖4 反演結(jié)果

圖5 迭代擬合曲線

4 結(jié)論

本文提出了快速擬線性近似方法，為多源電磁法正、反演的深入研究奠定了良好的基礎(chǔ)；實(shí)現(xiàn)了Cole-Cole模型參數(shù)范圍約束的三維頻譜激電反演.取得以下幾條結(jié)論：

（1） 利用并矢格林函數(shù)及一次場(chǎng)的空間性質(zhì)可以有效減少擬線性近似方法的計(jì)算量，縮短正演的計(jì)算時(shí)間.掌握格林函數(shù)的物理實(shí)質(zhì)，挖掘格林函數(shù)更多的性質(zhì)并靈活運(yùn)用，有利于進(jìn)一步加快擬線性近似方法的計(jì)算速度.

（2） 頻譜激電方法的參數(shù)種類較多，各個(gè)參數(shù)變化范圍不一，這增加了三維頻譜激電反演的復(fù)雜性和多解性.Cole-Cole參數(shù)范圍約束函數(shù)能剔除不合實(shí)際的解，有效改善反演的多解性.Cole-Cole模型各參數(shù)的靈敏度有高有低，其中時(shí)間常數(shù)的靈敏度最低，因此該參數(shù)的反演結(jié)果不佳，但其他參數(shù)的反演效果較好.

（3） Cole-Cole參數(shù)范圍約束的三維頻譜激電反演方法對(duì)理論模型的反演效果良好，計(jì)算速度快、穩(wěn)定性高、占用內(nèi)存小，完全能夠在普通計(jì)算機(jī)上處理實(shí)際數(shù)據(jù).

致謝 感謝賁放、伍亮、孫博、李鶴、蘇曉波以及關(guān)振瑋、張镕哲、陳帥在作者編寫(xiě)論文過(guò)程中的幫助.

Born M. 1933. Optic. New York: Springer. Commer M, Newman G A. 2008. New advances in three-dimensional controlled-source electromagnetic inversion.Geophys.J.Int., 172（2）: 513-535.

de Lugao P P, Wannamaker P E. 1996. Calculating the two-dimensional magnetotelluric Jacobian in finite elements using reciprocity.Geophys.J.Int., 127（3）:806-810

Eaton P A. 1989. 3D electromagnetic inversion using integral equations.Geophys.Prosp., 37（4）: 407-426.

Kim H J, Song Y, Lee K H. 1999. Inequality constraint in least-squares inversion of geophysical data.Earth,PlanetsandSpace, 51（4）: 255-259.

Lee K H, Xie G Q. 1993. A new approach to imaging with low-frequency electromagnetic fields.Geophysics, 58（6）: 780-796.

Li J P. 2008. Research on electromagnetic modeling of 3D complex resistivity with topography [Ph. D. thesis] （in Chinese）. Changchun: Jilin University.

Liu Y H, Yin C C. 2013. 3D inversion for frequency-domain HEM data.ChineseJ.Geophys. （in Chinese）, 56（12）: 4278-4287. Madden T M, Mackie R L. 1989. Three-dimensional magnetotelluric modelling and inversion.Proc.IEEE, 77（2）: 318-332.

Nekut A G. 1994. Electromagnetic ray-trace tomography.Geophysics, 59（3）: 371-377. Newman G A, Alumbaugh D L. 1996. Three-dimensional electromagnetic modeling and inversion on massively parallel computers. Sandia Report SAND, 96-0852.

Oldenburg D W, Li Y, Farquharson C G, et al. 1998. Applications of geophysical inversions in mineral exploration.TheLeadingEdge, 17:461-465.

Oldenburg D W, Haber E, Shekhtman R. 2013. Three dimensional inversion of multisource time domain electromagnetic data.Geophysics, 78（1）: E47-E57.

Oristaglio M L, Wang T, Hohmann G, et al. 1993. Resistivity imaging of transient EM data by conjugate-gradient method.∥ 63rd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 347-350.

Pellerin L, Johnston J, Hohmann G. 1993. 3-D inversion of EM data.∥ 63rd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 360-363. Pelton W, Ward S, Hallof P, et al. 1978. Mineral discrimination and removal of inductive coupling with multrfrequency IP.Geophysics, 43（3）:588-609.

Smith J T, Booker J R. 1991. Rapid inversion of two- and three-dimensional magnetotelluric data.J.Geophys.Res., 96（B3）: 3905-3922.

Takasugi S, Keisaku T, Noriaki K, Shigeki M. 1992. High spatial resolution of the resistivity structure revealed by a dense network MT measurement—A case study in the Minamikayabe area, Hokkaido, Japan: J. Geomag.Geoelectr., 44, 289-308.

Torres-Verdin C, Habashy T M. 1995. An overview of the extended Born approximation as a nonlinear scattering approach, and its application to cross-well resistivity imaging.∥ Progress in Electromagnetic Research Symposium Proc., 323.

Wang B N, Li G F, Cao W W. 2003. On the extension of Einstein′s summation convention.JournalofUniversityofShanghaiforScienceandTechnology（in Chinese）, 25（1）: 5-7.

Weng A H, Liu Y H, Jia D Y, et al. 2012. Three-dimensional controlled source electromagnetic inversion using non-linear conjugate gradients.ChineseJ.Geophys. （in Chinese）, 55（10）: 3506-3515, doi: 10.6038/j.issn.0001-5733.2012.10.034.

Xie G, Lee K H. 1995. Nonlinear inversion of 3-D electromagnetic data.∥ Progress in Electromagnetic Research Symposium Proc., 323.

Xu K J. 2007. Study on 2. 5-D complex resistivity electromagnetic forward and inversion [Ph. D. thesis] （in Chinese）. Changchun: Jilin University.

Yoshioka K, Zhdanov M S. 2005. Three-dimensional nonlinear regularized inversion of the induced polarization data based on the Cole-Cole model.PhysicsoftheEarthandPlanetaryInteriors,

150（1-3）: 29-43. Zhang H. 2006. Research of complex resistivity 3-D electromagnetic forward and inversion [Ph. D. thesis] （in Chinese）. Changchun: Jilin University.Zhdanov M S, Fang S. 1996a. Quasi-linear approximation in 3-D electromagnetic modeling.Geophysics, 61（3）: 646-665.

Zhdanov M S, Fang S. 1996b. Three-dimensional quasi-linear electromagnetic inversion.RadioSci., 31（4）: 741-754.

Zhdanov M S, Fang S, Hursan G. 2000. Electromagnetic inversion using quasi-linear approximation.Geophysics, 65（5）: 1501-1513.

Zhdanov M S, Tartaras E. 2002. Three-dimensional inversion of multitransmitter electromagnetic data based on the localized quasi-linear approximation.Geophys.J.Int., 148（3）: 506-519.

附中文參考文獻(xiàn)

李建平. 2008. 帶地形的三維復(fù)電阻率電磁場(chǎng)正反演研究[博士論文]. 長(zhǎng)春: 吉林大學(xué).

劉云鶴, 殷長(zhǎng)春. 2013. 三維頻率域航空電磁反演研究. 地球物理學(xué)報(bào), 56（12）: 4278-4287.

王伯年, 李過(guò)房, 曹偉武. 2003. 愛(ài)因斯坦求和約定的推廣. 上海理工大學(xué)學(xué)報(bào), 25（1）: 5-7.

翁愛(ài)華, 劉云鶴, 賈定宇等. 2012. 地面可控源頻率測(cè)深三維非線性共軛梯度反演. 地球物理學(xué)報(bào), 55（10）: 3506-3515, doi: 10.6038/j.issn.0001-5733.2012.10.034.

徐凱軍. 2007. 2. 5維復(fù)電阻率電磁場(chǎng)正反演研究[博士論文]. 長(zhǎng)春: 吉林大學(xué).

張輝. 2006. 復(fù)電阻率三維電磁場(chǎng)正反演研究[博士論文]. 長(zhǎng)春: 吉林大學(xué).

（本文編輯 胡素芳）

Fast quasi-linear approximation and the three-dimensional spectrum of induced polarization inversion study

LIU Yong-Liang1, LI Tong-Lin1*, HU Ying-Cai2, GAO Wen1, ZHU Cheng1, CHAI Lun-Wei1

1CollegeofGeo-explorationSciencesandTechnology,JilinUniversity,Changchun130026,China2NuclearIndustryGeologicalResearchInstituteinBeijing,Beijing100029,China

The quasi-linear （QL） approximation method proposed by Zhdanov effectively tackles the problem of slow speed as well as high memory usage of integral equation method. However, as the location of transmitters of spectrum induced polarization method （SIP） needs to be changed in different modeling computation, we need to recalculate the green′s function and primary field each time, so the speed of three-dimensional SIP modeling based on QL is still slow. In this paper, we proposed a fast QL approximation method, which is based on the spatial character of primary field and Green′s tensor. After that, the Cole-Cole model was introduced to our modeling computation algorithm and then we inversed three-dimensional SIP theoretical data with conjugate gradient method and parameter constraints. The results of model trial shows that: on one hand, the fast QL approximation method has higher speed than QL approximation method and is more appropriate for three-dimensional SIP modeling. On the other hand, the inversion method in this paper can obtain good inversion results for SIP parameters such as direct current resistivity, intrinsic chargeability and so on. Normally, the computing time of single inversion iteration is about 0.2 seconds and it takes up about 140 MB memory when the inhomogeneity domain is divided into two thousand blocks in the computing example.

Fast; Quasi-linear; Cole-Cole model; Conjugate gradient; Scope of constraint functions

深部探測(cè)技術(shù)與實(shí)驗(yàn)研究專項(xiàng)SinoProbe-03-05和國(guó)家重大科學(xué)儀器設(shè)備開(kāi)發(fā)專項(xiàng)（2011YQ05006009）聯(lián)合資助.

劉永亮，男，1986年生，博士，主要從事電磁法理論及三維正反演研究.E-mall:liuyongliang198718@163.com

10.6038/cjg20151231.

10.6038/cjg20151231

P631

2015-05-16，2015-11-11收修定稿

*通訊作者 李桐林,男，1962年生，教授，主要從事電磁法理論與應(yīng)用教學(xué)和研究.E-mall:lilaoshizh@163.com

劉永亮， 李桐林， 胡英才等. 2015. 快速擬線性近似方法及三維頻譜激電反演研究.地球物理學(xué)報(bào)，58（12）:4709-4717,

Liu Y L, Li T L, Hu Y C, et al. 2015. Fast quasi-linear approximation and the three-dimensional spectrum of induced polarization inversion study.ChineseJ.Geophys. （in Chinese）,58（12）:4709-4717,doi:10.6038/cjg20151231.