大地電磁阻抗計(jì)算方法進(jìn)展

馬海榮 ,王世彪，郭榮文,柳建新

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙 410083;2.中南大學(xué) 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410083；3.廣東省有色地質(zhì)環(huán)境中心，廣東 廣州 510080)

?

大地電磁阻抗計(jì)算方法進(jìn)展

馬海榮1,2,王世彪3，郭榮文1,2,柳建新1,2

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙 410083;2.中南大學(xué) 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410083；3.廣東省有色地質(zhì)環(huán)境中心，廣東 廣州 510080)

大地電磁法主要依據(jù)視電阻率及相位曲線(xiàn)的形態(tài)特征來(lái)分析大地電性分布，因此阻抗張量的計(jì)算十分重要。在實(shí)際探測(cè)中，噪聲的存在使得計(jì)算得到的阻抗張量往往存在一定的偏差，干擾了數(shù)據(jù)解釋的準(zhǔn)確性。本文概述了大地電磁阻抗張量計(jì)算方法的發(fā)展，首先介紹了現(xiàn)在常用的幾種大地電磁阻抗計(jì)算方法及其不足之處，然后介紹了如何利用兩種Jackknife方法來(lái)判斷阻抗估計(jì)的穩(wěn)定性。

大地電磁法；阻抗張量；Jackknife

1　引　言

大地電磁測(cè)深法(Magnetotelluric, MT)誕生于20世紀(jì)50年代，隨著技術(shù)及儀器的發(fā)展，已成為油氣田、煤炭普查、巖石圈結(jié)構(gòu)探測(cè)、地下水與地?zé)豳Y源尋找、工程勘探等領(lǐng)域的一種重要手段，是一種應(yīng)用較廣的地球物理方法。在實(shí)際探測(cè)中，大地電磁往往受到較強(qiáng)的電磁噪聲干擾，且現(xiàn)在的勘探環(huán)境相對(duì)復(fù)雜，勘探難度相對(duì)加大，如何獲得質(zhì)量較好的原始數(shù)據(jù)及對(duì)原始數(shù)據(jù)進(jìn)行有效的處理，對(duì)后續(xù)反演等工作有著積極的作用。在工程物探中，大地電磁阻抗計(jì)算采用傳統(tǒng)的最小二乘法，但是因?yàn)榇蟮仉姶艤y(cè)深應(yīng)用的是天然場(chǎng)源，易受到外界噪聲的干擾，為了提高信噪比，獲得無(wú)偏的阻抗估計(jì)，提出了兩種改進(jìn)方法：一是Robust穩(wěn)健阻抗估計(jì)法[1-3]，它的核心是減少較大異常值對(duì)整體數(shù)據(jù)的影響；二是遠(yuǎn)參考法[4]，它的核心是設(shè)置參考點(diǎn)，消除輸入場(chǎng)(磁場(chǎng))的噪聲帶來(lái)的偏移。隨后，Chave等評(píng)述了大地電磁阻抗估計(jì)中的許多細(xì)節(jié)，并提出了Jackknife估計(jì)阻抗誤差的方法[5]。

2　大地電磁阻抗計(jì)算方法

大地電磁阻抗反映了地下電性結(jié)構(gòu)特征，它的計(jì)算原理為：大地電磁野外采集數(shù)據(jù)為電場(chǎng)與磁場(chǎng)x、y方向的分量即Ex、Ey、Hx、Hy的時(shí)間序列，首先將它們進(jìn)行Fourier變換，得到頻率域數(shù)據(jù)，然后用頻率域數(shù)據(jù)計(jì)算大地電磁阻抗Z。下面分別介紹幾種常用的阻抗張量計(jì)算方法。

2.1最小二乘處理

最小二乘法(Least Square Estimator)是最早的大地電磁阻抗張量的計(jì)算方法，它需滿(mǎn)足3個(gè)假設(shè)前提：①誤差方差相等；②誤差為高斯分布；③信號(hào)與噪聲不相關(guān)[6]。它是最傳統(tǒng)的阻抗估計(jì)方法[7]，滿(mǎn)足公式

(1)

(2)

應(yīng)用最小二乘法，使其殘差和最小，即

(3)

則最小二乘解為

(4)

其中，*表示轉(zhuǎn)置矩陣的共軛矩陣。

一般得到的阻抗形式為

Zxx=

(5)

Zxy=

(6)

Zyx=

(7)

Zyy=

(8)

判斷阻抗值是否穩(wěn)定及可靠，可通過(guò)計(jì)算其協(xié)方差來(lái)判斷。最小二乘法協(xié)方差可由下式計(jì)算：

(9)

(10)

2.2Robust 處理

Robust是近幾十年來(lái)提出的一種統(tǒng)計(jì)方法，其實(shí)質(zhì)是在實(shí)際情況與理想模型有微小偏離時(shí)，對(duì)這種偏離不敏感的一種統(tǒng)計(jì)方法。傳統(tǒng)方法假定數(shù)據(jù)服從高斯分布，但實(shí)際上往往不是，原來(lái)的計(jì)算方法就會(huì)造成阻抗張量計(jì)算的分散或偏離。Robust穩(wěn)健阻抗估計(jì)是一種主流的提高阻抗估計(jì)質(zhì)量的方法。穩(wěn)健估計(jì)的成功依賴(lài)于整個(gè)數(shù)據(jù)段中性態(tài)良好的數(shù)據(jù)占主要地位，當(dāng)噪聲貫穿大部分甚至整個(gè)觀測(cè)時(shí)段時(shí)，該方法就失效了[8-12]。

2.2.1M回歸估計(jì)

M回歸估計(jì)[13]是一種加權(quán)最小二乘法，主要研究不同數(shù)據(jù)長(zhǎng)度的Fourier變換，它根據(jù)觀測(cè)誤差和剩余功率譜的大小對(duì)數(shù)據(jù)進(jìn)行加權(quán)，注重未受干擾的數(shù)據(jù)，降低“飛點(diǎn)”的權(quán)，使之對(duì)阻抗函數(shù)的估算影響最小。該方法的關(guān)鍵在于尺度估計(jì)和權(quán)系數(shù)的處理。

Huber提出了最大似然的估算方法，即M回歸估計(jì)法，其基本思路就是不允許少量異常數(shù)據(jù)在阻抗函數(shù)的估算中起主要控制作用，因此對(duì)于這類(lèi)線(xiàn)性方程，就是使下式最?。?/p>

(11)

式中：ρ(r)稱(chēng)為損失函數(shù)；x為誤差；β為尺度函數(shù)。根據(jù)Huber的定義，M回歸的損失函數(shù)為

(12)

式中：一般取k=1.5,k為調(diào)整量。當(dāng)r小于1.5時(shí)，認(rèn)為該點(diǎn)為正常值，損失函數(shù)與最小二乘一致；當(dāng)r大于1.5時(shí)，用不同的損失函數(shù)來(lái)降低對(duì)應(yīng)數(shù)據(jù)的權(quán)重。

定義影響函數(shù)ψ(r)為損失函數(shù)的一階偏導(dǎo)數(shù)，即ψ(r)=ρ′(r)，令w(r)=ψ(r)/r，滿(mǎn)足

(13)

這里，W為權(quán)重的對(duì)角矩陣。上式可用迭代方法求解，是阻抗Z的Robust回歸迭代估計(jì)。

Robust法協(xié)方差計(jì)算相對(duì)復(fù)雜，因?yàn)榇藭r(shí)的傳遞函數(shù)估計(jì)值不能簡(jiǎn)單的用線(xiàn)性公式來(lái)準(zhǔn)確表示，協(xié)方差可由下式計(jì)算：

(14)

(15)

(16)

這里，wi表示權(quán)重。

2.2.2有界影響估計(jì)

從2.2.1的論述可知，M回歸估計(jì)降低奇異值的權(quán)重均是在殘差r的基礎(chǔ)上討論的，這表明僅在電道存在噪聲時(shí)，M回歸估計(jì)可以給出穩(wěn)定的結(jié)果。在實(shí)際工作中，盡管磁道信號(hào)抗干擾能力強(qiáng)于電道，但磁道存在噪聲的情況非常普遍，在這種情況下，M回歸估計(jì)效果變差。因此人們需要一種在電道和磁道均含噪聲的情況下仍然穩(wěn)健的估計(jì)方法。為了解決這一問(wèn)題，提出了有界影響估計(jì)[14]。

在統(tǒng)計(jì)學(xué)中，自變量(磁道)中的噪聲稱(chēng)為杠桿點(diǎn)(Leveragepoint)。在回歸理論中，帽子矩陣 (Hatmatrix)是一個(gè)重要的輔助量，用于判斷自變量中是否存在奇異值。

HAT=H(H*H)-1H*

帽子矩陣H是一個(gè)投影矩陣，且由公式可以看出，它僅與輸入道磁道有關(guān)，而與輸出道電道無(wú)關(guān)。帽子矩陣滿(mǎn)足如下性質(zhì)：

1)HAT為對(duì)稱(chēng)矩陣，冪等矩陣；

2)HAT的特征值或者為0或者為1，非0個(gè)數(shù)等于矩陣的秩；

3)HAT對(duì)角線(xiàn)元素hii滿(mǎn)足0

4)HAT對(duì)角線(xiàn)元素hii的期望值為p/N，其中p為待估計(jì)參數(shù)個(gè)數(shù)，在大地電磁阻抗估計(jì)中為2，N為參與估計(jì)的數(shù)據(jù)段數(shù)，在大地電磁阻抗估計(jì)中通常為分段處理的段數(shù)n，通常N?P;

5)—般而言，當(dāng)HAT對(duì)角線(xiàn)元素hii滿(mǎn)足ii?2p/N時(shí)，認(rèn)為對(duì)應(yīng)的第i道輸入存在奇異值。類(lèi)似于M回歸估計(jì)迭代格式，有界影響估計(jì)阻抗的迭代格式為

(17)

其中，v與M回歸估計(jì)相同，可釆用Huber權(quán)重函數(shù)或者Thomson權(quán)重函數(shù)；W為用于降低杠桿點(diǎn)的矩陣，其對(duì)角線(xiàn)上元素為

(18)

有界影響估計(jì)通過(guò)矩陣v降低電道異常數(shù)據(jù)在阻抗估計(jì)中的權(quán)重，通過(guò)矩陣W降低磁道噪聲數(shù)據(jù)在回歸問(wèn)題中的權(quán)重。

2.2.3重復(fù)中值估計(jì)

重復(fù)中值估計(jì)最早由Siegel提出，SmimovPai將它應(yīng)用到了大地電磁阻抗估計(jì)中。重復(fù)中值估計(jì)的思路與前面幾種估計(jì)方法不盡相同，巧妙地利用了中值遠(yuǎn)比均值穩(wěn)健這一特性，它是一種簡(jiǎn)潔高效的估計(jì)方法。

對(duì)于參與估計(jì)的第i、j組數(shù)據(jù)，大地電磁阻抗估計(jì)模型可寫(xiě)成如下形式：

(19)

(20)

2.3遠(yuǎn)參考處理

遠(yuǎn)參考大地電磁測(cè)量，是將一遠(yuǎn)參考點(diǎn)處的磁信號(hào)作為測(cè)點(diǎn)處的磁分量來(lái)進(jìn)行阻抗張量的估算，求取阻抗張量Z的估計(jì)值Zr。和傳統(tǒng)的Z估計(jì)值不同，只要參考場(chǎng)與電、磁道噪音不相關(guān)的話(huà)，Zr不受任一場(chǎng)中噪音功率偏倚的影響[15,16]。

實(shí)際觀測(cè)值是實(shí)際信號(hào)和噪聲干擾之和，由以下式子表示：

Ex=Exs+Exn

(21)

式中，Ex表示電場(chǎng)實(shí)際觀測(cè)值；Exn表示噪聲干擾；Exs表示真實(shí)信號(hào)。而大地電磁場(chǎng)的線(xiàn)性關(guān)系只對(duì)信號(hào)項(xiàng)才成立，即

通常，對(duì)于二維介質(zhì)，遠(yuǎn)參考法得到的阻抗張量表達(dá)式為

Zxx=

(22)

Zxy=

(23)

Zyx=

(24)

Zyy=

(25)

其中

(26)

式中，下標(biāo)r表示參考道測(cè)得的磁場(chǎng)分量。

遠(yuǎn)參考法中，協(xié)方差可由下式計(jì)算：

(27)

(28)

×(R*H)-1(R*R)(H*R)-1

(29)

2.4其他

除此之外，多道數(shù)據(jù)的利用也是學(xué)者們努力的方向。常見(jiàn)的一種思路是利用測(cè)站間的電磁場(chǎng)分量相互關(guān)系對(duì)含噪數(shù)據(jù)進(jìn)行預(yù)測(cè)，如大地電流-大地電磁(Telluric-Magnetotelluric, T-MT)方案[17]，“偽遠(yuǎn)參考”和組合站間轉(zhuǎn)換函數(shù)[18]等一批基于測(cè)站間“轉(zhuǎn)換函數(shù)”的采集處理方案，這些方案不僅改善了MT的數(shù)據(jù)質(zhì)量，而且拓展了數(shù)據(jù)采集的思路。然而，參考道的數(shù)據(jù)質(zhì)量決定了此類(lèi)方法的適用范圍。另一些學(xué)者提出了算法方面的改進(jìn)，如Kappler[19]提出了IARWR(Intersite Activity Ratio Wiener-filter Replacement)時(shí)間域去噪方法，王輝等[20]提出了利用多道同步依賴(lài)關(guān)系進(jìn)行阻抗估計(jì)，Cui等[21]提出引入獨(dú)立分量分析(Independent Component Analysis，ICA)的方法，Varentsov等[22]提出利用多道遠(yuǎn)參考聯(lián)合使用的方法。另外，Egbert等系統(tǒng)地闡述了利用陣列數(shù)據(jù)進(jìn)行信噪分離的阻抗估計(jì)方法[2,23,24]。隨后，在其基礎(chǔ)上，Smirnov等[25]發(fā)展了陣列數(shù)據(jù)的主成分分析穩(wěn)健估計(jì)方法(Principal Component Analysis，PCA)。

3　阻抗估計(jì)穩(wěn)定性

用標(biāo)準(zhǔn)方差來(lái)估計(jì)傳遞函數(shù)誤差棒的穩(wěn)定性常常是不理想的，因?yàn)閭鹘y(tǒng)的估計(jì)方法都是參數(shù)估計(jì)法，如果統(tǒng)計(jì)參數(shù)呈非正態(tài)分布，則會(huì)使估計(jì)值有偏差，不能給出正確的區(qū)間估計(jì)。而Jackknife法(刀切法)是一種非參數(shù)估計(jì)法，它不受參數(shù)統(tǒng)計(jì)分布類(lèi)型的限制，能夠減小估計(jì)偏差，給出近似的置信區(qū)間[26]。

(30)

計(jì)算偽值的算數(shù)平均值，得到Z的刀切法估計(jì)量

(31)

則刀切法誤差協(xié)方差計(jì)算公式為

(32)

當(dāng)測(cè)量數(shù)據(jù)為非線(xiàn)性時(shí)，使用刀切法來(lái)計(jì)算誤差棒更加準(zhǔn)確。刀切法提供了一個(gè)相對(duì)簡(jiǎn)單的方式來(lái)計(jì)算估計(jì)值的誤差棒。它也可用于張量分解的誤差分析。

3.1固定權(quán)重刀切法

(33)

權(quán)重hi由帽子矩陣定義。

以遠(yuǎn)參考法為例，有

(34)

3.2刪除子集刀切法

刀切法能夠計(jì)算復(fù)雜非線(xiàn)性預(yù)測(cè)值的誤差棒，為了更好地運(yùn)用這一點(diǎn)，在計(jì)算每一個(gè)偽值時(shí)，刪除某一固定區(qū)域的數(shù)值，每個(gè)子集上都進(jìn)行完全Robust計(jì)算。偽值、傳遞函數(shù)、方差可由公式(30)～(32)計(jì)算得到，其中I為子集數(shù)，i表示每個(gè)子集中的頻點(diǎn)數(shù)[7]。

3.3不足之處

對(duì)于Jackknife計(jì)算，有以下總結(jié)：

1)固定權(quán)重刀切法計(jì)算得出的誤差估計(jì)有類(lèi)似的計(jì)算公式，且允許磁道的信號(hào)與噪聲有一定的相關(guān)性，如果相關(guān)性超出范圍，則會(huì)計(jì)算出錯(cuò)誤的誤差估計(jì)值。

2)刪除子集刀切法在偏移帶的計(jì)算效果較好，特別是刪除了較大偏移的數(shù)據(jù)。但是不足的是計(jì)算量相對(duì)較大。

4　結(jié)　語(yǔ)

大地電磁法中，阻抗計(jì)算尤為重要。傳統(tǒng)最小二乘法有較多實(shí)際意義上的限制，Robust法或遠(yuǎn)參考法也由于具體情況而難分好壞。Jackknife刀切法可以利用非參數(shù)方法計(jì)算阻抗張量，極大地增加了傳遞函數(shù)的穩(wěn)定性及可靠性。

此外，由于天然電磁場(chǎng)強(qiáng)度低，大地電磁數(shù)據(jù)的信噪比往往不高，為了獲得有用信號(hào)，可采用一些新的信噪分離技術(shù)，如數(shù)學(xué)形態(tài)濾波技術(shù)[28]、Hilbert-Huang變換等，張量阻抗也從引入單參考道到多道數(shù)據(jù)綜合利用。同時(shí)，局部畸變?nèi)匀皇乾F(xiàn)階段的難題，目前的主要解決方法還是阻抗張量分解、空間濾波等方法。

[1]Larsen J C. Electromagnetic response functions from interrupted and noisy data[J]. Journal of Geomagnetism and Geoelectricity,1980,32(supplement 1): SI89-SI103.

[2]Egbert G D, Booker J R. Robust estimation of geomagnetic transfer functions[J]. Geophysical Journal of the Royal Astronomical Society ,1986,87(1): 173-194.

[3]Chave A D, Thomson D J, Ander M. On the robustestimation of power spectra, coherences and transfer functions[J]. Journal of Geophysical research ,1987, 92(B1): 633-648.

[4]Gamble T, Goubau W, Clarke J. Error analysis for remote reference magnetotellurics[J]. Geophysics, 1976,44(5):959-968.

[5]Chave A D, Thomson D J. Some comments on magnetotelluric response function estimation[J]. Journal of Geophysical Resarch,1989, 94(B10):14 215-14 226.

[6]Sims W, Bostick F, Smith H. The estimation of magntotelluric impedance tensor elements from measured data[J]. Geophysics,1971, 36(5): 938-942.

[7]Markus Eisel, Egbert G D. On the stability of magnetotelluric transfer function estimates and the realiability of their variances[J]. Geophysical Journal International, 2000, 144(1):65-82.

[8]程遠(yuǎn)志，胡祥云，王程.45°N86°E大地電磁參數(shù)標(biāo)準(zhǔn)點(diǎn)建設(shè)方法研究[J].工程地球物理學(xué)報(bào)，2011,8(5)：515-520.

[9]Egbert G D, Livelybrooks D W. Single Station magnetotellric impedance estimation: coherence weighting and the regression M-estimate[J]. Geophysics ,1996,61(4):964-970.

[10]Smirnow M Y. Magnetotellric data processing with a robust statistical procedure having a high breakdown point[J]. Geophysical Journal International,2003,152(1):1-7.

[11]Sutarno D. Development of robust magnetotelluric impedance estimation: a review[J]. Indonesian Journal of Physical. 2008,16(3):79-89.

[12]Jones A G, Chave A D, Egbert G D, et al. A comparison of techniques for magnetotelluric response function estimation[J]. Journal of Geophysical Research: Solid Earth (1978-2012), 1989, 94(B10):14 201-14 213.

[13]Egbert G D, Livelybrooks D W. Single Station magnetotellric impedance estimation: coherence weighting and the regression M-estimate[J]. Geophysics ,1996,61(4):964-970.

[14]湯井田，張弛，肖曉,等.大地電磁阻抗估計(jì)方法對(duì)比[J].中國(guó)有色金屬學(xué)報(bào)，2013,23(9)：2 351-2 358.

[15]陳高，金祖發(fā)，馬永生,等.大地電磁測(cè)深遠(yuǎn)參考技術(shù)及應(yīng)用效果[J].石油物探，2001，40(3)：112-117.

[16]仇根根，張小博，白大為,等.大地電磁法遠(yuǎn)參考處理技術(shù)壓制噪聲干擾的應(yīng)用效果分析[J].工程地球物理學(xué)報(bào)，2014,11(3)：305-310.

[17]Garcia X, Jones A G. A new methodology for the acquisition and processing of Audio-Magnetotelluric(AMT) data in the AMT dead band[J]. Geophysics, 2005, 70(5):G119-G126.

[18]Campanya J, Ledo J, Queralt P, et al. A new methodology to estimate Magnetotelluric(MT) Tensor relationships: Estimation of Local transfer-function by Combining Interstation Transfer-functions(ELICIT)[J]. Geophysical Journal International, 2014,198(1):484-494.

[19]Kappler K N, A data variance technique for automated despiking of magnetotelluric data with a remote reference[J].Geophysical prospecting,2012,60(1):179-191.

[20]王輝，魏文博，金勝,等.基于同步大地電磁時(shí)間序列依賴(lài)關(guān)系的噪聲處理[J].地球物理學(xué)報(bào)，2014,57(2)：531-545.

[21]Cui J L, Deng M, Jing J E, et al. Using independent component analysis to process magnetotelluric data[J]. Applied Mechanics and Materials,2013：295-298，2 795-2 798.

[22]Varentsov I M, Sokolova E Y, Martanus E R, et al. System of electromagnetic field transfer oprators for the BEAR array of simultaneous soundings; methods and results[J]. Izvestiya-Physics of the Solid Earth, 2003, 39(2):118-148.

[23]Egbert G D. Robust Multiple-station magnetotelluric data processing[J]. Geophysical Journal International.1997,130(2):475-496.

[24]Egbert G D. Processing and interpretation of electromagnetic induction array data[J]. Surveys in Geophysics,2002,23(2-3):207-249.

[25]Smirnov M Y, Egbert G D. Robust principal component analysis of electromagnetic arrays with missing data[J]. Geophysical Journal International,2012,190(3):1 423-1 438.

[26]Efron B, Gong G. A leisurely look at the bootstrap, the jackknife, and cross-validation[J]. The American Statistician, 1983,37(1): 36-48.

[27]Hinkley D V. Jackknifing in unbalanced situations[J]. Technometrics,1977,19(3):285-292.

[28]湯磊，湯洪志，韓雪平,等.數(shù)字形態(tài)濾波在音頻大地電磁法去噪中的應(yīng)用[J].工程地球物理學(xué)報(bào)，2013,10(4)：533-538.

The Research of Magnetotelluric Impedance Estimation

Ma Hairong1,2,Wang Shibiao3,Guo Rongwen1,2,Liu Jianxin1,2

(1.SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,ChangshaHunan410083,China； 2.HunanKeyLaboratoryofNon-ferrousResourcesandGeologicalHazardDetection,CentralSouthUniversity,ChangshaHunan410083,China； 3.GuangdongNonferrousGeologicalEnvironmentCenter,GuangzhouGuangdong510080,China)

Magnetotelluric method is mainly used to analyze the electrical distribution based on the morphological characteristics of apparent resistivity and phase curve, so the impedance tensor calculation is very important. In the actual detection, noise always disturbs the result of this transfer function in detection. In this paper we outline the recent progress of magnetotelluric impedance estimation. At first we studied several kinds of impedance estimation methods now commonly used, and then introduced two kinds of Jackknife methods which could calculate the stability of error bars of the transfer function.

magnetotelluric; impedance tensor; Jackknife

1672—7940(2016)04—0423—06

10.3969/j.issn.1672-7940.2016.04.004

國(guó)家自然科學(xué)基金(編號(hào)：41174103)

馬海榮(1990-)，女，碩士研究生，主要研究方向?yàn)榈刭|(zhì)工程。E-mail:1203966758@qq.com

柳建新(1962-)，男，博士，教授，博士生導(dǎo)師，主要從事礦產(chǎn)資源勘探、工程勘察領(lǐng)域的理論與應(yīng)用研究。

E-mail: ljx6666@126.com

P631.3

A

2016-03-04