齊 鵬，孫建國

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

近年來，隨著計算機硬件水平的不斷提高，地震波數(shù)值模擬進入高精度時代，三維數(shù)值模擬已逐漸成為常規(guī)的地震波場數(shù)值模擬技術(shù)。但是，在處理大尺度數(shù)值模擬或是快速計算時，由于計算效率的限制，三維數(shù)值模擬的高效實現(xiàn)仍存在一些困難。而傳統(tǒng)的二維數(shù)值模擬由于未考慮三維幾何擴散，導(dǎo)致模擬的地震波振幅信息存在較大誤差，只能較準(zhǔn)確地反映地震波的運動學(xué)特性，模擬精度較低。為了在保證計算效率的同時又能保證模擬的準(zhǔn)確性，2.5維數(shù)值模擬得到了很大的發(fā)展[1]。2.5維數(shù)值模擬假設(shè)三維介質(zhì)沿走向方向不變，使用三維震源在計算平面內(nèi)研究地震波的傳播，可以得到較精確的地震波動力學(xué)信息，特別是針對具有沿走向方向穩(wěn)定變化特點的沉積巖層褶皺、斷層或其它含油氣地質(zhì)構(gòu)造，作為三維問題的一種特殊情況，2.5維模擬是一種很好的折中手段[3]。

鑒于聲波模型是反射地震偏移成像理論和應(yīng)用研究中的基本數(shù)據(jù)模型，所以對2.5維聲波方程數(shù)值模擬的研究具有重要意義。目前，2.5維聲波方程數(shù)值模擬主要采用兩種方法：Fourier變換法[2]和近似方程法[3-6]。其中，F(xiàn)ourier變換法是對于三維聲波方程的沿走向方向做Fourier變換，得到關(guān)于波數(shù)k的二維聲波方程，對不同波數(shù)k的二維方程進行數(shù)值模擬，再通過疊加得到2.5維數(shù)值模擬的結(jié)果，相當(dāng)于進行多次二維數(shù)值模擬，計算量較大；而近似方程法是通過某些近似條件推導(dǎo)出不同形式的2.5維的聲波近似方程，再求解近似方程，相當(dāng)于一次“二維數(shù)值模擬”，相對于Fourier變換法具有更高的計算效率?？紤]到近似方程法在計算速度上的優(yōu)越性，本文中采用近似方程法。

針對2.5維聲波近似方程的推導(dǎo)，人們利用不同的近似條件得到了不同形式的近似方程。Liner[3]將勻速介質(zhì)中的三維Green函數(shù)帶入二維聲波方程，得到常密度條件下L型2.5維聲波近似方程。Williamson等[4]對二維波場施加一個濾波器，獲得與L型方程末項不同的W型2.5維聲波近似方程。Narayan[5]通過對比L型方程和W型方程的數(shù)值模擬結(jié)果，得出兩個方程末項可缺省的結(jié)論，得到DL型方程。我們選用Liner提出的L型近似方程進行起伏地表條件下2.5維聲波數(shù)值模擬。

此外，在實際的地震勘探工作中，經(jīng)常會遇到復(fù)雜條件的起伏地表問題，起伏地表下地震波傳播問題的研究在地震學(xué)和勘探地震學(xué)領(lǐng)域具有重要意義[7-9]。起伏地表的地震波數(shù)值模擬方法很多，其中處理效果較好的有限元法[10]存在計算效率不高、網(wǎng)格生成復(fù)雜的缺點，邊界元法不適合模擬復(fù)雜地質(zhì)模型；而Tessmer等[11]提出了坐標(biāo)變換法可以將有限差分法應(yīng)用于起伏地表的數(shù)值模擬，具有很高的計算效率，而且能保證精度；Hestholm等[12-13]利用坐標(biāo)變換法實現(xiàn)了基于應(yīng)力-速度一階彈性波方程的起伏地表二維和三維數(shù)值模擬；唐文等[14]實現(xiàn)了基于二階彈性波方程的起伏地表數(shù)值模擬；王祥春等[15]采用類似的坐標(biāo)變換實現(xiàn)了二維聲波方程的起伏地表數(shù)值模擬。

相對于應(yīng)力-速度一階彈性波方程和二階彈性波方程可以模擬起伏地表產(chǎn)生的散射波、轉(zhuǎn)換波等各種復(fù)雜現(xiàn)象，起伏地表問題對聲波數(shù)值模擬的動力學(xué)特征會產(chǎn)生很大影響，值得研究探索。同時，作為現(xiàn)階段地震數(shù)據(jù)采集和處理的基本數(shù)據(jù)模型，針對起伏地表問題的聲波方程數(shù)值模擬同樣具有現(xiàn)實意義?？紤]到2.5維數(shù)值模擬具有兼顧模擬精度和計算效率的優(yōu)勢，我們將坐標(biāo)變換法引入到起伏地表2.5維聲波方程數(shù)值模擬中，推導(dǎo)出映射坐標(biāo)系下的2.5維L型聲波近似方程；采用有限差分法，實現(xiàn)起伏地表條件下2.5維聲波方程的數(shù)值模擬；并通過與起伏地表二維聲波方程數(shù)值模擬結(jié)果的對比，來驗證起伏地表2.5維聲波方程有限差分?jǐn)?shù)值模擬方法的有效性。

1 時間域2.5維聲波近似方程

如圖1所示，2.5維問題是三維問題的一種特殊情況。2.5維問題假設(shè)介質(zhì)條件沿走向(圖1中y方向)不變，震源使用三維震源，且與檢波點位于垂直于走向的同一條直線上(圖1中y=0的x-z平面內(nèi))時，就可以把三維數(shù)值模擬問題近似為一個2.5維數(shù)值模擬問題。這樣，2.5維地震波數(shù)值模擬既保留了三維地震波的傳播特點，又具有二維地震波數(shù)值模擬有較高計算效率的優(yōu)點，可以說是三維數(shù)值模擬的一種很好的折中手段。

圖1 2.5維地質(zhì)模型示意圖解

2.5維波動方程數(shù)值模擬的基礎(chǔ)是采用什么形式的2.5維波動方程，關(guān)于該問題Liner[3]假設(shè)在常密度介質(zhì)條件下，2.5維Green函數(shù)是某一空間方向延伸為0(例如y軸方向)的三維Green函數(shù)，即2.5維Green函數(shù)相當(dāng)于是三維Green函數(shù)的一個空間“切片”，表達式為

(1)

(2)

(3)

從左到右各項依次表示平面內(nèi)的二維波動方程、平面外振幅隨時間衰減的散射方程及平面內(nèi)的保幅諧振子方程。方程(3)即為本文要求解的L型2.5維聲波方程。

2 坐標(biāo)變換法

坐標(biāo)變換法的基本思想如圖2所示，通過坐標(biāo)變換將物理空間(x-z)域內(nèi)的起伏地表問題，轉(zhuǎn)化為計算空間(ξ-η)域內(nèi)的水平地表問題[15]。即對起伏地表條件下的L型2.5維聲波近似方程進行坐標(biāo)變換，利用有限差分法求解坐標(biāo)變換后的方程，再將計算結(jié)果映射回起伏地表條件，實現(xiàn)起伏地表條件下2.5維聲波方程有限差分法數(shù)值模擬。

圖2 物理空間(x-z)域內(nèi)的起伏地表(a)及計算空間(ξ-η)域內(nèi)的水平地表(b)

坐標(biāo)變換法的關(guān)鍵是選取合理的映射函數(shù)，我們選取的映射函數(shù)為[12]

(4)

式中：(x,z)為物理空間內(nèi)的起伏地表坐標(biāo)系；(ξ,η)為計算空間內(nèi)的水平地表坐標(biāo)系；z0(ξ)為描述地表起伏的函數(shù)；ηmax為(ξ,η)坐標(biāo)系內(nèi)垂直方向的最大坐標(biāo)；zmax為(x,z)坐標(biāo)系內(nèi)垂直方向的最大坐標(biāo)。

為推導(dǎo)計算空間內(nèi)的2.5維聲波方程，首先對物理空間內(nèi)起伏地表條件下的波場p關(guān)于x(或z)的一、二階偏導(dǎo)作如下變換。

求取波場p對x方向的一、二階偏導(dǎo)：

(5)

(6)

求取波場p對z方向的一、二階偏導(dǎo)：

(7)

(8)

式中：ξx，ηx，ξz和ηz是方程變換過程中的變換系數(shù)。由坐標(biāo)映射式(4)得

(9)

將(6)式、(8)式和(9)式代入(x-z)坐標(biāo)系下的2.5維L型方程(3)，得到映射坐標(biāo)系(ξ-η)下的2.5維L型聲波近似方程：

(10)

此外，由于2.5維聲波近似方程法是通過對二維聲波方程添加校正項，修改振幅、波形和波前等來模擬三維地震波場的，所以本文中采用二維起伏地表聲波方程的單程波動方程邊界條件[15]。由于4個方向的邊界條件及4個角點邊界條件相似，這里只給出左邊界條件及左上角點邊界條件。

左邊界條件：

(11)

左上角邊界條件：

(12)

其中，c1=1.05086，c2=0.75147。

至此，給出了映射坐標(biāo)系下的2.5維L型聲波近似方程(10)，邊界條件(11)和(12)，下面我們將詳細闡述如何采用有限差分格式來離散這些方程和邊界條件。

3 數(shù)值實現(xiàn)

通過以上的坐標(biāo)變換得到了計算空間內(nèi)的2.5維L型聲波近似方程及邊界條件，即水平地表條件下的規(guī)則網(wǎng)格2.5維問題。本文中采用有限差分法，以“差分代替微分”，建立并求解差分方程，來實現(xiàn)數(shù)值模擬。如圖3所示，根據(jù)網(wǎng)格點位置建立差分方程。

圖3 2.5維聲波方程及邊界條件的離散示意圖解

在網(wǎng)格內(nèi)部非邊界上的網(wǎng)格點(如圖3中的點1)，對時間和空間的二階導(dǎo)數(shù)采用二階中心差分格式，對空間的一階導(dǎo)數(shù)采用一階中心差分格式，對時間的一階導(dǎo)數(shù)采用一階向前差分格式即：

(13)

式中：x，z方向的網(wǎng)格間距均為h；Δt為時間步長；l，m是水平和垂直方向的網(wǎng)格指數(shù)；n是時間指數(shù)。然后將(13)式代入方程(10)，得到內(nèi)部網(wǎng)格點的差分方程：

(14)

式中：Ex,z=vx,z·Δt/h；vx,z代表聲波速度，是網(wǎng)格位置x，z的函數(shù)；T=Δt/t，t為旅行時。

在網(wǎng)格邊界及角點處(如圖3中的點2或點3)，同理利用邊界條件和角點條件得到差分方程[15](此處只給出左邊界條件及左上角點邊界條件)。

左邊界條件:

(15)

左上角點條件：

[p(l,m+1,n)-p(l,m+1,n-1)-

p(l,m,n)+p(l,m,n-1)]+E(l,m)·

n-1)-p(l,m,n)+p(l,m,n-1)]-

p(l,m+1,n)-p(l+1,m,n)+p(l,m,n)]-

[p(l,m,n-1)-2p(l,m,n)]

(16)

通過以上的差分方程可以得到計算空間內(nèi)的各網(wǎng)格點的聲壓值。最后，通過映射函數(shù)((4)式)將計算空間內(nèi)(ξ-η)坐標(biāo)系下得到的聲壓值逆映射到物理空間內(nèi)(x-z)坐標(biāo)系下的聲壓值，完成起伏地表條件下2.5維聲波方程的有限差分?jǐn)?shù)值模擬。

4 數(shù)值算例及分析

4.1 均勻介質(zhì)模型

為了驗證起伏地表對直達波的影響，選取一個起伏地表的均勻介質(zhì)模型，如圖4所示。地表函數(shù)為正弦函數(shù)z0=50sin(πx/125+π/2)+50。模型大小為1000m×1100m，網(wǎng)格間距為5m，采樣時間為0.4s，時間步長為0.1ms，震源為主頻30Hz的雷克子波，震源位于模型中間地表以下50m處。

圖4 起伏地表均勻介質(zhì)模型

對起伏地表均勻介質(zhì)模型應(yīng)用2.5維聲波方程有限差分?jǐn)?shù)值模擬給出的直達波示于圖5，可以看出，起伏地表對直達波場具有很大影響，檢波器接收到的直達波場呈現(xiàn)出與地表條件吻合的形狀。

4.2 層狀介質(zhì)模型

為進一步考察起伏地表條件下2.5維反射波場的傳播，建立了一個3層的起伏地表層狀介質(zhì)模型(圖6)，與起伏地表均勻介質(zhì)模型具有相同的地表函數(shù)，3層介質(zhì)的速度參數(shù)分別為1500，2000，2500m/s。除了采樣時間為1s外，其它模擬參數(shù)均與上述均勻介質(zhì)模型的數(shù)值模擬相同。

圖5 起伏地表均勻模型模擬的直達波

應(yīng)用起伏地表條件下2.5維聲波方程有限差分?jǐn)?shù)值模擬方法對起伏地表層狀介質(zhì)模型模擬出的反射波場如圖7所示。從圖7中可以發(fā)現(xiàn)，檢波器在0.4s左右接收到第1個反射界面的反射波，在0.7s左右接收到第2個反射界面的反射波；反射波場的形狀與起伏地表正弦函數(shù)的形狀吻合；存在少量頻散。

圖6 起伏地表層狀介質(zhì)模型

圖7 起伏地表層狀模型模擬的直達波

4.3 與二維模擬結(jié)果的對比

為了驗證2.5維數(shù)值模擬較二維數(shù)值模擬在模擬地震波動力學(xué)特征(振幅)方面的優(yōu)勢，對起伏地表條件下2.5維聲波近似方程數(shù)值模擬與起伏地表條件下二維聲波方程數(shù)值模擬的結(jié)果進行對比分析。選取起伏地表雙層介質(zhì)模型如圖8所示。

起伏地表條件下2.5維聲波方程和二維聲波方程模擬出的反射地震波場如圖9所示。對比可見，2.5維(圖9a)與二維(圖9b)模擬地震記錄的反射波接收時間基本一致。為了更精確的對比，選取兩張記錄中的第50道(橫坐標(biāo)250m處，炮檢距為250m)疊合放大如圖10a所示。分析圖10a可以得到幾點認(rèn)識：①2.5維與二維單道模擬記錄只是在振幅上存在差異，在時間上完全吻合，這就證明了本文中2.5維聲波方程數(shù)值模擬方法的正確性；②二維模擬的反射波振幅(虛線)明顯大于2.5維模擬的反射波振幅(實線)，這是因為二維聲波方程數(shù)值模擬沒有考慮三維球面幾何擴散，使得二維數(shù)值模擬的地震波動力學(xué)特征(振幅)數(shù)值偏大，偏離了實際的地震-地質(zhì)情況；③2.5維數(shù)值模擬得到地震波傳播時間特性與二維模擬結(jié)果具有相同的精度，模擬的地震波振幅特性比二維數(shù)值模擬更接近地震波的空間傳播特性，可見2.5維數(shù)值模擬相較于二維數(shù)值模擬更具優(yōu)勢。

圖8 起伏地表雙層介質(zhì)模型

為進一步驗證起伏地表條件下2.5維聲波方程數(shù)值模擬的衰減特性，選取位于炮點右側(cè)第40道(橫坐標(biāo)200m處，炮檢距為300m)的2.5維與二維模擬單道記錄疊合放大示于圖10b。與圖10a

圖9 起伏地表雙層模型的2.5維(a)和二維(b)模擬地震記錄

圖10 起伏地表雙層模型2.5維(實線)與二維(虛線)模擬結(jié)果的第50道(a)和第40道(b)記錄

對比可以發(fā)現(xiàn)，隨著炮檢距增大，2.5維與二維模擬響應(yīng)的振幅都逐漸減小，但2.5維聲波模擬響應(yīng)振幅的降低幅度明顯大于2維模擬響應(yīng)振幅降低的幅度，這就驗證了2.5維聲波方程數(shù)值模擬在時間和空間上的衰減特性。

5 結(jié)論與認(rèn)識

基于Liner提出的2.5維L型聲波方程，采用坐標(biāo)變換的方法，實現(xiàn)了起伏地表條件下2.5維聲波近似方程的有限差分法數(shù)值模擬。通過理論模型試算，將該方法與二維起伏地表條件下聲波方程數(shù)值模擬方法作對比，發(fā)現(xiàn)兩者的模擬結(jié)果在地震波運動學(xué)特征上基本一致，而在地震波動力學(xué)特征上存在差異，證實了本文給出的起伏地表2.5維聲波方程數(shù)值模擬方法相較于二維聲波方程數(shù)值模擬方法在模擬地震波動力學(xué)特征上更有優(yōu)勢。同時，通過不同炮檢距模擬記錄的振幅對比，驗證了2.5維聲波方程數(shù)值模擬在時間和空間上的衰減特性。研究結(jié)果表明，2.5維聲波方程數(shù)值模擬可以很好地處理起伏地表問題，在模擬精度和計算效率方面都是三維數(shù)值模擬的一種很好的折中方法。

參 考 文 獻

[1] 孫建國.2.5維地震波數(shù)值模擬評述:聲波模型[J].地球物理學(xué)進展,2009,24(1):20-34

Sun J G.A review of 2.5 dimensional seismic modeling methods:acoustic case[J].Progress in Geophysics,2009,24(1):20-34

[2] Novais A,Santos L T.2.5-D finite-difference solution of the acoustic wave equation[J].Geophysical Prospecting,2005,53(4):523-531

[3] Liner C L.Theory of a 2.5-D acoustic wave equation for constant density media[J].Geophysics,1991,56(12):2114-2117

[4] Williamson P R，Pratt P R.A critical review of acoustic wave modeling procedures in 2.5 dimensions[J].Geophysics,1995,60(2):591-595

[5] Narayan J P.2.5-D numerical simulation of acoustic wave propagation[J].Pure and Applied Geophysics,1998,151:47-81

[6] 李芳，孫建國，王雪秋,等.2.5維聲波數(shù)值模擬[J].吉林大學(xué)學(xué)報(地球科學(xué)版),2006,36(增刊):177-181

Li F,Sun J G,Wang X Q.Numerical simulation of 2.5D acoustic waves[J].Journal of Jilin University(Earth Science Edition),2006,36(S1):177-181

[7] 孫建國.復(fù)雜地表條件下地球物理場數(shù)值模擬方法評述[J].世界地質(zhì),2007,26(3):345-362

Sun J G.Methods for numerical modeling of geophysical fields under complex topographical conditions:a critical review[J].Global Geology,2007,26(3):345-362

[8] 楊海生.起伏地表波動方程波場延拓起始面研究[J].石油物探,2009,48(1)：66-71

Yang H S.Study on staring datum of wave equation wave field continuation for rolling surfaces[J].Geophysical Prospecting for Petroleum,2009,48(1):66-71

[9] 孫章慶,孫建國.2.5維起伏地表條件下坐標(biāo)變換法直流電場數(shù)值模擬[J].吉林大學(xué)學(xué)報(地球科學(xué)版),2010,40(2):425-431

Sun Z Q,Sun J G.2.5-D DC electric field numerical modeling including surface topography based on coordinate transformation method[J].Journal of Jilin University(Earth Science Edition),2010,40(2):425-431

[10] 郭宏偉,王尚旭，孫文博.基于非均質(zhì)體的波動方程有限元正演模擬[J].石油物探，2012,51(4):319-326

Guo H W,Wang S X,Sun W B.Wave equation modeling by finite-element method for heterogeneous body[J].Geophysical Prospecting for Petroleum,2012,51(4):319-216

[11] Tessmer E,Kosloff D.3-D elastic modeling with surface topography by a chebychev spectral method[J].Geophysics,1994,59(3):464-473

[12] Hestholm S,Ruud B.2-D finite-difference elastic wave modeling including surface topography[J].Geophysical Prospecting,1994,42(5):371-390

[13] Hestholm S,Ruud B.3-D finite-difference elastic wave modeling including surface topography[J].Geophysics,1998,63(2):613-622

[14] 唐文，王尚旭，袁三一.起伏地表二階彈性波方程差分策略穩(wěn)定性分析[J].石油物探，2013,52(5)：457-463

Tang W,Wang S X,Yuan S Y.Stability analysis of differential strategy for rugged by second-order elastic wave equation based on coordinate transformation[J].Geophysical Prospecting for Petroleum,2013,52(5):457-463

[15] 王祥春，劉學(xué)偉.起伏地表二維聲波方程地震波場模擬與分析[J].石油地球物理勘探,2007,42(3):268-276

Wang X C,Liu X W.Simulation an analysis of seismic wavefield in relief surface by 2D acoustic wave equation[J].Oil Geophysical Prospecting,2007,42(3):268-276