波动方程的变步长有限差分数值模拟

参数的最大特征值 Ａｘ 和 Ｚ 方向上的相应 特征值 Ａｚ 的和来决定。这里， Ｒ２ 可由（ ８） 式计算
! " ２ ２
Ｒ ＝ｃ
Ａｘ !!ｘ
"２ ＋
Ａｚ !!ｚ "２
ｍａｘ
（ ８）
２ 空间导数计算
通常情况下， 勘探区域的地表较为复杂且介质 中波的传播速度较低。为较好地压制数值频散并提 高模拟结果的分辨率， 就要在整个计算区域采用小 网格步长来做差分数值模拟。这样就造成了对高速 区域的过采样， 浪费了计算机资源。为减小内存需 求和节省计算时间， 可以用图 １ 所示方法进行模 拟。对于地层中有低速层的情形可以采用图 ２ 所示 方法来达到提高分辨率并且计算量增加不大的差 分算法。
１
"
（ ９）
对应的双步长网格定义如下
２
Ｎ
& ２
Ｄｚ
$Ｐ
% ｉ， ｊ
＝
"
Ｐ
２
"ｚ
＝
ｉ， ｊ
１ !２! ｚ
"２ ｋ ＝ １ $ｋ !Ｐｉ， ｊ＋２ｋ － ２Ｐｉ， ｊ ＋Ｐｉ， ｊ－ ２ｋ "
（ １０）
式中： ２Ｎ 是算子长度［８］。
第５卷 第３期
李胜军： 波动方程的变步长有限差分数值模拟
·７·
若在深度 ｚ０＝ｌ!ｚ 的步长是双倍的，（ ９） 式在长 度为２Ｎ的有限差分法可被使用到深度为（ ｌ－ Ｎ） !ｚ，
·８·
油气地球物理
２００７年 ７ 月
为压制频散， 提高模拟精度， 只有减小采样间隔， 如 用常网格 ４ｍ 步长算法进行模拟， 内存需求将会扩 大 １ 倍， 以上述（ １６００ｍ ×１６００ｍ ） 的模拟区域为例， 计算时间将会增大 １．２ 倍左右。用文中所述变步长 算法进行模拟， 低速地表用 ４ｍ 小步长， 地表以下用 ８ｍ 大步长， 变网格步长算法的内存需求仅是常网格
ｚ
３０００ｍ ／ｓ
图 ３ 二维模型
为了获得一个稳定的算法， 使用对称的有限差 分参数。记 Ｐ（ ｉ!ｘ， ｊ!ｚ） ＝Ｐｉ， ｊ， 一种 用对称有 限差分 参 数 近 似 的 二 阶 导 数 定 义 如 下［７］
２
Ｎ
& ２
Ｄｚ
$Ｐ
% ｉ， ｊ
＝
"
Ｐ
２
"ｚ
＝
ｉ， ｊ
１ !!ｚ
"２ ｋ
＝
#ｋ !Ｐｉ， ｊ＋ｋ － ２Ｐｉ， ｊ ＋Ｐｉ， ｊ－ ｋ
＝
１ #!ｚ $２ ｋ
＝
!
１
#Ｐ － ｋ ｉ， ｊ＋ｍ（ｋ）
$ ２Ｐｉ， ｊ ＋Ｐｉ， ｊ－ ｍ（ｋ）
（ １２）
差分因子
!
必须根据函数
ｋ
ｍ（
ｋ）
计算， 这可通过
Ｆｏｒｎｂｅｒｇ（ １９８８ａ） 提出的一种算法来实现［９］。图 ４ 表
示了方程（ １２） 在双倍步长区域差分算子 ２Ｎ＝ １０的
情况下有限差分参数的变化情况。显示了 ６ 个网格
２００７ 年 ７ 月
油气地球物理 ＰＥＴＲＯＬＥＵＭ ＧＥＯＰＨＹＳＩＣＳ
第５卷 第３期
波动方程的变步长有限差分数值模拟 ＊
李胜军 １， ２） 孙成禹 １） 张玉华 １） 倪长宽 １） １） 中国石油大学地球资源与信息学院； ２） 中石油勘探开发研究院西北分院
摘要： 有限差分算法是常用的正演模拟方法之一， 其包含的地震信息丰富， 且实现简单。传统的有限差分方法通常都 采用均匀网格步长， 在对含低速 ／ 高速介质、薄层 ／ 厚层介质的模型进行波场模拟时往往缺乏稳定性。文章介绍了一 种可以有效解决上述问题的变网格算法， 对常规有限差分法与变网格差分算法在内存需求、计算速率等方面的差别 进行了比较， 对变网格差分算法中的边界条件、时间积分的快速展开算法作了阐述， 进而总结了变网格算法的优点。 关键词： 变步长； 边界条件； 计算时间； 快速展开法； 数值模拟
０
０
０
０
０
ｍ（ １） １
１
１
１
１
２
ｍ（ ２） ２
２
２
２
３
４
ｍ（ ３） ３
３
３
４
５
６
ｍ（ ４） ４
４
５
６
７
８
ｍ（ ５） ５
６
７
８
９ １０
１８００ｍ ／ｓ ８００
１６００ （ ｍ）
２２００ｍ ／ｓ 图 ６ 低速地表地质模型
在 内 存 为 １．０Ｇ Ｂ ｙｔｅ、Ｃ ＰＵ 为 奔 腾 ＩＶ －２．８Ｇ Ｈ ｚ 的 微型计算机上模拟试算， 分别用常网格有限差分算 法和变网格有限差分进行模拟运算形成单炮记录。 在此基础上对 ４ｍ 常 网 格 步 长 、８ｍ 常 网 格 步 长 、变 网格步长算法从内存需求和计算时间进行比较分 析。图 ７（ ａ） 是常网格 ８ｍ 步长模拟的单炮记录， 可 以看出， 由于在低速地表的采样点太少， 图中的频散 现象比较严重， 甚至影响到了倾斜界面的反射记录。
$# % Ｐ !ｘ， ｚ， ｔ "＝ ｔ ｓｉｎ !!Ｌ "ｈ !ｔ－ ! "ｄ ! ｇ !ｘ， ｚ "（ ５） ０Ｌ 该常规解通过契比雪夫展开近似为
$ ! "% &Ｎ／２
Ｐ !ｘ，
ｚ，
ｔ
"
ｋ
＝
ｂ２ｋ＋１
０
Ｒ ｉＬ!
Ｑ２ｋ＋１
ｉＬ! Ｒ
ｇ !ｘ， ｚ " （ ６）
式中： ! 为变量； Ｌ! ２ 是空间差分算子 Ｌ２ 的数值 逼近；
在 地 震 资 料 采 集 、处 理 和 解 释 中 通 常 需 要 进 行 地震波场数值模拟： 假设已知地下的地质情况， 应 用地震波运动学和动力学的基本原理， 计算给定地 质模型的地震响应。这种做法对正确认识地震波的 运动学和动力学特征， 以及准确分析油气藏的反射 波场特征有着重要的指导意义。声波在介质中的正 演模拟研究为我们精确模拟地震波在复杂介质中 的 传 播 提 供 了 理 论 基 础［１］。
列表（ 黑圆点位置属于精细网格循环； 方框点属于 双倍的步长循环点； 五角星位置点的导数被计算） 。 在每一列用不同的差分模式来计算在不同深度的二
阶导数。针对十阶精度有限差分相应的 ｍ（ ｋ） 的值列 于表 １， 其他阶精度的 ｍ（ ｋ） 值可用相应的方法计算。
ｚ０ ｚ
图 ５ 必须内插压力的网格点
３ 算例分析
０ ０．０
８００
１６００ （ ｍ）
０
０．０
８００
１６００ （ ｍ）
０．５
０．５
１．０
１．０
１．５
（ ｓ）
（ ａ） 均匀网格
１．５
（ ｓ）
（ ｂ） 变网格
图 ７ 常网格与变网格得到的炮集记录
４ 结论
上述分析与理论模型试算结果表明， 变网格步 长差分算法是一种高效的正演模拟算法。它能较好 地实现对低速地表、高速层夹层、复杂界面的数值模 拟， 成功地解决小网格步长差分算法内存需求量大、 计算效率低等问题。在特殊地区用增大（ 减小） 步长 的模拟方法减小了内存需求， 提高了计算效率； 在特 定位置减小步长， 解决了对超薄层地质体、倾斜界面 的数值模拟效果不理想的问题： 解决了用常规高阶 有限差分或傅立叶变换法模拟存在的问题。
" $ ２
!
Ｐ
２
＝－
２
Ｌ
!ｔ
Ｐ !
－
１
２
２
Ｌ Ｐ＋ＰＬ
１
＋ｓ
!
!
（ ３）
% $ ２
－
２
Ｌ＝
!ｃ
２
２
２
!
２＋
!
２
!ｘ !ｚ
（ ４）
收稿日期： ２００７－０３－２３； 修订日期： ２００７－０４－２７ 作者简介： 李胜军， 男， 在读硕士研究生， 研究方向为地震 波 传 播 理 论 。 联 系 电 话 ：（ ０５４６） ８３９２０５５， Ｅ－ｍａｉｌ： ｈｄｐｕｌｉｓ＠１２６．ｃｏｍ， 通 讯 地 址 ： （ ２５７０６１） 中国石油大学（ 华东） 地球资信与信息学院。 ＊ 中国石油大学（ 华东） 研究生创新基金资助， 编号： Ｓ２００６—０６。
４ｍ 步长算法的 ６５％ 左右， 计算时间也缩短 ４５％ 左 右。与常网格 ８ｍ 步长算法相比， 计算时间是常网格 ８ｍ 步长的 １．２ 倍， 当然， 内存需求也相应有所扩大。 图 ７（ ｂ） 是变网格步长算法的模拟结果（ 为了对比 明显， 两图使用了相同的增益） ， 可以看出， 频散得 到了较好的压制， 分辨率也明显提高。
因而， 采用变网格算法将能改进有上覆低速层 情况模拟结果的有效性（ 对地层中间有超薄夹层的 情形， 必须用精细网格覆盖才能精确的对地层进行 模拟） 。应用这种变网格算法既能实现对夹层的模 拟， 又能保障计算量不增加。因此这种通过函数实 现在任意深度上网格步长变化的有限差分方法被
推广［４］。为了计算空间导数， 在 Ｘ 方向用傅立叶变换 法或有限差分算法， 在 Ｚ 方向使用高阶有限差分方 法。通过时间积分快速展开法（ Ｒ Ｅ Ｍ ） 来保障差分 方法的计算精度［６］。这种差分技 巧比二阶 时间差分 有较高的精确度且计算用时短。
（ ７） 式可在深度 ｌ!ｚ 以下使用。对剩余的（ Ｎ－ １） 个
网格点， 可以结合（ ９） 和（ １０） 式来计算其导数。在
深度（ ｌ－ Ｎ＋ｎ） !ｚ 的网格上的导数可写作
Ｎ－ ｎ
% ２
Ｄｚ
!Ｐ
" ｉ， ｊ
＝
１ #!ｚ $２ ｋ
＝
!
１
ｋ #Ｐｉ， ｊ＋ｋ －
２Ｐｉ， ｊ ＋Ｐｉ， ｊ－ ｋ
$
Ｎ
% ＋
Hale Waihona Puke Baidu ·６·
油气地球物理
２００７年 ７ 月
波动方程 的 时 间 积 分 可 通 过 Ｋ ｏｓｌｏｆｆ等 提 出 的 快速展开法（ Ｒ Ｅ Ｍ ） 来计算［６］。为 简便起见 ， 仅对基 本方程（ １） 的快速展开法做一些介绍。对变密度的 情况可用类似的方法进行推导。
对震源项 ｓ（ ｘ， ｚ， ｔ） ＝ｇ（ ｘ， ｚ） ｈ（ ｔ） ， 边界条件为 吸收边界时， 式（ １） 的常规解为
图 ３ 给出了简单的模型： 用固定小网格步长 !ｘ 和 "ｚ 采样， 可以很好地满足低速层的采样需求， 但 是这样造成了对高速层的过采样。为避免这种过采 样， 可以从某一深度 ｚ０（ 图 １） 开始采用双倍的步长 来采样。在修改的网格上， Ｘ 方向的导数用傅立叶法 或有限差分法计算， #ｘ 通过采样定理来确定。由于
傅立叶变换法和高阶有限差分法（ ＦＤ ） 已成为 计算声波方程空间导数的标准技术［２， ３］。虽然常网格 步长差分算法比较容易实现， 但是它们对大部分模 型都增大了不必要的计算量。例如， 对存在浅层低 速带的沉积盆地模型地面地震记录进行模拟时， 由 于低速地层阻抗小， 地震波传入其中会引起较大的 振幅和较长的延续时间（ 这与深层的高速层完全不 同） 。由于这些浅层低速层中地震波的波长较短、地 层厚度较小， 模拟时需要用小网格进行。这样， 常网 格步长算法就必须用小网格离散整个模型， 从而增 加了不必要的代价， 如内存、计算量的增大。
为验证上述算法的正确性、有效性， 设计了图 ６ 所示的低速地表地质模型。数值模拟时， 震源设置在 （ ８００ｍ ， ２４ｍ ） 处。
ａｂ
ｃ ｄｅｆ
０
８００
１６００ （ ｍ）
０
１０００ｍ ／ｓ
ｚ０
图 ４ 过渡带的差分格式
表 １ 图 ４ 显示的有限差分参数 ｍ（ ｋ） 值
ｍ（ ｋ） ａ
ｂ
ｃ
ｄ
ｅ
ｆ
ｍ（ ０） ０
在图 ２ 显示的网格中， 可使用这种变网格算法 来改变空间分辨率。在细网格区域不需要用傅立叶 变换法或高阶有限差分算子， 四阶或六阶精度的算 法就可以满足精度的要求。
ｘ
如果 Ｚ 方向上的步长通过其他参数来调整， 类 似的方程可被推导。一般通过函数 ｍ（ ｋ） 来产生， 有
Ｎ
% ２
Ｄｚ
!Ｐ
" ｉ， ｊ
１ 时间积分
均 匀 介 质 中 的 二 维 声 波 方 程 可 用 下 式 表 示［２］
２
!
Ｐ
２
２
＝－ Ｌ Ｐ＋ｓ
（ １）
!ｔ
" # ２
２
２２
－ Ｌ ＝ｃ
!２＋
!
２
（ ２）
!ｘ !ｚ
式中： Ｐ＝Ｐ（ ｘ， ｚ， ｔ） ， 代表压 力项； ｃ＝ｃ（ ｘ， ｚ） ， 代表速 度； ｓ＝ｓ（ ｓ， ｚ） ， 代表震源函数； Ｌ２ 为差分算子。 在 密 度 !＝!（ ｘ， ｚ） 变 化 的 情 况 下 ， 常 用 的 是 Ｖ ｉｄａｌｅ 给 出 的 公 式［５］
$ｘ 在 Ｘ 方向上是连续变化的， 用常规的有限差分 算法就可计算。但是在 Ｚ 方向上， %ｚ 是不连续的， 基 于连续变化的傅立叶变换法不再适用， 所以 Ｚ 方向上的导数要用高阶有限差分法近似。
ｘ
ｚ
ｚ０
图 １ 低速地表模型差分网格分布示意图 ｘ
ｚ
图 ２ 低速夹层模型差分网格分布示意图 ｘ １５００ｍ ／ｓ
Ｑｋ 指修改了的契比雪夫多项式； 系数 ｂｋ 为
# ｂｋ＝
ｔ
Ｊ２ｋ＋１
!!Ｒ
" ｈ
!ｔ－
!
"ｄ!
０Ｒ
（ ７）
式中： Ｊｋ 为 ｋ 阶精度的贝赛尔函数。 如果空间参数 － Ｌ! ２ 的特征值很靠近实轴、Ｒ２ 比
最大特征值大， 当 Ｎ＞Ｒ·ｔ 时， （ ６） 、（ ７） 式按指数规 律收敛。参数 － Ｌ! ２的最大特征值由 Ｘ 方向上的差分
１ #２!ｚ
$２ ｋ
＝
!
１
# $ Ｐ － ２Ｐ ＋Ｐ ｋ ｉ， ｊ＋２# ｋ＋ｎ $ ｉ， ｊ ｉ， ｊ－ ２# ｋ＋ｎ $ （ １１）
为保证这种差分算法的稳定性， 压力必须内插 在水平方向两倍步长边界上的一些点上。对参数长 度 ２Ｎ＝ １０ 的那些点在图 ５ 中已列出。内插计算时 可使用特别精确且没有错误发生的三角内插算法。