聲波方程數(shù)值模擬矩形網(wǎng)格有限差分系數(shù)確定法

梁文全 王彥飛 楊長春

（①龍巖學(xué)院資源工程學(xué)院，福建龍巖364000；②中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，北京100029）

聲波方程數(shù)值模擬矩形網(wǎng)格有限差分系數(shù)確定法

梁文全①王彥飛*②楊長春②

（①龍巖學(xué)院資源工程學(xué)院，福建龍巖364000；②中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，北京100029）

壓制數(shù)值頻散是有限差分方法的關(guān)鍵問題之一。目前壓制數(shù)值頻散的方法大多假設(shè)不同方向空間偏導(dǎo)數(shù)的空間步長相同，導(dǎo)致算法精度低，計(jì)算效率低。為此，提出使用線性方法壓制聲波方程矩形網(wǎng)格有限差分算子的數(shù)值頻散，并進(jìn)行了穩(wěn)定性分析、頻散分析和數(shù)值模擬。通過頻散分析和數(shù)值模擬，驗(yàn)證了本文方法能夠有效壓制矩形網(wǎng)格有限差分?jǐn)?shù)值頻散，相較于泰勒展開方法和最小二乘方法，線性方法計(jì)算有限差分系數(shù)的效率更高，可以替代傳統(tǒng)的正方形有限差分網(wǎng)格和相應(yīng)的系數(shù)用于聲波方程數(shù)值延拓。

聲波模擬 時(shí)間—空間域 有限差分格式 矩形網(wǎng)格 數(shù)值頻散

1 引言

地震波數(shù)值模擬延拓是地震學(xué)和勘探地震學(xué)的重要基礎(chǔ)，在理論和應(yīng)用中都發(fā)揮了重要的作用［1］。地震波數(shù)值模擬延拓的方法包括有限差分法、偽譜法、傳統(tǒng)有限元法、譜元法、有限體積法、間斷Galerkin法等［2，3］。有限差分方法因?yàn)橛?jì)算效率高、所需內(nèi)存小、實(shí)現(xiàn)簡單而廣泛應(yīng)用于地震正演研究，同時(shí)是疊前逆時(shí)深度偏移成像、全波形反演、速度分析技術(shù)迅速發(fā)展的基礎(chǔ)［4，5］。

數(shù)值頻散是對波動方程的時(shí)間和空間偏導(dǎo)數(shù)的離散化造成的，使相速度變成了網(wǎng)格間距的函數(shù)，降低了模擬的精度。在地震波數(shù)值模擬中，一般使用正方形（正方體）差分網(wǎng)格進(jìn)行數(shù)值模擬［6，7］。在實(shí)際地震資料處理中，需要使用矩形網(wǎng)格有限差分格式。Chen［8］基于平均導(dǎo)數(shù)法提出一種新的9點(diǎn)有限差分算子，可以解決頻率—空間域數(shù)值模擬中不同方向空間采樣間距不同的問題，并提高了頻率域數(shù)值模擬的精度。Tang等［9］使用平均導(dǎo)數(shù)法優(yōu)化了17點(diǎn)矩形網(wǎng)格有限差分格式在頻率—空間域聲波方程中的應(yīng)用。Chu等［10］討論了矩形網(wǎng)格（立方體）隱格式有限差分地震波數(shù)值模擬。但是鮮見對時(shí)間域矩形網(wǎng)格（立方體）有限差分算子系數(shù)如何確定的討論，并且一般采用泰勒展開法確定有限差分系數(shù)［11］。為此本文使用效率較高的線性方法確定聲波方程矩形網(wǎng)格的有限差分算子系數(shù)，以期提高地震波模擬的精度和效率。

通常在空間域確定波動方程空間偏導(dǎo)數(shù)的有限差分算子系數(shù)。但是，地震波在時(shí)間域和空間域同時(shí)傳播，因此在時(shí)間—空間域確定空間偏導(dǎo)數(shù)的有限差分算子系數(shù)能同時(shí)抑制時(shí)間頻散和空間頻散。Liu等［12］通過平面波理論和泰勒展開方法在時(shí)間—空間域確定有限差分算子系數(shù)，在壓制空間頻散的同時(shí)壓制了時(shí)間頻散。Yang等［13］提出和改進(jìn)了近似解析離散化方法以提高數(shù)值模擬的精度。Zhang等［14］提出使用模擬退火方法確定有限差分系數(shù)并指出選擇合適的頻散誤差范圍的重要性。Liu等［15］提出新的差分模板并用優(yōu)化方法確定相應(yīng)的系數(shù)，提高了數(shù)值模擬的精度。文中運(yùn)用線性方法確定了聲波方程矩形網(wǎng)格有限差分系數(shù)，并進(jìn)行穩(wěn)定性分析、頻散分析和數(shù)值模擬，驗(yàn)證了本文方法的有效性。

2 有限差分算子確定方法

對于二維聲波方程

在（x，z）兩個(gè)方向上的二階空間偏導(dǎo)數(shù)，使用不同的空間步長計(jì)算，可以得到

對于時(shí)間偏導(dǎo)數(shù)，采用二階有限差分近似，得到

將式（2）和式（3）代入式（1），可得

對式（4）進(jìn)行傅里葉變換，可得

式中：k x=kcosθ，k z=ksinθ；k為波數(shù)；θ是平面波的傳播角。

根據(jù)震源頻率、波傳播速度和網(wǎng)格間距確定需要滿足頻散關(guān)系的波數(shù)占總波數(shù)的比例［16］，即

式中：f是震源最高頻率；v表示地震波速度；h=max（h z，h x），假設(shè)h z＜h x；ku表示有限差分算子需要優(yōu)化的波數(shù)；ktotal表示總的波數(shù)。

在式（6）確定的波數(shù)范圍內(nèi)，假設(shè)M+1個(gè)均勻分布的波數(shù)點(diǎn)滿足式（5）所示的時(shí)間—空間域頻散關(guān)系，則得到線性方程組

式中：ki，l（i=1，2，…，M+1；l=x，z）均勻分布于0和之間，而R通過式（6）確定；由式（7）可以確定聲波方程矩形網(wǎng)格有限差分算子系數(shù)。

由式（5）可以得到一般來說，隨著波數(shù)（頻率）的增加，頻散誤差增大，令

將式（9）代入式（8），得到聲波方程矩形網(wǎng)格有限差分算子的穩(wěn)定性條件，即

圖1是聲波方程矩形網(wǎng)格有限差分算子的穩(wěn)定性條件。從圖1可以看出，在hz確定時(shí)，h x越大，穩(wěn)定的時(shí)間范圍越大。以Marmousi模型為例，其最高波速為4700m／s，假如hz=h x=4m，則其最大時(shí)間步長≤0.46ms；假如h z=4m，h x=12.5m，則其時(shí)間步長可以達(dá)到0.6ms。

3 頻散分析

二維聲波方程數(shù)值頻散定義為［17］

式中：vphase-grid表示離散相速度；vreal表示真實(shí)相速度；1-δ（θ）表示數(shù)值頻散，1-δ（θ）的絕對值越小，數(shù)值頻散越小。

圖1 矩形網(wǎng)格有限差分的穩(wěn)定性條件

有限差分方法在一個(gè)空間網(wǎng)格上傳播的時(shí)間誤差為［17］

圖2為泰勒展開方法和本文方法確定的矩形網(wǎng)格有限差分系數(shù)的數(shù)值頻散，以M=30，τ=0.5ms泰勒展開方法作為近似解析解。由圖2可以看出，在kh∈［0，2.5］范圍內(nèi)，本文方法的數(shù)值頻散（圖2b、圖2e）要小于泰勒展開方法的數(shù)值頻散（圖1a、圖1d），而接近于近似解析解數(shù)值頻散（圖1c、圖1f）。

圖2 不同方法矩形網(wǎng)格有限差分系數(shù)數(shù)值頻散誤差

4 數(shù)值模擬

4.1 均勻介質(zhì)模型

首先考慮均勻介質(zhì)模型，對于所有有限差分算子v=2000m／s，h z=10m，h x=20m，τ=1ms，M=7。模型網(wǎng)格數(shù)為350×700。震源在模型的中心位置，震源是高斯函數(shù)的二階導(dǎo)數(shù)［18］，即

式中：f0=55Hz；是主頻頻率；c是常數(shù)。

圖3為不同有限差分算法對應(yīng)的1500ms時(shí)的波場快照及t=1500ms時(shí)的橫向地震記錄。使用τ=0.5ms，M=30的傳統(tǒng)泰勒展開方法有限差分算子系數(shù)模擬結(jié)果作為近似解析解。由圖可見，對于各種方法確定的有限差分算子系數(shù)，z方向的網(wǎng)格間距小，數(shù)值頻散都不明顯；x方向網(wǎng)格間距較大，數(shù)值頻散明顯（圖3a、圖3b）。由圖3d、圖3e以及泰勒方法和本文方法與近似解析解方法的殘差（圖3f）可見，較泰勒展開方法，本文方法壓制了數(shù)值頻散，更適用于確定矩形網(wǎng)格有限差分算子系數(shù)。對比圖3d和圖3f可知直達(dá)波的誤差達(dá)到了10%左右，主要是由于本文方法時(shí)間步長較大引起的。如果將圖3b和圖3c的時(shí)間步長各縮小為原來時(shí)間步長的一半，則直達(dá)波的誤差可以限制在2.5%以內(nèi)。

圖3 不同有限差分方法對應(yīng)的1500ms處波場快照（上）及t=1500ms時(shí)的橫向地震記錄（下）

4.2 雙層介質(zhì)模型

不同的采樣間隔比dx／dz顯著提高了地震波數(shù)值模擬的效率，但是對地震記錄的走時(shí)有較為顯著的影響。例如一個(gè)雙層介質(zhì)模型，從0～1000m深度的波速是2000m／s，大于1000m深度的波速是2500m／s（圖4a）。深度z方向分別采用5、10、20m的網(wǎng)格劃分，其界面深度分別在1002.5m（200.5×5m），1005m（100.5×10m），1010m（50.5×20m圖4b），因此造成了地震記錄走時(shí)的不同。雙層介質(zhì)模型（圖4a）得到的地震記錄如圖4c（震源在地表）所示。三種采樣間隔比得到的單道地震記錄如圖4d（震源處接收）所示。由圖4e和圖4f可知，當(dāng)深度z方向分別采用5、10、20m的網(wǎng)格劃分時(shí)，其反射波第一個(gè)波峰時(shí)間分別是1076、1080、1086ms。dz=5m與dz=10m時(shí)網(wǎng)格劃分走時(shí)相差4ms（理論上是2.5ms），dz=10m與dz=20m網(wǎng)格劃分走時(shí)相差6ms（理論上是5ms）。因此網(wǎng)格間隔比的設(shè)計(jì)要根據(jù)實(shí)際情況綜合考慮計(jì)算效率和精度。

圖4 雙層介質(zhì)模型地震波數(shù)值模擬結(jié)果

4.3 Marmousi模型

圖5為Marmousi速度模型。震源函數(shù)采用式（13），其中震源頻率f0=70Hz。對于所有的有限差分算子hx=12.5m，hz=4m。使用單程波和雙程波混合邊界條件壓制邊界反射［19］。圖6為不同有限差分算子對應(yīng)的地震記錄及t=600ms時(shí)的橫向地震記錄。其中傳統(tǒng)泰勒展開方法時(shí)間步長τ=0.6ms，M=7；本文方法時(shí)間步長τ=0.6ms，M=7；以τ=0.3ms，M=30的泰勒展開法數(shù)值模擬作為解析解。從圖6a～圖6c可以看出，本文方法較傳統(tǒng)泰勒展開方法在相同算子長度和時(shí)間步長時(shí)有效壓制了空間頻散（圖中框內(nèi)），效果更接近于近似解析解。圖6d～圖6f比較了600ms時(shí)的橫向地震記錄（z=0），同樣可以看出，本文方法有效地壓制了數(shù)值頻散。

圖5 Marmousi速度模型

圖6 不同有限差分方法對應(yīng)的地震記錄（上）及t=600ms時(shí)的橫向地震記錄（下）

5 結(jié)論

（1）在相同算子長度和時(shí)間步長時(shí)本文方法較傳統(tǒng)泰勒展開方法有效地壓制了空間頻散，效果更接近于近似解析解。

（2）不同的網(wǎng)格間隔比dx／dz顯著提高了地震波數(shù)值模擬的效率，但是對地震記錄的走時(shí)有較顯著的影響，因此網(wǎng)格間隔比的設(shè)計(jì)要根據(jù)實(shí)際情況綜合考慮計(jì)算效率和精度。

（3）通過頻散分析和數(shù)值模擬，驗(yàn)證了本文方法能夠有效壓制矩形網(wǎng)格有限差分?jǐn)?shù)值頻散，相較于泰勒展開方法和最小二乘方法，線性方法計(jì)算有限差分系數(shù)的效率更高。因此本文方法可以替代傳統(tǒng)的正方形有限差分網(wǎng)格和相應(yīng)的系數(shù)用于聲波方程數(shù)值延拓。

［1］ Dablain M A.The application of high-order differencing to the scalar wave equation.Geophysics，1986，51（1）：54-66.

［2］ De Basabe J D and Sen M K.Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations.Geophysics，2007，72（6）：T81-T95.

［3］ De Basabe J D and Sen M K.A comparison of finitedifference and spectral-element methods for elastic wave propagation in media with a fluid-solid interface.Geophysical Journal International，2015，200（1）：278-298.

［4］ Liu H W，Li B，Liu H et al.The issues of prestack reverse time migration and solutions with graphic processing unit implementation.Geophysical Prospecting，2012，60（5）：906-918.

［5］ 劉洋.波動方程時(shí)空域有限差分?jǐn)?shù)值解及吸收邊界條件研究進(jìn)展.石油地球物理勘探，2014，49（1）：35-46.Liu Yang.Research progress on time-space domain finite-difference numerical solution and absorption boundary conditions of wave equation.OGP，2014，49（1）：35-46.

［6］ 梁文全，楊長春，王彥飛等.用于聲波方程數(shù)值模擬的時(shí)間—空間域有限差分系數(shù)確定新方法.地球物理學(xué)報(bào)，2013，56（10）：3497-3506.Liang Wenquan，Yang Changchun，Wang Yanfei et al.Acoustic wave equation modeling with new time-space domain finite difference operators.Chinese Journal of Geophysics，2013，56（10）：3497-3506.

［7］ 張志禹，譚顯波，黃璐瑤等.抗頻散有限差分波動方程數(shù)值模擬及逆時(shí)偏移.石油地球物理勘探，2014，49（6）：1115-1121.Zhang Zhiyu，Tan Xianbo，Huang Luyao et al.Antidispersion finite difference simulation and reversetime migration for wave equations.OGP，2014，49（6）：1115-1121.

［8］ Chen J B.An average-derivative optimal scheme for frequency-domain scalar wave equation.Geophysics，2012，77（6）：T201-T210.

［9］ Tang X，Liu H，Zhang H et al.An adaptable 17-point scheme for high-accuracy frequency-domain acoustic wave modeling in 2D constant density media.Geophysics，2015，80（6）：T211-T221.

［10］ Chu C and Stoffa P L.Implicit finite-difference simulations of seismic wave propagation.Geophysics，2012，77（2）：T57-T67.

［11］ Chen H，Zhou H and Li Y.Application of unsplit convolutional perfectly matched layer for scalar arbitrarily wide-angle wave equation.Geophysics，2014，79（6）：T313-T321.

［12］ Liu Y，Sen M K.Acoustic VTI modeling with a timespace domain dispersion-relation-based finite-difference scheme.Geophysics，2010，75（3）：A11-A17.

［13］ Yang D，Lu M，Wu R et al.An optimal nearly analytic discrete method for 2D acoustic and elastic wave equations.Bulletin of the Seismological Society of A-merica，2004，94（5）：1982-1992.

［14］ Zhang J H，Yao Z X.Optimized finite-difference operator for broadband seismic wave modeling.Geophysics，2012，78（1）：A13-A18.

［15］ Liu H F，Dai H X，Niu F L et al.An explicit time evolution method for acoustic wave propagation.Geophysics，2014，79（3）：T117-T124.

［16］ Wang Y F，Liang W Q，Nashed Z et al.Seismic modeling by optimizing regularized staggered-grid finitedifference operators using a time-space-domain dispersion-relationship-preserving method.Geophysics，2014，79（5）：T277-T285.

［17］ Liu Y，Sen M K.A new time-space domain high-order finite-difference method for the acoustic wave equation.Journal of Computational Physics，2009，228（23）：8779-8806.

［18］ 王彥飛，斯捷潘諾娃I E，提塔連科V N等.地球物理數(shù)值反演問題.北京：高等教育出版社，2011.

［19］ Liu Y，Sen M K.A hybrid scheme for absorbing edge reflections in numerical modeling of wave propagation.Geophysics，2010，75（2）：A1-A6.

P631

A

10.13810／j.cnki.issn.1000-7210.2017.01.009

梁文全，王彥飛，楊長春.聲波方程數(shù)值模擬矩形網(wǎng)格有限差分系數(shù)確定法.石油地球物理勘探，2017，52（1）：56-62.

1000-7210（2017）01-0056-07

*北京市朝陽區(qū)北士城西路19號，100029。Email：yfwang＠m(xù)ail.iggcas.ac.cn

本文于2015年11月3日收到，最終修改稿于2016年11月1日收到。

本項(xiàng)研究受國家自然科學(xué)基金項(xiàng)目（41325016，91630202）及福建省自然科學(xué)基金（2016J05104）和龍巖學(xué)院博士科研啟動基金（LB2014010）資助。

（本文編輯：金文昱）

梁文全 講師，1983年生；2006年獲南陽師范學(xué)院計(jì)算機(jī)系學(xué)士學(xué)位；2009年獲得中國地質(zhì)大學(xué)（北京）地理信息系統(tǒng)專業(yè)碩士學(xué)位；2012年獲得中國科學(xué)院地質(zhì)與地球物理研究所固體地球物理專業(yè)博士學(xué)位；目前在龍巖學(xué)院從事地震波數(shù)值模擬、偏移及反演研究及教學(xué)工作。