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

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

0 引言

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

1998年，Huang et al.[2]提出了Hilbert-Huang Transform(HHT)，用來處理分析非平穩(wěn)、非線性信號。經(jīng)驗?zāi)B(tài)分解(EMD)是HHT的關(guān)鍵步驟之一。本征模態(tài)函數(shù)(IMF)由EMD方法分解得到，對IMF進(jìn)行希爾波特變換可以得到具有明確物理意義的瞬時信號[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á)數(shù)據(jù)去噪方法，該方法可以通過獨立分量分析來篩選IMF分量。EMD方法雖然被不斷研究改進(jìn)，但仍存在一些缺點。例如EMD方法很容易受到噪聲干擾，致使得到錯誤的IMF分量；EMD方法會出現(xiàn)模態(tài)混疊現(xiàn)象，這種情況下，IMF分量就失去了物理意義；EMD方法分解得到的IMF分量，沒有嚴(yán)格的數(shù)學(xué)定義，依賴經(jīng)驗性。

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

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

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

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

1.1 VMD 基本原理

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

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

(1)

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

(2)

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

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

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

(3)

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

(4)

1.2 VMD主要步驟

(5)

(6)

(7)

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

1.3 含噪信號分析檢驗

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

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分量對應(yīng)的瞬時能譜圖Fig.2 Instantaneous energy spectra corresponding to IMF component of each method

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

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

正演方法選用時域有限差分(FDTD)。天線主頻選擇300 MHz，采樣間隔0.101 ns，時間窗選擇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個IMF分量，有效信息在圖4c和圖4f中，噪聲信息主要在圖4b、圖4d和圖4e中。

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

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

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

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

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

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

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

實測數(shù)據(jù)采用的是第32次中國南極科考期間，HiCARS機載雷達(dá)測量的F13R47a測線的一部分?jǐn)?shù)據(jù)。該實測數(shù)據(jù)主頻為60 MHz，每道采樣點數(shù)為2 000個點，總共有1 000道。對其進(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ù)(虛線表示的是測試道位置)Fig.6 Data processed by conventional ground penetrating radar data processing method

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

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

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

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

表2 原始實測數(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ù)據(jù)的剖面更加清晰，可以較好地識別等時層與基巖面(圖9)。我們使用峰值信噪比(PSNR)定量評價該方法的去噪能力。重組后數(shù)據(jù)的PSNR是20.3 dB，而原始數(shù)據(jù)PSNR只有15.4 dB，存在明顯的差距。因此，可以進(jìn)一步證實VMD方法在雷達(dá)數(shù)據(jù)去噪方面的有效性。

4 結(jié)論

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

(2)VMD方法相比其他3種方法，可以更好地識別探地雷達(dá)數(shù)據(jù)中等時層與基巖面。原始數(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)