張千祥,王德利,周進(jìn)舉

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

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

張千祥,王德利,周進(jìn)舉

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

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

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

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

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

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

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

1 方法原理

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

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

(1)

(2)

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

(3)

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

1.2 傅里葉有限差分法

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

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

(4)

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

(5)

(6)

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

(7)

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

(8)

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

(9)

(10)

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

(11)

1.4 快速展開法

定義偽-差分算子L為：

(12)

方程(11)可以表示為：

(13)

(14)

進(jìn)而得出p-t項(xiàng)為：

(15)

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

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

(16)

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

(17)

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

(18)

式中：I是單位矩陣。

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

2 數(shù)值計(jì)算實(shí)例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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

[1] 陳洪杰.基于聲波方程的數(shù)值模擬與逆時(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] 張永剛.地震波場(chǎng)數(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] 周進(jìn)舉,王德利.含傾斜裂隙頁巖介質(zhì)地震波場(chǎng)傳播特征研究[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

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

(A1)

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

(A2)

其中

(A3)

方程(A2)有如下的解：

(A4)

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

(A5)

(A6)

v1是對(duì)稱面上的P波相速度,v2是垂直于對(duì)稱面方向上的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ì)波場(chǎng)正、反演理論和高精度地震勘探研究工作。

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

P631

A

1000-1441(2015)05-0485-08

10.3969/j.issn.1000-1441.2015.05.001