η2
=
1 4
r21
+ 4 r1 s1 + r1 s1 ( r1 +
4 r2 s1 r2 )
,
(14)
η3 = -
1 4
r2 ( r1
+
r21 r2 ) ( r1
+
r2
+
s1 )
,
(15)
η4 = -
1 4
s1 ( r1
+
r21 s1 ) ( r1
+
r2
+
s1 )
,
(16)
同理 ,容易得到55 <y和55 <z . 如果是规则网格 ( s1 = r1 = r2 =Δx) ,则
图 1 交错网格有限差分示意图 Fig. 1 Staggered grids finite difference
笛卡儿坐标系斜方晶系各向异性模型弹性波动
方程速度 - 应力公式为[18]
ρv = D ·T ,
(1)
T = C ·DT ·v ,
(2)
其中 ρ( x) 为 介 质 密 度 , 弹 性 介 质 振 动 速 度 矢 量
第 47 卷 第 2 期
地球物理学报
Vol. 47 , No. 2
2004 年 3 月
CHINESE JOURNAL OF GEOPHYSICS
Mar. , 2 0 0 4
Sun W T , Yang H Z. A 32D finite difference method using irregular grids for elastic wave propagation in anisotropic media. Chinese J . Geo2 phys. (in Chinese) , 2004 , 47 (2) :332～337
r1 2
+
r2
3
+
54 < 5 x4
1 4!
r1 2
+
r2
4
+ …,
(11)

由上述方程组
,55
<表示为
x
<
的线性组合
5 5
<
x
= η1 <m +1
-
η2 <m
+ η3 <m+2
-
η4 <m- 1
,
(12)
η1
=
1 4
r21
+ 4 r1 r2 + 4 s1 r2 r1 r2 ( r1 + s1 )
,
(13)
1 2!
r1 2
+ s1
2
-
53 < 5 x3
1 3!
r1 2
+
s1
3
+
54 < 5 x4
1 4!
r1 2
+
s1
4
+ …,
(8)
<m
= <i
-
5 < r1 5x 2
+
52 < 5 x2
1 2!
r1 2
2
-
53 < 1 5 x3 3 !
r1 2
3
+
54 < 5 x4
1 4!
r1 2
4
+ …,
s1 = xm - xm- 1 ,
r1 = xm +1 - xm ,
r2 = xm +2 - xm+1 .
(7)
节点 m 和 m + 1 的中点命名为节点 i . 节点
m - 1 、m 、m + 1 、m + 2 处的波场值 < 写成级数形式
<m - 1 = <i
-
5< 5x
r1 2
+
s1
+
52 < 5 x2
(4)
5 5x
0
0
0
55 5z 5y
D=
0
5 5y
0
5 5z
0
5 5x
,
(5)
0
0
555 5z 5y 5x
0
c11 c12 c13 0 0 0
c12 c22 c23 0 0 0
c13 c23 c33 0 0 0
C=
. (6)
0 0 0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66 波动方程中包含速度和应力的一阶导数 ,下面 推导具有四阶空间精度和二阶时间精度的不规则网 格差分算子. 在交错网格 (见图 1) 上离散一阶偏微 分算子 ,笛卡儿坐标轴上四节点示意图见图 2 ,它们 之间的间距分别是 s1 、r1 、和 r2 .
基金项目 中国石油天然气集团公司基金资助 (2002CXKF24) . 作者简介 孙卫涛 ,男 ,1975 年生 ,1998 年毕业于大连理工大学工程力学系 ,2003 年在清华大学工程力学系获硕士 、博士学位 ,现在清华大学计
算机系做博士后研究 ,主要从事弹性波动力学 、高性能计算的理论和方法研究. E2mail : sunwt @mail . tsinghua. edu. cn
η1
= η2
=
9 8Δx
,
(17)
η3
= η4
=-
1 24Δx
,
(18)
5 5
<
x
=
9 8
( <m +1 Δx
<m )
-
1 24
( <m +2 - <m- 1 ) Δx
,
(19)
这是 Levander[5] 给出的规则网格四阶差分算子. 本
文提出具有更一般形式的不规则网格公式 ,模拟弹
(9)
<m +1
= <i
+
5 5
<
x
r1 2
+
52 < 5 x2
1 2!
r1 2
2
+
53 < 5 x3
1 3!
r1 2
3
+
54 < 5 x4
1 4!
r1 2
4
+ …, (10)
<m +2
= <i
+
5 5
<
x
r1 2
+
r2
+
52 < 5 x2
1 2!
r1 2
+
r2
2
+
53 < 5 x3
1 3!
2 期 孙卫涛等 :各向异性介质弹性波传播的三维不规则网格有限差分方法
333
1 引 言
20 世纪 70 年代以来 ,在地震波传播和地震地 表运动的数值模拟研究中 ,有限差分方法成为求解 波动方程的最有力工具之一. 地震波正演模拟中有 限差分方法的早期研究见 Alterman[1] , Kelly 等[2] 和 Virieux[3 ,4] 等学者的著作 ,Levander[5] 在 P2SV 波模拟 中引入了四阶空间差分算子 , Graves[6] 给出了等效物 质参数的三维四阶速度 - 应力有限差分方法 , Dablain[7] 提出了高阶差分算子的方法. 这些工作全 部基于笛卡儿坐标系中的规则网格 ,用普通网格模 拟曲线界面时出现“阶梯状”边界 ,在地形构造模型 中产生虚假衍射波. 另外 ,局部物理参数的变化也 会要求加密整个模型网格 ,导致计算量的增加. Shortley 等[8] 首先研究了 Laplace 方程中的不规则网 格有限差分方法 , Jastram 等[9] 和 Falk 等[10] 给出了 交错网格上的可变网格差分方法 , Tessmer 等[11] 和 Hestholm 等[12] 用变形的矩形网格模拟曲线边界 , Ivo Oprsal 等[13] 研究了非均匀介质中的二阶波动方程的 不规则网格差分方法 , Pitarka[14] 提出了各向同性介 质中不规则网格的有限差分方法 ,董良国等[15 ,16] 研 究了交错网格高阶差分方法 ,张剑锋[17] 通过坐标映 射得到非规则网格差分法. 但是 ,以往变形网格有 限差分方法需要波动方程在不同坐标系之间的正反 变换 ,正演模拟的计算量相当大 ;对于具有非平面界 面的三维非均匀各向异性介质模型 ,不规则网格有 限差分方法目前还鲜有相应的研究 ,普通有限差分 方法仍然面临着计算规模和计算精度的限制. 高精 度 、高计算速度正演成为三维复杂介质弹性波传播 问题的研究热点. 本文给出一种具有高阶精度 、非 均匀各向异性波动方程不规则网格有限差分方法. 这种方法无需坐标变换和网格间的插值 ,简单易行 而且占用内存少 、计算量小 ,适合于三维非均匀各向 异性介质模型的弹性波正演问题.
Abstract This paper presents a new 3D finite2difference ( FD) method using spatially irregular grids to sim2 ulate elastic wave propagation in heterogeneous anisotropic media with topographic structures. The method ap2 proximates the first2order elastic wave equations by the finite difference operators on irregular grids with sec2 ond2order time precise and fourth2order spatial precise . Unlike the multi2grid scheme , this method has no in2 terpolation between the fine and coarse grids. All grids are computed at the same spatial iteration. Complicat2 ed geometrical structures like rough submarine interfaces , faults and nonplanar interfaces are treated with fine irregular grids. Theoretical analysis and numerical simulations show that this method saves considerable memo2 ry and computing time , at the same time , has satisfactory stability and accuracy. The proposed scheme is more efficient than conventional methods in simulating seismic wave propagation in complex topographic struc2 tures. Key words Seismic wave , Irregular grid , Finite difference ,Anisotropy media.