基于變分模態(tài)分解的地震隨機噪聲壓制方法

方江雄 溫志平 顧華奇 劉 軍 張 華

(①東華理工大學(xué)核技術(shù)應(yīng)用教育部工程研究中心，江西南昌 330013；②東華理工大學(xué)地球物理與測控技術(shù)學(xué)院，江西南昌 330013；③江西省基礎(chǔ)地理信息中心，江西南昌 330209)

0 引言

地震信號隨機噪聲衰減是地震數(shù)據(jù)處理領(lǐng)域的熱點和難點問題。從觀測數(shù)據(jù)中有效去除隨機噪聲干擾、提高信號信噪比和分辨率，是正、反演計算和地質(zhì)解釋的前提。目前，學(xué)者們提出了大量的地震信號隨機噪聲壓制方法，比如基于濾波理論的方法[1-2]、基于小波域變換方法[3-6]、基于矩陣?yán)碚摰淖儞Q方法[7-10]、基于信號分解理論方法[11-16]等。近年來，由于在處理非平穩(wěn)及非線性數(shù)據(jù)上具有很高的信噪比，基于信號分解的經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)方法[17]已成為信號分解領(lǐng)域研究的熱點。

學(xué)者們在EMD方法基礎(chǔ)上，進一步提出兩種經(jīng)典的改進方法，即集合經(jīng)驗?zāi)B(tài)分解(Ensemble EMD，EEMD)[18]和完備集合經(jīng)驗?zāi)B(tài)分解(Complete EEMD，CEEMD)[19]。其中，EEMD避免了EMD方法中的模態(tài)混疊，而CEEMD改善了EEMD方法中因輔助噪聲造成的噪聲污染。但是，CEEMD方法中依然采用遞歸迭代式篩選分解過程，非平穩(wěn)地震信號極值點插值和包絡(luò)計算過程耗時過長，在處理多維度和多尺度地震數(shù)據(jù)時存在一定限制。為了解決這些問題，以進一步提升信號分解的精度，本文提出全新的非線性、非遞歸的完全自適應(yīng)時序序列分解算法——基于變分模態(tài)分解(Variational Modal Decomposition，VMD)的地震信號分解方法，并應(yīng)用于二維和三維地震數(shù)據(jù)隨機噪聲壓制中。

1 基本原理

VMD的基本思路是：在EMD的固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)分量基礎(chǔ)上，引入更嚴(yán)格的頻率域帶寬先驗約束(Band-Limited Intrinsic Mode Function，BIMF)信號，又稱帶限分量，解決EMD中模態(tài)混疊和高頻率信號丟失問題。IMF被定義為在信號極值點個數(shù)和過零點個數(shù)相等或相差一個時，任一時刻點由局部極值點所確定包絡(luò)均值為零的信號。

VMD基本原理是： 在IMF信號分量基礎(chǔ)上，通過引入調(diào)幅—調(diào)頻(AM-FM)信號

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

(1)

式中：Ak(t)和φk(t)分別為uk(t)的瞬時幅值和瞬時相位。φk(t)為非減函數(shù)，即瞬時角頻率

(2)

與φk(t)相比，Ak(t)和ωk(t)的變化較為緩慢，即在短時范圍內(nèi)可認(rèn)為是一個幅值和頻率不變的諧波信號。將滿足上述條件的BIMF賦予具體的物理意義，即VMD將輸入信號自適應(yīng)分解成有限個BIMF分量之和，分解過程滿足每個帶限BIMF分量的估計聚集帶寬之和為最小。與傳統(tǒng)EMD方法不同，在獲取BIMF分量時，VMD方法的遞歸迭代篩選模式將信號的分解過程轉(zhuǎn)換至變分問題的求解，在尋找變分模型最優(yōu)解的過程中完成信號的自適應(yīng)分解。

2 地震隨機噪聲壓制

2.1 地震信號的域表示

線性地震數(shù)據(jù)可表示為

(3)

式中：t和x為時間和炮檢距軸；w為地震子波函數(shù)，如Ricker子波；c為波在介質(zhì)中的傳播速度。將同相軸數(shù)據(jù)沿t(時間軸)作傅里葉變換，將t-x域轉(zhuǎn)換到f-x域，得到原始數(shù)據(jù)的f-x頻譜

(4)

將式(4)進行離散化，即可得到D(f-x)的離散表示

Df(m)≡D(f,mΔx)m=1,2,…,M

(5)

式中： Δx表示沿x軸方向的道間距；m為地震采樣道數(shù)，M為總道數(shù)。統(tǒng)計相鄰兩道可表示為

(6)

式(6)是一階遞歸微分方程，對應(yīng)每一個頻率切片Df(m)信號，由一個復(fù)調(diào)和諧波函數(shù)(調(diào)幅—調(diào)頻)組成。因此，通過在t-x域中的無噪聲線性數(shù)據(jù)中找到遞歸系數(shù)ei2πfΔx/c，使其在f-x域中完全可預(yù)測。類似地，t-x域中N個線性事件的疊加相當(dāng)于f-x域中的N個復(fù)數(shù)諧波的疊加。

2.2 VMD框架求解

基于上述地震數(shù)據(jù)域分析，在頻率域中估計BIMF分量頻率帶寬目標(biāo)函數(shù)[20]

(7)

通過將VMD方法的變分優(yōu)化框架(式(7))拓展至復(fù)數(shù)空間，而復(fù)數(shù)空間正弦信號的傅里葉頻譜是單邊的，因此式(7)改寫為

(8)

為求解式(8)中變分問題，采用增廣Lagrange函數(shù)[14]計算其最優(yōu)解。通過引入二次懲罰因子α和拉格朗日乘法算子λ，將約束性變分問題轉(zhuǎn)換為非約束性變分形式

(9)

式中：L(·)表示目標(biāo)函數(shù)；“∶”是一個函數(shù)的標(biāo)記；〈·,·〉表示內(nèi)積；二次懲罰因子α是控制數(shù)據(jù)保真度的均衡參數(shù)，用于平衡變分正則項和二次約束項，在含噪聲情形時可保證信號重構(gòu)精度；λ可以保證模型約束條件的嚴(yán)格性。

采用交替方向乘子算法 (Alternate Direction Method of Multipliers，ADMM)[22]求解式(9)，具體步驟如下：

(2)n=n+1，執(zhí)行主循環(huán)；

(3)Fork=1∶K-1，執(zhí)行第一個內(nèi)循環(huán)更新uk

(10)

(4)k=K，結(jié)束第一個內(nèi)循環(huán)；

(5)Fork=1∶K-1，執(zhí)行第二個內(nèi)循環(huán)更新ωk

(11)

(6)k=K，結(jié)束第二個內(nèi)循環(huán)；

(7)對于所有ωk>0,k∈[1,K]， 雙重對偶上升，更新λ

(12)

式中τ表示噪聲容限參數(shù)。在噪聲壓制任務(wù)中(而不是信號的重構(gòu))，更新參數(shù)τ=0以得到更好的去噪效果；

(8)給定判定精度ε>0,k∈[1,K]， 重復(fù)步驟(2)～(7)，直至滿足迭代停止條件

(13)

(9)結(jié)束迭代，即得到K個帶限BIMF分量。

從上述步驟可知，ADMM的求解過程包含VMD的模態(tài)和頻率中心更新。其中，ωk的更新由對應(yīng)模態(tài)的能量譜重新得到，uk模態(tài)更新對應(yīng)于1/βω2的Wiener濾波器結(jié)構(gòu)(β為白噪聲方差，1/ω2表示信號的能量譜為低通形式)。參數(shù)α控制著Wiener濾波器的寬度，本文稱為保真度均衡參數(shù)。增大α值，Wiener濾波器寬度變窄，可以濾除更多的噪聲，但也使得BIMF分量包含更少的真實信號峰值信息，同時算法趨于發(fā)散不收斂的幾率增加，反之亦然。式(7)的增廣Lagrange函數(shù)形式與Tikhonov正則一致，具有保真項和正則項。其中，VMD方法中分解的帶限BIMF分量具有明顯的稀疏特性，因而正則項可用稀疏先驗。α與正則參數(shù)類似，用于控制保真項和正則項的權(quán)重，使BIMF模態(tài)分量的更新公式(式(10)～式(12))具備Wiener濾波特性，這正是VMD方法具有較強抗噪性能的關(guān)鍵因素。

2.3 基于VMD的二維地震信號去噪方法

根據(jù)上述理論，基于VMD方法的二維地震信號隨機噪聲壓制處理算法流程如下：

(1)選擇一個時間窗口將原始含噪地震信號d(x,t)作傅里葉變換至f-x域；

(2)對每一個頻率切片數(shù)據(jù)進行復(fù)數(shù)VMD分解；

(3)將VMD分解得到的BIMF分量組合得到濾波后的信號；

(4)將信號作傅里葉逆變換回t-x域；

(5)對下一個時間窗口重復(fù)以上操作，地震記錄全部處理完畢后，即得到最終的二維地震去噪結(jié)果。

2.4 基于VMD的三維地震信號去噪方法

在三維地震數(shù)據(jù)中，平面波可表示為

(14)

式中：x=(m,h)為中心點和炮檢距坐標(biāo)；n為波的傳播單位方向。三維地震模型的頻率切片Df是由沿n方向的一個f-m-h域復(fù)諧波模式組成。將d(t,x)沿時間軸上作傅里葉變換得到f-m-h譜，每個頻率切片數(shù)據(jù)即為Df。相應(yīng)地t-m-h域中p個平面波的疊加與f-m-h域中p個復(fù)諧波疊加等效。本文將二維VMD拓展到復(fù)值，拓展后的目標(biāo)優(yōu)化函數(shù)為

(15)

式中vk為二維平面的波數(shù)向量。與一維VMD求解類似，通過引入二次懲罰和拉格朗日乘數(shù)(增廣Lagrange)重建約束變分框架，由ADMM優(yōu)化求解，二維VMD的三維地震隨機噪聲壓制算法步驟如下：

(1)選擇一個時間窗口將原始含噪地震信號d(m,h,t)作傅里葉變換至f-m-h域；

(2)對每一個頻率切片數(shù)據(jù)進行二維復(fù)數(shù)VMD分解；

(3)組合VMD分解的BIMF分量得到濾波后的信號；

(4)將信號作傅里葉逆變換回t-m-h域；

(5)對下一個時間窗口重復(fù)以上操作；

(6)地震記錄全部處理完畢后，即得到最終的三維地震去噪結(jié)果。

3 算例

3.1 VMD參數(shù)優(yōu)選

VMD變分模型目標(biāo)函數(shù)是非凸的，算法的收斂性與參數(shù)的設(shè)置有密切聯(lián)系，有三個重要參數(shù)，即帶限BIMF分量個數(shù)K、ωk和α。其中，K依據(jù)BIMF分量瞬時頻率均值的極大值點確定；ωk采用匹配追蹤算法(Matching Pursuits，MP)[23]確定；α與原始信號的噪聲水平有關(guān)，用來控制保真項和正則項的權(quán)重，本文采用“L”曲線法則確定α。為驗證VMD方法的有效性，本文將三個不同中心頻率余弦信號累加構(gòu)成合成信號s(t)

(16)

圖1a、圖1c、圖1e顯示了合成信號s(t)經(jīng)VMD分解得到的3個BIMF分量信號，圖1b、圖1d、圖1f分別為上述相應(yīng)的BIMF分量信號局部放大結(jié)果。從分解結(jié)果可以看出，各個BIMF分量信號與原始組分信號幾乎一致，僅在信號兩端點處出現(xiàn)微小誤差。

圖2為原始合成信號的頻譜分布與VMD分解后的3個BIMF分量信號的頻譜圖(雙對數(shù)坐標(biāo))，每個重構(gòu)BIMF模態(tài)分量的頻譜都有一個最高峰值，對應(yīng)其中心頻率，與原始信號的期望中心頻率2、24、288Hz高度一致(非等間距橫坐標(biāo))，VMD分解得到的3個BIMF分量均能成功捕捉對應(yīng)中心頻率。

圖1 各組分信號及VMD分解BIMF分量

圖3為分解各模態(tài)信號的頻率中心隨迭代次數(shù)的變化曲線(半對數(shù)坐標(biāo))。由圖可見，圖3a迭代更新過程較慢，在150次后才找到最優(yōu)的中心頻率，而圖3b經(jīng)3次迭代即可獲得最佳頻率中心分布。

圖2 VMD分解各模態(tài)頻譜分布

圖3 VMD分解各模態(tài)頻率中心變化

圖4 合成地震數(shù)據(jù)(左)及加噪數(shù)據(jù)(右)頻譜對比

圖4是一組由四個線性事件組成的合成地震數(shù)據(jù)、加噪數(shù)據(jù)及其對應(yīng)的f-x頻譜圖、f-k頻譜圖。對圖4b中第40個采樣點切片含噪數(shù)據(jù)作VMD分解，結(jié)果如圖5所示。設(shè)置參數(shù)α=2000，K=4，以匹配追蹤算法初始化ωk分解得到四個近似諧波信號，各個BIMF分量的頻譜與原始信號的頻譜峰值高度一致。

3.2 二維地震去噪實例

本文采用信噪比(SNR)、 局部相似度(Local Similarity，LS)[24]、 結(jié)構(gòu)相似度(Structural Simila-rity Index，SSIM)[25]評價算法的有效性。局部相似度定義為

(17)

式中： T表示轉(zhuǎn)置；c1、c2通過最小二乘問題最小化得到

圖5 含噪f-x譜第40個采樣點切片數(shù)據(jù)作VMD分解

(18)

(19)

式中：a和b為矩陣；A、B分別為a、b構(gòu)成的對角算子。LS值變化范圍為[0，1]，LS值越大表示去噪結(jié)果與噪聲剖面的局部相似程度越高，有效能量泄漏越嚴(yán)重，算法的幅值保持能力越差。

結(jié)構(gòu)相似度定義為

(20)

式中：μx、μy分別表示圖像x、y的均值；σx、σy分別表示圖像x、y的標(biāo)準(zhǔn)差；σx、σy分別表示圖像x、y的方差；σxy代表圖像x和y的協(xié)方差；C1、C2為常數(shù)，C1=(k1L)2，C2=(k2L)2，L為像素動態(tài)范圍，一般k1=0.01、k2=0.03、L=255；SSIM值變化范圍為[0，1]，值越大表示兩個圖像的結(jié)構(gòu)相似程度越高，去噪性能越強。

圖6表示一個含有多層傾斜地層、一個不整合面、一條斷層和多套正弦起伏形態(tài)地層組成的理論合成地震記錄Sigmoid模型[26]，含256道，每道有256個時間采樣。為驗證VMD方法的有效性，本文將VMD與Curvelet(曲波)、NLM(非局部均值)、BM 3D(三維塊匹配)、CEEMD四種方法進行對比。

圖6 合成地震記錄Sigmoid模型

在本實驗中，VMD方法中K=4，α=2000，ωk初始化采用匹配追蹤方法，時間軸窗口長度取64個采樣點，空間軸窗口長度取64道，時間軸和空間軸的滑動步長重疊度均設(shè)置為50%；CEEMD方法集總次數(shù)I取100；Curvelet、NLM和BM3D三種方法采用工具包中的默認(rèn)參數(shù)。算法運行平臺參數(shù)為MATLAB 2017a Parallel Pool×8、Win 10×64、i7-7700K CPU 4.2GHz。

圖7a～圖7f左圖為地震數(shù)據(jù)加噪后用不同方法去噪的結(jié)果，右圖為相應(yīng)的噪聲剖面局部相似度。由圖可知，五種方法均達到一定的去噪效果，但Curvelet方法去噪結(jié)果有一定的變形，在正弦狀起伏的強反射同相軸周圍發(fā)生明顯偽影現(xiàn)象，斷層現(xiàn)象不明顯，弱反射同相軸連續(xù)性差，對應(yīng)的局部相似度圖表明有效能量泄漏明顯(圖7b)；NLM去噪結(jié)果中仍保留有肉眼可視的噪聲信息，對應(yīng)的局部相似度圖表明線性和雙曲線有效能量均出現(xiàn)泄漏，對原圖像的結(jié)構(gòu)信息保護不夠，這與NLM算法中在計算圖像塊相似性時只考慮塊的平移特性有關(guān)(圖7c)；相對Curvelet和NLM方法，CEEMD方法在SNR數(shù)值上有些提升，但去噪結(jié)果中仍保留部分噪聲信息，斷層部分尚可辨識，對應(yīng)的局部相似度圖在正弦狀起伏的強反射同相軸有些許能量損失(圖7d)；BM 3D方法的去噪效果較好，但也可見部分低能量有效信號被當(dāng)做噪聲被去除(圖7e)；VMD的去噪結(jié)果中SNR提升最為明顯，優(yōu)于其余四種方法，合成模型中的傾斜、正弦狀起伏地層以及斷層位置都得到了很好的保留，噪聲得到有效壓制，視覺表現(xiàn)上與原始無噪模型最為接近。對比局部相似度圖可知，VMD方法有效能量泄漏程度最輕微，保真度性能最好(圖7f)。

圖8給出了五種方法的第175道信號去噪結(jié)果對比，可以看出VMD去噪結(jié)果與無噪信號最為接近。因此，VMD方法在去除強噪聲的同時，可最大程度的減少有效地震能量的損失。

Curvelet、NLM、BM 3D、CEEMD和VMD五種方法處理耗時分別為0.724、11.291、181.039、2.207和2.513s，Curvelet去噪處理效率最高，BM3D處理效率最低。

圖7 Sigmoid模型五種方法去噪結(jié)果(左)及相應(yīng)的局部相似度(右)

圖8 Sigmoid模型五種方法第175道去噪結(jié)果對比

3.3 三維地震去噪實例

處理效率是考慮算法是否適用的重要指標(biāo)，在三維甚至五維地震資料的隨機噪聲壓制處理中，尤其需考慮計算耗時。本文將Curvelet 3D、NLM 3D、BM 4D和二維VMD(VMD 2D)四種方法運用于三維實際含噪地震數(shù)據(jù)處理。取某地隨機噪聲分布較多的三維疊后實際資料，數(shù)據(jù)尺度為221×271×752(依次表示x、y和時間三個方向的采樣點數(shù))，采樣間隔分別為20m和4ms。取實際數(shù)據(jù)x、y方向的第100個切片，如圖9所示。

圖10為四種方法去噪結(jié)果，圖11為相應(yīng)的噪聲剖面局部相似度圖，去噪結(jié)果與噪聲剖面的局部相似程度越高，表明有效能量泄漏越嚴(yán)重。對比各圖可知，Curvelet 3D方法去噪結(jié)果過于平滑，同相軸邊緣細節(jié)信息模糊，較難辨別，對應(yīng)局部相似度表現(xiàn)較多的有效能量損失；NLM 3D方法去噪結(jié)果中局部塊狀細節(jié)消失，可見部分噪聲依然分布于去噪結(jié)果中，對應(yīng)局部相似度出現(xiàn)較多的局部高異常區(qū)域；BM 4D方法的去噪效果有所改善，剖面視覺效果較好，局部相似度表現(xiàn)為有效信息損失少；本文二維VMD 2D方法的去噪結(jié)果最優(yōu)，剖面上同相軸細節(jié)、斷裂走向清晰，局部相似度呈低能均勻分布，有效同相軸信息丟失最少。

圖9 三維疊后含噪實際地震剖面

圖11 三維實際資料去噪后噪聲剖面的局部相似度

Sigmoid模型和實際資料的去噪結(jié)果數(shù)值統(tǒng)計如表1所示，可見Curvelet方法的SNR數(shù)值提升最小，但時間效率最高；CEEMD方法處理耗時最長；與CEEMD遞歸式模態(tài)分解方法不同，VMD在變分框架約束下一次性求解得到所有帶限BIMF分量，在SNR和SSIM指標(biāo)上均達最優(yōu)，去噪結(jié)果與原始無噪信號在細節(jié)上最為接近。由于在匹配追蹤優(yōu)化下快速定位到頻率中心，VMD方法時間效率高于NLM、CEEMD和BM 3D方法。三維實際地震資料的去噪結(jié)果LS同樣驗證其有效性。因此，VMD方法具備優(yōu)異的噪聲壓制、幅值保持性能的同時，還有具有較高的計算效率，可以滿足實時處理要求。

表1 去噪結(jié)果統(tǒng)計

4 結(jié)束語

本文在EMD信號分解理論基礎(chǔ)上，提出了在頻率域內(nèi)復(fù)數(shù)VMD方法。通過將信號分解過程轉(zhuǎn)移到變分框架內(nèi)，引入增廣Lagrange函數(shù)迭代求解各分量的頻率中心及帶寬參數(shù)，生成窄帶約束模態(tài)分量，使VMD方法在噪聲壓制和有效能量保持間尋求最優(yōu)平衡。實驗結(jié)果表明，在二維地震數(shù)據(jù)的隨機噪聲壓制處理方面，與CEEMD、Curvelet、NLM和BM 3D方法相比，VMD方法在SNR、SSIM和視覺效果方面均保持最優(yōu)趨勢，處理后同相軸細節(jié)變得清晰連續(xù)，微弱信號被有效保留；在三維地震數(shù)據(jù)的隨機噪聲壓制處理方面，與Curvelet 3D、NLM 3D和BM 4D方法相比，VMD方法具備優(yōu)異的噪聲壓制性能且計算效率高。