静校正问题一直是制约复杂地区地震勘探的瓶颈, 处理不好会严重影响地震资料的成像效果^[1]。针对复杂地区的静校正问题, 人们研究了折射静校正、层析静校正、波动方程静校正等各种技术, 在近地表变化较为复杂的地区往往要联合两种以上的静校正方法才能取得较好的效果^[2]。沙漠地区是我国西部油气勘探的主战场, 沙漠对地震勘探的主要影响表现在起伏剧烈的地表引起严重的静校正问题。由于沙漠地区近地表近似为连续介质, 因此有人提出了基于折射波初至的连续介质速度模型反演静校正方法^[3], 但是其应用效果并不理想。杨贵明等^[4]根据微测井测量得到的初至时间与深度的关系, 拟合得到近地表时深关系量板(时深曲线), 根据该量板及高速层顶面计算出所有炮点和检波点的静校正量。该方法不需要拾取记录初至, 特别适用于当前高密度地震采集, 自推出以来凭借其经济高效的优势在沙漠地区地震勘探中得到了广泛应用, 成为沙漠地区主要的静校正技术。1996年, 许亚军等^[5]在塔里木沙漠成功应用了时深曲线静校正技术。2010年, 张恒超等^[6]在ZGE沙漠地区对时深曲线静校正、模型静校正和折射静校正的效果进行了对比, 认为时深曲线静校正方法能较好地解决沙漠区的静校正问题。2015年, 尚新民^[7]在准中地区沙层地球物理性质研究的基础上, 综合应用时深曲线量板解决了该区沙丘引起的静校正问题。但是, 在一些沙丘成因复杂、表层物性差别大的地区, 特别是大沙漠的边缘地区, 根据简单的时深曲线获得的静校正量应用效果不理想。针对这一问题, 学者们提出了很多对策：王增明等^[8]将模型约束技术与时深曲线结合计算静校正量, 在准噶尔盆地沙漠地区取得了一定的应用效果, 但该方法仍然需要拾取初至, 而且实现过程较为复杂；肖泽阳等^[9]在塔里木沙漠地区对常见的几种静校正方法应用效果进行了对比, 提出按照沙层厚度或分区域建立时深曲线库的静校正方法, 取得了明显效果, 但该方法应用的效果依赖于分区策略, 对各分区之间如何实现静校正量的闭合则需要进一步研究；秦亚玲等^[10-11]先后在2007年和2009年提出了利用多项式拟合速度和深度的关系进而计算静校正量的方法, 该方法是常规时深曲线静校正方法的另一种形式, 并没有考虑沙层物性空间上的差异。

虽然沙漠地区有些沙层物性空间上具有较大差异, 但是其纵向上仍然保持着连续介质的特性, 因此时深曲线静校正应用的前提条件依然满足。为了更好地应用时深曲线方法, 需要根据沙层物性空间差异对方法进行改进。受肖泽阳等^[9]的方法启发, 本文以微测井测量数据的地质统计学分析结果为依据, 提出了一种逐点时深曲线静校正方法。通过克里金插值获得每个炮点和检波点处不同厚度的地震波垂直传播时间, 利用最小二乘法拟合得到每个炮点和检波点处的时深曲线量板, 进而计算静校正量。该方法通过沙层物性空间上的相关性, 采用地质统计学分析和处理手段, 实现了真正意义上的逐点时深曲线静校正量计算。

借鉴CHAMBERS等^[12-13]应用地质统计学的经验, 本文首先对微测井数据进行地质统计学分析, 明确了样本数据的统计特征。

S1X工区大部分地表为沙漠, 北部为盐碱滩, 南部为梭梭林浮土区, 古尔图河从工区中东部穿过。工区西部主要为高大的垄状沙丘, 东北部为高大的蜂窝状沙丘, 东南部为相对较小的沙丘。

搜集S1X工区64口微测井(图 1)数据并拾取初至, 选取10m、30m两个激发深度, 分析同一深度地震波垂向传播时间的统计特征及空间分布特征, 以确定合适的空间数据插值方法。

图 1(Fig. 1)

图 1 S1X工区卫星图片及微测井分布(红色三角形为微测井位置)

克里金插值法以变异函数理论和结构分析为基础, 利用对待插值点有影响的距离范围内的采样点来估计待插值点的属性值, 区域化变量存在空间相关性是克里金插值方法的应用条件^[14-17]。在应用克里金插值方法之前, 需要对空间变量的变异函数进行分析, 根据变异函数分析区域化变量是否存在空间相关性。如果存在空间相关性, 则可以利用克里金方法进行插值。肖斌等^[18]和李黎等^[19]对克里金方法进行了较为系统的研究, 指出克里金插值方法要求样本数据符合正态分布, 不符合正态分布的数据应转化成为正态分布的形式。绘制正态分位数图(quantile-quantile plot)是评估待插值数据是否为正态分布的一种直观方式, 如果数据散点接近一条直线, 则表明它服从正态分布。

图 2分别为深度10m、30m处的垂直传播时间正态分位数图, 两个深度处的地震波垂直传播时间分位数散点均接近直线, 基本符合正态分布特征, 因此可用于后续的空间插值。

图 2(Fig. 2)

图 2 不同深度垂直传播时间正态分位数 a深度10m; b深度30m

在使用克里金方法进行地质统计学插值之前, 必须对样本数据进行变异分析。实际计算时, 为了满足随机变量内在平稳的假设, 将样本数据配对分组为条柱单元, 取半变异函数值的平均值进行分析。变异函数定义为:

式中:h为成对位置相隔距离; n为距离h的样点对个数; z为与x相关的变量。

根据半变异函数图, 可检验随机变量是否具有空间相关性, 确定采用哪种变异函数模型进行插值计算。

为便于分析, 选择各微测井20m深度测量的地震波垂直传播时间绘制全方向半变异函数图, 并利用正态分布模型进行拟合(图 3)。结果表明:垂直传播时间变程值在15000m之内变化小, 相关性大。利用正态分布模型得到的插值结果精度更高, 可将其作为本工区变异函数理论模型。正态分布模型表达式为:

图 3(Fig. 3)

图 3 全方向半变异函数与理论模型

式中, γ_m(h)为根据正态模型拟合的变异函数; h为距离; c₀、c、a为待求系数, 本例拟合结果分别为0.585ms、9.628ms、13666.458m。

由上述分析可知, 研究区近地表同一深度地震波垂直传播时间的半变异函数存在且平稳, 因此本文将根据微测井测量值, 应用克里金插值方法完成所有炮点和检波点不同深度垂直传播时间的插值计算。

利用克里金插值方法计算出每个炮点和检波点处不同深度的地震波初至时间, 再对各个炮点和检波点处的深度-初至散点进行最小二乘法多项式拟合, 可以得到炮点和检波点的逐点时深曲线。

沙丘底(高速顶界面)足够稳定是可以采用时深曲线进行静校正的前提。当这一条件满足时, 可以通过小折射和微测井测量得到的高速顶高程进行插值, 得到整个工区沙丘底。静校正量的计算采用两步法, 首先将炮点和检波点校正到沙丘底, 然后使用统一的替换速度将炮点和检波点从沙丘底校正到最终基准面(或中间基准面)上。在炮点静校正量的计算中, 应考虑实际炮点对应的井深。

沙漠区受地表条件影响, 在沙丘厚的地方一般用浅井多井组合激发, 此时炮点位于沙丘底界之上; 在沙丘薄的地方采用单井深井激发, 这种情况下炮点位于沙丘底界之下。因此, 炮点和检波点静校正量的计算公式可以分为两种情况。

当炮点或检波点在沙丘底界之上时, 为:

当炮点或检波点位于沙丘底界之下时, 为:

式中:f(h)为根据时深曲线量板计算的值, 单位为s; h₁为沙丘底界离地表面的深度, h₂为炮点或检波点的实际埋深, 单位为m; d为基准面高程, b为炮点或检波点处沙丘底界高程, 单位为m; v_R为沙丘底界到基准面之间的替换速度, 单位为m/s; e为炮点高程, 单位为m。

将S1X工区64口微测井测量得到的深度-时间散点叠合显示, 结果如图 4所示。可见沙丘深度-时间关系差异较大, 相同深度最大时差在10ms以上, 因此常规时深曲线静校正方法明显不适用。

图 4(Fig. 4)

图 4 S1X工区全部微测井测量的深度-时间散点

由于插值得到的每个炮点和检波点处的深度-时间散点是否准确决定了后续逐点时深曲线量板及静校正计算结果是否可靠, 因此本文采用交叉验证的方法, 对10m、30m两个深度点的克里金插值结果进行了交叉验证。方法如下:首先移除一个样本数据, 利用剩余样本数据进行克里金插值得到移除样本点的预测值, 然后将移除样本点的真实值与预测值比较; 对所有样本点进行同样处理, 将所有样本点预测值与真实值交会显示。

图 5为根据样本数据交叉验证的结果, 克里金插值得到的预测值与真实值高度逼近, 整体误差为[-2ms, 2ms], 满足静校正计算要求。即使存在剩余静校正量, 也能通过后续的反射波剩余静校正予以消除, 从而验证了方法的有效性。

图 5(Fig. 5)

图 5 不同深度地震波垂直传播时间克里金插值结果交叉验证误差分析 a深度10m; b深度30m

为了验证纵向上深度-时间散点插值结果的可靠性, 选取工区东部和西部两个样本点位置, 通过插值得到其深度-时间散点, 并将其与实测深度-时间散点进行对比, 结果如图 6所示。可以看出, 预测得到的深度-时间散点与实测值高度逼近, 只存在[-1ms, 2ms]的随机误差。由于后续求取逐点时深曲线时要进行多项式拟合, 因此深度-时间散点的随机误差并不会影响静校正量的计算。

图 6(Fig. 6)

图 6 工区不同位置微测井深度-时间散点插值结果与测量值的比较 a工区西部; b工区东部

S1X工区内沙丘厚度从几米到70m左右均有, 其中30~40m厚度的大沙丘最为常见。为了验证逐点时深曲线静校正方法的应用效果, 选择工区内具有代表性的大沙丘区, 对不同方法静校正处理后的共炮点道集和叠加剖面进行了对比。图 7a展示了沙丘厚度在40m左右的共炮点道集, 图 7b和图 7c分别为常规时深曲线方法和逐点时深曲线方法静校正后的共炮点道集。可以看出, 常规时深曲线静校正方法虽然一定程度上校正了沙丘起伏造成的时差, 但是其初至和反射波同相轴上仍然存在明显的静校正问题, 这表明在大沙丘区应用常规时深曲线静校正方法效果有限(图 7b)。与常规时深曲线静校正结果相比, 使用逐点时深曲线静校正方法处理的结果初至平滑, 沙丘起伏引起的同相轴扭曲得到了有效恢复, 中短波长静校正问题得到了较好解决。图 8对比了两种方法计算的检波点静校正量曲线, 可以看出, 两种方法计算的静校正量总体趋势相同, 但受地表沙丘变化的影响, 局部细节上差别较大。

图 7(Fig. 7)

图 7 大沙丘区共炮点道集静校正效果对比 a原始记录; b常规时深曲线静校正; c逐点时深曲线静校正

图 8(Fig. 8)

图 8 大沙丘区不同静校正方法计算的检波点静校正量

采用同样的叠加速度对两种方法静校正后的道集叠加成像(图 9), 可以看出, 常规时深曲线静校正(图 9a)只是消除了一部分地表起伏引起的时差, 在沙丘较大处同相轴连续性较差, 由沙丘引起的时差仍然没有得到很好的校正。而用本文方法静校正处理后(图 9b), 同相轴连续性明显变好, 不存在明显的中长波长静校正问题, 虽然短波长静校正问题仍有残留, 但通过后续的剩余静校正基本可以将其消除。采用折射静校正方法进行了对比处理, 结果如图 9c所示。可以看出, 大沙丘区使用折射静校正方法处理的效果仍不理想, 甚至无法达到常规时深曲线静校正方法处理的效果。

图 9(Fig. 9)

图 9 大沙丘区叠加剖面静校正效果对比 a常规时深曲线静校正; b逐点时深曲线静校正; c折射静校正

本文将地质统计学方法应用于沙漠区静校正量的计算, 实现了一种逐点时深曲线静校正方法。通过计算变异函数和对微测井样本进行克里金插值, 获得时深曲线量板, 对实际资料进行逐点时深曲线静校正值计算, 克服了常规时深曲线静校正技术在沙丘物性差异较大地区计算误差大的问题, 进一步提高了时深曲线静校正方法的适用性, 为解决我国西部沙漠、黄土塬等地区的静校正问题提供了一条新的途径。S1X工区实际资料应用表明, 改进后的逐点时深曲线静校正方法比常规的时深曲线静校正方法和折射静校正方法的应用效果有明显改善。