2.4 常用DEM插值算法

纵观DEM插值算法研究的发展历程，反距离加权插值算法、改进谢别德插值算法、径向基函数插值算法、克里格插值算法是常用的DEM插值算法。

2.4.1 反距离加权插值算法

反距离加权插值算法（Inverse Distance Weighted，IDW）最早是由气象学家和地质工作者提出的（王建等，2004），是空间数据插值最常见的算法之一。

反距离加权插值算法基于相近、相似的原理（卢华兴，2008），每个采样点都对插值点有一定的影响，即权重。权重随着采样点和插值点之间距离的增加而减弱，距离插值点越近的采样点的权重越大；当采样点在距离插值点一定范围以外时，权重可以忽略不计。任意插值点的值是各采样点的权重之和（王家华等，1999），如式（2.3）所示。

式中，zp为插值点的高程值，λi为第i个点的权重，di为第i个采样点到插值点的距离，d-u为距离衰减函数，幂指数u具有随着距离的增加减小其他位置影响的作用（de smith et al.，2007）。当u=0时，距离没有影响；当u=1时，距离的影响是线性的；当u≫1时，快速地减小遥远位置的影响。幂指数u通常取值为1或2（Lam，1983），但是大多数学者认为，幂指数采用2将取得更好的实验效果（Declercq，1996）。

反距离加权插值算法的计算易受采样点集群的影响，导致在采样点附近局部出现明显的隆起或凹陷的“牛眼”效应（Johns，1998）。

另外，由于反距离加权插值算法是一种精确性插值算法，因此插值生成的最大值和最小值只会出现在采样点处。这会直接导致出现山顶高程被降低、山谷高程被抬高的局部细节湮没。

基于反距离加权插值算法的缺点，许多学者提出了反距离加权插值算法的多种改进形式，用于克服上述缺点。

（1）给距离衰减函数增加一个平滑参数：一个微小的距离增量t，即，这样可能导致平滑，而不是精确插值。

（2）给距离衰减函数增加一个调和参数：基于最远点的距离R调整权重值，即

（3）使用其他的距离衰减函数，如高斯函数，即。

2.4.2 改进谢别德插值算法

谢别德插值算法（Modified Shepards Method，SPD）由南非地质学家Shepard最早提出，本质上是一种标准的导数距离加权过程（王金玲等，2010），权函数如式（2.5）所示。

式中，wi为权重，di为第i点距待插值点的距离，r为调整距离。

改进谢别德插值算法一般存在两种变化形式。

一是基于最远点的距离（在整个数据集中或在给定搜索半径范围内）调整权重。假设最远距离为r，那么修正的距离倒数加权函数如式（2.6）所示。

二是使用拟合的局部二次多项式调整权重，即参与距离倒数加权函数的高程值并不是原始采样点的高程值，而是使用拟合的局部二次多项式修正的高程值，如式（2.7）所示。

式中，zj为待插值点的高程值；，为插值点至采样点的距离；δ为平滑因子，当δ=0时为精确性插值，当δ≠0时为非精确性插值；Qi为二次多项式函数；u为权指数。

2.4.3 径向基函数插值算法

径向基函数插值（Radial Basis Functions，RBF）算法是一系列用于精确插值算子的统称（de smith et al.，2007）。它来源于Hardy的多面函数法，其插值原理是任何一个表面都可以使用多个曲面的线性组合逼近（卢华兴，2008）。在多数情况下，径向基函数插值算法与地统计插值算法相似，但具有不需要分析半变异函数模型的优点，而且不需要有关采样点的任何假设（除了非共线性）。

在通常情况下，径向基函数插值算法可以表述为两个部分之和（Mitasova and Mitas，1993），即

式中，z p为插值点的高程值；λi为第i个点的权重；di为第i个采样点到插值点的距离；φ（di）为径向基函数，它代表第j个核函数对多层叠加曲面的贡献；f j（x）为“趋势”函数，是次数小于m的基本多项式函数，由于fj（x）并不能提高插值的精度，因此在插值过程中不考虑“趋势”函数的影响（de smith et al.，2007）。径向基函数插值算法的解算过程可以使用矩阵符号表示为如下步骤。

步骤1：计算源数据中所有（x，y）点对的点间距离构成的n×n阶矩阵D；

步骤2：对D中的每个矩阵值应用选择的径向基函数φ（·），从而产生一个新的矩阵Φ；

步骤3：用单位列矢量和单位行矢量增大矩阵Φ，并且在位置（n+1，n+1）处插入零值，称这个增广矩阵为A；

步骤4：计算从格网点p到用来创建D的每个源数据点间的距离构成的列矢量r；

步骤5：将选择的径向基函数应用于r中的每个距离产生一个列矢量φ，然后生成一个n+1阶的列矢量c，它由φ加上元素1构成，即

步骤6：计算矩阵积b=A-1c。这样就给出了用于计算p点的估计值的n个权重，即

径向基函数插值算法可以选用许多不同的径向基函数（见表2-3）。

表中，c为光滑因子，一般由用户指定。c的值取决于对插值结果产生重要影响的采样点的数目、高程、空间分布等因素（Rippa，1999）。

对于如何确定c，没有普遍认可的方法，但也有一些学者提出了各种方法。Hardy（1971）使用c=0.815d，其中，，di为第i个点到其最近邻的距离。Franke（1982）使用替换d，其中，D是数据集最小外接圆的直径，于是建议使用。Foley（1987）做出了和Franke类似的建议，不过使用数据集的最小外接矩形的边长代替最小外接圆的直径。Rippa （1999）提出了使用递归算法寻找使得插值表面全局误差最小的参数c的方法。Aguilar等（2005）认为，在MQF插值算法和MLF插值算法中应当使用接近于零的光滑因子；在IMQF插值算法、NCSF插值算法、TPSF插值算法中则应当使用非常大的光滑因子，因为在IMQF插值算法、NCSF插值算法、TPSF插值算法中如果使用较小的光滑因子，将会产生显著的数值不稳定性。

径向基函数插值算法作为一种精确的插值算法，不同于局部多项式插值算法。局部多项式插值算法作为一种非精确的插值算法，并不要求表面经过所有的采样点。径向基函数插值算法和同为精确插值算法的反距离加权插值算法的不同之处在于，反距离加权插值算法不能计算出高于或低于采样点的插值点的值，而径向基函数插值算法则可以计算出高于或低于采样点的插值点的值。

2.4.4 克里格插值算法

克里格插值算法也称局部估计插值或空间局部插值，是地统计学的两大主要内容之一（张景雄，2008）。地统计学源于20世纪50年代Krige在地质和采矿业方面的工作，1963年法国学者Matheron发表了专著《应用地质统计学》，提出了区域化变量理论，并且给出了地统计学的概念：以区域化变量理论为基础，以变异函数为主要工具，研究在空间分布上既有随机性又有结构性的自然现象的科学（侯景儒，1998）。

1.区域化变量

区域化变量是以空间采样点x的三维直角坐标（xu，xv，xw）为自变量的随机场函数Z（xu，xv，xw）=Z（x），在对其进行一次观测后，就得到随机场函数Z（x）的一个具体实现z（x）。在空间的每个点取某个确定的数值后，当由一个点移动到下一个点时，函数具体实现z（x）是变化的。

区域化变量具有随机性和结构性的双重特征。随机性是指区域化变量在具体实现时表现出一定的不规则特征；结构性是指区域化变量在不同的空间方位具有某种程度的空间自相关性。

地形表面作为一个连续的随机场表面，符合区域化变量的双重特征，因此以区域化变量理论为基础的“地统计学”在地形建模、空间分析等方面的应用方兴未艾。

2.半变异函数

半变异函数是一种空间变量相关性的定量化描述模型。当空间采样点在一维轴x上变化时，区域化变量在x和x+h处的值为Z （x）和Z （x+h），两者之差的方差的一半定义为区域化变量在x轴上的半变异函数，记为γ（x，h）。

在二阶平稳假设下，有

那么γ（x，h）可以改写成

从式（2.17）可知，半变异函数依赖于两个自变量x和h，当半变异函数γ（x，h）与位置x无关时，它仅依赖于分隔两个采样点之间的距离，那么γ（x，h）可以改写成γ（h），即

在通常情况下，半变异函数值随着采样点间距的增加而增大，并在到达某一个间距值后趋于稳定（见图2-4）。半变异函数具有3个重要的参数，分别是块金值（Nugget）、基台值（Sill）和变程（Range），它们表示区域化变量在一定尺度上的空间变异性和相关性。

（1）块金值。根据半变异函数的定义，理论上当h=0时，半变异函数值应等于0。但是，由于采样误差等原因，即使两个采样点之间的距离h很小，其变量依然存在差异，这表示区域化变量在小于观测尺度时的非连续性变异。

（2）基台值。基台值表示半变异函数随着间距递增到一定程度时出现的平稳值，即C0+C。其中，C称为结构方差（或拱高），在数值上等于基台值与块金值之间的差值，代表由于样本数据中存在空间相关性而引起的方差变化的范围。

（3）变程。变程表示半变异函数达到基台值时的距离，反映了空间采样点的自相关距离尺度。在变程距离之内，空间上越近的点之间的相关性越大；当h大于变程时，空间采样点之间不具备自相关性，除非半变异函数具有周期性变化特征。更为重要的是，变程表示了空间插值的极限距离，选择在变程范围内的采样点参与插值才有意义（张仁铎，2005）。

在实际插值估计中，由于空间采样点是离散的，无法获取半变异函数γ（h）的理论值，所以需要通过实验方法获得实验半变异函数值γ*（h），即

式中，N （h）是近似相隔h的采样点对的数目。

之后，根据实验半变异函数值选择合适的理论半变异函数模型，并且拟合半变异函数模型的基本参数。

理论半变异函数模型包括几种简单的模型，如图2-5所示。

（1）线性模型（LINE）：

（2）球形模型（SPHERE）：

当C0=0、C=1时，球形模型称为标准球形模型。

（3）指数模型（EXP）：

当C0=0、C=1时，指数模型称为标准指数模型。

（4）高斯模型（GAUSS）：

当C0=0、C=1时，高斯模型称为标准高斯模型。

3.普通克里格插值算法

克里格插值算法包括简单克里格插值算法、普通克里格插值算法、通用克里格插值算法、指标克里格插值算法、概率克里格插值算法、分离克里格插值算法、分层克里格插值算法、联合克里格插值算法、因子克里格插值算法等20多种不同的变形形式（de Smith et al.，2007），但是所有的克里格插值算法都是基于式（2.24）的微小变异。

式中，m为整个区域内所有采样数据的均值，λi是克里格权重，n是以x0为中心的、指定搜索区域内的参与克里格插值的采样点个数，m（x0）是指定搜索区域内的采样点均值。

当m为已知参数时，克里格插值称为简单克里格插值；当m为未知参数时，克里格插值称为普通克里格插值。

从式（2.24）可以看出，克里格插值的关键在于求解克里格权重λi，并且克里格权重λi必须满足无偏条件，使估计方差最小。

其中，无偏条件的数学表达式为

估计方差表示为

式中，C（xi-xj）=Cov[Z（xi）-Z （xj）]为协方差函数，协方差函数和半变异函数之间具有如下关系：

式中，C（0）为区域化变量的Z （x）的方差。

要使估计方差在无偏条件下最小，则问题变为一个求解条件极值的方程，可以采用标准拉格朗日乘数法求解。依据拉格朗日原理构造函数F，即

式中，μ为拉格朗日乘数法。

分别对F求λi和μ的偏导，可得

式（2.29）是一个n+1阶线性方程组，有n个未知数λi和1个未知数μ。因此，可以建立n+1维线性方程组，即

将已知采样点数据代入式（2.30），即可解得λi。