TTI介質(zhì)qP波偽譜法正演模擬

張慶朝 朱國(guó)維 周俊杰 朱聰聰 劉衛(wèi)剛

(中國(guó)礦業(yè)大學(xué)(北京)煤炭資源與安全開采國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100083)

0 引言

為了簡(jiǎn)便起見，在地震勘探中通常把地下介質(zhì)視為各向同性介質(zhì)。但研究結(jié)果表明，介質(zhì)的各向異性性質(zhì)廣泛存在。人們根據(jù)介質(zhì)性質(zhì)，又將各向異性介質(zhì)細(xì)分為不同類型，常見的有橫向各向同性介質(zhì)(Transversely Isotropic Media,TI)。具有垂直對(duì)稱軸的TI介質(zhì)稱為VTI介質(zhì)，如果對(duì)稱軸是傾斜的，則稱為TTI介質(zhì)。正演模擬是認(rèn)識(shí)地震波傳播規(guī)律的有效手段，也是實(shí)現(xiàn)逆時(shí)偏移(Reverse Time Migration、RTM)、波動(dòng)方程層析成像以及全波形反演(Full Waveform Inversion，F(xiàn)WI)的核心內(nèi)容。各向異性全彈性波動(dòng)方程可準(zhǔn)確、全面地描述地震波在地層中的傳播情況，但是全彈性波正演、成像與反演還面臨諸多困難，目前TI介質(zhì)正演模擬、RTM以及參數(shù)反演主要采用簡(jiǎn)化的qP波方程[1]。

與各向同性介質(zhì)相比，TI介質(zhì)是有5個(gè)獨(dú)立的彈性系數(shù)，波動(dòng)方程形式更復(fù)雜。在垂直傳播面內(nèi)，平行對(duì)稱軸和垂直對(duì)稱軸兩個(gè)方向的縱波(P波)和橫波(SV波)是解耦的，其他傳播方向的P波和SV波是耦合的。另外，TI介質(zhì)的P波極化方向和傳播方向不一致，SV波極化方向不垂直于傳播方向，因此又分別稱為準(zhǔn)縱波(qP波)和準(zhǔn)橫波(qSV波)。由于qP波和qSV波耦合，在TI介質(zhì)波動(dòng)方程數(shù)值模擬中分離P、S波存在困難。為了得到獨(dú)立傳播的qP波和qSV波方程，一般需要對(duì)介質(zhì)做近似處理，如小傾角近似、橢圓近似、弱各向異性近似、聲學(xué)近似等[2]。

聲學(xué)近似簡(jiǎn)單、方便，計(jì)算量小，在RTM成像中發(fā)揮著重要作用。Alkhalifah[3-4]提出了聲學(xué)近似，將沿著對(duì)稱軸方向的S波速度設(shè)置為零，從Thomsen參數(shù)表示的頻散關(guān)系出發(fā)，推導(dǎo)了四階qP-qSV波耦合波動(dòng)方程，方程中的四階混合偏導(dǎo)數(shù)增加了數(shù)值模擬的復(fù)雜度。Zhou等[5]基于相同的頻散關(guān)系，通過引入輔助函數(shù)，導(dǎo)出了qP-qSV波耦合的二階波動(dòng)方程，數(shù)值模擬也更易實(shí)現(xiàn)。隨后Zhou等[6]把qP-qSV波耦合的二階波動(dòng)方程推廣到TTI介質(zhì)。Zhang等[7]基于聲學(xué)近似假設(shè)，推導(dǎo)了TTI介質(zhì)的四階P波波動(dòng)方程。Zhang等[8]、Zhang等[9]實(shí)現(xiàn)了三維TTI介質(zhì)RTM。Fletcher等[10]把VTI介質(zhì)的偽聲波波動(dòng)方程擴(kuò)展到TTI介質(zhì)。梁鍇等[11]導(dǎo)出了三維TTI介質(zhì)聲學(xué)近似qP波波動(dòng)方程。雖然基于聲學(xué)近似的偽P波波動(dòng)方程可以精確地保留P波運(yùn)動(dòng)學(xué)特征，但其P、S波并非完全解耦，容易產(chǎn)生兩個(gè)問題，一是偽橫波噪聲干擾，二是復(fù)雜介質(zhì)條件下方程的解不穩(wěn)定。在各向異性介質(zhì)中，地震波的傳播速度是傳播方向的函數(shù)，僅把對(duì)稱軸方向的S波速度人為設(shè)置為零，以這種方式推導(dǎo)出的偽聲波方程并不能真正地消除S波的影響，所以嚴(yán)格意義上不應(yīng)稱為聲波波動(dòng)方程。當(dāng)震源位于VTI或TTI介質(zhì)中時(shí)，在數(shù)值模擬時(shí)會(huì)出現(xiàn)S波，另外在界面處也會(huì)產(chǎn)生轉(zhuǎn)換S波。在橢圓型各向異性介質(zhì)中不會(huì)激發(fā)S波，實(shí)際介質(zhì)中橢圓型各向異性介質(zhì)很少見。各向異性參數(shù)匹配法是壓制偽橫波噪聲的一種有效手段，在聲學(xué)假設(shè)條件下，通過在震源附近設(shè)置各向同性層或者橢圓各向異性層消除偽橫波噪聲的影響，但在各向異性參數(shù)變化劇烈的區(qū)域仍然會(huì)產(chǎn)生S波干擾。Yoon等[12]在陡傾區(qū)域設(shè)置橢圓各向異性層減少S波干擾，將劇烈變化的模型各向異性參數(shù)簡(jiǎn)化為連續(xù)變化，也能改善模擬效果，但在改變模型局部特征的同時(shí)，也改變了波場(chǎng)的運(yùn)動(dòng)學(xué)特征。張巖等[13]借助于輔助場(chǎng)，給出了一種空間域的簡(jiǎn)單、高效的壓制偽橫波噪聲的方法。Fowler等[14]推導(dǎo)了對(duì)稱軸方向S波速度不為零的VTI介質(zhì)二階耦合波動(dòng)方程。Fletcher等[15]通過把對(duì)稱軸方向的S波速度設(shè)置為非零值，在一定程度上解決了不穩(wěn)定性問題，但會(huì)增強(qiáng)偽橫波噪聲干擾強(qiáng)度，隨著遞推時(shí)間增加，數(shù)值頻散也變得更嚴(yán)重。

一般通過直接推導(dǎo)qP-qSV波完全解耦的純P波方程解決由聲學(xué)近似帶來的數(shù)值不穩(wěn)定問題。Du等[16-17]基于弱各向異性近似和平方根近似，推導(dǎo)了TTI介質(zhì)時(shí)間—波數(shù)域純qP波解耦方程。Liu等[18]推導(dǎo)了VTI介質(zhì)時(shí)間—空間域純qP波波動(dòng)方程。Pestana等[19-20]基于快速展開法(Rapid Expansion Method，REM)實(shí)現(xiàn)了VTI介質(zhì)時(shí)間—波數(shù)域純qP波RTM成像。Chu等[21]利用泰勒展開公式，推導(dǎo)了TTI介質(zhì)解耦純qP波波動(dòng)方程。Zhan等[22]利用REM的混合有限差分與偽譜法實(shí)現(xiàn)了TTI介質(zhì)純qP波RTM。王偉國(guó)等[23]實(shí)現(xiàn)了TTI介質(zhì)qP波偽譜法RTM。黃建平等[24]基于偽譜法實(shí)現(xiàn)了TTI介質(zhì)一階qP波方程正演模擬。黃金強(qiáng)等[1]借助Low-rank分解，實(shí)現(xiàn)了一種間接的純qP波波場(chǎng)外推方案。郭成鋒等[25]利用偽譜法和有限差分法聯(lián)合實(shí)現(xiàn)了純qP波波場(chǎng)延拓，提高了計(jì)算效率。

本文從Tsvankin[26]的VTI介質(zhì)qP波精確相速度表達(dá)式出發(fā)，引入單位向量，利用坐標(biāo)變換的方法把qP波的精確相速度擴(kuò)展到三維TTI介質(zhì)?；谌醺飨虍愋越坪推椒礁?，推導(dǎo)了三維TTI介質(zhì)純qP波近似相速度和頻散關(guān)系。根據(jù)頻散關(guān)系構(gòu)建了純qP波時(shí)間—波數(shù)域的波場(chǎng)遞推格式，并進(jìn)行了相應(yīng)的穩(wěn)定性分析。最后利用偽譜法對(duì)不同的TI模型進(jìn)行了qP波數(shù)值模擬。

1 方法原理

1.1 TTI介質(zhì)的相速度

對(duì)于VTI介質(zhì)，Tsvankin[26]給出了qP波精確相速度表達(dá)式

(1)

當(dāng)介質(zhì)的ε和δ遠(yuǎn)小于1時(shí)，稱為弱各向異性介質(zhì)，實(shí)際介質(zhì)大多為弱各向異性介質(zhì)。因?yàn)槭?1)較復(fù)雜，將式中的根式進(jìn)行泰勒展開，根據(jù)弱各向異性假設(shè)，舍去ε和δ的二次及以上高次項(xiàng)，得到線性近似的相速度表達(dá)式[26]

(2)

式(2)為純qP波相速度表達(dá)式，無S波項(xiàng)。

對(duì)于TTI介質(zhì)，將式(1)、式(2)進(jìn)行坐標(biāo)旋轉(zhuǎn)可得到相應(yīng)的qP波相速度。郝重濤等[28]基于坐標(biāo)變換的方法給出了TI介質(zhì)qP、qSV和qSH波精確相速度表達(dá)式。姚振岸等[29]利用Bond變換導(dǎo)出了TI介質(zhì)彈性波速度解析式，但精確相速度表達(dá)式較復(fù)雜。本文基于弱各向異性近似，導(dǎo)出相應(yīng)的線性近似相速度公式。由于目前的TI介質(zhì)RTM主要基于qP波，故本文主要討論qP波相速度。

設(shè)TI介質(zhì)對(duì)稱軸與z軸的夾角為θ0,當(dāng)θ0=0°時(shí)為VTI介質(zhì)，當(dāng)θ0≠0°時(shí)則為TTI介質(zhì)。利用旋轉(zhuǎn)坐標(biāo)法，將TTI介質(zhì)的qP波相速度表示為

(3)

同理，可近似簡(jiǎn)化為

(4)

式中γ為qP波傳播方向與對(duì)稱軸的夾角，式(3)、式(4)中其他參數(shù)的含義同式(1)、式(2)。可以看到，除了θ換成γ外，式(3)、式(4)和式(1)、式(2)在形式上完全一樣。本文用θ表示地震波傳播方向與z軸的夾角，γ表示地震波傳播方向與對(duì)稱軸的夾角。很顯然，VTI介質(zhì)的對(duì)稱軸與z軸重合，即γ=θ； TTI介質(zhì)的對(duì)稱軸與z軸的夾角為θ0，即γ=θ+θ0。

郝重濤等[28]的研究結(jié)果表明，TI介質(zhì)中體波速度隨傳播方向變化的速度圖案取決于對(duì)稱軸方向，并與TI介質(zhì)的Thomsen參數(shù)相關(guān)。速度特征只依賴于傳播矢量與對(duì)稱軸的夾角，具有一定的對(duì)稱性、漸變性和重復(fù)性。在各向異性參數(shù)確定的條件下，式(3)和式(4)只與γ相關(guān)，這與文獻(xiàn)[28]的研究結(jié)果一致。通過合理的簡(jiǎn)化，將對(duì)稱軸的方位角設(shè)置為0，即對(duì)稱軸位于垂直面xoz內(nèi)，這樣只需考慮θ0、θ以及傳播方向的方位角φ等三個(gè)變量。

為了計(jì)算的方便，引入單位向量，其中s為TI介質(zhì)對(duì)稱軸方向的單位向量，r為地震波傳播方向的單位向量，i為x軸方向的單位向量，j為y軸方向的單位向量，k為z軸方向的單位向量。根據(jù)地震波傳播方向和坐標(biāo)軸的關(guān)系(圖1)，單位向量表示為

(5)

s和r的夾角為γ，通過s與r的矢量積及數(shù)量積得到

cosγ=sinθcosφsinθ0+cosθcosθ0

(6)

sin2γ= sin2θsin2φ+

(sinθcosφcosθ0-cosθsinθ0)2

(7)

將式(6)、式(7)代入式(3)和式(4)即可求出相速度。雖然計(jì)算涉及到θ、θ0和φ三個(gè)變量，但最終都?xì)w結(jié)到γ一個(gè)變量。

表1為TI介質(zhì)各向異性參數(shù)，根據(jù)對(duì)稱軸方向的變化，對(duì)比、分析TTI介質(zhì)的qP波線性近似相速度與精確相速度的逼近程度。

圖1 地震波傳播方向與坐標(biāo)軸的關(guān)系

VP0/(m·s-1)VS0/(m·s-1)εδ300015000.250.1

圖2為由式(3)、式(4)得到的相速度。由圖可見：①近似式(式(4))的均方差ΔVRMS較小，表明由式(4)能較好地逼近精確公式(式(3))；②在VTI介質(zhì)中，當(dāng)θ0=0°、φ=0°(圖2a)與θ0=0°、φ=90°(圖2b)時(shí)得到的曲線完全一樣，說明VTI介質(zhì)的速度特征與傳播方位角無關(guān)；③在TTI介質(zhì)中速度特征具有方位性(圖2c、圖2d)。

1.2 時(shí)間—波數(shù)域波動(dòng)方程的推導(dǎo)

偽譜法的基本思想是在波數(shù)域計(jì)算空間導(dǎo)數(shù)，采用有限差分法計(jì)算時(shí)間導(dǎo)數(shù)[24]。在三維(x,y,z)空間內(nèi)傳播的平面波有下列相角關(guān)系

(8)

式中：VP為相速度；kx、ky、kz為波數(shù)；ω為角頻率。

將式(8)代入式(1)，得到VTI介質(zhì)的qP-qSV波頻散關(guān)系。Fletcher等[15]利用旋轉(zhuǎn)波數(shù)導(dǎo)出了TTI介質(zhì)的頻散關(guān)系

圖2 由式(3)、式(4)得到的相速度

(9)

(10)

式中φ0為對(duì)稱軸的方位角。令VSz=0，式(9)變?yōu)門TI介質(zhì)的聲學(xué)近似qP波頻散關(guān)系[4]

(11)

令f=1，將式(6)、式(7)代入式(3)也可導(dǎo)出式(11)。對(duì)式(11)進(jìn)行傅里葉逆變換，得到時(shí)間—空間域的4階偏微分方程，方程的解較為復(fù)雜。通過引入輔助函數(shù)

p(kx,ky,kz,ω)

(12)

Fletcher等[15]將四階偏微分方程轉(zhuǎn)換為等價(jià)的2階耦合方程。式(11)兩邊乘以p(kx,ky,kz,ω)，化簡(jiǎn)后有

(13)

(14)

對(duì)式(13)、式(14)進(jìn)行傅里葉逆變換，得到聲學(xué)近似二階耦合波動(dòng)方程

(15)

(16)

式中p為壓力波場(chǎng)。式中的偏微分算子為

(17)

(18)

(19)

從式(4)出發(fā)，導(dǎo)出純qP波時(shí)間—波數(shù)域波動(dòng)方程。將式(6)和式(7)代入式(4)，等式兩邊同乘以時(shí)間—波數(shù)域的波場(chǎng)U(kx,ky,kz,t)；再利用式(8)，同時(shí)僅在頻域應(yīng)用傅里葉逆變換，便可得到TTI介質(zhì)純qP波時(shí)間—波數(shù)域波動(dòng)方程

(20)

式(20)等號(hào)左邊的二階時(shí)間偏導(dǎo)數(shù)可寫成有限差分的形式，等號(hào)右邊則用簡(jiǎn)單的波數(shù)計(jì)算替代復(fù)雜的空間導(dǎo)數(shù)計(jì)算。

波動(dòng)方程數(shù)值模擬的吸收邊界條件大致分為兩種：一類是傍軸近似的吸收邊界條件；另一類是阻尼層吸收邊界條件。當(dāng)前通用的吸收邊界條件是最佳匹配層(PML)吸收邊界條件，屬于第二種吸收邊界條件。PML吸收邊界條件是由一階波動(dòng)方程得到的，無論是分裂算法還是非分裂算法，推廣到二階系統(tǒng)時(shí)算法較為復(fù)雜，對(duì)內(nèi)存需求也較大，故本文采用Cerjan等[30]提出的阻尼吸收邊界條件。

在數(shù)值模擬計(jì)算時(shí)，首先必須考慮方程解的穩(wěn)定性條件，此處給出θ0=0°時(shí)(VTI介質(zhì))的穩(wěn)定性條件(詳細(xì)推導(dǎo)過程見附錄A)

(21)

式中： Δt為時(shí)間步長(zhǎng)； Δd為網(wǎng)格空間步長(zhǎng)；α=min[abs(ε),abs(δ)]。

2 數(shù)值計(jì)算

2.1 均勻模型

首先設(shè)計(jì)一個(gè)VTI均勻介質(zhì)模型，圖3為由式(15)、式(16)及式(20)模擬得到的VTI均勻介質(zhì)模型在139ms的波場(chǎng)快照。由圖可見：由聲學(xué)近似二階耦合方程(式(15)、式(16))、純qP波方程(式(20))得到的波場(chǎng)快照的波場(chǎng)傳播特征十分吻合，但前者存在明顯的qSV波干擾(圖3a上～圖3c上)，后者沒有qSV波干擾(圖3a下～圖3c下)；由于VTI介質(zhì)對(duì)稱軸垂直，在任意方向的垂直切片的地震波場(chǎng)特征一致，故垂直切片的波前形態(tài)一致(圖3b、圖3c)，而水平切片反映了各向同性層，波前呈圓形(圖3d)。

圖4為由式(15)、式(16)及式(20)得到的TTI均勻介質(zhì)模型在139ms的波場(chǎng)快照。由圖可見：由純qP波方程(式(20))得到的波場(chǎng)快照只有純qP傳播，沒有偽橫波噪聲干擾(圖4a下～圖4d下)；由于TTI介質(zhì)對(duì)稱軸在xoz面內(nèi)傾斜，故相對(duì)于VTI介質(zhì)xoz垂直切片的波場(chǎng)發(fā)生旋轉(zhuǎn)(圖4b)；由于對(duì)稱軸位于yoz面外，在yoz面內(nèi)地震波傳播方向與對(duì)稱軸方向的夾角缺少0°～45°范圍內(nèi)的值，故在xoz面(圖4b)與yoz面的垂直切片(圖4c)的波場(chǎng)特征不同；由于xoz面外的任意水平或垂直平面都與對(duì)稱軸相交，這些平面內(nèi)的地震波傳播方向與對(duì)稱軸方向的夾角缺少0°～β(小于對(duì)稱軸傾角的角度)范圍內(nèi)的值，故不同方位的垂直波場(chǎng)切片形態(tài)也不同， 因此xoy面內(nèi)的波場(chǎng)切片(圖4d)表征了TTI介質(zhì)的方位各向異性。

2.2 Hess VTI模型

圖5為Hess VTI模型?？紤]到計(jì)算穩(wěn)定性條件，選取網(wǎng)格間距dx=dz=6.096m，時(shí)間采樣間隔為0.5ms，采樣時(shí)長(zhǎng)為1.5s，震源坐標(biāo)為(x,z)=(11009m,3645m)，震源函數(shù)采用雷克子波，主頻為20Hz。圖6為Hess VTI模型在1.5s的波場(chǎng)快照。由圖可見，由式(15)、式(16)(圖6a)和式(20)(圖6b)得到的波場(chǎng)快照的波場(chǎng)特征基本一致，說明聲學(xué)近模型網(wǎng)格數(shù)為201×201×201，網(wǎng)格間距為dx=dy=dz=5m，各向異性參數(shù)：VPz=3000m/s；ε=0.25；δ=0.1，VTI介質(zhì)對(duì)稱軸傾斜角θ0=0°。震源函數(shù)采用雷克子波，主頻為100Hz，震源坐標(biāo)為(x,y,z)=(500m,500m,500m)。數(shù)值計(jì)算的時(shí)間步長(zhǎng)為0.1ms，總采樣時(shí)長(zhǎng)為150ms似方程與純qP波方程的等價(jià)性，但前者存在偽橫波，在震源附近出現(xiàn)頻散，后者沒有偽橫波干擾。

圖3 由式(15)、式(16)(上)及式(20)(下)得到的VTI均勻介質(zhì)模型在139ms的波場(chǎng)快照

圖4 由式(15)、式(16)(上)及式(20)(下)得到的TTI均勻介質(zhì)模型在139ms的波場(chǎng)快照

2.3 BP 2007 TTI模型

圖7為BP 2007 TTI模型。由圖可見，存在強(qiáng)各向異性地層，局部區(qū)域發(fā)生對(duì)稱軸傾斜，呈TTI介質(zhì)特征。由于原始模型較大，本文選擇部分典型區(qū)域進(jìn)行正演計(jì)算。圖8為BP 2007 TTI模型在1.5s的波場(chǎng)快照。由圖可見，由聲學(xué)近似方程(圖8a)和純qP波方程(圖8b)得到的波場(chǎng)快照的波場(chǎng)特征基本一致，但后者的波場(chǎng)穩(wěn)定，未出現(xiàn)頻散現(xiàn)象，說明本文方法對(duì)復(fù)雜TTI模型具有較好的適應(yīng)性。

圖5 Hess VTI模型

圖6 Hess VTI模型在1.5s的波場(chǎng)快照

圖7 BP 2007 TTI模型

圖8 BP 2007 TTI模型在1.5s的波場(chǎng)快照

偽譜法涉及到傅里葉變換，計(jì)算效率低于有限差分法，但在同樣精度條件下，偽譜法可選取更大的空間網(wǎng)格間距和時(shí)間采樣間隔[24]。進(jìn)行大區(qū)域或三維模型數(shù)值模擬時(shí)，偽譜法的計(jì)算時(shí)間較長(zhǎng)。如果計(jì)算機(jī)硬件條件有限，可以采用混合法計(jì)算，即用有限差分法計(jì)算簡(jiǎn)單的空間偏導(dǎo)數(shù)項(xiàng)，用偽譜法計(jì)算復(fù)雜的空間偏導(dǎo)數(shù)項(xiàng)，減少了傅里葉變換次數(shù)，提高了計(jì)算效率，可在一定程度上減少計(jì)算耗時(shí)。

3 結(jié)束語

本文利用坐標(biāo)變換的方法，把Tsvankin[26]的VTI介質(zhì)qP波精確相速度擴(kuò)展到三維空間，基于弱各向異性近似推導(dǎo)了三維線性近似qP波相速度表達(dá)式，構(gòu)建了時(shí)間—波數(shù)域的純qP波波動(dòng)方程，并給出了穩(wěn)定性條件。通過對(duì)比、分析數(shù)值計(jì)算結(jié)果，得到以下認(rèn)識(shí)。

(1)在各向異性參數(shù)確定的條件下，任意對(duì)稱軸取向的TI介質(zhì)中qP波的傳播速度取決于傳播方向與TI對(duì)稱軸的夾角。

(2)把TTI介質(zhì)對(duì)稱軸的方位角取為0°，簡(jiǎn)化了計(jì)算，文中的各向異性模型為2.5維，但不同方向的切片仍然可以表征TTI介質(zhì)的方位各向異性。

(3)本文的純qP波波動(dòng)方程可以較好地應(yīng)用于TTI介質(zhì)正演模擬，得到的波場(chǎng)快照沒有偽橫波干擾，對(duì)復(fù)雜TTI模型仍然保持較好的波場(chǎng)穩(wěn)定性。偽譜法計(jì)算量較大，可以和有限差分法聯(lián)合使用以減少計(jì)算時(shí)間。

感謝Hess公司提供了VTI模型，感謝BP公司提供了BP 2007 TTI模型。

附錄A TTI介質(zhì)純qP波時(shí)間—波數(shù)域波場(chǎng)遞推格式

UP(kx,ky,kz,t+Δt)-2UP(kx,ky,kz,t)+

(A-1)

把試驗(yàn)解ei(kxx+kyy+kzz-ω t)代入式(A-1)，得

(A-2)

(A-3)

式(A-3)形式仍然很復(fù)雜，將式(A-3)等號(hào)兩邊開平方，并且取θ0=0°，得到VTI介質(zhì)的穩(wěn)定性條件

2>ΔtVPz×

(A-4)

令α=min[abs(ε),abs(δ)]，則式(A-4)簡(jiǎn)化為

(A-5)

根據(jù)采樣定理，信號(hào)帶寬應(yīng)小于奈奎斯特頻率(采樣頻率的二分之一)，即有

式中Δd為網(wǎng)格尺寸(空間步長(zhǎng))。整理上式，得

(A-6)