梁鍇， 孫上饒， 曹丹平， 印興耀

中國石油大學(華東)地球科學與技術學院， 青島 266580

0 引言

TTI介質彈性波頻散關系方程可以通過求解Christoffel方程得到，利用不同的近似方法對其進行近似處理，并利用傅里葉逆變換將頻率-波數(shù)域算子變換到時空域，可以得到對應的解耦波動方程.Alkhalifah和Tsvankin(1995)、Alkhalifah(1998, 2000)開創(chuàng)性地提出了聲學假設近似方法，即將沿對稱軸的橫波速度近似視做零值，推導了一種簡化的VTI介質彈性波頻散關系方程，進而得到qP波方程.多種純qP方程相繼被提出(Klíe and Toro, 2001; Zhou et al., 2006; Hestholm, 2007; Du et al., 2008, 2010; Duveneck et al., 2008; Chu et al., 2011, 2013; Bloot et al., 2013; Schleicher and Costa, 2016; Li and Zhu, 2018; Xu et al., 2020)，但是聲學假設條件下的波動方程存在兩組共軛解，分別對應qP波和退化qSV波，而后者是不被期望的“干擾”波，在地震波數(shù)值模擬的過程中，隨著迭代次數(shù)的增加，退化qSV波會產(chǎn)生嚴重的數(shù)值頻散問題，特別在各向異性參數(shù)ε<δ的情況下，描述退化qSV波的解發(fā)散，誤差隨時間累計，導致數(shù)值解不穩(wěn)定.Jin和Stovas(2018)定義一組新的參數(shù)用于描述VTI介質中退化qSV波運動學特征，并分析了其傳播特征.為了能完全消除聲學近似中的qSV波，Klíe和Toro(2001)在此基礎上進行改進，消除了描述退化qSV波的一組解析解，得到純qP波波動方程.此外也可以通過引入輔助變量，實現(xiàn)對VTI介質qP波方程的降階處理，使其更容易實現(xiàn)(Zhou et al., 2006; Du et al., 2008; Chu et al., 2011).Liu等(2009)基于VTI介質聲學近似方程，解耦得到了一種穩(wěn)定的純qP波和退化qSV波的波動方程，并應用于逆時偏移中；Fletcher等(2009)基于坐標旋轉法推導了一般TTI介質qP波、qSV波的波動方程，并指出若人為設置TTI介質對稱軸方向的橫波速度足夠大，則可以移除qSV的波面三角區(qū)，進而使波場能穩(wěn)定傳播；Chu等(2011)對Fletcher推導的TTI介質擬聲波方程進行因式分解，利用一階和高階Taylor展開近似qP波方程，并指出采用一階Taylor展開近似的qP波方程足以準確運用于大部分實際應用中；梁鍇等(2009)從Thomsen弱各向異性近似和聲學假設近似出發(fā)，對TTI介質彈性波波動方程進行了分解，導出TTI介質弱各向異性條件近似qP波和qSV波波動方程以及聲學假設近似的qP波波動方程；黃翼堅等(2011)對VTI介質彈性波精確相速度表達式進行多項式近似，得到了VTI介質純qP波方程；楊鵬等(2017)利用近似展開法，通過分解偏微分算子得到了TI介質純qP波動方程；張慶朝等(2019)根據(jù)弱各向異性近似的假設，推導了任意空間取向TI介質彈性波相速度近似表達式，并進一步導出了qP波和qSV波的解耦波動方程，波場模擬結果顯示該波動方程有著較好的穩(wěn)定性；Chu等(2013)利用偽譜法對TTI介質擬聲波方程(Chu et al., 2011)進行了數(shù)值模擬，誤差分析表明一階Taylor近似適合較弱的各向異性介質，強各向異性介質中需要采用更高階的近似，以控制誤差在允許范圍內(nèi)；杜啟振等(2015)采用偽譜法和有限差分法相結合的混合法求解VTI介質純qP方程，極大地提升了運算效率.

本文在TTI介質彈性波精確頻散關系的基礎上，利用近似的配方法(梁鍇等，2018；孫上饒等，2021)，推導TTI介質qP波和qSV波近似頻散關系，并將頻率-波數(shù)域算子反變換到時-空域，導出了qP波和qSV波近似解耦的波動方程.使用兩組模型參數(shù)對近似頻散關系方程進行數(shù)值計算，得到了三維頻散關系曲面和二維曲線，驗證了近似頻散關系的有效性.隨后考察了近似式與精確式的差異項數(shù)值在不同各向異性強度下的變化情況，并分析了XOZ面內(nèi)近似頻散關系曲線的相對誤差分布.最后，使用有限差分方法求解TTI介質純qP波和純qSV波近似波動方程，模擬了qP波和qSV波在均勻、層狀及復雜TTI介質中的傳播.數(shù)值模擬結果表明，在η<0和各向異性傾角變化較大的介質中，純qP波和純qSV波近似波動方程依然可以保持穩(wěn)定.

1 方法原理

1.1 TTI介質彈性波近似頻散關系方程

TTI介質彈性波相速度和頻散關系可通過將平面波位移方程代入到Christoffel方程求解得到.TTI介質qP波和qSV波精確頻散關系方程可表示為：

(1)

圖1 TTI介質模型示意圖Fig.1 Schematic diagram of TTI media model

依照近似配方法的思想(梁鍇等, 2018)消除D項的一次根號，簡化頻散關系方程(1)，其中式(1)中D項可改寫為：

(2)

(3)

(4)

令式(4)中右邊第二項近似為零，即：

(5)

那么，D項可近似表示為D≈(a+b)2=Da.將D的近似值Da代入到精確qP波和qSV波頻散關系(1)中，即可得到相對應的基于近似配方法的頻散關系方程：

(6)

(7)

1.2 近似純qP波和純qSV波波動方程

(8)

(9)

解耦波動方程(8)和(9)既包含對時間和空間的混合偏導，也存在混合空間導數(shù)，這使得其計算量大于單一變量的空間導數(shù)(Fletcher et al., 2009)，不利于波動方程的數(shù)值求解.考慮波動方程(8)的二維形式，并引入輔助變量QP，變換得到等效的純qP波波動方程為：

(10)

(11)

同理引入輔助變量QSV，對方程(9)進行變換，得到等效qSV波波動方程為：

(12)

(13)

可以利用方程(10)—(13)分別進行純qP波和純qSV波的數(shù)值模擬.

2 數(shù)值示例

2.1 數(shù)值試算和誤差分析

Vernik和Liu(1997)對多組含油氣地層進行了地震波相速度和各向異性參數(shù)的測量，本文從測量數(shù)據(jù)中選取兩組中強各向異性頁巖作為TTI介質模型，參數(shù)如表1所示.

表1 TTI介質模型參數(shù)Table 1 Parameters of TTI models

通過對比qP波和qSV波的精確和近似頻散關系中的D項，定義差異項ΔD為：

(14)

利用表1的model 1的速度參數(shù)和傾角θ0對ΔD進行數(shù)值計算，在|σ|值固定不變的情況下，僅改變參數(shù)η和δ的值，結果如圖2所示，圖中空心箭頭處為沿對稱軸方向傳播，實心箭頭處為垂直對稱軸方向傳播.參數(shù)|η|越小，差異ΔD越小，并且在對稱軸傳播方向附近的誤差相對垂直對稱軸傳播方向小；在沿對稱軸和垂直對稱軸的傳播方向，近似值與精確值相等.特別地，當η=0時，介質呈現(xiàn)橢圓各向異性，近似頻散關系方程(5)等于精確式(1).

圖2 ΔD隨相角θ的變化特征(a) η>0； (b) η<0.Fig.2 Variation characteristics of the difference ΔD with phase angle θ

對TTI介質qP波和qSV波頻散關系進行三維數(shù)值計算，結果如圖3所示.圖3a、b為model 1的精確qP波和qSV波頻散關系曲面，圖3e、f為model 1的近似qP波和qSV波頻散關系曲面，圖3c、d為model 2的精確qP波和qSV波頻散關系曲面，圖3g、h為model 2的近似qP波和qSV波頻散關系曲面.兩組模型中的近似qP波和qSV波頻散關系曲面(圖3e—h)與精確頻散關系曲面(圖3a—d)方位特征較為相似，其中qP波的精度明顯高于qSV波.

圖3 TTI介質彈性波頻散關系曲面(a) model 1 qP波精確值; (b) model 1 qSV波精確值; (c) model 2 qP波精確值; (d) model 2 qSV波精確值; (e) model 1 qP波近似值; (f) model 1 qSV波近似值; (g) model 2 qP波近似值; (h) model 2 qSV波近似值.Fig.3 Surfaces of elastic wave dispersion relation in TTI media(a) Exact value of qP wave in model 1； (b) Exact value of qSV wave in model 1； (c) Exact value of qP wave in model 2； (d) Exact value of qSV wave in model 2； (e) Approximate value of qP wave in model 1； (f) Approximate value of qSV wave in model 1； (g) Approximate value of qP wave in model 2； (h) Approximate value of qSV wave in model 2.

為考察近似值與精確值的誤差，這里討論頻散曲面在XOZ平面內(nèi)的頻散關系曲線的相對誤差分布情況.在觀測坐標系XOZ平面內(nèi)，基于近似配方法的頻散關系方程(6)可表示為：

(15)

定義相對誤差(RE)為：

(16)

根據(jù)表1中兩組模型參數(shù)可得到二維TTI介質彈性波頻散關系曲線，如圖4所示，圖中短實線表示TTI介質對稱軸方向，qP波和qSV波近似頻散關系曲線與精確值在對稱軸和垂直對稱軸方向誤差為零，且誤差關于上述兩個傳播方向呈現(xiàn)對稱分布的形式.對比圖4a、b可知，qP波近似頻散關系精度整體上高于qSV波，原因在于qSV波近似頻散關系的精度在較大程度上受組合參數(shù)σ的控制，該參數(shù)中縱橫波速度比起到放大ε和δ之間差值的作用，從而導致了誤差的增大.圖5較為直觀地呈現(xiàn)了qP波和qSV波近似頻散關系曲線相對誤差在θ∈[0,180°]范圍內(nèi)的分布情況，本例中沿對稱軸(空心箭頭處)和垂直對稱軸(實心箭頭處)的兩個傳播方向對應的相對誤差為零，并且相對誤差在這兩個傳播方向附近一定角度范圍內(nèi)很小，這同時也契合差異項ΔD的變化趨勢(圖2).

圖4 TTI介質彈性波頻散關系曲線(a) model 1 qP波； (b) model 1 qSV波; (c) model 2 qP波； (d) model 2 qSV波; 灰色實線表示精確頻散關系曲線，黑色虛線表示近似頻散關系曲線.Fig.4 Dispersion curves of elastic waves in TTI media(a) qP wave in model 1； (b) qSV wave in model 1； (c) qP wave in model 2； (d) qSV wave in model 2. Gray solid curve represents exact values, and black represents approximate values.

圖5 qP波和qSV波近似頻散關系曲線的相對誤差(a) model 1相對誤差； (b) model 2 相對誤差.其中灰色實線表示qP波頻散關系相對誤差EP，黑色虛線表示qSV波頻散關系相對誤差ESV.Fig.5 Relative error of approximate dispersion relation curve for qP and qSV wave(a) Model 1； (b) Model 2. Gray solid curve represents relative error of dispersion relationship for qP wave (EP), and the black dashed curve represents relative error of dispersion relationship for qSV wave (ESV).

2.2 波場模擬

為驗證文中導出的TTI介質qP波和qSV波近似波動方程(10)—(13)的有效性，本文利用有限差分法對其進行均勻、層狀和復雜介質波場數(shù)值模擬.首先，設置均勻介質模型大小為401×401，空間間隔為Δx=Δz=10 m，時間采樣間隔為Δt=1 ms.模型參數(shù)為：vP0=3.00 km·s-1，vS0=1.73 km·s-1，η=0.143，δ=0.2，σ=0.601，θ0=π/4；震源子波是主頻為25 Hz的雷克子波，位于模型中心點.高階空間有限差分近似有利于壓制數(shù)值頻散，但會增加運算成本，合理選擇差分近似的階數(shù)可以在適當增加運算時間的情況下得到更好的數(shù)值模擬結果.考慮到泊松方程的求解效率問題，本文對方程(10)—(13)取空間十階、時間二階精度的差分近似進行波場模擬.圖6展示了TTI介質彈性波以及解耦的qP波和qSV波波場快照.對比圖6a、c可知純qP波的運動學特征基本與彈性波中qP波的特征相似，即使在大偏移距處，兩者的波場特征也很相似；對比圖6b、c可知，相比于彈性波qSV波場，純qSV波的波前面三角區(qū)不明顯，但其運動學特征整體保持較好的相似度.

圖6 均勻TTI介質波場快照(a) 純qP波； (b) 純qSV波; (c) 彈性波.Fig.6 Snapshots of homogeneous TTI media(a) Pure qP wave; (b) Pure qSV wave; (c) Elastic wave.

其次，對層狀各向異性介質模型進行數(shù)值模擬.模型大小為401×401，空間步長為Δx=Δz=10 m，時間采樣間隔為Δt=1 ms，模型參數(shù)如表2所示.注意到表2中l(wèi)ayer 1和layer 2為η>0的TTI介質模型，layer 3是η<0的TTI介質模型，這樣設置是為了測試純qP波和純qSV波波動方程的穩(wěn)定性.震源子波與均勻介質模型中相同，位于(2 km, 0.1 km)處，檢波點設置于深度為0.1 km的水平面上，如圖7a所示；圖7b—d為0.8 s時刻的TTI介質地震波波場快照.通過對比可知，純qP波(圖7b)和純qSV波(圖7c)入射到界面時只產(chǎn)生反射、透射的同類波，不產(chǎn)生轉換波.相對常規(guī)彈性波場(圖7d)，前兩者有著更加清晰的波場，利于研究同類型波的傳播特征.此外，數(shù)值模擬結果表明，純qP波方程在η<0和各向異性傾角變化較大情況下，仍然可以保持較好的穩(wěn)定性.

圖7 層狀TTI介質模型(a)、純qP波波場快照(b)、純qSV波波場快照(c)以及彈性波波場快照(d)Fig.7 Layered TTI media model (a), snapshots of pure qP wave (b), pure qSV wave (c) and elastic wave (d)

表2 層狀各向異性介質模型Table 2 Layered anisotropic medium model

最后，本文對BP鹽丘模型進行數(shù)值模擬；模型大小為601×601，空間間隔和時間采樣間隔與上述均勻介質模型一致，模型參數(shù)如圖8所示.震源子波與均勻介質模型中相同，位于(3 km, 3 km)處.數(shù)值模擬結果表明，純qP波波場中只可見透射和反射qP波的信息，其波前面與彈性波qP波的波前面位置幾乎一致；相應地，純qSV波中只可見透射和反射qSV波的信息；從相速度誤差分析中可見，qSV波的誤差大于qP波的誤差，因此純qSV波的波前面與彈性波SV波波前面存在一定的差異，但總體上在可接受的范圍內(nèi)(圖9).

圖8 BP鹽丘模型(a) vP0; (b) vS0; (c) η; (d) δ; (e) θ0.Fig.8 Salt dome model of BP

圖9 BP鹽丘模型波場快照(a) 純qP波; (b) 純qSV波; (c) 彈性波.Fig.9 Snapshots of BP salt dome model(a) Pure qP wave; (b) Pure qSV wave; (c) Elastic wave.

3 結論和認識

本文在TTI介質彈性波精確頻散關系的基礎上，利用近似配方法推導了qP波和qSV波近似頻散關系，通過傅里葉逆變換將其從頻率-波數(shù)域變換到時-空域，進而導出了純qP波和純qSV波近似解耦波動方程.

(1) 理論分析和數(shù)值計算表明，近似和精確頻散關系的差異項ΔD與|η|呈正相關，在橢圓各向異性介質中差異項為零，即近似式與精確式等效.

(2) 近似頻散關系曲線相對誤差整體較小，且在沿對稱軸和垂直對稱軸的傳播方向上相對誤差為零.

(3) 參數(shù)σ中縱橫波速度比使得qSV波近似頻散關系曲線相對誤差普遍較qP波大.

(4) 數(shù)值結果表明，基于近似配方法解耦的純qP波和純qSV波波動方程在η<0和對稱軸傾角變化較大的介質模型中依然可以保持穩(wěn)定；相對于彈性波波場，兩者均不產(chǎn)生轉換波，波場更加清晰.

本文推導的TTI介質解耦波動方程可以被進一步推廣到任意空間取向的TI介質中，此時需要從三維空間的角度討論解耦縱橫波的傳播特征.本文的解耦波動方程可以用于各向異性介質正演模擬和逆時偏移算法中，但因用常規(guī)有限差分方法求解泊松方程的效率較低，可以考慮采用偽譜法、低秩分解法、有限差分法和偽譜法的混合法等效率較高的方法求解純qP波和純qSV波方程，以提高正演模擬和偏移成像的效率.

附錄A

在qP波和qSV波近似頻散關系(6)兩邊同時乘以波場PP(ω,kx,ky,kz)和PSV(ω,kx,ky,kz)，再進行傅里葉逆變換(ω→i?/?t,kx→i?/?x,ky→i?/?y,kz→i?/?z)，化簡合并后得到三維TTI介質純qP波和純qSV波的近似波動方程：

(A1)

(A2)

(A3)

(A4)