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

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

0 引言

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

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

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

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

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

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

1 方法原理

1.1 三維各向異性彈性波動(dòng)方程特征分析

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

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

(1)

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

(2)

(3)

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

通過特征分析，式(3)的特征值及對(duì)應(yīng)的特征向量(詳細(xì)推導(dǎo)見附錄A)為

(4)

(5)

(6)

(7)

(8)

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

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

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

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

(9)

(10)

(11)

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

(12)

(13)

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

圖1 三維TTI介質(zhì)對(duì)稱軸與觀測(cè)坐標(biāo)系的關(guān)系示意圖

(14)

新坐標(biāo)系和觀測(cè)坐標(biāo)系之間空間偏導(dǎo)數(shù)的關(guān)系為

(15)

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

(16)

三維TTI分解算子在新坐標(biāo)系下可以近似為三維VTI分解算子。當(dāng)θ=0°時(shí)，三維TTI分解算子

1.3 三維TTI彈性波場(chǎng)分解的實(shí)現(xiàn)

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

(17)

式中w是一個(gè)輔助波場(chǎng)[24]。

(18)

則P波和S波波場(chǎng)可表示為

(19)

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

(20)

(21)

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

(22)

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

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

W=-f(ε,δ)U

(23)

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

r4[kz2cos2θ+ky2sin2θ-kykzsin2θ]

(24)

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

(25)

(26)

(27)

(28)

(29)

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

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

(30)

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

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

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

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

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

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

(5)利用式(19)計(jì)算P波和S波波場(chǎng)。

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

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

2 數(shù)值算例

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

圖3 三維弱各向異性均勻模型正演彈性波場(chǎng)

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

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

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

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

圖7 三維強(qiáng)各向異性均勻模型正演彈性波場(chǎng)

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

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

為了驗(yàn)證本文推導(dǎo)的三維TTI介質(zhì)分解方法在層狀介質(zhì)中的有效性，建立三維兩層TTI介質(zhì)模型，對(duì)于反射波和透射波，采用兩種分解算子，比較應(yīng)用效果。第一層彈性參數(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模型的正演波場(chǎng)，既有反射PP、PS、SP和SS波(黑色字母所示)，還有透射PP、PS、SP和SS波(紅色字母所示)。由于各向異性對(duì)稱軸并沒有沿著垂直方向，而且第一層和第二層的對(duì)稱軸傾角不同，所以反射和透射波場(chǎng)既不左右對(duì)稱，也不沿界面上下對(duì)稱。圖11為利用VTI分解算子得到的P波和S波波場(chǎng)，反射和透射彈性波場(chǎng)中，都存在能量較強(qiáng)的串?dāng)_噪聲(黑色箭頭所示)。圖12為三維TTI介質(zhì)分解算子得到的分解結(jié)果，由于考慮了各向異性對(duì)稱軸的傾角，P波波場(chǎng)只包含干凈的反射和透射PP、SP波，S波波場(chǎng)僅有干凈的反射和透射PS、SS波。通過以上的數(shù)值算例，說明本文推導(dǎo)的三維TTI分解算子可以有效應(yīng)用于層狀各向異性介質(zhì)彈性波場(chǎng)。

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

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

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

3 討論

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

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

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

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

4 結(jié)束語

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

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

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

(A-1)

(A-2)

(A-3)

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

(A-4)

(A-5)

(A-6)

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

為了求解特征值λ對(duì)應(yīng)的特征向量，將等式|λ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)

則對(duì)應(yīng)于特征值λ1的特征向量為

(A-13)

同理可以獲得對(duì)應(yīng)于特征值λ2和λ3的特征向量分別為

(A-14)

(A-15)