基于變分模態(tài)分解算法的探地雷達(dá)信號(hào)去噪研究

冉利民， 李健偉， 杜娟， 程思遠(yuǎn)， 宋二喬， 劉四新

1.中石化經(jīng)緯有限公司，鄭州 450000；2.吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

0 引言

探地雷達(dá)通過接收來自地下目標(biāo)的反射回波探測(cè)地下地層，測(cè)量信號(hào)是一種非線性、非平穩(wěn)信號(hào)，傳統(tǒng)數(shù)據(jù)處理方法已經(jīng)滿足不了人們對(duì)探測(cè)精度越來越高的要求[1]。前人[2-9]對(duì)EMD方法不斷進(jìn)行改進(jìn)，出現(xiàn)了總體經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)、完全總體經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMD)等方法，但仍然存在不足之處，例如容易受到噪聲干擾，存在模態(tài)混疊現(xiàn)象，依賴經(jīng)驗(yàn)性等問題。

1998年，Huang et al.[2]提出了Hilbert-Huang Transform(HHT)，用來處理分析非平穩(wěn)、非線性信號(hào)。經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)是HHT的關(guān)鍵步驟之一。本征模態(tài)函數(shù)(IMF)由EMD方法分解得到，對(duì)IMF進(jìn)行希爾波特變換可以得到具有明確物理意義的瞬時(shí)信號(hào)[3-4]。2009年，Wu et al.[5]提出了EEMD方法，在一定程度上解決了EMD分解過程中出現(xiàn)的模態(tài)混疊問題。王宏等[6]采用了EEMD這種改進(jìn)方法，并將之應(yīng)用到穿墻雷達(dá)數(shù)據(jù)去噪中，相比于EMD方法，模態(tài)混疊現(xiàn)象有了明顯改善。許軍才等[7]把EEMD方法應(yīng)用到探地雷達(dá)正演模擬數(shù)據(jù)去噪中，噪聲得到有效壓制，但仍存在加噪去除不干凈的問題。2011年，Torres et al.[8]通過改進(jìn)EEMD方法，提出了CEEMD方法，有效解決了EEMD加噪去除不干凈的問題。雷文太等[9]為解決CEEMD方法中需要人為篩選有效IMF分量的問題，提出一種基于自動(dòng)反相校正和峰度值比較的雷達(dá)數(shù)據(jù)去噪方法，該方法可以通過獨(dú)立分量分析來篩選IMF分量。EMD方法雖然被不斷研究改進(jìn)，但仍存在一些缺點(diǎn)。例如EMD方法很容易受到噪聲干擾，致使得到錯(cuò)誤的IMF分量；EMD方法會(huì)出現(xiàn)模態(tài)混疊現(xiàn)象，這種情況下，IMF分量就失去了物理意義；EMD方法分解得到的IMF分量，沒有嚴(yán)格的數(shù)學(xué)定義，依賴經(jīng)驗(yàn)性。

2014年，Dragomiretskiy et al.[10]提出VMD算法。從原理上看，EMD方法是基于時(shí)域遞歸分解的，而VMD方法是基于頻率完全非遞歸分解的，因此VMD方法能較好地解決EMD方法存在的問題[11-12]。VMD方法使用迭代方法來搜索變分模型的最優(yōu)解，從而確定IMF分量的帶寬和中心頻率。該方法能自適應(yīng)地剖分信號(hào)頻域，從而使IMF被有效地分離[13]，這樣得到的IMF具有很好的魯棒性[14]。VMD現(xiàn)在已經(jīng)是一種很成熟的方法，被應(yīng)用到眾多領(lǐng)域中，例如大壩信號(hào)監(jiān)測(cè)、零件故障檢測(cè)等[12-20]，在地球物理領(lǐng)域也有應(yīng)用[21-27]。在地震數(shù)據(jù)處理中，Xue et al.[21]使用VMD方法將地震信號(hào)分解成若干IMF分量，這些分量又是非負(fù)平滑變化瞬時(shí)頻率AM-FM信號(hào)，具有窄帶特征。經(jīng)過實(shí)際數(shù)據(jù)和模擬數(shù)據(jù)對(duì)比發(fā)現(xiàn)，相比于EMD方法，VMD方法具有更強(qiáng)的局部分解能力，且對(duì)噪聲更加敏感。相比于小波變換等變換方法，VMD方法可以得到更高空間和時(shí)間分辨率的瞬時(shí)頻譜，從而更有利于識(shí)別煤炭地震數(shù)據(jù)中的層位信息[21]。張杏莉等[22]把能量熵理論與VMD方法進(jìn)行了結(jié)合，提出了一種自適應(yīng)的微震噪聲去除方法。相比于EMD方法，該方法去噪效果得到了顯著的提高。

在探地雷達(dá)數(shù)據(jù)處理中，Zhang et al.[26]使用VMD方法把探地雷達(dá)信號(hào)分解成IMF分量，通過實(shí)驗(yàn)室模擬數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)驗(yàn)證，表明VMD方法可以有效地應(yīng)用在雷達(dá)數(shù)據(jù)噪聲處理中，從而可以更清楚地識(shí)別地下異常體。許軍才等[27]在VMD方法的基礎(chǔ)上，使用樣本熵來選擇高階模量，從而成功地去除了探地雷達(dá)中的白噪聲，降噪效果得到了明顯的改善。Cheng et al.[28]把VMD方法應(yīng)用到探冰雷達(dá)中，有效消除了噪音影響，獲得了更加清晰的冰層結(jié)構(gòu)。因此，筆者將VMD方法引入到機(jī)載探地雷達(dá)數(shù)據(jù)處理中。

本文首先介紹VMD方法的原理，然后通過對(duì)正余弦信號(hào)與高斯白噪聲合成的含噪聲信號(hào)分別用 EMD、EEMD、CEEMD和VMD算法進(jìn)行分解，分析對(duì)比證明VMD算法的優(yōu)越性。通過時(shí)域有限差分方法(FDTD)對(duì)一個(gè)水平層狀冰蓋模型進(jìn)行模擬得到合成數(shù)據(jù)，并使用VMD方法對(duì)其進(jìn)行去噪處理，最后使用VMD方法對(duì)一個(gè)機(jī)載雷達(dá)實(shí)測(cè)數(shù)據(jù)進(jìn)行處理，經(jīng)過處理的探地雷達(dá)信號(hào)峰值信噪比從15.4 dB增加到20.3 dB。

1 變分模態(tài)分解算法(VMD)

1.1 VMD 基本原理

VMD算法將一個(gè)復(fù)雜的解析信號(hào)f分解為預(yù)設(shè)尺度的k個(gè)具有特定稀疏特性的固有模態(tài)函數(shù)IMF[10-11]。IMF可以表示為：

uk(t)=Ak(t)cos(φk(t))

(1)

假設(shè)IMF在中心頻率ωk(t)附近是有限帶寬的，則通過估計(jì)每個(gè)IMF的帶寬，建立使IMF帶寬之和最小的變分約束模型：

(2)

式中：uk為分解得到的k個(gè)IMF分量，ωk為uk的中心頻率，?t為對(duì)函數(shù)求時(shí)間的偏導(dǎo)數(shù)，δ(t)為單位脈沖函數(shù)，j為虛數(shù)單位，*表示卷積。

引入二次罰函數(shù)項(xiàng)和拉格朗日乘法算子來求解該變分約束模型(式(2))，因而獲得了擴(kuò)展的Lagrange表達(dá)式：

L({uk},{ωk},λ)=

(3)

式中：α為帶寬參數(shù)，通過設(shè)定足夠大的正數(shù)，使信號(hào)的重構(gòu)精度得到保證，λ為L(zhǎng)agrange乘子。式(3)可以通過交替算子法來求解，頻率域上每個(gè)IMF的表達(dá)式為：

(4)

1.2 VMD主要步驟

(5)

(6)

(7)

式中：τ為噪聲容限參數(shù)。當(dāng)信號(hào)中含有強(qiáng)噪聲時(shí)，為了獲得良好的降噪效果，可以設(shè)定τ=0。VMD的詳細(xì)完整算法可在[10-11]中找到。基于模式更新的維納濾波使得VMD算法對(duì)噪聲具有較強(qiáng)的魯棒性。

1.3 含噪信號(hào)分析檢驗(yàn)

首先合成一個(gè)含噪信號(hào)，由2個(gè)頻率分別為5 Hz、20 Hz的正弦波信號(hào)和高斯白噪聲合成的信號(hào)組成。然后應(yīng)用 EMD、EEMD、CEEMD和VMD算法對(duì)含噪信號(hào)進(jìn)行分解得到IMF分量。最后，為更好地比較信噪分離效果，把分解得到的IMF分量用頻譜圖(圖1)和瞬時(shí)能譜圖(圖2)表示。由圖1可以看到，5 Hz的頻率成分被分解到IMF1分離中，20 Hz的頻率成分被分解到IMF2分量中，從低頻到高頻信號(hào)被依次分解，模態(tài)混疊現(xiàn)象在VMD方法中被有效地避免。從圖2可以看到，VMD方法比另外3種方法有更高的時(shí)間頻率精確度和分辨率。

a.EMD；b.EEMD；c.CEEMD；d.VMD。圖1 各方法的IMF分量頻譜圖Fig.1 IMF component spectrogram of each method

a.EMD；b.EEMD；c.CEEMD；d.VMD。圖2 各方法的IMF分量對(duì)應(yīng)的瞬時(shí)能譜圖Fig.2 Instantaneous energy spectra corresponding to IMF component of each method

2 探地雷達(dá)合成數(shù)據(jù)處理

我們?cè)O(shè)計(jì)一個(gè)由冰下基巖、冰蓋以及空氣組成的模型。模型的大小設(shè)定為 45 m×40 m，空氣和冰面分界線為0 m，-5～0 m為空氣層，0～35 m為冰層(圖3)。

正演方法選用時(shí)域有限差分(FDTD)。天線主頻選擇300 MHz，采樣間隔0.101 ns，時(shí)間窗選擇400 ns。網(wǎng)格尺寸為0.05 m×0.05 m，道間距為0.375 m，采集道數(shù)為121道。天線置于距地面 5 m 的空氣中，收發(fā)距為0.1 m。

使用模擬數(shù)據(jù)(圖4a)進(jìn)行VMD分解的處理數(shù)據(jù)。設(shè)定k=5，可以得到5個(gè)IMF分量，有效信息在圖4c和圖4f中，噪聲信息主要在圖4b、圖4d和圖4e中。

a.介電常數(shù)模型；b.電導(dǎo)率模型。圖3 合成數(shù)據(jù)模型(考慮實(shí)際情況增加了隨機(jī)介質(zhì)的干擾)Fig.3 Synthetic data model

圖4 試驗(yàn)數(shù)據(jù)圖和各IMF分量圖Fig.4 Test data diagram and each IMF component diagram

表1為實(shí)驗(yàn)數(shù)據(jù)與各個(gè)IMF分量的相關(guān)系數(shù)，用于定量分析各分量的有效信息。通過在模擬數(shù)據(jù)中加入隨機(jī)等效介質(zhì)來產(chǎn)生噪聲。當(dāng)有效信息占主體，噪聲信息比例較小時(shí)，定量分析以各分量與實(shí)驗(yàn)數(shù)據(jù)的相關(guān)系數(shù)為標(biāo)準(zhǔn)。如果IMF分量的有效信息占比較少，那么其相關(guān)系數(shù)也會(huì)越小，從而剔除該IMF分量。反之，如果IMF分析的有效信息占比較多，那么相關(guān)系數(shù)會(huì)比較大，該IMF分量需要保留下來。因此，筆者以相關(guān)系數(shù)>0.6為標(biāo)準(zhǔn)，對(duì)IMF分量進(jìn)行重構(gòu)。

表1 實(shí)驗(yàn)數(shù)據(jù)與IMF的相關(guān)系數(shù)

從表1可以看到，IMF2和IMF5的相關(guān)系數(shù)滿足標(biāo)準(zhǔn)，對(duì)這兩個(gè)分量進(jìn)行重組，如圖5所示。對(duì)比圖5和圖4a可以看到，剖面中的噪聲明顯減少，整體來看更加清晰。因此通過剔除噪聲占比較多的IMF分量，重組有效信息占比較多的IMF分量，可以明顯地提高信噪比。經(jīng)過上面的對(duì)比，可以證明VMD方法在探地雷達(dá)信號(hào)去噪中的有效性。

圖5 探地雷達(dá)模擬數(shù)據(jù)IMF2與IMF5分量組合圖Fig.5 Combination of IMF2 and IMF5 components of ground penetrating radar simulation data

3 探地雷達(dá)實(shí)測(cè)數(shù)據(jù)處理

實(shí)測(cè)數(shù)據(jù)采用的是第32次中國(guó)南極科考期間，HiCARS機(jī)載雷達(dá)測(cè)量的F13R47a測(cè)線的一部分?jǐn)?shù)據(jù)。該實(shí)測(cè)數(shù)據(jù)主頻為60 MHz，每道采樣點(diǎn)數(shù)為2 000個(gè)點(diǎn)，總共有1 000道。對(duì)其進(jìn)行一些常規(guī)探地雷達(dá)數(shù)據(jù)處理，包括去直流、減平均等，其中帶通濾波頻率分別為1 MHz、20 MHz、150 MHz和300 MHz，常速度偏移中速度為0.167 m/ns。之后，使用VMD方法分解處理完畢后的數(shù)據(jù)。

圖6 使用常規(guī)探地雷達(dá)數(shù)據(jù)處理方法處理后的數(shù)據(jù)(虛線表示的是測(cè)試道位置)Fig.6 Data processed by conventional ground penetrating radar data processing method

從圖6可以看到，即使經(jīng)過上述常規(guī)探地雷達(dá)數(shù)據(jù)處理，實(shí)測(cè)數(shù)據(jù)的信噪比還是比較低，存在較明顯的干擾信號(hào)，這對(duì)等時(shí)層判斷造成了嚴(yán)重的干擾。使用VMD方法對(duì)第405道實(shí)測(cè)數(shù)據(jù)進(jìn)行處理，設(shè)定k=5，可以得到5個(gè)IMF分量(圖7)。同樣，使用VMD方法對(duì)整個(gè)剖面數(shù)據(jù)進(jìn)行處理，k仍為5，得到5個(gè)IMF分量圖(圖8)。從圖8可以看到，有效信息在圖8c和圖8f中，噪聲信息主要在圖8b、圖8d和圖8e中。

采取與合成數(shù)據(jù)中一樣的評(píng)價(jià)方法，使用相關(guān)系數(shù)對(duì)各IMF分量中有效信息進(jìn)行評(píng)價(jià)?；谠搶?shí)測(cè)數(shù)據(jù)為有效數(shù)據(jù)的前提下，那么有效信號(hào)占比要明顯大于噪聲信號(hào)。因而把實(shí)測(cè)數(shù)據(jù)與IMF分量的相關(guān)系數(shù)作為評(píng)價(jià)標(biāo)準(zhǔn)是可行的。表2為實(shí)測(cè)數(shù)據(jù)與各個(gè)IMF分量的相關(guān)系數(shù)，計(jì)算獲得的相關(guān)系數(shù)與根據(jù)IMF分量得到的結(jié)論是一致的。

圖7 第405道數(shù)據(jù)和各IMF分量圖Fig.7 The 405th channel data and each IMF component

圖8 原始實(shí)測(cè)數(shù)據(jù)與各IMF分量Fig.8 Original measured data and each IMF component

表2 原始實(shí)測(cè)數(shù)據(jù)與IMF的相關(guān)系數(shù)

同樣以相關(guān)系數(shù)>0.6作為標(biāo)準(zhǔn)，選擇IMF1和IMF4分量進(jìn)行重組，重組后的效果見圖9。經(jīng)過VMD分解，再挑選合適的IMF分量進(jìn)行重組的結(jié)果剖面，要比原始實(shí)測(cè)數(shù)據(jù)的剖面更加清晰，可以較好地識(shí)別等時(shí)層與基巖面(圖9)。我們使用峰值信噪比(PSNR)定量評(píng)價(jià)該方法的去噪能力。重組后數(shù)據(jù)的PSNR是20.3 dB，而原始數(shù)據(jù)PSNR只有15.4 dB，存在明顯的差距。因此，可以進(jìn)一步證實(shí)VMD方法在雷達(dá)數(shù)據(jù)去噪方面的有效性。

4 結(jié)論

(1)VMD方法可以有效避免模態(tài)混疊現(xiàn)象，相比EMD、EEMD和CEEMD等方法有更高的時(shí)間頻率精確度和分辨率。

(2)VMD方法相比其他3種方法，可以更好地識(shí)別探地雷達(dá)數(shù)據(jù)中等時(shí)層與基巖面。原始數(shù)據(jù)只有15.4 dB的峰值信噪比，經(jīng)過VMD方法處理，峰值信噪比可以提高到20.3 dB。

圖9 原始數(shù)據(jù)剖面(a)與IMF1和IMF2重組后的剖面(b)Fig.9 Original data section (a) and reorganized section of IMF1 and IMF2(b)