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