基于偽地震數(shù)據(jù)模式學(xué)習(xí)的多次波自適應(yīng)相減方法

姜博午,劉金朋,陸文凱

(1.北京化工大學(xué)機(jī)電工程學(xué)院,北京100029;2.智能技術(shù)與系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100084;3.北京信息科學(xué)與技術(shù)國(guó)家研究中心,北京 100084;4.清華大學(xué)自動(dòng)化系,北京 100084;5.中海油田服務(wù)股份有限公司物探事業(yè)部特普公司,廣東湛江524057)

工業(yè)界常用的基于伴隨算子的地震偏移成像算法假設(shè)地震數(shù)據(jù)中只存在一次反射波(一次波),而多次反射波(多次波)的存在不符合傳統(tǒng)成像算法的假設(shè)條件,造成有效波振幅、頻率、相位畸變,影響地震解釋[1]。多次波廣泛存在于陸地和海洋地震數(shù)據(jù)中,特別是在海洋地震數(shù)據(jù)中,由于存在海面和海底兩個(gè)強(qiáng)反射界面,因而導(dǎo)致海洋地震數(shù)據(jù)的表面多次波干擾問題更為嚴(yán)重。對(duì)于陸上地震數(shù)據(jù)而言,當(dāng)?shù)叵麓嬖趶?qiáng)反射界面,諸如基巖面、火成巖侵入體、不整合面等,地震波在強(qiáng)反射界面間或其內(nèi)部發(fā)生多次反射,形成層間多次波[2]。由于多次波傳播路徑更加復(fù)雜,因而壓制層間多次波更有挑戰(zhàn)性。多次波壓制方法主要分為兩類:濾波法和預(yù)測(cè)相減法[3]。濾波類方法利用一次波和多次波在變換域內(nèi)的差異壓制多次波,包括F-K濾波[4]、Radon變換濾波方法[5-8]等。通常當(dāng)一次波和多次波速度差異較大時(shí),濾波類方法效果較理想;但當(dāng)?shù)叵陆橘|(zhì)復(fù)雜時(shí),濾波類方法很難在保護(hù)有效信號(hào)的同時(shí)壓制多次波。預(yù)測(cè)相減法首先基于波動(dòng)理論得到預(yù)測(cè)多次波模型,然后通過自適應(yīng)相減壓制多次波。這類算法更加適用于復(fù)雜介質(zhì)情況。BERKHOUT等[9]與VERSCHUUR等[10]基于反饋迭代模型,提出了表面多次波去除方法(SRME),該方法已經(jīng)成為工業(yè)上成熟應(yīng)用的表面多次波去除方法。針對(duì)層間多次波的復(fù)雜波場(chǎng)特征,諸多國(guó)內(nèi)外地球物理研究人員發(fā)展了相應(yīng)的層間多次波壓制方法,包括逆散射級(jí)數(shù)法[7,11-13]、共聚焦點(diǎn)技術(shù)[14-15]和Marchenko層間多次波壓制技術(shù)[16-17]等。

但是,在實(shí)際地震數(shù)據(jù)采集過程中,觀測(cè)系統(tǒng)往往存在缺陷,震源子波常常不一致等。這些因素使得預(yù)測(cè)多次波模型往往存在相位、頻帶、振幅的差異[9]。因此,自適應(yīng)相減方法作為預(yù)測(cè)相減法不可或缺的關(guān)鍵步驟,能夠校正預(yù)測(cè)多次波模型中的誤差。工業(yè)界最常用的基于最小二乘匹配濾波器的多次波自適應(yīng)相減方法[10]利用預(yù)測(cè)多次波直接匹配觀測(cè)數(shù)據(jù),以去除多次波后的數(shù)據(jù)的L2范數(shù)最小為優(yōu)化目標(biāo)。該方法效率很高,可以快速估計(jì)出觀測(cè)數(shù)據(jù)中的多次波。但是其僅在一次波和多次波正交假設(shè)下,才能取得較好的效果。隨后,基于L1范數(shù)的匹配濾波器方法被提出[18]。該方法可以更好地處理一次波和多次波不正交的數(shù)據(jù)。獨(dú)立成分分析方法也被應(yīng)用于多次波自適應(yīng)相減問題[19-20]。

基于模式學(xué)習(xí)的多次波自適應(yīng)相減方法相繼被提出[21-24]。模式學(xué)習(xí)方法首先從預(yù)測(cè)多次波中學(xué)習(xí)多次波特征,然后利用學(xué)習(xí)到的多次波特征去重構(gòu)原始數(shù)據(jù),這樣就可以得到原始數(shù)據(jù)中的多次波,進(jìn)而將其從原始數(shù)據(jù)中減去即可得到一次波數(shù)據(jù)。模式學(xué)習(xí)方法將多次波衰減分為多次波模式學(xué)習(xí)和自適應(yīng)相減兩個(gè)獨(dú)立過程,能夠更好地保護(hù)有效信號(hào)。本文基于混合L2-L1范數(shù)模式學(xué)習(xí)方法[24],衰減地震數(shù)據(jù)中的多次波。同時(shí),針對(duì)預(yù)測(cè)多次波模型中存在的相位和頻帶誤差,引入偽地震數(shù)據(jù)[25-28],在不損傷有效信號(hào)的同時(shí),進(jìn)一步減少多次波殘留。最后利用模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)驗(yàn)證本文方法的有效性。

1 技術(shù)方法

1.1 混合L2-L1范數(shù)模式學(xué)習(xí)方法

采用SRME預(yù)測(cè)多次波模型m∈h×w,其中,w表示地震道數(shù),h表示地震道采樣點(diǎn)數(shù)。我們利用主成分分析(Principle Component Analysis,PCA)方法來(lái)提取預(yù)測(cè)多次波的關(guān)鍵模式。首先將預(yù)測(cè)多次波模型m劃分為n個(gè)互相重疊的大小為d=k×l的圖像塊,再將這些圖像塊向量化之后合成一個(gè)矩陣M∈d×n,其中,n=(h-k+1)×(w-l+1)。對(duì)該矩陣進(jìn)行PCA特征編碼,有:

(1)

式中:B∈d×p表示學(xué)習(xí)到的p個(gè)d維特征向量構(gòu)成的特征矩陣;X∈p×n表示對(duì)應(yīng)的特征系數(shù)。特征矩陣B的構(gòu)建過程如下。

算法1:SVD求解公式。

1)初始化:通常λ=0.1。

2)對(duì)數(shù)據(jù)矩陣M作SVD分解,有[U,Σ,V]=SVD(M)。

3)設(shè)Σ=diag({σi}1≤i≤N),選擇前p個(gè)最大的特征值,使得σp+1<λσ1,且σp≥λσ1。

4)設(shè)U=[u1,u2,…,un],選擇前p個(gè)最大的特征值對(duì)應(yīng)的特征向量,構(gòu)成特征矩陣B=[u1,u2,…,up]。

考慮分頻編碼策略,用帶通濾波器將預(yù)測(cè)多次波模型分為低、中、高3個(gè)不同的頻段,即mLF,mMF,mHF。相應(yīng)的低、中、高頻數(shù)據(jù)矩陣記為MLF,MMF,MHF。利用算法1學(xué)習(xí)分頻段特征矩陣BLF,BMF,BHF。

最終的分頻特征矩陣為:

B=[BLF,BMF,BHF]

(2)

由此,我們可以得到預(yù)測(cè)多次波m的特征矩陣B[24]。在完成模式學(xué)習(xí)并學(xué)習(xí)到多次波的特征之后,我們便可以利用學(xué)習(xí)到的特征去重構(gòu)觀測(cè)數(shù)據(jù)中包含的多次波,進(jìn)而分離出觀測(cè)數(shù)據(jù)的一次波,即自適應(yīng)相減過程。本質(zhì)上是估計(jì)一組系數(shù),使得多次波關(guān)鍵模式線性組合后,與觀測(cè)數(shù)據(jù)的差異最小。

與模式學(xué)習(xí)過程相對(duì)應(yīng),我們將觀測(cè)數(shù)據(jù)s(包含多次波和一次波)劃分為若干個(gè)大小為k×l的重疊圖像塊,再將這些圖像塊向量化之后組合成一個(gè)矩陣S∈d×n。在數(shù)據(jù)矩陣S中,每個(gè)列向量對(duì)應(yīng)一個(gè)原數(shù)據(jù)s中的圖像塊。我們用L1范數(shù)來(lái)重構(gòu)多次波[24]:

(3)

令P=S-BY即一次波數(shù)據(jù)矩陣,其中P∈d×n,公式重寫為:

(4)

由拉格朗日乘數(shù)法,引入一個(gè)罰函數(shù),我們將公式(4)轉(zhuǎn)化為一個(gè)無(wú)約束的優(yōu)化問題:

(5)

其中,μ為懲罰系數(shù),且μ>0。當(dāng)μ較小時(shí),公式(5)求解不準(zhǔn)確;當(dāng)μ較大時(shí),公式(5)求解數(shù)值不穩(wěn)定。為避免μ調(diào)參復(fù)雜,我們利用裂步Bregman算法[29]求解公式(5)。引入中間變量Ak,且A1=0,則:

(6)

Ak+1=Ak+S-BYk+1-Pk+1

(7)

在裂步Bregman算法[29]中,公式(7)用來(lái)添加約束,避免選擇一個(gè)較大的μ帶來(lái)的數(shù)值不穩(wěn)定問題。通常,我們選擇μ=1。公式(6)兩個(gè)變量P和Y是解耦的,我們可以利用交替下降法近似求解公式(6):

(8)

(9)

公式(8)和公式(9)的解析解為:

Yk+1=(BTB)-1BT(S+Ak-Pk)

(10)

(11)

裂步Bregman算法[29]求解自適應(yīng)相減問題(公式(3))的步驟總結(jié)如下。

算法2:裂步Bregman算法[29]求解公式。

1)初始化:k=1,A1=0,P1=0,P0=S,D=(BTB)-1。

2)While‖Pk-Pk-1‖2>γdo。

3)Yk+1=DBT(S+Ak-Pk)。

5)Ak+1=Ak+(S-BYk+1-Pk+1)。

6)k=k+1。

7)End while。

1.2 基于偽地震數(shù)據(jù)的模式學(xué)習(xí)

(12)

在得到特征矩陣B后,利用算法2重構(gòu)觀測(cè)數(shù)據(jù)中的多次波,進(jìn)而分離一次波和多次波。圖1展示了基于偽地震數(shù)據(jù)的多次波模式學(xué)習(xí)和自適應(yīng)相減方法流程。

圖1 基于偽地震數(shù)據(jù)的多次波模式學(xué)習(xí)和自適應(yīng)相減方法流程

2 應(yīng)用實(shí)例

2.1 Sigsbee模擬數(shù)據(jù)

采用Sigsbee模擬數(shù)據(jù)驗(yàn)證算法的有效性[23]。圖2a,圖2b和圖2c分別展示了真實(shí)一次波、預(yù)測(cè)多次波模型和觀測(cè)數(shù)據(jù)。我們對(duì)真實(shí)多次波進(jìn)行如下處理生成預(yù)測(cè)多次波模型:縱向錯(cuò)動(dòng)3個(gè)網(wǎng)格,橫向錯(cuò)動(dòng)2個(gè)網(wǎng)格,然后與主頻20Hz,長(zhǎng)30個(gè)單位網(wǎng)格的雷克子波進(jìn)行褶積。圖3a,圖3b和圖3c分別展示了采用分頻模式學(xué)習(xí)方法估計(jì)的一次波、多次波和誤差。圖4a,圖4b和圖4c分別展示了采用偽地震數(shù)據(jù)模式學(xué)習(xí)方法估計(jì)的一次波、多次波和誤差。對(duì)比圖3c 和圖4c可以發(fā)現(xiàn),偽地震數(shù)據(jù)模式學(xué)習(xí)方法誤差更小,在未損傷一次波的基礎(chǔ)上,多次波殘留更少。定量來(lái)看,分頻模式學(xué)習(xí)方法和偽地震數(shù)據(jù)模式學(xué)習(xí)方法估計(jì)的一次波信噪比分別為19.03dB和20.06dB,偽地震數(shù)據(jù)模式學(xué)習(xí)方法的信噪比提升了1dB。圖2中箭頭所指位置強(qiáng)能量多次波造成分頻模式學(xué)習(xí)方法存在多次波殘留,而偽地震數(shù)據(jù)模式學(xué)習(xí)方法明顯殘留更少。為了更清晰地對(duì)比,我們將圖2 至圖4中紅色框部分放大顯示(圖5至圖7)。由于繞射多次波預(yù)測(cè)不準(zhǔn),分頻模式學(xué)習(xí)方法殘留繞射多次波嚴(yán)重。偽地震數(shù)據(jù)模式學(xué)習(xí)方法幾乎完全壓制了繞射多次波。因此,相比于分頻模式學(xué)習(xí)方法,偽地震數(shù)據(jù)模式學(xué)習(xí)方法在不損傷一次波的基礎(chǔ)上,能夠有效避免多次波殘留,特別是對(duì)繞射多次波壓制效果更佳。

圖2 Sigsbee模擬數(shù)據(jù)a 真實(shí)一次波;b 預(yù)測(cè)多次波;c 觀測(cè)數(shù)據(jù)

圖3 采用分頻模式學(xué)習(xí)方法估計(jì)的一次波(a)、多次波(b)和誤差(c)

圖4 采用偽地震數(shù)據(jù)模式學(xué)習(xí)方法估計(jì)的一次波(a)、多次波(b)和誤差(c)

圖5 圖2紅色框位置的放大顯示結(jié)果a 真實(shí)一次波;b 預(yù)測(cè)多次波;c 觀測(cè)數(shù)據(jù)

圖6 圖3紅色框位置的放大顯示結(jié)果a 分頻模式學(xué)習(xí)估計(jì)的一次波;b 分頻模式學(xué)習(xí)估計(jì)的多次波;c 分頻模式學(xué)習(xí)的誤差

圖7 圖4紅色框位置的放大顯示結(jié)果a 偽地震數(shù)據(jù)模式學(xué)習(xí)估計(jì)的一次波;b 偽地震數(shù)據(jù)模式學(xué)習(xí)估計(jì)的多次波;c 偽地震數(shù)據(jù)模式學(xué)習(xí)的誤差

2.2 海上實(shí)際數(shù)據(jù)

我們對(duì)海上實(shí)際地震資料進(jìn)行了多次波壓制。圖8a和圖8b分別顯示了觀測(cè)數(shù)據(jù)和預(yù)測(cè)多次波模型;圖9a和圖9b分別顯示了采用分頻模式學(xué)習(xí)方法估計(jì)的一次波和多次波;圖10a和圖10b分別顯示了采用偽地震數(shù)據(jù)模式學(xué)習(xí)方法估計(jì)的一次波和多次波。對(duì)比圖8 至圖10中箭頭所指處可以看出,強(qiáng)能量多次波導(dǎo)致分頻模式學(xué)習(xí)方法存在一個(gè)低頻殘留;而偽地震數(shù)據(jù)模式學(xué)習(xí)方法不存在這一問題。放大顯示圖8至圖10中的矩形框1和矩形框2,如圖11至圖16所示。圖11至圖13中箭頭所指位置尖滅構(gòu)造生成復(fù)雜多次波。由于多次波預(yù)測(cè)過程中多次褶積的影響,存在相位誤差。分頻模式學(xué)習(xí)方法存在一個(gè)明顯和多次波走勢(shì)相同的殘留。偽地震數(shù)據(jù)模式學(xué)習(xí)方法減少了多次波殘留。圖14至圖16中箭頭所指位置繞射多次波發(fā)育。但由于大角度繞射多次波預(yù)測(cè)誤差較大,導(dǎo)致分頻模式學(xué)習(xí)方法并沒有很好地壓制掉大角度的繞射多次波,而偽地震數(shù)據(jù)模式學(xué)習(xí)方法能夠緩解這一問題。但是對(duì)于更大角度處,預(yù)測(cè)多次波時(shí)移誤差更大,仍然存在部分殘留。

圖8 海上實(shí)際地震資料a 觀測(cè)數(shù)據(jù);b 預(yù)測(cè)多次波

圖9 采用分頻模式學(xué)習(xí)方法估計(jì)的一次波(a)和多次波(b)

圖10 采用偽地震數(shù)據(jù)模式學(xué)習(xí)方法估計(jì)的一次波(a)和多次波(b)

圖11 圖8中矩形框1的放大顯示a 觀測(cè)數(shù)據(jù);b 預(yù)測(cè)多次波

圖12 圖9中矩形框1的放大顯示a 分頻模式學(xué)習(xí)估計(jì)的一次波;b 分頻模式學(xué)習(xí)估計(jì)的多次波

圖13 圖10中矩形框1的放大顯示a 偽地震數(shù)據(jù)模式學(xué)習(xí)估計(jì)的一次波;b 偽地震數(shù)據(jù)模式學(xué)習(xí)估計(jì)的多次波

圖14 圖8中矩形框2的放大顯示結(jié)果a 觀測(cè)數(shù)據(jù);b 預(yù)測(cè)多次波

圖15 圖9中矩形框2的放大顯示結(jié)果a 分頻模式學(xué)習(xí)估計(jì)的一次波;b 分頻模式學(xué)習(xí)估計(jì)的多次波

圖16 圖10中矩形框2的放大顯示結(jié)果a 偽地震數(shù)據(jù)模式學(xué)習(xí)估計(jì)的一次波;b 偽地震數(shù)據(jù)模式學(xué)習(xí)估計(jì)的多次波

3 結(jié)論

本文提出了一種基于偽地震數(shù)據(jù)的多次波模式學(xué)習(xí)和自適應(yīng)相減方法。該方法一方面能夠繼承模式學(xué)習(xí)方法的優(yōu)勢(shì),更好地保護(hù)有效信號(hào);另一方面,利用了偽地震數(shù)據(jù)中額外的相位和高頻信息,避免了多次波殘留。Sigsbee模擬數(shù)據(jù)實(shí)驗(yàn)結(jié)果表明,相比于分頻模式學(xué)習(xí)方法,本文提出的偽地震數(shù)據(jù)模式學(xué)習(xí)方法一次波信噪比能夠提升1dB。海上實(shí)際數(shù)據(jù)應(yīng)用結(jié)果表明,分頻模式學(xué)習(xí)方法在多次波復(fù)雜情況下效果不佳,本文方法能夠更好地消除大角度的繞射多次波。但是當(dāng)繞射多次波預(yù)測(cè)時(shí)移誤差更大時(shí),仍然存在一定殘留。事實(shí)上,可以通過進(jìn)一步調(diào)參來(lái)消除更大時(shí)移誤差,但會(huì)損傷部分有效信號(hào)。因此,對(duì)于實(shí)際地震數(shù)據(jù),多次波預(yù)測(cè)誤差變化較大時(shí),必須要在保護(hù)一次波和壓制多次波之間進(jìn)行折衷處理。