2016, Vol. 59
2. 国土资源部应用地球物理综合解释理论开放实验室—波动理论与成像技术实验室, 长春 130026;
3. 东方地球物理公司物探技术研究中心, 河北 涿州 072751;
4. 中石化江苏油田分公司物探技术研究院, 南京 210046
2. Laboratory for Integrated Geophysical Interpretation Theory of the Ministry for Land and Resources of China-Laboratory for Wave Theory and Imaging Technology, Changchun 130026, China;
3. Research Center of BGP, Zhuozhou Hebei 072751, China;
4. Geophysical Research Institute of Jiangsu Oilfield Company, SINOPEC, Nanjing 210046, China
在复杂山地地震数据处理中，最重要的环节就是如何消除复杂地形和山地复杂介质所带来的各种影响.目前这些影响主要采用各种校正(静校正、低速层校正等)方法进行处理(Berryhill，1984；Yang et al.，1999)，但是即使是非常仔细的校正方法也很难很好地消除复杂地形和复杂山地介质所带来的影响(Jäger et al.，2003；孙建国，2007).因此，研发直接基于复杂山地介质的地震数据处理方法非常有必要，因为这样做的主要优点在于可以避免各种校正所造成的有效接收数据的有用信息的丢失.而对直接基于复杂山地介质的地震波正演模拟技术进行研究，则是对这些地震数据处理技术进行研究的前提和基础(孙建国，2007；孙建国和蒋丽丽，2009；Sun et al.，2011；Zhang et al.，2013).通常，在复杂山地地震波数值模拟和偏移成像研究中(Jäger et al.，2003；Sun et al.，2011；Zhang et al.，2013)，均使用确定性非均匀介质模型.该模型中，介质的非均匀性有以下两种实现形式：(1) 将物理参数表示为连续变化的空间点函数；(2) 将物理参数表示为离散变化的空间点函数，即只在计算网格节点上给出的离散函数. 然而，在实际的山地地质条件下，物理参数的空间分布具有很强的随机性(Ikelle et al.，1993).因此，在复杂山地地震波数值模拟和偏移成像研究中，更接近于实际地质条件的模型应该是能刻画物理参数随机空间分布的随机介质模型.然而到目前为止，涉及复杂山地浅表介质的随机性描述的相关研究极少.换句话说，与复杂山地有关的地震波数值模拟和偏移成像研究中，还没有考虑随机介质模型.显然，这不利于进一步深入研究地震波在复杂山地的传播特征和偏移成像效果.针对上述问题，本文以研究复杂山地随机介质模型中的射线追踪问题为基本目的.
关于随机介质模型的生成及该介质中地震波运动学特征研究方面，Ikelle等(1993)提出了随机介质的生成方法和流程；Ergintav和Canitez(1997)提出了基于指数型和高斯型自相关函数生成随机介质的方法；Iooss(1998)研究了二维随机介质中反射地震波走时问题；Baig等(2003)研究了三维随机介质中地震波的走时问题.国内，主要是姚姚和奚先研究随机介质的生成问题和随机介质中的波动方程数值模拟问题(姚姚和奚先，2004;奚先和姚姚，2001，2002)；李灿苹等(2006)采用随机介质的方法对非均匀地质体中的散射波进行了初步研究；朱生旺等(2008)研究了基于随机介质的空洞型油气储层问题.近些年，Guo等(2013)和Xu等(2013)研究了三维随机介质的生成方法.在以上研究工作中，关于复杂山地的随机介质的相关研究工作极少.考虑到复杂山地中介质小尺度的随机扰动普遍存在，并对复杂山地地震勘探，尤其是近地表地震勘探工作有重要的影响.因此，本文将对复杂山地随机介质的生成问题展开相应的研究工作.
关于射线追踪的问题，由于大量地震数据处理技术(如Kirchhoff偏移、层析成像、走时反演等)的广泛需求(瞿辰等，2007)，在过去几十年里各种射线追踪方法分别被提出：两点射线追踪法(Julian and Gubbins，1977)、波前构建法(Vinje et al.，1993)、最短路径射线追踪法(Moser，1991；张美根等，2006)、走时插值法(Asakawa and Kawanaka，1993；张东等，2009)，以及各种有限差分法：Vidale 方法(Vidale，1988)、波前扩展法(Qin et al.，1992)、快速推进法(Fast Marching Method，FMM)(Sethian andPopovici，1999)、群推进法(Grouping Marching Method，GMM)(Kim and Folie，2000)、 快速扫描法(Fast Sweeping Method，FSM)(Zhao，2004)、快速迭代法(Fast Iterative Method，FIM)(Jeong and Whitaker，2008)等. 与常规的射线追踪方法不同，复杂山地随机介质中的射线追踪方法需要具备以下特性：①能够灵活地处理复杂地形引起的不规则计算边界问题；②能够稳定高效地适应随机介质的介质参数分布特征；③能够尽可能好地兼顾精度和效率.
综合考虑以上因素及现有各种射线追踪方法的优势，提出一种新的射线追踪算法：GMM-ULTI方法(Grouping Marching Method with Upwind LinearTraveltime Interpolation，GMM-ULTI).该方法采用GMM的波前扩展方式作为走时算法的整体实现步骤，采用迎风走时插值法(ULTI)作为局部走时和射线路径计算策略，同时GMM与ULTI能通过迎风思想无条件稳定地兼容为一个整体，因此新方法同时兼顾如下优势：①迎风差分思想在面对任意复杂介质时的无条件稳定性；②走时插值法近似于二阶有限差分法的计算精度；③GMM方法的O(N)计算效率，以及对随机介质的特定适应性；④能灵活处理复杂地形引起的不规则计算边界的能力.综上所述，新方法能无条件稳定、灵活、兼顾精度和效率地处理复杂山地随机介质中的射线追踪问题.接下来，首先阐述复杂山地随机介质模型的生成方法及步骤，然后再详细阐述GMM-ULTI的具体内容.
2 复杂山地随机介质模型的生成随机介质模型的核心思想是将复杂介质视为一个大尺度的背景速度和一个小尺度的随机扰动的组合.常规随机介质的大尺度的背景介质一般采用均匀介质(整个模型的平均值)，而小尺度的非均质则是背景介质的随机扰动(Ikelle et al.，1993)，其可以描述为：
|
(1) |
其中，x 代表空间坐标向量；v0代表大尺度的非均质，常规随机介质假设其为整个随机介质参数平均值的均匀介质；δv(x)代表小尺度的随机非均质扰动，其一般定义为加在背景介质上的扰动δv(x)=σ(x)v0；σ(x)为均值为零、具有一定方差和自相关函数的随机过程；v(x)为最终要构建的随机介质模型.
由式(1)可知，生成随机介质的关键是计算随机过程σ(x).根据随机过程理论，随机过程σ(x)的功率谱其实是其自相关函数φ(x)的傅氏变换Φ(kx).因此，可由已知的功率谱函数Φ(kx)来模拟产生由它描述的随机过程σ(x)，并最终获得随机介质模型.基于以上的阐述，选取自相关函数的类型，是决定生成随机介质特点的关键因素，为了灵活控制该因素在此选择混合型自相关函数(奚先和姚姚，2002)：
|
(2) |
其中，a、b分别为x、z方向上的自相关长度；r为粗糙度因子，若其分别取0和1时混合型自相关函数将变为高斯型和指数型椭圆自相关函数.参数a、b直接决定着随机介质模型的特点.关于该问题，将在4.1节中进行详细的实例讨论.
与常规随机介质模型的生成不同，为了生成更接近实际介质特点的复杂山地随机介质模型，在此创新性地引入地形修正和背景介质修正的策略(如图 1所示).同时，因为采用混合型自相关函数，从而能更加灵活地控制随机介质的特点.下面详细阐述复杂山地随机介质模型的生成过程：
|
图 1 复杂山地随机介质模型的生成 Fig. 1 The generation of complex mountainous random media model |
(1) 采用快速傅氏变换对自相关函数φ(x)进行二维傅氏变换获得Φ(k x)；
(2) 生成 0，2π 内满足均匀分布的独立的二维随机数θ(k x)；
(3) 在FFT域，计算
(4) 对Ψ(k x)进行二维傅氏反变换回空间域，进而获得空间的随机扰动ψ(x)；
(5) 对空间中的随机扰动进行地形修正，获得带地形的空间随机扰动T(x)；
(6) 计算复杂山地随机扰动T(x)的均值μ和方差d2；
(7) 规范化参数均值为0、方差为ε2，进而获得随机过程σ(x)；
(8) 综合利用(1)式和随机过程σ(x)生成最终的复杂山地随机介质模型.
上述步骤(5)中，地形修正的具体方法为：首先，读入地形高程数据，根据地形的起伏程度确定在深度方向上地形修正的范围，即地形最高点和最低点间的高程差；其次，以地形高程数据为判断条件，将地表以上的随机扰动ψ(x)的值均赋值为零(即地表以上为空气，没有随机扰动)，将地表以下的随机扰动ψ(x)保持原值不变.当然，上述地形修正会破坏随机扰动整体的规范化，所以在此特将地形修正放在模型规范化(步骤(6)、(7))之前进行.因此，最后生成的带地形的随机介质模型是满足随机介质的规范化要求的.此外，要想生成的复杂山地随机介质模型更加符合实际情况，则最好在修正过程中读入来自于实际大地测量的地形高程数据.不过，不管是采用数学模拟的高程数据，还是采用实际测量的高程数据，只要将修正步骤放在模型规范化之前进行即可保证生成随机介质的规范性，也就是说上述如此操作的地形修正是合理的.
上述步骤(8)中，关于(1)式中的背景介质的修正，本文的修正思想实际上是对常规随机介质模型生成思想的广义化.其目的是让新生成的随机介质不仅能包含介质参数小尺度上的随机扰动，还能包括大尺度上的非均匀性.其中大尺度上的非均匀性是通过非均匀的背景介质来表述.为了不失一般性，本文采用随深度线性变化的梯度介质(即v(z)=v0+az km·s-1，其中v0为地表速度，a为变化因子，z为深度)作为背景介质.当然，该梯度介质肯定不能描述所有复杂介质的分布特点.所以除了梯度介质以外，还可以根据介质的复杂程度和特点采用其他复杂的背景介质(如包含局部介质参数突变的背景介质等).不过，无论采用怎样的背景介质，均是可以实施上述随机介质模型的生成步骤的.因为背景介质的选择和加入是在步骤(8)被实施的，因此其不影响空间随机扰动的生成(步骤(1)—(7)).基于此，严格意义上讲，在此背景介质是可以选择任意复杂介质的，也即本文的算法是可以模拟任意复杂的随机介质的.
综上所述，与常规的随机介质模型生成不同，上述复杂山地随机介质模型的生成步骤主要引入两方面的核心修正内容，即地形修正和背景介质修正.这两方面的修正内容的共同作用，使得最终生成的复杂山地随机介质不仅包含了复杂地表形态，还能同时更接近实际的描述复杂介质的大尺度非均值性和小尺度的随机扰动.
3 基于GMM-ULTI的新射线追踪算法 3.1 算法的提出如图 2所示，粗实线圈定的即为复杂山地随机介质模型.在该介质中进行射线追踪时，与常规的射线追踪算法不同，算法需要同时处理因地表引起的不规则计算边界、随机介质的参数分布的复杂性问题、算法的稳定性、计算精度、计算效率等问题.为了同时很好地处理如上问题：首先，采用三角网格和正方形网格的混合网格剖分复杂地形，该网格能保证对地表形态的准确刻画，以及算法的绝大部分计算工作在规则的正方形网格中进行；其次，综合FMM中的迎风思想、LTI的局部走时计算公式的建立思想、GMM波前扩展的基本思想、以及LTI射线路径计算的基本方法，提出了一种新的射线追踪算法：GMM-ULTI算法.
|
图 2 新算法的提出 Fig. 2 The introduction of the new algorithm |
首先从整体上讲，新射线追踪算法采用的是一个两步法，即正向的走时计算和逆向的射线路径计算.在正向的走时计算环节，综合利用了GMM、迎风思想和LTI算法.其中，GMM用于整体高效的波前扩展实现，该扩展方式能很好地适应随机介质特征；LTI算法则用于构建局部走时计算公式，该算法具有近似于二阶迎风有限差分法的精度；迎风思想用于保证：算法的无条件稳定性，LTI算法与GMM的无缝对接，以及走时计算的局部实现策略等问题.而逆向的射线路径计算环节是在走时计算完成后开始的，其基本原理是基于地震波波前与射线路径相互垂直的物理事实.具体实施时实际上追踪的是走时场的负梯度方向，通过追踪一个个子震源点来获得从震源点出发的某一个运动质点的运动轨迹，即射线路径.下面将分别阐述正向走时计算和逆向射线路径计算的具体算法和实现策略.
3.2 正向走时计算正向的走时计算实际上包括两个方面的核心问题： ① 算法的整体实现步骤，即波前扩展方式； ② 算法的局部实现策略和计算公式.在此前者采用改进后的GMM算法，后者采用ULTI算法，二者综合起来被称为GMM-ULTI算法.
GMM(Group Marching method)是一种高效的波前扩展方式，其由Kim和Folie(2000)提出.它是FMM的一种更高级的形式，与FMM通过查找波前上的全局最小值作为扩展点不同，GMM则是更加贴近波前传播实际情况地一次性同时选取满足一定条件的一组局部极小值点作为扩展点，这实际上是Huygens原理的更准确体现和实现.此外，相关研究还表明GMM对随机介质有独特的适应性(Jeong and Whitaker，2008)，其能高效地处理介质速度参数的随机扰动问题.如图 3e所示，GMM在实现时将网格节点的属性分别标注为idTT=0、1、2，它们分别表示走时计算未开始、波前上、走时计算已完成的网格节点.GMM的波前推进则是通过不断更新网格节点的属性变换来实现的.波前扩展时，波前上idTT=1的满足要求的一系列网格节点将作为波前上的一系列子震源点，即扩展点.更新迭代的过程包括： ① 将未开始计算的其他网格节点的属性由0变为1； ② 将作为扩展点的子震源点的属性由1变为2.如此循环迭代下的结果是计算区域内的网格节点的属性分别由0变为1，再由1变为2，最后当所有网格节点的属性均为2时则计算完成.
|
图 3 GMM-ULTI算法的基本原理示意图 Fig. 3 The fundamental principle of GMM-ULTI method |
常规GMM一般是在规则的矩形计算区域被实现的，为了让其适应复杂山地随机介质，只需要在初始化时将地表以上的网格节点类型定义为idTT=2(由如下的初始化和迭代过程可知该属性能保证波前扩展不能穿透地表进入地表以上，这也是满足地震波的传播规律的).总体上讲，如图 3e所示，复杂山地随机介质中GMM的实现步骤为：
初始化过程：
(I1) 将震源点(si，sj)及地表以上网格节点的走时值TT(si，sj)设为零，其属性设为：
idTT(si，sj)=2；
(I2) 将与震源点直接相邻的网格节点的走时值设为解析值，属性设为：
idTT(，)=1，它们构成初始波前；
(I3) 将其他的网格节点走时值设为无穷大(e.g.，TT(，)=1.0×106)，并将它们的属性设置为：
idTT(，)=0；
(I4) 设置波前内扩展点选取的控制条件 delTAU=
波前迭代扩展过程：
(M1) 设置TM=TM+delTAU，即波前向前推进一步;
(M2) 对波前上的所有节点(i，j)进行判断，如果TT(i，j)≤TM，则重新计算节点属性idTT≤1的邻近网格节点的走时值；
(M3) 对波前上属性满足TT(i，j)≤TM的扩展点邻近的网格节点作如下处理：
(a) 如果邻近点的属性为idTT≤1，则计算其走时值；
(b) 如果邻近点的属性为idTT=0，则将其属性改为idTT=1，即将其纳入波前；
(c) 将波前上作为扩展点的网格节点移出波前，其属性改为idTT=2；
(M4) 判断波前上点数是否为零，不是则跳回M1继续计算，否则计算完成.
与常规的GMM算法相比，此处复杂地形条件下的GMM仅改进了初始化的环节：将地表以上的网格节点的属性设为idTT=2.然而正是这一小小的改进，却将波前的扩展传播限定在了地表以下，而不能穿越地空界面，进而满足了地震波在复杂地形区域的传播规律，同时也让修正后的GMM算法能够灵活地处理由复杂地形引起的不规则计算边界问题.
针对复杂山地随机介质地震波走时计算的GMM-ULTI算法，除了如图 3e所示的修正后的GMM外，还需要另外一个核心算法，即ULTI算法.该算法要解决的核心问题是：如图 3e所示在修正后的GMM循环递推实施过程(步骤M2、M3)中，均需要计算当前扩展点周围的网格节点的走时值.而解决该问题的主要难点为：如图 3e虚线矩形框圈定的网格单元所示，由于复杂地表的存在使得被算点所处网格形态不同，同时被算点周围网格节点的走时值的已知情况也是不同的(周围网格节点的属性不同，且分布情况也不同).针对被算点所处网格形态不同的问题，如图 3a、3c、3f所示将被算点的类型分为：地表上的点(如图 3a中的点C)、地表以下与地表直接相邻的规则网格节点(如图 3c中的点C)以及与地表不相邻完全规则的网格节点(如图 3f中的点C).关于被算点周围邻近点走时值已知情况不同的问题，如图 3a、3c、3f所示，一般采用在包围点C的所有线段(图 3a上的1~4条，图 3c上的1~8条、图 3f上的1~8条)上分别进行试算，然后根据Fermat原理，选取最小的走时计算值作为最终的计算结果.但是该算法需要消耗很大的计算量.为了无条件稳定地处理该问题，在此引入FMM中的迎风思想(Sethian and Popovici，1999).该思想的核心是：先基于Fermat原理确定走时值分布最小的网格线，然后再以此为已知条件进行走时计算.该原理认为：地震波总是从走时值更小的区域向走时值未知的区域传播.所以，应该选取走时值分布最小的网格线作为计算图 3a、3c、3f中被算点C的走时值的已知条件.
下面以规则网格节点为例(其他两个类型的网格与此类似，只是参与筛选的网格节点更少)，假设该网格线为如图 3f所示中的线段AminBmin，则点Amin和Bmin的走时值TAmin、TBmin应分别满足：TAmin= min(TA1，TA2，TA3，TA4)、TBmin=min(TA+，TA-)，A+和A- 则为直接与Amin相邻的网格节点.如此一来，直接以线段AminBmin为已知条件，即可一次性的计算被算点C的走时值.而该计算策略满足地震波传播的熵守恒定律和Fermat原理(Sethian and Popovici，1999)，因此是无条件稳定的.确定AminBmin为已知条件后，接下来最重要的问题是如何建立被算点C的走时TC的计算公式.如图 3b、3d、3g所示，整体上可以将公式的建立分为两种情况，即三角形网格中的和正方形网格中的.具体计算时引入LTI算法，该算法具有精度高且不受网格形态限制的优点(Asakawa and Kawanaka，1993).下面首先讨论任意三角形网格中的公式推导.
如图 3b、3d所示，为了公式推导方便，在此建立局部直角坐标系xDy；假定到达点C的射线由线段AminBmin上的点E传入，设点E在局部直角坐标系xDy上的坐标为(x，0)，则可以建立TC关于x的表达式f(x).根据Fermat原理，求取f(x)的最小值，即为TC的值；而此时对应的坐标x则指示E点的位置，将该坐标x带回表达式f(x)即可获得计算TC的公式：
|
(3) |
其中，TC为被算点C的走时值，TAmin、TBmin分别为走时值已知的点Amin、Bmin的走时值，ΔT=TAmin-TBmin；d1、d2、d3分别是线段AminBmin、AminD、DC的长度，s为局部计算网格单元的慢度(速度的倒数)的平均值.(3)式中分别对应的E点的坐标为：
|
(4) |
其中，ΔT、d1、d2、d3、s的意义和公式(3)相同.公式(3)便是局部走时计算公式.而公式(4)则表明的是C点的上一级子震源点E在各种走时值已知情况下的坐标，该公式可以用于解决在不同走时值分布情况下的上一级子震源的定位问题，其为射线路径计算的核心公式.同理可以确定在规则的正方形网格(与图 3b、3d类似，不过此时图 3g中的规则网格中，D点与Bmin点重合)中，计算TC和E点的坐标的公式分别为：
|
(5) |
|
(6) |
其中，h为网格间距.同样，公式(6)表明的是C点的上一级子震源点E在各种走时值已知情况下的坐标，该公式可以用于解决在不同走时值分布情况下的上一级子震源的定位问题，其为射线路径计算的核心公式.
综上所述，正向走时计算的核心算法为GMM-ULTI算法.其中，GMM用于整体波前扩展计算，而邻近点的走时值的更新与计算则通过LTI算法来实现.此外，迎风思想用于GMM与LTI算法的无条件稳定的衔接，一方面其能处理GMM迭代计算过程中被算点周围的走时值的任意已知情况；另一方面其能让LTI算法无条件稳定且简洁地处理用于走时计算时的线段的选取问题.最后，ULTI在局部走时计算公式推导时，实际就隐含着将射线路径计算的局部公式建立起来了，下面将结合公式(4)、(6)来阐述逆向的射线路径计算方法.
3.3 逆向射线路径计算如图 4所示，射线路径计算时采用的是从接收点到震源的逆向追踪策略，即从接收点出发一步一步追踪上一级震源点的过程.下面以一条射线路径的计算为例，阐述射线路径的计算步骤：
|
图 4 射线路径的计算策略 Fig. 4 The computational strategy of ray path |
(T1) 确定接收点G的位置，其作为射线路径计算的起点；
(T2) 对G点周围的走时值分布情况做判断，根据Fermat原理走时值最小的分布区域为G点上一级子震源点R1所在的区域，并根据公式(4)或(6)计算R1的坐标；
(T3) 将R1作为类似于G点的新的接收点，然后采用类似于(T2)的步骤计算其上一级子震源点 R2的坐标；
(T4) 重复上述计算步骤，直到追踪到震源点S为止，则运动轨迹(S，Rn，Rn-1，…，R2，R1，G)即为追踪的射线路径.
上述射线路径计算过程实际上是根据Fermat原理追踪上一级子震源点的过程.该射线追踪方法最显著的优势在于，其可以灵活地追踪从震源点出发到达指定位置的射线路径，因为该算法是一种面向目标的算法，更贴切地说其是直接始于目标点的算法.与此同时，其无需做震源点处的初始角度的迭代试射，也无需做震源点和接收点固定的最小走时射线路径的迭代优化计算，其直接追踪出来的即为最小走时路径.所以进一步说，该射线追踪算法是一种直接面向目标无需繁琐低效迭代更新与优化的一种射线追踪算法.此外，该射线追踪方法还可以通过加密空白区内的接收点的密度，来处理传统射线追踪算法在介质区域内覆盖不均匀的问题.
4 算法分析 4.1 复杂山地随机介质模型的参数分析为了分析各模型参数对随机介质生成的影响，以及研究采用本文方法生成的复杂山地随机介质的特点，在此列出了几个模型生成的实例.如图 5所示，模型的大小均为：10.0 km×7.0 km，随机介质模型在x和z方向上的自相关长度、地表形态、背景介质速度分别如图题所示.其中，背景为均匀介质的速度均值v=2.0 km·s-1；背景为梯度介质的速度表达式为：v(z)=1000.0+4.0×z km·s-1.图 5a—5f分别显示了基于不同参数生成的随机介质模型.下面将对比分析讨论这些随机介质模型的特点.
|
图 5 复杂山地随机介质的参数分析 (a)a=5.0 m，b=5.0 m，常规随机介质；(b)a=5.0 m，b=5.0 m，复杂山地随机介质；(c)a=20.0 m，b=20.0 m 常规随机介质；(d)a=20.0 m，b=20.0 m 复杂山地随机介质；(e)a=50.0 m，b=20.0 m，复杂山地随机介质；(f)a=200.0 m，b=20.0 m，复杂山地随机介质. Fig. 5 The parameter analysis of complex mountainous random media (a)a=5.0 m，b=5.0 m，the conventional random media;(b)a=5.0 m，b=5.0 m，the complex mountainous random media;(c)a=20.0 m，b=20.0 m，the conventional random media;(d)a=20.0 m，b=20.0 m，the complex mountainous random media;(e)a=50.0 m， b=20.0 m，the complex mountainous random media;(f)a=200.0 m，b=20.0 m，the complex mountainous random media. |
分析图 5取不同参数时生成的随机介质模型可以发现： ① 从整体上看，随着自相关长度a、b的逐渐增大，介质随机扰动的尺度逐渐增大(对比图 5a、5c可知)； ② 自相关长度a、b分别影响着x、z方向上的扰动尺度大小，a越大x方向上的扰动尺度越大(如图 5e所示)，当a比b大很多时整个介质趋向于薄互层状介质(如图 5f所示)； ③ 复杂山地随机介质是对常规随机介质经过地形修正后而得到的，它有着比常规随机介质更复杂的地表边界，该复杂边界会给各类地震波数值模拟带来一些难题(对比分析图 5a、5b)； ④ 采用梯度介质作为背景介质，既能保证随机介质扰动的特点，又能更加接近实际介质的分布特点，因为实际介质的地震波速度一般情况下是随着地下深度的增加而增加(对比分析图 5a、5c与图 5b、5d—5f可知).
综上所述，本文提出的随机介质的生成方法同时处理了地形存在和地下实际介质大体为梯度介质的问题.进而，使得最终生成的随机介质能够带地形，同时介质速度更贴近实际介质参数分布的特点，并且还很好地保留了随机介质随机扰动的特点.因此，笔者认为本文生成的复杂山地随机介质模型，能够作为更贴近实际的、更能准确描述介质参数分布特点的一种地下介质的模型刻画方式.
4.2 走时算法分析由以上阐述可知，新射线追踪GMM-ULTI算法能兼顾精度和效率，下面将通过具体的计算结果的精度和效率分析来验证该观点.精度和效率分析时，选取一个均匀介质模型，模型大小为2.0 km×2.0 km，介质的速度为1.0 km·s-1，震源点置于(0.0 km，0.0 km)处.在此，将引入快速推进一阶迎风差分法(FMM-1stUFD)、快速推进二阶迎风差分法FMM-2ndUFD、快速推进迎风线性插值法(FMM-ULTI)作对比.
对比分析表 1、表 2可以发现：在计算精度方面，新算法具有和FMM-ULTI算法一样的计算精度，它们均有着和二阶迎风差分法相当的计算精度，却有着比一阶迎风差分法高很多的计算精度；在计算效率方面(计算的硬件设备参数为：Intel Core i5-3230 2.60 GHz CPU，4.00 G内存)，新算法有着最高的计算效率，其CPU耗时仅为快速推进一阶迎风差分法的61.93%，快速推进二阶迎风差分法的45.74%，快速推进迎风插值法的56.92%.综合分析数据可以得出结论：新算法与快速推进一阶迎风差分法、二阶迎风差分法、迎风线性插值法相比，具有最高的计算精度和效率.因此，新算法是兼顾精度和效率的.
|
表 1 均匀介质中的计算效率对比(CPU耗时，s) Table 1 The efficiency(CPU running time，s) contrast in uniform media |
|
|
表 2 均匀介质中的计算精度对比(平均相对误差，%) Table 2 The accuracy(the average errors，%) contrast in uniform media |
得出上述结论的理由：新算法与FMM-ULTI均采用迎风线性插值法作为局部走时计算公式，所以它们有着相同的计算精度.此外，因为新算法采用了比其他3种算法的FMM更高级的波前扩展方式(GMM)，所以它有着更高的计算效率.与FMM相比，GMM在波前扩展时选取满足条件的一组扩展点，同时进行扩展计算，而不是FMM的一个扩展点.这不仅提高了效率，而且是波前传播过程更加贴切的模拟过程，因为波前是同时向外传播的，而并非单点向外扩展的.
此外，为了分析新算法对随机介质的适应性.在上述同样大小和计算参数的随机介质模型(随机介质的自相关长度为：a=b=100.0 m)中，对上述四种算法进行计算效率分析(因为不易获得随机介质中的解析走时值，所以在此不做计算精度分析).表 3给出了四种算法的CPU耗时，分析表中的数据可知在随机介质中新算法仍然是四种算法中计算效率最高的.同时对比表 3和表 1可知，当介质变复杂后(由均匀介质变为复杂随机介质)新算法的计算效率并没有明显的降低(CPU耗时平均仅增加2.23%)，但是与新算法计算精度相当的二阶迎风差分法的计算效率却因为介质的复杂而计算效率明显降低(CPU耗时平均增加38.35%)，这说明新算法在面对随机介质时有很好的适应性和稳定性.
|
|
表 3 随机介质中的计算效率对比(CPU耗时，s) Table 3 The efficiency(CPU running time，s) contrast in random media |
为了分析本文算法的射线路径计算的精度，以及对比分析均匀介质与随机介质中射线路径的差别，采用与图 5a大小一致的均匀介质模型，介质的速度参数v=2000.0 m·s-1.图 6a—6d分别给出了网格间距为100.0 m、50.0 m、25.0 m、10.0 m时，数值计算的射线路径与解析值的对比.通过分析图 6中的变化规律可以发现：在相同的计算区域内网格间距越小，数值计算出来的射线路径与理论解析的射线路径吻合得越好(尤其是当网格间距取10.0 m时，如图 6d所示数值计算的射线路径基本与实际射线路径重合，因此在相同模型大小中后续计算实例采用10.0 m作为网格间距)；同时离震源较远的区域计算精度相对较高，离震源近的区域计算精度相对低一些.总体上讲数值计算出来的射线路径和实际射线路径吻合得很好，说明算法能够保证很好的计算精度.
|
图 6 射线路径计算的精度分析 图中实线为理论解析计算的射线路径，虚线为本文算法计算得出的射线路径；(a—d)为网格间距分别取100.0 m、50.0 m、25.0 m、10.0 m时的计算结果. Fig. 6 The accuracy analysis of the ray path computation In the figure，the solid lines are the analytical solutions and the dotted lines are the computational solutions;(a—d)are the computational results when the grid spacing are 100.0 m，50.0 m，25.0 m，10.0 m，respectively. |
为了验证本文射线追踪算法在面对复杂山地随机介质时的稳定性和有效性，分别在如图 5所示的各个随机介质模型中采用本文算法进行射线追踪.分析计算结果可得出如下结论与现象：
① 对比分析图 6d和图 7a可知，与常规均匀介质相比，介质中速度参数的随机扰动，对地震波传播的波前面和射线路径均有着不同程度的随机扰动，地震波的等时面也从均匀的同心圆弧变为了局部曲率剧烈变化的曲线，射线路径则由直线变为了随机扰动的曲线；
|
图 7 随机介质中的等时线与射线路径分布 (a)a=5.0 m，b=5.0 m，常规随机介质；(b)a=5.0 m，b=5.0 m，复杂山地随机介质；(c)a=20.0 m，b=20.0 m，常规随机介质；(d)a=20.0 m，b=20.0 m 复杂山地随机介质；(e)a=50.0 m，b=20.0 m，复杂山地随机介质；(f)a=200.0 m，b=20.0 m，复杂山地随机介质. Fig. 7 The distribution of traveltime contour and ray path in random media (a)a=5.0 m，b=5.0 m，the conventional random media;(b)a=5.0 m，b=5.0 m，the complex mountainous random media;(c)a=20.0 m，b=20.0 m，the conventional random media;(d)a=20.0 m，b=20.0 m，the complex mountainous random media;(e)a=50.0 m，b=20.0 m，the complex mountainous random media;(f)a=200.0 m，b=20.0 m，the complex mountainous random media. |
② 对比分析图 7b、7d、7f可知，随着随机介质的随机扰动尺度的增加，地震波的波前面和射线路径被改造的程度越大，其中等时面的波前曲率变化程度越大，而射线路径弯曲和分布不均匀的程度越严重；
③ 对比分析图 7c、7e可知，复杂山地随机介质中地震波传播特征更为复杂，在地表附近等时面和射线路径除了受随机介质的局部扰动影响外，还受地形起伏的影响，在很多地形起伏剧烈的区域还产生了绕射波等.
6 结论与讨论为了准确刻画复杂山地介质的非均质特点和研究复杂山地介质中的射线追踪问题.首先，通过在常规随机介质生成过程中，加入地形修正和背景介质修正来建立复杂山地随机介质模型.通过数值实例验证，新提出的复杂山地随机介质模型具有如下特点： ① 新模型带有复杂的地表形态； ② 新模型采用梯度介质为背景介质，使得模型更贴近实际介质速度参数的分布特点.其次，通过融合多种射线追踪算法的优势，提出了一种新的基于GMM-ULTI法的射线追踪算法.通过数值分析和计算实例验证，新提出的射线追踪算法具有如下特点： ① 算法在走时计算时兼顾精度和效率； ② 算法能保证射线路径的计算精度； ③ 算法能无条件稳定地处理复杂地形引起的不规则计算边界条件的问题； ④ 算法能无条件稳定地处理介质的随机扰动问题.最后，通过计算实例分析可得出复杂山地随机介质中地震波的一些传播规律： ① 复杂地表附近地震波的传播特点非常复杂，地形的剧烈起伏会引起明显的绕射波，进而引起地震波等时面和射线路径的严重扰动； ② 随机介质中，地震波的等时面会因为介质不同程度的随机扰动而引起波前曲率的随机扰动，而地震波的射线路径也由直线变成了曲率不同的随机扰动曲线，同时地下射线路径的分布也因为介质的随机扰动而分布不均匀.这些模型建立方法、射线追踪算法、地震波传播规律对于提高复杂山地区域地震勘探的质量具有一定的理论意义和实际价值.
| Asakawa E, Kawanaka T. 1993. Seismic ray tracing using linear traveltime interpolation. Geophysical Prospecting , 41(1): 99–111. | |
| Baig A M, Dahlen F A, Hung S H. 2003. Traveltimes of waves in three-dimensional random media. Geophysical Journal International , 153(2): 467–482. | |
| Berryhill J R. 1984. Wave-equation datuming before stack. Geophysics , 49(11): 2064–2066. | |
| Ergintav S, Canitez N. 1997. Modeling of multi-scale media in discrete form. Journal of Seismic Exploration , 6(1): 77–96. | |
| Guo N C, Shen Z H, Peng G, et al. 2013. An improved method to establish 3D random media with flexible preferred orientation. // 83rd Annual International Meeting, SEG, Expanded Abstracts, 3559-3563. | |
| Ikelle L T, Yung S K, Daube F. 1993. 2-D random media with ellipsoidal autocorrelation functions. Geophysics , 58(9): 1359–1372. | |
| Iooss B. 1998. Seismic reflection traveltimes in two-dimensional statistically anisotropic random media. Geophysical Journal International , 135(3): 999–1010. | |
| Jäger C, Hertweck T, Spiner M. 2003. True-amplitude Kirchhoff migration from topography. //73rd Annual International Meeting, SEG, Expanded Abstracts, 909-912. | |
| Jeong W K, Whitaker R T. 2008. A fast iterative method for eikonal equations. SIAM Journal on Scientific Computing , 30(5): 2512–2534. | |
| Jing X L, Yang C C, Wang S Q. 2007. A improved seismic reflection tomographic method. Chinese J. Geophys. (in Chinese) , 50(6): 1831–1836. | |
| Julian B R, Gubbins D. 1977. Three-dimensional seismic ray tracing. J. Geophys. , 43(1-2): 95–113. | |
| Kim S, Folie D. 2000. The group marching method: An O(N) level set eikonal solver. // 70th Annual International Meeting, SEG, Expanded Abstracts, 2297-2300. | |
| Li C P, Liu X W, Li M F, et al. 2006. Preliminary study on scattering wave characteristics of heterogeneous geologic bodies. Geophysical Prospecting for Petroleum (in Chinese) (in Chinese) , 45(2): 134–140. | |
| Moser T J. 1991. Shortest path calculation of seismic rays. Geophysics , 56(1): 59–67. | |
| Qin F H, Luo Y, Olsen K B, et al. 1992. Finite-difference solution of the eikonal equation along expanding wavefronts. Geophysics , 57(3): 478–487. | |
| Qu C, Zhou H L, Zhao D P. 2007. Deep structure beneath the west margin of Philippine Sea Plate and South China Sea from P and S wave travel time tomography. Chinese J. Geophys. (in Chinese) , 50(6): 1757–1768. | |
| Sethian J A, Popovici A M. 1999. 3-D traveltime computation using the fast marching method. Geophysics , 64(2): 516–523. | |
| Sun J G. 2007. Methods for numerical modeling of geophysical fields under complex topographical conditions: a critical review. Global Geology (in Chinese) (in Chinese) , 26(3): 345–362. | |
| Sun J G, Jiang L L. 2009. Orthogonal curvilinear grid generation technique used for numeric simulation of geophysical fields in undulating surface condition. Oil Geophysical Prospecting (in Chinese) (in Chinese) , 44(4): 494–500. | |
| Sun J G, Sun Z Q, Han F X. 2011. A finite difference scheme for solving the eikonal equation including surface topography. Geophysics , 76(4): T53–T63. | |
| Vidale J E. 1988. Finite-difference calculation of traveltimes. Bull. Seismol. Soc. Am. , 78(6): 2062–2076. | |
| Vinje V, Iversen E, Gjøystdal H. 1993. Traveltime and amplitude estimation using wavefront construction. Geophysics , 58(8): 1157–1166. | |
| Xi X, Yao Y. 2001. 2-D random media and wave equation forward modeling. Oil Geophysical Prospecting (in Chinese) (in Chinese) , 36(5): 546–552. | |
| Xi X, Yao Y. 2002. Simulations of random medium model and intermixed random medium. Earth Science-Journal of China University of Geosciences (in Chinese) (in Chinese) , 27(1): 67–71. | |
| Xu Z W, An Y, Zhi L, et al. 2013. Simulation of 3D random medium in seismic prospecting using fractal method. // 83rd Annual International Meeting, SEG, Expanded Abstracts, 4672-4676. | |
| Yang K, Wang H Z, Ma Z T. 1999. Wave equation datuming from irregular surfaces using finite difference scheme. //69th Annual International Meeting, SEG, Expanded Abstracts, 1465-1468. | |
| Yao Y, Xi X. 2004. Regionalized multi-scale random medium model and its wavefield analysis. Geophysical Prospecting for Petroleum (in Chinese) (in Chinese) , 43(1): 1–7. | |
| Zhang D, Fu X R, Yang Y, et al. 2009. 3-D seismic ray tracing algorithm based on LTI and partition of grid interface. Chinese J. Geophys. (in Chinese) , 52(9): 2370–2376. doi: 10.3969/j.issn.0001-5733.2009.09.023. | |
| Zhang D L, Zhan G, Schuster G T. 2013. Multi-source least-squares reverse time migration with topography. //83rd Annual International Meeting, SEG, Expanded Abstracts, 3736-3740. | |
| Zhang M G, Jia Y G, Wang M Y, et al. 2006. A global minimum traveltime raytracing algorithm of wavefront expanding with interface points as secondary sources. Chinese J. Geophys. (in Chinese) , 49(4): 1169–1175. | |
| Zhao H K. 2004. A fast sweeping method for eikonal equations. Mathematics of Computation , 74(250): 603–628. | |
| Zhu S W, Wei X C, Qu S L, et al. 2008. Description of the carbonate karst reservoir with random media model. Acta Geologica Sinica (in Chinese) (in Chinese) , 82(3): 420–427. | |
| 井西利, 杨长春, 王世清. 2007. 一种改进的地震反射层析成像方法. 地球物理学报 , 50(6): 1831–1836. | |
| 李灿苹, 刘学伟, 李敏锋, 等. 2006. 非均匀地质体散射波特征初探. 石油物探 , 45(2): 134–140. | |
| 瞿辰, 周蕙兰, 赵大鹏. 2007. 使用纵波和横波走时层析成像研究菲律宾海板块西边缘带和南海地区的深部结构. 地球物理学报 , 50(6): 1757–1768. | |
| 孙建国. 2007. 复杂地表条件下地球物理场数值模拟方法评述. 世界地质 , 26(3): 345–362. | |
| 孙建国, 蒋丽丽. 2009. 用于起伏地表条件下地球物理场数值模拟的正交曲网格生成技术. 石油地球物理勘探 , 44(4): 494–500. | |
| 奚先, 姚姚. 2001. 二维随机介质及波动方程正演模拟. 石油地球物理勘探 , 36(5): 546–552. | |
| 奚先, 姚姚. 2002. 随机介质模型的模拟与混合型随机介质. 地球科学——中国地质大学学报 , 27(1): 67–71. | |
| 姚姚, 奚先. 2004. 区域多尺度随机介质模型及其波场分析. 石油物探 , 43(1): 1–7. | |
| 张东, 傅相如, 杨艳, 等. 2009. 基于LTI和网格界面剖分的三维地震射线追踪算法. 地球物理学报 , 52(9): 2370–2376. | |
| 张美根, 贾豫葛, 王妙月, 等. 2006. 界面二次源波前扩展法全局最小走时射线追踪技术. 地球物理学报 , 49(4): 1169–1175. | |
| 朱生旺, 魏修成, 曲寿利, 等. 2008. 用随机介质模型方法描述孔洞型油气储层. 地质学报 , 82(3): 420–427. | |