黃建平, 張入化, 國運東, 雍 鵬, 李 闖

(1.中國石油大學(華東)地球科學與技術(shù)學院,山東青島 266580;2.西安交通大學電子信息工程學院,陜西西安 710000)

近年來,油氣勘探的地質(zhì)目標越來越精細化,對勘探、開發(fā)的要求進一步提高,發(fā)展更加精確的偏移成像技術(shù)成為必然趨勢[1-2]?；诜囱菟枷氲淖钚《似?LSM)和最小二乘逆時偏移(LSRTM)雖然具有高精度保幅成像的特點,但計算量較大。多炮編碼策略如相位編碼、振幅編碼、極性編碼、平面波編碼等[3-9]能一定程度提高計算效率,但多炮數(shù)據(jù)混疊也會引入串擾噪聲。李闖等[10]將隨機最優(yōu)化相位編碼引入多炮編碼中進一步減少了串擾噪聲。PLSRTM也可以通過疊加不同射線參數(shù)的平面波成像結(jié)果來壓制串擾噪聲[11-12],但平面波道集數(shù)量增加會增大計算量。Xue等[13]在目標函數(shù)后加上了整形正則化約束項,能一定程度地改善多震源數(shù)據(jù)的成像質(zhì)量。Li等[14]則將奇異譜值分析(SSA)引入到PLSRTM中作為正則化約束。除此之外,根據(jù)地震數(shù)據(jù)在曲波域稀疏分布的特點,Herrmann等[15]將曲波變換用作偏移的預(yù)處理。相較于來自圖像算法的曲波變換,Seislet變換更加適用于地震數(shù)據(jù),Xue[16]首先將Seislet約束用于全波形反演中,Dutta[17]將Seislet變換引入最小二乘逆時偏移中,有效地提高了成像質(zhì)量。在前人研究的基礎(chǔ)上,筆者利用Riemann-Liouville分數(shù)階積分理論推導得到分數(shù)階閾值函數(shù),并與Seislet變換相結(jié)合得到Seislet分數(shù)階閾值算法,再將其作為稀疏約束算子,實現(xiàn)Seislet分數(shù)階閾值算法約束的平面波最小二乘逆時偏移。

1 方法原理

1.1 平面波編碼

多炮編碼的方式主要有相位編碼、振幅編碼、平面波編碼等方式。在二維地震數(shù)據(jù)中對共接收點道集(CRP)進行平面波編碼。

(1)

式中,d(xr,xs,ω)為共接收點道集,由位于橫坐標xs處的震源產(chǎn)生,再被橫坐標位于xr處的檢波點所接收的記錄;eiωpxs為代表時間域的時移量,ω為角頻率,而最主要的時移參數(shù)則是pxs,p為射線參數(shù),p=sinθ/v,θ為地層傾角,v為地層速度。等式(1)主要是將道集中參與時移的道疊加在一起,形成一個平面波道集。

平面波編碼的關(guān)鍵是射線參數(shù)的選擇,最大射線參數(shù)pmax應(yīng)該滿足不等式

(2)

式中,fmax為最大頻率;Δxs為炮間隔。射線參數(shù)的間隔Δp應(yīng)滿足

(3)

式中,θ1為最小傾角;θ2為最大傾角;Np為平面波道集的數(shù)量。射線參數(shù)p不同,角度范圍便不同,照明的區(qū)域便存在差異,增加平面波道集的數(shù)量能夠加強對地下復(fù)雜區(qū)域的有效照明。但過多的平面波道集也會增大計算量,降低計算效率。

1.2 平面波編碼最小二乘逆時偏移

正演的方式主要分為有限差分正演和Born線性正演。其中Born線性正演過程可以表示為

d=Lm.

(4)

式中,L為Born線性正演算子;d為觀測數(shù)據(jù);m為反射系數(shù)。

偏移過程則可以表達為

mmig=LHd.

(5)

式中,mmig為偏移結(jié)果;H為共軛轉(zhuǎn)置;LH為偏移算子,在本文中為逆時偏移算子[18]。

假設(shè)有n個平面波道集,則反射系數(shù)模型M可表示為

M=[m1,m2,m3,…,mn]T.

(6)

式中,mn為第n個平面波道集所對應(yīng)的偏移成像結(jié)果。平面波道集D可以表示為

(7)

式中,pi為平面波的Born正演算子。

平面波最小二乘逆時偏移的誤差函數(shù)可以定義為

(8)

其中誤差函數(shù)的梯度解為

(9)

其中

1.3 約束的平面波最小二乘逆時偏移

模型M的更新等式為

M(i+1)=M(i)+α(i)z(i),

(10)

z(i)=-f(M)(i)+β(i)z(i-1).

(11)

式中,i為第i次迭代;α(i)為搜索步長;z(i)為共軛梯度;β(i)為共軛梯度修正因子。一個n維函數(shù)的共軛方向最多只有n個,大于n次迭代時共軛梯度的計算將失效[19]。為避免迭代誤差的積累,本文中在連續(xù)計算4次共軛梯度后將使用一次最速下降算法再重新使用共軛梯度法。

對模型的更新等式(10)加上約束得

M(i+1)=S-1TS(M(i)+α(i)z(i)).

(12)

式中,S為Seislet正變換;S-1為Seislet反變換;T為閾值函數(shù),本文中T為分數(shù)階閾值函數(shù)。約束算子S-1TS的目的主要在于消除偏移假象和因編碼時炮數(shù)據(jù)混疊產(chǎn)生的串擾噪聲。

如果將某個稀疏變換作為一個整體的算子,那么等式(12)就等同于用線性Bregman迭代法使得下列的正則化目標函數(shù)最小化[20]。

(13)

(14)

Seislet分數(shù)階閾值算法約束的PLSRTM流程如圖1,其中對初始模型M(0)賦零值。

圖1 Seislet分數(shù)階閾值算法約束的PLSRTM流程Fig.1 Workflow of PLSRTM with Seislet fractional threshold algorithm constraint

1.4 分數(shù)階閾值函數(shù)

傳統(tǒng)的閾值函數(shù)主要有硬閾值函數(shù)和軟閾值函數(shù)[21-22]。硬閾值函數(shù)在閾值附近存在斷點,重構(gòu)信號會產(chǎn)生偽吉布斯現(xiàn)象,一般用得比較少。常用的主要是軟閾值函數(shù),但軟閾值函數(shù)處理后的系數(shù)與原系數(shù)相比有恒定的偏差,偏差在數(shù)值上等于閾值,這會使得重構(gòu)信號與原信號之間產(chǎn)生較大的誤差。有學者也提出過半軟閾值,就是在硬閾值函數(shù)和軟閾值函數(shù)中加了一個權(quán)系數(shù)[23],但這個權(quán)系數(shù)需根據(jù)不同的數(shù)據(jù)測試后再確定,且與真實系數(shù)的相比仍存在一定偏差。

(15)

式中,s為Seislet稀疏變換后的系數(shù);λ為閾值。

設(shè)閾值為2,軟閾值函數(shù)如圖2所示,經(jīng)過軟閾值函數(shù)處理后的系數(shù)與原系數(shù)之間存在偏差。

圖2 軟閾值函數(shù)圖Fig.2 Soft threshold function

分數(shù)階閾值函數(shù)應(yīng)用了分數(shù)階微積分理論(Riemann-Liouville分數(shù)階積分公式)[24],用該分數(shù)階積分公式對函數(shù)進行積分處理積分后,能夠增強函數(shù)的連續(xù)性且保留原函數(shù)特征。

本文中提出的分數(shù)階閾值函數(shù)如下:

(16)

圖3 分數(shù)階閾值函數(shù)圖Fig.3 Fractional threshold function

由圖3可見,分數(shù)階閾值函數(shù)在系數(shù)偏差上得到了很好的控制,端點處也得到了一定的平滑。

2 Seislet變換測試

常見的稀疏變換有Fourier變換、小波變換、曲波變換、Seislet變換等。其中Seislet變換是專門針對地震數(shù)據(jù)的一種信號變換,能將數(shù)據(jù)變換到Seislet域[25],但在進行Seislet變換前需要采用平面波解構(gòu)濾波器(PWD)對地震數(shù)據(jù)進行局部傾角估計。Seislet變換需要局部傾角數(shù)據(jù)作為輸入,再沿著傾角方向?qū)Φ卣饠?shù)據(jù)作Seislet變換。Fourier變換、小波變換、曲波變換等均以采樣點為變換的單位,而Seislet變換則是以地震道為單位進行變換。

選擇某一復(fù)雜模型作Seislet變換測試和小波變換測試,該模型數(shù)據(jù)的橫向長度為5 120 m,縱向長度為1 500 m,道數(shù)為512,縱向采樣點數(shù)為150,間距均為10 m。首先對復(fù)雜模型分別進行Seislet變換和小波變換,得到對應(yīng)的系數(shù),再通過選取同樣百分比的Seislet系數(shù)和小波系數(shù)進行反變換得到重構(gòu)模型,再進行對比分析。

圖4(a)所示為某一復(fù)雜模型的真實反射系數(shù),再對其進行Seislet變換,得到Seislet系數(shù)(圖4(b))。由圖4(b)可知,Seislet系數(shù)主要集中在尺度小于50的位置,且整體分布較為集中。隨后分別將Seislet變換和小波變換的系數(shù)經(jīng)過從大到小排序后,選取前1%的系數(shù)得到重構(gòu)后的結(jié)果(圖4(c)、(d)),前者基本重構(gòu)了真實反射系數(shù)的輪廓,但卻沒有重構(gòu)出繞射點(圖4(c)),而后者的右側(cè)反射系數(shù)很雜亂沒有得到有效重構(gòu),但能重構(gòu)下部繞射點(圖4(d))。圖4(e)、(f)則是分別將Seislet變換和小波變換的系數(shù)經(jīng)過從大到小排序后,選取前7%的系數(shù)進行重構(gòu)后的結(jié)果,前者基本與真實反射系數(shù)一致,而后者中部斷層處的反射系數(shù)仍比較模糊。

圖4 Seislet變換測試Fig.4 Seislet transform test

將同樣閾值截取的數(shù)據(jù)投影到原有的尺度矩陣中(橫坐標為尺度)。圖5(a)、(b)是分別將Seislet變換和小波變換的系數(shù)經(jīng)過從大到小排序后,選取前1%的系數(shù)投影到尺度矩陣的結(jié)果,其中黑色數(shù)據(jù)為前1%的系數(shù)。從圖5 (a)中可以看出系數(shù)主要集中在尺度小于15的位置且分布較為集中, 圖5 (b)中的系數(shù)主要在尺度小于50且分布相對圖5 (a)較為分散。同理,再選取前7%的系數(shù)投影到原來的尺度矩陣(圖5(c)、(d)),黑色數(shù)據(jù)為前7%的系數(shù)。從圖5(c)中得知系數(shù)主要集中在尺度小于100的位置,在尺度小于256的范圍都有一定分布;由圖5 (d),可以發(fā)現(xiàn)系數(shù)在全局范圍內(nèi)都有分布,且分布較為分散。

最終,可知Seislet變換比小波變換能夠更好地壓縮地震數(shù)據(jù),更能表示地震數(shù)據(jù)的稀疏性。

圖5 系數(shù)在原來尺度矩陣的分布Fig.5 Distribution of coefficients in original scale matrices

3 數(shù)值測試

3.1 簡單模型測試

選擇一個4層的層狀模型進行測試。該模型的尺寸為1 500 m×5 120 m,橫向采樣點數(shù)為512,縱向采樣點數(shù)為150,橫縱向采樣點距均為10 m。共激發(fā)120炮,炮間隔為40 m,檢波點為512個,檢波點距10 m,最大采樣時間為8 000 ms,震源子波采用雷克子波,震源主頻率為20 Hz。先使用時間2階,空間8階的有限差分算法及PML邊界條件進行正演模擬[26-27]得到單炮記錄,再用平面波靜態(tài)編碼技術(shù)將多炮數(shù)據(jù)合成具有不同射線參數(shù)的平面波道集。

圖6(a)為4層水平模型,第一層速度為1 600 m/s,第二層速度為1 900 m/s,第三層速度為2 300 m/s,第四層速度為2 600 m/s,模型長度為5 120 m,深度為1 500 m,圖6(b)為偏移速度場。根據(jù)公式(1)合成了不同平面波射線參數(shù)p的平面波道集,p的取值范圍為-0.3～0.3 ms/m,且呈線性變化,共計11個平面波道集。

圖7(a)、(c)、(e)分別為簡單層狀模型經(jīng)過平面波最小二乘逆時偏移迭代5、10、25次的結(jié)果?？梢钥闯鰣D7(a)、(c)、(e)都有不同程度的串擾噪音和少量采集腳印,但圖7 (c)的噪音明顯小于圖7 (a),圖7 (e)的噪音也小于圖7 (c),這是因為最小二乘逆時偏移本身隨著迭代次數(shù)的增加能夠一定程度地壓制噪音。

圖7(b)、(d)、(f)分別為PLSRTM迭代5、10、25次結(jié)果對應(yīng)的局部傾角估計。根據(jù)顏色條,顏色越深傾角越大,且經(jīng)過測試發(fā)現(xiàn)層狀模型本身的傾角值為0,表示為最淺的顏色。由圖7(b)、(d)、(f)可知,隨著迭代次數(shù)增加,局部傾角估計圖中顏色團的顏色逐漸變淺,而對應(yīng)偏移結(jié)果的串擾噪音逐漸減弱(圖7(a)、(c)、(e))。因此層狀模型偏移結(jié)果中的非零局部傾角值是由噪音引起的。

圖6 層狀模型速度場和偏移速度場Fig.6 Velocity field of layered model and migration velocity field

圖7 層狀模型PLSRTM結(jié)果及其局部傾角估計Fig.7 The PLSRTM result of layered model and its local dip estimation

圖8(a)為簡單層狀模型經(jīng)Seislet分數(shù)階閾值函數(shù)約束的PLSRTM迭代10次后的結(jié)果,圖8(b)為簡單層狀模型經(jīng)Seislet軟閾值函數(shù)約束的PLSRTM迭代10次后的結(jié)果,圖8(a)、(b)都一定程度上壓制了噪音,但前者的壓制效果更好,見圖中的紅圈部分。圖8(c)、(d)分別為圖8(a)、(b)所對應(yīng)的局部傾角估計。圖8(c)和圖8(d)相比,不同顏色的顏色團幾乎消失,大部分為局部傾角值為0的顏色,則表明經(jīng)過Seislet分數(shù)階閾值函數(shù)約束和軟閾值函數(shù)約束的PLSRTM都能夠壓制串擾噪音和采集腳印。但圖8(d)在上部還存在的顏色團較大,而圖8(c)在上部附近顏色團稍小,因此經(jīng)Seislet分數(shù)階閾值函數(shù)約束PLSRTM的效果優(yōu)于經(jīng)Seislet軟閾值函數(shù)約束PLSRTM。

圖8 層狀模型下不同Seislet閾值算法約束的PLSRTM及其局部傾角估計對比Fig.8 Comparison of PLSRTM with different Seislet threshold function constraints and comparison of dip estimation in layered model

3.2 復(fù)雜模型測試

為進一步驗證本文方法的可行性,選用某一復(fù)雜模型進行成像測試。該模型尺寸為1 500 m×5 120 m,橫向采樣點數(shù)為512,縱向采樣點數(shù)為150,橫縱向采樣點距均為10 m。共激發(fā)120炮,炮間隔為40 m,檢波點為512個,檢波點距為10 m,最大采樣時間為8 000 ms,震源采用雷克子波,震源主頻率為20 Hz。

先使用時間2階,空間8階的有限差分算法及PML邊界條件進行正演模擬得到單炮記錄,再使用平面波靜態(tài)編碼技術(shù)將多炮數(shù)據(jù)合成具有不同射線參數(shù)的平面波道集,射線參數(shù)p呈線性變化,平面波道集數(shù)為11個。圖9(a)為該復(fù)雜模型的真實速度模型,圖9(b)為偏移速度場。

圖9 復(fù)雜模型速度場和偏移速度場Fig.9 Velocity Field of complex model and migration velocity field

先分別對復(fù)雜模型進行PLSRTM迭代15次和35次,得到成像結(jié)果(圖10(a)、(b))。由圖10(a)、(b),經(jīng)過15次和35次迭代的平面波最小二乘逆時偏移后,偏移結(jié)果中的構(gòu)造成像效果較好、斷層邊界較為清楚,但存在較為嚴重的串擾噪聲,這是由于多炮數(shù)據(jù)混疊在一個道集中所帶來的噪音,且淺層部分存在一定的采集腳印。迭代35次的成像結(jié)果比迭代15次的成像結(jié)果更為清晰且串擾噪音更少。圖10(c)、(d)分別是復(fù)雜模型經(jīng)過Seislet軟閾值算法約束和Seislet分數(shù)階閾值函數(shù)約束的PLSRTM迭代15次后的成像結(jié)果。由圖10(c)、(d)可知,在Seislet約束的平面波最小二乘逆時偏移后的結(jié)果中,串擾噪聲都得到了一定的壓制,對斷層邊界的成像效果都很好,但Seislet分數(shù)階閾值函數(shù)的約束效果優(yōu)于Seislet軟閾值函數(shù)的約束效果,具體見紅圈。圖10(e)、(f)分別是復(fù)雜模型經(jīng)過Seislet軟閾值函數(shù)約束和Seislet分數(shù)階閾值函數(shù)約束的PLSRTM迭代35次后的成像結(jié)果。圖10(e)在部分反射軸附近仍存在部分串擾噪音,圖10(f)的串擾噪音幾乎被壓制,淺層的采集腳印也得到一定的壓制,具體位置見紅圈。Seislet分數(shù)階閾值函數(shù)的約束效果明顯好于Seislet軟閾值函數(shù)的約束效果。

圖10 復(fù)雜模型下不同Seislet閾值函數(shù)約束的PLSRTM對比Fig.10 Comparison of PLSRTM with different Seislet threshold function constraints in complex model

雖然增加平面波道集的數(shù)量可以一定程度地提高PLSRTM的成像質(zhì)量,壓制串擾噪音,但會增加計算量。本文中先分別使用30個平面波道集的PLSRTM迭代15和35次后,得到偏移結(jié)果(圖11(a)、(c))。圖11(b)、(d)為使用11個平面波道集且經(jīng)過Seislet分數(shù)階閾值函數(shù)約束的PLSRTM迭代15和35次后的結(jié)果。由圖11(a)、(b)、(c)、(d)可知,使用30個平面波道集的PLSRTM結(jié)果(圖11(a)、(c))比未經(jīng)過約束使用11個平面波道集的PLSRTM結(jié)果(圖10(a)、(b))噪音更少,成像質(zhì)量更高,但使用11個平面波道集且采用Seislet分數(shù)階閾值函數(shù)約束的PLSRTM(圖11(b)、(d)),其串擾噪音被壓制得更為徹底,具有和使用較多平面波道集的傳統(tǒng)PLSRTM相同的成像效果。再將使用11個道集的約束PLSRTM、11個道集的傳統(tǒng)方法,以及30個道集的傳統(tǒng)方法進行運算時間對比。

由表1可知,經(jīng)過約束且使用11個平面波道集的PLSRTM時間略少于使用11個平面波道集的傳統(tǒng)PLSRTM,更少于使用30個平面波道集的傳統(tǒng)PLSRTM。因此該方法能夠一定程度地提高計算效率。

表1 運算時間對比Table 1 Calculation time comparison s

圖11 使用不同平面波道集數(shù)目下的PLSRTM與約束PLSRTM對比Fig.11 Comparison between PLSRTM and constrained PLSRTM using different number of plane-wave gather

對復(fù)雜模型進行偏移成像后,分別抽取其第150道、350道和200道作單道振幅對比(圖12)。在圖12(a)、(b)中,藍色虛線為經(jīng)過Seislet分數(shù)階閾值算法約束的PLSRTM迭代20次后的振幅,橘色實線為經(jīng)過PLSRTM迭代20次后的振幅,黃色為真實反射系數(shù)。從圖12(a)、(b)中可以看出,經(jīng)過Seislet分數(shù)階閾值算法約束的20次迭代PLSRTM結(jié)果保幅性更好,更加收斂于真實的反射系數(shù),成像更為準確,且振幅曲線波動更小,相對于未約束的PLSRTM抗噪能力更強。由圖12(c),對11個道集的Seislet分數(shù)階閾值函數(shù)約束的PLSRTM、11個道集的傳統(tǒng)方法,以及30個道集的傳統(tǒng)方法迭代第20次結(jié)果抽取第200道對比后,經(jīng)過約束的PLSRTM結(jié)果更加收斂于真實的反射系數(shù),具有更好的保幅性。

圖12 復(fù)雜模型單道振幅對比Fig.12 Single-trace amplitude comparison of complex model

為了進一步驗證本文方法對噪聲的適應(yīng)性,對所有炮集添加隨機噪聲,使其信噪比為2 dB。再對炮集進行平面波編碼,得到11個平面波道集, (圖13)。計算信噪比RSN的方法為能量疊加法[28],公式如下:

(17)

式中,M為時間采樣點;N為道數(shù);xij為第j道第i點的振幅。再分別進行Seislet分數(shù)階閾值算法約束的PLSRTM和傳統(tǒng)方法PLSRTM,并分別抽取第15次和第35次迭代結(jié)果進行分析(圖14)。由圖14可知,傳統(tǒng)方法PLSRTM結(jié)果的噪音較為嚴重,尤其在構(gòu)造附近,而經(jīng)過Seislet分數(shù)階閾值算法約束的PLSRTM后,噪音得到有效壓制,且仍然能夠較好地還原構(gòu)造。

圖13 射線參數(shù)p為0的含噪聲平面波道Fig.13 Plane-wave gather for p=0 ms/m encoded by shot data with noise

圖14 含噪聲下的PLSRTM與約束PLSRTM對比Fig.14 Comparison between PLSRTM and constrained PLSRTM using plane-wave gather with noise

4 結(jié) 論

(1)Seislet變換測試表明,地震數(shù)據(jù)在Seislet域具有較好的稀疏性,但在使用閾值法進行地震數(shù)據(jù)重構(gòu)時閾值不宜太小,以免丟失部分有效信息。

(2)數(shù)值測試中的模型測試和單道振幅測試表明,Seislet分數(shù)階閾值算法約束的PLSRTM能夠在傳統(tǒng)PLSRTM的串擾噪音較強和保幅性等方面得到改善,且Seislet分數(shù)階閾值算法比Seislet軟閾值算法性能更優(yōu)越。含噪數(shù)據(jù)測試結(jié)果表明,本文算法約束的PLSRTM能夠在炮記錄含噪的情況下得到較好的成像結(jié)果。

(3)基于Seislet分數(shù)階閾值算法約束的PLSRTM可以使用較少的平面波道集達到與傳統(tǒng)PLSRTM相同的成像效果,且運行時間也少于傳統(tǒng)PLSRTM,提高了計算效率。

(4)Seislet變換對傾角過于依賴,傾角預(yù)測不精確時容易對成像造成影響,以后可以在如何減弱Seislet對傾角依賴的方面進行改進。