郭 萌,張會(huì)星,劉明珠

(1.中國(guó)海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,山東青島266100;2.青島海洋科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室,山東青島266071)

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

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

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

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

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

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

d=Lm

(1)

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

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

d=Km

(2)

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

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

m=Φ-1x

(3)

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

d=KΦ-1x=Ax

(4)

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

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

(5)

式中:‖·‖0為L(zhǎng)0范數(shù)。

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

(6)

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

(7)

或

(8)

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

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

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

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

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

(9)

(10)

(11)

i=1,…,n

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

(12)

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

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

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

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

(13)

2.2 分裂Bregman迭代去噪算法

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

(14)

(15)

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

(16)

(17)

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

(18)

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

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

(19)

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

(20)

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

(21)

(22)

分裂Bregman迭代的迭代停止準(zhǔn)則為:

(23)

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

圖1 雙重Bregman迭代流程

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

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

圖2 Marmousi部分模型

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

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

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

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

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

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

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

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

圖7 圖6藍(lán)色方框的放大顯示

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

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

5 結(jié)論

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

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