張 澎， 余小東， 許 洋， 呂仁斌， 陸從德,c

(1.成都理工大學(xué) a.信息科技與技術(shù)學(xué)院,b.地球物理學(xué)院,c. 地球勘探與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室, 成都 610059;2. 四川省蜀通巖土工程公司, 成都 610081)

半航空時間域電磁數(shù)據(jù)一維自適應(yīng)正則化反演

張 澎1a， 余小東1b， 許 洋2， 呂仁斌1b， 陸從德1b,c

(1.成都理工大學(xué) a.信息科技與技術(shù)學(xué)院,b.地球物理學(xué)院,c. 地球勘探與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室, 成都 610059;2. 四川省蜀通巖土工程公司, 成都 610081)

半航空時間域電磁法集成了地面和航空電磁法各自的優(yōu)勢，不僅容易實(shí)現(xiàn)大發(fā)射磁矩，增大勘探深度，同時還可以提高勘探效率。但是由于半航空時間域電磁法在測量過程中，收、發(fā)距會發(fā)生變化，因此與地面和航空電磁法相比，其數(shù)據(jù)處理與反演解釋更加困難。以長導(dǎo)線電性源半航空電磁法為例，將自適應(yīng)正則化反演算法應(yīng)用于半航空時間域電磁數(shù)據(jù)反演中。首先給出最平緩模型約束條件下的半航空時間域電磁數(shù)據(jù)自適應(yīng)正則化反演算法，并引入并行技術(shù)實(shí)現(xiàn)了該反演算法的并行計(jì)算；然后對兩種三層地電模型和一種六層地電模型的理論電磁響應(yīng)加高斯白噪聲后的數(shù)據(jù)進(jìn)行反演計(jì)算；最后對其反演結(jié)果和并行計(jì)算的效率進(jìn)行了詳細(xì)地分析與討論。實(shí)驗(yàn)結(jié)果分析表明：自適應(yīng)正則化反演算法是一種有效的半航空時間域電磁數(shù)據(jù)反演解釋方法，該方法具有很好的收斂性和穩(wěn)定性；自適應(yīng)正則化反演算法并行計(jì)算能大大減少反演計(jì)算時間，提高反演效率。

半航空時間域電磁法； 一維反演； 自適應(yīng)正則化； 并行計(jì)算

0 引言

航空電磁法(Airborne Electromagnetic Method, AEM)具有速度快、成本低、可大面積勘探、能適應(yīng)復(fù)雜地形條件等優(yōu)點(diǎn)，已廣泛應(yīng)用于礦產(chǎn)資源勘查、地質(zhì)填圖、水資源勘查和環(huán)境監(jiān)測等領(lǐng)域[1-2]。然而，AEM系統(tǒng)是將發(fā)射和接收線圈安裝于飛機(jī)或者吊艙上。發(fā)射磁矩容易受飛機(jī)負(fù)載的限制，極大地影響了航空電磁法的勘探深度。為了兼顧勘探深度和效率，Elliott等[3]提出了半航空電磁法(Semi-Airborne Electromagnetic Method, SAEM)。半航空電磁法通常是在地面鋪設(shè)長達(dá)幾公里的長導(dǎo)線電性源或大回線磁性源作為發(fā)射系統(tǒng)，并將接收安置在直升機(jī)、無人機(jī)等飛行器上。SAEM集成了地面電磁法和航空電磁法各自的優(yōu)勢，它不僅具有AEM系統(tǒng)速度快、效率高、空間分辨率高等優(yōu)點(diǎn)，還具有地面電磁法發(fā)射磁矩大、信噪比高、勘探深度大等優(yōu)勢。

20世紀(jì)50年代，國外就開始了半航空電磁勘查技術(shù)的研究。但由于相關(guān)技術(shù)的限制，尚未得到商用。直到20世紀(jì)90年代，Elliott等[4]研制了使用大回線源作為發(fā)射源的半航空電磁系統(tǒng)(Fixed Loop Airborne Transient Electromagnetics Method, FLAITTEM)，該系統(tǒng)首先在澳大利亞進(jìn)行了有效性測試，并在巴布亞新幾內(nèi)亞成功找到了一個硫化物礦藏。Fugro公司1997年研發(fā)了一套半航空時間域電磁測量系統(tǒng)(TerraAir系統(tǒng))，該系統(tǒng)實(shí)現(xiàn)了全波形數(shù)據(jù)采集；Smith[5]利用該系統(tǒng)在加拿大進(jìn)行了一次試驗(yàn)，并對比分析了地面、航空和半航空三種電磁法的數(shù)據(jù)，結(jié)果表明半航空電磁數(shù)據(jù)的信噪比要明顯高于航空電磁數(shù)據(jù)；Mogi等[7]基于長偏移距瞬變電磁法(Long Offset Transient Electromagnetic Method, LOTEM)原理[6]，設(shè)計(jì)了一套長導(dǎo)線電性源地空系統(tǒng)裝置，即GREATEM系統(tǒng)。近年來，該系統(tǒng)在日本已經(jīng)進(jìn)行較多的試驗(yàn)，成功應(yīng)用于火山區(qū)域和沿海區(qū)域的電阻率結(jié)構(gòu)探測[8-10]。國內(nèi)半航空電磁系統(tǒng)探測的研究起步較晚，嵇艷鞠等[10]描述了吉林大學(xué)團(tuán)隊(duì)自主研制的無人飛艇長導(dǎo)線源時域地空電磁勘探系統(tǒng)及其試驗(yàn)結(jié)果，證明了地空電磁探測方法具有一定的可行性，但是無人飛艇容易受到風(fēng)的影響[11]。

由于半航空電磁法測量裝置的特殊性，發(fā)射源與接收線圈之間的距離(即收發(fā)距)在飛行測量過程中會發(fā)生變化，這就使得半航空電磁數(shù)據(jù)的處理與解釋相對于航空和地面電磁數(shù)據(jù)而言更加困難，為此，半航空電磁數(shù)據(jù)的處理與解釋方法是亟待解決的關(guān)鍵問題。李肅義等[12]針對基于飛艇的電性源半航空時間域系統(tǒng)電磁數(shù)據(jù)中的主要噪聲，提出了一種綜合小波消除法，結(jié)果表明這種綜合消除法對多種噪聲均有很好的去噪效果，是一種實(shí)用有效的半航空時間域電磁數(shù)據(jù)去噪方法；嵇艷菊等[13]利用神經(jīng)網(wǎng)絡(luò)技術(shù)，對電性源半航空電磁數(shù)據(jù)進(jìn)行視電阻率反演研究，用于解決電性源時域半航空電磁法視電阻率反演效率低的問題。

由于反演問題一般是不確定的，反演結(jié)果容易出現(xiàn)不穩(wěn)定現(xiàn)象。要獲得更加穩(wěn)定準(zhǔn)確的反演結(jié)果，不僅需要從儀器上提高觀測數(shù)據(jù)的精度，還需要在反演過程中施加約束條件。正則化反演方法是常用的約束反演方法之一，它通過加入先驗(yàn)約束條件來增強(qiáng)反演過程的穩(wěn)定性，減少反演結(jié)果的非唯一性，是電磁數(shù)據(jù)反演中一種非常有效的反演方法。在正則化反演方法中，正則化因子反演結(jié)果的影響很大，因此正則化因子選取方案成為正則化反演方法研究的主要內(nèi)容之一。OCCAM反演[14]是正則化反演方法中的一種，它采用線性搜索來求取每次迭代過程中的最佳正則化因子，該算法已經(jīng)在航空電磁數(shù)據(jù)反演中得到了應(yīng)用[15]。但是OCCAM反演每次迭代需要額外求解多次反演方程和正演計(jì)算，使得計(jì)算量過大，限制了其在電磁數(shù)據(jù)反演中的應(yīng)用。為了解決這一問題，陳小斌[16]針對大地電磁數(shù)據(jù)反演，提出了新的正則化因子自適應(yīng)選取解決方案，并實(shí)現(xiàn)了最平緩模型約束下的大地電磁數(shù)據(jù)一維連續(xù)介質(zhì)反演；毛立峰[17]將自適應(yīng)正則化反演算法(Adaptive Regularized Inversion Algorithm, ARIA)應(yīng)用于直升機(jī)航空瞬變電磁數(shù)據(jù)一維反演解釋中，結(jié)果表明該方法具有良好的穩(wěn)定性和可靠性。

我們將自適應(yīng)正則化反演算法，應(yīng)用于長導(dǎo)線電性源半航空時間域電磁(Semi-Airborne Time-Domain Electromagnetics Method, SATEM)數(shù)據(jù)反演中，實(shí)現(xiàn)了最平緩模型約束下的半航空時間域電磁數(shù)據(jù)自適應(yīng)正則化反演算法，并對幾種典型地電模型的正演模擬數(shù)據(jù)進(jìn)行反演計(jì)算。反演結(jié)果分析表明，ARIA方法應(yīng)用于SATEM數(shù)據(jù)反演中具有較好的收斂性和穩(wěn)定性，是一種有效的SATEM數(shù)據(jù)反演方法。同時，還引入并行技術(shù)，實(shí)現(xiàn)了SATEM數(shù)據(jù)并行反演算法，以減少反演計(jì)算時間，提高反演效率。

1 SATEM一維正演模擬

長導(dǎo)線電性源SATEM一維層狀介質(zhì)正演模型如圖1所示，測量系統(tǒng)采用2 km～3 km的長接地導(dǎo)線作為發(fā)射源，將接收系統(tǒng)固定在直升機(jī)或無人機(jī)上進(jìn)行測量。水平收發(fā)距為R，接收線圈距地面高度為z．第n層的電阻率為ρn，厚度為dn．當(dāng)向接地導(dǎo)線供電流I時，可以得到頻率域磁場的垂直分量Hz的表達(dá)式如式(1)所示[18]

圖1 半航空電磁法正演模型Fig.1 Forward model of SAEM

(1)

其中：y為觀測點(diǎn)的y坐標(biāo)；R為偏移距；rTE為反射系數(shù)；λ為積分變量。

我們首先采用線性數(shù)字濾波方法[19]和Gauss數(shù)值積分方法，分別計(jì)算式(1)中的外層和內(nèi)層積分，從而獲得頻率域的半航空電磁響應(yīng)，然后采用G-S變換[20]將頻率域的電磁響應(yīng)轉(zhuǎn)換到時間域，為后面的反演提供數(shù)據(jù)基礎(chǔ)。

2 自適應(yīng)正則化反演算法

2.1 正則化反演的基本原理

正則化反演目標(biāo)函數(shù)可表示為式(2)。

Φ(m)=Φ1(m)+λΦ2(m)→min

(2)

其中：Ф是總體目標(biāo)函數(shù)；Ф1和Ф2分別是觀測數(shù)據(jù)和模型約束的目標(biāo)函數(shù)；m為模型參數(shù)向量；λ是正則化因子。

觀測數(shù)據(jù)目標(biāo)函數(shù)Ф1可表示為

Φ1=ΔdTσdΔd

(3)

其中：σd為數(shù)據(jù)方差相關(guān)矩陣；Δd為模型理論響應(yīng)與觀測數(shù)據(jù)之差向量。

最小模型約束、最平緩模型約束、最光滑模型約束，是三種常見的模型約束目標(biāo)函數(shù)(即模型粗糙度函數(shù))構(gòu)建方法。我們采用模型參數(shù)導(dǎo)數(shù)的平方和最小的最平緩模型約束條件來構(gòu)建模型粗糙度函數(shù)，假設(shè)模型可由一連續(xù)函數(shù)m(s)表示，則模型約束目標(biāo)函數(shù)可表示為式(4)。

Φ2=∫[?m(s)σ(s)?m(s)]ds

(4)

其中，σ(s)是模型方差相關(guān)函數(shù)，通常令其σ(s)≡1，則式(4)可改寫為式(5)。

Φ2=∫[?m(s)]2ds

(5)

2.2 粗糙度核矩陣的構(gòu)建

對于式(5)中?m(s)，假設(shè)m(s)是N個模型參數(shù)求解域已知的某種插值函數(shù)，即

(6)

其中，C=(C1(s)、C2(s)、…、CN(s))T，那么Ф2可表示為式(7)。

Φ2=mT[∫Ω?C?CTds]m=mTRm

(7)

(8)

式(8)滿足Guass積分的條件，可采用Guass數(shù)值積分進(jìn)行快速計(jì)算。

2.3 自適應(yīng)正則化因子調(diào)整方案

正則化反演算法常用的正則化因子自適應(yīng)調(diào)節(jié)方案有2種：①M(fèi)D方案；②CMD方案[15]。

(9)

(10)

這種正則化因子自適應(yīng)調(diào)整方案的反演結(jié)果，能同時滿足總體目標(biāo)函數(shù)Ф和觀測數(shù)據(jù)目標(biāo)函數(shù)Ф1以及模型約束目標(biāo)函數(shù)Ф2的極小。式(9)正則化因子調(diào)整方案在反演開始之前，需要設(shè)置一個初始的正則化因子值，而式(10)這種方案不需要。這兩種正則化因子的自適應(yīng)調(diào)節(jié)方法，不需要額外的正演計(jì)算，相對于傳統(tǒng)正則化方法減少了反演計(jì)算。

2.4 模型參數(shù)修正量的求取

根據(jù)總體目標(biāo)函數(shù)(2)極小原則，反演方程可表示為式(11)。

(11)

其中：m0為當(dāng)前已知模型參數(shù)向量；Δm是模型參數(shù)的修正量；G0是雅可比矩陣。因此，求解式(11)可得模型參數(shù)修正量為

(12)

從而可得到新模型參數(shù)向量為式(13)。

m=m0+Δm

(13)

2.5 反演迭代終止條件

由式(12)和式(13)計(jì)算得到新模型參數(shù)向量后，便可計(jì)算新模型的正演電磁響應(yīng)。如何判斷新模型是否滿足反演結(jié)果的要求、反演迭代過程何時結(jié)束，這就需要設(shè)定相關(guān)的反演迭代終止條件。目前主要采用相對擬合誤差rms(式(14))來進(jìn)行判斷。

(14)

通常而言，當(dāng)rms到達(dá)預(yù)先設(shè)定的誤差期望值時，反演迭代結(jié)束。但實(shí)際情況由于各種噪聲的影響，這一條件可能難以得到滿足，因此需要設(shè)定多個反演終止條件，從而避免反演迭代陷入死循環(huán)。筆者除了設(shè)定預(yù)期的rms值作為終止條件外，還通過設(shè)定最大迭代次數(shù)等作為終止條件。

3 反演算例分析

為了驗(yàn)證自適應(yīng)正則化反演算法應(yīng)用于半航空時間域電磁數(shù)據(jù)反演的正確性和穩(wěn)定性，對兩種典型的三層地電模型和六層復(fù)雜地電模型的理論電磁響應(yīng)加10%的高斯白噪聲后的數(shù)據(jù)進(jìn)行反演計(jì)算，并對結(jié)果進(jìn)行分析和討論。理論響應(yīng)計(jì)算時，系統(tǒng)的參數(shù)統(tǒng)一設(shè)置為：長接地線源長度為1 km；發(fā)射電流為1 A；偏移距為200 m；飛行高度為100 m。

3.1 兩種三層地電模型的反演結(jié)果與分析

分別對三層H型和K型地電模型理論電磁響應(yīng)加噪聲的數(shù)據(jù)進(jìn)行反演計(jì)算。三層H型地電模型的參數(shù)設(shè)置為：電阻率分別為300 Ω·m、50 Ω·m、300 Ω·m；層厚度分別為：100 m、100 m、∞。三層K型地電模型的參數(shù)設(shè)置為：電阻率分別為50 Ω·m、300 Ω·m、50 Ω·m；層厚度分別為：100 m、100 m、∞。反演初始模型設(shè)置為100 Ω·m的均勻半空間模型，并劃分為15層，每層厚度設(shè)置為20 m。反演計(jì)算過程中，正則化因子初始值λ0=1.0，誤差期望值設(shè)置為1.0e-4，迭代次數(shù)最大為20次。反演結(jié)果分別如圖2和圖3所示，圖2、圖3中給出了迭代過程中早期、中期和晚期的反演結(jié)果模型。

從圖2和圖3兩種三層地電模型反演結(jié)果可以得出：反演迭代計(jì)算早期反演結(jié)果模型都與真實(shí)模型相差較大，三層地電模型中間層為低阻時最后的反演結(jié)果模型更接近真實(shí)模型，反演迭代10次就達(dá)到收斂終止條件；而中間層為高阻時，反演結(jié)果不理想，中間層電阻率和真實(shí)模型相差較大，并且反演迭代20次后，相對擬合誤差都沒達(dá)到預(yù)期設(shè)置的誤差值，這與時間域電磁法本身對高阻不敏感有關(guān)。通過兩種三層地電模型的反演結(jié)果分析，表明自適應(yīng)正則化反演算法，能有效地應(yīng)用于半航空時間域電磁數(shù)據(jù)的反演解釋。

圖2 H型地電模型反演結(jié)果Fig.2 Inversion results of H model

圖3 K型地電模型反演結(jié)果Fig.3 Inversion results of K model

3.2 六層地電模型的反演結(jié)果與分析

為了進(jìn)一步說明自適應(yīng)正則化反演算法的穩(wěn)定性和有效性，采用復(fù)雜的六層地電模型的理論電磁響應(yīng)加噪聲數(shù)據(jù)進(jìn)行反演計(jì)算，其模型參數(shù)設(shè)置為：電阻率分別為10 Ω·m、200 Ω·m、20 Ω·m、100 Ω·m、5 Ω·m、10 Ω·m；層厚度分別為：100 m、500 m、200 m、1 000 m、1 000 m、∞。該模擬數(shù)據(jù)的反演使用了較多的反演層數(shù)作為初始模型，即電阻率為 100 Ω·m的均勻半空間模型，并劃分為 50層，第一層厚度為 10 m，此后每層的厚度以等比因子1.037等比增加。反演計(jì)算過程中，正則化因子初始值λ0=1.0，誤差期望值設(shè)置為0.000 5，最大迭代次數(shù)20次。反演結(jié)果如圖4所示，分別給出了第1、5、10、15、17次反演迭代后的模型結(jié)果。另外，圖5至圖7還分別給出了該六層地電模型反演過程中自適應(yīng)正則化因子、目標(biāo)函數(shù)和相對擬合誤差隨迭代次數(shù)的變化曲線。

圖4 六層地電模型反演結(jié)果Fig.4 Inversion results of six layer model

圖5 正則化因子隨迭代次數(shù)變化曲線Fig.5 The changing curve of regularization factor

圖6 目標(biāo)函數(shù)隨迭代次數(shù)變化曲線Fig.6 The changing curve of objective function

圖7 相對擬合誤差隨迭代次數(shù)變化曲線Fig.7 The changing curve of RMS

從圖4可知，反演結(jié)果是逐漸擬合真實(shí)模型的，迭代5次后，第二層高阻層顯現(xiàn)出來，總體成三層模型結(jié)構(gòu)；迭代10次后，第三層低阻層也開始顯現(xiàn)，但是各層的電阻率和真實(shí)模型相差還較大；第15次迭代完成后，反演結(jié)果已經(jīng)非常接近真實(shí)模型，結(jié)果模型也十分光滑；但是此時反演計(jì)算并沒有終止，直到迭代17次后才終止，從最后第17次迭代的反演模型結(jié)果可以看出，高阻層的電阻率與真實(shí)模型之間有一定的差距，因此在反演過程中，并不是迭代次數(shù)越多越好，也不是相對誤差越小越好，要綜合考慮各方面因素，才能得出更接近真實(shí)模型的反演結(jié)果。該六層模型反演結(jié)果表明，自適應(yīng)正則化反演算法應(yīng)用于半航空時間域電磁數(shù)據(jù)的反演解釋具有較好的穩(wěn)定性。

從圖5至圖7中可以看出，在前5次迭代過程中，它們的值下降得十分迅速，說明反演過程收斂很快。大約迭代5次左右后，這些值的曲線變得相對平緩，說明反演結(jié)果逐漸趨于穩(wěn)定。

4 自適應(yīng)正則化反演算法的并行計(jì)算與分析

由于基于長導(dǎo)線電性源的半航空時間域電磁響應(yīng)的正演數(shù)值計(jì)算涉及大量的積分計(jì)算，反演過程中往往又包含多次正演計(jì)算，從而使得半航空時間域電磁數(shù)據(jù)的反演處理是非常耗時的。對自適應(yīng)正則化反演算法計(jì)算過程進(jìn)行分析，可以發(fā)現(xiàn)其中正演數(shù)值計(jì)算和雅可比矩陣計(jì)算是相對獨(dú)立的計(jì)算，滿足并行計(jì)算的伯恩斯坦準(zhǔn)則。因此，為了解決反演計(jì)算量的問題，引入了并行技術(shù)。

為了更好地實(shí)現(xiàn)SATEM數(shù)據(jù)自適應(yīng)正則化反演算法的并行計(jì)算，對三層H型地電模型數(shù)據(jù)反演的常規(guī)串行計(jì)算程序中，各部分的計(jì)算時間進(jìn)行了統(tǒng)計(jì)。表1給出了在華碩X42J PC機(jī)上該模型數(shù)據(jù)反演的串行計(jì)算所需時間。

該反演串行計(jì)算耗時測試中，反演總共進(jìn)行了8次迭代計(jì)算。從表1中可以看出，雅可比矩陣的計(jì)算時間占了總計(jì)算時間的87.05%，這是由于雅可比矩陣計(jì)算過程中涉及到多次正演計(jì)算，并且計(jì)算時間隨反演參數(shù)個數(shù)的增加而增多。因此，在編寫本文中的反演算法并行計(jì)算程序時，只需將正演計(jì)算和雅可比矩陣計(jì)算進(jìn)行并行處理即可。SATEM數(shù)據(jù)一維自適應(yīng)正則化反演算法并行計(jì)算的流程圖如圖8所示。

表1 三層H型模型數(shù)據(jù)反演串行計(jì)算耗時Tab. 1 Serial algorithm calculation time of H model data

圖8 SATEM數(shù)據(jù)一維自適應(yīng)正則化反演并行算法流程圖Fig.8 Flow chart of ARIA for SATEM data

為了驗(yàn)證該并行算法的計(jì)算效率，三層H型地電模型的電磁響應(yīng)數(shù)據(jù)在三臺不同配置的計(jì)算機(jī)上進(jìn)行并行反演計(jì)算。三臺計(jì)算機(jī)的配置如表2所示，并行計(jì)算的時間及效率對比結(jié)果如表3至表5所示，其中并行計(jì)算加速比等于串行算法計(jì)算耗時除以并行算法計(jì)算耗時，并行計(jì)算效率等于并行計(jì)算加速比除以并行計(jì)算的進(jìn)程數(shù)。

表2 計(jì)算機(jī)配置Tab. 2 Computer Configuration

并行計(jì)算結(jié)果統(tǒng)計(jì)分析表明,該反演的并行計(jì)算效率不僅與計(jì)算機(jī)CPU核心數(shù)有關(guān)，還與CPU的線程數(shù)有直接關(guān)系。由表3可以得出，串行計(jì)算耗時隨計(jì)算機(jī)配置的提高不會有較明顯地改善，但并行計(jì)算耗時隨線程數(shù)的增加而減少，也隨計(jì)算機(jī)配置的提高而減少。由表4和表5可知，雖然并行計(jì)算加速比隨線程數(shù)的增加而增大，并且隨計(jì)算機(jī)配置的提高也有所增大，但是在相同CPU的情況下，并行計(jì)算的效率隨進(jìn)程數(shù)的增加反而減小。這是由于在相同CPU的情況下，并行算法進(jìn)程數(shù)目越多，各進(jìn)程之間相互通信需要的時間就會增加。另外如果所有進(jìn)程中有一個進(jìn)程的耗時相對其他進(jìn)程較多時，對其他進(jìn)程的計(jì)算就會產(chǎn)生影響，從而使得并行計(jì)算的總時間成倍增加，降低了并行計(jì)算的效率。但是，從總體而言，并行反演算法的計(jì)算效率還是明顯高于傳統(tǒng)串行反演算法的計(jì)算效率，大大減少了反演的運(yùn)行時間。

表3 不同配置計(jì)算機(jī)并行計(jì)算耗時Tab. 3 Parallel algorithm calculation time of different computer

表4 不同配置計(jì)算機(jī)并行計(jì)算加速比Tab. 4 Parallel algorithm calculation speed-up ratio of different computer

表5 不同配置計(jì)算機(jī)并行計(jì)算效率Tab. 5 Parallel algorithm calculation efficiency of different computer

5 結(jié)論

將自適應(yīng)正則化反演方法應(yīng)用于半航空時間域電磁數(shù)據(jù)反演解釋中，并實(shí)現(xiàn)了該反演算法的并行計(jì)算。通過對三種地電模型模擬數(shù)據(jù)反演結(jié)果的分析，可得出以下結(jié)論：

1)自適應(yīng)正則化反演算法是一種有效的半航空時間域電磁數(shù)據(jù)反演解釋方法，該方法具有很好的收斂性和穩(wěn)定性。

2)自適應(yīng)正則化反演算法對低阻層有更好的反演效果，高阻層的反演結(jié)果與真實(shí)模型相差較大，這與時間域電磁法本身對高阻不敏感有關(guān)。

3)自適應(yīng)正則化反演算法并行計(jì)算能大大減少反演計(jì)算時間，提高反演解釋效率，但并行計(jì)算時并不是進(jìn)程數(shù)越多越好。

[1] 雷棟, 胡祥云, 張素芬. 航空電磁法的發(fā)展現(xiàn)狀[J]. 地質(zhì)找礦論叢, 2006, 21(1): 40-53. LEI D, HU X Y, ZHANG S F. Development status quo of airborne electromagnetic [J]. Contributions to Geology and Mineral Resources Research, 2006, 21(1): 40-53.(In Chinese)

[2] 余小東. 時間域直升機(jī)航空電磁法2.5維反演[D]. 成都：成都理工大學(xué), 2014. YU X D.Inversion of 2.5 Dimensionol tine domain helicopter electromagnetic method[D].Chengdu:Chengdu University of Techndogy,2014.(In Chinese)

[3] ELLIOTT P.New airborne electromagnetic method provides fast deep-target data turnaround [J].The Leading Edge,1996,15(4):309-310.

[4] ELLIOTT P. The Principles and Practice of FLAIRTEM [J]. Exploration Geophysics, 1998, 29(2): 58-60.

[5] SMITH R S . A comparison of data from airborne, semi-airborne, and ground [J]. Geophysics, 2001, 66(5): 1379-1385.

[6] STACK K. M. Exploration with deep transient electromagnetics[M]. Amsterdam: Elsevier Science, 1992.

[7] MOGI T, TANAKA Y, KUSUNOKI K, et al. Eevelopment of grounded electrical source airborne transient EM (GREATEM) [J]. Exploration Geophysics, 1998, 29(2): 61-64.

[8] MOGI T. Grounded electrical-source airborne transient electromagnetic (GREATEM) survey of Mount Bandai, north-eastern Japan [J]. Exploration Geophysics, 2009, 40(1):1-7.

[9] ITO H, KAIEDA H, MOGI T, et al. Grounded electrical-source airborne transient electromagnetics (GREATEM) survey of Aso Volcano, Japan [J]. Exploration Geophysics, 2014, 45(1): 43-48.

[10]ALLAH S. A, MOGI T, ITO H, et al. Three-dimensional resistivity modelling of grounded electrical-source airborne transient electromagnetic (GREATEM) survey data from the Nojima Fault, Awaji Island, south-east Japan [J]. Exploration Geophysics, 2014, 45(1): 49-61.

[11]嵇艷菊, 王遠(yuǎn), 徐江,等. 無人飛艇長導(dǎo)線源時域地空電磁勘探系統(tǒng)及其應(yīng)用[J]. 地球物理學(xué)報, 2013, 56(11): 3640-3650. JI Y J, WAMG Y, XU J, et al. Development and application of the grounded long wire source airborne electromagnetic exploration system based on an unmanned airship [J]. Chinese J. Geophys, 2013, 56(11): 3640-3650. (In Chinese)

[12]李肅義, 林君, 陽貴紅,等. 電性源時域地空電磁數(shù)據(jù)小波去噪方法研究[J]. 地球物理學(xué)報, 2013, 56(9): 3145-3152. LI X Y, LIN J, YANG G H, et al. Ground-airborne electromagnetic signals de-noising using a combined wavelet transform algorithm [J]. Chinese J. Geophys, 2013, 56(9): 3145-3152. (In Chinese)

[13]嵇艷鞠, 徐江, 吳瓊等. 基于神經(jīng)網(wǎng)絡(luò)電性源半航空視電阻率反演研究[J]. 電波科學(xué)學(xué)報, 2014, 29(5): 973-980. JI Y J, XU J, WU Q, et al. Apparent resistivity inversion of electrical source semi-airborne electromagnetic data based on neural network [J]. Chinese Journal of radio Science, 2014, 29(5): 973-980. (In Chinese)

[14]CONSTABLE S C,PARLER R L,CONSTABLE C G. Occam’s inversion: a practical algorithm for generating smooth models from electromagnetic sounding data [J]. Geophysics, 1987, 52(3): 289-300.

[15]VALLéE M.A.,SMITH R.S.Application of Occam's inversion to airborne time-domain electromagnetics [J].The Leading Edge,2009,28(3):284-287.

[16]陳小斌,趙國澤,湯吉,等.大地電磁自適應(yīng)正則化反演算法[J]. 地球物理學(xué)報, 2005, 48(4): 937-946. CHEN X B, ZHAO G Z, TANG J, et al. An adaptive regularized inversion algorithm for magnetotelluric data [J]. Chinese J. Geophys, 2005, 48(4): 937-946. (In Chinese)

[17]毛立峰, 王緒本, 陳斌. 直升機(jī)航空瞬變電磁自適應(yīng)正則化一維反演方法研究[J]. 地球物理學(xué)進(jìn)展, 2011, 26(1): 300-305. MAO L F，WANG X B，CHEN B．Study on an adaptive regularized 1 D inversion method of helicopter TEM data [J]．Progress in Geophys, 2011, 26(1): 300-305. (In Chinese)

[18]NABIGHIAN M. N. Electromagnetic Methods in Applied Geophysics [M]. Tulsa, Oklahoma: Society of Exploration Geophysicists, 1988.

[19]GUPTASARMA D, SINGH B. New digital linear filter for Hankel J0 and J1 transforms [J]. Geophysical Prospecting, 1997, 45(5): 745-762.

[20]KNIGHT J H, RAICHE A P. Transient electromagnetic calculations using the Gaver-Stehfest inverse Laplace transform method [J]. Geophysics, 1982, 47(1): 47-50.

[21]FOUNTAIN D.Airborne electromagnetic systems-50 years of devecopment[J].Exploration Geophysics,1998,299(1/2):1-11.

An adaptive regularized inversion of 1D semi-airborne time-domain electromagnetic data

ZHANG Peng1a, YU Xiaodong1b, XU Yang2, LV Renbin1b, LU Congde1b,c,

(1. Chengdu University of Technology a. College of Information Science & Technology,b. College of Geophysics, c. Key Laboratory of Earth Exploration and Information Technology of Ministry of Education, Sichuan, Chengdu 610059, China;2. Sichuan Shutong Geotechnical Engineering Company, Sichuan, Chengdu 610081, China)

Semi-airborne time-domain electromagnetic method (SATEM) combined with the advantages of ground and airborne electromagnetic method not only has large magnetic moment, but also has deep exploration depth. At the same time, it can improve the efficiency of exploration. However, in the measurement process of SATEM, the distance between transmitter and receiver may change. So, the SATEM data processing and inversion are harder than ground and airborne EM. Few of studies on the inversion method for SATEM data has been conducted currently. In this paper, based on long wire grounded source SAEM, the adaptive regularized inversion algorithm (ARIA) is applied to the inversion of SATEM data. Firstly, this paper describes the ARIA of SATEM data under the most gentle model constraints, and introduces the parallel technology to achieve the parallel computation of this inversion algorithm. Then, we get the inversion results of the theoretical response data with Gauss white noise of two three layers model and one six layers model. At last, we analyze the inversion results and the efficiency of parallel computation. The test results show that ARIA is an effective SATEM data inversion method, it has good convergence and stability. The parallel computation of ARIA can greatly reduce the inversion calculation time and improve the efficiency of inversion.

semi-airboune time domain electromagnetic method; 1D inversion; adaptive regularized inversion algorithm; parallel computation

2016-08-03 改回日期：2016-10-19

國土資源行業(yè)專項(xiàng)(201311037)

張澎(1975-)，男，博士，講師，研究方向?yàn)閿?shù)字信號處理與計(jì)算機(jī)應(yīng)用技術(shù)，E-mail：uzhangpeng@qq.com。

陸從德(1972-)，男，博士，副教授，研究方向?yàn)殡姶判盘柼幚砼c正反演解釋，E-mail: cdlu@163.com。

1001-1749(2017)01-0001-08

P 631.3

A

10.3969/j.issn.1001-1749.2017.01.01