基于ICEEMD-RAKSVD的地震弱信號去噪方法

朱劍兵

(中國石化 勝利油田分公司物探研究院，東營 257022)

0 引言

眾所周知，地震勘探由于噪音來源多，一些弱信號往往會淹沒在噪音里，常規方法很難去識別弱信號中的構造、儲層信息。為了有效識別弱信號，前人引入了很多方法，其中包括：基于非線性動力系統的混沌理論[1]、基于獨立分量分析的盲源分離技術[2]、基于奇異值分解(SVD)的方法[3]、基于經驗模式分解(EMD)[4]和基于高階統計量[5]等識別技術。

Aharon[6]提出了KSVD算法；韓文功等[7]提出的SVD分解方法可以較好進行弱信號的識別和去噪。KSVD算法由于需要進行多次SVD分解的操作,因此，耗時過長, Rubinstein[8]在KSVD的基礎上提出了AKSVD算法(Approximate KSVD)，但是其效果比較差；Dumitrescud等[9]在KSVD基礎上引入了正則化，得到RKSVD與RAKSVD(正則化近似KSVD)方法,并進行了對比，RAKSVD在效果幾乎與KSVD一致的情況下，大大提高了效率。在模態分解去噪方面，Huang等[10]在1998年首次提出了經驗模態分解理論(EMD)。然而EMD也并非完美，仍存有一定的缺陷性，諸如部分模態之間的混疊現象、端點效應以及篩分準則等；Wu等[11]引入了總體經驗模態分解(EEMD)，對EMD算法進行了的改進，但又出現了新的問題，恢復信號中仍含有一定程度的殘余噪音成分，無法將其消除干凈；Torres[12]對EEMD算法進行了改進，引入的完備總體經驗模態分解(CEEMD)擁有高精度去噪的優勢，但也不能盲目去除高頻模態，其中還保存著很多有效信息；楊濤[13]結合KSVD和CEEMD的特點，利用KSVD對弱信號的識別優勢及CEEMD的高效去噪，設計了一個基于KSVD的CEEMD去噪識別方法，去噪效果較好，但計算效率相對較低。筆者將改進的CEEMD(ICEEMD)與正則化近似KSVD(RAKSVD)有機結合起來，形成了基于ICEEMD-RAKVSD的去噪方法，該方法能夠對地震弱信號有效識別，取得了較好的弱信號去噪效果，同時計算效率也得到了提升。

1 RAKSVD方法原理

由于KSVD中每次原子更新都要進行一次SVD分解，從而導致了計算效率較低。Rubinstein等[8]提出了一種近似KSVD算法(AKSVD)，該算法能夠明顯提高計算效率，但由于引入了近似，SVD分解結果存在一定的誤差，因而影響了分解效果。針對這一問題，Irofti等[14]提出了正則化近似KSVD(RAKSVD)算法，建立嶺估計模型，解決稀疏分解中的病態問題，從而使SVD分解結果的誤差得到很大程度上的控制，提高了稀疏分解效果。

1.1 KSVD字典

KSVD字典是學習型字典之一，具有更適合反映數據結構特點的優勢，能獲得較好的去噪效果。KSVD算法通過對初始化字典進行逐列更新，用來尋找一個最佳字典的過程，假定訓練的信號為Y∈Rm×n，其中m為信號的行數，n為信號的列數，每列對應于一個訓練樣本。學習型字典模型可以描述為式(1)。

(1)

式中：D∈Rm×p為待學習字典，p為字典的列數；X∈Rp×n為稀疏表示系數；‖·‖0表示l0范數；‖·‖2表示l2范數；T為每列個數的最大值。

首先固定字典D，利用正交匹配追蹤算法(OMP)對已知樣本數據進行分解，獲取初始稀疏分解系數矩陣。

字典更新實施步驟如下：

1)記當前更新的第k列字典原子為dk，X中使用到dk原子的索引集合為Ik。

2)計算殘差：

(2)

式中：xj為x的第j列；Ek為去除dk外的殘差。

(3)

式中：U表示左奇異矩陣；V表示右奇異矩陣；Λ是奇異值合集；r表示非零奇異值的數量；ui表示左奇異矩陣的第i行；σi表示第i個奇異值；νi表示右奇異矩陣的第i列。

4)原子和稀疏系數更新為：

xk=σ1ν1，dk=u1

(4)

1.2 近似KSVD

Rubinstein等[8]在KSVD的基礎上提出了AKSVD算法,主要是在上述步驟3)中采用近似計算更新原子,AKSVD算法更新后的原子和相關系數分別為式(5)和式(6)。

(5)

(6)

1.3 正則化AKSVD

稀疏分解是一種不適定問題,而正則化方法能有效提高分解效果。正則化字典習模型可描述為式(7)。

(7)

式中，μ>0為正則化參數。同樣的，記當前更新的字典原子為dk，X中使用到原子的索引集合為Ik，則公式(2)可以變為式(8)。

μ‖Xk,Ik‖2

(8)

由于相關系數里除了第k行外，都為固定項。

由公式(2)則可得到：

(9)

(10)

根據式(10)，從KSVD算法的正則化[9]可以得出，正則化約束后的更新原子方法和原算法一致，而正則化參數是應用到了稀疏系數的更新，使系數變為1/(1+μ)倍，通過給稀疏系數一個權值，每次通過這個權值影響后續的原子更新，最終完成字典的更新。所以RAKSVD的原子和稀疏系數更新[14]為式(11)和式(12)。

(11)

(12)

1.4 效果對比

圖1為KSVD、AKSVD、RAKSVD三種訓練字典的誤差對比圖，其中訓練集Y的稀疏度為5，字典維數為64*256，正則化參數μ為0.1，整個分解過程重復迭代120次，取平均值得到均方根誤差RMSE(單位無量綱)隨迭代次數的關系。

圖1 誤差對比圖Fig.1 Error comparison

由圖1可以看出，這三種KSVD訓練字典在迭代20次后RMSE基本趨于穩定，其中AKSVD訓練效果里相對最差，RAKSVD與KSVD訓練效果接近，經前人研究AKSVD相比于KSVD效率提高,故而RAKSVD在提高效率的同時，也能達到較好的分解效果。

2 ICEEMD方法原理

首先，EEMD算法在解決EMD存在的模態混疊、端點效應等缺陷的同時，也導致了分解的IMF分量存在噪音殘余，而CEEMD通過在不同階段加入不同高斯白噪音進行調整，比例大小根據待分解信號自適應計算得到，同時每個符合條件的模態分量對應的計算殘差值是固定唯一的。不僅解決了模態混疊問題，同時使原始信號復原的精度變得更高。

按照Torres[12]所設計的CEEMD方法步驟對信號做處理，原始信號中各個頻率成分的信號，基本上都能夠得到較好的分解，但該算法仍存在可能出現虛假模態的問題，同時由于高斯白噪聲的加入，分解的模態中會混入一定程度的噪聲。為了解決上述問題，改進的CEEMD(Improved CEEMD,ICEEMD)算法重新定義了各個模態局部均值的求取。

假設：Ej(·)為CEEEMD分解得到的第j階模態,s(t)為原始地震信號，ωi為加入的第i個高斯白噪聲，白噪聲加入的比例為εk。則改進后的ICEEMD實現步驟如下：

1)將不同正負成對的白噪聲分別加入原始信號中，合成得到T個新信號(T為偶數)，計算集合平均得到IMF1(t)。

(13)

2)令r0(t)=s(t)，對k=1、…、K，計算k階殘差rk(t)為式(14)。

rk(t)=rk-1(t)-IMFk(t)

(14)

3)對k階殘差添加不同比例系數的白噪聲：rk(t)+εkEk+1(ωi(t))，經EMD分解得到IMF1(t)，計算總體平均并將其作為IMF(k+1)(t)。

(15)

4)重復上述步驟，直到達到迭代條件，得到最終的殘差E(t)為式(16)。

(16)

CEEMD去噪往往是直接去除IMF1或IMF1和IMF2分量，但其中往往包含有效信息，所以，雖能有效去除高頻噪音，但對于弱信號去噪的問題，容易丟掉真正的有效信號，使同相軸不連續。

相比于CEEMD算法，改進后的CEEMD(ICEEMD)算法可以加入更少的白噪聲個數，就能使信號得到很好的模態分離，同時，在相同參數的條件下，ICEEMD分解的效果更好，計算效率更高。

3 去噪方法

KSVD字典稀疏去噪方法是屬于對壓縮感知原理的改進算法，通過訓練樣本求取冗余字典，較好地刻畫出信號特征，從而更加精確地重構原信號，實現噪聲去除和保留有效信號的目標。KSVD字典雖然能將噪音基本壓制和消除，但在信號的波峰波谷處仍會存有殘留一些毛刺。ICEEMD算法不需要提前設置基函數，它是以信號本身的特點作為基礎來自動分解和重構信號，打破了測不準原理的限制，達到了高精度識別弱信號的能力。

具體實現步驟為：

1)利用ICEEMD算法對地震信號進行分解得到固有模態函數集，根據分解效果，將其劃分為噪聲主導模態、信號主導模態和過渡模態。

2)對于噪聲主導模態，則直接舍棄。

3)對于信號主導模態，將其與原始信號疊加，形成主導信號集，通過RAKSVD方法進行主導信號集的分解和去噪。

4)對于過渡模態，先將其疊加，再進行ICEEMD分解；舍棄噪音模態后，與原始信號疊加，重構過渡信號集，通過RAKSVD方法進行過渡信號集的分解和去噪。

5)將步驟3)和步驟4)中的兩個分解結果重構并平均，得到最終的去噪結果。

4 模型測試

含有弱信號模型如圖2(a)所示，模型共有25道，每道300個采樣點，Ricker子波主頻25 Hz。在200 ms處有一弱信號。圖2(b)為含噪模型，加入噪聲后的峰值信噪比(PSNR)為1.2，圖2中強軸能夠明顯看到，但弱軸在噪聲的埋沒中已經不明顯了。對圖2(b)分別采用五種算法進行去噪處理，得到的結果如圖3所示，圖4為模型的15道的單道方法去噪結果，表1是各方法的去噪效果與時間對比。

表1 五種方法針對模型數據的去噪效果對比Tab.1 Comparison of the effects of five methods for denoising the model data

圖3 五種方法針對模型數據的去噪結果Fig.3 The result of five methods for denoising the model data(a)KSVD去噪；(b)RAKSVD去噪;(c)CEEMD去噪；(d)CEEMD-KSVD去噪；(e)ICEEMD-RAKSVD去噪

圖4 單道去噪方法結果對比Fig.4 Comparison of single channel denoising results

從圖3、圖4和表1可以看到，如果CEEMD單獨使用，則對弱信號的識別和去噪效果最差；KSVD和RAKSVD對模型中弱信號的識別和去噪效果較好，說明該類方法相比于CEEMD，對弱信號識別有著獨特的優勢，但仍有一些毛刺的殘留，在實際應用中，當地下介質不均勻或由薄互層形成復波時，其缺點會被放大；CEEMD-KSVD和ICEEMD-RAKSVD是將模態分解和字典訓練聯合在一起進行去噪，將模態分解的信號分解為主導模態和過渡模態，利用KSVD和RAKSVD對弱信號識別的優勢，最終使含噪弱信號中的噪聲得到了較好的壓制，且保留了有效弱信號，ICEEMD-RAKSVD方法能夠達到最好的去噪效果。由表1可以看到，由于CEEMD算法需要多次迭代，每次迭代需要加入大量成對的高斯白噪聲，因此，其計算效率低。KSVD和RAKSVD計算時間少，RAKSVD較KSVD效率提高了60%。由于CEEMD的引入，CEEMD-KSVD和ICEEMD-RAKSVD算法明顯增加了計算時間，ICEEMD-RAKSVD 相對于CEEMD-KSVD效率提高6%，對于簡單的稀疏模型，CEEMD的計算時間占主導，因此，計算效率提升幅度不是很大。

5 實際資料處理

針對實際資料，對本方法的去噪效果進行了進一步測試。圖5(a)為原始局部疊前地震道集，圖5中隨機噪聲十分發育，信噪比非常低，其中弱同相軸不清晰、不連續。上述五種方法的去噪結果如圖5(b)～圖5(f)所示，圖6為單道去噪結果對比展示，圖7為去噪后殘差剖面效果對比，去噪效果與時間對比見表2。

表2 五種方法針對實際資料的去噪效果對比Tab.2 Comparison of the effects of five methods for denoising the actual data

圖5 五種方法針對實際資料的去噪結果Fig.5 The result of five methods for denoising the actual data(a)原始疊前地震道集；(b)KSVD去噪；(c)RAKSVD去噪；(d)CEEMD去噪；(e)CEEMD-KSVD去噪；(f)ICEEMD-RAKSVD去噪

圖6 單道去噪方法結果對比Fig.6 Comparison of single channel denoising results

從圖5和表2可以看出，對于弱地震信號，實際資料的處理結果與模型測試得到的結論相一致：CEEMD去噪效果較差，去噪后PSNR最低，將部分弱同相軸的有效信號充當噪音去除，弱同相軸損失嚴重；RAKSVD與KSVD稀疏去噪較明顯，但仍存留了許多噪音。CEEMD-KSVD和ICEEMD-RAKSVD算法都能夠很好地壓制噪聲，紅方框標記可以看出ICEEMD-RAKSVD有效弱信號得以保留。從圖5(a)可以更為清晰看出五種方法對比，而圖6可以得出 CEEMD-KSVD和ICEEMD-RAKSVD算法去噪效果幾乎一致，能有效去除噪音，圖7藍圓圈中可以看出，CEEMD-KSVD去除的殘余還仍然保留些有效信號，而本文方法去除的噪音沒有保留有效信號，去噪效果更好，同時由于實際地震數據非常復雜，訓練字典的時間大大增加，從表2可以看出，RAKSVD較KSVD計算時間縮短了65%，下降幅度較為明顯，ICEEMD-RAKSVD相對于CEEMD-KSVD時間縮短了40%。對實際地震資料的去噪處理，更能突顯本文方法在提高計算效率方面的優勢。

圖7 去噪后殘差剖面效果對比Fig.7 Comparison of residual profile effects after denoising(a)CEEMD_KSVD;(b)ICEEMD_RAKSVD

6 總結

由于弱地震信息往往受噪音影響更大，通過去噪能更好地識別弱信號。CEEMD去噪方法直接舍棄高頻模態，容易造成高頻中有效信息的丟失，而完整地保留了其他含有噪音模態，KSVD去噪方法和RAKSVD去噪方法，對于弱信號的識別效果較好，但去噪后殘留毛刺，其中RAKSVD計算效率得到明顯提高，CEEMD-KSVD去噪方法能達到較好的去噪效果，但計算效率相對較低，這里指出的基于ICEEMD-RAKSVD去噪方法能在保證較好的去噪效果情況下，還能有效提升了計算效率，在實際地震處理中提升效率尤為重要，該方法具有一定的應用價值。