郭 萌,張會星,劉明珠

(1.中國海洋大學海底科學與探測技術(shù)教育部重點實驗室,山東青島266100;2.青島海洋科學與技術(shù)國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室,山東青島266071)

陸上地震勘探受地形、障礙物、禁采區(qū)等采集環(huán)境以及采集成本、采集設備等的影響,采集到的地震數(shù)據(jù)常常出現(xiàn)道間距過大、壞道、部分偏移距數(shù)據(jù)缺失等問題,使得數(shù)據(jù)連續(xù)性變差,用這種不連續(xù)的地震數(shù)據(jù)進行成像處理時往往會產(chǎn)生空間假頻或虛假繞射,從而降低了地震剖面的分辨率[1]。海上拖纜和海底地震勘探受成本、設備與洋流等影響,采集到的數(shù)據(jù)也常常出現(xiàn)空間方向的不規(guī)則、道距稀疏等現(xiàn)象,不規(guī)則的稀疏數(shù)據(jù)同樣會在偏移剖面中產(chǎn)生空間假頻、繞射噪聲或虛假同相軸等問題[2],造成這些問題的原因在于現(xiàn)有的許多處理技術(shù)都假設采集到的數(shù)據(jù)是完整連續(xù)的,當實際輸入數(shù)據(jù)不滿足這一假設條件時必然會產(chǎn)生誤差,因此需要在成像處理之前采用合適的技術(shù)對這種缺道地震數(shù)據(jù)進行重構(gòu),使其滿足后續(xù)成像處理所需的假設條件。

地震數(shù)據(jù)插值重構(gòu)方法的研究始于20世紀80年代,LARNER等[3]基于同相軸傾角一致的假設,給出了一種地震數(shù)據(jù)分時窗傾角掃描并內(nèi)插的方法。之后業(yè)界先后發(fā)展出一系列地震數(shù)據(jù)插值技術(shù),主要分為以下4類,第一類為基于預測濾波的方法,如SPITZ[4]提出的f-x域缺失地震數(shù)據(jù)插值恢復方法,該方法能夠一定程度上克服空間假頻,但只對規(guī)則采樣且當?shù)卣鹩涗浀耐噍S為線性時的數(shù)據(jù)有用。國九英等[5]在f-x域預測濾波的基礎上將其拓展到f-k域,提高了插值計算效率。隨后GULUNAY等[6]提出抗假頻f-k域缺失地震道插值方法,并對該方法進行擴展,提出適用性更廣泛的廣義f-k域缺失地震道插值方法[7]?；陬A測濾波的重構(gòu)方法可以加密空間采樣率,防止空間假頻,但要求輸入數(shù)據(jù)必須規(guī)則采樣。第二類為基于稀疏變換的方法,如傅里葉變換[2]、Radon變換[8]、曲波變換[9-10]等,這類方法對于非規(guī)則缺失數(shù)據(jù)插值效果較好,但該類方法大都需要人為設定閾值,插值效果依賴于處理人員的經(jīng)驗。第三類為矩陣降秩的方法[11],該類方法將高維數(shù)據(jù)通過低秩矩陣進行表示然后重建,計算簡單快速,并發(fā)展到了五維地震數(shù)據(jù)重構(gòu)[12],但該方法的計算量會隨著地震數(shù)據(jù)維度的增加而成倍地增長,對硬件要求很高。第四類為基于波動方程的方法,FOMEL[13]基于有限差分法進行波動方程缺失數(shù)據(jù)重構(gòu),LIU等[14]基于快速傅里葉變換的共軛梯度并結(jié)合快速矢量矩陣乘法的數(shù)值算法進行了空間帶限信號的波場重建。這類方法需要速度信息,并且計算量大,因而沒有得到廣泛應用,而且該類方法主要是對滿足傳統(tǒng)Nyquist采樣定律的地震數(shù)據(jù)進行插值加密處理,對于不滿足傳統(tǒng)Nyquist采樣定理的地震數(shù)據(jù)插值效果不甚理想。

壓縮感知理論的提出,給信號恢復和重構(gòu)帶來新的思路,地球物理學者將其引入地震數(shù)據(jù)恢復領域[15-17]。該理論以信號的稀疏性為前提,可以對遠低于Nyquist標準的欠采樣數(shù)據(jù)進行準確恢復,打破了Nyquist采樣定理的約束。FOMEL等[18]和LIU等[19]基于該理論提出了適用于地震數(shù)據(jù)的Seislet變換,實現(xiàn)了地震數(shù)據(jù)插值重構(gòu)。此外基于壓縮感知的Bregman迭代[20]、迭代閾值法[21]等地震數(shù)據(jù)重構(gòu)方法得到不斷研究和應用。

由于采集環(huán)境等因素的影響,野外地震數(shù)據(jù)往往存在隨機噪聲干擾,這種噪聲使得缺失數(shù)據(jù)得不到充分重建,或重構(gòu)精度較低,因此在重構(gòu)過程中需要對隨機噪聲進行適當壓制[22-23]。ZHANG等[24]利用二維曲波變換對地震數(shù)據(jù)同時進行重構(gòu)和去噪,取得了較好的效果。本文基于壓縮感知理論,提出雙重Bregman迭代算法,即將分裂Bregman迭代作為內(nèi)部迭代用于壓制隨機噪聲,將線性Bregman迭代作為外部迭代用于數(shù)據(jù)重構(gòu),更高效地實現(xiàn)含噪缺失地震數(shù)據(jù)的同步重構(gòu)與去噪。

1 地震數(shù)據(jù)重構(gòu)與去噪的數(shù)學反問題

地震數(shù)據(jù)重建可以歸納為一類線性估計問題:假設需重建的地震數(shù)據(jù)為m,實際觀測到的數(shù)據(jù)為d,則d和m之間的關(guān)系可表示為:

d=Lm

(1)

式中:L為空間采樣算子。由觀測到的已知信號d求解需重建數(shù)據(jù)m的過程就是地震數(shù)據(jù)重構(gòu)過程。

FOMEL[25]指出:對缺失地震數(shù)據(jù)進行重構(gòu),可設置一個掩膜矩陣K代替空間采樣算子L,則公式(1)可表示為:

d=Km

(2)

掩膜矩陣K為對角矩陣,K中待重構(gòu)位置處的元素為0,其它位置處的元素為1。

本文目標是利用壓縮感知理論由已知的觀測數(shù)據(jù)d求取未知的完整數(shù)據(jù)m。由于壓縮感知理論假設信號m具有稀疏特性,因此設計一個稀疏基Φ對m進行稀疏表示:

m=Φ-1x

(3)

式中:x為m的變換域系數(shù)。將公式(3)代入公式(2)可得:

d=KΦ-1x=Ax

(4)

式中:A=KΦ-1為感知矩陣,掩膜矩陣K與稀疏矩陣Φ滿足等距約束性準則[26]。

由此,求解完整數(shù)據(jù)m的過程轉(zhuǎn)換為求解稀疏信號x的過程,根據(jù)壓縮感知理論,將公式(4)轉(zhuǎn)化成L0范數(shù)問題進行求解:

(5)

式中:‖·‖0為L0范數(shù)。

L0范數(shù)的求解較困難,利用L1范數(shù)代替L0范數(shù),將L0范數(shù)最優(yōu)化問題轉(zhuǎn)化為L1范數(shù)的凸優(yōu)化問題:

(6)

構(gòu)造懲罰函數(shù)將公式(6)統(tǒng)一轉(zhuǎn)化成下列形式:

(7)

或

(8)

式中:μ和λ為罰參數(shù),用于平衡正則項和保真項。

對(7)式或(8)式進行求解可獲得稀疏信號x,進而通過公式(3)得到完整地震數(shù)據(jù)m。如果利用單位矩陣E代替掩膜矩陣K,則地震數(shù)據(jù)重構(gòu)問題便轉(zhuǎn)換成了地震數(shù)據(jù)去噪問題。問題的關(guān)鍵是如何準確求解公式(7)或公式(8)。

2 基于雙重Bregman迭代的地震數(shù)據(jù)重構(gòu)與去噪

2.1 Bregman迭代重構(gòu)算法

對于公式(7)或公式(8)的求解,YIN等[27]以Bregman迭代正則化算法為基礎,結(jié)合不動點延續(xù)(fixed-point continuation,FPC)算法形成了一種新的Bregman迭代框架,即:

(9)

(10)

(11)

i=1,…,n

其中,閾值算子shrink的定義為:設實數(shù)空間R中,y∈R,α∈R+,則

(12)

為提高公式(10)的求解效率,使用FPC方法求解。根據(jù)文獻[28]可知,FPC算法具體迭代流程如表1所示。

表1 FPC算法實現(xiàn)步驟

YIN等[27]的研究結(jié)果表明,由公式(9)迭代算法重構(gòu)出的數(shù)據(jù)更精確,圖像對比度更高,且由于使用了FPC算法,計算速度更快,因此本文采用公式(9)進行缺失地震數(shù)據(jù)插值處理。

表1迭代過程中,迭代次數(shù)k是決定能否得到最優(yōu)值的關(guān)鍵參數(shù),k過小無法收斂到最優(yōu)解,k過大則會降低計算效率,因此需要通過設定迭代停止準則來控制迭代次數(shù),本文采用的迭代停止準則為:

(13)

2.2 分裂Bregman迭代去噪算法

分裂Bregman迭代算法由GOLDSTEIN等[29]提出,本文基于Rudin-Osher-Fatemi(ROF)模型[29]推導分裂Bregman迭代去噪框架。ROF模型利用數(shù)據(jù)的全變分將圖像降噪問題轉(zhuǎn)化為求解最優(yōu)化問題,對于高斯噪聲,具體的ROF模型表示形式為:

(14)

(15)

參照公式(8)加入罰函數(shù)λ1,可得:

(16)

(17)

可以看出公式(17)與公式(8)形式相同。參照上述Bregman的迭代步驟,利用公式(9)即可將公式(16)轉(zhuǎn)化為以下多變量的Bregman迭代形式:

(18)

上述公式含有兩個變量x,D,根據(jù)文獻[29],為獲得公式(18)的最優(yōu)解,解耦公式(18)中的L1部分和L2部分,分別對D和x進行迭代最小化,迭代過程如下。

第一步:固定Dk,更新xk。公式為:

(19)

第二步:固定xk,更新Dk。公式為:

(20)

第三步:更新bk。公式為:

(21)

(22)

分裂Bregman迭代的迭代停止準則為:

(23)

本文采用傅里葉變換對地震數(shù)據(jù)進行稀疏表示,采用掩膜算子作為觀測矩陣,以YIN等[27]提出的Bregman迭代作為外部迭代,GOLDSTEIN等[29]提出的分裂Bregman迭代作為內(nèi)部迭代,形成雙重Bregman迭代,在滿足迭代停止準則的條件下完成對含噪聲的缺失地震記錄進行重構(gòu)和去噪。具體實現(xiàn)流程見圖1。

圖1 雙重Bregman迭代流程

3 模型數(shù)據(jù)測試

為驗證雙重Bregman算法的迭代效果,選擇Marmousi部分模型(圖2)正演數(shù)據(jù)的縱波分量進行算法測試。正演記錄由彈性波方程有限差分模擬得到,模擬所用的參數(shù)如下:地表激發(fā)和接收,激發(fā)震源為主頻30Hz的Ricker子波,接收道為601道,道間距為5m,采樣間隔為0.35ms,記錄長度為2100ms。圖3a為模擬單炮記錄,圖3b為加入信噪比為3.46dB的高斯隨機噪聲后的單炮記錄,隨機抽取30%的地震道充零(圖3c),對充零后的地震記錄進行雙重Bregman迭代。根據(jù)迭代準則(13)式和(23)式,設置內(nèi)部迭代允許最大迭代誤差為10-3,外部迭代允許最大誤差為10-5;為避免出現(xiàn)程序無限迭代不能正常退出的情況,經(jīng)多次試驗,同時限定了內(nèi)部和外部迭代的最大迭代次數(shù)為20次。計算得到的重構(gòu)結(jié)果如圖3d所示,圖3e所示為算法誤差剖面,是圖3a和圖3d相減后的結(jié)果。結(jié)合圖3d和圖3e可以看出,對含隨機噪聲的缺失地震數(shù)據(jù),采用本文算法實現(xiàn)了地震數(shù)據(jù)的重構(gòu)與去噪同步進行,可以有效壓制隨機噪聲,并同時獲得完整的重構(gòu)數(shù)據(jù),提高了信噪比。

圖2 Marmousi部分模型

利用本文方法對不同缺失程度下地震數(shù)據(jù)的重構(gòu)與去噪效果以及該方法計算誤差剖面如圖4所示。其中,圖4a和圖4d分別為對信噪比為3.46dB的Marmousi地震數(shù)據(jù)進行50%地震道充零和70%地震道充零,然后進行與30%地震道充零相同的雙重Bregman迭代重構(gòu)與去噪處理,分別得到如圖4b和圖4e所示的重構(gòu)與去噪結(jié)果,以及圖4c和圖4f的誤差剖面。圖5為3種情況下的誤差曲線。另外,為對比本文方法與線性Bregman迭代的重構(gòu)效果與計算效率,對上述3種地震數(shù)據(jù)同時進行了線性Bregman迭代計算。計算環(huán)境參數(shù)為:CPU 2.60GHz,RAM 4.0GB,64位Windows 7系統(tǒng)。重構(gòu)和去噪時間以及處理后的信噪比量化評價結(jié)果見表2。綜合圖3、圖4、圖5以及表2,可得如下結(jié)論:當欠采樣數(shù)據(jù)為完整數(shù)據(jù)的70%時,利用本文方法能夠得到相對完整的重構(gòu)數(shù)據(jù),處理后的信噪比為8.73,得到了大幅度提升;當欠采樣數(shù)據(jù)為完整數(shù)據(jù)的50%時,本文方法基本上實現(xiàn)了對缺失數(shù)據(jù)的重構(gòu),并且重構(gòu)與處理后的信噪比由原來的3.46提升到4.46;當欠采樣數(shù)據(jù)為完整數(shù)據(jù)的30%時,本文方法能夠較完整地重建出缺失數(shù)據(jù),但仍有部分數(shù)據(jù)未能得到準確重構(gòu),而且可能由于缺失程度較大、迭代過程中引入重構(gòu)噪聲等因素,迭代處理后的數(shù)據(jù)信噪比有所降低。經(jīng)測試發(fā)現(xiàn),當欠采樣數(shù)據(jù)為完整數(shù)據(jù)的一半及以上時,能夠利用本文方法得到相對完整的重構(gòu)與去噪成果,并且數(shù)據(jù)缺失越少,重構(gòu)和去噪所需計算時間越少,誤差曲線下降越快,信噪比提升也越高;在不同缺失程度下,雙重Bregman迭代與線性Bregman迭代的運行時間相差不大。這是因為分裂Bregman迭代的速度相當快。GOLDSTEIN等[29]曾指出,分裂Bregman迭代速度很快的主要原因是由其計算方法決定的,分裂Bregman迭代是將L1和L2部分進行解耦,然后分別用快速收斂的Gauss-Seidel方法和閾值算子方法求解,這種算法加快了分裂Bregman迭代的速度;另外,分裂Bregman的迭代速度與參數(shù)選擇有很大關(guān)系,參數(shù)選擇得當,該算法非常快速。綜上,雙重Bregman迭代不僅有效去除了隨機噪聲,其結(jié)果的信噪比明顯高于線性Bregman迭代結(jié)果,而且運算時間并沒有大量增加,因此,對于含有隨機噪聲的地震數(shù)據(jù),本文方法優(yōu)于線性Bregman迭代方法,具有一定的研究意義及應用價值。

圖3 Marmousi模型含噪聲的缺失地震記錄重構(gòu)與去噪

圖4 不同缺失程度下雙重Bregman迭代對含噪聲地震數(shù)據(jù)的重構(gòu)與去噪結(jié)果和誤差剖面

表2 雙重Bregman迭代和線性Bregman迭代在不同缺失程度下的重構(gòu)和去噪時間及信噪比

圖5 不同缺失程度下含噪聲地震數(shù)據(jù)的重構(gòu)與去噪誤差曲線

4 實際數(shù)據(jù)測試

利用海上某工區(qū)實際地震數(shù)據(jù)進行測試,原始地震數(shù)據(jù)共301道,每道10001個采樣點,重構(gòu)工作在CDP道集(圖6a)上進行,對該數(shù)據(jù)隨機抽取30%的地震道充零后的結(jié)果如圖6b,對抽稀后的地震數(shù)據(jù)進行雙重Bregman迭代,迭代停止準則采用公式(13)和公式(23),經(jīng)多次試驗,同時設置內(nèi)部迭代允許最大迭代次數(shù)為50次,外部迭代允許最大迭代次數(shù)為100次,得到的結(jié)果如圖6c。為說明本文算法的重構(gòu)與去噪效果,對圖6b進行線性Bregman迭代計算,迭代停止準則為公式(13),允許最大迭代次數(shù)為100次,迭代結(jié)果如圖6d。結(jié)合圖6中藍色方框?qū)木植糠糯蟛ㄐ物@示(圖7)可以看出,兩種迭代算法都能夠?qū)θ笔?shù)據(jù)進行較好恢復,但本文迭代算法與原始數(shù)據(jù)的重合度更高。將兩種算法迭代后的地震數(shù)據(jù)與原始地震數(shù)據(jù)求差,其中1.5～2.5s時間段第100～150道(即藍色方框?qū)秶?的誤差如圖8所示,整體地震數(shù)據(jù)的重構(gòu)誤差如圖9所示,對比可以看出,雙重Bregman算法迭代誤差遠遠小于線性Bregman迭代,進一步顯現(xiàn)出本文迭代算法的優(yōu)越性。

圖6 實測地震記錄的重構(gòu)與去噪

圖7 圖6藍色方框的放大顯示

圖8 兩種Bregman迭代局部結(jié)果的精度比較

圖9 兩種Bregman迭代整體結(jié)果的精度比較

5 結(jié)論

本文結(jié)合Bregman迭代重構(gòu)算法與分裂Bregman迭代去噪算法,提出了一種基于雙重Bregman迭代的地震數(shù)據(jù)重構(gòu)與去噪方法。Marmousi模型數(shù)值模擬實驗和實際數(shù)據(jù)應用結(jié)果表明:本文算法沒有因兩種獨立算法的結(jié)合而減弱對數(shù)據(jù)的處理效果,且只需較少次數(shù)的迭代就能高效實現(xiàn)含噪聲缺道地震數(shù)據(jù)的同時重構(gòu)與去噪,為地震數(shù)據(jù)恢復提供了一種新的缺失地震數(shù)據(jù)處理方法。

基于壓縮感知的地震數(shù)據(jù)處理方法需要先將地震數(shù)據(jù)通過數(shù)學變換進行稀疏表示,本文采用了最常用的傅里葉稀疏表示方法,而該方法只適用于同相軸近似線性或者平穩(wěn)變化的地震數(shù)據(jù),因此本文進一步的工作是尋找復雜地質(zhì)構(gòu)造下地震數(shù)據(jù)的最佳稀疏表示方法,結(jié)合雙重Bregman迭代,進一步研究復雜地質(zhì)構(gòu)造下含噪聲缺失地震數(shù)據(jù)的高精度重構(gòu)與去噪方法。