数学物理方程有限差分法

数学物理方法课程报告题目：声波有限差分法数值模拟
学生姓名：ｘxx
学号：xｘx
学院:地球科学与技术学院
专业班级：ｘxxｘ
教师：xｘｘ
２016年 4月1２日
声波有限差分法数值模拟
Xxx
（地球科学与技术学院研15级 学号：x ｘx ）
摘要:数值模拟是最常用的正演模拟的方法。

它通过给出的结构模型和物理参数，模拟
地震波的传播轨迹,了解其规律以及过程，然后通过计算来推断观测点的地震记录。

根据求解方法，地震波方程数值解法可分为有限元法、伪谱法、有限差分法。

根据本门课程的要求，并且有限差分法具有内存占用较小，精度较高等优点,本文主要采用这种方法进
行模拟。

关键词：数值模拟,声波,有限差分
正文
1、 引言
在勘探过程中,数值模拟的作用很大.例如:１、采集上,可用于设计或者优化野外观测系统；２、处理上,可以通过数值模拟来检验是否采用了正确的反演方法。

将正演反演不断的逼近,从而使结果更加准确;3、解释上,还可以检测一下解释的资料是否正确。

而有限差分法是数值模拟最常用的方法，本文利用有限差分法，通过对声波进行正演模拟，来了解其在地下的传播规律及特点。

2、 二维各向同性介质声波方程数值模拟 使用规则网格差分对二阶方程进行求解。

具体过程:
在x 方向上，关于0x 对称分布的2N 个网格节点的坐标分别为x q x N ∆-0， x q x N ∆--10,……,x q x ∆-10，x q x ∆+10,……x q x N ∆+-10,x q x N ∆+0.其
中，x ∆表示节点间的最小间距；i q 表示任意正整数.2N 个网格节点所对应的函
数值已知，分别为()x q x f N ∆-0，()x q x f N ∆--10，……，()x q x f ∆-10， ()x q x f ∆+10……,()x q x f N ∆+-10，()x q x f N ∆+0.利用Ｔaylor 级数展开求解
()x f 在点0x 处的一阶导数近似值。

()()()()()()()()()()()()()[]
120220220100!
21
!
21
+∆+∆+
+∆+
∆+=∆+N i N N i i i i x q O x f x q N x f x q x f x q x f x q x f
()()()()()()()()()()()()()[]
120220220100!
21
!
21
+∆+∆+
+∆+
∆-=∆-N i N N i i i i x q O x f x q N x f x q x f x q x f x q x f
其中,i=1，2，…，N
将上述两式相加,省略式中的误差项，得到
()()()[]()()()()()()()()()()022*********!
21
!41!21221
x f x q N x f x q x f x q x q x f x f x q x f N N i i i i i ∆+
+∆+∆=∆-+-∆+
(１）
将相减后得到的式子整理成矩阵形式,有
()()()()()()()()()()()()()()()()()()()⎥⎥⎥⎥
⎦
⎤⎢⎢⎢⎢⎣⎡∆-+-∆+∆-+-∆+∆-+-∆+∆=⎥⎥
⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∆∆⨯⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣
⎡-x q x f x f x q x f x q x f x f x q x f x q x f x f x q x f x x f x N x f x x f q q q q q q q q q N N N N N N N
N N N
000200201001020222042
0224
2224
2222141
2
122221!21!41!
21
(２)
为了简化矩阵，可以记作
⎥⎥⎥
⎥⎥⎦
⎤
⎢⎢⎢⎢⎢⎣⎡=N N N
N
N N
q q q q q q q q q A 242224222
214
1
21 ，()()()()()()()()()()⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎣⎡∆-+-∆+∆-+-∆+∆-+-∆+∆=x q x f x f x q x f x q x f x f x q x f x q x f x f x q x f x D N N 00020020100102
22221
同时，构造两个简单矩阵,辅助计算
N N I ⨯⎥
⎥
⎥⎥⎦⎤⎢⎢⎢
⎢⎣
⎡=111 整理的, 1001⨯⎥⎥⎥⎥
⎦⎤⎢⎢⎢⎢⎣⎡=N E 假设存在1-A ，使得I AA =-1,也可得()
I A A T T
=-1
;即()T
A 1-为T A 的逆，得到
()
I A A T
T =-1
.式子两边右乘向量E 就可得
()
E E A A T
T =-1 （3)
由式（2）可得
()
()D A E x f T 1022
1-=ﻩﻩ(4） 同时，假设
()
()T
N T T
c c c C E A ,,,211 ==-
(5）
将()N c c c C ,,,21 =带入式(4)，得
()()()
()()()[]x q x f x f x q x
f c x x f n n N
n n ∆-+-∆+∆=∑=000
1
20222121
(6）
整理得 ()()
()()()()[]x q x f x q x f c x f c x f
x n n N
n n ∆-+∆++=∆∑=001
00022
可结合式（３）和式（5）,可得到矩阵计算式：
⎥⎥
⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡=
⎥⎥⎥⎥⎦⎤⎢
⎢⎢⎢
⎣⎡⨯⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡00121222
21
4424
12
2
2
2
1 N N N N N N N c c c q q q q q q q q q ﻩ
（７)
∑=-=N
i i c c 1
02
当i q 的值确定后，可根据式(7)来求解n c 的值,从而计算出()()01x f 的值。

利用式(７）可以求得对称任意节点间距的一阶导数差分系数.其中，当i q 取值为
),2,1(N n n =,则式(７）可表示为
()()⎥⎥
⎥⎥
⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣
⎡0012121212122244422
2
N N N
N c c c N N N ﻩﻩ(８) 此时,所求得的()N n c n ,,2,1 =就是等节点间距的一阶导数的规则网格不同差分精度的差分系数（表1所示）。

二维声波方程的形式可表示为：
()22222
221z
u
x u t u v p ∂∂+∂∂=∂∂ﻩﻩ(9） 时间导数采用２阶，空间导数采用２Ｎ阶近似,即
()())(2)(222
t t u t u t t u t
u
t ∆-+-∆+=∂∂∆
()()()()()[]x n x u x n x u c x u c x x u
x N
n n ∆-+∆++=∂∂∆∑=001
000222
带入式(9)中，可得到在固定网格下,差分格式为
()()()()[]()()()[]⎭
⎬
⎫⎩⎨⎧∆-+∆++
⎪⎪⎭
⎫
⎝⎛∆∆+⎭⎬⎫⎩
⎨⎧∆-+∆++⎪⎪⎭⎫ ⎝⎛∆∆+∆--=∆+∑∑==N
n n p N
n n p z n z u z n z u a z u a z t v x n x u x n x u a x u a x t v t t u t u t t u 102
10
2
)(2)( ﻩ(10）
3、模型测试：
震源选取：
正演模拟过程中采用雷克子波作为震源子波,雷克子波的表达式为
Ｓｏurce （i ｔ） =((1-２π f ｍ （t －t ０）2
)e
—2π f ｍ （t —ｔ0
)2
模型建立:
建立了一个两层介质模拟,其上层纵波速度为v=200０m/s ，下层纵波速度为v ＝3000m/s.模型大小为2０0×２00,空间采样间隔为dx=dz=10m 。

采用30Hz 的雷克子波作为震源子波，震源位于模型（7０,10０）处,时间采样间隔为1ms 。

结果分析：
it=50 it=１0０ it=1５0
it＝200 it＝２５0
it=300
iｔ=35０ｉt=400
图２不同时刻波场快照
图中可以看出,在未遇到界面前，地震波在均匀介质中的波前面一个圆。

当遇到地层界面之后,在界面处发生了反射、透射和折射现象。

沿测线方向的地震记录如图２所示。

记录中存在两条直线状的同相轴和两条近似双曲线的同相轴。

由于直达波的时距曲线是直线，因此两条直线同相轴对应直达波;由于反射波的时距曲线是近似双曲线,因此近似双曲线同相轴对应的是反射波。

参考文献
［1］刘庆敏,高阶差分数值模拟方法研究与应用，中国石油大学（华东）硕士论文，20０4年9月
[2］孙成禹、李振春，地震波动力学基础，石油工业出版社，2011年4月
［3］王元名，数学物理方程与特殊函数,高等教育出版社，20１２年１2月。