0 引言

天然气水合物不仅是一种资源储量巨大的潜在能源，而且对海底地质灾害、气候变化及环境等方面有重要影响 (张尧声等，2015)，因此，确定海洋天然气水合物在海底地层中的分布意义重大.地震勘探是海洋天然气水合物勘探的主要技术之一，BSR (Bottom Simulating Reflector) 作为指示天然气水合物的主要标志被广泛使用.但是，BSR与天然气水合物之间并没有一一对应的关系，有时不存在BSR时也会有天然气水合物存在 (刘培栋等，2016)，所以仅以BSR作为水合物存在与否的评判标准不一定可靠.另一方面，含天然气水合物地层的地震波速度会明显增加，而水合物下伏地层因含游离气或含水，地震速度则会有明显下降 (刘培栋等，2016).正是由于水合物层与下伏地层之间的这种阻抗差异，才使得地震剖面上出现了BSR.因此，精确的海底地层速度结构能为天然气水合物的存在提供更有力的证据，而且，地层速度剖面还可以反映天然气水合物储层的顶底界面，为圈定天然气水合物体积提供依据.

目前，地震勘探中建立速度模型的方法主要有叠加速度分析、偏移速度分析和层析成像三大类 (张志明等，2016)，相对于前两类速度建模方法，层析成像更能适应复杂地质构造，因而成为近年来重点关注的复杂介质速度建模方法 (金昌昆和张建中，2014).在国外天然气水合物勘探中，层析成像技术已得到了广泛应用.Singh等 (1993)应用单道全波形反演方法研究了温哥华近海BSR上下的速度结构，并被大洋钻探所证实；Tinivella和Accaino (2000)应用地震反射层析成像研究了南极洲南Shetlan群岛岸外、布莱克海台的纵波速度结构，并推断了水合物和游离气的分布；Thakur等 (2007)利用多道地震反射走时层析成像法反演了印度西南部Kerala-Konkan海域的二维速度结构，预测了该海域天然气水合物的存在，并估算出了水合物层和游离气层的厚度；Cheng等 (2010)利用反射走时层析成像反演了台湾西南部海域的三维纵波速度结构，并估计了水合物的饱和度.在国内应用层析成像进行水合物预测起步较晚.宋海斌等和徐宁应用一维波形层析反演，分别建立了日本东南海海槽双BSR的对应地层、东海冲绳海槽南部水合物沉积层的一维速度结构，并根据反演的速度结构，对相应区域的天然气水合物储层进行了预测 (宋海斌等，2003；徐宁，2007).但是，沿深度变化的一维速度结构，不能反映地层速度的横向变化.二维波形层析反演方法虽然理论上精度高，但因影响波形因素多，子波估计不准、反演多解性严重、计算量大等问题，实际应用效果受到限制.与之相比，反射走时层析成像使用比较可靠的走时数据，能够建立地下速度模型.另一方面，相对于石油勘探，海洋水合物储层位于海底数百米的浅层，目标区地震同相轴往往比较清楚，可以拾取不同层位的反射走时，使地震走时层析成像具有应用条件.但是，反射走时层析成像依然存在诸如反射同相轴的识别和走时的拾取、同时反演速度和反射界面深度不同物理量的处理等诸多问题 (刘玉等，2014).本文应用人机交互多域拾取方法 (马飞等，2015) 拾取反射走时；采用归一化技术解决了同时反演速度和反射界面深度时不同物理量量级不统一的问题；在反演中引入正则化和平滑处理，有效提高了反演过程的稳定性；并利用本文研发的基于不规则单元的二维反射走时层析成像方法，对南海北部陆坡神狐海域的地震数据进行了相应处理，建立了沿地震测线地层的2D速度模型，结合钻井和测井资料，确定了天然气水合物储层.

1 反射走时层析成像方法

根据地震射线理论，走时t与慢度s有如下关系：

(1)

其中，t为地震波走时，L(s) 为射线路径，s(x, z) 为地层慢度，dl为沿射线路径的距离增量.由于射线路径与慢度有关，上式是一个非线性关系式.

1.1 模型离散化

为了对 (1) 式进行正反演数值计算，需要对地层速度模型进行离散化.对于反射层析成像，既要考虑每个地层速度的纵向和横向变化，还要考虑地层边界 (即反射界面) 的起伏.本文具体做法是：首先在x方向上，用平行于z轴的直线等间距地把模型进行剖分，在x方向的各个剖分网格内，地层界面用直线近似，再在每一个地层内根据需要设置一定层数的网格，将该层在z方向上均分.这样，就把每个地层划分成x方向长度相等、z方向厚度各异的不规则单元.各个单元左右边界平行，上下边界倾斜 (张建中和丁兴号，2007)，如图 1所示.用各单元角点的速度值表示模型的速度分布，用地层界面上单元角点的深度表示地层界面的变化.这种离散化不仅能够很好地拟合起伏的反射界面，较好地考虑每个地层内速度的纵横向变化，而且易于计算机自动实现.

1.2 层析反演

本文采用高斯-牛顿法，迭代地求解反射走时层析反演问题.在每步迭代中，在初始模型的基础上，对式 (1) 进行离散化，同时考虑地层速度和反射界面深度量纲和量级不同，对反演方程进行了归一化处理，这里，对系数矩阵中与深度有关的元素乘一系数，使各元素量级相近，从而得到反射走时残差与地层慢度修正量和反射界面深度修正量满足的线性层析反演方程为

(2)

其中，Δt=[Δt₁, …, Δt_m]^T，为m维走时残差向量；Δx=[Δs₁, …, Δs_n, Δd₁, …, Δd_r]^T为模型校正向量，Δs_j(j=1, 2, …, n) 为第j个单元的慢度校正量，Δd_k(k=1, 2, …, r) 为第k个反射界面节点的深度校正量，m是观测到的反射走时数目，n是模型离散单元数目，r是所有反射界面上节点的总数，它等于反射界面数与模型在x方向节点数的乘积.C_refl为m×(m+n) 阶矩阵，C_refl中元素(i=1, 2, …, m，j=1, 2, …, n) 为第i个反射走时对第j个单元慢度的微分，等于第i条射线在第j个网格单元中的长度；表示第i个反射走时对第k个反射界面节点深度的微分，且(i=1, 2, …, m；k=1, 2, …, r)(Ditmar等，1999)，s_R为反射点上部地层的慢度，β为反射波入射角，θ为反射面与水平面的夹角，x_R为反射点的横坐标，x_k-1和x_k分别为反射点左右两侧节点的横坐标，如图 2所示.ε是归一化系数，其作用是使走时对慢度的微分与走时对深度的微分的数量级相同，由此方程求出的深度为(k=1, 2, …, r).

同时反演速度和深度的反演问题病态严重，为了使层析反演稳定收敛，得到可靠的反演结果，需要对慢度和反射界面深度分别施加光滑约束以及先验速度约束.

慢度光滑约束.根据最平坦解的思想，对模型横向和纵向节点的慢度取一阶导数正则化，建立的慢度正则化方程为

(3)

其中，为 (n-1)×(n+m) 阶横向慢度正则化矩阵；B₂=为 (n-1)×(n+m) 阶纵向慢度正则化矩阵，设模型横向网格数目为c，则B₂每行元素中-1和1所在的列数相差c；，Δs₁=[s₁⁰-s₂⁰, …, s_n-1⁰-s_n⁰]^T，Δs₂=[s₁⁰-s_c+1⁰, …, s_n-c⁰-s_n⁰]^T，Δs₁和Δs₂都为n-1维列向量，它们的元素s_j⁰(j=1, 2, …, n) 是第j个节点上的初始慢度或前次迭代的慢度.

反射界面深度光滑约束.类似地，根据最平坦解的思想，对模型反射界面的深度取一阶导数正则化，建立模型反射界面深度正则化方程，公式为

(4)

其中，为 (m-1)×(n+m) 阶矩阵；Δd=[(d₁⁰-d₂⁰)/ε, …, (d_m-1⁰-d_m⁰)/ε]^T，为m-1维向量，它的元素d_k⁰(k=1, 2, …, r) 是第k个反射界面节点上的初始深度或前次迭代的深度.对于不属于反射界面的节点，只需设矩阵C_depth中双对角阵对应元素为0即可.

先验速度约束.先验速度信息可归纳为线性方程组形式公式为

(5)

其中，C_prior=[D₁, D₂]为n×(n+m) 阶矩阵，D₁为n×n阶单位矩阵，D₂为n×m阶零矩阵，Δs_prior为n维列向量，其元素为s′_i-s_i⁰(i=1, 2, …, n)，s′_i为先验慢度值，s_i⁰为初始慢度值或前次迭代的慢度.当节点i处有先验慢度时，矩阵D₁中节点i对应的元素赋为1，否则赋为0.

联立 (2) 到 (5) 式就构成了每步迭代的约束层析反演方程，为了调整各个方程的作用大小，增加了各方程的权重系数，最终形成的层析反演方程为

(6)

其中λ_refl、λ_s、λ_d、λ_prior分别表示反射波走时、模型慢度、反射界面深度和先验速度的权重系数，它们也起到了对方程之间进行归一化的作用.

计算当前模型的反射射线路径和走时，建立上述反演方程，并用LSQR求解反演方程，得到模型修正量Δx，再对模型进行修改，如此迭代进行，直到计算走时与观测走时之间的误差达到设定的条件为止.

1.3 正演计算

在反射走时层析反演的每一次迭代中，都要确定反射射线路径，计算反射走时，并与实测走时进行对比.本文采用Huang等 (2011)提出的方法，建立了地层速度单元和反射界面深度完全耦合模型的二维地震反射射线追踪和反射走时计算方法.射线追踪通过先确定波前走时、再根据波前走时确定射线路径两个步骤来实现 (Zhang et al., 2004).第一步，计算波前走时.在前述用不规则单元对模型离散化的基础上，采用不规则单元双线性走时插值和GMM波前扩展算法 (Zhang et al., 2011) 计算所有网格节点的波前走时.对于反射波，先从源点开始，向下计算波前走时至反射界面，再将反射界面上波前走时最小的节点作为新震源，向上计算各节点上的反射波前时间.第二步，根据波前走时确定射线路径.先利用上行波波前时间，从接收点开始，反向确定满足Fermat原理的射线，直至进行到反射界面为止，界面上的该点，即为对应的反射点；再利用下行波时间，从反射点开始，反向追踪反射点到震源点的射线，这样，连接震源点至反射点以及反射点到接收点的射线，便得到从震源点到接收点的反射射线路径.该方法具有很高的计算精度和对复杂介质的适应性.

2 理论模型数据测试

为了检验本文反射层析成像方法建立含天然气水合物地层速度模型的有效性，根据南海神狐地区钻探证实的天然气水合物储层情况，设计了理论模型，合成了该理论模型的地震数据，把层析反演结果与理论模型进行对比.设计的模型如图 3a所示.模型在x方向长9000 m，z方向长1600 m，由五层介质组成，第一层为水层，速度1500 m/s，平均厚度约600 m；第二层为含水合物的沉积层，速度1800 m/s，厚度约300 m，其中部的局部高速体为水合物，速度2100 m/s，厚度约80 m；第三层为含游离气层，速度1700 m/s，厚度约100 m；第四层为不含水合物和游离气的沉积层，速度2000 m/s，厚度约200 m；第五层为速度2500 m/s的高速地层.观测系统模拟海上二维拖缆观测系统，炮间距25 m，道间距12.5 m，炮检距125 m，总共328炮，每炮55道，2 ms采样，道长2 s.用声波方程模拟得到了地震数据，图 3b是偏移距为487.5 m的共偏移距记录.从图 3b中可以看到，在150~250道之间，在高速体底界对应的反射同相轴出现了明显的极性反转，即BSR.这是由于第三层速度低于其上覆含水合物地层的速度造成的.因此，在走时拾取的过程中要注意拾取极性的变化.

层析反演的初始模型为水平层状、横向速度均匀的五层速度模型，前四层速度为1500 m/s，第五层速度为2500 m/s，如图 3c所示.反演中取海水速度为1500 m/s，将其作为先验信息加入反演方程，第7次迭代的反演结果如图 3d所示.图 3e为层析反演模型与理论模型的速度差值剖面，图 3f为层析反演模型中各界面深度相对于理论模型各界面深度的误差分布图.可以看出，反演结果中各层的速度和厚度与理论模型吻合良好，水合物对应的位置有明显的高速异常，且位置和展布范围也与理论模型基本一致，说明该反射走时层析成像方法稳定性好，反演精度高，能够有效地给出海域含水合物地层的速度结构.

3 实际资料应用 3.1 研究区概况

研究区为南海北部陆坡神狐海域，位于珠江口盆地珠Ⅱ凹陷，水深600~3000 m.晚上新世以来，该区陆源碎屑供给充足，沉积物快速堆积，具有很好的成烃条件.受新构造活动影响，区域内断层、气烟囱、泥底辟等与水合物相关的特殊地质体极为发育，为水合物形成提供了运移通道 (苏明等，2015).2006年，中国地质调查局在神狐海域获取了高分辨率地震数据，地震剖面上有明显的BSR，BSR延伸长度一般为1.5~2.5 km，垂向分布深度一般处于海底以下150~350 m之间.2007年，在该区实施了天然气水合物钻探，并在三个站位获得了水合物样品，含天然气水合物层段位于海底以下153~225 m，厚度为18~34 m (苏正等，2014).自从神狐海域钻获水合物样品以来，该区的地震资料得到广泛研究 (李军红等，2015；张建中等，2015).我们把SH2井附近的含水合物沉积层作为目标区，利用走时层析成像方法建立速度模型，进行水合物储层预测.

3.2 地震资料处理及走时拾取

地震资料为神狐海域的一条过SH2井的二维测线数据，道间距12.5 m，炮间距25 m，共692炮，每炮192道，采样长度5 s，采样率2 ms.由于神狐海域海水深度较深、海底起伏大、地质构造复杂、特殊地质体发育，再加上在地震采集时由于海流影响电缆羽角较大，这些因素导致地震剖面上绕射波异常发育、反射同相轴连续性差，这给走时拾取造成了困难.为了准确地拾取反射走时，我们对原始地震资料进行了精细的去噪处理，并采用蒋波等 (2014)基于反射波层位拉平的绕射波分离技术对原始地震记录中的绕射波进行了压制.图 5a是研究区的偏移叠加剖面.

采用“多域联合拾取走时”的人机交互软件 (马飞等，2015) 拾取了反射波走时.由于地震数据在共炮点域、共中心点域和共偏移距域等不同域的同相轴特征及连续性不同，当在单一域下同相轴不够清晰，连续性不好时，在多个域交互拾取，相互补充、相互校正，大大提高了走时拾取的准确度.同一反射界面同相轴在各域的形态可能不同，但同一道上的同相轴位置在各域一定相同，以此可以验证在各域拾取走时的正确性.根据偏移叠加剖面水合物目标层以上同相轴的展布，我们选取了四个主要的同相轴进行走时拾取，图 4a和b分别是一个共炮域和一个共偏移距域拾取的四个同相轴.

3.3 初始模型建立

根据拾取的四个主要同相轴，建立了如图 5b所示五层介质的初始模型.各反射界面的深度这样获得：提取每炮第一道来自四个反射界面的双程走时，当作来自四个反射界面的自激自收时间，把该自激自收时间的二分之一乘以从海面到某个反射界面的平均速度，作为该界面的深度.第一层为海水层，速度设为1500 m/s，第二层到第五层的初始速度按SH2井声波测井速度的平均值设置，依次为1700 m/s、1900 m/s、1800 m/s和2500 m/s.这种做法近似地将叠加剖面的时间信息转化为地层反射界面的深度信息，在一定程度上比较接近实际反射界面.

3.4 层析反演

在x方向上等间距离散模型，在z方向上，同一地层按等单元数划分地层，各地层在z方向上的单元数可以不同.反演过程中模型离散化的具体参数：x方向网格宽度为20 m，在z方向上，第一层 (海水层) 分为76个单元，单元平均厚度约为20 m，第二层分为13个单元，单元平均厚度约为6 m，第三层 (含水合物层) 分为104个单元，单元平均厚度约为2 m，第四层 (含游离气层) 分为20个单元，单元平均厚度约为6 m，第五层由于未参与反演，未进行剖分.把海水速度作为约束，取1500 m/s，对反演方程进行正则化和归一化处理，第十次迭代反演结果如图 5c所示.从图中可以看出，在水平方向2000 m到3200 m的区间内，从海底向下，速度逐渐增大，到第三层下部速度达到最大，再向下速度减小，在第三层下部存在一个厚度约50 m左右的高速体，其速度为2000~2050 m/s，其下方存在一个低速层，速度为1830~1890 m/s之间，与BSR对应.在SH2井位置层析反演速度与测井声波速度曲线的对比如图 6所示，层析反演速度和测井声波速度的基本变化趋势一致.SH2井在海底以下186~229 m钻遇水合物层 (苏正等，2014)，结合钻探结果推断第三层下部的高速体为天然气水合物储层.

4 结束语

本文研发了一种基于不规则单元、同时反演地层速度和分界面深度的反射走时层析成像方法，采用平滑和先验信息约束、正则化和归一化技术，使反演过程更加稳定，反演效果得到了改善.利用多域人机交互方法拾取了南海北部陆坡神狐海域某地震测线的反射走时，应用反射走时层析成像获得了沿该测线的海底地层展布及其速度结构，与声波测井速度的变化趋势一致，并与钻探结果吻合.虽然层析反演速度不如测井速度的垂向分辨率高，但却能刻画整个剖面上的地层及其速度的变化特征，能为天然气水合物储层识别提供依据.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！