基于瞬變電磁矩變換的快速三維反演方法

饒麗婷，武欣，吳超，張曉娟，方廣有

1 中國科學院電磁輻射與探測技術(shù)重點實驗室，北京 1001902 中國科學院大學，北京 100049

?

基于瞬變電磁矩變換的快速三維反演方法

饒麗婷1,2，武欣1,2，吳超1,2，張曉娟1，方廣有1

1 中國科學院電磁輻射與探測技術(shù)重點實驗室，北京 1001902 中國科學院大學，北京 100049

瞬變電磁法的嚴格三維反演計算復雜、占用資源多，在普通計算機上難以實現(xiàn).本文引入瞬變電磁矩變換的概念，提出一種快速三維反演方法.該方法基于阻性限制(resistive limit)特性，建立包含異常體的三維大地的一階矩響應正演算法，根據(jù)不同約束條件，選擇優(yōu)化的最速下降法實現(xiàn)瞬變電磁快速三維反演.文中通過含異常體的三維大地正演一階矩與仿真數(shù)據(jù)一階矩的對比，驗證了快速三維正演算法的有效性，之后在不同約束條件下，利用優(yōu)化的最速下降法實現(xiàn)了對含噪聲的仿真瞬變電磁數(shù)據(jù)的快速三維反演.結(jié)果表明，該方法能夠在普通計算機上短時間內(nèi)較為準確地反演出地下異常體的體積和位置，在瞬變電磁數(shù)據(jù)的實時解釋工作中具有良好的應用前景.

瞬變電磁法；瞬變電磁矩變換；阻性限制；快速三維反演

1 引言

瞬變電磁法是一種建立在電磁感應原理上的時間域人工源電磁探測方法(Nabighian，1988)，廣泛應用于礦產(chǎn)資源勘探、地下水探查、地質(zhì)調(diào)查與地質(zhì)填圖等領(lǐng)域(Fitterman and Stewart，1986；Meju et al.，2002；He et al.，2012；嵇艷鞠等，2013；Xue et al.,2013;殷長春等，2015).瞬變電磁反演建立在正演算法基礎(chǔ)上，由于高維正演算法的復雜性,高維反演問題尚未妥善解決(薛國強等，2008；Yin et al.,2015).Cox等(2012)采用積分方程法實現(xiàn)了航空瞬變電磁的正反演，Oldenburg等(2013)利用有限元法實現(xiàn)了瞬變電磁的三維正反演，然而這些嚴格三維反演方法受限于復雜的三維正演算法，數(shù)據(jù)量巨大，占用資源過多，在普通計算機上幾乎無法運行，故以其為基礎(chǔ)的嚴格三維反演運算速度緩慢，不能真正投入實際應用.

為了使三維反演易于實現(xiàn)，必須盡量壓縮數(shù)據(jù)量并簡化正演算法，因此本文引入瞬變電磁矩變換.瞬變電磁矩變換是Smith和Lee (2002)在阻性限制和感性限制的基礎(chǔ)上擴展而來，其中，阻性限制定義為當頻率趨近于零時頻域瞬變磁場斜率的幅度值(Grant and West,1965)，該限制等于時域中理想負階躍響應(磁場)的面積積分即瞬變電磁一階矩.

Smith (2000)利用阻性限制對航空瞬變電磁數(shù)據(jù)成像，發(fā)現(xiàn)圖像中包含off-time觀測方式無法探測的地電特征；Reid和Macnae (2002)利用阻性限制對數(shù)據(jù)量龐大的航空電磁勘探數(shù)據(jù)建模，簡化了耗時巨大的航空數(shù)據(jù)解釋處理工作.電磁場在傳播過程中，當滿足阻性限制條件時，其已充分穿透目標體，目標體內(nèi)部感應作用可忽略.故可將目標體網(wǎng)格化，通過幾何耦合因子與時間常數(shù)乘積計算各目標體微元的響應，對所有微元響應進行線性疊加可得到目標體的總響應(King,1997).

由以上可知，利用瞬變電磁一階矩變換不僅可以大幅壓縮處理數(shù)據(jù)量，而且能夠簡化目標體響應計算，因此本文基于瞬變電磁矩變換，提出一種快速三維反演方法.首先，以瞬變電磁感性限制和阻性限制概念為基礎(chǔ)，引入瞬變電磁矩變換定義；然后以矩形回線源為例，推導了點導電球體以及均勻半空間的一階矩響應，并在此基礎(chǔ)上建立了含異常體大地的一階矩響應的近似三維正演算法；之后，采用含松弛因子的優(yōu)化最速下降法作為迭代算法，并結(jié)合約束條件，闡述快速三維反演的原理；最后，在不同約束條件下，對含噪聲的仿真瞬變電磁數(shù)據(jù)實施了快速三維反演，結(jié)果表明，該方法能夠在普通計算機條件下短時間內(nèi)較為準確地反演出地下異常體的體積和位置，本文的研究成果可應用于瞬變電磁數(shù)據(jù)的實時解釋工作中.

2 瞬變電磁矩變換

對于理想負階躍響應，可由一系列指數(shù)函數(shù)求和表示(Stolz and Macnae,1998)：

(1)

其中Ai和τi分別表示第i個指數(shù)函數(shù)的幅度和時間常數(shù).在頻域，負階躍響應表示為

(2)

當頻率趨于無窮大時，負階躍響應達到感性限制IL(Inductive Limit):

(3)

感性限制時，感應電流僅存在于目標體表面以抵抗一次場法向分量，此時，電流沒有穿透到目標體內(nèi)部，感性限制僅與發(fā)射接收的幾何信息有關(guān)，與目標體的電導率無關(guān)，目標體感性限制可表示如下：

(4)

上式中，Gk為幾何耦合因子，僅與發(fā)射接收的幾何信息有關(guān).

當頻域負階躍響應的斜率趨近于零頻時，達到阻性限制RL(Resistive Limit):

(5)

阻性限制時，電磁場已經(jīng)充分穿透目標體，當ω→0，磁場變化近乎為零，目標體內(nèi)部感應作用可忽略.此時，可將目標體網(wǎng)格化，通過幾何耦合因子與時間常數(shù)的乘積計算單個微元響應，將所有微元響應線性疊加得到目標體響應，表示如下：

(6)

上式中，τk為時間常數(shù)，與微元的電導率有關(guān).

基于阻性限制和感應限制的概念，Smith和Lee (2002)提出了瞬變電磁矩變換，定義如下:

(7)

即脈沖響應與時間加權(quán)的積分，其中g(shù)(t)表示脈沖響應，n表示時間加權(quán)的階數(shù).脈沖響應可由負階躍響應表示如下：

(8)

將矩變換采用磁場時間導數(shù)定義，可表示為

(9)

當n=0時，零階矩等于感性限制，此時零階矩響應與導體的電導率無關(guān)，僅僅與發(fā)射接收的幾何信息有關(guān)；當n=1時，一階矩等于阻性限制，此時目標體一階矩響應可由其所有微元響應線性疊加，每個微元響應分解成幾何耦合因子與時間常數(shù)的乘積，時間常數(shù)與導體的電導率相關(guān)；當n取更高階時，高階矩更加強調(diào)晚期信號的加權(quán).目前，還沒有明確的物理意義，從定義來看，高階矩適合探測更深部的異常體，然而對晚期信號加權(quán)也會放大噪聲，并且目標體高階矩響應計算更加復雜.綜上，本文將采用瞬變電磁一階矩，首先建立含異常體的大地的三維正演理論.

3 近似三維正演理論

3.1 點導電球體的一階矩

圖1 導電球體幾何示意圖Fig.1 The geometric configuration of the conducting sphere

Grant和West(1965)、Kaufman(1978)等先后推導了導電球體的全時域磁場響應，對時域磁場作一階矩變換后，導電球體響應的表達式更加簡化(Smith and Lee,2001)，如圖1所示，經(jīng)過矩形回線源激勵后，自由空間中點導電球體的一階矩為

(10)

在自由空間中，半徑為a和電導率為σ的圓球的時間常數(shù)定義為

(11)

將(10)式中時間常數(shù)外的部分作為幾何耦合因子，表示為

(12)

根據(jù)畢奧薩伐爾定律，一次場B0可表示為

(13)

(14)

(15)

沿x方向線電流源積分得到一次場的y分量和z分量，表達式分別為

(16)

(17)

沿回線各個方向的線電流定積分并求和，可得到總一次場B0.

3.2 均勻半空間的一階矩

設一個矩形發(fā)射線框置于均勻半空間表面，如圖2所示，直角坐標系原點位于線框中心，均勻半空間任意接收點處二次場的垂直磁場表達式為(Schaa and Fullagar,2012)：

圖2 矩形回線源裝置示意圖Fig.2 The geometry of a rectangular loop on the surface of half space

(18)

上式中，G(x,y,t)為單位電流元激勵的垂直磁場響應，

x1=XE-XR，x2=XW-XR，y1=YN-YR，y2=YS-YR.G(x,y,t)表達式為

(19)

(20)

(21)

上式中，當t→0時，中括號部分的表達式等于零，同時，尾部的積分可進一步化簡，因此，沿y方向的線電流激勵形成的一階矩響應表達式為

(22)

上式表示從端點(x,y1)到(x,y2)定積分.在矩形線框源激勵下，均勻半空間的總一階矩等于四條邊框方向的一階矩之和.

3.3 三維正演理論

(23)

微元的幾何耦合因子可通過等體積導電球體的幾何耦合因子近似，表示如下：

(24)

微元的時間常數(shù)τk同樣可由相同電導率和體積的導電球體時間常數(shù)近似，可得

(25)

其中，L為立方體微元的邊長.

4 三維快速反演原理

4.1 優(yōu)化最速下降法

TEM測量數(shù)據(jù)中包含背景響應和異常體響應，為了簡化反演問題，首先從N個測量數(shù)據(jù)一階矩中剔除背景一階矩得到異常體響應的參考一階矩d=(d1,d2,…,dN)T，令實測數(shù)據(jù)一階矩的誤差為q=(q1,q2,…,qN)T，將可能存在異常的區(qū)域劃分為K個立方體微元，微元時間常數(shù)為τ=(τ1,τ2,…,τK)T，令異常區(qū)域的理論一階矩為c=(c1,c2,…,cN)T.根據(jù)三維近似正演模型，第n個接收點處異常區(qū)域的理論一階矩cn可以表示為

(26)

上式中，Gnk為第n個接收點處第k個微元的幾何耦合因子.由于數(shù)據(jù)誤差的存在，根據(jù)加權(quán)最小二乘法的原理，定義異常區(qū)域的理論一階矩與參考一階矩的擬合差目標函數(shù)為

(27)

對于這種線性反演問題，選擇最速下降法進行迭代優(yōu)化.在迭代過程中，令τi代表第i次迭代起始時最優(yōu)時間常數(shù)矢量，經(jīng)過迭代后，第i+1次起始時最優(yōu)化時間常數(shù)矢量τi+1=τi+δτi，δτi為時間常數(shù)的擾動參數(shù)矢量.在反演迭代過程中，最速下降法以負梯度方向作為下降方向，因此，δτi可表示為

(28)

下面推導每次迭代中步長αi的大小，第i+1次迭代目標函數(shù)的表達式如下：

(29)

(30)

將式(30)代入式(29)，這樣，(29)式可寫為

(31)

其中，

(32)

(33)

(34)

為了使L(τi+1)取得最小值，有

(35)

將式(31)和式(33)代入到式(35)，求導可得

(36)

由上式可求得

(37)

最速下降法雖然簡單，計算快，無矩陣求逆，但是收斂速度慢，尤其是當越接近最優(yōu)解時，收斂速度進一步放慢.最速下降法收斂性慢的根本原因是最優(yōu)化步長的選擇，而不是負梯度下降方向(Raydan and Svaiter,2002).為了加快收斂速度，本文引入松弛因子ri對步長做調(diào)整，每次迭代中線性搜索當前最優(yōu)松弛因子，此時δτi為

(38)

4.2 約束條件

由于三維反演問題的欠定性，直接反演一般無法得到合理結(jié)果，為了促使反演結(jié)果貼近真實地下結(jié)構(gòu)，本文給出兩種約束條件：

(1) 電導率加權(quán)：起始時間常數(shù)為零，根據(jù)視電導率深度圖(conductivity-depth-imaging,CDI)對相應位置處微元賦予加權(quán)值，反演過程中，所有微元的時間常數(shù)限制為正值.

加權(quán)值由CDI決定，分布于0和1之間，第k微元相應的加權(quán)值為

(39)

其中，σk為第k微元相應的視電導率值，σmax為整個CDI中最大的視電導率.電導率加權(quán)使CDI中大電導率區(qū)域占主導作用，減弱小電導率區(qū)域的幾何耦合因子影響.需要注意的是，電導率加權(quán)值依賴CDI，如果CDI存在虛假異常，可能反演不出合理結(jié)果.

(2) 深度加權(quán)：起始時間常數(shù)為零，對不同深度的微元賦予加權(quán)值，反演過程中，所有微元的時間常數(shù)限制為正值.

由于淺層幾何耦合因子Gk較大，為了減弱淺層耦合因子的影響并促使反演朝著較深層變化(Oldenburg and Li,2005)，根據(jù)經(jīng)驗函數(shù)計算深度加權(quán)值，函數(shù)定義如下：

(40)

上式中，zk為第k微元的深度，s0為斜率因子，d0為參考深度.s0影響加權(quán)值的增長速度，斜率因子越小，加權(quán)值隨著深度變化增長越慢；設定d0相當于引入了表面層，在參考深度以上的加權(quán)值為零.不同d0和s0對應的加權(quán)值變化曲線如圖3所示.

圖3 不同斜率因子s0和參考深度d0的加權(quán)值變化曲線Fig.3 Depth weights for various s0 and d0factors

實際反演中，斜率因子s0和參考深度d0的選擇不是固定的.如果異常體分布在較深處，則應選擇較小的斜率因子，參考深度d0的引入是為了減弱發(fā)射線圈下方淺層的較強幾何耦合因子的影響，d0可通過CDI中低電導率逐漸變化到高電導率的分界位置決定.經(jīng)驗表明，CDI中電導率最高處對應的加權(quán)值為0.5時，是斜率因子的最優(yōu)值.

4.3 快速三維反演流程圖

圖4給出了快速三維反演的總體流程圖.

圖4 快速三維反演總體流程圖Fig.4 The flowchart of rapid 3D inversion procedure

5 三維正演驗證與反演實現(xiàn)

5.1 瞬變電磁一階矩驗證

從瞬變電磁矩變換的定義可以看出，一階矩的積分時間范圍是從零到無窮大，而實測數(shù)據(jù)是有限時間內(nèi).為了得到實測數(shù)據(jù)的一階矩，需補齊實測數(shù)據(jù)時間外積分，因此，實測數(shù)據(jù)一階矩可通過下式求解：

(41)

其中，t1和tn分別表示實測數(shù)據(jù)的起始時間窗和截止時間窗，σ1和σn分別表示t1和tn時間處的視電導率，頭部表示從0到t1時間內(nèi)磁場積分，中部表示實測數(shù)據(jù)的積分，尾部表示從tn到∞時間內(nèi)的磁場積分.頭部通過(22)式均勻半空間的完整一階矩減去(21)式中從t1到∞的數(shù)值積分求得；尾部直接對(21)式從tn到∞數(shù)值積分求得；中部通過對實測數(shù)據(jù)B(t)三次樣條插值得到擬合函數(shù)，然后對擬合函數(shù)在t1和tn時間內(nèi)數(shù)值積分.

圖5給出了電導率為1 mS·m-1和200 mS·m-1的均勻半空間上垂直磁場一階矩的計算結(jié)果，圖中發(fā)射源為矩形線框回線，中心位于(0E,0N)，其中，E，N分別代表東西向、南北向，線框大小為250E×250N，發(fā)射電流為1A，圖中測線在0N方位從-500E到500E，測線上接收點間距為50 m，仿真數(shù)據(jù)的時間窗從0.001 ms到100 ms.其中理論一階矩采用(21)式計算，理論中部是利用(21)式和(22)式計算出時間窗范圍內(nèi)的磁場積分，仿真中部即對仿真數(shù)據(jù)三次樣條插值后求定積分，拼接一階矩是仿真中部加上理論計算的頭尾一階矩，即(41)式的計算方法.這里，仿真數(shù)據(jù)是采用Gaver-Stehfest變換對任意形狀發(fā)射回線頻域響應求時間域響應得到(齊有政等，2013).

由圖5可以看出，理論一階矩與拼接一階矩完全吻合，理論中部與仿真中部同樣也在每個數(shù)據(jù)點上顯著地一致.圖5驗證了(21)式中均勻半空間的理論一階矩的正確性，同時也說明(41)式對實測數(shù)據(jù)進行拼接頭部和尾部得到完整一階矩的可行性.

5.2 近似三維正演算法驗證

EMIT Maxwell軟件的Marco算法模塊基于三維積分方程，可計算層狀大地中含多個棱柱體異常目標的瞬變電磁響應(Xiong and Tripp,1995).本文利用該算法模塊的結(jié)果驗證文中所提三維正演算法的有效性.

根據(jù)Marco算法模塊的特性，建立驗證近似三維正演算法的計算模型，該模型三維示意圖如圖6所示，在電導率為1 mS·m-1的均勻大地背景中放置電導率為1 S·m-1的高導異常體，該異常體大小為

圖5 均勻大地上垂直磁場的一階矩.(a)電導率為1 mS·m-1；(b)電導率為200 mS·m-1Fig.5 First order vertical TEM moment of a half space.(a) Conductivity of 1 mS·m-1;(b) Conductivity of 200 mS·m-1

圖6 仿真計算模型的三維幾何示意圖Fig.6 3D geometry of model for simulated calculation

圖7 仿真模型俯視圖與發(fā)射源位置及測線布置Fig.7 Top view of the simulated model and position of the transmitting source and the survey line arrangement

600E×600N×300Z，上表面中心的坐標為(-100E,0N,-300Z)，其中，Z代表深度向，矩形線框發(fā)射源的中心位于(0E,0N,0Z)，線框大小為450N×450E.圖7是發(fā)射線框和測線分布平面圖，如圖所示，在南北方位上從-500N到500N分布11條測線，測線間距為100 m，每條測線的走向從-500E到500E，均勻分布有21個接收點.

利用Marco算法模塊仿真磁場垂直分量的數(shù)據(jù)，發(fā)射電流波形為雙極性方波，發(fā)射電流1A，共30個接收時間窗口，時間窗口范圍從0.1 ms到53 ms.

仿真計算完成后，將時間域磁場數(shù)據(jù)通過5.1節(jié)中的方法轉(zhuǎn)換為一階矩，這里是正演驗證，計算頭部和尾部的積分時，直接采用背景介質(zhì)的設定電導率.由于異常體處于傳導性環(huán)境中，無法直接采用自由空間中時間常數(shù)公式計算其時間常數(shù)，因此，本文采用經(jīng)驗的時間常數(shù)分析方法，該方法通過對磁場晚期信號進行指數(shù)函數(shù)擬合，得出異常體的時間常數(shù)估值為1.6 ms.將平板剖分成108000個邊長為10 m的立方體微元，經(jīng)過正演計算后，在各條測線上垂直分量一階矩曲線如圖8所示，由于對稱性，只顯示了0N到-500N測線上的結(jié)果，可以看出，Marco仿真數(shù)據(jù)一階矩與近似三維正演的理論一階矩在每條測線上都有明顯的一致性，同時，RMS也顯示了仿真數(shù)據(jù)一階矩與理論一階矩之間只有較小的誤差.

5.3 快速三維反演實現(xiàn)

5.3.1 預處理工作

采用5.2節(jié)中的計算模型進行反演，參考圖4的反演流程圖，首先進行預處理工作，設定異常區(qū)域范圍為-500E到500E，-500N到500N，-1500Z到-100Z，將異常區(qū)域劃分成89600個邊長為25 m立方體單元，根據(jù)發(fā)射接收及微元的幾何參數(shù)計算幾何耦合因子矩陣G，矩陣大小為231×89600，需要注意的是，該幾何耦合因子矩陣與測量數(shù)據(jù)無關(guān)，因此，該矩陣可在測量之前完成計算，在反演中作為一個加載數(shù)據(jù)項.

對Marco仿真磁場數(shù)據(jù)形成CDI，圖9是在0N方位上經(jīng)過插值后的CDI切面圖，如圖所示，高電導率區(qū)域頂部區(qū)域一定程度反應了異常體的位置，然而高電導率區(qū)域的底部位置遠遠超出了實際異常體，該插值CDI可以用于計算電導率加權(quán)值，而且在實地測量數(shù)據(jù)反演中，CDI可作為參考，驗證反演結(jié)果的可靠性；從中心接收點的CDI中讀取t1和tn時間處的視電導率σ1和σn，利用(41)式計算出仿真數(shù)據(jù)一階矩，將仿真數(shù)據(jù)一階矩中加入5%高斯白噪聲作為實測一階矩，引入的數(shù)據(jù)誤差為0.0035pTs/A，從CDI中估算出背景介質(zhì)的電導率，之后計算出背景一階矩，從實測一階矩中減去背景一階矩得到異常區(qū)域的參考一階矩.

圖8 仿真數(shù)據(jù)一階矩、背景一階矩、正演一階矩在各條測線上曲線圖Fig.8 TEM Bz-moments of synthetic data,background,forward theoretical response at a range of Northings

圖9 CDI切面圖Fig.9 Conductivity-depth sections

5.3.2 無約束反演

時間常數(shù)的初值設定為零，在反演過程中，時間常數(shù)無任何約束，初始擬合差為69.6，耗時24 s后達到收斂，反演的時間常數(shù)切面圖如圖10所示，淺層較強的幾何耦合因子使反演的異常區(qū)域貼近地面，同時，時間常數(shù)出現(xiàn)負值，可見反演結(jié)果并不能反映異常體的真實地電結(jié)構(gòu)，因此，為了使反演貼近實際地下結(jié)構(gòu)，反演過程中，設定約束條件是非常必要的.

5.3.3 電導率加權(quán)反演

利用插值CDI得到相應的電導率加權(quán)值，加權(quán)值切面圖如圖11a所示，根據(jù)電導率約束條件，反演出最優(yōu)時間常數(shù)，其切面圖如圖11b所示，與無限制條件的反演結(jié)果不同，電導率加權(quán)的反演結(jié)果中，高時間常數(shù)區(qū)域局限于實際異常體的體積范圍內(nèi)，更加準確反映地下目標的位置和體積.

圖10 無約束反演時間常數(shù)切面圖Fig.10 Recovered time constant model based on the unconstrained starting model

圖11 (a) 電導率加權(quán)值切面圖;(b) 反演時間常數(shù)切面圖Fig.11 (a) The sections of conductivity weights;(b) The same sections of the recovered time constant model

5.3.4 深度加權(quán)反演

采用兩種深度加權(quán)參數(shù)進行反演，第一種參數(shù)為s0=0.004，d0=250，對應的深度加權(quán)值切面圖和反演時間常數(shù)切面圖為圖12a和12b.第二種參數(shù)為s0=0.004，d0=150，對應的深度加權(quán)切面圖和反演時間常數(shù)切片圖為圖13a和13b.其高時間常數(shù)區(qū)域與異常體的位置和體積較好的吻合，第二種參數(shù)反演出高時間常數(shù)區(qū)域較第一種情況擴大一些，且高時間常數(shù)區(qū)域的頂部比實際平板位置高一些，這是由于第二種參數(shù)的參考深度d0小一些，使異常體頂部的淺層耦合因子影響比較大.因此，在反演中合理設置深度加權(quán)參數(shù)是非常重要的.

5.3.5 快速三維反演效率

表1給出了分別采用傳統(tǒng)最速下降法和優(yōu)化最速下降法的反演效率對比.觀察表1，相比嚴格三維反演方法，本文中三維反演方法無論采用傳統(tǒng)還是優(yōu)化最速下降法，都非常高效地得到反演結(jié)果.同時，經(jīng)過松弛因子優(yōu)化的最速下降法，收斂速度可提高幾十倍甚至百倍，反演經(jīng)過幾秒到幾十秒就可得到收斂結(jié)果.

表1 傳統(tǒng)最速下降法和改進最速下降法反演效率對比

6 總結(jié)

本文提出了一種基于瞬變電磁矩變換的快速三維反演方法.文中應用瞬變電磁一階矩變換，將任意接收點處的一道瞬變電磁數(shù)據(jù)壓縮成了點數(shù)據(jù)，數(shù)據(jù)量的大幅壓縮不僅加速了反演的處理速度，同時也使三維反演能夠在普通計算機上實現(xiàn)；在阻性限制約束下，目標體內(nèi)部感應作用可以忽略，可將目標體網(wǎng)格化，通過幾何耦合因子與時間常數(shù)乘積計算各目標體微元的響應，對所有微元響應進行線性疊加可得到目標體的總響應，基于此特性，建立了一種近似簡化的三維正演算法，將此正演算法應用到三維反演中，加速了三維反演的速度；在反演問題的優(yōu)化迭代算法中，引入松弛因子，改進了傳統(tǒng)的最速下降法，使收斂速度進一步提高；為了促使反演結(jié)果貼近真實地下結(jié)構(gòu)，文中給出了兩種約束條件，實際應用中可靈活運用.最后需要強調(diào)的是，本文雖然以矩形回線源為例推導了三維正演算法，但是本文所提方法對瞬變電磁的發(fā)射源沒有限制，既適用于電性源，也適用于磁性源，同時，該方法不僅可用于地面，也可用于航空或者半航空.

圖12 深度加權(quán)參數(shù)s0=0.004，d0=250(a) 深度加權(quán)值切面圖;(b) 反演時間常數(shù)切面圖.Fig.12 Weight with parameters s0=0.004 and d0=250(a) The sections of depth weights;(b) The same sections of the recovered time constant model.

圖13 深度加權(quán)參數(shù)s0=0.004,d0=150(a) 深度加權(quán)值切面圖;(b) 反演時間常數(shù)切面圖.Fig.13 Weight with parameters s0=0.004 and d0=150(a) The sections of depth weights;(b) The same sections of the recovered time constant model.

三維瞬變電磁反演方法是國內(nèi)外地球物理方向的重要研究熱點之一，復雜的嚴格三維反演具有各種技術(shù)條件限制，使其無法投入實際應用中.本文在普通計算機上實現(xiàn)了一種近似的快速三維反演，結(jié)果表明，該方法在瞬變電磁實時解釋工作中具有良好的應用前景.

致謝 作者向中國科學院電子所李光博士、王友成博士在論文準備過程中提供的幫助致謝.

Cox L H,Wilson G A,Zhdanov M S.2012.3D inversion of airborne electromagnetic data.Geophysics,77(4):WB59-WB69.

Fitterman D V,Stewart M T.1986.Transient electromagnetic sounding for groundwater.Geophysics,51(4):995-1005,doi:10.1190/1.1442158.

Grant F S,West G F.1965.Interpretation theory in applied geophysics.New York:McGraw-Hill Book Co.

He Z X,Zhao Z,Liu H Y,et al.2012.TFEM for oil detection:Case studies.The Leading Edge,31(5):518-521.

Ji Y J,Wang Y,Xu J,et al.2013.Development and application of the grounded long wire source airborne electromagnetic exploration system based on an unmanned airship.Chinese Journal of Geophysics (in Chinese),56(11):3640-3650,doi:10.6038/cjg20131105.

Kaufman A.1978.Frequency and transient responses of electromagnetic fields created by currents in confined conductors.Geophysics,43(5):1002-1010.King A.1997.Inversion of localized electromagnetic anomalies [Ph.D.thesis].Australia:Macquarie University.

Meju M A.2002.Geoelectromagnetic exploration for natural resources:Models,case studies and challenges.Surveys in Geophysics,23(2-3):133-206.

Nabighian M N.1988.Electromagnetic methods in applied geophysics-theory volume I.Society of Exploration Geophysicists,313-503.

Oldenburg D W,Li Y G.2005.Inversion for applied geophysics:A tutorial.∥Butler D K ed.Near-Surface Geophysics.SEG,89-150.

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

Qi Y Z,Huang L,Zhang J G,et al.2013.Study on the electromagnetic fields of an extended source over layered models.Acta Physica Sinica (in Chinese),62(23):234201.

Raydan M,Svaiter B F.2002.Relaxed steepest descent and cauchy-barzilai-borwein method.Computational Optimization and Applications,21(2):155-167.

Reid J E,Macnae J C.2002.Resistive limit modeling of airborne electromagnetic data.Geophysics,67(2):492-500.

Schaa R,Fullagar P K.2012.Vertical and horizontal resistive limit formulas for a rectangular-loop source on a conductive half-space.Geophysics,77(1):E91-E99.

Smith R S.2000.The realizable resistive limit:A new concept for mapping geological features spanning a broad range of conductances.Geophysics,65(4):1124-1127.

Smith R S,Lee T J.2002.The moments of the impulse response:A new paradigm for the interpretation of transient electromagnetic data.Geophysics,67(4):1095-1103.

Smith R S,Lee T J.2001.The impulse-response moments of a conductive sphere in a uniform field,a versatile and efficient electromagnetic model.Exploration Geophysics,32(2):113-118.

Stolz E M,Macnae J.1998.Evaluating EM waveforms by singular-value decomposition of exponential basis functions.Geophysics,63(1):64-74.

Xiong Z H,Tripp A C.1995.A block iterative algorithm for 3-D electromagnetic modeling using integral equations with symmetrized substructures.Geophysics,60(1):291-295.

Xue G Q,Li X,Di Q Y.2008.Research progress in TEM forward modeling and inversion calculation.Progress in Geophysics (in Chinese),23(4):1165-1172.

Xue G Q,Cheng J L,Zhou N N,et al.2013.Detection and monitoring of water-filled voids using transient electromagnetic method:A case study in Shanxi,China.Environmental Earth Sciences,70(5):2263-2270,doi 10.1007/s12665-013-2375-2.

Yin C C,Ren X Y,Liu Y H,et al.2015.Review on airborne electromagnetic inverse theory and applications.Geophysics,80(4):W17-W31.

Yin C C,Zhang B,Liu Y H,et al.2015.Review on airborne EM technology and developments.Chinese Journal of Geophysics (in Chinese),58(8):2637-2653,doi:10.6038/cjg20150804.

附中文參考文獻

齊有政,黃玲,張建國等.2013.層狀介質(zhì)上時空展源瞬變電磁響應的計算方法研究.物理學報,62(23):234201.

薛國強,李貅,底青云.2008.瞬變電磁法正反演問題研究進展.地球物理學進展,23(4):1165-1172.

殷長春,張博,劉云鶴等.2015.航空電磁勘查技術(shù)發(fā)展現(xiàn)狀及展望.地球物理學報,58(8):2637-2653,doi:10.6038/cjg20150804.

嵇艷鞠,王遠,徐江等.2013.無人飛艇長導線源時域地空電磁勘探系統(tǒng)及其應用.地球物理學報,56(11):3640-3650,doi:10.6038/cjg20131105.

(本文編輯 胡素芳)

A rapid 3D inversion based on TEM moment transformation

RAO Li-Ting1,2,WU Xin1,2,WU Chao1,2,ZHANG Xiao-Juan1,FANG Guang-You1

1 Key Laboratory of Electromagnetic Radiation and Sensing Technology, Chinese Academy of Sciences,Beijing 100190,China2 University of Chinese Academy of Sciences,Beijing 100039,China

A full EM solution of the 3D inversion problem is difficult to carry out on common computer due to its complex forward calculation and high demand for user resources.These restrictions make rigorous 3D inversion still impractical in the interpretation of TEM data sets.A rapid 3D inversion scheme for interpreting TEM data is presented,utilizing the concept of TEM moment.Approximate forward modeling of the 3D inversion scheme has been accomplished by adopting the first order moment transform which is the resistive limit.The modified steepest decent method is employed for inversion of the underground target response.The inherent non-uniqueness associated with the underdetermined inversion problem is balanced by including constraints in the inversion such as depth weights or conductivity weights.The fully 3D electromagnetic modeling software is employed to calculate the synthetic time-domain data to verify the proposed 3D inversion scheme.A comparison of the forward theoretical results and the moment transformed results of synthetic time-domain data shows that forward calculation gives reliable results.In order to test the inversion algorithm and examine the variability of the inversion results due to non-uniqueness,the synthetic inversion examples are explored based on different constraints.The inversion was successfully completed in dozens of seconds to recover time constant models which approximately indicated the location of the true conductive body.In conclusion,the proposed rapid 3D inversion scheme can be efficiently implemented on the normal computer which indicates a promising application prospect in real-time interpretation of TEM data.

Transient electromagnetic method (TEM)；TEM moment；Resistive limit；Rapid 3D inversion

饒麗婷，武欣，吳超等.2016.基于瞬變電磁矩變換的快速三維反演方法.地球物理學報，59(11):4338-4348,

10.6038/cjg20161133.

Rao L T,Wu X,Wu C,et al.2016.A rapid 3D inversion based on TEM moment transformation.Chinese J.Geophys.(in Chinese),59(11):4338-4348,doi:10.6038/cjg20161133.

國家重大科研裝備研制項目(ZDYZ2012-1-03-01)資助.

饒麗婷，女，1989年生，博士，主要從事瞬變電磁正反演理論與方法技術(shù)研究.E-mail:raoliting11@mails.ucas.ac.cn

*通訊作者 武欣，男，1982年生，助理研究員，主要從事時間域電磁法理論與系統(tǒng)研發(fā).E-mail:wu_xin18@mail.ie.ac.cn

10.6038/cjg20161133

P631

2015-11-11，2016-07-06收修定稿