王輝，魏文博，3，金勝，3，葉高峰，3，景建恩，張樂天，董浩，李波，謝成良

1 中國(guó)地質(zhì)大學(xué)（北京）地球物理與信息技術(shù)學(xué)院，北京 100083

2 地質(zhì)過程與礦產(chǎn)資源國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100083

3 地下信息探測(cè)技術(shù)與儀器教育部重點(diǎn)實(shí)驗(yàn)室，北京 100083

1 引言

大地電磁測(cè)深（Magnetotellurics，簡(jiǎn)稱 MT）是一種通過研究天然電磁場(chǎng)的傳播規(guī)律來探測(cè)地下介質(zhì)電性結(jié)構(gòu)的勘探方法（Cagniard，1953）.方法提出之初，由于受微弱天然場(chǎng)信號(hào)高精度采集技術(shù)的限制，發(fā)展緩慢.20世紀(jì)90年代隨著科技的快速發(fā)展，MT儀器設(shè)備有了很大發(fā)展（魏文博，2002），如今，MT已經(jīng)成為研究區(qū)域斷裂構(gòu)造（謝成良等，2012）、流體運(yùn)動(dòng)（Wei et al.，2001；Bai et al.，2010），以及板塊構(gòu)造和大陸動(dòng)力學(xué)（金勝等，2012；張樂天等，2012）等問題的重要地球物理方法之一.

通常認(rèn)為大地電磁測(cè)深的高頻信號(hào)主要來自于雷電等中高空天氣活動(dòng)，低頻則來自太陽風(fēng)與地球磁層的相互作用，這些場(chǎng)源距地球很遠(yuǎn)；當(dāng)它傳播到達(dá)地球表面時(shí)，可以認(rèn)為基本滿足平面波場(chǎng)的理論條件，而其他不滿足平面波場(chǎng)條件的電磁場(chǎng)信號(hào)，則是MT的干擾噪聲（Junge，1996）.實(shí)際上，大地電磁測(cè)深（MT）的實(shí)測(cè)數(shù)據(jù)必然包含了大地電磁場(chǎng)信號(hào)和噪聲.因此，如何從受干擾的信號(hào)中求取精確的MT阻抗張量是通過大地電磁測(cè)深數(shù)據(jù)反演獲取較準(zhǔn)確的地下介質(zhì)電性結(jié)構(gòu)模型的關(guān)鍵.盡管最小二乘法可以有效去除高斯隨機(jī)噪聲；但對(duì)于實(shí)測(cè)數(shù)據(jù)而言，噪聲卻往往是非高斯分布的情況，因而利用常規(guī)的最小二乘算法通常很難得到理想的結(jié)果.為此，1986年Egbert和Booker在最小二乘法的基礎(chǔ)上提出了Robust估計(jì).Robust的實(shí)質(zhì)是離差加權(quán)估計(jì)，在統(tǒng)計(jì)過程中對(duì)離差大的數(shù)據(jù)給予較小的權(quán)，離差小的數(shù)據(jù)給予較大的權(quán)，從而使曲線連續(xù)圓滑，剔除“飛點(diǎn)”，但對(duì)輸入端的磁場(chǎng)噪聲無能為力.1978年Gamble等提出遠(yuǎn)參考技術(shù)，基于一定距離內(nèi)同步信號(hào)相關(guān)而噪聲不相關(guān)的原理，抑制不相關(guān)噪聲，得到阻抗的無偏估計(jì).隨著衛(wèi)星同步采集技術(shù)的成熟，遠(yuǎn)參考技術(shù)得到了廣泛應(yīng)用，該方法對(duì)不相關(guān)的磁場(chǎng)噪聲抑制效果明顯，但會(huì)增大誤差棒（Gamble et al.，1979b）；實(shí)際應(yīng)用中我們發(fā)現(xiàn)，對(duì)于電場(chǎng)噪聲和低頻噪聲的處理效果不佳，還存在參考距的影響等問題（Shalivahan and Bhattacharya，2002）.加拿大PHEONIX公司的處理系統(tǒng)，通過人工挑選，刪除受干擾時(shí)段的功率譜，提高M(jìn)T響應(yīng)數(shù)據(jù)的質(zhì)量，但處理過程非常耗時(shí)，處理人員須具有一定經(jīng)驗(yàn).隨著我國(guó)工業(yè)化進(jìn)程日益加快，人文電磁噪聲對(duì)大地電磁場(chǎng)的干擾日趨嚴(yán)重；在某些經(jīng)濟(jì)發(fā)達(dá)地區(qū)，微弱的MT信號(hào)幾乎被強(qiáng)干擾噪聲淹沒（湯井田等，2012b），尤其當(dāng)MT電場(chǎng)信號(hào)受到強(qiáng)相關(guān)噪聲的干擾時(shí)，傳統(tǒng)的MT數(shù)據(jù)處理方法得不到令人滿意的處理結(jié)果.

在傳統(tǒng)MT信號(hào)處理方法中，時(shí)頻轉(zhuǎn)換是以傅里葉變換（FFT）為基礎(chǔ)的技術(shù)，而FFT主要用于對(duì)“平穩(wěn)信號(hào)”的變換；但MT信號(hào)具有較強(qiáng)的“非平穩(wěn)性”，針對(duì)MT信號(hào)非平穩(wěn)性的特點(diǎn)（王書明和王家映，2004），近些年來，更多處理“非平穩(wěn)”信號(hào)的方法被引入到MT數(shù)據(jù)處理中，如小波變換（Wavelet）、希爾伯特－黃變換（HHT）和廣義S變換.小波變換的優(yōu)點(diǎn)在于對(duì)信號(hào)局部特征的分析，可以有效抑制隨機(jī)噪聲和局部相關(guān)噪聲，但小波變換的處理效果依賴于母波的選取（徐義賢和王家映，2000；Daniel and Travassos，2000）.HHT的優(yōu)點(diǎn)是自適應(yīng)性和時(shí)頻分析能力，其最小二乘法處理結(jié)果與FFT經(jīng)Robust估計(jì)的結(jié)果相當(dāng)，顯示出了優(yōu)越的噪聲處理潛力（Chen et al.，2012），缺點(diǎn)在于對(duì)信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解時(shí)，窗口的首尾端會(huì)引入誤差，而且分解速度慢，實(shí)用性不強(qiáng)（蔡劍華，2010）.廣義S變換是一種優(yōu)于小波變換的時(shí)頻分析方法，在時(shí)頻譜中分析噪聲，可以改善低頻段數(shù)據(jù)質(zhì)量，但窗函數(shù)易引起邊緣效應(yīng)（景建恩等，2012）.

在時(shí)域中對(duì)時(shí)間序列進(jìn)行預(yù)處理，直接濾除強(qiáng)干擾噪聲，可以防止短時(shí)間的強(qiáng)干擾噪聲影響到頻域中的多個(gè)頻點(diǎn)甚至整個(gè)頻帶的信號(hào)質(zhì)量.時(shí)域大地電磁場(chǎng)的噪聲處理方法主要有：人工神經(jīng)網(wǎng)絡(luò)濾波、Kalman濾波、經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）、形態(tài)濾波以及IARWR.人工神經(jīng)網(wǎng)絡(luò)濾波利用神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)和自適應(yīng)功能，對(duì)數(shù)據(jù)進(jìn)行分類，區(qū)分噪聲與信號(hào)，分類的正確性與訓(xùn)練網(wǎng)絡(luò)時(shí)的輸入信號(hào)和噪聲特征息息相關(guān)，該方法只能去除特征固定的噪聲（Manoj and Nagarajan，2003）.Kalman濾波的優(yōu)點(diǎn)是可以去除不相關(guān)的脈沖噪聲（Andrzej et al.，2009），實(shí)踐過程中我們發(fā)現(xiàn)，過程噪聲和測(cè)量噪聲的協(xié)方差難以確定.EMD可以按頻率高低自適應(yīng)地分解復(fù)雜的非平穩(wěn)信號(hào)，去除頻率固定的噪聲，如工頻干擾，但對(duì)于頻率成分復(fù)雜的噪聲，閥值的設(shè)定往往難以把握，分解速度慢，實(shí)用性不強(qiáng)（蔡劍華，2010）.形態(tài)濾波可以有效提取MT信號(hào)中的基線漂移和大尺度強(qiáng)干擾噪聲，計(jì)算速度快，但依賴于結(jié)構(gòu)元素形態(tài)和大小的選取（湯井田等，2012a）.IARWR方法通過方差比識(shí)別噪聲，利用維納濾波去除突變?cè)肼?，提高低頻數(shù)據(jù)質(zhì)量，但對(duì)濾波后數(shù)據(jù)精確度缺乏有效驗(yàn)證，還會(huì)引入階躍噪聲（Kappler，2012）.有學(xué)者把時(shí)域中的MT噪聲處理比喻為給原始觀測(cè)數(shù)據(jù)“動(dòng)手術(shù)”，用以說明時(shí)域去噪的必要性和需要注意的問題.一則，對(duì)于強(qiáng)干擾噪聲，如礦區(qū)的方波噪聲和脈沖噪聲，目前頻域的去噪方法很難消除，只能在時(shí)域中進(jìn)行預(yù)處理；二則，在時(shí)域去除噪聲必須十分謹(jǐn)慎，確保有效信號(hào)不被濾除，否則即使得到了光滑的MT測(cè)深曲線，其正確性也無法保證.總而言之，對(duì)強(qiáng)干擾噪聲，尤其是電場(chǎng)噪聲，目前還缺少有效的抑制方法，從時(shí)域或時(shí)頻域去噪將是MT噪聲處理的主要發(fā)展方向.

隨著我國(guó)“深部探測(cè)技術(shù)與實(shí)驗(yàn)研究（SinoProbe）”專項(xiàng)的啟動(dòng)，在全國(guó)范圍內(nèi)開展了陣列式大地電磁測(cè)深數(shù)據(jù)的觀測(cè)（董樹文等，2012；葉高峰等，2010），在對(duì)實(shí)測(cè)資料的處理中我們發(fā)現(xiàn)，山西和華南某些地區(qū)的數(shù)據(jù)受到嚴(yán)重的人文噪聲干擾，干擾源主要有：高壓線、鐵路、公路、發(fā)電站、采礦區(qū)等；這些噪聲最主要的特點(diǎn)是能級(jí)高、頻帶寬、噪聲波形復(fù)雜多變（湯井田，2012b），采用傳統(tǒng)處理方法得到的測(cè)深曲線具有明顯的近場(chǎng)源效應(yīng)，甚至全頻段的數(shù)據(jù)都無法使用；如何在強(qiáng)干擾地區(qū)提高M(jìn)T資料的信噪比是一項(xiàng)具有重要實(shí)際意義的工作.

本文從多道同步的MT時(shí)間序列出發(fā)，討論了多道信號(hào)之間的依賴關(guān)系，通過對(duì)多站同步實(shí)測(cè)數(shù)據(jù)的分析驗(yàn)證了依賴關(guān)系的穩(wěn)定性；同時(shí)，通過壓制仿真噪聲的試算，驗(yàn)證了基于多道時(shí)間序列依賴關(guān)系噪聲處理的可行性和正確性；最后，處理了實(shí)測(cè)含噪聲數(shù)據(jù).結(jié)果表明，該方法可以有效去除脈沖、方波和充放電等強(qiáng)干擾噪聲，抑制近場(chǎng)源噪聲對(duì)MT測(cè)深曲線的影響，提高數(shù)據(jù)質(zhì)量.

2 基本原理

可以把MT工作的基本原理看成一個(gè)線性系統(tǒng)（Jones et al.，1989），磁場(chǎng)為輸入信號(hào)，通過地球這個(gè)特殊的系統(tǒng)，輸出電場(chǎng)，通過測(cè)量電磁場(chǎng)就可以估算系統(tǒng)的傳遞函數(shù)，即阻抗，阻抗是穩(wěn)定的，通常情況下與時(shí)間無關(guān)，這是我們利用MT方法探測(cè)地下介質(zhì)電性結(jié)構(gòu)的物理基礎(chǔ).因此，雖然電磁場(chǎng)信號(hào)的隨機(jī)性很強(qiáng)（Luciano et al.，2011），但同一測(cè)點(diǎn)磁場(chǎng)與電場(chǎng)是互相依賴的，不同測(cè)點(diǎn)同步的磁場(chǎng)信號(hào)具有較強(qiáng)的相關(guān)性（圖1），電場(chǎng)信號(hào)之間的差異主要是由測(cè)點(diǎn)間地下介質(zhì)的電性結(jié)構(gòu)差異引起的，但這種差異是穩(wěn)定的，如同兩個(gè)具有強(qiáng)相關(guān)的輸入信號(hào)，分別通過穩(wěn)定的系統(tǒng)后，得到的輸出信號(hào)應(yīng)具有某種相關(guān)性，因此不同測(cè)點(diǎn)的電場(chǎng)信號(hào)之間也應(yīng)具有相對(duì)穩(wěn)定的關(guān)系，總而言之，同步多道電磁場(chǎng)信號(hào)之間應(yīng)具有相對(duì)穩(wěn)定的依賴關(guān)系.當(dāng)本地點(diǎn)電場(chǎng)或磁場(chǎng)信號(hào)在某些時(shí)刻受到噪聲干擾時(shí)，利用這種依賴關(guān)系和參考點(diǎn)的數(shù)據(jù)，合成本地道受干擾時(shí)段的數(shù)據(jù)，用合成數(shù)據(jù)代替噪聲段數(shù)據(jù)，組成新數(shù)據(jù)，就可以達(dá)到去除噪聲的目的.

圖1 同步大地電磁時(shí)間序列信號(hào)Fig.1 Time series of synchronous MT signals

實(shí)現(xiàn)這種方法需要解決兩個(gè)基本問題，一是依賴關(guān)系的求取方法，二是依賴關(guān)系的穩(wěn)定性.我們選取四個(gè)同步測(cè)點(diǎn)的高信噪比數(shù)據(jù)作為理想信號(hào)，用來討論這兩個(gè)基本問題，測(cè)點(diǎn)的點(diǎn)位分布如圖2所示.測(cè)點(diǎn)位于西藏?zé)o人區(qū)，遠(yuǎn)離人文噪聲干擾，其測(cè)深曲線如圖3所示，曲線連續(xù)光滑，誤差棒很小，忽略儀器噪聲，可以將此數(shù)據(jù)作為理想信號(hào).

2.1 依賴關(guān)系的求取

不失一般性，我們選取08232E10（參考點(diǎn)）的數(shù)據(jù)合成一段08232（本地點(diǎn)）的數(shù)據(jù)，具體計(jì)算方法如下，首先在本地道（如Ex）t0時(shí)刻，取出長(zhǎng)度為T的數(shù)據(jù)段，記為d，M為此時(shí)參考點(diǎn)K道（RREx，RREy，RRHx，RRHy，K＝4）數(shù)據(jù)組成的矩陣，假定d與M之間關(guān)系為φ，則有：

展開有：

將（2）式改寫成矩陣表示為

圖2 MT測(cè)點(diǎn)點(diǎn)位分布圖Fig.2 Topography map showing MT stations layout

圖3 高信噪比實(shí)測(cè)數(shù)據(jù)測(cè)深曲線Fig.3 MT sounding curves of practical data with high SNR

其中：mi，j表示第i個(gè)道的第j個(gè)數(shù)據(jù)，φi，j是mi，j對(duì)應(yīng)的濾波系數(shù).

由（3）式可以看出，d中的第n個(gè)值是由K個(gè)道的第（n－h(huán)Q）～（n＋hQ）數(shù)據(jù)加權(quán)求和得到的，φ的本質(zhì)是加權(quán)系數(shù)或?yàn)V波系數(shù).

為了求解濾波系數(shù)矩陣φ，假定φ不隨時(shí)間變化，因此可以選取另一時(shí)段t1，長(zhǎng)度為Tp的數(shù)據(jù)重新組成M和d，求解t1時(shí)段的φ，注意此時(shí)M和d數(shù)據(jù)必須是高信噪比的，代表的是天然大地電磁場(chǎng)信號(hào)的特征，稱為先驗(yàn)數(shù)據(jù)，這相當(dāng)于利用理想的輸入輸出信號(hào)求解系統(tǒng)響應(yīng)，基于最小均方誤差準(zhǔn)則，類似于最小二乘法估算阻抗，容易得到：

其中：

由于先驗(yàn)數(shù)據(jù)噪聲水平低，因此用（4）式求得的濾波系數(shù)φ表示的是信號(hào)之間的依賴關(guān)系，再由（1）式就可以求得t0時(shí)段本地道（Ex）的合成數(shù)據(jù)，注意（1）式中的M是t0時(shí)段的.

2.2 依賴關(guān)系穩(wěn)定性測(cè)試

從圖3中不難看出，08232S50和08232E10測(cè)點(diǎn)的兩個(gè)模式視電阻率曲線基本重合，說明兩者一維性較強(qiáng)，但兩者視電阻率值相差較大，08232和08533W30測(cè)點(diǎn)高頻部分XY模式和YX模式電阻率曲線重合，兩種模式電阻率曲線在中低頻分開，說明淺部介質(zhì)一維性較強(qiáng)，深部可能存在二維或三維的電性構(gòu)造，但08232視電阻率值較小，08533W30視電阻率值較大.表1是四個(gè)測(cè)點(diǎn)視電阻率值和相位值的統(tǒng)計(jì)結(jié)果，從表中也不難發(fā)現(xiàn)四個(gè)測(cè)點(diǎn)視電阻率值具有明顯差異，因此，四個(gè)測(cè)點(diǎn)地下介質(zhì)的電性構(gòu)造存在明顯差異，從圖1中也可以看出，四者的電場(chǎng)信號(hào)存在一定差異，但磁場(chǎng)信號(hào)仍具有很強(qiáng)的相關(guān)性.

測(cè)試方法如下：選擇與本地點(diǎn)具有不同參考距離，且存在明顯電性結(jié)構(gòu)差異的三個(gè)參考點(diǎn)，再從中選取不同時(shí)段、不同長(zhǎng)度的先驗(yàn)數(shù)據(jù)求取信號(hào)間的依賴關(guān)系，最后合成本地點(diǎn)（08232）的同一段數(shù)據(jù).用合成數(shù)據(jù)與原始理想數(shù)據(jù)的相干度，衡量合成數(shù)據(jù)的精確度，以此驗(yàn)證φ的穩(wěn)定性.如圖4所示，從08232選取一道長(zhǎng)度為30的原始數(shù)據(jù)（如Ey），記為dori，在dori的后方分別選取30段長(zhǎng)度不同的先驗(yàn)數(shù)據(jù)dpriori，先驗(yàn)數(shù)據(jù)的長(zhǎng)度按對(duì)數(shù)等分，原始數(shù)據(jù)和先驗(yàn)數(shù)據(jù)的時(shí)段分別為t0、tp，另外對(duì)不同參考點(diǎn)而言，選取的先驗(yàn)數(shù)據(jù)時(shí)段是不同的，如圖中tps，最后合成t0時(shí)段的Ey道數(shù)據(jù)記為dsyn，合成數(shù)據(jù)和原始數(shù)據(jù)相干度與先驗(yàn)數(shù)據(jù)長(zhǎng)度的關(guān)系如圖5所示.

表1 圖2中測(cè)點(diǎn)視電阻率值和相位值的統(tǒng)計(jì)Table 1 Statistics of apparent resistivity and phase values calculated from the sounding sites in Fig.2

圖4 先驗(yàn)數(shù)據(jù)的選取示意圖Fig.4 Selecting scheme of priori data

圖5 參考點(diǎn)和先驗(yàn)數(shù)據(jù)長(zhǎng)度對(duì)合成數(shù)據(jù)精度的影響Fig.5 Influence on accuracy of synthetic data by altering reference sites and length of priori data

圖6 不同參考點(diǎn)和先驗(yàn)數(shù)據(jù)長(zhǎng)度的合成數(shù)據(jù)Fig.6 Synthetic data of different reference sites and lengths of priori data

由圖5可知：①利用不同參考點(diǎn)求得的依賴關(guān)系即使是不同的，卻都可以精確地合成同一段本地點(diǎn)的數(shù)據(jù)；②無論磁場(chǎng)或電場(chǎng)信號(hào)，都可以通過參考點(diǎn)數(shù)據(jù)合成，且精度都較高；③在對(duì)數(shù)坐標(biāo)下，當(dāng)先驗(yàn)數(shù)據(jù)長(zhǎng)度與合成數(shù)據(jù)長(zhǎng)度的比值小于0.5時(shí)，合成數(shù)據(jù)與原始數(shù)據(jù)的相干度小于0.8，且不穩(wěn)定，當(dāng)兩者比值大于0.5（先驗(yàn)數(shù)據(jù)長(zhǎng)度約為合成數(shù)據(jù)長(zhǎng)度的3倍）時(shí)，合成數(shù)據(jù)與原始數(shù)據(jù)的相干度接近于1，兩者基本一致.圖6是不同先驗(yàn)數(shù)據(jù)長(zhǎng)度下，合成數(shù)據(jù)與原始數(shù)據(jù)的對(duì)比圖.

由圖6可知，當(dāng)先驗(yàn)數(shù)據(jù)長(zhǎng)度為91和153時(shí)，合成數(shù)據(jù)與原始數(shù)據(jù)基本一致，兩者的相干度在0.95以上，當(dāng)先驗(yàn)數(shù)據(jù)長(zhǎng)度為32時(shí)，兩者相干度小于0.8，說明選擇較長(zhǎng)時(shí)段信號(hào)求得的依賴關(guān)系包含了短時(shí)段信號(hào)之間的關(guān)系，從大地電磁信號(hào)特征的角度也不難理解，即短時(shí)段信號(hào)的頻譜信息包含在長(zhǎng)時(shí)段信號(hào)當(dāng)中.

上述試算結(jié)果表明，選擇不同參考距離、不同電性結(jié)構(gòu)參考點(diǎn)的不同時(shí)段的先驗(yàn)數(shù)據(jù)，都可以精確地合成同一段本地點(diǎn)的數(shù)據(jù)，說明多道同步大地電磁信號(hào)之間的依賴關(guān)系是穩(wěn)定的，其一，對(duì)同一參考點(diǎn)而言，依賴關(guān)系與先驗(yàn)數(shù)據(jù)的時(shí)段無關(guān)，只與先驗(yàn)數(shù)據(jù)的長(zhǎng)度有關(guān)；其二，即使不同參考點(diǎn)信號(hào)存在差異，依賴關(guān)系不同，但只要選擇足夠長(zhǎng)高信噪比的先驗(yàn)數(shù)據(jù)，就可以精確合成同一段本地點(diǎn)的數(shù)據(jù).因此，選擇較長(zhǎng)的、低噪聲的先驗(yàn)數(shù)據(jù)（大于3倍噪聲段數(shù)據(jù)長(zhǎng)度）尤為重要，在實(shí)際應(yīng)用過程中，需要人為挑選先驗(yàn)數(shù)據(jù).

3 仿真試算

天然大地電磁信號(hào)微弱，頻帶寬，與復(fù)雜的人文噪聲耦合在一起，難以區(qū)分.時(shí)域中處理噪聲的難點(diǎn)在于，去除強(qiáng)干擾噪聲的同時(shí)要保留微弱的有效信號(hào)成分，過分追求連續(xù)光滑的測(cè)深曲線，而忽略微弱的信號(hào)，可能會(huì)導(dǎo)致對(duì)地下介質(zhì)電性結(jié)構(gòu)的錯(cuò)誤認(rèn)識(shí).為了試驗(yàn)本文方法的有效性和正確性，進(jìn)行如下仿真試算.

（1）仿真噪聲

對(duì)08232（本地點(diǎn)）測(cè)點(diǎn)四個(gè)水平分量加入仿真噪聲，以08232E10測(cè)點(diǎn)作為參考點(diǎn)，兩點(diǎn)相距10.1km.仿真噪聲由方波、三角波和脈沖這三種常見的強(qiáng)干擾噪聲（湯井田，2012a）隨機(jī)組合而成，其形態(tài)和幅值也是隨機(jī)的，對(duì)每個(gè)信號(hào)道，隨機(jī)選擇20%的窗口加入噪聲；加入噪聲后的數(shù)據(jù)如圖8中dnoise所示.

（2）噪聲識(shí)別

利用本地道與參考道時(shí)間序列的方差比識(shí)別噪聲（Kappler，2012）.首先對(duì)同步后的數(shù)據(jù)加窗，用wk，n表示 第k道 第n個(gè) 窗 口 的 數(shù) 據(jù)，var（wk，n）和rrvar（wk，n）分別表示為本地道數(shù)據(jù)和參考道數(shù)據(jù)的方差，當(dāng)var（wk，n）／rrvar（wk，n）＞th時(shí)，th為閥值，則認(rèn)為本地點(diǎn)第k道第n個(gè)窗口的數(shù)據(jù)含有噪聲.這種方法基于強(qiáng)噪聲具有突變性從而引起方差變化的特點(diǎn)，從時(shí)域中定位本地點(diǎn)噪聲所在的通道和時(shí)段.圖7是對(duì)仿真噪聲的識(shí)別結(jié)果，識(shí)別率達(dá)到99.8%.

圖7 方差比對(duì)仿真噪聲的識(shí)別Fig.7 Identification of simulation noise by variance ratio method

（3）合成數(shù)據(jù)

選擇一段未加入噪聲的本地點(diǎn)數(shù)據(jù)和參考點(diǎn)數(shù)據(jù)作為先驗(yàn)數(shù)據(jù)，為了保證先驗(yàn)數(shù)據(jù)長(zhǎng)度大于噪聲段數(shù)據(jù)長(zhǎng)度的3倍，先驗(yàn)數(shù)據(jù)長(zhǎng)度為2min（time＝2012－06－1114∶06∶02—14∶08∶02），共1800個(gè)采樣點(diǎn)，窗口長(zhǎng)度為30個(gè)采樣點(diǎn)，按2.1節(jié)所述，用本文方法對(duì)識(shí)別后的噪聲進(jìn)行時(shí)域處理，用合成數(shù)據(jù)替換噪聲數(shù)據(jù)得到新數(shù)據(jù).圖8是原始數(shù)據(jù)、噪聲數(shù)據(jù)和合成數(shù)據(jù)的對(duì)比圖.

從圖8中可以看出合成數(shù)據(jù)成功去除了仿真噪聲，與原始信號(hào)數(shù)據(jù)一致，兩者的相干度高達(dá)0.99.

（4）測(cè)深曲線對(duì)比

圖8 原始數(shù)據(jù)、人工噪聲數(shù)據(jù)和合成數(shù)據(jù)的對(duì)比Fig.8 Comparison of original signals，artificial noisy data and synthetic data

圖9 人工噪聲數(shù)據(jù)和新數(shù)據(jù)測(cè)深曲線的對(duì)比Fig.9 Comparison of MT sounding curves between artificial noisy data and new data

仿真噪聲數(shù)據(jù)和去噪后數(shù)據(jù)的測(cè)深曲線如圖9所示，可以看出，加入噪聲后，視電阻率和相位曲線出現(xiàn)了很多“飛點(diǎn)”，難以分辨曲線形態(tài)，無法用于后續(xù)反演計(jì)算，新數(shù)據(jù)的測(cè)深曲線與原始曲線（圖3中08232）一致.仿真噪聲試驗(yàn)結(jié)果表明，用方差比方法可以精確識(shí)別方波、脈沖和三角波等典型的干擾噪聲，利用本文方法得到的合成數(shù)據(jù)去除了噪聲，保留了微弱的有效信號(hào)，提高了大地電磁測(cè)深資料的數(shù)據(jù)質(zhì)量和可信度.

4 參考點(diǎn)噪聲測(cè)試

為了研究參考點(diǎn)數(shù)據(jù)噪聲水平對(duì)合成數(shù)據(jù)精度的影響，進(jìn)行如下試算，對(duì)08232E10號(hào)測(cè)點(diǎn)（參考點(diǎn)）四道所有時(shí)段的數(shù)據(jù)，加入不同信噪比的高斯噪聲，再對(duì)上述仿真噪聲數(shù)據(jù)進(jìn)行處理，處理步驟與第3節(jié)相同.選擇圖8中的Ey數(shù)據(jù)，考察合成數(shù)據(jù)精度與參考點(diǎn)信噪比之間的關(guān)系，如圖10所示，隨著參考點(diǎn)數(shù)據(jù)信噪比的提高，合成數(shù)據(jù)與原始數(shù)據(jù)的相干度也隨之增大，當(dāng)信噪比大于等于40dB時(shí)，兩者的相干度大于0.94，當(dāng)信噪比小于30dB時(shí)，相干度急劇下降，說明當(dāng)參考點(diǎn)噪聲較強(qiáng)時(shí)，合成數(shù)據(jù)的誤差較大，這是因?yàn)閰⒖键c(diǎn)噪聲較強(qiáng)時(shí)，求得的本地點(diǎn)與參考點(diǎn)數(shù)據(jù)之間的關(guān)系并不完全表征信號(hào)間的依賴關(guān)系，而且用于計(jì)算合成數(shù)據(jù)的參考點(diǎn)數(shù)據(jù)含有噪聲，也會(huì)增大合成數(shù)據(jù)的誤差.

圖10 合成數(shù)據(jù)精度與參考點(diǎn)數(shù)據(jù)信噪比的關(guān)系Fig.10 Relationship between the accuracy of synthetic data and SNR of reference data

圖11 是參考點(diǎn)加入信噪比分別為10dB、20 dB、40dB和80dB的噪聲后，合成數(shù)據(jù)與原始數(shù)據(jù)的對(duì)比圖.從圖中容易看出，當(dāng)參考點(diǎn)信噪比增大時(shí)，合成數(shù)據(jù)越接近于原始數(shù)據(jù)，當(dāng)信噪比較低時(shí)，合成數(shù)據(jù)與原始數(shù)據(jù)的變化趨勢(shì)基本一致，幅值相差較大，合成數(shù)據(jù)的精度較低，但與仿真噪聲數(shù)據(jù)（圖8）相比，合成數(shù)據(jù)的誤差比噪聲要小得多.

圖11 參考點(diǎn)信噪比不同時(shí)的合成數(shù)據(jù)Fig.11 Variation of synthetic data with reference data for different SNR

圖12 參考點(diǎn)信噪比不同時(shí)08232測(cè)點(diǎn)的測(cè)深曲線Fig.12 MT sounding curves of site 08232with reference data for different SNR

圖13 3005（本地點(diǎn)）和1605（參考點(diǎn)）測(cè)點(diǎn)的先驗(yàn)數(shù)據(jù)Fig.13 Priori data of 3005（local）and 1605（reference）sites

圖12 是參考點(diǎn)信噪比分別為20dB和40dB時(shí)，用本文方法處理仿真噪聲數(shù)據(jù)后得到的測(cè)深曲線圖.與圖9對(duì)比易知，即使參考點(diǎn)數(shù)據(jù)含有噪聲，也可以得到正確且較光滑的測(cè)深曲線，當(dāng)參考點(diǎn)數(shù)據(jù)信噪比增大時(shí)，測(cè)深曲線越接近于原始測(cè)深曲線.在基本理論部分，已經(jīng)論證了合成數(shù)據(jù)的精度與參考點(diǎn)的距離及電性結(jié)構(gòu)無關(guān)，因此實(shí)際應(yīng)用中，首先要考慮將參考點(diǎn)布設(shè)在沒有人文噪聲干擾的“安靜”環(huán)境中，從而保證合成數(shù)據(jù)的精度.

5 實(shí)測(cè)數(shù)據(jù)試算

3005（本地點(diǎn)）和1605（參考點(diǎn)）大地電磁測(cè)點(diǎn)位于江西贛州，兩點(diǎn)相距112km，1605號(hào)點(diǎn)位于稻田中，周圍1km內(nèi)無明顯干擾源，數(shù)據(jù)質(zhì)量較高，3005號(hào)測(cè)點(diǎn)距一個(gè)大型火力發(fā)電站約2km，周圍有很多高壓線.從其測(cè)深曲線（圖15a）上可以發(fā)現(xiàn)，中高頻段數(shù)據(jù)受到明顯的近場(chǎng)源噪聲干擾，視電阻率曲線呈約45°直線上升，相位接近于0°或180°.從時(shí)間序列信號(hào)中，我們發(fā)現(xiàn)很多類似方波、脈沖、三角波等突變?cè)肼?，噪聲?qiáng)度大，分布廣，部分時(shí)間序列如圖14中dori所示.經(jīng)本文方法處理后得到的合成數(shù)據(jù)如圖14中dsyn，新數(shù)據(jù)測(cè)深曲線如圖15b所示.本例中選擇的先驗(yàn)數(shù)據(jù)如圖13所示，長(zhǎng)度為2min，共1800個(gè)采樣點(diǎn)，窗口長(zhǎng)度為30個(gè)采樣點(diǎn)，先驗(yàn)數(shù)據(jù)不含明顯噪聲，其中Hx和Hy相關(guān)性較好，Ex和Ey相差較大，但這并不影響去噪的結(jié)果.

圖14 3005測(cè)點(diǎn)原始數(shù)據(jù)與合成數(shù)據(jù)的對(duì)比Fig.14 Comparison of original data and synthetic signal of site 3005

從圖14中容易看出，合成數(shù)據(jù)去除了原始數(shù)據(jù)中的強(qiáng)干擾噪聲，保留了微弱變化的有效信號(hào).從圖15中不難看出，新數(shù)據(jù)有效抑制了中高頻段的近場(chǎng)源效應(yīng)，提高了視電阻率和相位曲線數(shù)據(jù)質(zhì)量，只是在1Hz附近幾個(gè)頻點(diǎn)上還存在“跳點(diǎn)”.

6 存在的問題及改進(jìn)方法

當(dāng)時(shí)間序列中含有階躍噪聲時(shí)，單用方差比方法是無法識(shí)別的，如圖16中橢圓標(biāo)出區(qū)域所示，方差代表的是數(shù)據(jù)的波動(dòng)情況，整體抬升或下降的噪聲數(shù)據(jù)，其方差是不變的，因此單用方差比無法識(shí)別階躍噪聲，所以新數(shù)據(jù)中仍保留了階躍噪聲，如圖16中dnew所示.

圖15 3005號(hào)點(diǎn)測(cè)深曲線（a）原始數(shù)據(jù)測(cè)深曲線；（b）新數(shù)據(jù)測(cè)深曲線.Fig.15 MT sounding curves of site 3005（a）Original data；（b）New data.

（1）平移法去除方波噪聲

方波噪聲可以看成一種特殊的階躍噪聲，識(shí)別方波噪聲最直接的方法是延長(zhǎng)窗口長(zhǎng)度，如圖16a所示，當(dāng)窗口長(zhǎng)度大于c2時(shí)，仍用方差比方法就能識(shí)別c2段的噪聲數(shù)據(jù).但是在試算的過程中我們發(fā)現(xiàn)，為了保證合成數(shù)據(jù)的精度，延長(zhǎng)窗口長(zhǎng)度的同時(shí)必須增大先驗(yàn)數(shù)據(jù)的長(zhǎng)度，而在強(qiáng)干擾地區(qū)，很難找到持續(xù)時(shí)間較長(zhǎng)的高信噪比數(shù)據(jù)作為先驗(yàn)數(shù)據(jù)，因?yàn)榇翱陂L(zhǎng)度太長(zhǎng)，先驗(yàn)數(shù)據(jù)含有明顯噪聲或長(zhǎng)度不夠，將導(dǎo)致合成數(shù)據(jù)的精度大大降低，甚至引入“假信號(hào)”，正如圖6中先驗(yàn)數(shù)據(jù)長(zhǎng)度為32時(shí)的情況，因此延長(zhǎng)窗口長(zhǎng)度的方法不適用于強(qiáng)干擾地區(qū)去除方波噪聲.

根據(jù)階躍噪聲具有突變性和突變前后信號(hào)變化不大的特點(diǎn)，如圖16a中c1至c2有突變，但兩者本身變化不大，因此可以用（5）式進(jìn)行識(shí)別.

std為求標(biāo)準(zhǔn)差，th為閥值，c2（1）為c2段第一個(gè)數(shù)據(jù)，c1（end）為c1段最后一個(gè)數(shù)據(jù)，為了避免將信號(hào)誤判為噪聲，th的值不易太小，試算發(fā)現(xiàn)電場(chǎng)閥值th＝5～10，磁場(chǎng)閥值th＝10～20時(shí)，可以識(shí)別幅值較大的階躍噪聲.定位噪聲后，對(duì)c2段數(shù)據(jù)按（6）式平移，就可以去除方波噪聲，記為ddesquare.

試算結(jié)果如圖16所示.

圖16表明，結(jié)合（5）式和（6）式可以有效識(shí)別和去除幅值較大的方波噪聲，保留了微弱的有效信號(hào).圖17是3005號(hào)測(cè)點(diǎn)去除方波噪聲后，單點(diǎn)處理和遠(yuǎn)參考處理得到的測(cè)深曲線圖，對(duì)比圖17a與圖15b可知，去除方波噪聲后的測(cè)深曲線基本消除了中高頻近場(chǎng)源效應(yīng)的影響和個(gè)別頻點(diǎn)“跳點(diǎn)”的情況，再經(jīng)遠(yuǎn)參考處理后，測(cè)深曲線更加光滑.因此，通過本文時(shí)域噪聲處理方法和頻域遠(yuǎn)參考技術(shù)的綜合處理，可以有效抑制近場(chǎng)源干擾，提高強(qiáng)干擾地區(qū)MT資料的信噪比.

（2）人為噪聲

平移法只能去除較理想的方波噪聲，真實(shí)噪聲是復(fù)雜的，平移后的數(shù)據(jù)看似消除了方波噪聲，實(shí)則引入了人為噪聲，但人為噪聲相對(duì)于原始方波噪聲是較小的，通過Robust估計(jì)和遠(yuǎn)參考處理可以抑制人為噪聲的影響，如圖17所示.為了確保去噪后數(shù)據(jù)的正確性，閥值的選擇不易太大，以免去除有效信號(hào)，另外，（5）式也可以識(shí)別階躍噪聲，但（6）式不能去除階躍噪聲，其難點(diǎn)在于信號(hào)“基準(zhǔn)線”難以確定.

7 結(jié)論

本文從理論上探討了多道同步大地電磁時(shí)間序列信號(hào)之間的依賴關(guān)系φ，選擇一段低噪聲的數(shù)據(jù)作為先驗(yàn)數(shù)據(jù)，用最小二乘法估算φ，再結(jié)合參考道的數(shù)據(jù)，合成本地道含噪聲時(shí)段的數(shù)據(jù)，用合成數(shù)據(jù)代替噪聲數(shù)據(jù)得到新數(shù)據(jù)，就可以在時(shí)域處理大地電磁噪聲.該方法避免了傳統(tǒng)方法中信噪分離的思路，直接利用參考道的信號(hào)合成本地道信號(hào)，無需考慮復(fù)雜的人文噪聲，且合成數(shù)據(jù)的精度較高.

對(duì)西藏?zé)o人區(qū)高信噪比實(shí)測(cè)數(shù)據(jù)的試算結(jié)果表明，φ是相對(duì)穩(wěn)定的，對(duì)同一參考點(diǎn)而言，與選擇的先驗(yàn)數(shù)據(jù)時(shí)段無關(guān)，只與先驗(yàn)數(shù)據(jù)的長(zhǎng)度有關(guān)；對(duì)不同參考點(diǎn)而言，即使φ不相同，但結(jié)合各參考點(diǎn)的數(shù)據(jù)，都可以精確合成本地點(diǎn)的同一段信號(hào)，與參考點(diǎn)地下電性結(jié)構(gòu)和參考距離無關(guān).

仿真噪聲的試算結(jié)果表明，方差比方法可以有效定位方波、脈沖和三角波等強(qiáng)干擾噪聲，合成數(shù)據(jù)與原始數(shù)據(jù)一致，對(duì)應(yīng)的測(cè)深曲線也一致；對(duì)參考點(diǎn)數(shù)據(jù)加入隨機(jī)噪聲的試算結(jié)果表明，合成數(shù)據(jù)的精度會(huì)隨著參考點(diǎn)數(shù)據(jù)信噪比的減小而降低，但與強(qiáng)噪聲相比，合成數(shù)據(jù)的影響要小得多.實(shí)測(cè)數(shù)據(jù)的處理結(jié)果表明，該方法可以有效抑制近場(chǎng)源強(qiáng)干擾噪聲，保留微弱的有效信號(hào)，提高測(cè)深曲線光滑性和連續(xù)性的同時(shí)保證了其正確性.針對(duì)方差比無法識(shí)別持續(xù)時(shí)間較長(zhǎng)的方波噪聲的問題，提出了一種簡(jiǎn)單的改進(jìn)方法，對(duì)實(shí)測(cè)數(shù)據(jù)的處理結(jié)果表明，該方法可以有效識(shí)別和去除方波噪聲.將本文方法與遠(yuǎn)參考技術(shù)相結(jié)合，可以有效抑制MT資料的近場(chǎng)源效應(yīng)，提高強(qiáng)干擾地區(qū)MT資料的利用率.

實(shí)際應(yīng)用本文方法時(shí)，要盡可能選擇沒有人文噪聲干擾的環(huán)境布設(shè)參考點(diǎn)，提高參考點(diǎn)的信噪比；先驗(yàn)數(shù)據(jù)的選擇非常關(guān)鍵，需要人為手動(dòng)挑選沒有明顯噪聲干擾時(shí)段的數(shù)據(jù)，保證先驗(yàn)數(shù)據(jù)求得的φ，表征的是天然場(chǎng)信號(hào)間的依賴關(guān)系，先驗(yàn)數(shù)據(jù)的長(zhǎng)度至少要大于噪聲段數(shù)據(jù)長(zhǎng)度的3倍.

圖16 方波噪聲的識(shí)別及去除Fig.16 Identification and removal of square wave noise

圖17 3005測(cè)點(diǎn)測(cè)深曲線（a）去除方波噪聲＋單點(diǎn)處理；（b）去除方波噪聲＋遠(yuǎn)參考處理.Fig.17 MT sounding curves of site 3005（a）Square wave removed with standard procedure；（b）Square wave removed with remote reference technique.

雖然強(qiáng)干擾地區(qū)噪聲強(qiáng)度大、成分復(fù)雜，但只要能夠找到一段沒有明顯噪聲的數(shù)據(jù)，就可以利用該方法去噪，而且處理過程中沒有改變時(shí)間序列數(shù)據(jù)格式，后續(xù)數(shù)據(jù)處理采用成熟的商業(yè)軟件即可完成，實(shí)用性較強(qiáng).

致謝 非常感謝魏文博教授在研究過程中細(xì)心的教導(dǎo)和幫助，以及對(duì)文章初稿的審閱和給出的寶貴意見，還要感謝張樂天師兄提供最新的文獻(xiàn)資料和董浩師兄在編寫代碼時(shí)給予的幫助.

Andrzej L，Danek T，Wojdyla M.2009.Application of kalman filter to noise reduction in multichannel data.SchedaeInformaticae，17（1）：63－73.

Bai D H，Unsworth M J，Meju M A，et al.2010.Crustal deformation of the eastern Tibetan plateau revealed by magnetotelluric imaging.NatureGeoscience，3（5）：358－362.

Cagniard L.1953.Basic theory of the magnetotelluric method of geophysical prospecting.Geophysics，18（3）：605－635.

Cai J H.2010.Study on processing of method of magnetotelluric signal and its application based on Hilbert－Huang transform［Ph.D.thesis］（in Chinese）.Changsha：Central South University.Chen J，Heincke B，Jegen M，et al.2012.Using empirical mode decomposition to process marine magnetotelluric data.Geophysical JournalInternational，190（1）：293－309.

Daniel O T，Travassos J M.2000.Wavelet filtering of magnetotelluric data.Geophysics，65（2）：482－491.

Dong S W，Li T D，Chen X H，et al.2012.Progress of deep exploration in mainland China：A review.ChineseJournalof Geophysics（in Chinese），55（12）：3884－3901.

Egbert G D，Booker J R.1986.Robust estimation of geomagnetic transfer functions.GeophysicalJournalInternational，87（1）：173－194.

Egbert G D.1997.Robust multiple－station magnetotelluric data processing.GeophysicalJournalInternational，130（2）：475－496.

Gamble T D，Goubau W M，Clarke J.1979a.Magnetotellurics with a remote magnetic reference.Geophysics，44（1）：53－68.

Gamble T D，Goubau W M，Clarke J.1979b.Error analysis for remote reference magnetotellurics.Geophysics，44（5）：959－968.

Jin S，Zhang L T，Jin Y J，et al.2012.Crustal electrical structure along the Hezuo－Dajing profile across the northeastern margin of the Tibetan Plateau.ChineseJournalofGeophysics（in Chinese），55（12）：3979－3990.

Jing J E，Wei W B，Chen H Y，et al.2012.Magnetotelluric sounding data processing based on generalized S transformation.ChineseJournalofGeophysics（in Chinese），55（12）：4015－4022.

Jones A G，Chave A D，Egbert G，et al.1989.A comparison of techniques for magnetotelluric response function estimation.JournalofGeophysicalResearch，94（10）：201－213.

Junge A.1996.Characterization of and correction for cultural noise.SurveysinGeophysics，17（4）：361－391.

Kappler K N.2012.A data variance technique for automated despiking of magnetotelluric data with a remote reference.GeophysicalProspecting，60（1）：179－191.

Luciano T，Lovallo M，Han－Lun H，et al.2011.Analysis of dynamics in magnetotelluric data by using the Fisher－Shannon method.PhysicaA：StatisticalMechanicsandItsApplications，390（7）：1350－1355.

Manoj C，Nagarajan N.2003.The application of artificial neural networks to magnetotelluric time－series analysis.Geophysical JournalInternational，153（2）：409－423.

Shalivahan，Bhattacharya B B.2002.How remote can the far remote reference site for magnetotelluric measurements be？JournalofGeophysicalResearch，107（6）：119－126.

Tang J T，Li J，Xiao X，et al.2012a.Mathematical morphology filtering and noise suppression of magnetotelluric sounding data.ChineseJournalofGeophysics（in Chinese），55（5）：1784－1793.

Tang J T，Xu Z M，Xiao X，et al.2012b.Effect rules of strong noise on magnetotelluric sounding in the Luzong ore cluster area.ChineseJournalofGeophysics（in Chinese），55（12）：4147－4159.

Wang S M，Wang G Y.2004.Statistical analysis of magnetotelluric signal.ActaSeismologicaSinica（in Chinese），26（6）：669－674.

Wei W B.2002.Progress and prospects of magnetotelluric sounding in China.ProgressinGeophysics（in Chinese），17（2）：245－254.

Wei W B，Unsworth M，Jones A G，et al.2001.Detection of widespread fluids in the Tibetan crust by magnetotelluric studies.Science，292（5517）：716－719.

Xie C L，Ye G F，Wei W B，et al.2012.Electrical features of the main faults beneath northern Tibetan Plateau.ChineseJournal ofGeophysics（in Chinese），55（12）：3991－4002.

Xu Y X，Wang G Y.2000.Power spectrum estimation for magnetotelluric signal based on continuous wavelet transform.ChineseJournalofGeophysics（in Chinese），43（5）：677－683.Ye G F，Wei W B，Deng M，et al.2010.Construction methods and experiments for magnetotelluric standard network at the Tibetan Plateau and North China.ActaGeologicaSinica（in Chinese），84（6）：801－807.

Zhang L T，Jin S，Wei W B，et al.2012.Electrical structure of crust and upper mantle beneath the eastern margin of Tibetan plateau and the Sichuan basin.ChineseJournalofGeophysics（in Chinese），55（12）：4126－4137.

附中文參考文獻(xiàn)

蔡劍華.2010.基于Hilbert－Huang變換的大地電磁信號(hào)處理方法與應(yīng)用研究［博士論文］.長(zhǎng)沙：中南大學(xué).

董樹文，李廷棟，陳宣華等.2012.我國(guó)深部探測(cè)技術(shù)與實(shí)驗(yàn)研究進(jìn)展綜述.地球物理學(xué)報(bào)，55（12）：3884－3901.

金勝，張樂天，金永吉等.2012.青藏高原東北緣合作—大井剖面地殼電性結(jié)構(gòu)研究.地球物理學(xué)報(bào)，55（12）：3979－3990.

景建恩，魏文博，陳海燕等.2012.基于廣義S變換的大地電磁測(cè)深數(shù)據(jù)處理.地球物理學(xué)報(bào)，55（12）：4015－4022.

湯井田，李晉，肖曉等.2012a.數(shù)學(xué)形態(tài)濾波與大地電磁噪聲壓制.地球物理學(xué)報(bào)，55（5）：1784－1793.

湯井田，徐志敏，肖曉等.2012b.廬樅礦集區(qū)大地電磁測(cè)深強(qiáng)噪聲的影響規(guī)律.地球物理學(xué)報(bào)，55（12）：4147－4159.

王書明，王家映.2004.大地電磁信號(hào)統(tǒng)計(jì)特征分析.地震學(xué)報(bào)，26（6）：669－674.

魏文博.2002.我國(guó)大地電磁測(cè)深新進(jìn)展及瞻望.地球物理學(xué)進(jìn)展，17（2）：245－254.

謝成良，葉高峰，魏文博等.2012.藏北高原主要斷裂帶電性結(jié)構(gòu)特征.地球物理學(xué)報(bào)，55（12）：3991－4002.

徐義賢，王家映.2000.基于連續(xù)小波變換的大地電磁信號(hào)譜估計(jì)方法.地球物理學(xué)報(bào)，43（5）：677－683.

葉高峰，魏文博，鄧明等.2010.青藏及華北陣列式區(qū)域大地電磁場(chǎng)標(biāo)準(zhǔn)觀測(cè)網(wǎng)建設(shè)方法與實(shí)驗(yàn).地質(zhì)學(xué)報(bào)，84（6）：801－807.

張樂天，金勝，魏文博等.2012.青藏高原東緣及四川盆地的殼幔導(dǎo)電性結(jié)構(gòu)研究.地球物理學(xué)報(bào)，55（12）：4126－4137.