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

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

0 引言

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

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

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

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

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

1 理論方法

1.1 傾角域RTM

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

(1)

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

(2)

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

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

1.2 傾角域Kirchhoff積分偏移方法

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

(6)

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

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

(7)

(8)

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

1.3 中值濾波

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

(9)

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

2 傾角域波特征分析

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

x0=xR+zRtanα0

(10)

(11)

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

zR(xR)=z1+xRtanα0

(12)

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

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

(13)

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

(14)

偏移過(guò)程是將地表數(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)

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

(20)

令式(20)等于0，可得

(21)

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

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

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

x0=xD+zDtanβ0

(22)

(23)

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

(24)

(25)

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

(26)

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

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

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

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

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

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

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

3 數(shù)值實(shí)驗(yàn)

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

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

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

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

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

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

圖9 繞射體模型

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

圖11 部分Sigbee2A速度模型

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

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

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

4 結(jié)論

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