Radon域下的K--SVD字典的混采分離

李慧, 韓立國, 張良,賈帥

吉林大學地球探測科學與技術學院，長春130026

0 引言

傳統的海上或者陸地地震數據采集中，為避免鄰近炮源之間的干擾常常使用連續震源，但是連續震源激發需要震源間隔大，這就產生了時間長、資金消耗多的缺點。針對該問題，在近十年，提出了同步源采集的思想，它在減少采集費用及改善地震數據質量上有很大的突破[1]。同步源采集是多個炮點同時激發的方法，它實現了用更少的時間間隔來密集采樣的要求,在同一時間內收集到多個炮點的混合數據,減少采集時間，增加了空間采樣，但地震數據分析與處理目前需要傳統采集的地震單炮數據，當前已相繼用Seislet變換、多級中值濾波和單純字典學習等來研究去除混采噪聲，并取得一定成就。稀疏表示是一種高效的信號處理方法，為此2008年Berkhout[2]在此基礎上進行擴展，提出混合采集的概念，闡述一種在鄰近炮點之間、同一時間內的連續重疊方法，這在減少采集時間和增加空間采樣方面有巨大的優勢。

目前大部分的地震處理工作和疊前屬性分析仍然依靠傳統的單炮采樣數據，需要將同步源采集的數據進行分離，便于地震數據分析及處理，即混采分離,以降低采集時間和資金花費。混合數據的處理是在混合數據中直接移除和分離數據：在相同時間內，相對于某一個震源，其他震源在隨機延遲時間內激發，此時其他震源激發的能量類似于一個不同步噪聲，幾個震源的能量被記錄在一個常規接收的時間域內，由于延遲的隨機性，兩個源之間的串擾能量類似于不同步的相干噪聲，那么該問題也就轉化為去噪的問題。Ikelle[3]在2007年指出，混采分離可以視為盲源分離問題，可以用獨立成分分析算法解決。隨后一些方法相繼提出來去除典型噪聲：Doulgeris和Huo et al.[4,5]提出基于相干濾波的算法，該方法計算效率較高而且得到了預期的結果；Mansour et al.[6]將混采分離作為一個反演問題來處理，然后優化這個方案的方法；韓立國等[7]提出多級中值濾波與Curvelet變換相結合的方法來去除混采噪聲。

在壓制相干噪聲的方法中，二維面內所有點的函數值都與原始函數線積分有映射關系，Radon變換根據這個原理，沿同相軸進行積分求和使形狀規則的同相軸的能量在Radon域中收斂，根據同相軸的截距和曲率不同使其在Radon域內分離，以此來去除噪聲。

在地震應用中，字典學習主要被用來去除隨機噪聲，基于稀疏表示的字典在地震信號處理中效果顯著，可以應用于地震去噪和一維信號重建[8]。然而由于初始字典的選擇困難，用同步源采集和訓練數據在實際處理中難以執行。除此之外，在假設字典有特殊結構的前提下，一個經過訓練的緊框架能用于地震插值。字典的學習能用來設計數據驅動字典，可以更好地訓練數據、改善稀疏信號[9]。字典學習的具體方法是：假設的參考信號通過字典原子線性地稀疏表示，用設計的最優化方法來學習字典。相對于固定變換基的字典，一個學習的字典更符合信號，也能更稀疏地表示信號。

筆者結合Radon域下的字典學習原理，用K--SVD的稀疏字典的學習方法，將數據轉換到Radon域進行字典學習來分離混合數據，用合成數據和實際地震資料進行測試，相較于中值濾波、單純字典和小波變換，該方法具有更好的分離效果。

1 方法原理

1.1 Radon變換

時間域Radon變換組順著同相軸進行線積分，轉換到Radon域收斂成一個點或一個線(圖1)，這一過程有明顯的物理意義—將時間--空間域t-x的一條直線映射到τ-p域。一般地，Radon變換可以表示為:

(1)

該公式描述了一個二維無限連續方程d(t,h)的Radon正變換。

若從Radon域返回到t-h域則需進行反Radon變換，其定義為：

(2)

需要注意的是，這個方程并不是正演方程的精確反演，兩個方程合成正反變換對。正反變換是由Radon變換對的算法實現的。

在地震勘探中，時間域和空間域的采樣是離散的，不能使用連續的函數方程，故用離散的累加來替代連續的積分方程。同時Radon變換有拋物、線性和雙曲線3種主要變換類型。雙曲Radon變換(Hyperbolic Radon Transform)在地震資料的去噪、插值重建、速度估計和多次波衰減等方面應用廣泛。雙曲Radon變換可收斂類似于雙曲線的反射波同相軸，從而使反射波與其他波分離，傳統的雙曲Radon變換是時變的，Bickel[10]改進離散雙曲Radon變換方程為：

(3)

時不變的雙曲Radon變換和傳統的時變雙曲Radon變換相比，具有可以在頻率--空間域進行計算的優勢，極大提高了雙曲Radon正變換的計算效率，本文采用時不變的雙曲積分路徑。

1.2 混合數據的偽分離

在隨機時間延遲假設下的同步源采集中，野外采集到的混合地震信號b用矩陣形式可以表示為：

b=Γf

(4)

式中：Γ為混合矩陣；f為常規采集的單炮記錄，在已知b中恢復未知的f。這個式子是欠定的，所以我們可以求它的最小平方解：

Γ-1=[ΓHΓB]-1ΓH

(5)

在稀疏編碼中進行線性組合，那么偽分離信號f′為：

f′=ΓHb

(6)

偽分離數據在不同時域呈現不同的分布[7]，從圖1中可以看出，偽分離數據在共檢波點道集中，混疊噪聲為隨機分布，所以筆者將偽分離數據轉換到共檢波點域進行分離去噪。

圖1 不同時域的偽分離數據Fig.1 Pseudo-separation data in different domains

1.3 字典分塊的炮分離

基于共檢波點道集的炮分離可以視為一個反演問題。通過字典學習的稀疏表示主要有兩步：

第一步稀疏編碼。信號f的恢復是復雜的正則化過程，在隨機同步源采集的猜想下，f的恢復可以利用壓縮感知和稀疏反演的公式(7)：

(7)

式中：D為預定義的字典；m為對應的稀疏系數向量；T為非零系數的數量。D與m按列對應，字典中的每一列按照Mi(稀疏系第i行)為系數進行線性組合就會得到單炮記錄F。

第二步更新字典。思想是用獲得的系數稀疏向量m，更新字典D，經過K次迭代完成字典的一次更新：

(8)

公式(7)和(8)是用來學習最優化字典和最稀疏的表示的迭代方法[11]。

同步源地震數據的字典首先從混合數據b中學習，但是由于混合記錄中有很大的噪聲，所以本文的字典從單炮記錄中學習。訓練字典時，單炮記錄中每一塊的數據被選擇，等式變為

(9)

(10)

式(9)是NP-hard問題[12]，只能用窮舉法且不容易得到，本文用的是正交匹配追蹤(OMP)算法解決，是貪婪的逐個窮舉來找到最精確的解。迭代更新字典D的(10)式本文用了K奇異值分解(K--SVD)算法。

1.4 用K--SVD算法學習字典

K--SVD方法是解決等式(14)的一個有效地算法。在已知F和M的情況下，字典中的第K列在固定稀疏系數M的情況下一次次被更新，為更新第K列，等式變為:

(11)

E=U∑VT

(12)

式中：這里U中的第一列為dk，V的第一列和Σ的第一個元素乘積為稀疏系數向量mk。在K列全部被更新之后，再用OMP算法解決第一個等式求出稀疏矩陣M，即可求出分離后的單炮數據：

f=dm

(13)

本文采用K--SVD與Radon變換相結合的雙稀疏字典，在有效去除混疊噪聲的基礎上提升速度。

2 測試與應用

2.1 合成數據處理

資料處理時采用最小偏移距為20 m，道間距和炮間距均為2 m,每道采樣點數226個，采樣間隔1 ms的模型數據，常規采集得到90個單炮記錄，之后對單炮數值混合成45個超級炮，混合度為2的混疊炮記錄(圖2)。

圖2 模擬數據的混采Fig.2 Blended acquisition of simulation data

圖3是混疊數據和偽分離數據在不同時域的偽分離記錄，可以明顯看出，在共檢波點域將偽分離數據變換為有效同相軸和隨機噪聲，為之后的去噪工作提供條件；圖4是原炮記錄和偽分離記錄分別在Radon域中的表示，可見炮記錄在Radon域中能量收斂，而噪聲則分散分布在Radon域中，在這里用閾值濾波進行分離去噪。分離后的數據(圖5)，雖然仍有一些殘余噪聲，但是分離效果較明顯，接下來將用字典進行分離殘余的噪聲。可以看出，分離后的數據會稍有損失，這是使用Radon濾波的缺點。

圖3 不同時域的偽分離數據Fig.3 Pseudo-separation data in different domains

圖6a為Radon域的炮記錄訓練成的字典，字典原子總數為64個，圖像塊尺寸大小為8×8。圖6b是經過字典學習后的數據，將Radon域的數據轉換到時間域，得到圖7a，將Radon閾值未去除的噪聲去除得更加干凈。圖7b為分離后的數據和單炮記錄的差值，可以看出有效波有損失的，但是完成了混采分離的目標。圖8是單純使用中值濾波或者字典、二層小波進行處理，可以看出中值濾波的有效波損耗比較少，但是濾波效果較差，能量較高的隨機噪聲殘留太多；字典處理仍有部分損耗，而且去除效果差，比中值濾波更差；小波處理效果最差，基本沒有去掉較大的能量。相比較而言，用本文方法的處理結果比較干凈，有一定優勢。

2.2 實際數據處理

繼續處理混合度為2、常規采集為90個單炮、90個道的實際地震記錄(圖9)，以進一步驗證方法的有效性。數據道間距25 m,炮間距25 m,每道采樣點數900個，采樣間隔為1 ms。圖10為Radon域的數據，用雙曲Radon變換，將單炮數據變換到Radon域。

圖4 Radon域中的偽分離記錄Fig.4 Pseudo-separation data in Radon domain

圖5 Radon域分離去噪Fig.5 Denoising in Radon domain

圖6 分塊的字典和字典分離后的結果Fig.6 Blocking dictionary and result of dictionary separation

圖7 最終的分離結果及去除的噪聲Fig.7 Results of separation and noises attenuated

圖8 單純中值濾波或字典或小波的處理結果Fig.8 Results of median filter or dictionary or wavelets

圖9 實際數據的混采Fig.9 Field data by blended acquisition

圖10 Radon域數據Fig.10 Data of Radon domain

圖11 Radon域分離去噪Fig.11 Denoising in Radon domain

圖12 訓練字典進行分離Fig.12 Separation by training dictionary

圖13 單純中值濾波處理和字典、小波處理結果Fig.13 Results of median filter or dictionary or wavelets

將數據轉換到共檢波點域使混疊數據變為隨機分布的噪聲(圖9)，對偽分離信號進行Radon變換(圖10d)，可以看出由于噪聲較大，壓制了有效波的顯示。在圖11中可以看出Radon去噪的有效性，濾掉的部分是圖11d，同時也去掉了一些Radon域內能量較小的有效波，這個方面還需后續的改進；再進行K--SVD字典去噪，將Radon域內有效波進行分塊字典處理，字典原子依然是64個，圖像塊尺寸大小為8×8，完成最終分離過程。圖12c是最終的分離結果，雖然有效波損失但是證明了方法的有效性。圖13是用中值或字典、二層小波來進行數據處理，可以看出中值濾波效果最差，噪聲殘留明顯，字典處理較好，二層小波去掉的主要是小能量噪聲，實際效果不大，都沒有本文方法的處理結果干凈。

3 結論

(1)同步源采集具有省時、省資金和有效的特點，是未來采集方法的發展方向。在Radon域中進行閾值濾波和K--SVD字典學習，滿足了偽分離數據分離為常規單炮數據的要求。

(2)從模擬數據和實際數據的處理結果中可知，該方法分離效果明顯，同中值濾波、單純字典和小波變換相比，雖然會有少量振幅損耗但能更有效地去噪，達到分離的目的，在實際處理中具有實用價值。