基于自適應(yīng)K-SVD的能量泄漏恢復(fù)研究

何真,曹思遠(yuǎn),郝婳婕,段緒遠(yuǎn),葉力夏提,張韻佩

(1.中國石油新疆油田分公司 勘探開發(fā)研究院，新疆 克拉瑪依 834000;2.中國石油大學(xué)(北京) 地球物理學(xué)院，北京 102249;3.中國石油大學(xué)(北京) CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京 102249;4.中國石油新疆油田分公司 數(shù)據(jù)公司信息研究所，新疆 克拉瑪依 834000;5.中國石油新疆油田分公司 新港公司，新疆 克拉瑪依 834000)

0 引言

在地震資料處理中，隨機(jī)噪聲的壓制研究一直都是熱點(diǎn)與難點(diǎn)，處理結(jié)果的好壞直接影響著后續(xù)解釋與儲層預(yù)測等工作。目前生產(chǎn)實(shí)際中常用的去噪手段有小波變換、奇異值分解(SVD)、F-X域去噪、Radon變換、三維塊匹配算法(BM3D)等一系列方法。其中，小波變換[1]作為一種經(jīng)典的數(shù)學(xué)變換，對于不同特征的帶噪信號采用軟閾值函數(shù)或硬閾值函數(shù)進(jìn)行降噪處理，雖有一定效果但仍然有待提高，該方法在提出以后也不斷有專家學(xué)者對其進(jìn)行討論與改進(jìn)[2-3]。奇異值分解算法(SVD)[4]借助地震信號之間的相關(guān)性與隨機(jī)噪聲之間的不相關(guān)來進(jìn)行信噪分離，根據(jù)能量的不同分為若干奇異值，大奇異值對應(yīng)有效信號，小奇異值對應(yīng)噪聲，但地震同相軸斜率較大時則不適用該方法。后來劉志鵬、趙偉等[5]提出局部頻率域SVD對地震波進(jìn)行去噪，基于有效頻帶頻譜存在較好的相關(guān)性來進(jìn)行噪聲壓制。f-x域去噪[6-7]方法假定相干信號在f-x域可預(yù)測而隨機(jī)噪聲不可預(yù)測的原理對信號和噪聲加以區(qū)分來達(dá)到去噪的目的。Radon變換去噪主要有線性、拋物線和雙曲線拉東變換3種去噪方法[8-10]，根據(jù)同相軸的不同時差特征來選取合適的變換類型進(jìn)行去噪。三維塊匹配算法[11-12](BM3D)是近些年熱度較高的圖像去噪算法之一，它將噪聲圖像沿垂直與水平方向搜索相似塊并分組，再進(jìn)行協(xié)同濾波，該方法解決了圖像去噪后模糊的問題，保留了圖像的邊緣細(xì)節(jié)，在一定程度上優(yōu)于傳統(tǒng)的空間域降噪技術(shù)和頻率域降噪技術(shù)。

盡管去噪方法和理論在不斷改進(jìn)和完善，但在常規(guī)處理中由于參數(shù)選擇的不合理或者基于一些不合理假設(shè)導(dǎo)致在對信號進(jìn)行降噪的同時難免會“傷及”一部分有效信號，即殘差剖面中存在有效信號的情形，稱其為能量泄漏現(xiàn)象。對于噪聲中含有的有效信號部分很少有人進(jìn)行研究。本文提出基于優(yōu)質(zhì)樣本的學(xué)習(xí)字典進(jìn)行能量恢復(fù)研究，優(yōu)質(zhì)樣本構(gòu)建的字典包含較多的原地震信息以及較少的隨機(jī)噪聲，并且與殘差剖面中殘留的有效信號有較強(qiáng)的相關(guān)性，可以有效地將其提取并回收利用。

1 自適應(yīng)K-SVD算法能量泄漏恢復(fù)原理

1.1 K-SVD算法

實(shí)際地震數(shù)據(jù)可表示為：

f=s+σ*n,

(1)

其中：f為含噪地震數(shù)據(jù)；s為有效信號；n為隨機(jī)噪聲；σ為噪聲強(qiáng)度，一般不可知。壓制隨機(jī)噪聲的過程就是將σ*n這一部分從地震數(shù)據(jù)f中去除，現(xiàn)用任意手段將信號f進(jìn)行去噪處理得到：

D_noise=Method(f) ,

(2)

其中：D_noise表示去噪后信號，Method表示任一種去噪方法，可以是小波閾值去噪、SVD分解、中值濾波等。

一次去噪后信號D_noise即本文研究的部分，將其寫成如下形式：

D_noise=D*X,

(3)

i=1,2,…,n

(4)

T0為稀疏表示中非零分量個數(shù)的最大值。字典D的更新是逐列進(jìn)行，假設(shè)字典D和稀疏矩陣X都是固定的，當(dāng)更新字典D中第h列原子dh時，則有：

(5)

(6)

(7)

(8)

(9)

Σ=diag(σ1,σ2,σ3,…,σr) 。

(10)

σr(r=1,2,3,…,r)表示矩陣的所有非零奇異值，r為矩陣的秩。

(11)

(12)

1.2 正交匹配追蹤(OMP)

OMP算法是根據(jù)匹配度選擇最優(yōu)的原子進(jìn)入原子集，它在更新信號殘差時，減去的是當(dāng)前殘差在所有已被選擇的原子上的正交投影，這樣就使得該算法在迭代過程中避免對已選的原子的重復(fù)挑選。并且通過求解一個最小二乘問題找到原信號的最優(yōu)近似解，再對殘差進(jìn)行更新迭代，最后用原子的線性組合來表征信號。具體算法如下：

1)輸入m×n的字典D以及需要稀疏重構(gòu)的信號y，給定預(yù)設(shè)誤差ε,稀疏度K；

2)初始化殘差ε0=y，索引集Λ0=?，迭代次數(shù)t=1；

3)將字典D每列與殘差εt做內(nèi)積，選擇內(nèi)積最大的一個即信號在原子垂直投影長度中最長的一個，并記錄其索引φt；

4)令Λt=Λt-1∪{φt};

7)判斷是否滿足εt<ε，若滿足，停止迭代，否則返回3)。

1.3 相關(guān)系數(shù)匹配準(zhǔn)則

(13)

(14)

此時停止OMP迭代，輸出結(jié)果。

1.4 本文算法實(shí)現(xiàn)流程

圖1 本文算法流程圖Fig.1 Algorithm flow chart

2 實(shí)驗(yàn)結(jié)果與分析

實(shí)驗(yàn)分為對模型資料的處理以及實(shí)際資料的處理，首先通過楔形模型和傾斜模型的處理證明本文方法的正確性和有效性，再對實(shí)際資料進(jìn)行分析、去噪、能量恢復(fù)獲得高信噪比高保真度的地震剖面。

實(shí)驗(yàn)中為統(tǒng)一標(biāo)準(zhǔn)，設(shè)置地震數(shù)據(jù)塊大小為8×8，字典中的原子數(shù)為256，初始化字典采用DCT字典。在評價去噪效果時采用能量疊加法[15]計(jì)算信噪比，公式為：

(15)

(16)

式中：M、N分別為采樣點(diǎn)數(shù)和地震道數(shù)，SNR為信噪比，E為地震道總能量，EN為噪聲總能量，ES為有效信號能量。

2.1 楔形模型去噪分析

圖2a是楔形模型，道數(shù)30道，采樣點(diǎn)數(shù)500，子波主頻35 Hz；圖2b是加入了50%的隨機(jī)噪聲，其信噪比為2.18 dB?，F(xiàn)在利用KSVD算法分別對無噪和帶噪情況下的模型求取特征字典，得到圖2c(帶噪字典)，圖2d(無噪字典)。不難發(fā)現(xiàn)在無噪情況下，特征字典提取的都是有效信號的成分，而噪聲情況下的字典則較為“模糊”?，F(xiàn)對加噪信號分別采用帶噪和無噪的字典進(jìn)行稀疏求解，得到去噪結(jié)果圖2e和圖2f。信噪比分別為18.34 dB和20.05 dB，兩種方法對噪聲均有較好的壓制作用，但對比殘差剖面圖2g、h和單道15道信號圖2i、j、k、l來看，采用無噪字典去噪的剖面在信號的去噪效果和保幅性上更勝一籌。

2.2 傾斜模型能量泄漏恢復(fù)效果分析

上文2.1中通過楔形模型證明利用一個優(yōu)質(zhì)的學(xué)習(xí)樣本(即無噪聲或噪聲盡可能少的情況下的剖面)進(jìn)行字典學(xué)習(xí)對噪聲壓制以及信號重構(gòu)的效果要優(yōu)于傳統(tǒng)對帶噪信號進(jìn)行學(xué)習(xí)(即普通樣本)的效果。接下來采用傾斜模型，對其加入一定比例隨機(jī)噪聲，先利用傳統(tǒng)去噪手段對其進(jìn)行去噪處理，再分別對去噪后剖面和原始帶噪剖面作為樣本進(jìn)行字典學(xué)習(xí)，將學(xué)習(xí)到的字典對殘差剖面進(jìn)行稀疏分解，實(shí)現(xiàn)信號和噪聲的二次分離，最后將分離出的有效信號補(bǔ)回去噪后的剖面并分析討論兩種方法的恢復(fù)結(jié)果。

a—楔形模型;b—加噪楔形模型;c—噪聲情況下構(gòu)建的字典;d—無噪聲情況下構(gòu)建的字典;e—噪聲字典去噪剖面;f—無噪聲字典去噪剖面;g—噪聲字典去噪后殘差;h—無噪聲字典去噪后殘差;i—楔形模型第15道信號;j—加噪楔形模型第15道;k—噪聲字典去噪后第15道;l—無噪字典去噪后第15道a—wedge model;b—noise wedge model;c—dictionary constructed under noise;d—noise dictionary denoising profile;e—noise dictionary denoising profile;f—noiseless dictionary denoising profile;g—noise dictionary denoising residuals;h—noiseless dictionary denoising residuals;i—wedge model 15th signal;j—noise wedge model 15th;k—noise dictionary denoising after 15th;l—No15 after noiseless dictionary denoising圖2 楔形模型試驗(yàn)結(jié)果Fig.2 Wedge model test results

圖3a為傾斜模型，道數(shù)81道，采樣點(diǎn)數(shù)801，子波主頻55 Hz；圖3b加入65%的隨機(jī)噪聲，信噪比為3.31 dB。利用傳統(tǒng)小波變換對其進(jìn)行去噪處理，得到圖3c，信噪比為19.62dB，去噪后剖面信噪比有了大幅度的提升，但從殘差剖面觀察(圖3d)，去除的噪聲中存在大量有效信號。以普通樣本和優(yōu)質(zhì)樣本為基礎(chǔ)進(jìn)行字典學(xué)習(xí)對殘差剖面進(jìn)行處理，得到圖3e和圖3f的信噪分離結(jié)果。從分離結(jié)果可以看出，以優(yōu)質(zhì)樣本學(xué)習(xí)的字典在殘差剖面中信噪分離的效果更加徹底，將其分別補(bǔ)回去噪后剖面得到能量恢復(fù)結(jié)果(圖3g、h)，信噪比分別為18.13dB和21.25 dB?？梢姳疚奶岢龅姆椒▽τ行盘柕奶崛「泳珳?zhǔn)和可靠，下文的單道信號分析中更能清晰地體現(xiàn)這一結(jié)論。

a—傾斜模型;b—加噪傾斜模型;c—小波閾值去噪剖面;d—小波閾值去噪殘差;e—普通樣本分離的有效信號;f—優(yōu)質(zhì)樣本分離的有效信號;g—普通樣本能量恢復(fù)結(jié)果;h—優(yōu)勢樣本能量恢復(fù)結(jié)果a—tilt model;b—charging tilt model;c—wavelet threshold denoising profile;d—wavelet threshold denoising residual;e—effective signal for normal sample separation;f—effective signal for good sample separation;g—ordinary sample energy recovery results;h—high quality sample energy recovery results圖3 傾斜模型剖面試驗(yàn)結(jié)果Fig.3 Tilt model section test results

對圖3傾斜模型時間域第75道，觀察普通樣本的字典學(xué)習(xí)恢復(fù)效果(圖4d)和優(yōu)質(zhì)樣本的字典學(xué)習(xí)恢復(fù)效果(圖4e)，優(yōu)質(zhì)樣本學(xué)習(xí)下恢復(fù)的結(jié)果形態(tài)更加接近原始信號，而普通樣本學(xué)習(xí)下恢復(fù)的結(jié)果中含有較多的隨機(jī)噪聲，信噪分離程度相對不夠理想。

a—傾斜模型;b—傾斜模型加噪;c—小波閾值去噪結(jié)果;d—普通樣本能量恢復(fù)結(jié)果;e—優(yōu)質(zhì)樣本能量恢復(fù)結(jié)果a—tilt model;b—tilt model noise;c—wavelet threshold denoising result;d—normal sample energy recovery result;e—good sample energy recovery result圖4 傾斜模型第75道信號試驗(yàn)結(jié)果Fig.4 Test results of the 75th signal of tilt model

為更好地體現(xiàn)本文方法的實(shí)用性，現(xiàn)將帶噪傾斜剖面分別采用f-x域去噪、改進(jìn)小波閾值[16]、ASVD[17]3種方法對其進(jìn)行一次去噪，再分別對3個去噪剖面進(jìn)行字典構(gòu)造進(jìn)而提取殘差剖面的有效信號，并對結(jié)果進(jìn)行分析。

由于補(bǔ)償恢復(fù)后的結(jié)果從時間域剖面中很難看出差異，所以對以上3種恢復(fù)結(jié)果只從單道信號進(jìn)行觀察討論。從去噪結(jié)果看，3種方法均有較好的效果，噪聲基本被去除，但f-x域去噪相對于改進(jìn)小波閾值和ASVD方法，在去除噪聲的同時也傷及了大量有效信號。仍然以恢復(fù)后的第75道為例，見圖5g、h、i，對f-x域去噪、改進(jìn)小波閾值以及ASVD三種去噪結(jié)果進(jìn)行能量恢復(fù)后隨機(jī)噪聲含量均比較小，但在有效信號振幅能量上后兩者恢復(fù)得更好，更加接近原始模型資料。經(jīng)對比也不難發(fā)現(xiàn)：原始剖面降噪效果越好，有效信號保留越多，在構(gòu)建字典原子時信號的特征越能被更好地提取和表達(dá)，從而實(shí)現(xiàn)更加理想的恢復(fù)效果。

a—f-x域去噪剖面;b—f-x域去噪殘差;c—改進(jìn)小波閾值去噪剖面;d—改進(jìn)小波閾值去噪殘差;e—ASVD去噪剖面;f—ASVD去噪殘差;g—f-x域去噪后能量恢復(fù)結(jié)果(第75道);h—改進(jìn)小波閾值去噪后能量恢復(fù)結(jié)果(第75道);i—ASVD去噪后能量恢復(fù)結(jié)果(第75道)a—f-x domain denoising profile;b—f-x domain denoising residual;c—improved wavelet threshold denoising profile;d—improved wavelet threshold denoising residual;e—ASVD denoising proflie;f—ASVD denoising residual;g—f-x domain denoising energy recovery result(75th channel);h—improved wavelet threshold denoising energy recovery result(75th channel);i—ASVD denoising energy recovery result(75th channel)圖5 本文方法實(shí)用性試驗(yàn)結(jié)果Fig.5 Result of the practical test of the method

2.3 實(shí)際資料處理及效果分析

選取某實(shí)際二維疊加數(shù)據(jù)進(jìn)行處理分析。圖6a為某區(qū)塊實(shí)際地震資料，剖面上隨機(jī)噪聲較為發(fā)育，能量較弱的信號基本被湮沒，給接下來的小層解釋、儲層預(yù)測等工作帶來了較大難度。實(shí)際資料的一次去噪選取fx域和ASVD法，進(jìn)行f-x域去噪，得到去噪結(jié)果圖6b和殘差剖面圖6c，采用本文方法對殘差剖面處理實(shí)現(xiàn)噪聲信號二次分離，將分離的有效信號加回去噪后剖面得到能量恢復(fù)結(jié)果(圖6d)，同樣地，采用上述方法得到ASVD的處理結(jié)果(圖6e、f、g)，對于一次去噪結(jié)果(圖6b和圖6e)，兩種方法均去除大部分隨機(jī)噪聲，使得有效波變得清晰，但殘差剖面中仍能發(fā)現(xiàn)較多的有效信號。而恢復(fù)后的剖面圖6d、g中同相軸信息被較好地提取(見紅圈及紅框部分)，剖面變得更加自然、連續(xù)。不難發(fā)現(xiàn)，ASVD處理后再進(jìn)行能量恢復(fù)的剖面在細(xì)節(jié)表現(xiàn)上要優(yōu)于f-x去噪后再恢復(fù)的結(jié)果。

a—某地區(qū)二維實(shí)際資料;b—f-x域去噪剖面;c—f-x域去噪殘差;d—f-x域去噪后恢復(fù)剖面;e—ASVD去噪剖面;f—ASVD去噪殘差;g—ASVD去噪后恢復(fù)剖面a—2D actual data in a region;b—f-x domain denoising profile;c—f-x domain denoising residual;d—f-x domain denoising recovery profile;e—ASVD denoising proflie;f—ASVD denoising residual;g—ASVD domain denoising recovery profile圖6 實(shí)際資料處理Fig.6 Actual data processing

3 結(jié)論與建議

本文基于K-SVD字典學(xué)習(xí)方法進(jìn)行地震能量泄露恢復(fù)研究，通過楔形、傾斜模型以及實(shí)際資料處理分析得到以下結(jié)論:

1)以無噪或者噪聲盡可能小的剖面作為樣本進(jìn)行字典學(xué)習(xí)、稀疏表示得到的去噪結(jié)果優(yōu)于傳統(tǒng)的字典學(xué)習(xí)方法。

2)利用一次去噪后的剖面作為樣本進(jìn)行學(xué)習(xí)，原始剖面的降噪性、保真性越好，在構(gòu)建字典原子時信號的特征越能被更好的提取和表達(dá)，從而實(shí)現(xiàn)更加理想的恢復(fù)效果。

3)根據(jù)不同的地震信號可以構(gòu)建不同的特征字典，具有較好的自適應(yīng)性，并且在適用范圍上有一定應(yīng)用前景。