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

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

0 引言

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

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

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

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

1 算法原理

1.1 奇異值分解

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

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

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的奇異值.

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

（2）

式中1

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

1.2 多分辨率奇異值分解

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

1.2.1 三分遞推SVD算法

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

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

（3）

（2）對(duì)矩陣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）分別為對(duì)應(yīng)的左正交矩陣和右正交矩陣，相應(yīng)的對(duì)角矩陣為Σ1=（diag（α1,1,α1,2,α1,3）,0）（Σ1∈R3×3）.其中，α1,1、α1,2、α1,3分別為信號(hào)第1次分解得到的奇異值.

（3）通過(guò)式（4）得到三個(gè)奇異值α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對(duì)應(yīng)的奇異值大，反映信號(hào)的主體成分，為近似矩陣；HD1和Hd1對(duì)應(yīng)的奇異值小，反映信號(hào)的細(xì)節(jié)成分，為細(xì)節(jié)矩陣；從這三個(gè)矩陣可以分別恢復(fù)出近似信號(hào)A1和細(xì)節(jié)信號(hào)D1、d1.

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

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

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

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

1.2.3 三分遞推SVD的分解性能

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

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

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

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

（7）

2 模擬仿真實(shí)驗(yàn)

2.1 大地電磁信號(hào)分布特性

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

圖2 原始信號(hào)的時(shí)域波形和頻譜Fig.2 Time domain waveform and frequency spectrum of original signal

（續(xù)圖3）

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

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

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

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

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

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

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

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

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

圖5 測(cè)試樣本奇異值分布特性（a） 無(wú)干擾；（b） 類(lèi)充放電三角波；（c） 類(lèi)脈沖；（d） 類(lèi)方波.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 測(cè)點(diǎn)QH401504中數(shù)據(jù)段的奇異值分布特性Fig.6 Singular value distribution characteristics of data segments in site QH401504

圖7 經(jīng)MRSVD分解1層得到的近似信號(hào)和細(xì)節(jié)信號(hào)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）分別為近似信號(hào)標(biāo)準(zhǔn)差和細(xì)節(jié)信號(hào)標(biāo)準(zhǔn)差.

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

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

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

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

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

圖9 不同干擾的辨識(shí)效果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，對(duì)大地電磁數(shù)據(jù)均勻分段；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖18 測(cè)點(diǎn)BL22200J的視電阻率-相位曲線對(duì)比圖Fig.18 Comparison of apparent resistivity-phase curve for site BL22200J

圖19 測(cè)點(diǎn)EL22189A的視電阻率-相位曲線對(duì)比圖Fig.19 Comparison of apparent resistivity-phase curve for site EL22189A

圖20 測(cè)點(diǎn)EL22195A的視電阻率-相位曲線對(duì)比圖Fig.20 Comparison of apparent resistivity-phase curve for site EL22195A

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

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

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

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

3.1 時(shí)間序列分析

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

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

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

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

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

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

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

4 結(jié)論

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

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