声波方程有限差分数值模拟的变网格步长算法

2 时间积分
均匀介质中的二维声波方程可用下式表示[ 3]
2P t2
=
-
L2 P +
s
( 1a)
2
2
- L2 = c2 [ x 2 + z 2 ]
( 1b)
这里 P = P ( x , z , t) 代表压力项, c = c( x , z )
是速度, s = s( s, z ) 是震源函数。在密度 = ( x ,
209
图 1 低速地表差分网格分布示 意图 F ig . 1 Sketch map of finite difference g rid distr ibutio n for w eathering lay er
图 3 二维模型 F ig . 3 2D mo del
图 2 低速夹层差分网 格分布示意图 F ig . 2 Sketch map of finite difference
算[ 5] 。为简便起见, 仅对基本方程( 1) 的快速展开
法做一些介绍, 但对变密度的情况可用类似的方
法进行推导。
对震源项 s = ( x , z , t ) = g( x , z ) h( t ) 边界条
件为吸收边界时, 式( 1) 的常规解为
P(x , z, t) =
[
sin( L) h( t 0L
基金项目: 中国石油大学( 华东) 研究生创新基金( 编号: S2006- 06) 资助。 作者简介: 李胜军( 1979 ) , 男, 中国石油大学( 华东) 资信学院硕士研究生, 研究方向为地震波传播理论。E_mail: hdpul is@ 126. com
孙成禹( 1968 ) , 男, 博士, 副教授, 主要从事地震勘探理论和方法的教学与科研工作。
R
R) h( t -
)d
( 4b)
Qk 指修改了的契比雪夫多项式, Jk 为 k 阶精
#
度的贝赛尔函数。如果空间参数- L2 的特征值 很靠近实轴, R 2 选择的比最大特征值大, 当 N >
#
Rt 时式( 4a, b) 按指数规律衰减。参数- L2 的最
大特征值由 x 方向上的差分参数的最大特征值
针对图 3 展示的简 单模型, 用固定 步长 x 和 z 采样, 可以很好的满足低速层的采样需求, 但是这样造成了对高速层的过采样。为了避免过 采样可以从某一深度 z 0( 图 1) 开始采用双倍的步 长来采样。在修改的网格上, x 方向的导数用傅
第3期
李胜军等: 声波方程 有限差分数值模拟的变网格步长算法
N + n) z 的网格上的导数可写作
N- n
∀ D
2 z
[
P
]
i,
j
=
(
1 z)
2
k=
1
∀#[
Pi, j+
k
-
2P i, j
N
∀ +
Pi, j- k ] +
(
2
1 z
)
2
k=
1
∀#[
P
i,
j+
2( k+
n)
- 2 Pi, j + P ] i, j- 2( k+ n)
( 8)
21 0
工 程 地 球 物 理 学 报 ( Chinese Jour nal of Engineer ing Geo phy sics)
适用, 所以 z 方向上的导数用高阶有限差分法近
似。
为了获得一个稳定的算法, 使用对称的有限
差分参数。记 P ( i x , j z ) = P i, j , 一种用对称有
限差分参数近似的二阶导数定义如下[ 6]
∀ D
2 z
[
P
]
i,
j
=
!2 P !z 2
=
i, j
(
N
1 z )2
k=
∀#
I
! [ Pi, j+ k - 2P i, j + P i, j- k ]
( 7)
这里 2N 是算 子长度[ 7] 。如 果在深 度 z 0 =
l z 的步长是双倍的, 式( 6) 在长度为 2N 的有限
差分法可被使用到深度为( l - N ) z , 式( 7) 可在
深度 l z 以下使用。在剩余的( N - 1) 个网格点,
可结合式( 6) 和式( 7) 来计算其导数。在深度( l -
傅立叶变换法和高阶有限差分法( F D, F init e Diff erence) 已成为计算声波方程空间导数的标准 技术[ 2, 3] 。两种方 法用的都是固 定的网格 步长。 由于采用了均匀网格, 这些方法在一些实例中的 有效性受到限制。然而, 采用变网格算法将能改 进有上覆低速层情况模拟结果的有效性。对地层 中间有超薄夹层的情形, 必须用精细网格覆盖才 能精确的对地层进行模拟, 应用这种变网格算法 既能实现对夹层的模拟, 又能保障计算量不增加。 因此这种通过整数 m 实现在任意给定深度上网 格步长变化的有限差分方法被推广[ 4] 。为了计算 空间导数, 在 x 方向用傅立叶变换法或有限差分 算法, 在 z 方向使用高阶有限差分方法。通过用 对时间积分 的快速展开法( REM , Rapid Ex pan sion M et hod) 来 保障 差分 方 法的 计算 精 度[ 5] 。 这种差分技巧比二阶时间差分有较高的精确度, 且计算用时短。
文献标识码: A
收稿日期: 2007 03 19
Acoustic Equation Numerical Modeling on a Grid of Varying Spacing
L i Shengjun1, 2 , Sun Chengyu1, N i Changkuan1, Zhang Yuhua1
g rid distr ibution for embedded layer
立叶法或有限差分法计算, x 通过采样定理来确
定。如, 在最短波长上至少有 2 个采样点。由于
x 在 x 方向上是连续变 化的, 用传 统的差分算
法进行计算是没什么问题的。但是在 z 方向上,
是不连续的, 基于连续变化的傅立叶变换法不再
Abstract: It is an import ant met hod t o study the feat ure of seismic w ave propag at es in the me dium by t he forw ard modeling of seismic w av e f ield. T radit ional F ourier spect ral met hod and finit e diff erence use unifo rm grids in numerical simulatio n of low - speed/ hig h- speed m edi um , t hin layer/ t hick lay er . So t he t wo methods are all lack of f lexibilit y. In t his paper, an al g orit hm o f mov ing g rid is recomm ended, w hich has solved the pro blem descr ibed above. T he diff erence o f sto rage requirement and comput at ional ef ficiency betw een t radit ional finit e diff erence and t he algorit hm o f mov ing grid is discussed. And the boundary condit ion, rapid ex pansion method ( REM) f or t he tim e int eg rat io n of t he algo rithm of mo ving gr id are int ro duced and t he advant ages of t his alg orit hm are sum marized. Key words: f orw ard modeling; mov ing gr id; boundary condit io n; grid spacing
z ) 变化的情况下, 常用的是 Vidale 给出的公式
2P t2
=-
L2
P
-1 L2 P + P源自L2 P + s ( 2a)
- L2 =
c2 2[
2
x2 +
2
z2 ]
( 2b)
波动方程的时间积分可通过 Ko slof f 等提出
的快速 展开法 ( Rapid Expansion M et hod) 来计
第4卷
如果 z 方 向上的步长通过 其它的参 数来调
整, 方程 可用 类似 的方 法推 导。一 般 通过 函数
20 8
工 程 地 球 物 理 学 报 ( Chinese Jour nal of Engineer ing Geo phy sics)
第4卷
1引 言
随着油气勘探工作的深入, 现代油气勘探工 作正面临着勘探目标越来越复杂和勘探精度越来 越高的挑战, 如此艰巨的勘探任务要求我们必须 提高对复杂勘探目标中地震波传播规律的认识, 并在此基础上, 采用一种特殊的、有针对性的方法 和技术, 对薄层和小的断块进行精确模拟[ 1] 。声 波在介质中传播的正演模拟研究, 能为我们精确 模拟地震波在复杂介质中的传播、做好精确勘探 提供理论基础。
( 1. College of Geo - Res our ces and I nf or mation , China Univer sity of P etr oleum , D ongy ing Shandong 257061, China;