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

汪天池 劉少勇* 顧漢明 唐永杰 嚴 哲

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

0 引言

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

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

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

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

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

1 理論方法

1.1 傾角域RTM

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

(1)

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

(2)

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

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

1.2 傾角域Kirchhoff積分偏移方法

遠場近似下Kirchhoff偏移公式[25]為

(6)

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

通常，Kirchhoff積分深度偏移需要進行旅行時計算。基于旅行時梯度場，可以方便地得到炮點波場和檢波點的波場傳播方向矢量Ps和Pr[26-27]

(7)

(8)

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

1.3 中值濾波

對輸出傾角道集進行中值濾波，達到分離繞射波的目的[15，28]

(9)

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

2 傾角域波特征分析

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

x0=xR+zRtanα0

(10)

(11)

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

zR(xR)=z1+xRtanα0

(12)

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

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

(13)

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

(14)

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

(15)

(16)

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

(17)

(18)

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

(19)

對式(19)求關于α的偏導數，有

(20)

令式(20)等于0，可得

(21)

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

均勻介質中，在幾何地震學假設下，針對地下繞射點(圖2b)，其自激自收記錄與地下繞射點關系為

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

x0=xD+zDtanβ0

(22)

(23)

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

(24)

(25)

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

(26)

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

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

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

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

圖4 兩層介質速度模型

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

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

3 數值實驗

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

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

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

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

圖8 RTM繞射波成像結果

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

圖9 繞射體模型

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

圖11 部分Sigbee2A速度模型

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

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

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

4 結論

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