朱鶴文,韓立國(guó)，陳瑞鼎

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

0 引言

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

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

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

1 方法原理

1.1 主成分分析算法

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

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

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

(1)

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

(2)

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

C=UΛUT

(3)

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

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

(4)

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

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

(5)

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

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

(6)

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

(1)準(zhǔn)備階段：給定的原始訓(xùn)練集為Y。給定初始化字典為D，初始化字典可以通過(guò)兩種方式獲得：①?gòu)臉颖局刑崛〕跏蓟值?；②給定一個(gè)字典作為初始化字典。

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

(7)

(8)

(9)

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

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

這一過(guò)程表示為：

(10)

2 數(shù)值模擬實(shí)驗(yàn)

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

(11)

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

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

(12)

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

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

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

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

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

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

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

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

圖4 實(shí)際地震記錄Fig.4 Actual seismogram

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

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

4 結(jié)論

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

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