射线追踪数值模拟方法研究及效果分析
第１５卷第４期
２００８年７月
断块油气田
收稿日期：２００７－１１－２０；改回日期：２００８－０４－２９。
作者简介：王向前，男，１９７８年生，助理工程师，２００３年毕业于
江汉石油学院，现从事地震采集应用软件研发及应用工作。
射线追踪数值模拟方法研究及效果分析
王向前
王
昀
（胜利石油管理局物探公司，山东东营２５７１００）
摘要地震波数值模拟方法是和地震波传播理论紧密结合在一起的。由于可以用射线理论、积分方程、微分方程来
描叙地震波的传播，模拟方法也相应地有射线追踪方法、
积分方程数值求解方法及微分方程数值求解方法。射线追踪方法通过求解程函方程计算地震波旅行时，通过求解传输方程计算地震波振幅。是以高频近似为前提，适合于物性缓变模型中地震波传播模拟，但不适合物性变化较大的模型中。其突出优点是模型简单，计算速度快。
关键词
地震模型；射线追踪；高斯束；波前扩展；反射波旅行时
文章编号：１００５－８９０７（２００８）０４－０４４－０４
中图分类号：ＴＥ１３２．１文献标识码：Ａ
Ｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎｏｆｒａｙ－ｔｒａｃｉｎｇａｎｄｔｈｅａｎａｌｙｓｉｓｏｆｉｔｓｅｆｆｅｃｔ
ＷａｎｇＸｉａｎｇｑｉａｎ（ＧｅｏｐｈｙｓｉｃａｌＥｘｐｌｏｒａｔｉｏｎＣｏｍｐａｎｙ，ＳｈｅｎｇｌｉＰｅｔｒｏｌｅｕｍＡｄｍｉｎｉｓｔｒａｔｉｏｎＢｕｒｅａｕ，ＳＩＮＯＰＥＣ，Ｄｏｎｇｙｉｎｇ２５７１００，Ｃｈｉｎａ），
ＷａｎｇＹｕｎ．
Ｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎｏｆｓｅｉｓｍｉｃｗａｖｅａｎｄｓｅｉｓｍｉｃｗａｖｅｐｒｏｐａｇａｔｉｏｎｉｓａｃｌｏｓｅｌｙｃｏｍｂｉｎｅｄｔｈｅｏｒｙ．Ａｓｓｅｉｓｍｉｃｗａｖｅｐｒｏｐａｇａｔｉｏｎｃａｎｂｅｄｅｓｃｒｉｂｅｄｂｙｒａｙｔｈｅｏｒｙ，ｉｎｔｅｇｒａｌｅｑｕａｔｉｏｎｓａｎｄｄｉｆｆｅｒｅｎｔｉａｌｅｑｕａｔｉｏｎｓ，ｓｏｔｈｅｓｉｍｕｌａｔｉｏｎｍｅｔｈｏｄｓａｒｅｃｏｒｒｅｓｐｏｎｄｉｎｇｌｙｏｆｒａｙ－ｔｒａｃｉｎｇｍｅｔｈｏｄ，ｎｕｍｅｒｉｃａｌｓｏｌｕｔｉｏｎｏｆｉｎｔｅｇｒａｌｅｑｕａｔｉｏｎｍｅｔｈｏｄａｎｄｎｕｍｅｒｉｃａｌｓｏｌｕｔｉｏｎｏｆｄｉｆｆｅｒｅｎｔｉａｌｅｑｕａｔｉｏｎｓｍｅｔｈｏｄ．Ｒａｙ－ｔｒａｃｉｎｇｍｅｔｈｏｄｃａｌｃｕｌａｔｅｓｔｈｅｔｉｍｅｏｆｓｅｉｓｍｉｃｗａｖｅｓｔｒａｖｅｌｂｙｓｏｌｖｉｎｇｔｈｅｅｉｋｏｎａｌｅｑｕａｔｉｏｎａｎｄｃａｌｃｕｌａｔｅｓｓｅｉｓｍｉｃｗａｖｅａｍｐｌｉｔｕｄｅｂｙｓｏｌｖｉｎｇｅｑｕａｔｉｏｎｏｆｔｒａｎｓｆｅｒｅｑｕａｔｉｏｎ．Ｉｔｔａｋｅｓｓｉｍｉｌａｒｉｔｙｏｆｈｉｇｈｆｒｅｑｕｅｎｃｙａｓｔｈｅｐｒｅｍｉｓｅａｎｄｆｉｔｓｆｏｒｅａｓｅｏｆｃｈａｎｇｅｏｆｔｈｅｍｏｄｅｌｓｅｉｓｍｉｃｗａｖｅｐｒｏｐａｇａｔｉｏｎｓｉｍｕｌａｔｉｏｎ，ｂｕｔｎｏｔｆｏｒｔｈｅｍｏｄｅｌｗｉｔｈｌａｒｇｅｒｃｈａｎｇｅｓｉｎｐｈｙｓｉｃａｌｃｈａｒａｃｔｅｒｉｓｔｉｃｓ．Ｔｈｅｃｏｎｓｐｉｃｕｏｕｓａｄｖａｎｔａｇｅｉｓｔｈａｔｃｏｍｐｕｔａｔｉｏｎｉｓｆａｓｔｗｈｅｎｔｈｅｍｏｄｅｌｉｓｓｉｍｐｌｅ．
Ｋｅｙｗｏｒｄｓ：ｓｅｉｓｍｉｃｍｏｄｅｌ，ｒａｙ－ｔｒａｃｉｎｇ，Ｇａｕｓｓｉａｎｂｅａｍ，ｗａｖｅｆｒｏｎｔｓｐｒｅａｄ，ｒｅｆｌｅｃｔｅｄｗａｖｅｔｒａｖｅｌｔｉｍｅ．
地震模型与正演是利用已知资料建立地震模型，或建立复杂二维和三维地质模型，进行地震模拟和面元分析，得到该地质模型的人工合成地震记录和面元分析结果，为观测系统的设计提供参考依据，如有效覆盖次数分析、ＣＲＰ面元照明度分析，并对设计方案进行评价，从而实现基于地质模型的地震采集技术，同时为地震资料解释提供工具，为反演方法提供基础。
为了达到对设计方案指导的目的，需要建立符合实际地质属性的地质模型和高精度的模拟数值计算方法。该文主要针对射线追踪数值模拟记算方法进行研究，并通过《ＳｅｉｓＷａｙ!１．０地震模型与正演》模块验证模拟计算的实际应用效果。
１
模拟方法
１．１
原理
射线追踪根据费马原理，计算出其稳定值，求解出
特征方程，再根据射线追踪方法计算出模拟数据。费马原理：
ｔ＝
Ａ′
Ａ
"ｄｌｖ（ｘ１
，ｘ２
，ｘ３
）
式中：ｌ为地震波在介质中传播的几何路程；ｘ１，ｘ２，ｘ３为反射界面点水平坐标；ｖ（ｘ１，ｘ２，ｘ３）为此点地震波速度；ｔ为反射界面点总旅行时。
１．２射线追踪方法
射线追踪所用的方法有打靶法、迭代追踪法、程函
方程射线追踪法及高斯束方法等。
１．２．１程函方程射线追踪法
１．２．１．１波前扩展法计算初至波旅行时
１）计算震源周围邻近点初至波旅行时
在三维介质中，震源点Ａ周围共有２６个邻点（ｉ＝１，２，３，…，２６），点Ｂ处的初至波旅行时间可用距离乘以平均慢度算出。计算２６个邻点初至波走时的公式为
断块油气田ＦＡＵＬＴ－ＢＬＯＣＫＯＩＬ＆ＧＡＳＦＩＥＬＤ
第１５卷第４期２００８年７月
Ｔｂｉ＝Ｓ!
ｉ
ＡＢ
ｉ
（ｉ＝１，２，…，２６）
Ｓ!＝１
２（Ｓ
Ａ
＋Ｓ
Ｂ
）
式中：Ｓ!为平均慢度；ＳＡ，ＳＢ分别为点Ａ，Ｂ的速度；ＡＢ
ｉ为震源点到不同邻点距离；Ｔｂ为邻点初至波走时。
作出反射界面点上的总旅行时与反射界面点的水平坐标关系曲线，按费马原理，震源到接收点的旅行时是稳定的。
２）波前扩展
在包含各旅行时计算点的点周边（存储于周边数列）找到最小旅行时点。该点是波前最先到达此周边面的点，解域最先从该点开始扩展。用有限差分法计算该点邻点的旅行时，将新算的点及其旅行时存入周边数列中，若该点成为内部点，则将其从周边数列中删去，即对周边进行更新。这样，周边就从这个点向它的邻区扩展，形成一个新的解域。重复该过程直到计算了模型内全部点上的时间为止。
１．２．１．２反射波旅行时计算
基于地震波传播互易性原理利用程函方程计算反射波旅行时的的方法、步骤为：
１）用有限差分法解程函方程，求出震源点到反射界面上任意点处的初至时间。
２）将接受点作为震源，通过解程函方程，求出接受点到反射界面上任意点处的初至时间。
３）将上述两步计算出的反射界面点的初至时间相加，得到界面点处的总旅行时。
４）作出反射界面点上的总旅行时与反射界面点的水平坐标关系曲线。
曲线上的极值点对应的旅行时即为反射波的旅行时，相对应的界面水平坐标值即为反射点的水平坐标。当总旅行时与界面水平坐标关系曲线存在多个极小值，说明接收点能够收到多个界面反射波。因此，只要在计算时，将界面离散点取得足够密，就可以把所有的反射波旅行时求出。
１．２．１．３射线路径追踪
利用差分法解程函方程得到的只是地震波旅行时，如何利用旅行时计算射线路径是程函方程射线追踪方法需要解决的另外一个重要问题。在程函方程有限差分法中，现在有２种射线路径追踪方法：走时最大梯度回追法和走时方程路径搜索法。这２种方法均基于费马原理。走时最大梯度反向追踪射线法是基于地震波垂直于波前面传播的事实提出的。根据地震波传播理论可知，地震波射线路径垂直于波前面（等时面），即平行于走时最大梯度方向。射线路径是从检波点向震源方向追，走时梯度最大的方向即为射线方向。射线路径之所以反向追踪是因为若从震源向检波点追，有些波的射线路径追不出来。路径分支于某条射线的折射或绕射波的分叉点，所起的作用就相当于几何射线族的一个源点。
走时方程路径搜索法的基础是走时方程。根据费马原理可知，地震波从震源Ｓ到检波点Ｒ的射线路径是所有可能路径中旅行时最小的一条。设地震波沿路径传播的旅行时为ｔ（Ｓ，Ｒ），点Ｐ为曲线上的一点，则有如下等式：
ｔ（Ｓ，Ｐ）＋ｔ（Ｐ，Ｒ）＝ｔ（Ｓ，Ｒ）
式中：ｔ（Ｓ，Ｐ）为地震波从震源Ｓ到点Ｐ的传播时间；ｔ（Ｐ，Ｒ）为地震波从点Ｐ到检波点Ｒ的传播时间；ｔ（Ｓ，Ｒ）为地震波从震源Ｓ到检波点Ｒ的传播时间。
称之为走时方程。在追踪路径之前，首先根据问题的复杂程度确定路径搜索区域Ｄ（Ｐ，ｎ）的大小，Ｐ为曲线ＳＲ上任意点；ｎ为搜索区域半径；在搜索区域Ｄ（Ｒ，ｎ）中找到点Ａ，使之满足：
ｔ（Ｓ，Ａ１）＋ｔ（Ａ１，Ｒ）－ｔ（Ｓ，Ｒ）＝ｍｉｎ
Ｐ"Ｄ（ｒ，Ｎ）
ｔ（Ｓ，Ｐ）＋ｔ（Ｐ，Ｒ）－ｔ（Ｓ，Ｒ
#$
）
式中：ｔ（Ｓ，Ａ１）代表地震波沿路径ＳＡ１传播旅行时；ｔ（Ａ１，Ｒ）代表地震波沿路径ＲＡ１传播旅行时，
这样就找到了一小段路径ＬＡ
１
Ｒ
再从Ａ１点开始继续寻找下一个点，依此做法进行直到一个点为源点Ｓ为止，即求出了整条射线路径。按照这种做法可求出任意点到源点的射线路径。
走时最大梯度回追法和走时方程路径搜索追踪射线路径方法都很简单。相对而言，第１种方法因为在射线追踪过程中不再需要计算走时，仅算走时梯度，所以速度要快一些；而第２种方法由于采用了择优策略进行射线追踪，因此稳定性更强。
１．２．２高斯束方法
为了研究沿某个主导方向传播的波，采用抛物型波动方程方法是很方便的。在射线方法中，只需要动力学射线追踪系统的实数解。但在高斯束方法中，波动方程在射线附近的解需要用到动力学射线追踪系统的复数解。在射线附近波动方程最简单的解就是高斯束。当速度结构有较大的横向变化时，克希霍夫方法可以通过射线追踪方法求出格林函数，但是当速度结构很复杂时，这种方法会遇到困难。在焦散点附近，标准的射线方法所得的振幅趋于无穷大［１］。尽管有时可以通过对不规则的振幅进行空间的平滑来得到比较满意的
王向前，等：射线追踪数值模拟方法研究及效果分析
第１５卷第４期２００８年７月断块油气田
结果，但是对多值走时的处理依然很困难。
高斯束作为波动方程集中在射线附近的高频渐近
解，即使在介质中的不规则区域也能给出较准确的结
果。其公式的推导是依据Ｃｅｒｖｅｎｙ１９８２年提出的［２］。
实现步骤是用高斯束进行正演模拟。首先将震源
处的波场分解为一系列的高斯束，也就是在震源处以
初射角!发出一系列射线，用动力学射线追踪计算其
路径和高斯束参数ｖ（ｓ），ｑ（ｓ），ｋ（ｓ），Ｐ（ｓ）和Ｌ（ｓ），射线
（射线坐标为）相对应的高斯束如下所示
ｕ!（ｓ，ｎ）＝
ｖ（ｓ）
ｕ（ｓ）
!"１／２ｅｘｐｉ"#（Ｓ）＋１
２
Ｐ（ｓ）
ｑ（ｓ）
ｎ２
!"
#$
式中：ｕ!（ｓ，ｎ）为空间点波场；ｓ，ｎ为空间点在与射线相联系的中心射线坐标系中的坐标；!为ｓ处射线的初始方向；"为频率，表示点在该射线上的位置（沿射线的弧长）。
这样就可得到以这些射线为中心射线的一系列高斯束。对特定检波器位置将所有高斯束的贡献求和（通过对积分），得到该处的波场值。
２应用效果分析
２．１二维模型射线追踪模拟应用
通过建立包含有透镜体、断层等复杂地质体，验证程函数射线追踪方法模拟数值计算的精度和正确性。观测系统参数：
炮点起始坐标：（５０００，０），ｍ；炮数：１；
检波点起始坐标：（１００，０），ｍ；检波器个数：９８；
道距：１００ｍ。
２．１．１模型建立
模型建立主要通过两个过程来完成，首先通过输入和编辑组成模型的曲线，来进行地下地质构造的描述，即地层编辑（见图１）；其次通过给闭合曲线填充不同颜色来显示不同地质体，鼠标点击不同地质体，弹出介质定义对话框，定义其纵波、横波速度以及密度。
图１地层编辑２．１．２射线追踪
根据观测系统参数，进行炮检点布设以后，得到的其中一个单炮的在地质模型中追踪的射线路径（见图２），以及创建模型的地震剖面合成记录（见图３）。
图２射线路径
图３模拟记录
２．２三维模型射线追踪模拟应用
２．２．１模型建立
通过设置模型范围、测线、地层及地层线编辑，进行Ｄｅｌａｕｎａｙ剖分网格化处理方法［３］，生成三维模型，然后通过“分析－加载观测系统”菜单加载进来（见图４）。
图４加载观测系统
２．２．２射线路径显示
模型建立和加载观测系统以后，在三维模型上可显示出来单炮点在地下构造中的射线路径（见图５）。
第１５卷第４期２００８年７月
图５单炮射线路径
２．２．３模拟结果分析（ＣＲＰ分析）
通过射线追踪，可以进行ＣＲＰ（共反射点）面元属性分析。从覆盖次数分布可以看出，平层的覆盖次数分布均匀，与观测系统设计理论上达到的覆盖次数一致，由此可以看出，射线追踪数值模拟算法计算结果精确，能达到预期的效果．第２层为起伏地层，可以看出覆盖次数分布不均匀（见图６）。
图６第２层覆盖次数
用不同颜色代表地层能量分布的情况，可以将地震波照射在界面上的最强能量的位置和能量覆盖的范围反映出来，如果对一个地质目标用一个设计的观测系统进行连续放炮和向下照射，即可观察地震波能否到达地质目标的每一个面元，从而分析评价实际采集中目的层的有效性和资料的信噪比大小（见图７）。图８为单炮模拟的地层情况，可以看出模拟的精度合理。２．３结果分析
通过《ＳｅｉｓＷａｙ!１．０地震模型与正演模块》对射线追踪数值模拟方法理论进行实例应用．通过应用结果分析，该方法理论计算结果正确合理，达到预期效果；为实际地震解释提供正确的炮集数据；通过覆盖次数、能量、方位角分析，达到指导观测系统设计的目的。
图７能量分析
图８单炮模拟记录
３结束语
通过射线追踪数值模拟计算方法的研究和实例效果分析，对其方法原理计算数值的正确性和精度进行了验证；同时射线追踪数值模拟通过高斯束计算方法很好地解决了复杂不规则性区域数值计算问题，实现了较好的模拟效果。由于实例是利用模拟工区进行分析，可能存在一定的局限性，在实际应用中通过工区实际属性使用合适的稳定性条件和选择合适的参数。
地质模型的人工合成地震记录和面元分析结果，为合理的进行观测系统的设计提供参考依据，并对设计方案进行评价，从而实现基于地质模型的地震采集技术。
参考文献
［１］ＡｋｉＫ，ＲｉｃｈａｒｄｓＰＧ．ＱｕａｎｔｉｔａｔｉｖｅＳｅｉｓｍｏｌｏｇ［Ｍ］．ＳａｕｓａｌｉｔｏＣＡ：ＵｎｉｖｅｒｓｉｔｙＳｃｉｅｎｃｅＢｏｏｋｓ，２００２．
［２］ＣｅｒｖｅｎｙＶ，ＰｒｏｐｏｖＭＭ，ＰｓｅｎｃｉｋＩ．Ｃｏｍｐｕｔａｔｉｏｎｏｆｗａｖｅｆｉｅｌｄｓｉｎｉｎｈｏｍｏｇｅｎｅｏｕｓｍｅｄｉａ－Ｇａｕｓｓｉａｎｂｅａｍａｐｐｒｏａｃ［Ｊ］．ＧｅｏｐｈｙｓＪＲ，１９８２：５０－８０．
［３］孟小红．地质模型计算机辅助设计原理与应用［Ｍ］．北京：地质出版社，２００１：１５－３０．
（编辑杨会朋）
王向前，等：射线追踪数值模拟方法研究及效果分析