克里金估值法中各向异性土体参数的空间变异函数套合研究

适用于普通克里金法（ordinary Kriging，OK）的随机变量需满足以下2个假设：1）随机变量的均值存在并且与空间距离h无关；2）随机变量的方差存在并且与点x位置无关。当已知观测数据满足上述假设条件时，空间距离h对应的变异函数值γ（h）可以表示为^［4］

式中：Z（x_i）为点x_i位置处的观测值；N_h表示空间距离为h的观测数据点对的个数。需要注意的是，对于各向异性土体而言，变异函数值不仅依赖于2点之间的空间距离，还依赖于2点之间的方向。

变异函数中一个非常重要的参数就是变程，通过某方向上变异函数值与空间距离的对应关系可以确定该方向的变程。如图1所示，空间距离h越大，变异函数值越大，表明空间距离越近的点具有更相近的特征。当空间距离h达到某个值时，变异函数值将趋于稳定，表明此后的点对不再具有相关性，该空间距离临界值即为当前方向的变程。

图1  变异函数值随空间距离的变化

Fig.1  Variation of variogram value with spatial distance

采用澳大利亚国家软土现场试验中心在新南威尔士州巴利纳附近约100 m×100 m的区域内通过孔压静力触探试验（piezocone penetration test， CPTU）所获取的土体锥尖阻力值^［17］，计算竖直剖面不同方向上锥尖阻力的变程。该试验一共采集了26个钻孔的数据，每个钻孔竖向测量点的间距为0.01 m。综合考虑钻孔间距、土层分布和剖面选取等因素，本文在原有26个钻孔中选取21个钻孔的观测数据进行后续计算及分析，如图2所示。

图2  CPTU钻孔位置

Fig.2  Location of CPTU boreholes

将图2中剖面1作为典型竖直剖面来计算其锥尖阻力值在不同方向上的变程，剖面1包含了编号为NS9、NS14、NS20、NS26和NS12等5个钻孔，通过CPTU试验所获取的钻孔土体锥尖阻力值如图3所示。

图3  钻孔锥尖阻力观测值

Fig.3  Observed cone-tip resistance of boreholes

竖直剖面1包含了有机土层和其他土层等多个土层，不同土层的锥尖阻力值存在明显差异。为了方便分析，本文统一选取有机土层展开研究。如图4所示，能够利用包含有机土层的5个钻孔的观测数据点对进行竖直剖面变程计算的方向为-14°~9°。其中，沿着水平方向逆时针旋转形成的角度为正，顺时针旋转形成的角度为负。

图4  竖直剖面1变程计算方向

Fig.4  Calculation direction of range for vertical profile 1

为了对比竖直剖面多个方向的变程，在-14°~9°范围中依次选取-14°、-9°、-3°、0°、3°和9°等6个方向。为了充分利用观测数据并减少不必要的运算，将竖直方向上的锥尖阻力值每隔0.2 m取一个局部平均值作为观测数据。现以3°方向为例，计算有机土层锥尖阻力在该方向的变程，具体步骤如下：

1）应用最小二乘法对3°方向上的锥尖阻力观测值进行拟合，获取该方向上锥尖阻力的线趋势。

2）运用下式计算3°方向上的锥尖阻力的残差值：

式中：R（x， y）表示点（x， y）处锥尖阻力的残差值；

Z（x， y）表示点（x， y）处锥尖阻力的观测值；T（x， y）表示点（x， y）处锥尖阻力的趋势值。

3）基于公式（1），利用锥尖阻力残差值计算不同空间距离对应的变异函数值。将空间距离h作为横轴，将变异函数值γ（h）作为纵轴，绘制散点图。

4）利用经典的变异函数模型诸如球型模型、指数模型、立方体模型和高斯模型等分别对上述散点进行拟合，选择拟合效果最佳的模型作为实际变异函数模型，最终确定有机土层锥尖阻力在3°方向的变程，如图5所示。

图5  竖直剖面1在3°方向的变异函数

Fig.5  Variogram in 3° direction for vertical profile 1

从图5可以看出，变异函数值首先随空间距离的增大而增大，然后趋于稳定，这与变异函数的理想变化情况一致。高斯模型能够较好地对散点进行拟合（拟合优度R²=0.93），继而得出有机土层锥尖阻力在3°方向的变程为22.8 m。

利用前述步骤，计算有机土层锥尖阻力在-14°、-9°、-3°、0°和9°等其余5个方向的变程。此外，分别计算5个钻孔中的有机土层锥尖阻力在竖直方向上的变程，取5个结果的平均值作为竖直剖面1在90°方向上的变程。已计算的竖直剖面1有机土层锥尖阻力在多个方向上的变程列于表1。由表1可以看出，本文场地土竖直剖面在除90°方向之外的其他方向的变程差异不大，然而，0°方向变程达到90°方向变程的15倍，显示出二者存在显著差异。

将表1所示竖直剖面1有机土层锥尖阻力在7个方向上的变程绘制于直角坐标系中，变程散点与原点的连线长度即为变程，连线与横轴的夹角即为变程方向，横轴数值表示变程在水平方向上的分量，纵轴数值表示变程在竖直方向上的分量。利用椭圆函数对变程散点进行拟合，绘制于直角坐标系中的变程散点以及椭圆拟合结果如图6所示。

图6  竖直剖面1变程散点及椭圆拟合

Fig.6  Range scattered points and ellipse fitting for vertical profile 1

从图6可以看出，竖直剖面1在0°方向附近的变程变化幅度较小，同时与90°方向的变程差异较大。椭圆函数难以描述这种变程变化特征，因此对上述场地竖直剖面1不同方向变程散点的拟合效果不佳。

卡西尼卵形线定义为平面上到2个定点距离之积为常数的所有点组成的图形。2个定点即为焦点，坐标分别为（-c，0）和（c，0）（c>0），常数为a²（a>0）。卡西尼卵形线的数学表达式如下：

假设卡西尼卵形线与直角坐标系的横轴交点即水平方向变程为r_x，与纵轴交点即竖直方向变程为r_y，根据公式（3）可以推导出：

已知r_x>r_y>0，卡西尼卵形线的形状由参数a和c的相对关系决定。当c<a<2–√c时，卡西尼卵形线是1条光滑的曲线，曲线中部有凹进的“细腰”；当a=2–√c时，与前种情况一样，但曲线中部变平；当a>2–√c时，与前种情况一样，但曲线中部凸起。

当r_x≫r_y，即c<a<2–√c时，利用卡西尼卵形线函数对已知变程散点进行拟合，其结果如图7所示。

图7  竖直剖面1变程散点及卡西尼卵形线拟合

Fig.7  Range scattered points and Cassini oval line fitting for vertical profile 1

从图7可以看出，使用卡西尼卵形线函数可以较好地对上述场地竖直剖面1不同方向变程散点进行拟合（R²=0.96）。图中卡西尼卵形线在左右两端呈现出光滑的弧线，在中部则呈现出凹进的“细腰”，比较符合竖直剖面1不同方向变程所表现出来的变化特征。

为了检验本文所涉场地其他竖直剖面锥尖阻力不同方向的变程是否同样适用卡西尼卵形线函数，在场地上选取6个新的竖直剖面，如图8所示。

图8  场地竖直剖面划分

Fig.8  Vertical profiles division of the site

参照竖直剖面1不同方向变程的计算方法，对图8所示的竖直剖面2~7在不同方向上的变程进行计算，并利用卡西尼卵形线函数进行拟合，结果如图9所示。由图9可以看出，使用卡西尼卵形线函数可以较好地对场地多个竖直剖面上的变程散点进行拟合，再次验证了卡西尼卵形线函数用于拟合各向异性土体竖直剖面参数变程的可行性。

图9  不同竖直剖面变程散点及卡西尼卵形线拟合

Fig.9  Range scatter and Cassini oval line fitting in different vertical profiles

考虑到卡西尼卵形线相较于椭圆更加适合描述本文所涉场地竖直剖面不同方向的变程，有必要进一步提出等效应卡西尼卵形线套合模型。

与等效应椭圆套合模型类似^［9］，本文所提出的等效应卡西尼卵形线套合模型满足以下3个假设条件：

1）竖直剖面上有2个方向的变异函数值是已知的，而且这2个方向的变异函数类型一致。

2）竖直剖面不同方向上的变程满足卡西尼卵形线函数。

3）竖直剖面不同方向上的变异函数值满足卡西尼卵形线函数。

在上述条件成立的前提下，提出等效应卡西尼卵形线套合模型，用于计算竖直剖面待求方向上的变异函数值，具体步骤如下：

1）将竖直剖面上的待求方向记为L，将方向L与水平方向的夹角记为θ，将方向L上的空间距离记为h_L，将方向L上的变程记为r_L。

2）分别计算竖直剖面水平方向和竖直方向上的变异函数值，利用同一类型的拟合函数求出该竖直剖面水平方向变程r_x和竖直方向变程r_y。

3）根据假设可知，依据式（4），基于r_x和r_y可以绘制一条等效应变程卡西尼卵形线。由式（3）可知，待求方向L上的变程r_L为

4）计算待求方向L上的空间距离h_L在水平方向和竖直方向上的等效应空间距离：

5）将上述水平方向和竖直方向上的等效应空间距离代入各自的变异函数，求得对应的变异函数值γ_x（h_x）和γ_y（h_y）。

6）根据假设可知，参照式（4），基于γ_x（h_x）和γ_y（h_y）可以绘制以a_γ和c_γ为参数的等效应变异函数卡西尼卵形线：

7）由式（3）可知，待求方向L上的变异函数值γ_L（h_L）为

克里金法作为一种最优线性无偏估计方法，它要求预测值的平均误差接近0，同时误差的方差最小，据此可推出如下线性方程组^［4］：

式中：N为观测点的个数；γ_ij（i， j = 1， 2， …， N）为观测点x_i与观测点x_j之间的变异函数值；γ_ik为观测点x_i与预测点x_k之间的变异函数值；λ_i为分配给观测点x_i处数据的权重；μ为拉格朗日乘数。

若土体竖直剖面的变异函数套合模型已知，即竖直剖面上任意2点之间的变异函数值已知，就可以利用式（9）反推权重矩阵λ，利用普通克里金法对预测点x_k处的参数进行估算^［4］：

式中：Z（x_k）为预测点x_k处的估计值。此外，可以得到普通克里金法估算的方差σOK，以衡量普通克里金法的插值精度^［4］：

以竖直剖面1为例：将编号为NS9、NS20、NS26和NS12的钻孔作为已测位置，其钻孔数据作为观测值；将编号NS14的钻孔作为待估位置，其钻孔数据作为待估值。在已知位置的9.5 ~13.5 m深度范围内（全部位于有机土层）每隔0.1 m取1个点作为观测点，在待估位置取类似点作为预测点，如图10所示。

图10  竖直剖面1观测点与预测点

Fig.10  Observation and prediction points in vertical profile 1

基于普通克里金法，分别利用等效应卡西尼卵形线套合模型和等效应椭圆套合模型，根据竖直剖面1观测点的锥尖阻力值估算预测点的锥尖阻力值。不同深度处的观测值及基于等效应卡西尼卵形线套合模型和等效应椭圆套合模型的预测值如图11所示，统计参数如表2所示。

图11  基于2种套合模型的预测值对比

Fig.11  Comparison of prediction values based on two nested overlap models

从图11可以看出，相比基于等效应椭圆套合模型的预测结果，基于等效应卡西尼卵形线套合模型的预测结果更接近观测值。基于等效应椭圆套合模型的预测结果的误差在20%以上的点数占总预测点数的13%，而基于等效应卡西尼卵形线套合模型的预测结果的最大误差为9%，表明等效应卡西尼卵形线套合模型的应用效果明显优于等效应椭圆套合模型。从表2也可以进一步看出，基于等效应椭圆套合模型的预测值的平均值、极差和标准差都更加接近观测值所对应的统计参数。尽管基于等效应椭圆套合模型预测值的平均值也接近观测值的平均值，但是较大的极差与标准差表明，基于等效应椭圆套合模型的预测值的离散性较大。这也使得等效应椭圆套合模型预测值的误差更大，可以从更大的OK标准差的均值得到印证。基于等效应椭圆套合模型的预测值之所以产生这么大的误差，主要源于其对于竖直剖面不同方向变程的拟合效果欠佳，变程效果拟合的好坏显著影响了最终的预测结果。

图8所示的不同竖直剖面在其0°方向的变程可以组合成本文所涉场地土水平剖面在不同方向的变程。由于剖面6与剖面7的朝向是平行的，取两者在其0°方向变程的平均值作为水平剖面在该朝向的变程。将其0°方向具有最大变程的竖直剖面朝向作为水平剖面的0°方向，其他竖直剖面朝向的方位角以逆时针旋转为正而顺时针旋转为负的原则确定，最终确定的水平剖面在不同方向的变程如表3所示。

已知r_y<r_x<3–√r_y，即a>2–√c，此时的卡西尼卵形线是一条中部凸起的光滑曲线，其形态与椭圆相似。分别利用卡西尼卵形线函数和椭圆函数对水平剖面不同方向的变程进行拟合，结果如图12所示。可以看出，椭圆函数与卡西尼卵形线函数都可以很好地对本文所涉场地水平剖面锥尖阻力不同方向的变程散点进行拟合。

图12  水平剖面变程散点及两种函数拟合

Fig.12  Range scatter and ellipse and Cassini oval line fitting for horizontal profile

为了检验等效应卡西尼卵形线套合模型在水平剖面参数预测方面的适用性，利用普通克里金法对本文所涉各向异性场地土水平剖面锥尖阻力进行预测，并与等效应椭圆套合模型的预测结果进行对比。将深度为11 m处的有机土层水平剖面作为研究对象，在该水平剖面上取14个钻孔作为观测点，而另外7个钻孔作为预测点，如图13所示。

图13  水平剖面观测点与预测点（深度为11 m）

Fig.13  Observation and prediction position for horizontal profile (11 m depth)

基于普通克里金法，分别利用等效应卡西尼卵形线套合模型和等效应椭圆套合模型，根据水平剖面观测点的锥尖阻力值估算预测点的锥尖阻力值。不同位置处的观测值与基于等效应卡西尼卵形线套合模型和等效应椭圆套合模型的预测值的对比如图14所示，统计参数如表4所示。

图14  水平剖面观测值与预测值对比

Fig.14  Comparison of observation and prediction values for horizontal profile

从图14可以看出，基于等效应卡西尼卵形线套合模型的预测结果与等效应椭圆套合模型的预测结果接近，两者与观测值的误差大都小于12%。编号为NS28的钻孔预测值与该钻孔观测值的误差为20%，主要是因为该钻孔与其他观测点的位置较远，与观测点之间的相关性较弱，从而造成该预测点的预测误差较大。从图14可看出，基于等效应卡西尼卵形线套合模型的预测结果更接近于观测值。从表4也可以进一步看出，相较于等效应椭圆套合模型预测值，等效应卡西尼卵形线套合模型预测值的平均值、极差和标准差都更加接近观测值所对应的统计参数。此外，对应于等效应卡西尼卵形线套合模型的OK标准差的均值也小于等效应椭圆套合模型。综上可见，等效应卡西尼卵形线套合模型同样适用于水平剖面参数的预测，预测性能稍好于等效应椭圆套合模型。