地面與半航空瞬變電磁法三維聯合反演

劉明宏， 蔡紅柱,2*， 楊浩， 熊詠春, 胡祥云,2

1 中國地質大學(武漢)地球物理與空間信息學院， 武漢 430074 2 地質過程與礦產資源國家重點實驗室， 武漢 430074

0 引言

瞬變電磁法(Transient electromagnetic method, TEM)是一種電磁感應類地球物理探測方法，通過觀測激勵源斷電后的二次場響應，對觀測信號進行成像或反演分析可獲得地下空間的電性分布信息.目前，瞬變電磁法已經被廣泛應用于礦產勘查(武欣等, 2016; Yang and Oldenburg, 2012)、油氣勘探(Li and Constable, 2010)和地下水調查(Auken et al., 2017; Trabelsi et al., 2013)等領域.

瞬變電磁法通常采用磁性回線源或電性接地導線源作為激勵，可在地面或空中等平臺進行觀測.傳統的地面中心回線源瞬變電磁觀測裝置具有探測深度大、分辨率高、抗干擾能力強等特點，應用非常廣泛，但在山地、叢林、沼澤、冰川等艱險環境，地面工作難度極大.航空電磁法是探測地面環境復雜區域的有效手段，然而航空電磁測量對儀器設備技術要求高，成本昂貴、飛行風險大.將發射源設置在地面，使用無人機在空中觀測的SATEM系統成為了一種新興的地球物理探測技術(Meng et al., 2020; Wu et al., 2019).SATEM系統相對于地面系統，采用空中接收工作方式，具有效率高、不受地表環境限制的優點；相對于航空系統，采用地面發射方式，具有發射功率大、數據信噪比高、勘探深度大的優勢(Ito et al., 2014; Smith et al., 2001).圖1展示了地面與半航空瞬變電磁裝置系統的示意圖.地面中心線源裝置通常在回線框內的一定范圍內采集數據，半航空瞬變電磁一般采用接地長導線作為發射源，在偏離源的區域飛行采集數據.

圖1 地面瞬變電磁與半航空瞬變電磁裝置示意圖Fig.1 Schematic diagram of the ground-based TEM (GTEM) and semi-airborne TEM

半航空電磁方法的開發旨在于增加航空電磁探測的深度.Nabighian(1988)提出了SATEM的概念.Elliott(1996, 1998)開發了基于回線源的FLAIRTEM系統，并成功應用于巴布亞新幾內亞的礦床調查.加拿大Fugro公司開發了回線源TerraAir系統(Smith et al., 2001).Mogi等(1998)開發了首套采用接地導線發射源SATEM系統(GREATEM)，成功應用于地熱調查、火山構造探測和海水入侵等多個領域，展示了電性激勵源SATEM系統具有高信噪比、探測深度大的優勢(Abd Allah et al., 2013; Ito et al., 2014; Mogi et al., 2009).我國在SATEM系統的研發與應用上也取得了不錯的進展.嵇艷鞠等(2013)和Wang等(2013)研制了一套以無人飛艇為搭載平臺的時域電磁接收系統，并開展了海水入侵調查和水資源探測試驗.劉富波等(2017)基于無人直升機搭載平臺，在山東昌邑蓮花山鐵礦區開展了電磁探測飛行試驗.隨著無人機技術的迅速發展，以無人機為搭載平臺的SATEM成為了近年的研究熱點.方濤等(2015)、Xue等(2018)、吳啟龍(2019)和謝小國等(2021)將無人機SATEM系統應用于地下采空區、地下巷道、隧道工程調查、古河道探測等領域.

在解釋方面，一些學者采用視電阻率成像的方法進行瞬變電磁數據解釋(Christiansen et al., 2006; 馬振軍等, 2021; 李貅等, 2015; 張瑩瑩等, 2015; 趙越等, 2015).此外，采用一維反演獲得層狀電導率模型是瞬變電磁數據解釋的主要方法(Farquharson and Oldenburg, 1993; 李永興等, 2010; 李展輝和黃清華, 2018)，但根據一維反演結果構建的擬二維剖面通常會產生電阻率突變以至于難以進行地質解釋.隨著橫向約束和空間約束方法的引入，瞬變電磁一維反演結果的空間連續性得到了增強(Auken et al., 2008; Viezzoli et al., 2008).然而在復雜的地質情況下，傳統的視電阻率成像方法和一維反演很難恢復精確的地下結構.

三維反演是對TEM數據精確解釋的有效方法.正演是反演的基礎，前人在TEM三維正演上進行了大量研究.積分方程法(Wannamaker et al., 1984; Zhdanov et al., 2006; Zhu et al., 2019; 殷長春和劉斌, 1994)、有限差分法(Commer et al., 2015; Commer and Newman, 2004; Maa?, 2007; Wang and Hohmann, 1993; Yu et al., 2017; 孫懷鳳等, 2013; 余翔等, 2017)被應用于電磁三維正演模擬.隨著計算硬件的快速發展，尤其是高性能集群的應用，TEM三維反演已成為可能.基于積分方程法，結合時頻轉換技術和移動足跡(“moving footprint”)方法，Cox等(2012)開發了TEM三維反演方法，Wang等(2021)開發了適用于多道瞬變電磁數據的三維并行反演算法.基于有限差分法和時頻轉換方法，Liu和Yin(2016)實現了航空TEM三維反演.隨后，Su等(2021)將其推廣到TEM各向異性三維反演.基于規則網格的有限體積法的TEM三維反演也有少量報道(Oldenburg et al., 2013; Ren et al., 2018; Yang et al., 2014).

在先前報道的基于積分方程、有限差分或有限體積法的TEM三維反演，采用了規則網格對模型進行離散，難以考慮復雜的地質構造(如地形)的影響.近年來，基于非結構化四面體網格的有限元法在地球物理電磁模擬中的應用受到廣泛關注(Badea et al., 2001; Cai et al., 2017; Fu et al., 2015; Jahandari et al., 2017; Li et al., 2018; Um et al., 2010, 2012; Yin et al., 2019).Um等(2010)使用基于非結構化四面體網格的有限元方法來計算時域電磁場.此外，他們還提出了一種時間步長加倍的方法減少所需計算的時間步進(Um et al., 2012)，隨后引入了并行策略加速計算(Fu et al., 2015).Cai等(2017)基于非結構化四面體網格的時域有限元方法，結合自適應時間步進法和混合邊界條件策略，提高了正演的準確性和計算效率.殷長春等(2019)提出了基于自適應有限元方法的時間域三維正演算法，提高計算效率.Jiang等(2019)采用全場控制方程的有限元法進行了地-井時域電磁模擬，并利用煤礦含水層模型進行了驗證.Li等(2020)開發了考慮地形影響的接地線源瞬變電磁法有限元三維正演算法，并將其應用于Ovoid地區硫化物礦床的數值模擬.

近年來，基于矢量有限元的三維TEM反演的研究也取得進展.Liu等(2019)基于無約束四面體網格矢量有限元法實現了三維TEM反演，反演采用BFGS優化方法.Qi等(2022)開發了一種結合局部網格技術和解耦正反演網格技術的航空TEM數據三維反演方法.Zhang等(2021)提出了一種三重網格方法的TEM三維反演策略，正演計算使用細的非結構化網格以保證精度，靈敏度計算使用粗的非結構化網格以減少計算，反演采用規則的結構化網格以便于約束反演參數.基于有限元的TEM三維反演多報道于地面裝置和航空裝置，但少有關于SATEM的三維反演報道.

可靠性和分辨率是地球物理反演最關注的問題.單一的電磁探測裝置對地下結構的靈敏度具有局限性，難以獲得地質結構的全面信息.不同TEM裝置對地下的電性結構的靈敏度不同，聯合地面與半航空TEM數據將形成互補優勢，比單獨觀測方式對地下的電性結構具有更高的靈敏度.另一方面，地球物理反演普遍存在多解性問題，多種地球物理方法的數據聯合反演是提高反演結果可靠性的有效方法.聯合反演已被廣泛應用于地球物理數據解釋(Gallardo and Meju, 2004; Lelièvre et al., 2012; Zhdanov et al., 2012)，但關于地面TEM與SATEM的三維聯合反演未見報道.

為了實現可靠的瞬變電磁三維反演，本文開發了并行的瞬變電磁三維正演與反演算法.正演采用基于非結構化四面體網格的矢量有限元法，可以精確模擬復雜結構.為了提高計算效率，在求解正演和靈敏度計算時，本文采用自適應時間步長的方法減少時間步進，并利用Intel MKL PARDISO并行直接求解器求解大型線性方程組.在反演中，為了提高反演收斂速度，采用具有擬二階收斂性的高斯-牛頓法反演策略，通過最小二乘QR分解(LSQR)方法獲得模型更新方向.為結合不同裝置的探測優勢以及降低反演的多解性影響，基于構建統一目標函數和數據權重調整的反演策略，實現了地面TEM和SATEM數據的聯合反演.最后通過理論模型研究，分析了地面TEM和SATEM的探測能力，并對比單獨反演與聯合反演結果，驗證了聯合反演方法的優勢.

1 瞬變電磁三維正演

三維正演模擬是認識電磁規律與開展電磁反演的基礎.從時域的準靜態Maxwell方程組出發，在忽略位移電流和考慮真空磁導率的情況下得到(Ward and Hohmann, 1988)：

(1)

(2)

(3)

本文采用基于非結構化四面體網格的一階矢量有限元方法求解方程(3)(Cai et al., 2017; Jin, 2014; Liu et al., 2019; Um et al., 2012).在不同的時刻，電場被賦值到單元邊緣上，四面體單元內的電場可用沿四面體邊界的電場的線性組合表示：

(4)

方程(3)的時間導數項需要使用有限差分近似表示.在時域有限元方法中，后退歐拉法因其無條件穩定而常被采用(Haber, 2014; Jin, 2014; Um et al., 2012).在相同步進次數條件下，二階后退歐拉法的時間差分比一階方法具有更高精度.殷長春等(2019)和Liu等(2019)使用了具有相同步長的二階后退歐拉方法近似時間導數項.為了減少所需模擬的時間步進次數，本文采用了靈活的自適應時間步進方法，在一定的步數后，嘗試增加時間步長(Um et al., 2012; Cai et al., 2017).

采用二階后退歐拉方法，在tn時刻的電場時間導數可以表示為

(5)

通過伽遼金有限元分析獲得如下的有限元線性方程組(Jin, 2014)：

AnEn=bn，

(6)

其中

(7)

(8)

式中An∈(Ned×Ned)，是一個大型的稀疏矩陣，bn∈(Ned×1)，K和M是剛度矩陣(Jin, 2014).

求解方程6需要給定適當邊界條件和初始條件.本文采用Dirichlet邊界條件，假設電場在建模區域的外邊界上趨近于零.為了匹配Dirichlet邊界條件，將模擬區域拓展至一個很大的邊界.在測點和發射源的附近加密網格以確保計算精度.在核心區域以外逐漸增大網格的尺寸，以減少網格數量.對于初始條件的選擇，通常可以采用兩種方式.當計算全波形響應時采用電場為零的初始條件.當忽略上升電流與持續電流只考慮關斷電流的響應時，通常求解直流電場的拉普拉斯方程獲得初始電場.本文采用了如圖2所示的梯形發射波，其脈寬時間較短，不能忽略上升電流(I)、持續電流的影響，因此需要進行全波形計算，考慮電場為零的初始條件：E0=0.本文通過將離散的發射電流密度加入公式(8)實現激勵源的加載，可以模擬任意發射波形.電流關斷后啟用自適應時間步進，初始時間步長大小取決于發射波形的離散.

圖2 梯形波Fig.2 Trapezoidal waveform

確定邊界條件和初始條件后，可通過迭代法或直接法求解方程(6)得到某時刻的電場En.由于每個時間步進都必須求解有限元方程組，考慮到時間步長數量很大，采用迭代法并不實用(Jin, 2014; Um et al., 2012).采用直接法求解更具優勢，當時間步長相同時，剛度矩陣保持不變，An矩陣分解的結果可以被重復使用(Jin, 2014; Oldenburg et al., 2013).在自適應步進方法中，矩陣An只需改變幾次，可以顯著減少矩陣An的分解次數.為了提高計算效率，采用 Intel MKL PARDISO并行直接求解器求解方程(6)(Schenk et al., 2001).

通過正向步進的方式逐步求解方程(6)得到電場，測點的數據可以通過以下插值矩陣獲得：

d=QE，

(9)

其中，Q是與時間和空間相關的插值矩陣，E為所有時刻的電場.

2 高斯-牛頓法反演

TEM反演是不適定問題.為了克服這個問題以得到具有地球物理意義的解，Tikhonov正則化被引入到反演目標泛函中(Tikhonov and Arsenin, 1977)：

φ(m)=φd(m)+λφm(m)，

(10)

其中，φd是數據擬合項，φm是模型穩定項，λ是用于平衡數據擬合和模型約束正則化參數.φd和φm分別定義為

(11)

(12)

其中，m是模型參數，為了保證反演中模型參數具有物理意義，取對數域的電導率作為模型參數，m=log(σ).Wd為數據加權矩陣，瞬變電磁數據跨越多個數量級，為了平衡不同時間通道的數據權重，采用觀測數據的倒數進行構建，并在數據變號處數據權重調整為零，以避免無窮大值，類似的方法已成功應用于瞬變電磁反演(Cox et al., 2012).F(m)表示模型m的正演響應，dobs為觀測數據，mref為含先驗信息的參考模型.Wm為模型粗糙度矩陣，由反演單元與鄰近單元的體積和距離的權重構建(Cai et al., 2021).

對目標函數進行二階泰勒級數展開，并忽略高階項可以獲得以下法方程：

式中δmn為模型更新量，數據誤差rn-1=F(mn-1)-dobs，n為迭代次數.法方程通常采用共軛梯度法(Conjugate gradient，CG)求解，該方法需要求解近似海森矩陣以獲得有效的預條件因子.在數值計算上，法方程等價于如下最小二乘形式(Hansen, 1998; Aster et al., 2018)：

(14)

LSQR算法是利用Lanczos子空間投影方法求解最小二乘問題的有效方法(Paige and Saunders, 1982)，相對于PCG算法，LSQR算法具有更高的穩定性(Grayver et al., 2013).求解方程組(14)前需要獲得靈敏度矩陣Jd，Jd的準確性影響著反演的穩定性和準確性.為了求取數據的靈敏度矩陣，需要先求取電場對模型參數的導數Je,n，式(6)對m取導數和整理得

(15)

第n步的數據對模型參數的靈敏度(Jacobian)可以表示為

(16)

式中，Nd和Ne分別表示數據的數量和模型空間四面體單元數量.通常不需要求解整個模型空間的靈敏度矩陣，只考慮反演區域的靈敏度矩陣，Jd,n的大小變為Nd×Nm，Nm為反演區域的模型參數的數量.數據的數量Nd通常遠少于模型參數的數量Nm，因此本文求取靈敏度矩陣的轉置.因為系數矩陣An是對稱矩陣，則有

(17)

獲得模型更新方向δm后，通過線搜索策略找到最優的模型更新步長l，模型更新可以表示為(Nocedal and Wright, 2006)：

mn=mn-1+lδm.

(18)

正則化參數在權衡數據擬合項和模型約束項中起著重要作用，關系著反演過程的穩定.在初始反演迭代中，本文采用近似譜分析來獲得正則化參數(Long et al., 2020)，在后續迭代中使用冷卻方法減小正則化參數.此外，觀測數據與預測數據之間的歸一化擬合差可用以下公式表示(Zhdanov, 2009)：

(19)

當歸一化擬合差達到設定閾值或收斂停止時，反演終止迭代.

3 聯合反演

聯合反演可以結合不同探測方法的優勢以降低反演的多解性.對于不同物性的聯合反演方法，通常需要構建不同物性之間的相關性約束.在地球物理反演中構建相關性約束的方法主要有：交叉梯度(Gallardo and Meju, 2004)、Gramian約束(Zhdanov et al., 2012)、模糊c均值聚類法(Lelièvre et al., 2012)等.地面TEM與SATEM是對同一物性(電導率)進行反演，因此，在聯合反演策略中，直接采用與單獨反演相同的目標函數，而不需要構建相關性約束關系(Commer and Newman, 2009).

在聯合反演中，觀測數據dobs為地面TEM的數據和SATEM的數據組合成的列向量：

(20)

式中，dGTEM和dSATEM分別表示地面TEM的數據子集和SATEM的數據子集.

與單獨反演不同，在聯合反演中數據的權重需要進行適當調整.地面TEM與SATEM的數據采集密度不同，在相同測量面積上，半航空的數據量可能遠大于地面數據.如果簡單采用與單獨反演相同的數據加權矩陣，當地面TEM數據量遠小于SATEM數據時，地面數據在反演中難以發揮作用.我們希望在聯合反演中，不同數據子集能發揮相應的作用，以結合兩種裝置的探測能力.因此，本文采用改變不同數據子集的數據加權矩陣的方式，來調整地面TEM數據和SATEM數據的權重(Commer and Newman, 2009).

本文以相同觀測區域上不同地球物理數據的數量為依據確定權重系數α1和α2.假設子集dGTEM和dSATEM的數據量分別為N1和N2，則

α1×N1=α2×N2，

(21)

當指定α1=1時，α2=N1/N2.相應的Wd則為

(22)

上式中WGTEM、WSATEM分別表示地面TEM和SATEM的數據加權矩陣.在聯合反演中，只需通過改變α2，就可以調節不同數據子集的權重，確保不同數據子集權重的平衡.聯合反演中的靈敏度計算、法方程求解等流程與單獨反演完全相同.

4 理論模型研究

為了驗證所提出的反演算法的有效性，本文設計了兩個理論模型.在數值實驗中，地面TEM與SATEM均采用如圖2所示的梯形發射波.觀測時間取等對數間隔10-4.5～10-1.5s，共25個時間道.在所有的模型案例中，使用添加2%的高斯噪聲的垂直分量dHz/dt合成數據進行反演.

4.1 模型1：無地形模型

模型1電性結構和觀測裝置如圖3所示.導線源長1500 m，第1號源在X=-1500 m處，第2號源在X=1500 m處.設置6個大小均為500 m×500 m的大回線源，僅在回線圈內觀測，采用兩種觀測點距，50 m點距的數據量為6×81×25=12150個，200 m點距的數據量為6×9×25=1350個.SATEM測區范圍在X:-700～700 m，Y:-450～450 m之間，點距為50 m，數據量為551×25=13775個.淺部兩個異常體大小為300 m×600 m×60 m，中心點深度為60 m，電導率分別為0.1 S·m-1和0.001 S·m-1；深部異常體大小為400 m×600 m×200 m，中心點深度為350 m，電導率為0.2 S·m-1；背景電導率為0.01 S·m-1.

圖3 模型1的切片圖(a) Z=60 m切片，紅色線段、白色線框分別為SATEM發射源、地面TEM回線框發射源，黑色、白色測點的距離分別為50 m、200 m； (b) Y =0 m切片.Fig.3 Slice view of the model 1(a) Slice at Z=60 m, red line and white loops represent the SATEM electrical sources and the ground TEM loop sources respectively, the spacing of the black and white dots are 50m and 200 m, respectively; (b) Slice at Y = 0.

在環境復雜地區難以進行密集數據采集.為了測試地面TEM的中心回線源裝置的探測能力與不同測點數量對地面瞬變電磁三維反演結果的影響，對不同點距的數據進行反演.正演與反演使用不同的網格，反演網格取消對異常體處加密，在模型1的所有反演中均使用同一套網格，整個模型區域大小為20 km×20 km×20 km.反演區域大小為3 km×2 km×1 km，反演區域四面體單元數量為393940，整個模型的四面體單元數量為1037619，反演初始模型電導率為0.01 S·m-1的均勻半空間.

本文的數值實驗在中國地質大學(武漢)高性能計算集群上完成，每個節點包含2個Intel(R) Xeon(R) 2.5 GHz的CPU，每個CPU包含20個核，單節點內存為384GB.模型1的反演使用1個節點計算，50 m點距和200 m點距的每次反演迭代分別需要342 min和346 min，測點數量對反演時間影響較小.反演收斂情況如圖7a所示，最終反演結果的數據歸一化擬合差在2%左右.

圖4為地面TEM單獨反演的結果.淺層兩個薄板異常體的位置和電導率恢復較好，但深部異常體的反演結果類似一個“環狀”結構，只能恢復異常體的邊界，在異常體中心形成了一個空洞.隨著測點數量的減少，反演結果分辨率降低，假異常增加.對于實際探測，難以進行合理的地質解釋.

圖4 模型1的地面TEM反演結果(a,b) 50 m點距數據反演結果; (c,d) 200 m點距數據反演結果； (a,c) 三維視圖； (b,d) Z=350 m切片.Fig.4 The ground TEM inversion results for model 1(a,b) Inversion results with station spacing of 50 m; (c,d) Inversion results with station spacing of 200 m; (a,c) Three-dimensional view; (b,d) Slice view at Z=350 m.

為了測試SATEM裝置的探測能力，對兩個不同發射源的模擬數據分別進行反演.網格、初始模型和其他參數與上相同，每次反演迭代約4 h.反演收斂情況如圖7b所示，最終反演結果的數據歸一化擬合差在5%左右.因為在較小的偏移距內進行觀測，衰減曲線會出現“變號”，數據擬合差稍大于噪聲水平.

SATEM單獨反演結果如圖5所示.兩個反演結果均恢復了大致的地下電性分布情況.但與地面TEM裝置反演結果相比，SATEM對異常體邊界的分辨能力較差.SATEM單獨反演的結果也存在明顯的假異常，且在不同發射源的情況下，反演結果的假異常不同，這給數據解釋帶來巨大挑戰.

圖5 模型1的SATEM反演結果(a,b) 第1號源數據反演結果； (c,d) 第2號源數據反演結果； (a,c) 三維視圖； (b,d) Z=350切片.Fig.5 SATEM inversion results for model 1(a,b) Inversion results with the data excited by the first source; (c,d) Inversion results with the data excited by the first source; (a,c) Three-dimensional view; (b,d) Slice view at Z=350 m.

單一觀測方式探測能力的局限性和反演的多解性嚴重影響反演結果的準確性.聯合多種測量裝置進行反演是降低多解性的有效手段.我們分別將不同采集密度的地面TEM數據與SATEM第2號源的數據進行聯合反演.聯合反演采用與單獨反演相同的網格、初始模型和反演參數.聯合反演耗時與地面TEM單獨反演相當，每次反演迭代約341 min.反演收斂情況如圖7c所示，地面TEM點距50 m、200 m與SATEM的聯合反演最終結果的數據歸一化擬合差分別為5.6%和3.9%.

圖6為SATEM與地面TEM聯合反演結果.聯合反演對三個異常體的邊界和電導率均得到有效重建.地面TEM單獨反演的深部異常體的“環狀”結構以及SATEM單獨反演的假異常，在聯合反演中均得到有效壓制.尤其采用50 m點距的地面TEM數據時，聯合反演結果基本上沒有明顯假異常.當地面TEM數據為200 m點距時，雖然數據量只有50 m點距的1/9，但聯合反演結果也明顯好于單獨反演的結果.

圖6 模型1的地面TEM與SATEM聯合反演結果(a,b) 地面TEM點距50 m數據與SATEM數據聯合反演； (c,d) 地面TEM點距200 m數據與SATEM數據聯合反演； (a,c) 三維視圖； (b,d) Z=350切片.Fig.6 Results of ground TEM and SATEM joint inversion for model 1(a,b) Joint inversion between the ground TEM data with 50 m station spacing and the SATEM data; (c,d) Joint inversion between the ground TEM data with 200 m station spacing and the SATEM data; (a,c) Three-dimensional view; (b,d) Slice view at Z=350 m.

圖7 模型1反演的收斂曲線(a) 地面TEM反演； (b) SATEM反演； (c) 地面TEM與SATEM聯合反演.Fig.7 Convergence plot for different inversions of model 1(a) The ground TEM inversion; (b) SATEM inversion; (c) Joint inversion between the ground TEM data and the SATEM data.

4.2 模型2：帶地形模型

為了更接近真實地質情況，建立了一個帶有地形和多個異常體的模型.模型2的地形、電性結構和觀測裝置如圖8所示.該模型地形起伏高差約100 m，X軸負方向較為平坦，X軸正方向地形起伏，分別模擬喀斯特地貌的平原與峰叢.導線源長1000 m，在X=-1500 m處.設置8個大小為500 m×500 m的大回線源，僅在回線內觀測，點距50 m，觀測數據量為8×81×25=16200個.SATEM測區范圍在X:-950～950 m，Y:-450～450 m之間，點距50 m，數據量為741×25=18525個.異常體1的大小為300 m×600 m×100 m，中心點埋深約100 m，電導率為0.001 S·m-1；異常體2的大小為600 m×600 m×100 m，傾斜放置，電導率為0.2 S·m-1；異常體3和異常體4的大小均為600 m×200 m×100 m，中心點埋深約100～200 m，電導率分別為0.2 S·m-1和0.001 S·m-1；背景電導率為0.01 S·m-1.

反演網格取消對異常體處加密，模型2的所有反演中使用同一套網格.反演區域X方向范圍：-2000～2000 m，Y方向范圍：-1000～1000 m，Z方向厚度約1000 m.反演區域四面體單元數量為 367077，整個模型的四面體單元數量為713753，反演初始模型為0.01 S·m-1的均勻半空間.反演收斂情況如圖12所示，地面TEM、SATEM單獨反演和聯合反演最終結果的數據歸一化擬合差分別為8.3%、5.4%和10.2%.

圖9、10和11為帶地形模型反演結果的不同方向切片圖.在圖9b、10b和11b的地面TEM單獨反演結果中，地面中心回線觀測方式對該模型的反演結果很不理想，兩個高導異常體連在一起，無法恢復高導體的準確信息.兩個低導異常體的恢復效果較差，只能看到模糊的大致情況.我們認為是因為中心回線裝置的觀測數據中缺少不同偏移距離的信息導致難以區分該模型的兩個高導異常體.在圖9c、10c和11c的SATEM單獨反演結果中，四個異常體的恢復結果較好，對于底部埋深約有600 m的傾斜高導體也恢復了異常體的總體形態，但在異常體附近存在一些小的假異常.圖9d、10d和11d的聯合反演結果中，四個異常體的形態恢復較好，SATEM單獨反演的局部假異常也被有效弱化.在X=-1500 m附近出現了新的假異常，其產生原因可能是由于此處不在觀測數據的覆蓋范圍內，數據對該區域靈敏度較小.此外，該假異常較之單獨反演的假異常更為微弱，對數據解釋影響較小.

圖9 模型2在Z=100 m的模型與反演結果切片(a) 理論模型； (b) 地面TEM反演； (c) SATEM反演； (d) 聯合反演.Fig.9 Slice view of the true model and the inversion results at Z=100 m for model 2(a) Synthetic model; (b) The ground TEM inversion; (c) SATEM inversion; (d) Joint inversion.

圖10 模型2在Y=-200 m的模型與反演結果切片(a) 理論模型； (b) 地面TEM反演； (c) SATEM反演； (d) 聯合反演.Fig.10 Slice view of the true model and the inversion results at Y=-200 m for model 2(a) Synthetic model; (b) The ground TEM inversion; (c) SATEM inversion; (d) Joint inversion.

圖11 模型2在Y =200 m的模型與反演結果切片(a) 理論模型； (b) 地面TEM反演； (c) SATEM反演； (d) 聯合反演.Fig.11 Slice view of the true model and the inversion results at Y = 200 m for model 2(a) Synthetic model; (b) The ground TEM inversion; (c) SATEM inversion; (d) Joint inversion.

圖12 模型2反演的收斂曲線Fig.12 Convergence plot for different inversions of model 2

5 結論

本文提出了一套三維TEM反演算法.正演算法采用基于非結構化四面體網格的矢量有限元法，可以對起伏地形、復雜地質結構、任意形狀發射源和任意波形進行精確模擬.在三維正演算法的基礎上開發出了基于高斯-牛頓法的高效并行反演算法；構建了地面TEM和SATEM的聯合反演方案.該算法可適用于不同TEM系統的單獨反演及不同系統的聯合反演.通過理論模型研究驗證了算法的有效性.在模型1中，對比了地面TEM點距50 m數據和點距200 m數據的單獨反演，點距200 m數據反演結果分辨率較高、假異常較少，說明了增加TEM數據量可以一定程度上提高反演結果質量.對比兩個模型的地面TEM與SATEM的單獨反演可以發現，單一的TEM探測方式的探測能力具有局限性，單獨反演通常伴隨一定的假異常，且不同的探測裝置的單獨反演結果產生不同的假異常，給解釋帶來巨大的挑戰，為獲得可靠結果有必要開展聯合反演.聯合地面TEM和SATEM數據聯合反演，可以結合不同觀測方式的探測優勢，以及壓制由于反演多解性而形成的假異常，反演結果的可靠性和分辨率得到提高.本文所提出的方法算法可以為復雜地質條件下開展地面TEM與SATEM探測提供理論支撐.