数学物理方程--有限差分法-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN数学物理方法课程报告题目：声波有限差分法数值模拟学生姓名：xxx学号：xxx学院：地球科学与技术学院专业班级：xxxx教师：xxx2016年 4月12日声波有限差分法数值模拟Xxx（地球科学与技术学院研15级 学号：xxx ）摘要：数值模拟是最常用的正演模拟的方法。

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

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

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

关键词：数值模拟，声波，有限差分正文1、 引言在勘探过程中，数值模拟的作用很大。

例如：1、采集上，可用于设计或者优化野外观测系统；2、处理上，可以通过数值模拟来检验是否采用了正确的反演方法。

将正演反演不断的逼近，从而使结果更加准确；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。

利用Taylor 级数展开求解()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!21221x 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 ∆++∆+∆=∆-+-∆+(1)将相减后得到的式子整理成矩阵形式，有 ()()()()()()()()()()()()()()()()()()()⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∆-+-∆+∆-+-∆+∆-+-∆+∆=⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∆∆⨯⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡-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 NNN N0002002010010202220420224222422221412122221!21!41!21(2)为了简化矩阵，可以记作⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=N N NNN Nq q q q q q q q q A 242224222214121 ，()()()()()()()()()()⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∆-+-∆+∆-+-∆+∆-+-∆+∆=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 0002002010010222221同时，构造两个简单矩阵，辅助计算N N I ⨯⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=111 整理的， 1001⨯⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=N E 假设存在1-A ，使得I AA =-1，也可得()I A A T T=-1；即()TA 1-为T A 的逆，得到()I A A TT =-1。

式子两边右乘向量E 就可得()E E A A TT =-1 (3)由式（2）可得()()D A E x f T 10221-= (4)同时，假设()()TN T Tc c c C E A ,,,211 ==-（5）将()N c c c C ,,,21 =带入式（4），得()()()()()()[]x q x f x f x q xf c x x f n n Nn n ∆-+-∆+∆=∑=000120222121（6）整理得 ()()()()()()[]x q x f x q x f c x f c x fx n n Nn n ∆-+∆++=∆∑=00100022可结合式（3）和式（5），可得到矩阵计算式：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡00121222214424122221 N N N N N N N c c c q q q q q q q q q（7）∑=-=Ni i c c 102当i q 的值确定后，可根据式（7）来求解n c 的值，从而计算出()()01x f 的值。

利用式（7）可以求得对称任意节点间距的一阶导数差分系数。

其中，当i q 取值为),2,1(N n n =,则式（7）可表示为()()⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡00121212121222444222N N NN c c c N N N （8）此时，所求得的()N n c n ,,2,1 =就是等节点间距的一阶导数的规则网格不同差分精度的差分系数（表1所示）。

二维声波方程的形式可表示为：()22222221zux u t u v p ∂∂+∂∂=∂∂ （9）时间导数采用2阶，空间导数采用2N 阶近似，即()())(2)(222t t u t u t t u tut ∆-+-∆+=∂∂∆()()()()()[]x n x u x n x u c x u c x x ux Nn n ∆-+∆++=∂∂∆∑=001000222带入式（9）中，可得到在固定网格下，差分格式为 ()()()()[]()()()[]⎭⎬⎫⎩⎨⎧∆-+∆++⎪⎪⎭⎫⎝⎛∆∆+⎭⎬⎫⎩⎨⎧∆-+∆++⎪⎪⎭⎫ ⎝⎛∆∆+∆--=∆+∑∑==Nn n p Nn 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 102102)(2)(（10）3、模型测试：震源选取：正演模拟过程中采用雷克子波作为震源子波，雷克子波的表达式为Source (it) =((1-2π f m (t-t0)2 )e-2π fm (t-t0)2模型建立：建立了一个两层介质模拟，其上层纵波速度为v=2000m/s，下层纵波速度为v=3000m/s。

模型大小为200×200，空间采样间隔为dx=dz=10m。

采用30Hz的雷克子波作为震源子波，震源位于模型(70,100)处，时间采样间隔为1ms。

结果分析：it=50 it=100 it=150it=200 it=250 it=300it=350 it=400图2 不同时刻波场快照图中可以看出，在未遇到界面前，地震波在均匀介质中的波前面一个圆。

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

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

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

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

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