唐國(guó)維， 程鑫華， 張巖

(東北石油大學(xué) 計(jì)算機(jī)與信息技術(shù)學(xué)院， 大慶 163318)

0 引言

現(xiàn)在的高分辨率地震勘探技術(shù)對(duì)地震數(shù)據(jù)的規(guī)則性以及完整性提出了更高的要求。然而由于客觀地質(zhì)條件、勘探成本以及其他環(huán)境因素，在地震勘探時(shí)，無(wú)效的地震道采樣、壞道、以及各種不規(guī)則的干擾數(shù)據(jù)通常都會(huì)對(duì)地震數(shù)據(jù)進(jìn)行一定程度的干擾，使得野外勘探得到的數(shù)據(jù)往往是不規(guī)則和不完整的。若不對(duì)缺失的地震數(shù)據(jù)進(jìn)行處理，會(huì)影響多次波壓制、偏移等處理結(jié)果，并最終影響地震資料的解釋和油氣藏的評(píng)價(jià)。在實(shí)際勘探過(guò)程中，處理不規(guī)則采樣數(shù)據(jù)的一些相對(duì)簡(jiǎn)單的方法，如拷貝或線性插值出鄰近道、忽略空缺道等做法均有不足之處。因此，尋求一個(gè)快速而有效的地震數(shù)據(jù)重建方法，對(duì)于地震數(shù)據(jù)的分析以及后續(xù)工作具有重大意義。

地震數(shù)據(jù)重建，即利用一定的方法和手段對(duì)已有的道缺失地震數(shù)據(jù)進(jìn)行插值重建處理，生成完整的或者采樣率更高的地震數(shù)據(jù)。地震數(shù)據(jù)的重建方法大致一共有3類：第一類是基于濾波器的策略：將欠采樣地震數(shù)據(jù)與某種插值濾波器做卷積，常用的一種是預(yù)測(cè)誤差濾波器，但這類方法通常將隨機(jī)采樣數(shù)據(jù)當(dāng)作規(guī)則數(shù)據(jù)來(lái)處理，并使用高斯窗進(jìn)行插值操作，這種操作有著非常大的不確定性，會(huì)造成很大的誤差；第二類是基于波動(dòng)方程的方法：通過(guò)DMO或AMO正、反演算子迭代進(jìn)行插值處理，此種方法能夠處理隨機(jī)采樣的地震數(shù)據(jù)，但是計(jì)算量巨大，并且需要有地下結(jié)構(gòu)的先驗(yàn)信息，對(duì)于使用較粗網(wǎng)格的欠采樣地震數(shù)據(jù)的重建效果并不理想；第三類是基于變換函數(shù)的方法：利用特定的變換函數(shù)對(duì)欠采樣地震數(shù)據(jù)進(jìn)行變換，在變換域內(nèi)進(jìn)行插值等處理之后再進(jìn)行逆變換，最終得到重建后的地震數(shù)據(jù)。

在過(guò)去幾十年里，奈奎斯特(Nyquist)采樣定理在采樣、存儲(chǔ)、傳輸和處理等信號(hào)處理的傳統(tǒng)領(lǐng)域中起到了非常關(guān)鍵的作用。最近發(fā)展出的壓縮感知(Compressed Sensing,CS)理論表明[1,2]，若待處理的數(shù)據(jù)是稀疏的，通過(guò)一定的采樣方法，就可以突破奈奎斯特定理，利用更少的數(shù)據(jù)來(lái)重建出滿足一定精度要求的重建數(shù)據(jù)。基于以上特點(diǎn)，將壓縮感知理論應(yīng)用于地震數(shù)據(jù)重建[3,4]:國(guó)內(nèi)已有白蘭淑提出了一種基于曲波變換的聯(lián)合迭代重建算法[5]，周亞同提出了一種基于K-奇異值分解的字典訓(xùn)練的地震數(shù)據(jù)重建方法[6]，唐剛則使用泊松碟采樣的方式對(duì)地震數(shù)據(jù)進(jìn)行重建等。因此壓縮感知的稀疏表示部分通常是先將欠采樣的地震數(shù)據(jù)進(jìn)行某種變換來(lái)滿足壓縮感知理論中對(duì)于數(shù)據(jù)稀疏性的要求。通常使用的變換方法有離散余弦變換(DCT)、傅立葉變換、離散小波變換等。離散小波變換作為信號(hào)和數(shù)據(jù)處理工具取得了巨大的成果，離散小波變換由于其時(shí)頻局部化特性、多分辨率特征性、邊緣檢測(cè)特性等使得其在地震數(shù)據(jù)重建方面有著重要應(yīng)用。但是傳統(tǒng)的離散小波具有的自身局限性，主要體現(xiàn)在：缺乏平移不變性；缺乏方向選擇性；震蕩性；頻譜混疊性。普通的復(fù)小波變換可以解決上述局限性的問(wèn)題。但是由于超過(guò)一層分解的復(fù)小波變換的輸入都是復(fù)數(shù)形式，所以很難找到與之對(duì)應(yīng)的完全重構(gòu)的濾波器。為了解決這個(gè)問(wèn)題，Nick Kingsbury等提出了雙樹復(fù)小波變換(DTCWT)[7]，它按照一定的規(guī)則采用雙樹濾波的形式設(shè)計(jì)，既保留了一般復(fù)小波的優(yōu)點(diǎn)，又可以完全重建，能夠更有效的處理地震數(shù)據(jù)。

2 地震數(shù)據(jù)稀疏表示

2.1 壓縮感知理論

壓縮感知(Compressed Sensing,CS)理論表明，如果采集到的原始數(shù)據(jù)具有稀疏屬性，同時(shí)使用一定的采樣方法，就可以通過(guò)少量的欠采樣數(shù)據(jù)來(lái)重建出具有一定精度的接近原數(shù)據(jù)的地震數(shù)據(jù)。從以上描述可以得知，壓縮感知理論已在地震數(shù)據(jù)重建領(lǐng)域占有了一席之地。但對(duì)于絕大部分的地震圖像來(lái)說(shuō)，稀疏性不是一個(gè)共通的屬性，所以需要對(duì)地震圖像的稀疏表示，即進(jìn)行某種正交變換，使其滿足壓縮感知理論中的稀疏特性。欠采樣的地震數(shù)據(jù)重建過(guò)程可以表示為模型為式(1)。

y=Rf

(1)

來(lái)表述。式中:N維向量f表示原始的完整地震數(shù)據(jù)；R是一個(gè)M×N階的采樣矩陣，未采集的點(diǎn)的位置為0，其余采集點(diǎn)的值為1；N維向量y表示向量化后的不完整數(shù)據(jù)。由于采集數(shù)據(jù)的不完整性，R是一個(gè)欠定矩陣，整個(gè)地震數(shù)據(jù)重建的過(guò)程，就是由采樣矩陣R和不完整的地震數(shù)據(jù)y重建出完整的、達(dá)到一定精度的地震數(shù)據(jù)f的過(guò)程。

根據(jù)壓縮感知理論，如使用某種變換D∈RNXL使得向量a=DHf是稀疏的，則稀疏過(guò)程可以表示為式(2)。

y=RDa=Aa

(2)

其中A=RD為感知矩陣。因?yàn)橄蛄縜本身具有稀疏性，則原始重建問(wèn)題轉(zhuǎn)化為求解最小l0范式問(wèn)題為式(3)。

min‖a‖0y=RDa

(3)

由于實(shí)際操作中允許誤差的存在，則求解式 的最優(yōu)化問(wèn)題可以用一個(gè)最簡(jiǎn)單的近似求解形式代替為式(4)。

min‖a‖0‖y-RDa‖≤ε

(4)

其中ε為一個(gè)非常小的常量。l0范數(shù)最小使得稀疏非零元素的個(gè)數(shù)最少，而K

2.2 雙樹復(fù)小波變換

普通的復(fù)小波變換可以解決離散小波的局限性問(wèn)題。但是由于一層分解以上的復(fù)小波變換的輸入都是復(fù)數(shù)形式，所以很難找到與之對(duì)應(yīng)的完全重建的濾波器。Nick Kingsbury等提出了雙樹復(fù)小波(DTCWT)解決了這個(gè)問(wèn)題，它使用了雙樹濾波的形式進(jìn)行設(shè)計(jì)，既保留了一般復(fù)小波的優(yōu)點(diǎn)，又可以完全重建[6]。復(fù)小波可以表示為式(5)。

ψ(t) =ψr(t) +jψi(t)

(5)

ψr(t)，ψi(t)分別表示復(fù)小波的實(shí)部和虛部，他們都是實(shí)函數(shù)，這樣雙樹復(fù)小波變換可以表示為兩個(gè)獨(dú)立的實(shí)小波變換，一個(gè)給出實(shí)部，一個(gè)給出虛部，如圖1所示。

圖1 雙樹復(fù)小波變換示意圖

樹a和樹b為一對(duì)平行樹，樹a給出雙樹復(fù)小波變換的實(shí)部，樹b給出雙樹復(fù)小波變換的虛部。雙樹復(fù)小波變換的優(yōu)良性能來(lái)源于其解析性，為使ψ(t)滿足解析性，則需要ψr(t)與ψi(t)組成Hilbert變換對(duì)，兩個(gè)正小波函數(shù)組成Hilbert對(duì)的充要條件是：兩低通濾波器滿足半采樣延遲條件，確保樹b中的第一層向下采樣取到樹a中因隔點(diǎn)采樣而舍棄的，未保留采樣值。所以雙樹復(fù)小波具備了頻譜單邊性的優(yōu)良特性，同時(shí)在二抽樣條件下，讓雙樹復(fù)小波具有其自身優(yōu)勢(shì)：

1) 平移不變性：雙樹復(fù)小波變換具有平移不變性，即信號(hào)的微小平移不會(huì)導(dǎo)致在各尺度上能量的變化。

2) 多方向選擇性：雙樹復(fù)小波不單融合了離散小波所具有的良好視頻特性，同時(shí)還有更好的方向分析手段。

3) 具有完全重構(gòu)特性，地震圖像數(shù)據(jù)在分解后可以保證完全重建，保證了數(shù)據(jù)的重建效果。

4) 數(shù)據(jù)冗余較為有限。對(duì)于一維信號(hào)其冗余為2∶1；對(duì)于二維信號(hào)冗余為4∶1。

5) 較少的計(jì)算量。分解重建過(guò)程的計(jì)算量比非抽象離散小波變換都要少很多。

基于以上所述的雙樹復(fù)小波所具有的優(yōu)良特性，能夠在地震數(shù)據(jù)重建領(lǐng)域有著很好的應(yīng)用。

3 重建算法

雙樹復(fù)小波變換后，利用閾值軟迭代算法(IST)[7]，進(jìn)行地震數(shù)據(jù)的重建，同時(shí)考慮到相同尺度內(nèi)的子帶維數(shù)均相同，同位置的小波系數(shù)在同一尺度內(nèi)的各個(gè)子帶的位置是固定的，這樣有利于分析同一位置小波系數(shù)間的統(tǒng)計(jì)特性。用合適的分布擬合小波系數(shù)分布，并通過(guò)分析系數(shù)間的相關(guān)性特征，建立合適的統(tǒng)計(jì)模型，可提高重建算法的精度。本文利用地震信號(hào)的雙樹復(fù)小波變換系數(shù)中同一方向子帶與當(dāng)前子帶的父帶小波系數(shù)之間的相關(guān)性，利用雙閾值軟迭代法[8,11]進(jìn)行地震數(shù)據(jù)的最后重建步驟。

用a1、a2表示對(duì)采樣前完整地震數(shù)據(jù)進(jìn)行稀疏變換得來(lái)的小波系數(shù)及其父系數(shù)，u1、u2表示有缺失道地震數(shù)據(jù)稀疏變換后的小波系數(shù)及其父系數(shù)，ε1、ε2表示由于道缺失對(duì)地震數(shù)據(jù)造成的假頻影響的小波變換系數(shù)及其父系數(shù)，則有式(6)。

u1=a1+ε1

(6)

雙閾值法考慮的是當(dāng)前系數(shù)與其父系數(shù)之間的關(guān)系，則有式(7)。

u2=a2+ε2

(7)

則上述式子可以改寫為式(8)。

u=a+ε

(8)

(9)

根據(jù)條件概率公式進(jìn)行推導(dǎo)，上式可變?yōu)槭?10)。

(10)

在迭代重建過(guò)程中，注意到地震數(shù)據(jù)的道缺失具有隨機(jī)性，則對(duì)于道缺失對(duì)于雙樹復(fù)小波變換后地震數(shù)據(jù)的影響ε，此處設(shè)ε服從均值為0，方差為σε高斯隨機(jī)分布，即式(11)。

(11)

設(shè)原地震數(shù)據(jù)在同一方向的小波系數(shù)與其父系數(shù)的雙變量聯(lián)合概率密度函數(shù)為式(12)。

(12)

其中ɑ是待定參數(shù)，σ2表示被估計(jì)信號(hào)小波系數(shù)的方差。

則將上述式子代入，通過(guò)計(jì)算得到當(dāng)前系數(shù)ɑ1的MAP估計(jì)為式(13)。

(13)

在以上計(jì)算推導(dǎo)過(guò)程中，為了計(jì)算得到當(dāng)前系數(shù)的估計(jì)值 ，必須先得到模型中的未知參數(shù) 與被估計(jì)數(shù)據(jù)的方差 。在實(shí)際應(yīng)用中，由于地震數(shù)據(jù)的卻是造成的信號(hào)干擾的方差一般是未知的，通常需要進(jìn)行估計(jì)。由于道缺失對(duì)于原始地震數(shù)據(jù)造成的影響主要集中在高頻子帶，且在各個(gè)高頻子帶中不盡相同且無(wú)法準(zhǔn)確估計(jì)，此處，我們可以利用小波系數(shù)來(lái)進(jìn)行中值估計(jì) (Donoho,1994)[12]為式(14)。

(14)

其中表示第i個(gè)高頻子帶的雙樹復(fù)小波分解系數(shù)。

對(duì)于道缺失的地震數(shù)據(jù)的小波系數(shù)的方差 、被估計(jì)小波系數(shù)的方差 和假頻小波系數(shù)方差 ，假設(shè)這三者之間滿足關(guān)系式為式(15)。

(15)

對(duì)于方差 的估計(jì)值計(jì)算公式為式(16)

(16)

其中N是鄰域U的大小。利用上式可以得到 的估計(jì)值為式(17)。

(17)

其中max表示原地震圖像數(shù)據(jù)中的最大值，對(duì)于灰度圖像的地震數(shù)據(jù)即255、MSE表示原圖像與處理圖像之間的均方誤差。PSNR值越大，表示與原地震數(shù)據(jù)的差異越小。

雙閾值收縮法插值重建由以下幾個(gè)步驟組成：

1) 對(duì)欠采樣地震數(shù)據(jù)進(jìn)行雙樹復(fù)小波變換。利用正交小波變換的快速算法獲得低分辨率下的尺度系數(shù)以及各個(gè)分辨率下的小波系數(shù)，其中尺度系數(shù)和小波系數(shù)共N個(gè)。

2) 對(duì)小波系數(shù)進(jìn)行非線性閾值處理。此處利用不同頻率間子帶與父帶的相關(guān)性，使用以上描述的雙閾值算法進(jìn)行重建，對(duì)分解過(guò)程中的低頻系數(shù)不做處理。硬閾值雖然可以保留更多的信號(hào)特征，但是在平滑方面還有所欠缺。根據(jù)本文情況選用雙閾值軟收縮法。

3) 進(jìn)行逆向雙樹復(fù)小波變換。由所有低頻尺度系數(shù)、以及經(jīng)由閾值處理后的小波系數(shù)做逆變換進(jìn)行重建，提取欠采樣的地震道，用插值法與欠采樣地震數(shù)據(jù)插值重建為較為接近原地震數(shù)據(jù)的結(jié)果。

本文使用6層分解的雙樹復(fù)小波變換，需要逐層處理高頻分解系數(shù)，每次插值對(duì)每層小波系數(shù)都進(jìn)行雙閾值處理，則整體的算法復(fù)雜度為o(N)。

4 實(shí)驗(yàn)結(jié)果及分析

實(shí)驗(yàn)的硬件平臺(tái)采用雙核CPU主頻3.3G的Intel I5微機(jī)，內(nèi)存容量為4G。系統(tǒng)軟件為32位Windows7操作系統(tǒng)，仿真實(shí)驗(yàn)軟件使用Matlab R2013b。實(shí)驗(yàn)數(shù)據(jù)采用廣泛使用的marmousi模型，此模型由法國(guó)石油研究院于1998年做出，這種模型及其聲波有限差分合成數(shù)據(jù)已被全世界成百上千的科研人員用于許多地球物理科研項(xiàng)目，直到今天仍然是出版最多的地球物理資料集之一。地震數(shù)據(jù)重建效果的衡量指標(biāo)采用峰值信噪比(PSNR)，如式(18)。

(18)

其中MSE為不含噪聲原始地震數(shù)據(jù)與去噪后地震數(shù)據(jù)的均方誤差，定義為式(19)。

(19)

為了對(duì)比本文算法的性能，對(duì)采樣地震數(shù)據(jù)分別進(jìn)行傳統(tǒng)的迭代閾值法與文中上述雙閾值插值重建法進(jìn)行重建。整個(gè)過(guò)程經(jīng)過(guò)200次迭代，如圖2—圖7所示。

圖2 迭代算法示意

圖3 原始地震數(shù)據(jù)

圖4 添加隨機(jī)道缺失之后的地震數(shù)據(jù)，采樣率為60%：(PSNR=20.829 8)

5 總結(jié)

由上述結(jié)果可知：本文算法相對(duì)于單閾值軟迭代算法更為優(yōu)秀，PSNR提升約為1.01，相對(duì)于欠采樣的地震數(shù)據(jù)PSNR提升為3.6525。且從上述結(jié)果圖中可以較為直接的看出，本文算法在地震數(shù)據(jù)重建中能夠更好的還原地震數(shù)據(jù)的紋理與細(xì)節(jié)。需要指出的是，由于使用的雙樹復(fù)小波變換，相對(duì)傳統(tǒng)離散小波變換其濾波器的構(gòu)造更為復(fù)雜，且雙閾值迭代每次迭代都需要與當(dāng)前子層的父層一起進(jìn)行迭代運(yùn)算，需要更多的運(yùn)算時(shí)間。

圖6 使用雙樹復(fù)小波雙閾值軟迭代進(jìn)行重建后的地震數(shù)據(jù)：(PSNR=24.482 3)

圖7 雙閾值迭代過(guò)程中PSNR的變化曲線