凌振寶， 王沛元， 萬云霞*， 王言章， 程德福， 李桐林

1 吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院， 長春 130021 2 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長春 130021

?

強(qiáng)人文干擾環(huán)境的電磁數(shù)據(jù)小波去噪方法研究

凌振寶1， 王沛元1， 萬云霞1*， 王言章1， 程德福1， 李桐林2

1 吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院， 長春 130021 2 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長春 130021

大地電磁法(MT)以成本低廉，探測(cè)深度大、水平方向分辨能力高等優(yōu)點(diǎn)，在礦產(chǎn)資源勘探方面得到廣泛應(yīng)用.然而，在老礦區(qū)或者礦區(qū)周圍進(jìn)行二次探礦時(shí)，強(qiáng)人文干擾嚴(yán)重影響觀測(cè)數(shù)據(jù)的質(zhì)量，導(dǎo)致反演結(jié)果出現(xiàn)偏差，甚至出現(xiàn)錯(cuò)誤的解釋結(jié)果.因此，需要對(duì)觀測(cè)數(shù)據(jù)進(jìn)行降噪處理.結(jié)合多分辨率分析算法和小波閾值算法的特點(diǎn)，本文提出了綜合小波算法：采用db3小波基；基于多分辨率分析算法，去除長周期噪聲；基于小波閾值算法，將Bayes估計(jì)配合改進(jìn)型閾值函數(shù)去除短周期噪聲干擾.對(duì)實(shí)測(cè)數(shù)據(jù)的處理結(jié)果顯示，處理后的數(shù)據(jù)的時(shí)間序列以及視電阻率曲線質(zhì)量都有了明顯的改善，近源效應(yīng)得到有效的抑制.

人文干擾； 多分辨率分析； 小波閾值； 視電阻率； 近源效應(yīng)

1 引言

大地電磁法(MT)是一種探測(cè)地下電性結(jié)構(gòu)的重要方法，觀測(cè)數(shù)據(jù)的質(zhì)量對(duì)后續(xù)資料的處理解釋至關(guān)重要.隨著人類社會(huì)的發(fā)展，工業(yè)化程度越來越高，人文噪聲充斥在人類生活的方方面面.大地電磁數(shù)據(jù)由于頻帶寬、信號(hào)弱，極易受到周圍環(huán)境噪聲的影響，使得采集到的信號(hào)摻雜著多種電磁干擾噪聲，反演結(jié)果存在很多不確定性.

通過觀察分析九瑞礦集區(qū)電磁數(shù)據(jù)的時(shí)頻曲線，可將影響MT數(shù)據(jù)質(zhì)量的人文干擾歸結(jié)為兩類：長周期干擾(方波噪聲、階躍噪聲等)，短周期干擾(脈沖噪聲，似充放電衰減模式的三角波噪聲等)(范翠松，2009).噪音主要影響MT信號(hào)的中低頻段，且主要存在于電場(chǎng)信號(hào)中，由于干擾源距離測(cè)點(diǎn)較近，其產(chǎn)生的干擾信號(hào)不符合平面波的要求，因而被稱之為近場(chǎng)或近源干擾(王剛等，2015).該類干擾在時(shí)間域信號(hào)中通常表現(xiàn)為信號(hào)振幅很大(幅值可能比正常MT信號(hào)的幅值大幾個(gè)數(shù)量級(jí))；在視電阻率曲線上表現(xiàn)為45°上升；而相應(yīng)的相位曲線則接近于0°.通過分析九瑞測(cè)區(qū)MT實(shí)測(cè)數(shù)據(jù)的視電阻率曲線后發(fā)現(xiàn)，其形態(tài)與上述近源干擾的視電阻率形態(tài)相近，判斷為受人文干擾而造成的近源效應(yīng).

近些年來，國內(nèi)外的一些學(xué)者在處理大地電磁信號(hào)的噪聲問題上做了很多努力.Gamble T D，Goubau W M等人先后在1978年和1979年提出了互功率譜法(Goubau et al., 1978)和遠(yuǎn)參考大地電磁測(cè)深法(Gamble et al., 1979)，用于消除不相關(guān)的電磁噪聲以及同源相關(guān)的電、磁噪聲；1998年，中國的嚴(yán)良俊、胡文寶利用遠(yuǎn)參考法處理南方碳酸鹽地區(qū)的MT強(qiáng)干擾噪聲后認(rèn)為，遠(yuǎn)參考法在較強(qiáng)的噪聲環(huán)境下收效甚微(嚴(yán)良俊等，1998；李桐林等，2009)；后來?xiàng)钌热私?jīng)過研究后也認(rèn)為遠(yuǎn)參考法無法壓制電場(chǎng)中的近源干擾(楊生等，2002)；Egbert等在1986年提出了基于Robust統(tǒng)計(jì)處理方法，用于消除非高斯正態(tài)分布的噪聲，此方法無法剔除輸入端噪聲，對(duì)受到較強(qiáng)近源干擾的電磁數(shù)據(jù)效果不明顯(Egbert and Booker, 1986)； Huang等在1998年提出了一種HHT方法(Huang et al., 1998)，可以有效抑制大地電磁信號(hào)中的工頻干擾，但由于EMD分解的自適應(yīng)性，該方式無法反映每時(shí)段數(shù)據(jù)特征的微妙變化(湯井田等，2008)；湯井田等人提出的數(shù)學(xué)形態(tài)濾波法可以很好地還原大地電磁信號(hào)的原始特征，但是由于結(jié)構(gòu)元素的選取及其尺寸必須通過反復(fù)實(shí)驗(yàn)獲得，其應(yīng)用推廣還需要進(jìn)一步的深入研究(湯井田等，2015，2012)；此外，吉林大學(xué)的范翠松等人提出的人機(jī)聯(lián)合去噪法在壓制強(qiáng)人文干擾方面效果明顯，但是操作過程中人為的經(jīng)驗(yàn)因素太多，不但耗時(shí)間耗力，也難以推廣應(yīng)用(范翠松，2009).

小波變換被稱為“數(shù)學(xué)顯微鏡”，在時(shí)間域和頻率域都具有良好的局部化特性，對(duì)于不同頻段的信號(hào)都可以更好地分析其中噪聲的特點(diǎn)，非常適合處理非平穩(wěn)信號(hào)和提取信號(hào)的局部特征.本文根據(jù)大地電磁信號(hào)特點(diǎn)，把多分辨分析算法和基于經(jīng)驗(yàn)貝葉斯估計(jì)的改進(jìn)閾值函數(shù)的小波閾值算法相結(jié)合，用于壓制觀測(cè)數(shù)據(jù)中的方波、脈沖等人文干擾，改善了數(shù)據(jù)質(zhì)量，獲得了較為理想的效果.

2 綜合小波去噪法原理

2.1 小波變換原理

(1)

其逆變換(重構(gòu)公式)為

(2)

(3)

對(duì)應(yīng)的逆變換為

(4)

小波系數(shù)Wf(j,k)是小波變換后用來表征信號(hào)特征信息的無量綱結(jié)果.對(duì)信號(hào)進(jìn)行小波分析，實(shí)質(zhì)上就是對(duì)變換后的小波系數(shù)進(jìn)行分析處理.

從(3)式中可以看出，在對(duì)信號(hào)f(t)進(jìn)行小波變換時(shí)，選擇不同的母小波函數(shù)ψ(t)，就會(huì)得到不同的小波系數(shù)組Wf(j,k).再由(4)式對(duì)處理之后的小波系數(shù)進(jìn)行逆變換(即重構(gòu))，從而得到處理后的信號(hào).

2.2 母小波函數(shù)選取

在對(duì)信號(hào)進(jìn)行小波變換時(shí)，由2.1節(jié)可知，無論是多分辨分析法還是小波閾值降噪法，選擇合適的母小波函數(shù)是進(jìn)行小波分析的基礎(chǔ).

在實(shí)際應(yīng)用中，分別用Daubechies(dbN)小波系、Coiflet(coifN)小波系和SymletsA(symN)小波系對(duì)大地電磁數(shù)據(jù)進(jìn)行多次分解重構(gòu)實(shí)驗(yàn)(其中N為消失矩).一般來說消失矩N越大對(duì)應(yīng)濾波器也就越平坦，但是小波分解高頻分量的零點(diǎn)越多，去除的有用信號(hào)就越多，所以消失矩也并不是越大越好.此外通過分析處理后信號(hào)的視電阻率曲線圖，又結(jié)合大地電磁數(shù)據(jù)信號(hào)的特點(diǎn)和小波函數(shù)的特性，本文最終選擇db3作為綜合小波去噪算法的小波函數(shù).2.3 基于多分辨分析去除低頻干擾

在礦集區(qū)進(jìn)行電磁探測(cè)時(shí)，由于周邊電力通訊網(wǎng)絡(luò)發(fā)達(dá)，電話線傳輸?shù)皖l模擬信號(hào)時(shí)會(huì)產(chǎn)生長周期方波噪聲，對(duì)電磁數(shù)據(jù)造成嚴(yán)重干擾(范翠松等，2008).而多分辨分析法(胡昌華等，2008)在小波分解的基礎(chǔ)上，將低頻部分進(jìn)一步分解，高頻部分不予考慮.分解層數(shù)越多，頻率的分辨率越高，得到的低頻小波系數(shù)更能反映低頻噪聲的性質(zhì)，適合于提煉與壓制低頻(如基線漂移、長周期方波等)噪聲.多分辨分析法分解原理如圖1所示.

本文選取的實(shí)驗(yàn)數(shù)據(jù)采樣率為24 Hz，根據(jù)奈奎斯特定理，有效電磁數(shù)據(jù)頻譜的頻帶范圍為0～12 Hz.通過觀察方波干擾、階躍噪聲等長周期干擾的時(shí)間域波形及相應(yīng)的頻譜數(shù)據(jù)，發(fā)現(xiàn)該類干擾的頻率主要集中在0.04～0.2 Hz.根據(jù)信號(hào)有效頻帶與采樣率之間的關(guān)系，確定分級(jí)層數(shù)取10層，可將低頻噪聲與有用信號(hào)分離.即采用db3小波對(duì)信號(hào)進(jìn)行10層分解，所得到的低頻層作為對(duì)長周期噪聲的估計(jì)，除去低頻層后將所有高頻部分相加得到校正后的信號(hào).

2.4 小波閾值降噪

除了2.3節(jié)中提到的長周期噪聲外，在數(shù)據(jù)采集過程中，由于周邊大型機(jī)械的突然開關(guān)也會(huì)產(chǎn)生脈沖等短周期高頻干擾，嚴(yán)重影響數(shù)據(jù)質(zhì)量.

圖1 多分辨分析原理圖cA代表信號(hào)的低頻小波系數(shù)；cD代表信號(hào)的高頻小波系數(shù)，LPF表示低通濾波器，HPF表示高通濾波器，Li(i=1,2,3…)表示第i層分解后的低頻分量，Hi(i=1,2,3…)表示第i層分解后的高頻分量.Fig.1 Elementary diagram of multiresolution analysiscA is approximation coefficient, cD is the detailed coefficients, LPF is low pass filters and HPF is high pass filters. Li is the low frequency component after decomposition at level i. Hi is the high frequency component after decomposition at level i.

對(duì)于短周期噪聲，高頻小波系數(shù)能夠更好地表征噪聲特性，采用小波閾值法對(duì)小波高頻噪聲進(jìn)行門限閾值處理， 抑制信號(hào)中的噪聲部分，達(dá)到降噪效果.

一個(gè)含有噪聲的一維信號(hào)的模型可以用以下式子進(jìn)行表示(Daubechies, 1990):

(5)

其中，f(i)是真實(shí)信號(hào)；e(i)是噪聲信號(hào)；s(i)是含噪信號(hào).

一般來說，一維信號(hào)的降噪過程大都可以分為3個(gè)步驟來進(jìn)行：

(1) 一維信號(hào)的小波分解.選擇一個(gè)母小波函數(shù)并確定分解層數(shù)對(duì)含噪信號(hào)s(i)進(jìn)行小波分解.

(2) 對(duì)小波分解的高頻系數(shù)進(jìn)行閾值量化.選擇一種閾值規(guī)則對(duì)各層小波系數(shù)進(jìn)行量化處理.

(3) 一維信號(hào)的小波重構(gòu).將處理后的小波系數(shù)進(jìn)行小波重構(gòu).

從降噪過程可以看出，保證小波閾值降噪法去噪效果的關(guān)鍵就是選擇適當(dāng)?shù)姆纸鈱訑?shù)和小波閾值規(guī)則.

2.4.1 分解層數(shù)選取

不同于多分辨分析法，小波閾值降噪法的分解層數(shù)并不是分解的層數(shù)越高，閾值去噪效果越好.因?yàn)樾〔ㄩ撝捣治鰧?shí)質(zhì)上是將信號(hào)通過多個(gè)濾波器，層數(shù)越多意味著濾波器越多，會(huì)造成信號(hào)的移位增多；同時(shí)也會(huì)引起信號(hào)的失真和能量的損失(李肅義等，2013)，從而會(huì)導(dǎo)致信號(hào)重構(gòu)時(shí)信息量的缺失.在實(shí)際操作時(shí)，可根據(jù)觀測(cè)數(shù)據(jù)的采樣率及噪聲所在頻段確定分解層數(shù)，將分解后的噪聲所在層的系數(shù)進(jìn)行重點(diǎn)處理，一般選擇尺度為3～6層.本文處理的數(shù)據(jù)采樣率為24 Hz，低于0.2 Hz的低頻干擾已經(jīng)通過多分辨率方法進(jìn)行了處理，因此層數(shù)確定為5層，將分解后高于0.375 Hz的高頻分量系數(shù)進(jìn)行閾值量化處理，以便壓制高頻干擾.

2.4.2 小波閾值函數(shù)選取

常用的小波閾值函數(shù)分為硬閾值和軟閾值(Donoho, 1995)，其公式分別如下:

硬閾值函數(shù):

(6)

軟閾值函數(shù):

(7)

本文采用一種改進(jìn)型的閾值方法(柯熙政等，2008)，公式如下:

(8)

改進(jìn)型的閾值函數(shù)減小了使用軟閾值所產(chǎn)生的恒定誤差，同時(shí)對(duì)較大的小波系數(shù)進(jìn)行緩慢壓縮，盡可能地保留了有用信號(hào)，從而減小了重構(gòu)之后信號(hào)與真實(shí)信號(hào)之間的偏差.

2.4.3 小波閾值規(guī)則選擇

小波閾值規(guī)則的選用，是確定各層閾值的關(guān)鍵.閾值過大則容易將有用信號(hào)去除，使噪聲數(shù)據(jù)失真，如若過小，則去噪效果不明顯.一般的閾值選擇規(guī)則有(周祥鑫等，2014)：基于Stein的無偏似然估計(jì)的SURE閾值，啟發(fā)式閾值，以及固定式閾值和極大極小原理的閾值等四種.由于Ebayes閾值規(guī)則是以經(jīng)驗(yàn)貝葉斯方法為原理(Jansen and Bultheel, 1999; Johnstone and Silverman, 2005; 劉偉等，2013)，可以更好地計(jì)算出平均平方差的最小值，所得閾值更接近最優(yōu)閾值，因此本文選用該閾值規(guī)則.

經(jīng)驗(yàn)貝葉斯估計(jì)Ebayes計(jì)算閾值T的步驟如下：

(1) 小波分解得到各層的小波系數(shù)ω(j,k).

(2) 根據(jù)以下公式估計(jì)各高頻層的噪聲方差：

(9)

(10)

(4)根據(jù)下式計(jì)算各層的閾值T：

(11)

2.5 綜合小波法降噪步驟

礦集區(qū)大地電磁信號(hào)中受到的人文噪聲既含有方波、基線漂移等長周期噪聲，又含有脈沖噪聲等短周期噪聲.結(jié)合2.3節(jié)和2.4節(jié)的分析結(jié)果，將多分辨分析法和小波閾值降噪法相結(jié)合對(duì)電磁數(shù)據(jù)中的人文干擾進(jìn)行壓制.首先載入實(shí)測(cè)信號(hào)f1(t)，采用多分辨分析法將低頻噪聲去除，得到信號(hào)f2(t)，根據(jù)選好的小波基函數(shù)以及確定的分解層數(shù)N對(duì)信號(hào)f2(t)再進(jìn)行小波分解，得到低頻系數(shù)和每一層的高頻系數(shù).采用經(jīng)驗(yàn)貝葉斯估計(jì)配合改進(jìn)的閾值函數(shù)進(jìn)行閾值的量化處理，得到第一層到第N層的高頻系數(shù)，最后進(jìn)行小波重構(gòu)，得到降噪后的信號(hào).如圖2所示為小波綜合降噪法的處理流程圖.

3 仿真數(shù)據(jù)處理

為了驗(yàn)證綜合小波降噪法在處理含人文干擾的MT數(shù)據(jù)方面的有效性，采用數(shù)值模擬的方法產(chǎn)生了含噪電磁數(shù)據(jù)作為被處理對(duì)象，在0.01～0.8 Hz的頻帶內(nèi)選取0.01 Hz、0.06 Hz、0.3 Hz、0.7 Hz四個(gè)頻率產(chǎn)生正弦波進(jìn)行疊加，作為MT信號(hào)的模擬.在MT模擬信號(hào)的不同時(shí)間段分別加入不同寬度的方波和脈沖，同時(shí)加入基線漂移，產(chǎn)生含有噪聲的MT仿真數(shù)據(jù).利用本文提出的方法對(duì)仿真含噪數(shù)據(jù)進(jìn)行降噪處理，通過處理前后數(shù)據(jù)的信噪比(SNR, Signal Noise Rate)與均方根誤差(MSE，Mean Square Error)的改變判斷處理效果.

圖2 綜合小波降噪法流程圖Fig.2 Flow chart of combined wavelet transform algorithm

圖3為模擬仿真數(shù)據(jù)的時(shí)間序列降噪前后的對(duì)比.其中，圖3a為未加干擾的MT模擬數(shù)據(jù)，圖3b為添加了方波和脈沖干擾的MT模擬數(shù)據(jù)，圖3c為經(jīng)過綜合小波降噪法處理后得到的MT時(shí)間序列.對(duì)比圖3b和3c可以看出，加入的方波和脈沖干擾得到明顯抑制，時(shí)間序列變得平穩(wěn)，連續(xù)性更好；對(duì)比圖3a和3c可以看出，含噪數(shù)據(jù)處理后，有用信號(hào)被很好地還原保留，波形無明顯變化，且無相位偏移，信號(hào)整體趨勢(shì)未發(fā)生明顯改變.證明該方法在處理MT數(shù)據(jù)時(shí)既可消除噪聲干擾，又可保證不損失有用信號(hào).表1分別給出了降噪前后模擬數(shù)據(jù)的SNR和MSE的對(duì)比.經(jīng)過降噪處理后，模擬MT數(shù)據(jù)的信噪比從-13.7535提升至2.6096，表明噪聲得到有效抑制；相應(yīng)的均方根誤差從48.7864降至1.0125，表明去噪后信號(hào)與原信號(hào)相似程度較高.

圖3 仿真信號(hào)時(shí)間序列去噪前后對(duì)比(a) 原始模擬信號(hào)；(b) 含噪模擬信號(hào)；(c) 降噪后模擬信號(hào).Fig.3 Comparison of simulated data time series de-noising before and after(a) Original simulated data; (b) Noisy simulated data; (c) De-noised simulated data.

表1 仿真數(shù)據(jù)降噪前后SNR和MSE對(duì)比Table 1 Comparison of SNR and MSE of simulated date

4 實(shí)測(cè)數(shù)據(jù)處理

本文處理的實(shí)測(cè)數(shù)據(jù)是來自九瑞礦集區(qū)實(shí)測(cè)MT數(shù)據(jù).選取一個(gè)距礦山較近的測(cè)點(diǎn)的數(shù)據(jù)作為處理對(duì)象，以此驗(yàn)證綜合小波降噪法的去噪有效性.如圖4所示，分別為處理前后的實(shí)測(cè)電磁數(shù)據(jù)的時(shí)間序列，采樣率為24 Hz(朱威等，2011)，顯示采樣時(shí)長為1000 s.

從圖4a中可以看出，處理前磁道和電道都存在零點(diǎn)漂移現(xiàn)象，同時(shí)，電道的信號(hào)存在嚴(yán)重的方波和大量脈沖噪聲，導(dǎo)致數(shù)據(jù)整體連續(xù)性差，有用信號(hào)被噪聲淹沒.處理后的電磁信號(hào)時(shí)間序列變得平穩(wěn)，零點(diǎn)漂移現(xiàn)象消失，方波、脈沖等人文噪聲得到了較好的壓制，電場(chǎng)和磁場(chǎng)的方差均有明顯降低，降噪效果顯著.去噪前后電場(chǎng)和磁場(chǎng)數(shù)據(jù)的最大/最小值及方差對(duì)比結(jié)果見表2.

圖4 實(shí)測(cè)信號(hào)時(shí)間序列去噪前后對(duì)比(a) 處理前的時(shí)間序列； (b) 處理后的時(shí)間序列.Fig.4 Comparison of measured data time series de-noising before and after(a) Before de-noising; (b) After de-noising.

圖5為y方向電場(chǎng)數(shù)據(jù)Ey處理前后的細(xì)節(jié)對(duì)比(所選時(shí)間段在550～590 s之間).

從圖5a中可以看出，在562～580 s時(shí)間段內(nèi)，信號(hào)受到一個(gè)較大能量的方波干擾，同時(shí)，此時(shí)間段內(nèi)信號(hào)一直受到強(qiáng)脈沖干擾，有用信號(hào)被淹沒，無法辨別.從頻率域分析得知，噪聲頻率有0.06 Hz，0.25 Hz，0.83 Hz，0.9 Hz，0.95 Hz，2.3 Hz等，主要分布在0.01～1 Hz之間.這些噪聲導(dǎo)致視電阻率曲線發(fā)生畸變，同時(shí)相位曲線接近0°.利用綜合小波降噪法對(duì)信號(hào)處理之后，時(shí)間序列如圖5b所示，方波干擾被剔除，強(qiáng)脈沖干擾得到了明顯壓制，信號(hào)波形變得平穩(wěn).

圖6為處理前后的視電阻率和相位曲線對(duì)比.

從圖6a—6b可以看出，在0.1 Hz左右，處理前視電阻率曲線圖呈45°上升，相位曲線接近于0°，這是典型的近源干擾現(xiàn)象(徐志敏等，2012).對(duì)于1 Hz以上的數(shù)據(jù)為中頻和高頻采樣率得到的數(shù)據(jù)，暫時(shí)不列為本文的研究對(duì)象，本文重點(diǎn)對(duì)頻率在0.8～0.01 Hz段內(nèi)的噪聲進(jìn)行降噪處理.從處理之后的視電阻率和相位曲線可以看出，該頻段內(nèi)的數(shù)據(jù)質(zhì)量得到明顯提升，近源干擾得到有效抑制.本文處理的數(shù)據(jù)采樣率為24 Hz，其參與視電阻率運(yùn)算的頻率成分最高到0.8 Hz左右，因此，圖6只給出了頻率低于0.8 Hz的頻點(diǎn).

為了檢驗(yàn)綜合小波降噪法在降噪過程中未對(duì)有用信號(hào)的影響，我們從同一測(cè)線上選擇了數(shù)據(jù)質(zhì)量較好的第273號(hào)測(cè)點(diǎn)進(jìn)行降噪處理，該測(cè)點(diǎn)離礦區(qū)較遠(yuǎn)，近源干擾現(xiàn)象不明顯.圖7(a,c，e)中所示為該點(diǎn)的時(shí)間序列及視電阻率和相位曲線.分析電磁場(chǎng)時(shí)間序列，發(fā)現(xiàn)存在較多的似充放電三角波波形，并且幅值較正常的大地電磁信號(hào)高.除此之外，整體的時(shí)間序列中無其他的非平穩(wěn)信號(hào)，可以表征天然大地電磁場(chǎng)的特征.從視電阻率曲線上看，除了受到似充放電三角波干擾而影響的0.1～0.6 Hz頻段的曲線有所跳變外，其他頻段的電阻率曲線形態(tài)光滑，誤差棒較小，無明顯近源現(xiàn)象存在.圖7(b,d，f)為采用本文方法處理后的273測(cè)點(diǎn)的時(shí)間序列及視電阻率和相位曲線.可以看出，1～0.1 Hz頻段內(nèi)方向的電阻率和相位曲線形態(tài)有明顯的改善，其他頻段的曲線形態(tài)無較大改變，說明本文方法可在損失較少有用信號(hào)的條件下剔除大部分的噪聲干擾.

表2 實(shí)測(cè)信號(hào)降噪前后的最大/最小值及方差比較Table 2 Compared of minimum/maximum and variance of measured data

圖5 電場(chǎng)數(shù)據(jù)Ey處理前后對(duì)比圖(時(shí)間段為550～590 s)(a) 處理前的時(shí)間序列； (b) 處理后的時(shí)間序列.Fig.5 Comparison of measured data in Ey de-noising before and after(a) Before de-noising; (b) After de-noising.

圖6 第208號(hào)測(cè)點(diǎn)視電阻率曲線變化(a) 處理前幅值曲線； (b) 處理后幅值曲線； (c) 處理前相位曲線；(d) 處理后相位曲線.Fig.6 Comparison of the curves of 208# point apparent resistivity(a) Magnitude before de-noising; (b) Magnitude after de-noising； (c) Phase before de-noising; (d) Phase after de-noising.

由于天然的大地電磁場(chǎng)是由不同的源激勵(lì)生成的，因此其極化方向是隨機(jī)變化的，隨著時(shí)間表現(xiàn)出無序性.如果區(qū)域中存在強(qiáng)烈的電磁干擾，電磁場(chǎng)方向則可能被這些主動(dòng)源所統(tǒng)治，密集地集中在某個(gè)方向上(張弛，2013)，表現(xiàn)出很強(qiáng)的方向性.Weckmann等提出用電磁場(chǎng)極化方向作為電磁場(chǎng)信噪識(shí)別的重要指標(biāo)(Weckmann et al., 2005)，其極化方向定義如下：

(12)

(13)

對(duì)208測(cè)點(diǎn)處理前后的電場(chǎng)數(shù)據(jù)進(jìn)行了極化圖的繪制.分別選取0.03 Hz和0.6 Hz頻率處的數(shù)據(jù)進(jìn)行計(jì)算，對(duì)比處理前后的電場(chǎng)數(shù)據(jù)極化方向的變化，如圖8、9所示.從圖中可以看出，處理前的電場(chǎng)極化方向角度非常集中，具有明顯的方向性，表示受近源干擾嚴(yán)重.處理后的電場(chǎng)數(shù)據(jù)極化方向相對(duì)分散，隨機(jī)性較大，符合天然場(chǎng)數(shù)據(jù)的極化特征，表示數(shù)據(jù)中的近源干擾得到了有效的抑制.

圖7 第273號(hào)測(cè)點(diǎn)的時(shí)間域波形以及視電阻率和相位曲線(a)、(c)、(e) 降噪前時(shí)間序列及視電阻率和相位曲線； (b)、(d)、(f) 降噪后時(shí)間序列及視電阻率和相位曲線.Fig.7 The curves of time series and apparent resistivity for 273# point(a)、(c)、(e) The curves of time series and apparent resistivity before de-nosing; (b)、(d)、(f) The curves of time series and apparent resistivity after de-nosing.

圖8 頻率為0.03 Hz時(shí)電場(chǎng)數(shù)據(jù)處理前后的極化方向(a) 處理前; (b) 處理后.Fig.8 The comparison chart of polarization direction at 0.03 Hz(a) Before de-noising; (b) After de-noising.

圖9 頻率為0.6 Hz時(shí)電場(chǎng)數(shù)據(jù)處理前后的極化方向(a) 處理前; (b) 處理后.Fig.9 The comparison chart of polarization direction at 0.6 Hz(a) Before de-noising; (b) After de-noising.

5 結(jié)論

天然的大地電磁信號(hào)頻率成分豐富，觀測(cè)時(shí)間長，其長期觀測(cè)的時(shí)間序列可以看成是一種低能量的平穩(wěn)隨機(jī)信號(hào)，很容易受到強(qiáng)電磁場(chǎng)源(如電氣化鐵路、電站、工廠、油田等人工設(shè)施)的干擾.本文根據(jù)小波變換的時(shí)頻高分辨特性，將小波變換中的多分辨率分析算法和基于經(jīng)驗(yàn)貝葉斯閾值函數(shù)的小波閾值算法相結(jié)合，提出了綜合小波降噪法.采用db3小波基，通過對(duì)模擬信號(hào)和實(shí)測(cè)電磁信號(hào)進(jìn)行降噪處理，結(jié)果顯示該方法可有效地抑制MT資料中近源效應(yīng)的影響.本文目前的研究?jī)H僅針對(duì)中低頻段的MT信號(hào)，對(duì)高頻MT信號(hào)需要采用不同的參數(shù)，將在下一步研究中探討.同時(shí)，采用多分辨率分析法，對(duì)于0.01 Hz以下的信號(hào)損失較大，可能造成有用信號(hào)的缺失.因此，采用該方法處理后的MT數(shù)據(jù)進(jìn)行視電阻率計(jì)算及后期資料解釋時(shí)，0.01 Hz以下的數(shù)據(jù)不參與計(jì)算.

致謝 感謝吉林大學(xué)地探學(xué)院范翠松對(duì)本研究提供的技術(shù)幫助.

Fan C S, Li T L, Wang D Y. 2008. Treatment of wavelet transform for square wave noise in MT data.JournalofJilinUniversity(EarthScienceEdition) (in Chinese), 38(S1): 61-63.

Fan C S. 2009. The strong noise characteristics of MT in ore concentration area and research of Denoise method (in Chinese)[Master′s thesis]. Changchun: Jilin University.

Gamble T D, Goubau W M, Clarke J. 1979. Error analysis for remote reference magnetotellurics.Geophysics, 44(5): 959-968.

Goubau W M, Gamble T D, Clarke J. 1978. Magnetotelluric data analysis: Removal of bias.Geophysics, 43(6): 1157-1166.

Hu C H, Li G H, Zhou T. 2008. System Analysis and Design Based on MATLAB 7.x-Wavelet Analysis (in Chinese). 3nd ed. Xi′an: Xidian University Publisher.

Huang N E, Shen Z, Long S R, et al. 1998. The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis.ProceedingsoftheRoyalSocietyA:Mathematical,PhysicalandEngineeringSciences, 454(1971): 903-995. Jansen M, Bultheel A. 1999. Multiple wavelet threshold estimation by generalized cross validation for images with correlated noise.IEEETransactionsonImageProcessing, 8(7): 947-953.

Johnstone I M, Silverman B W. 2005. Ebayes thresh: R programs for empirical Bayes thresholding.JournalofStatisticalSoftware, 12(8): 1-38. Ke X Z, Wang L, Ni G R. 2008. Application of improved threshold Denoising based on wavelet transform to pulsar signal processing.JournalofXi′anUniversityofTechnology(in Chinese), 24(1): 18-21.Li S Y, Lin J, Yang G H, et al. 2013. Ground-Airborne electromagnetic signals de-noising using a combined wavelet transform algorithm.ChineseJ.Geophys. (in Chinese), 56(9): 3145-3152, doi: 10.6038/cjg20130927.

Li T L, Fan C S, Wang D Y. 2009. The characteristic of strong EM noises of magnetotelluric data in ore deposit concentration area and their suppression techniques (in Chinese).∥Proceedings of the 25th Annual Conference of the Chinese Geophysical Society. Hefei: Chinese Geophysical Society, 259.

Liu W, Cao S Y, Wang Z, et al. 2013. The self-adaptive random noise attenuation in curvelet domain based on Bayes estimation.GeophysicalProspectingforPetroleum(in Chinese), 52(2): 115-120.

Tang J T, Hua X R, Cao Z M, et al. 2008. Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data.ChineseJ.Geophys. (in Chinese), 51(2): 603-610, doi: 10.3321/j.issn:0001-5733.2008.02.034.

Tang J T, Li J, Xiao X, et al. 2012. Magnetotelluric sounding data strong interference separation method based on mathematical morphology filtering.JournalofCentralSouthUniversity(ScienceandTechnology) (in Chinese), 43(6): 2215-2221.

Tang J T, Li J, Xiao X, et al. 2015. Mathematical morphology filtering and noise suppression of magnetotelluric sounding data.ChineseJ.Geophys. (in Chinese), 55(5): 1784-1793, doi: 10.6038/j.issn.0001-5733.2012.05.036.

Wang G, Wang S M, Zhu W, et al. 2015. Method of suppressing the adjacent source interference in AMT.GeophysicalandGeochemicalExploration(in Chinese), 39(2): 371-375. Weckmann U, Magunia A, Ritter O. 2005. Effective noise separation for magnetotelluric single site data processing using a frequency domain selection scheme.GeophysicalJournalInternational, 161(3): 635-652.

Xu Z M, Tang J T, Qiang J K. 2012. An analysis of the magnetotelluric strong interference types in ore concentration areas.GeophysicalandGeochemicalExploration(in Chinese), 36(2): 214-219.

Yan L J, Hu W B, Chen Q L, et al. 1998. Application of remote reference MT to noisy area.JournalofJianghanPetroleumInstitute(in Chinese), 20(4): 34-38.

Yang S, Bao G S, Zhang Q S. 2002. A study on the application of remote reference magnetotelluric sounding technique.GeophysicalandGeochemicalExploration(in Chinese), 26(1): 27-31, 49.Zhang C. 2013. Magnetotelluric data quality assessment and impedance estimation (in Chinese)[Master′s thesis]. Changsha: Central South University.Zhou X X, Wang X M, Yang Y, et al. 2014. De-noising of high-speed turnout vibration signals based on wavelet threshold.JournalofVibrationandShock(in Chinese), 33(23): 200-206.

Zhu W, Fan C S, Yan D W, et al. 2011. Noise source analysis and noise characteristics study of MT in an ore concentration area.GeophysicalandGeochemicalExploration(in Chinese), 35(5): 658-662.

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

范翠松, 李桐林, 王大勇. 2008. 小波變換對(duì)MT數(shù)據(jù)中方波噪聲的處理. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 38(S1): 61-63.

范翠松. 2009. 礦集區(qū)強(qiáng)干擾大地電磁噪聲特點(diǎn)及去噪方法研究[碩士論文]. 長春: 吉林大學(xué).

胡昌華, 李國華, 周濤. 2008. 基于MATLAB 7.x的系統(tǒng)分析與設(shè)計(jì)——小波分析. 3版. 西安: 西安電子科技大學(xué)出版社.

柯熙政, 汪麗, 倪廣仁. 2008. 改進(jìn)的小波閾值消噪法應(yīng)用于脈沖星弱信號(hào)處理. 西安理工大學(xué)學(xué)報(bào), 24(1): 18-21.

李肅義, 林君, 陽貴紅等. 2013. 電性源時(shí)域地空電磁數(shù)據(jù)小波去噪方法研究. 地球物理學(xué)報(bào), 56(9): 3145-3152, doi: 10.6038/cjg20130927.

李桐林, 范翠松, 王大勇. 2009. 礦集區(qū)大地電磁強(qiáng)干擾噪聲的特點(diǎn)與壓制技術(shù)研究. ∥中國地球物理學(xué)會(huì)第二十五屆年會(huì)論文集. 合肥: 中國地球物理學(xué)會(huì), 259.

劉偉, 曹思遠(yuǎn), 王征等. 2013. 基于貝葉斯閾值估計(jì)的曲波域自適應(yīng)隨機(jī)噪聲衰減. 石油物探, 52(2): 115-120.

湯井田, 化希瑞, 曹哲民等. 2008. Hilbert-Huang變換與大地電磁噪聲壓制. 地球物理學(xué)報(bào), 51(2): 603-610, doi: 10.3321/j.issn:0001-5733.2008.02.034.

湯井田, 李晉, 肖曉等. 2012. 基于數(shù)學(xué)形態(tài)濾波的大地電磁強(qiáng)干擾分離方法. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 43(6): 2215-2221.

湯井田, 李晉, 肖曉等. 2015. 數(shù)學(xué)形態(tài)濾波與大地電磁噪聲壓制. 地球物理學(xué)報(bào), 55(5): 1784-1793, doi: 10.6038/j.issn.0001-5733.2012.05.036.

王剛, 王書民, 朱威等. 2015. AMT近源干擾壓制方法. 物探與化探, 39(2): 371-375.

徐志敏, 湯井田, 強(qiáng)建科. 2012. 礦集區(qū)大地電磁強(qiáng)干擾類型分析. 物探與化探, 36(2): 214-219.

嚴(yán)良俊, 胡文寶, 陳清禮等. 1998. 遠(yuǎn)參考MT方法及其在南方強(qiáng)干擾地區(qū)的應(yīng)用. 江漢石油學(xué)院學(xué)報(bào), 20(4): 34-38.

楊生, 鮑光淑, 張全勝. 2002. 遠(yuǎn)參考大地電磁測(cè)深法應(yīng)用研究. 物探與化探, 26(1): 27-31, 49.

張弛. 2013. 大地電磁數(shù)據(jù)質(zhì)量評(píng)價(jià)與阻抗估計(jì)[碩士論文]. 長沙: 中南大學(xué).

周祥鑫, 王小敏, 楊揚(yáng)等. 2014. 基于小波閾值的高速道岔振動(dòng)信號(hào)降噪. 振動(dòng)與沖擊, 33(23): 200-206.

朱威, 范翠松, 姚大為等. 2011. 礦集區(qū)大地電磁噪聲場(chǎng)源分析及噪聲特點(diǎn). 物探與化探, 35(5): 658-662.

(本文編輯 胡素芳)

A combined wavelet transform algorithm used for de-noising magnetotellurics data in the strong human noise

LING Zhen-Bao1, WANG Pei-Yuan1, WAN Yun-Xia1*, WANG Yan-Zhang1, CHENG De-Fu1，LI Tong-Lin2

1TheCollegeofInstrumentationandElectricalEngineering,JilinUniversity，Changchun130021，China2TheCollegeofGeoexplorationScienceandTechnology,JilinUniversity，Changchun130021，China

Magnetotellurics (MT) is normally applied for mineral exploration, because of its low cost and deeper penetration, as well as high resolving power in horizontally. However, the quality of MT data is affected by the human noise seriously around the mine area, which result in a unreliable result of inversion and even a wrong explanation. Therefore, it is important to remove the various disturbance. The multi-resolution algorithm (MRA) has domination in terms of making a high resolution for the frequency domain of MT data, and the wavelet thresholding algorithm (WTA) is at an advantage when removing the high frequency noise. A combined wavelet transform algorithm based on MRA and WTA is proposed. db3 wavelet is adopted in processing. Baseline drift and periodic noise of square wave can be removed by multi-resolution algorithm. Bayes estimation combined with improved threshold function which based on algorithm of wavelet threshold value is used to eliminate the impulse noise, triangular waveform as well as other forms of high frequency noises. Processing case indicated that the quality of the processed data, including the time series data and apparent resistivity curves, is improved significantly. Meanwhile near-source effect, which is caused by human interference, is suppressed effectively.

Human disturbance; Multi-resolution analysis; Wavelet threshold; Apparent resistivity; Source effect

Daubechies I. 1990. The wavelet transform, time-frequency localization and signal analysis.IEEETransactionsonInformationTheory, 36(5): 961-1005. Donoho D L. 1995. De-noising by soft-thresholding.IEEETransactionsonInformationTheory, 41(3): 613-627.Egbert G D, Booker J R. 1986. Robust estimation of geomagnetic transfer functions.GeophysicsJournalInternational, 87(1): 173-194.

10.6038/cjg20160926.

國家自然科學(xué)基金項(xiàng)目(41404094)資助.

凌振寶，男，1966年生，教授，1998年碩士畢業(yè)于長春地質(zhì)學(xué)院電磁測(cè)量技術(shù)及儀器專業(yè)，主要從事方向傳感器技術(shù)應(yīng)用.

E-mail:lingzhenbao@jlu.edu.cn

10.6038/cjg20160926

P631

2016-02-01，2016-06-06收修定稿

凌振寶， 王沛元， 萬云霞等. 2016. 強(qiáng)人文干擾環(huán)境的電磁數(shù)據(jù)小波去噪方法研究. 地球物理學(xué)報(bào)，59(9):3436-3447,

Ling Z B, Wang P Y, Wan Y X, et al. 2016. A combined wavelet transform algorithm used for de-noising magnetotellurics data in the strong human noise.ChineseJ.Geophys. (in Chinese),59(9):3436-3447,doi:10.6038/cjg20160926.

*通訊作者 萬云霞，女，1980年生，工程師，博士，2013年畢業(yè)于吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院.主要研究領(lǐng)域是頻率域電磁探測(cè)和數(shù)字信號(hào)處理.E-mail:wanyx@jlu.edu.cn