子空間增強在大地電磁信噪分離中的應用研究

燕 歡， 李 晉,2， 李春陽

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

子空間增強在大地電磁信噪分離中的應用研究

燕 歡1， 李 晉1,2， 李春陽1

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

為了提高強干擾地區大地電磁測深數據質量、保留大地電磁信號低頻段的有用信息，壓制大尺度強干擾已刻不容緩。針對礦集區典型的強干擾類型，引入子空間增強方法對大地電磁信號與強干擾進行信噪分離研究。通過計算機模擬常見的大尺度方波和充放電三角波干擾，討論了子空間增強算法中幀長與特征值判別閾值的最優選取范圍，給出了算法流程，并對礦集區實測大地電磁數據進行了信噪分離處理。實驗結果表明：該方法是切實可行的，能有效剔除時間域序列中的大尺度強干擾，近源效應得到了有效壓制；經處理后，重構的大地電磁信號中包含了更為豐富的細節成分，視電阻率曲線更為光滑、連續，低頻段的大地電磁數據質量得到了明顯改善。

大地電磁； 信噪分離； 子空間增強； 強干擾

0 引言

大地電磁測深(Magnetotelluric，簡稱MT)法是由Tikhonov和Cagniard[1-2]提出的以天然交變電磁場為場源的電磁勘探方法。由于天然電磁場信號弱、頻帶寬，因此極易受到各種噪聲的干擾，是典型的非線性、非平穩信號[3-5]。隨著人類文明的不斷進步和經濟的高速發展，人文電磁干擾日益嚴重，尤其在礦集區，電話線路、無線通訊基站、采礦用的大功率直流電機車等引起強能量的電磁干擾，導致原本微弱的大地電磁信號幾乎完全被湮沒在這些“大尺度”干擾中。諸多現代信號處理方法(如Hilbert-Huang變換[6]、廣義S變換[7]、方差比維納濾波[8]、數學形態濾波[9]、同步時間序列[10]等),被應用到大地電磁噪聲壓制，均在一定程度上改善了大地電磁測深數據質量，推動了大地電磁測深法在各種噪聲中的發展。

子空間增強是一種非線性、非平穩信號分析方法，上世紀90年代該方法已被應用于語音增強領域。由于語音信號的噪聲方差可以估計，根據子空間增強算法原理，帶噪語音信號的向量空間可分解為信號子空間和噪聲子空間，通過剔除噪聲子空間，再利用信號子空間即可恢復較為純凈的語音信號。子空間增強法現已被逐步推廣到故障診斷、地震資料處理等領域[11-12]。由于大地電磁信號是典型的非線性、非平穩信號，為此借助語音增強、地震資料處理等領域廣泛應用的子空間增強法，提出分離礦集區微弱大地電磁信號和典型強干擾的研究思路，分析了關鍵參數的最優選取方案，給出了具體的算法處理流程。

1 子空間增強基本原理

子空間分解的降噪方法早期由Dendrinos[13]提出，該算法利用帶噪信號組成具有Toeplitz結構的矩陣，隨后對該矩陣進行奇異值分解；根據最小二乘估計，忽略較小的奇異值后可將降維得到的數據根據Toeplitz矩陣結構恢復出原信號。然而，子空間增強法具有一定的局限性，只能處理白噪聲、不能處理有色噪聲。在此基礎上，Ephraim等[14]做了進一步完善，利用信號協方差矩陣的特征值分解(Eigen value decomposition，EVD)，將帶噪信號的向量空間分解為噪聲子空間和信號子空間兩個相互正交的子空間再進行語音增強處理。隨后，Jensen等[15]將基于SVD(Singular value decomposition，SVD)的子空間算法擴展到有色噪聲領域。

子空間增強具有控制信號失真和噪聲殘留的特征，其基本原理是將受干擾的帶噪信號映射到兩個子空間①信號子空間；②噪聲子空間，通過數學方法將噪聲子空間全部置零，濾除信號子空間中所包含的噪聲干擾，再利用信號子空間恢復出較為純凈的原始信號[16]。

以一維離散信號為例，假定噪聲是加性的，且與純凈信號相互獨立，則帶噪信號可表示為：

y=x+n

(1)

式中：y、x和n分別表示帶噪信號、純凈信號和加性噪聲的向量表示,對其進行分幀，可得到相應的K維向量，其帶噪信號、純凈信號和噪聲的協方差矩陣Ry、Rx和Rn關系如下：

Ry=E[yyT]=E[xxT]+E[nnT]=

Rx+Rn

(2)

式中:E[·]表示均值運算。

假設純凈信號估計為：

(3)

式中:H表示K×K維的線性估計器矩陣，則該估計器的真實值與估計值的誤差ε為：

εx+εn

(4)

式中：εx和εn分別表示信號的失真向量和殘余噪聲向量。子空間增強即通過線性估計投影到信號子空間，從而得到原始信號的較好估計。

tr(HRxHT-HRx-RxHT+Rx)

(5)

tr(HRnHT)

(6)

通過求解以下時域約束條件方程，可獲得優化的線性估計器：

(7)

式中:σ2為正常數。在可接受的噪聲殘差水平σ2下，該估計器能最小化信號失真。

上述方程的解為：

Hopt=Rx(Rx+μRn)-1

(8)

式中:μ為拉格朗日乘子。

由Rx的特征值分解為：

Rx=UΛxUT

(9)

式中：U表示Rx的歸一化特征向量矩陣；Λx表示由Rx的特征值組成的對角矩陣：

(10)

Hopt=UΛx(Λx+μUTRnU)-1UT

(11)

當噪聲n為方差σn的加性白噪聲時，Rn可表示為：

Rn=σnI

(12)

當噪聲n為有色噪聲時，往往用對角矩陣Λn來近似UTRnU：

(13)

為此，Hopt可優化為式(14)。

Hopt=UΛx(Λx+μΛn)-1UT

(14)

2 仿真信號分析

2.1 子空間增強關鍵參數分析

觀測礦集區海量大地電磁數據的時間域波形可知，大尺度方波和充放電三角波干擾通常出現在各種采樣數據中，是最為典型且嚴重的噪聲干擾類型[17]。為了驗證文中所提方法的可行性，計算機模擬大尺度方波和充放電三角波進行仿真實驗。大尺度方波通過選取Matlab自帶的block信號疊加4dB的白噪聲來模擬(圖1)。大尺度充放電三角波通過選取幅值為5%、峰谷距為10、間距為100的充放電三角波疊加隨機噪聲來模擬(圖2)。

圖1 模擬大尺度方波干擾Fig.1 Simulated large-scale square interferences

圖2 模擬大尺度充放電三角波干擾Fig.2 Simulated large-scale charging and discharging triangular interferences

由于子空間增強算法中分幀幀長和協方差矩陣特征值的判別閾值對濾波精度尤為關鍵，為此引入波形相似度對該關鍵參數進行性能評價。相似度的定義為式(15)。

(15)

式中:f(n)、g(n)分別表示兩個離散序列。NCC的值在-1～1之間，-1表示變換前后波形反向，“0”表示正交，“1”表示完全相同；NCC越接近“1”，則表示波形之間越相似。

從圖3可知，大尺度方波的幀長取值在0～8范圍內相似度逐漸上升；當幀長取值在8～18時相似度趨于平穩，且達到最大值0.9；幀長取值在18以后，相似度逐漸下降。大尺度充放電三角波的幀長取值在0～10范圍時，相似度逐漸上升；在10～18時，相似度接近于0.8；幀長為18以后，相似度逐漸下降。分析圖4可知，大尺度方波的判別閾值在10～100時，相似度在0.75～0.8之間；大尺度充放電三角波具有同樣的變化趨勢，其相似度在0.9～0.95之間。由于模擬的大尺度強干擾相比實測大地電磁信號要單一得多，而判別閾值本身即為噪聲殘留與信號失真的折中，且每一幀噪聲的特征值變化沒有實測信號復雜。綜上分析，當幀長選取10～18、判別閾值選取40～60時，濾波效果最優。

圖3 不同幀長濾波性能對比Fig.3 Filtering performance comparison of different frames

圖4 不同判別閾值濾波性能對比Fig.4 Filtering performance comparison of different criterion thresholds

2.2 算法流程

根據子空間增強基本原理，將該算法應用到模擬大尺度信噪分離，其具體步驟如下：

1)將模擬帶噪信號分段，每段為一幀，幀長在最優范圍內選取。

2)構造協方差矩陣，對其進行特征值分解。

3)假設大尺度噪聲與微弱信號互不相關，設置最優判別閾值，分離出噪聲子空間和信號子空間。

4)將噪聲子空間的特征值置零，重構原始信號。

圖5所示為模擬大尺度干擾經子空間增強后的去噪效果圖。文中幀長為12，判別閾值為40。

從圖5可知，疊加在大尺度干擾上的信號極其微弱，經文中所提方法處理后，提取的大尺度干擾輪廓較為光滑、連續，重構信號中較好地保留了更為豐富的細節成分。

圖5 模擬大尺度干擾去噪效果圖Fig.5 Denoising effect of simulated large-scale interferences(a)模擬大尺度干擾(方波);(b)子空間分離出的大尺度干擾(方波);(c)重構信號(方波);(d)模擬大尺度干擾(充放電三角波);(e)子空間分離出的大尺度干擾(充放電三角波);(f)重構信號(充放電三角波)

3 實測數據處理

3.1 時間域去噪效果

為了驗證子空間增強法的實用性，對礦集區受強干擾嚴重的大地電磁數據進行噪聲壓制。由于礦集區大地電磁數據量龐大，文中給出某測點在同一時刻電道(EX、EY)和磁道(HY、HX)中四道時間域片段的去噪效果，如圖6所示。

從圖6可知，由于該測點所處地區噪聲干擾源眾多，導致時間域波形中出現大量能量強、頻帶寬的大尺度類方波、類充放電三角波等噪聲干擾；波形的整體形態呈不連續狀態、跳躍性強，且電道EX/EY和磁道HY/HX出現干擾的時刻具有一定的相關性，微弱的大地電磁有用信號幾乎被完全湮沒。經子空間增強法處理后，大尺度干擾的輪廓特征幾乎被完整地提取出來，分離出的大地電磁微弱信號在剔除基線漂移的同時，保留了有用信號更多的細節信息。由于實測大地電磁信號所含噪聲種類眾多，即使是同種形態的噪聲其幅值、間距也是多種多樣的。在如此復雜的干擾環境下，子空間增強算法仍能較清晰地提取出疊加在原始大地電磁信號上的大尺輪廓，說明該方法在時間域對大尺度干擾的分離是可行的。

3.2 視電阻率曲線

選用加拿大鳳凰公司生產的V5-2000大地電磁測深儀采集的數據進行分析，該儀器采樣頻率通常為高頻(2 560 Hz)、中頻(320 Hz)和低頻(24 Hz)，數據的存儲格式為TSH、TSL，其中TSH文件記錄了2 560 Hz和320 Hz采樣率的中高頻數據，TSL文件記錄了24 Hz采樣率的低頻數據[18]。由于高頻和中頻是交替采集，每5 min采集1次高頻或中頻，其中只有1 s的高頻數據和連續8 s的中頻數據，所含信息量較少；低頻數據為全時間段采集，能記錄豐富的大地電磁數據，且更能反映測點深部的電性結構信息，為此我們重點對低頻信號(TSL文件)進行去噪處理。

圖7所示為廬縱礦集區B40993測點的原始時間域波形(EX、EY、HX、HY)和視電阻率曲線。

從圖7(a)可知，原始數據電道EX、EY中出現大量大尺度類方波干擾，磁道HX、HY中出現大量大尺度類充放電三角波干擾。分析圖7(b)可知，原始數據的視電阻率曲線整體形態連續性差。在大于40 Hz時，視電阻率曲線比較平穩，且變化趨勢一致；在0.05 Hz～50 Hz，視電阻率曲線呈45°左右漸近線快速上升，0.05 Hz左右視電阻率值達到最大值接近于100 000 Ω·m，表現為典型的近源效應；在0.005 Hz～0.05 Hz，視電阻率曲線突然下降至100 Ω·m左右，誤差棒增大，并出現不同程度的突跳與畸變。分析原始數據的時間域波形和視電阻率曲線可知，該測點受到了各種復雜干擾源的影響，其數據已不能客觀反映真實的地電信息。雖然功率譜篩選具有較強的適用性，但對礦集區高強度近源干擾也無能為力。

圖6 實測大地電磁數據去噪效果Fig.6 Denoising effect of measured magnetotelluric data(a)原始實測信號EX;(b)子空間分離出的大尺度干擾EX;(c)重構大地電磁信號EX;(d)原始實測信號EY;(e)子空間分離出的大尺度干擾EY;(f)重構大地電磁信號EY;(g)原始實測信號HY;(h)子空間分離出的大尺度干擾HY;(i)重構大地電磁信號HY;(j)原始實測信號HX;(k)子空間分離出的大尺度干擾HX;(l)重構大地電磁信號HX

圖8為子空間增強法得到的時間域波形和視電阻率曲線。其中，圖8(a)為經子空間增強法處理后，重新還原至V5 System 2000觀測的時間域波形。

從圖8(a)可知，充斥在時間域中的大尺度類方波和類充放電三角波已被基本剔除、基線漂移得到了有效壓制，處理后的時間域波形基本連續、無明顯跳變。分析圖8(b)可知，經子空間增強法處理后，0.05 Hz～50 Hz呈近45°上升的近源趨勢完全消失，曲線變得光滑、連續、同步；在0.05 Hz～5 Hz，視電阻率最大值下降了近3個數量級，近15個頻點的數據質量得到了明顯改善。分析對比圖7可知，時、頻域的處理更加真實、合理、置信，其結果已逐步接近于測點本身所固有的地下介質電性結構。

圖9所示為在子空間增強法基礎上做簡單的功率譜篩選得到的視電阻率曲線。

分析圖9可知，經簡單功率譜篩選后，視電阻率曲線在0.000 5 Hz～0.05 Hz的突跳頻點得到了較好恢復，曲線下掉的趨勢抬升顯著，且誤差棒減小；視電阻率曲線的整體形態更為光滑、連續，視電阻率曲線相對穩定。實驗結果表明，經子空間增強法處理后，受近源干擾嚴重的測點其低頻段的整體數據質量得到了明顯改善。

圖7 測點B40993原始大地電磁數據Fig.7 Original MT data at site B40993(a)時間域波形;(b)視電阻率曲線

圖8 子空間增強方法Fig.8 The subspace enhancement method(a)時間域波形;(b)視電阻率曲線

圖9 經簡單功率譜篩選的視電阻率曲線Fig.9 Apparent resistivity curve through simple power spectrum filtering

4 結論

1)將子空間增強算法應用于礦集區大地電磁信噪分離，通過模擬典型的強干擾類型，分析了算法中影響濾波精度的幀長及特征值判別閾值的最優選值范圍，給出了算法流程。

2)通過模擬信號仿真和實測數據處理，驗證了方法的可行性。受近源干擾嚴重的測點經子空間增強法處理后，時間域的大尺度干擾和基線漂移得到了有效壓制，視電阻率曲線的整體形態更為光滑、連續，低頻段的大地電磁數據質量得到了明顯改善。子空間增強法對提升礦集區大地電磁測深深部勘探能力，具有重要的實際應用價值。

[1] TIKHONOV A N. On determining electrical characteristics of the deep layers of the Earth’s crust[J]. Dok1. Akad. Nauk. SSSR, 1950, 73(2): 295-297.

[2] CAGNIARD L.Basic theory of the magnetotelluric method of geophysical prospecting[J].Geophysics,1953,18(3):605-635.

[3] 朱威, 范翠松, 姚大為, 等. 礦集區大地電磁噪聲場源分析及噪聲特點[J]. 物探與化探, 2011, 35(5):658-662. ZHU W, FAN C S , YAO D W,et al. Noise source analysis and noise characteristic study of MT in an ore concentration area[J]. Geophysical and Geochemical Exploration, 2011, 35(5): 658-662. (In Chinese)

[4] 王書明, 王家映. 大地電磁信號統計特征分析[J]. 地震學報, 2004, 26(6): 669-674. WANG S M, WANG J Y. Analysis on statistic characteristics of magnetotelluric signal[J]. Acta Seismologica Sinca, 2004, 26(6): 669-674. (In Chinese)

[5] 王書明, 王家映. 大地電磁信號非最小相位性的討論[J]. 地球物理學進展, 2004, 19(2): 216-221. WANG S M, WANG J Y. Discussion on the non-minimum phase of magnetotelluric signals[J]. Progress in Geophysics, 2004, 19(2): 216-221. (In Chinese)

[6] 湯井田, 化希瑞, 曹哲民, 等. Hilbert-Huang變換與大地電磁噪聲壓制[J]. 地球物理學報, 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[J]. Chinese Journal Geophysics, 2008, 51(2):603-610. (In Chinese)

[7] 景建恩, 魏文博, 陳海燕, 等. 基于廣義S變換的大地電磁測深數據處理[J]. 地球物理學報, 2012, 55 (12): 4015-4022. JING J E, WEI W B, CHEN H Y, et al. Magnetotelluric sounding data processing based on generalized S transformation[J]. Chinese Journal Geophysics，2012, 55(12): 4015-4022. (In Chinese)

[8] 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.

[9] 湯井田, 李晉, 肖曉, 等. 數學形態濾波與大地電磁噪聲壓制[J]. 地球物理學報, 2012, 55(5): 1784-1793. TANG J T, LI J, XIAO X, et al. Mathematical morphology filtering and noise suppression of magnetotelluric sounding data[J]. Chinese Journal Geophysics, 2012, 55(5): 1784-1793. (In Chinese)

[10]王輝, 魏文博, 金勝, 等. 基于同步大地電磁時間序列依賴關系的噪聲處理[J]. 地球物理學報, 2014, 57(2): 531-545. WANG H, WEI W B, JIN S, et al. Removal of magnetotelluric noise based on synchronous time series relationship[J]. Chinese Journal Geophysics, 2014, 57(2): 531-545. (In Chinese)

[11]唐貴基，龐彬，劉尚坤. 基于奇異差分譜和平穩子空間分析的滾動軸承故障診斷[J]. 振動與沖擊, 2015, 34 (11): 83-87. TANG G J, PANG B, LIU S K. Fault diagnosis of rolling bearings based on difference spectrum of singular value and stationary subspace analysis[J]. Journal of Vibration and Shock, 2015, 34 (11): 83-87. (In Chinese)

[12]陸文凱, 丁文龍, 張善文, 等. 基于信號子空間分解的三維地震資料高分辨率處理方法[J]. 地球物理學報, 2009, 52 (8): 2174-2181. LU W K, DING W L, ZHANG S W, et al. A high resolution processing technique for 3_D seismic data based on signal subspace decomposition[J]. Chinese Journal Geophysics, 2009, 52(8): 2174-2181. (In Chinese)

[13]DENDRINOS M, BAKAMIDIS S, CARAYANNI G. Speech enhancement from noise: A regenerative approach[J]. Speech Communication, 1991, 10(2): 45-57.

[14]EPHRAIM Y, VAN TREES H L. A signal subspace approach for speech enhancement[J]. IEEE Transactions on Speech and Audio Processing, 1995, 3(4): 251-266.

[15]JENSEN S H, HANSEN P C, HANSEN S D, et al. Reduction of broad-band noise in speech by truncated QSVD[J]. IEEE Transactions on Speech and Audio Processing, 1995, 3(6): 439-448.

[16]李晉, 湯井田, 王玲, 等. 基于信號子空間增強和端點檢測的大地電磁噪聲壓制[J]. 物理學報, 2014,63(1): 019101. LI J, TANG J T, WANG L, et al. Noise suppression for magnetotelluric sounding data based on signal subspace enhancement and endpoint detection[J]. Acta Physica Sinica, 2014, 63(1): 019101. (In Chinese)

[17]湯井田, 劉子杰, 劉峰屹, 等. 音頻大地電磁法強干擾壓制試驗研究[J]. 地球物理學報, 2015, 58(12): 4636-4647. TANG J T, LIU Z J, LIU F Y, et al. The denoising of the audio magnetotelluric data set with strong interferences[J]. Chinese Journal Geophysics, 2015, 58(12): 4636-4647. (In Chinese)

[18]代小威. 基于V5-2000格式MT時間序列處理與功率譜估計及軟件開發[D]. 武漢：長江大學, 2015. DAI X W. Processing and power spectrum estimation of MT time serizes in V5-2000 format and software development[D]. Wuhan: Yangtze University, 2015. (In Chinese)

Application research of subspace enhancement for magnetotelluric signal to noise separation

YAN Huan1, LI Jin1,2, LI Chunyang1

(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)

In order to improve the quality of magnetotelluric sounding data and retain the useful information of low frequency band for magnetotelluric, suppressing the large-scale strong interference is imperative. Aimed at typical strong interference types of ore concentration area, we introduced the subspace enhancement method to study signal to noise separation of magnetotelluric data and strong interference. Through simulated large-scale square wave and charge and discharge triangle wave by computer, we discussed the optimal range of frame and criterion threshold of eigenvalue in subspace enhancement, and the concrete steps of the algorithm are presented. Then, the measured magnetotelluric data of ore concentration area were processed by signal to noise separation. Experiment results show that the method is feasible, which can effectively eliminate the large-scale strong noise interference in the time domain and the near source interference of ore concentration area has been effectively suppressed. At the same time, the apparent resistivity curve is more smoothly and continuously. Moreover, the quality of magnetotelluric sounding data is improved significantly in the low frequency band.

magnetotelluric sounding data; signal to noise separation; subspace enhancement; strong interference

2016-06-14 改回日期：2016-12-15

國家自然科學基金資助項目(41404111)；湖南省自然科學基金資助項目(2015JJ3088)；中國博士后科學基金資助項目(2015M570687) ;湖南省研究行科研創新項目(CX2017B224)

燕歡(1991-)，女，碩士，主要從事礦集區大地電磁強干擾壓制研究，E-mail：aviva163@126.com。

李晉(1981-)，男，博士后，副教授，碩士生導師，主要從事信號處理及電磁勘探研究， E-mail：geologylj@163.com。

1001-1749(2017)03-0346-08

P 631.2

A

10.3969/j.issn.1001-1749.2017.03.08