基于同步大地電磁時間序列依賴關系的噪聲處理

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

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

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

3 地下信息探測技術與儀器教育部重點實驗室，北京 100083

1 引言

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

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

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

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

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

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

2 基本原理

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

圖1 同步大地電磁時間序列信號Fig.1 Time series of synchronous MT signals

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

2.1 依賴關系的求取

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

展開有：

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

圖2 MT測點點位分布圖Fig.2 Topography map showing MT stations layout

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

其中：mi，j表示第i個道的第j個數(shù)據(jù)，φi，j是mi，j對應的濾波系數(shù).

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

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

其中：

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

2.2 依賴關系穩(wěn)定性測試

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

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

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

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

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

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

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

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

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

3 仿真試算

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

（1）仿真噪聲

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

（2）噪聲識別

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

圖7 方差比對仿真噪聲的識別Fig.7 Identification of simulation noise by variance ratio method

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

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

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

（4）測深曲線對比

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

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

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

4 參考點噪聲測試

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

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

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

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

圖12 參考點信噪比不同時08232測點的測深曲線Fig.12 MT sounding curves of site 08232with reference data for different SNR

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

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

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

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

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

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

6 存在的問題及改進方法

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

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

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

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

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

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

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

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

（2）人為噪聲

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

7 結(jié)論

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

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

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

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

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

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

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

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

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.

附中文參考文獻

蔡劍華.2010.基于Hilbert－Huang變換的大地電磁信號處理方法與應用研究［博士論文］.長沙：中南大學.

董樹文，李廷棟，陳宣華等.2012.我國深部探測技術與實驗研究進展綜述.地球物理學報，55（12）：3884－3901.

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

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

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

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

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

魏文博.2002.我國大地電磁測深新進展及瞻望.地球物理學進展，17（2）：245－254.

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

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

葉高峰，魏文博，鄧明等.2010.青藏及華北陣列式區(qū)域大地電磁場標準觀測網(wǎng)建設方法與實驗.地質(zhì)學報，84（6）：801－807.

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