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

1 中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院， 武漢 430074 2 地質(zhì)過(guò)程與礦產(chǎn)資源國(guó)家重點(diǎn)實(shí)驗(yàn)室， 武漢 430074

0 引言

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

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

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

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

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

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

在先前報(bào)道的基于積分方程、有限差分或有限體積法的TEM三維反演，采用了規(guī)則網(wǎng)格對(duì)模型進(jìn)行離散，難以考慮復(fù)雜的地質(zhì)構(gòu)造(如地形)的影響.近年來(lái)，基于非結(jié)構(gòu)化四面體網(wǎng)格的有限元法在地球物理電磁模擬中的應(yīng)用受到廣泛關(guān)注(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)使用基于非結(jié)構(gòu)化四面體網(wǎng)格的有限元方法來(lái)計(jì)算時(shí)域電磁場(chǎng).此外，他們還提出了一種時(shí)間步長(zhǎng)加倍的方法減少所需計(jì)算的時(shí)間步進(jìn)(Um et al., 2012)，隨后引入了并行策略加速計(jì)算(Fu et al., 2015).Cai等(2017)基于非結(jié)構(gòu)化四面體網(wǎng)格的時(shí)域有限元方法，結(jié)合自適應(yīng)時(shí)間步進(jìn)法和混合邊界條件策略，提高了正演的準(zhǔn)確性和計(jì)算效率.殷長(zhǎng)春等(2019)提出了基于自適應(yīng)有限元方法的時(shí)間域三維正演算法，提高計(jì)算效率.Jiang等(2019)采用全場(chǎng)控制方程的有限元法進(jìn)行了地-井時(shí)域電磁模擬，并利用煤礦含水層模型進(jìn)行了驗(yàn)證.Li等(2020)開發(fā)了考慮地形影響的接地線源瞬變電磁法有限元三維正演算法，并將其應(yīng)用于Ovoid地區(qū)硫化物礦床的數(shù)值模擬.

近年來(lái)，基于矢量有限元的三維TEM反演的研究也取得進(jìn)展.Liu等(2019)基于無(wú)約束四面體網(wǎng)格矢量有限元法實(shí)現(xiàn)了三維TEM反演，反演采用BFGS優(yōu)化方法.Qi等(2022)開發(fā)了一種結(jié)合局部網(wǎng)格技術(shù)和解耦正反演網(wǎng)格技術(shù)的航空TEM數(shù)據(jù)三維反演方法.Zhang等(2021)提出了一種三重網(wǎng)格方法的TEM三維反演策略，正演計(jì)算使用細(xì)的非結(jié)構(gòu)化網(wǎng)格以保證精度，靈敏度計(jì)算使用粗的非結(jié)構(gòu)化網(wǎng)格以減少計(jì)算，反演采用規(guī)則的結(jié)構(gòu)化網(wǎng)格以便于約束反演參數(shù).基于有限元的TEM三維反演多報(bào)道于地面裝置和航空裝置，但少有關(guān)于SATEM的三維反演報(bào)道.

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

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

1 瞬變電磁三維正演

三維正演模擬是認(rèn)識(shí)電磁規(guī)律與開展電磁反演的基礎(chǔ).從時(shí)域的準(zhǔn)靜態(tài)Maxwell方程組出發(fā)，在忽略位移電流和考慮真空磁導(dǎo)率的情況下得到(Ward and Hohmann, 1988)：

(1)

(2)

(3)

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

(4)

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

采用二階后退歐拉方法，在tn時(shí)刻的電場(chǎng)時(shí)間導(dǎo)數(shù)可以表示為

(5)

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

AnEn=bn，

(6)

其中

(7)

(8)

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

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

圖2 梯形波Fig.2 Trapezoidal waveform

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

通過(guò)正向步進(jìn)的方式逐步求解方程(6)得到電場(chǎng)，測(cè)點(diǎn)的數(shù)據(jù)可以通過(guò)以下插值矩陣獲得：

d=QE，

(9)

其中，Q是與時(shí)間和空間相關(guān)的插值矩陣，E為所有時(shí)刻的電場(chǎng).

2 高斯-牛頓法反演

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

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

(10)

其中，φd是數(shù)據(jù)擬合項(xiàng)，φm是模型穩(wěn)定項(xiàng)，λ是用于平衡數(shù)據(jù)擬合和模型約束正則化參數(shù).φd和φm分別定義為

(11)

(12)

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

對(duì)目標(biāo)函數(shù)進(jìn)行二階泰勒級(jí)數(shù)展開，并忽略高階項(xiàng)可以獲得以下法方程：

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

(14)

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

(15)

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

(16)

式中，Nd和Ne分別表示數(shù)據(jù)的數(shù)量和模型空間四面體單元數(shù)量.通常不需要求解整個(gè)模型空間的靈敏度矩陣，只考慮反演區(qū)域的靈敏度矩陣，Jd,n的大小變?yōu)镹d×Nm，Nm為反演區(qū)域的模型參數(shù)的數(shù)量.數(shù)據(jù)的數(shù)量Nd通常遠(yuǎn)少于模型參數(shù)的數(shù)量Nm，因此本文求取靈敏度矩陣的轉(zhuǎn)置.因?yàn)橄禂?shù)矩陣An是對(duì)稱矩陣，則有

(17)

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

mn=mn-1+lδm.

(18)

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

(19)

當(dāng)歸一化擬合差達(dá)到設(shè)定閾值或收斂停止時(shí)，反演終止迭代.

3 聯(lián)合反演

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

在聯(lián)合反演中，觀測(cè)數(shù)據(jù)dobs為地面TEM的數(shù)據(jù)和SATEM的數(shù)據(jù)組合成的列向量：

(20)

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

與單獨(dú)反演不同，在聯(lián)合反演中數(shù)據(jù)的權(quán)重需要進(jìn)行適當(dāng)調(diào)整.地面TEM與SATEM的數(shù)據(jù)采集密度不同，在相同測(cè)量面積上，半航空的數(shù)據(jù)量可能遠(yuǎn)大于地面數(shù)據(jù).如果簡(jiǎn)單采用與單獨(dú)反演相同的數(shù)據(jù)加權(quán)矩陣，當(dāng)?shù)孛鎀EM數(shù)據(jù)量遠(yuǎn)小于SATEM數(shù)據(jù)時(shí)，地面數(shù)據(jù)在反演中難以發(fā)揮作用.我們希望在聯(lián)合反演中，不同數(shù)據(jù)子集能發(fā)揮相應(yīng)的作用，以結(jié)合兩種裝置的探測(cè)能力.因此，本文采用改變不同數(shù)據(jù)子集的數(shù)據(jù)加權(quán)矩陣的方式，來(lái)調(diào)整地面TEM數(shù)據(jù)和SATEM數(shù)據(jù)的權(quán)重(Commer and Newman, 2009).

本文以相同觀測(cè)區(qū)域上不同地球物理數(shù)據(jù)的數(shù)量為依據(jù)確定權(quán)重系數(shù)α1和α2.假設(shè)子集dGTEM和dSATEM的數(shù)據(jù)量分別為N1和N2，則

α1×N1=α2×N2，

(21)

當(dāng)指定α1=1時(shí)，α2=N1/N2.相應(yīng)的Wd則為

(22)

上式中WGTEM、WSATEM分別表示地面TEM和SATEM的數(shù)據(jù)加權(quán)矩陣.在聯(lián)合反演中，只需通過(guò)改變?chǔ)?，就可以調(diào)節(jié)不同數(shù)據(jù)子集的權(quán)重，確保不同數(shù)據(jù)子集權(quán)重的平衡.聯(lián)合反演中的靈敏度計(jì)算、法方程求解等流程與單獨(dú)反演完全相同.

4 理論模型研究

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

4.1 模型1：無(wú)地形模型

模型1電性結(jié)構(gòu)和觀測(cè)裝置如圖3所示.導(dǎo)線源長(zhǎng)1500 m，第1號(hào)源在X=-1500 m處，第2號(hào)源在X=1500 m處.設(shè)置6個(gè)大小均為500 m×500 m的大回線源，僅在回線圈內(nèi)觀測(cè)，采用兩種觀測(cè)點(diǎn)距，50 m點(diǎn)距的數(shù)據(jù)量為6×81×25=12150個(gè)，200 m點(diǎn)距的數(shù)據(jù)量為6×9×25=1350個(gè).SATEM測(cè)區(qū)范圍在X:-700～700 m，Y:-450～450 m之間，點(diǎn)距為50 m，數(shù)據(jù)量為551×25=13775個(gè).淺部?jī)蓚€(gè)異常體大小為300 m×600 m×60 m，中心點(diǎn)深度為60 m，電導(dǎo)率分別為0.1 S·m-1和0.001 S·m-1；深部異常體大小為400 m×600 m×200 m，中心點(diǎn)深度為350 m，電導(dǎo)率為0.2 S·m-1；背景電導(dǎo)率為0.01 S·m-1.

圖3 模型1的切片圖(a) Z=60 m切片，紅色線段、白色線框分別為SATEM發(fā)射源、地面TEM回線框發(fā)射源，黑色、白色測(cè)點(diǎn)的距離分別為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.

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

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

圖4為地面TEM單獨(dú)反演的結(jié)果.淺層兩個(gè)薄板異常體的位置和電導(dǎo)率恢復(fù)較好，但深部異常體的反演結(jié)果類似一個(gè)“環(huán)狀”結(jié)構(gòu)，只能恢復(fù)異常體的邊界，在異常體中心形成了一個(gè)空洞.隨著測(cè)點(diǎn)數(shù)量的減少，反演結(jié)果分辨率降低，假異常增加.對(duì)于實(shí)際探測(cè)，難以進(jìn)行合理的地質(zhì)解釋.

圖4 模型1的地面TEM反演結(jié)果(a,b) 50 m點(diǎn)距數(shù)據(jù)反演結(jié)果; (c,d) 200 m點(diǎn)距數(shù)據(jù)反演結(jié)果； (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.

為了測(cè)試SATEM裝置的探測(cè)能力，對(duì)兩個(gè)不同發(fā)射源的模擬數(shù)據(jù)分別進(jìn)行反演.網(wǎng)格、初始模型和其他參數(shù)與上相同，每次反演迭代約4 h.反演收斂情況如圖7b所示，最終反演結(jié)果的數(shù)據(jù)歸一化擬合差在5%左右.因?yàn)樵谳^小的偏移距內(nèi)進(jìn)行觀測(cè)，衰減曲線會(huì)出現(xiàn)“變號(hào)”，數(shù)據(jù)擬合差稍大于噪聲水平.

SATEM單獨(dú)反演結(jié)果如圖5所示.兩個(gè)反演結(jié)果均恢復(fù)了大致的地下電性分布情況.但與地面TEM裝置反演結(jié)果相比，SATEM對(duì)異常體邊界的分辨能力較差.SATEM單獨(dú)反演的結(jié)果也存在明顯的假異常，且在不同發(fā)射源的情況下，反演結(jié)果的假異常不同，這給數(shù)據(jù)解釋帶來(lái)巨大挑戰(zhàn).

圖5 模型1的SATEM反演結(jié)果(a,b) 第1號(hào)源數(shù)據(jù)反演結(jié)果； (c,d) 第2號(hào)源數(shù)據(jù)反演結(jié)果； (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.

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

圖6為SATEM與地面TEM聯(lián)合反演結(jié)果.聯(lián)合反演對(duì)三個(gè)異常體的邊界和電導(dǎo)率均得到有效重建.地面TEM單獨(dú)反演的深部異常體的“環(huán)狀”結(jié)構(gòu)以及SATEM單獨(dú)反演的假異常，在聯(lián)合反演中均得到有效壓制.尤其采用50 m點(diǎn)距的地面TEM數(shù)據(jù)時(shí)，聯(lián)合反演結(jié)果基本上沒(méi)有明顯假異常.當(dāng)?shù)孛鎀EM數(shù)據(jù)為200 m點(diǎn)距時(shí)，雖然數(shù)據(jù)量只有50 m點(diǎn)距的1/9，但聯(lián)合反演結(jié)果也明顯好于單獨(dú)反演的結(jié)果.

圖6 模型1的地面TEM與SATEM聯(lián)合反演結(jié)果(a,b) 地面TEM點(diǎn)距50 m數(shù)據(jù)與SATEM數(shù)據(jù)聯(lián)合反演； (c,d) 地面TEM點(diǎn)距200 m數(shù)據(jù)與SATEM數(shù)據(jù)聯(lián)合反演； (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聯(lián)合反演.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：帶地形模型

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

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

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

圖9 模型2在Z=100 m的模型與反演結(jié)果切片(a) 理論模型； (b) 地面TEM反演； (c) SATEM反演； (d) 聯(lián)合反演.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的模型與反演結(jié)果切片(a) 理論模型； (b) 地面TEM反演； (c) SATEM反演； (d) 聯(lián)合反演.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的模型與反演結(jié)果切片(a) 理論模型； (b) 地面TEM反演； (c) SATEM反演； (d) 聯(lián)合反演.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 結(jié)論

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