半变异函数及取样间距对克里金法在海洋地层分析中的影响研究*

杨 阳 李 飒 何福耀 施展玲

(①天津大学建筑工程学院 天津 300350)(②中海油上海分公司 上海 200032)

0 引 言

在海洋工程建设前，对施工场地进行工程地质勘察，了解海底地形、地貌以及海底地层分布特征尤为重要。工程钻探是一种常用的勘察技术手段，取样获得的钻孔数据可以很好地反映地层垂向分布特征(Lemon et al., 2003)。陆地施工方便，获得钻孔数据较容易且成本低。而对于海洋工程，由于勘察条件的限制及海底地理环境的复杂，获取钻孔数据极其困难，且成本较高。相比陆地，海上钻孔取样密度大幅度降低，仅依靠少量的钻孔数据显然不能完整地获取地层厚度分布特征(李一保等， 2007； 李安龙等， 2016)。浅剖数据在测线方向上具有良好的连续性，若将极有限的钻孔数据和连续的浅剖数据有机结合起来，将有利于提高地层厚度判别精度(周晓， 2014)。

地层厚度在空间分布上具有一定的随机性和结构性。因此，可将地层厚度视为区域化变量，采用克里金法进行空间插值估计地层厚度(Cressie， 1991； Chilès et al.，1999； 周小文等， 2005； 许国等， 2013； 刘晓磊等， 2018)。在利用克里金法估计地层厚度方面，国内外学者展开了广泛的研究。李晓军等(2014)采用IK-OK方法较准确地估计了崇明越江隧道工程段地层厚度； 杜文凤等(2010)将地层厚度和地震振幅作为源数据，利用协同克里金法预测煤层厚度； Hengl et al. (2004)则采样回归克里金法预测西伯利亚某土层厚度； Jang et al. (2016)采用指示克里金多阈值组合的方法在估计冲积扇地区低渗性土材料厚度中取得良好的效果。但大多数学者在利用克里金法估计地层厚度方面，研究区域多局限在陆地，针对海底地层厚度估计的研究较少。正确且准确地识别海底地层厚度分布对海洋工程非常重要。因此，将克里金法应用在海底地层厚度识别中具有重要的意义。

对地层厚度进行克里金插值前，首先需要对已有数据进行变异结构分析，通过变异结构分析构建一个半变异函数模型； 再利用该模型计算未知区域地层厚度，从而得到地层厚度分布特征(张征等， 1995； 崔洪庆等， 2003)。目前国内外针对地层厚度插值的研究多在分析比较不同克里金法插值的优劣性(Penizek et al.，2006； Kuriakose et al.，2009； 易湘生等， 2012； 张琳娜等， 2016)，如：Penizek et al. (2006)采用了不同的克里金法对土壤厚度进行插值，发现协同克里金优于普通克里金、线性克里金和回归克里金； Kuriakose et al. (2009)则认为回归克里金优于普通克里金。对半变异函数模型选取的研究较少，而半变异函数模型选取适当与否将直接影响插值精度。因此，有必要对最优半变异函数模型作进一步探讨(Kambhammettu et al.，2011)。同时，对于半变异函数模型最优性的判别，通常采取决定系数越大越佳，残差越小越佳的原则确定(刘爱利等， 2012)。但是，很多时候会出现各模型的决定系数和残差的差异很小而无法择优(施小清等， 2009； 张涛等， 2016)，甚至出现决定系数最大的模型，其插值结果却不是误差均值最小的情况(巫振富等， 2012)。若仅依据上述原则选择模型可能无法保证所选模型为最优，造成预测结果的偏差，所以对于最优的模型需进一步验证其可靠性。

此外，插值精度不仅受半变异函数模型的影响，也受样本数量、取样间距、原始数据等的影响(朱会义等， 2004)。通过目前的研究资料发现，对插值精度的研究主要涉及取样空间尺度(Robeson， 1993)、样本容量(Journel et al., 1978； Bourennane et al.， 2000)等的影响。因此，了解取样间距对插值精度的影响尤为重要。

为了研究适合海洋地层厚度的最佳半变异函数模型，本文基于GS+8.0软件，以某井场3个反射层组的浅剖数据为源数据，钻孔数据为实测验证数据，采用克里金法，利用不同的半变异函数模型分别对地层厚度变异函数进行拟合对比分析。一方面探讨更适合该井场地层厚度插值的模型，同时研究取样间距对预测误差的影响，为该井场地层厚度估计提供依据； 另一方面给出海底地层厚度估计中半变异函数模型选择的一般思路，为其他海域地层厚度插值提供参考。

1 地统计学方法

克里金插值法属于地质统计法的一种，它是以无偏为约束，寻求估计方差最小的一种估计方法(孙洪泉， 1990； 王鹏等， 2013)。根据地质统计理论，半变异函数既能考虑空间统计特征，又能同时描述区域化变量在空间分布的随机性和结构性，完整反映区域化变量的空间变异规律。由于地层厚度符合区域化变量的特征，因此可以用地层厚度取样值建立半变异理论函数，利用半变异理论函数求出研究区域内任意一点处的地层厚度。

1.1 变异结构分析

半变异函数是计算研究区域内地层厚度变化的关键。关于地层厚度的理论半变异函数定义为地层厚度Z(x)在点x和x+h处数值之差的方差之半，其表达式为：

(1)

在实际计算γ(x，h)时，采用地层厚度实验半变异函数计算。测线方向上相距h的xi和xi+h的取样点地层厚度值为Z(xi)和Z(xi+h)，测线方向上相距h的取样点的数据对数目为N。对于任意的h，实验半变异函数γ*(x，h)可以表示为：

(2)

假定任意一点x0处的地层厚度待估值Z*(x0)是周围实测值Z(xi)的线性组合，即：

(3)

式中，λi为权重系数；Z*(x0)为点x0处的地层待估值；Z(xi)为已知点i=1，2，3，…，n处的地层厚度实测值。根据无偏约束条件可推导出：

(4)

根据最小方差条件，可得：

(5)

式中，γ(xi，xj)为地层厚度采样点xi和xj的理论变异函数。

(6)

1.2 半变异函数理论模型

半变异函数是半变异值γ(x，h)与距离之间的函数。其中半变异值表示空间相关性，等于所有取样点中距离为h的厚度值之差的期望。实验半变异函数通常由3个参数表征：(1)块金值c0，主要是由于空间变异和测量误差导致； (2)基台值c0+c，一般来说表征数据对间的最大方差； (3)变程a，代表土性参数存在相关性的最远距离。实验变异函数是根据采集的数据获得的非连续分布函数，需采用理论变异函数对实验变异函数进行拟合以获得连续的变异性评价，岩土工程中常用的半变异函数的理论模型见表 1。

表 1 岩土工程常用变异函数模型Table 1 Typical theoretical semi-variogram models in geotechnical engineering

2 地层的厚度

地层厚度是识别地层结构和构造的重要依据之一。对于海洋地质勘查，通常通过浅地层剖面仪所观测到的海底地层声波接收反射特征得到地层分布特征，反映地层厚度变化规律。

浅地层剖面仪是一种基于水声学的走航式水下浅部地层结构和构造的地球物理探测方法，其工作原理是利用声波在不同介质中传播速度的差异性。浅地层剖面仪的发射基阵(震源)按一定时间间隔向下发射声波，声波穿透海水抵达海底后，一部分声波在分界面处反射并返回，另一部分则被接受基阵(水听器)接收。

产生反射的条件是分界面两侧存在波阻抗差(反射系数R)。在浅地层剖面调查中，通常认为声波是垂直入射的，此时反射系数R有：

(7)

图 1 井场过井剖面Fig. 1 Shallow profile section of over well

式中，ρ1V1、ρ2V2分别是上下两层介质的声阻力率，ρ和V分别表示介质的密度和声速。当相邻地层介质存在较大的声阻力差，在界面就会有较强的声强，接收到较强的反射信号，在浅地层剖面仪终端显示器上反映灰度较强的分界线。通过观测分析海底地层的声波反射时间可计算出地层厚度。

3 地层数据分析

为尽可能满足采样点数据对的要求，选定采样点的距离间隔为50im，每组地层厚度数据点数目为78个，符合数据配对数目应大于30的原则(Robeson， 1993)。

3.1 半变异函数拟合

GS+是目前专业的常用地统计学软件之一，在对区域化变量的变异结构分析方面具有较明显的优势。因此，本文利用GS+8.0软件，选取高斯函数、球状函数和指数函数进行地层厚度变异函数拟合。有效滞后距为1950im。拟合结果如图 2所示，模型具体参数见表 2。

图 2 拟合半变异函数图Fig. 2 Fitted semivariograma. 第Ⅰ层； b. 第Ⅱ层； c. 第Ⅲ层

表中，r2为决定系数，RSS为残差。块金值c0与基台值c0+c的比值称为空间相关度，可以反映影响因素中结构性因素和随机性因素的作用。根据Cambardella et al.(1994)的划分原则，空间相关度小于25%，表明变量空间相关性强； 大于75%，表明空间相关性较弱； 处于两者中间，说明空间相关性中等。从3种变异函数模型参数可以看出，3组地层厚度在3种模型下的块金值c0与基台值c0+c的比值均小于10%，说明研究区域地层厚度分布以结构性变异为主，具有强烈的自相关性。

表 2 变异函数参数统计Table 2 Variance function parameter statistics

海洋地层在沉积过程中大致受岩土体形成及后期改造过程中各种宏观规律控制，使得不同点的地层厚度之间存在极强的空间结构性。当c0和c0+c发生变化时，变程a的变化幅度较大，其中变程大小依次为指数模型、球状模型和高斯模型，但均未超过研究区域范围。通过决定系数越大，残差越小的原则确定模型最优性(刘爱利等， 2012)。对比模型拟合参数，同组地层各模型的决定系数较为接近，高斯模型略高于球状模型和指数模型。

表 3 交叉验证结果分析Table 3 Analysis of cross validation results

对比残差，指数模型最大，球状模型和高斯模型次之。仅从拟合效果角度分析，3种函数模型中最优的是高斯模型，其次是球状模型。

3.2 交叉验证

在对3组地层厚度进行克里金插值前，应对实测采样数据模型进行整体的精度评价。交叉验证是一种常用的精度验证方法，即每次移除一个或多个实测采样数据，然后使用其他位置的数据来预测与其相关的数据。重复操作，最后将每个数据点的真实值和估计值进行统计回归分析，进而评价模型的精度。表 3为交叉验证误差统计结果。

图 3 第Ⅰ层插值结果Fig. 3 Ordinary Kriging of layer Ⅰa. 指数模型； b. 球状模型； c. 高斯模型

图 4 第Ⅱ层插值结果Fig. 4 Ordinary Kriging of layer Ⅱa. 指数模型； b. 球状模型； c. 高斯模型

表3中，SE是回归系数的标准误差，其值越小说明精度越高；r2是相关系数的平方，其值越大说明回归效果越好；SE预测值是标准误差预测值，其值越大说明误差越大。通过这3个指标可以综合评价不同半变异理论模型的预测精度。从表3中可以看到，在第Ⅰ、Ⅱ层中，3种模型的r2接近，较难看出3种模型回归效果的优劣，但从SE预测值中可以看出球状模型的误差明显小于其余两种模型。在模型数据交叉验证环节，球状模型表现出的整体精度最高，高斯模型次之，指数函数的误差最大。可见，结构分析最佳的模型不一定是误差最小的模型。在出现结构分析和误差分析的择优结果不一致的情况下，很难说明球状模型和高斯模型哪一个是绝对最优。因此，需对插值效果进一步对比分析以选择最适合的模型。

图 5 第Ⅲ层插值结果Fig. 5 Ordinary Kriging of layer Ⅲa. 指数模型； b. 球状模型； c. 高斯模型

3.3 地层厚度插值效果分析

在对数据完成结构分析的基础上，经交叉验证的方法对不同理论模型精度评定后，采用普通克里金法对3组地层厚度(沿测线方向)进行插值。插值结果(图 3～图 5)。对比3组地层厚度插值云图，可以看到，球状模型和高斯模型在第Ⅰ反射层组0、690im、1300im和1600im处，第Ⅱ反射层组910im、1560im、1950im和2340im处，第Ⅲ反射层组300im、520im、1410im和3380im处表现出较高的相似性，与结构分析和交叉验证分析结果相吻合。

表 4 K1插值结果对比Table 4 The actual thickness of K1 and prediction

利用研究区域内的浅钻孔K1(非采样点)的实测厚度数据和利用普通克里金法插值得出的地层厚度进一步检验(表 4)。从图中可发现，不同理论模型插值的结果在浅钻实测点均存在不同程度的误差。其中，球状模型的最小，在3组地层厚度误差分别是1.9%、4.4%和2.3%。在地层厚度变化幅度较大的第Ⅰ、Ⅱ层中尤为明显。由以上分析可知，球状模型在该井场地层厚度离散数据插值中具有较高的实用性。

从K1插值结果数据对比可以看到，结构分析最佳的模型，即决定系数最大，残差最小的模型，不一定是误差最小的模型。因此，在进行半变异函数的选择时，对于其可靠性验证需要通过以下步骤：首先进行结构分析，其次通过交叉验证评定模型整体精度，最后进一步验证克里金插值效果。通过上述方法综合分析来选择最为适宜的半变异函数模型。仅根据结构分析，就有可能造成判断结果的偏差。图 6为最优半变异函数模型选择流程图。

图 6 最优半变异函数模型选择流程图Fig. 6 Flow diagram of selection of optimal Semi-variogram function

4 取样间距的影响

图 7 第Ⅲ层交叉验证Fig. 7 Cross validation of layer Ⅲa. 50 m； b. 75im； c. 100im； d. 125im； e. 150im

求解半变异函数值时，若取样间距增大，取样点空间分布离散稀疏，无法反应小范围内的空间变异性，即无法准确反映地层厚度变异规律(Cambardella et al.， 1994)。本文通过设置不同的取样间距来研究其对插值结果的影响。为满足数据对数目大于30的要求，故将取样间距最大值控制在150im以内，分别为50im、75im、100im、125im、150im。半变异理论模型选取上文中得到的最优模型球状模型。将不同取样间距获得的地层厚度数据进行交叉验证。图 7分别为第Ⅲ组地层不同取样间距时计算出的预测值和实际值的回归曲线，其中横坐标为预测值，纵坐标为实际值，虚线为理想回归曲线。从3组地层厚度数据交叉验证误差统计结果(图 8)可以看出，取样间距越大，SE值越大，拟合精度越低。第Ⅰ、Ⅱ层的相关系数r2受采样间距影响较小。第Ⅲ层的相关系数有一定程度的降低，当取样间距为150im时，相关系数为0.801。通常认为r2>0.8，即可认为具有较好的回归效果，因此该条件下地层厚度仍具有较强的自相关性。同时也说明对于球状模型，增大取样间距对地层厚度变化较大的地层回归效果影响较小，而对于地层厚度变化较平缓的地层回归效果影响较大。

图 8 交叉验证Fig. 8 Cross validationa. SE； b. 相关系数r2； c. SE预测值

SE预测值随着取样间距的增大而增大，说明预测误差逐渐增大。定义：

SE预测值变化率=

(8)

第Ⅰ层SE预测值变化率分别为77.38%、54.7%、45.77%和24.48%； 第Ⅱ层SE预测值变化率分别为119.11%、25.45%、16.25%和12.56%； 第Ⅲ层SE预测值变化率分别为25.32%、20.01%、7.33%和7.23%。同一取样间距下，地层厚度变化较剧烈的地层SE预测值变化率大于地层厚度变化不大的地层。对于地层厚度变化剧烈的地层，减小取样间距可以有效的减少插值误差，而对于地层厚度变化不大的地层，减小取样间距对插值精度提高的意义不大。不同取样间距的插值结果表明，取样间距在一定程度上会影响地层厚度的插值精度。通过交叉验证数据比较分析，可以看到，合理选择取样间距有利于提高预测结果的可靠性。取样间距的选取取决于使用者对精度的要求。对于本算例，当选择r2不小于0.8且SE预测值不大于1作为精度要求，对于土层厚度变化较剧烈的土层(第Ⅰ、Ⅱ层)，取样间距控制在100im以内； 对于土层变化不大的土层(第Ⅲ层)，本文的最大取样间距150im仍可以满足要求。

5 结 论

本文基于GS+8.0软件，以浅剖数据为源数据，钻孔实测数据为验证数据，采用克里金法对海洋地层厚度进行空间插值。通过比较不同变异函数模型和不同采样间距的结果，得到以下结论：

(1)海洋地层厚度在空间分布上具有一定结构性与局部的变异性。因此，可将浅剖数据得到的地层厚度作为区域化变量，采用普通克里金进行空间插值。插值结果表明，普通克里金是一种有效的海底地层厚度预测方法。

(2)利用不同变异函数理论模型对实际半变异函数拟合，在结构变异分析的基础上，进行交叉验证获得模型精度评价，并将钻孔K1实测数据作为验证。结果表明，决定系数最大，残差最小的模型不一定是误差最小的模型，应对不同模型下的插值结果进行综合分析来选择最合适的模型。在本文的井场地层厚度估计中，球状模型为最优模型，高斯模型次之。

(3)取样间距对插值精度有一定的影响。随着取样间距增大，插值拟合精度越低，预测误差越大。对于球状模型，增大取样间距对地层厚度变化剧烈的地层回归效果影响较小，对于地层厚度变化不大的地层回归效果影响较大。同时，减小取样间距可有效减小地层厚度变化剧烈的地层的插值误差，而对地层厚度变化不大的地层插值精度的提高意义不大。