TTI介質(zhì)純準(zhǔn)P波一階壓力-速度方程及求解方法

李佳珂，張會星，白 冰，張建敏

(1.中國海洋大學(xué) 山東 青島 266100； 2.海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東 青島 266100)

0 引言

地球是一個(gè)非均勻、廣泛存在各向異性的介質(zhì)體。其中具有傾斜對稱軸的橫向各向同性(Titled Transverse Isotropic，TTI)介質(zhì)作為描述地層中各向異性介質(zhì)最具一般代表性的模型[1]，受到了地球物理勘探界的廣泛關(guān)注。研究TTI介質(zhì)中地震波的傳播規(guī)律對于提高裂縫性地層的地震勘探精度具有重要意義[1-10]。

前人已對于聲學(xué)各向異性理論進(jìn)行了大量的研究。Alkhalifah[11-12]首先提出了著名的聲學(xué)假設(shè)：即將沿對稱軸方向的剪切波速度設(shè)置為零，并基于該理論推導(dǎo)了具有垂直對稱軸的橫向各向同性(Vertical Transverse Isotropic，VTI)介質(zhì)四階偽聲波波動(dòng)方程；吳國忱等[13]在頻率—空間域試算了該四階偽聲波方程的數(shù)值解；Du等[14]對方程進(jìn)行化簡，引入輔助函數(shù)，推導(dǎo)出VTI介質(zhì)二階偽聲波波動(dòng)方程；Duveneck[15]通過變換彈性張量中的系數(shù)，使得輔助波場具有實(shí)際的物理意義，也獲得了相應(yīng)的VTI耦合準(zhǔn)P波方程；Hestholm等[16]通過引入輔助參數(shù)，推導(dǎo)出了VTI介質(zhì)下耦合準(zhǔn)P波方程的一階形式；韓令賀和何兵壽[17]針對該方程，利用高階交錯(cuò)網(wǎng)格有限差分進(jìn)行了求解；Zhou[5]、Fletcher[6-7]將VTI介質(zhì)二階方程推廣至TTI介質(zhì)；Flowler等[18]研究發(fā)現(xiàn)該類耦合波動(dòng)方程都是等價(jià)的，并給出了TI介質(zhì)二階耦合準(zhǔn)P波方程的一般形式。然而，實(shí)際模型試算結(jié)果表明，該類TTI介質(zhì)耦合波動(dòng)方程存在菱形偽橫波干擾噪聲，對各向異性介質(zhì)中準(zhǔn)P波傳播規(guī)律產(chǎn)生影響；在Thomsen參數(shù)δ>ε時(shí)，方程難以穩(wěn)定求解；此外當(dāng)TTI介質(zhì)空間各點(diǎn)的對稱軸傾角陡然驟變時(shí)，方程會產(chǎn)生數(shù)值不穩(wěn)定的現(xiàn)象[8]。

近年來，針對消除聲波方程中的偽SV波，國內(nèi)外學(xué)者們給出了幾種不同的各向異性介質(zhì)準(zhǔn)P波方程。Xu等[19]通過把擬微分算子分解為標(biāo)量算子和微分算子，實(shí)現(xiàn)對各向異性介質(zhì)聲波方程相速度的控制，由此推導(dǎo)出VTI介質(zhì)的二階純準(zhǔn)P波控制方程；Zhan等[9]提出旋轉(zhuǎn)波數(shù)的概念，科學(xué)地將VTI介質(zhì)解耦波動(dòng)方程旋轉(zhuǎn)到了TTI介質(zhì)中，并且推演出TTI介質(zhì)下的純P波方程形式，利用偽譜法和有限差分法聯(lián)合求解了該波動(dòng)方程；Zhang等[10]在Xu等[19]的基礎(chǔ)上，通過提出新的旋轉(zhuǎn)方式推導(dǎo)出TTI介質(zhì)三維二階純準(zhǔn)P波方程，該方程相比較Zhan等[9]的方程，不存在偽微分算子，且結(jié)構(gòu)簡單，便于利用有限差分法快速求解。

本文以Zhang等[10]的研究成果為基礎(chǔ)，通過引入輔助波場推導(dǎo)出TTI介質(zhì)純準(zhǔn)P波一階壓力-速度方程；利用旋轉(zhuǎn)交錯(cuò)網(wǎng)格技術(shù)實(shí)現(xiàn)了數(shù)值模擬；并且引入了完全匹配層(Perfectly Matched Layers，PML)吸收邊界條件處理人工邊界問題。

1 TTI介質(zhì)二階純準(zhǔn)P波方程的建立

Xu[20]根據(jù)Alkhalifah[11]所提出的聲學(xué)假設(shè)下VTI介質(zhì)頻散關(guān)系，通過引入橢圓微分輔助算子R，推導(dǎo)并改進(jìn)了VTI介質(zhì)中的純準(zhǔn)P波方程：

式中橢圓微分輔助算子

二階偏導(dǎo)數(shù)

TTI介質(zhì)可以認(rèn)為是VTI介質(zhì)的一般形式，可以通過投影變換將VTI介質(zhì)轉(zhuǎn)換為TTI介質(zhì)。故要得到TTI介質(zhì)下的純準(zhǔn)P波方程，可將傳播方向從原始坐標(biāo)系(XYZ)進(jìn)行兩次旋轉(zhuǎn)變換投影到TTI介質(zhì)本構(gòu)坐標(biāo)系(X′Y′Z′)中，設(shè)第一次旋轉(zhuǎn)產(chǎn)生的過渡坐標(biāo)系為(X0Y0Z0)。圖1給出了兩次坐標(biāo)系旋轉(zhuǎn)之間相互關(guān)系。

圖1 坐標(biāo)系空間關(guān)系示意圖Figure 1 A schematic diagram of coordinate system spatial relation

(2)

第二次坐標(biāo)系旋轉(zhuǎn)，將(X0Y0Z0)坐標(biāo)系以Z0軸為對稱軸逆時(shí)針旋轉(zhuǎn)φ角度，得到(X′Y′Z′)坐標(biāo)系。對應(yīng)的旋轉(zhuǎn)矩陣為：

(3)

將兩次旋轉(zhuǎn)過程合并為一個(gè)旋轉(zhuǎn)關(guān)系：

(4)

進(jìn)而得到波場傳播方向向量的坐標(biāo)轉(zhuǎn)換關(guān)系：

(5)

接下來求取TTI介質(zhì)本構(gòu)坐標(biāo)系(X′Y′Z′)下波場傳播的方向向量表達(dá)式。根據(jù)式(1)可以很容易得到該坐標(biāo)系下TTI介質(zhì)純準(zhǔn)P波方程形式：

(6)

(7)

將式(5)代入式(7)并對其作傅里葉逆變換至?xí)r間-空間域，即可得到TTI介質(zhì)純準(zhǔn)P波二階方程：

+C5uxy+C6uyz)

(8)

式中R的求取同式(6)；系數(shù)Ci如下：

2 TTI介質(zhì)純準(zhǔn)P波一階壓力-速度方程的建立及求解方法

2.1 方程的建立

TTI介質(zhì)純準(zhǔn)P波一階壓力-速度方程方程的推導(dǎo)采用類似各向同性介質(zhì)中二階聲波方程降階的方法。首先引入輔助波場p、vx、vy、vz，設(shè)：

(9)

將式(9)帶入式(8)可以得到TTI介質(zhì)純準(zhǔn)P波一階壓力-速度方程組：

2.2 數(shù)值求解方法

對于本文建立的一階壓力-速度方程，可以采用旋轉(zhuǎn)交錯(cuò)網(wǎng)格[22]高階有限差分方法來求解。旋轉(zhuǎn)交錯(cuò)網(wǎng)格發(fā)展于交錯(cuò)網(wǎng)格[23-25]的基礎(chǔ)之上，它通過網(wǎng)格旋轉(zhuǎn)差分求解TTI介質(zhì)中的地震波方程，避免了常規(guī)網(wǎng)格差分求解中的插值和擬合問題[26]。旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分具體實(shí)現(xiàn)過程[22]是：首先計(jì)算出沿對角線方向的波場差分值，然后用這些對角線方向差分的線性組合求解出沿坐標(biāo)軸的差分?jǐn)?shù)值，其網(wǎng)格剖分方法如圖2所示。

圖2 三維旋轉(zhuǎn)交錯(cuò)網(wǎng)格示意圖Figure 2 A schematic diagram of 3D rotating staggered grid

(11)

(12)

式中en為差分系數(shù)，可通過Taylor展開求解[27]，其同階具體數(shù)值與交錯(cuò)網(wǎng)格[28]相同。經(jīng)過線性迭代組合便可得到沿坐標(biāo)軸方向的旋轉(zhuǎn)交錯(cuò)網(wǎng)格差分算子計(jì)算公式：

(13)

2.3 人工邊界條件

邊界條件是地震波動(dòng)方程求解的重要研究內(nèi)容[23-24,29-30]，本文應(yīng)用PML吸收邊界條件處理人工邊界問題。依據(jù)完全匹配層的分裂思路[29]，在式(10)中引入沿x、y、z方向的邊界衰減因子dx(x)、dy(y)、dz(z)，以p分量為例，給出了TTI介質(zhì)純準(zhǔn)P波一階壓力-速度方程的PML吸收邊界條件：

(14)

3 數(shù)值模擬

由于考慮到計(jì)算效率，本文將在二維情況下對上述純準(zhǔn)P波一階壓力-速度方程及其他學(xué)者給出的方程進(jìn)行數(shù)值模擬，并對數(shù)值模擬結(jié)果進(jìn)行對比分析。

3.1 均勻TTI介質(zhì)模型

給定一個(gè)單層均勻的TTI介質(zhì)模型，模型參數(shù)見表1；其網(wǎng)格尺寸4 000 m×4 000 m，空間步長Δx=Δz=10 m，時(shí)間采樣步長Δt=1 ms；激發(fā)震源選用Ricker子波置于介質(zhì)中心，主頻16Hz。對上述模型分別采用Zhang提出的純準(zhǔn)P波二階方程[8]及本文提出的純準(zhǔn)P波一階壓力-速度方程進(jìn)行數(shù)值模擬，獲得0.6s的波前快照對比如圖3所示。

表1 均勻TTI介質(zhì)模型參數(shù)

由圖3可知兩種方程的波前快照傳播規(guī)律相同，但純準(zhǔn)P波一階壓力-速度方程的計(jì)算精度略高于純準(zhǔn)P波二階方程；此外，PML邊界在兩種方程下都能夠很好地吸收能量，解決人工邊界問題。在相同條件下，一階方程計(jì)算效率明顯高于二階方程：二階方程耗時(shí)344s，而一階方程只耗時(shí)190s。

3.2 水平層狀模型

給定一個(gè)兩層水平層狀模型，其中第一層為各向同性介質(zhì)，第二層為TTI介質(zhì)。模型參數(shù)見表2；模型尺寸規(guī)模為2 000 m×2 000 m，空間步長Δx=Δz=5 m，時(shí)間采樣步長Δt=0.2ms；震源同樣選用Ricker子波，主頻24.0Hz，并置于(1 000m，0m)處；接收點(diǎn)位于地表，界面位于模型中間，記錄時(shí)間層數(shù)為4 000，400道接收。對上述模型分別采用Zhou提出的二階耦合準(zhǔn)P波方程[5]及本文提出的純準(zhǔn)P波一階壓力-速度方程進(jìn)行數(shù)值模擬，對比如圖4所示。

表2 水平層狀介質(zhì)模型參數(shù)

由圖4可知，本文給出的純準(zhǔn)P波一階壓力-速度方程沒有偽橫波噪音干擾，而二階耦合準(zhǔn)P方程在模擬TTI介質(zhì)中地震波傳播的過程中存在圖4(a)所示的低頻低速的偽橫波。這種不可消除的波型在逆時(shí)偏移時(shí)會產(chǎn)生成像假象，從而影響我們對地層結(jié)構(gòu)的判斷。

(a)純準(zhǔn)P波二階方程；(b)純準(zhǔn)P波一階壓力-速度方程圖3 均勻TTI介質(zhì)波前快照對比Figure 3 Homogeneous TTI medium wavefront snapshots comparison

(a)二階耦合準(zhǔn)P波方程波前快照；(b)純準(zhǔn)P波一階壓力-速度方程波前快照；(c)純準(zhǔn)P波一階壓力-速度方程單炮地震記錄圖4 水平層狀介質(zhì)波前快照和單炮地震記錄Figure 4 Horizontal layered medium wavefront snapshot and single-shot seismic record

3.3 TTI介質(zhì)Wedge模型試算

通過均勻TTI介質(zhì)模型試算，可知本文推導(dǎo)的純準(zhǔn)P波一階壓力-速度方程與二階形式方程的數(shù)值模擬結(jié)果在計(jì)算精度上并無明顯區(qū)別，但計(jì)算效率更高。但當(dāng)模型對稱軸方向變化較為劇烈時(shí)，是否會產(chǎn)生不穩(wěn)定問題，并不能獲得很好的驗(yàn)證。

為說明該問題，本文給出了一個(gè)TTI介質(zhì)楔形模型。模型尺寸為3 000 m×3 000 m；空間步長Δx=Δz=10 m，時(shí)間采樣步長Δt=0.5 ms空間采樣間隔為10m；震源為Ricker子波，主頻15Hz，埋深10m；檢波器位于地表300道接收；模型介質(zhì)參數(shù)如圖5、圖6所示。對上述模型分別采用Zhou提出的二階耦合準(zhǔn)P波方程[5]及本文提出的純準(zhǔn)P波一階壓力-速度方程進(jìn)行數(shù)值模擬。

純準(zhǔn)P波一階壓力-速度方程數(shù)值模擬結(jié)果如圖7所示；二階耦合準(zhǔn)P波方程數(shù)值模擬結(jié)果如圖8所示。由波場快照結(jié)果對比可知：論文推導(dǎo)的一階純準(zhǔn)P波方程，在各向異性對稱軸角度變化幅度最大的地方，波場并沒有發(fā)生不穩(wěn)定的現(xiàn)象，仍然能夠在介質(zhì)中保持穩(wěn)定傳播；而二階耦合準(zhǔn)P波方程，在t>1.5s時(shí)，已經(jīng)出現(xiàn)了不穩(wěn)定現(xiàn)象。超過3.5s后，結(jié)果已經(jīng)基本無法顯示。

4 結(jié)論

①推導(dǎo)了TTI介質(zhì)純準(zhǔn)P波一階壓力-速度方程，并給出了數(shù)值求解的方法。相比純準(zhǔn)P波二階方程，本文給出的方程在計(jì)算精度相同的情況下，具有更高的計(jì)算效率；

②本文給出的一階方程可以更好地適用PML吸收邊界條件，同時(shí)能夠結(jié)合旋轉(zhuǎn)交錯(cuò)網(wǎng)格技術(shù)，避免了進(jìn)行各向異性介質(zhì)數(shù)值模擬時(shí)由于網(wǎng)格剖分產(chǎn)生的特殊誤差，提高了數(shù)值模擬結(jié)果的精度；

③設(shè)計(jì)了楔形TTI介質(zhì)模型， 驗(yàn)證了本文的方程能夠穩(wěn)定地模擬對稱軸參數(shù)變化劇烈的情況。由此可認(rèn)為該純準(zhǔn)P波方程對于改善各向異性逆時(shí)偏移成像的質(zhì)量以及提高計(jì)算效率有一定推動(dòng)作用。

(a)縱波速度Vp0；(b)極化角θ圖5 楔形模型速度和角度參數(shù)Figure 5 Wedge model velocity and angle parameters

(a)Thomsen參數(shù)δ；(b) Thomsen參數(shù)ε圖6 楔形模型各向異性參數(shù)Figure 6 Wedge model anisotropic parameters

(a)t=1.6s時(shí)刻波前快照；(b)單炮地震記錄圖7 純準(zhǔn)P波一階壓力-速度方程數(shù)值模擬結(jié)果Figure 7 Pure quasi-P wave first-order pressure-velocity equation numerical simulation results

(a)t=1.6s時(shí)刻波前快照；(b)t=3.6s時(shí)刻波前快照圖8 二階耦合準(zhǔn)P波方程數(shù)值模擬結(jié)果Figure 8 Second-order coupled quasi-P wave equation numerical simulation results