徐 智 唐 剛* 劉 偉 李鐘曉

(1.北京化工大學 機電工程學院， 北京 100029； 2.青島大學 電子信息學院， 青島 266071)

引 言

由于人工、環(huán)境和地質(zhì)條件等因素的影響，采集的地震數(shù)據(jù)常含有大量隨機噪聲，這會嚴重影響地震數(shù)據(jù)的后續(xù)處理和解釋，因此有效壓制隨機噪聲對提高地震數(shù)據(jù)的信噪比具有重要意義。

小波分析因其多分辨的特點，在隨機噪聲壓制方面具有較好的效果[1-3]。陳曉玉[4]使用小波閾值去噪降低了地震數(shù)據(jù)中的隨機噪聲。崔少華等[5-6]研究了分解層數(shù)、重構(gòu)個數(shù)、去噪閾值以及閾值處理方式對地震數(shù)據(jù)去噪效果的影響。但受基函數(shù)固定、多分辨率恒定等限制，小波分析缺乏自適應性。

經(jīng)驗模態(tài)分解(empirical mode decomposition, EMD)可以自適應壓制地震數(shù)據(jù)中的噪聲，提升信噪比[7-8]，但EMD方法存在模態(tài)混疊和端點效應等不足，導致其分解結(jié)果不穩(wěn)定、不唯一。集合經(jīng)驗模態(tài)分解[9]和完備集合經(jīng)驗模態(tài)分解[10]能在一定程度上緩解模態(tài)混疊問題，但同時也極大增加了計算量，且端點效應問題仍然存在。為解決這一問題，何元等[11]和康佳星等[12]使用變分模態(tài)分解(variational mode decomposition, VMD)[13]進行地震數(shù)據(jù)處理，結(jié)果都表明VMD方法抗噪性更好，運算效率更高，且去噪效果僅依賴于分解參數(shù)的選取。然而，目前VMD參數(shù)大都依賴經(jīng)驗選取，結(jié)果具有較大隨機性和偶然性。本文提出使用估計信噪比作為評價參數(shù)，自適應地優(yōu)化VMD的參數(shù)分解和有效模態(tài)選取過程，以改善地震數(shù)據(jù)去噪效果。

1 理論方法

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

變分模態(tài)分解是將信號f(t)分解為有限個具有限定帶寬的模態(tài)函數(shù)uk(t)，每個模態(tài)函數(shù)圍繞在該模態(tài)的中心頻率ωk附近，變分問題的解為最小帶寬之和。

1.1.1變分模態(tài)分解模型構(gòu)建

首先對每個模態(tài)分量進行Hilbert變換。

(1)

式(1)中，t為時間，uk(t)為第k個模態(tài)分量。

其次對每個模態(tài)預設(shè)一個中心頻率，將各模態(tài)得到的解析信號頻譜通過移頻方式移至基帶上。

(2)

式(2)中，ωk為第k個模態(tài)預設(shè)的中心頻率。

最后計算解調(diào)信號梯度的范數(shù)，估計出各頻帶的帶寬。

(3)

1.1.2變分模態(tài)分解模型求解

首先引入二次懲罰因子α和拉格朗日乘法算子λ(t)，將約束變分問題轉(zhuǎn)化為無約束變分問題，其表達式為

(t)]e-jωkt2+f(t)-∑kuk(t)2+λ(t),f(t)-

∑kuk(t)

(4)

然后利用乘子交替方向算法(ADMM)求取無約束變分問題的鞍點uk和ωk，即為該變分問題的解。uk和ωk的具體更新步驟如下。

1) 更新uk將式(4)轉(zhuǎn)化為等價最小化問題

(5)

利用二范數(shù)下的Parseval/Plancherel傅里葉等距，通過式(6)在傅里葉域求解式(5)

(6)

使用ω-ωk代替式(6)中ω，結(jié)果如式(7)

(7)

利用重構(gòu)保真項中實數(shù)信號的埃爾米特對稱性將式(7)寫成非負頻率上的半空間積分

(8)

利用式(9)求解式(8)的二次優(yōu)化問題

(9)

2) 更新ωk由于重構(gòu)保真項中不含有ωk，所以式(4)可以簡化為

(10)

將式(10)轉(zhuǎn)換到傅里葉域

(11)

對公式(11)求解如下

(12)

最后迭代拉格朗日算子λ

(13)

1.1.3VMD算法步驟

綜合1.1.1和1.1.2兩節(jié)可得VMD算法的步驟如下：

(2)根據(jù)公式(9)和公式(12)迭代計算uk和ωk；

(3)根據(jù)公式(13)迭代計算λ；

1.2 敏感參數(shù)優(yōu)化

1.2.1敏感參數(shù)

分數(shù)高斯噪聲(fractional Gaussian noise，F(xiàn)GN)是由分數(shù)布朗運動(fractional Brown motion，F(xiàn)BM)表示的時間離散序列。本文采用FGN的數(shù)值模擬實驗對VMD等效濾波特性進行分析[14]。設(shè)xH(m)(m=…，-1,0,1，…)為FGN序列，自相關(guān)函數(shù)pH(l)=E{xH(m)xH(m+l)}滿足式(14)[15]

(14)

式中，H為影響FGN特性的Hurst指數(shù)，σ為FGN的標準差。0

使用式(14)對VMD方法進行實驗可知，VMD分解等價于一系列的帶通濾波器，帶通濾波器主要由分解參數(shù)K進行控制，K越大，帶通濾波器的個數(shù)越多，帶寬越小。由VMD實驗結(jié)果還可以得出，VMD地震數(shù)據(jù)去噪嚴重依賴于模態(tài)分量個數(shù)K的大小。由于K值主要根據(jù)經(jīng)驗設(shè)定，具有偶然性和隨機性；并且與經(jīng)驗模態(tài)分解類似，分解后的模態(tài)中有若干個模態(tài)可分為有效模態(tài)(包含有效信息)和噪聲模態(tài)(包含噪聲信息)，說明對有效模態(tài)的選取也會影響去噪效果。故K的選擇及有效模態(tài)的選取成為VMD地震數(shù)據(jù)降噪處理的關(guān)鍵步驟。

1.2.2估計信噪比

通常采用信噪比來定量評價地震數(shù)據(jù)的去噪效果，但實際地震數(shù)據(jù)由于沒有真實信號可作為參考而無法直接計算其信噪比，所以需要對地震數(shù)據(jù)進行信噪比估計。劉洋等[16]比較研究了疊加法、時域奇異值分解法和頻域奇異值分解法對地震數(shù)據(jù)信噪比的估計效果，結(jié)果表明頻域奇異值分解法在相鄰道的地震數(shù)據(jù)間同時存在振幅變化和相位變化時都可以很好地估計信噪比，對于不同地震數(shù)據(jù)的估計能力優(yōu)于其他兩種方法。故本文選取頻域奇異值分解法作為信噪比估計的方法。

在同時考慮振幅變化和相位變化時，地震信號可以表示為

(15)

(16)

任意兩道的頻域互相關(guān)系數(shù)為

(17)

(18)

對矩陣R進行奇異值分解后可以寫成R=UΛVT，其中U是RRT的特征值對應的特征向量矩陣，Λ=diag(σ1,σ2,…,σN)是RRT特征值的非負平方根按遞減順序組成的對角矩陣；V是RTR的特征值對應的特征向量矩陣。可以證明，矩陣R的奇異值為

(19)

(20)

(21)

故信號的能量可由式(22)計算得出

(22)

估計信噪比可由式(23)計算得出

(23)

式中，Es為信號能量，En為噪聲能量。

1.2.3基于估計信噪比的VMD敏感參數(shù)優(yōu)化

本文建立了基于估計信噪比的變分模態(tài)分解敏感參數(shù)優(yōu)化方法，具體步驟如下。

(1)設(shè)置參數(shù)。設(shè)對照分解模態(tài)數(shù)目K0=0，對照含噪信號估計信噪比為xSNR0。

(2)確定K的范圍。VMD的分量個數(shù)至少應為2，故K的下限取為2；EMD分解時產(chǎn)生模態(tài)混疊且存在殘余分量，其上限值應為EMD分解模態(tài)分量個數(shù)KEMD。故K的取值范圍為[2,KEMD]。

(3)在取值范圍內(nèi)選取K值。使用VMD方法分解地震數(shù)據(jù)得到K個模態(tài)，使用頻域奇異值分解方法求取模態(tài)估計信噪比xSNRm。xSNRm>0時該模態(tài)中的有效信號能量強于噪聲能量，故將該模態(tài)作為有效模態(tài)；xSNRm<0時該模態(tài)中的噪聲能量強于有效信號能量，故將該模態(tài)作為噪聲模態(tài)。將有效模態(tài)疊加得到重構(gòu)信號，使用頻域奇異值分解方法求取重構(gòu)信號估計信噪比xSNRr，并與xSNR0進行比較。若xSNRr>xSNR0，則xSNR0=xSNRr，K0=K；否則K0和xSNR0均保持不變。

(4)遍歷取值范圍內(nèi)的K值后迭代更新K0。

(5)使用K0對地震數(shù)據(jù)進行VMD分解重構(gòu)，重構(gòu)結(jié)果即為基于估計信噪比的變分模態(tài)分解參數(shù)優(yōu)化方法的去噪結(jié)果。

基于VMD參數(shù)優(yōu)化的地震數(shù)據(jù)去噪流程如圖1所示。

圖1 VMD參數(shù)優(yōu)化流程圖Fig.1 VMD parameter optimization flow chart

1.3 去噪質(zhì)量評價指標

為了對比各種方法對模型數(shù)據(jù)去噪處理的效果，以信噪比、互相關(guān)系數(shù)和均方根誤差[17]作為地震去噪質(zhì)量的評價指標。

信噪比(signal- noise ratio，SNR)是真實信號與隨機噪聲能量比值的函數(shù)(下文稱為實際信噪比)，其計算公式為

(24)

式(24)中，Ps和Pn分別代表信號和噪聲的能量，f(i)為原始信號，(i)為重構(gòu)信號，N′為信號長度。xSNR值越大越好。

均方根誤差(root-mean-square error，RMSE)代表重構(gòu)信號同原始信號之間的偏差，其計算公式為

(25)

xRMSE越小越好。

互相關(guān)系數(shù)(mutual relationship，MR)表示原始信號與重構(gòu)信號線性相關(guān)的程度，其計算公式為

(26)

式(25)中，σ(f(t))為原始信號標準差，σ((t)為重構(gòu)信號標準差，cov (f(t),(t))為原始信號和重構(gòu)信號的協(xié)方差。xMR越接近1，相關(guān)程度越高，去噪效果越好。

2 數(shù)值實驗

2.1 估計信噪比

使用合成的二維地震剖面(圖2)對信噪比估計方法進行實驗驗證。分別加入強度為[1 dBW,10 dBW]的隨機噪聲，得到估計信噪比與實際信噪比如表1所示。

圖2 合成的二維地震剖面Fig.2 Synthesized 2D seismic profile intensity

噪聲強度/dBW實際信噪比/dB估計信噪比/dB15.84445.539524.84624.748133.74193.872342.83602.951751.84792.003360.84290.97377-0.1575-0.03358-1.1568-1.13919-2.1484-2.219110-3.1591-3.3638

由表1可以看出，估計信噪比隨噪聲強度的增加而減小，能有效反映地震數(shù)據(jù)的噪聲程度；并且估計信噪比與實際信噪比誤差很小，可以將估計信噪比作為有效模態(tài)選取和參數(shù)優(yōu)化的評價參數(shù)。

2.2 仿真模型數(shù)據(jù)

2.2.1二維剖面

為了驗證本文方法的有效性，首先進行了仿真模型數(shù)據(jù)實驗。該模型數(shù)據(jù)包含1 000道，每道1 000個采樣點，采樣間隔2 ms。原始數(shù)據(jù)如圖3(a)所示，加入隨機噪聲后數(shù)據(jù)如圖3(b)所示。

圖3 二維數(shù)據(jù)各方法去噪結(jié)果對比Fig.3 Comparison of the results of methods for obtaining 2D data

EMD分解模態(tài)數(shù)目為12個，故參數(shù)搜索范圍為[2,12]。使用本文方法進行搜索篩選得到最優(yōu)分解模態(tài)數(shù)K=4，再使用K=4對地震數(shù)據(jù)進行VMD分解重構(gòu)，結(jié)果如圖3(c)所示，其差剖面如圖3(d)所示。比較圖3(c)、(b)可以看出，處理后噪聲明顯減小，有效信號得以凸顯；由圖3(d)看出，差剖面中僅包含噪聲，不含同相軸信息，說明本文方法可以在有效壓制地震信號噪聲的同時完整保留有效地質(zhì)信息。

為了驗證本文方法的優(yōu)越性，同時使用小波閾值去噪方法及EMD去噪方法進行對比實驗。小波閾值去噪使用symmlet 6小波進行信號分解，分解層數(shù)為3[6],利用“3σ”方法選取閾值后通過硬閾值方法處理小波系數(shù)[5]，再對處理后小波系數(shù)進行重構(gòu)，結(jié)果如圖3(e)所示，其差剖面如圖3(f)所示。EMD去噪通過對信號的EMD分解選取包含有效信息的模態(tài)進行重構(gòu)[8]，得到結(jié)果如圖3(g)所示，其差剖面如圖3(h)所示。可以看出，小波分析和EMD均有一定的去噪效果，但仍有部分噪聲無法去除。

綜上可得，參數(shù)優(yōu)化的VMD相較小波變換及EMD方法能夠在進一步壓制噪聲的同時完整保留有效信息，提高地震數(shù)據(jù)的信噪比。

圖4 單道仿真模型數(shù)據(jù)各方法處理結(jié)果對比Fig.4 Comparison of experimental single channel simulation model data processing results

為避免主觀判斷去噪效果可能出現(xiàn)的偏差，通過計算信噪比、互相關(guān)系數(shù)和均方根誤差對去噪效果進行量化分析，結(jié)果如表2所示。由表2可以得出，與含噪地震剖面相比，EMD方法信噪比提升7.886 0 db，小波方法信噪比提升8.518 7 db，而本文方法信噪比提升16.674 1 db；EMD方法均方根誤差降低了1.885 4，占總誤差的59.66%，小波方法均方根誤差降低了1.974 9，占總誤差的62.50%，而本文方法均方根誤差降低了2.696 6，占總誤差的85.34%；EMD方法互相關(guān)系數(shù)提升了0.293 3，占含噪數(shù)據(jù)互相關(guān)系數(shù)誤差的68.40%，小波方法互相關(guān)系數(shù)提升了0.309 1，占含噪數(shù)據(jù)互相關(guān)系數(shù)誤差的72.08%，而本文方法相關(guān)系數(shù)提升了0.407 1，占含噪數(shù)據(jù)互相關(guān)系數(shù)誤差的94.94%。評價指標定量分析結(jié)果表明，小波分析、EMD和參數(shù)優(yōu)化的VMD方法均能壓制地震信號中的噪聲，但參數(shù)優(yōu)化的VMD方法壓制噪聲效果最好。

表2 二維數(shù)據(jù)處理結(jié)果的評價指標

2.2.2單道數(shù)據(jù)

從2.2.1節(jié)仿真模型數(shù)據(jù)中抽取第100道數(shù)據(jù)對本文方法進行驗證。該道數(shù)據(jù)包含1 000個采樣點，采樣間隔2 ms。原始信號如圖4(a)所示，加入隨機噪聲后數(shù)據(jù)如圖4(b)所示。利用參數(shù)優(yōu)化的VMD方法去噪后的單道數(shù)據(jù)如圖4(c)所示，去除的噪聲如圖4(d)所示。由圖4(a)、(b)對比可以看出，含噪信號噪聲很大，有效信號被噪聲淹沒；由圖4(c)、(b)對比可知，經(jīng)參數(shù)優(yōu)化VMD方法處理后的干擾噪聲得到有效壓制，子波信號被較好地恢復，從而顯著增強了信號的信噪比和分辨率。

利用小波變換和EMD進行對比實驗，小波去噪的單道處理結(jié)果和去除的噪聲分別如圖4(e)、(f)所示，EMD去噪的單道處理結(jié)果和去除的噪聲分別如圖4(g)、(h)所示。由圖4各方法處理結(jié)果對比可以看出，參數(shù)優(yōu)化的VMD方法比小波分析和EMD方法的子波恢復能力更強，壓制噪聲的效果更明顯。

圖5 實際地震數(shù)據(jù)處理結(jié)果Fig.5 Actual seismic data processing results

由表3可以得出，EMD方法信噪比提升8.494 5 db，小波閾值去噪方法信噪比提升9.273 0 db，而本文方法信噪比提升16.988 1 db；EMD方法均方根誤差降低了1.967 4，占總誤差的61.25%，小波閾值去噪方法均方根誤差降低了2.074 1，占總誤差的64.57%，而本文方法均方根誤差降低了2.743 9，占總誤差的85.43%；EMD方法互相關(guān)系數(shù)提升了0.291 7，占含噪數(shù)據(jù)互相關(guān)系數(shù)誤差的71.16%，小波閾值去噪方法相關(guān)系數(shù)提升了0.312 7，占含噪數(shù)據(jù)互相關(guān)系數(shù)誤差的76.29%，而本文方法相關(guān)系數(shù)提升了0.39，占含噪數(shù)據(jù)互相關(guān)系數(shù)誤差的95.15%。以上評價指標計算結(jié)果表明參數(shù)優(yōu)化VMD的去噪效果明顯優(yōu)于小波去噪和EMD去噪。

表3 單道模型數(shù)據(jù)處理結(jié)果的評價參數(shù)

2.3 實際數(shù)據(jù)

本文所用實際地震數(shù)據(jù)包含800道，每道數(shù)據(jù)包含800個采樣點，采樣點時間間隔4 ms，原始實際地震數(shù)據(jù)如圖5(a)所示。

設(shè)EMD分解的K值為10，利用本文方法在[2,10]范圍內(nèi)對K進行篩選搜索，求得最優(yōu)K值為5，使用K=5對地震數(shù)據(jù)進行VMD分解重構(gòu)如圖5(b)所示，差剖面如圖5(c)所示。對比圖5(a)、(b)可以看出，實際地震數(shù)據(jù)中的噪聲被去除，有效信號得以保留，去噪效果明顯，顯著提高了信號的信噪比。

通過小波閾值去噪方法處理的結(jié)果如圖5(d)所示，差剖面如圖5(e)所示；通過EMD方法處理的結(jié)果如圖5(f)所示，差剖面如圖5(g)所示。

由圖5(b)、(d)和(f)對比可以看出，小波分析和EMD分解有一定的去噪效果，但是在0～200道處噪聲干擾大，噪聲壓制效果不佳；參數(shù)優(yōu)化的VMD方法的整體噪聲壓制效果明顯優(yōu)于上述兩種方法。通過圖5(c)可以看出，參數(shù)優(yōu)化的VMD方法去除的噪聲中不含有效地質(zhì)信息，同樣證明本文方法可以在保留有效地震信息的同時顯著壓制噪聲，是一種較優(yōu)的地震信號噪聲去除方法。

3 結(jié)論

VMD分解可以有效去除地震數(shù)據(jù)中的噪聲,相較小波分析和EMD方法，VMD壓制噪聲效果更明顯。利用仿真模型數(shù)據(jù)和實際數(shù)據(jù)對多種方法進行實驗驗證的結(jié)果表明，基于估計信噪比的參數(shù)優(yōu)化的VMD方法去噪效果良好，能夠有效克服人為選取參數(shù)的主觀隨機性，提升了VMD分解壓制噪聲效果的穩(wěn)定性。

致謝

在本文完成過程中，哈爾濱工業(yè)大學馬堅偉教授和于四偉博士提供了有益討論和幫助，東方地球物理公司提供了實際地震數(shù)據(jù)，在此一并表示感謝。