声波方程隐式高阶有限差分数值模拟_智敏.pdf

第 44 卷 第 2 期 煤田地质与勘探 Vol. 44 No.2 2016 年 4 月 COAL GEOLOGY numerical simulation; implicit difference; numerical frequency dispersion 声波方程数值模拟已广泛应用于地震正演、 逆时 偏移、全波形反演等各个方面，如何进一步提高数值 模拟精度、 减小频散和计算量一直是有限差分算法研 究的核心问题[1-2]。目前，相关文献给出的空间二阶 导数计算大多采用显式差分格式， 这种方法通常具有 较强的频散现象[3-4]。而为了有效地降低数值频散， 通常采用减小网格尺寸、提高差分阶数的方法，这种 办法虽然提高了数值模拟的精度， 但也大幅增加了计 算量和内存空间[5]。 Boris 和 Book 将求解流体力学方 程的通量校正传输方法FCT应用到波动方程数值模 拟中， 该方法通过差分运算降低高波数端能量压制了 频散噪声，但不利于提高模拟结果的分辨率[3-7]。 采用隐式差分格式求解波场的空间二阶导数较 显式法具有更高的精度。前人研究表明隐式差分格 式的2N2阶差分算子的精度与显式差分算子4N2 阶的计算精度相当[8]。因此，笔者研究了空间二阶导 数的隐式求解方法，理论分析及模拟结果均表明该 方法较显式求解方法具有更小的空间频散及更高的 模拟精度。 1 方法原理 1.1 基本方程 二维声波方程如式1所示 222 2222 1ppp xyvt ∂∂∂ ∂∂∂ 1 式中 p为地震波场；v为纵波速度。 如果只考虑 x 方 向，根据波场传播理论可有（）（ ） i e k x p xxp x Δ Δ， 因此根据欧拉公式可推出下式 （ ）（）（）（ ）（） （ ） 22 222 24sin/2p xp xxp xxp xk x p x xxx δ δ Δ-Δ--Δ ΔΔ 2 式中 2 2 p x x δ δ 为差分算子；k为波数；xΔ为空间采样 间隔。 为了提高差分精度，可以取 224 2 224 x ppp qb x xxx δδ δδ ∂ ≈- Δ ∂ 3 式中 b为一待定常数。由于 1 1,1 1 x x x ≈ -＜＜ ，所以 ChaoXing 第2期 智敏 声波方程隐式高阶有限差分数值模拟 107 式3可以转化为式4 [6] 2 22 2 2 222 2 2 1 1 x p x qb xp xx b x x δ δδ δ δδδ δ ■■ ≈- Δ≈ ■■ ■■ ■■ Δ 4 由式4可得 22 2 22 x x q qb xp xx δδ δδ Δ 5 将式5中的差分算子利用式2展开， 可得式6 （） （）（ ）（） 2 2 12 xxx bqxxb qxbqxxp x δ δ -Δ-Δ 6 如果将差分算子写成形如式7的差分格式， （ ）（）（） 2 0 22 1 1 M m m p c p xcp xm xp xm x xx δ δ ■■ ■■ Δ--Δ■■ ■■ ■■ Δ■■ ■■ Σ 7 则在式4可以写成 （） （） （） 0 2 1 2 2cos 14 sin2 M m m ccmk x k x bk x Δ -Δ≈ -Δ Σ 8 由此可以确定一组系数， 让差分算子与微分算子 高阶逼近。 1.2 差分系数的确定 令Xk x Δ，则式8化简可得 （）（） 2 0 1 122 cos2cos M m m XbbXccmX --≈■■ ■■Σ 9 将式9两边都对余弦函数进行taylor展开， 可得 如下关系 （） （） （） （） （） 2 2 2 0 101 11 122 22 nn n M n m nnm X XbccmX nn ∞∞ ■■ -- ■■-≈ ■■ ■■ ΣΣΣ 10 移项得 （） （） （） （） 22 0 11 1 22 21 21 11 20 221 MM mm mm nn M nn m nm ccc mX bc mX nn - ∞ ■■■■ - ■■■■ ■■■■ ■■■■ -- ■■■■ ≈ ■■■■ -■■ ■■■■■■ ■■■■ ΣΣ ΣΣ 11 若要式11对任意X成立，则X的前M次方系 数均为0，即 0 1 20 M m m cc Σ 12 2 1 10 M m m c m - Σ 13 （）（） 2 1 11 0 221 M n m m bc m nn ■■ -■■ -■■ ■■ ■■■■ Σ 14 式13式14可以写成形如式15的矩阵形式 （）（） （）（） 0000 2222 2222 2222 1230 1236 1232221 1232221 ■■ ■■ -■■ ■■ ■■ ■■ - ■■ ■■ ■■ - ■■ ■■ ⋮⋮⋮⋮⋮⋮ ⋮⋮⋮⋮⋮⋮ nnnn MMMM M M Mnn MMM 2 1 2 2 2 3 2 1 1 20 0 3 0 0 ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■ ■ ■■ ⋮ ⋮ m c c c M c b 15 通过求解线性方程组式15可以得到差分系数 （）1,2,, m cmM及待定系数b，再结合式12即可 得到差分系数 0 c。得到差分系数后，就可以应用追 赶法隐式求解式6得到波场在X方向上的空间二阶 导数。 用同样的方法可以求得波场在Z方向上的空间 二阶导数 z q，再利用式16对波场进行更新，即可得 到新时刻的波场值。 （） 2 11 2 2 kkk xz t PqqPP v - Δ - 16 式中 tΔ为时间上的采样间隔。 从上述分析可知所确定的差分系数可以保证 Xk x Δ的前2M次方的系数均为零。因此，该方法 的差分精度可以达到2M2阶。 表1是不同阶数的差 分系数表。 表 1 隐式差分格式差分系数表 Table 1 Coefficients of implicit difference M b c0 c1 c2 c3 c4 c5 c6 1 0.083 333 –2 1 2 0.133 333 –1.7 0.8 0.05 3 0.160 714 –1.490 08 0.656 250 0.091 071 –0.002 282 4 0.177 778 –1.345 99 0.556 049 0.121 975 –0.005 2200.000 190 5 0.189 394 –1.242 67 0.483 481 0.145 503 –0.008 1540.000 526 –0.000 021 6 0.197 802 –1.165 42 0.428 838 0.163 854 –0.010 8660.000 946 –0.000 068 0.000 003 ChaoXing 108 煤田地质与勘探 第44卷 1.3 边界条件 数值模拟是在有限的空间内进行，不可避免地 会产生边界反射，需对边界进行处理，减弱边界反 射，提高模拟精度[9]。这里采用的是海绵吸收边界 与二阶Clayton-Engquist-majda吸收边界相结合的方 法进行边界处理。 1.3.1 海绵吸收边界条件 在模型边界设置一定厚度的吸收边界层来消除 回绕和反射能量；在每一步的计算过程中，吸收区 域的每一个网格节点的振幅均乘上一个递减的高斯 函数[9-11] （ ）（） 22 exp,0,1,2,,G iaiNbiNb-- 17 式中 Nb为吸收边界的宽度网格节点数；i为计算 节点距边界的网格点数；a为衰减系数。一般 20Nb ， 0.015a 。 1.3.2 二阶Clayton-Engquist-majda吸收边界条件 二阶Clayton-Engquist-majda边界条件的微分 表达式如下 222 222 12 0 PPP v t nvtn ∂∂∂ ∂ ∂∂∂ →→ 18 式中 n → 为模型边界的外法线方向。其具体差分格式为[12] （） （） （） 12211 0,0,1,2,0,1, 12211 ,,1,2,,1, 122 ,0,0,1, 2221212 2221212 2221 kkkkkk jxxjxxjxjxjxj kkkkkk Nx jxxNx jxxNxjxNxjxNx jxNxj kkk izzizzizi PAAPAAPA PAPA P PAAPAAPA PAPA P PAAPAAPA P -- -- --- ----- ----- --- （） 11 2,0,1 12211 ,,,1,2,,1 212 2221212 kkk zizi kkkkkk i Nzzzi Nzzzj Nzzi Nzzi Nzz i Nz APA P PAAPAAPA PAPA P -- -- --- ■ ■ ■ ■ ■ --■ ■ ----- ■ ■ 19 式中 / x Av tx ΔΔ；/ z Av tz ΔΔ。 2 基本流程 声波方程高阶隐式有限差分数值模拟的基本流 程如图1所示。 图 1 声波方程数值模拟流程 Fig.1 Flowchart of numerical simulation of acoustic wave equation 3 频散分析 利用式20来衡量隐式有限差分格式的频散特 征，（）f k xΔ越接近1，数值频散越小，反之数值频 散越大。 （） （） （）（）（） 0 1 2 2 2cos 14 sin2 M m m ccmk x f k x bk xk x Δ Δ- -ΔΔ Σ 20 图2是显式差分格式与隐式差分格式的数值频散 曲线对比图。可以看到随着差分阶数的增加，显式 与隐式差分算子的数值频散均在减小；当差分阶数相 同时， 隐式差分格式的频散要明显小于显式差分格式。 图2中的实线隐式与虚线显式所示， 隐式差分格式 的4阶差分精度可以达到显式差分8阶差分精度相当 的频散特征， 而隐式差分格式的12阶频散甚至要比显 式差分格式的22阶要小。 因此， 隐式差分格式要比显 式差分格式具有更小的数值频散。 图 2 显式与隐式差分格式频散曲线对比 Fig.2 Dispersion curve comparison of explicit and implicit difference scheme ChaoXing 第2期 智敏 声波方程隐式高阶有限差分数值模拟 109 4 算 例 为了验证方法的有效性， 分别利用隐式差分格式 和显式差分格式对均匀介质模型及Marmousi模型进 行了数值模拟。 4.1 均匀介质模型 均匀介质模型大小为2 200 m2 200 m，空间步 长为5 m，介质速度为2 000 m/s；采用中心频率为 60 Hz的Ricker子波作为震源，时间步长为0.5 ms， 震源位于模型中心1 100 m，1 100 m。分别采用显 式差分格式和隐式差分格式进行数值模拟， 时间均为 4阶精度，空间均为8阶精度。 图3为显式差分格式与隐式差分格式的波场对比 图。可以看出，在显式差分格式模拟的波场上的首波 后面出现了明显的尾波，这是由空间频散所致。而应 用隐式差分格式，首波后的尾波得到了很好的压制， 说明隐式方法较显式方法具有更高的计算精度。 图 3 显式与隐式差分方案波场对比 Fig.3 Comparison of wave field of the explicit and implicit difference scheme 图4是单个检波器接收到的声波波形， 从波形可 以看出，应用显式差分格式的波形频散较为严重，首 波后面的尾波能量较强，而应用隐式差分格式后，尾 波的能量明显减弱， 这说明隐式差分方法有效地压制 了频散，提高了数值模拟的精度。 图 4 显式与隐式波形曲线图对比 Fig.4 Comparison of wave curve of explicit and implicit difference scheme 图5是边界条件处理前后，t700 ms的波场快照 对比。可以看到，应用边界条件后，来自边界的反射 波得到很好地压制， 说明本文所述的边界处理方法具 好良好的效果。 图6是显式差分格式与隐式差分格式 的计算机时随差分阶数的变化曲线。 从该曲线可以看 出，随着差分阶数的增加，显式格式及隐式格式的计 算机时基本都呈线性增加。 图 5 边界条件处理前后波场快照对比 Fig.5 Comparison of wave field before and after snapshot of boundary conditions ChaoXing 110 煤田地质与勘探 第44卷 从前文分析得知，4阶隐式格式的频散特征与 8阶显式格式相当，12阶隐式差分格式的频散比22 阶显式差分格式小。8阶显式格式的计算机时为 55.462 s，4阶隐式格式的计算机时为58.256 s，22 阶显式格式的计算机时为143.376 s，12阶隐式格 式的计算机时为124.014 s，因此在频散特征相同 时，两种差分方案的低阶差分格式计算效率基本相 当，当差分阶数较高时，隐式格式的计算效率要高 于显式格式。 图 6 显式差分与隐式差分格式计算效率 Fig.6 Calculation efficiency of explicit and implicit difference schemes 4.2 Marmo usi 模型 图7所示为Marmousi速度模型。该模型横向长7 900 m，纵向长3 500 m，空间网格大小5 m，介质速度 1 5004 800 m/s； 采用主频为60 Hz的Ricker子波作为 震源进行模拟， 时间步长为0.5 ms， 记录长度3 000 ms， 震源位于4 250 m，40 m处，检波点置于Z40 m处， 边界采用海绵吸收与二阶Clayton-Engquist-majda吸收 边界相结合的方法进行处理， 其中海绵吸收边界的厚度 为20个节点。分别采用显式和隐式差分格式对该模型 进行数值模拟， 为了对比方便， 两种差分格式均采用时 间4阶、空间8阶的差分精度。 图 7 Marmousi 模型 Fig.7 Marmousi model 图8为840 ms时两种差分方法的波场快照。可 以看到两种差分方案在同样网格参数及差分阶数 时，隐式差分格式比显式差分格式的频散小，波场精 度高。 图 8 显式与隐式差分方案波场对比 Fig.8 Comparison of wave field of explicit and implicit difference schemes 图 9 显式与隐式差分方案单炮记录对比 Fig.9 Comparison of single shot record of explicit and implicit difference scheme ChaoXing 第2期 智敏 声波方程隐式高阶有限差分数值模拟 111 图9是隐式和显式差分的单炮记录对比， 可见显 式差分格式模拟结果在700850道的800 ms左右可 以看到明显的反射波拖尾现象图中椭圆区域，在隐 式差分格式模拟记录上，这种现象得到了压制，模拟 地震记录的信噪比和分辨率均得到了大幅改善。 图 10 显式与隐式差分方案自激自收剖面对比 Fig.10 Comparison of zero offset profile of explicit and im- plicit difference scheme 图10是隐式和显式差分的自激自收剖面对比局 部放大图，和单炮记录对比相似，在显式差分格式 模拟结果上可以看到较明显的频散现象图中椭圆区 域，而隐式差分格式模拟结果上，数值频散得到很 好的压制。 5 结 语 应用隐式差分格式求取波场的二阶空间导数较 传统的显式求解式具有更高的精度， 可以有效地压制 空间频散，提高数值模拟的信噪比和分辨率。 该方法可以较好地模拟叠前单炮和叠后水平叠 加剖面， 并方便地推广到三维声波方程数值模拟和逆 时深度偏移成像等方面。 参考文献 [1] 陈可洋．标量声波波动方程高阶交错网格有限差分法[J]．中国 海上油气，2009，214232-236． CHEN Keyang. High-order staggered-grid finite difference scheme for scalar acoustic wave equation[J]. China Offshore Oil And Gas，2009，214232-236. [2] 曹丹平．声波数值模拟中的多核并行方法研究[J]．计算机工程 与应用，2012，48369-13． CAO Danping. Multicore parallel programming of the acoustic equation numerical modeling[J]. Computer Engineering and Ap- plications，2012，48369-13. [3] 朱生旺，魏修成．波动方程数值模拟的隐式差分法[J]．石油物 探，2006，452151-156． ZHU Shengwang，WEI Xiucheng. Wave equation modeling with implicit difference technique[J]. Geophysical Prospecting for Pe- troleum，2006，452151-156. [4] 蔡其新，何佩军，秦广胜，等．有限差分数值模拟的最小频散 算法及其应用[J]．石油地球物理勘探，2003，383247-248． CAI Qixin，HE Peijun，QIN Guangsheng，et al. Minimum dis- persion algorithm by finite-difference numeric modeling and its applications[J]. Oil Geophysical Prospecting，2003，383 247-248. [5] 成景旺，顾汉明，刘琳，等．海上粘弹介质 FCT 有限差分正演 模拟[J]．石油天然气学报，2011，331283-85． CHENG Jingwang，GU Hanming，LIU Lin，et al. Flus correction transport FCT finite difference forward modeling for marine viscoelastic media[J]. Journal of Oil and Gas Technology，2011， 331283-85. [6] 李文杰，张改兰，姜大建，等．利用 FCT 方法压制弹性波数值 模拟中的数值频散[J]． 物探化探计算技术， 2011， 333 248-251. LI Wenjie，ZHANG Gailan，JIANG Dajian，et al. Suppressing numerical dispersion in finite difference modeling based on acous- tic and elastic equation using FCT scheme[J]. Computing Tech- niques for Geophysical and Geochemical Exploration，2011， 333248-251. [7] 吴国忱， 王华忠. 波场模拟中的数值频散分析与校正策略[J]． 地 球物理学进展，2005，20158-65. WU Guochen，WANG Huazhong. Analysis of numerical disper- sion in wave-feild simulation[J]. Progress in Geophysics，2005， 20158-65. [8] LIU Yan，SEN M K. A practical implicit finite-difference examples from seismic modeling[J]. Journal of Geophysics and Engineering，2009，63231-249. [9] 任志明，刘洋．一阶弹性波方程数值模拟中的混合吸收边界条 件[J]．地球物理学报，2014，572595-605. REN Zhiming，LIU Yang. Numerical modeling of the first-order elastic equations with the hybrid absorbing boundary condition[J]. Chinese Journal of Geophysics，2014，572595-605. [10] 张文生．波动方程成像方法及其计算[M]．北京科学出版社， 200976. [11] CERJAN C，KOSLOFF D，KOSLOFF R，et al. A nonreflecting boundary condition for discrete acoustic and elastic wave equa- tions[J]. Geophysics，1985，504705-708. [12] 邵秀民，刘臻．带吸收边界条件的声波方程显式差分格式的稳 定性分析[J]．计算数学，2001，232164-182． SHAO Xiumin，LIU Zhen. Stability analysis of explicit finite difference schemes for the acoustic wave equation with absorbing boundary conditions[J]. Mathematica Numerica Sinica，2001， 232164-182. 责任编辑 聂爱兰 ChaoXing