孔麗云， 胡志方*， 王一博

1 中國地質(zhì)調(diào)查局油氣資源調(diào)查中心， 北京 100083 2 中國科學院地質(zhì)與地球物理研究所， 北京 100029

0 引言

地球介質(zhì)中廣泛存在著各向異性，其縱橫波耦合的特性，使得基于各向同性假設的正演模擬方法不足以準確地描述實際地球介質(zhì)中地震波的傳播規(guī)律.在以往的研究中，考慮各向異性的彈性波動方程正演模擬取得了一定的進展.該方法是用三分量矢量場來描述彈性波波場，其正演結(jié)果能夠獲得比標量波場更豐富的波場信息.但其存在一定的缺陷(何兵壽和張會星，2006)：①各向異性參數(shù)復雜，物理意義不明確，在實際生產(chǎn)中難以獲得；②多分量數(shù)據(jù)占用資源大，計算效率低；③由于多分量地震勘探成本高，目前野外采集仍然以縱波地震資料為主，多分量彈性波理論的應用受到資料不足的限制.因此，需要對各向異性介質(zhì)的縱波分量進行提取，以適應目前階段地震數(shù)據(jù)采集和處理的需求.目前主要有兩種思路，一種是對介質(zhì)的彈性波場進行縱橫波波場的分離，另外一種就是通過頻散關系式對qP波近似方程進行推導.

縱橫波分離方面，常用的方法是Helmholtz分解(Yan and Sava，2008).其基本原理是基于縱橫波偏振方向的不同(縱波偏振方向平行于傳播方向，橫波偏振方向垂直于傳播方向)，利用波場的散度和旋度來分別對縱波和橫波進行求取.該算法在各向同性情況下適用，但在各向異性介質(zhì)中，縱、橫波的傳播方向發(fā)生了偏離，與極化方向不再是平行和垂直的關系，所以Helmholtz分解對于各向異性介質(zhì)不適用.Yan和Sava(2009)提出了基于空間域波場分離算子的各向異性介質(zhì)縱橫波分離算法，但該算法的計算效率太低，在3D情況下幾乎無法實現(xiàn).為了提高計算效率，Yan 和Sava(2011)提出利用不同權(quán)重的均勻反射模型對非均質(zhì)模型進行近似，這種方法實質(zhì)上是在計算效率和精度之間進行了一個權(quán)衡，但反射模型的選擇和計算代價都依賴于模型本身，難以應用于實際情況.Cheng 和Fomel(2014)轉(zhuǎn)換思路，將波場傳播方向和極化方向之間的投影關系重新改寫為空間波數(shù)域的積分算子，通過低階近似的方法對其進行數(shù)值求解.計算效率得到了部分提高，但仍然會隨著模型的尺寸和復雜度而降低.Zhou 和Wang(2017)通過將波場傳播方向局部旋轉(zhuǎn)到qP波的極化方向，發(fā)展了VTI介質(zhì)波場分離的新算子，有效提高了計算效率.另外，由于在分離過程中采用散度場代表縱波，旋度場代表橫波，這就意味著，縱波僅代表整體位移的輸入輸出，而橫波則代表旋轉(zhuǎn)變換的角速度，其物理意義與原波場并不相同，且其振幅和相位信息都發(fā)生了改變.在各向同性介質(zhì)中，Duan和Sava(2015)提出，可以通過在時間上對震源子波積分后乘以縱、橫波相速度的外推波場對分離后的波場進行矯正.在各向異性介質(zhì)中，由于波形改變變得更加復雜，沒有跟各向同性相似的矯正方法.Sun等(2001)通過Hilbert變換矯正相位的改變；Zhang和McMechan(2010)對VTI介質(zhì)下波數(shù)域二維和三維矢量波波場分離方法進行了研究；Sun等(2011)用縱橫波速度恢復 P波和S波的相對振幅比，通過地表P波速度和頻率恢復P和S波的真正振幅值.該方法雖然能恢復振幅與相位，但存儲空間耗費量巨大.

qP波近似方程方面，Alkhalifah從VTI介質(zhì)的精確qP-qSV波頻散關系出發(fā)，假設垂向剪切波速度Vs0=0，推導了VTI介質(zhì)四階擬聲波方程.為降低四階偏微分方程計算的復雜性(Alkhalifah，2000)，周熙焱等(2018)、Zhou等(2019)和Du等(2008)從聲學近似的頻散關系出發(fā)，通過引入不同的輔助波場函數(shù)，分別推導了兩種VTI介質(zhì)二階耦合擬聲波方程.為結(jié)合交錯網(wǎng)格技術(shù)實現(xiàn)精確的有限差分數(shù)值模擬，Hestholm推導了一階偏微分形式的VTI介質(zhì)擬聲波方程(Hestholm，2009).Duveneck等(2008)從聲學近似的胡克定律出發(fā)，也推導了二階形式和一階應力-速度形式的VTI介質(zhì)擬聲波方程.與Alkhalifah方程一樣，幾種擬聲波方程均存在低速度、低振幅的qSV人為干擾波問題(Grechka et al.，2004)，而這種干擾波的存在會影響正演模擬和偏移成像的最終結(jié)果.同時，由于無法考慮到轉(zhuǎn)換波的信息，利用qP波近似方程得到的各向異性介質(zhì)縱波數(shù)據(jù)并不完整.

為得到振幅和相位信息相對準確且波場信息相對完整的各向異性介質(zhì)縱波數(shù)據(jù)，本文采用波場分離的算法進行各向異性介質(zhì)的縱橫波分離，其基本思想是源于Zhang等(2007)提出的“方程解耦”的方法.該方法是直接將彈性剛度系數(shù)分量分解為縱波和橫波相關的系數(shù)，進而達到縱橫波方程分離的效果，因此不需要額外的計算.在各向同性介質(zhì)中，Zhang等(2007)得到精確的縱橫波分離的一階速度-應力方程；Zhou等(2016)將這種方法進一步拓展到三維各向同性模型的情況.在各向異性介質(zhì)中，縱橫波是耦合在一起的，這就使得其彈性波場的縱橫波分離需要在一定的近似條件下才能完成.Bloot等(2013)從Thomsen(1986)弱各向異性參數(shù)出發(fā)，通過假設δ=0(近垂直P波傳播中關鍵的各向異性參數(shù))，從彈性剛度系數(shù)矩陣中的耦合項C13中分離出縱波模量.

本文借助于Bloot等(2013)的近似方法，將彈性剛度系數(shù)分量中的耦合項C13近似為縱波模量和橫波模量的線性函數(shù)，推導得到了各向異性介質(zhì)的縱橫波分離的一階速度-應力方程，從而能夠?qū)崿F(xiàn)各向異性介質(zhì)的縱橫波波場分離.

1 傳統(tǒng)的波場分離算法

1.1 Helmholtz分解

對于縱橫波波場分離技術(shù)，本文所使用的是基于旋度場和散度場的交錯網(wǎng)格高階有限差分的縱橫波波場分離方法.首先，根據(jù)場的理論，任何一個矢量場都可定義為一個標量場的散度和一個矢量場的旋度相加的結(jié)果，因此，質(zhì)點速度v=(vx,vz)可以寫為如下表達形式：

(1)

其中，ψp和ψs分別代表速度場的標量位和矢量位，即縱波和橫波的位移位.其次，對于上述速度場分別求其散度和旋度，可以得到速度矢量的無旋場與無散場，也即縱波波場vp與橫波波場vs.其表達式分別為：

(2)

(3)

在二維情況下，上述方程可以簡化為：

(4)

(5)

從式(4)和式(5)中可以看出，Helmholtz分解得到的縱波是標量，在二維情況下，橫波也只有平行于y軸的一個分量，可視為標量.

1.2 聲波近似方程

首先，對VTI介質(zhì)的虎克定律進行經(jīng)典的聲波近似假設(Alkhalifah，2000)，令垂向剪切波速度vs0=0，得到用Thomsen參數(shù)表示的VTI介質(zhì)聲波近似胡克定律：

(6)

式中，σij(i,j=x,y,z)為應力分量，eij(i,j=x,y,z)為應變分量，ρ為密度，Vp0為VTI介質(zhì)對稱軸方向上的縱波速度，ε和δ為Thomsen各向異性參數(shù).

表達式(6)兩端對時間求導，并結(jié)合運動方程

ρ?vi/?t=?σij/?xii,j=x,y,z

(7)

可以得到一階應力-速度形式的VTI介質(zhì)qP波方程

(8)

式中，vx,vy,vz分別為沿x,y,z方向的速度分量，p=σ11=σ22和q=σ33為水平應力分量和垂直應力分量.

在二維情況下，方程(8)簡化為：

(9)

2 縱橫波方程解耦算法

2.1 基本思想

聲波和彈性波的區(qū)別在于其傳播介質(zhì)的不同，前者是在聲學介質(zhì)中傳播，只能傳播縱波；而后者則是在彈性空間中傳播，除了縱波還會產(chǎn)生橫波.而傳播介質(zhì)的不同主要體現(xiàn)在剛度系數(shù)矩陣上，因此，如果能將彈性剛度系數(shù)矩陣分離為分別只與縱波和橫波相關的兩部分，那么就可以將整個彈性波方程分解為縱波方程和橫波方程，從而達到縱橫波分離的目的.基于這一思想，給出適用于交錯網(wǎng)格高階有限差分算法的各向異性縱橫波分離方程的具體推導過程.

2.1.1 二階波動方程

三維彈性波動方程通常表示為：

(10)

式中，ρ為介質(zhì)密度；u=(ux,uy,uz)T為位移矢量；C為剛度系數(shù)矩陣；G為偏微分算子矩陣.

根據(jù)Zhang等(2007)的思想，假設完整的彈性系數(shù)矩陣C可以分為與縱波和橫波分別相關的彈性系數(shù)矩陣Cp和Cs兩部分，那么彈性波動方程就可以拆分為：

(11)

其中，up和us分別為縱波波場和橫波波場，全波場u=up+us.

2.1.2 一階速度-應力方程

胡克定律表達式為：

T=CE=CGTu,(12)

其中，T=(σxx,σyy,σzz,σyz,σxz,σxy)T為應力向量，E=(exx,eyy,ezz,eyz,exz,exy)T為應變向量.將上述表達式代入方程(11)中，可以得到縱橫波分離的一階速度-應力方程表達式：

(13)

從上述推導過程中可以看出，如何將完整的彈性系數(shù)矩陣C分解為分別只與縱波和橫波相關的Cp和Cs兩部分，是該方法的核心問題.對于各向同性介質(zhì)，由于縱波和橫波能夠完全解耦，所以能夠得到完全精確的縱波和橫波方程(Zhang et al.，2007).但是在各向異性介質(zhì)中，由于縱橫波無法解耦，我們需要對相關系數(shù)進行一定的假設，才能將彈性系數(shù)矩陣C進行分解，下面給出各向異性VTI介質(zhì)彈性系數(shù)矩陣分解的過程.

2.2 彈性系數(shù)矩陣的分解

2.2.1C13的線性近似表達

對于VTI介質(zhì)，剛度系數(shù)不為零分量用Thomsen參數(shù)可以表達為：

其中，Vp0、Vs0分別為qP波和qSV波沿VTI介質(zhì)對稱軸方向的相速度；ε、δ和γ是表示VTI介質(zhì)各向異性強度的三個無量綱因子：ε是度量qP波各向異性強度參數(shù)，ε越大，介質(zhì)的縱波各向異性越大，ε=0，縱波無各向異性；δ是連接Vp0和Vp90之間的一種過渡性參數(shù)；γ可以看成是度量qS波各向異性強度或橫波分裂強度的參數(shù)，γ越大，介質(zhì)的橫波各向異性越大，γ=0時，橫波無各向異性.

(15)

從(15)式可以看出，系數(shù)A、B的求解是C13和C23分解的關鍵.由于C13和C23的表達式相同，以C13為例，給出其系數(shù)A和B的近似求解過程.

2.2.2 系數(shù)的求取

對式(14)中C13的表達式進行重新整理得：

(16)

為求取橫波分量系數(shù)B，首先假設δ=0，得到

進一步，假設Vp0=0，可得：

B=-2.

(17)

將(17)式代入C13近似表達式(15)中，就可以計算得到系數(shù)A為：

(18)

2.2.3Cp和Cs近似表達式

將系數(shù)表達式(17)和(18)代入C13的近似表達式(15)中，結(jié)合式(14)中其他分量的表達式，可以得到VTI介質(zhì)中的縱波彈性模量Cp和橫波模量Cs的分量表達式分別為：

(19)

(20)

將表達式(19)、(20)代入方程(13)中，可以得到VTI介質(zhì)在三維情況下的縱橫波分離的一階速度-應力方程.

2.3 二維分解公式

在二維情況下，VTI介質(zhì)中的一階速度-應力方程表達式為：

(21)

其中，

D11m=(C13C13m-C33C11m)M-1,

D12m=(C13C11m-C11C13m)M-1,

D13m=C55m/C55,

D21m=(C13C33m-C33C13m)M-1(m=p,s)

D22m=(C13C13m-C11C33m)M-1,

D23m=C55m/C55,

(22)

3 理論模型測試

為說明“方程解耦”算法的精確性和可靠性，基于x和z兩個方向的地震波場，利用二維公式(22)對三層水平層狀介質(zhì)模型進行測試，并將波場分離結(jié)果與同樣考慮x、z兩個方向地震波場的Helmholtz分解(式(4)、(5))和擬聲波近似(式(9))等傳統(tǒng)算法進行比較.

3.1 理論模型及測試

設計三層水平層狀介質(zhì)模型如圖1所示，中間層為VTI介質(zhì)，其頂層和底層為各向同性彈性介質(zhì)模型，介質(zhì)參數(shù)如表1 所示.模型大小為704 m×704 m×528 m，空間網(wǎng)格間距為8 m；模擬的物理時間為0.32 s，時間步長為0.4 ms.震源采用爆炸震源，主頻25 Hz，波源函數(shù)表達式為：

(23)

其中，h(t)為雷克子波，(x0,z0)為震源中心點，α決定了震源子波的寬度，f0為主頻，t0為延遲時間.

表1 三層模型物理參數(shù)表Table 1 Physical parameters for three-layer model

圖1 三層模型示意圖Fig.1 Diagram of three-layer model

圖2 爆炸震源激發(fā)的縱波路徑示意圖(去除直達波)，包含入射縱波-反射縱波(P-P，r1，橙色)，入射縱波-透射轉(zhuǎn)換橫波-反射橫波-透射轉(zhuǎn)換縱波(P-S-S-P，r2，綠色)，入射縱波-透射轉(zhuǎn)換橫波-反射轉(zhuǎn)換縱波-透射縱波(P-S-P-P，r3，紅色)，入射縱波-透射縱波-反射轉(zhuǎn)換橫波-透射轉(zhuǎn)換縱波(P-P-S-P，r3，紅色)，入射縱波-透射縱波-反射縱波-透射縱波(P-P-P-P，r4，藍色)Fig.2 Diagram of P-wave path excited by explosion source (removing direct wave), including incident P-refeccted P (P-P, r1, orange), incident P-transmitted converted S-reflected S-transmitted converted P(P-S-S-P, r2, green), incident P-transmitted converted S-reflected converted P-transmitted P (P-S-P-P, r3, red), incident P-transmitted P-reflected converted S-transmitted converted P (P-P-S-P, r3, red), incident P-transmitted P-reflected P-transmitted P (P-P-P-P, r4, blue) (P is longitudinal wave and S is transverse wave)

3.2 波場的完整性比較

圖2和圖3給出了爆炸震源激發(fā)后，在接收點接收到的來自I、II兩個地層界面的縱橫波射線路徑示意圖.從圖中可以看出，縱波和橫波都存在四條反射路徑(其中紅色路徑是兩條路徑的重合，由于到達時間相同，在這里用一條路徑代替).所以，在分離后的地震剖面上，應該存在四條同相軸.圖4和圖5分別為基于方程解耦算法在x方向和z方向的縱橫波結(jié)果，為了說明該算法的有效性，將其與Helmholtz分解(圖6)和各向異性擬聲波方程(圖7)的縱橫波分離結(jié)果進行對比.由于Helmholtz分解算法在二維情況下得到的縱波和橫波都是標量，圖6中只顯示了縱波和橫波的綜合分離結(jié)果，沒有方向性.另外，由于擬聲波近似算法的結(jié)果中只有縱波信息，所以圖9中只給出了縱波分別在x方向和z方向的分離結(jié)果.

圖3 爆炸震源激發(fā)的橫波路徑示意圖(去除直達波)，包含入射縱波-反射轉(zhuǎn)換橫波(P-S，r′1，橙色)，入射縱波-透射轉(zhuǎn)換橫波-反射橫波-透射橫波(P-S-S-S，r′2，綠色)，入射縱波-透射轉(zhuǎn)換橫波-反射轉(zhuǎn)換縱波-透射轉(zhuǎn)換橫波(P-S-P-S，r′3，紅色)，入射縱波-透射縱波-反射轉(zhuǎn)換橫波-透射橫波(P-P-S-S，r′3，紅色)，入射縱波-透射縱波-反射縱波-透射轉(zhuǎn)換橫波(P-P-P-S，r′4，藍色)Fig.3 Diagram of S-wave path excited by explosion source (removing direct wave), including incident P-reflected converted S (P-S, r′1, orange), incident P-transmitted converted S-reflected S-transmitted S (P-S-S-S, r′2, green), incident P-transmitted converted S-reflected converted P-transmitted converted S (P-S-P-S, r′3, red), incident P-transmitted P-reflected converted S-transmitted S (P-P-S-S, r′3, red), incident P-transmitted P-reflected P-transmitted convered S (P-P-P-S, r′4, blue)(P is longitudinal wave and S is transverse wave)

圖4 基于方程解耦算法的數(shù)值模擬結(jié)果(x分量)(a) 完整的波場； (b) 縱波波場； (c) 橫波波場.Fig.4 Numerical simulation results based on equation decoupling algorithm (x component)(a) Complete wave field; (b) P-wave; (c) S-wave.

圖5 基于方程解耦算法的數(shù)值模擬結(jié)果(z分量)(a) 完整的波場； (b) 縱波波場； (c) 橫波波場.Fig.5 Numerical simulation results based on equation decoupling algorithm (z component)(a) Complete wave field; (b) P-wave; (c) S-wave.

圖6 基于Helmholtz分解的數(shù)值模擬結(jié)果(a) 縱波波場； (b) 橫波波場.Fig.6 Numerical simulation results based on Helmholtz decomposition(a) P-wave; (b) S-wave.

圖7 基于擬聲波近似的縱波模擬結(jié)果(a) x分量； (b) z分量.Fig.7 Numerical simulation results of longitudinal wave based on acoustic approximation(a) x component; (b) z component.

通過對比可以看出，基于方程解耦算法(圖4和圖5)和基于Helmholtz分解算法得到的波場(圖6)，都存在四條反射波同相軸，而基于擬聲波近似方法得到的波場(圖7)中，只存在兩條反射波同相軸.因此，基于方程解耦算法和Helmholtz分解算法得到的縱波數(shù)據(jù)，都能夠得到完整的反射波場.

3.3 振幅和相位的精確性比較

為了說明方程解耦算法在地震波振幅和相位方面的精確性，分別取方程解耦算法、Helmholtz分解、擬聲波近似等三種方法數(shù)值模擬結(jié)果中的第100道地震數(shù)據(jù)進行對比.圖8—11分別是x方向縱波分量、z方向縱波分量、x方向橫波分量、z方向橫波分量與相對應全波場信息的比較.為方便對比，將Helmholtz分解的標量結(jié)果同樣顯示在x方向和z方向.

對比縱波(圖8和圖9)可以看出，三種方法在地震波走時上與全波場都保持一致.方程解耦算法(紅色實線)和擬聲波近似方法(紫色實線)的振幅和相位信息與全波場(黑色虛線)基本保持一致，但是Helmholtz分解算法(藍色實線)得到的縱波信息，無論在振幅還是相位方面，都與全波場(黑色虛線)存在很大誤差.在橫波的對比結(jié)果中(圖10和圖11)同樣可以看出，相對于Helmholtz分解算法(藍色實線)，方程解耦算法在橫波振幅和相位上與原始波場保持高度一致性.因此，基于方程解耦算法得到的縱波數(shù)據(jù)，能夠得到準確的振幅和相位信息.

圖8 不同方法(藍實線：Helmholtz分解；紫實線：擬聲波方程；紅實線：方程解耦)得到的縱波信息與全波場信息(黑色虛線)的對比(x分量)(a) 0.1～0.38 s； (b) 0.38～0.64 s.Fig.8 Comparison of full wavefield (black dotted line) and P-wave obtained by different methods (blue solid line: Helmholtz decomposition; purple solid line: acoustic approximation; red solid line: equations decoupling)(x component)

圖9 不同方法(藍實線：Helmholtz分解；紫實線：擬聲波方程；紅實線：介質(zhì)分離)得到的縱波信息與全波場信息(黑色虛線)的對比(z分量)(a) 0.1～0.38 s； (b) 0.38～0.64 s.Fig.9 Comparison of full wavefield (black dotted line) and P-wave obtained by different methods (blue solid line: Helmholtz decomposition; purple solid line: acoustic approximation; red solid line: equations decoupling) (z component)

圖10 不同方法(藍實線：Helmholtz分解；紅實線：介質(zhì)分離) 得到的橫波信息與全波場信息(黑色虛線)的對比(x分量)(a) 0.1～0.38 s； (b) 0.38～0.64 s.Fig.10 Comparison of full wavefield (black dotted line) and S-wave obtained by different methods (blue solid line: Helmholtz decomposition; red solid line: equations decoupling) (x component)

圖11 不同方法(藍實線：Helmholtz分解；紅實線：介質(zhì)分離) 得到的橫波信息與全波場信息 (黑色虛線)的對比(z分量)(a) 0.1～0.38 s； (b) 0.38～0.64 s.Fig.11 Comparison of full wavefield (black dotted line) and S-wave obtained by different methods (blue solid line: Helmholtz decomposition; red solid line: equations decoupling)(z component)

4 結(jié)論

基于方程解耦的思想，通過將彈性剛度系數(shù)分量中的耦合項C13近似表示為縱波模量和橫波模量的線性函數(shù)，推導得到了各向異性介質(zhì)的縱橫波分離的一階速度-應力方程，從而能夠?qū)崿F(xiàn)各向異性介質(zhì)的縱橫波波場分離.與擬聲波方程算法相比，方程解耦算法得到的縱波波場信息較為完整；與Helmhotlz分解算法相比，方程解耦算法得到的縱波波場的振幅和相位信息都相對準確.因此，對于各向異性介質(zhì)而言，方程解耦算法既能保證分離后波場的完整性，又能保持波場振幅、相位信息的準確性，為各向異性非常規(guī)儲層地震數(shù)據(jù)處理、解釋的可靠性提供了重要的數(shù)據(jù)保障.