應用結(jié)構(gòu)聚類字典學習壓制地震數(shù)據(jù)隨機噪聲

張 巖 任偉建 唐國維

(①東北石油大學計算機與信息技術(shù)學院，黑龍江大慶 163318； ②東北石油大學電氣信息工程學院，黑龍江大慶 163318； ③東北石油大學現(xiàn)代教育技術(shù)中心，黑龍江大慶 163318)

1 引言

隨機噪聲壓制一直是地震資料處理中研究的熱點與難點。當前地震數(shù)據(jù)隨機噪聲壓制技術(shù)的思路是首先將地震數(shù)據(jù)稀疏表示，利用地震數(shù)據(jù)與噪聲信號在稀疏域內(nèi)的不同特征分離噪聲，要求所采用的稀疏表示基函數(shù)能夠準確地捕捉到地震波前信息，用較少的系數(shù)就可以描述這些主要特征。

地震數(shù)據(jù)稀疏表示的方法有正交變換法、多尺度幾何分析法、字典學習法等。正交變換法包括Fourier變換[1]、 Radon變換[2]、Wavelet變換[3]等；多尺度幾何分析法包括Ridgelet變換[4]、Curvelet變換[5]、Seislet變換[6]等。考慮到實際地震數(shù)據(jù)通常由多種元素組成，單一的固定變換基函數(shù)難以獲得最優(yōu)的稀疏表示效果。基于字典學習法的稀疏表示能根據(jù)地震數(shù)據(jù)本身的特點，通過對數(shù)據(jù)的學習和訓練，自適應地調(diào)整變換基函數(shù)以適應特定數(shù)據(jù)本身，因此具有更好的稀疏表示能力。近年來，基于學習型字典稀疏表示的地震數(shù)據(jù)去噪技術(shù)在業(yè)內(nèi)得到了廣泛關注。Beckouche等[7]提出一種根據(jù)包含噪聲地震數(shù)據(jù)自適應學習得到稀疏字典的去噪方法；Chen等[8]結(jié)合自適應學習字典和固定基函數(shù)的兩種稀疏表示方法，提出雙稀疏字典，能更好地處理去噪問題；Chen[9]以K-SVD法[10]為基礎，提出多維地震數(shù)據(jù)的快速字典學習稀疏表示去噪方法；Turquais等[11]提出了基于一致性約束的字典學習去噪方法。

上述學習型字典稀疏表示方法的基本思路是：首先將地震數(shù)據(jù)分塊，每一塊包含多個地震記錄道在一定的采樣時間段內(nèi)波形的信息，以地震數(shù)據(jù)塊為訓練樣本，利用超完備字典學習技術(shù)，根據(jù)地震數(shù)據(jù)本身的特點，自適應構(gòu)造超完備字典，稀疏表示地震數(shù)據(jù)，從而恢復數(shù)據(jù)的主要特征[12，13]。此類方法的共同點在于利用一個通用的字典稀疏表示整個地震數(shù)據(jù)，考慮到地震數(shù)據(jù)中不同空間位置的波形變化差異較大，通用字典不足以對整個數(shù)據(jù)的每一個局部內(nèi)容都能最優(yōu)地稀疏表示。在圖像處理領域中，為提高稀疏表示對整個圖像中局部結(jié)構(gòu)的適應性，Dong等[14，15]、Xu等[16]提出自適應局部字典學習方法，利用自然圖像的局部相似性，自適應得到局部字典，能夠更好地刻畫圖像局部結(jié)構(gòu)，獲得更好的重構(gòu)圖像效果。地震數(shù)據(jù)中隨機噪聲普遍滿足平穩(wěn)性、高斯性，噪聲壓制算法通常可采用高斯隨機噪聲進行模擬仿真，更重要的是地震數(shù)據(jù)中同一反射波在相鄰地震道上到達時間、振幅相近，波形相似。為利用局部數(shù)據(jù)道中同相軸振幅與波形的相似性，借鑒文獻[14-16]的思想，本文提出基于結(jié)構(gòu)聚類的局部超完備字典學習方法，壓制地震數(shù)據(jù)中的隨機噪聲。

2 基于結(jié)構(gòu)聚類的字典學習

2.1 全局字典稀疏表示分析

地震數(shù)據(jù)中相同顏色的分塊具有相似的結(jié)構(gòu)(圖1)，并且在地震數(shù)據(jù)中重復多次出現(xiàn)，說明地震數(shù)據(jù)塊結(jié)構(gòu)具有較強的自相似性。全局字典稀疏表示的一般模型均假設信號的稀疏表示系數(shù)之間是獨立同分布。對圖1地震數(shù)據(jù)利用K-SVD全局字典稀疏表示后，得到字典如圖2所示，字典中原子大小為8×8，原子的數(shù)目為256。圖3為圖1數(shù)據(jù)對應稀疏域的系數(shù)分布結(jié)構(gòu)，可見其系數(shù)分布不是完全隨機的，而是具有一定的規(guī)律性與冗余性，以此作為先驗知識應用到地震數(shù)據(jù)的隨機噪聲壓制處理中，可以對地震數(shù)據(jù)進行更穩(wěn)定、更稀疏的描述。

圖1 地震數(shù)據(jù)塊相似性(相同顏色子塊具有相似的結(jié)構(gòu))

圖2 K-SVD學習字典

圖3 K-SVD學習字典對應的系數(shù)分布白色像素代表幅值較大的系數(shù)，黑色代表幅值較小的系數(shù)

2.2 局部字典學習模型

借鑒文獻[15]的思路，利用地震數(shù)據(jù)塊結(jié)構(gòu)的自相似特點，提出地震數(shù)據(jù)結(jié)構(gòu)聚類字典學習模型的目標函數(shù)為

(1)

式中：μk表示第k個聚類Ck的質(zhì)心；Φk為第k個聚類的字典；y為無噪地震數(shù)據(jù)x的含噪結(jié)果；D為地震數(shù)據(jù)y的整體字典；α為地震數(shù)據(jù)在字典D表示下對應的系數(shù)；K為聚類的類別數(shù)目；λ1和λ2為正則化參數(shù)。稀疏系數(shù)α具有自身結(jié)構(gòu)冗余的特點，對聚類質(zhì)心進行更新與調(diào)整，并用各個分類的聚類質(zhì)心重新編碼，以得到更稀疏的系數(shù)，相當于對原始地震數(shù)據(jù)更稀疏的表示和描述。對各個聚類質(zhì)心也可以采用與原始數(shù)據(jù)一樣的字典進行表示，即μk=Φkβk，其中βk代表μk在字典Φk下的系數(shù)。受L1范數(shù)下的壓縮感知理論啟發(fā)，可用L1范數(shù)代替L2范數(shù)，將式(1)優(yōu)化為

(2)

上式利用地震數(shù)據(jù)的自相似性與稀疏系數(shù)結(jié)構(gòu)冗余性，將字典學習和結(jié)構(gòu)分類統(tǒng)一起來。βk可以看成是通過結(jié)構(gòu)分類所獲得聚類中心而得到的新編碼，依照先驗知識不斷促進系數(shù)的稀疏程度。

2.3 結(jié)構(gòu)聚類方法

地震數(shù)據(jù)塊結(jié)構(gòu)聚類的過程無需關心某一類的特征，僅需將相似的地震數(shù)據(jù)塊聚為一類，因此不需要對訓練數(shù)據(jù)進行學習，即為無監(jiān)督學習的過程。典型的聚類算法主要包括：屬于劃分法的K均值聚類算法[17]、K中心點算法[18]、Clarans算法[19]；屬于層次法的Birch算法[20]、Cure算法[21]、Chameleon算法[22]等；基于簇密度的Dbscan算法[23]、Optics算法[24]、Denclue算法[25]等；基于網(wǎng)格的Sting算法[26]、Clique算法[27]、Wave-Cluster算法[28]等。由于K均值聚類算法計算效率高，在對大規(guī)模數(shù)據(jù)聚類時被廣泛應用。利用K均值聚類算法對地震數(shù)據(jù)分塊結(jié)構(gòu)聚類的關鍵環(huán)節(jié)有3個。①初始質(zhì)心的選取：適當?shù)某跏假|(zhì)心是基本K均值聚類算法的關鍵。為提高計算效率，本文利用隨機選擇地震數(shù)據(jù)塊的方式作為初始質(zhì)心。每次運行過程中,在地震數(shù)據(jù)的分塊中使用一組不同的隨機初始質(zhì)心，然后選取具有最小誤差平方和的簇集作為質(zhì)心。該策略計算效率高，效果相對好。②度量距離的選擇：聚類采用的距離度量方法(相似性度量方法)有很多種，常用的有歐氏距離、閔可夫斯基距離、曼哈頓距離、皮爾遜相關系數(shù)、余弦相似度等。余弦相似度更多地是從方向上區(qū)分差異，對絕對的數(shù)值不敏感。本文重點關注地震數(shù)據(jù)各個分塊的結(jié)構(gòu)相似情況，因此采用余弦相似度作為距離的度量。③質(zhì)心的計算：在K均值算法中，將簇中所有地震數(shù)據(jù)樣本塊的均值作為該簇的質(zhì)心，即向量各維取平均作為新的質(zhì)心。

2.4 字典學習方法

文獻[10]中字典學習采用主成分(PCA)分解提取主成分，雖然在構(gòu)造的字典中各個原子滿足正交條件，計算效率相對高，但是直接將該算法應用于地震數(shù)據(jù)的去噪存在一定的局限性，原因在于地震數(shù)據(jù)中相當一部分中、高頻信息比較重要，去噪過程應該注意保留，同時對地震數(shù)據(jù)PCA分解需要計算協(xié)方差矩陣，對于某些矩陣求協(xié)方差時很可能會損失精度，更重要的是冗余度較大的超完備字典更有利于地震數(shù)據(jù)細節(jié)的保留，因此本文對每一類數(shù)據(jù)塊采用奇異值分解構(gòu)造超完備字典，稀疏編碼該類地震數(shù)據(jù)。

2.5 模型求解算法

針對式(2)目標函數(shù)的最大后驗概率估計問題為

(3)

根據(jù)貝葉斯公式將式(3)改寫為

(4)

式中：前一項對應于原目標函數(shù)中的誤差項；后一項為先驗概率項。由于觀測模型采用的是加性高斯噪聲模型，并且用獨立同分布的拉普拉斯分布函數(shù)描述α，最終可得

(5)

式中：σ1和σ2分別為式(4)中兩項的標準差；σn為噪聲標準差。則有

(6)

從上式可以發(fā)現(xiàn)，實際上對各個稀疏系數(shù)局部方差的估計，可以從輸入數(shù)據(jù)自適應得到正則項系數(shù)的估計值。通過估計每一類系數(shù)的正則化參數(shù)，獲得實際地震數(shù)據(jù)的逼近，然后再利用逼近的地震數(shù)據(jù)更新正則化參數(shù)的估計，依此不斷迭代，最終收斂獲得最優(yōu)的地震數(shù)據(jù)逼近值。在字典確定的情況下，對于上述目標函數(shù)需要求解模型中雙L1范數(shù)的優(yōu)化問題，即需要聯(lián)合求解α和β，然而這是一個非凸的優(yōu)化問題，無法通過線性規(guī)劃方法求解，可采用交替迭代的優(yōu)化方法求解，即先固定待求量中的一個，對另一個待求量進行優(yōu)化； 再固定己經(jīng)優(yōu)化了的待求量轉(zhuǎn)而對另一個待求量進行優(yōu)化，兩步交替迭代進行，直到目標函數(shù)收斂。在稀疏系數(shù)α的聚類中心β確定的情況下，本文采用雙變量迭代閾值算法，交替更新α和β。借鑒替代函數(shù)的思想[29]，更新α

(7)

式中：Sτ1,τ2為迭代收縮算子；αj,t和βj,t分別代表向量αj和向量βj中的任意一個元素；

(8)

(9)

其中c1與c2是兩個給定的常數(shù)。

受文獻[32]的啟發(fā)， 通過下式更新恢復地震數(shù)據(jù)

(10)

其中

(11)

上式是實現(xiàn)迭代正則思想的算子，ω是一個控制迭代反饋的噪聲級小正數(shù)。

本文提出的基于結(jié)構(gòu)聚類超完備字典地震數(shù)據(jù)去噪算法，通過K均值聚類算法和奇異值分解得到字典Φ與稀疏系數(shù)α；然后迭代更新正則化參數(shù)，更新質(zhì)心估計β和地震數(shù)據(jù)估計值，最后輸出去噪后的地震數(shù)據(jù)，具體算法流程如圖4所示。

圖4 本文算法流程

3 實驗結(jié)果與分析

實驗過程中，設置每一個地震數(shù)據(jù)塊大小為8×8，相鄰數(shù)據(jù)塊之間的間隔設為2，初始化采用離散余弦基構(gòu)造超完備字典，結(jié)構(gòu)聚類數(shù)目為40，字典學習次數(shù)為5，結(jié)構(gòu)聚類次數(shù)為5。地震數(shù)據(jù)去噪效果的衡量指標采用峰值信噪比(RSNR)，其表達式為

(12)

式中RMSE為不含噪聲原始數(shù)據(jù)與去噪后數(shù)據(jù)的均方誤差，定義為

(13)

式中：N為時間樣點數(shù)；M為地震記錄道數(shù)。為了模擬地震數(shù)據(jù)中包含的隨機噪聲，本文在原始地震數(shù)據(jù)中加入零均值正態(tài)分布的高斯隨機噪聲，噪聲的標準差與原始地震數(shù)據(jù)的標準差呈正相關。噪聲標準差定義為

(14)

3.1 合成地震數(shù)據(jù)實驗分析

圖5a是合成原始單炮地震記錄，道數(shù)為140，樣點數(shù)為287；圖5b為加入噪聲強度比例因子為0.03的高斯隨機噪聲單炮記錄，其RSNR為30.42dB。圖6a為本文算法去噪后結(jié)果，RSNR為39.53dB，約提高了9dB，且地震數(shù)據(jù)的同向軸方向性與連續(xù)性保持較好；圖6b為K-SVD算法的去噪結(jié)果，RSNR為35.99dB，可見本文算法結(jié)果具有較高的RSNR，而且壓制隨機噪聲的效果更好。圖6c是本文算法去噪結(jié)果與原始數(shù)據(jù)差值，圖6d為本文算法結(jié)果的RSNR隨迭代次數(shù)的變化曲線，說明本文算法具有較高的穩(wěn)定性和有效性。

圖5 合成原始單炮記錄(a)及其加噪記錄(b)

圖6 合成地震數(shù)據(jù)實驗(a)本文算法隨機噪聲壓制結(jié)果; (b)K-SVD算法隨機噪聲壓制結(jié)果; (c)本文算法去噪結(jié)果與原始數(shù)據(jù)的差值; (d)本文算法結(jié)果的RSNR隨迭代次數(shù)的變化曲線

圖7a是5次聚類過程中質(zhì)心變化的圖示，按由上至下的順序分別給出聚類1～5次的質(zhì)心塊集合，由于聚類數(shù)目為40，因此每一行包括40個質(zhì)心塊，隨著聚類次數(shù)的增加，特征類似的質(zhì)心塊逐漸減少，不同類別質(zhì)心之間的差異逐漸明顯，表明聚類算法設計合理有效。圖7b為40個分類的地震數(shù)據(jù)塊，經(jīng)過5次迭代后訓練得到的字典，字典中的每列為一類數(shù)據(jù)塊的局部字典，每類的字典包括64個原子，每個原子大小為8×8，可以看出每一類數(shù)據(jù)塊對應的局部字典中原子存在一定的冗余度，保證了對地震數(shù)據(jù)更完備的描述。

圖7 聚類質(zhì)心變化(a)和超完備局部字典(b)

為分析聚類數(shù)目對算法效果的影響，圖8a給出了運行時間隨聚類數(shù)目的變化曲線，可見當聚類數(shù)目的增加，算法運算量逐漸增加，導致運行時間增加。圖8b給出了去噪結(jié)果的RSNR隨聚類數(shù)目的變化曲線，表明對于該合成地震數(shù)據(jù)，在聚類數(shù)目為5時，算法具有最高的去噪效果，當聚類數(shù)目超過5時，隨著聚類數(shù)目的增加，RSNR有下降趨勢，說明不同的聚類數(shù)目參數(shù)會影響算法整體的運行時間與效果。

3.2 Marmousi模型地震數(shù)據(jù)實驗分析

在Marmousi數(shù)據(jù)模型中隨機抽取單炮數(shù)據(jù)，截取原始地震數(shù)據(jù)207道，每道700個采樣點(圖9a)，加入噪聲強度比例因子為0.03的高斯隨機噪聲數(shù)據(jù)(圖8b)，RSNR為30.46dB。使用該數(shù)據(jù)進行本文算法與小波變換、曲波變換及K-SVD算法的效果對比。圖10a為小波變換去噪結(jié)果，小波基采用db5，分解級數(shù)為5，去噪之后噪聲得到一定的壓制，其RSNR為32.74dB；圖10b為Curvelet變換去噪結(jié)果，Curvelet分解級數(shù)為5，由于Curvelet具有較好的方向性，去噪效果相對較好，RSNR為33.66dB。圖10c為K-SVD字典法去噪結(jié)果，初始化采用離散余弦基構(gòu)造超完備字典，以可重疊的8×8固定大小的數(shù)據(jù)塊為單位，構(gòu)造的超完備字典可捕捉整體目標數(shù)據(jù)的主要特征，逼近原始數(shù)據(jù)，RSNR為33.84dB；圖10d為本文算法去噪結(jié)果，由于采用了結(jié)構(gòu)聚類字典學習的思想，對原始數(shù)據(jù)進行分類稀疏表示后，提高了去噪效果，RSNR為38.43dB。

圖8 運行時間(a)及RSNR(b)隨聚類數(shù)目的變化曲線

圖9 Marmousi模型合成單炮數(shù)據(jù)(a)及其加噪結(jié)果(b)

圖10 Marmousi模型數(shù)據(jù)不同方法隨機噪聲壓制效果(a)小波變換法; (b)Curvelet變換法; (c)K-SVD字典法; (d)本文算法

為測試本文算法對不同強度高斯隨機噪聲的敏感程度，表1分別給出了小波變換法、Curvelet變換法、K-SVD自適應學習字典法及本文算法在不同強度隨機噪聲下的去噪效果(RSNR)對比。從表1可以看出，隨著加入高斯隨機噪聲的標準差從地震數(shù)據(jù)標準差的0.01增加到0.08，去噪的難度逐漸增加，各個算法的去噪效果變差，相比之下本文算法在各個噪聲強度下均能取得最好的去噪效果。

表1 不同算法在加入不同強度高斯隨機噪聲下噪聲壓制后的RSNR(單位：dB)對比

3.3 實際地震數(shù)據(jù)實驗分析

圖11為M區(qū)原始疊后海洋地震數(shù)據(jù)片段，其中道數(shù)為300、樣點數(shù)為500。圖12a為小波變換去噪結(jié)果，RSNR=23.16dB； 圖12b為Curvelet變換去噪結(jié)果，RSNR=24.42dB； 圖12c為K-SVD字典學習法的去噪結(jié)果，RSNR=24.53dB； 圖12d為本文算法的去噪結(jié)果，RSNR=25.25dB。由此可見本文算法利用地震數(shù)據(jù)結(jié)構(gòu)聚類字典學習的思想，在地震數(shù)據(jù)同向軸連續(xù)區(qū)域，能夠較好地提取出地震數(shù)據(jù)的主要特征，稀疏表示地震數(shù)據(jù)，壓制隨機噪聲(如矩形框所示)。但在地震波能量較弱的區(qū)域，由于原始地震數(shù)據(jù)波前曲線特征不明顯，地震數(shù)據(jù)與噪聲數(shù)據(jù)振幅相當?shù)那闆r下，地震數(shù)據(jù)稀疏表示的能力受到影響，存在細節(jié)丟失的現(xiàn)象(如橢圓區(qū)域所示)。

圖11 實際地震剖面

圖12 實際數(shù)據(jù)不同方法噪聲壓制結(jié)果對比(a)小波變換法； (b)Curvelet變換法； (c)K-SVD字典法； (d)本文算法

4 結(jié)論

為提高超完備字典學習方法對地震數(shù)據(jù)的稀疏表示能力與去噪效果，本文提出結(jié)構(gòu)聚類字典學習稀疏表示的地震數(shù)據(jù)去噪方法。通過理論分析與模型實驗可以得出以下結(jié)論：

(1)利用結(jié)構(gòu)聚類字典學習可以提高地震數(shù)據(jù)稀疏表示的能力，減少全局字典稀疏表示系數(shù)的冗余程度，有利于地震數(shù)據(jù)去噪處理；

(2)基于結(jié)構(gòu)聚類字典學習的去噪算法效果和運行時間受聚類數(shù)目參數(shù)影響較大，當聚類數(shù)目與實際地震數(shù)據(jù)分塊類別接近時，去噪效果較好，而運行時間隨聚類數(shù)目增加而增加；

(3)本文算法對地震數(shù)據(jù)中能量弱的局部區(qū)域去噪會出現(xiàn)細節(jié)丟失的現(xiàn)象，該問題是今后重點改進的研究方向。