張千祥,王德利,周進舉

(吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，吉林長春 130026)

二維TTI介質(zhì)的純P波波動方程數(shù)值模擬

張千祥,王德利,周進舉

(吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，吉林長春 130026)

聲波各向異性數(shù)值模擬對地震數(shù)據(jù)處理和解釋起著重要的作用。基于Tsvankin提出的精確色散關(guān)系,通過平方根近似,在時間-波數(shù)域中推導(dǎo)出二維TTI介質(zhì)純P波聲波波動方程,并利用快速展開法(Rapid Expansion Method,REM)進行了數(shù)值模擬。與傳統(tǒng)的有限差分法求解二維TTI介質(zhì)耦合方程和傅里葉有限差分法在時間上進行波場外推相比,該方法的模擬結(jié)果精度更高,計算速度更快,并且成功去除橫波分量。

聲波各向異性數(shù)值模擬;純P波聲波方程;快速展開法;有限差分法;傅里葉有限差分法

地震數(shù)值模擬是地震勘探方法研究的前提和基礎(chǔ),在地震勘探和地震學(xué)的各項研究及生產(chǎn)工作中都扮演著重要的角色[1]。常用的地震波場數(shù)值模擬方法主要有幾何射線法、波動方程法和積分方程法[2]。波動方程模擬方法中的有限差分法由于計算速度快、占用計算機內(nèi)存小而被廣泛應(yīng)用。很早時候各國地球物理學(xué)家就對各向同性介質(zhì)和各向異性介質(zhì)彈性波地震數(shù)值模擬進行了深入研究。近年來,周進舉等[3]利用高階旋轉(zhuǎn)網(wǎng)格有限差分法研究了復(fù)雜介質(zhì)下彈性波數(shù)值模擬。

在對地下的各向異性介質(zhì)進行彈性波數(shù)值模擬時,由于彈性波方程復(fù)雜,各向異性參數(shù)多,導(dǎo)致模擬計算量大,耗時長,增加了彈性波偏移和干涉的難度。為了解決這些問題,我們采用聲波各向同性近似理論,通過設(shè)定彈性波中的橫波速度為零來簡化計算參量,在保證模擬精度的條件下提高計算效率。然而,由于地下介質(zhì)的不均勻性,這種各向同性的聲學(xué)假設(shè)常常是不恰當(dāng)?shù)摹Ｒ虼?在各向異性介質(zhì)的情況下,需要提出一個合理的正演模擬方法來提高聲波各向異性數(shù)值模擬的精度,并顯著縮短計算時間。

多年來,人們對于各向異性介質(zhì)發(fā)展了多種P波正演算子[4-5],但是這些正演算子大多針對具有垂直對稱軸的各向同性介質(zhì)(Transversely isotropic media with a vertical symmetry axis,VTI介質(zhì))的。VTI介質(zhì)在沉積盆地中廣泛存在[6],然而,這種假設(shè)只對一些簡單的地質(zhì)形態(tài)有效,例如介質(zhì)中的巖層面平行于記錄面。當(dāng)?shù)叵碌膸r層面傾斜時,各向同性介質(zhì)(Transversely isotropic media,TI介質(zhì))的對稱軸可能不垂直,這種介質(zhì)通常被稱為有傾斜對稱軸的各向同性介質(zhì)(Transversely isotropic media with a tilt symmetry axis,TTI介質(zhì)),通常在有背斜構(gòu)造或者逆沖推覆體的區(qū)域廣泛發(fā)育。自從Alkhalifah等[7]在VTI介質(zhì)研究中引入了“聲學(xué)”近似,許多地球物理學(xué)家在VTI介質(zhì)的建模和正演中進行了更深入研究,并逐漸向TTI介質(zhì)擴展,發(fā)展了許多形式不同的二維TTI介質(zhì)聲波各異性波動方程和各種高精度、低頻散的數(shù)值模擬方法。2001年Klie等[8],2006年Zhou等[9],2007年Hestholm[10]都基于Alkhalifah聲學(xué)假設(shè)對VTI介質(zhì)進行了聲波數(shù)值模擬。之后,2009年Fletcher等[11]通過旋轉(zhuǎn)對稱軸導(dǎo)出了TTI介質(zhì)聲波波動方程。近年來,Barrera等[12]提出了不含有橫波的純P波波動方程并成功的應(yīng)用到逆時偏移中。Kim等[13]在前人研究的基礎(chǔ)上對純P波波動方程在數(shù)值模擬實現(xiàn)上進行了改進,采用圖形處理器(Graphic processing unit,GPU)提高了計算速度。

我們從Tsvankin[14]提出的精確色散關(guān)系出發(fā),得到近似P波色散關(guān)系,并構(gòu)建一個二維TTI介質(zhì)的純P波模式的擬微分聲波方程,然后應(yīng)用快速展開法進行數(shù)值實現(xiàn)。與傳統(tǒng)的有限差分法(Finite difference method,FDM)[15]和傅里葉有限差分法(Fourier finite differences,FFD)[16]相比,該方法不僅能很好去除頻散，提高計算效率,還能去除SV波分量。

1 方法原理

1.1 二維TTI介質(zhì)耦合方程

首先從Tsvankin[17]提出的精確相速度表達(dá)式開始推導(dǎo):

(1)

(2)

其中微分算子H1和H2被定義為：

(3)

對耦合方程(2)中的時間導(dǎo)數(shù)和空間導(dǎo)數(shù)分別用有限差分法進行求解,但這種方法會隨著有限差分項的階數(shù)變高而變得繁瑣。因此,我們結(jié)合傅里葉變換和有限差分法來簡化方程的推導(dǎo)和求解。

1.2 傅里葉有限差分法

我們提出了一個波場延拓的新方法,該方法結(jié)合了傅里葉變換和有限差分法。具體推導(dǎo)過程見附錄A。對于復(fù)雜速度模型的正演模擬,該方法具有較高的精確性和穩(wěn)定性;對于二維TTI介質(zhì)的數(shù)值模擬,該方法能很好的去除橫波分量,但容易出現(xiàn)頻散效應(yīng),模擬效果不是最佳。我們從Tsvankin提出的精確相速度開始推導(dǎo),得出純P波聲波方程,并用快速展開法進行求解,取得了很好的效果。

1.3 二維TTI介質(zhì)的純P波波動方程

(4)

方程(4)對于P波色散關(guān)系式是一個很好的近似[19-20],當(dāng)滿足：

(5)

(6)

其中F=1+2ε/f,令ε=0,方程(6)簡化為：

(7)

方程(7)對于VTI介質(zhì)也是成立的。對于對稱軸有任意取向的TTI介質(zhì),可以由方程(7)通過一個描述繞Z軸逆時針旋轉(zhuǎn)的變量推導(dǎo)出來。在旋轉(zhuǎn)的坐標(biāo)系下,波速算子為：

(8)

其中，θ和φ是互余的關(guān)系,從方程(8)我們可以得出：

(9)

(10)

把方程(10)的兩邊在傅里葉域中同時乘以波場函數(shù)p(ω,kx,kz),然后,對方程兩邊同時用傅里葉逆變換,之后由關(guān)系式iω??/?t,最終得出二維TTI介質(zhì)在時間-波數(shù)域中純P波波動方程：

(11)

1.4 快速展開法

定義偽-差分算子L為：

(12)

方程(11)可以表示為：

(13)

(14)

進而得出p-t項為：

(15)

將方程(14)和方程(15)相加,應(yīng)用快速展開法進行數(shù)值求解,求解方程為：

p(t+Δt)=-p(t-Δt)+2cos(LΔt)p(t)

(16)

在方程(16)中余弦函數(shù)的正交多項式展開式可表示為[21]：

(17)

其中,當(dāng)k=0時,C2k=1;當(dāng)k>0時,C2k=2。R是L2的最大特征值。J2k為第一類貝塞爾函數(shù),Q2k為切比雪夫多項式,滿足如下的循環(huán)關(guān)系：

(18)

式中：I是單位矩陣。

當(dāng)M>RΔt時,方程(17)以指數(shù)形式收斂。因此,當(dāng)M的取值略微大于RΔt時,方程的總求和項可以被大大的縮短。Pestana等[19]已經(jīng)證實,當(dāng)M=1時,也就是說求和中只有兩項,這種使用切比雪夫多項式的余弦函數(shù)的近似相當(dāng)于在時間上的二階有限差分格式。當(dāng)M=2時,包含了L4算子項,這種近似等價于Dablain[22]提出的時間上的四階有限差分格式。

2 數(shù)值計算實例

為了驗證本文提出的純P波聲波方程的有效性,我們設(shè)計了一個模型并對其進行數(shù)值模擬。模型大小為4000m×4000m,震源是雷克子波,主頻30Hz,網(wǎng)格間距為Δx=Δz=10m,時間采樣間隔為0.1ms,vP0=3000m/s,ε=0.24,δ=0.1。

圖1a到圖1d分別代表θ為0,30°,60°,90°時所對應(yīng)的波場快照。從圖1可以看出數(shù)值模擬效果很好,成功去除了橫波,證明了本文方程的有效性。

接下來通過與傳統(tǒng)有限差分法求解聲波耦合方程和傅里葉有限差分法波場數(shù)值模擬的比較來證明快速展開法的優(yōu)勢。模擬時,震源子波采用主頻30Hz的雷克子波,傳播速度vP0=3000m/s,Thomsen參數(shù)ε=0.20,δ=0.05,模型大小同樣為4000m×4000m。

圖2a對應(yīng)有限差分法求解TTI耦合方程(2)得到的波場快照，空間網(wǎng)格間距10m×10m;圖2b 對應(yīng)傅里葉有限差分法進行時間上的波場外推得到的波場快照，空間網(wǎng)格間距5m×5m;圖2c 為應(yīng)用快速展開法在時間—波數(shù)域中求解二維TTI介質(zhì)純P波聲波方程(11)得到的波場快照，空間網(wǎng)格間距10m×10m。

圖1 二維TTI介質(zhì)在傳播時間為0.5s的波場快照

圖2 不同方法求解出的二維TTI介質(zhì)在傳播時間為0.6s的波場快照

有限差分法求解TTI耦合方程時，時間上采用二階、空間上采用六階有限差分形式，計算時間為273s，而應(yīng)用快速展開法的計算時間只需要147s。對比圖2a和圖2c可以看出，應(yīng)用快速展開法不但很好地去除了橫波分量，還大大縮短了計算時間。對比圖2b與圖2c發(fā)現(xiàn),雖然傅里葉有限差分波場外推也能很好地去除橫波分量,但計算時采用的網(wǎng)格間距小,圖2b上頻散效應(yīng)還很嚴(yán)重。當(dāng)兩種方法參數(shù)相同時,傅里葉有限差分法計算時間為156s。兩種方法都是在頻率域計算,模擬時間相差不多。綜上所述,本文采用的快速展開方法進行數(shù)值模擬具有很強的優(yōu)越性。

利用本文方法對三層均勻介質(zhì)模型進行正演模擬。速度模型如圖3所示。

震源子波采用主頻為40Hz的雷克子波,時間步長0.2ms,網(wǎng)格間距給定為15m×15m,傳播時間為0.5s時的波場快照如圖4所示。

圖5a和圖5b分別給出了合成的BP速度模型以及應(yīng)用快速展開法進行數(shù)值模擬的結(jié)果。模擬時采用峰值頻率為50Hz的雷克子波,水平網(wǎng)格間距為20m,垂直網(wǎng)格間距為10m,時間采樣間隔為0.5ms,傳播時間為1.5s。

圖3 三層均勻介質(zhì)速度模型

圖4 采用本文方法對三層均勻介質(zhì)模型正演模擬波場快照

由圖4可看出，在地下2.0km和3.1km處存在明顯的反射和透射，并且沒有橫波波前面，與圖3所示的速度模型相吻合。在圖5b中也顯示出了地下介質(zhì)中高速體的存在。由此可知，本文所提出的方法可以很好地應(yīng)用于復(fù)雜模型的正演模擬。

圖5 BP速度模型(a)和本文提出的快速展開法數(shù)值模擬結(jié)果(b)

3 結(jié)論

1) 利用傳統(tǒng)的有限差分法進行數(shù)值模擬時,我們只是假定了沿對稱軸的剪切波速度為零,而不是在任意方向剪切波的速度都為零。所以快照圖中還會產(chǎn)生橫波分量,且出現(xiàn)一定的頻散效應(yīng),加大有限差分的階數(shù)減小頻散時,又會大大地增加計算量和計算時間。

2) 通過傅里葉有限差分法對地震波場在時間上進行外推時,雖然橫波分量被很好地去除掉,但在網(wǎng)格間距很小的情況下頻散效應(yīng)還是很嚴(yán)重,效果不理想,不適合應(yīng)用到實際情況下的數(shù)值模擬。

3) 通過快速展開法在時間-波數(shù)域?qū)ΧSTTI介質(zhì)純P波聲波方程進行數(shù)值模擬時,有效去除了橫波分量,縮短了計算時間,提高了計算效率；網(wǎng)格間距很大時也能取得良好的模擬效果,沒有頻散效應(yīng)。所以當(dāng)模擬區(qū)域比較大的時候,在保證相同的模擬效果情況下能大大地縮短模擬時間。

4) 本文的研究成果也能很好地應(yīng)用于復(fù)雜介質(zhì)的正演模擬,為以后TTI介質(zhì)的逆時偏移(RTM)、聲波干涉和全波形反演(FWI)等提供了可靠的技術(shù)手段。

[1] 陳洪杰.基于聲波方程的數(shù)值模擬與逆時偏移方法研究[D].吉林大學(xué),2010 Chen H J.Numerical simulation and reverse-time migration method research base on acoustic wave equation[D].Jilin University,2010

[2] 張永剛.地震波場數(shù)值模擬方法[J].石油物探,2003,42(2)：143-148 Zhang Y G.On numerical simulations of seismic wavefield[J].Geophysical Prospecting for Petroleum,2003,42(2):143-148

[3] 周進舉,王德利.含傾斜裂隙頁巖介質(zhì)地震波場傳播特征研究[J].石油物探,2014,53(5):501-508 Zhou J J,Wang D L.Research of wave field propagation characteristics of shale medium containing tilt fracture[J].Geophysical Prospecting for Petroleum,2014,53(5):501-508

[4] Han Q,Wu R S.A one-way dual-domain propagator for scalar qP-waves in VTI media[J].Geophysics,2005,70(2):D9-D17

[5] Zhang L,Hua B,Calandra H.3D Fourier finite-difference anisotropic depth migration[J].Expanded Abstracts of 75thAnnual Internat SEG Mtg,2005,1914-1917

[6] Crampin S,Chesnokov E M,Hipkin R A.Seismic anisotropy—the state of the art[J].First Break,1984,20(3):9-18

[7] Alkhalifah T.An acoustic wave equation for anisotropic media[J].Geophysics,2000,65(11):1239-1250

[8] Klie H,Toro W.A new acoustic wave equation for modeling in anisotropic media[J].Expanded Abstracts of 71stAnnual Internat SEG Mtg,2001,1171-1174

[9] Zhou H,Zhang G,Bloor R.An anisotropic acoustic wave equation for VTI media[J].Expanded Abstracts of 68thEAGE Annual Conference,2006,H033

[10] Hestholm S.Acoustic VTI modeling using high-order finite-differences[J].Expanded Abstracts of 77thAnnual Internat SEG Mtg,2007,139-143

[11] Fletcher R P,Du X,Fowler P J.Reverse-time migration in tilted transversely isotropic(TTI) media[J].Geophysics,2009,74( 6):WCA179-WCA187

[12] Barrera D F,Pestana R C.Pseudo-acoustic wave equation for modeling and reverse time migration in TTI media[J].Journal of Seismic Exploration.2013,22(1):33-48

[13] Kim Y S.Acceleration of stable TTI P-wave reverse-time migration with GPUs[J].Computers & Geosciences,2013,52(1):204-217

[14] Tsvankin I.Seismic signatures and analysis of reflection data in anisotropic media[M].Tulsa:SEG,2012,1-433

[15] Zhou H B,Zhang G Q,Bloor R.An anisotropic acoustic wave equation for modeling and migration in 2D TTI media[J].Expanded Abstracts of 76thAnnual Internat SEG Mtg,2006,194-198

[16] Song X L,Fomel S.Fourier finite-difference wave propagation[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010,3204-3209

[17] Tsvankin I.P-wave signatures and notation for transversely isotropic media:an overview[J].Geophysics,1996,61(2):467-483

[18] Thomsen L.Weak elastic anisotropy[J]:Geophysics,1986,51(11):1954-1966

[19] Pestana R C,Stoffa P L.Time evolution of the wave equation using rapid expansion method[J].Geophysics,2010,75( 4):T121-T131

[20] Pestana R C,Ursin B,Stoffa P L.Separate P- and SV-wave equations for VTI media[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011,163-167

[21] Tal-Ezer H,Kosloff D,Koren Z.An accurate scheme for seismic forward modeling[J].Geophysical Prospecting,1987,35( 5):479-490

[22] Dablain M A.The application of high-order differencing to the scalar wave equation[J].Geophysics,1986,51(1):54-66

附錄A

對各向同性介質(zhì)有：

(A1)

其中p(x,t)是地震壓力場,v(x)是傳播速度。假定一個常速度v,經(jīng)過空間中的傅里葉變換,我們得到如下的表達(dá)式：

(A2)

其中

(A3)

方程(A2)有如下的解：

(A4)

對時間應(yīng)用二階有限差分和逆傅里葉變換，得到：

(A5)

(A6)

v1是對稱面上的P波相速度,v2是垂直于對稱面方向上的P波相速度,η是各向異性彈性參數(shù),它與Thomsen彈性參數(shù)ε和δ有關(guān)：

(A7)

4) 將通過泰勒展開得到的系數(shù)應(yīng)用到q(x,)來得到q(x,t+Δt);

5)q(x,t+Δt)+2p(x,t)-p(x,t-Δt)→p(x,t+Δt)。

(編輯：朱文杰)

Acoustic wave equation numerical simulation for pure P-wave in 2D TTI medium

Zhang Qianxiang,Wang Deli,Zhou Jinju

(JilinUniversity,Changchun130026,China)

Anisotropic numerical simulation of acoustic wave plays an important role in seismic data processing and interpretation.Starting with the exact dispersion relationship derived by Tsvankin and a square root approximation，we proposed a pure P-wave acoustic equation for 2D TTI medium in time-wavenumber domain,then the rapid expansion method (REM) is applied for numerical simulation.Compared with the conventional finite difference method to solve the 2D TTI medium coupled equations and Fourier finite difference to implement wavefield extrapolation in time,the synthetic data test results of our method is advanced in high precision and fast computing speed; what’s more,the SV-wave component has been removed successfully.

anisotropic numerical simulation of acoustic wave,pure P-wave acoustic equation,rapid expansion method,finite difference method,Fourier finite difference method

2015-01-14;改回日期：2015-03-30。

張千祥(1992—),男,碩士,現(xiàn)從事聲波各向異性正演數(shù)值模擬研究工作。

王德利(1973—),男,教授,博士生導(dǎo)師,主要從事各向異性介質(zhì)波場正、反演理論和高精度地震勘探研究工作。

國家自然科學(xué)基金項目(41374108)和國家科技重大專項(2011ZX05023-005-008)聯(lián)合資助。

P631

A

1000-1441(2015)05-0485-08

10.3969/j.issn.1000-1441.2015.05.001