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

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

0 引言

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

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

1 基本原理

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

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

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

(1)

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

(2)

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

2 地震隨機(jī)噪聲壓制

2.1 地震信號(hào)的域表示

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

(3)

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

(4)

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

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

(5)

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

(6)

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

2.2 VMD框架求解

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

(7)

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

(8)

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

(9)

式中：L(·)表示目標(biāo)函數(shù)；“∶”是一個(gè)函數(shù)的標(biāo)記；〈·,·〉表示內(nèi)積；二次懲罰因子α是控制數(shù)據(jù)保真度的均衡參數(shù)，用于平衡變分正則項(xiàng)和二次約束項(xiàng)，在含噪聲情形時(shí)可保證信號(hào)重構(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í)行第一個(gè)內(nèi)循環(huán)更新uk

(10)

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

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

(11)

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

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

(12)

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

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

(13)

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

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

2.3 基于VMD的二維地震信號(hào)去噪方法

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

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

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

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

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

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

2.4 基于VMD的三維地震信號(hào)去噪方法

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

(14)

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

(15)

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

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

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

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

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

(5)對(duì)下一個(gè)時(shí)間窗口重復(fù)以上操作；

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

3 算例

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

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

(16)

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

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

圖1 各組分信號(hào)及VMD分解BIMF分量

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

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

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

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

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

3.2 二維地震去噪實(shí)例

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

(17)

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

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

(18)

(19)

式中：a和b為矩陣；A、B分別為a、b構(gòu)成的對(duì)角算子。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為像素動(dòng)態(tài)范圍，一般k1=0.01、k2=0.03、L=255；SSIM值變化范圍為[0，1]，值越大表示兩個(gè)圖像的結(jié)構(gòu)相似程度越高，去噪性能越強(qiáng)。

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

圖6 合成地震記錄Sigmoid模型

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

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

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

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

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

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

3.3 三維地震去噪實(shí)例

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

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

圖9 三維疊后含噪實(shí)際地震剖面

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

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

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

4 結(jié)束語(yǔ)

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