TTI介質qP波偽譜法正演模擬

張慶朝 朱國維 周俊杰 朱聰聰 劉衛剛

(中國礦業大學(北京)煤炭資源與安全開采國家重點實驗室，北京 100083)

0 引言

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

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

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

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

本文從Tsvankin[26]的VTI介質qP波精確相速度表達式出發，引入單位向量，利用坐標變換的方法把qP波的精確相速度擴展到三維TTI介質。基于弱各向異性近似和平方根近似，推導了三維TTI介質純qP波近似相速度和頻散關系。根據頻散關系構建了純qP波時間—波數域的波場遞推格式，并進行了相應的穩定性分析。最后利用偽譜法對不同的TI模型進行了qP波數值模擬。

1 方法原理

1.1 TTI介質的相速度

對于VTI介質，Tsvankin[26]給出了qP波精確相速度表達式

(1)

當介質的ε和δ遠小于1時，稱為弱各向異性介質，實際介質大多為弱各向異性介質。因為式(1)較復雜，將式中的根式進行泰勒展開，根據弱各向異性假設，舍去ε和δ的二次及以上高次項，得到線性近似的相速度表達式[26]

(2)

式(2)為純qP波相速度表達式，無S波項。

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

設TI介質對稱軸與z軸的夾角為θ0,當θ0=0°時為VTI介質，當θ0≠0°時則為TTI介質。利用旋轉坐標法，將TTI介質的qP波相速度表示為

(3)

同理，可近似簡化為

(4)

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

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

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

(5)

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

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)即可求出相速度。雖然計算涉及到θ、θ0和φ三個變量，但最終都歸結到γ一個變量。

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

圖1 地震波傳播方向與坐標軸的關系

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

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

1.2 時間—波數域波動方程的推導

偽譜法的基本思想是在波數域計算空間導數，采用有限差分法計算時間導數[24]。在三維(x,y,z)空間內傳播的平面波有下列相角關系

(8)

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

將式(8)代入式(1)，得到VTI介質的qP-qSV波頻散關系。Fletcher等[15]利用旋轉波數導出了TTI介質的頻散關系

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

(9)

(10)

式中φ0為對稱軸的方位角。令VSz=0，式(9)變為TTI介質的聲學近似qP波頻散關系[4]

(11)

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

p(kx,ky,kz,ω)

(12)

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

(13)

(14)

對式(13)、式(14)進行傅里葉逆變換，得到聲學近似二階耦合波動方程

(15)

(16)

式中p為壓力波場。式中的偏微分算子為

(17)

(18)

(19)

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

(20)

式(20)等號左邊的二階時間偏導數可寫成有限差分的形式，等號右邊則用簡單的波數計算替代復雜的空間導數計算。

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

在數值模擬計算時，首先必須考慮方程解的穩定性條件，此處給出θ0=0°時(VTI介質)的穩定性條件(詳細推導過程見附錄A)

(21)

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

2 數值計算

2.1 均勻模型

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

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

2.2 Hess VTI模型

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

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

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

2.3 BP 2007 TTI模型

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

圖5 Hess VTI模型

圖6 Hess VTI模型在1.5s的波場快照

圖7 BP 2007 TTI模型

圖8 BP 2007 TTI模型在1.5s的波場快照

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

3 結束語

本文利用坐標變換的方法，把Tsvankin[26]的VTI介質qP波精確相速度擴展到三維空間，基于弱各向異性近似推導了三維線性近似qP波相速度表達式，構建了時間—波數域的純qP波波動方程，并給出了穩定性條件。通過對比、分析數值計算結果，得到以下認識。

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

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

(3)本文的純qP波波動方程可以較好地應用于TTI介質正演模擬，得到的波場快照沒有偽橫波干擾，對復雜TTI模型仍然保持較好的波場穩定性。偽譜法計算量較大，可以和有限差分法聯合使用以減少計算時間。

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

附錄A TTI介質純qP波時間—波數域波場遞推格式

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

(A-1)

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

(A-2)

(A-3)

式(A-3)形式仍然很復雜，將式(A-3)等號兩邊開平方，并且取θ0=0°，得到VTI介質的穩定性條件

2>ΔtVPz×

(A-4)

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

(A-5)

根據采樣定理，信號帶寬應小于奈奎斯特頻率(采樣頻率的二分之一)，即有

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

(A-6)