主成分分析和字典學(xué)習(xí)聯(lián)合地震數(shù)據(jù)去噪方法

朱鶴文,韓立國，陳瑞鼎

吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，長春 130026

0 引言

在地震資料采集過程中，由于受到周圍環(huán)境和地下因素的影響，不可避免地會采集到隨機噪聲。這會對有效信號產(chǎn)生嚴重干擾甚至導(dǎo)致有效信號變形，從而使地震數(shù)據(jù)的質(zhì)量降低。為獲得高信噪比的地震數(shù)據(jù)，使地震資料能更真實地反映地下的地質(zhì)情況，如何突出有效波的同時抑制干擾波成為地震數(shù)據(jù)處理的主要研究內(nèi)容之一。為解決這個問題，前人嘗試了許多方法來壓制地震數(shù)據(jù)中的隨機噪聲[1-2]，目前最常用的方法有傅里葉變換[3]、小波變換[4-6]和多項式擬合[7]等，雖然這些方法在一定程度上可以抑制隨機噪聲的影響，但各自均存在一些不足。例如，傅里葉變換是一種全時域變換，不能很好地區(qū)分地震數(shù)據(jù)的局部特征，小波變換對地震數(shù)據(jù)邊緣的處理能力較弱等。

近年來，隨著稀疏表示理論[8]的興起，開始利用地震資料在某個變換域內(nèi)的可稀疏性和可分離性去除隨機噪聲的影響[9-10],例如使用一系列的固定基字典方法消除噪聲。由于這些方法的基是提前選定的，在數(shù)據(jù)處理過程中并不會隨數(shù)據(jù)特征的變化而變化，導(dǎo)致其只對某些信號有較好的處理效果，不具備自適應(yīng)性。后來提出了根據(jù)數(shù)據(jù)特點自適應(yīng)的改變基函數(shù)的方法，如Olshausen et al.[11]在完備字典的基礎(chǔ)上提出了自適應(yīng)學(xué)習(xí)型超完備字典，該字典通過不斷的學(xué)習(xí)訓(xùn)練更新字典獲得更稀疏的信號表示。唐剛等[12]將學(xué)習(xí)型超完備字典用于地震數(shù)據(jù)去噪并取得了很好的效果。

超完備字典的精確性決定了地震數(shù)據(jù)稀疏表示效果的好壞，傳統(tǒng)的K-SVD字典中使用奇異值分解來更新字典原子中的誤差項，其重構(gòu)的誤差項與真實誤差項之間存在一定殘差，因此得到的字典不是最優(yōu)的。本文將K-SVD字典與PCA算法相結(jié)合，利用PCA代替SVD對誤差項進行特征提取，用提取到的最大特征向量來更新字典原子，PCA重構(gòu)的誤差項與真實誤差項之間殘差更小，訓(xùn)練后的字典對數(shù)據(jù)稀疏表示效果更好。筆者通過對復(fù)雜模型合成地震記錄以及實際地震記錄進行去噪實驗，對比PCA、K-SVD以及本文算法對地震數(shù)據(jù)的去噪效果，驗證了本文方法在地震數(shù)據(jù)去噪中更有優(yōu)勢。

1 方法原理

1.1 主成分分析算法

主成分分析法(principal component analysis,PCA)是一種基于統(tǒng)計理論的字典學(xué)習(xí)算法，被廣泛應(yīng)用于數(shù)據(jù)壓縮[13]、人臉識別[14]和圖像去噪[15-16]等領(lǐng)域。PCA算法主要思想是將方差的大小視為對信息衡量的標準，方差越大，代表其所提供的信息就越多，反之，則越少。PCA算法的實現(xiàn)是通過對數(shù)據(jù)的協(xié)方差矩陣進行特征值分解，分解后的特征向量對應(yīng)數(shù)據(jù)的主成分，與特征向量相對應(yīng)的特征值代表數(shù)據(jù)在各個主成分上的權(quán)重。PCA是一種線性變換，變換后各個主成分分量之間彼此線性無關(guān)。隨著主成分編號的增加,該分量所包含的信息量就越小。筆者認為除主要成分外，其他主成分分量為噪聲分量,這使得這種變換在損失數(shù)據(jù)較小的情況下去除了噪聲。根據(jù)PCA的原理，結(jié)合地震數(shù)據(jù)實際情況，給出PCA算法的過程：

對于由n個檢波點與時間上m個采樣點得到的地震數(shù)據(jù)矩陣X(X∈Rm×n),PCA算法主要分為以下4步：

(1)將數(shù)據(jù)中心化處理：

(1)

(2)計算協(xié)方差矩陣

(2)

(3)將協(xié)方差矩陣特征值分解

C=UΛUT

(3)

式中:C為數(shù)據(jù)的協(xié)方差矩陣；Λ(Λ∈Rm×m)表示特征值矩陣，是一個對角矩陣，Λ中元素數(shù)值越大，代表其對應(yīng)的特征向量所包含的信息越多。U(U∈Rm×m)表示特征值矩陣Λ對應(yīng)的特征向量矩陣，UT(UT∈Rm×m)是U的轉(zhuǎn)置，被稱為旋轉(zhuǎn)矩陣，用來將數(shù)據(jù)轉(zhuǎn)換為主成分。

(4)去除含有噪聲的低階主成分，完成數(shù)據(jù)重構(gòu)：

(4)

1.2 K-SVD字典學(xué)習(xí)

K-SVD[16]是K-means算法的一種泛化形式。其主要思想是利用奇異值分解方法對誤差矩陣進行更新，從而得到更新后的字典原子以及相應(yīng)的稀疏系數(shù)。K-SVD字典學(xué)習(xí)算法求解如下優(yōu)化問題：

(5)

式中:Y表示給定訓(xùn)練數(shù)據(jù)；D表示字典矩陣；X表示稀疏系數(shù)矩陣；F表示F范數(shù)；xi是稀疏編碼矩陣X中的行向量，代表字典矩陣的系數(shù)；T0表示稀疏系數(shù)中非零元素最大值。

式(1)是一個有約束問題，利用拉格朗日乘子將其轉(zhuǎn)化為無約束問題：

(6)

K-SVD字典學(xué)習(xí)過程可分為3個步驟：

(1)準備階段：給定的原始訓(xùn)練集為Y。給定初始化字典為D，初始化字典可以通過兩種方式獲得：①從樣本中提取初始化字典；②給定一個字典作為初始化字典。

(2)稀疏編碼階段：根據(jù)給出的初始化字典，根據(jù)訓(xùn)練集Y與初始化字典矩陣D求解稀疏系數(shù)矩陣X。求解方法有很多，如匹配追蹤算法[17]、正交匹配追蹤算法[18-19]、快速迭代收縮閾值算法[20]等，本文采用正交匹配追蹤算法。

(7)

(8)

(9)

1.3 主成分分析和字典學(xué)習(xí)聯(lián)合方法的去噪原理

在K-SVD算法中，通過對誤差項進行奇異值分解來更新字典原子，這個過程屬于對誤差項的近似重構(gòu)，與真實誤差項存在一定誤差。因此在K-SVD算法中，利用奇異值分解更新字典原子所構(gòu)成的字典矩陣并不是最優(yōu)的。本文結(jié)合K-SVD字典理論與PCA算法強大的特征提取能力，提出使用PCA算法代替SVD算法對誤差項進行特征提取，用PCA算法提取的最大特征向量來更新字典原子。

這一過程表示為：

(10)

2 數(shù)值模擬實驗

為驗證基于主成分字典學(xué)習(xí)方法在地震數(shù)據(jù)去噪方面的優(yōu)越性，本文對復(fù)雜模型進行數(shù)值模擬，分別用小波變換、PCA、K-SVD字典學(xué)習(xí)以及本文方法進行去噪處理比較。定義信噪比參數(shù)SNR(Signal Noise Ratio)來衡量去噪效果的好壞。

(11)

式中:s代表去除隨機噪聲后的地震數(shù)據(jù)；s0代表原始未加隨機噪聲的地震數(shù)據(jù)。

對合成記錄加入的高斯噪聲的概率密度服從高斯分布：

(12)

式中:均值μ=0；標準差σ=1；n為隨機噪聲的量級，即產(chǎn)生n個滿足標準正態(tài)分布的隨機噪聲。

為驗證本文方法的去噪性能，對復(fù)雜模型進行合成地震記錄(圖1a)，共有100道，道間距10 m，每道700個采樣點，采樣時間間隔0.001 s。對模擬地震記錄加入0.6量級的隨機噪聲(圖1b)，加入后信噪比為7.79 dB，對比圖1a和圖1b中紅色箭頭和黃色箭頭所指部分，可以看出隨機噪聲對有效信號的干擾很嚴重，造成有效信號能量較弱的部分模糊形變。

a.合成地震記錄；b.加噪聲記錄。圖1 復(fù)雜模型合成地震記錄和加噪聲地震記錄Fig.1 Complex model synthetic seismogram and seismic records with added noise

經(jīng)過筆者多次實驗，使用PCA對地震數(shù)據(jù)進行去噪時，選取前65～70個主成分(占比約0.945)，能取得最高信噪比。這里選取前67個主成分。對比圖2中4種方法的去噪效果可以看出，這4種算法均能對地震記錄中的隨機噪聲進行壓制，但小波變換(圖2a)對同相軸的表示能力不足，會使同相軸產(chǎn)生形變。PCA算法(圖2b)紅色箭頭處的同相軸依然模糊，表明其對能量薄弱部分的有效信號去噪效果明顯不如小波變換、K-SVD字典學(xué)習(xí)和本文方法，這是由于PCA算法會將弱能量的有效信號視為噪聲并將其去除，使能量弱的同相軸變得模糊。K-SVD算法(圖2c)和本文方法(圖2d)去噪效果明顯好于前兩種算法，對噪聲的壓制比較明顯，但對比之下可以發(fā)現(xiàn)，本文方法較K-SVD算法對隨機噪聲剔除的更干凈，從黃色箭頭部分也可以看出，本文方法對弱能量有效信號保護好于K-SVD算法。同時根據(jù)表1看出，本文算法去噪后信噪比較其他三者都有提升，其中較K-SVD算法提升約1.2 dB，表明本文方法在保護有效信號和壓制噪聲兩方面都好于K-SVD算法及其他兩種方法。

表1 復(fù)雜模型合成地震記錄去噪信噪比Table 1 Denoise SNR of complex model synthetic seismogram

a.小波去噪；b.PCA去噪；c.K-SVD去噪；d.本文方法。圖2 復(fù)雜模型合成地震記錄去噪效果Fig.2 Denoising effect of complex model synthetic seismogram

為驗證本文算法在不同噪聲強度下的去噪效果，對合成地震記錄加入不同量級的隨機噪聲，分別使用以上4種方法進行去噪處理并繪制去噪后的信噪比曲線(圖3)，可以看出隨著噪聲量級的增加，4種方法去噪后的信噪比都在減小，但通過對比可以發(fā)現(xiàn)，本文方法去噪后的信噪比明顯高于其他3種方法,相較于傳統(tǒng)的K-SVD算法高出1～1.5dB，說明在不同噪聲量級下，本文方法對K-SVD算法去噪效果有提升，這種提升在低量級隨機噪聲中表現(xiàn)尤為明顯。

圖3 復(fù)雜模型合成地震記錄在不同隨機噪聲量級下的信噪比曲線Fig.3 SNR curves of complex model synthetic seismograms under different random noise magnitude

3 實際地震數(shù)據(jù)及分析

為檢測本文方法在實際地震數(shù)據(jù)中的去噪能力，分別用小波變換、PCA、K-SVD以及本文方法對實際含噪聲地震數(shù)據(jù)進行去噪實驗對比。采用實際地震數(shù)據(jù)如圖4所示，地震記錄共61道，每道3 000個采樣點，時間采樣間隔為0.001 s。由紅色箭頭可以看出隨機噪聲對反射波同相軸產(chǎn)生嚴重畸變，導(dǎo)致反射波同相軸模糊，黑色箭頭處可以看到背景噪聲比較明顯。

圖4 實際地震記錄Fig.4 Actual seismogram

圖5為實際地震數(shù)據(jù)去噪效果圖，可以看出對于實際地震數(shù)據(jù)去噪處理中，4種方法都取得了不錯的效果。但小波變換(圖5a)對于紅色箭頭處反射波同相軸部分的去噪效果不理想，PCA算法(圖5b)同小波變換一樣，對黑色箭頭處的背景噪聲壓制有一定效果，紅色箭頭處有效信號部分去噪效果不明顯，同時PCA算法在去除隨機噪聲的同時引入了條帶狀噪聲。相比之下K-SVD算法(圖5c)與本文方法(圖5d)對于背景噪聲的壓制明顯好于小波變換和PCA算法，但K-SVD算法在去噪過程中，對有效信號的損失比較嚴重，導(dǎo)致同相軸模糊，而本文算法無論是在去噪效果還是在對有效信號的保護上都要好其他3種算法。

a.小波去噪；b.PCA去噪；c.K-SVD去噪；d.本文方法。圖5 實際地震數(shù)據(jù)去噪效果Fig.5 Denoising effect of actual seismogram

4 結(jié)論

(1)對復(fù)雜模型合成地震記錄的實驗可知，在不同強度的隨機噪聲下，本文算法都能獲得比其他3種算法更高的信噪比，在噪聲能量較弱時表現(xiàn)尤為明顯。

(2)對實際地震記錄的實驗可知，相較于K-SVD及其他2種算法，本文算法對有效信號的損害較少，噪聲去除的更干凈，同相軸更清晰。