波动方程的变步长有限差分数值模拟
5第五讲 典型模型方程-波动方程的差分格式
修正方程：
u u c x 2u c x2 3u 1 2 2 1( 1) 3 c x x t x 2 6
a t x
稳定条件 ：
1
耗散及频散特性. G 1 c c cos ic sin
Euler差分方法（如下）不稳定
为了使其为稳定差分方法，修正如下
：
因此，此差分方程不一定总相容于PDE，当 则相容，为条件相容。
t 0, x 0 x
t 0, x 0
2
t
0
?
时，v保持为常数，
Lax差分方法的耗散性和色散性分析
当v不为1时，耗散性很强。
理解下图？？？
Euler差分方法（如下）不稳定
初值问题的精确解
由数学物理方程得：
（3） （4）
由（4）分别求t、x导数，代入（2）式可以验证解的正确性， 即： cF cF 0
欧拉显格式:Euler Explicit Method
都为一步显格式
1、 2、
： ：
1 un j
时间、空间前差
时间前差、空间中心差
两种差分方程的时间项都为一阶精度。 每个方程中未知量只有 ，都为一步显格式 显格式: 递推;而不同求解方程组
蛙跳格式的耗散性和色散性分析
蛙跳格式的缺点
1、初始值需要给两个时间层上的值。解决：可以 在初始计算时，仅给n＝0时刻的初始值，采用两层 差分格式来得出n＝1层上的值，然后在采用蛙跳格 式循环计算。
n u 不依赖于 j ，导致了 2 n 1 3 1 2n 4 2 两个不相关的解。 u j ...u j u j ，u j ...u j u j
一阶迎风格式的放大因子：
地震波波动方程数值模拟方法(严选优质)
地震波波动方程数值模拟方法地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶变换法、有限元法和有限差分法等。
克希霍夫积分法引入射线追踪过程，本质上是波动方程积分解的一个数值计算，在某种程度上相当于绕射叠加。
该方法计算速度较快，但由于射线追踪中存在着诸如焦散、多重路径等问题，故其一般只能适合于较简单的模型，难以模拟复杂地层的波场信息。
傅里叶变换法是利用空间的全部信息对波场函数进行三角函数插值，能更加精确地模拟地震波的传播规律，同时，利用快速傅里叶变换(FFT)进行计算，还可以提高运算效率，其主要优点是精度高，占用内存小，但缺点是计算速度较慢，对模型的适用性差，尤其是不适应于速度横向变化剧烈的模型.波动方程有限元法的做法是：将变分法用于单元分析，得到单元矩阵，然后将单元矩阵总体求和得到总体矩阵，最后求解总体矩阵得到波动方程的数值解；其主要优点是理论上可适宜于任意地质体形态的模型，保证复杂地层形态模拟的逼真性，达到很高的计算精度，但有限元法的主要问题是占用内存和运算量均较大，不适用于大规模模拟，因此该方法在地震波勘探中尚未得到广泛地应用。
相对于上述几种方法，有限差分法是一种更为快速有效的方法。
虽然其精度比不上有限元法，但因其具有计算速度快，占用内存较小的优点，在地震学界受到广泛的重视与应用。
声波方程的有限差分法数值模拟对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：)()(2222222t S zu x u v t u +∂∂+∂∂=∂∂ （4-1） (,)v x z 是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，)(t s 为震源函数。
为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。
为此，先把空间模型网格化（如图4-1所示）。
设x 、z 方向的网格间隔长度为h ∆，t ∆为时间采样步长，则有：z∆,i j1,i j +2,i j+1,i j-h i x ∆= （i 为正整数）h j z ∆= (j 为正整数)t n t =∆ (n 为正整数)k j i u , 表示在(i,j)点，k 时刻的波场值。
CSC频率—空间域波动方程数值模拟
CSC频率—空间域波动方程数值模拟吕晓春;顾汉明;成景旺;周丽【摘要】针对频率空间域波动方程数值模拟需要巨大内存空间的现状,提出了利用列索引压缩存储(CSC)技术存储大型稀疏非对称复数型的矩阵系数.CSC技术将系数矩阵转化为三个一维数组来存储,分别存储系数非零元素、非零元素对应所在的行以及每列起始非零元素所在位置.经CSC技术压缩存储后显著减少了内存空间及计算量,在计算时只有少许的非零元素参加计算,且根据三个一维数组可以简便地找到对应的非零元素,进而采用LU分解快速而精确地求解.本文基于Jo等提出的最优化9点差分方法,首次应用CSC技术在频率空间域进行二维声波方程数值模拟.通过对Corner-edge模型和二维Marmousi模型进行试算,可以显著节省内存需求,明显提高计算速度,进而得到精度较高的正演结果.【期刊名称】《石油地球物理勘探》【年(卷),期】2014(049)002【总页数】7页(P288-294)【关键词】频率—空间域;CSC;系数矩阵;一维数组;LU分解;数值模拟【作者】吕晓春;顾汉明;成景旺;周丽【作者单位】中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074【正文语种】中文【中图分类】P6311 引言为了反演地下介质参数或研究地震波在各种复杂介质中的传播机制，需要进行波场数值模拟。
波动方程的数值模拟方法包括有限差分法、有限元法、伪谱法等。
在时间—空间域的数值模拟技术已较成熟，并广泛应用于复杂介质的正演模拟中。
一维波动方程的有限差分法
学生实验报告实验课程名称偏微分方程数值解开课实验室数统学院学院数统年级2013 专业班信计02班学生姓名学号开课时间2015 至2016 学年第 2 学期数学与统计学院制开课学院、实验室：数统学院实验时间：2016年6月20日五．实验结果及实例分析1、u(x,t)在t=0.5,1.0,1.5,2.0时刻的数值解、精确解以及绝对误差表1 u(x,t)在t=0.5,1.0,1.5,2.0时刻的数值解时刻tt=0.5,1.0,1.5,2.0时刻的数值解t=0.5 0 -0.0059 -0.0113 -0.0155 -0.0182 -0.0192 -0.0182 -0.0155 -0.0113 -0.0059 0 t=1.0 0 -0.3090 -0.5877 -0.8090 -0.9510 -0.9999 -0.9510 -0.8090 -0.5877 -0.3090 0 t=1.5 0 0.0020 0.0038 0.0052 0.0061 0.0064 0.0061 0.0052 0.0038 0.0020 0 t=2.0 0 0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0表2 u(x,t)在t=0.5,1.0,1.5,2.0时刻的精确解时刻tt=0.5,1.0,1.5,2.0时刻的精确解t=0.5 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 t=1.0 0 -0.3090 -0.5878 -0.8090 -0.9511 -1.0000 -0.9511 -0.8090 -0.5878 -0.3090 0 t=1.5 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 t=2.0 0 0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0表3 u(x,t)在t=0.5,1.0,1.5,2.0时刻的绝对误差时刻tt=0.5,1.0,1.5,2.0时刻的绝对误差t=0.5 0 0.0059 0.0113 0.0155 0.0182 0.0192 0.0182 0.0155 0.0113 0.0059 0 t=1.0 0 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0001 0.0000 0.0000 0 t=1.5 0 0.0020 0.0038 0.0052 0.0061 0.0064 0.0061 0.0052 0.0038 0.0020 0 t=2.0 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0说明：在t=0.5时刻的绝对误差最大，t=1.5时刻次之，t=1与t=2时刻的绝对误差均较小，由于0.11r hτ==<，该格式稳定，由数值计算得到的矩阵不难看出，数值解符合理论解。
有限差分法地震波传播数值模拟
=
kΔx
=
2π λ
Δx ≤ 1
，即只要一个波长包含几个空间步
长，随着差分精度2M的提高，上述高阶差分解法产生
的数值频散会逐渐减小。
不同差分精度空间频散曲线
不同差分精度时间频散曲线
五、边界问题
自由边界条件
内部边界条件 吸收边界条件
设计吸收边界条件的目标：
z 方程+边界条件数学上是非病态的
连续
问题 z 方程+边界条件可以近似描述无限介质中的物理过程 z 边界条件和内部点的计算方式是相容、不冲突的
-----------J.M. Carcione
地震波传播数值模拟应用领域
地震波传播理论
数据采集
理论指导 物性参数
研究传播规律
正演模拟
指导设计 观测系统
验证
地震解释
提供理论数据 试验处理流程
数据处理
提供正演方法
岩石物理
参数反演
断层下覆界面反射能量强
炮点
T=2000ms
炮点 T=2300ms
炮点位于11km处的单炮记录
?21?4?1?6?1?m?2m?1222m246333m24lllol622m32m?m4?m6??m?m?2m?m?2?c1m??m??c2?m??c3???m??cm??m??1??0????0???m????0??35?1?133335?555?135?mm?m2n?12n?12n?1?35?1llloln?2n?1?c1??1??n???3?2n?1??c2??0?n?5??0?2n?1?c??3???m??m??m?n?2n?1????2n?10???cn????ox2ox数值频散试验dxdz10mdt1ms10高阶差分为何会消除数值频散
基于波动方程有限差分数值模拟及FCT消频散分析
中 国 锰 业
２０１８年 ４月
ＣＨＩＮＡ′ＳＭＡＮＧＡＮＥＳＥＩＮＤＵＳＴＲＹ
Ｖｏｌ．３６Ｎｏ．２ Ａｐｒ．２０１８
基于波动方程有限差分数值模拟及 ＦＣＴ消频散分析
赵 威，李伟波，桂志先，周 游，于晓东
（长江大学 油气资源与勘探技术教育部重点实验室，湖北 武汉 ４３０１００）
Δｘ＝Δｚ＝ｈ
（４）
１．２ 一阶速度—应力弹性波方程
在二维各向同性介质中，不考虑体力的影响下，
弹性波方程可以化为一阶速度—应力的形式：
!"
$ ""
% "&
#
"
&
!& '
$ "& "
% && &
""
'
$()*
!" "
高，占用内存小等特点，但由于其本身时间和空间都 采用差分 离 散 形 式，往 往 会 出 现 数 值 频 散 问 题［２］。 怎样减少频散，提高模拟精度始终是有限差分数值 模拟的关键问题。为此，学者们做了大量工作，当信 号在一个波长范围内，具有多个离散节点是，频散现 象基本消失。一般随着网格间距增大，频散问题越 严重；反之若减小网格间距相应会增加网格点数，将 会大大增加计算量，造成计算时间过长。既能控制 模拟精度，又能节省计算时间的方法显得尤为重要。 Ｂｏｒｉｓ和 Ｂｏｏｋ等人 首 先 提 出 了 通 量 校 正 传 输 的 方 法，它最早起源于求解流体力学连续方程中，随后才 开始被运用于声波方程数值模拟的求解；国内杨顶 辉等［３］将其 应 用 于 弹 性 波 动 方 程 正 演 数 值 模 拟 频 散处理，取得一定效果。为了更好的对比分析未加 ＦＣＴ数值模拟和加了 ＦＣＴ方法的区别，本文基于声 波和弹性波两方面，借助于交错网格的模拟优势，对 各向同性介质波场进行了有限差分数值模拟，并对 其相应的波场快照和计算效率进行对比分析。
基于GID有限元前处理的波动方程数值模拟
基于GID有限元前处理的波动方程数值模拟刘静;文山师;黄晶晶【摘要】在地震波数值模拟计算过程中,缺乏简单易行的有限元前处理方法,使得复杂构造模型较难建立和分析.本文以二维声波方程为例结合GID软件,网格剖分部分采用三角形单元模拟速度界面,把单元内的场和波速均看作单元上的线性函数;GID 软件可以方便地进行网格剖分和设置网格控制节点,通过编写用户自定义”问题类型”,建立并输出已有的有限元计算程序的初始模型.将GID软件前处理与有限元计算程序整合,提高了方法的效率,简单易行.【期刊名称】《工程地球物理学报》【年(卷),期】2014(011)002【总页数】7页(P243-249)【关键词】数值模拟;有限元;GID;声波方程;三角形单元【作者】刘静;文山师;黄晶晶【作者单位】山西省煤炭地质115勘查院,山西大同037003;中石化西北油田分公司勘探开发研究院,新疆乌鲁木齐830011;中石化石油工程地球物理有限公司河南分公司,河南南阳473000【正文语种】中文【中图分类】P631.41 引言地震波场的数值模拟技术是在已知地下介质结构和参数的情况下，利用理论计算的方法研究地震波在地下介质中的传播规律，合成地震记录的一种技术。
地震勘探中的数值模拟方法主要以射线理论和波动方程理论为基础，有射线追踪法，柯西霍夫积分法，有限元法，有限差分法和伪谱法[1～6]。
有限差分法直接用差分代替微分，因其方法简单、精度高，在地震模拟中而得到了广泛的研究和应用。
但其固有缺陷是不能准确模拟具有复杂几何形态的物性界面，有限元法则是求解原问题等价泛函的变分或原问题的等效积分方程的弱解（当等价泛函不存在时），因而能够适应较有限差分更为剧烈的物性变化，加之种类繁多的插值形函数，使其能够模拟很复杂的几何界面。
有限元法的主要缺点是计算和存储量都很大，效率相对较低。
建立有限元分析模型比较复杂且存在困难，因此可以用一些成形的软件作为有限元网格剖分的工具，建立并输出可用于已有有限元计算程序的初始模型，将大大提高方法的效率[7]。
声波方程有限差分数值模拟的变网格步长算法
声波方程有限差分数值模拟的变网格步长算法声波方程有限差分数值模拟是一种常用的声波传播模拟方法，可以在计算机上通过数值计算求解声波传播的过程。
在进行这种数值模拟时，常常需要选择合适的网格步长，以保证计算结果的准确性和计算效率。
本文将介绍一种变网格步长算法，用于优化声波方程有限差分数值模拟的计算。
声波方程可以用下面的形式表示：∂^2p/∂t^2=c^2∇^2p其中p是声场变量，t是时间，c是声速，∇^2是Laplace算子。
为了将声波方程用有限差分方法进行离散化计算，我们需要将空间和时间分别离散化。
首先，将空间离散化为网格，在每个网格点上计算声场的值。
其次，将时间离散化为离散的时间步长，通过迭代计算不同时间步长上的声场分布。
为了保证计算结果的准确性，网格步长应当满足Nyquist采样定理的要求。
即网格步长应小于声波的最小波长的一半。
根据声波方程的性质，我们可以通过声速和最高频率来估计声波的最小波长。
然后，我们可以根据最小波长来选择合适的网格步长。
然而，在实际的声波传播计算中，声场的变化往往不是均匀的。
有些区域的声场变化较大，而其他区域的声场变化较小。
如果我们在整个计算区域都采用较小的网格步长，将会造成计算资源的浪费。
因此，需要一种方法能够根据声场的变化情况来自适应地调整网格步长。
变网格步长算法就是一种能够根据声场变化情况自动调整网格步长的算法。
其基本思想是根据声场在不同网格上的变化率来决定每个网格上的网格步长。
具体的算法步骤如下：1.初始化：选择一个合适的初始网格步长。
通常可以选择根据声波的最小波长来确定。
2.计算网格步长：在每个时间步长上，对于每个网格点，计算其周围网格点上的声场变化率。
常用的方法是计算声场在三个相邻时间步长上的差分值，然后取绝对值并求平均。
根据声场变化率，调整当前网格点上的网格步长。
变化率大的网格点应该有更小的网格步长，而变化率小的网格点则可以有更大的网格步长。
3.更新声场：根据调整后的网格步长，更新所有网格点上的声场值。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
·８·
油气地球物理
２００７年 ７ 月
为压制频散， 提高模拟精度， 只有减小采样间隔， 如 用常网格 ４ｍ 步长算法进行模拟， 内存需求将会扩 大 １ 倍， 以上述（ １６００ｍ ×１６００ｍ ） 的模拟区域为例， 计算时间将会增大 １．２ 倍左右。用文中所述变步长 算法进行模拟， 低速地表用 ４ｍ 小步长， 地表以下用 ８ｍ 大步长， 变网格步长算法的内存需求仅是常网格
" $ ２
!
Ｐ
２
＝－
２
Ｌ
!ｔ
Ｐ !
－
１
２
２
Ｌ Ｐ＋ＰＬ
１
＋ｓ
!
!
（ ３）
% $ ２
－
２
Ｌ＝
!ｃ
２
２
２
!
２＋
!
２
!ｘ !ｚ
（ ４）
收稿日期： ２００７－０３－２３； 修订日期： ２００７－０４－２７ 作者简介： 李胜军， 男， 在读硕士研究生， 研究方向为地震 波 传 播 理 论 。 联 系 电 话 ：（ ０５４６） ８３９２０５５， Ｅ－ｍａｉｌ： ｈｄｐｕｌｉｓ＠１２６．ｃｏｍ， 通 讯 地 址 ： （ ２５７０６１） 中国石油大学（ 华东） 地球资信与信息学院。 ＊ 中国石油大学（ 华东） 研究生创新基金资助， 编号： Ｓ２００６—０６。
０ ０．０
８００
１６００ （ ｍ）
０
０．０
８００
１６００ （ ｍ）
０．５
０．５
１．０
１．０
１．５
（ ｓ）
（ ａ） 均匀网格
１．５
（ ｓ）
（ ｂ） 变网格
图 ７ 常网格与变网格得到的炮集记录
４ 结论
上述分析与理论模型试算结果表明， 变网格步 长差分算法是一种高效的正演模拟算法。它能较好 地实现对低速地表、高速层夹层、复杂界面的数值模 拟， 成功地解决小网格步长差分算法内存需求量大、 计算效率低等问题。在特殊地区用增大（ 减小） 步长 的模拟方法减小了内存需求， 提高了计算效率； 在特 定位置减小步长， 解决了对超薄层地质体、倾斜界面 的数值模拟效果不理想的问题： 解决了用常规高阶 有限差分或傅立叶变换法模拟存在的问题。
１ 时间积分
均 匀 介 质 中 的 二 维 声 波 方 程 可 用 下 式 表 示［２］
２
!
Ｐ
２
２
＝－ Ｌ Ｐ＋ｓ
（ １）
!ｔ
" # ２
２
２２
－ Ｌ ＝ｃ
!２＋
!
２（ ｘ， ｚ， ｔ） ， 代表压 力项； ｃ＝ｃ（ ｘ， ｚ） ， 代表速 度； ｓ＝ｓ（ ｓ， ｚ） ， 代表震源函数； Ｌ２ 为差分算子。 在 密 度 !＝!（ ｘ， ｚ） 变 化 的 情 况 下 ， 常 用 的 是 Ｖ ｉｄａｌｅ 给 出 的 公 式［５］
$# % Ｐ !ｘ， ｚ， ｔ "＝ ｔ ｓｉｎ !!Ｌ "ｈ !ｔ－ ! "ｄ ! ｇ !ｘ， ｚ "（ ５） ０Ｌ 该常规解通过契比雪夫展开近似为
$ ! "% &Ｎ／２
Ｐ !ｘ，
ｚ，
ｔ
"
ｋ
＝
ｂ２ｋ＋１
０
Ｒ ｉＬ!
Ｑ２ｋ＋１
ｉＬ! Ｒ
ｇ !ｘ， ｚ " （ ６）
式中： ! 为变量； Ｌ! ２ 是空间差分算子 Ｌ２ 的数值 逼近；
２００７ 年 ７ 月
油气地球物理 ＰＥＴＲＯＬＥＵＭ ＧＥＯＰＨＹＳＩＣＳ
第５卷 第３期
波动方程的变步长有限差分数值模拟 ＊
李胜军 １， ２） 孙成禹 １） 张玉华 １） 倪长宽 １） １） 中国石油大学地球资源与信息学院； ２） 中石油勘探开发研究院西北分院
摘要： 有限差分算法是常用的正演模拟方法之一， 其包含的地震信息丰富， 且实现简单。传统的有限差分方法通常都 采用均匀网格步长， 在对含低速 ／ 高速介质、薄层 ／ 厚层介质的模型进行波场模拟时往往缺乏稳定性。文章介绍了一 种可以有效解决上述问题的变网格算法， 对常规有限差分法与变网格差分算法在内存需求、计算速率等方面的差别 进行了比较， 对变网格差分算法中的边界条件、时间积分的快速展开算法作了阐述， 进而总结了变网格算法的优点。 关键词： 变步长； 边界条件； 计算时间； 快速展开法； 数值模拟
·６·
油气地球物理
２００７年 ７ 月
波动方程 的 时 间 积 分 可 通 过 Ｋ ｏｓｌｏｆｆ等 提 出 的 快速展开法（ Ｒ Ｅ Ｍ ） 来计算［６］。为 简便起见 ， 仅对基 本方程（ １） 的快速展开法做一些介绍。对变密度的 情况可用类似的方法进行推导。
对震源项 ｓ（ ｘ， ｚ， ｔ） ＝ｇ（ ｘ， ｚ） ｈ（ ｔ） ， 边界条件为 吸收边界时， 式（ １） 的常规解为
４ｍ 步长算法的 ６５％ 左右， 计算时间也缩短 ４５％ 左 右。与常网格 ８ｍ 步长算法相比， 计算时间是常网格 ８ｍ 步长的 １．２ 倍， 当然， 内存需求也相应有所扩大。 图 ７（ ｂ） 是变网格步长算法的模拟结果（ 为了对比 明显， 两图使用了相同的增益） ， 可以看出， 频散得 到了较好的压制， 分辨率也明显提高。
在图 ２ 显示的网格中， 可使用这种变网格算法 来改变空间分辨率。在细网格区域不需要用傅立叶 变换法或高阶有限差分算子， 四阶或六阶精度的算 法就可以满足精度的要求。
ｘ
如果 Ｚ 方向上的步长通过其他参数来调整， 类 似的方程可被推导。一般通过函数 ｍ（ ｋ） 来产生， 有
Ｎ
% ２
Ｄｚ
!Ｐ
" ｉ， ｊ
１ #２!ｚ
$２ ｋ
＝
!
１
# $ Ｐ － ２Ｐ ＋Ｐ ｋ ｉ， ｊ＋２# ｋ＋ｎ $ ｉ， ｊ ｉ， ｊ－ ２# ｋ＋ｎ $ （ １１）
为保证这种差分算法的稳定性， 压力必须内插 在水平方向两倍步长边界上的一些点上。对参数长 度 ２Ｎ＝ １０ 的那些点在图 ５ 中已列出。内插计算时 可使用特别精确且没有错误发生的三角内插算法。
$ｘ 在 Ｘ 方向上是连续变化的， 用常规的有限差分 算法就可计算。但是在 Ｚ 方向上， %ｚ 是不连续的， 基 于连续变化的傅立叶变换法不再适用， 所以 Ｚ 方向上的导数要用高阶有限差分法近似。
ｘ
ｚ
ｚ０
图 １ 低速地表模型差分网格分布示意图 ｘ
ｚ
图 ２ 低速夹层模型差分网格分布示意图 ｘ １５００ｍ ／ｓ
图 ３ 给出了简单的模型： 用固定小网格步长 !ｘ 和 "ｚ 采样， 可以很好地满足低速层的采样需求， 但 是这样造成了对高速层的过采样。为避免这种过采 样， 可以从某一深度 ｚ０（ 图 １） 开始采用双倍的步长 来采样。在修改的网格上， Ｘ 方向的导数用傅立叶法 或有限差分法计算， #ｘ 通过采样定理来确定。由于
列表（ 黑圆点位置属于精细网格循环； 方框点属于 双倍的步长循环点； 五角星位置点的导数被计算） 。 在每一列用不同的差分模式来计算在不同深度的二
阶导数。针对十阶精度有限差分相应的 ｍ（ ｋ） 的值列 于表 １， 其他阶精度的 ｍ（ ｋ） 值可用相应的方法计算。
ｚ０ ｚ
图 ５ 必须内插压力的网格点
３ 算例分析
Ｑｋ 指修改了的契比雪夫多项式； 系数 ｂｋ 为
# ｂｋ＝
ｔ
Ｊ２ｋ＋１
!!Ｒ
" ｈ
!ｔ－
!
"ｄ!
０Ｒ
（ ７）
式中： Ｊｋ 为 ｋ 阶精度的贝赛尔函数。 如果空间参数 － Ｌ! ２ 的特征值很靠近实轴、Ｒ２ 比
最大特征值大， 当 Ｎ＞Ｒ·ｔ 时， （ ６） 、（ ７） 式按指数规 律收敛。参数 － Ｌ! ２的最大特征值由 Ｘ 方向上的差分
＝
１ #!ｚ $２ ｋ
＝
!
１
#Ｐ － ｋ ｉ， ｊ＋ｍ（ｋ）
$ ２Ｐｉ， ｊ ＋Ｐｉ， ｊ－ ｍ（ｋ）
（ １２）
差分因子
!
必须根据函数
ｋ
ｍ（
ｋ）
计算， 这可通过
Ｆｏｒｎｂｅｒｇ（ １９８８ａ） 提出的一种算法来实现［９］。图 ４ 表
示了方程（ １２） 在双倍步长区域差分算子 ２Ｎ＝ １０的
情况下有限差分参数的变化情况。显示了 ６ 个网格
参数的最大特征值 Ａｘ 和 Ｚ 方向上的相应 特征值 Ａｚ 的和来决定。这里， Ｒ２ 可由（ ８） 式计算
! " ２ ２
Ｒ ＝ｃ
Ａｘ !!ｘ
"２ ＋
Ａｚ !!ｚ "２
ｍａｘ
（ ８）
２ 空间导数计算
通常情况下， 勘探区域的地表较为复杂且介质 中波的传播速度较低。为较好地压制数值频散并提 高模拟结果的分辨率， 就要在整个计算区域采用小 网格步长来做差分数值模拟。这样就造成了对高速 区域的过采样， 浪费了计算机资源。为减小内存需 求和节省计算时间， 可以用图 １ 所示方法进行模 拟。对于地层中有低速层的情形可以采用图 ２ 所示 方法来达到提高分辨率并且计算量增加不大的差 分算法。
为验证上述算法的正确性、有效性， 设计了图 ６ 所示的低速地表地质模型。数值模拟时， 震源设置在 （ ８００ｍ ， ２４ｍ ） 处。
ａｂ
ｃ ｄｅｆ
０
８００
１６００ （ ｍ）
０
１０００ｍ ／ｓ
ｚ０
图 ４ 过渡带的差分格式
表 １ 图 ４ 显示的有限差分参数 ｍ（ ｋ） 值
ｍ（ ｋ） ａ
ｂ
ｃ
ｄ
ｅ
ｆ
ｍ（ ０） ０
１
"
（ ９）
对应的双步长网格定义如下
２
Ｎ
& ２
Ｄｚ
$Ｐ
% ｉ， ｊ
＝
"
Ｐ