基于變分模態分解的地震隨機噪聲壓制方法

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

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

0 引言

地震信號隨機噪聲衰減是地震數據處理領域的熱點和難點問題。從觀測數據中有效去除隨機噪聲干擾、提高信號信噪比和分辨率，是正、反演計算和地質解釋的前提。目前，學者們提出了大量的地震信號隨機噪聲壓制方法，比如基于濾波理論的方法[1-2]、基于小波域變換方法[3-6]、基于矩陣理論的變換方法[7-10]、基于信號分解理論方法[11-16]等。近年來，由于在處理非平穩及非線性數據上具有很高的信噪比，基于信號分解的經驗模態分解(Empirical Mode Decomposition，EMD)方法[17]已成為信號分解領域研究的熱點。

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

1 基本原理

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

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

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

(1)

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

(2)

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

2 地震隨機噪聲壓制

2.1 地震信號的域表示

線性地震數據可表示為

(3)

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

(4)

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

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

(5)

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

(6)

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

2.2 VMD框架求解

基于上述地震數據域分析，在頻率域中估計BIMF分量頻率帶寬目標函數[20]

(7)

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

(8)

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

(9)

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

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

(2)n=n+1，執行主循環；

(3)Fork=1∶K-1，執行第一個內循環更新uk

(10)

(4)k=K，結束第一個內循環；

(5)Fork=1∶K-1，執行第二個內循環更新ωk

(11)

(6)k=K，結束第二個內循環；

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

(12)

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

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

(13)

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

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

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

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

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

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

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

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

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

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

在三維地震數據中，平面波可表示為

(14)

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

(15)

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

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

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

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

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

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

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

3 算例

3.1 VMD參數優選

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

(16)

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

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

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

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

圖2 VMD分解各模態頻譜分布

圖3 VMD分解各模態頻率中心變化

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

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

3.2 二維地震去噪實例

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

(17)

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

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

(18)

(19)

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

結構相似度定義為

(20)

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

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

圖6 合成地震記錄Sigmoid模型

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

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

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

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

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

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

3.3 三維地震去噪實例

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

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

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

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

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

表1 去噪結果統計

4 結束語

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