五种插值法的对比研究论文

nearest：执行速度最快，输出结果为直角转折；linear：默认值，在样本点上斜率变化很大；spline：最花时间，但输出结果也最平滑；cubic：最占内存，输出结果与spline差不多。

地层面的拟合是多源地质建模中最为重要的步骤。无论是通过Delaunay细分方法增加节点还是通过网格分级增加的节点，都需要进一步求取其高程值。因此，必须借助插值方法来对有限的数据点信息所形成的初始地层面进行精确、光滑处理。插值是指在根据已知的数据计算未知值的过程，结果是形成一个连续分布的数据场(Spragueetal.，2005)。目前存在的插值方法众多，但适合三维数据场表达的插值方法还是比较有限的(胡小红等，2007)，常用的有距离加权反比法(Inverse-Distance Weighting，IDW)和普通Kriging法，其中距离加权反比法属于一种确定性差值方法，Lu et al.(2008)对该方法进行扩展，根据样本的数量和分布密度特征，使其能够根据样本的特征来确定其参数的取值；而Kriging法属于一种不确定性差值方法，Jessell(2001)对其进行了深入研究，基于该方法提出一种势场(potential-field)的插值方法，能够处理存在断层的不连续数据场。这里，仅对这两种方法进行讨论。

5.3.3.1 距离加权反比法

距离反比加权法是最常用的地质数据插值方法之一。它首先由气象学家及地质学工作者提出，后由D.Shepard进行改进，故该方法被称为Shepard方法。

距离反比加权法的基本思想是将插值函数f(P)定义为各数据点函数fk的加权平均，它认为与待插值点距离最近的若干个已知采样点对待插值点的值贡献最大，其贡献与距离的某次幂成反比。

距离反比加权法的基本原理可用下式表示：

数字地下空间与工程三维地质建模及应用研究

式中：f(P)是待插值点P的估计值；fi是第i(i=1，…，n)个已知采样点Pi的样本值；di是第i个样本点Pi与待插值点P的广义距离；v是距离的幂，它显著影响内插的结果，它的选择标准是最小平均绝对误差。相关研究结果表明，幂越高，内插结果越具有平滑的效果(Lu et al.，2008)。在幂指数为2时，不仅能得出较满意的内插结果，而且具有容易计算的优点。在实际应用中，一般用距离平方反比法来求待估算值。

传统的距离反比加权法是非常简单和自然的，但应用于实际的地质特征插值时，却存在以下明显的缺陷：

(1)当数据点的数目非常庞大时，f(P)的计算量将变得十分巨大，计算量的庞大甚至可能导致该方法变得无法实现。

(2)该方法只考虑了从Pi到P的距离，而没有考虑其方向。事实上，只考虑距离的大小是不充分的，有的已知离散点虽然距离待插值点较近，但它对待插值点的影响可能会被其他点屏蔽掉。

(3)在已知采样点Pi的邻域内，由于di≈0计算的误差将变得非常敏感，尤其是当两个项形式占优而又符号相反时，计算的误差更是如此。

因此，在实际应用时，必须对传统的距离反比加权法加以改进。考虑到地质特征数据的空间相关性，对距离反比加权法可附加地质体结构和影响距离两个限制条件，以提高空间几何或属性数据插值的合理性和精度。具体改进如下：

(1)地质体几何结构限制条件。在空间特征插值过程中，在选择影响待估计值的原始样本数据时只选取同一地层岩性地质内的样本；不是同一地层岩性地质体内的样本，即便距离很近，也不采用。即按照地质构造划分空间单元，并只在同一地层岩性地质体内选取样本点。

(2)邻近样本点的选择条件。选择待插值点的邻近点时，可考虑三个原则：一是距离原则，即根据地质属性数据的特征给出一距离r，在该距离之外的样本点对待插值点的估算无影响；二是点数原则，即给定一数据m，以距离待插值点最近m个样本点进行估算；三是利用Voronoi图求取待插值点的邻近点。

(3)建立数据点索引表，提高待插值点周围样本点的搜索效率，从而大大减少大数据量队的计算量。

5.3.3.2 普通Kriging插值方法

设研究区域为A，区域化变量(即欲研究的物理属性变量)为 {Z(x)∈A}，x表示空间位置。Z(x)在采样点xi(i=1，2，…，n)处的属性值(或称为区域化变量的一次实现)为Z(xi)(i=1，2，…，n)，则根据普通Kriging插值原理，未采样点x0处的属性值Z(x0)估计值是n个已知采样点属性值的加权和，即：

数字地下空间与工程三维地质建模及应用研究

其中，λi(i=1，2，…，n)为待求权系数。

假设区域化变量Z(x)在整个研究区域内满足二阶平稳假设：

(1)Z(x)的数学期望存在且等于常数：E[Z(x)]=m(常数)。

(2)Z(x)的协方差Cov(xi，xj)存在且只与两点之间的相对位置有关。或满足本征假设：

(3)E[Z(xi)-Z(xj)]=0。

(4)增量的方差存在且平稳：Var[Z(xi)-Z(xj)]=E[Z(xi)-Z(xj)]2。

依据无偏性要求：E[Z*(x0)]=E[Z(x0)]。

推导可得： 。

在无偏条件下使估计方差达到最小，即：

min{ Var[Z*(x0)-2μ( (λi-1))}，其中μ为拉格朗日乘子。

可求得求解权系数λii(=1，2，…，n)的方程组：

数字地下空间与工程三维地质建模及应用研究

求出诸权系数λii(=1，2，…，n)后，就可以求出采样点x0处的属性值Z*(x0)。

上述求解权系数λii(=1，2，…，n)的方程组中协方差Cov(xi，xj)若用变异函数γ(xi，xj)表示时，形式为：

数字地下空间与工程三维地质建模及应用研究

变异函数的定义为：

数字地下空间与工程三维地质建模及应用研究

由Kringing插值所得到的方差为：

数字地下空间与工程三维地质建模及应用研究

或

数字地下空间与工程三维地质建模及应用研究

A.Kringing插值法中的变异函数

变异函数是Kringing插值法插值的基础。插值中需要首先确定所研究的区域化变量的变异函数。假设研究的区域为A，区域A中有一区域化变量Z(x)，它在位置xi(i=1，2，…，N)上的一次采样为Z(xi)(i=1，2，…，N)，则Z(x)的变异函数的定义为：

数字地下空间与工程三维地质建模及应用研究

一个空间变量的空间变异性是指这个变量在空间中如何随着位置的不同而变化的性质。变异函数通过其自身的结构及其各项参数从不同的角度反映空间变异性，确定变异函数的过程就是一个对空间变异性进行结构分析的过程。

设h是一个模为r=|h|，方向为a的向量，如果存在着被向量h所隔开的Nh对观测数据点，则在a方向上相应于向量h的实验变异函数γ*(h)可表示为如下形式：

数字地下空间与工程三维地质建模及应用研究

其中，z(xi+h)和z(xi)分别位于点xi+h和xi(i=1，2，…，Nh)上的观测数据。

B.变异函数理论模型

当获取实验变异函数值后，需要先选择变异函数理论模型，然后对所选择的变异函数理论模型进行参数拟合，这一过程被称为“结构分析”。

变异函数理论模型参数一般包括：变程(range，一般用a表示)、基台(sill，一般用C(0)表示)、拱高(一般用C表示)、块金常数(Nugget，一般用C0表示)，如图5.9所示。

图5.9 理论变差函数

变程a表示了从空间相关性状态(|h|＜a)向不存在相关性状态(|h|＞a)转变的分界线；变异函数在原点处的间断性称为块金效应，相应的常数C0= (h)称为块金常数；基台C(0)具有协方差函数：C(h)=C(0)-γ(h)的二阶平稳区域化变量Z(x)的先验方差：Var{Z(x)}=C(0)=γ(h)；拱高C为变异函数中基台C(0)与块金常数C0之差：C=C(0)-C0。

C.变异函数理论模型分类

变异函数理论模型一般分为有基台值和无基台值两大类。有基台值的变异函数理论模型包括球状模型、指数模型、高斯模型等(图5.10)。最常用的是球状模型。球状模型公式为：

数字地下空间与工程三维地质建模及应用研究

图5.10 变异函数中的球状模型、指数模型、高斯模型

D.变异函数对空间变异性结构的反映

变异函数作为定量描述空间变异性的一种统计学工具，通过其自身的结构及其各项参数，从不同角度反映了空间变异性结构。利用变异函数可以对空间变量的连续性、相关性、变量的影响范围、尺度效应、原点处的间断性、各向异性等要素进行描述。

E.变异函数理论模型参数拟合

变异函数理论模型参数拟合就是利用原始采样点数据或实验变异函数取值对所选定的理论模型参数以特定的方法进行估计。拟合方法一般采用手工拟合法。

手工拟合就是依据实验变异函数的取值，一方面通过观察实验变异函数图；另一方面对所研究的区域化变量进行必要的分析，采用肉眼观察来确定变异函数模型参数，并对参数反复进行交叉验证，最终确定模型参数。其拟合的大致过程如下：

(1)首先对所研究的区域化变量进行必要的结构、背景等方面的分析，结合专家经验，确定变异函数理论模型。

(2)利用实验变异函数散点图确定变异函数参数中的块金常数、基台值、变程、各向异性角度以及各向异性比值。

(3)交叉验证。

三种插值方法的比较如下：

（1）拉格朗日插值评述

拉格朗日插值法无谓就是利用已知的个插值节点及其所在节点处的函数值，在每个插值节点处构造相应的插值基函数，再根据特定的线性关系将这个插值基函数进行线性组合，即得拉格朗日插值函数。

用几何的语言来描述这种方法就是将有限个点通过一条光滑的且与高度契合的次数不超过的函数来表示，其方法简洁明了，但是拉格朗日插值多项式在实际应用的过程中也暴露了本身存在的问题。如果对数据点的个数进行增加，那么原来我们所得的拉格朗日插值函数就毫无用处，必须从基函数构造重新开始整个过程。

在实际应用中节点的增减是特别普遍常见的，面临这种情况拉格朗日插值法就难免会面临较大的局限性，不仅会浪费时间，也会造成先前劳动力的浪费，这样就会极大的抑制大机器的生产，更加体现不出函数插值法的优化作用。

（2）牛顿插值评述

牛顿插值很好地解决了上述拉格朗日插值中的局限，即当增加节点时已得成果无法被利用的问题。牛顿插值法仅需在已有的多项式的基础上添加一项即可，这就很好的解决了上述拉格朗日插值方法所遇到的当增加节点时已得成果全部作废无法被继续使用的问题。

在日常实际问题解决过程中，利用有限个插值节点所构造的插值函数可能并不能达到我们所要求的插值精度。对于这个问题如果我们仅从增加插值节点个数这一方面来入手很可能会起到相反的作用。

（3）埃尔米特插值评述

通过对前面拉格朗日插值法和牛顿插值法的分析，我们可以很明显的观察到这两种插值方法的构造仅仅与插值节点以及插值节点处的函数值有关，并没有涉及到其它约束条件。但是如果插值条件不仅含有对节点处的函数值的约束，而且还增加对节点处的导数的限制，解决这一类问题的方法就要利用埃尔米特插值多项式。

对比上述拉格朗日插值方法和牛顿插值方法，埃尔米特插值具有较高的精度可以应用的领域更加宽泛，更加适应于实际问题的解决，所以在现实生活中也就凸现出高度的灵活性和适应性。

三种插值方法的比较

最近点插值

在一维空间中，最近点插值就相当于四舍五入取整。在二维图像中，像素点的坐标都是整数，该方法就是选取离目标点最近的点。计算方式如下：假设原图为A[aw,ah]，宽度为aw，高度为ah。目标图为B[bw,bh]，宽度为bw，高度为bh。已知A[aw,ah]的宽度，高度及其中每个点的颜色值，B[bw,bh]中每个点像素值的计算方式如下：

线性插值

线性插值是以距离为权重的一种插值方式。在一维空间中，假设有点A，B，其距离为LAB。A，B之间任意一点C的值为A*LBC/LAB+B*LAC/LAB。在二维空间中，需要在两个方向上做插值。如下图所示：

已知Q11，Q21，Q12，Q22，计算P点的值时，需要先由Q11和Q21插值得到R1，由Q12和Q22插值得到R2，再由R1和R2插值得到P。

兰索斯插值(lanczos)

一维的线性插值，是在目标点的左边和右边各取一个点做插值，这两个点的权重是由线性函数计算得到。而一维的兰索斯插值是在目标点的左边和右边各取四个点做插值，这八个点的权重是由高阶函数计算得到。