2007射线追踪与波动方程正演模拟方法对比研究
４７
科技资讯
科技资讯 ＳＣＩＥＮＣＥ ＆ ＴＥＣＨＮＯＬＯＧＹ ＩＮＦＯＲＭＡＴＩＯＮ２００７ ＮＯ．１２
ＳＣＩＥＮＣＥ ＆ ＴＥＣＨＮＯＬＯＧＹ ＩＮＦＯＲＭＡＴＩＯＮ
工 业 技 术
地震正演模拟作为反演解释的反过程，是验证解释成果的有效手段，进行必要可靠的正演模拟可以有效的监控反演解释。地震学一般可以分为几何地震学和物理地震学，在几何地震学中进行的正演模拟方法就是我们通常所说的射线追踪法，射线追踪法是在合成记录时用地震子波和界面或地质体的反射系数进行反褶积运算，即。运算的最大特点是说明了地震波传播的运动学特征。而在物理地震学中应用波动方程法合成的地震记录是通过求解波动方程的数值解来模拟地震波场的。在波动方程合成的地震记录中不单保持了地震波传播运动学特征，还说明了地震波传播的动力学特征。本文将分别用射线追踪和波动方程的方法合成地震记录。
１ 基于射线追踪的合成地震响应
射线追踪法的主要理论基础是，在高频近
射线追踪与波动方程正演模拟方法对比研究
王志美 畅永刚
（长江大学油气资源与勘探技术教育部重点实验室 湖北荆州 ４３４０２３）
摘 要：地震学一般可以分为几何地震学和物理地震学，几何地震学中进行正演模拟方法就是射线追踪法，射追踪法是在合成记录时用地震子波和界面或地质体的反射系数进行反褶积运算，即。运算的最大特点就是说明了地震波传播的运动学特征。而在物理地震学中的波动方程法合成的地震记录是通过求解波动方程的数值解来模拟地震波场。在波动方程合成的地震记录中不单保持了地震波传播
运动学特征外，还说明了地震波传播的动力学特征。本文将分别用射线追踪和波动方程的方法合成地震记录。关键词：射线追踪 波动方程 正演模拟
中图分类号：Ｐ３１５文献标识码：
Ａ
文章编号：１６７２－３７９１（２００７）０４（ｃ）－００４７－００
图１
射线追踪正演模拟（１）
图
２
逐段迭长示意图
图
３ 射线追踪正演模拟（２）
图４ 波动方程正演模拟结果
似条件下，地震波的主能量沿射线轨迹传播。基于这种认识，运用惠更斯原理和费马原理来重建射线路径，并利用程函方程来计算射线的旅行时。在旅行时计算中应用有限差分等方法，以获得快速的解。射线法的主要优点是概念明确，显示直观，运算方便，适应性强；其缺陷是应用有一定限制条件，计算结果在一定程度上是近似的，对于复杂构造进行两点三维射线追踪往往比较麻烦。为了计算波沿射线的旅行时和波的传播路径，叙述如下。
如图１所示，首先给出连接Ｓ（激发点）和Ｒ（接收点）之间的初始射线路径射线的振幅变化，首先必须知道地震波在实际地层中传播的射线路径。
由于地震波在整条路径上满足同一个射线参数，因此射线路径上任意连续三点也将满足同一个参数，而三点间的射线表现形式为Ｓｎｅｌｌ定律。按照Ｓｎｅｌｌ 定律，可导出一个求
取中间点的一阶近似公式。当前后两点位于界面两边时，中间点为透射点，所求路径为透射路径；当前后两点位于界面的同一边时，中间点为反射点，所求路径为反射路径。为此，可以从任一端点出发，连续地选取三点，通过一阶近似公式进行逐段迭代取中间点，利用新求出的点代替原来的点，然后以一点的跨跃作为步长，顺序地逐段迭代下去，直到另一端点。这样，新计算出的中间点和两个端点就构成了一次迭代射线路径，如图２中所示。如果整条射线路径上校正量的范数之和满足一定的精度要求，则认为射线追踪过程结束，否则从追踪出的射线路径开始，继续重复上述过程，直到满足精度要求为止。最后一次追踪到的中间点和两个端点，构成整条射线路径。图３基于多层倾斜界面模型通过射线追踪正演模拟地震响应。从模拟结果可以直观的看出基于几何地震学的原理正演模拟结果只能反映地震波的几何传播路径。在实际的工程设计中通过正演模拟可以在地表确定地下观测范围，节约设备提高工程效率，但不能反映
物理地震学中的地震属性，例如振幅，频率和相位等。更不能反映地震波的动力学特征。
２ 波动方程的合成地震响应
２．１ 波动方程的建立
非均匀介质的声波方程：
（１） （２）
可由对连续介质方程（１）式的两端对时间求导，并利用欧拉方程推得：
（３）
其中：Ｐ是波数，Ｖ是质点振动的速度向量，ρ是密度，ｃ是波速，ρ和ｃ是随着空间参数χ和ｚ变化的，这里ρ给定为常数，只有ｃ
是地质模型的控制参数。χ和Ｚ分别是在地面水平距离和深度。这样（３）式就可以变为：
（４）
其中：ｃ＝ν （χ，ｚ）；（４）式即是所求的弹性波动方程。
２．２ 数值计算及稳定性
求解弹性波动方程的方法有多种，付立叶变换法是对弹性波动方程的波场进行付立叶变换，优点是运算速度快。克希霍夫积分法是基于均匀模型，利用格林函数公式计算曲面积分，求出空间波场值，但这种方法不能适应
４８科技资讯
科技资讯 ＳＣＩＥＮＣＥ ＆ ＴＥＣＨＮＯＬＯＧＹ ＩＮＦＯＲＭＡＴＩＯＮ
２００７ ＮＯ．１２
ＳＣＩＥＮＣＥ ＆ ＴＥＣＨＮＯＬＯＧＹ ＩＮＦＯＲＭＡＴＩＯＮ
工 业 技 术
１ 工程概况：
惠阳区柠檬酸厂位于广汕公路８号，该厂房由上海轻工业设计院设计，建于１９９０～１９９１年，框架结构、采用天然基础、层高５ｍ、总高度２５ｍ、长４７．７５ｍ、宽１４．６ｍ、建筑面积约３４００ｍ２、共五层。由于该厂房生产的是有重度腐蚀性的柠檬酸，建成后使用不久，柱混凝土的表层发生严重破坏，混凝土柱表层出现裂缝，部分钢筋表层已有被腐蚀的现象。经业主委托广东中科华大工程技术检测有限公司的检测评估结论：一层柱强度普通
偏低和出现了大裂缝，缝宽１－４ｍｍ，其中柱网×、×、×、×和×共五条柱截面出现裂缝，×轴柱柱头表层混凝土已经蚀去５０ｍｍ，呈松散状，仅用手指便可以取下，较严重的是Ｚ２底层柱，混凝土强度偏低，原设计Ｃ２０的混凝土，检测强度仅有Ｃ１５。其中柱编号为：Ｚ１、Ｚ２、Ｚ３、Ｚ４、Ｚ５ 、柱截面尺寸ｂ×ｈ＝４００×６００，具体分布详见附图（一）。
２ 设计依据：
１）根据《混凝土结构加固技术规范》
柱外包钢的加固设计应用实例
孔越东
（惠州市正富工程建设监理有限公司 广东惠州 ５１６００１）
摘 要：通过实例分析解决柱外包钢在加固技术中的应用关键词：结构设计 外包钢 加固中图分类号： ＴＵ７４文献标识码：Ａ
文章编号：１６７２－３７９１（２００７）０４（ｃ）－００４８－０２
ＣＥＣＳ１４６：２００３，设计图纸和该根据工程检测
报告编号ＢＯｂＬＯＪＧ０３３，本工程采用加大截面加固法、外粘钢加固法等，其工作程序如下：
可靠性稳定→加固方案→加固设计→施工→验收
２）材料：外包钢采用Ｑ２３５材料Ｌ８０×８０×５的Ｂ型角钢，Ａｓ＝７９１ｍｍ２、ｆｙ＝２１５Ｎ／ｍｍ２加固箍筋用扁铁４０×４，外包混凝土用Ｃ２５（ｆＣ＝１１．９Ｎ／ｍｍ２），外包钢加固后的尺寸ｂ×ｈ＝５００×７００，角钢与扁钢的连接采用焊接，ｈｆ＝５ｍｍ焊缝饱满，焊条Ｅ４３、Ｅα＝２．１×１０５Ｎ／ｍｍ２，粘钢采用改性环氧树脂胶粘剂。
３）单根柱加固前的承载力计算为：Ｎ１＝４．５×３．５×２０×５＝１５７５ＫＮ
Ｎ２＝６×３．５×２０×５＝２１００ＫＮＮ３＝６×３．５×２０×５＝２１００ＫＮＮ４＝６×３．５×２０×５＝２１００ＫＮＮ５＝６×３．５×２０×５＝２１００ＫＮ
原混凝土强度等级按Ｃ２０的标准取值ｆｙ＝３００Ｎ／ｍｍ２、Ａｓ＝８×４９０＝３９２０ｍｍ２。外包混凝土强度Ｃ２５、ｆＣ＝１１．９Ｎ／ｍｍ２，α＝０．８
４）根据混凝土结构加固技术规范，当用加大截面加固钢筋混凝土轴心受压构件时，其正截面承载力应按下列公式计算：
Ｎ＜Ψ〔ｆｃｏＡｃｏ＋ｆ／ｙｏＡ／ｓｏ＋（ｆｃＡｃ＋ｆ／ｙＡ／ｓ）〕（４．１．２）
式中Ｎ
为构件加固后的轴向力设计值’Ψ为构件的稳定系数，
以加固后截面为
准，
按现行国家标准《混凝土结构设计规范》速度变化情况下的介质。有限元素法能对复杂地质构造模型进行地震波场数值模拟，通常计算速度较慢，有限差分法是一种求解弹性波
动方程的有效方法，而且地质模型的复杂程度并不影响运算速度。２．３ 波动方的差分解法
取网格：
其中ｈ＞０是χ方向的步长，ι＞０是ｔ的步长，ｈ＝ι／Ｎ，ι＝ｔ／ｊ，空间域和时间域的离散化处理，如果我们选择空间在横向和深度的二维波动方程，求波函数Ｐ（ｘ，ｚ，ｔ）的数值解，取空间采样间隔ｄｘ＝ｄｚ＝１０ｍ，时间采样间隔ｄｔ＝ｌｍｓ，在空间采样点数在横向ｘ上为ｍ，在深度ｚ上为ｎ，在时间ｔ上为ｊ，则ｘ＝ｍｄｘ，ｚ＝ｎｄｚ，
ｔ＝ｊｄｔ；；这样（４）式就可以变为：
（５）
其中：
Ａ＝Ｐ（ｍ＋１，ｎ，ι）－２Ｐ（Ｍ，ｎ，ι）＋Ｐ（ｍ－１，ｎ，ι）
Ｂ＝Ｐ（ｍ，ｎ＋１，ι）－２Ｐ（ｍ，ｎ，ι）＋Ｐ（ｍ－１，ｎ，ι）
图４分别是处理边界效应和没有处理边界效应的基于物理地震学的波动方程正演模拟结果，在处理边界问题中采用了多道接收，取其中近激发点道记录，这样可以消除边界效应。波动方程模拟结果明显的说明了波动的客观时距关系同时也可以表征能量等物理地震学的意义。
３ 结语
基于射线追踪的正演模拟，形象的说明了地震波传播的运动学特征。用波动方程进行的数值正演模拟可以把地下构造在每个时刻振动都记录下来构造（介质）的振动情况。数值模拟又分为波动方程的解析解和数值解（逼近值）。特别复杂地质模型解析解很难获得，这只能靠波动方程的数值解才可以说明构造
对地震波传播的影响。其实解析解只是数字解的检验。
参考文献
［１］ ＣｅｒｖｅｎｙＶ，Ｍｏｌｏｔｋｏｖ ＩＡ，Ｐｓｅｎｃｉｋ Ｉ１ 地
震学中的射线方法［Ｍ］．刘福田译１北京：地质出版社，１９８６．
［２］ 何樵登．地震勘探原理和方法［Ｍ］．北京：地
质出版社，１９８６．
［３］ 朱介寿．地震学中的射线方法［Ｍ］．北京： 地
震出版社，１９８８．
［４］ 高尔根，徐果明．任意界面下的整体迭代射
线追踪方法的研究［Ｊ］．学学报，２００２，２７（３）：２８２～２８５．
［５］ 高尔根，徐果明，赵炎．一种任意界面的逐
段迭代射线追踪方法［Ｊ］．石油地球物理勘探，１９９８，３３（１）：５４～６０．
［６］ 朱海波．基于射线追踪法的地震数据采集
参数优化［Ｊ］．江汉石油学院学报，２００４，２６（１）５７～５９．
一维热传导方程的差分格式
《微分方程数值解》 课程论文 学生姓名1：许慧卿学号：20144329 学生姓名2：向裕学号：20144327学生姓名3：邱文林学号：20144349学生姓名4：高俊学号：20144305学生姓名5：赵禹恒学号：20144359学生姓名6：刘志刚学号： 20144346 学院：理学院 专业：14级信息与计算科学 指导教师：陈红斌 2017年6 月25日
《偏微分方程数值解》课程论文 《一维热传导方程的差分格式》论文 一、《微分方程数值解》课程论文的格式 1）引言：介绍研究问题的意义和现状 2）格式：给出数值格式 3）截断误差：给出数值格式的截断误差 4）数值例子：按所给数值格式给出数值例子 5）参考文献：论文所涉及的文献和教材 二、《微分方程数值解》课程论文的评分标准 1）文献综述：10分； 2）课题研究方案可行性：10分； 3）数值格式：20分； 4）数值格式的算法、流程图：10分； 5）数值格式的程序：10分； 6）论文撰写的条理性和完整性：10分； 7）论文工作量的大小及课题的难度：10分； 8）课程设计态度：10分； 9）独立性和创新性：10分。 评阅人： - 2 -
一维热传导方程的差分格式 1 引言 考虑如下一维非齐次热传导方程Dirichlet 初边值问题 22(,),u u a f x t t x ??=+?? ,c x d << 0,t T <≤ (1.1) (,0)(),u x x ?= ,c x d ≤≤ (1.2) (,)(),u c t t α= (,)(),u d t t β= 0t T <≤ (1.3) 的有限差分方法, 其中a 为正常数,(,),(),(), ()f x t x t t ?αβ为已知常数, ()(0),c ?α= ()(0).d ?β= 称(1.2)为初值条件, (1.3)为边值条件. 本文将给出(1.1) (1.3)的向前Euler 格式, 向后Euler 格式和Crank Nicolson -格式, 并给出其截断误差和数值例子. 经对比发现, Crank Nicolson -格式误差最小, 向前 Euler 格式次之, 向后Euler 格式误差最大. 2 差分格式的建立 2.1 向前Euler 格式 将区间[,]c d 作M 等分, 将[]0,T 作N 等分, 并记 ()/h d c M =-, /T N τ=, j x c jh =+,0j M ≤≤, k t k τ=,0k N ≤≤. 分别称h 和τ为空间步长和时间步长.用 两组平行直线 j x x =, 0j M ≤≤, k t t =, 0k N ≤≤ 将Ω分割成矩形网格.记{} |0h j x j M Ω=≤≤, {}|0k t k N τΩ=≤≤, h h ττΩ=Ω?Ω. 称() ,j k x t 为结点[1] . 定义h τΩ上的网格函数 {}|0,0k j U j M k N Ω=≤≤≤≤, 其中() ,k j j k U u x t =. 在结点() ,j k x t 处考虑方程(1.1),有
【文献综述】热传导方程差分格式的收敛性和稳定性
文献综述 信息与计算科学 热传导方程差分格式的收敛性和稳定性在实际研究物理问题过程中, 往往能给出问题相应的数学表达式, 但是由于实际物理问题的复杂性, 它的解却一般不容易求出. 由此计算物理应运而生, 计算物理是以计算机为工具, 应用数学的方法解决物理问题的一门应用性学科, 是物理、数学和计算机三者结合的交叉性学科. 它产生于二战期间美国对核武器的研究, 伴随着计算机的发展而发展. 计算物理的目的不仅仅是计算, 而是要通过计算来解释和发现新的物理规律. 这一点它与传统的实验物理和理论物理并无差别, 所不同的只是使用的工具和方法. 计算物理早已与实验物理和理论物理形成三足鼎立之势, 甚至有人提出它将成为现代物理大厦的“栋梁”. 在一个物理问题中一个数值解往往比一个式子更直观, 更有价值. 在实际求解方程时, 除了一些特殊的情况下可以方便地求得其精确解外, 在一般情况下, 当方程或定解条件具有比较复杂的形式, 或求解区域具有比较复杂的形状时, 往往求不到, 或不易求到其精确解. 这就需要我们去寻找方程的近似解, 特别是数值近似解, 简称数值解. 这里主要研究的是热传导方程. 有限差分法是微分方程和积分微分方程数值解的方法. 其基本思想是把连续的定解区域用有限个离散点构成的网格来代替, 这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似, 积分用积分和来近似, 于是原微分方程和定解条件就近似地代之以代数方程组, 即有限差分方程组, 解此方程组就可以得到原问题在离散点上的近似解. 然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解. 热传导的差分法是求解热传导方程的重要方法之一. 对于差分格式的的求解, 我们首先要关注差分格式的收敛性和稳定性. 对于一个微分方程建立的各种差分格式, 为了有实用意义, 一个基本要求是它们能够任意逼近微分方程, 即相容性要求. 一个差分格式是否有用, 就要看差分方程的精确解能否任意逼近微分方程的解, 即收敛性的概念. 此外, 还有一个重要的概念必须考虑, 即差分格式的稳定性. 因为差分格式的计
叠加地震记录的相移波动方程正演模拟数值模拟实验共22页
《地震数值模拟》实验报告 一、实验题目 叠加地震记录的相移波动方程正演模拟
二、实验目的 1.掌握各向同性介质任意构造、水平层状速度结构地质模型的相移波动方程正演模拟基本理论 2.实现方法与程序编制 3.由正演记录初步分析地震信号的分辨率。 三、实验原理 1、地震波传播的波动方程 设（x,z）为空间坐标，t为时间，地震波传播速度为v(x,z)，则二位介质中任意位置、任意时刻的地震波场为p(z,x,t)：压缩波——纵波。则二维各向同性均匀介质中地震波传播的遵循声波方程为 2、傅里叶变换的微分性质 p(t)与其傅里叶变换的P(w)的关系： 3、地震波传播的相移外推公式 令速度v不随x变化，只随z变化，则利用傅里叶变换微分性质把波动方程(变换到频率-波数域，得： 4、初始条件和边界条件 按照爆炸界面理论，反射界面震源在t=0时刻同时起爆，此时刻的波场就是震源。根据不同情况，可直接使用反射系数脉冲或子波作震源。如果直接使用反射系数作震源脉冲，则初始条件可表示为： 5、边界处理
（1）边界反射问题 把实际无穷空间区域中求解波场的问题化为有穷区域求解时，左右两边使用零边界条件。物理上假设探区距Xmin与Xmax两个端点很远，在两个端点上收到的反射波很弱。但是，上述条件在实际中不能成立，造成零边界条件反而成为绝对阻止波通过的强反射面。在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。 （2）边界强反射的处理 镶边法、削波法、吸收边界都能有效消除边界强反射。 削波法就是在波场延拓过程中，没延拓一次，在其两侧均匀衰减到零，从而消除边界强反射的影响。假设横向总长度为NX，以两边Lx道吸波为例，有以下吸波公式： 四、实验内容
一维导热方程 有限差分法 matlab实现
第五次作业（前三题写在作业纸上） 一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.pdf 文件，热扩散系数α=const ， 22T T t x α??=?? 1. 用Tylaor 展开法推导出FTCS 格式的差分方程 2. 讨论该方程的相容性和稳定性，并说明稳定性要求对求解差分方程的影响。 3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。 4. 编写M 文件求解上述方程，并用适当的文字对程序做出说明。（部分由网络搜索得到，添加，修改后得到。） function rechuandaopde %以下所用数据，除了t 的范围我根据题目要求取到了20000，其余均从pdf 中得来 a=0.00001;%a 的取值 xspan=[0 1];%x 的取值范围 tspan=[0 20000];%t 的取值范围 ngrid=[100 10];%分割的份数，前面的是t 轴的，后面的是x 轴的 f=@(x)0;%初值 g1=@(t)100;%边界条件一 g2=@(t)100;%边界条件二 [T,x,t]=pdesolution(a,f,g1,g2,xspan,tspan,ngrid);%计算所调用的函数 [x,t]=meshgrid(x,t); mesh(x,t,T);%画图，并且把坐标轴名称改为x ，t ，T xlabel('x') ylabel('t') zlabel('T') T%输出温度矩阵 dt=tspan(2)/ngrid(1);%t 步长 h3000=3000/dt;
h9000=9000/dt; h15000=15000/dt;%3000,9000,15000下，温度分别在T矩阵的哪些行T3000=T(h3000,:) T9000=T(h9000,:) T15000=T(h15000,:)%输出三个时间下的温度分布 %不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面 %稳定性讨论,傅里叶级数法 dx=xspan(2)/ngrid(2);%x步长 sta=4*a*dt/(dx^2)*(sin(pi/2))^2; if sta>0,sta<2 fprintf('\n%s\n','有稳定性') else fprintf('\n%s\n','没有稳定性') error end %真实值计算 [xe,te,Te]=truesolution(a,f,g1,g2,xspan,tspan,ngrid); [xe,te]=meshgrid(xe,te); mesh(xe,te,Te);%画图，并且把坐标轴名称改为xe，te，Te xlabel('xe') ylabel('te') zlabel('Te') Te%输出温度矩阵 %误差计算 jmax=1/dx+1;%网格点数 [rms]=wuchajisuan(T,Te,jmax) rms%输出误差
波动方程正演模型的研究与应用
波动方程正演模型的研究与应用 郑鸿明* 娄 兵 蒋 立 （新疆油田公司勘探开发研究院地物所） 摘要野外采集的地震数据是经过大地滤波后的畸变信号，处理的地震剖面只是间接地反映了地下构造和地质体的特征，虽然目前有很多方法和手段可以分析并提取相关的地质信息，但由于处理对波场的改造和噪声的存在以及方法本身的多解性问题降低了识别地质信息的可靠性。处理中每一步对有效信息的影响有多大，对地震属性解释的影响有多大，没有一个定量的标准，只能凭经验和认识来定性地判断。正演模型在弹性波理论指导下，遵循严格的数学公式，可以最佳模拟地下各种情况。各种处理方法和不同的处理流程所得到的结果能否符合或最佳逼近波动方程建立的数学模型，正演模型是判断处理工作合理性的良好准则。 主题词地质模型波动方程正演模型地震响应模块测试 1 引 言 随着地震勘探的不断深入，地震勘探也由构造型油气藏勘探进入精细的岩性勘探阶段，要求地震勘探能够反映地下地质体岩性变化，以及识别含油、气、水的地震响应特征，分辨薄互层、低幅度构造的能力。地球物理学家们在长期的实践中已经研究开发了很多相关的技术，虽然理论上这些方法都能够成立，这些技术应用成功的实例也很多，但也不乏有失败的教训，往往产生多解性，或与钻探的结论不符。这里除了复杂地表和复杂地下构造形成的复杂地震波场而不满足建立在简单地质模型处理理论的因素外，与处理过程对地震波场的改造也有很大关系。从地震数据的采集到最终处理的地震剖面，整个过程是一个系统工程，地下地质结构、地质体的岩性变化以及含流体的性质，对处理人员来说是看不见、摸不着的“黑匣子”，我们所看到的只是经过大地滤波后产生畸变的地震波场，如何从这个畸变的地震波场中去伪存真、恢复真实的构造形态、提取储层的相关地震属性信息，这是岩性处理的最终目标。处理中的每一步环环相扣、相互影响、相互制约，而我们对处理中的每一步产生的中间结果所应达到的标准只是凭经验、感觉进行定性判定，加入了很多人为因素，这些因素或多或少影响着我们对解释成果的正确认识。另外，处理技术发展很快，相应的地震处理软件越来越多，应用这些模块之前对各模块所起的作用以及它们所产生的结果都需要有一个定量的认识，以及验证处理流程的合理性是当前迫切需要解决的问题。究竟什么样的结果满足岩性解释的要求、什么样的结果反映的是真正地下地质体的响应、什么样的处理方法满足保振幅处理和地震属性分析的应用等等一系列问题，这都是当前岩性处理中迫切需要解决的主要问题。它直接关联着处理成果的真伪及后续解释的可靠性，关联着勘探的投资风险。 随着计算机运算能力发展迅猛，特别是微机群的出现，为波动方程算法提供了硬件环境，开展此项技术的研究与应用已成为可能。此次模型的设计全面考虑了地表和地下的典型地质特征并将这些特征容入到模型中，真实模拟了实际地质结构。应用该地质模型正演叠前炮集的地震响应。 2 模型的建立 模型分物理模型和数学模型两种，目前的物理模型只能做非常简单的模拟，只有用数学模型才能模拟各种复杂的地质现象。20世纪70年代，美国哥伦比亚大学在郭宗汾
波动方程的变步长有限差分数值模拟
收稿日期：２００７－０３－２３；修订日期：２００７－０４－２７ 作者简介：李胜军，男，在读硕士研究生，研究方向为地震波传播理论。联系电话：（０５４６）８３９２０５５，Ｅ－ｍａｉｌ：ｈｄｐｕｌｉｓ＠１２６．ｃｏｍ，通讯地址：（２５７０６１）中国石油大学（华东）地球资信与信息学院。 ＊中国石油大学（华东）研究生创新基金资助，编号：Ｓ２００６—０６。 油气地球物理 ２００７年７月 ＰＥＴＲＯＬＥＵＭＧＥＯＰＨＹＳＩＣＳ 第５卷第３期 在地震资料采集、处理和解释中通常需要进行地震波场数值模拟：假设已知地下的地质情况，应用地震波运动学和动力学的基本原理，计算给定地质模型的地震响应。这种做法对正确认识地震波的运动学和动力学特征，以及准确分析油气藏的反射波场特征有着重要的指导意义。声波在介质中的正演模拟研究为我们精确模拟地震波在复杂介质中的传播提供了理论基础［１］。 傅立叶变换法和高阶有限差分法（ＦＤ）已成为计算声波方程空间导数的标准技术［２，３］。虽然常网格步长差分算法比较容易实现，但是它们对大部分模型都增大了不必要的计算量。例如，对存在浅层低速带的沉积盆地模型地面地震记录进行模拟时，由于低速地层阻抗小，地震波传入其中会引起较大的振幅和较长的延续时间（这与深层的高速层完全不同）。由于这些浅层低速层中地震波的波长较短、地层厚度较小，模拟时需要用小网格进行。这样，常网格步长算法就必须用小网格离散整个模型，从而增加了不必要的代价，如内存、计算量的增大。 因而，采用变网格算法将能改进有上覆低速层情况模拟结果的有效性（对地层中间有超薄夹层的情形，必须用精细网格覆盖才能精确的对地层进行模拟）。应用这种变网格算法既能实现对夹层的模拟，又能保障计算量不增加。因此这种通过函数实现在任意深度上网格步长变化的有限差分方法被 推广［４］。为了计算空间导数，在Ｘ方向用傅立叶变换法或有限差分算法，在Ｚ方向使用高阶有限差分方法。通过时间积分快速展开法（ＲＥＭ）来保障差分方法的计算精度［６］。这种差分技巧比二阶时间差分有较高的精确度且计算用时短。 １时间积分 均匀介质中的二维声波方程可用下式表示［２］ 式中：Ｐ＝Ｐ（ｘ，ｚ，ｔ），代表压力项；ｃ＝ｃ（ｘ，ｚ），代表速度；ｓ＝ｓ（ｓ，ｚ），代表震源函数；Ｌ２为差分算子。在密度!＝!（ｘ，ｚ）变化的情况下，常用的是Ｖｉｄａｌｅ给出的公式［５］ 波动方程的变步长有限差分数值模拟＊ 李胜军１，２） 孙成禹１） 张玉华１） 倪长宽１） １）中国石油大学地球资源与信息学院；２）中石油勘探开发研究院西北分院 摘要：有限差分算法是常用的正演模拟方法之一，其包含的地震信息丰富，且实现简单。传统的有限差分方法通常都采用均匀网格步长，在对含低速／高速介质、 薄层／厚层介质的模型进行波场模拟时往往缺乏稳定性。文章介绍了一种可以有效解决上述问题的变网格算法，对常规有限差分法与变网格差分算法在内存需求、计算速率等方面的差别进行了比较，对变网格差分算法中的边界条件、 时间积分的快速展开算法作了阐述，进而总结了变网格算法的优点。关键词：变步长；边界条件；计算时间；快速展开法；数值模拟 !２ Ｐ!ｔ２＝－Ｌ２Ｐ＋ｓ （１） （２） －Ｌ２ ＝ｃ ２ !２!ｘ２＋!２ !ｚ ２" # （３） （４） !２ Ｐ!ｔ ２＝－Ｌ２Ｐ!"$ －１!Ｌ２Ｐ＋ＰＬ２１!＋ｓ －Ｌ２ ＝!ｃ ２ ２ !２!ｘ２＋!２ !ｚ ２% $
打靶法进行射线追踪实验
《地震走时成像方法及技术》实验一打靶法进行射线追踪实验 专业名称：地球物理学 学生姓名： 学生学号： 指导老师：白超英 提交日期：2016-12-6
一．试验目的； 一维线性增加模型中会用打靶法进行射线追踪； 二．试验设备； 计算机，射线路径转换软件，正演模型软件 三．试验内容及要求； 1、构建包含速度异常体的井间模型，计算并绘制： 1）单炮走时场及射线路径； 2）模型和射线分布图。 注意：由于是井间模型，震源和检波器应位于井内(检波器可位于地表，即VSP)； 模型速度范围2-4km/s 为宜； 模型尺度：横向100-800m左右，纵向400m-2000m左右； 2、构建一速度随深度线性增加的模型，同样绘制： 1）单炮走时场及射线路径； 2）模型和射线分布图。 注意：震源和检波器均位于地表； 模型速度范围4-8km/s 为宜； 模型尺度：横向100-200km左右，纵向40-90km左右； 3、构建均匀模型，分析对比解析解和数值解的绝对误差和相对误差（图4）。 解析解直接由距离除以速度得到； 绝对误差= 解析解- 数值解 相对误差= 绝对误差/ 解析解* 100% 四．试验结果； 1.所选模型：水平层状，线性速度增加深地震测深模型 2.模型参数：以地面作为X轴，井间为Z轴，震源为（0，0），检波器在x=40的井间；第一层速度为v1= 3.8km/h,增加速度为0.2km/h.
图1 均匀线性增加多层介质单炮走时场及射线路径（真实情况z 轴应倒转） 计算走时 理论走时 绝对误差 相对误差 图2 走时误差分析（真实情况z 轴应倒转） 五．结果分析； 经测试，对于水平层状，线性速度增加深地震测深模型，打靶法能较好地追踪射线得到一个较准确的走时场模型。 X/distance （km ） Z /d e p t h (k m )
声波方程数值模拟实验报告
声波方程数值模拟实验报告 一．基础理论知识 需要的已知条件包括： 1．1）震源函数 2）地层速度（波速） 3）边界条件 2．弹性波方程：?????????+??=??+??+??=??) ()()(222222 22222 222z w x w v t w t S z u x u v t u s p 声波方程的有限差分法数值模拟 对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述： )()(2222 222t S z u x u v t u +??+??=?? （4-1） (,)v x z 是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，)(t s 为震 源函数。 为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。为此，先把空间模型网格化（如图4-1所示）。 设x 、z 方向的网格间隔长度为h ?，t ?为时间采样步长，则有： h i x ?= （i 为正整数） h j z ?= (j 为正整数)t n t =? (n 为正整数) k j i u , 表示在(i,j)点，k 时刻的波场值。 将1 ,+k j i u 在(i,j)点k 时刻用Taylor 展式展开： z ?,i j 1,i j +2,i j +1,i j -2,i j -,2 i j -,2 i j +,1i j +,1 i j -1,1i j -+1,2 i j -+2,1i j -+2,2 i j -+1,2 i j ++2,2 i j ++1,1 i j +-2,1i j +-1,1i j ++2,1i j ++1,1i j --1,2i j +-2,2i j +-2, 2i j --2,1 i j --1,2i j --x ?
2007射线追踪与波动方程正演模拟方法对比研究
４７ 科技资讯 科技资讯 ＳＣＩＥＮＣＥ ＆ ＴＥＣＨＮＯＬＯＧＹ ＩＮＦＯＲＭＡＴＩＯＮ２００７ ＮＯ．１２ ＳＣＩＥＮＣＥ ＆ ＴＥＣＨＮＯＬＯＧＹ ＩＮＦＯＲＭＡＴＩＯＮ 工 业 技 术 地震正演模拟作为反演解释的反过程，是验证解释成果的有效手段，进行必要可靠的正演模拟可以有效的监控反演解释。地震学一般可以分为几何地震学和物理地震学，在几何地震学中进行的正演模拟方法就是我们通常所说的射线追踪法，射线追踪法是在合成记录时用地震子波和界面或地质体的反射系数进行反褶积运算，即。运算的最大特点是说明了地震波传播的运动学特征。而在物理地震学中应用波动方程法合成的地震记录是通过求解波动方程的数值解来模拟地震波场的。在波动方程合成的地震记录中不单保持了地震波传播运动学特征，还说明了地震波传播的动力学特征。本文将分别用射线追踪和波动方程的方法合成地震记录。 １ 基于射线追踪的合成地震响应 射线追踪法的主要理论基础是，在高频近 射线追踪与波动方程正演模拟方法对比研究 王志美 畅永刚 （长江大学油气资源与勘探技术教育部重点实验室 湖北荆州 ４３４０２３） 摘 要：地震学一般可以分为几何地震学和物理地震学，几何地震学中进行正演模拟方法就是射线追踪法，射追踪法是在合成记录时用地震子波和界面或地质体的反射系数进行反褶积运算，即。运算的最大特点就是说明了地震波传播的运动学特征。而在物理地震学中的波动方程法合成的地震记录是通过求解波动方程的数值解来模拟地震波场。在波动方程合成的地震记录中不单保持了地震波传播 运动学特征外，还说明了地震波传播的动力学特征。本文将分别用射线追踪和波动方程的方法合成地震记录。关键词：射线追踪 波动方程 正演模拟 中图分类号：Ｐ３１５文献标识码： Ａ 文章编号：１６７２－３７９１（２００７）０４（ｃ）－００４７－００ 图１ 射线追踪正演模拟（１） 图 ２ 逐段迭长示意图 图 ３ 射线追踪正演模拟（２） 图４ 波动方程正演模拟结果 似条件下，地震波的主能量沿射线轨迹传播。基于这种认识，运用惠更斯原理和费马原理来重建射线路径，并利用程函方程来计算射线的旅行时。在旅行时计算中应用有限差分等方法，以获得快速的解。射线法的主要优点是概念明确，显示直观，运算方便，适应性强；其缺陷是应用有一定限制条件，计算结果在一定程度上是近似的，对于复杂构造进行两点三维射线追踪往往比较麻烦。为了计算波沿射线的旅行时和波的传播路径，叙述如下。 如图１所示，首先给出连接Ｓ（激发点）和Ｒ（接收点）之间的初始射线路径射线的振幅变化，首先必须知道地震波在实际地层中传播的射线路径。 由于地震波在整条路径上满足同一个射线参数，因此射线路径上任意连续三点也将满足同一个参数，而三点间的射线表现形式为Ｓｎｅｌｌ定律。按照Ｓｎｅｌｌ 定律，可导出一个求 取中间点的一阶近似公式。当前后两点位于界面两边时，中间点为透射点，所求路径为透射路径；当前后两点位于界面的同一边时，中间点为反射点，所求路径为反射路径。为此，可以从任一端点出发，连续地选取三点，通过一阶近似公式进行逐段迭代取中间点，利用新求出的点代替原来的点，然后以一点的跨跃作为步长，顺序地逐段迭代下去，直到另一端点。这样，新计算出的中间点和两个端点就构成了一次迭代射线路径，如图２中所示。如果整条射线路径上校正量的范数之和满足一定的精度要求，则认为射线追踪过程结束，否则从追踪出的射线路径开始，继续重复上述过程，直到满足精度要求为止。最后一次追踪到的中间点和两个端点，构成整条射线路径。图３基于多层倾斜界面模型通过射线追踪正演模拟地震响应。从模拟结果可以直观的看出基于几何地震学的原理正演模拟结果只能反映地震波的几何传播路径。在实际的工程设计中通过正演模拟可以在地表确定地下观测范围，节约设备提高工程效率，但不能反映 物理地震学中的地震属性，例如振幅，频率和相位等。更不能反映地震波的动力学特征。 ２ 波动方程的合成地震响应 ２．１ 波动方程的建立 非均匀介质的声波方程： （１） （２） 可由对连续介质方程（１）式的两端对时间求导，并利用欧拉方程推得： （３） 其中：Ｐ是波数，Ｖ是质点振动的速度向量，ρ是密度，ｃ是波速，ρ和ｃ是随着空间参数χ和ｚ变化的，这里ρ给定为常数，只有ｃ 是地质模型的控制参数。χ和Ｚ分别是在地面水平距离和深度。这样（３）式就可以变为： （４） 其中：ｃ＝ν （χ，ｚ）；（４）式即是所求的弹性波动方程。 ２．２ 数值计算及稳定性 求解弹性波动方程的方法有多种，付立叶变换法是对弹性波动方程的波场进行付立叶变换，优点是运算速度快。克希霍夫积分法是基于均匀模型，利用格林函数公式计算曲面积分，求出空间波场值，但这种方法不能适应
地震波数值模拟方法研究综述.
地震波数值模拟方法研究综述 在地学领域，对于许多地球物理问题，人们已经得到了它应遵循的基本方程(常微分方程或偏微分方程)和相应的定解条件，但能用解析方法求得精确解的只是少数方程性质比较简单，且几何形状相当规则的问题。对于大多数问题，由于方程的非线性性质，或由于求解区域的几何形状比较复杂，则不能得到解析解。这类问题的解决通常有两种途径。一是引入简化假设，将方程和几何边界简化为能够处理的情况，从而得到问题在简化状态下的解答。但这种方法只是在有限的情况下是可行的，过多的简化可能导致很大的误差甚至错误的解答。因此人们多年来寻找和发展了另一种求解方法——数值模拟方法。 地震数值模拟(SeismicNumericalModeling)是地震勘探和地震学的基础，同时也是地震反演的基础。所谓地震数值模拟，就是在假定地下介质结构模型和相应的物理参数已知的情况下，模拟研究地震波在地下各种介质中的传播规律，并计算在地面或地下各观测点所观测到的数值地震记录的一种地震模拟方法。地震波场数值模拟是研究复杂地区地震资料采集、处理和解释的有效辅助手段，这种地震数值模拟方法已经在地震勘探和天然地震领域中得到广泛应用。 地震数值模拟的发展非常迅速，现在已经有各种各样的地震数值模拟方法在地震勘探和地震学中得到广泛而有效
的应用。这些地震波场数值模拟方法可以归纳为三大类，即几何射线法、积分方程法和波动方程法。波动方程数值模拟方法实质上是求解地震波动方程，因此模拟的地震波场包含了地震波传播的所有信息，但其计算速度相对于几何射线法要慢。几何射线法也就是射线追踪法，属于几何地震学方法，由于它将地震波波动理论简化为射线理论，主要考虑的是地震波传播的运动学特征，缺少地震波的动力学信息，因此该方法计算速度快。因为波动方程模拟包含了丰富的波动信息，为研究地震波的传播机理和复杂地层的解释提供了更多的佐证，所以波动方程数值模拟方法一直在地震模拟中占有重要地位。 1地震波数值模拟的理论基础 地震波数值模拟是在已知地下介质结构的情况下，研究地震波在地下各种介质中传播规律的一种地震模拟方法，其理论基础就是表征地震波在地下各种介质中传播的地震波传播理论。上述三类地震波数值模拟方法相应的地震波传播理论的数学物理表达方式不尽相同。射线追踪法是建立在以射线理论为基础的波动方程高频近似理论基础上的，其数学表形式为程函方程和传输方程。积分方程法是建立在以惠更斯原理为基础的波叠加原理基础上的，其数学表达形式为波动方程的格林函数域积分方程表达式和边界积分方程表达式。波
基于GPU的波动方程正演模拟的实现
基于GPU的波动方程正演模拟的实现 袁崇鑫;邓飞 【期刊名称】《电脑知识与技术》 【年(卷),期】2014(000)018 【摘要】随着计算机技术的发展，使得波动方程正演由理论研究应用到实际地震勘探中成为了可能。而有限差分技术作为地震波场模拟的一种有效数值方法，它具有实现简单，速度快，从而被广泛应用正演计算密集的波形正反演中。地震波正演的计算量大，通过CPU来计算地震波正演模拟严重影响整体运算效率，GPU通用计算技术的产生及其在内的数据并行性有望改变这一状况。该文主要研究波动方程正演在GPU上的模拟实现。%With the development of computer technology, the wave equation forward by the application of theory to real seismic exploration as possible. The finite-difference seismic wave field simulation technology as an effective numerical methods, it has a simple, fast, and thus is widely used computationally intensive forward modeling and inversion of the waveform. Computationally intensive seismic forward modeling of seismic waves through the CPU to calculate the forward modeling seriously affect the over-all operational efficiency, GPU general computing technologies, including the generation and data parallelism is expected to change this situation. This paper studies the wave equation forward simulation on the GPU. 【总页数】6页(4333-4337,4340)
声波波动方程正演模拟程序总结
声波波动方程正演模拟程序 程序介绍： 第一部分：加载震源，此处选用雷克子波当作震源。 编写震源程序后，我将输出的数据复制，然后我用excel做成了图片，以检验程序编写是否正确。以下为雷克子波公式部分的程序： for(it=0;it 模型构建与试算： 1、我首先建立了一个均匀介质模型，首先利用不同时间，进行了数值模拟，得到波场快照如图所示： 100ms 200ms 300ms 此处，纵波速度为v=3000m/s。模型大小为200×200，空间采样间隔为dx=dz=10m。采用30Hz的雷克子波作为震源子波，时间采样间隔为1ms，图中可以看出，波场快照中的同相轴是圆形的，说明在均匀各向同性介质中，点源激发的波前面是一个圆，这与理论也是吻合的。并且随着时间的增大，波前面的面积逐渐增大，说明地震波从震源中心向外传播。 2、我在建立的均匀模型的基础上，改变差分算子的精度，分别采用2阶、6阶、12阶精度进行试算。时间统一采用300ms的时候。得到的波长快照如下： 2阶精度6阶精度12阶精度 16 工程地质计算机应用 2004年 第 4 期 总 36期 最短路径法射线追踪的MATLAB 实现 李志辉 刘争平（西南交通大学土木工程学院 成都 610031） 【摘要】MATLAB 是一种广泛应用的优秀的科学计算语言和工具。本文简要阐述了最短路径法的基 本原理，探讨了在MATLAB 环境中实现最短路径射线追踪的方法和步骤，并通过数值模拟演示了所 编程序在射线追踪正演计算中的应用。 【关键词】最短路径法 射线追踪 MATLAB 数值模拟 利用地震初至波确定近地表介质结构，在矿产资源的勘探开发及工程建设中有重要作 用。地震射线追踪方法是研究地震波传播的有效工具，目前常用的方法主要有有限差分解方 程法和最小路径法。最短路径方法起源于网络理论，首次由Nakanishi 和Yamaguchi 应用于 地震射线追踪中。Moser 以及Klimes 和Kvasnicha 对最短路径方法进行了详细研究。通过科 技人员的不断研究，最短路径方法目前已发展较为成熟，其基本算法的计算程序也较为固定。 被称作是第四代计算机语言的MATLAB 语言，利用其丰富的函数资源把编程人员从繁 琐的程序代码中解放出来。MATLAB 语言用更直观的、符合人们思维习惯的代码，为用户提 供了直观、简洁的程序开发环境。本文介绍运用Matlab 实现最短路径法的方法和步骤，便于 科研院校教学中讲授、演示和理解最短路径方法及其应用。 1 最短路径法射线追踪方法原理 最短路径法的基础是Fermat 原理及图论中的最短 路径理论。其基本思路是，对实际介质进行离散化，将 这个介质剖分成一系列小单元，在单元边界上设置若干 节点，并将彼此向量的节点相连构成一个网络。网络中， 速度场分布在离散的节点上。相邻节点之间的旅行时为 他们之间欧氏距离与其平均慢度之积。将波阵面看成是 由有限个离散点次级源组成，对于某个次级源（即某个 网格节点），选取与其所有相邻的点（邻域点）组成计 算网格点；由一个源点出发，计算出从源点到计算网格 点的透射走时、射线路径和射线长度；然后把除震源之 外的所有网格点相继当作次级源，选取该节点相应的计算网格点，计算出从次级源点到计算网格点的透射走时、射线路径和射线长度；将每次计算 出来走时加上从震源到次级源走时，作为震源点到该网格节点的走时，记录下相应射线路径 位置及射线长度（见图1）。 根据Fermat 原理逐步计算最小走时及射线方向。设Ω为已知走时点q 的集合，p 为与其 相邻的未知走时点，tq 分别和p 点的最小走时，tqp 为q 至p 最小走时。r 为p 的次级源位置， 图1 离散化模型（星点表示震源或次级震源，空心点为对应计算网格点） 三维波动方程正演及模型应用研究 1熊晓军，贺振华，黄德济 成都理工大学油气藏地质与开发工程国家重点实验室（610059） E-mail:xiongxiaojun@https://www.wendangku.net/doc/6812252449.html, 摘 要：为了真实准确地反映三维地质体的波场特征，在频率-波数域将二维波场延拓算子推广到三维空间，采用三维波动方程延拓方法实现了三维地质模型的快速叠后正演。该方法采用傅立叶变换进行计算不仅计算迅速，算法稳定，而且可以采用相位移加插值方法处理一定的横向变速的情况，为更加灵活方便地模拟地下复杂的三维地质体提供了一种有效的工具。作为实例，首先进行了三维French 模型的数值模拟，得到了和实际物理模型实验结果相一致的正演记录，并对比分析了三维偏移剖面和二维偏移剖面的偏移效果，验证了该方法的正确性和适用性。然后进行了三维缝洞地质模型的正演计算，得到了高信噪比的正演记录，对认识和解释三维缝洞地质体的地震波场特征很有帮助。 关键词：三维地震正演，波动方程，波场延拓，French 模型，缝洞模拟 1. 引 言 在地震资料采集、处理和解释中通常要进行地震波场的数值模拟，即假设已知地下的地质情况，应用地震波的运动学和动力学的基本原理，计算给定地质模型的地震响应，其对正确认识地震波的运动学和动力学特征，以及准确地分析油气藏的反射波场特征有着重要的指导意义。常规的正演方法主要有波动方程法和几何射线法两大类[1]。几何射线法[2]将地震波波动理论简化为射线理论，主要考虑地震波传播的运动学特征，所得的地震波的传播时间比较准确，但是计算结果很难保持地震波的动力学特征，而且对复杂的地质构造会出现盲区。波动方程法通过求解地震波波动方程，建立地下地震波运动的波场，它包含了地震波传播的所有信息，在地震模拟中占有重要地位。目前，常规的波动方程正演方法，如有限差分法[3]、积分法[4]以及F-K 域[5]等方法，主要进行二维地质体的数值模拟，而实际的地下构造往往是三维的，因此，十分有必要研究基于三维地质体的数值模拟方法。本文在二维正演方法的基础上，将二维波场延拓算子推广到三维空间，在频率—波数域进行三维叠后正演模拟。该方法不仅算法稳定，计算速度快，具有处理一定横向变速的能力，而且可以得到高信噪比的零偏移距正演记录。 2. 方法原理 利用波动方程进行地震波场数值模拟的核心是波场延拓，对于垂向变速介质，利用二维 标量波动方程，在频率—波数域可以得到各个深度间隔内的相位移延拓的正演和偏移公式[6]， i zi Z ik i x i x e z k P z k P ?+=),,(),,(1ωω (1) i zi Z ik i x i x e z k P z k P ??+=),,(),,(1ωω (2) 1 本课题得到高等学校博士学科点专项科研基金资助课题（SRFDP ）资助，编号：20040616001。 - 1 - 地震波波动方程数值模拟方法 地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶变换法、有限元法和有限差分法等。 克希霍夫积分法引入射线追踪过程，本质上是波动方程积分解的一个数值计算，在某种程度上相当于绕射叠加。该方法计算速度较快，但由于射线追踪中存在着诸如焦散、多重路径等问题，故其一般只能适合于较简单的模型，难以模拟复杂地层的波场信息。 傅里叶变换法是利用空间的全部信息对波场函数进行三角函数插值，能更加精确地模拟地震波的传播规律，同时，利用快速傅里叶变换(FFT)进行计算，还可以提高运算效率，其主要优点是精度高，占用内存小，但缺点是计算速度较慢，对模型的适用性差，尤其是不适应于速度横向变化剧烈的模型. 波动方程有限元法的做法是：将变分法用于单元分析，得到单元矩阵，然后将单元矩阵总体求和得到总体矩阵，最后求解总体矩阵得到波动方程的数值解；其主要优点是理论上可适宜于任意地质体形态的模型，保证复杂地层形态模拟的逼真性，达到很高的计算精度，但有限元法的主要问题是占用内存和运算量均较大，不适用于大规模模拟，因此该方法在地震波勘探中尚未得到广泛地应用。。 相对于上述几种方法，有限差分法是一种更为快速有效的方法。虽然其精度比不上有限元法，但因其具有计算速度快，占用内存较小的优点，在地震学界受到广泛的重视与应用。 声波方程的有限差分法数值模拟 对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述： )()(2222222t S z u 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 j 1,i j +2,i j +1,i j - 王红落,常 旭,陈传仁.基于波动方程有限差分算法的接收函数正演与偏移.地球物理学报,2005,48(2):415～422 Wang H L ,Chang X ,Chen C R.Receiver function forward m odeling and migration based on wave 2equation finite difference method.Chinese J .G eophys .(in Chinese ),2005,48(2):415～422 基于波动方程有限差分算法的接收函数正演与偏移 王红落1 ,常 旭1 ,陈传仁 2 1中国科学院地质与地球物理研究所,北京 1000292长江大学地球物理与石油资源学院,荆州 434023 摘 要 针对接收函数正演与偏移,本文采用波动方程有限差分算法.借鉴成熟的勘探地震学方法,引入等效速度概念,建立接收函数转换波与地震勘探反射波的等效走时方程,实现了基于波动方程有限差分算法的接收函数正演与偏移.数值计算表明,波动方程有限差分叠后偏移方法可以对点绕射和穹隆构造模型实现高精度成像.本文利用数值计算讨论了波动方程有限差分叠后偏移与K irchhoff 叠后偏移对于接收函数偏移的适用性,还对偏移过程中速度模型的误差进行了分析. 关键词 有限差分正演,接收函数阵列,共转换点(CCP )叠加,有限差分偏移 文章编号 0001-5733(2005)02-0415-08 中图分类号 P631 收稿日期 2004-04-13,2004-12-08收修定稿 基金项目 国家自然科学基金项目(40174017,40235055)资助. 作者简介 王红落,男,1975年生,1998年毕业于大庆石油学院,现为中国科学院地质与地球物理研究所博士研究生,主要从事勘探地震数 字处理方面的工作.E 2mail :whl @https://www.wendangku.net/doc/6812252449.html, R eceiver function for w ard modeling and migration based on w ave 2equation finite difference method W ANG H ong 2Luo 1,CHANG Xu 1,CHE N Chuan 2Ren 2 1Institute o f G eology &G eophysics ,Chinese Academy o f Sciences ,Beijing 100029,China 2The College o f G eophysics and P etroleum Resources ,Yangtze Univer sity ,Jingzhou 434023,China Abstract The wave 2equation finite difference alg orithm is applied here as a s olution for forward m odeling and migration of receiver function.With the help of exploration seism ology concepts ,the equivalent velocity and traveltime equation are established ,therefore the forward calculation and migration of receiver function are im plemented by wave 2equation finite difference method.The numerical calculation indicates that imaging of point diffractor and dome can be accom plished with high accuracy em ploying poststack wave 2equation finite difference migration.Then the applicabilities for receiver function imaging between poststack finite difference migration and K irchhoff migration are discussed according to numerical results.Finally ,the effects on velocity error for imaging are analyzed. K eyw ords Finite difference m odeling ,Receiver function array ,C omm on converted point stack ,Finite difference migration 1 引 言 地震台站记录是震源时间函数、地下介质结构 以及仪器特性耦合作用的结果.接收函数是消除震 源效应以及仪器特性后,台站下方介质结构的径向响应.早期的接收函数方法是用单个固定台站远震体波资料来反演台站下方地壳的一维速度结 第48卷第2期2005年3月 地 球 物 理 学 报 CHI NESE JOURNA L OF GE OPHY SICS V ol.48,N o.2 Mar.,2005最短路径法射线追踪的MATLAB实现
三维波动方程正演及模型应用研究
地震波波动方程数值模拟方法
基于波动方程有限差分算法的接收函数正演与偏移