三维复杂黄土塬区变加密不等距网格初至走时计算与特征分析

韩令贺， 孙章庆， 胡自多， 刘威， 韩复兴， 田彦灿

1 中国石油勘探开发研究院西北分院， 兰州 730020 2 中国石油大学(北京)地球物理学院， 北京 102249 3 吉林大学地球探测科学与技术学院， 长春 130026

0 引言

复杂黄土塬区在我国有较大的分布范围，尤其是在我国西北部的鄂尔多斯盆地和塔里木盆地西南部.这两个区域均蕴藏着丰富的油气资源.但是，由于经历了长期的风化、侵蚀、冲刷、切割，黄土塬区形成了塬、梁、坡、沟等复杂地貌，地表高程变化剧烈，黄土层厚度剧烈变化且黄土疏松、弹性差、速度低，表层结构复杂多变，因此黄土塬区一直被视为地震勘探的“禁区”(阎世信和吕其鹏，2002).这些复杂的地震地质条件使得复杂黄土塬区也面临着地震激发接收困难、地震资料信噪比低、噪声异常发育、相干和次生干扰严重、近地表建模困难、静校正问题严重、中深层成像困难或基本不成像等诸多挑战，因此对复杂黄土塬区地震勘探的相关问题进行研究具有重要的理论和实际意义(刘宝国等，2017).

黄土塬区的复杂地震地质条件通常会引起地震波场结构异常复杂，如：初至扭曲严重、面波散射发育、体波地形散射发育、黄土层内吸收衰减严重、黄土层内多次波发育、中深层反射能量很弱且非常不连续等(阎世信和吕其鹏，2002；刘宝国等，2017).对这些复杂波场结构进行深入解析是解决上述复杂黄土塬区诸多挑战的基础.然而，目前专门针对复杂黄土塬区地震波传播规律与低信噪比数据的形成机制的相关研究并不多见，这几乎也是一直以来制约复杂黄土塬区地震数据处理效果欠佳的一个瓶颈(阎世信和吕其鹏，2002；刘宝国等，2017).面对该难题，很难一蹴而就的全面解析该区域的复杂地震波场结构.不过，在它们中初至波是相对最重要、最易、最应先解析的地震波场成分，因为它对近地表速度建模和静校正，这两项位于相对处理前端又极其重要的地震数据处理环节具有重要的意义.鉴于此，本文将紧密围绕“三维复杂黄土塬分区初至走时计算与特征分析”这一主题展开，意在从地震波场运动学角度，分析复杂黄土塬区初至的基本特征，并为后续该区近地表建模和静校正等建立稳定、灵活且兼顾精度和效率的初至走时算法.

综上所述，为了精细刻画三维复杂黄土塬的复杂近地表介质特征，建立稳定、高效、高精度、灵活的初至波走时算法，并深入分析该介质条件下初至波的各方面特征，本文提出一种三维复杂黄土塬分区变加密不等距网格初至波走时计算方法，该方法采用三维分区局部变加密不等距网格精细剖分复杂黄土塬模型，采用一种三维迎风变网格双线性插值法精细计算初至波走时，最后基于此得出三维复杂黄土塬区初至波的一些运动学特征.下面将紧密围绕这几个方面的内容展开.

1 三维复杂黄土塬模型的网格剖分

三维黄土塬模型通常非常复杂，采用怎样的网格剖分才能精细地刻画该模型，并便于后续初至波走时计算非常重要.为了同时兼顾计算精度和效率，本文提出一种三维分区局部变加密不等距网格剖分三维复杂黄土塬区模型，具体步骤如下：①读入三维地表高程数据、计算空间范围、计算网格间距等信息；②根据步骤①读入的信息，在计算范围内，采用网格间距较大的均匀的正方体粗网格剖分整个计算空间；③对地表不做任何近似，让它根据高程直接落在步骤②的均匀的正方体网格线上，并去掉地表以上非必要网格，形成地表附近的不等距网格；④从地表处所在的网格单元开始向地下由浅至深，按照(2n,2n-1,2n-2,…,2)指数衰减的加密程度对步骤②的单个正方体粗网格，进行网格加密,其中n为加密因子，它用于控制网格的加密程度；⑤以初至走时计算时的震源点为中心点由内向外，按照(2n,2n-1,2n-2,…,2)指数衰减的加密程度，对步骤④的网格，实施震源附近的网格加密.

在上述复杂网格生成的过程中，步骤④和⑤最为关键，它们均涉及加密因子n的选取和加密范围的确定两个核心问题.其中，在实施步骤④的地表附近的网格加密时，加密因子n的选取是由最粗的网格间距的大小、地表处需要加密的程度和计算精度共同决定的，n越大地表附近加密程度越大，计算精度越高；关于加密范围，由于采用的是变加密策略，网格的加密程度是渐变的，数值实验表明：由地表向下，上一级加密程度的网格的加密范围不超过两个下一级网格的网格间距，即可达到加密效果和数值计算的稳定，这也能最大程度的保证计算效率.此外，在实施步骤⑤的震源附近的网格加密时，加密因子n的选取是由最粗的网格间距的大小、震源附近需要加密的程度和计算精度共同决定的；而加密范围则是根据震源附近需要控制的计算误差来确定的，较低误差要求的范围越大则加密的范围越大，达到上一级加密程度的误差控制范围后，再过渡到下一级加密网格，直到过渡到常规均匀粗网格.

以上实现步骤最终获得了如图1所示的精细剖分三维复杂黄土塬区的分区局部变加密不等距网格.与其他网格相比，此处采用的网格有如下优缺点：①地表附近的不等距网格，不对地表高程做任何近似，完全准确的保留地表高程点的实际位置；②从地表开始向地下由浅入深的变加密网格，能更加细致的描述复杂黄土塬的地表形态和近地表复杂介质分布(尤其是薄风化层)，同时还有利于保证复杂地表处初至波走时计算的稳定性和精度；当然，该加密网格一定程度上会增加计算时的存储空间和计算量，但是由于实施的仅是地表附近的变加密，故仅会增加非常有限的存储空间和计算量；③由震源处出发由内而外的变加密网格，能够在不过多增加计算量的前提下，很好的解决常规初至波走时算法在震源附近误差较大的问题；当然，这种误差的控制是通过大幅度减少网格间距来实现震源附近截断误差降低的，虽然它不能完完全全从本质上解决震源奇异性的问题，但这种变加密策略是一种仅花费很少计算量的增加来大幅度提高计算精度的高性价比算法；④与曲线网格相比，上述网格剖分方法无需单独花费网格生成的计算量而实现对不规则边界不做任何近似处理；⑤在编程实现时，只要合理的处理好各个加密级别之间的调用和重新初始化，仅需在地表附近单独开辟用于存储加密节点的数组空间，而在震源加密时也仅需在一个数组中来回采样和反复初始化赋值多次，即可完成不同加密程度的变加密运算；⑥因为采用局部加密，所以算法还能保证后续初至波走时计算的绝大部分计算工作量在规则均匀的正方体网格中进行.

图1 三维复杂黄土塬模型的网格剖分Fig.1 Mesh generation of 3D complex loess plateau model

2 面向三维复杂黄土塬区的初至走时计算

上述采用了三维分区局部变加密不等距网格剖分三维复杂黄土塬模型，在该复杂网格中实现面向三维复杂黄土塬区的初至波走时计算，核心需要解决复杂网格中的局部走时算法和复杂网格中的整体波前扩展方式这两个问题，下面将详细阐述这两方面问题.

2.1 复杂网格中的局部走时算法

如图2a所示，图1中的三维分区变加密不等距网格，实际上包括了如图2b—e所示的四种由简单到复杂的情况：①网格间距没有变化的均匀立方体网格(图2b)；②网格间距有变化的不均匀过渡区立方体网格(图2c)；③邻近地表附近的复杂网格(图2d)；④地表处的复杂网格(图2e).这4种不同的网格情况，需要采用不同的局部走时计算公式.目前，关于这方面的三维走时计算问题，有学者曾提出过一种迎风混合法(孙章庆等，2017)，该方法采用不等距差分公式来简洁的处理地表附近的走时计算问题，但是相比于统一的迎风双线性插值法(孙章庆等，2015)，该方法在适应上述变加密不等距网格时相对困难，且地表附近的有限差分法精度相对有限.针对该问题，在此提出适应上述复杂网格的迎风变网格双线性插值法.该方法的核心思想是：首先基于Fermat原理，通过采用无条件稳定的迎风思想，在封闭被算点的邻近所有网格单元中，挑选走时分布最小的网格单元，然后再在该网格单元上采用双线性插值法进行局部走时计算公式的建立.在具体局部走时计算公式建立时，迎风变网格双线性插值法可以将上述提及的4种情况简化为如图2f—h所示的3种情况：

(1)非地表附近的规则立方体网格(图2f)中的局部走时计算公式

如图2所示，图2b、c对应的局部网格形态最终都有可能通过基于Fermat原理的迎风思想的挑选，简化为如图2f所示的规则立方体网格.网格中G点或N点的走时TG or N为需要计算的走时值.阴影部分网格单元MACB为基于Fermat原理的迎风思想，从图2b、c中封闭G点或N点的24个正方形网格单元中挑选出的，走时分布为最小的网格单元，其端点M点、A点、C点、B点的走时值TM、TA、TC、TB为已知.在这种情况下，可以直接基于双线性插值法(孙章庆等，2015)建立局部走时计算公式：

(1)

其中，s为当前计算网格单元的平均慢度(速度的倒数)，h为正方体网格的网格间距.公式(1)给出了在一个网格单元，其网格端点走时值均为已知且网格单元内走时为线性分布假设的情况下，计算其邻近点走时值的局部计算公式.但是，与常规的双线性插值法(孙章庆等，2015)不同，为了无条件稳定的适应如图1所示的复杂变化的网格，在此提出相应的迎风变网格双线性插值法，其核心是在复杂变化的网格中如何确定公式(1)建立时的已知条件，即如何确定如图2f所示的网格单元MACB及其中的端点M点、A点、C点、B点的走时值TM、TA、TC、TB.确定图2f所示的网格单元实际上涉及两种情况：①如图2b所示，当G点位于均匀立方体网格中时，根据Fermat原理M点应为与G点直接相邻的6个网格节点中走时值最小的网格节点，A点和B点则分别为与M点直接相邻的另外两个坐标轴方向上走时值较小者，即TM、TA、TB需满足：

图2 迎风变网格双线性插值法(a) 三维复杂地表附近的局部网格； (b—e) 不同位置的网格形态；(f—h)简化后的计算网格.Fig.2 The method of upwind variable grid bilinear interpolation(a) Local grid nearby the 3D complex earth′s surface; (b—e) The shapes of grid at different locations; (f—h) Simplified computational grid.

TM=min(TMi)i=1,2,…,6,

(2a)

TA=min(TMix+,TMix-),TB=min(TMiy+,TMiy-),i=1或6,

(2b)

TA=min(TMiy+,TMiy-),TB=min(TMiz+,TMiz-),i=2或4,

(2c)

TA=min(TMix+,TMix-),TB=min(TMiz+,TMiz-),i=3或5.

(2d)

公式(2a)表述了M点为与G点直接相邻的6个网格节点中走时值最小的网格节点的含义，而公式(2b)—(2d)则表述了A点和B点分别为与M点直接相邻的另外两个坐标轴方向上走时值较小者(例如：如图2b所示，若M=M6时，TA和TB则分别可以表示为：TA=min(TM6x+,TM6x-),TB=min(TM6y+,TM6y-)，其中：TM6x+,TM6x-分别为x轴方向上M6邻近的点M6x+和M6x-的走时值，TM6y+,TM6y-则为y轴方向上M6邻近的点M6y+和M6y-的走时值)；②如图2c所示，当N点位于网格间距有变化的不均匀立方体网格时，同样根据Fermat原理来确定网格单元MACB及其中的端点M点、A点、C点、B点的走时值TM、TA、TC、TB，但是此时情况变得复杂了很多，其中M点应为与N点直接相邻的10个网格节点中走时值最小的网格节点，而A点和B点的确定方法已有所变化，TM、TA、TB需满足：

TM=min(TMi)i=1,2,…,10,

(3a)

TA=min(TMix+,TMix-),TB=min(TMiy+,TMiy-),i=1或10,

(3b)

TA=min(TMiy+,TMiy-),TB=TM1x+orTM1x-,i=2或4,

(3c)

TA=min(TMix+,TMix-),TB=TM1y+orTM1y-,i=3或5,

(3d)

TA=min(TMiy+,TMiy-),TB=TM10x+orTM10x-,i=6或8,

(3e)

TA=min(TMix+,TMix-),TB=TM10y+orTM10y-,i=7或9.

(3f)

公式(3a)表述了M点为与G点直接相邻的10个网格节点中走时值最小的网格节点的含义；而公式(3b)与公式(2b)类似，它表述了当M=M1or M10时，A点和B点分别为与M点直接相邻的另外两个坐标轴方向上走时值较小者；但是，公式(3c)—(3f)则不同，此时A点和B点分别为与M点直接相邻的另外两个坐标轴方向上的邻近点，但是其中一个方向为两个邻近点的走时值较小者，而另外一个方向上则仅有一个选择(例如：如图2c所示，若M=M6时，TA和TB则分别可以表示为：TA=min(TM6y+,TM6y-),TB=TM10x+，其中：TM6y+,TM6y-分别为y轴方向上M6邻近的点M6y+和M6y-的走时值，而在y轴方向上TB只能取TM10x+).

至此就建立了规则立方体网格(图2f)中的局部走时计算公式.该公式看似和常规双线性插值走时计算的公式区别不大，但是该公式在其已知条件的确定时，则隐含着基于Fermat原理的迎风思想；同时，该思想还能灵活的将网格间距变化的复杂网格处的局部走时计算策略简化到规则立方体网格中，并保证局部算法的无条件稳定性.

(2)邻近地表的不规则网格(图2g)中的局部走时计算公式

如图2d所示，当计算地表邻近点H点的走时值时，其也是处于一种复杂的不规则网格中.此时同样首先采用基于Fermat原理的迎风思想，挑选封闭H点的所有网格单元中走时分布最小的网格单元.但是，这种情况比图2b、c两种情况要复杂一些，因为这里除了考虑与H点直接相邻的网格点外，还需要考虑邻近地表点的影响，所以此处先给出确定M点走时值TM的表达式，再根据TM的取值不同，给出该复杂局部网格中的走时计算公式.在该复杂局部网格中，M点的走时值TM应该满足：

TM=min[(TMi)i=1,2,…,6, (TSi)i=1,2,…,8]

(4)

公式(4)实际上表述了M点应为所有与其相邻的网格点或地表点中走时值最小的节点.如图2d所示，当M点为不同情况时，TH的计算公式将不一样：①当M=M1时或M=(Mi)i=2,3,4,5且A=M1x+,M1y-,M1x-orM1y+，则图2d的复杂网格可以简化成图2f的规则立方体网格中的局部走时计算公式，此时根据双线性插值方法，直接采用公式(1)即可；②当M=(Mi)i=2,3,4 or 5且A=S1,S3,S5orS7时，此时图2d的复杂网格可以简化成图2g所示的不规则网格，同样采用双线性插值法，可得TH的计算公式：

(5)

其中：s为当前计算网格单元的平均慢度(速度的倒数)，h为正方体网格的网格间距，d为M点到A点的距离.③当M=M6or (Si)i=1,2,…,8时，即M点为地表上的点时，则直接根据Huygens原理，将M点视为子震源点，则有：

TH=TM+sd，

(6)

其中：d为M点到H点的距离.

至此就建立了邻近地表的不规则网格(图2d)处的局部走时计算公式，建立时首先需要基于Fermat原理的迎风思想给出确定M点的公式(4)，然后分别根据公式(4)的挑选结果，将此时的复杂网格简化为图2f的规则立方体网格、图2g的不规则网格或地表上点的直达波这3种情况中的一种，进而对应的采用公式(1)、(5)或(6)作为计算地表附近邻近H点的局部走时计算公式，该公式的建立也是综合利用迎风思想、Fermat原理和Huygens原理的结果，其能灵活地适应地表附近的复杂网格并能保证局部算法的无条件稳定性.

(3)地表处的不规则网格(图2h)中的局部走时计算公式

如图2e所示，当计算地表上的S点的走时值时，其也是处于一种复杂的不规则网格中.此时同样首先采用基于Fermat原理的迎风思想，挑选封闭S点的所有网格单元中走时分布最小的网格单元.不过，此时情况变得略微简单点，因为如图2e所示，与S点直接相邻的网格点仅有M1和(Si)i=1,2,…,8点.因此首先根据基于Fermat原理的迎风思想，M点的走时值TM应该满足：

TM=min[TM1,(TSi)i=1,2,…,8].

(7)

同样当M点为不同情况时，TS的计算公式将不一样：①当M=(Si)i=1,2,…,8时，则可以直接根据Huygens原理将M点视为子震源点，采用公式(6)的形式计算TS，不同之处仅为此时的d为M点到S点的距离；②当M=M1时，则图2e的复杂网格可以简化成图2h的不规则网格，同样采用双线性插值法，可得TS的计算公式：

(8)

其中：s为当前计算网格单元的平均慢度(速度的倒数)，h为正方体网格的网格间距，d为M点到S点的距离.

至此就建立了地表处不规则网格(图2e)中的局部走时计算公式，建立时首先需要基于Fermat原理的迎风思想给出确定M点的公式(7)，然后分别根据公式(7)的挑选结果，将此时的复杂网格简化为图2h的不规则网格或地表上点的直达波这2种情况中的一种，进而对应的采用公式(6)或(8)作为计算地表附近邻近S点的局部走时计算公式，该公式的建立也是综合利用迎风思想、Fermat原理和Huygens原理的结果，其能灵活地适应地表处的复杂网格并能保证局部算法的无条件稳定性.

综上所述，为了计算复杂黄土塬区的复杂网格中的初至波走时，提出了一种迎风变网格双线性插值法，该方法首先通过基于Fermat原理的迎风思想，将复杂网格简化为图2f—h所示的3种网格，并综合利用双线性插值法和Huygens原理建立了这些网格中的局部走时计算公式.该方法具有如下优点：①采用复杂网格剖分复杂黄土塬模型，对地表高程没做任何近似处理，并在地表附近采用加密的网格用于精细刻画地表形态，进而提高对不规则计算边界刻画的精细程度；②上述局部走时计算公式的建立方法，能够灵活的将复杂网格中的走时计算问题，在迎风思想、Fermat原理和Huygens原理的综合应用下简化为3种网格中的双线性插值问题；③算法能保证在复杂网格中的任意网格位置均采用计算精度相对较高的双线性插值法，进而保证整个计算的精度；④迎风思想、Fermat原理和Huygens原理的综合应用能够保证算法的灵活性和无条件稳定性；⑤上述变加密网格仅在地表附近和震源附近进行网格加密，因此能够在不增加过多计算量的前提下很好地解决震源附近和不规则计算边界处误差较大的问题.

2.2 复杂网格中的整体波前扩展方式

以上给出了局部走时算法，它可以解决在复杂网格的不同情况下的局部走时计算问题.但是，要想完成整个计算区域的初至走时计算，还需要一种整体的波前扩展方式.为了很好的解决该问题，在此引入群推进法(Kim and Folie, 2000).该方法在三维情况下有很高的计算效率，同时只需稍加改进就能灵活的适应复杂地表的问题.如图3所示，常规群推进法将计算区域内的网格节点的属性(用Gid的数值表示)分为三类，即Gid=0、1、2三种(图3中分别用白色、灰色、黑色充填的原点表示).其中，Gid=0表示未开始走时计算的网格节点，Gid=1表示当前波前上的网格节点，Gid=2表示已经完成走时计算的网格节点.

图3 三维复杂黄土塬条件下修正后的群推进法Fig.3 Modified group marching method in 3D complex loess plateau condition

以上群推进法的实现过程，实际上是一个通过不断改变网格节点属性，来近似模拟波前传播的过程，不断纳入波前点的运算表征着波前不断向前传播到达新的网格节点的过程，不断移除波前点的过程表征着波前到达过的区域为走时计算完成的区域.上述波前扩展过程是同时选取波前上的很多网格节点作为子震源点，这与Huygens原理是完全吻合的.此外，为了适应复杂地表问题，只需永久性的将地表以上的网格节点的属性设为Gid=2，走时值设为无穷大，即可实现地震波无法穿越地空界面进入空气中物理事实.

综上所述，综合采用上述复杂网格中的局部走时算法和整体波前扩展方式，即可无条件稳定、在不过多增加计算量的前提下高精度的计算复杂黄土塬区的初至波走时.当然，该算法能否被用于分析复杂黄土塬区初至波特征，还需要对算法进行精度和效率评估后才能决定，因此接下来首先对算法的精度与效率进行分析.

3 计算精度与效率分析

为了对比分析本文提出的新算法的计算精度与效率，采用如图4a所示的黄土塬模型.模型的大小为1.6 km×0.8 km×1.2 km，模型的地表部分由两个半径为0.4 km的半球形黄土塬(便于获得精确的地表上任意一点的解析走时值)构成，地下为均匀介质，地震波在介质中的传播速度为1.0 km·s-1.计算时，震源点置于(0.8 km，0.4 km，0.4 km)处.图4b给出了网格加密策略：在离震源最近的半径为80.0 m的球体内网格间距为2.5 m，依次向外进行球形扩展，在径向距离分别为80.0～160.0 m和160.0～320.0 m的球环内，计算的网格间距分别为5.0 m和10.0 m，然后在320.0～1600.0 m的非加密区域，网格间距过渡为均匀的20.0 m.上述震源附近的加密策略，相当于1.2节中阐述的(2n,2n-1,2n-2,…,2)指数衰减的加密策略取加密因子n=3的情况.此外，地表附近则直接采用1.2节阐述的加密策略.

图4 计算精度与效率(a) 模型； (b) 网格加密策略.Fig.4 The computational accuracy and efficiency(a) Model; (b) The strategy for increasing the density of grid.

图5a—c和图5d分别给出了当单独采用网格间距20.0 m、10.0 m、5.0 m和采用上述变加密网格剖分整个模型时，地表上所有网格节点走时计算相对误差的分布，对比分析可知：随着网格间距的减小，走时计算的相对误差逐渐减小(图5a—c)，但是即使当网格间距减少到5.0 m时，整个计算区域计算误差均存在震源点附近误差较大的问题；不过，通过采用本文的网格加密策略，却能很好的解决震源附近误差较大的问题(图5d).

图5 地表上的计算精度(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h为变加密.Fig.5 The computational accuracy of the earth′s surface(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h is variable density.

图6和图7分别给出了Y=0.4 km和X=0.8 km剖面上，各个网格点的走时计算误差.与地表上分析的结果类似，在地表以下的剖面上，同样存在：随着网格间距的减小，走时计算的相对误差逐渐减小(图6a—c和图7a—c)；而本文采用的变加密网格方法(图6d和图7d)，可以很好的解决震源附近误差较大的问题.此外，同时分析图6和图7的计算误差分布，可以发现一个有趣的现象：走时计算误差的分布具有明显的方向性.其中，在震源向外辐射的0°、45°、90°、135°、180°、225°、270°和315°的八个方向及附近区域误差相对较小，而在这八个方向中的相邻的两个方向的平分线方向及附近区域误差相对较大.产生这种现象的原因，初步分析为计算局部走时的计算公式在均匀介质特定方向为解析解无误差.

图6 Y=0.4 km剖面上的计算精度(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h为变加密.Fig.6 The computational accuracy of Y=0.4 km section(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h is variable density.

图7 X=0.8 km剖面上的计算精度(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h为变加密.Fig.7 The computational accuracy of X=0.8 km section(a) h=20.0 m； (b) h=10.0 m； (c) h=5.0 m； (d) h is variable density.

表1给出了计算精度和计算效率对比的定量化统计数据，对比分析可以发现：①无论在哪个统计区域(地表上、Y=0.4 km剖面或X=0.8 km剖面上)，当网格间距减小时平均相对计算误差均减小；②采用上述变加密网格策略获得的计算精度(平均相对误差)，介于网格间距单独取10.0 m和5.0 m之间；③变网格策略的计算耗时(计算的硬件设备参数为：Intel(R) Xeon(R) CPU E5-2620 0 @2.00 GHz，32 G运行内存)仅为h=10.0 m和h=5.0 m时的19.22%和1.53%.综上所述，本文的变加密网格策略能在仅增加少量计算量的前提下，大幅提高计算精度，并能很好的解决震源附近误差较大的问题.

4 三维复杂黄土塬区初至特征分析

黄土塬区地震地质条件非常复杂，地表高程剧烈变化(沟、坡、梁、塬等复杂地貌并存)、表层黄土疏松且风化严重、黄土层厚度较大且变化剧烈、表层结构复杂多变、巨厚黄土层内存在低速层、降速层和高速层、黄土层内随机非均质性严重、部分区域存在明显的砾石层散射.这些复杂地震地质条件给黄土塬区地震勘探带来了巨大的挑战.为获得三维复杂黄土塬区地震初至波走时特征，利用本文所提出的方法计算了三维水平地表层状模型和三维复杂黄土塬模型的地震波走时场.

图8a给出了三维水平地表层状介质模型的形态和震源点的位置，其中：模型的大小为1.2 km×1.2 km×2.55 km，由浅及深三层介质中地表、界面1和界面2的深度分别为0.15 km、0.37 km和1.4 km；三层介质中地震波的传播速度依次为0.8 km·s-1、2.0 km·s-1和3.0 km·s-1，计算时震源位于(0.0 km，0.0 km，0.15 km)处，震源和地表附近的网格加密策略采用第3节“计算精度与效率分析”中的网格加密方式.图8b—e分别给出了地表上、Y=0.0 km剖面上、X=0.0 km剖面上和X-Y=0.0 km剖面上的走时分布.

分析图8b—e可以得出在三维层状介质模型中初至波的一些特征：①在地表上接收的初至波可以分为直达波和折射波两个区域(图8b中半径约为0.6 km的1/4红色虚线圆弧分开的两个区域)，其中直达波等时线间隔距离远小于折射波等时线间隔距离(图8b)；②由于模型为层状介质，其在横向上是均匀分布的，因此初至波等时线在地表上是以同心圆形式存在的(图8b)；③在地下介质中，由于由浅至深地震波速度是逐层增加的，因此在地表以下两个界面处都会产生折射波，不过由于各层间的速度比不一样，其产生折射波的临界角也有所差别，第1个界面处临界角较小，其先产生折射波.

图9a给出了三维复杂黄土塬模型的形态和震源点的位置，其中：模型的大小为1.2 km×1.2 km×2.55 km，模型的地表起伏形态是大致根据一个实际黄土塬工区近地表调查的黄土层厚度生成的(图9b)；近地表速度结构复杂，由浅及深，分别由风化层、砾石散射层、低速层、降速层和高速层构成，其中：风化层、低速层和降速层中地震波的传播速度分别为：0.5 km·s-1、0.9 km·s-1和1.5 km·s-1，砾石散射层由一个随机介质构成，高速层参考了一个塔西南黄土塬区深层的复杂构造模型.整体上分析，该三维复杂黄土塬模型在大尺度上考虑了复杂构造的非均匀性，在几何形态上考虑了各种复杂地形，在小尺度上考虑了砾石层的随机扰动，因此其全面细致的描述了三维复杂黄土塬区的复杂地震地质条件.计算时，震源位于(0.0 km，0.0 km，0.35 km)处.图9c—f分别给出了地表上、Y=0.0 km剖面上、X=0.0 km剖面上和X-Y=0.0 km剖面上的走时分布.

图9 三维复杂黄土塬模型中的初至走时分布(a) 模型； (b) 地表高程的等高线； (c) 地表上走时的等值线在平面上投影； (d) Y=0.0 km剖面上的走时分布； (e) X=0.0 km剖面上的走时分布； (f) X-Y=0.0 km剖面上的走时分布.Fig.9 Travel-time distribution of first arrival in 3D complex loess plateau model(a) Model； (b) The contour of the earth surface′s elevation； (c—f) Travel-time distribution on the earth′s surface，on the section of Y=0.0 km，on the section of X=0.0 km，on the section of X-Y=0.0 km.

分析图9c—f可以得出复杂黄土塬模型中初至波的一些特征：①同样，在地表上接收的初至波可以分为直达波和折射波两个区域(图9c中半径约为0.06 km的1/4红色虚线圆弧分开的两个区域)，其中直达波等时线间隔距离远小于折射波等时线间隔距离(图9c)；但是，由于极浅风化低速层的存在，地表上直达波覆盖的区域很小，很快地表上就接收到了折射波；②由于复杂黄土塬区地形的剧烈起伏，地表高程变化引起了较大的初至走时时差，因此地表上的初至走时等时线出现了很多极值点(图9c中的一系列近似同心圆的圆心)；③由于近地表极浅层低速风化层的存在，该薄层中折射波发育，等时线的间距很小(图9d—f地表附近的等时线)；④砾石层中介质速度的随机扰动，对初至波的等时线有明显的随机改造，等时线有明显的局部随机扰动(图9d—f中随机扰动层的等时线)；⑤同样，在地下介质中，由于至浅而深地震波速度是逐层增加的，因此在地表以下速度剧变的界面处都会产生折射波(图9d—f)；⑥由于大体上地下介质速度是随深度的增加而增加的，因此初至波等时线的间隔，由地表向下是逐渐增加的(图9d—f).

对比分析图8和图9可以发现：①与三维水平地表层状介质模型相比，三维复杂黄土塬模型地表接收到的初至走时的直达波等时线形态在近偏移距同样为同心圆，但是随着偏移距增加其形态受近地表黄土层厚度横向变化快和地下速度结构复杂的影响较大，其等时线分布非常复杂，到达中远偏移距折射波的等时线形态与地形的等高线有一个很好的对应关系：地形高程的高值闭合区(图9b)对应着等时线的高值闭合区(图9c)，这与地震波传到黄土塬的顶部需要花更多走时的物理事实是相吻合的；②与三维水平地表层状介质模型相比，由于近地表薄低速风化层的存在，复杂黄土塬模型中地表初至成分里直达波覆盖的区域比水平层状模型的要小很多，前者的覆盖范围仅为半径约为0.06 km的1/4红色虚线圆弧圈定的扇形区(图9c)，而后者的覆盖半径约为0.6 km(图8b)；③即使是在如图9a所示的三维复杂黄土塬模型中，本文算法也能保证很好的稳定性，能很好的适应地形剧烈起伏、风化层、砾石散射层、低速层、降速层和高速层等近地表及地下复杂介质(图9c—f).

如图10所示，为了进一步对比分析三维复杂黄土塬模型中的初至特征，分别提取了三维水平地表层状介质模型各个剖面地表上的时距曲线，以及三维复杂黄土塬模型各个剖面地表上的时距曲线，对比二者可以发现：①三维水平地表模型的时距曲线是典型的折线形态，每段折线的斜率与各层介质的速度相关，其中直达波和折射波明显分开，但是三维复杂黄土塬模型的时距曲线则是不规则的曲线，其中也能明显看出直达波的直线段，但是其形态也一定程度受到了地表形态和近地表复杂介质的影响；②三维复杂黄土塬模型地表上的时距曲线形态虽然为不规则的曲线，但是它们的形态和地表高程线的形态之间有一定的对应关系，而三维水平地表层状介质模型则没有；③同样，图10中红色的虚线非常明显的划分开了直达波和折射波的时距曲线的区间，三维复杂黄土塬模型的直达波覆盖的范围远远小于三维水平地表层状介质模型，这主要是三维复杂黄土塬模型近地表薄风化层带来的影响；④二者形态的对比可以预料到，与常规水平地表地下分层很强的区域相比，三维复杂黄土塬区的初至拾取会因为时距曲线的严重不规则而面临着很大的挑战，同时利用该初至时距曲线进行初至层析获取近地表速度模型时也将面临着很大的挑战.

图10 地表上时距曲线的对比(a)、(c)、(e)三维水平地表层状介质模型的Y=0.0 km、X=0.0 km、X-Y=0.0 km剖面的时距曲线；(b)、(d)、(f)三维复杂黄土塬模型的Y=0.0 km、X=0.0 km、X-Y=0.0 km剖面的时距曲线.Fig.10 The comparison of time-distance curves on the earth′s surface(a),(c),(e) The time-distance curves of the 3D model with horizontal earth′s surface and layered medium on the section of Y=0.0 km，X=0.0 km and X-Y=0.0 km；(b),(d),(f) The time-distance curves of the 3D complex loess plateau model on the section of Y=0.0 km，X=0.0 km and X-Y=0.0 km.

5 结论与讨论

为了快速稳定的计算三维复杂黄土塬区初至波走时并分析其特征，提出一种群推进迎风变网格双线性插值法；通过对该方法的理论分析、计算精度和效率以及计算实例的分析，可以得出如下结论：

(1)采用三维分区局部变加密不等距网格剖分三维复杂黄土塬模型，具有精细刻画地表形态、提高地表及震源附近计算精度、在增加很少计算量的前提下大幅提高整个模型区域的走时计算精度、保证整个计算区域稳定性等优势；

(2)计算精度和计算效率对比分析表明：为了计算三维复杂黄土塬区复杂网格中的初至波走时而提出的一种迎风变网格双线性插值法，具有能无条件稳定适应复杂地形和近地表及地下复杂介质、灵活稳定适应复杂网格、保证计算精度、在不过多增加计算量的前提下解决震源附近误差较大的问题等优势；

(3)计算实例分析表明：本文算法能稳定适应三维复杂黄土塬地震地质条件，同时还得出了三维复杂黄土塬区初至成分组成及分布特征，对于充分拾取和利用复杂黄土塬区的初至走时信息具有一定的理论意义和实际价值.