基于自適應(yīng)多分辨率奇異值分解的大地電磁數(shù)據(jù)處理

李晉，馬翻紅，湯井田，李勇

1 湖南師范大學(xué)信息科學(xué)與工程學(xué)院，長沙 410081 2 中南大學(xué)有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室，長沙 410083 3 中國地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所，自然資源部地球物理電磁法探測技術(shù)重點實驗室，河北廊坊 065000

0 引言

20世紀(jì)50年代，Tikhonov和Cagniard提出了大地電磁測深法（Magnetotelluric,MT）（Tikhonov,1950;Cagniard,1953）；該方法是一種探測地殼深部結(jié)構(gòu)的地球物理方法（趙國澤等,2004），具有探測深度大、成本低、施工方便、垂向分辨能力和水平分辨能力高等優(yōu)點，在地震預(yù)測、油氣田勘探與普查、礦產(chǎn)資源勘查等領(lǐng)域已得到廣泛應(yīng)用（湯吉等,2005;魏文博等,2009;金勝等,2010;胡祥云等,2012;董樹文等,2012）.然而，頻帶范圍極寬的天然大地電磁信號非常微弱，且不可避免地會受到各類噪聲的干擾，尤其是在礦集區(qū)，大地電磁數(shù)據(jù)質(zhì)量急劇下降，導(dǎo)致視電阻率-相位曲線發(fā)生畸變、阻抗估計過度失真（白登海等,2002）.此時的數(shù)據(jù)已不能客觀反映地下電性結(jié)構(gòu)，嚴(yán)重影響了后續(xù)電磁反演成像的可靠性和可解釋性，給地球物理勘探帶來極大困擾（湯井田等,2012b;呂慶田等,2015）.因此，有效壓制強干擾，提升數(shù)據(jù)質(zhì)量是大地電磁測深領(lǐng)域最具挑戰(zhàn)的任務(wù)之一，對精細(xì)勘探地下深部結(jié)構(gòu)、勘查深部礦產(chǎn)資源具有重要的理論意義和實際應(yīng)用價值.

國內(nèi)外學(xué)者針對不同的噪聲提出了許多有針對性的去噪方法.Gamble等（1979）提出了遠(yuǎn)參考大地電磁測深法，該方法用于去除同源相關(guān)電磁噪聲；Egbert和Booker（1986）提出了消除非高斯分布噪聲的Robust統(tǒng)計方法；Kao和Rankin（1977）提出了利用最小二乘法中的循環(huán)計算來消除非相關(guān)噪聲；由于小波變換的多分辨特性，Trad和Travassos（2000）將其用于非平穩(wěn)的大地電磁數(shù)據(jù)處理，隨后該方法在阻抗估計和長（短）周期噪聲壓制等方面均得到廣泛應(yīng)用（徐義賢和王家映,2000;范翠松等,2008;凌振寶等,2016;Ling et al.,2019）；王書明和王家映（2004）將高階統(tǒng)計量運用于大地電磁數(shù)據(jù)處理，有效抑制了高斯有色噪聲；為了壓制工頻干擾，湯井田等（2008）和Cai（2017）將Hilbert-Huang變換引入大地電磁數(shù)據(jù)處理；緊接著，互補集成經(jīng)驗?zāi)B(tài)分解等系列方法在大地電磁數(shù)據(jù)處理中得以應(yīng)用，該方法解決了集成經(jīng)驗?zāi)B(tài)分解（史恒等,2011）出現(xiàn)的分解誤差，以及經(jīng)驗?zāi)B(tài)分解（Empirical Mode Decomposition,EMD）出現(xiàn)的模態(tài)混疊和端點效應(yīng)（Li et al.,2019）；湯井田等利用數(shù)學(xué)形態(tài)濾波有效抑制了具有明顯形態(tài)特征的大尺度干擾（湯井田等,2012a），并結(jié)合稀疏表示和字典學(xué)習(xí)對大地電磁數(shù)據(jù)進(jìn)行去噪（湯井田等,2018）；李晉等研究了匹配追蹤和變分模態(tài)分解（Variational Mode Decomposition,VMD）的大地電磁信噪分離方法（李晉等,2018,2019;Zhang et al.,2021）.

分析國內(nèi)外相關(guān)文獻(xiàn)可知，大地電磁噪聲強度大且環(huán)境復(fù)雜，尤其是在強干擾區(qū)，現(xiàn)有方法具有一定的優(yōu)勢，同時也存在一定的局限性.

近年來，奇異值分解（Singular Value Decomposition,SVD）在信號處理領(lǐng)域迅猛發(fā)展.與小波分析相比，該方法具有穩(wěn)健性強、信噪比高、零相移、波形失真小等優(yōu)點，在故障診斷、電能質(zhì)量檢測、信號識別和消噪，以及數(shù)據(jù)挖掘等諸多領(lǐng)域得到廣泛應(yīng)用.Phillips等（2009）將SVD應(yīng)用于衛(wèi)星遙感數(shù)據(jù)降維，獲得了比主成分分析（王文波等,2013）更精確的特征，分類速度更快且更準(zhǔn)確；Lehtola等（2008）將數(shù)學(xué)形態(tài)學(xué)和SVD相結(jié)合，對心電信號進(jìn)行消噪處理，提高了心電信號的準(zhǔn)確性.趙學(xué)智等（2010,2017）在故障診斷領(lǐng)域，結(jié)合SVD對矩陣分解及相空間矩陣構(gòu)造等進(jìn)行探索，并將小波的多分辨思想應(yīng)用于SVD，提出了多分辨率奇異值分解（Multi-Resolution Singular Value Decomposition,MRSVD），同時分析了SVD和小波變換在信號處理方面的優(yōu)勢.Luo和Zhang（2019）將MRSVD應(yīng)用于滾動軸承周期脈沖特征提取，實現(xiàn)了初期的故障檢測.He等（2020）采用MRSVD抑制振動噪聲，準(zhǔn)確地檢測了振動信號中的故障點，并利用長短時記憶神經(jīng)網(wǎng)絡(luò)預(yù)測了軸承性能；Zhang等（2020）研究了各種多分結(jié)構(gòu)的奇異值去噪模型.鑒于MRSVD在去噪方面的優(yōu)越性，本文結(jié)合MRSVD和標(biāo)準(zhǔn)差差值提出基于自適應(yīng)多分辨率奇異值分解（Adaptive Multi-Resolution Singular Value Decomposition,AMRSVD）的大地電磁數(shù)據(jù)處理方法，并與EMD、VMD等方法進(jìn)行分析對比.結(jié)果表明，該方法不受噪聲先驗信息的影響、處理效率高，且能有效分離相關(guān)性較好的噪聲，保留了大地電磁數(shù)據(jù)的原始特征.

1 算法原理

1.1 奇異值分解

Kaleman（1996）提出奇異值分解，該方法原理描述如下：

對任意一個矩陣H∈Rm×n，都有正交矩陣U=（u1,u2,…,um）∈Rm×m及正交矩陣V=（v1,v2,…,vn）∈Rn×n，使得矩陣H滿足

H=UΣVT，

（1）

式中，Σ=（diag（α1,α2,…,αr）,0）（Σ∈Rm×n），r=min（m,n），其中α1≥α2≥…≥αr>0，αi（i=1,2,…,r）為矩陣H的奇異值.

對信號進(jìn)行SVD處理時，首先將信號構(gòu)造Hankel矩陣H.設(shè)待處理信號x（t）=[x1,x2,x3,…,xN]，矩陣H的形式如下：

（2）

式中1

一般地，只對信號構(gòu)造一次Hankel矩陣，然后運用SVD進(jìn)行處理.從本質(zhì)上講，這種對信號采用單一的SVD處理，其結(jié)果是利用正交矩陣U（V）的行（列）矢量得到的，都屬于同層次的矢量空間.

1.2 多分辨率奇異值分解

為了將SVD拓展到不同層次的矢量空間，趙學(xué)智等（2010）借鑒小波的多分辨率特征，在SVD的基礎(chǔ)上提出一種二分遞推矩陣構(gòu)造方法;Zhang等（2020）在二分遞推矩陣的基礎(chǔ)上提出多分遞推矩陣構(gòu)造方法，用于將復(fù)雜信號分解到不同層次的子空間.

1.2.1 三分遞推SVD算法

三分遞推SVD算法步驟如下：

（1）對待處理數(shù)據(jù)x（t）用三分遞推構(gòu)造一個Hankel矩陣H1：

（3）

（2）對矩陣H1進(jìn)行SVD處理：

（4）

式中，U1=|_u1,1,u1,2,u1,3_|（U1∈R3×3）、V1=|_v1,1,v1,2,v1,3_|（V1∈R（N-2）×3）分別為對應(yīng)的左正交矩陣和右正交矩陣，相應(yīng)的對角矩陣為Σ1=（diag（α1,1,α1,2,α1,3）,0）（Σ1∈R3×3）.其中，α1,1、α1,2、α1,3分別為信號第1次分解得到的奇異值.

（3）通過式（4）得到三個奇異值α1,1、α1,2、α1,3，并進(jìn)行SVD重構(gòu)得到：

（5）

式中，u1,1、u1,2和u1,3為3×1的矩陣，分別表示u1,1=（u1,1,1,u1,1,2,u1,1,3）T、u1,2=（u1,2,1,u1,2,2,u1,2,3）T及u1,3=（u1,3,1,u1,3,2,u1,3,3）T，將其代入式（5）得到：

（6）

同理，得到HD1和Hd1.HA1對應(yīng)的奇異值大，反映信號的主體成分，為近似矩陣；HD1和Hd1對應(yīng)的奇異值小，反映信號的細(xì)節(jié)成分，為細(xì)節(jié)矩陣；從這三個矩陣可以分別恢復(fù)出近似信號A1和細(xì)節(jié)信號D1、d1.

（4）通過上述步驟得到第1次分解的結(jié)果，再利用近似信號A1構(gòu)造Hankel矩陣H2返回步驟（2）進(jìn)行同樣的處理，可將原始信號分解為不同分辨率的近似信號和細(xì)節(jié)信號，如圖1所示.

圖1 三分遞推SVD分解流程Fig.1 The decomposition process of three recursive SVD

1.2.2 三分遞推SVD的多分辨率特性

多分辨率指的是L2（R）的子空間Vj中每個子空間Vj+1都是它前一級Vj的精細(xì)化.從式（6）可知，在三分遞推SVD中，近似信號Aj是利用vj,1得到的，而細(xì)節(jié)信號Dj和dj分別是由vj,2和vj,3得到；vj,1、vj,2及vj,3稱為近似基矢量和細(xì)節(jié)基矢量，通過這些基矢量可得到不同分辨率的近似信號和細(xì)節(jié)信號.從三分遞推SVD的分解過程可知，下一層分解的基矢量是由上一層的近似基矢量得到，說明三分遞推SVD在分解過程中實現(xiàn)了類似于小波的多分辨率分解，即多分辨率奇異值分解（MRSVD）.

1.2.3 三分遞推SVD的分解性能

為了測試MRSVD的優(yōu)越性，構(gòu)造如圖2所示的原始信號做模擬仿真實驗，其中f1=sin（80πt），f2=sin（300πt）×（1+2sin（30πt）），f3=f1+f2；左邊為時域波形圖，右邊為其對應(yīng)的頻譜.

將f3（合成信號）作為分解對象，采用EMD、VMD、SVD和MRSVD四種方法進(jìn)行處理，如圖3所示；其中左邊為時域波形圖，右邊為其對應(yīng)的頻譜.分析圖3可知，（a）中EMD分解得到的IMF1和IMF2出現(xiàn)端點效應(yīng)；（b）中VMD雖克服了EMD分解存在的缺陷，但單個分量的頻率成分不連續(xù)；（c）中SVD將原始信號分成了三個分量，每個分量雖都包含原信號信息，但仍無法得出精確的分解效果；（d）中MRSVD將原始信號分解為A1～A7七個近似信號，最后一個分量為細(xì)節(jié)信號的疊加；原始信號經(jīng)MRSVD分解得到的近似信號A6的時間域波形與f1的相似度達(dá)到0.9921，A7與f1的相似度為0.9927，且（d）中最后一個分量（f3-A7）與f2的時間域波形相似度達(dá)到0.9866，說明信號分解到A7時對應(yīng)的細(xì)節(jié)信號的頻譜即f3-A7與f2最為接近.分析不同方法的分解效果可知，MRSVD能將原始信號中的復(fù)雜信號以細(xì)節(jié)信號的形式分解出來，凸顯了隱藏的原始信號特征.

1.3 標(biāo)準(zhǔn)差

標(biāo)準(zhǔn)差在概率統(tǒng)計中通常用來作為統(tǒng)計分布程度上的測量，反映數(shù)據(jù)集個體間的離散程度.標(biāo)準(zhǔn)差越大，表示數(shù)據(jù)和平均值之間的差異較大；反之，則表示越接近平均值.

（7）

2 模擬仿真實驗

2.1 大地電磁信號分布特性

大地電磁場指的是由地球外部場源引起天然電磁場短周期變化形成的電場和磁場，其頻率范圍極

圖2 原始信號的時域波形和頻譜Fig.2 Time domain waveform and frequency spectrum of original signal

（續(xù)圖3）

圖3 不同分解方法的時域波形和頻譜對比圖（a） EMD;（b） VMD;（c） SVD;（d） MRSVD.Fig.3 Comparison of time domain waveform and frequency spectrum of different decomposition methods

寬、信號極其微弱，且場源極化方向隨機；同時具有復(fù)雜的時變特性，野外采集的大地電磁信號極易受到各種電磁噪聲的污染.由于天然大地電磁信號是典型的非線性、非平穩(wěn)信號，實測數(shù)據(jù)受到的大尺度噪聲和低頻有用信號在統(tǒng)計分布上具有明顯區(qū)別.

圖4所示為一組測試樣本，分別為青海地區(qū)某測點中一段幾乎未受電磁干擾的大地電磁數(shù)據(jù)段和礦集區(qū)測點中含類充放電三角波干擾、類脈沖干擾、類方波干擾的數(shù)據(jù)段.

圖5所示為圖4中的測試樣本經(jīng)式（3）和（4）得到的三個奇異值α1,1、α1,2、α1,3的分布特性，其中每200個采樣點為1段，共分為15段.分析圖5可知，無干擾信號的三個奇異值在數(shù)值上沒有明顯差異；類充放電三角波干擾α1,1明顯大于α1,2和α1,3，且α1,2和α1,3非常接近；類方波干擾和類脈沖干擾在有干擾的部分，α1,1也明顯大于α1,2和α1,3；同樣地，圖6所示為青海點QH401504中一段電道Ex和磁道Hx數(shù)據(jù)經(jīng)式（3）和（4）得到的三個奇異值α1,1、α1,2、α1,3的分布特性，其中左邊為時域波形圖，右邊為其對應(yīng)的奇異值分布特性.從圖6可知，電道和磁道數(shù)據(jù)在受干擾的情況下同圖5相比，其奇異值分布具有相似性.分析圖5和圖6可知，奇異值大的可表征干擾數(shù)據(jù)，奇異值小的則可代表低頻信號成分.由于α1,2和α1,3分別為Dj和dj的奇異值，且相差不大，為此后續(xù)將細(xì)節(jié)信號定義為Dj和dj的疊加.

圖7所示為圖4測試樣本中的無干擾信號、類充放電三角波干擾、類脈沖干擾和類方波干擾數(shù)據(jù)分別采用MRSVD分解1層得到的近似信號和細(xì)節(jié)信號.

分析圖7可知，圖4中的無干擾信號經(jīng)MRSVD得到的細(xì)節(jié)信號和近似信號沒有明顯區(qū)別；類充放電三角波干擾分解得到的近似信號主要表現(xiàn)為強干擾，細(xì)節(jié)信號則主要由微弱的有用信號構(gòu)成；類脈沖和類方波干擾分解之后的細(xì)節(jié)信號中出現(xiàn)大量的脈沖干擾，分解效果不佳.

為了區(qū)分強干擾和微弱的大地電磁有用信號，分別計算圖7中的近似信號標(biāo)準(zhǔn)差與細(xì)節(jié)信號標(biāo)準(zhǔn)差的差值（簡稱標(biāo)準(zhǔn)差差值），如圖8所示.

分析圖8可知，經(jīng)MRSVD分解1層后，無干擾信號的標(biāo)準(zhǔn)差差值在基線附近，類充放電三角波干擾基本都在0.5以上；類脈沖干擾在采樣點為500時有明顯干擾，但差值仍在基線附近，干擾數(shù)據(jù)區(qū)分不明顯；類方波干擾在采樣點0～500處有明顯干擾，其標(biāo)準(zhǔn)差差值也在0.5以上.為此，文中采用MRSVD分解1層得到的標(biāo)準(zhǔn)差差值對類充放電三角波干擾和類諧波干擾進(jìn)行信噪辨識，其辨識參數(shù)定義為：

圖4 測試樣本數(shù)據(jù)Fig.4 Test sample data

圖5 測試樣本奇異值分布特性（a） 無干擾；（b） 類充放電三角波；（c） 類脈沖；（d） 類方波.Fig.5 Singular value distribution characteristics of test sample data（a） No interference;（b） Charge-discharge-like triangle wave;（c） Pulse-like wave;（d） Square-like wave.

圖6 測點QH401504中數(shù)據(jù)段的奇異值分布特性Fig.6 Singular value distribution characteristics of data segments in site QH401504

圖7 經(jīng)MRSVD分解1層得到的近似信號和細(xì)節(jié)信號Fig.7 The approximate signal and detail signal obtained by MRSVD decomposition layer 1

|Δψ|=|ψA1-ψ（D1+d1）|，

（8）

約束條件定義為|Δψ|<θ（θ為0.4～0.9），其中ψA1、ψ（D1+d1）分別為近似信號標(biāo)準(zhǔn)差和細(xì)節(jié)信號標(biāo)準(zhǔn)差.

圖9所示為兩段實測數(shù)據(jù)（含類充放電三角波干擾和類諧波干擾）采用上述方法和約束條件得到的信噪辨識效果.分析圖9可知，有用信號段和強干擾數(shù)據(jù)段可以有效區(qū)分.

圖8 近似信號標(biāo)準(zhǔn)差和細(xì)節(jié)信號標(biāo)準(zhǔn)差的差值對比Fig.8 Comparison of difference between approximate signal standard deviation and detail signal standard deviation

圖10所示為測試樣本（圖4）中含類充放電三角波干擾的原始數(shù)據(jù)經(jīng)MRSVD分解11層得到的近似信號.分析圖10可知，類充放電三角波干擾的噪聲輪廓隨著分解層數(shù)的增大逐漸分解到近似信號部分.分析A1～A8可知，噪聲輪廓逐漸光滑，低頻信號被分解到細(xì)節(jié)信號部分，A8～A11的噪聲輪廓則沒有發(fā)生明顯變化，說明一定的分解層數(shù)對MRSVD的去噪效果敏感.

圖11所示為圖9中兩段實測數(shù)據(jù)經(jīng)MRSVD分解50層得到的相鄰細(xì)節(jié)信號標(biāo)準(zhǔn)差差值的變化情況.

分析圖11可知，兩段含強干擾數(shù)據(jù)的相鄰細(xì)節(jié)信號標(biāo)準(zhǔn)差差值隨分解層數(shù)的增加，差值增大隨后逐漸減小直到0.001～0.01趨于穩(wěn)定狀態(tài).結(jié)合圖10分析可知，噪聲輪廓隨著分解層數(shù)的增加逐步變得光滑，但當(dāng)分解到一定層數(shù)后，近似信號不再變化，對應(yīng)的細(xì)節(jié)信號也不再變化.考慮到大地電磁信號的主要特征集中在細(xì)節(jié)信號（微弱的低頻成分），為此文中將相鄰細(xì)節(jié)信號標(biāo)準(zhǔn)差差值|Δψ（Dj+dj）|作為最佳分解層數(shù)的約束條件（劉嫣和湯偉,2016），即若存在分解層數(shù)j使得：

圖9 不同干擾的辨識效果Fig.9 Identification effect of different interference

|Δψ（Dj+dj）|<ω，

（9）

則MRSVD停止分解，即完成自適應(yīng)多分辨率奇異值分解（AMRSVD），式中ω為0.001～0.01.

2.2 算法流程

基于AMRSVD的大地電磁數(shù)據(jù)處理算法流程如圖12所示.

具體步驟如下：

（1）初始化參數(shù)：分解層數(shù)j=1、ω取0.001～0.01、θ取0.4～0.9，對大地電磁數(shù)據(jù)均勻分段；

（2）對大地電磁數(shù)據(jù)構(gòu)建Hankel矩陣H，對其進(jìn)行SVD分解得到近似信號和細(xì)節(jié)信號；

（3）計算近似信號和細(xì)節(jié)信號的標(biāo)準(zhǔn)差差值|Δψ|，判斷是否滿足|Δψ|<θ，若滿足則認(rèn)為是大地電磁有用信號，否則進(jìn)入步驟（4）進(jìn)行信噪分離處理；

圖10 原始信號經(jīng)MRSVD分解得到的近似信號Fig.10 The approximate signal obtained by MRSVD decomposition of the original signal

圖11 相鄰細(xì)節(jié)信號標(biāo)準(zhǔn)差差值變化Fig.11 Variation of standard deviations of adjacent detail signal

（4）對分解出的近似信號重建Hankel矩陣且j=j+1，利用SVD分解得到下一層的近似信號和細(xì)節(jié)信號；

（5）計算相鄰細(xì)節(jié)信號的標(biāo)準(zhǔn)差差值|Δψ（Dj+dj）|，判斷是否滿足|Δψ（Dj+dj）|<ω，若滿足則停止分解，得到最終的近似信號即噪聲輪廓，然后用原始數(shù)據(jù)減去近似信號得到低頻信號，并與步驟（3）中的有用信號進(jìn)行拼接、重構(gòu)，否則返回步驟（4）.

2.3 EMTF數(shù)據(jù)模擬

將EMTF開源代碼包中的純凈電道（Ex、Ey）數(shù)據(jù)和磁道（Hx、Hy）數(shù)據(jù)分別加入相關(guān)性較強的大尺度噪聲進(jìn)行分析，引入均方誤差（MSE）、信噪比（SNR）、歸一化互相關(guān)性（NCC）和運行效率（Runtime）進(jìn)行評判（Li et al,2020a,2020b,2021）.

圖13和圖14所示分別為EMTF含噪數(shù)據(jù)（Ex、Ey）采用EMD、VMD、SVD和AMRSVD的時頻域去噪效果，其中左邊為時域波形圖，右邊為其對應(yīng)的頻譜.分析可知，經(jīng)EMD處理后，時域波形中大部分有用信號丟失嚴(yán)重、頻譜失真；經(jīng)VMD和SVD處理后殘留了大量的尖脈沖，去噪性能較低；AMRSVD可以較好地實現(xiàn)大尺度強噪聲和低頻信號的分離.

圖15所示為采用不同方法處理圖13和圖14的含噪EMTF數(shù)據(jù)得到的視電阻率曲線對比圖.分析圖15可知，原始EMTF數(shù)據(jù)的視電阻率曲線呈平滑直線狀態(tài)，當(dāng)Ex和Ey分別加入大尺度噪聲后視電阻率曲線發(fā)生了顯著跳變.經(jīng)EMD和SVD處理后，視電阻率曲線呈下降趨勢，低頻有用信號丟失；VMD雖然比EMD和SVD有較大提升，但中低頻段處理不佳；相比而言，AMRSVD得到的視電阻率曲線改善明顯.

圖12 AMRSVD數(shù)據(jù)處理流程圖Fig.12 Flow chart of AMRSVD data processing

圖13 Ex道數(shù)據(jù)經(jīng)不同方法處理的時頻域去噪效果Fig.13 Denoising effect of different methods in time-frequency domain of Ex data

圖15 EMTF數(shù)據(jù)的視電阻率曲線對比圖Fig.15 Comparison of apparent resistivity curve for EMTF data

圖16 實測數(shù)據(jù)時間域去噪效果圖（a） 類充放電三角波；（b） 類諧波.Fig.16 Time domain denoising effect of measured data（a） Charge-discharge-like triangle wave;（b） Harmonic-like wave.

圖17 四道實測數(shù)據(jù)段去噪前后對比（a） 去噪前；（b） 去噪后.Fig.17 Comparison of four measured data segments before and after denoising（a） Before denoising；（b） After denoising.

圖18 測點BL22200J的視電阻率-相位曲線對比圖Fig.18 Comparison of apparent resistivity-phase curve for site BL22200J

圖19 測點EL22189A的視電阻率-相位曲線對比圖Fig.19 Comparison of apparent resistivity-phase curve for site EL22189A

圖20 測點EL22195A的視電阻率-相位曲線對比圖Fig.20 Comparison of apparent resistivity-phase curve for site EL22195A

表1和表2分別為Ex道和Ey道數(shù)據(jù)經(jīng)不同方法的去噪性能對比.分析表1、表2可知，AMRSVD在NCC、MSE、SNR均優(yōu)于EMD、VMD和SVD的去噪效果，尤其是AMRSVD的去噪速度更快、效率更高.

表1 Ex道數(shù)據(jù)不同方法的去噪性能對比Table 1 Comparison of denoising performance of different methods of Ex data

表2 Ey道數(shù)據(jù)不同方法的去噪性能對比Table 2 Comparison of denoising performance of different methods of Ey data

3 實測數(shù)據(jù)處理

3.1 時間序列分析

為了驗證方法的實用性，將其應(yīng)用于廬樅礦集區(qū)相關(guān)性較強的大尺度類充放電三角波干擾和類諧波干擾的大地電磁時間序列，如圖16所示.

分析圖16a可知，AMRSVD可以有效去除大地電磁數(shù)據(jù)中的類充放電三角波干擾、保留更多的有用信號.分析圖16b可知，AMRSVD對一些類諧波干擾也可以提取出光滑的噪聲輪廓曲線.

圖17所示為廬樅礦集區(qū)某測點中同一時段的電道和磁道數(shù)據(jù)經(jīng)AMRSVD法去噪前后的時域波形對比圖.

分析圖17a可知，Ex/Hy和Ey/Hx表現(xiàn)出很強的噪聲相關(guān)性；分析圖17b可知，四道數(shù)據(jù)中的強干擾均被有效去除，低頻信號得到了更好保留.

3.2 視電阻率-相位曲線分析

圖18—20所示分別為廬樅礦集區(qū)測點BL22200J、EL22189A和EL22195A經(jīng)遠(yuǎn)參考（RR）、EMD、VMD、SVD和AMRSVD處理的視電阻率-相位曲線對比圖.

分析圖18—20可知，3個測點的原始數(shù)據(jù)視電阻率曲線在低頻段急劇上升，30～0.5 Hz的視電阻率曲線基本呈45°，視電阻率值從102Ωm上升到了106Ωm，提升了4個數(shù)量級，對應(yīng)的相位在0°或-180°附近且分布紊亂，這些測點的數(shù)據(jù)真實性急劇下降，表現(xiàn)為典型的近源效應(yīng).經(jīng)RR處理后，視電阻率曲線在30～5 Hz變得連續(xù)、光滑，但在5～0.3 Hz處的視電阻率曲線跳躍大且分布凌亂，表明遠(yuǎn)參考對近源干擾無能為力.經(jīng)EMD處理后，視電阻率曲線在30～5 Hz較為連續(xù)，但小于5 Hz的曲線部分頻點的值急劇下降且相位較為紊亂，低頻信息丟失嚴(yán)重.經(jīng)VMD處理后，30～5 Hz視電阻曲線仍呈45°上升，3～0.3 Hz視電阻率曲線呈下降趨勢，去噪效果不佳.SVD由于分解矩陣難以確定，導(dǎo)致SVD的去噪效果較差，視電阻率曲線紊亂.分析廬樅礦集區(qū)大量原始大地電磁數(shù)據(jù)可知，電道中含有大量的類充放電三角波噪聲，磁道中含有不同程度的類充放電三角波噪聲和類諧波噪聲，且相關(guān)性強，這些噪聲可能是引起近源效應(yīng)的主要原因，經(jīng)AMRSVD處理后視電阻率曲線45°上升的趨勢得到明顯緩解.然而由于該礦集區(qū)采集的原始大地電磁數(shù)據(jù)中廣泛分布脈沖干擾，導(dǎo)致1～10 Hz左右的頻段數(shù)據(jù)嚴(yán)重受損，而文中方法對脈沖噪聲分離效果不佳，重構(gòu)的信號中殘留有尖脈沖干擾.

4 結(jié)論

針對如何從強干擾中有效提取微弱的大地電磁信號，文中以類充放電三角波干擾和類諧波干擾作為研究對象，提出一種基于自適應(yīng)多分辨率奇異值分解的大地電磁數(shù)據(jù)處理方法.研究了大地電磁數(shù)據(jù)的奇異值分布特性，發(fā)現(xiàn)奇異值的大小能分別表征近似信號和細(xì)節(jié)信號；通過將大地電磁數(shù)據(jù)分解成不同分辨率的近似信號和細(xì)節(jié)信號，并引入標(biāo)準(zhǔn)差差值對大地電磁數(shù)據(jù)進(jìn)行信噪辨識；同時，結(jié)合相鄰標(biāo)準(zhǔn)差差值和MRSVD，提出利用AMRSVD對強干擾數(shù)據(jù)進(jìn)行去噪處理.仿真實驗和實測數(shù)據(jù)處理結(jié)果表明，本文方法能有效去除時間序列中形狀規(guī)則、相關(guān)性較強的干擾，視電阻率-相位曲線得到改善，去噪后的信號可靠性提升；與EMD、VMD和SVD等方法相比，AMRSVD的處理效果更優(yōu)、分解效率更高，為強干擾下開展大地電磁數(shù)據(jù)處理提供了一種新的研究思路.

由于本文方法的處理效果與噪聲的相關(guān)性及出現(xiàn)的頻率等有關(guān)，若噪聲的能量、幅值等特征不明顯，辨識度將降低、信噪分離會有一定的偏差.同時，該方法對尖脈沖干擾處理效果不佳.為此，接下來將重點研究方法對不同噪聲類型的穩(wěn)健性，以及引入智能算法優(yōu)化MRSVD的自適應(yīng)分解過程，提升方法的抗噪性能.