傾角域逆時偏移繞射波成像方法

汪天池 劉少勇* 顧漢明 唐永杰 嚴(yán) 哲

(①中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，湖北武漢 430074；②地球內(nèi)部多尺度成像湖北省重點實驗室，湖北武漢 430074)

0 引言

根據(jù)惠更斯原理，震源激發(fā)的一次球面波遇到不連續(xù)體時，該不連續(xù)體會作為二次震源產(chǎn)生繞射波場。繞射波可以照明一次反射波難以照明的縫洞、斷層、尖滅點等不規(guī)則地質(zhì)目標(biāo)體，對該類探區(qū)地下儲層探測有重要意義[1-2]。然而，繞射波相對于反射波能量往往較低，在溶洞發(fā)育的碳酸鹽巖地區(qū)，若存在較強反射界面，常規(guī)成像剖面中繞射波能量往往會被掩蓋[3]。因此，在特定的地震地質(zhì)條件下，開發(fā)合適的繞射波成像方法，突出繞射能量有重要的應(yīng)用價值[4-6]。

通常，實現(xiàn)繞射波成像可以在數(shù)據(jù)域先進(jìn)行繞射波和反射波的分離再成像，也可以在全波場數(shù)據(jù)成像過程中對反射波能量進(jìn)行壓制而實現(xiàn)繞射波成像。

先分離后成像的方法大多基于反射波與繞射波的運動學(xué)和動力學(xué)特征差異，在數(shù)據(jù)域進(jìn)行反射波和繞射波分離，然后分別進(jìn)行處理。Landa等[7]基于繞射波與反射波在共炮檢距域的運動學(xué)特征差異提出繞射波剖面法(D-Section)，通過描述沿繞射波時距曲線分布的波場值偏離其平均值的估計標(biāo)準(zhǔn)差以增強繞射點處的振幅。Nowak[8]利用繞射波和反射波同相軸軌跡差異，通過Radon變換從地震記錄中分離出繞射波數(shù)據(jù)。劉太臣等[9]利用奇異值譜分析法壓制共炮檢距域中的反射波強能量線性信號，重構(gòu)Hankel矩陣，提取繞射波信息。

先分離后成像的方法重點在于數(shù)據(jù)域的繞射波和反射波分離，而在成像過程中進(jìn)行反射波壓制實現(xiàn)繞射波成像則更多地依賴成像算子本身。Kirchhoff積分偏移方法高效靈活，可以用來實現(xiàn)繞射波成像。Moser等[10]根據(jù)第一菲涅耳帶反射能量集中、繞射能量分散的特點，修改Kirchhoff積分偏移公式中的加權(quán)函數(shù)，從而提取繞射能量。Klokov等[11]利用平面波分解和Radon變換的方法實現(xiàn)傾角道集內(nèi)繞射波和反射波能量的分離。為達(dá)到更好的反射能量壓制效果，Klokov等[12]又在前期研究基礎(chǔ)上提出在傾角域使用相似譜分析，以拾取更準(zhǔn)確的反射同相軸頂點。Hou等[13]則將傾角道集按角度重新編碼，使擬雙曲線形態(tài)的反射波同相軸變成了噪點，然后利用混合濾波方法去除反射波能量。隨著計算技術(shù)的發(fā)展，基于波動方程的成像方法也被用于繞射波成像。Nakata等[14]提出幾何平均成像條件下的逆時偏移(Reverse Time Migration，RTM)成像，該方法通過修改互相關(guān)成像條件，將檢波器反傳波場的算術(shù)平均項更改為幾何平均項，從而提取繞射信息。Liu等[15]則根據(jù)傾角道集內(nèi)表現(xiàn)為噪點的反射波和表現(xiàn)為直線的繞射波同相軸之間的形態(tài)差異濾除反射信號，實現(xiàn)繞射波的偏移成像。

實際上，傾角道集中反射波和繞射波的運動學(xué)特征依賴于偏移算子[16-19]。研究不同偏移算子情況下繞射波在傾角域的特征，然后基于合適的成像算子進(jìn)行傾角域濾波實現(xiàn)繞射波成像。本文從傾角域中繞射波和反射波共成像點道集的形態(tài)差異入手，利用RTM和Kirchhoff積分偏移分別對特定的模型數(shù)據(jù)進(jìn)行傾角域角度道集輸出，對比兩種成像方法在繞射波成像能力方面的差異并分析原因，最后基于RTM傾角道集，利用中值濾波提取繞射波信息以達(dá)到對繞射波偏移成像的目的。

1 理論方法

1.1 傾角域RTM

RTM成像過程簡要概括為：由震源端外推波場記錄地下每個成像點不同時刻的波場值；然后對地震記錄在時間方向上進(jìn)行反向外推并存儲地下每個成像點不同時刻的波場值[20-21]；最后利用成像條件進(jìn)行成像。經(jīng)典RTM的互相關(guān)成像條件[22]可以表示為

(1)

式中：us(x,t)和ur(x,t)分別為炮點和檢波點波場；x為地下成像點的位置；tmax為最大記錄時間；IRTM(x)為逆時偏移成像結(jié)果。在正傳源端波場和反推檢波端波場時，通過計算波印廷矢量可以獲得波場的傳播方向。波印廷矢量可表示為[23-24]

(2)

式中：θs,r為入射波與反射波的夾角，即張角；θd為反射界面法向矢量與垂直矢量Z的夾角，即需要計算的界面傾角；θs，z為入射波與Z的夾角。

圖1 簡單介質(zhì)情況下地下反射界面(a)和繞射點(b)與炮檢點的幾何關(guān)系示意圖

1.2 傾角域Kirchhoff積分偏移方法

遠(yuǎn)場近似下Kirchhoff偏移公式[25]為

(6)

式中：ξ為地面點坐標(biāo)；θ為成像點處入射射線與法線的夾角；v為成像點處地震波傳播速度；r為成像點與檢波器之間的距離。

通常，Kirchhoff積分深度偏移需要進(jìn)行旅行時計算?；诼眯袝r梯度場，可以方便地得到炮點波場和檢波點的波場傳播方向矢量Ps和Pr[26-27]

(7)

(8)

式中Ts和Tr分別為炮點和檢波點旅行時場。進(jìn)而通過式(3)～式(5)可計算傾角。

1.3 中值濾波

對輸出傾角道集進(jìn)行中值濾波，達(dá)到分離繞射波的目的[15，28]

(9)

式中：Idiff為繞射波場成像結(jié)果；n為滑動窗口的總數(shù)量；mi為濾波窗口的中值；Mθd表示沿傾角θd方向的中值濾波算子，其滑動窗口范圍為[mi-n,mi+n]，通過調(diào)控n值不同程度地抑制反射波能量。將濾波后的傾角道集疊加可得到繞射波成像結(jié)果。

2 傾角域波特征分析

均勻介質(zhì)情況下反射和繞射點的自激自收示意圖如圖2所示。在幾何地震學(xué)假設(shè)下，針對地下反射界面(圖2a)，自激自收記錄與地下反射點關(guān)系為

x0=xR+zRtanα0

(10)

(11)

式中：x0為地表自激自收點M的橫坐標(biāo)；t為M點自激自收的旅行時；α0為地下界面真實傾角； (xR,zR)為地下任意反射點R的坐標(biāo)。在單反射界面情況下，反射界面深度可以表示為

zR(xR)=z1+xRtanα0

(12)

式中z1為反射界面與z軸交點的深度。由式(10)和式(12)可得

zR(x0)=(z1cosα0+x0sinα0)cosα0

(13)

將式(13)代入式(11)，可得到任意檢波點的自激自收時間

(14)

偏移過程是將地表數(shù)據(jù)剖面u(x0,t)轉(zhuǎn)換成地下成像剖面R(xm,zm),即

(15)

(16)

式中：α為偏移傾角，也即傾角道集橫軸坐標(biāo)；vm為偏移速度。將式(14)代入式(15)和式(16)，可得

(17)

(18)

聯(lián)立式(17)和式(18)消去x0，并將xm替換為x，zm替換為zα(x,α)，可得傾角道集的幾何形態(tài)曲線

(19)

對式(19)求關(guān)于α的偏導(dǎo)數(shù)，有

(20)

令式(20)等于0，可得

(21)

由式(21)可知，在偏移速度正確，即vm=v的情況下，α=α0為式(19)的一個極小值點。因此，對于當(dāng)偏移速度正確時，傾角道集中反射波是一個有極小值點的擬雙曲線，且該曲線極小值點對應(yīng)著反射界面真實的界面深度和真實的界面傾角。

均勻介質(zhì)中，在幾何地震學(xué)假設(shè)下，針對地下繞射點(圖2b)，其自激自收記錄與地下繞射點關(guān)系為

圖2 自激自收觀測下反射(a)和繞射(b)時深關(guān)系示意圖

x0=xD+zDtanβ0

(22)

(23)

式中：β0為地表自激自收點M和繞射點D(xD,zD)連線與Z的夾角。同樣將地表數(shù)據(jù)剖面u(x0,t)轉(zhuǎn)換成地下成像剖面D(xm,zm)，即

(24)

(25)

式中：β為M點和Dm點連線關(guān)于Z的夾角。聯(lián)立式(24)和式(25)，消去β0，并將xm替換為x，zm替換為zβ(x,β)，可得

(26)

由式(26)可見，在偏移速度正確且記錄點位于繞射點正上方的情況下，即vm=v、x=xD時，該式可簡寫為zβ(x,β)=zD，此時關(guān)于繞射點的傾角道集是一條水平的直線。在x≠xD情況下，繞射波傾角道集表現(xiàn)為一條沒有極值點的曲線[29]。圖3展示了一個繞射點和一個水平反射界面在傾角道集的幾何形態(tài)，其中藍(lán)線表示水平反射界面，星號表示繞射點。在圖3b展示的地表1000m處的傾角道集中，繞射波被拉平，位于真實深度處，反射波深度—傾角關(guān)系表現(xiàn)為擬雙曲線形態(tài)，曲線頂點對應(yīng)反射界面真實的傾角和深度。

圖3 典型的反射和繞射模型(a)及其對應(yīng)的傾角道集形態(tài)(b)示意圖

為驗證以上理論分析的正確性，選取一個兩層模型，網(wǎng)格數(shù)為400×400，網(wǎng)格間距為10m×10m。模型速度如圖4所示，在下層介質(zhì)中包含一高速繞射點。利用有限差分法進(jìn)行模擬，正演參數(shù)包括：子波為主頻20Hz的雷克子波，炮間距為50m，共81炮，最大炮檢距為1500m，采用雙邊接收，檢波器間距為10m。Kirchhoff積分偏移和RTM的成像結(jié)果如圖5所示，其輸出的傾角道集如圖6所示。

在圖6a和圖6b中，在繞射點正上方(CDP200)繞射點傾角道集響應(yīng)均為直線。對于反射波同相軸，Kirchhoff方法為一條擬雙曲線且能量聚焦在擬雙曲線的頂點，而RTM法則聚焦為一個點，兩者反射波能量主要聚集在真實的界面深度和傾角處。圖6c和圖6d展示了繞射點右側(cè)(CDP250)的傾角道集，其中Kirchhoff方法得到繞射波的響應(yīng)為一條沒有駐點的曲線，RTM方法則看不到繞射波的響應(yīng)。造成圖6c與圖6d存在差異原因可以歸結(jié)為二者成像的實現(xiàn)方式不同。經(jīng)典Kirchhoff積分法偏移的實現(xiàn)方式是利用輸出道的觀點，依照旅行時關(guān)系在成像域進(jìn)行數(shù)據(jù)到成像結(jié)果的投影；而RTM則基于全波波動方程，采用有限差分進(jìn)行波場外推，在成像點進(jìn)行相關(guān)成像。相比而言，后者對反射波場有更好的聚焦能力。

圖4 兩層介質(zhì)速度模型

圖5 兩層模型Kirchhoff積分偏移(a)和RTM(b)的成像結(jié)果對比

圖6 兩層模型Kirchhoff積分偏移和RTM在不同CDP處傾角道集對比a)CDP200，Kirchhoff積分偏移； (b)CDP200，RTM； (c)CDP250，Kirchhoff積分偏移； (d)CDP250，RTM

3 數(shù)值實驗

由于RTM偏移算子對反射波聚焦能力更好，本文采用RTM進(jìn)行傾角道集輸出，并在傾角域進(jìn)行中值濾波處理，以便能有效濾除傾角道集中位于真實傾角和深度位置處的反射波能量，保留繞射波能量，實現(xiàn)繞射波成像。

首先，基于上文兩層模型CDP200和CDP250處傾角道集，進(jìn)行中值濾波處理，結(jié)果如圖7所示。對比圖7與圖6b和圖6d可以看出，兩處傾角道集中反射波的點狀響應(yīng)被視作“噪點”濾除而繞射波同相軸得到有效的保留，其最后疊加成像結(jié)果如圖8所示。對比圖5和圖8，可以看出圖5中反射界面能量被濾除而繞射點能量得到保留。

然后，采用稍復(fù)雜的繞射體模型(圖9)進(jìn)一步測試本文方法效果。如圖9所示，該模型四個繞射體位于第二層介質(zhì)中，界面角點也可產(chǎn)生繞射波。模型網(wǎng)格數(shù)為600×600，網(wǎng)格間距為10m×10m。采用主頻為20Hz的雷克子波用有限差分法進(jìn)行波場模擬，雙邊接收，共121炮，炮間距為50m，最大炮檢距為1500m，檢波器間距為10m。傳統(tǒng)RTM成像和繞射波RTM成像結(jié)果如圖10所示，可以看出，在繞射波RTM成像剖面中，反射界面能量被充分壓制，藍(lán)色方框所示的繞射點的能量則得到了有效的保留。

圖7 中值濾波后CDP200(a)和CDP250(b)處傾角道集

圖8 RTM繞射波成像結(jié)果

最后，采用含繞射點的部分Sigbee2A模型驗證本文方法的有效性，其速度模型如圖11所示。該模型網(wǎng)格數(shù)為500×700，采樣間隔為25ft×25ft。對該模型模擬數(shù)據(jù)采用經(jīng)典RTM成像和RTM繞射波成像，結(jié)果如圖12所示，中值濾波前、后CDP200處傾角道集如圖13所示。由圖12可以看出，常規(guī)RTM成像繞射波波能量淹沒在反射同相軸中，在成像剖面中不夠明顯。RTM繞射波成像則使繞射波的能量得到保留，反射波的能量得到極大的削弱,突出了繞射點所在的位置。由于該模型較為復(fù)雜，在該情況下采用波印廷矢量，存在嚴(yán)重的波場交叉，可能得到不準(zhǔn)確的波場傳播方向，使得計算出來的傾角不夠準(zhǔn)確，反射波的響應(yīng)很難聚焦為“點”，而是表現(xiàn)為“短線”(圖13a)，但在該情況下反射波還是被較徹底地去除了(圖13b)。

圖9 繞射體模型

圖10 繞射體模型常規(guī)RTM(a)和繞射波RTM(b)成像結(jié)果對比藍(lán)色方框所示為繞射點位置

圖11 部分Sigbee2A速度模型

圖12 部分Sigbee2A模型常規(guī)RTM(a)和繞射波RTM(b)成像結(jié)果對比藍(lán)色方框所示為繞射點位置

當(dāng)繞射點位于模型深部，RTM低頻噪聲對繞射波成像影響有限。但對于某些構(gòu)造模型，在低頻噪聲與繞射點重疊的部分時，為消除低頻噪聲對成像剖面質(zhì)量的影響，可以根據(jù)低頻噪聲的角度特性，在互相關(guān)成像條件中添加角度加權(quán)因子[30-34]。

圖13 部分Sigbee 2A模型CDP200處中值濾波處理前(a)、后(b)傾角道集對比藍(lán)色方框內(nèi)為繞射波在傾角道集的響應(yīng)

4 結(jié)論

本文對傾角域反射波和繞射波的幾何形態(tài)進(jìn)行了理論分析，利用Kirchhoff積分偏移和RTM算子分別進(jìn)行傾角道集輸出，研究了不同成像算子在傾角域?qū)@射波的成像能力。Kirchhoff積分偏移計算的傾角道集中，反射波同相軸形態(tài)表現(xiàn)為下凹的擬雙曲線形，駐點位于實際傾角位置，繞射波同相軸則為一條拉平的直線； RTM得到的傾角道集中反射波同相軸表現(xiàn)為孤立的“噪點”，繞射波同相軸為一條直線。數(shù)值實驗結(jié)果表明，在傾角域，RTM算子對反射波有著更好的聚焦能力。在RTM輸出的傾角道集中，利用中值濾波可將反射波能量濾除而保留繞射波能量，進(jìn)而實現(xiàn)繞射波的成像。