滿開峰， 殷長(zhǎng)春， 劉云鶴， 孫思源， 熊彬

1 中國(guó)科學(xué)院深圳先進(jìn)技術(shù)研究院， 深圳 518055 2 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 130026 3 深圳市自然資源和不動(dòng)產(chǎn)評(píng)估發(fā)展研究中心(深圳市地質(zhì)環(huán)境監(jiān)測(cè)中心)， 深圳 518040 4 中國(guó)自然資源航空物探遙感中心， 北京 100083 5 桂林理工大學(xué)地球科學(xué)學(xué)院， 桂林 541006

0 引言

時(shí)間域航空電磁法較傳統(tǒng)地面電磁法具有勘探效率高、成本低、無(wú)需地面人員接近等明顯技術(shù)優(yōu)勢(shì)，特別適合于地形復(fù)雜以及人員難以接近的地區(qū)(如高山、沙漠、湖泊沼澤及森林覆蓋區(qū))進(jìn)行大面積勘查，已在油氣、礦產(chǎn)、地下水和地?zé)豳Y源等領(lǐng)域獲得廣泛應(yīng)用 (殷長(zhǎng)春等，2015).基于實(shí)電阻率模型的時(shí)間域航空電磁中心回線(或重疊回線)裝置得到的觀測(cè)數(shù)據(jù)理論上是一個(gè)不變號(hào)的衰減曲線.然而，實(shí)際勘探中時(shí)間域航空電磁中心回線(或重疊回線)裝置晚期道響應(yīng)中常出現(xiàn)符號(hào)反轉(zhuǎn)現(xiàn)象(如圖1所示).眾多學(xué)者對(duì)這種現(xiàn)象開展了深入的研究，認(rèn)為這種符號(hào)反轉(zhuǎn)現(xiàn)象是由地下介質(zhì)存在激發(fā)極化效應(yīng)(Induced polarization，IP效應(yīng))引起的，并一致同意基于傳統(tǒng)實(shí)電阻率模型的成像和反演方法無(wú)法對(duì)此類符號(hào)反轉(zhuǎn)數(shù)據(jù)進(jìn)行有效解釋(Lee，1975, 1981；Raiche，1983；Lewis and Lee，1984；Hohmann and Newman，1990；Flores and Peralta-Ortega，2009).因此，對(duì)此類電磁數(shù)據(jù)進(jìn)行反演時(shí)需要考慮激電效應(yīng)的影響.為描述方便，下文提到的航空電磁系統(tǒng)均為中心回線(或重疊回線)裝置.

圖1 時(shí)間域航空電磁響應(yīng)(存在和不存在IP效應(yīng))Fig.1 Time-domain airborne EM responses (with IP and without IP)

電磁法勘探中激發(fā)極化參數(shù)反演的研究主要是近10年.Kratzer和Macnae(2012)提出了一種利用基函數(shù)擬合時(shí)間域航空電磁響應(yīng)的方法.基于該方法時(shí)間域航空電磁數(shù)據(jù)可被分解成電磁感應(yīng)和激發(fā)極化兩部分，進(jìn)而可以從激發(fā)極化數(shù)據(jù)中計(jì)算出視極化率.Chen等(2015)進(jìn)一步提出了電磁響應(yīng)分步提取技術(shù)：先利用基函數(shù)擬合提取出電磁感應(yīng)產(chǎn)生的基本響應(yīng)，然后根據(jù)基本響應(yīng)和極化響應(yīng)的關(guān)系再利用基函數(shù)擬合方法計(jì)算出極化響應(yīng)，最后計(jì)算視極化率.隨后，Macnae(2016)基于薄板模型設(shè)定基函數(shù)并提取出IP信息.此外，對(duì)于帶激電效應(yīng)的時(shí)間域電磁數(shù)據(jù)一維反演也取得了很好的進(jìn)展.Yu等(2013)基于SVD方法對(duì)含IP效應(yīng)的實(shí)測(cè)TEM數(shù)據(jù)進(jìn)行反演，并從實(shí)測(cè)數(shù)據(jù)中分離出IP信息.Viezzoli等(2017)使用一維橫向約束方法反演時(shí)間域航空電磁IP參數(shù)，結(jié)果表明當(dāng)時(shí)間常數(shù)在一定約束范圍內(nèi)，該方法對(duì)電阻率和極化率有較好的求解能力.Lin等(2019)在一維橫向約束基礎(chǔ)上對(duì)反演做進(jìn)一步改進(jìn)，通過(guò)采用模型空間轉(zhuǎn)換、在數(shù)據(jù)符號(hào)反轉(zhuǎn)處增大誤差、建立穩(wěn)定的初始模型、設(shè)置合理的正則化因子等技術(shù)，獲得了穩(wěn)定的IP參數(shù)反演結(jié)果.繆佳佳等(2020)嘗試?yán)肙ccam方法提取時(shí)間域航空電磁數(shù)據(jù)中心的激電信息，得到與真實(shí)模型較為接近的電阻率和極化率值.劉衛(wèi)強(qiáng)等(2021)利用粒子群算法對(duì)各向異性激發(fā)極化進(jìn)行反演，得到各地層在平行層理與垂直層理兩個(gè)方向上的電阻率、充電率及時(shí)間常數(shù).路俊濤等(2022)采用橫向約束反演同時(shí)計(jì)算激電參數(shù)及層厚，增加約束條件從而改善多解性問題.此外，對(duì)于三維反演也有部分學(xué)者進(jìn)行了嘗試.Kang和Oldenburg(2016)對(duì)時(shí)間域航空電磁激發(fā)極化數(shù)據(jù)進(jìn)行了三維反演.該方法將早期響應(yīng)作為背景電阻率模型響應(yīng)，并用總電磁響應(yīng)減去背景電阻率模型響應(yīng)近似得到極化響應(yīng)后，再反演得到極化率.然而，該方法的反演精度受限于對(duì)背景電導(dǎo)率模型的準(zhǔn)確估計(jì)和對(duì)極化響應(yīng)的有效校正，實(shí)用性受到較大限制.

時(shí)間域航空電磁激發(fā)極化效應(yīng)涉及電阻率，極化率，時(shí)間常數(shù)，頻率相關(guān)系數(shù)等多個(gè)參數(shù).由于激電參數(shù)之間靈敏度差異較大，時(shí)間域電磁法中激電參數(shù)反演存在嚴(yán)重的多解性，給反演結(jié)果的穩(wěn)定性帶來(lái)了極大的考驗(yàn)(繆佳佳等, 2020).另一方面，地下極化體的電阻率和極化率之間存在一定的相關(guān)關(guān)系，本文電磁數(shù)據(jù)反演中我們將這種相關(guān)性作為正則化項(xiàng)對(duì)反演進(jìn)行約束，從而減小激電參數(shù)反演的多解性.統(tǒng)計(jì)中的Pearson相關(guān)系數(shù)可以描述兩個(gè)變量的線性相關(guān)性，本文利用Pearson相關(guān)系數(shù)構(gòu)建電阻率和極化率的相關(guān)關(guān)系，并在反演目標(biāo)函數(shù)中引進(jìn)由其構(gòu)建的相關(guān)性約束項(xiàng)，從而實(shí)現(xiàn)基于Pearson相關(guān)性約束的時(shí)間域航空電磁數(shù)據(jù)反演.此外，實(shí)測(cè)數(shù)據(jù)反演中還需要考慮時(shí)間常數(shù)和頻率相關(guān)系數(shù)，由于兩者的靈敏度很小，與電阻率和極化率同步反演很難得到真實(shí)的參數(shù)值，本文采用深度學(xué)習(xí)方法預(yù)測(cè)時(shí)間常數(shù)和頻率相關(guān)系數(shù)，通過(guò)給定時(shí)間常數(shù)和頻率相關(guān)系數(shù)較小的約束范圍，進(jìn)一步減少反演多解性.我們將通過(guò)理論和實(shí)測(cè)數(shù)據(jù)驗(yàn)證本文反演算法的有效性.

1 三維正演計(jì)算

目前，用于描述巖礦石激電特性的模型有很多，使用較為廣泛的是Pelton等(1978)提出的Cole-Cole模型，其表達(dá)形式為：

(1)

其中，ρ0表示零頻電阻率；η為極化率，τ為時(shí)間常數(shù)，c為頻率相關(guān)系數(shù)，i為虛數(shù)單位，ω為角頻率.

假設(shè)時(shí)諧因子為eiω t，將總場(chǎng)分離為背景場(chǎng)和散射場(chǎng)，則由電磁場(chǎng)雙旋度方程可得散射場(chǎng)滿足如下偏微分方程(Newman and Alumbaugh, 1995)：

(2)

其中，μ0=4π×10-7H·m-1是自由空間磁導(dǎo)率，Ep是背景電場(chǎng)，Es是散射電場(chǎng)，σ(ω)是與頻率相關(guān)的復(fù)電導(dǎo)率.

利用基于交錯(cuò)網(wǎng)格的有限差分方法(Alumbaugh et al., 1996)對(duì)式(2)進(jìn)行離散，得到的線性方程組可表示為：

KEs=s,

(3)

其中，K為系數(shù)矩陣，s是右端項(xiàng).我們采用擬殘差法(QMR)(Siripunvaraporn et al.，2002)求解該線性方程組，得到電場(chǎng)值Es，并與一次場(chǎng)Ep相加得到總電場(chǎng)值E，進(jìn)而由法拉第感應(yīng)定律(式(4))可得到頻率域的磁感應(yīng)強(qiáng)度B.之后，我們利用殷長(zhǎng)春等(2013)提出的反余弦變換法并與發(fā)射波形褶積，即可得到任意發(fā)射波形條件下時(shí)間域航空電磁響應(yīng).式(4)為：

(4)

2 基于Pearson相關(guān)性約束的三維電磁反演方法

統(tǒng)計(jì)學(xué)中的Pearson相關(guān)系數(shù)主要用來(lái)衡量?jī)蓚€(gè)變量之間的線性相關(guān)程度.其值介于-1和1之間.當(dāng)兩個(gè)變量的線性相關(guān)增強(qiáng)時(shí)，相關(guān)系數(shù)趨于1或-1.當(dāng)兩個(gè)變量的線性相關(guān)越弱時(shí)，相關(guān)系數(shù)越趨于0；若相關(guān)系數(shù)等于0，則表明兩變量之間不相關(guān).若存在兩組長(zhǎng)度相同的變量X={x1,x2,…,xN}、Y={y1,y2,…,yN}，則它們的Pearson相關(guān)系數(shù)可寫為：

(5)

(6)

(7)

綜合式(6)、(7)可以得到如下基于子域Pearson相關(guān)系數(shù)的相關(guān)性約束項(xiàng)，即：

(8)

在使用梯度類反演算法對(duì)目標(biāo)函數(shù)進(jìn)行求解的過(guò)程中，需要計(jì)算目標(biāo)函數(shù)的梯度.我們采用復(fù)合函數(shù)求導(dǎo)法則計(jì)算φs(m1,m2)對(duì)地下物性參數(shù)的導(dǎo)數(shù)，可得：

(9)

則每個(gè)子域內(nèi)pi對(duì)物性參數(shù)m的導(dǎo)數(shù)可表示為：

(10)

(11)

(12)

基于上面討論的Pearson相關(guān)性約束，我們可以進(jìn)一步討論三維時(shí)間域航空電磁數(shù)據(jù)反演問題.為此，引入公式(8)的Pearson相關(guān)性約束項(xiàng)作為穩(wěn)定泛函，構(gòu)建如下正則化反演目標(biāo)函數(shù)：

(13)

其中，d為觀測(cè)數(shù)據(jù)，m0為初始模型參數(shù)，Cd為數(shù)據(jù)協(xié)方差矩陣，Cm為模型協(xié)方差矩陣.F(m)為正演算子，是一個(gè)非線性函數(shù)，需要對(duì)其進(jìn)行線性化處理.為此，對(duì)F(m)進(jìn)行泰勒展開，并忽略二階及以上的高階項(xiàng)可得：

F(m)=F(mk)+Jk(m-mk).

(14)

將式(14)代入式(13)，對(duì)m求導(dǎo)并化簡(jiǎn)后得到：

進(jìn)而，令gk=0得到如下反演求解方程：

(16)

式(16)可利用共軛梯度法進(jìn)行求解，得到模型更新量Δm.該向量指向目標(biāo)函數(shù)的下降方向.我們采用Haber(2014)提出的線性搜索方式得到步長(zhǎng)s，再利用公式mk+1=mk+sΔm獲得每次迭代后新的模型參數(shù).

3 激電參數(shù)預(yù)測(cè)

如上所述，激電參數(shù)之間靈敏度差異很大，特別是時(shí)間常數(shù)和頻率相關(guān)系數(shù)的靈敏度遠(yuǎn)遠(yuǎn)小于電阻率和極化率，因此同步直接反演四個(gè)參數(shù)難度很大.本文先利用深度學(xué)習(xí)對(duì)激電參數(shù)進(jìn)行估計(jì)，以此為基礎(chǔ)給定時(shí)間常數(shù)和頻率相關(guān)系數(shù)一個(gè)較小的約束范圍后再反演電阻率和極化率，從而解決四個(gè)參數(shù)同步反演的多解性問題.

4)對(duì)于海底光纜頻繁受損現(xiàn)狀，在加強(qiáng)與海事、港監(jiān)、邊防等政府部門單位聯(lián)絡(luò)協(xié)調(diào)的基礎(chǔ)上，利用現(xiàn)有的海纜監(jiān)測(cè)平臺(tái)增強(qiáng)實(shí)時(shí)監(jiān)控能力，對(duì)禁錨海域展開全天候巡視，一旦發(fā)現(xiàn)情況，立即警告驅(qū)離，防止船只錨損事件重復(fù)發(fā)生。及時(shí)修復(fù)受損海底光纜，對(duì)承載有繼電保護(hù)、圖像監(jiān)控和視頻會(huì)議等業(yè)務(wù)的海底光纜提前做好相關(guān)的應(yīng)急預(yù)案。

神經(jīng)網(wǎng)絡(luò)可以近似為如下非線性映射關(guān)系，即：

fθ∶X,H→Y，

(17)

其中，X是電磁響應(yīng)數(shù)據(jù)，Y是激電參數(shù)模型.θ代表神經(jīng)網(wǎng)絡(luò)的參數(shù)，主要是指權(quán)重和偏置.f通常是一個(gè)復(fù)雜的非線性函數(shù).我們可以通過(guò)神經(jīng)網(wǎng)絡(luò)來(lái)近似表示這種映射關(guān)系.

神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)的構(gòu)建通常依據(jù)經(jīng)驗(yàn)設(shè)置.由于長(zhǎng)短期記憶神經(jīng)網(wǎng)絡(luò)(LSTM)擅長(zhǎng)于處理時(shí)間序列的數(shù)據(jù)，我們選擇使用該神經(jīng)網(wǎng)絡(luò).同時(shí)，為了提取數(shù)據(jù)特征，我們使用CNN+LSTM的網(wǎng)絡(luò)模型，即在LSTM模型之前使用卷積神經(jīng)網(wǎng)絡(luò)CNN提取數(shù)據(jù)特征.在卷積層中，激活函數(shù)采用ReLU函數(shù)，在LSTM層中則采用tanh函數(shù).網(wǎng)絡(luò)收斂條件采用Early stopping方法.

我們可以通過(guò)訓(xùn)練集訓(xùn)練神經(jīng)網(wǎng)絡(luò)獲得式(17)的映射關(guān)系.訓(xùn)練集主要包含神經(jīng)網(wǎng)絡(luò)的輸入數(shù)據(jù)和輸出數(shù)據(jù)(標(biāo)簽).輸入數(shù)據(jù)為電磁響應(yīng)，而輸出數(shù)據(jù)對(duì)應(yīng)地下激電參數(shù).顯然，我們的目標(biāo)在于構(gòu)建訓(xùn)練集來(lái)訓(xùn)練一個(gè)神經(jīng)網(wǎng)絡(luò)，使得向該神經(jīng)網(wǎng)絡(luò)輸入一個(gè)電磁響應(yīng)后，輸出與該電磁相應(yīng)對(duì)應(yīng)的地下激電參數(shù)模型.通過(guò)在合理的參數(shù)范圍內(nèi)設(shè)計(jì)隨機(jī)激電參數(shù)，進(jìn)而使用一維正演軟件 (Man et al., 2020)生成相應(yīng)的電磁響應(yīng)數(shù)據(jù).在準(zhǔn)備訓(xùn)練集過(guò)程中，需要對(duì)樣本數(shù)據(jù)進(jìn)行預(yù)處理.激電參數(shù)設(shè)計(jì)為均勻半空間模型，其中零頻電阻率ρ0的取值范圍為1～10000 Ωm，并按照對(duì)數(shù)域設(shè)計(jì)隨機(jī)值.首先取[0,4]范圍內(nèi)的隨機(jī)數(shù)R1，即R1∈[0,4]，則電阻率的取值為10R1.極化率范圍設(shè)定為[0,1]，同樣按照對(duì)數(shù)域取值，即假設(shè)對(duì)數(shù)極化率取值為10R2，R2∈[-4,0].同理，時(shí)間常數(shù)取值為10R3，R3∈[-5,-2].頻率相關(guān)系數(shù)取值為10R4，R4∈[-2,0].

設(shè)計(jì)大小為200000的樣本數(shù)據(jù)，包括訓(xùn)練集195000個(gè)和測(cè)試集5000個(gè).在訓(xùn)練神經(jīng)網(wǎng)絡(luò)的過(guò)程中，我們使用minbatch gradient-descent方法對(duì)神經(jīng)網(wǎng)絡(luò)參數(shù)進(jìn)行更新，主要過(guò)程是從訓(xùn)練集中選取一小部分樣本子集(本文中batch的大小選取64)，通過(guò)梯度下降法使得該樣本構(gòu)成的代價(jià)函數(shù)值最小.本文使用的代價(jià)函數(shù)為Mean Absolute Error (MAE), 定義為預(yù)測(cè)值與期望值之間的平均絕對(duì)誤差，即：

(18)

其中，K是矢量y的長(zhǎng)度.式(18)可以通過(guò)反向傳播BP算法求解.

4 理論模型反演

為了驗(yàn)證本文算法的有效性，我們首先設(shè)計(jì)理論模型進(jìn)行反演測(cè)試.理論模型采用VTEM航空電磁系統(tǒng).發(fā)射波形如圖2所示.數(shù)值仿真時(shí)，我們假設(shè)發(fā)射線圈半徑為17.5 m，匝數(shù)為2匝，發(fā)射電流峰值1000 A，發(fā)射磁矩峰值為9621100 Am2.接收機(jī)位于發(fā)射線框的中心，發(fā)射線框和接收機(jī)距離地表高度均為30 m.

圖2 VTEM系統(tǒng)發(fā)射波形Fig.2 Transmitting wave of VTEM system

4.1 電阻率和極化率的反演

電阻率和極化率是反映地下電性結(jié)構(gòu)的重要信息，因此我們討論如何從觀測(cè)的電磁數(shù)據(jù)中反演得到地下介質(zhì)的電阻率和極化率.由Cole-Cole模型可知，地下極化體包含電阻率、極化率、時(shí)間常數(shù)和頻率相關(guān)系數(shù)四個(gè)參數(shù)，而得到的含有激電效應(yīng)的電磁響應(yīng)數(shù)據(jù)是由這四個(gè)參數(shù)共同作用的結(jié)果，因此實(shí)際上四個(gè)參數(shù)應(yīng)該同步參與反演.考慮到四個(gè)激電參數(shù)的靈敏度差別很大，極化率的靈敏度小于電阻率，而時(shí)間常數(shù)和頻率相關(guān)系數(shù)的靈敏度相對(duì)更小，因此本文先討論賦予時(shí)間常數(shù)和頻率相關(guān)系數(shù)真實(shí)值，并將其固定不變，再同步反演電阻率和極化率.由于電阻率和極化率靈敏度差異也較大，同步反演多解性也很嚴(yán)重，我們采用基于Pearson相關(guān)約束方法對(duì)電阻率和極化率施加相關(guān)性約束，以改善電阻率和極化率反演效果，減少多解性.

4.1.1 雙棱柱異常體模型

下面首先設(shè)計(jì)一個(gè)雙棱柱異常體模型，其xy和xz截面如圖3所示.兩個(gè)棱柱體的大小分別為150 m×180 m×46 m和150 m×180 m×55 m，埋深分別為30 m和37 m，兩者的水平間距為120 m.兩個(gè)棱柱體的電阻率均為10 Ωm，極化率為0.5，時(shí)間常數(shù)為0.001 s，頻率相關(guān)系數(shù)為0.5.圍巖的電阻率設(shè)計(jì)為100 Ωm，極化率為0.1，時(shí)間常數(shù)為0.0001 s，頻率相關(guān)系數(shù)為0.1.我們采用30個(gè)斷電時(shí)間道(off-time)計(jì)算時(shí)間域電磁響應(yīng)并加入5%的隨機(jī)噪聲，將最終的合成數(shù)據(jù)作為待反演的理論數(shù)據(jù).

圖3 雙棱柱體模型示意圖(a) z=54 m處xy剖面圖； (b) y=705 m處xz剖面圖.Fig.3 Schematic diagram of double prism model(a) xy section at z=54 m; (b) xz section at y=705 m.

首先對(duì)圖3中的雙棱柱異常體模型在固定時(shí)間常數(shù)和頻率相關(guān)系數(shù)的條件下進(jìn)行電阻率和極化率同步反演.圖4和圖5分別給出高斯-牛頓和基于Pearson相關(guān)約束的反演結(jié)果.從圖4中可以看出，高斯-牛頓反演得到的兩個(gè)異常體電阻率橫向分布范圍與真實(shí)模型基本一致，電阻率值與真實(shí)值較為接近，但兩個(gè)異常體的電阻率縱向分布范圍較大，兩個(gè)異常體的深部邊界均向下的延伸，左側(cè)淺部異常體的地面位置出現(xiàn)了部分低阻假異常.對(duì)于極化率，兩個(gè)異常體極化率的橫向分布范圍比真實(shí)模型稍大，而兩個(gè)異常體在深部均有較大的延伸，在深部?jī)蓚€(gè)異常體連接在一起，難以區(qū)分.從圖5中可以看出，基于Pearson相關(guān)約束反演得到的兩個(gè)異常體電阻率橫向分布范圍和電阻率值均與真實(shí)模型基本一致，電阻率縱向分布范圍也與真實(shí)模型較為接近，電阻率整體反演效果良好.對(duì)于極化率，兩個(gè)異常體極化率橫向分布范圍與真實(shí)模型基本一致，而縱向分布范圍稍大，有一定縱向延伸.

圖4 高斯-牛頓法雙棱柱模型電阻率和極化率反演結(jié)果Fig.4 Inversion results of resistivity and chargeability of double prism model by Gauss Newton method

圖5 Pearson相關(guān)約束雙棱柱模型電阻率和極化率反演結(jié)果Fig.5 Inversion results of resistivity and chargeability of double prism model by Pearson correlation constrained method

對(duì)比兩種反演方法得到的電阻率和極化率反演結(jié)果可以看出，兩種方法得到的電阻率異常形態(tài)和電阻率與真實(shí)模型較為接近，電阻率反演結(jié)果總體比極化率的反演結(jié)果好.這主要是由于電阻率的靈敏度比極化率的靈敏度大很多，從而使得同步反演兩個(gè)參數(shù)時(shí)電阻率較為靈敏，容易得到好的反演結(jié)果.對(duì)比高斯-牛頓和基于Pearson相關(guān)約束的電阻率和極化率的反演結(jié)果可以看出，基于Pearson相關(guān)約束的反演結(jié)果明顯好于高斯-牛頓的反演結(jié)果，尤其是極化率的反演，高斯-牛頓反演不能在深部區(qū)分兩個(gè)異常體，而基于Pearson相關(guān)約束反演得到的異常體在深部能夠得到區(qū)分.這主要是由于基于Pearson相關(guān)約束反演考慮電阻率和極化率的相關(guān)性，減少了極化率反演的多解性.

4.1.2 拱形異常體模型

為進(jìn)一步驗(yàn)證Pearson相關(guān)約束反演的效果，我們對(duì)圖6中的拱形異常體模型進(jìn)行反演試算.拱形異常體上頂面的水平方向大小為90 m×240 m，上頂面埋深為30 m，中部下低面的深度為60 m，左右兩端下底面深度為120 m.異常體沿y方向的延伸長(zhǎng)度均為240 m.拱形異常體的電阻率為10 Ωm，極化率為0.5，時(shí)間常數(shù)為0.001 s，頻率相關(guān)系數(shù)為0.5.圍巖的電阻率為100 Ωm，極化率為0.1，時(shí)間常數(shù)為0.0001 s，頻率相關(guān)系數(shù)為0.1.對(duì)拱形模型進(jìn)行正演計(jì)算，將生成的正演響應(yīng)加入5%的隨機(jī)噪聲，最終合成的數(shù)據(jù)作為待反演的理論響應(yīng)數(shù)據(jù).

圖6 拱形模型示意圖(a) z=54 m位置xy切片； (b) y=705 m位置xz切片.Fig.6 Schematic diagram of arch model(a) xy section at z=54 m; (b) xz section at y=705 m.

圖7和圖8分別給出高斯-牛頓和基于Pearson相關(guān)約束的電阻率和極化率反演結(jié)果，其中黑色實(shí)線框?yàn)檎鎸?shí)模型的位置.從圖7中可以看出，高斯-牛頓反演得到的拱形異常體電阻率橫向分布范圍與真實(shí)模型基本一致，而拱形異常體的電阻率縱向分布與真實(shí)模型有一定的差異，左右兩端下底面相對(duì)真實(shí)模型均呈現(xiàn)向下延伸，而拱形異常體上頂面上方出現(xiàn)低阻假異常.對(duì)于極化率，拱形異常體極化率的橫向分布范圍比真實(shí)模型大，中間和兩端的下底面相比真實(shí)模型均有較大的延伸，左右兩只連成一體.從圖8可以看出，基于Pearson相關(guān)約束反演得到的拱形異常體電阻率和極化率橫向分布范圍與真實(shí)模型相當(dāng)一致.特別是極化率呈現(xiàn)左右兩只分開的形態(tài)，更接近真實(shí)模型分布特征.由此得出結(jié)論，基于Pearson相關(guān)約束的反演效果明顯好于高斯-牛頓方法.

圖7 高斯-牛頓法拱形模型電阻率和極化率反演結(jié)果Fig.7 Inversion results of resistivity and chargeability of arch model by Gauss Newton method

圖8 Pearson相關(guān)約束拱形模型電阻率和極化率反演結(jié)果Fig.8 Inversion results of resistivity and chargeability of arch model by Pearson correlation constrained method

4.2 考慮四個(gè)激電參數(shù)的反演

為了實(shí)現(xiàn)四個(gè)激電參數(shù)的三維反演，本文提出深度學(xué)習(xí)估計(jì)激電參數(shù)和基于Pearson相關(guān)約束相結(jié)合的聯(lián)合反演策略.為驗(yàn)證其有效性，我們利用深度學(xué)習(xí)對(duì)圖6中拱形模型的時(shí)間常數(shù)和頻率相關(guān)系數(shù)進(jìn)行估算后，再進(jìn)行電阻率和極化率反演.從圖9給出的結(jié)果可以看出，利用該聯(lián)合反演策略，拱形異常體電阻率和極化率形態(tài)與真實(shí)模型基本一致，但縱向分布范圍有所增大.這是由于時(shí)間常數(shù)和頻率相關(guān)系數(shù)估算結(jié)果不精確造成的.將圖9中的反演結(jié)果與圖8進(jìn)行對(duì)比發(fā)現(xiàn)，利用深度學(xué)習(xí)估算時(shí)間常數(shù)和頻率相關(guān)系數(shù)后，電阻率和極化率反演結(jié)果與給定真實(shí)時(shí)間常數(shù)和頻率相關(guān)系數(shù)的電阻率和極化率反演效果相當(dāng)，僅在縱向分布上向下拓展.這種現(xiàn)象不難理解.事實(shí)上，由于機(jī)器學(xué)習(xí)無(wú)法對(duì)時(shí)間常數(shù)和頻率相關(guān)系數(shù)進(jìn)行精確估算，因此聯(lián)合反演無(wú)法獲得與使用真實(shí)模型參數(shù)時(shí)相同的反演結(jié)果.

圖9 利用深度學(xué)習(xí)估計(jì)時(shí)間常數(shù)和頻率相關(guān)系數(shù)后拱形異常體電阻率和極化率反演結(jié)果Fig.9 Inversion results of resistivity and chargeability of arch model after estimating time constant and frequency exponent by depth learning method

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

為了進(jìn)一步測(cè)試本文反演策略的有效性，我們對(duì)澳大利亞馬斯格雷夫省北部實(shí)測(cè)數(shù)據(jù)進(jìn)行反演.航空電磁采用SkyTEM系統(tǒng).發(fā)射波形如圖10所示，on-time為5 ms.發(fā)射線圈面積為337 m2，匝數(shù)12匝，最大發(fā)射磁矩為473000 Am2，發(fā)射基頻25 Hz.發(fā)射和接收線圈高度約為45 m.圖11中紅線給出航空電磁測(cè)線位置.本文中我們選取測(cè)線L502201、L502301、L502501、L502601、L502701中Northing 702739.4～702939.4范圍內(nèi)的實(shí)測(cè)數(shù)據(jù)作為測(cè)試區(qū)，測(cè)試區(qū)位于圖12中藍(lán)色矩形位置，每條測(cè)線點(diǎn)距為100 m，測(cè)點(diǎn)間距為250 m，包含有25個(gè)off-time時(shí)間道.為簡(jiǎn)單起見，下面討論中將它們簡(jiǎn)單標(biāo)記為L(zhǎng)1、L2、L3、L4和L5.

圖10 SkyTEM系統(tǒng)發(fā)射波形Fig.10 Transmitting wave of SkyTEM system

圖11 澳大利亞馬斯格雷夫省航空電磁測(cè)區(qū)示意圖圖片來(lái)自:https:∥ ecat.ga.gov.au/geonetwork/srv/chi/catalog.search#/metadata/110284.Fig.11 Schematic diagram of airborne electromagnetic survey area in Masgrave Province, Australia Picture from https:∥ecat.ga.gov.au/geonetwork/srv/chi/catalog.search#/metadata/110284.

圖12 澳大利亞馬斯格雷夫省實(shí)測(cè)數(shù)據(jù)反演結(jié)果(a) 不考慮IP效應(yīng)反演的電阻率結(jié)果(不包含負(fù)響應(yīng)的19個(gè)時(shí)間道數(shù)據(jù))； (b) 不考慮IP效應(yīng)反演的電阻率結(jié)果(包含負(fù)響應(yīng)的25個(gè)時(shí)間道數(shù)據(jù))； (c) 考慮IP效應(yīng)反演的電阻率結(jié)果(包含負(fù)響應(yīng)的25個(gè)時(shí)間道數(shù)據(jù)); (d) 考慮IP效應(yīng)反演的極化率結(jié)果(包含負(fù)響應(yīng)的25個(gè)時(shí)間道數(shù)據(jù)).Fig.12 Inversion results of measured data in Musgrave Province, Australia (a) Resistivity inversion results without considering IP effect (excluding 19 time channel data with negative response); (b) Resistivity inversion results without considering IP effect (including 25 time channel data with negative response); (c) Resistivity inversion results considering IP effect (including 25 time channel data with negative response); (d) Chargeability inversion results considering IP effect (including 25 time channel data with negative response).

下面分別對(duì)不包含負(fù)響應(yīng)的數(shù)據(jù)(19個(gè)時(shí)間道)進(jìn)行不考慮IP效應(yīng)的反演，對(duì)包含負(fù)響應(yīng)的數(shù)據(jù)(取25個(gè)時(shí)間道)進(jìn)行不考慮IP效應(yīng)的反演和考慮IP效應(yīng)的反演，反演結(jié)果如圖12所示.圖12a給出對(duì)19個(gè)時(shí)間道的正響應(yīng)數(shù)據(jù)進(jìn)行不考慮IP效應(yīng)的反演結(jié)果，圖中沿Y方向依次為測(cè)線L1～L5的電阻率反演剖面.從圖可以看出，測(cè)線L2、L3、L4地下大致分為兩層，淺部為低阻層，電阻率在4～10 Ωm，厚度5～100 m，深部為相對(duì)高阻層，電阻率為50～200 Ωm.測(cè)線L1中X=1000～1500 m處淺部低阻層消失，而地下存在一個(gè)高阻異常體.測(cè)線L5地下則表現(xiàn)為三層電性分布特征.圖12b給出對(duì)包含負(fù)響應(yīng)的25時(shí)間道數(shù)據(jù)在不考慮IP效應(yīng)條件下的反演結(jié)果.與圖12a相比，L4和L5測(cè)線的低阻層變厚，而測(cè)線L1、L2、L3均出現(xiàn)了電阻率為500 Ωm的高阻體.圖12c給出對(duì)包含負(fù)響應(yīng)的25個(gè)時(shí)間道數(shù)據(jù)，在考慮IP效應(yīng)條件下的電阻率反演結(jié)果.從圖可以看出，相比于圖12a、b，測(cè)線L1、L2、L3、L4中低阻層更加連續(xù)，界面更加清晰，而測(cè)線L5的反演結(jié)果分為三層，與圖12a反演結(jié)果一致.圖12d給出對(duì)包含負(fù)響應(yīng)的25個(gè)時(shí)間道數(shù)據(jù)進(jìn)行考慮IP效應(yīng)的極化率反演結(jié)果.從圖可以看出，極化率表現(xiàn)為與電阻率分布特征類似的層狀結(jié)構(gòu)，高極化與低電阻率存在對(duì)應(yīng)關(guān)系.綜合圖12c、d可以看出，測(cè)線L1、L2和L3均表現(xiàn)為淺部低阻高極化層，而深部為高阻低極化基底.測(cè)線L4和L5電性分為三層，中間層為低阻高極化層.

圖13給出不考慮IP效應(yīng)19個(gè)時(shí)間道數(shù)據(jù)反演、不考慮IP效應(yīng)25個(gè)時(shí)間道數(shù)據(jù)反演和考慮IP效應(yīng)25個(gè)時(shí)間道數(shù)據(jù)反演結(jié)果與實(shí)測(cè)數(shù)據(jù)的對(duì)比結(jié)果(以測(cè)線L3為例).從圖可以看出，對(duì)不考慮IP效應(yīng)19個(gè)時(shí)間道數(shù)據(jù)反演，理論和實(shí)測(cè)數(shù)據(jù)在19個(gè)時(shí)間上均擬合較好，而對(duì)不考慮IP效應(yīng)25個(gè)時(shí)間道數(shù)據(jù)反演，理論和實(shí)測(cè)數(shù)據(jù)在x=1800～2000 m區(qū)間早期道沒能很好地?cái)M合，而晚期20～25時(shí)間道的數(shù)據(jù)無(wú)法得到擬合.對(duì)于考慮IP效應(yīng)25個(gè)時(shí)間道數(shù)據(jù)反演，理論實(shí)測(cè)數(shù)據(jù)在早期道有較好的擬合，在包含負(fù)響應(yīng)的晚期道也有一定程度的擬合.綜上，對(duì)正響應(yīng)數(shù)據(jù)進(jìn)行不考慮IP效應(yīng)的反演和對(duì)包含負(fù)響應(yīng)數(shù)據(jù)進(jìn)行不考慮IP效應(yīng)反演均無(wú)法很好地?cái)M合全時(shí)間道數(shù)據(jù)，而考慮IP效應(yīng)的反演能夠相對(duì)較好地?cái)M合全時(shí)間道數(shù)據(jù)，說(shuō)明考慮IP效應(yīng)進(jìn)行電磁數(shù)據(jù)反演的必要性.

圖13 (a) 不考慮IP效應(yīng)19個(gè)時(shí)間道數(shù)據(jù)反演結(jié)果與實(shí)測(cè)數(shù)據(jù)對(duì)比； (b) 不考慮IP效應(yīng)25個(gè)時(shí)間道數(shù)據(jù)反演結(jié)果與實(shí)測(cè)數(shù)據(jù)對(duì)比； (c) 考慮IP效應(yīng)25個(gè)時(shí)間道數(shù)據(jù)反演結(jié)果與實(shí)測(cè)數(shù)據(jù)對(duì)比Fig.13 (a) Comparison of 19 time channels measured data and the results without considering IP effect; (b) Comparison of 25 time channels measured data and the results without considering IP effect; (c) Comparison of 25 time channels measured data and the results with considering IP effect

6 結(jié)論

本文將電阻率和極化率的Pearson相關(guān)性約束和深度學(xué)習(xí)預(yù)測(cè)激電參數(shù)成功應(yīng)用于帶激電效應(yīng)的時(shí)間域航空電磁數(shù)據(jù)反演中，通過(guò)對(duì)理論和實(shí)測(cè)數(shù)據(jù)反演得出如下結(jié)論：

(1)高斯-牛頓法反演得到的極化率分布較為發(fā)散，而基于Pearson相關(guān)約束反演，由于考慮電阻率和極化率的相關(guān)性特征，得到的異常體極化率橫向和縱向分布相對(duì)聚集，相比沒有施加約束的高斯-牛頓法的反演結(jié)果更接近真實(shí)模型.電阻率反演結(jié)果也有一定的改善.

(2)本文提出的時(shí)間域航空電磁多激電參數(shù)聯(lián)合反演是較為成功的反演策略.利用深度學(xué)習(xí)估計(jì)時(shí)間常數(shù)和頻率相關(guān)系數(shù)后，基于估計(jì)的參數(shù)對(duì)電阻率和極化率進(jìn)行相關(guān)約束反演，可以改善極化率反演效果.

(3)對(duì)不包含負(fù)響應(yīng)數(shù)據(jù)進(jìn)行不考慮IP效應(yīng)的反演，對(duì)包含負(fù)響應(yīng)數(shù)據(jù)進(jìn)行不考慮IP效應(yīng)的反演和考慮IP效應(yīng)的反演結(jié)果表明，考慮IP效應(yīng)反演的低阻層更加連續(xù)，界面更加清晰，數(shù)據(jù)擬合程度高，說(shuō)明反演中考慮IP效應(yīng)的必要性.

致謝感謝吉林大學(xué)電磁研究團(tuán)隊(duì)成員在本文準(zhǔn)備過(guò)程中給予的幫助.感謝Geoscience Australia提供的實(shí)測(cè)航空電磁數(shù)據(jù) (https:∥ecat.ga.gov.au/geonetwork/srv/chi/catalog. search#/metadata/110284).最后，對(duì)兩位審稿人提出的建設(shè)性意見表示衷心感謝！