基于頻率域小波去噪的大地電磁信號工頻干擾處理

蔡劍華,李 晉

(1. 湖南文理學院信息研究所,湖南常德 415000； 2. 湖南師范大學物理與信息科學學院, 湖南長沙 410081)

?

基于頻率域小波去噪的大地電磁信號工頻干擾處理

蔡劍華1,李 晉2

(1. 湖南文理學院信息研究所,湖南常德 415000； 2. 湖南師范大學物理與信息科學學院, 湖南長沙 410081)

大地電磁測深(Magnetotelluric，MT)在油氣勘查中得到越來越多的應用，針對MT中日益嚴重的工頻干擾，從頻率域著手，結合小波閾值去噪方法，提出了基于頻率域小波去噪的MT信號工頻干擾處理方法。先對受噪的MT信號進行傅里葉變換，得到其實部和虛部，再用小波閾值去噪的方法對實部序列和虛部序列分別進行去噪處理，最后將去噪后的實部和虛部聯(lián)合，進行反傅里葉變換得到去噪后的信號。給出了去噪方法的原理、步驟，并用仿真信號和實測大地電磁信號驗證了其有效性。結果表明：頻率域小波閾值去噪的大地電磁信號工頻干擾處理方法是正確、有效的，能有效且自適應地壓制大地電磁信號中的工頻干擾，突出被工頻干擾淹沒了的有用信號的信息。

頻率域 小波閾值 大地電磁信號 去噪

Cai Jian-Hua, Li Jin. Suppression of power line interference on MT signals based on the frequency domain wavelet method[J]. Geology and Exploration, 2015, 51(2):0353-0359.

近年來大地電磁測深得以迅速地發(fā)展，廣泛應用于石油天然氣勘查(孫衛(wèi)斌等，2003)。在復雜山區(qū)的油氣勘探應用效果日益突出，特別是在地震勘探困難區(qū)，油氣電磁法勘探技術已成為首選的非地震勘探技術之一。油氣勘探中利用電磁測深法，可以查明斷裂，斷裂破碎帶，地質構造填圖，巖層孔隙率調查及火成巖地質體等方面(衛(wèi)平生等，2006; 王斌等，2014)。而大地電磁信號極其微弱，易受到噪聲的污染，尤其是近些年來，隨著國民經濟和現(xiàn)代化建設的發(fā)展，電磁背景噪聲日益嚴重，縱橫交錯的高壓線、鐵路電網(wǎng)、大型工業(yè)設備等造成的工頻干擾成了大地電磁信號檢測中難以避開的噪聲來源(王書明等，2004;郭曉東等，2009)。小波去噪、經驗模態(tài)分解、形態(tài)濾波等方法為抑制工頻干擾提供了有利的手段，取得了不錯的效果(嚴家斌等，2008；湯井田等, 2008)。尤其是小波方法，自上世紀80年代應用于地球物理領域來得到了很大的發(fā)展，被成功應用在大地電磁信號的去噪、譜估計和靜態(tài)效應抑制等方面(Tradetal, 2000; 嚴家斌等，2008)。常規(guī)的小波去噪，通常是在時間域對信號進行分解，再對分解得到的系數(shù)進行閾值去噪，最后再重構達到去噪的目的(徐宏斌等，2012;陳濤等,2013)。本文從頻率域著手，提出了基于頻率域小波去噪的大地電磁信號工頻干擾處理方法，先對受噪的MT信號進行傅里葉變換，得到其實部和虛部，再用小波閾值去噪的方法對實部序列和虛部序列分別進行去噪處理，最后將去噪后的實部和虛部聯(lián)合，進行反傅里葉變換得到去噪后的信號。試圖為大地電磁信號去噪探索新的技術手段。

1 基本原理

1.1 傅里葉變換和反變換

給定實自變量t的非周期函數(shù)f(t)，可以是實函數(shù)，也可以是復函數(shù)，做積分

(1)

若上式對參數(shù)ω的任何實數(shù)值都存在， 則稱F(ω)為f(t)的傅里葉變換， 或稱傅里葉積分。f(t)可用F(ω)的積分表示，稱為傅里葉逆變換。

(2)

離散時間信號x(n)的連續(xù)傅立葉變換定義為

(3)

式中X(ejω)是一個連續(xù)函數(shù)，不能直接在計算機上做數(shù)字運算。為了在計算機上實現(xiàn)頻譜分析，必須對x(n)的頻譜作離散近似。有限長離散信號的離散傅立葉變換定義為

(4)

式中ωN=exp(-j2π/N)，其反變換定義為

(5)

工程上都通過構造離散傅里葉變換的快速算法，即快速離散傅里葉變換來實現(xiàn)傅氏變換。傅里葉變換的振幅譜、相位譜和能量譜分別為(孫麗穎等，2005)

(6)

(7)

(8)

式中R(ω)，I(ω)分別表示傅里葉變換的實部和虛部。

1.2 小波閾值去噪

對于連續(xù)的情況，小波為(Mallatetal,1992; Donohoetal,1994)

(9)

式中：Ψa,b(t)為一個基本小波函數(shù)或母小波函數(shù)，a為伸縮(尺度) 因子；b為平移因子。任給函數(shù)f(t) ∈L2(R)，其連續(xù)小波變換及重構公式(逆變換) 為

(10)

(11)

對于離散情況，小波為(Donohoetal,1994,1995)

(12)

與函數(shù)f(t) 相應的離散小波變換系數(shù)及其重構公式為

(13)

(14)

式中：C是一個與信號無關的常數(shù)。根據(jù)小波變換的Mallat算法，對信號其進行分解，得到不同尺度上的小波系數(shù)Cj,k，再對Cj,k進行閾值處理，然后再進行重構得到新信號。本文采用常用的軟閾值算法如下(Donohoetal,1995; Johnstoneetal,1997):

(15)

(16)

其中，N為信號長度，σ為噪聲的方差(Donohoetal,1995)

σ=MADj/0.6745

(17)

MADj為小波分解后第j層高子帶系數(shù)的中值。

2 去噪流程

基于頻率域小波去噪的MT信號工頻干擾處理流程如圖1所示：(1) 通過Fourier變換將原始含噪的大地電磁信號變換到頻率域中，得到信號傅里葉變換的實部和虛部；(2) 對實部和虛部分別進行小波變換，并對各層系數(shù)進行小波軟閾值處理；(3) 將處理后的實部和虛部分量組合成新的復數(shù)序列；(4) 最后通過反Fourier變換將此序列轉換至時間域，達到在頻率域用小波閾值方法抑制大地電磁信號噪聲的目的。

3 仿真處理

仿真信號用蒙特卡洛方法產生互不相干的一組隨機信號，噪聲為類似實際信號采集中常見工頻干擾的正弦信號，幅值為仿真信號最大值的1.2倍，頻率為50 Hz。原始仿真信號和加噪后仿真信號的時域圖如圖2(a)和圖2(b)所示。顯然，由于正弦噪聲干擾的加入，仿真信號完全被覆蓋，整體表現(xiàn)為正弦波的形態(tài)，從時域無法看出信號的細節(jié)特征。加噪前后仿真信號的統(tǒng)計參數(shù)為：加噪前信號的方差為480.7422，能量為4.8069×105mV2；加噪聲后信號方差為3.630×103，能量為3.6269×106mV2，信噪比為-8.9075 dB。圖3(a)、圖3(b)、圖3(c)分別給出了仿真信號、所加的噪聲信號和加噪后信號的時頻譜圖和邊際譜圖。從時頻譜圖可以看到，噪聲的能量非常強，在圖3(c)中，加入的噪聲能量幾

圖1 頻率域小波去噪的流程圖Fig.1 Flow chart of de-noising based on frequency domain Wavelet

乎淹沒了其他頻率成分的信息。

圖2 去噪前后的仿真信號Fig.2 Simulated signal of before and after de-noisinga-無噪仿真信號; b-加噪后仿真信號; c-去噪后仿真信號a-original signal; b-noised signal; c-de-noised signal

按照本文提出的方法，對含噪信號進行去噪處理，先對加噪后的信號進行快速傅里葉變換，得到其實部和虛部，再用小波閾值去噪的方法對實部序列和虛部序列分別進行去噪處理(本文采用小波軟閾值去噪方法，閾值函數(shù)如公式(15)，閾值如公式(16)所示)，最后將去噪后的實部和虛部聯(lián)合，進行反傅里葉變換得到去噪后的信號。去噪聲后信號的方差為482.0253，能量為4.8150×105mV2，信噪比為16.5859 dB?？梢?，噪聲能量得到了有效的抑制，信號能量由去噪前的3.6269×106mV2降到了4.8150×105mV2，而去噪后信號的能量基本與原始信號能量4.8069×105mV2，保持相當。圖2(c)為信號去噪后的時域圖，圖3(d)為信號去噪后的時頻圖，我們分別將他與其對應的圖2(a)和圖3(a)比較，可以看出信號和頻譜形態(tài)基本恢復了原樣。顯然，從圖譜和統(tǒng)計參數(shù)上都可以看出，本文提出的基于頻率域小波去噪的方法能有效去除50 Hz的干擾信號。

4 實測信號處理

周期正弦噪聲往往由觀測點附近高壓線或其他大型工業(yè)電器設施產生，由于強度大，在時間序列主要表現(xiàn)為正弦波的形式，它幾乎完全掩蓋了電磁信號，形成擬等振幅的50 Hz 噪聲干擾。具有如下表達形式

(18)

其Fourier 變換形式為：

(19)

實測數(shù)據(jù)是來自安徽某礦區(qū)某測點的EH-4數(shù)據(jù)，作業(yè)區(qū)離礦區(qū)很近，有多條高壓電線跨過。圖4所示為典型的受工頻干擾大地電磁信號的時間序列，從圖4可以看出：在強正弦噪聲干擾下，天然電磁場信號完全被掩蓋，在時間序列上很難識別出天然電磁信號的存在。為作去噪前后的對比，先給出信號去噪前的統(tǒng)計參量為：Hy分量最大值為1.4676×104nT，最小值-1.4152×104nT，方差為3.3915×107，能量為1.3911×1011nT2；Ex分量最大值為 631 mV，最小值-770 mV，方差為5.1922×104，能量為2.2192×108mV2；Hx分量最大值為1.2560 ×104nT，最小值-1.2753×104nT，方差為2.7146×107，能量為1.1123×1011nT2；Ey分量最大值為4.998 ×103mV，最小值-5.236×103mV，方差為3.5392×106，能量為1.4493×1010mV2。由于極化方向的原因，這類干擾可能在某個方向上強度較大，在另一個方向上強度較小，對單一主頻的噪聲往往只對其主頻和諧波頻率上的阻抗計算產生較大的影響。

圖3 去噪前后仿真信號的時頻譜圖Fig.3 Time-frequency spectrum of simulated signal of before and after de-noisinga-無噪仿真信號；b-所加噪聲；c-加噪后仿真信號；d-去噪后仿真信號a-Original signal；b-Noise；c-Noised signal；d-De-noised signal

圖4 實測大地電磁信號時域圖Fig.4 Time-frequency spectrum of measured MT signal

圖5所示為4道信號去噪前后所對應的功率譜。功率譜能很好地區(qū)分出能量隨頻率的細微變化，能反應大地電磁信號每頻段的頻率特性和能量差異。由圖5(左)去噪前的功率譜和時頻譜可以清晰的看到，50 Hz頻點及其諧波成分處的信號能量極大，是其他頻率成分的好幾倍，幾乎掩蓋了其他信號的信息。為作對比，圖5(右)給出了傳統(tǒng)小波閾值去噪的結果，從圖中可以看出，50 Hz工頻干擾的幅值很大程度上得到了有效抑制，但還能看到此頻點出突出的幅值。采用本文提出的去噪方法，對各道信號進行去噪處理。圖5(中)所示為4道信號去噪后所對應的功率譜。與圖5左右兩邊的譜圖對比，可以看出：本文方法去噪后50 Hz的工頻干擾及其低頻段部分工頻干擾的諧波成分也得到了有效的抑制，且本文提出的去噪方法也較傳統(tǒng)的小波閾值去噪效果好。

圖6所示為4道信號去噪后信號的時域圖，消噪后信號的統(tǒng)計參量：Hy分量最大值為2.8016×103nT，最小值為-2.6927×103nT，方差為8.6367×105，能量為3.5368×109nT2；Ex分量最大值為173.7995 mV，最小值為-156.8419 mV，方差為1.4768×103，能量為6.0478×106mV2；Hx分量最大值為2.9489×103nT，最小值為-2.5824×103nT，方差為8.2590×105，能量為3.3821×109nT2；Ey分量最大值1.7922×103mV，最小值為-1.3696×103mV，方差為 2.1610×105，能量為 8.8495×108mV2。對比去噪前后的統(tǒng)計參數(shù)，4道信號的方差和能量都減小了，磁道信號方差和能量強度下降了2個數(shù)量級，電道信號方差和能量強度下降了1個數(shù)量級，這也與工頻干擾主要影響磁道信號的事實相符。由去噪前后的統(tǒng)計參數(shù)，結合前面分析的功率譜，可以看到：消噪后信號顯得較為平穩(wěn)，看不到周期波形，信號形態(tài)符合大地電磁信號的規(guī)律，突出了被噪聲掩蓋的有用信息，這就為后續(xù)的功率譜估計和阻抗估算提供了條件。

圖5 去噪前后大地電磁信號的功率譜(左為去噪前的，中為本文方法去噪后的，右為傳統(tǒng)小波方法去噪后的)Fig.5 Power spectrum of measured signal of before and after de-noising(the left is before de-noising, the middle is after de-noising based on presented method, the right is after de-noising based on traditional Wavelet method)a-Hy分量；b-Ex分量；c-Hx分量；d-Ey分量；a-Hy component；b-Ex component；c-Hx component；d-Ey component

進一步，用去噪前后的信號計算其功率譜，并估算了視電阻率曲線和相位曲線，結果如圖7所示。

去噪前后響應曲線的整體趨勢是一致的，但細節(jié)上發(fā)現(xiàn)了明顯的變化。由去噪前的曲線可以看到，由于受到強烈的工頻干擾，在50 Hz頻點處，視電阻率曲線和相位曲線參數(shù)都存在飛點，而去噪后響應曲線變得平滑，曲線突變點得到有效抑制，MT測深資料的質量得到了提高了。

圖6 消噪后實測大地電磁信號的時域圖Fig.6 Time domain wave of measured MT signal after de-noising

5 結論

日益嚴重的工頻干擾給大地電磁信號檢測帶來了麻煩，這類干擾可能在某個方向上強度較大，在另一個方向上強度較小，對視電阻率和相位等MT響應參數(shù)的計算產生較大影響。本文從頻率域著手，結合小波閾值去噪方法，提出了基于頻率域小波去噪的大地電磁信號工頻干擾的抑制方法，仿真和實測信號的處理驗證了方法的有效性。結果表明本文方法突出了被工頻干擾淹沒了的有用信號的信息，這就為后續(xù)的功率譜估計和阻抗估算提供了條件，為大地電磁信號去噪提供了新的技術手段。

圖7 消噪前后的大地電磁響應曲線對比Fig.7 Comparison of MT response curve before and after de-noisinga-ρxy曲線；b-ρyx曲線；c-Фxy曲線；d-Фyx曲線a-ρxy curve；b-ρyx curve；c-Фxy curve；d-Фyx curve

Cai Jian-hua, Wang Xian-chun, Hu Wei-wen.2013.De-noising of MT signal based on empirical mode decomposition and wavelet threshold method[J]. Geophysical Prospecting, 2013, 48(2):303-307 (in Chinese with English abstract)

Chen Tao, Wang Wei, Ji-qing. 2013. De-noising method of fiber gyro signal based on wavelet transform[J] . Ordnance Industry Automation, 32(3):70-71,82( in Chinese with English abstract)

Donoho D L.1995. De-noising by soft-thresholding[J]. IEEE Transaction on Information Theory, 41(3):613-627

Donoho D L, Johnstone I M. 1994. Ideal spatial adaptation via Wavelet shrinkage[J]. Biometrika, 81:425-455

Johnstone L M and Silverman B W.1997.Wavelet threshold estimators for data with correlated noise[J]. J R Statist Soc,59(B):319-351

Guo Xiao-dong, Chen Xiao-qiang, Wang Zhi-hua, Cheng Rui-lin, Fan Zhan-jun, Xu De-li. 2009. Application of EH4 continuous conductivity measure in Baoxingchang district[J]. Geology and Prospecting, 45(1):52-58(in Chinese with English abstract)

Sun Li-yin, Qu D, Yan T. 2005. Application of Fourier Transform and Wavelet Transform to Signal Fault Diagnosis[J]. Journal of Liaoning Institute of Technology, 25(3): 155-157( in Chinese with English abstract)

Mallat S, Hwang W L.1992.Singularity detection and processing with Wavelet[J]. IEEE Transation on Information Theory,38(2):617-643

Sun Wei-bin, Song Qun-hui, Zheng Li. 2003. The development of MT sounding technique and its application to the oil and gas exploration[J].Geology and prospecting, 39(10): 1-7 ( in Chinese with English abstract)

Trad D O, Travassos J.M. 2000. Wavelet filtering of magnetotellurie data[J]. Geophysics, 65: 482-491

Tang Jin-tian, Hua Xi-rui, Cao Zhe-min. 2008. Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data[J]. Chinese J. Geophys, 51(2): 603-610(in Chinese with English abstract)

Wei Pin-sheng, Tan Kai-jun, Zhang Hu-quan. 2006. Application of magnetotelluric sounding to oil and gas exploration: Take the Minghe basin as an example[J]. Natural gas industry, 26 (8):37-40 (in Chinese with English abstract)

Wang Bin,Chen Xiao-qiang, Fan Zhan-jun, Chen Rui-lin. 2014. Application of EH4 continuous conductivity meter to the study of ore-controlling structures in the Jinlongshan gold deposit of Zhen’an country,Shaanxi Province [J]. 50(3): 0564-0571 (in Chinese with English abstract)

Wang Shu-min, Wang Jia-yin. 2004. Analysis on statistic characteristics of Magnetotelluric signal[J]. Acta Seismologica Sinica, 26 (6):669-674( in Chinese with English abstract)

Xu Hong-bin, Li Shu-lin, Chen Ji-jing. 2012. A study on method of signal de-noising based on wavelet transform for micro-seismicity monitoring in large-scale rockmass structures[J]. Acta Seismologica Sinica, 34(1): 81-96( in Chinese with English abstract)

Yan Jia-bin, Liu Gui-zhong, Liu Jian-xin. 2008. Application of wavelet transform in processing nature electromagnetic field time series[J]. Geology and Prospecting, 44(3): 75-78( in Chinese with English abstract)

[附中文參考文獻]

蔡劍華, 王先春, 胡惟文. 2013. 基于經驗模態(tài)分解與小波閾值的MT信號去噪方法[J]. 石油地球物理勘探,48(2):303-3070

陳 濤,王 偉,吉 清. 2013. 基于小波變換的光纖陀螺信號去噪方法[J]. 兵工自動化, 32(3):70-71,82

郭曉東,陳孝強,王治華,陳瑞林,樊占軍,徐德利. 2009. EH4 連續(xù)電導率測量在寶興廠礦區(qū)的應用[J]. 地質與勘探,45( 1) : 52-58

徐宏斌, 李庶林, 陳際經. 2012. 基于小波變換的大尺度巖體結構微震監(jiān)測信號去噪方法研究[J]．地震學報,34(1):81-96

孫衛(wèi)斌, 宋群會, 鄭 莉. 2003.大地電磁測深技術發(fā)展及在油氣勘探的應用[J]. 地質與勘探, 39(10):1-7

孫麗穎, 屈 丹, 閆 鈿. 2005. 傅里葉變換與小波變換在信號故障診斷中的應用[J]. 遼寧工學院學報, 25 (3):155-157

湯井田, 化希瑞, 曹哲民. 2008. Hilbert-Huang變換與大地電磁噪聲壓制[J].地球物理學報, 51(2):603-610

衛(wèi)平生,譚開俊,張虎權. 2006.大地電磁測深在油氣勘探中的應用—以民和盆地為例[J]. 天然氣工業(yè), 26 (8):37-40

王 斌, 陳孝強, 樊戰(zhàn)軍,陳瑞林. 2014. EH4 連續(xù)電導率儀在陜西省鎮(zhèn)安縣金龍山礦區(qū)控礦構造研究中的應用[J]. 地質與勘探, (50)3: 0564-0571

王書明,王家映. 2004. 大地電磁信號統(tǒng)計特征分析[J].地震學報, 26(6) :669-674

嚴家斌, 劉貴忠, 柳建新. 2008.小波變換在天然電磁場信號時間序列處理中的應用[J]. 地質與勘探, 44(3): 75-78

Suppression of Power Line Interference on MT Signals Based on the Frequency Domain Wavelet Method

CAI Jian-Hua1, LI Jin2

(1.InstituteofInformation,HunanUniversityofArtsandScience,Changde,Hunan415000; 2.InstituteofPhysicsandInformationScience,HunanNormalUniversity,Changsha,Hunan410081)

The magnetotelluric(MT)method is widely applied in oil and gas exploration. While power line interference is an increasingly serious problem with the development of society. In the frequency domain, combined with the wavelet threshold method, an approach to restrain power line interference is proposed for MT signal processing. Firstly, the Fourier transform is applied to the noised MT signal and its real and imaginary parts can be obtained. Then, wavelet threshold method is used to de-noise for the imaginary and the real part sequences. Finally, the de-noised real and the de-noised imaginary are assembled, and the anti-Fourier transform is conducted and the de-noised signal is obtained. The principle and steps of de-noising methods are presented. The simulated signals and the measured MT data are processed to verify its validity. The results show that the proposed method is correct and effective. It can suppress the power line interference on MT signal effectively and adaptively, and the useful information is highlighted which is concealed by the power line interference.

Frequency domain, Wavelet threshold, Magnetotelluric signal, De-noising

2014-09-26；[修改日期]2014-12-28；[責任編輯]陳偉軍。

國家自然科學基金項目(41304098，41404111), 湖南省自然科學基金項目(12JJ4034), 湖南省教育廳青年項目(13B076), 湖南省重點建設學科—光學基金, 湖南省重點實驗室“光電信息集成與光學制造技術”,“湖南省光電信息技術校企聯(lián)合人才培養(yǎng)基地”共同資助。

蔡劍華(1979年-), 男, 2010年畢業(yè)于中南大學，獲博士學位, 副教授, 從事大地電磁數(shù)據(jù)處理的研究。E-mail:cjh1021cjh@163.com。

P361.4

A

0495-5331(2015)02-0353-07