三維TTI介質(zhì)中的彈性波場分解

左佳卉 張樂樂 帥 達 朱成宏 徐蔚亞 趙 楊*

(①頁巖油氣富集機理與有效開發(fā)國家重點實驗室，北京 100083；②中國石化彈性波理論與探測技術(shù)重點實驗室，北京 100083；③中國石油大學(北京)非常規(guī)油氣科學技術(shù)研究院，北京 102249；④中國石化石油勘探開發(fā)研究院，北京 100083)

0 引言

逆時偏移不受傾角和偏移孔徑的限制，在復雜地層結(jié)構(gòu)中，可比較精確地保存地震波振幅和相位等信息。隨著三分量寬方位數(shù)據(jù)采集技術(shù)的不斷發(fā)展，彈性逆時偏移成像技術(shù)也逐漸應用于實際數(shù)據(jù)。該方法考慮了轉(zhuǎn)換波的特征，能夠獲得反映特殊結(jié)構(gòu)(例如氣藏)等更豐富的有效信息[1-4]；通過求解彈性波方程，重建P波和S波波場，然后運用成像條件分別得到PP波、PS波、SP波和SS波成像結(jié)果[5-6]。但是，地震波在實際地層中傳播時，耦合的P波和S波波場會產(chǎn)生串擾噪聲，嚴重降低成像剖面的質(zhì)量，所以在應用成像條件之前，需要準確分離P波和S波波場。

地面多分量地震記錄可通過矢量合成法計算P波、S波的矢量方向[7]，把地震波振幅旋轉(zhuǎn)到真正的縱、橫波矢量軸，實現(xiàn)P波、S波波場分離[8-10]。另外，也可以在逆時偏移波場延拓過程中進行波場分離，獲得純P波、純S波波場。基于Helmholtz理論，通過散度和旋度算子得到標量和矢量勢，然后計算P波、S波波場[11]；也可以直接對散度和旋度算子計算的勢進行成像[12]。另外，采用P波、S波解耦的彈性波方程，可以得到三維介質(zhì)中分離的P波、S波波場[13]。通過Poynting矢量求取的質(zhì)點振動方向并進行標量化，可得到標量PP波和PS波成像結(jié)果，但是該方法只適用于各向同性介質(zhì)。

實際地下介質(zhì)普遍呈現(xiàn)各向異性特征，通常由頁巖層、定向排列的垂直裂縫以及具有周期性的薄互層引起[14]，用橫向各向同性即TI介質(zhì)描述。若各向異性對稱軸沿垂直方向，則稱為VTI介質(zhì)。當?shù)貙咏?jīng)受了強烈的構(gòu)造運動，其對稱軸與垂直方向形成一定夾角時，TTI介質(zhì)是更精確的近似[15]。面對各向異性介質(zhì)，以上分離方法不再適用，因為此時彈性波的偏振方向不再與傳播方向平行或垂直，旋度和散度算子會造成不正確的投影，在彈性波場中產(chǎn)生不同波型能量殘留[16]，導致錯誤的成像結(jié)果。所以需要發(fā)展適合三維TTI介質(zhì)的高效精確的彈性波場解耦方法，以提高彈性波逆時偏移的成像精度。

通過求解 Christoffel 方程，可以計算各向異性介質(zhì)中波的偏振方向，進而對彈性波場進行分離[17]。很多學者都沿著這個思路對各向異性介質(zhì)中的波場分離方法進行了研究，如利用非穩(wěn)態(tài)濾波構(gòu)建歸一化的偏振方向，可以在非均勻VTI介質(zhì)中實現(xiàn)P波和S波波場的分離[18]。考慮到地震波的矢量特征，可以在波數(shù)域?qū)崿F(xiàn)VTI介質(zhì)的彈性波場分解，且分解前后的波場具有相同的振幅、相位和量綱[19]。引入低秩近似思想后，矢量波場分解法的計算過程可以簡化為混合域(空間波數(shù)域)的矩陣乘法形式[20]。另外，在波數(shù)域?qū)TI介質(zhì)Christoffel方程進行特征分析，可以計算不同波模式的偏振方向，然后構(gòu)建合適的矢量彈性波場分解算子[21]。但是，這些適用于VTI模型的分解算子無法考慮傾斜對稱軸的情況，不能準確估計TTI介質(zhì)中波的偏振方向。直接求解TTI介質(zhì)的 Christoffel 方程，可以計算與各向異性對稱軸傾角有關的偏振方向，實現(xiàn)TTI介質(zhì)的彈性波場分解[22]。但是，分解前、后的波場具有不同的振幅、相位和量綱。依據(jù)矢量彈性波場分解，有學者提出在波數(shù)域?qū)ΧS Christoffel 方程進行特征分析，通過構(gòu)建VTI分解算子以及坐標旋轉(zhuǎn)，得到TTI分解算子，最終實現(xiàn)二維TTI介質(zhì)彈性波場分解[23]。該方法應用矩陣分解算法提高了計算效率，但是分解前、后的波場振幅不一致，當對稱軸傾角以及各向異性參數(shù)都為零時，分解算子無法退化到各向同性的形式。

目前針對各向異性介質(zhì)的彈性波場分解方法，主要都應用于二維介質(zhì)，面對更復雜的三維TTI介質(zhì)，分解方法主要還是基于散度和旋度算子，通過投影得到分離的標量P波場和矢量S波場。但是，在應用成像條件時，分離后的標量和矢量波場將產(chǎn)生三個PS波成像結(jié)果。而且在PS波成像界面中會出現(xiàn)極性反轉(zhuǎn)的問題，如果不進行適當校正，會導致非相干疊加，產(chǎn)生假象。

為了解決三維TTI介質(zhì)彈性波場分解面臨的問題，本文首先從VTI彈性波動方程出發(fā)，對 Christoffel 方程進行特征分析，再推導出三維VTI分解算子；然后根據(jù)各向異性對稱軸與觀測坐標系之間的關系，利用坐標旋轉(zhuǎn)導出了三維TTI介質(zhì)波場分解算子；最后構(gòu)建泊松方程，實現(xiàn)三維TTI介質(zhì)彈性波場分解。為了避免直接求解泊松方程，引入了改進的快速算法計算輔助波場。該三維TTI波場分解算子考慮了隨空間變化的彈性參數(shù)和對稱軸傾角，適用于復雜的各向異性介質(zhì)。對均勻和非均勻的三維TTI介質(zhì)模型進行彈性波波場分解的數(shù)值模擬結(jié)果表明，該分解算子能夠有效分離三維TTI介質(zhì)中的彈性波場。

1 方法原理

1.1 三維各向異性彈性波動方程特征分析

對各向異性介質(zhì)的彈性波場進行特征分析，其特征值以及對應的特征向量分別表示不同的波模式。通過求解VTI介質(zhì) Christoffel 方程的特征向量，可以構(gòu)建分解算子，實現(xiàn)VTI介質(zhì)的彈性波場分解。針對TTI介質(zhì)，將觀測坐標系進行旋轉(zhuǎn)，使垂向坐標軸沿著傾斜對稱軸方向，此時可近似為VTI介質(zhì)，于是適用于VTI介質(zhì)的彈性波場分解方法就可以推廣到TTI介質(zhì)。

根據(jù)牛頓第二定理和廣義胡克定律，三維VTI二階彈性波動方程為

(1)

式中：ux、uy和uz表示位移三分量；C11、C13、C33、C44和C66為剛度矩陣中的獨立元素；ρ為密度。在頻率—波數(shù)域中，將式(1)寫成矩陣形式

(2)

(3)

式中：ρω2是矩陣A的特征值，可以計算得到不同波型的相速度ω/k；U是矩陣A的特征向量，表示不同波的偏振方向。

通過特征分析，式(3)的特征值及對應的特征向量(詳細推導見附錄A)為

(4)

(5)

(6)

(7)

(8)

式中vP和vS分別是沿各向異性對稱軸方向的P波和S波速度。

在三維VTI介質(zhì)中，包含各向異性對稱軸的平面都是對稱面，所以對任意方向入射的地震波，P波和SV波的偏振方向都在對稱面內(nèi)，而SH波的偏振方向在各向同性面內(nèi)，與彈性參數(shù)無關。

1.2 構(gòu)建三維TTI分解算子

根據(jù)傅里葉變換的性質(zhì)，將式(5)～式(7)從波數(shù)域轉(zhuǎn)換到空間域，有

(9)

(10)

(11)

根據(jù)計算的偏振方向，將對應的波模式投影到各自的偏振方向上，實現(xiàn)彈性波場分解。在三維VTI介質(zhì)中，P波、SV波和SH波的偏振方向互相垂直，于是可以利用P波的偏振方向，構(gòu)建空間域的分解算子為

(12)

(13)

在此基礎上，根據(jù)三維TTI介質(zhì)各向異性對稱軸與觀測坐標系之間的幾何關系，構(gòu)建一個旋轉(zhuǎn)矩陣，使垂向坐標軸沿著傾斜各向異性對稱軸方向。定義θ為三維TTI介質(zhì)對稱軸在yOz平面內(nèi)沿z軸順時針方向的夾角(圖1)，三維坐標旋轉(zhuǎn)公式為

圖1 三維TTI介質(zhì)對稱軸與觀測坐標系的關系示意圖

(14)

新坐標系和觀測坐標系之間空間偏導數(shù)的關系為

(15)

代入式(12)，可得三維TTI分解算子

(16)

三維TTI分解算子在新坐標系下可以近似為三維VTI分解算子。當θ=0°時，三維TTI分解算子

1.3 三維TTI彈性波場分解的實現(xiàn)

根據(jù)Helmholtz定理，任何一個彈性波場u都可以分解為一個無旋場uP和一個無散場uS。矢量分解將彈性波場分別投影到各自的偏振方向上，實現(xiàn)波場分離。通過構(gòu)建泊松方程u=2w，可以在空間域?qū)崿F(xiàn)各向同性介質(zhì)波場分離

(17)

式中w是一個輔助波場[24]。

(18)

則P波和S波波場可表示為

(19)

將式(18)展開并整理，有

(20)

(21)

式(20)在頻率—波數(shù)域可表示為

(22)

式中W分別為w的頻率—波數(shù)域形式。

利用式(22)，通過Gaussian消元法可以直接在頻率域內(nèi)計算輔助波場，但是系數(shù)r1、r2和r3包含了隨空間變化的彈性參數(shù)，所以波場分解只適用于平滑的各向異性模型[21]。為了在復雜各向異性介質(zhì)里實現(xiàn)彈性波場分解，Zhang等[25]提出一種在二維VTI介質(zhì)中快速求解輔助波場的方法，可以在空間域計算系數(shù)r1、r2和r3，以及構(gòu)建w，本文將該快速算法擴展到三維TTI介質(zhì)。式(22)可寫為

W=-f(ε,δ)U

(23)

基于弱各向異性介質(zhì)的假設，將f在ε=δ=0處做Taylor展開，有

r4[kz2cos2θ+ky2sin2θ-kykzsin2θ]

(24)

W=-U1-r4U2cos2θ-r4U3sin2θ+r4U4sin2θ

(25)

(26)

(27)

(28)

(29)

對式(25)做反傅里葉變換，有

w=-u1-r4u2cos2θ-r4u3sin2θ+r4u4sin2θ

(30)

式中u1、u2、u3和u4分別是U1、U2、U3和U4在時間—空間域的形式。式(30)即為空間域計算輔助波場的公式，適用于三維復雜TTI介質(zhì)。

三維TTI介質(zhì)中彈性波場分解的具體計算步驟為：

(1)對三分量原始波場u求三維傅立葉變換得到U；

(2)在頻率波數(shù)域計算U1、U2、U3和U4；

(3)做反傅立葉變換得到u1、u2、u3和u4；

(4)在空間域計算含傾角θ的三角函數(shù)和參數(shù)r4，用式(30)計算空間域的三維輔助波場；

(5)利用式(19)計算P波和S波波場。

VTI介質(zhì)是特殊的TTI介質(zhì)，不用考慮對稱軸傾角，輔助波場的求解過程大大簡化(圖2)。

圖2 三維TTI與VTI介質(zhì)彈性波場分解過程的對比

2 數(shù)值算例

首先建立均勻TTI介質(zhì)模型，對比各向同性和各向異性分解方法的結(jié)果，驗證本文三維TTI彈性波場分解算子的有效性。模型的網(wǎng)格數(shù)為200×200×200，空間采樣間隔均為10m。在模型中心設置一個主頻為25Hz的Ricker子波作為爆炸源。為了滿足弱各向異性介質(zhì)的假設，設置三維TTI介質(zhì)的彈性參數(shù)為：vP=3.25km/s，vS=1.90km/s，ρ=2.00g/cm3，ε=0.22，δ=0.22，γ=0，θ=30°。

圖3 三維弱各向異性均勻模型正演彈性波場

圖4 三維弱各向異性均勻模型正演波場的各向同性算子分解結(jié)果

圖5 三維弱各向異性均勻模型正演波場的VTI算子分解結(jié)果

圖6 三維弱各向異性均勻模型正演波場的TTI算子分解結(jié)果

當?shù)卣鸩ㄔ趯嶋H三維各向異性介質(zhì)中傳播時，除了P波和SV波，還存在SH波。為了體現(xiàn)強橫波各向異性特征，令上述三維TTI均勻介質(zhì)的γ=0.5。在正演的三分量彈性波場中，藍色箭頭所示為SH波場(圖7)。VTI分解算子不能將彈性波場投影到實際的偏振方向(圖8)，分解的P波和S波波場中都存在串擾噪聲(黑色箭頭所示)。只有同時考慮了各向異性和傾斜對稱軸的影響，反映了真實彈性波偏振方向的分解算子才能壓制串擾噪聲(圖9)。比較不同算子的分解結(jié)果可見，對于三維TTI均勻介質(zhì)的彈性波場，本文推導的三維TTI介質(zhì)分解算子能更有效地分離P波和S波波場。

圖7 三維強各向異性均勻模型正演彈性波場

圖8 三維強各向異性均勻模型正演波場的VTI算子分解結(jié)果

圖9 三維強各向異性均勻模型正演波場的TTI算子分解結(jié)果

為了驗證本文推導的三維TTI介質(zhì)分解方法在層狀介質(zhì)中的有效性，建立三維兩層TTI介質(zhì)模型，對于反射波和透射波，采用兩種分解算子，比較應用效果。第一層彈性參數(shù)為：vP=3.25km/s，vS=1.90km/s，ρ=2.00g/cm3，ε=0.22，δ=0.22，γ=0，θ=10°；第二層彈性參數(shù)為：vP=3.65km/s，vS=2.50km/s，ρ=2.35g/cm3，ε=0.24，δ=0.23，γ=0，θ=15°。界面深度為1000m，三分量位移震源的深度為200m。

圖10為三維兩層TTI模型的正演波場，既有反射PP、PS、SP和SS波(黑色字母所示)，還有透射PP、PS、SP和SS波(紅色字母所示)。由于各向異性對稱軸并沒有沿著垂直方向，而且第一層和第二層的對稱軸傾角不同，所以反射和透射波場既不左右對稱，也不沿界面上下對稱。圖11為利用VTI分解算子得到的P波和S波波場，反射和透射彈性波場中，都存在能量較強的串擾噪聲(黑色箭頭所示)。圖12為三維TTI介質(zhì)分解算子得到的分解結(jié)果，由于考慮了各向異性對稱軸的傾角，P波波場只包含干凈的反射和透射PP、SP波，S波波場僅有干凈的反射和透射PS、SS波。通過以上的數(shù)值算例，說明本文推導的三維TTI分解算子可以有效應用于層狀各向異性介質(zhì)彈性波場。

圖10 三維兩層TTI介質(zhì)模型正演彈性波場

圖11 三維兩層TTI介質(zhì)模型正演波場的VTI算子分解結(jié)果

圖12 三維兩層TTI介質(zhì)模型正演波場的TTI算子分解結(jié)果

3 討論

為了解決三維TTI介質(zhì)彈性波場的分解問題，本文利用三維VTI分解算子和傾斜對稱軸與觀測坐標系之間的幾何關系，推導了適用于三維TTI介質(zhì)的分解算子。相比于前人提出的波場分解方法，本文的方法在物理意義和數(shù)值穩(wěn)定性方面均具有一定優(yōu)勢。

(1)目前已實現(xiàn)的二維TTI各向異性波場分解方法中，已有研究學者通過求解 Christoffel 方程的特征向量計算各向異性分解算子[23]。但是當θ=0°和ε=δ=0時，無法退化為各向同性介質(zhì)的形式。而本文推導的三維TTI分解算子，當介質(zhì)退化為各向同性時，其表達式與彈性參數(shù)無關，只與傳播方向有關，更符合實際地震波傳播特征。

(2)在適用于三維TTI介質(zhì)的彈性波場分解方法中，所有基于散度和旋度的波場分離方法只能得到標量P波波場和矢量S波波場，這會導致PS波在成像界面出現(xiàn)極性反轉(zhuǎn)[22]，需要進一步對波場的振幅和相位進行校正。本文推導的三維TTI分解算子能同時計算矢量P波和S波波場，可以避免這一的問題。

本文的波場分離方法只針對波場傳播過程中，三維各向異性介質(zhì)P波、SV波和SH波極化方向垂直的情況，面對非正交快、慢橫波的分離，還具有一定的局限性。可以將本文的分離方法引入到彈性波逆時偏移，但由于三維TTI介質(zhì)分解算子包含了對稱軸傾角參數(shù)，其波場分解的數(shù)值計算比VTI介質(zhì)算子更復雜。在每一時刻求解輔助波場時，需要進行兩次反傅里葉變換，導致基于三維TTI分解算子的彈性波逆時偏移需要更多的計算時間、計算量和內(nèi)存。為解決這一問題，今后考慮將GPU并行引入三維TTI介質(zhì)彈性波逆時偏移，可極大地提高計算效率，實現(xiàn)適用于復雜三維各向異性介質(zhì)模型的、基于彈性波場分解的逆時偏移。

4 結(jié)束語

本文通過對波數(shù)域的 Christoffel 方程進行特征分析，計算了不同波的偏振方向，構(gòu)建適用于三維VTI介質(zhì)的分解算子。在此基礎上，根據(jù)傾斜對稱軸與觀測坐標系之間的幾何關系，利用坐標旋轉(zhuǎn)推導了適用于三維TTI介質(zhì)的分解算子。為了提高計算效率，避免直接求解泊松方程，引入了快速算法計算輔助波場，最終實現(xiàn)三維TTI介質(zhì)中矢量彈性波場的分解。數(shù)值算例的結(jié)果說明，本文推導的空間域三維TTI分解算子能有效實現(xiàn)三維復雜TTI介質(zhì)中的彈性波場的縱、橫波分離。

附錄A 三維VTI介質(zhì)特征值和特征向量的導出

為了計算矩陣A的特征值λ，構(gòu)建等式|λE-A|=0(其中E為單位矩陣)，其三個解為

(A-1)

(A-2)

(A-3)

根據(jù)剛度矩陣和Thomsen參數(shù)之間的關系(式(8))，式(A-1)、式(A-2)可寫為

(A-4)

(A-5)

(A-6)

基于弱各向異性假設，將根號項用泰勒展開，可得式(4)的特征值。

為了求解特征值λ對應的特征向量，將等式|λE-A|=0展開，有

(C11-C66)kxkyUy+(C13+C44)kxkzUz

(A-7)

(A-8)

λUz=(C13+C44)kxkzUx+(C13+C44)kykzUy+

(A-9)

將λ1代入式(A-8)方程，且令Ux=Uy，可得

(A-10)

將式(A-10)代入式(A-9)，整理后有

(A-11)

(A-12)

則對應于特征值λ1的特征向量為

(A-13)

同理可以獲得對應于特征值λ2和λ3的特征向量分別為

(A-14)

(A-15)