基于壓縮感知的小波域地震數據實時壓縮與高精度重構

陳祖斌 王麗芝 宋 楊 龍 云*

(①吉林大學儀器科學與電氣工程學院,吉林長春 130026； ②吉林大學地球信息探測儀器教育部重點實驗室,吉林長春130026)

1 引言

數據實時傳輸是制約無纜地震儀發展的最主要因素。隨著石油地球物理勘探技術的飛速發展，地震勘探數據量隨之呈指數增長，給無纜地震儀的無線傳輸、存儲器的容量以及計算機的處理速度帶來了極大的壓力。通過對地震數據進行合理的壓縮處理，可以提高實時處理速度，改善地震儀無線通訊數據傳輸性能[1]。因此，縮減壓縮時間、提高重構精度的需求變得很迫切。

國內外學者在此領域開展了卓有成效的研究。Wang等[2]根據地震資料的特殊性，提出了基于整數小波變換的無損壓縮算法，研究發現CDF(2，n)系列雙正交小波比其他CDF(Cohen Daubechies Feauveau)小波能達到更好的壓縮比。石水勇[3]使用不同提升小波處理地震數據，再分別用EZW(嵌入式零樹小波編碼，Embedded Zero-tree Wavelet)、SPIHT(多級樹集合分裂算法，Set Partitioning in Hierarchical Trees)、EBCOT(優化截斷點的嵌入塊編碼，Embedded Block Coding with Optimized Truncation)三種算法進行壓縮，實驗結果表明db4和bior6.8提升小波變換與SPIHT算法結合的效果最好。Rubin等[4]研究地震數據有損壓縮算法，將CS算法與KRLE(K Run Length Encoding)、LTC(Lightweight Temporal Compression)、WQTR(Wavelet Quantization Thresholding and Run-length)、FFT(Fast Fourier Transform)四種壓縮算法進行多方位比較。Xie等[5]提出HESPIHT(High Efficiency SPIHT)算法，為了驗證其有效性，選擇10個地震剖面做實驗，結果表明,與傳統的SPIHT算法相比，該算法數據恢復精度更高，編解碼效率也更高。徐鋒濤等[6]對地震信號進行二維小波變換，采用EZW算法進行編碼，試驗結果表明，對疊前地震數據壓縮4倍時，信噪比可達50dB以上，而選用適用小波基對疊后地震數據壓縮16倍時，信噪比仍可達30dB以上。Fajardo等[7]在研究二維提升小波壓縮地震數據過程中，實驗發現Bior小波消失距選擇5～8時效果更好，分解層次越多信噪比越低；由于非均勻量化在大幅值時容易出錯，所以在地震數據應用中使用均勻量化獲得最小熵，在實現信噪比高于40dB時，用Huffman編碼比算術編碼更好；因此驗證了壓縮效果與濾波器類型及長度、分解層數、量化方式以及編碼方式有關。以上壓縮方案大多需要對完整采集的地震數據進行操作，而不能對一個或少量幾個采樣數據處理，無法對數據流實時壓縮；并且解碼純粹是編碼的逆運算，壓縮倍數是以損失信號精度為代價，無法做到降低壓縮比的同時，減少重構誤差。

壓縮感知CS(Compressive Sensing)理論認為，若信號是稀疏的或者是在某個變換域內可壓縮，則原始高維信號就能通過測量矩陣轉變成低維信號[8]，有需要時再利用重構算法恢復完整原始信號。目前在地震數據缺道重構領域應用廣泛[9,10]，它與常規壓縮方式最大的不同在于，不是先采集再壓縮，而是邊采集邊壓縮，解碼不受編碼的制約。

為了滿足實時壓縮與高精度重構的需求，本文運用壓縮感知理論[11,12]對無纜地震儀數據采集做壓縮重構處理。通過構造混沌伯努利測量矩陣(CBMM)，對地震數據小波變換后的稀疏系數觀測，得到壓縮數據；再采用貝葉斯小波樹結構壓縮感知(BTSWCS)算法，對觀測值重構，從而恢復完整數據。實驗表明，本文CBMM壓縮所用時間比SPIHT編碼算法縮短很多，特別是待重構信號低信噪比的情況下，BTSWCS算法的重構精度遠高于常見貪婪算法，并且相比SPIHT解碼算法，該算法可以通過改變算法迭代次數而靈活地提高重構精度。

2 地震數據壓縮感知模型

地震數據采集可以表示為如下數學模型[13]

y=Φθ+n

(1)

式中：θ∈RN表示原始地震信號x∈RN在變換域Ψ∈RN×N的稀疏系數，即θ=Ψ∈x，θ中有K個元素不為零， ‖θ‖0=K；Φ∈RM×N表示測量矩陣；n表示測量噪聲，n～N(0,σ2)；y∈RM表示采集的地震數據，由于其中K

式(1)的最稀疏解可以描述為如下問題的解

(2)

由式(1)和式(2)可知，本文基于壓縮感知理論的地震數據壓縮重構技術分三步。

(1)地震信號在時間域不是稀疏的，在某個合適的變換域Ψ對它稀疏表示，得到稀疏系數θ；

(2)構造測量矩陣Φ，用測量矩陣對稀疏系數θ觀測，得到壓縮數據y；

3 方法概述

3.1 地震信號的小波變換

地震信號在時域不是稀疏的，根據CS理論，地震信號首先需要稀疏表示。地震信號是非平穩、非線性信號，而小波變換在分析非平穩信號時優勢很明顯，能自適應地分析各種頻率成分。同時，地震信號在二維空間上相關性很高，具有分頻帶性強和幅度變化大等特征。小波基形態上同地震信號具有相似性，因此對地震信號進行小波變換后，得到的小波系數具有很好的稀疏性(可壓縮性)。無纜地震儀一般是單通道或三通道，對一維信號(單道地震信號)的離散小波多分辨率分析，用Mallet算法描述公式為

(3)

式中：Cj,k、Dj,k分別表示原始信號在尺度2-j下分解得到的低頻分量和高頻分量；h(m-2k)、g(m-2k)分別是低通濾波器和高通濾波器的系數。相應的重建公式為

(4)

對二維地震信號進行分解與重構可以分別對行和列進行一維小波變換得到，其中行號和列號分別對應地震記錄的時間和道號。小波基的選擇直接關系到數據的壓縮效果。對常見小波基的正交性、對稱性、支撐區間和光滑度等性質進行分析，結合文獻[3]得出結論：Biorthogonal系列和Daubechies系列適合用于地震勘探數據的稀疏表示。

3.2 基于Logistic混沌—伯努利測量矩陣的構造

高性能的測量矩陣Φ是數據壓縮測量的關鍵，在測量矩陣的設計方面，矩陣必須滿足RIP(有限等距性質，Restricted Isometry Property)條件[14]，即測量矩陣與稀疏基之間必須是不相關的。常見的隨機測量矩陣很容易滿足此條件，但是除此之外，還應該考慮到測量矩陣的實用性。為推動CS從理論向實際應用的發展，2011年CS理論創始人Candes等[15]提出RIPless理論，放寬了對壓縮測量矩陣的要求。因此，本文對伯努利隨機矩陣做了改進，用Logistic混沌[16]構造滿足RIPless特性的混沌伯努利測量矩陣(Chaotic Bernoulli Measurement Matrix，CBMM)。

混沌序列是根據混沌系統方程迭代產生的，所以兼具確定性和隨機性。混沌映射能夠通過電路實現，因此，本文用混沌序列構成的CBMM具有硬件可實現的潛在優勢。

已知Logistic混沌系統方程如下

(5)

當上式中參數λ∈[1.872,2.0]，初值u0取0.23、0.37和0.7時， Logistic系統進入混沌狀態。用該映射生成混沌序列{un}， 再利用符號函數生成序列{an}，即

(6)

文獻[17]指出當λ=2.0時，Logistic映射產生的混沌序列{an}滿足Bernoulli分布，則可以利用該序列構造測量矩陣。本文構造混沌伯努利測量矩陣的具體步驟為：

(1)根據Logistic混沌映射性質，取參數λ=2.0，初值u0=0.37，生成混沌序列{un}；

(2)為使序列{un}處于穩定狀態，選取{u101}為采樣的起點，并對序列進行等間隔采樣(每隔15個間隔采樣)，確保混沌序列的元素相互獨立。加上初值，得到序列的總長度為M×N；

(3)將第(2)步得到的混沌序列通過式(6)生成序列{an}；

(4)把第(3)步生成的序列{an}取長度為N截斷形成M×N維的矩陣Φ′；

(5)最后對矩陣Φ′作正交化處理，得到Logistic混沌伯努利測量矩陣Φ。

構造完測量矩陣后，需要驗證CBMM是否滿足RIPless特性，只有滿足該特性，才能從觀測值y中準確地恢復出原始信號。RIPless理論認為，測量矩陣Φ中的行向量φk(k=1,2,…,m)是從概率分布F中隨機獨立地抽取的向量，且概率分布F要滿足迷向性和不相干性。

由于本文首先用小波變換對地震信號稀疏表示，即有公式θ=Ψx；再對得到的稀疏系數θ壓縮采樣，即y=Φθ+n。由此可知，本文構造的CBMM用于把高維向量θ轉變為低維向量y。因此，迷向性和不相干性可以通過以下方法證明。

迷向性驗證：因為θ本身具有可壓縮性，所以只需考慮本文的混沌伯努利測量矩陣Φ的行向量φ=(φ[1],φ[2],…,φ[N])的迷向性。也就是檢驗

(7)

CBMM中M個行向量φk∝F，就有E(φφ*)=I，I為單位向量。所以測量矩陣Φ中φ∝F元素的分布滿足迷向特性。

不相干性驗證：假定相干性參數用μ(F)表示，同樣因為θ本身是可壓縮的，那么不相干性表示的就是測量矩陣Φ的行向量的所有元素滿足

MAX1≤t≤N|φ[t]|2≤μ(F)

(8)

因為矩陣元素的取值范圍是[-1,1]，所以相干參數μ(F)=1，則Φ不相干。

綜上所述，本文構造的CMBB滿足RIPless特性，符合CS理論的測量矩陣構造原則，適用于數據壓縮測量。

3.3 貝葉斯小波樹結構壓縮感知重構算法的設計

3.3.1 貝葉斯先驗模型

實際地震勘探數據通常是含噪聲的信號，既含有高斯噪聲，又含有非高斯噪聲。應用地質統計學相關知識，根據式(1)可得到高斯似然函數，即樣本集y的條件概率密度函數。這個高斯分布的均值和方差會隨Φ和θ的變化而變化，是一種分層混合高斯分布

(9)

式中αn表示方差。在CS框架下，已知似然函數表達式，稀疏系數θ的重構問題就能用Bayesian理論進行求解。首先，本文設定θ中每個元素都服從釘鉚先驗分布[22](Spike-and-slab Prior)，即

(10a)

(10b)

(10c)

(10d)

3.3.2 馬爾科夫鏈蒙特卡洛推理

用基于Gibbs采樣的MCMC方法[24]進行式(10)中參數后驗計算。該方法多使用足夠多的樣本數估計后驗分布，通過迭代每個隨機變量(模型參數和中間變量)收集這些樣本，迭代的隨機變量是從它的條件后驗分布給出的所有其他隨機變量的最新值。隨機變量的先驗獨立設為

(11)

因為給出的隨機變量先驗分布與高斯似然函數共軛，所以可以推導出樣本集的條件后驗概率。在每次MCMC迭代，可得出樣本集滿足以下條件后驗分布

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

式中：θs,i是θ中的第i個元素；Φ(j)表示測量矩陣Φ∈M×N的第j列；L(V)是指示函數，如果V是真的，L(V)=1， 若是假的，L(V)=0。為了得到完整后驗的精確表示，MCMC通常要迭代成千上萬次，考慮算法快速收斂，所以每次MCMC迭代，θ按順序抽樣。

4 實驗結果及分析

4.1 壓縮實驗及結果分析

無纜地震儀采樣率為0.25-1ms，所以本文選取采樣點總數分別為214、212、210、28的數據做壓縮實驗。數據來自中國西北地區某次地震勘探。首先對實驗數據進行小波變換，得到稀疏系數，再分別用SPIHT編碼算法和本文構造的CBMM進行壓縮比為0.20、0.25、0.30、0.35、0.40、0.45、0.50、0.55的壓縮處理，記錄壓縮所用時間，因為測量矩陣具有隨機性，每次實驗矩陣元素都不相同，所以結果取50次實驗的平均值，如圖1所示。

圖1a～圖1d很直觀地表明，壓縮相同數據，達到相同壓縮比情況下，本文CBMM要比SPIHT編碼算法所用時間更短。特別地，隨著待壓縮數據量的減少，CBMM壓縮的優勢更加突出。對長度為28的數據用CBMM觀測，所用時間更是縮短至10-5s級；也就是說，若地震儀采樣率是1ms，本文CBMM對無纜地震儀采集的0.25s的數據基本可以實現實時壓縮。

圖1 CBMM和SPIHT算法壓縮所用時間對比

4.2 重構實驗及結果分析

(22)

(23)

圖2 不同壓縮比重構所需的MCMC迭代次數曲線

當選用實際數據采樣點總數為214時，在壓縮比分別為0.20、0.25、0.30、0.35、0.40、0.45、0.50和0.55的情況下，BTSWCS算法精確重構所需的MCMC迭代次數如圖2所示。由圖2可知，當壓縮比為0.2時，200次迭代足以恢復總長度為214的實際數據。

對圖3a所示的實際地震數據，用CBMM進行壓縮比為0.2的觀測，然后分別用貪婪算法和BTSWCS算法重構出完整數據。圖3b是貪婪算法重構結果，圖3c是BTSWCS算法迭代200次重構結果?？梢钥闯?，用BTSWCS算法重構效果明顯比貪婪算法好。貪婪算法對本身含有噪聲和測量噪聲的數據進行重構處理時，重構性能明顯下降，不能達到重構精度要求。

圖4和圖5是對圖3c處理結果的細節評價。圖4a是實際數據小波變換的小波系數圖像，圖4b是BTSWCS算法恢復出的小波系數圖像。圖5a是實際數據第2道地震記錄和與之對應的重構結果的時域圖和頻譜圖，圖5b是第15道地震記錄和與之對應的重構結果的時域圖和頻譜圖。從這些細節圖不難看出，即使是原數據的20%，BTSWCS算法依然可以對其精確重構。

圖3 OMP算法和BTSWCS算法重構實際數據(壓縮比為0.2)

圖4 BTSWCS算法重構與原始數據小波系數比較

圖5 原始數據(上)與重構數據(下)的時域比較(左)和頻域比較(右)

圖6 OMP五種算法性能比較

對圖3a所示地震數據分別用SPIHT編碼算法、CBMM做壓縮比為0.20、0.25、0.30、0.35、0.40、0.45、0.50和0.55的實驗，再通過SPIHT解碼算法、OMP算法和BTSWCS算法恢復出完整數據。不同的壓縮比下，BTSWCS算法統一選擇迭代30次，所有結果取30次實驗結果的平均值。各種算法重構性能如圖6所示。

圖6a是PSNR值曲線，可見對實際數據的處理，OMP算法整體峰值信噪比小于65dB，SPIHT解碼算法和BTSWCS算法重構的PSNR值均大于OMP算法。尤其是BTSWCS算法，當壓縮比大于等于0.25的時候，PSNR值均大于70dB，它重構數據的峰值信噪比比OMP算法至少高出5dB，隨著壓縮比或迭代次數的增加，PSNR值會更高。圖6b是RMSE值曲線，貪婪算法在重構實際數據時，重構誤差均大于0.1。SPIHT解碼算法和BTSWCS算法要比貪婪算法重構誤差小很多。當壓縮比大于等于0.25時，BTSWCS算法的RMSE值均小于0.05，特別是相對于SPIHT解碼算法，BTSWCS算法可以通過改變迭代次數，靈活地改變重構精度。隨著迭代次數的增加，誤差基本趨近于0。

5 結論

地震數據壓縮重構是目前緩解海量地震數據高速傳輸壓力的一項關鍵技術。常規的變換編碼壓縮算法由于本身的局限性，很難達到實時壓縮與高精度解壓。本文根據壓縮感知理論，提出了一種更適用于地震數據壓縮重構的方法。通過對野外實際數據處理得到以下結論和認識。

(1)鑒于地震勘探采集站對實時壓縮的需求，本文通過構造基于Logistic混沌序列—伯努利測量矩陣，對地震數據小波變換后的稀疏系數采集和壓縮同時進行。并且，通過實際數據壓縮發現，CBMM矩陣測量要比SPIHT編碼算法壓縮時間縮短很多。對采樣點為28的數據用CBMM矩陣觀測，所用壓縮時間更是縮短至1.0×10-5s。也就是說，若地震儀采樣率是1ms，CBMM測量矩陣對它采集的0.25s的數據基本可以實現實時壓縮。

(2)鑒于常用壓縮感知重構算法在低信噪比時失效，本文提出BTSWCS算法提高重構精度。該算法根據小波系數內部結構統計特性，先設置小波系數服從釘鉚先驗模型，再用MCMC方法對模型參數進行后驗估計，從而恢復完整數據。對壓縮的實際數據恢復發現，BTSWCS算法重構效果明顯優于常用貪婪算法，重構結果峰值信噪比能達到70dB以上，重構誤差能小于0.05，PSNR值至少提升5dB。