楊海濤, 李 瓊,1b, 錢浩東, 宋 鑫, 雍 鵬

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

0 引言

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

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

1 方法原理

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

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

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

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

y=Cf

(1)

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

f=Ψx

(2)

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

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

y=Cf=CΨx

(3)

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

y=Bx

(4)

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

(5)

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

當滿足RIP條件時，信號重構(gòu)最直接的方法是L0范數(shù)的最優(yōu)化問題：

(6)

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

argmin‖x‖1s.ty=Bx

(7)

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

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

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

1.2 信號質(zhì)量客觀評價方法

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

信噪比：

(8)

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

峰值信噪比：

(9)

其中，MSE為均方差誤差。

1.3 正交匹配追蹤算法

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

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

1.4 DCT字典

DCT 字典是通過離散余弦變換的正交變換來得到的，能較好地對周期信號進行分辨，其特點為能量集中，低頻部分為其主要集中部分，同時可以去掉高頻信號中的無用的一部分，因此適用于對地震資料局部的特征進行處理反變換為式(11)。

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

(10)

(11)

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

Yc=CNy

(12)

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

圖3 DCT字典Fig.3 DCT dictionary

1.5 K-SVD字典學習

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

(13)

圖4 K-SVD學習字典Fig.4 K-SVD learning dictionary

圖5 K-SVD學習字典流程Fig.5 Algorithm flowchart of K-SVD

3 模型測試

為了證明K-SVD與DCT變換矩陣下的OMP算法的去噪效果，筆者選用主頻30 Hz的雷克子波，設(shè)計了如圖6(a)的三層模型的彈炮合成記錄與圖6(b)所示的六層地層模型，圖6(a)設(shè)計為總長為2 560 m，道間距為5 m，波速分別為2 000 m/s、2 400 m/s、2 300 m/s，因此該模型513道接收。圖6(b)設(shè)計為總長為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)是為加噪模型的單道地震合成記錄，對比原始模型可以明顯看到噪音的存在。以兩種字典學習型進行去噪，從而得到如圖7所示的去噪效果圖。對于不同地層，就峰值信噪比來說，表1的PSNR值是由Matlab計算出來的，從表1中可以看出，兩種學習字典的峰值信噪比PSNR都隨著地層的增加而減??；為了更清晰對比，由圖7(e)與7(f)、7(g)與7(h)可知，對同一地層來說，從圖7中紅線框中可以看DCT字典對同相軸的損害要更高一點，可以看到K-SVD字典明顯比DCT字典的效果好，但K-SVD較之DCT，由于在字典的稀疏編碼階段通過奇異值分解不斷地更新字典，在噪聲強度較大的情況下，去噪效果有了顯著的提高，且具有更強的穩(wěn)定性；缺點是由于迭代的次數(shù)比較多，導致運行速度較慢。K-SVD算法中求解系數(shù)矩陣時采用OMP算法，是在MP的基礎(chǔ)上進行施密特正交化，這樣就不用再多次選擇原子，故減少了運算時間。

表1 降噪方法的評價指標表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 實際資料處理效果評價

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

圖8(b)為加上隨機噪聲后的剖面圖，分別應(yīng)用DCT與K-SVD字典進行去噪處理，得到圖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)對比分析，再綜合表2的PNSR(峰值信噪比)可知，K-SVD和DCT去噪效果能力都較好，都能很好地去除大量的隨機噪聲。殘差圖中只顯示噪聲點，從圖9(c)、圖9(d)中對比可看到，K-SVD去噪精度更高一些，而且DCT字典對有用信息有所損害(圖中紅色虛線框)，學習字典則對地震資料的紋理細節(jié)無損害，沒有對原始信號品質(zhì)產(chǎn)生破壞。

表2 降噪方法的評價指標表Tab.2 Noise reduction evaluation index table

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

5 結(jié)論

筆者通過兩種學習字典的去噪與重建效果對比研究：

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

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

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