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

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

?

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

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

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

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

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

1 引言

隨著找礦難度的加大以及電磁法儀器和計算機技術(shù)的快速發(fā)展，人們對電磁法解決地質(zhì)問題的期待已越來越高，傳統(tǒng)的一維、二維解釋模型已不能滿足對復雜地質(zhì)模型的刻畫，三維電磁法正、反演技術(shù)已提到日程，對此，國內(nèi)外研究學者作了大量研究工作（如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; 翁愛華等,2012等）.但是，三維電磁法內(nèi)存消耗大、計算速度慢，目前還很難投入到實際生產(chǎn)應(yīng)用當中.Zhdanov和Fang（1996a,1996b）提出擬線性（QL）近似法，通過對比積分方程法、QL近似方法和Born近似方法對不同模型的計算結(jié)果， 驗證了QL近似方法的計算精度，有效解決了積分方程計算速度慢、占用內(nèi)存多的缺陷.Zhdanov等（2000）實現(xiàn)了三維擬線性反演，對日本Minamikayabe地區(qū)7個不同頻率（96、48、24、12、6、3和1.5）161個觀測點的MT數(shù)據(jù)進行反演，其結(jié)果與Takasugi等（1992）的二維反演結(jié)果相近，但是該三維反演結(jié)果提供了更為詳細的地質(zhì)結(jié)構(gòu)，為解釋該地區(qū)高導有利含水區(qū)提供了依據(jù).

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

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

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

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

2 正演

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

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

E=En+Ea,

（1）

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

×En（r）dv,

（2）

（3）

（4）

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

（5）

然后可以得到

（6）

（7）

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

（8）

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

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

（9）

其中

（10）

（11）

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

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

（12）

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

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

（13）

因為場具有方向性，因此上式兩端加上絕對值才成立.

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

（14）

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

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

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

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

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

（16）

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

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

2.2 正演模擬

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

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

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

圖1 三維模型

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

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

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

3.1 反演理論

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

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

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

（17）

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

（18）

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

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

（20）

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

擬合精度定義為

（21）

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

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

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

（4）hpi=p·tdi；

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

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

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

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

在反演過程中拉格朗日因子λL可以采取自適應(yīng)方式進行選擇（劉云鶴和殷長春，2013），具體如下：開始給定一個相對較大的λL值進行反演迭代，當反演步長小到一定程度時，將λL值變?yōu)樵瓉淼氖种唬敠薒值或擬合差小于閥值時終止迭代.

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

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

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

（22）

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

（23）

這里σ（ω）=σn+Δσ（ω），σn為背景電導率.

對于（23）式，通過四個頻率列出非線性方程組，利用共軛梯度迭代求解出Cole-Cole模型參數(shù).考慮到Cole-Cole模型參數(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個觀測點，將反演區(qū)域均勻剖分成1920塊，采用偶極-偶極觀測方式如圖3所示，在13條線上進行觀測，每次發(fā)射的接收點數(shù)為12個，每條線換源22次.對該模型通過0.1 Hz、1 Hz、4 Hz及16 Hz四個頻率進行正演模擬，得到的地面觀測總場作為觀測數(shù)據(jù)，加入1%的高斯噪聲進行反演.通過反演試算，測試該反演的速度、占用內(nèi)存大小、穩(wěn)定性及對高阻體和低阻體的反演效果.

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

反演物質(zhì)性質(zhì)系數(shù)m時，給定的初始模型為均勻半空間模型，m（r）=0，r=1,2,…,n；反演Cole-Cole模型參數(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，隨著反演的進行，λL隨著迭代步長變小而逐漸減小，當達到給定精度或步長低于閥值時終止迭代.設(shè)定的精度為10-5，步長閥值設(shè)定為10-4.

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

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

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

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

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

圖4 反演結(jié)果

圖5 迭代擬合曲線

4 結(jié)論

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

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

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

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

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

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.

附中文參考文獻

李建平. 2008. 帶地形的三維復電阻率電磁場正反演研究[博士論文]. 長春: 吉林大學.

劉云鶴, 殷長春. 2013. 三維頻率域航空電磁反演研究. 地球物理學報, 56（12）: 4278-4287.

王伯年, 李過房, 曹偉武. 2003. 愛因斯坦求和約定的推廣. 上海理工大學學報, 25（1）: 5-7.

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

徐凱軍. 2007. 2. 5維復電阻率電磁場正反演研究[博士論文]. 長春: 吉林大學.

張輝. 2006. 復電阻率三維電磁場正反演研究[博士論文]. 長春: 吉林大學.

（本文編輯 胡素芳）

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

深部探測技術(shù)與實驗研究專項SinoProbe-03-05和國家重大科學儀器設(shè)備開發(fā)專項（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)用教學和研究.E-mall:lilaoshizh@163.com

劉永亮， 李桐林， 胡英才等. 2015. 快速擬線性近似方法及三維頻譜激電反演研究.地球物理學報，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.