湯井田，李 晉＊，肖 曉，張林成，呂慶田

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083

2 中國(guó)地質(zhì)科學(xué)院礦產(chǎn)資源研究所，北京 100037

數(shù)學(xué)形態(tài)濾波與大地電磁噪聲壓制

湯井田1，李 晉1＊，肖 曉1，張林成1，呂慶田2

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083

2 中國(guó)地質(zhì)科學(xué)院礦產(chǎn)資源研究所，北京 100037

數(shù)學(xué)形態(tài)濾波是一種新型的非線性濾波方法.介紹了形態(tài)學(xué)濾波的基本原理，針對(duì)大地電磁信號(hào)表現(xiàn)出的非線性、非平穩(wěn)性和非最小相位特性，綜合結(jié)構(gòu)元素特征我們提出一種基于數(shù)學(xué)形態(tài)濾波的大地電磁噪聲壓制方法.為了有效抑制目標(biāo)信號(hào)中的噪聲干擾及修正統(tǒng)計(jì)偏倚現(xiàn)象，通過選用合理的結(jié)構(gòu)元素及形態(tài)開－閉、閉－開組合，將正、負(fù)結(jié)構(gòu)元素級(jí)聯(lián)構(gòu)造組合廣義形態(tài)濾波器對(duì)實(shí)測(cè)大地電磁信號(hào)進(jìn)行噪聲壓制.實(shí)驗(yàn)結(jié)果表明，該方法是切實(shí)可行的，有效地剔除了大尺度干擾及基線漂移，較好地還原了大地電磁信號(hào)的原始特征，修正了標(biāo)準(zhǔn)形態(tài)算子所產(chǎn)生的統(tǒng)計(jì)偏倚現(xiàn)象，去噪精度高.該方法計(jì)算速度快，具有潛在優(yōu)勢(shì)，為礦集區(qū)海量大地電磁信號(hào)與強(qiáng)干擾的分離提供了一條新的途徑，應(yīng)用前景廣闊.

數(shù)學(xué)形態(tài)濾波，結(jié)構(gòu)元素，非線性濾波，大地電磁，噪聲壓制

1 引 言

大地電磁測(cè)深法（Magnetotelluric，MT）自20世紀(jì)50年代誕生起，以野外施工簡(jiǎn)便、探測(cè)深度大、垂向分辨能力和水平分辨能力高等優(yōu)點(diǎn)，在探測(cè)地殼深部結(jié)構(gòu)方面已得到廣泛應(yīng)用，并已成為礦產(chǎn)資源勘查、地下水和地?zé)峥碧?、油氣普查、地震預(yù)報(bào)、巖石圈深部結(jié)構(gòu)探測(cè)等領(lǐng)域不可或缺的重要地球物理方法之一［1］.然而，由于天然電磁場(chǎng)信號(hào)弱、頻帶寬，礦區(qū)附近重工業(yè)密集等因素造成的環(huán)境噪聲以及人類活動(dòng)等人文電磁噪聲日益嚴(yán)重，導(dǎo)致大地電磁測(cè)量數(shù)據(jù)受到嚴(yán)重污染，這給地球物理勘探工作帶來巨大困難［2－3］.因此，對(duì)大地電磁信號(hào)中的強(qiáng)噪聲干擾進(jìn)行有效壓制，從而來改善大地電磁測(cè)量數(shù)據(jù)質(zhì)量，獲得無偏的阻抗估計(jì)是取得良好勘探結(jié)果的關(guān)鍵，也一直是大地電磁測(cè)深領(lǐng)域的研究前沿和熱點(diǎn)之一.通過研究證明，實(shí)際觀測(cè)的大地電磁信號(hào)是典型的非線性、非平穩(wěn)信號(hào)［4－5］.

分析國(guó)內(nèi)外相關(guān)文獻(xiàn)可知，目前的大地電磁去噪方法其能力都存在一定的局限性.互功率譜和加權(quán)平均功率譜法雖對(duì)不相關(guān)噪聲具有一定的抑制能力，但電磁噪聲往往是同時(shí)作用在各道電磁信號(hào)上，屬于相關(guān)噪聲，導(dǎo)致該方法失效［6－7］；遠(yuǎn)參考道法雖能消除同一測(cè)點(diǎn)各道之間不相關(guān)的電磁噪聲，并對(duì)各道之間同源電磁噪聲也有較好的效果，但經(jīng)遠(yuǎn)參考處理后，單點(diǎn)數(shù)據(jù)的誤差棒不同程度地變大了，特別是在受電磁干擾污染嚴(yán)重、校正量較大的數(shù)據(jù)段，該現(xiàn)象尤為明顯，且參考道距離的選擇是一個(gè)比較實(shí)際的問題［8］；Robust處理法對(duì)于MT資料中的非高斯正態(tài)分布噪聲的消除具有優(yōu)越性，能有效減少視電阻率和相位的分散度，但Robust法無法消除輸入端的噪聲，且無法剔除噪聲較多、能量較強(qiáng)時(shí)的近源電磁相關(guān)噪聲對(duì)數(shù)據(jù)的干擾［9－10］；小波變換雖能有效壓制大地電磁局部相關(guān)噪聲，但小波變換過分依賴于小波基函數(shù)的選取，有時(shí)隨著尺度增大，相應(yīng)正交基函數(shù)的頻譜局部性變差，使其對(duì)大地電磁信號(hào)更精細(xì)分解受到限制［11－12］；近年來，Hilbert－Huang變換應(yīng)用到電法勘探，并能有效抑制大地電磁信號(hào)中的工頻干擾，與小波變換相比不需要選擇基函數(shù)，具有更強(qiáng)的時(shí)頻刻畫能力.但因經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）是自適應(yīng)的，無法揭示每時(shí)段的頻率特性和能量差異所具有的細(xì)微性變化，分解得到的固有模態(tài)函數(shù)（IMF）具有多分辨性，對(duì)于每階IMF分量在大地電磁信號(hào)中的物理意義還有待進(jìn)一步探究，且該算法占用大量運(yùn)算時(shí)間，不適合實(shí)測(cè)大地電磁信號(hào)處理［13－14］；人機(jī)聯(lián)作去噪法是基于可視化技術(shù)的思想，將MT原始數(shù)據(jù)通過計(jì)算機(jī)圖形界面顯示處理，并使用人機(jī)聯(lián)作的方式去除噪聲.該方法雖能較好地降低噪聲能量，改善測(cè)量數(shù)據(jù)質(zhì)量，但操作時(shí)涉及了太多的人為因素，耗費(fèi)了很多時(shí)間和精力，不適合處理觀測(cè)周期長(zhǎng)和噪聲較多的數(shù)據(jù)，且操作者必須具備豐富的噪聲識(shí)別經(jīng)驗(yàn)，否則效果會(huì)適得其反［15］.

2008年以來，我國(guó)開展了深部探測(cè)技術(shù)與實(shí)驗(yàn)研究，需要在礦集區(qū)進(jìn)行大地電磁探測(cè)工作，尤其是我國(guó)東部和華南的礦集區(qū)，經(jīng)濟(jì)發(fā)達(dá)、礦山密布、人煙稠密、礦山開采的大功率直流電機(jī)車、高壓電網(wǎng)、電視塔、各種金屬管網(wǎng)、廣播電臺(tái)、雷達(dá)、通訊電纜及信號(hào)發(fā)射塔等造成的電磁干擾，嚴(yán)重污染了實(shí)際大地電磁信號(hào)，導(dǎo)致MT數(shù)據(jù)采集與處理相當(dāng)困難［16］.由于礦集區(qū)電磁噪聲類型復(fù)雜、能量大，相關(guān)性強(qiáng)，且頻率覆蓋范圍極寬，現(xiàn)有在頻率域的處理方法對(duì)該類強(qiáng)噪聲干擾無能為力.因此，如何從礦集區(qū)微弱大地電磁信號(hào)中壓制或分離強(qiáng)干擾是一項(xiàng)極具挑戰(zhàn)性的工作.針對(duì)這一系列不利因素，鑒于國(guó)內(nèi)對(duì)資源的強(qiáng)勁需求及快速處理礦集區(qū)海量MT數(shù)據(jù)，我們另辟蹊徑，考慮從時(shí)間域?qū)ふ矣嗅槍?duì)性的方法在礦集區(qū)獲取有效、優(yōu)質(zhì)的MT測(cè)量數(shù)據(jù).這項(xiàng)工作的開展對(duì)礦集區(qū)電磁法探測(cè)結(jié)果的處理和解釋具有重要意義，同時(shí)對(duì)探測(cè)地殼精細(xì)結(jié)構(gòu)，尋找深部控礦構(gòu)造具有非常重要的實(shí)際價(jià)值.

數(shù)學(xué)形態(tài)學(xué)（Mathematical Morphology，MM）是基于積分幾何、隨機(jī)集合論等數(shù)學(xué)理論建立起來的一種非線性信號(hào)處理方法［17］，形態(tài)濾波器是從數(shù)學(xué)形態(tài)學(xué)發(fā)展起來的一種新型的非線性濾波技術(shù).近年來，隨著形態(tài)學(xué)理論的飛速發(fā)展，形態(tài)學(xué)濾波被逐步推廣到一維信號(hào)處理領(lǐng)域［18－19］.本文針對(duì)大地電磁信號(hào)的特點(diǎn)，嘗試性地引入形態(tài)濾波技術(shù)對(duì)礦集區(qū)的大地電磁強(qiáng)干擾在時(shí)間域進(jìn)行信噪分離，采用形態(tài)學(xué)中的開、閉運(yùn)算，選擇合理的結(jié)構(gòu)元素，構(gòu)造適合大地電磁信號(hào)的組合廣義形態(tài)濾波器.實(shí)驗(yàn)結(jié)果表明，該方法能有效抑制大地電磁信號(hào)中的強(qiáng)噪聲干擾，較好地還原了大地電磁信號(hào)的原始特征，改善了MT測(cè)量數(shù)據(jù)質(zhì)量.

2 形態(tài)學(xué)濾波

數(shù)學(xué)形態(tài)學(xué)是1964年由法國(guó)數(shù)學(xué)家Matheron G和Serra J共同創(chuàng)立的一種信號(hào)分析方法.此后，兩人通過詳細(xì)闡述拓?fù)鋵W(xué)、隨機(jī)論、遞增映射和凸性分析等問題，逐步建立并完善了數(shù)學(xué)形態(tài)學(xué)的理論體系［20－21］.該方法最早是以圖像的形態(tài)特征作為研究對(duì)象，現(xiàn)已成功應(yīng)用于圖像處理、圖形分析、計(jì)算機(jī)視覺以及電能擾動(dòng)等工程實(shí)踐領(lǐng)域［22－25］.方法的基本思想是通過集合來描述目標(biāo)信號(hào)，集合各部分之間的關(guān)系則說明目標(biāo)信號(hào)的結(jié)構(gòu)特點(diǎn)，即設(shè)計(jì)一個(gè)稱為結(jié)構(gòu)元素的“探針”，通過探針在信號(hào)中不斷移動(dòng)來考察信號(hào)各部分之間的關(guān)聯(lián)［26－28］.該方法僅取決于待處理信號(hào)的局部特征，利用結(jié)構(gòu)元素對(duì)信號(hào)的幾何特征進(jìn)行局部匹配或修正，同時(shí)保留目標(biāo)信號(hào)主要的形狀特征，以達(dá)到提取有用信息、抑制噪聲的目的［29－30］.

2.1 基本原理

數(shù)學(xué)形態(tài)學(xué)中的形態(tài)變換包括腐蝕、膨脹、形態(tài)開、形態(tài)閉以及形態(tài)開、閉線性組合.以一維離散信號(hào)為例，數(shù)學(xué)描述如下［21］：

設(shè)輸入信號(hào)f（n）為定義在F＝｛0，1，…，N－1｝上的離散函數(shù)，結(jié)構(gòu)元素g（n）為定義在G＝｛0，1，…，M－1｝上的離散函數(shù)，且N?M，則f（n）關(guān)于g（n）的腐蝕和膨脹運(yùn)算分別定義為式中，符號(hào)Θ和⊕分別表示腐蝕和膨脹運(yùn)算.從以上定義可知，腐蝕和膨脹運(yùn)算等價(jià)于離散函數(shù)在滑動(dòng)濾波窗（相當(dāng)于結(jié)構(gòu)元素）內(nèi)的最小值和最大值濾波.其中，腐蝕運(yùn)算表示一個(gè)收縮過程，用來剔除邊界不平滑的凸起部分，算法減少了峰值，加寬了谷域，使目標(biāo)收縮，孔洞擴(kuò)張；膨脹運(yùn)算表示一個(gè)擴(kuò)張過程，用來填平邊界不平滑的凹陷部分，算法增大了谷值，擴(kuò)展了峰頂，使目標(biāo)擴(kuò)張，孔洞收縮.

形態(tài)開、閉運(yùn)算是在腐蝕和膨脹級(jí)聯(lián)的基礎(chǔ)上衍生而來的復(fù)合運(yùn)算，從而組成數(shù)學(xué)形態(tài)學(xué)中最基本的形態(tài)濾波器.f（n）關(guān)于g（n）的形態(tài)開、閉運(yùn)算分別定義為

式中，符號(hào)°和·分別表示形態(tài)開和形態(tài)閉運(yùn)算.其中，形態(tài)開運(yùn)算是對(duì)同一結(jié)構(gòu)元素進(jìn)行先腐蝕后膨脹，目的是消除目標(biāo)信號(hào)中的細(xì)節(jié)和毛刺，使目標(biāo)信號(hào)的輪廓光滑，從而剔除尖峰，抑制正脈沖噪聲；形態(tài)閉運(yùn)算是對(duì)同一結(jié)構(gòu)元素進(jìn)行先膨脹后腐蝕，目的是填平目標(biāo)信號(hào)中的小洞及裂縫，濾除低谷噪聲，從而補(bǔ)償谷底，抑制負(fù)脈沖噪聲.因此，實(shí)際應(yīng)用中常采用形態(tài)開、閉運(yùn)算相結(jié)合來抑制正、負(fù)脈沖噪聲.

Maragos［31－32］采用相同尺寸形狀的結(jié)構(gòu)元素，通過形態(tài)開、閉運(yùn)算級(jí)聯(lián)，定義了形態(tài)開－閉（OC）和閉－開（CO）濾波器：

由以上定義可知，形態(tài)開－閉（OC）和形態(tài)閉－開（CO）濾波器都能同時(shí)濾除目標(biāo)信號(hào)中的正、負(fù)脈沖噪聲，但由于形態(tài)開運(yùn)算的收縮性導(dǎo)致形態(tài)開－閉濾波器輸出偏小，而形態(tài)閉運(yùn)算的擴(kuò)張性導(dǎo)致形態(tài)閉－開濾波器輸出偏大，以至濾波結(jié)果存在統(tǒng)計(jì)偏倚現(xiàn)象［33］.因此，單獨(dú)使用它們并不能取得較好的濾波效果.為了有效抑制各種噪聲干擾和克服統(tǒng)計(jì)偏倚現(xiàn)象，常采用兩者的平均值來逼近原信號(hào)，構(gòu)造形態(tài)開－閉和形態(tài)閉－開的組合濾波器（OC－CO）如下：

式中，y（n）表示形態(tài)濾波器的輸出結(jié)果，ΨOC（CO）（g）表示形態(tài)濾波器的基本濾波運(yùn)算單元.

2.2 結(jié)構(gòu)元素的選取

形態(tài)濾波的質(zhì)量取決于所選擇的形態(tài)變換和結(jié)構(gòu)元素.其中，形態(tài)變換的選擇必須滿足一些基本的約束條件，而結(jié)構(gòu)元素的選取則需要根據(jù)具體情況來確定.

結(jié)構(gòu)元素在形態(tài)運(yùn)算中的作用類似于一般信號(hào)處理時(shí)的濾波窗口或參考模板，其尺寸和形狀都將對(duì)形態(tài)學(xué)運(yùn)算產(chǎn)生很大影響.采用不同的結(jié)構(gòu)元素可以提取出目標(biāo)信號(hào)不同的形狀特征，常見的結(jié)構(gòu)元素有直線型、三角型、圓盤型、正弦型、拋物線型以及其他多邊形組合.結(jié)構(gòu)元素的設(shè)計(jì)通常取決于待處理信號(hào)的形狀，一般一種結(jié)構(gòu)元素對(duì)一種噪聲有較好的濾除效果.待處理信號(hào)的形狀不同，則應(yīng)選用不同的結(jié)構(gòu)元素才能達(dá)到噪聲壓制的目的.相對(duì)而言，結(jié)構(gòu)元素越復(fù)雜，其濾除噪聲的能力就越強(qiáng)，但所耗費(fèi)的時(shí)間也就越長(zhǎng).

2.3 組合廣義形態(tài)濾波器的構(gòu)造

由于大地電磁噪聲類型復(fù)雜多樣，為了能靈活應(yīng)用不同形狀和尺寸的結(jié)構(gòu)元素以及修正標(biāo)準(zhǔn)形態(tài)算子所產(chǎn)生的統(tǒng)計(jì)偏倚現(xiàn)象，文中引入廣義形態(tài)濾波器對(duì)大地電磁信號(hào)進(jìn)行噪聲壓制［34］.廣義形態(tài)開－閉和形態(tài)閉－開濾波器的數(shù)學(xué)描述分別定義為

式中，g1、g2分別表示不同的結(jié)構(gòu)元素.

廣義形態(tài)濾波器的基本濾波運(yùn)算單元ΨGOC（GCO）（g1，g2）定義為

采用廣義形態(tài)開－閉和閉－開組合，能更好地保持目標(biāo)信號(hào)的幾何結(jié)構(gòu)特征，修正標(biāo)準(zhǔn)形態(tài)算子所產(chǎn)生的統(tǒng)計(jì)偏倚現(xiàn)象［35］.

由于圓盤型結(jié)構(gòu)元素具有旋轉(zhuǎn)不變性，避免了直線結(jié)構(gòu)平滑程度不夠的缺點(diǎn)，而拋物線型結(jié)構(gòu)元素能有效抑制脈沖噪聲干擾［36］.因此，文中選用圓盤型和拋物線型兩種結(jié)構(gòu)元素相結(jié)合來設(shè)計(jì)廣義形態(tài)濾波器.結(jié)構(gòu)元素形狀如圖1所示.

圖1 兩種結(jié)構(gòu)元素（a）圓盤型結(jié)構(gòu)元素；（b）拋物線型結(jié)構(gòu)元素.Fig.1 Two structuring elements（a）Disc type structuring element；（b）Parabolic structuring element.

圓盤型結(jié)構(gòu)元素定義為

拋物線型結(jié)構(gòu)元素定義為

以上兩種結(jié)構(gòu)元素只需改變參數(shù)k和L就能改變其幅度和寬度，從而控制結(jié)構(gòu)元素的大小.

考慮到大地電磁信號(hào)的準(zhǔn)對(duì)稱性和克服基線漂移，文中將正、負(fù)結(jié)構(gòu)元素級(jí)聯(lián)組成如圖2所示的組合廣義形態(tài)濾波器，其目的是進(jìn)一步抑制目標(biāo)信號(hào)中的各種噪聲干擾和消除統(tǒng)計(jì)偏倚現(xiàn)象［37－38］.

圖2 組合廣義形態(tài)濾波器Fig.2 The combination generalized morphological filter

圖2中，ΨGOC（GCO）（－g1，－g2）表示采用負(fù)的結(jié)構(gòu)元素組成的廣義形態(tài)基本濾波單元.

經(jīng)組合廣義形態(tài)濾波后，重構(gòu)的大地電磁有用信號(hào)定義為

3 實(shí)際資料應(yīng)用

圖3所示為安徽廬樅礦集區(qū)某測(cè)點(diǎn)的一段實(shí)測(cè)MT原始數(shù)據(jù).對(duì)該段數(shù)據(jù)的時(shí)間序列進(jìn)行分析可知，電道和磁道均不同程度受到了周期性的突跳、波動(dòng)等信號(hào)干擾，這些信號(hào)與穩(wěn)定的天然電磁場(chǎng)信號(hào)相比，具有振幅大、能量強(qiáng)、周期性明顯等特征.分析大量MT信號(hào)的時(shí)間序列可知，礦集區(qū)中的強(qiáng)噪聲干擾通常以方波噪聲、脈沖噪聲、三角波噪聲和似充放電三角波噪聲為主，這些噪聲都是由于礦集區(qū)中復(fù)雜的外界因素和人文因素所造成的.

為了驗(yàn)證本文方法的實(shí)用性，對(duì)該礦集區(qū)中受強(qiáng)干擾嚴(yán)重的MT數(shù)據(jù)進(jìn)行數(shù)學(xué)形態(tài)濾波研究.鑒于大地電磁信號(hào)的數(shù)據(jù)量龐大、噪聲類型極其復(fù)雜，文中僅選用具有以上典型干擾特征的電場(chǎng)分量EX和EY運(yùn)用形態(tài)濾波法進(jìn)行討論.

3.1 不同類型結(jié)構(gòu)元素濾波效果對(duì)比

圖4所示為采用不同類型的結(jié)構(gòu)元素對(duì)實(shí)測(cè)大地電磁信號(hào)EX進(jìn)行形態(tài)濾波的仿真效果圖.

從圖4可知，含方波噪聲的大地電磁信號(hào)分別經(jīng)三種不同類型的結(jié)構(gòu)元素的形態(tài)濾波器后，將小于或等于結(jié)構(gòu)元素的信號(hào)進(jìn)行了濾除，只保留了比結(jié)構(gòu)元素大的信號(hào)單元.結(jié)果表明：數(shù)學(xué)形態(tài)濾波對(duì)大地電磁強(qiáng)干擾具有較好的去噪能力.比較圖4可知，圓盤型和拋物線型結(jié)構(gòu)元素較直線型結(jié)構(gòu)元素濾波效果明顯，提取的形態(tài)輪廓清晰、平滑，重構(gòu)的大地電磁信號(hào)有效地剔除了大尺度干擾和基線漂移，突出了MT有用信號(hào)的相關(guān)局部特征.

3.2 同一類型不同尺寸結(jié)構(gòu)元素濾波效果對(duì)比

圖5所示為采用不同尺寸的圓盤型結(jié)構(gòu)元素對(duì)實(shí)測(cè)大地電磁信號(hào)EX的去噪效果圖.其中，結(jié)構(gòu)元素分別采用3點(diǎn)結(jié)構(gòu)元和5點(diǎn)結(jié)構(gòu)元.

分析圖5可知，5點(diǎn)圓盤型結(jié)構(gòu)元素的濾波效果比3點(diǎn)結(jié)構(gòu)元素的效果明顯.去噪前，含大尺度方波噪聲的能量達(dá)到正常信號(hào)的幾十倍，完全把MT正常信號(hào)湮沒.經(jīng)過5點(diǎn)形態(tài)濾波后，提取出大尺度方波噪聲的形態(tài)曲線，且曲線輪廓較3點(diǎn)結(jié)構(gòu)元自然、光滑，較好地保持了原始信號(hào)的特征.重構(gòu)后的MT信號(hào)基本濾除了由噪聲引起的突跳波形，強(qiáng)噪聲干擾與MT正常信號(hào)得到了有效分離.

圖6所示為采用不同尺寸的拋物線型結(jié)構(gòu)元素對(duì)實(shí)測(cè)大地電磁信號(hào)EY的去噪效果圖.

對(duì)比圖6可知，5點(diǎn)拋物線型結(jié)構(gòu)元素的濾波效果比3點(diǎn)結(jié)構(gòu)元素好，通過形態(tài)濾波提取出的含大尺度似充放電三角波噪聲的輪廓較3點(diǎn)結(jié)構(gòu)元素光滑，重構(gòu)的大地電磁信號(hào)波形較為平穩(wěn)，且保留了豐富的細(xì)節(jié)成分.

通過以上分析可知，選擇合適的結(jié)構(gòu)元素的尺寸能較好地獲取疊加在MT有用信號(hào)上的噪聲形態(tài)輪廓，重構(gòu)后的MT信號(hào)則基本還原了MT有用信號(hào)的原始特征.

3.3 傳統(tǒng)與組合廣義形態(tài)濾波效果對(duì)比

圖7所示為傳統(tǒng)形態(tài)濾波和根據(jù)圖2設(shè)計(jì)的組合廣義形態(tài)濾波的效果對(duì)比圖.

圖7 兩種算法的濾波效果對(duì)比（a）EX分量；（b）EY分量.Fig.7 Comparison of filtering effect chart for two algorithms（a）EXComponent；（b）EYComponent.

分析圖7可知，傳統(tǒng)形態(tài)濾波在獲取噪聲輪廓上出現(xiàn)很嚴(yán)重的毛刺現(xiàn)象，曲線不光滑、連續(xù)性差，且在部分曲率最大處造成了信號(hào)的失真，而組合廣義形態(tài)濾波幾乎能完整地勾勒出整段大尺度噪聲輪廓，曲線自然、光滑，重構(gòu)的MT信號(hào)則較好地保留了有用信號(hào)的細(xì)節(jié)信息，重現(xiàn)了原始MT信號(hào)的基本特征，從而保證了MT有用信號(hào)的可靠性與準(zhǔn)確性.圖8所示為該礦集區(qū)某測(cè)點(diǎn)的MT原始數(shù)據(jù)在TM模式下的視電阻率曲線圖.

從圖8可知，原始數(shù)據(jù)的視電阻率曲線整體連續(xù)性較差.在大于5.5Hz時(shí)曲線形態(tài)較為平穩(wěn)，5.5～0.055Hz處曲線呈45°左右漸近線快速上升，表現(xiàn)為典型的近源效應(yīng).在0.055Hz左右時(shí)，視電阻率值超過10000Ωm，在0.0055Hz時(shí)，視電阻率值快速下降至100Ωm，且低頻段誤差棒增大，并出現(xiàn)不同程度的突跳畸變，這些現(xiàn)象表明該測(cè)點(diǎn)數(shù)據(jù)受到了礦集區(qū)強(qiáng)噪聲干擾.

圖9所示為組合廣義形態(tài)濾波后重構(gòu)的MT數(shù)據(jù)的視電阻率曲線圖.

對(duì)比圖8可知，5.5～0.055Hz處的近源干擾基本消除，低頻段突跳頻點(diǎn)得到了有效恢復(fù)，0.055～0.0055Hz的低頻誤差棒明顯減小，曲線形態(tài)光滑、平穩(wěn)，整體連續(xù)性大為提高.實(shí)驗(yàn)結(jié)果表明，該測(cè)點(diǎn)經(jīng)組合廣義形態(tài)濾波后的數(shù)據(jù)質(zhì)量較原始數(shù)據(jù)有明顯改善，得到的視電阻率曲線為地下電性結(jié)構(gòu)提供了資料的可解釋性.

4 結(jié) 論

本文將數(shù)學(xué)形態(tài)濾波引入到大地電磁噪聲壓制領(lǐng)域，提出了一種基于數(shù)學(xué)形態(tài)濾波的大地電磁噪聲壓制方法，并應(yīng)用到礦集區(qū)實(shí)測(cè)大地電磁信號(hào)與強(qiáng)干擾分離中.介紹了數(shù)學(xué)形態(tài)學(xué)的基本原理，分析了結(jié)構(gòu)元素的選取方法，采用正、負(fù)結(jié)構(gòu)元素級(jí)聯(lián)構(gòu)造了組合廣義形態(tài)濾波器，比較了傳統(tǒng)形態(tài)濾波和組合廣義形態(tài)濾波的去噪效果.實(shí)驗(yàn)結(jié)果表明，形態(tài)學(xué)濾波能有效抑制大地電磁信號(hào)中的大尺度干擾和基線漂移，較好地保持了目標(biāo)信號(hào)的局部特征；設(shè)計(jì)的組合廣義形態(tài)濾波器在大地電磁噪聲壓制中的效果明顯，進(jìn)一步修正了由標(biāo)準(zhǔn)形態(tài)算子所產(chǎn)生的統(tǒng)計(jì)偏倚現(xiàn)象.方法的優(yōu)勢(shì)在于，不需考慮噪聲干擾是何種類型，只要選擇與目標(biāo)信號(hào)相匹配的結(jié)構(gòu)元素，設(shè)計(jì)合適的形態(tài)學(xué)濾波器，則能較好地還原大地電磁信號(hào)的原始特征.該方法較好地改善了MT測(cè)量數(shù)據(jù)品質(zhì)，且計(jì)算速度快，適合礦集區(qū)海量MT數(shù)據(jù)處理.可以預(yù)測(cè)：形態(tài)學(xué)濾波在該領(lǐng)域?qū)⒕哂袕V闊的應(yīng)用前景，為大地電磁信號(hào)與強(qiáng)干擾的有效分離以及開展大地電磁深部探測(cè)與深部找礦提供了新的解決途徑.

由于結(jié)構(gòu)元素的選取對(duì)形態(tài)濾波器的濾波效果至關(guān)重要，迄今為止，選取何種結(jié)構(gòu)元素及其尺寸只能通過反復(fù)實(shí)驗(yàn)獲得.因此，如何根據(jù)背景噪聲和待處理信號(hào)的形狀自適應(yīng)地選取結(jié)構(gòu)元素的類型和尺寸，以及合理運(yùn)用形態(tài)變換的線性組合將有待進(jìn)一步深入研究.

（References）

［1］ 湯井田，何繼善.可控源音頻大地電磁法及其應(yīng)用.長(zhǎng)沙：中南大學(xué)出版社，2005.Tang J T，He J S.Controllable Source Audio Magnetotelluric Method and Its Application（in Chinese）.Changsha：Central South University Press，2005.

［2］ Cai J H，Tang J T，Hua R H，et al.An analysis method for magnetotelluric data based on the Hilbert－Huang transform.Exploration Geophysics，2009，40（2）：197－205.

［3］ 肖曉，湯井田，周聰?shù)?廬樅礦集區(qū)大地電磁探測(cè)及電性結(jié)構(gòu)初探.地質(zhì)學(xué)報(bào)，2011，85（5）：873－886.Xiao X，Tang J T，Zhou C，et al.Magnetotelluric Sounding in the Lujiang－Zongyang ore－district and preliminary study of electrical structure.Acta Geologica Sinica（in Chinese），2011，85（5）：873－886.

［4］ 王書明，王家映.關(guān)于大地電磁信號(hào)非最小相位性的討論.地球物理學(xué)進(jìn)展，2004，19（2）：216－221.Wang S M，Wang J Y.Discussion on the non－minimum phase of magnetotelluric signals.Progress in Geophysics（in Chinese），2004，19（2）：216－221.

［5］ 王書明，王家映.大地電磁信號(hào)統(tǒng)計(jì)特征分析.地震學(xué)報(bào)，2004，26（6）：669－674.Wang S M，Wang J Y.Analysis on statistic characteristics of magnetotelluric signal.Acta Seismologica Sinica（in Chinese），2004，26（6）：669－674.

［6］ Kao D W，Rankin D.Enhancement of signal－to－noise ratio in magnetotelluric data.Geophysics，1977，42（1）：103－110.

［7］ Goubau W M，Gamble T D，Clarke J.Magnetotelluric data analysis：removal of Bias.Geophysics，1978，43（6）：1157－1169.

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

［9］ Egbert G D.Robust multiple－station magnetotelluric data processing.Geophysical Journal International，1997，130（2）：475－496.

［10］ 柳建新，嚴(yán)家斌，何繼善等.基于相關(guān)系數(shù)的海底大地電磁阻抗Robust估算方法.地球物理學(xué)報(bào)，2003，46（2）：241－245.Liu J X，Yan J B，He J S，et al.Robust estimation method of sea magnetotelluric impedance based on correlative coefficient.Chinese J.Geophys.（in Chinese），2003，46（2）：241－245.

［11］ Trad D O，Travassos J M.Wavelet filtering of magnetotelluric data.Geophysics，2000，65（2）：482－491.

［12］ 范翠松，李桐林，王大勇.小波變換對(duì)MT數(shù)據(jù)中方波噪聲的處理.吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），2008，38（增刊）：61－63.Fan C S，Li T L，Wang D Y.Treatment of wavelet transform for square wave noise in MT date.Journal of Jilin University（Earth Science Edition）（in Chinese），2008，38（Suppl.）：61－63.

［13］ 湯井田，化希瑞，曹哲民等.Hilbert－Huang變換與大地電磁噪聲壓制.地球物理學(xué)報(bào)，2008，51（2）：603－610.Tang J T，Hua X R，Cao Z M，et al.Hilbert－Huang transformation and noise suppression of magnetotelluric sounding data.Chinese J.Geophys.（in Chinese），2008，51（2）：603－610.

［14］ 湯井田，蔡劍華，任政勇等.Hilbert－Huang變換與大地電磁信號(hào)的時(shí)頻分析.中南大學(xué)學(xué)報(bào)（自然科學(xué)版），2009，40（5）：1399－1405.Tang J T，Cai J H，Ren Z Y，et al.Hilbert－Huang transform and time frequency analysis of magnetotelluric signal.Journal of Central South University（Science and Technology）（in Chinese），2009，40（5）：1399－1405.

［15］ 王大勇.長(zhǎng)江中下游礦集區(qū)綜合地質(zhì)地球物理研究——以九瑞、銅陵礦集區(qū)為例［博士論文］.吉林：吉林大學(xué)，2010.Wang D Y.The integrated geophysical and geological study in the ore belt of the middle and lower reach of the Yangtze river－The cases study of Tongling and Jiurui ore district［Ph.D.thesis］（in Chinese）.Jilin：Jilin University，2010.

［16］ 呂慶田，史大年，湯井田等.長(zhǎng)江中下游成礦帶及典型礦集區(qū)深部結(jié)構(gòu)探測(cè)——SinoProbe－03年度進(jìn)展綜述.地球?qū)W報(bào)，2011，32（3）：257－268.LüQ T，Shi D N，Tang J T，et al.Probing on deep structure of middle and lower reaches of the Yangtze Metallogenic Belt and typical ore concentration area：a review of annual progress of SinoProbe－03.Acta Geoscientica Sinica（in Chinese），2011，32（3）：257－268.

［17］ 龔煒，石青云.數(shù)字空間中的數(shù)學(xué)形態(tài)學(xué)——理論及應(yīng)用.北京：科學(xué)出版社，1997.Gong W，Shi Q Y.Mathematical Morphology in Digital Space（in Chinese）.Beijing：Science Press，1997.

［18］ 胡愛軍，唐貴基，安連鎖.基于數(shù)學(xué)形態(tài)學(xué)的旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)降噪方法.機(jī)械工程學(xué)報(bào)，2006，42（4）：127－130.Hu A J，Tang G J，An L S.De－noising technique for vibration signals of rotating machinery based on mathematical morphology filter.Chinese Journal of Mechanical Engineering（in Chinese），2006，42（4）：127－130.

［19］ 張文斌，楊辰龍，周曉軍.形態(tài)濾波方法在振動(dòng)信號(hào)降噪中的應(yīng)用.浙江大學(xué)學(xué)報(bào)（工學(xué)版），2009，43（11）：2096－2099.Zhang W B，Yang C L，Zhou X J.Application of morphology filtering method in vibration signal de－noising.Journal of Zhejiang University（Engineering Science）（in Chinese），2009，43（11）：2096－2099.

［20］ Matheron G.Random Sets and Integral Geometry.New York：Wiley Press，1975.

［21］ Serra J.Image Analysis and Mathematical Morphology.New York：Academic Press，1982.

［22］ 岳蔚，劉沛.基于數(shù)學(xué)形態(tài)學(xué)消噪的電能質(zhì)量擾動(dòng)檢測(cè)方法.電力系統(tǒng)自動(dòng)化，2002，26（7）：13－17.Yue W，Liu P.Detection of power quality disturbances based on mathematical morphology filter.Automation of Electric Power Systems（in Chinese），2002，26（7）：13－17.

［23］ 蔣立輝，耿蒙，趙春暉.基于廣義形態(tài)濾波和模糊邏輯的散斑噪聲抑制.紅外與激光工程，2005，34（1）：80－83.Jiang L H，Geng M，Zhao C H.Speckle suppressing based on generalized morphological filter and fuzzy logic.Infrared and Laser Engineering（in Chinese），2005，34（1）：80－83.

［24］ 張建成，吳新杰.形態(tài)濾波在實(shí)時(shí)信號(hào)處理中應(yīng)用的研究.傳感技術(shù)學(xué)報(bào)，2007，20（4）：828－831.Zhang J C，Wu X J.Research on applications of morphological filtering in real－time signal processing.Chinese Journal of Sensors and Actuators（in Chinese），2007，20（4）：828－831.

［25］ 李兵，張培林，任國(guó)全等.基于數(shù)學(xué)形態(tài)學(xué)的分形維數(shù)計(jì)算及在軸承故障診斷中的應(yīng)用.振動(dòng)與沖擊，2010，29（5）：191－194.Li B，Zhang P L，Ren G Q，et al.Mathematic morphologybased fractal dimension calculation and its application in fault diagnosis of roller bearings.Journal of Vibration and Shock（in Chinese），2010，29（5）：191－194.

［26］ 舒泓，王毅.基于數(shù)學(xué)形態(tài)濾波和Hilbert變換的電壓閃變測(cè)量.中國(guó)電機(jī)工程學(xué)報(bào)，2008，28（1）：111－114.Shu H，Wang Y.Voltage flicker detection method based on mathematical morphology filter and Hilbert transform.Proceedings of the Chinese Society for Electrical Engineering（in Chinese），2008，28（1）：111－114.

［27］ 李春枝，何榮建，田光明.數(shù)學(xué)形態(tài)濾波在振動(dòng)信號(hào)分析中的應(yīng)用研究.計(jì)算機(jī)工程與科學(xué)，2008，30（9）：126－128.Li C Z，He R J，Tian G M.Research on the application of the mathematical morphology filtering in vibration signal analysis.Computer Engineering ＆Science（in Chinese），2008，30（9）：126－128.

［28］ 程揚(yáng)軍，黃純，何朝暉等.基于自適應(yīng)順序形態(tài)濾波的電能質(zhì)量去噪算法.計(jì)算機(jī)仿真，2009，26（12）：218－220，314.Cheng Y J，Huang C，He Z H，et al.A de－nosing method of power quality signal based on adaptive order morphology.Computer Simulation（in Chinese），2009，26（12）：218－220，314.

［29］ 柏林，劉小峰，秦樹人.小波－形態(tài)－EMD綜合分析法及其應(yīng)用.振動(dòng)與沖擊，2008，27（5）：1－4，35.Bo L，Liu X F，Qin S R.Hybrid wavelet－morphology－EMD analysis and its application.Journal of Vibration and Shock（in Chinese），2008，27（5）：1－4，35.

［30］ 陳輝，郭科，胡英.數(shù)學(xué)形態(tài)學(xué)在地震信號(hào)處理中的應(yīng)用研究.地球物理學(xué)進(jìn)展，2009，24（6）：1995－2002.Chen H，Guo K，Hu Y.A study on application of mathematical morphology to seismic signal processing.Progress in Geophysics（in Chinese），2009，24（6）：1995－2002.

［31］ Maragos P，Schafer R W.Morphological filters－Part I：Their set－theoretic analysis and relations to linear shift－invariant filters.IEEE Trans.on ASSP，1987，35（8）：1153－1169.

［32］ Maragos P，Schafer R W.Morphological filters－Part II：Their relations to median，order－statistic，and stack filters.IEEE Trans.on ASSP，1987，35（8）：1170－1184.

［33］ 杜必強(qiáng)，唐貴基，石俊杰.旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)形態(tài)濾波器的設(shè)計(jì)與分析.振動(dòng)與沖擊，2009，28（9）：79－81.Du B Q，Tang G J，Shi J J.Design and analysis of morphological filter for vibration signals of a rotating machinery.Journal of Vibration and Shock（in Chinese），2009，28（9）：79－81.

［34］ 沈路，周曉軍，張文斌等.廣義數(shù)學(xué)形態(tài)濾波器的旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)降噪.振動(dòng)與沖擊，2009，28（9）：70－73.Shen L，Zhou X J，Zhang W B，et al.De－noising for vibration signals of a rotating machinery based on generalized mathematical morphological filter.Journal of Vibration and Shock（in Chinese），2009，28（9）：70－73.

［35］ 趙靜，何正友，錢清泉.利用廣義形態(tài)濾波與差分熵的電能質(zhì)量擾動(dòng)檢測(cè).中國(guó)電機(jī)工程學(xué)報(bào)，2009，29（7）：121－126.Zhao J，He Z Y，Qian Q Q.Detection of power quality disturbances utilizing generalized morphological filter and difference－entropy.Proceedings of the CSEE（in Chinese），2009，29（7）：121－126.

［36］ 趙曉群，王津.一種基于形態(tài)學(xué)的語(yǔ)音增強(qiáng)方法.同濟(jì)大學(xué)學(xué)報(bào)（自然科學(xué)版），2006，34（10）：1394－1397.Zhao X Q，Wang J.A speech enhancement method based on morphology.Journal of Tongji University（Natural Science）（in Chinese），2006，34（10）：1394－1397.

［37］ Trahanias P E.An approach to QRS complex detection using mathematical morphology.IEEE Trans.on Biomedical Engineering，1993，40（2）：201－205.

［38］ Wang J，Xu G H，Zhang Q，et al.Application of improved morphological filter to the extraction of impulsive attenuation signals.Mechanical Systems and Signal Processing，2009，23（1）：236－245.

Mathematical morphology filtering and noise suppression of magnetotelluric sounding data

TANG Jing－Tian1，LI Jin1＊，XIAO Xiao1，ZHANG Lin－Cheng1，LüQing－Tian2
1 School of Geosciences and Info－Physics，Central South University，Changsha410083，China
2 Institute of Mineral Resources，Chinese Academy of Geological Sciences，Beijing100037，China

Mathematical morphology filtering is a new non－linear filtering method.The basic principle of morphology filtering was introduced in this paper.Considering that magnetotelluric signals are non－linear，non－stationary，and non－minimum phase and taking into account the properties of structuring elements，we proposed a new method based on mathematical morphology filtering for noise suppression of magnetotelluric sounding data.In order to effectively suppress the noise interference of target signal and correct for statistical dominance phenomenon，we chose reasonable structuring elements and the combination of morphological open－closing and close－opening to construct the combination generalized morphology filtering.The filtering is composed of cascade positive and negative structuring elements，and its purpose is to suppress noise of actual magnetotelluric sounding data.Experimental results indicate that the proposed method is feasible and can effectively eliminate larger scale disturbance and baseline driftof magnetotelluric.In addition，the proposed method can better restore the original features of magnetotelluric，amend the statistical dominance phenomenon produced by standard morphological operators，and has better precision.The method is faster and has potential advantages，which provides a new way for massive magnetotelluric signal and strong interference separation in ore district，moreover，it has broad application prospects.

Mathematical morphology filtering，Structuring element，Non－linear filtering，Magnetotelluric sounding，Noise suppression

10.6038／j.issn.0001－5733.2012.05.036

P631

2011－12－01，2012－04－12收修定稿

國(guó)家科技專項(xiàng)“深部探測(cè)技術(shù)與實(shí)驗(yàn)研究專項(xiàng)（SinoProbe－03）”、國(guó)家自然科學(xué)基金資助項(xiàng)目（40930418、41104071）和湖南省教育廳資助科研項(xiàng)目（11B074）資助.

湯井田，男，1965年生，教授，博士生導(dǎo)師，主要從事信號(hào)處理及電磁法數(shù)值模擬研究.E－mail：jttang＠m(xù)ail.csu.edu.cn

＊通訊作者李晉，男，1981年生，博士研究生，主要從事大地電磁強(qiáng)噪聲壓制及信號(hào)處理研究.E－mail：geologylj＠163.com

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

10.6038／j.issn.0001－5733.2012.05.036.

Tang J T，Li J，Xiao X，et al.Mathematical morphology filtering and noise suppression of magnetotelluric sounding data.Chinese J.Geophys.（in Chinese），2012，55（5）：1784－1793，doi：10.6038／j.issn.0001－5733.2012.05.036.

（本文編輯 何 燕）