三維各向異性TI介質中的P/S波快速解耦技術及在彈性波逆時偏移中的應用

張輝,尹國慶,徐珂,王志民,王海應,梁景瑞,來姝君,左佳卉,鮮成鋼,申潁浩,趙楊

1 中國石油塔里木油田公司,新疆庫爾勒 841000 2 中國石油大學(北京)非常規科學技術研究院,北京 102249

0 引言

基于聲波方程的逆時偏移,假設橫波速度為零,其實現過程簡單,被廣泛用于實際地震勘探中(Duveneck et al.,2008; Fletcher et al.,2009; Zhan et al.,2012; Liu et al.,2019).但是在實際模擬各向異性介質的波場傳播過程時,會產生菱形偽橫波,導致計算不穩定,嚴重降低逆時偏移的成像精度(Grechka et al.,2004; Zhang et al.,2009; Yoon et al.,2010).隨著多分量數據采集技術的發展,基于彈性波方程的逆時偏移方法被逐漸應用于石油天然氣勘探中.基于彈性波方程的偏移技術考慮到了PS轉換波,其包含了更多信息,例如地層的各向異性特征(Jin et al.,2010; Weibull and Arntsen,2014).綜合應用P和S波能對復雜地層特征進行更準確的成像(Caldwell,1999; Zhao and Li,2018).

盡管彈性波逆時偏移有著較大優勢,但在實際地震波傳播過程中,其耦合的P波和S波將產生串擾噪聲,尤其在各向異性介質中較為嚴重,影響了偏移的成像質量.針對這一問題,可通過求解Christoffel方程得到彈性波的偏振方向,實現P和S波場解耦(Dellinger and Etgen,1990; Tsvankin,2012; Zhao et al.,2018).基于這一思想,Yan和Sava(2009)在非均勻VTI/TTI介質中引入非穩態濾波技術,成功構建了歸一化的偏振算子.但其分解后的標量P波場和矢量S波場將產生多個PS成像剖面,使得彈性波逆時偏移無法應用在三維介質中.而另外一套是基于矢量彈性波場分解的方法(Zhang and McMechan,2010),但考慮到分解過程中介質速度變化劇烈問題,該方法只能應用在均勻介質中.為了提高復雜介質適用性,可以在空間波數域引入低秩近似策略,將矢量分解法簡化為矩陣相乘的形式(Cheng and Fomel,2014).不過需要根據低秩近似的系數來進行多次傅里葉變換,龐大的計算量使矢量彈性波場分解方法無法在三維介質中應用.

為了提高P/S波解耦技術的計算效率,可引入LU分解算法(Yang et al.,2019).但用該方法構建的偽亥姆霍茲算子,未對與彈性系數有關的加權系數進行歸一化,因此分解后的P/S波場振幅與原始波場不一致.Zuo等(2022)對上述提到的加權系數進行了歸一化修正,并將修正后的亥姆霍茲算子應用到了三維VTI介質中,成功得到具有正確相位和振幅的P/S波矢量波場.針對更為復雜的各向異性地層,例如具有傾斜對稱軸的TI介質,一些學者同樣進行了波場傳播特征以及波場分離方面的研究(郝重濤和姚陳,2007;張慶朝等,2019;左佳卉等,2022),但在三維空間域計算效率較低,難以應用于彈性波逆時偏移中.

針對上述提到的問題,本文在改進的三維VTI分解算子基礎上,利用矩陣旋轉推導適合三維TI介質的彈性波場分解算子.在實現彈性波場分解時,應用一種快速計算方法,有效提高了計算效率.最后將算子引入到彈性波逆時偏移中,利用點積成像條件,獲得了高精度的成像界面.

1 方法

1.1 三維VTI介質中波的偏振方向和分解算子

圖1b展示了在觀測坐標系(x,y,z)下,具有傾斜對稱軸的三維TI介質示意圖,其中θ代表了對稱軸與觀測坐標系垂直方向的傾角,φ代表了對稱軸在水平面投影的方位角.對比具有垂直對稱軸的三維VTI介質(圖1a),如果將觀測坐標系進行旋轉,使垂直方向沿著傾斜對稱軸方向,此時三維TI介質就可以看成是在新坐標系(x′,y′,z′)下的三維VTI介質,適用于三維VTI介質的彈性波場分解方法就可以直接擴展到三維TI介質中.

圖1 三維VTI介質和具有傾斜各向異性對稱軸的三維TI介質

首先回顧改進的三維VTI彈性波場分解方法(Zuo et al.,2022).通過求解三維VTI介質下的Christoffel方程,可以得到三個特征值:

(1)

(2)

(3)

其對應的特征向量為:

(4)

(5)

公式(4)(5)(6)分別代表了P、SH和SV波的偏振方向.其中,kx、ky和kz是x軸、y軸和z軸三個方向上的波數,vP和vS是沿著各向異性對稱軸上的縱波和橫波速度,ρ是密度,ε、γ和δ是Thomson參數,ε描述了縱波速度在垂直與水平方向的差別,γ描述了SH波和SV波在水平方向的差別,δ描述了縱波垂直方向各向異性變化的快慢.在三維VTI介質中,三個波模式的偏振方向相互垂直,所以分解算子可以統一用P波的偏振方向來表示.因此空間域的三維VTI分解算子公式可表示為

(7)

(8)

1.2 三維TI介質分解算子公式推導

根據圖1中的幾何關系,通過三維旋轉矩陣,得到新的坐標系

于是新坐標系下的空間偏導為

(9)

將式(9)代入式(7),得到適用于三維TI介質的彈性波場分解算子公式:

(10)

當θ=φ=0°時,公式(10)可以簡化為三維VTI分解公式.

1.3 三維TI介質中的矢量彈性波場分解

(11)

將式(11)代入式(4)—(6),得到三維TI介質中對應的P、SH和SV波的偏振方向:

(12)

(13)

(14)

通過波場投影公式(Zhang et al.,2023),關于P波場有a′1·U=a′1·UP、a′2·UP=0和a′3·UP=0,將偏振方向代入,化簡得到三分量P波場的表達式:

(15)

(16)

(17)

(18)

(19)

(20)

其中k′x=k′sinθ0cosφ0,k′y=k′sinθ0sinφ0和k′z=k′cosθ0,θ0是波傳播方向與坐標z′軸的夾角,φ0是波傳播方向在x′-y′面投影與坐標軸x′的夾角.直接在波數下求解式(18)—(20),需要在每一個波場傳播的時間步長中做多次三維傅里葉變換,其計算量非常巨大.因為我們將式(18)—(20)轉換到時間空間域,三分量P波場為

(21)

(22)

(23)

其中,u″x、u″y和u″z分別是Ux/ω2、Uy/ω2和Uz/ω2的逆傅里葉變換,可以在彈性波場分解前,對子波進行兩次時間積分后延拓波場得到.

SH波場投影公式為a′2·U=a′2·USH、a′1·USH=0和a′3·USH=0,將偏振方向代入,化簡得到三分量SH波場的表達式:

(24)

(25)

(26)

(27)

于是,空間域的SH波場計算公式為

(29)

(30)

(31)

SV波場的投影公式為a′3·U=a′3·USV、a′1·USV=0和a′2·USV=0,將偏振方向代入,并化簡得到三分量SV波場的表達式:

(32)

(33)

(34)

(35)

(36)

(37)

此時空間域SV波場的計算公式為

(38)

(39)

(40)

相對于傳統波數域內求取解耦算子的方法,本文的方法在空間域內可以直接求解,極大程度地提高了三維各向異性介質中的波場分解計算效率,因此使其在三維成像中的可操作性大大增強.

1.4 成像條件

(41)

(42)

(43)

(44)

在此基礎上,省去與相角有關的參數,分解后的三分量波場計算公式可以簡化為

(45)

(46)

(47)

(48)

(49)

(50)

(51)

(52)

(53)

將新坐標系下的空間偏導公式(9)代入,最終我們利用公式(45)—(53)來計算分解的P和S波場.

2 理論測試

首先構建一個均勻的三維TI模型,大小為2000 m×2000 m×2000 m,采樣間隔為10 m.彈性參數為:vP=3.15 km·s-1、vS=1.78 km·s-1、ρ=2.30 g·cm-3、ε=0.20、δ=0.20、γ=0.00、θ=28°和φ=45°.設置頻率為25 Hz的雷克子波作為爆炸源,位于模型中心,波場延拓采用時間二階和空間八階有限差分.圖2展示了原始三分量波場,箭頭分別指示了P波和S波.圖3為利用各向同性分解算子(公式(2))的結果,在P和S波場中都存在很強的串擾噪聲,如黑色箭頭所示.三維VTI分解算子考慮了各向異性的影響,得到的結果如圖4所示,依然存在串擾噪聲.圖5展示了三維TI分解算子的結果,說明只有同時考慮各向異性參數以及傾斜對稱軸,才能將矢量P和S波場分離干凈.

圖2 三維均勻TI模型中,原始三分量彈性波場

圖3 利用三維各向同性分解算子得到的P/S波場

圖5 利用三維TI分解算子得到的P/S波場

3 數值算例

將不同三維分解算子引入到彈性波逆時偏移中,比較不同算子提高偏移成像精度的有效性.建立一個兩層的各向異性模型,大小為3000 m×3000 m×2000 m,采樣間隔為10 m.為了測試在強剪切各向異性情況下,該彈性波場分解方法的有效性,設置較高的各向異性參數γ,第一層的彈性參數為:vP=3.15 km·s-1、vS=1.78 km·s-1、ρ=2.30 g·cm-3、ε=0.33、δ=0.33、γ=0.40、θ=28°和φ=20°;第二層的彈性參數為:vP=3.40 km·s-1、vS=2.00 km·s-1、ρ=2.35 g·cm-3、ε=0.35、δ=0.36、γ=0.42、θ=30°和φ=24°.設置頻率為20 Hz的雷克子波作為爆炸源.一共有85炮和8100個檢波點,深度都為100 m,且在水平面均勻分布.利用準確的兩層彈性模型生成多分量地震記錄,而與第一層彈性參數一樣的均勻介質作為偏移模型.

圖6模擬了一炮的原始波場,由于存在強剪切各向異性,SH波場呈現為橢圓,在界面處,存在反射和透射的P波.圖7展示了利用三維VTI分解算子得到的彈性波場,都存在較明顯的串擾噪聲(黑色箭頭),無法獲得分離干凈的P和S波場,特別是在SV波場中,還存在部分SH波的能量.圖8是利用三維TI分解算子得到的結果,三種彈性波場都被分離得很干凈,相比于SH波,SV波場的能量很強.

圖6 三維兩層TI模型中,原始三分量彈性波場

圖7 利用三維VTI分解算子得到的P/S波場

將推導的三維各向異性分解算子引入彈性波逆時偏移框架中,首先疊加了16炮(圖9).利用VTI波場延拓以及VTI分解算子的偏移成像結果,由于分解算子沒有考慮到傾斜各向異性對稱軸,剖面上存在明顯的偏移假象(黑色箭頭),導致水平界面的成像精度低,特別是PS剖面.經過了85炮的疊加(圖10),VTI逆時偏移可以有效壓制串擾噪聲帶來的假象,但是三維TI分解算子能獲得更精確且干凈的成像結果,且成像界面的能量更聚焦在正確的水平界面位置(黑色箭頭).

圖9 經過不同分解方法的三維16炮成像結果

圖10 經過不同分解方法的三維85炮成像結果

將彈性波場分解方法應用到更加復雜的三維TI介質彈性波逆時偏移中,針對公開的SEAM 3D Arid模型(HTI),首先建立一個大小為2000 m×2000 m×2000 m,采樣間隔同樣為10 m的部分模型.然后將該HTI模型轉換成TI,彈性參數如圖11所示,為了加強各向異性特征,分別增大了參數ε、δ和γ的值,其中各向異性對稱軸傾角和方位角通過vP和vS估計得到(圖12).最后設置頻率為20 Hz的雷克子波作為爆炸源.一共有85炮和8100個檢波點,深度都為100 m,且在水平面均勻分布.

圖11 通過SEAM 3D Arid部分模型構建的彈性參數

圖12 通過vP和vS構建的角度參數

在實際逆時偏移計算中,應用準確的模型參數獲得三維地震記錄和進行波場正傳模擬,應用平滑的模型參數進行波場反傳模擬.圖13展示了利用三維各向同性、VTI和TI分解算子得到的PP和PS結果.通過各向同性算子,偏移剖面同相軸無法很好的聚焦(黑色箭頭).相比較于三維VTI分解算子,利用三維TI算子的偏移能得到更清晰,且與原始模型結構更一致的結果,特別是在高陡界面處,能進行準確成像,并壓制了更多的假象.另外,比較PP成像界面,PS因波長更短,具有較高的成像分辨率,能對更細小的地質結構進行成像.但是,PS剖面在深層的成像效果沒有PP剖面好,這是由于P波和S波在地下傳播過程中,照明角度不一樣導致的.

圖13 經過不同三維分解方法的85炮成像結果

4 討論

本文推導的三維TI算子能有效壓制串擾噪聲,提高三維TI彈性波逆時偏移的成像精度.與以往的彈性波場分解方法相比,本文工作主要為:

(1)在三維VTI分解算子的基礎上,利用坐標旋轉,推導三維TI分解算子公式的過程簡便.另一種方法是通過Bond變換(Winterstein,1990)計算三維TI介質的剛度矩陣,求解三維TI介質的Christoffel方程.Bond矩陣旋轉增加了公式推導的難度.

(2)本文實現了空間域的彈性波場分解,在三維TI介質中,計算效率高.而基于求解泊松方程來進行波場分解的方法,例如在三維VTI介質中通過LU分解(Zuo et al.,2022),或者在三維TTI介質中通過快速算法(左佳卉等,2022),需要更多的計算量.另一種通過低秩近似的矢量彈性波場分解方法,計算量在三維各向異性介質中同樣較大(Cheng and Fomel,2014).

本文引用了空間域分解P和S波場的快速方法,且采用的成像條件公式能得到正確的相位和走時.在彈性波逆時偏移中,高效進行真振幅成像,一直是研究難點.在未來的進一步研究中,將把真振幅恢復融入到空間域彈性波分解中.

5 結論

基于對各向異性Christoffel方程進行特征分析,本文提出了一種適用于三維TI介質的矢量彈性波場分解方法,并將其應用于彈性波逆時偏移.通過旋轉觀測坐標,推導了與各向異性對稱軸傾角和方位角有關的三維TI分解算子的表達式,并在空間域實現快速的矢量彈性波場分解.本文推導的波場分解算子,能有效提高三維各向異性彈性波逆時偏移的成像精度.