FK和Shearlet域聯合壓縮感知數據重構技術

閆海洋 周 輝 劉海波 徐朝紅 孫贊東 劉 昭

(①中國石油大學(北京)油氣資源與探測國家重點實驗室，北京 102249； ②中國石油大學(北京)CNPC物探重點實驗室，北京 102249； ③中國石油大學(北京)地球物理學院，北京 102249； ④東方地球物理公司海洋物探處，天津 300457)

0 引言

稀疏炮檢點采集或野外采集因素造成地震數據的不規則影響后續地震資料成像質量?；趬嚎s感知理論的重構方法能夠有效恢復地震數據。壓縮感知數據重構技術是將重構問題轉化為稀疏優化問題。由于空間隨機缺失的地震道引起的空間假頻在稀疏域呈白噪分布，有效信號在稀疏域集中分布，可采用反演迭代的方法逐步恢復缺失的地震數據。

Ma等[1]實現了基于Fourier變換的快速迭代閾值收縮算法，在Fourier域重構地震數據； Zhang等[2]在FK域數據重構過程中引入了凸集投影(POCS)算法，在迭代過程中引入了指數衰減閾值參數。Fourier變換是全局變換，不能描述地震數據的方向特征，且稀疏度不夠，數據重構精度有待提高，于是能夠描述地震信號方向特征的曲波變換應用到壓縮感知數據重構中。劉國昌等[3]和Yang等[4]實現了曲波域的POCS方法地震數據重構，采用指數閾值提高了迭代收斂的速度； Zhang等[5]提出了一種基于曲波變換的三維數據同步重構和噪聲壓制方法； Zhang等[6]提出了一種基于非等距快速離散曲波變換的三維重構方法，采用線性Bregman迭代提高反演效率； Cai等[7]在曲波域基于L1范數的譜梯度投影(SPG-L1)算法實現了數據重構； Liu等[8]提出了局部隨機采樣的曲波域數據重構； Wang等[9]提出了一種基于曲波變換和高分辨率Radon變換的多域稀疏約束高精度地震數據恢復方法。

曲波變換采用網格旋轉方式實現方向表征，無法實現對線奇異在連續系統和離散系統的統一，而Shearlet變換克服了這個缺點，采用剪切方式代替曲波變換的旋轉方式，能更加稀疏表達地震信號。馮飛等[10]提出了采用Jitter欠采樣的Shearlet變換稀疏約束地震數據重構方法； 張良等[11]實現了基于壓縮感知的Shearlet變換地震數據重構，其重構效果優于Fourier變換、離散余弦變換、小波變換和曲波變換； 劉成明等[12]采用Shearlet變換與基于Landweber加速下降迭代法進行插值，在提高計算效率的同時保證了插值精度； 楊冠雨等[13]采用基于Shearlet變換的雙正則化方法進行地震數據重構，同時兼顧了信號的稀疏性和地質結構的復雜性； 王常波[14]在Shearlet域利用快速凸集投影(FPOCS)算法和指數閾值進行壓縮感知地震數據重構，與Fourier變換、曲波變換的重構結果對比表明，Shearlet變換重構精度最高。

與固定稀疏基不同，形態分量分析(Morphological Component Analysis，MCA)框架采用不同基函數線性組合方式對地震數據稀疏表達，通過迭代方式逐一求解并重建各形態分量，最后通過合并重建后的形態分量完成數據的重構。張凱等[15-16]提出了一種Shearlet字典和離散余弦變換(DCT)字典組合的地震數據重構方法； 周亞同等[17]提出在MCA框架下的基于DCT與曲波字典組合的地震信號重構方法。MCA框架的壓縮感知數據重構需要迭代求解、重建各形態分量，計算效率有待提升。

本文采用Shearlet變換對地震數據稀疏表示，與常規對有缺失的地震數據進行Sheartlet變換不同，是將時空域的地震道重構轉化為FK域的隨機噪聲壓制問題。首先對有缺失的地震數據做FK變換，然后對FK域數據做Shearlet變換，通過閾值迭代法消除FK域由于地震道缺失引起的空間假頻，從而實現地震道的時空域重構。與形態分量分析(MCA)框架下的地震數據重構不同，MCA認為信號由多個不同形態的分量線性組合而成，用不同的稀疏字典組合方式進行地震數據重構，本文采用的 FK變換加Shearlet變換，可以看作一種新的稀疏基變換，不需要迭代求解各形態分量。

稀疏采樣方法是基于壓縮感知觀測系統設計的重要內容。全局隨機采樣空間假頻呈白噪分布，分段隨機采樣空間假頻具有藍色噪聲的頻譜特征，在保持采樣隨機性的同時，有效解決了重采樣結果疏密不勻問題，提高了地震數據的重構效果。最后通過Shearlet域和全局隨機采樣、Shearlet域和分段隨機采樣、FK+Shearlet域和全局隨機采樣、FK+Shearlet域和分段隨機采樣四種數據重構結果對比，表明FK+Shearlet域和分段隨機采樣結合重構效果最好，信噪比和保真度最高。

1 方法原理

1.1 FK+Shearlet域聯合重構方法原理

考慮到傳統采樣定理的不足及信號的稀疏特性，Candes等[18]、Donoho[19]提出了壓縮感知理論。壓縮感知數據重構包括三個方面：①信號的稀疏表述，信號是稀疏的或在變換域內是稀疏的； ②測量矩陣的設計，測量矩陣要與稀疏基完全不相干； ③重構算法的構造，信號的有效恢復依賴高效高精度重構優化算法。

基于壓縮感知的地震數據重構可表示為求解方程

d=Mf

(1)

式中：f∈Rp為完整數據；d∈Rq為觀測數據；M為對角采樣矩陣，對角線上的元素與d對應，采樣前、后數據的維數分別為p和q，且p>q。當d中地震道數據未缺失時，對應的M對角線上的元素為1； 地震道數據缺失時，對應的M對角線上的元素為0。M一般選取高斯隨機矩陣。

時空域的地震數據不是稀疏的，但可以通過數學變換滿足稀疏特性。張良等[11]通過實驗認為基于Shearlet變換的數據重構精度高于Fourier變換、DCT變換、小波變換和曲波變換。

本文利用FK+Shearlet變換對數據進行稀疏表示，將時空域不連續的地震數據重構問題轉為連續FK域數據隨機噪聲壓制問題，稀疏基為Shearlet基函數

x=φΨf

(2)

式中：Ψ為FK正變換；φ為Shearlet正變換；x為f的FK+Shearlet變換結果。

將式(2)代入式(1)，可得

d=MΨ-1φ-1x

(3)

式中：Ψ-1為FK逆變換；φ-1為Shearlet逆變換。令

θ=MΨ-1φ-1

(4)

為感知矩陣。

式(3)為欠定問題，無法直接求解。感知矩陣θ只有滿足有限等距性質(RIP)時，則可以求取唯一解[18-20]。高斯隨機矩陣與稀疏基不相干，θ滿足RIP條件。地震數據重構問題可以表示為

(5)

式(5)在一定的條件下，L0范數最小化求解問題可以轉化為L1范數最小化問題[22]，即

(6)

設ε為重構誤差，式(6)可寫為

(7)

關于高精度重構方法，國內外學者對Bregman迭代[6,21]、POCS[14，22-23]、SPG-L1[7，24-25]、正交匹配追蹤法[26-27]、交替乘子方向算法(ADMM)[13，28]等方法進行研究或改進。溫睿等[29]通過實驗對比分析了POCS、IHT(Iterative Hard Thresholding)、Bregman和JRSI(Joint Reconstruction by Sparsity-promoting Inversion)四種迭代算法，認為POCS和IHT算法較為穩定，經過多次迭代能夠達到較高的信噪比。本文采用POCS進行壓縮感知數據重構，算法如下

(8)

(9)

1.2 稀疏采樣

稀疏采樣方法是壓縮感知地震采集研究的重要內容，決定了由于地震數據缺失引起的空間假頻在稀疏域的隨機分布特征。相比全局隨機采樣，Jitter稀疏采樣在保持隨機性的同時，分布更均勻，減少了假頻和有效信號的混疊，重建效果優于全局隨機采樣[10，12，30]。本文采用王漢闖等[31]提出的分段隨機采樣方法，該方法不僅具有更為優越的藍譜特征，還突破了Jitter采樣對段內樣點為奇數的限制，具有更強的實用性。

以一道為例分析重采樣因子對數據頻譜(波數譜)的改造。對完整數據進行稀疏采樣，就是對重采樣因子與完整數據進行點乘，即

d(t)=r(t)f(t)

(10)

時間(空間)域的乘積對應頻率(波數)域的褶積，即

D(ω)=R(ω)*F(ω)

(11)

式中：f為單道完整地震數據；r為重采樣因子；d為采樣結果；F、R和D分別為f、r和d的Fourier變換。完全采樣時，r(t)≡1。

圖1為規則、全局隨機、分段隨機1/2采樣示意圖。規則采樣時采樣因子r中0和1規則交替出現； 1/2全局隨機采樣的采樣因子r中0和1全局隨機出現； 1/2分段隨機采樣中2個采樣點為1段，每段采樣因子r中0和1隨機出現。假定原始數據采樣間隔為2ms，采樣點數為3000，則奈奎斯特頻率為250Hz。圖2a為完全采樣因子及其頻譜，0頻為有效信號頻譜，頻譜值為3000，不存在假頻。圖2b為1/2規則采樣因子及其頻譜，0頻處為有效信號頻譜，頻譜值降為原來的1/2(1500)，假頻在-250Hz處，與0頻處有效信號頻譜值相等。圖2c為1/2全局隨機采樣因子及其頻譜，0頻處譜值降為原來的1/2(1500)，假頻的頻譜呈白噪特征，假頻和有效信號混疊。圖2d為1/2分段隨機采樣因子及其頻譜，0頻處譜值降為原來的1/2(1500)，假頻的頻譜呈藍譜特征，假頻和有效信號混疊。

圖1 規則(a)、全局隨機(b)、分段隨機(c)1/2采樣示意圖實心表示采樣； 空心表示未采樣

圖3a和圖3b分別為完全采樣和1/2規則采樣的數據及其頻譜。規則采樣后，如果信號頻率和假頻混疊，則無法保真恢復原始信號。圖3c和圖3d分別為全局1/2隨機采樣和分段1/2隨機采樣的數據和頻譜。隨機采樣后，假頻振幅值相對較小，有效數據的頻譜大于假頻頻譜，方便后續通過反演迭代進行全頻帶數據恢復，理論上可以突破奈奎斯特頻率的限制。但分段隨機采樣數據的藍譜特征明顯，由于地震數據為帶限數據，分段隨機更有利于后期數據重構。

圖2 不同稀疏采樣因子(左)及其頻譜(右)(a)完全采樣因子； (b)1/2規則采樣因子； (c)1/2全局隨機采樣因子； (d)1/2分段隨機采樣因子

1.3 FK+Shearlet變換

信號的稀疏表示就是將信號投影到正交變換基時，有效信號對應的變換系數集中分布，系數相對較大，噪聲(對地震數據重構來說，主要指由于數據缺失造成的空間假頻)對應的變換系數的絕對值很小。Shearlet變換[32-34]可以有效地表示地震數據各向異性特征。由剪切方式代替了曲波變換的旋轉方式，能夠在離散系統和連續系統做統一處理。

對于任意二維平方可積函數b(τ)，其Shearlet變換可表示為

SHφb{j,l,m}=〈b,φj,l,m〉

(12)

式中

(13)

(14)

其中：A為各向異性膨脹矩陣，控制尺度；B為剪切矩陣，控制方向，A和B均為2階可逆矩陣；j為尺度參數；l為剪切參數；m為平移參數。

Shearlet變換是一類多尺度幾何分析方法，通過對基函數的縮放、剪切和平移等仿射變換生成具有

圖3 不同稀疏采樣因子采樣后的數據(左)及其頻譜(右)(a)完全采樣因子； (b)1/2規則采樣因子； (c)全局1/2隨機采樣因子； (d)分段1/2隨機采樣因子

不同特征的剪切波函數，具有較好的地震數據的方向性描述和稀疏表達能力。

應用模擬單炮數據(采樣間隔為2ms，512個樣點，512道)對比Shearlet變換與FK+Shearlet變換效果。采用1尺度的Shearlet分解，將模擬單炮分解成1個不帶方向特征的低頻分量和8個帶方向特征的高頻分量。圖4為模擬單炮及FK譜，圖5為對模擬單炮進行Shearlet分解的9個系數，由于分解的8個高頻系數(圖5a～圖5h)能量較弱，對其振幅放大10倍顯示。

圖6為圖5的9個Shearlet域系數對應的時空域數據。圖7為對模擬單炮進行FK+Shearlet分解的9個系數譜，圖8為對模擬單炮進行FK+Shear-let分解的9個系數對應的時空域數據。表1對比了模擬單炮不同稀疏基變換分解系數能量占比。對模擬單炮直接進行Shearlet變換，不帶方向特征的低頻分量系數能量占比為98.75%，8個帶方向特征的高頻分量系數能量占比僅為1.25%； 對模擬單炮進行FK+Shearlet變換，不帶方向特征的低頻能量占比為8.81%，8個帶方向特征的高頻分量系數能量占比為91.19%，FK+Shearlet域的帶有方向特征的高頻系數相對Shearlet域能量更強，能更好地觀測地震數據的方向(線奇異)特征。

圖4 模擬單炮(a)及其FK譜(b)

圖5 模擬單炮Shearlet分解的9個系數(a)～(h)帶方向特征的8個高頻系數譜，振幅放大10倍顯示； (i)不帶方向特征的低頻系數譜

圖6 圖5的9個系數對應的時空域數據(a)～(h)帶方向特征的8個高頻系數譜，振幅放大10倍顯示； (i)不帶方向特征的低頻系數譜

圖7 模擬單炮分解的FK+Shearlet域 9個系數譜(a)～(h)帶方向特征的8個高頻系數譜； (i)不帶方向特征的低頻系數譜

圖8 圖7的9個系數譜對應的時空域數據(a)～(h)帶方向特征的8個高頻系數譜； (i)不帶方向特征的低頻系數譜

表1 模擬單炮不同變換分解的系數能量占比

2 模擬數據實驗

采用重構信噪比評價重構效果

(15)

式中f*表示重構的地震數據。

應用上述模擬單炮記錄對比全局隨機和分段隨機采樣250道(采樣率48.8%)在Shearlet域與FK+Shearlet域的數據重構效果。表2 為模擬單炮重構信噪比統計。

圖9為模擬單炮記錄的全局隨機采樣結果及其FK譜。圖10為Shearlet域和FK+Shearlet域數據重構結果及誤差。Shearlet域重構數據(圖10a)的信噪比為13.39dB，FK+Shearlet域重構數據(圖10c)的信噪比為21.05 dB； Shearlet域數據重構的誤差(圖10b)明顯大于FK+Shearlet域數據重構誤差(圖10d)，表明FK+Shearlet域重構效果優于Shearlet域。全局隨機采樣數據FK譜(圖9b)的空間假頻呈白噪分布，FK+Shearlet域重構數據有效信號的FK譜(圖11b)能量的恢復及空間假頻的壓制效果明顯優于Shearlet域重構數據的FK譜(圖11a)。

圖9 模擬單炮全局隨機采樣結果(a)及其FK譜(b)

圖10 模擬單炮全局隨機采樣Shearlet域與FK+Shearlet域模擬單炮重構效果分析(a)Shearlet域重構結果； (b)Shearlet域重構誤差； (c)FK+Shearlet域重構結果； (d)FK+Shearlet域重構誤差

圖11 模擬單炮全局隨機采樣不同域重構結果的FK譜(a)Shearlet域； (b)FK+Shearlet域

圖12為模擬單炮記錄的分段隨機采樣結果及其FK譜。圖13為Shearlet域和FK+Shearlet域數據重構結果及誤差。Shearlet域重構數據(圖13a)的信噪比為13.89dB，FK+Shearlet域重構數據(圖13c)的信噪比為23.52dB，Shearlet域數據重構的誤差(圖13b)大于FK+Shearlet域數據重構誤差(圖13d)，FK+Shearlet域重構效果優于Shearlet域。分段隨機采樣數據的FK譜(圖12b)的空間假頻呈藍譜特征，空間假頻與有效信號混疊相對全局隨機(圖9b)要弱。FK+Shearlet域重構數據有效信號的FK譜(圖14b)的能量恢復及空間假頻的壓制效果優于Shearlet域重構數據的FK譜(圖14a)。從數據重構信噪比(表2)可見，Shearlet域分段隨機重構效果(圖13a)優于全局隨機重構(圖10a)，FK+Shearlet域分段隨機重構效果(圖13c)優于全局隨機重構(圖10c)。

通過不同稀疏采樣比例進一步驗證FK+Shearlet的重構效果。圖15為分段隨機采樣下不同稀疏采樣比例和不同稀疏基變換模擬單炮重構信噪比對比。整體上，FK+Shearlet比Shearlet域重構信噪比高。當稀疏采樣比例低于60%時，FK+Shearlet比Shearlet域重構信噪比高6dB以上； 當采樣比例高于70%時，隨著采樣比例的增加，兩種稀疏基重構效果逐漸接近。

圖12 模擬單炮分段隨機采樣結果(a)及其FK譜(b)

圖13 模擬單炮分段隨機采樣Shearlet域與FK+Shearlet域重構效果分析(a)Shearlet域重構結果； (b)Shearlet域重構誤差； (c)FK+Shearlet域重構結果； (d)FK+Shearlet域重構誤差

圖14 模擬單炮分段隨機采樣不同域重構結果的FK譜(a)Shearlet域； (b)FK+Shearlet域

表2 模擬單炮不同域重構信噪比統計

圖15 不同稀疏采樣比例FK+Shearlet(紅色)與Shearlet域(藍色)重構數據信噪比對比

3 實際數據應用

選擇海域地震資料(圖16a)進行壓縮感知數據重構。由于海洋資料頻帶較寬，同時原始數據也存在空間采樣不足引起的空間假頻，給數據重構帶來一定難度，以70%全局隨機采樣(圖16b)和分段隨機采樣(圖16c)對比Shearlet域與FK+Shearlet域數據重構效果。圖17為原始、全局隨機采樣、分段隨機采樣數據的FK譜，分段隨機采樣數據FK譜(圖17c)的空間假頻呈藍譜特征，空間假頻與有效信號混疊相對全局隨機(圖17b)要弱。表3為不同稀疏采樣和稀疏基變換實際數據重構信噪比對比。

圖16 實際原始地震數據(a)及其全局(b)和分段隨機采樣結果(c)

圖17 實際原始地震數據(a)及其全局(b)和分段隨機采樣結果的FK譜(c)

在全局隨機采樣情況下，Shearlet域重構數據(圖18a)的信噪比為10.08dB，FK+Shearlet域重構數據(圖18b)的信噪比為15.66dB，Shearlet域數據重構的誤差(圖18c)明顯大于FK+Shearlet域數據重構誤差(圖18d)，FK+Shearlet域重構效果優于Shearlet域重構效果。在分段隨機采樣情況下，Shearlet域重構數據(圖19a)的信噪比為10.47dB，FK+Shearlet域重構數據(圖19b)的信噪比為16.52dB，Shearlet域數據重構的誤差(圖19c)大于FK+Shearlet域數據重構誤差(圖19d)，FK+Shearlet域重構效果優于Shearlet域重構。

在全局隨機采樣情況下，FK+Shearlet域重構數據有效信號的FK譜(圖20b)能量的恢復及空間假頻的壓制效果明顯優于Shearlet域重構數據的FK譜(圖20a)。在分段隨機采樣情況下，FK+Shearlet域重構數據有效信號的FK譜(圖20d)能量恢復及空間假頻的壓制效果優于Shearlet域重構數據(圖20c)。

圖18 實際數據全局隨機采樣下不同稀疏域重構結果對比(a)Shearlet域重構結果； (b)FK+Shearlet域重構結果； (c)Shearlet域重構誤差； (d)FK+Shearlet域重構誤差

圖19 實際數據分段隨機采樣下不同稀疏域重構結果對比(a)Shearlet域重構結果； (b)FK+Shearlet域重構結果； (c)Shearlet域重構誤差； (d)FK+Shearlet域重構誤差

圖20 實際數據兩種采樣方式、不同域重構結果的FK譜對比(a)全局隨機，Shearlet域； (b)全局隨機，FK+Shearlet域； (c)分段隨機，Shearlet域； (d)分段隨機，FK+Shearlet域

由表3的數據重構信噪比可見，Shearlet域分段隨機重構效果(圖19a)優于全局隨機重構(圖18a)，FK+Shearlet域分段隨機重構效果(圖19b)優于全局隨機重構(圖18b)。

表3 實際數據的重構信噪比統計

4 結論

壓縮感知數據重構技術要求地震數據隨機稀疏分布，由于數據缺失引起的假頻在稀疏域表現為噪聲，有效信號集中分布，假頻和有效信號混疊，但是假頻相對較弱，通過迭代的方法能進行數據重構。分段隨機采樣方法具有更為優越的“藍色噪聲”頻譜分布，分段隨機采樣在保持采樣隨機性的同時，有效解決了重采樣結果疏密不均問題。本文采用FK+Shearlet變換作為稀疏基進行分段隨機采樣數據重構。實驗表明，FK+Shearlet域比Shearlet域能更好地描述地震數據的方向(線奇異)特征，分段隨機采樣FK+Shearlet域重構精度高于全局隨機采樣Shearlet域重構、分段隨機采樣Shearlet域重構和全局隨機采樣FK+Shearlet域重構。