逯宇佳,曹俊興,田仁飛,呂雪松,何 沂

(成都理工大學地球物理學院,四川成都610059)

獨立分量分析(ICA)是20世紀末出現(xiàn)的一種新型信號處理方法,該方法可在信源未知的條件下,僅依賴混合信號將混疊在一起的不同信號分離開來[1]。其基本思想是假設原始信號為相互獨立的非高斯信號,利用信號的高階統(tǒng)計量信息,提取觀測信號獨立分量成分,最終得到相互獨立的源信號的近似估計[2-7]。1994年,COMON[8]首次提出了基于特征值分解的獨立分量分析方法。1998年,CICHOCKI等[9]利用獨立分量分析對含有加性噪聲的混合信號進行解混,取得較好的效果。在地球物理領域中,ICA主要應用于地震數(shù)據(jù)去噪處理。2003年,劉喜武等[10]假設地震記錄中有效信號和隨機干擾在統(tǒng)計上獨立且服從非高斯分布,對相鄰兩道地震記錄進行ICA處理實現(xiàn)地震資料的去噪;LU等[11]提出基于獨立分量分析的最大峰度自適應衰減算法實現(xiàn)了模型數(shù)據(jù)的多次波壓制;2006年,MIRKO[12]將獨立分量分析應用于P波和S波的波場分離;2011年,DONNO[13]將獨立分量分析與最小二乘法相結合提出了改進的多次波消除算法;呂文彪等[14-16]先后利用改進的ICA算法對疊后數(shù)據(jù)、地震屬性以及疊前數(shù)據(jù)進行了去噪處理;2012年,王維強等[17]將經(jīng)驗模態(tài)分解(EMD)技術與ICA技術結合應用于地震信號隨機噪聲壓制;2016年,孫成禹等[18]通過構造度量數(shù)據(jù)非高斯性的目標函數(shù),將地震數(shù)據(jù)轉(zhuǎn)換到ICA域,結合閾值法實現(xiàn)了復雜地震資料的去噪處理。

傳統(tǒng)的ICA算法如基于負熵的快速ICA算法只適用于地震同相軸較平緩的疊后數(shù)據(jù),且算法不夠穩(wěn)定,容易出現(xiàn)解混失敗現(xiàn)象,計算效率較低。針對這個問題,本文運用兩步奇異值分解法對傳統(tǒng)ICA算法的白化過程進行改進,將語音識別中的動態(tài)時間規(guī)整(DTW)與ICA算法相結合,利用動態(tài)時間規(guī)整進行地震道間的模式匹配,度量時間序列的相似程度,將地震數(shù)據(jù)的同相軸一一匹配,再利用ICA算法進行去噪處理。

1 方法原理

1.1 ICA算法

獨立分量分析處理需要對觀測信號矩陣進行去均值、預白化、ICA迭代處理。其中傳統(tǒng)的白化過程需要利用類似于主成分分析(PCA)方法,如零時滯的協(xié)方差矩陣的特征值分解來完成。利用兩步奇異值分解改進白化過程。

當觀察信號個數(shù)(m)大于源的個數(shù)(n)時,ICA算法中白化過程可利用兩步奇異值分解法(圖像處理中又稱AMUSE算法,即多個未知信號提取算法)[19]完成。這一過程能夠?qū)崿F(xiàn)源信號與噪聲信號的分離。主要實現(xiàn)步驟如下。

1) 觀測信號X減去它的均值,得到零均值觀測信號,實現(xiàn)X的中心化。

2) 估計觀測信號X的相關矩陣,即為協(xié)方差矩陣RX(0):

(1)

3) 利用奇異值分解(SVD)算法對RX(0)進行分解:

(2)

式中:VS為前n個特征值ΛS(從大到小排序)相對應的特征向量;VN為ΛN中后面m-n個不重要特征值ΛN=diag{λn+1,λn+2,…,λm}相對應的特征向量,λ為特征值。

5) 進行穩(wěn)健的預白化變換:

(3)

(4)

(5)

(6)

經(jīng)過以上兩步奇異值分解,即可完成解混過程,避免了傳統(tǒng)方法中復雜的迭代過程,在提高計算效率的同時也增加了算法的穩(wěn)定性。

1.2 動態(tài)時間規(guī)整(DTW)

設時間序列X=(x1,x2,…,xn),Y=(y1,y2,…,yn),則X,Y的DTW距離[20-22]DDTW(xi,yj)定義為[23]:

(7)

式中:Dbase(xi,yj)表示向量點xi和yj之間的基距離,可以根據(jù)情況選擇不同的距離度量。不失一般性,本文使用歐氏距離作為度量。根據(jù)公式(7),得到如圖1所示的一個鄰接矩陣,然后采用回溯法在鄰接矩陣上對DDTW(xi,yj)值進行遞歸回溯搜索DTW彎曲路徑。

圖1 DTW彎曲路徑示意

DTW距離允許序列點自我復制后再進行對齊匹配,能夠很好地支持時間軸彎曲,并且它可以對非等長時間序列進行度量,也支持時間軸伸縮。DTW距離實際上就是確定序列X與Y上每個點之間的對齊匹配關系(點對匹配)。

彎曲路徑必須滿足以下3個基本條件[24]。

1) 邊界條件:路徑起始于點(x1,y1)、終止于點(xn,yn),它表示兩個序列的起始點和結束點對應匹配。

2) 連續(xù)性:路徑上的任意兩個相鄰點(xi1,yj1)和(xi2,yj2)滿足條件0≤|i1-i2|≤1,0≤|j1-j2|≤1。

3) 單調(diào)性:若(xi1,yj1)和(xi2,yj2)為路徑上前后兩個點,則須滿足i2-i1≥0,j2-j1≥0。

滿足上述條件的彎曲路徑有很多,每條彎曲路徑都代表一種點對匹配關系。設彎曲路徑為F=(f1,f2,…,fk,…,fK),fk=(i,j)k是彎曲路徑上第k個元素,它表示xi和yj建立的匹配關系,路徑長度滿足max(n,m)≤K≤n+m-1。

點對匹配關系中,點對基距離之和的最小值即為DTW距離,對應的彎曲路徑為最佳路徑。DTW距離表示為:

(8)

求解最佳路徑需要構造一個m行n列的累積距離矩陣Mm×n,矩陣中的每個元素γi,j定義為:

γi,j=Dbase(xi,yj)+min{γi-1,j,γi-1,j-1,γi,j-1}

(9)

γi,j為序列X[1∶i]與序列Y[1∶j]的DTW距離,因此,DDTW(X,Y)=γm,n,γm,n可以用動態(tài)時間規(guī)整法求解[25]。

1.3 基于動態(tài)時間規(guī)整的ICA算法實現(xiàn)步驟

基于動態(tài)時間規(guī)整的ICA算法去噪過程主要分為以下3個步驟:

1) 選取第1道地震記錄作為模型道(模型道可以任選),利用動態(tài)時間規(guī)整算法對正演記錄的同相軸進行特征匹配(具體實現(xiàn)過程如圖1所示),在保證上覆水平層不變的情況下對下層傾斜地層同相軸進行校正拉平;

2) 利用兩步奇異值分解法進行去噪處理;

3) 利用動態(tài)時間規(guī)整提取的道間時差將拉平且去噪后的地震數(shù)據(jù)還原為原始形態(tài),完成整個去噪過程。

2 應用實例

2.1 模型測試

為了驗證基于動態(tài)時間規(guī)整的ICA算法在地震隨機噪聲壓制中的可行性,建立如圖2a所示的速度模型,模型包含水平界面、傾斜界面、斷層以及彎曲界面。通過正演模擬得到自激自收地震記錄如圖2b所示。在生成的正演記錄中加入信噪比為0.5的高斯隨機噪聲見圖3a。利用ICA算法直接(未使用DTW)對加噪數(shù)據(jù)進行去噪處理,圖3b為去噪后的地震剖面,圖3c為去除的噪聲。由圖3b可以看出,ICA算法在同相軸水平的位置去噪效果較好,而當同相軸傾斜或彎曲時,能夠明顯看到一部分有效波的信息被當成噪聲去除了(圖3c),有效信號損失嚴重,斷層位置處尤為明顯。

利用動態(tài)時間規(guī)整ICA算法對加噪數(shù)據(jù)進行去噪處理。圖4a為動態(tài)時間規(guī)整特征匹配后的拉平剖面,圖4b為兩步奇異值分解法對拉平剖面的去噪結果,圖5a為基于動態(tài)時間規(guī)整ICA算法去噪的結果剖面,與圖3b對比可以看出,經(jīng)過動態(tài)時間規(guī)整處理的去噪結果中隨機噪聲得到壓制,且在復雜界面如斷層、傾斜和彎曲界面的有效波信號保留較好,差剖面中沒有殘留有效信號(圖5b)。

2.2 疊前去噪應用

高信噪比的地震數(shù)據(jù)是進行地震解釋的重要基礎。解釋工作開始前,往往需要對疊前CRP道集數(shù)據(jù)進行層拉平或去噪處理,但這一工作非常費時,因此在保證去噪效果的同時,算法的計算效率顯得尤為重要。為了驗證兩步奇異值分解法的實際效果,本文選用疊前CRP道集進行驗證。疊前CRP道集是來自地下同一反射點的一系列地震響應特征,其每一道數(shù)據(jù)之間有極大的相似性,因此可利用ICA算法對其進行去噪處理。

圖2 速度模型(a)和正演記錄(b)

圖3 利用ICA算法處理的去噪效果(未使用DTW)a 加噪剖面;b 去噪后剖面;c 去除的噪聲

去噪效果如圖6所示。其中,圖6a為原始CRP道集,圖6b為基于負熵的快速ICA算法去噪結果,圖6c為兩步奇異值分解法的去噪剖面,圖6d為兩步奇異值分解法去除的噪聲。相比于原始記錄,兩種去噪方法都能夠有效去除原始剖面中的隨機噪聲,增強同相軸連續(xù)性。但是基于負熵的快速ICA算法處理結果中出現(xiàn)了由于迭代不收斂產(chǎn)生的壞道(圖6b紅色箭頭處),破壞了原有的有效波信息。而兩步奇異值分解法速度快,穩(wěn)定性高,不需要復雜的迭代過程,可以直接對矩陣進行分解,最終得到的剖面同相軸連續(xù),去噪效果有明顯的改善。圖7為圖6 中疊前CRP數(shù)據(jù)對應的疊加剖面。對比圖7a,圖7b和圖7c可以看出,利用兩步奇異值分解法對疊前資料進行去噪處理后所得疊后剖面同相軸連續(xù)性增強(紅色矩形框內(nèi)較為明顯),信噪比有所提高。

圖4 基于動態(tài)時間規(guī)整ICA算法去噪效果(Ⅰ)a 動態(tài)時間規(guī)整特征匹配拉平剖面;b 兩步奇異值分解法對拉平剖面去噪結果

圖5 基于動態(tài)時間規(guī)整ICA算法的去噪效果(Ⅱ)a 基于動態(tài)時間規(guī)整還原后的去噪結果剖面;b 去除的噪聲

2.3 疊后應用實例

將基于動態(tài)時間規(guī)整ICA算法應用于疊后地震資料隨機噪聲壓制,結果如圖8和圖9所示。其中,圖8a 是處理前的原始地震剖面;圖8b是經(jīng)過DTW處理后的層拉平剖面,由于原始地震資料的品質(zhì)限制,中心位置拉平效果一般;圖8c是拉平后兩步奇異值算法去噪結果。圖9b是基于負熵的快速ICA算法去噪結果,圖9c是基于動態(tài)時間規(guī)整的ICA算法去噪結果。對比圖9a,圖9b和圖9c可以看出,針對地層傾角較大且有斷層存在的復雜剖面,基于負熵的快速ICA算法雖然能去除一部分的隨機噪聲,但是處理后的有效波同相軸能量變?nèi)?連續(xù)性反而變差,這是因為在處理過中相鄰地震道存在的時差使得地震道之間的相似性降低,從而導致算法不收斂產(chǎn)生壞道。采用基于動態(tài)時間規(guī)整的ICA算法去噪后,同相軸連續(xù)性明顯增強,弱振幅信號得到放大,整體能量更為均一,斷層信息更為突出,斷點清晰可見(圖9c紅色橢圓標注)。圖10對比了基于動態(tài)時間規(guī)整ICA算法去噪前、后的頻譜,可以看出,本文方法去噪前、后頻譜基本一致,未改變原始資料的頻譜特征。相比于傳統(tǒng)的方法,基于動態(tài)時間規(guī)整的ICA算法更適合復雜地震資料的隨機噪聲壓制,有較好的實用價值。

圖6 疊前CRP道集去噪效果a 原始記錄;b 基于負熵的快速ICA算法去噪結果;c 兩步奇異值分解法去噪結果;d 兩步奇異值分解法去除的噪聲

圖9 疊后剖面去噪效果a 原始剖面;b 基于負熵的快速ICA算法去噪結果;c 基于動態(tài)時間規(guī)整ICA算法去噪結果;d 基于動態(tài)時間規(guī)整ICA算法去除的噪聲

圖10 基于動態(tài)時間規(guī)整ICA算法去噪前、后的頻譜對比

3 結論

本文將ICA算法與動態(tài)時間規(guī)整相結合,提出了基于動態(tài)時間規(guī)整ICA算法,實現(xiàn)了疊前疊后數(shù)據(jù)的隨機噪聲壓制處理,獲得以下結論和認識:

1) 利用兩步奇異值分解法對矩陣進行穩(wěn)健白化,可將有效波和噪聲分離。提高了ICA算法的穩(wěn)定性和計算效率,隨機噪聲壓制效果明顯;

2) 將動態(tài)時間規(guī)整與ICA算法相結合,解決了僅利用ICA算法對非平緩同相軸的地震資料去噪不適應性問題;

3) 基于動態(tài)時間規(guī)整的ICA算法在實際資料應用中取得較好效果。既可以應用于疊前資料的隨機噪聲壓制,也可以應用于疊后地震資料,具有普適性。