基于變分模態(tài)分解和奇異譜分析的GPR信號去噪

戴前偉，丁 浩，張 華，張 豪

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

0 引言

探地雷達(ground-penetrating radar, GPR)作為一種快速、高效的地球物理探測方法，通過接收并分析反射回波中攜帶的目標體信息來研究其形態(tài)大小、空間分布等物理特征[1-2]。因其具有操作簡便、采集快速、適應性強和無損探測等優(yōu)點，已在水文現(xiàn)象監(jiān)測[3]、城市管線測量[4]、隧道檢測[5-6]和地質災害調查[7-8]等領域廣泛應用。

在GPR的采集過程中，混入的噪聲信號具有隨機、非平穩(wěn)的特點，對此，國內(nèi)外學者已經(jīng)采取了許多去噪技術來針對性地提高GPR信號的分辨率[9]。小波閾值去噪是一種常用的表征非平穩(wěn)信號頻率隨時間變化的去噪分析方法。王勇等[10]對GPR信號的小波基進行了尋優(yōu)，戴前偉等[11]提出一種基于小波分解頻帶劃分與系數(shù)相關性相結合的閾值計算方法，鄒根等[12]對閾值選取方式和閾值函數(shù)進行了優(yōu)化，這些都取得了比常規(guī)小波閾值處理更好的效果。張夢殊等[13]在小波閾值去噪的基礎上，針對低頻小波系數(shù)進行維納濾波，聯(lián)合小波閾值和維納濾波重構得到去噪信號，更接近原始信號。然而，使用小波函數(shù)降噪過程中，小波基、分解層數(shù)和小波閾值受人為選擇影響容易使結果出現(xiàn)偏差。王憲楠等[14]發(fā)現(xiàn)在Shearlet域，GPR數(shù)據(jù)能夠得到更加稀疏的表示，再通過閾值去噪的方法將隨機噪聲去除。Liu等[15]提出了一種主成分分析-奇異值分解(PCA-SVD)混合降噪方法，通過對信號子空間進行奇異值分解和主成分分析，將子空間的雜波噪音用低秩表示并去除，從而達到去噪的效果。劉雷等[16]在時域、頻域分別使用主成分分析去噪算法對合成雷達剖面分別進行去噪分析，提高了GPR數(shù)據(jù)的信噪比。但基于主成分分析的方法在選擇前K個特征值的問題上，對去噪的效果影響很大。K越小，目標信號越能夠完整保留，抑制效果越差；K越大，抑制效果越好但會損失部分目標信號。K值的選擇目前很大程度上仍依賴于經(jīng)驗計算或者根據(jù)數(shù)值模擬的結果進行人工選擇。經(jīng)驗模態(tài)分解(empirical mode decomposition, EMD)[17]基于信號自身時間尺度特點，無需設定基函數(shù)或特征值參數(shù)，具有自適應時頻處理信號的優(yōu)勢，適合于處理非線性非平穩(wěn)信號。楊建軍等[18]論證了EMD方法在GPR信號處理上更優(yōu)于小波閾值去噪方法；馮德山等[19]使用EMD方法實現(xiàn)了在低信噪比GPR信號中突出雷達剖面中異常體特征，提高了信噪比和數(shù)據(jù)解釋精度。然而，對于復雜的信號，采用EMD方法進行信號分解時往往會出現(xiàn)模態(tài)混疊及端點飛翼的現(xiàn)象。為了解決這一問題，許軍才等[20]應用了集成經(jīng)驗模態(tài)分解(ensemble empirical mode decomposition, EEMD)對GPR信號進行去噪處理，得到了優(yōu)于EMD的去噪效果。EEMD是人為添加白噪聲，再進行平均得到各分量，解決了EMD的模態(tài)混疊問題，但在重構信號中仍存在殘留白噪聲且循環(huán)過程中增大了算法的計算量[21]。變分模態(tài)分解(variational mode decomposition, VMD)[22]是近年提出的一種時頻分析方法，在故障診斷[23]、特征提取[24]、微震信號[25]及GPR[26]數(shù)據(jù)處理中都得到了成功應用。但該方法在模態(tài)數(shù)的選擇和有效信號的提取上仍依賴于去噪效果對比和人工經(jīng)驗分析選擇，且去噪后的信號中仍含有部分中低頻噪聲，使得信號整體存在振蕩現(xiàn)象。

因此，本文提出一種自適應的VMD和奇異譜分析(singular spectrum analysis, SSA)[27]相聯(lián)合的去噪方法，進行高效準確的GPR信號去噪。通過合成Ricker子波實驗和合成雷達剖面模擬實驗，對比EEMD、VMD和本文方法，檢驗該方法的正確性和有效性。最后，將該方法應用于實際工程的GPR信號去噪處理，以去除噪聲，突出目標信號。

1 方法理論

1.1 變分模態(tài)分解

VMD是一種時頻分析方法。該算法假定所有分量都是集中在各自中心頻率附近的窄帶信號，因此，根據(jù)分量窄帶條件建立約束優(yōu)化問題，從而估計信號分量的中心頻率，能夠將多分量信號一次性分解成多個單分量調幅調頻信號，避免迭代過程中遇到的端點效應和模態(tài)混疊問題[22]。變分問題的構造步驟[28]如下。

1)采用Hilbert變換計算每個模態(tài)分量的解析信號，并獲取其對應的單邊頻譜。

2)通過在各模態(tài)中加入指數(shù)項來調整各自估計的中心頻率，將每個模態(tài)分量頻譜調整至其對應的基帶。

3)構建約束變分模型的表達式：

(1)

式中：uk為模態(tài)分量；ωk為uk的中心頻率；K為模態(tài)數(shù)；δ(t)為脈沖函數(shù)；t為時間域變量；j為虛數(shù)單位；s為原信號；*為卷積運算。

4)在變分問題的求解中，引入二次懲罰因子α和Lagrange算子λ(t)，將式(1)轉化為無約束方程：

(2)

式中，〈,〉為點積。

5)利用交替方向乘子迭代算法，結合傅里葉等距變換，優(yōu)化得到相應的模態(tài)分量和中心頻率，并搜索增廣拉格朗日表達式的鞍點。式(2)的解為：

(3)

6)重復步驟5)直到滿足迭代條件，并輸出K個分量，判別式為

(4)

式中，ε為判斷精度，且ε>0。

1.1.1 自適應模態(tài)數(shù)選擇

進行變分模態(tài)分解需要預先設置分解的模態(tài)數(shù)，K值過大會導致過量分解，使得有效信號頻率體現(xiàn)在多個分量中；過小會導致分解不完全，有效信號混疊在噪音信號中。所以選擇最合適的模態(tài)數(shù)，是運用VMD方法的關鍵。

經(jīng)變分模態(tài)分解后各分量是正交關系，因此各分量的能量和應與原信號的能量相等。但在實際操作過程中往往會因無法正交分解而產(chǎn)生能量的損失，損失的能量越少說明分解的正交性越好。因此，可以通過對信號損失能量進行定量分析選取分解的模態(tài)數(shù)K[29]。

設采集到有限時長L的GPR單道信號為si(i=1,2,3，…，L)，經(jīng)變分模態(tài)分解后其各分量為uk(k=1,2,3,…,K)，則定義該信號總能量為

(5)

分量信號能量之和為

(6)

選取K=2,3,4,…,N，N為K的最大取值，綜合統(tǒng)計規(guī)律和計算量，本文中N=15。此時分析其能量損失比η=(E-EIMF)/E，當η達到最小時，K為最合適的模態(tài)數(shù)。

1.1.2 皮爾遜相關系數(shù)法

相關系數(shù)法是通過計算各分量與原信號之間的皮爾遜相關系數(shù)，來判定有效信號與噪聲信號。為了避免幅值較小的有效信號被去除，需先將所有的分量與原信號進行歸一化處理。其歸一化后的相關系數(shù)為[30]

(7)

以相關系數(shù)的標準差為判斷閾值，其計算式為

(8)

若ri>ξ，則認為第i個分量為有效信號；否則為噪聲信號。

1.2 奇異譜分析

經(jīng)變分模態(tài)分解得到的有效信號中，往往會殘留一部分噪聲，導致信號振蕩。為了進一步去除這部分噪聲信號，改善振蕩現(xiàn)象，引入SSA算法。

SSA是一種處理非線性時間序列數(shù)據(jù)的方法，通過對所要研究的時間序列的軌跡矩陣進行分解，計算時間序列矩陣的特征值，將較小特征值認定為噪聲子空間，而保留較大特征值重構信號，從而實現(xiàn)對噪聲的抑制[27]。

SSA可以概述為以下兩個步驟[31]。

1)分解：對有效信號分量進行SSA去噪。首先需要選擇合適的窗口長度m，通常情況下取m

(9)

式中，g為任一有效信號分量。計算其協(xié)方差矩陣S：

(10)

式中，0≤τ≤m-1。對S進行奇異值分解，得到m個特征值并排序：λ1≥λ2≥…≥λm≥0。其對應的特征向量為U1,U2,…,Um，且有：

(11)

(12)

式中：Uk,j為第k個特征值對應的特征向量中的第j個元素；0

2)重構：制定合適的閾值。選取前v個特征值，對應的{λv}構成新的子集，則可將m×n的矩陣g重構成長度為L的時間序列。具體重構式為：

(13)

式中，xk,i為第k個特征值對應的重構信號中的第i個元素。

1.3 算法步驟

運用VMD和SSA進行信號去噪的算法步驟如下：

1)預先設置K的取值范圍(K=2, 3, …,N，N=15)，計算其對應的η，并選取最小η對應的K值。

2)對含噪信號進行變分模態(tài)分解得到k個本征模函數(shù)分量。

3)將含噪信號及各分量進行歸一化，根據(jù)式(7)計算各分量的皮爾遜相關系數(shù)，并根據(jù)式(8)計算劃分有效信號和噪聲信號的閾值。

4)對有效信號分量進行奇異譜分析，進一步去除中低頻噪聲。

5)將經(jīng)步驟4)處理后的分量進行重構，得到去噪后的信號。

綜上所述，本文提出的基于VMD和SSA的GPR去噪方法的算法流程如圖1所示。

附偽代碼如下。

基于VMD和SSA的GPR去噪方法輸入:GPR信號(M道),VMD參數(shù),SSA參數(shù)輸出:去噪信號for i=1: Mfor K=2: N 變分模態(tài)分解; 計算原信號能量及分量總能量;end計算損失能量比最小時的K值;變分模態(tài)分解;歸一化:原信號及分量;根據(jù)式(7)和式(8)計算相關系數(shù)及閾值;識別并提取有效信號分量;SSA去噪;重構;end

2 數(shù)值實驗

2.1 合成Ricker子波實驗

為驗證本文方法對低頻含噪信號和高頻含噪信號的去噪性能，分別用EEMD方法、VMD方法和本文方法(VMD-SSA)對合成Ricker子波信號進行去噪對比。合成信號由在20、50、80 ns處的Ricker子波信號(100和900 MHz)和高斯白噪聲組成。合成信號表達式s(t)為：

(14)

式中：f0為Ricker子波中心頻率；RN(t)為高斯白噪聲。

為了更好地對比去噪效果，通過計算去噪前后的均方根誤差(root mean squared error,ERMS)和信噪比(signal to noise ratio,RSN)進行定量分析：

(15)

式中，s′為降噪后的信號。RSN越大，且ERMS越小，表明信號中含有效成分越多，去噪效果越好。

圖1 基于VMD-SSA的GPR去噪算法流程圖

為了探究VMD模態(tài)數(shù)的自適應性選擇，在100和900 MHz的Ricker子波合成信號中添加同一噪聲信號，計算其能量損失比隨模態(tài)數(shù)的變化曲線如圖2所示。由圖2可知，隨著K增大，η先減小到最小值，其后總體呈現(xiàn)上升趨勢。其中:對于100 MHz的Ricker子波信號，η值在K為5時達到最??；對于900 MHz的Ricker子波信號，η值在K為3時達到最小。表明對于含相同噪聲的不同信號，依據(jù)能量損失比最小可以得到最適合該信號的K值。

向100 MHz的Ricker子波合成信號中添加不同信噪比的噪聲信號，統(tǒng)計其規(guī)律如表1所示。隨著信號含噪增多，K值呈現(xiàn)上升趨勢，且其對應的η值也逐漸增大。

圖3和圖4分別為100和900 MHz Ricker子波信號去噪效果以及對應的時頻譜圖。對于加噪的Ricker子波信號，EEMD、VMD和本文方法都能有效壓制噪聲，恢復信號。從圖3a、b、c、d和圖4a、b、c、d中可知：針對不同頻率的信號，EEMD和VMD兩種方法去噪效果不穩(wěn)定，表現(xiàn)為對于低頻信號，前者相比于后者振蕩現(xiàn)象更小，更貼合無噪信號，對于高頻信號則相反，且此時EEMD重構信號中含噪更多；而本文方法對于不同主頻Ricker子波信號都能最大程度壓制噪聲，使得信號整體更加光滑，更能凸顯出有效信號。從對應的時頻譜中圖3e--h、圖4e--h可以看出，經(jīng)3種方法去噪后，雜波區(qū)域都有改善，且信號的中心頻率都有恢復。在100 MHz Ricker子波去噪后的時頻譜中，EEMD方法相比于VMD方法噪點更少，在900 MHz Ricker子波去噪后的時頻譜中，結論恰好與之相反；而本文方法無論對于高頻還是低頻信號，雜波區(qū)域噪點都是最少的，中心頻率所對應的能量更加突出，去噪效果最好。這與時域的結論是一致的。結合表2中信噪比分析可知，與EEMD和傳統(tǒng)VMD方法相比，經(jīng)本文方法處理后的信噪比最大提升13.587 8 dB。本文方法在高頻和低頻情況下都得到了最高的信噪比和最低的均方根誤差，表明本文方法去噪效果更好，去噪后的信號更接近原始信號，具備有效性和穩(wěn)定性。

a. 100 MHz Ricker子波；b. 900 MHz Ricker子波。

表1 100 MHz Ricker子波合成信號不同信噪比下的K值及其對應的η值

2.2 合成雷達剖面模擬實驗

建立背景為均勻介質的空間，相對介電常數(shù)設置為6，電導率為0.005 S/m，相對磁導率為1。在介質層中設有3個異常體，分別為空氣管線、金屬管線及含水管線，其半徑為0.1 m，管頂距空氣層底0.2 m。每個異常體模型參數(shù)如表3所示，位置如圖5所示。

設定計算區(qū)域大小為5 m×1 m，頂部空氣層厚度0.1 m，計算網(wǎng)格大小為0.005 m，時間步長為0.01 ns，時間步為2 500步，時窗長度為25 ns。發(fā)射天線和接收天線均置于空氣層，采取自激自收的收發(fā)方式。以布萊克曼-哈里斯脈沖作為激勵源，主頻為900 MHz，模擬采集200道，道間距為0.025 m，通過時域有限差分(finite-difference time-domain, FDTD)方法正演得到GPR剖面、并添加高斯白噪聲后的合成GPR信號如圖6a所示，該剖面背景噪聲強烈，大部分反射回波都被淹沒在噪聲中，難以獲取有效信息，分辨異常體特征。

圖6中b、c、d分別為EEMD、VMD和本文方法去噪結果。經(jīng)EEMD及VMD方法去噪后，可以從剖面中分辨出空氣管線、金屬管線、含水管線的上反射界面回波信息及含水管線的多次回波，但背景噪聲依然十分強烈；本文方法壓制了背景噪聲，使得剖面整體更加平滑，突出異常體信息特征，可以清晰地獲取到空氣管線和金屬管線的微弱回波信息，且含水管線的多次回波信息也更加豐富。

EEMD、VMD和本文方法對合成GPR剖面去噪后的信噪比如表4所示，本文方法處理后剖面的信噪比較EEMD和傳統(tǒng)VMD方法分別提高了3.765 9 dB和2.655 7 dB，獲得了最高信噪比，最大程度壓制噪聲，突出有效回波信息，這與合成Ricker子波實驗結果是一致的。

取本文方法去噪后的GPR剖面第100道展示了更直觀的去噪效果，計算η隨K的變化曲線如圖7所示。當K取3時，η為最小值。

圖8中a、b、c、d分別表示第100道合成GPR信號及3種方法去噪結果，且將3～6 ns的波形進行放大顯示。由圖8b、c可知，EEMD方法和VMD方法在一定程度上都能恢復出該位置的反射回波，但波形振蕩依然嚴重，重構信號中依然含噪較多。由8d可知，經(jīng)本文方法去噪后的信號在反射回波處與無噪信號波形幾乎重合，波形振蕩得到了很好的改善，整體波形更加光滑，去噪效果最好。

a. 加噪信號；b. EEMD去噪信號；c. VMD去噪信號；d. 本文方法去噪信號；e. 加噪信號時頻譜；f. EEMD去噪信號時頻譜；g. VMD去噪信號時頻譜；h. 本文方法去噪信號時頻譜。

a. 加噪信號；b. EEMD去噪信號；c. VMD去噪信號；d. 本文方法去噪信號；e. 加噪信號時頻譜；f. EEMD去噪信號時頻譜；g. VMD去噪信號時頻譜；h. 本文方法去噪信號時頻譜。

表2 三種方法去噪效果對比結果

表3 異常體模型參數(shù)

圖5 正演模型

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

采用中心頻率為400 MHz的GPR天線對廣州市某地地下電力套管的目標區(qū)域進行掃描探測，得到的GPR剖面如圖9a所示。選用點測模式進行測量，道間距為0.05 m，時窗為100 ns。圖9b為使用本文方法去噪后的剖面，圖10為單道去噪效果。由圖9可知，50道左右，40～100 ns處有反射回波信息，可以判斷此處存在套管。自60 ns起，存在較強噪聲干擾。經(jīng)本文方法去噪后，60～80 ns處的噪聲基本去除，80～100 ns處噪聲也明顯受到壓制，圖像中噪聲區(qū)域整體變得平滑，更加突出有效反射回波信息。為了更好地顯示背景噪聲的去除效果，截取第100道和第200道單道數(shù)據(jù)處理結果，VMDK值均為4。由圖10可知：前50 ns無噪?yún)^(qū)域波形基本保持一致，直達波及部分反射回波的信號得到保留；從約60 ns起，雜亂的噪聲波形得到了很好的壓制，達到了去噪的效果。

a. 加噪剖面；b. EEMD去噪后剖面；c .VMD去噪后剖面；d. 本文方法去噪后剖面。

表4 剖面去噪后的信噪比結果

圖7 第100道信號不同模態(tài)數(shù)對應的能量損失比

a. 加噪信號；b. EEMD去噪信號；c. VMD去噪信號；d. 本文方法去噪信號。

a. 實測剖面；b. 本文方法去噪剖面。

a. 第100道；b. 第200道。

4 結論

1)在變分模態(tài)分解過程中，分析了基于能量損失比的分解方法。對于不同中心頻率的Ricker子波信號，隨模態(tài)數(shù)增加，能量損失比呈現(xiàn)先減小后增大的規(guī)律；對于添加不同強度噪聲的同一信號，模態(tài)數(shù)及最小能量損失比均隨信噪比減小而增大。結果表明：針對不同的含噪信號，可以實現(xiàn)模態(tài)數(shù)的自適應選擇，避免了人為選擇的主觀性。

2)本文利用SSA方法對變分模態(tài)分解后的信號進行二次濾波，結果表明，本文方法優(yōu)于EEMD和傳統(tǒng)VMD方法的去噪效果。其中，合成Ricker子波實驗結果中，與EEMD和傳統(tǒng)VMD方法相比，經(jīng)本文方法處理后信噪比最大提高了13.587 8 dB，突出有效信號的同時使得信號整體更加平滑；合成雷達模擬實驗結果中，本文方法處理后信噪比較EEMD和傳統(tǒng)VMD方法分別提高3.765 9 dB和2.655 7 dB，剖面中能夠分辨出更微弱的反射回波信息，且背景噪聲得到更好的壓制。

3)本文方法在實測資料的應用中，能夠有效壓制背景噪聲和雜波干擾，突出有效信號，進一步證實了該方法用于處理GPR信號的適用性和魯棒性，為GPR實際資料的數(shù)據(jù)處理提供了有價值的參考。