王正海，耿 欣，姚卓森，胡光道，方 臣

1.中山大學(xué)地球科學(xué)系，廣州 510275

2.中國(guó)科學(xué)院地質(zhì)與地球物理研究所，北京 100029

3.中國(guó)地質(zhì)大學(xué)資源學(xué)院，武漢 430074

0 引言

近年來，以天然電磁場(chǎng)為場(chǎng)源的大地電磁（MT）觀測(cè)在深部礦產(chǎn)勘查和工程物探中應(yīng)用廣泛。但在應(yīng)用中，由于噪聲干擾，有用信息往往被噪聲所掩蓋，嚴(yán)重影響后續(xù)視電阻率和阻抗的計(jì)算及目標(biāo)信息的提取[1-2]。如何利用采集得到的時(shí)間序列數(shù)據(jù)進(jìn)行測(cè)區(qū)的噪聲分析，從而進(jìn)一步實(shí)現(xiàn)相應(yīng)噪聲的壓制，是當(dāng)前MT數(shù)據(jù)預(yù)處理的主要內(nèi)容之一。脈沖類噪聲是一種能引起電磁波的電脈沖，其強(qiáng)度很大，但持續(xù)時(shí)間較短，頻帶很寬，存在于MT信號(hào)低頻段、中頻段和高頻段。脈沖干擾在頻率域中沒有任何明顯特征，但是在時(shí)間域上，由于其能量很大，振幅是正常信號(hào)的幾倍，非常容易識(shí)別[3-4]。因此，需要在時(shí)間域上尋找合適的方法進(jìn)行脈沖類噪聲抑制處理。

傳統(tǒng)大地電磁法數(shù)據(jù)處理以常規(guī)傅里葉變換和功率譜分析為基礎(chǔ) ，通常假定信號(hào)為平穩(wěn)或分段平穩(wěn)，采用適當(dāng)?shù)姆治龇椒?，如傅里葉變換、小波變換等對(duì) MT信號(hào)進(jìn)行線性化處理[5-6]。近年來，很多學(xué)者已經(jīng)證明大地電磁信號(hào)具有非平穩(wěn)、非線性、非高斯的特征。上述基于線性平穩(wěn)假設(shè)的處理方法不能很好地進(jìn)行MT信號(hào)噪聲處理。1998年，華裔科學(xué)家黃鍔[7]提出了經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）算法，將信號(hào)分解成多個(gè)固有模態(tài)函數(shù)（IMF），且IMF是單分量函數(shù)；然后對(duì)每一個(gè)IMF進(jìn)行Hibert變換，并定義瞬時(shí)頻率，分析其時(shí)-頻譜。EMD分解基于數(shù)據(jù)自身的特點(diǎn)，具有自適應(yīng)的特征，而IMF又具有多尺度濾波的特點(diǎn)[7-8]。因此，EMD是處理非平穩(wěn)、非線性信號(hào)的有效方法，被成功運(yùn)用到流體力學(xué)、地震信號(hào)分析、基礎(chǔ)結(jié)構(gòu)監(jiān)測(cè)等多個(gè)領(lǐng)域[9-10]。筆者主要討論EMD在大地電磁測(cè)深系統(tǒng)信號(hào)分析和脈沖類噪聲處理中的應(yīng)用效果。

1 經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）

經(jīng)驗(yàn)?zāi)B(tài)分解又叫Huang變換，為滿足IMF為固有模態(tài)函數(shù)，它必須滿足以下2個(gè)假設(shè)[9-10]：1）在整個(gè)數(shù)據(jù)序列中，極值點(diǎn)的數(shù)量與過零點(diǎn)的數(shù)量必須相等，或者最多相差不能多于一個(gè)。2）在任一時(shí)間點(diǎn)上，由局部最大值和局部最小值分別確定的2條包絡(luò)線的平均值為0。

EMD分解通過三次樣條方法多次迭代擬合實(shí)現(xiàn)：1）找出原始數(shù)據(jù)序列x（t）的局部最大值和局部最小值，采用三次樣條方法擬合，得到信號(hào)的emax（t）和emin（t）。2）求每個(gè)時(shí)刻下emax（t）和emin（t）的平均值，即瞬時(shí)平均值m（t），代表著每一時(shí)刻信號(hào)的“瞬時(shí)平衡位置”：

3）用原始數(shù)據(jù)x（t）減去瞬時(shí)平均值 m（t），得到c1（t）：

式中，c1（t）被稱為用原始數(shù)據(jù)的第一個(gè)本征模態(tài)函數(shù)分量（IMF）。然后，用原始數(shù)據(jù)減去c1（t）即可得到剩余數(shù)據(jù)r1（t）：

重復(fù)上述相同的迭代過程，通過EMD的一次次篩選，就可以獲得信號(hào)的多個(gè)IMF分量和一個(gè)逼近分量rn（t），從而信號(hào)可表示為

經(jīng)驗(yàn)?zāi)B(tài)分解其實(shí)是對(duì)非線性非平穩(wěn)信號(hào)的一種平穩(wěn)化處理，將信號(hào)中不同尺度的信息分解開來，使其產(chǎn)生的一系列不同尺度的信號(hào)均能滿足后續(xù)Hilbert變換要求的窄帶特性[11-12]。經(jīng)驗(yàn)?zāi)B(tài)分解是完全可逆的，將分解所得的每一個(gè)序列信息逐一相加即可得到最初的原始信息。因此，經(jīng)驗(yàn)?zāi)B(tài)分解可以在不同尺度下對(duì)數(shù)據(jù)進(jìn)行處理后，再重新構(gòu)建原數(shù)據(jù)，達(dá)到對(duì)原數(shù)據(jù)處理的目的。利用EMD變換進(jìn)行信號(hào)處理時(shí)，針對(duì)不同目的，對(duì)每一階IMF分量頻譜進(jìn)行相應(yīng)處理，再分別利用式（4）進(jìn)行信號(hào)重構(gòu)，獲得相應(yīng)的信號(hào)處理效果。

2 脈沖干擾信號(hào)的EMD去噪

2.1 脈沖噪聲分析

脈沖干擾多出現(xiàn)在電場(chǎng)中，電場(chǎng)中的脈沖干擾主要來源于測(cè)點(diǎn)附近介質(zhì)層中的脈沖型游散電流，多為測(cè)點(diǎn)周圍接地用電動(dòng)力設(shè)備的開關(guān)或其負(fù)荷的突然改變?cè)斐傻摹ｋ妶?chǎng)和磁場(chǎng)中同時(shí)出現(xiàn)的脈沖干擾則多來自雷暴，或者是人類的電磁感應(yīng)信號(hào)。

圖1 EH-4實(shí)測(cè)數(shù)據(jù)中各頻段的電磁脈沖干擾Fig.1 The electromagnetic pulse interference of the band of EH-4data

圖1為西藏白容礦區(qū)利用美國(guó)EMI公司和Geo-metrics公司聯(lián)合研制的EH-4大地電磁測(cè)深系統(tǒng)實(shí)測(cè)剖面28號(hào)測(cè)點(diǎn)的含脈沖干擾的MT信號(hào)，單純的電脈沖較少出現(xiàn)，但是在幾乎所有的測(cè)點(diǎn)內(nèi)，都不同程度地出現(xiàn)了或多或少的電磁脈沖干擾。圖1中低頻段和高頻段的電磁脈沖干擾在時(shí)間域中很明顯地區(qū)分出來。由圖1可知，該地區(qū)電磁脈沖能量極強(qiáng)，基本都達(dá)到了EH-4當(dāng)前增益系數(shù)下的飽和測(cè)量值。在這種情況下，很多有用的信號(hào)都被強(qiáng)的電磁脈沖干擾所覆蓋，不利于目標(biāo)信息的提取及后續(xù)視電阻率和阻抗的計(jì)算。這類干擾若出現(xiàn)在電場(chǎng)中，將使阻抗的幅值嚴(yán)重向上偏移；而若出現(xiàn)在磁場(chǎng)中，則幅值向下偏移。

2.2 EMD去噪算法與步驟

EH-4大地電磁測(cè)深系統(tǒng)信號(hào)中脈沖干擾處理的關(guān)鍵是將脈沖噪聲從MT信號(hào)中分離出來，進(jìn)而對(duì)其進(jìn)行相應(yīng)處理，抑制脈沖干擾對(duì)MT信號(hào)的影響。

利用EMD變換進(jìn)行脈沖干擾處理分三步進(jìn)行：1）將受脈沖干擾影響的 MT信號(hào)x（t）經(jīng)EMD分解后得到N個(gè)固有模態(tài)函數(shù)IMF；2）對(duì)每一階的IMF選擇一個(gè)合適的瞬時(shí)頻率閥值，截?cái)喑鲈撻y值的部分，進(jìn)行噪聲處理；3）重建處理后的各IMF分量，得到去噪后信號(hào)。

利用EMD變換進(jìn)行脈沖干擾處理的關(guān)鍵是如何確定消除脈沖干擾噪聲的瞬時(shí)頻率閥值。

對(duì)于第j階固有模態(tài)函數(shù)，消除脈沖干擾噪聲的瞬時(shí)頻率閥值為

其中：N為IMF序列中所包含的數(shù)據(jù)數(shù)；k＝1或2，分別對(duì)應(yīng)于數(shù)據(jù)的Gaussian分布和Rayleigh分布；是第j階IMF的總體信號(hào)水平，定義為

其中：IMFj（i）為第j階IMF數(shù)據(jù)中第i個(gè)數(shù)據(jù)點(diǎn)的值，在整個(gè)時(shí)間域內(nèi)取值；median[IMFj（i）]為第j 階 IMF 的 中 位 差；median｛｜IMFj（i）-median[IMFj（i）]｜｝是第j 階IMF 的絕對(duì)中位 差（MAD），也稱為絕對(duì)中值偏差，相對(duì)于均方差來說，它對(duì)于脈沖干擾這種突起點(diǎn)，能夠進(jìn)行很好的抑制。σk是一個(gè)具有統(tǒng)計(jì)意義的固定常數(shù)：如果信號(hào)屬于Gaussian分布，則σ＝0.6745；如果信號(hào)屬于Rayleigh分布，σ＝0.4485。

式中：IMF′j（i）為閥值去噪后的IMFj（i）；w（i）為超出閥值部分的脈沖干擾進(jìn)行加權(quán)處理的權(quán)函數(shù)。當(dāng)?shù)趈階IMF數(shù)據(jù)的絕對(duì)值超過閥值時(shí)，乘以權(quán)函數(shù)w（i），使得其值迅速減?。划?dāng)數(shù)據(jù)的絕對(duì)值小于等于閥值時(shí)，保持不變。

Thomson函數(shù)具有嚴(yán)格的自適應(yīng)性的穩(wěn)健性，筆者選擇Thomson函數(shù)作為權(quán)函數(shù)：

3 實(shí)例

圖2 28號(hào)測(cè)點(diǎn)中頻段受脈沖干擾信號(hào)EMD分解各階IMF圖Fig.2 The intrinsic mode functions image of the medial signal interferenced by pulse at point 28

圖3 受脈沖干擾信號(hào)EMD除噪處理前后效果對(duì)比圖Fig.3 Comparison diagram of the interference signal processed by EMD transform

圖4 EMD除噪處理前（a）后（b）所得視電阻率（ρa(bǔ)）對(duì)比圖Fig.4 Comparison diagram of the resistivity before（a）and after（b）processing by EMD transform

將上述除噪方法應(yīng)用于白容礦區(qū)EH-4實(shí)測(cè)剖面第28號(hào)測(cè)點(diǎn)中頻段的Ex數(shù)據(jù)脈沖干擾處理，以驗(yàn)證該方法的有效性和合理性。將該數(shù)據(jù)進(jìn)行EMD分解，共得到12階IMF以及一階Res（圖2）。由于IMF1頻率過高，表現(xiàn)為信號(hào)所包含的白噪聲；Res為低頻趨勢(shì)項(xiàng)。采用式（5）分別計(jì)算各階IMF序列的閥值，并用此閥值對(duì)各階數(shù)據(jù)值進(jìn)行修正，對(duì)超出的脈沖部分采用式（7）進(jìn)行加權(quán)處理，最后得到消噪后的信號(hào)（圖3）。可以看出，原始數(shù)據(jù)在時(shí)間域明顯受到了脈沖干擾（圖3a）；經(jīng)過消噪處理后的時(shí)間序列，很顯然對(duì)脈沖干擾的抑制和消除效果顯著（圖3b）；圖3c是將消噪處理后的數(shù)據(jù)選擇一個(gè)合適的比例尺來表現(xiàn)，這樣可以將被脈沖干擾壓制的信號(hào)細(xì)節(jié)信息展示出來。經(jīng)過仔細(xì)對(duì)比，可以發(fā)現(xiàn)，這種除噪方法不僅消除了脈沖干擾的影響，而且保留甚至更細(xì)致地刻畫出信號(hào)的局部信息。

將經(jīng)過EMD除噪后的信號(hào)重新輸入EH-4數(shù)據(jù)處理系統(tǒng)，經(jīng)過計(jì)算所得的視電阻率曲線與未處理數(shù)據(jù)形成了鮮明的對(duì)比（圖4）。未處理的信號(hào)由于脈沖干擾的影響，在相位和視電阻率上經(jīng)常出現(xiàn)較大間斷和數(shù)值的突然增加；而經(jīng)過除噪后，整個(gè)信號(hào)顯得更加連續(xù)和可靠，有效地提高了數(shù)據(jù)的可信度和質(zhì)量水平。

4 結(jié)論

1）EMD變換可以分解出多階IMF序列，屬于一種從高頻到低頻的層層篩分過程，可以將信號(hào)分為若干個(gè)不同的具有各自物理意義的序列，是一種多尺度、多分辨分析的濾波過程。

2）利用EMD變換進(jìn)行大地電磁測(cè)深（MT）信號(hào)脈沖類干擾去除的關(guān)鍵是如何選擇一個(gè)合適的瞬時(shí)頻率閥值，對(duì)MT信號(hào)經(jīng)EMD分解后得到N個(gè)固有模態(tài)函數(shù)IMF進(jìn)行噪聲篩選，然后進(jìn)行EMD重構(gòu)。筆者提出的利用絕對(duì)中位差結(jié)合Thomson權(quán)函數(shù)對(duì)各階IMF進(jìn)行噪聲篩選是可行、有效的。

3）實(shí)測(cè)數(shù)據(jù)處理表明，利用EMD變換進(jìn)行大地電磁測(cè)深（MT）信號(hào)脈沖類干擾處理，改正前后信號(hào)的相關(guān)性高，能量損失小，能有效消除MT信號(hào)中脈沖類干擾。

（References）：

[1]化希瑞，曹哲明，劉鐵，等.高頻大地電磁系統(tǒng)數(shù)據(jù)處理方法研究[J].工程地球物理學(xué)報(bào)，2009，6（增刊）：25-32.Hua Xirui，Cao Zheming，Liu Tie，et al.The Research on High-Frequency Magnetotelluric System’s Data Processing Technique[J].Chinese Journal of Engineering Geophysics，2009，6（Sup.）：25-32.

[2]孫潔，晉文光，白登海，等.大地電磁測(cè)深資料的噪聲干擾[J].物探與化探，2000，22（2）：120-127.Sun Jie，Jin Wenguang，Bai Denghai，et al.The Noise Interference of Magnetotelluric Sounding Data[J].Geophysical and Geochemical Exploration，2000，22（2）：120-127.

[3]席振銖，龍霞，董晨，等.EH-4系統(tǒng)中工頻干擾的處理與改進(jìn)[J].地球物理學(xué)進(jìn)展，2010，25（3）：1105-1109.Xi Zhenzhu，Long Xia，Dong Chen，et al.An Improved Method to Eliminate the Power-Line Interference in the EH-4System[J].Progress in Geophysics，2010，25（3）：1105-1109.

[4]蔡劍華，湯井田，王先春.基于經(jīng)驗(yàn)?zāi)B(tài)分解的大地電磁資料人文噪聲處理[J].中南大學(xué)學(xué)報(bào)：自然科學(xué)版，2011，42（6）：1786-1790.Cai Jianhua，Tang Jingtian，Wang Xianchun.Human Noise Elimination for Magnetotelluric Data Based on Empirical Mode Decomposition[J].Journal of Central South University：Science and Technology，2011，42（6）：1786-1790.

[5]嚴(yán)家斌，劉貴忠.基于小波變換的脈沖類電磁噪聲處理[J].煤田地質(zhì)與勘探，2007，35（5）：61-65.Yan Jiabin，Liu Guizhong.Like-Impulse Electromagnetic Noise Processing Based Wavelet Transform[J].Coal Geology ＆ Exploration，2007，35（5）：61-65.

[6]Trad O D，Travassos M J.Wavelet Filtering of Magnetotelluric Data[J].Geophysics，2000，65（2）：482-491.

[7]Huang E N.The Empirical Mode Decomposition and Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis[M].London：Proc R Soc，1998.

[8]Huang E N，Wu Zhaohua.A Review on Hilbert-Huang Transform：The Method and Its Applications on Geophysical Studies[J].Rev Geophys，2008，46：1-23.

[9]Huang E N，Shen S S P.The Hilbert-Huang Transform and Its Applications[M].[S.l.]：World Scientific Publishing Co Pte Ltd，2005.

[10]湯井田，化希瑞，曹哲民，等.Hilbert-Huang變換與大地電磁噪聲壓制[J].地球物理學(xué)報(bào)，2008，51（2）：603-610.Tang Jingtian，Hua Xirui，Cao Zhemin，et al.Hilbert-Huang Transformation and Noise Suppression of Magnetotelluric Sounding Data[J].Chinese Journal of Geophysics，2008，51（2）：603-610.

[11]董烈乾，李振春，劉磊，等.基于經(jīng)驗(yàn)?zāi)B(tài)分解的曲波閾值去噪方法[J].吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2012，42（3）：838-844.Dong Lieqian，Li Zhenchun，Liu Lei，et al.A Method of Curvelet Threshold Denoiosing Based on Empirical Mode Decomposition[J].Journal of Jilin University：Earth Science Edition，2012，42（3）：838-844.

[12]陳凱.基于經(jīng)驗(yàn)?zāi)Ｊ椒纸獾娜ピ敕椒╗J].石油地球物理勘探，2009，44（5）：603-608.Chen Kai.A New Denoising Method Based on Empirical Mode Decomposition（EMD）[J].Oil Geophysical Prospecting，2009，44（5）：603-608.