一種新的K-SVD 字典學習地震數據去噪方法

周子翔，吳娟，袁成，白敏*，桂志先

（1. 油氣資源與勘探技術教育部重點實驗室（長江大學），湖北武漢 430100；2. 長江大學地球物理與石油資源學院，湖北武漢 430100；3. 中國石油勘探開發(fā)研究院西北分院，甘肅蘭州 730020）

0 引言

隨著地震勘探開發(fā)難度的不斷增大，對地震資料質量的要求也越來越高。在現場采集地震數據時，由于地表等因素的影響，不可避免地會引入噪聲，降低了地震資料信噪比，因此地震數據去噪就顯得非常重要。地震數據去噪處理的眾多方法可分為幾大類：基于濾波的方法［1-2］、基于變換的方法［3-8］、基于降秩類的方法［9-10］和基于學習的方法等。

近年來，基于學習的方法發(fā)展迅速，它是學習含噪圖像到“干凈”圖像的一種映射關系，可細分為基于字典學習的方法［11-14］和基于深度學習的方法［15-16］?；谧值鋵W習的方法旨在得到含原始數據信息的自適應字典，需對數據進行多次學習訓練，根據有效信號與隨機噪聲在字典上表現的不同特點，突顯數據本質特征，壓制隨機噪聲，達到去噪目的?；谏疃葘W習的方法是一種數據驅動類方法，去噪時需從大量數據樣本學習得到有效信號與噪聲的區(qū)別，搭建適用的深度神經網絡以壓制噪聲，包括準備數據、搭建網絡、訓練網絡、最優(yōu)模型等環(huán)節(jié)，各環(huán)節(jié)的處理經驗對最終訓練效果有顯著影響。本文主要探討字典學習方法在地震數據去噪中的應用。

字典學習源于壓縮感知理論［17-18］。在稀疏表示中，若待處理的數據是稀疏的，就可通過某種稀疏變換稀疏表示原始地震信號，從而實現數據去噪。但由于地震信號的稀疏表示能力因采用不同的變換方法表達效果會有明顯差異，或因實際采集的地震數據具有不同的局部特性，這時若僅使用某一種變換所組成的稀疏字典往往不能充分、有效地將其全部表達，將對去噪結果產生一定影響。因此，基于數據本身特征，自適應地構造一種新的適用學習型方法［19］，更有效地應用于數據去噪顯得尤為重要。

Olshausen 等［20］率先提出學習型過完備字典概念，并將該理論應用于圖像去噪領域；在此基礎上，Aharon 等［21］提出一種新的自適應字典算法，實現信號的稀疏表示。將K 奇異值分解（K-SVD）作為字典更新方法，自適應地獲取表達信號本身特征的原子，訓練出學習型過完備字典，將學習訓練得到的字典作為圖像稀疏表示的稀疏字典，以有效壓制隨機噪聲。之后，基于K-SVD 算法的字典學習廣泛地應用于地震數據去噪［22］。但這種傳統的K-SVD 字典學習方法在去噪精度及計算效率等方面仍有進一步提升的空間。

針對以上問題，本文提出一種新的用于地震數據去噪的K-SVD 字典學習方法：采用一種將任意位置的塊作為學習的對象，并對這種不定的稀疏表示模型進行學習算法的改進；然后對稀疏編碼算法關鍵步驟進行優(yōu)化，使新的K-SVD 字典學習方法在信噪比、計算時間和數據局部特征保持等方面更具優(yōu)勢。文中首先闡述了字典學習方法基本理論；然后基于傳統KSVD 字典學習方法，進行改進并給出新方法的去噪原理及處理流程；最后，通過數值模擬試驗，分析、評估本文方法的去噪性能。

1 理論方法

先介紹字典學習的基本理論及去噪模型，再重點闡述新方法的K-SVD 字典更新理論。在此基礎上，使用新的字典學習方法對地震數據進行去噪。

1.1 字典學習

字典學習應用于地震去噪領域有兩個關鍵步驟：一是對地震數據進行稀疏編碼；二是學習更新字典。

對于含噪地震數據模型，可將其近似看作

式中Q0是待估測的無噪地震數據，通常所采集的是含有隨機噪聲ε的有噪數據Q。去除隨機噪聲相當于從采集的含噪地震數據Q中近似得到無噪數據Q0，簡言之，即是利用一種方法壓制噪聲ε。

1.1.1 稀疏編碼

基于K-SVD 算法的字典學習，是設計出一種基于小塊的過完備字典算法，對數據樣本進行研究；通過對每一小塊進行去噪，使其更接近理想化結果。字典學習包含的主要參數有可稀疏線性表示［23］的原始信號樣本Q（m×n，m為樣本個數，n為每個樣本的特征維度）、包含k個原子的字典矩陣D（m×k）、及稀疏矩陣X（k×n），這樣就可理想化表示為Q=DX。該思想過程可表示為

式中：xi(i=1，2，???，k)是稀疏矩陣X的行向量，其元素是字典矩陣中對應原子的系數；T為稀疏度；F 是指Frobenius范數。

實際上，在壓縮感知與稀疏表示理論框架中，利用拉格朗日乘子法就可將上述去噪問題轉變成求解無約束的正則化問題

式中：第一項是數據真實值與誤差之間的控制項，即控制重建后地震信號與含噪信號之間的誤差；第二項是xi的P范數，它是為了保證系數足夠稀疏。由此得到理想化的既稀疏又能準確地表達地震數據的稀疏解。當P=0 時，即需找出系數的0 范數，通常這是一個NP 問題［24］，相對來講較難解決，但若利用其默認滿足的有限等距條件，即可將0 范數替換為1 范數，從而簡化計算以便于求解

式中正則化參數λ可確定目標函數中數據精確性和系數稀疏性的權重，通過設置和調節(jié)λ值可獲得更佳方案。

針對稀疏編碼系數，通?？捎谜黄ヅ渥粉櫍∣rthogonal Matching Pursuit，OMP）這一貪婪算法近似求解。當然還有很多類似算法可求取近似解，如基追蹤（Basis Pursuit，BP）［25］、Bregman 算法等。

1.1.2 字典更新

通過K-SVD 算法進行字典更新時，因難以同時求解D和X，故K-SVD 算法采用迭代思想，在每次迭代中交替求解D和X。字典更新是固定稀疏編碼階段的X，從而求解出字典D，在多次迭代后得到最優(yōu)稀疏系數X及對應的字典D。此階段的約束條件［26-27］為

K-SVD 算法采用逐列更新字典方式，如對字典D的第t列做更新時，記D的第t列列向量為dt，并稱其為原子，記X的第t行行向量為xt，更新過程如下

利用奇異值分解方法分解誤差矩陣，使得

由此求解出最優(yōu)的dt、xt′。dt即取U的第1 列更新字典的第t列，xt′即取Σ的對角線上的第1個元素與V的第1 列的乘積，由此進行更新，并替換xt，這樣便使式（6）的值為最小。其中

上述式中：U是的特征向量構成的矩陣，稱作左奇異矩陣（m×m）；V是的特征向量構成的矩陣，也稱作右奇異矩陣（n×n）；Σ是Λ1（m×m）和Λ2（n×n）兩個矩陣的對角線上的非零元素由大到小構成的對角矩陣（m×n）。

1.2 地震數據去噪

新的K-SVD 字典學習方法對地震數據進行去噪包括兩個階段：首先對樣本數據進行訓練，得到學習字典；然后利用學習到的字典對含噪數據做降噪處理，得到較理想的最終去噪結果。整個過程通過多次循環(huán)實現。

在訓練集足夠大的情況下，本文從含噪數據中訓練提取任意位置的塊作為訓練集，使得用K-SVD 算法得到的字典更能較準確地捕獲數據中的特征，提高分塊后含噪數據的去噪質量。該算法對字典進行訓練時需從樣本數據中選取大量圖像塊作為訓練集，隨著訓練樣本的增加，求取系數時的計算量也會增加。雖然訓練集的數量是保證稀疏表示后圖像精度的前提，但數量過多的訓練集會加重運算負擔，降低算法效率。因此，根據所采集地震數據選擇合適的塊的大小尤為重要。在后面數值試驗部分，本文方法設置原子的尺寸為4×4 的塊，學習字典的尺寸［28］為16×64的塊，即字典大小設置為原子的4倍，即一個字典包含維度為16的64個原子。

在得到學習字典之后，重要的是對從含噪數據中提取的重疊部分的每個塊進行去噪。利用稀疏編碼參數中的約束條件，通過學習的字典計算每個塊的稀疏逼近進行去噪，然后對去噪后的塊進行平均以重建整個數據，完成對地震數據的去噪。這里使用OMP進行稀疏編碼，塊向量x與其稀疏逼近?之間的誤差需滿足下列約束條件［29］

式中：通過絕對中位差（Median Absolute Deviation，MAD）假設噪聲標準差為σ，在處理噪聲數據時，也可根據噪聲方差估計正則化參數λ值；C為增益因子，取1.15［21］。增益因子是據經驗估計的，目的是使稀疏編碼步長更靈活；同時，也是數據精確度和稀疏性的權重參數。增益因子越大，字典中塊的表示就越稀疏，精度就越低。

2 處理流程

圖1a 是傳統K-SVD 字典學習算法［22］去噪流程圖，本文主要從初始化字典、求稀疏系數的OMP 算法、字典的學習過程等方面進行優(yōu)化。其中從含噪數據中選取任意位置的塊作為字典學習的對象，這種隨機性使得能更準確地捕獲地震數據特征，進而實現不定的稀疏表示，以此提高K-SVD 字典學習的去噪效果。同時，在OMP 算法中增加了稀疏表示每個信號允許的最大表示誤差的約束，使得OMP 算法的效率大幅度提升等。

圖1b 是本文所提K-SVD 字典學習算法實現地震數據去噪的流程圖。具體去噪過程是：首先設置適用參數，如K-SVD 算法的迭代次數、字典學習過程的迭代次數、稀疏度等；然后從樣本訓練數據中提取隨機位置的塊，并訓練刪除空白塊，初始化字典D，根據設置的迭代次數進行字典學習，包括用OMP 算法稀疏編碼得到稀疏系數X，用K-SVD 算法更新字典直至遍歷所有原子，在達到設定的全局迭代次數時完成字典學習，同時得到最新字典D和最新稀疏編碼系數X；最后，在含噪數據中構建新塊，利用學習到的最新字典D對構建的新塊使用OMP 算法得到稀疏系數X，由字典D和新稀疏系數X求出近似塊，對該近似塊插入平均值，再重構圖像塊得到去噪后結果，完成地震數據去噪。

3 數值試驗

為了驗證本文方法去噪性能，分別利用理論數據和現場采集數據做數值試驗，將其去噪結果與f-x反褶積方法和傳統K-SVD 字典學習方法［22］進行比較。

為定量評價去噪效果，將信噪比定義為原始信號的功率（Power of Signal）與噪聲的功率（Power of Noise）之比值取對數，再乘以10［1］

式中Ps、Pn分別為不含噪聲原始地震數據和噪聲數據的功率。

3.1 理論數據

采用圖2 所示雙曲同相軸合成理論數據進行試驗。圖2a 是原始無噪數據，大小為256×256，即共有256 道，每道256 個采樣點，時間采樣間隔為4 ms。對圖2a加入標準差為0.2的隨機噪聲（圖2c），得到圖2b所示含噪數據。在本文理論數據研究部分，隨機噪聲均通過randn 函數產生。因為randn 函數產生的隨機數服從均值為0、標準差為1 的標準正態(tài)分布（高斯分布），所以加入的隨機噪聲也是高斯噪聲。

圖2 原始數據（a）、含噪數據（b）及加入的隨機噪聲（c）

分別采用f-x反褶積、傳統K-SVD 字典學習及本文方法對含噪數據做去噪處理（圖3）。從f-x反褶積方法去噪結果（圖3a）中，可見噪聲殘留仍相對較多；而去除的噪聲中仍包含較多有效信號（圖3b）。傳統K-SVD 字典學習方法去噪效果（圖3c）比f-x反褶積方法的好，但仍會損失部分原始信號（圖3d）。與前面兩種方法相比，本文方法獲得更好去噪結果（圖3e），在去除的噪聲剖面（圖3f）上基本看不到殘留有效信號，表明它更好地保存了原始數據特征。原始含噪數據信噪比為-2.59 dB，利用上述三種方法去噪之后的信噪比分別是3.10（f-x反褶積）、9.16（傳統K-SVD 字典學習）和10.46 dB（本文方法）。對比可知，傳統K-SVD 字典學習算法對信噪比有顯著提升，而本文方法相對又進一步提高了信噪比。試驗過程中，三種方法的耗時分別是2.76（f-x反褶積），271.67（傳統K-SVD 字典學習）和34.89 s（本文方法）。由此可知，本文方法既具有更高信噪比，計算效率也比傳統K-SVD 字典學習方法顯著提升。

字典學習即是先給定一個初始字典，然后通過學習更新字典。本文方法的字典是通過在含噪數據中隨機選擇64 個訓練塊開展隨機初始化。圖4a、圖4c分別是兩種字典學習方法的初始化字典原子，得到相應的更新后64 個字典原子（圖4b、圖4d）?？梢姼虑暗淖值湓虞^單一，存在表現特征完全一樣的字典原子，圖4a 中的原子形狀是剛性的，不能最優(yōu)地表示高度非平穩(wěn)的地震數據，圖4c的原子則較具柔性。兩種方法更新后的字典原子均顯得更豐富，且字典原子的表現特征也不同，字典中的每個原子都被重塑，相比于傳統K-SVD 字典學習方法，本文方法更新的字典在每個原子的色域上更寬泛，說明本文方法可更充分地突顯地震數據局部特性。

數值試驗中，另一較重要的參數是迭代次數。學習過程中迭代次數通常設置為100 次，但實際上大多不需迭代100 次即能獲得較滿意結果，因此本次試驗設定的學習迭代次數為50，可相應地縮短整個學習和去噪過程的時長。

為了對比不同噪聲水平下前述三種方法的去噪效果，利用randn 函數對原始理論數據分別加入標準差為0.1、0.2、0.3、0.4、0.5 的隨機噪聲，測試結果如表1、圖5a和圖5b所示。

表1 不同噪聲水平下去噪效果和時間對比

圖5 三種方法去噪后的信噪比（a）及去噪處理時間（b）的對比

從表1 對比數據和圖5a 可看出，在去噪效果上，不同方法去噪后的信噪比隨噪聲強度的逐漸增加都在降低，本文方法的去噪效果好于傳統K-SVD 字典學習方法和f-x反褶積方法。

從三種方法的去噪時間（表1 和圖5b）看，f-x反褶積方法的去噪時間隨加入噪聲的強弱變化較小，整體趨于穩(wěn)定。兩種K-SVD 字典學習算法的去噪時間在加入較弱～較強的隨機噪聲的過程中，在開始階段有一個跳躍的變化，然后隨著加入的噪聲越來越強，去噪時間的增幅并不大。與傳統K-SVD 字典學習方法相比，本文方法在提高精度的同時也大幅度節(jié)省了時間。

圖6 為原始數據的f-k譜。圖7 為圖3 數據的f-k譜對比，其中圖7a、圖7c 和圖7e 分別是上述三種方法去噪結果對應的f-k譜：可見f-x反褶積方法在去噪的同時也去除了很多有效信號，有效波損傷較大（圖7b）；而兩種K-SVD 方法去噪結果的fk譜均與原始數據f-k譜差別不大，進一步對比這兩種去噪結果頻譜與原始信號頻譜之差（圖7d 與圖7f），可知本文方法的去噪殘差中有效信號損失相對更少。

圖6 原始數據的f-k 譜

圖7 三種方法去噪結果的f-k 譜對比分析

3.2 實際采集數據

圖8 是實際采集的地震數據，受采集條件及現場環(huán)境因素影響，與理論數據相比，實際采集數據更復雜，噪聲的存在淹沒了部分有效信號。

圖8 實際地震數據

應用上述三種方法進行去噪，從所得處理結果（圖9a、圖9c、圖9e）可見：f-x反褶積方法雖可壓制部分噪聲能量，但它損失了較多有效信號（圖9b）；傳統K-SVD 字典學習方法去噪效果（圖9c）雖比f-x反褶積方法的略好，但噪聲剖面（圖9d）中也有部分有效信號殘留，去噪后數據同相軸連續(xù)性較差；相比而言，本文方法有效壓制了絕大部分噪聲能量，去噪相對徹底（圖9e），去噪后數據同相軸較連續(xù)，幾乎未損失有效信號（圖9f），去噪效果更好。

圖9 三種方法去噪結果（左）及其去除的噪聲（右）

利用這三種方法對實際地震數據去噪后的信噪比分別是6.04（f-x反褶積）、8.79（傳統K-SVD 字典學習）和9.54 dB（本文方法）?？梢姳疚姆椒ㄈピ牒蟮男旁氡茸罡?，表明該方法在較完整保留有效信號的同時，還顯著提升了信噪比。三種方法對本次試驗數據實施去噪分別耗時為2.53（f-x反褶積）、47.26（傳統K-SVD 字典學習）和4.81 s（本文方法）。從該數據可看出，本文方法在提高信噪比的同時，相比于傳統K-SVD 字典學習去噪方法還大幅度地提高了計算效率。

信噪比是針對全局的度量，它不能展示去噪性能較差的區(qū)域；而局部相似圖［22］是針對去噪性能的一個局部度量，能更直觀地顯示去噪結果與去除噪聲之間的局部相似性，可更準確地評估信號泄漏的程度。

為進一步展示去噪性能，下面采用局部相似圖分析這三種方法對有效信號的保護效果。若出現一些高相似性的異常區(qū)域，即表明去噪結果在相應位置上存在較多有效信號泄漏。

圖10 展示三種方法去噪后結果數據與去除的噪聲之間的局部相似圖。由局部相似圖原理可知，能量強的部分指示對原始信號的損失較大。從該圖中可見：f-x反褶積方法的局部相似圖（圖10a）中存在較強能量，表明去噪后損失了很多原始信號；傳統K-SVD方法的局部相似圖（圖10b）中能量弱于f-x反褶積方法，但有效信號仍殘留較多；本文方法的局部相似圖（圖10c）中能量最弱，表明原始信號損失最少，更完整地保留了有效信號局部特征。

圖10 三種方法去噪結果與去除噪聲之間的局部相似性

4 結論與展望

本文提出一種用于地震數據去噪的新的K-SVD字典學習方法，數值試驗顯示它具有良好的去噪性能。該方法是通過從含噪數據中提取隨機位置的塊實現字典初始化；并用OMP算法稀疏編碼和K-SVD算法字典更新進行學習，構建出自適應的字典，使原子特征逐漸接近樣本數據；再用字典對含噪地震數據分塊進行去噪。與f-x反褶積和傳統K-SVD 字典學習方法相比，本文方法去噪效果更好，信噪比更高，在提升計算效率和保持地震數據局部特征方面更具優(yōu)勢。

雖然K-SVD 是字典學習算法中很有效的算法之一，但不論哪一種方法，總會有其欠完善之處。如計算成本還較高，需數以千計的奇異值分解，在三維或五維地震數據中可能有更突出的體現，這也阻礙了其在地震數據處理中的廣泛應用。對此，擬從三方面做后續(xù)深入研究：①探究適合求解大型稀疏矩陣的0范數或1范數問題，使稀疏表示能力更強，更能突出地震數據的信息；②在字典學習階段可嘗試利用在線字典學習算法中的最小角回歸（Least Angle Regression，LARS）算法或改進的OMP 算法，使學習得到的字典符合地震數據的結構特征和反映地震數據的局部信息變化規(guī)律；③一些初始參數的設定可考慮根據經驗或實驗數據的不同稍做優(yōu)化處理，可能會有更好結果。這樣可望改進字典學習算法在地震去噪領域的應用效果，也為實現工業(yè)化應用提供更可靠的理論基礎。

文中未與深度學習去噪方法對比。當前深度學習面臨的問題是如何加入物理約束，并賦予其更明確的物理意義。因此，后續(xù)研究方向將是發(fā)展字典學習約束的深度學習去噪方法，以獲得更佳去噪效果。