局域均值分解和小波閾值在大地電磁噪聲壓制中的應用

李 晉，彭 沖，湯井田，燕 歡，蔡劍華

(1.湖南師范大學 物理與信息科學學院, 長沙 410081;2.中南大學 地球科學與信息物理學院, 有色金屬成礦預測與地質環境監測教育部重點實驗室, 長沙 410083;3.湖南文理學院 物理與電子科學學院, 湖南 常德 415000)

局域均值分解和小波閾值在大地電磁噪聲壓制中的應用

李 晉1,2，彭 沖1，湯井田2，燕 歡1，蔡劍華3

(1.湖南師范大學 物理與信息科學學院, 長沙 410081;2.中南大學 地球科學與信息物理學院, 有色金屬成礦預測與地質環境監測教育部重點實驗室, 長沙 410083;3.湖南文理學院 物理與電子科學學院, 湖南 常德 415000)

大地電磁測深法是基于電磁感應原理，利用天然交變電磁場來研究地下巖層的電學性質及其分布特征。然而，天然電磁場頻帶范圍寬、信號微弱，在實際測量中大地電磁信號極易受到各種電磁噪聲干擾，嚴重影響了后續的電磁法反演解釋水平。針對這一難題，將局域均值分解(LMD)的自適應性和小波分析的多分辨性相結合，提出基于局域均值分解和小波閾值的大地電磁噪聲壓制方法。將含噪信號進行LMD分解得到若干階乘積函數(PF)分量；根據大地電磁信噪特征保留PF1分量，僅對其余各階PF分量選取合適的小波閾值進行降噪處理；疊加重構獲得大地電磁有用信號。通過計算機模擬典型強干擾，研究不同小波函數、分解層數及閾值方式下算法的去噪性能，并將其應用于礦集區實測大地電磁數據處理。實驗結果表明，所提方法能較好地提取出疊加在微弱大地電磁信號上的大尺度強干擾的輪廓特征，視電阻率曲線更為光滑、連續，低頻段的大地電磁數據質量得到了明顯改善。

局域均值分解； 小波閾值； 大地電磁； 噪聲壓制

大地電磁測深法(Magnetotelluric，MT)誕生于20世紀50年代，是一種以天然交變電磁場為場源，通過測量地表相互正交的電場和磁場，獲得地下電性結構信息的地球物理方法。與有源的電磁勘探方法相比，天然大地電磁場頻帶范圍寬且本身信號極其微弱，野外觀測到的大地電磁信號不可避免地會受到各種噪聲的污染[1]。尤其是在礦集區，隨處可見的高壓電網、廣播電臺、通訊電纜、信號發射塔、各種金屬管網以及用于礦山開采的大功率直流電機車等嚴重影響了實測大地電磁信號的采集，大地電磁測深數據質量極具下降，阻抗估算偏差嚴重及測量獲得的視電阻率值過度失真等狀況導致不能客觀反映地下電性分布[2]。為此，如何消除大地電磁信號中的噪聲干擾、提高大地電磁測深數據質量是國內外長期矚目并不斷取得進展的研究課題。諸多現代信號處理方法，如Hilbert-Huang變換[3]、廣義S變換[4]、方差比維納濾波[5]、同步時間序列依賴[6]、數學形態濾波[7]、子空間增強[8]等均被應用到該領域，并在一定程度上對大地電磁測深數據質量起到了積極的改善作用，推動了MT法在與各種噪聲競爭中不斷成長。

局域均值分解(Local Mean Decomposition, LMD)是一種分析非線性、非平穩信號的時頻分析方法，能自適應地將待處理信號分解為若干階乘積函數(Production Functiion, PF)的線性組合[9]；該方法已成功應用于機械故障盲源分離、滾動軸承故障診斷等領域[10-12]。小波分析法具有良好的時頻分析和多分辨率分析能力，能較好地分析非線性、非平穩信號。本文根據大地電磁信號的特征，借助于機械振動中廣泛使用的LMD法和小波分析法，將LMD的自適應性和小波分析的多分辨性相結合，提出一種基于局域均值分解和小波閾值的大地電磁噪聲壓制方法。通過對模擬信號和礦集區實測大地電磁數據進行仿真實驗，結果表明該方法可用于壓制礦集區典型大尺度方波和充放電三角波干擾，大地電磁測深數據質量得到了明顯改善。

1 局域均值分解

(1)

(2)

(3)

(4)

h1(t)=x(t)-m1(t)

h2(t)=s1(t)-m2(t)

?

hn(t)=sn-1(t)-mn(t)

(5)

式中:

s1(t)=h1(t)/a1(t)

s2(t)=h2(t)/a2(t)

?

sn(t)=hn(t)/an(t)

(6)

其中:迭代終止的條件為

(7)

(6) 將迭代過程中產生的所有局域包絡估計函數ai(t)相乘，得到PF分量的包絡信號，即瞬時幅值函數a1(t)：

(8)

(7) 將包絡信號a1(t)和純調頻信號sn(t)相乘，得到原始信號的第一個PF分量：

PF1(t)=a1(t)sn(t)

(9)

式中：PF1(t)包含了原始信號x(t)中最高的頻率成分，且是一個單分量的調幅-調頻信號;瞬時頻率f1(t)可由純調頻信號sn(t)求出，即：

(10)

(8) 將原始信號x(t)減去第一個PF分量PF1(t)，得到一個新的余量信號r1(t)，并將r1(t)作為新的輸入數據重復上述步驟、循環k次，直到rk(t)為單調函數：

r1(t)=x(t)-PF1(t)

r2(t)=r1(t)-PF2(t)

?

rk(t)=rk-1(t)-PFk(t)

(11)

經過上述步驟，原始信號x(t)被分解為k個PF分量與rk(t)之和：

(12)

圖1所示為仿真信號x(t)的時域波形及LMD分解效果圖，x(t)的表達式如下：

(13)

圖1 x(t)的LMD分解效果Fig.1 LMD decomposition effect of x(t)

2 小波閾值選取及閾值去噪

小波變換的實質是信號在小波基空間進行分解，小波基空間則由選定的小波母函數經過伸縮和平移構成。

(14)

將ψ(t)平移和伸縮可得到一組小波序列，對于連續情況小波序列表示為

(15)

式中:a為伸縮因子;b為平移因子。

(16)

其逆變換表示為

(17)

小波閾值的選取方式主要有選擇啟發式閾值(Heursure)、極小極大方差閾值(Minimaxi)、基于Stein無偏似然估計閾值(Rigrsure)和固定閾值(Sqtwolog)四種方式，其閾值函數可分為軟閾值函數和硬閾值函數兩種情況。

小波閾值去噪的基本步驟如下：

(1) 選擇與噪聲形狀相似的小波函數和合適的分解層數，將待處理信號進行小波變換得到分解后的小波系數；

(2) 利用不同的閾值選取方式確定閾值，選取合適的閾值函數(軟、硬閾值)對各層小波系數進行量化處理；

(3) 對新得到的小波系數進行重構，獲得去噪后的信號估計。

3 仿真信號分析

觀測礦集區海量實測大地電磁數據可知，時間域波形中通常出現大尺度方波和充放電三角波干擾[16-17]。為此，計算機模擬典型的大尺度方波和充放電三角波干擾進行仿真實驗。降噪性能采用曲線相似度NCC來衡量，定義如下：

(18)

圖2和圖3分別所示為含不同寬度和幅值的大尺度方波和充放電三角波模擬信號經LMD分解后得到的各階PF分量效果圖。

圖2 大尺度方波模擬信號LMD分解圖

Fig.2 LMD decomposition of large-scale square simulation signal

圖3 大尺度充放電三角波模擬信號LMD分解圖

Fig.3 LMD decomposition of large-scale charging and discharging triangular simulation signal

從圖2和圖3可知，局域均值分解將大尺度方波和充放電三角波模擬信號自適應地分解為7階PF分量，每階PF分量都與一定的物理過程相對應，將各階PF分量的瞬時頻率和瞬時幅值組合則可得到原始信號完整的時頻分布。根據大地電磁信噪特征分析大尺度方波和充放電三角波模擬信號的LMD分解結果可以發現，分解得到的PF1分量通常由大量類似于大地電磁信號的隨機信息構成。為此，基于局域均值分解和小波閾值的大地電磁噪聲壓制方法其具體步驟如下：首先，利用局域均值分解的自適應性將大地電磁信號分解為若干階PF分量；然后，保留PF1分量，僅對其余各階PF分量選擇合適的閾值進行小波閾值去噪；最后，將處理后的小波系數疊加重構獲得大地電磁有用信號。由于硬閾值函數在閾值處不連續，易導致重構信號出現偽吉布斯現象，因此文中均采用軟閾值函數進行處理，目的是使其盡可能多地保留有用信號，減少重構信號與真實信號之間的偏差。

圖4所示為大尺度方波和充放電三角波模擬信號采用本文所提方法的噪聲壓制效果圖。

圖4 大尺度模擬信號去噪效果圖
Fig.4 Denoising effect of large-scale simulation signal

從圖4可知，文中所提方法在壓制大尺度方波和充放電三角波的同時，原始信號中豐富的細節成分得到了有效保留。

為進一步分析算法性能，表1所示為模擬信號經LMD分解后，再利用不同小波函數、分解層數及閾值方式獲得的曲線相似度對比情況。

表1 曲線相似度對比Tab.1 The contrastion of curve similarity

分析表1可知，對于大尺度方波模擬信號在Heursure 和Rigrsure方式下采用sym6小波函數進行5層分解得到的曲線相似度最高(NCC=0.995 4)；對于大尺度充放電三角波模擬信號在Sqtwolog方式下采用sym6小波函數進行5層分解得到的曲線相似度最高(NCC=0.975 8)。上述仿真實驗為有針對性地選取濾波參數處理實測大地電磁強干擾類型提供了一定的依據。

4 實際資料處理

4.1 時間域去噪效果

為了驗證本文方法的實用性，對礦集區采集到的受強干擾嚴重的大地電磁數據進行噪聲壓制。鑒于大地電磁數據量龐大且噪聲類型復雜多樣，文中給出一段某測點在同一時刻采集的含強干擾的電、磁場分量去噪效果圖，如圖5所示。

從圖5可知，由于礦集區人煙稠密、重工業密集，導致電道和磁道中采集的大地電磁數據均不同程度地受到了類似于周期性突跳和波動等大尺度方波及充放電三角波干擾，且電道Ex/Ey和磁道Hy/Hx出現干擾的時刻具有一定的相關性，微弱的大地電磁有用信號幾乎被完全湮沒。經本文所提方法處理后，大尺度強干擾的輪廓特征被較好地提取，且曲線光滑、連續；重構的大地電磁信號則在剔除基線漂移的同時，保留了有用信號更多的細節信息，原始大地電磁信號的基本特征得到了較好的還原。

圖5 實測MT信號去噪效果圖
Fig.5 Denoising effect of measured MT signal

4.2 實測點處理效果

文中選用的實測大地電磁數據由加拿大鳳凰公司生產的V5-2000大地電磁測深儀采集，該儀器的數據采集方式是1-8-5模式。以TSL、TSH文件為例，低頻段的數據為全時間段采集(采樣率為24 Hz)，數據存儲為TSL文件；高頻和中頻是交替采集，每5 min采集1次高頻或中頻，其中有1 s的高頻數據(采樣率2 560 Hz)和連續8 s的中頻數據(采樣率為320 Hz)，采樣起止時間相同，數據存儲為TSH文件。因此，大部分反映測點信息的數據集中在低頻段，對礦集區大地電磁數據進行噪聲壓制主要是針對低頻段的典型強干擾異常波形進行處理。為此，文中重點對TSL文件中的大地電磁數據進行降噪處理。分析V5-2000數據存儲格式可知，參與計算的視電阻率值其頻率成分最多至0.8 Hz左右。

圖6和圖7分別所示為廬樅礦集區實測點B40492和B44249原始數據和本文所提方法在時間域波形及視電阻率曲線兩方面的去噪效果對比圖。兩測點數據采集時間均在20 h以上，其中時間域波形的觀測環境為V5 System 2000。

分析圖6(a)和圖7(a)可知，由于兩測點面臨復雜的噪聲干擾環境，原始數據的電道(Ex、Ey)和磁道(Hx、Hy)中均出現大尺度方波和充放電三角波干擾，且能量幅值遠大于大地電磁有用信號幾個數量級，正常的大地電磁信號幾乎被完全湮沒。圖6(b)和圖7(b)為經文中方法處理后，將大地電磁數據重新還原至V5 System 2000觀測的時間域波形。對比分析圖6(a)和圖7(a)可知，大尺度方波和充放電三角波干擾均得到了有效剔除，微弱的大地電磁有用信號基本集中在基線附近且無明顯跳變。

觀測實測點B40492原始數據的視電阻率曲線(圖6(c))可知，大約在50～0.05 Hz頻段視電阻率曲線呈近45°漸近線快速上升，且曲線光滑、幾乎沒有誤差棒，該頻段內近13個頻點的數據完全失真；在0.05 Hz左右，視電阻率值達到最大，接近于1 000 000 Ω·m；0.05～0.005 Hz頻段的視電阻率值突然下降至1 000 Ω·m左右，視電阻率曲線的誤差棒增大，且曲線跳變劇烈，該現象屬于典型的近源效應。因此，由原始測點數據得到的視電阻率曲線儼然已不能客觀反映測點本身所固有的地下介質電性結構。分析圖6(d)可知，經本文所提方法處理后，0.8～0.05 Hz頻段的視電阻率曲線近45°上升的近源效應有明顯緩解，視電阻率最大值下降了近2個數量級，該頻段內挽救了近9個頻點的數據；0.05～0.000 5 Hz超低頻段，視電阻率曲線更為光滑、連續，視電阻率值相對穩定，且誤差棒減小。分析實測點B44249原始數據的視電阻率曲線(圖7(c))可知，曲線形態同樣具有典型的近源效應，且在低頻段曲線震蕩劇烈、誤差棒增大。經文中所提方法(圖7(d))處理后，視電阻率曲線近45°上升的近源效應得到了有效壓制，視電阻率曲線形態更為光滑、連續，誤差棒減小、視電阻率值相對穩定，低頻段的大地電磁數據質量得到了明顯改善。

(a) 原始數據時間域波形

(b) 本文所提方法時間域波形

(c) 原始數據視電阻率曲線

(d) 本文所提方法視電阻率曲線

(a) 原始數據時間域波形

(b) 本文所提方法時間域波形

(c) 原始數據視電阻率曲線

(d) 本文所提方法視電阻率曲線

5 結 論

(1) 提出了一種基于局域均值分解和小波閾值的大地電磁噪聲壓制方法。將局域均值分解的自適應性和小波分析的多分辨性相結合，應用于礦集區大地電磁微弱信號和大尺度強干擾的信噪分離，達到了噪聲壓制目的，大尺度異常波形和基線漂移均得到了有效壓制。

(2) 通過模擬信號仿真分析和實測資料數據處理，證明了方法的實用性。受近源干擾嚴重的測點經文中所提方法處理后，視電阻率曲線的整體形態更為光滑、連續，近源干擾得到了有效壓制，其結果更加真實地反映了測點本身所固有的大地電磁深部構造信息。方法為今后在礦集區開展大地電磁測深提供了一條新的解決途徑，應用前景廣闊。

[1] CHEN J, HEINCKE B, JEGEN M, et al.Using empirical mode decomposition to process marine magnetotelluric data[J].Geophysical Journal International, 2012, 190(1): 293-309.

[2] 湯井田, 徐志敏, 肖曉, 等.廬樅礦集區大地電磁測深強噪聲的影響規律[J].地球物理學報, 2012, 55(12): 4147-4159.

TANG Jingtian, XU Zhimin, XIAO Xiao, et al.Effect rules of strong noise on magnetotelluric (MT) sounding in the Luzong ore cluster area[J].Chinese J.Geophys., 2012, 55(12): 4147-4159.

[3] CAI J H, TANG J T, HUA X R, et al.An analysis method for magnetotelluric data based on the Hilbert-Huang transform[J].Exploration Geophysics, 2009, 40(2): 197-205.

[4] 景建恩, 魏文博, 陳海燕, 等.基于廣義S變換的大地電磁測深數據處理[J].地球物理學報, 2012, 55 (12): 4015-4022.

JING Jianen, WEI Wenbo, CHEN Haiyan, et al.Magnetotelluric sounding data processing based on generalized S transformation[J].Chinese J.Geophys., 2012, 55(12): 4015-4022.

[5] KAPPLE K N.A data variance technique for automated despiking of magnetotelluric data with a remote reference[J].Geophysical Prospecting, 2012, 60(1): 179-191.

[6] 王輝, 魏文博, 金勝, 等.基于同步大地電磁時間序列依賴關系的噪聲處理[J].地球物理學報, 2014, 57(2): 531-545.

WANG Hui, WEI Weibo, JIN Sheng, et al.Removal of magnetotelluric noise based on synchronous time series relationship[J].Chinese J.Geophys., 2014, 57(2): 531-545.

[7] 湯井田, 李晉, 肖曉, 等.數學形態濾波與大地電磁噪聲壓制[J].地球物理學報, 2012,55(5): 1784-1793.

TANG Jingtian, LI Jin, XIAO Xiao, et al.Mathematical morphology filtering an noise suppression of magnetotelluric sounding data[J].Chinese J.Geophys., 2012, 55(5): 1784-1793.

[8] 李晉, 湯井田, 王玲, 等.基于信號子空間增強和端點檢測的大地電磁噪聲壓制[J].物理學報, 2014, 63(1): 019101.

LI Jin, TANG Jingtian, WANG Ling, et al.Noise suppression for magnetotelluric sounding data based on signal subspace enhancement and endpoint detection[J].Acta Phys.Sin., 2014, 63(1): 019101.

[9] 程軍圣, 楊宇, 于德介.局部均值分解方法及其在齒輪故障診斷中的應用[J].振動工程學報, 2009, 22(1): 76-84.

CHENG Junsheng, YANG Yu, YU Dejie.The local mean decomposition method and its application to gear fault diagnosis[J].Journal of Vibration Engineering, 2009, 22(1): 76-84.

[10] 李志農, 劉衛兵, 易小兵.基于局域均值分解的機械故障欠定盲源分離方法研究[J].機械工程學報, 2011, 47(7): 97-102.

LI Zhinong, LIU Weibing, YI Xiaobing.Underdetermined blind source separation method of machine faults based on local mean decomposition[J].Journal of Mechanical Engineering, 2011, 47(7): 97-102.

[11] 武哲, 楊紹普, 張建超.基于LMD自適應多尺度形態學和Teager能量算子方法在軸承故障診斷中的應用[J].振動與沖擊, 2016, 35(3): 7-13.

WU Zhe, YANG Shaopu, ZHANG Jianchao.Bearing fault feature extraction method based on LMD adaptive multiscale morphology and energy operator demodulating[J].Journal of Vibration and Shock, 2016, 35(3): 7-13.

[12] 張焱, 湯寶平, 鄧蕾, 等.基于局域均值分解的自適應濾波滾動軸承故障特征提取[J].振動與沖擊, 2015, 34(23): 25-30.

ZHANG Yan, TANG Baoping, DENG Lei, et al.Fault feature extraction for rolling bearing based on adaptive wavelet filtering and LMD[J].Journal of Vibration and Shock, 2015, 34(23): 25-30.

[13] SMITH J S.The local mean decomposition and its application to EEG perception data[J].Journal of the Royal Society Interface, 2005, 2(5): 443-454.

[14] 侯高燕, 呂勇, 肖涵, 等.基于LMD的多尺度形態學在齒輪故障診斷中的應用[J].振動與沖擊, 2014, 33(19): 69-73.

HOU Gaoyan, Lü Yong, XIAO Han, et al.Based on the LMD and multi-scale morphology used in gear fault diagnosis[J].Journal of Vibration and Shock, 2014, 33(19): 69-73.

[15] 張亢, 程軍圣, 楊宇.基于局部均值分解與形態譜的旋轉機械故障診斷方法[J].振動與沖擊, 2013, 32(9): 135-140.

ZHANG Kang, CHENG Junsheng, YANG Yu.Rotating machinery fault diagnosis based on local mean decomposition and pattern spectrum[J].Journal of Vibration and Shock, 2013, 32(9): 135-140.

[16] 湯井田, 張弛, 肖曉, 等.大地電磁阻抗估計方法對比[J].中國有色金屬學報, 2013, 23(9): 2351-2358.

TANG Jingtian, ZHANG Chi, XIAO Xiao, et al.Comparison of methods for magnetotelluric impedance estimation[J].The Chinese Journal of Nonferrous Metals, 2013, 23(9): 2351-2358.

[17] 李晉, 湯井田, 肖曉, 等.基于組合廣義形態濾波的大地電磁資料處理[J].中南大學學報(自然科學版), 2014, 45(1): 173-185.

LI Jin, TANG Jingtian, XIAO Xiao, et al.Magnetotelluric data processing based on combined generalized morphological filter[J].Journal of Central South University (Science and Technology), 2014, 45(1): 173-185.

Application of local mean decomposition and wavelet threshold in magnetotelluric noise suppression

LI Jin1,2, PENG Chong1, TANG Jingtian2, YAN Huan1, CAI Jianhua3

(1.Institute of Physics and Information Science, Hunan Normal University, Changsha 410081, China;2.School of Geosciences and Info-Physics, Key Laboratory of Metallogenic Prediction of Non-Ferrous Metals and Geological Environment Monitor, Ministry of Education, Central South University, Changsha 410083, China;3.Department of Physics and Electronics, Hunan University of Arts and Science, Changde 415000, China)

The magnetotelluric sounding method is a method based on the principle of electromagnetic induction, to detect electrical properties and distribution characteristics of underground rock layers by use of the natural alternating electromagnetic field.However, the natural electromagnetic field has a wider frequency band range and weak signals, and these signals are easy to be disturbed by all kinds of electromagnetic noise in the actual measurement, so the interpretation level of the subsequent electromagnetic inversion method is seriously effected.In order to solve this problem, combining the adaptability of local mean decomposition (LMD) with the multi-resolution of wavelet, here the magnetotelluric noise suppression method based on local mean decomposition and wavelet threshold was proposed.First of all, LMD was used to divide a noisy signal into a number of product functiion (PF) components.Then, according to magnetotelluric signal-to-noise characters, PF1component was retained, the appropriate wavelet threshold was chosen to denoise the resteach PF component.Finally, the magnetotelluric useful signal was obtained with superposition and reconstruction.Using a computer to simulate typical strong interferences, the denoising performances of the proposed method were studied under the conditions of different wavelet functions, decomposition layers and threshold modes, and the method was applied to process the magnetotelluric data measured in ore concentration areas.The results showed that the proposed method can better extract the outline features of large-scale strong interferences superimposed on weak magnetotelluric signals, and the apparent resistivity curve is more smooth and continuous; the quality of magnetotelluric data is improved significantly in lower frequency bands.

local mean decomposition; wavelet threshold; magnetotelluric; noise suppression

國家自然科學基金(41404111; 41304098)；國家863計劃(2014AA06A602)；湖南省自然科學基金(2015JJ3088);中國博士后科學基金(2015M570687)

2016-04-14 修改稿收到日期：2016-08-05

李晉 男，博士后，副教授，1981年生

彭沖 男，碩士生，1989年生

TN911.4； P631

A

10.13465/j.cnki.jvs.2017.05.021