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

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

?

大地電磁阻抗計算方法進展

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

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

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

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

1　引　言

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

2　大地電磁阻抗計算方法

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

2.1最小二乘處理

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

(1)

(2)

應用最小二乘法，使其殘差和最小，即

(3)

則最小二乘解為

(4)

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

一般得到的阻抗形式為

Zxx=

(5)

Zxy=

(6)

Zyx=

(7)

Zyy=

(8)

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

(9)

(10)

2.2Robust 處理

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

2.2.1M回歸估計

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

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

(11)

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

(12)

式中：一般取k=1.5,k為調(diào)整量。當r小于1.5時，認為該點為正常值，損失函數(shù)與最小二乘一致；當r大于1.5時，用不同的損失函數(shù)來降低對應數(shù)據(jù)的權重。

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

(13)

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

Robust法協(xié)方差計算相對復雜，因為此時的傳遞函數(shù)估計值不能簡單的用線性公式來準確表示，協(xié)方差可由下式計算：

(14)

(15)

(16)

這里，wi表示權重。

2.2.2有界影響估計

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

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

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

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

1)HAT為對稱矩陣，冪等矩陣；

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

3)HAT對角線元素hii滿足0

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

5)—般而言，當HAT對角線元素hii滿足ii?2p/N時，認為對應的第i道輸入存在奇異值。類似于M回歸估計迭代格式，有界影響估計阻抗的迭代格式為

(17)

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

(18)

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

2.2.3重復中值估計

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

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

(19)

(20)

2.3遠參考處理

遠參考大地電磁測量，是將一遠參考點處的磁信號作為測點處的磁分量來進行阻抗張量的估算，求取阻抗張量Z的估計值Zr。和傳統(tǒng)的Z估計值不同，只要參考場與電、磁道噪音不相關的話，Zr不受任一場中噪音功率偏倚的影響[15,16]。

實際觀測值是實際信號和噪聲干擾之和，由以下式子表示：

Ex=Exs+Exn

(21)

式中，Ex表示電場實際觀測值；Exn表示噪聲干擾；Exs表示真實信號。而大地電磁場的線性關系只對信號項才成立，即

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

Zxx=

(22)

Zxy=

(23)

Zyx=

(24)

Zyy=

(25)

其中

(26)

式中，下標r表示參考道測得的磁場分量。

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

(27)

(28)

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

(29)

2.4其他

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

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

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

(30)

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

(31)

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

(32)

當測量數(shù)據(jù)為非線性時，使用刀切法來計算誤差棒更加準確。刀切法提供了一個相對簡單的方式來計算估計值的誤差棒。它也可用于張量分解的誤差分析。

3.1固定權重刀切法

(33)

權重hi由帽子矩陣定義。

以遠參考法為例，有

(34)

3.2刪除子集刀切法

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

3.3不足之處

對于Jackknife計算，有以下總結：

1)固定權重刀切法計算得出的誤差估計有類似的計算公式，且允許磁道的信號與噪聲有一定的相關性，如果相關性超出范圍，則會計算出錯誤的誤差估計值。

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

4　結　語

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

此外，由于天然電磁場強度低，大地電磁數(shù)據(jù)的信噪比往往不高，為了獲得有用信號，可采用一些新的信噪分離技術，如數(shù)學形態(tài)濾波技術[28]、Hilbert-Huang變換等，張量阻抗也從引入單參考道到多道數(shù)據(jù)綜合利用。同時，局部畸變?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]程遠志，胡祥云，王程.45°N86°E大地電磁參數(shù)標準點建設方法研究[J].工程地球物理學報，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].中國有色金屬學報，2013,23(9)：2 351-2 358.

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

[16]仇根根，張小博，白大為,等.大地電磁法遠參考處理技術壓制噪聲干擾的應用效果分析[J].工程地球物理學報，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]王輝，魏文博，金勝,等.基于同步大地電磁時間序列依賴關系的噪聲處理[J].地球物理學報，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)濾波在音頻大地電磁法去噪中的應用[J].工程地球物理學報，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

國家自然科學基金(編號：41174103)

馬海榮(1990-)，女，碩士研究生，主要研究方向為地質(zhì)工程。E-mail:1203966758@qq.com

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

E-mail: ljx6666@126.com

P631.3

A

2016-03-04