基于字典學(xué)習(xí)算法的地震資料去噪處理應(yīng)用研究

楊海濤, 李 瓊,1b, 錢(qián)浩東, 宋 鑫, 雍 鵬

(1.成都理工大學(xué) a.地球物理學(xué)院;b.地球勘探與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，成都 610059;2.中國(guó)石油集團(tuán) 川慶鉆探工程有限公司鉆采工程技術(shù)研究院，廣漢 618300)

0 引言

信噪比、分辨率、保真度是地震數(shù)據(jù)質(zhì)量的三大評(píng)價(jià)標(biāo)準(zhǔn)，其中高信噪比是三者(高信噪比、高分辨率、高保真度)的基礎(chǔ)[1]。由于野外采集數(shù)據(jù)的時(shí)候總會(huì)受到周?chē)h(huán)境的影響以及多次波的存在，使去噪成為解釋作業(yè)之前的重中之重。相較于已有的去噪方法諸如：離散余弦變換[2]、F-K 域?yàn)V波[3]、小波變換[4]、Radon變換[5]等，壓縮感知理論[6-7]的出現(xiàn)，打破了之前采樣間隔的限制，在稀疏及重建提供了新方向。諸多學(xué)者提出了一系列稀疏重建算法和字典學(xué)習(xí)方法(ISTA算法、MP算法、BP算法以及DCT字典等)，后續(xù)又有不少學(xué)者進(jìn)行優(yōu)化得到了OMP算法、FISTA算法以及K-svd字典。

基于壓縮感知的去噪，就是將實(shí)際數(shù)據(jù)在稀疏域中采用對(duì)應(yīng)的稀疏基變成系數(shù)矩陣，再通過(guò)系數(shù)矩陣來(lái)重建地震道，再對(duì)系數(shù)矩陣進(jìn)行處理，最終達(dá)到預(yù)期去噪效果。筆者應(yīng)用OMP算法解決K-SVD字典的迭代問(wèn)題，并與DCT字典對(duì)比研究。

1 方法原理

1.1 壓縮感知理論(Compression sensing，CS)

壓縮感知改變了以往抽樣定理采樣的前提——信號(hào)x(t)是有限帶寬的，其頻譜的最高頻為fc，壓縮感知采樣的前提是 “信號(hào)x(t)可以在某變換域稀疏表示，其稀疏度為K”。與之前相比，CS是邊采樣邊壓縮，將模擬信息轉(zhuǎn)化為AIC，大大降低了地震數(shù)據(jù)的采集數(shù)據(jù)量。壓縮感知可分為三步：稀疏、觀測(cè)以及重建(圖1(b))。

圖1 信號(hào)采樣壓縮流程Fig.1 Signal sampling compression process(a)傳統(tǒng)信號(hào)壓縮過(guò)程;(b)壓縮感知理論框架

壓縮感知的數(shù)學(xué)模型可以表示成式(1)。

y=Cf

(1)

其中：y為長(zhǎng)度為m的采樣結(jié)果即觀測(cè)信號(hào)；C∈Rm×n(m?n)為m×n維的采樣矩陣；f為長(zhǎng)度為n的待采樣信號(hào)。由此可知，壓縮感知需要的采樣數(shù)據(jù)量與之前相比大大降低了，但所采集的地震數(shù)據(jù)一般都是時(shí)間域內(nèi)的非稀疏信號(hào)，因此信號(hào)是需要稀疏表示來(lái)滿(mǎn)足CS采樣的前提——信號(hào)在某變換域內(nèi)是K稀疏的。常用來(lái)做稀疏基的有小波、曲波(Cuverlet)、Seislet、Radon、離散余弦等。根據(jù)信號(hào)的稀疏化表示，絕大多數(shù)地震信號(hào)都能在某一個(gè)變換域內(nèi)對(duì)其稀疏化表示，因此可表示為[8]式(2)。

f=Ψx

(2)

其中：Ψ∈Rn×n是稀疏矩陣；f是被信號(hào)x的稀疏表示。

將式(1)代入式(2)式可得式(3)。

y=Cf=CΨx

(3)

可以將CΨ合并成一個(gè)矩陣B∈Rm×n，可得：

y=Bx

(4)

由于矩陣B是一個(gè)m×n維矩陣，且m?n，因此方程組是一個(gè)欠定方程組。Candes等[6]提出有限等距性質(zhì)(Restricted Isometry Principle, RIP)即：零空間特性、約束等距性和相干程度低。由上述的理論可知(常數(shù)矩陣)δk∈(0,1)應(yīng)需滿(mǎn)足：

(5)

觀測(cè)矩陣一般使用Gauss矩陣，因?yàn)殡S機(jī)高斯矩陣與大部分固定稀疏基矩陣是不相干的。對(duì)地震數(shù)據(jù)進(jìn)行觀測(cè)，故選用正交基作為變換基時(shí)，B是滿(mǎn)足RIP性質(zhì)的[9]。

當(dāng)滿(mǎn)足RIP條件時(shí)，信號(hào)重構(gòu)最直接的方法是L0范數(shù)的最優(yōu)化問(wèn)題：

(6)

由于式(6)的求解是個(gè)NP(欠定)問(wèn)題。Candès[6]指出，可以把復(fù)雜性和不穩(wěn)定性較高的L0范數(shù)最優(yōu)化問(wèn)題，轉(zhuǎn)化為等價(jià)的L1范數(shù)最優(yōu)化問(wèn)題，因此式(6)可轉(zhuǎn)化為：

argmin‖x‖1s.ty=Bx

(7)

而關(guān)于信號(hào)重構(gòu)，目前可大致分為兩大類(lèi)：

1)凸規(guī)劃算法。在稀疏重建模型中，Candès[5]指出可以把復(fù)雜性和不穩(wěn)定性較高的L0范數(shù)最優(yōu)化問(wèn)題，轉(zhuǎn)化為等價(jià)的L1范數(shù)最優(yōu)化問(wèn)題。諸多學(xué)者相繼提出基追蹤(BP)、梯度下降法、凸集交替投影法(POCS)以及迭代閾值法等方法。

2)貪婪迭代算法。這類(lèi)方法通過(guò)從感知矩陣列中選出一個(gè)和觀測(cè)矩陣最匹配的原子，經(jīng)過(guò)不停地迭代重構(gòu)出稀疏信號(hào)，主要包括匹配跟蹤算法以及其優(yōu)化算法——正交匹配追蹤算法等。故OMP算法被經(jīng)常使用，是基于MP算法的基礎(chǔ)做斯密特正交化，使迭代速率明顯加快了。

1.2 信號(hào)質(zhì)量客觀評(píng)價(jià)方法

信號(hào)質(zhì)量客觀評(píng)價(jià)方法大致以下幾種：①信噪比(Signal to Noise Ratio, SNR);②均方誤差(Mean Square Error, MSE);③峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)[10]。

信噪比：

(8)

其中：x、y分別為去噪信號(hào)和原始信號(hào)，且均為M*N階矩陣。

峰值信噪比：

(9)

其中，MSE為均方差誤差。

1.3 正交匹配追蹤算法

正交匹配跟蹤算法(Orthogonal Matching Pursuit，OMP)的原理是基于MP算法原理，是利用篩選與誤差相關(guān)性最大的原子，利用原子的線性組合去重構(gòu)出最優(yōu)近似信號(hào)[11]，OMP算法具體流程如圖2所示。

圖2 OMP算法流程圖Fig.2 Algorithm flowchart of OMP

1.4 DCT字典

DCT 字典是通過(guò)離散余弦變換的正交變換來(lái)得到的，能較好地對(duì)周期信號(hào)進(jìn)行分辨，其特點(diǎn)為能量集中，低頻部分為其主要集中部分，同時(shí)可以去掉高頻信號(hào)中的無(wú)用的一部分，因此適用于對(duì)地震資料局部的特征進(jìn)行處理反變換為式(11)。

給定一組序列y(n)，經(jīng)過(guò)DCT變換得到式(10)。

(10)

(11)

式中0≤n≤N-1，其矩陣形式為：

Yc=CNy

(12)

其中CN經(jīng)過(guò)DCT變換就可以得到字典D(圖3)。

圖3 DCT字典Fig.3 DCT dictionary

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

K-SVD是一種字典學(xué)習(xí)的構(gòu)建算法，從矩陣角度分析，其目標(biāo)就職是把輸入的訓(xùn)練樣本F矩陣分解成兩部分：(字典)D和(稀疏矩陣)S。K-SVD字典學(xué)習(xí)其實(shí)是K-means算法的一種改進(jìn)算法，利用稀疏約束條件來(lái)進(jìn)行制約，然后利用逼近信號(hào)與原始信號(hào)差值的奇異值分解來(lái)不斷迭代更新，求解有兩個(gè)階段，即稀疏編碼和字典學(xué)習(xí)[12]。與固定字典DCT相比，K-SVD字典多了更新訓(xùn)練的過(guò)程，在對(duì)地震資料進(jìn)行去噪處理時(shí)，是通過(guò)對(duì)實(shí)際信號(hào)所得出的自適應(yīng)字典，比固定字典更具穩(wěn)定性。首先假設(shè)Y為M*N階訓(xùn)練樣本矩陣，K-SVD算法中求解系數(shù)矩陣時(shí)采用OMP算法。在更新字典階段：

(13)

圖4 K-SVD學(xué)習(xí)字典Fig.4 K-SVD learning dictionary

圖5 K-SVD學(xué)習(xí)字典流程Fig.5 Algorithm flowchart of K-SVD

3 模型測(cè)試

為了證明K-SVD與DCT變換矩陣下的OMP算法的去噪效果，筆者選用主頻30 Hz的雷克子波，設(shè)計(jì)了如圖6(a)的三層模型的彈炮合成記錄與圖6(b)所示的六層地層模型，圖6(a)設(shè)計(jì)為總長(zhǎng)為2 560 m，道間距為5 m，波速分別為2 000 m/s、2 400 m/s、2 300 m/s，因此該模型513道接收。圖6(b)設(shè)計(jì)為總長(zhǎng)為5 120 m，道間距為10 m，波速分別為2 300 m/s、3 200 m/s、5 500 m/s、3 700 m/s、6 000 m/s，因此該模型也為513道接收。

圖6 原始模型及加噪單炮合成記錄Fig.6 The original model and noisy single shot synthesis record(a)三層單炮記錄；(b)六層單炮記錄；(c)三層加噪單炮記錄；(d)六層加噪單炮記錄

圖6(c)和6(d)是為加噪模型的單道地震合成記錄，對(duì)比原始模型可以明顯看到噪音的存在。以?xún)煞N字典學(xué)習(xí)型進(jìn)行去噪，從而得到如圖7所示的去噪效果圖。對(duì)于不同地層，就峰值信噪比來(lái)說(shuō)，表1的PSNR值是由Matlab計(jì)算出來(lái)的，從表1中可以看出，兩種學(xué)習(xí)字典的峰值信噪比PSNR都隨著地層的增加而減小；為了更清晰對(duì)比，由圖7(e)與7(f)、7(g)與7(h)可知，對(duì)同一地層來(lái)說(shuō)，從圖7中紅線框中可以看DCT字典對(duì)同相軸的損害要更高一點(diǎn)，可以看到K-SVD字典明顯比DCT字典的效果好，但K-SVD較之DCT，由于在字典的稀疏編碼階段通過(guò)奇異值分解不斷地更新字典，在噪聲強(qiáng)度較大的情況下，去噪效果有了顯著的提高，且具有更強(qiáng)的穩(wěn)定性；缺點(diǎn)是由于迭代的次數(shù)比較多，導(dǎo)致運(yùn)行速度較慢。K-SVD算法中求解系數(shù)矩陣時(shí)采用OMP算法，是在MP的基礎(chǔ)上進(jìn)行施密特正交化，這樣就不用再多次選擇原子，故減少了運(yùn)算時(shí)間。

表1 降噪方法的評(píng)價(jià)指標(biāo)表Tab.1 Noise reduction evaluation index table

圖7 去噪圖及殘差圖Fig.7 Result and residual graph(a)三層模型DCT效果圖；(b)三層模型K-SVD效果圖；(c)六層模型DCT效果圖；(d)六層模型K-SVD效果圖；(e)三層模型DCT殘余圖；(f)三層模型K-SVD殘余圖；(g)六層模型DCT殘余圖；(h)六層模型K-SVD殘余圖

4 實(shí)際資料處理效果評(píng)價(jià)

實(shí)際地震資料為某工區(qū)疊后資料，根據(jù)已有的地質(zhì)、測(cè)井等資料的查證，可以知道該工區(qū)地層的目的層頻帶范圍大體位于8 Hz至70 Hz之間，主頻分布在20 Hz至28 Hz之間。由于研究區(qū)的侏羅系三工河組儲(chǔ)層砂體構(gòu)造整體復(fù)雜，砂層橫向變化大，且整體較薄，同時(shí)受煤層強(qiáng)反射的影響，造成了砂體識(shí)別困難。

圖8(b)為加上隨機(jī)噪聲后的剖面圖，分別應(yīng)用DCT與K-SVD字典進(jìn)行去噪處理，得到圖9 (a) 與圖9 (b)，為了更清晰地看到去噪效果，特輸出如圖9(c)與圖9(d)的殘差圖。

圖8 原始剖面與含噪疊后地震剖面Fig.8 The original profile noisy post-stack seismic section(a)原始剖面；(b)加噪剖面

圖9 去噪圖及殘差圖Fig.9 Result and residual graph(a)DCT字典去噪圖；(b)K-SVD去噪效果圖；(c)DCT殘余圖；(d)K-SVD殘余圖

圖9(a)、圖9(b)與圖8(b)對(duì)比分析，再綜合表2的PNSR(峰值信噪比)可知，K-SVD和DCT去噪效果能力都較好，都能很好地去除大量的隨機(jī)噪聲。殘差圖中只顯示噪聲點(diǎn)，從圖9(c)、圖9(d)中對(duì)比可看到，K-SVD去噪精度更高一些，而且DCT字典對(duì)有用信息有所損害(圖中紅色虛線框)，學(xué)習(xí)字典則對(duì)地震資料的紋理細(xì)節(jié)無(wú)損害，沒(méi)有對(duì)原始信號(hào)品質(zhì)產(chǎn)生破壞。

表2 降噪方法的評(píng)價(jià)指標(biāo)表Tab.2 Noise reduction evaluation index table

從重構(gòu)信息部分來(lái)看，就地震能量破壞程度來(lái)說(shuō)，K-SVD較DCT輕，同時(shí)圖9(b)包含的信息量更為豐富，紅色線框處同相軸對(duì)比更明顯，有效信息破壞程度較之9(a) 更輕，幾乎沒(méi)有破壞原始信息，信噪比高很多。

5 結(jié)論

筆者通過(guò)兩種學(xué)習(xí)字典的去噪與重建效果對(duì)比研究：

1)表明兩種字典都可以有效地對(duì)地震資料進(jìn)行去噪，對(duì)比分析可知，K-SVD結(jié)合了壓縮感知下的匹配追蹤算法的去噪效果要優(yōu)于DCT字典，DCT字典對(duì)地震資料的紋理細(xì)節(jié)有所損害，而K-SVD學(xué)習(xí)字典則幾乎沒(méi)有損害。

2)隨著地層復(fù)雜性的提高，K-SVD學(xué)習(xí)字典能夠根據(jù)實(shí)際資料訓(xùn)練出自適應(yīng)的學(xué)習(xí)字典，更具穩(wěn)定性，去噪效果之后剖面也更加清晰，能得到更多有用的分層信息。

3)字典學(xué)習(xí)可以提高構(gòu)造解釋與儲(chǔ)層預(yù)測(cè)的準(zhǔn)確性，對(duì)于海量地震數(shù)據(jù)的解釋分析工作具有極其重要的意義。