時(shí)間域航空電磁2.5維非線性共軛梯度反演

強(qiáng)建科， 滿開(kāi)峰*， 龍劍波， 魯凱， ZHU Yue， 陳龍偉， 李俊營(yíng)， 毛先成

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院， 長(zhǎng)沙 410083 2 中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)驗(yàn)室， 長(zhǎng)沙 410083 3 Department of Geology and Geophysics, University of Utah, Salt Lake City, Utah, USA

?

時(shí)間域航空電磁2.5維非線性共軛梯度反演

強(qiáng)建科1,2， 滿開(kāi)峰1,2*， 龍劍波1,2， 魯凱1,2， ZHU Yue3， 陳龍偉1,2， 李俊營(yíng)1,2， 毛先成1,2

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院， 長(zhǎng)沙 410083 2 中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)驗(yàn)室， 長(zhǎng)沙 410083 3 Department of Geology and Geophysics, University of Utah, Salt Lake City, Utah, USA

對(duì)于時(shí)間域航空電磁法二維和三維反演來(lái)說(shuō)，最大的困難在于有效的算法和大的計(jì)算量需求. 本文利用非線性共軛梯度法實(shí)現(xiàn)了時(shí)間域航空電磁法2.5維反演方法，著重解決了迭代反演過(guò)程中靈敏度矩陣計(jì)算、最佳迭代步長(zhǎng)計(jì)算、初始模型選取等問(wèn)題．在正演計(jì)算中，我們采用有限元法求解拉式傅氏域中的電磁場(chǎng)偏微分方程，再通過(guò)逆拉氏和逆傅氏變換高精度數(shù)值算法得到時(shí)間域電磁響應(yīng)．在靈敏度矩陣計(jì)算中，采用了基于拉式傅氏雙變換的伴隨方程法，時(shí)間消耗只需計(jì)算兩次正演，從而節(jié)約了大量計(jì)算時(shí)間．對(duì)于最佳步長(zhǎng)計(jì)算，二次插值向后追蹤法能夠保證反演迭代的穩(wěn)定性．設(shè)計(jì)兩個(gè)理論模型，檢驗(yàn)反演算法的有效性，并討論了選擇不同初始模型對(duì)反演結(jié)果的影響．模型算例表明：非線性共軛梯度方法應(yīng)用于時(shí)間域航空電磁2.5維反演中穩(wěn)定可靠，反演結(jié)果能夠有效地反映地下真實(shí)電性結(jié)構(gòu)．當(dāng)選擇的初始模型電阻率值與真實(shí)背景電阻率值接近時(shí)，能得到較好的反演結(jié)果，當(dāng)初始模型電阻率遠(yuǎn)大于或遠(yuǎn)小于真實(shí)背景電阻率值時(shí)反演效果就會(huì)變差．

時(shí)間域航空電磁法； 2.5D瞬變電磁反演； 伴隨方程法； 非線性共軛梯度； 靈敏度矩陣

1 引言

航空電磁法(AEM)分為頻率域航空電磁(AFEM)和時(shí)間域航空電磁或航空瞬變電磁(ATEM)兩種，其中頻率域方法發(fā)展較早，研究成果相對(duì)較多(Alumbaugh and Alumbaugh，1995；Cox and Zhdanov，2008；Yin et al.，2014)．但AFEM效率較低、探測(cè)深度較小，而ATEM探測(cè)效率較高，能夠獲取較多的地質(zhì)信息．早在20世紀(jì)六七十年代，Wait(1969)、Nabighian (1970)、Ghosh(1972)、Fraser(1972)、Singh(1973)等研究了時(shí)間域電磁(TEM)電磁場(chǎng)的傳播特征．但由于TEM探測(cè)方法的理論相對(duì)復(fù)雜，在二維和三維的反演方面研究進(jìn)展緩慢，主要工作集中在二維、三維正演模擬(王華軍和羅延鐘，2003; Commer and Newman，2004；熊彬和羅延鐘，2006；李建慧，2011；周俊杰，2011；王宇航，2013；強(qiáng)建科等，2015；殷長(zhǎng)春等，2015)和一維反演與定性解釋方面(Yin and Fraser，2004；Yin and Hodges，2007； Viezzoli et al.，2008；陳斌等，2014)．然而地下真實(shí)情況往往是二維和三維結(jié)構(gòu)，一維ATEM的反演和解釋不能滿足實(shí)際工作需要，因此有必要繼續(xù)研究高維ATEM的反演問(wèn)題．

三維時(shí)間域電磁正反演最大的障礙是內(nèi)存和時(shí)間消耗巨大，目前仍然處于理論研究階段(Haber et al.，2007；Cox et al.，2012；Zaslavsky et al.，2012；Oldenburg et al.，2013)，而2.5維ATEM正反演相對(duì)于三維而言，對(duì)計(jì)算資源和時(shí)間的消耗要小許多，也是最有可能實(shí)現(xiàn)突破的．截至目前，對(duì)于航空瞬變電磁法2D、2.5D的反演研究相對(duì)較少．Farquharson(1995)提出了頻率域中基于伴隨場(chǎng)的靈敏度矩陣的計(jì)算方法，為高維電磁法反演中靈敏度矩陣的求解提供理論依據(jù)．熊彬等(2004)利用原始場(chǎng)和輔助場(chǎng)之間的點(diǎn)乘和面積積分求得電磁響應(yīng)值對(duì)地下空間電導(dǎo)率的偏導(dǎo)數(shù)矩陣，但其反演算法仍在研究中．Wilson(2006)實(shí)現(xiàn)了基于阻尼特征參數(shù)法的頻率域2.5維航空電磁數(shù)據(jù)反演．Guillemoteau(2012)利用玻恩近似和經(jīng)驗(yàn)判斷的方法進(jìn)行快速2D航空瞬變電磁法的反演研究．余小東(2014)分別使用阻尼特征參數(shù)反演法與阻尼最小二乘反演法實(shí)現(xiàn)了2.5D航空瞬變電磁反演算法，并對(duì)比研究了兩種方法的反演結(jié)果．

非線性共軛梯度法(Non-linear Conjugate Gradient，NLCG)不需要對(duì)最優(yōu)化問(wèn)題進(jìn)行線性化，在反演時(shí)間和內(nèi)存的耗費(fèi)上相對(duì)較低，非常適合解決高維電磁反演問(wèn)題．自Newman和Alumbaugh(2000)提出采用非線性共軛梯度方法解決電磁反演問(wèn)題以來(lái)，該方法已經(jīng)應(yīng)用在大地電磁法、可控源電磁法的反演解釋中(Rodi and Mackie，2001；Newman and Boggs，2004；Kelbert et al.，2008；Commer and Newman，2009；翁愛(ài)華等，2012)．本文將非線性共軛梯度方法應(yīng)用到2.5D時(shí)間域航空電磁反演中．在正演過(guò)程中，我們利用有限單元法求解基于拉氏傅氏域中異常電磁場(chǎng)的微分方程，再通過(guò)逆拉氏逆傅氏雙變換得到時(shí)間域的瞬變電磁響應(yīng)．反演過(guò)程中，采用伴隨方程法求解靈敏度矩陣，通過(guò)二次插值向后追蹤法求解最佳步長(zhǎng)．通過(guò)對(duì)理論模型數(shù)據(jù)進(jìn)行反演，驗(yàn)證了時(shí)間域航空電磁2.5維非線性共軛梯度反演的有效性．此外，討論了噪聲、不同初始模型對(duì)反演結(jié)果的影響．

2 正演理論

由Maxwell方程組出發(fā)，同時(shí)將總感應(yīng)場(chǎng)劃分為背景場(chǎng)和異常場(chǎng)：

e=eb+ea, h=hb+ha，

(1)

式中e和h分別代表時(shí)間域中電場(chǎng)和磁場(chǎng)，下標(biāo)a代表異常場(chǎng)，下標(biāo)b代表背景場(chǎng)．將時(shí)間域Maxwell方程組經(jīng)過(guò)拉氏變換得到相應(yīng)背景場(chǎng)的偏微分方程:

(2)

(3)

式中E和H分別代表拉氏域中的電場(chǎng)和磁場(chǎng)，ε代表介電常數(shù)，σ代表介質(zhì)電導(dǎo)率，μ代表磁導(dǎo)率，s代表拉氏域變量，J代表場(chǎng)源電流密度矢量．

由總場(chǎng)等于背景場(chǎng)和異常場(chǎng)之和可得到異常場(chǎng)的偏微分方程：

(4)

(5)

其中Ja=σaEb.對(duì)式(4)、(5)的偏微分方程進(jìn)行y方向的一維傅里葉變換便可得到拉氏傅氏域中的偏微分方程，并利用有限單元法求解得到異常場(chǎng)．本文中的有限單元法采用三角形網(wǎng)格剖分技術(shù)；對(duì)于最終形成的大型線性方程組，采用Pardiso_64位直接解法求解器進(jìn)行求解，該方法的優(yōu)點(diǎn)在于求解精度高，多源情況下能實(shí)現(xiàn)并行回代．

經(jīng)過(guò)以上有限元法求解得到的是拉氏傅氏域中的異常場(chǎng)，需要再經(jīng)過(guò)逆傅氏變換和逆拉氏變換才能得到時(shí)間域中的異常場(chǎng)值．在逆傅氏變換中，由于采用漢克爾變換的濾波系數(shù)過(guò)多，計(jì)算量過(guò)大，本文采用三次樣條插值技術(shù)，從而在保證計(jì)算精度的情況下大大減少計(jì)算量．在逆拉氏變換中，采用Gaver-Stehfest變換技術(shù)，其優(yōu)點(diǎn)是計(jì)算量小，純實(shí)數(shù)運(yùn)算(羅延鐘和昌彥君，2000)．最后將得到時(shí)間域中的異常場(chǎng)，與解析法計(jì)算得到的時(shí)間域背景場(chǎng)相加就得到了總感應(yīng)場(chǎng)．

為了驗(yàn)證該正演算法的有效性與準(zhǔn)確性，將計(jì)算結(jié)果與層狀模型的半解析解(羅延鐘等，2003) 進(jìn)行了比較．收發(fā)裝置采用水平共面裝置(如圖1)，飛行高度為30 m，收發(fā)距為10 m，發(fā)射波形為階躍波，層狀模型的第一層、第三層電導(dǎo)率均為0.01 S·m-1，第二層電導(dǎo)率為0.05 S·m-1，第一層、第二層的厚度均為100 m，圖2a所示為層狀模型得到的數(shù)值解與半解析解對(duì)比情況．從圖2b可以看出，最大相對(duì)誤差2.40%，平均相對(duì)誤差1.49%．由此可知，在時(shí)間為1 μs到10 ms的范圍內(nèi)數(shù)值解與半解析解吻合較好，精度較高．

圖1 模型示意圖Fig.1 Sketch of model

圖2 層狀模型的數(shù)值解與解析解的對(duì)比結(jié)果(a) 垂直磁場(chǎng)數(shù)值解與解析解對(duì)比； (b) 數(shù)值解與解析解的相對(duì)誤差.Fig.2 Comparison of numerical and analytic solutions on the layering model(a) Comparison of numerical and analytic solutions for vertical magnetic field; (b) Relative errors of numerical and analytic solutions.

3 反演理論

3.1 NLCG方法的最優(yōu)化問(wèn)題

反演問(wèn)題的目標(biāo)函數(shù)可以寫(xiě)成如下形式：

Ψ(m)= (d-F(m))TWd(d-F(m))

+λ(m-mpre)TWm(m-mpre)，

(6)

式中，d為N維觀測(cè)數(shù)據(jù)向量，m為M維的模型空間向量，而F為正演算子，Ψ表示目標(biāo)函數(shù)，Wm為模型加權(quán)矩陣，Wd為數(shù)據(jù)加權(quán)矩陣，λ為正則化因子，上標(biāo)T表示轉(zhuǎn)置.目標(biāo)函數(shù)的梯度為

(7)

式中，Am為靈敏度矩陣．目標(biāo)函數(shù)Ψ(m)的最小值通過(guò)NLCG進(jìn)行求解，每次反演迭代可寫(xiě)成以下形式：

mk+1=mk+αkpk, k=0,1,…,K-1，

(8)

式中，k為迭代次數(shù)，αk為迭代步長(zhǎng)，pk是迭代方向矢量．

正則化因子λ在反演中約束步長(zhǎng)的大小，每減小一次λ導(dǎo)致正則化部分在目標(biāo)函數(shù)的比重變小，模型修正空間變大，可以搜索較大的步長(zhǎng)．這里正則化因子λ的選擇方法是借鑒Zhdanov(2009)的自動(dòng)化選取λ值的方法．首先給定一個(gè)初始值λ0，計(jì)算λk=λ0qk-1，q>0，k=1，2，3，…，K代表反演的迭代次數(shù)．

3.2 靈敏度矩陣的分析

靈敏度是指電磁響應(yīng)對(duì)模型參數(shù)的一階導(dǎo)數(shù)，即構(gòu)成靈敏度矩陣A的元素，其計(jì)算方法對(duì)于梯度類(lèi)的最優(yōu)化反演算法有著相當(dāng)重要的作用，本文采用伴隨方程法求解垂直磁場(chǎng)的靈敏度．拉氏域中垂直磁場(chǎng)的靈敏度為(龍劍波，2014)

|r=r0=∫Vφj(r)·E(r)·E+(r)dV，

(9)

式中的E為人工源區(qū)域V內(nèi)的電場(chǎng)空間分布函數(shù)，而E+為伴隨單位磁偶源激發(fā)的伴隨電場(chǎng)函數(shù)．因?yàn)轶w積分的被積函數(shù)中φj只在對(duì)應(yīng)的第j個(gè)網(wǎng)格內(nèi)不為0，所以每一個(gè)靈敏度矩陣的元素只需在相應(yīng)的單元內(nèi)進(jìn)行積分即可．反演中只需計(jì)算一次靈敏度矩陣，靈敏度矩陣求解只需要計(jì)算2次正演，從而節(jié)約了反演計(jì)算時(shí)間．

公式(9)可以進(jìn)一步簡(jiǎn)化為:

·E+(x,-ky,z)dxdzdky，

(10)

該式就是空間域中走向方向y=0處，場(chǎng)值Hz對(duì)電導(dǎo)率的偏導(dǎo)數(shù)計(jì)算式．為研究靈敏度的分布規(guī)律，設(shè)計(jì)如下均勻半空間模型，地下電導(dǎo)率為0.01 S·m-1，發(fā)射線框位于地面，邊長(zhǎng)為20 m，電流強(qiáng)度為10 A，選取采樣時(shí)間點(diǎn)分別為t1= 0.0020 ms，t2=0.0056 ms，t3=0.0155 ms，t4=0.0431 ms，t5=7.1876 ms，t6=20.0000 ms．圖3給出了各個(gè)時(shí)間點(diǎn)場(chǎng)值Hz對(duì)空間每個(gè)單元電導(dǎo)率的靈敏度分布特征(圖中數(shù)據(jù)為靈敏度幅值的對(duì)數(shù))．在時(shí)間上，靈敏度的幅值從早期(t1)到晚期(t6)不斷衰減，這種變化與場(chǎng)值在時(shí)間上的衰減類(lèi)似；在空間上，靈敏度幅值則從場(chǎng)源附近向四周迅速衰減，在早期t1～t4，靈敏度從淺部到深部的衰減較快，而在較晚期時(shí)間點(diǎn)t5～t6，靈敏度衰減速率相對(duì)變慢，在200 m內(nèi)已沒(méi)有較明顯的“輪廓線”．這表明，相對(duì)于晚期，早中期時(shí)間點(diǎn)的靈敏度對(duì)介質(zhì)更敏感．

3.3 步長(zhǎng)搜索

在NLCG主循環(huán)迭代過(guò)程中，由于目標(biāo)函數(shù)Ψ(mk+αkpk)不僅僅是迭代步長(zhǎng)α的二次函數(shù)(至少高于二次)，通常需要用迭代法獲得最佳步長(zhǎng)αk,best，這樣會(huì)增加時(shí)間消耗.本文結(jié)合前人的研究成果(劉云鶴和殷長(zhǎng)春，2013；Newman and Alumbaugh，2000)，采用二次插值向后追蹤方法(Backtracking)計(jì)算最佳步長(zhǎng)．在Wolfe準(zhǔn)則中，目標(biāo)函數(shù)充分下降需要滿足兩個(gè)不等式

(11)

(12)

式(11)稱(chēng)為充分下降條件，式(12)為曲率條件，其中的常數(shù)c1、c2為0到1之間的正數(shù)．本文采用的步長(zhǎng)搜索過(guò)程中最關(guān)鍵的是求取NLCG最佳步長(zhǎng)αk,best：

1) 求得Ψ(mk,1)、g(mk,1)、p(mk,1)，并計(jì)算：

(13)

然后取αt=min[1.0,1.01×temp_α]，min表示最小值；并計(jì)算f1=Ψ(mk,0+αtpk,0).

2) 若αt滿足下降條件：l=f1- f0-gk,0pk,0αt≤0，其中f0=Ψ(mk,0)，則αk,best= αt，結(jié)束搜索，令

(14)

圖3 異常體正上方(Offset=0)垂直磁場(chǎng)的靈敏度(絕對(duì)值)對(duì)數(shù)值的空間、時(shí)間變化Fig.3 Temporal and spatial variations of logarithm values in sensitivity of vertical magnetic field (absolute value) above the anomaly body (offset=0)

3) 若αt不滿足下降條件，進(jìn)行向后追蹤，令

(15)

驗(yàn)證αk,best是否滿足不等式(11)，若滿足，則結(jié)束搜索；若不滿足，則令αt=αk,best，并重復(fù)循環(huán)，直到滿足設(shè)計(jì)要求．

4 理論數(shù)據(jù)反演

4.1 反演算法穩(wěn)定性

為了測(cè)試反演算法的穩(wěn)定性與有效性，對(duì)理論模型正演的數(shù)據(jù)加入5%的高斯噪聲再進(jìn)行反演．模型如圖4a所示，在背景電阻率為100 Ωm的均勻介質(zhì)中，嵌入一個(gè)大小為100 m×20 m、電阻率為50 Ωm、埋深為20 m的低阻異常體．采用水平共面裝置(如圖1)，飛行高度為30 m，收發(fā)距為10 m，發(fā)射波形為階躍波．測(cè)點(diǎn)為15個(gè)，測(cè)點(diǎn)間距為30 m，采樣時(shí)間范圍為1.2 μs～10 ms，采樣時(shí)間點(diǎn)為10個(gè)．正演測(cè)道曲線見(jiàn)圖4b，早期測(cè)道平緩，晚期道電磁響應(yīng)變化明顯．

由反演結(jié)果(圖4c)看出，低阻異常體得到很好的還原，低阻體的埋深、橫向和縱向尺寸以及位置都能得到較準(zhǔn)確的反映，電阻率值也大致與理論模型的電阻率值相一致，可以較為清晰地反映地下真實(shí)電性結(jié)構(gòu)．

圖4d 給出了迭代過(guò)程中均方相對(duì)誤差(RMS)的變化情況．從圖中可以看出，① 加入5%噪聲后對(duì)反演有影響，但不是很大，說(shuō)明算法比較穩(wěn)定；② 迭代開(kāi)始時(shí)，均方相對(duì)誤差較大，但下降較快，從第9次迭代往后，誤差下降速率減慢，有一段較長(zhǎng)的平緩趨勢(shì)，這反映了非線性共軛梯度反演收斂比較穩(wěn)定．收斂誤差約20%左右，相比較其他的物探方法來(lái)說(shuō)，這個(gè)誤差是很大的，主要原因在于瞬變電磁的原始數(shù)據(jù)是二次衰變感應(yīng)場(chǎng)，最大數(shù)據(jù)和最小數(shù)據(jù)相差5～6個(gè)數(shù)量級(jí)，這樣的數(shù)據(jù)體計(jì)算的擬合誤差一般較大，因此說(shuō)，評(píng)價(jià)瞬變電磁反演效果時(shí)，均方相對(duì)誤差僅作參考．

4.2 初始模型對(duì)反演結(jié)果的影響

為了測(cè)試反演算法的適應(yīng)性，我們?cè)O(shè)計(jì)了一個(gè)低阻與高阻混合的較復(fù)雜模型，目的是檢驗(yàn)反演算法對(duì)電性差異較大異常體的分辨能力，以及不同初始模型對(duì)反演結(jié)果的影響．復(fù)雜模型如圖5a所示，背景電阻率為100 Ωm，低阻異常體大小為100 m×30 m、電阻率為20 Ωm、埋深為20 m，高阻異常體大小為100 m×20 m、電阻率為300 Ωm、埋深為20 m．采用水平共面裝置(如圖1)，飛行高度為30 m，收發(fā)距為10 m，發(fā)射波形為階躍波，橫向測(cè)點(diǎn)為17個(gè)，點(diǎn)距為30 m，初始模型電阻率分別設(shè)為150 Ωm、100 Ωm、50 Ωm．

圖4 反演算法穩(wěn)定性測(cè)試(a) 低阻體模型； (b) 低阻體正演響應(yīng)測(cè)道曲線圖； (c) 低阻體模型的反演結(jié)果； (d) 反演迭代過(guò)程中RMS變化曲線.Fig.4 Tests of inversion algorithm stability(a) Low resistivity model; (b) TEM response of low resistivity model; (c) Inversion result of the low resistivity model; (d) RMS variation curve during iteration inversion.

圖5 不同初始模型反演結(jié)果(a) 復(fù)雜模型； (b) 初始模型為150 Ωm； (c) 初始模型為100 Ωm； (d) 初始模型為50 Ωm.Fig.5 Inversion results of different initial models(a) Complicated model; (b) Initial resistivity model of 150 Ωm; (c) Initial resistivity model of 100 Ωm; (d) Initial resistivity model of 50 Ωm.

在反演中，當(dāng)初始模型選擇為150 Ωm時(shí)(大于實(shí)際模型的背景電阻率100 Ωm)，反演結(jié)果見(jiàn)圖5b，低阻與高阻異?；镜玫竭€原，但低阻異常體的暈圈范圍較大，低阻電阻率值偏高于真實(shí)低阻異常體；高阻異常體的電阻率值偏低于真實(shí)值．

當(dāng)初始模型設(shè)為100 Ωm時(shí)(等于實(shí)際模型的背景電阻率100 Ωm)，反演結(jié)果見(jiàn)圖5c，低阻與高阻異?；镜玫竭€原，而且低阻異常體的暈圈范圍變小一點(diǎn)，低阻電阻率值和高阻電阻率值都比較接近于真實(shí)的模型電阻率值，反演結(jié)果較好．

當(dāng)初始模型選擇為50 Ωm時(shí)(小于實(shí)際模型的背景電阻率100 Ωm)，反演結(jié)果見(jiàn)圖5d，大致能夠反映出低阻與高阻異常的中心位置，但低阻異常體的暈圈范圍變大了許多，對(duì)低阻體的反演結(jié)果變差，對(duì)高阻體反演結(jié)果影響不大，而且淺部出現(xiàn)了局部小異常．

通過(guò)以上比較可以看出，當(dāng)選擇的初始模型電阻率值與真實(shí)背景電阻率值差別不大時(shí)，能夠得到穩(wěn)定的反演結(jié)果．

4.3 正演與反演時(shí)間消耗

由于時(shí)間域電磁法數(shù)值計(jì)算需要多次變換，導(dǎo)致計(jì)算時(shí)間較長(zhǎng)．本次研究在PC機(jī)上完成，計(jì)算機(jī)處理器為Intel Core i7-2600，主頻3.4 GHz，內(nèi)存16 GB．采用Fortran語(yǔ)言實(shí)現(xiàn)正演與反演算法程序．

模型網(wǎng)格剖分為72×79，拉氏域波數(shù)為12個(gè)，漢克爾變換系數(shù)采用250個(gè)，時(shí)間采樣點(diǎn)為10個(gè)，測(cè)點(diǎn)數(shù)為17個(gè).通過(guò)分析，程序各模塊的耗時(shí)差異較大：數(shù)據(jù)文件輸入與輸出和背景垂直磁場(chǎng)計(jì)算環(huán)節(jié)時(shí)間消耗很少，可忽略不計(jì)．在反演流程中，正演程序是關(guān)鍵，每調(diào)用一次正演程序需要約272 s，在形成靈敏度矩陣時(shí)需要約632 s，這部分比較耗時(shí)；進(jìn)入正常反演迭代后，平均每次迭代耗時(shí)約614 s，總體反演時(shí)間隨迭代次數(shù)線性增加，一般迭代約10次即可(表1)，大約需要1.97 h．由此看出，2.5D ATEM反演非常耗時(shí)，投入實(shí)用化前必須解決并行運(yùn)算問(wèn)題．

表1 ATEM反演算法各模塊耗時(shí)統(tǒng)計(jì)

5 實(shí)測(cè)數(shù)據(jù)反演

ATEM的實(shí)測(cè)數(shù)據(jù)來(lái)自于Cox等(2012)論文，是關(guān)于Reid-Mahaffy實(shí)驗(yàn)場(chǎng)采用MEGATEM II系統(tǒng)獲取的L40測(cè)線dB/dt垂直分量數(shù)據(jù)，由于本反演算法速度所限，對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行了抽稀處理，即抽取的測(cè)點(diǎn)盡量落在異常密集部位(圖6a)．MEGATEM II系統(tǒng)設(shè)置為90 Hz基頻、42%占空比的半正弦發(fā)射波形，系統(tǒng)記錄20道信號(hào)，包括沿x和z方向的dB/dt和感應(yīng)磁場(chǎng)．我們反演時(shí)僅采用了斷電后的dBz/dt數(shù)據(jù)，最后一道的時(shí)間是5.5 ms．沿測(cè)線方向網(wǎng)格間距30 m，深度方向網(wǎng)格間距10 m到50 m不等，網(wǎng)格大小為82×91，迭代20次，擬合誤差從87%下降到34.3%．從反演結(jié)果看，得到了三個(gè)明顯的低阻異常體，但異常中心收斂不集中，反演深度總體偏淺，基本與Cox等(2012)采用積分方程的方法反演的結(jié)果一致，證明本算法是正確的，但在計(jì)算時(shí)間和精度方面還需要繼續(xù)優(yōu)化改進(jìn)．

圖6 來(lái)自MEGATEM儀器的L40測(cè)線dB/dt的反演結(jié)果(a) dB/dt的垂直分量曲線圖； (b) 反演電阻率結(jié)果.Fig.6 Inversion results of dB/dt data from survey line L40 using instrument MEGATEM(a) Vertical component of observed data; (b) Inversion results of resistivity.

6 結(jié)論

本文采用NLCG方法，實(shí)現(xiàn)了2.5D ATEM反演算法．采用將背景場(chǎng)和異常場(chǎng)分離的辦法，用解析法計(jì)算背景場(chǎng)，用有限單元法計(jì)算異常場(chǎng)，從而消除了場(chǎng)源處的奇異性問(wèn)題．此外，采用快速、高精度的直接解法求解線性方程組，以及Gaver-Stehfest變換技術(shù)求解逆拉氏變換和三次樣條插值方法求解逆傅氏變換，從而保證正演算法的計(jì)算精度和計(jì)算效率.在處理靈敏度矩陣時(shí)，采用伴隨方程法進(jìn)行求解，從而提高反演計(jì)算效率.在計(jì)算最佳迭代步長(zhǎng)時(shí)，采用二次插值向后追蹤方法(Backtracking)，為反演迭代的穩(wěn)定性提供了保障．理論模型計(jì)算表明，反演結(jié)果能夠基本反映地下真實(shí)電性結(jié)構(gòu)．反演得到的低阻異常體較真實(shí)異常體的范圍略大，可能與正演算法的計(jì)算精度、反演過(guò)程中靈敏度的計(jì)算以及迭代步長(zhǎng)選取等多種因素有關(guān)，還需要進(jìn)一步優(yōu)化算法．同時(shí)，研究了選擇不同初始模型情況下得到的反演結(jié)果，當(dāng)選擇的初始模型電阻率值與真實(shí)背景電阻率值相吻合時(shí)，得到的反演結(jié)果最好．此外，目前反演存在計(jì)算耗時(shí)較大的問(wèn)題，還需要通過(guò)并行算法提高計(jì)算效率．

致謝 感謝外審專(zhuān)家給予了詳細(xì)的點(diǎn)評(píng)和中肯的修改意見(jiàn)．

Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences.GeophysicalProspecting, 43(8): 1021-1042.

Chen B, Mao L F, Liu G D. 2014. The estimated prospecting depth of CHTEM-I system by the method of diffusion electric field.ChineseJ.Geophys. (in Chinese), 57(1): 303-309, doi: 10.6038/cjg20140126. Commer M, Newman G A. 2004. A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanic sources.Geophysics, 69(5): 1192-1202. Commer M, Newman G A. 2009. Three-dimensional controlled-source electromagnetic and magnetotelluric joint inversion.GeophysicalJournalInternational, 178(3): 1305-1316.

Cox L H, Zhdanov M S. 2008. Advanced computational methods of rapid and rigorous 3-D inversion of airborne electromagnetic data.CommunicationsinComputationalPhysics, 3(1): 160-179. Cox L H, Wilson G A, Zhdanov M S. 2012. 3D inversion of airborne electromagnetic data.Geophysics, 77(4): WB59-WB69. Farquharson C G. 1995. Approximate sensitivities for the multi-dimensional electromagnetic inverse problem [Ph. D. thesis]. Vancouver: University of British Columbia. Fraser D C. 1972. A new multicoil aerial electromagnetic prospecting system.Geophysics, 37(3): 518-537. Ghosh M K. 1972. Interpretation of airborne EM measurements based on thin sheet models [Ph. D. thesis]. Toronto: University of Toronto. Guillemoteau J, Sailhac P, Behaegel M. 2012. Fast approximate 2D inversion of airborne TEM data: Born approximation and empirical approach.Geophysics, 77(4): WB89-WB97.

Haber E, Oldenburg D W, Shekhtman R. 2007. Inversion of time domain three-dimensional electromagnetic data.GeophysicalJournalInternational, 171(2): 550-564.

Kelbert A, Egbert G D, Schultz A. 2008. Non-linear conjugate gradient inversion for global EM induction: resolution studies.GeophysicalJournalInternational, 173(2): 365-381.

Li J H. 2011. 3D numerical simulation for transient electromagnetic field excited by large source loop based on vector finite element method [Ph. D. thesis] (in Chinese). Changsha: Central South University. Liu Y H, Yin C C. 2013. 3D inversion for frequency-domain HEM data.ChineseJ.Geophys. (in Chinese), 56(12): 4278-4287, doi: 10.6038/cjg20131230.

Long J B. 2014. Non-linear conjugate gradient inversion for 2.5D airborne transient electromagnetic data [Master′s thesis] (in Chinese). Changsha: Central South University.

Luo Y Z, Chang Y J. 2000. A rapid algorithm for G-S transform.ChineseJ.Geophys. (in Chinese), 43(5): 684-690, doi: 10.3321/j.issn:0001-5733.2000.05.012.

Luo Y Z, Zhang S Y, Wang W P. 2003. A research on one-dimension forward for aerial electromagnetic method in time domain.ChineseJ.Geophys. (in Chinese), 46(5): 719-724.

Nabighian M N. 1970. Quasi-static transient response of a conducting permeable sphere in a dipole field.Geophysics, 35(2): 303-309. Newman G A, Alumbaugh D L. 2000. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients.GeophysicalJournalInternational, 140(2): 410-424.Newman G A, Boggs P T. 2004. Solution accelerators for large-scale three-dimensional electromagnetic inverse problems.InverseProblems, 20(6): S151-S170. Oldenburg D W, Haber E, Shekhtman R. 2013. Three dimensional inversion of multisource time domain electromagnetic data.Geophysics, 78(1): E47-E57.

Qiang J K, Zhou J J, Man K F. 2015. Synthetic study of 2.5-D ATEM based on finite element method.GeophysicalandGeochemicalExploration(in Chinese), 39(5): 1059-1062.

Rodi W L, Mackie R L. 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion.Geophysics, 66(1): 174-187.

Singh S K. 1973. Electromagnetic transient response of a conducting sphere embedded in a conductive medium.Geophysics, 38(5): 864-893.

Viezzoli A, Christiansen A V, Auken E, et al. 2008. Quasi-3D modeling of airborne TEM data by spatially constrained inversion.Geophysics, 73(3): F105-F113. Wait J R. 1969. Electromagnetic induction in a solid conducting sphere enclosed by a thin conducting spherical shell.Geophysics, 34(5): 753-759. Wang H J, Luo Y Z. 2003. Algorithm of a 2.5 dimensional finite element method for transient electromagnetic with a central loop.ChineseJ.Geophys. (in Chinese), 46(6): 855-862.

Wang Y H. 2013. The research on HTEM 2.5D forward modeling and curve analysis [Master′s thesis] (in Chinese). Chengdu: Chengdu University of Technology.

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.

Wilson G A, Raiche A P, Sugeng F. 2006. 2.5D inversion of airborne electromagnetic date.ExplorationGeophysics, 37(4): 363-371.

Xiong B, Luo Y Z, Qiang J K. 2004. Methods for calculating sensitivities for 2.5D transient electromagnetic inversion (Ⅰ).ProgressinGeophysics(in Chinese), 19(3): 616-620.

Xiong B, Luo Y Z. 2006. Finite element modeling of 2.5-D TEM with block homogeneous conductivity.ChineseJ.Geophys. (in Chinese), 49(2): 590-597.

Yin C C, Fraser D C. 2004. The effect of the electrical anisotropy on the response of helicopter-borne frequency-domain electromagnetic systems.GeophysicalProspecting, 52(5): 399-416.

Yin C C, Hodges G. 2007. Simulated annealing for airborne EM inversion.Geophysics, 72(4): F189-F195.

Yin C C, Huang X, Liu Y H, et al. 2014. Footprint for frequency-domain airborne electromagnetic systems.Geophysics, 79(6): E243-E254.

Yin C C, Zhang B, Liu Y H, et al. 2015. 2.5-D forward modeling of the time-domain airborne EM system in areas with topographic relief.ChineseJ.Geophys. (in Chinese), 58(4): 1411-1424, doi:

10.6038/cjg20150427. Yu X D. 2014. Inversion of 2.5 dimensional time domain helicopter airborne electromagnetic method [Master′s thesis] (in Chinese). Chengdu: Chengdu University of Technology. Zaslavsky M, Druskin V, Abubakar A, et al. 2012. Large-scale Gauss-Newton inversion of transient controlled-source electromagnetic data using the model reduction approach.∥ SEG Technical Program Expanded Abstracts 2012. SEG: 1-6.

Zhdanov M S. 2009. Geophysical Electromagnetic Theory and Methods. New York: Elsevier.

Zhou J J. 2011. Research on airborne transient electromagnetic 2.5-D forward modeling [Master′s thesis] (in Chinese). Changsha: Central South University.

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

陳斌, 毛立峰, 劉光鼎. 2014. 用擴(kuò)散電場(chǎng)法估算CHTEM-I系統(tǒng)的探測(cè)深度. 地球物理學(xué)報(bào), 57(1): 303-309, doi: 10.6038/cjg20140126.

李建慧. 2011. 基于矢量有限單元法的大回線源瞬變電磁法三維數(shù)值模擬[博士論文]. 長(zhǎng)沙: 中南大學(xué).

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

龍劍波. 2014. 2.5維航空瞬變電磁數(shù)據(jù)的非線性共軛梯度反演[碩士論文]. 長(zhǎng)沙: 中南大學(xué).

羅延鐘, 昌彥君. 2000. G-S變換的快速算法. 地球物理學(xué)報(bào), 43(5): 684-690, doi: 10.3321/j.issn:0001-5733.2000.05.012.

羅延鐘, 張勝業(yè), 王衛(wèi)平. 2003. 時(shí)間域航空電磁法一維正演研究. 地球物理學(xué)報(bào), 46(5): 719-724.

強(qiáng)建科, 周俊杰, 滿開(kāi)峰. 2015. 時(shí)間域航空電磁法2.5維有限元模擬. 物探與化探, 39(5): 1059-1062.

王華軍, 羅延鐘. 2003. 中心回線瞬變電磁法2.5維有限單元算法. 地球物理學(xué)報(bào), 46(6): 855-862.

王宇航. 2013. 吊艙式時(shí)間域直升機(jī)航空電磁2.5維正演及響應(yīng)曲線分析[碩士論文]. 成都: 成都理工大學(xué).

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

熊彬, 羅延鐘, 強(qiáng)建科. 2004. 瞬變電磁2.5維反演中靈敏度矩陣計(jì)算方法(Ⅰ). 地球物理學(xué)進(jìn)展, 19(3): 616-620.

熊彬, 羅延鐘. 2006. 電導(dǎo)率分塊均勻的瞬變電磁2.5維有限元數(shù)值模擬. 地球物理學(xué)報(bào), 49(2): 590-597.

殷長(zhǎng)春, 張博, 劉云鶴等. 2015. 2.5維起伏地表?xiàng)l件下時(shí)間域航空電磁正演模擬. 地球物理學(xué)報(bào), 58(4): 1411-1424, doi: 10.6038/cjg20150427.

余小東. 2014. 時(shí)間域直升機(jī)航空電磁法2.5維反演[碩士論文]. 成都: 成都理工大學(xué).

周俊杰. 2011. 航空瞬變電磁法2.5維正演模擬研究[碩士論文]. 長(zhǎng)沙: 中南大學(xué).

(本文編輯 何燕)

2.5D inversion of time domain airborne electromagnetic data using nonlinear conjugate gradients

QIANG Jian-Ke1,2, MAN Kai-Feng1,2*, LONG Jian-Bo1,2, LU Kai1,2, ZHU Yue3, CHEN Long-Wei1,2, LI Jun-Ying1,2, MAO Xian-Cheng1,2

1SchoolofGeosciencesandInfo-physicsofCentralSouthUniversity,Changsha410083,China2KeyLaboratoryofMetallogenicPredictionofNonferrousMetalsofMinistryofEducation,CentralSouthUniversity,Changsha410083,China3DepartmentofGeologyandGeophysics,UniversityofUtah,SaltLakeCity,Utah,USA

Inversion of time domain airborne electromagnetic (AEM) data are known for the difficulties in large computational requirement and effective algorithms, especially for two and three-dimensional problems. We have developed a 2.5 dimensional (2.5D) inverse algorithm for the time domain AEM data using the nonlinear conjugate gradient method with improved accuracy and efficiency. This paper focuses on solving the computation of the sensitivity matrix, the optimal step length and the initial model selection in this algorithm. In the forward modeling, we employ the finite element method (FEM) to solve the Maxwell′s equations in the Laplace and Fourier domains. The time domain responses are then obtained by the high-precision inverse Laplace and Fourier transforms. The sensitivity matrix is calculated by using the adjoint equation method with the Laplace and Fourier transforms, which requires only two forward modeling per iteration and reduces the time cost significantly. The backtracking method in the optimal step length computation ensures the stability of this iterative inverse algorithm. Then, we present two model studies and discuss the effects of different initial models. The synthetic studies demonstrate our inversion algorithm is stable and can yield reliable results, which reflect the underground electrical structure reasonably. It also turns out that the inversion result is good if the initial model is close to the true background. The inversion model becomes worse if the initial model is several times larger or smaller than the true values of the background resistivity.

Time domain airborne electromagnetic method; 2.5D transient electromagnetic inversion; Adjoint equation method; Nonlinear conjugate gradient method; Sensitivity matrix

10.6038/cjg20161229.

國(guó)家自然科學(xué)基金項(xiàng)目(41472301，41174104)；中南大學(xué)“創(chuàng)新驅(qū)動(dòng)計(jì)劃”項(xiàng)目(2015CX008)資助.

強(qiáng)建科，男，1967年生，副教授，從事地球物理電磁法正、反演研究.E-mail：qiangjianke@163.com

*通訊作者 滿開(kāi)峰，男，1990年生，碩士研究生，從事地球物理電磁法研究.E-mail：mankaifeng@csu.edu.cn

10.6038/cjg20161229

P631

2015-12-24，2016-11-17收修定稿

強(qiáng)建科， 滿開(kāi)峰， 龍劍波等. 2016. 時(shí)間域航空電磁2.5維非線性共軛梯度反演. 地球物理學(xué)報(bào)，59(12):4701-4709,

Qiang J K, Man K F, Long J B, et al. 2016. 2.5D inversion of time domain airborne electromagnetic data using nonlinear conjugate gradients.ChineseJ.Geophys. (in Chinese),59(12):4701-4709,doi:10.6038/cjg20161229.