TTI介質(zhì)彈性波近似解耦波動(dòng)方程

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

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

0 引言

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

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

1 方法原理

1.1 TTI介質(zhì)彈性波近似頻散關(guān)系方程

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

(1)

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

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

(2)

(3)

(4)

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

(5)

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

(6)

(7)

1.2 近似純qP波和純qSV波波動(dòng)方程

(8)

(9)

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

(10)

(11)

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

(12)

(13)

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

2 數(shù)值示例

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

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

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

通過對(duì)比qP波和qSV波的精確和近似頻散關(guān)系中的D項(xiàng)，定義差異項(xiàng)ΔD為：

(14)

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

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

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

圖3 TTI介質(zhì)彈性波頻散關(guān)系曲面(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)的頻散關(guān)系曲線的相對(duì)誤差分布情況.在觀測(cè)坐標(biāo)系XOZ平面內(nèi)，基于近似配方法的頻散關(guān)系方程(6)可表示為：

(15)

定義相對(duì)誤差(RE)為：

(16)

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

圖4 TTI介質(zhì)彈性波頻散關(guān)系曲線(a) model 1 qP波； (b) model 1 qSV波; (c) model 2 qP波； (d) model 2 qSV波; 灰色實(shí)線表示精確頻散關(guān)系曲線，黑色虛線表示近似頻散關(guān)系曲線.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波近似頻散關(guān)系曲線的相對(duì)誤差(a) model 1相對(duì)誤差； (b) model 2 相對(duì)誤差.其中灰色實(shí)線表示qP波頻散關(guān)系相對(duì)誤差EP，黑色虛線表示qSV波頻散關(guān)系相對(duì)誤差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 波場(chǎng)模擬

為驗(yàn)證文中導(dǎo)出的TTI介質(zhì)qP波和qSV波近似波動(dòng)方程(10)—(13)的有效性，本文利用有限差分法對(duì)其進(jìn)行均勻、層狀和復(fù)雜介質(zhì)波場(chǎng)數(shù)值模擬.首先，設(shè)置均勻介質(zhì)模型大小為401×401，空間間隔為Δx=Δz=10 m，時(shí)間采樣間隔為Δ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的雷克子波，位于模型中心點(diǎn).高階空間有限差分近似有利于壓制數(shù)值頻散，但會(huì)增加運(yùn)算成本，合理選擇差分近似的階數(shù)可以在適當(dāng)增加運(yùn)算時(shí)間的情況下得到更好的數(shù)值模擬結(jié)果.考慮到泊松方程的求解效率問題，本文對(duì)方程(10)—(13)取空間十階、時(shí)間二階精度的差分近似進(jìn)行波場(chǎng)模擬.圖6展示了TTI介質(zhì)彈性波以及解耦的qP波和qSV波波場(chǎng)快照.對(duì)比圖6a、c可知純qP波的運(yùn)動(dòng)學(xué)特征基本與彈性波中qP波的特征相似，即使在大偏移距處，兩者的波場(chǎng)特征也很相似；對(duì)比圖6b、c可知，相比于彈性波qSV波場(chǎng)，純qSV波的波前面三角區(qū)不明顯，但其運(yùn)動(dòng)學(xué)特征整體保持較好的相似度.

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

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

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

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

最后，本文對(duì)BP鹽丘模型進(jìn)行數(shù)值模擬；模型大小為601×601，空間間隔和時(shí)間采樣間隔與上述均勻介質(zhì)模型一致，模型參數(shù)如圖8所示.震源子波與均勻介質(zhì)模型中相同，位于(3 km, 3 km)處.數(shù)值模擬結(jié)果表明，純qP波波場(chǎng)中只可見透射和反射qP波的信息，其波前面與彈性波qP波的波前面位置幾乎一致；相應(yīng)地，純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鹽丘模型波場(chǎng)快照(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 結(jié)論和認(rèn)識(shí)

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

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

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

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

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

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

附錄A

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

(A1)

(A2)

(A3)

(A4)