8.7. 空间插值

整体插值方法 #

1）边界内插方法

边界内插方法假设任何重要的变化发生在边界上，边界内的变化是均匀的，同质的，即在各方向都是相同的。 这种概念模型经常用于土壤和景观制图，可以通过定义“均质的”土壤单元、景观图斑， 来表达其它的土壤、景观特征属性。

边界内插方法最简单的统计模型是标准方差分析（ANOVAR）模型：

式中，z是在 x_0 位置的属性值，μ是总体平均值， α_k 是k类平均值与μ的差， ε为类间平均误差（噪声）。

该模型假设每一类别k的属性值是正态分布； 每类k的平均值（ μ+α_k ）由一个独立样品集估计，并假设它们是与空间无关的； 类间平均误差ε假设所有类间都是相同的。

评价分类效果的指标是，为类间方差，为总体方差，比值越小分类效果越好。 分类效果的显著性检验可以用F检验。

实质上，边界内插方法的理论假设是：

• 属性值z在“图斑”或景观单元内是随机变化的，不是有规律的；

• 同一类别的所有“图斑”存在同样的类方差（噪声）；

• 所有的属性值都呈正态分布；

• 所有的空间变化发生在边界上，是突变而不是渐变。

在使用边界内插时，应仔细考虑数据源是否符合这些理论假设。

2）趋势面分析

某种地理属性在空间的连续变化，可以用一个平滑的数学平面加以描述。 思路是先用已知采样点数据拟合出一个平滑的数学平面方程，再根据该方程计算无测量值的点上的数据。 这种只根据采样点的属性数据与地理坐标的关系，进行多元回归分析得到平滑数学平面方程的方法， 称为趋势面分析。 它的理论假设是地理坐标 （x,y） 是独立变量，属性值Z也是独立变量且是正态分布的， 同样回归误差也是与位置无关的独立变量。

多项式回归分析是描述长距离渐变特征的最简单方法。 多项式回归的基本思想是用多项式表示线、面，按最小二乘法原理对数据点进行拟合。 线或面多项式的选择取决于数据是一维的还是二维的。

用一个简单的示例来说明，地理或环境调查中特征值z沿一个断面在 x_1 ， x_2 … x_n 处采样，若z值随x值增加而线性增大， 则该特征值的长期变化可以用下面一个回归方程进行计算：

其中， b_0 ， b_1 为回归系数，ε为独立于x的正态分布残差（噪声）。

然而许多情况下，不是以线性函数，而是以更为复杂的方式变化，则需用二次多项式进行拟合：

对于二维的情况，XY坐标的多元回归分析得到的曲面多项式，形式如下：

前三种形式分别是：

平面

斜平面

二次曲面

其中，p是趋势面方程的次数。

是趋势面多项式正常情况下的最少项数个数。 零次多项式是平面，有1个项数； 一次多项式是斜平面，有3个项数； 二次曲面有6个项数，三次趋势面有10个项数。

计算系数 b_i 是一个标准的多元回归问题。 趋势面分析的优点是非常容易理解，至少是在计算方面。 另外大多数情况下可用低次多项式进行拟合，但给复杂的多项式赋与明确的物理意义比较困难。

趋势面是个平滑函数，很难正好通过原始数据点，除非是数据点少且趋势面次数高才能是曲面正好通过原始数据点， 所以趋势面分析是一个近似插值方法。 实际上趋势面最有成效的应用是揭示区域中不同于总趋势的最大偏离部分，所以趋势面分析的主要用途是， 在使用某种局部插值方法之前，可用趋势面分析从数据中去掉一些宏观特征，不直接用它进行空间插值。

趋势面拟合程度的检验，同多元回归分析一样，可用F分布进行检验，其检验统计量为：

其中，U为回归平方和，Q为残差平方和（剩余平方和），p为多项式项数（不包括常数项 b_0 ）， n为使用数据点数目。 当 F>F_a 时，趋势面显著，否则不显著。

3）变换函数插值

根据一个或多个空间参量的经验方程进行整体空间插值，也是经常使用的空间插值方法， 这种经验方程称为变换函数。 下面以一个研究实例进行说明。

冲积平原的土壤重金属污染与几个重要因子有关，其中距污染源（河流）的距离，和高程两个因子最重要。 一般情况，携带重金属的粗粒泥沙沉积在河滩上，携带重金属的细粒泥沙沉淀在低洼的在洪水期容易被淹没的地方， 而那些洪水频率低的地方，由于携带重金属污染泥沙颗粒比较少，受到污染轻。 由于距河流的距离和高程是比较容易得到的空间变量，可以用各种重金属含量与它们的经验方程进行空间插值， 以改进对重金属污染的预测。 本例回归方程的形式如下：

式中是z(x)某种重金属含量（ppm）， b_0 … b_n 是回归系数， p_1 … p_n 是独立空间变量，本例 p_1 是距河流的距离因子， p_2 是高程因子。

这种回归模型通常叫做转换函数，大多数GIS软件都可以计算。 转换函数可以应用于其他独立变量， 如温度、高程、降雨量和距海、植被的距离关系可以组合为一个超剩含水量的函数。 地理位置及其属性可以尽可能多的信息组合成需要的回归模型，然后进行空间插值。 但应该注意的一点是，必须清楚回归模型的物理意义。 还要指出的是所有的回归转换函数都属于近似的空间插值。

整体插值方法通常使用方差分析和回归方程等标准的统计方法，计算比较简单。 其它的许多方法也可用于整体空间插值，如傅立叶级数和小波变换，特别是遥感影象分析方面， 但它们需要的数据量大。

局部插值方法 #

局部插值方法只使用邻近的数据点来估计未知点的值，包括几个步骤：

1. 定义一个邻域或搜索范围；

2. 搜索落在此邻域范围的数据点；

3. 选择表达这有限个点的空间变化的数学函数；

4. 为落在规则格网单元上的数据点赋值。重复这个步骤直到格网上的所有点赋值完毕。

使用局部插值方法需要注意的几个方面是：所使用的插值函数； 邻域的大小、形状和方向； 数据点的个数； 数据点的分布方式是规则的还是不规则的。

1）最近邻点法： 泰森多边形 方法

泰森多边形 （Thiessen， 又叫Dirichlet或Voronoi多边形）采用了一种极端的边界内插方法， 只用最近的单个点进行区域插值。 泰森多边形 按数据点位置将区域分割成子区域，每个子区域包含一个数据点，各子区域到其内数据点的距离小于任何到其它数据点的距离，并用其内数据点进行赋值。 连接所有数据点的连线形成Delaunay三角形，与不规则三角网TIN具有相同的拓扑结构。

GIS和地理分析中经常采用 泰森多边形 进行快速的赋值，实际上 泰森多边形 的一个隐含的假设是任何地点的气象数据均使用距它最近的气象站的数据。而实际上，除非是有足够多的气象站，否则这个假设是不恰当的，因为降水、气压、温度等现象是连续变化的，用 泰森多边形 插值方法得到的结果图变化只发生在边界上，在边界内都是均质的和无变化的。

2）移动平均插值方法：距离倒数插值

距离倒数插值方法综合了 泰森多边形 的邻近点方法和趋势面分析的渐变方法的长处， 它假设未知点 x_0 处属性值是在局部邻域内中所有数据点的距离加权平均值。 距离倒数插值方法是加权移动平均方法的一种。 加权移动平均方法的计算公式如下：

式中，权重系数由函数计算，要求当时，一般取倒数或负指数形式。 其中最常见的形式是距离倒数加权函数，形式如下：

式中， x_j 为未知点， x_i 为已知数据点。

加权移动平均公式最简单的形式称为线性插值，公式如下：

距离倒数插值方法是GIS软件根据点数据生成栅格图层的最常见方法。 距离倒数法计算值易受数据点集群的影响， 计算结果经常出现一种孤立点数据明显高于周围数据点的“鸭蛋”分布模式， 可以在插值过程中通过动态修改搜索准则进行一定程度的改进。

3）样条函数插值方法

在计算机用于曲线与数据点拟合以前，绘图员是使用一种灵活的曲线规逐段的拟合出平滑的曲线。 这种灵活的曲线规绘出的分段曲线称为样条。 与样条匹配的那些数据点称为桩点，绘制曲线时桩点控制曲线的位置。 曲线规绘出的曲线在数学上用分段的三次多项式函数来描述这种曲线，其连接处有连续的一阶和二阶连续导数。

样条函数是数学上与灵活曲线规对等的一个数学等式，是一个分段函数，进行一次拟合只有与少数点拟合， 同时保证曲线段连接处连续。 这就意味着样条函数可以修改少数数据点配准而不必重新计算整条曲线，趋势面分析方法做不到这一点， （图16）。

（a：当二次样条曲线的一个点位置变化时，只需要重新计算四段曲线；

b：而一次样条曲线的一个点位置变化时，只需要重新计算两段曲线）

一般的分段多项式 p(x) 定义为：

p(x)=p_i (x)

x_i < x < x_i+1 ( i=1, 2, 3…, k-1)

( j=0, 1, 2, ..., r-1; i=1, 2, ... , k-1 )

x_1 , ...,x_k-1 将区间 x_0 ，x_k 分成 k 个子区间，这些分割点称“断点”，曲线上具有这些 x 值的点称为“节”。函数 p_i(x) 为小于等于m次的多项式。 r 项用来表示样条函数的约束条件：

r=0 ，无约束；

r=1 ，函数连续且对它的导数无任何约束；

r=m-1 ，区间 [x_x_k] 可用一个多项式表示；

r=m，约束条件最多。

m=1,2,3时的样条分别为一次、二次、三次样条函数，其导数分别是0阶、1阶、2阶导数，二次样条函数的每个节点处必须有一阶连续导数，三次样条函数的每个节点初必须有二阶连续导数。 r=m 的简单样条只有 k+m 个自由度， r=m=3 有着特殊的意义，因为它是三次多项式，该函数首次被人们称为样条函数。术语“三次样条”用于三维情况，此时进行曲面内插，而不是曲线内插。

由于离散子区间的范围较宽，可能是一条数字化的曲线，在这个范围内计算简单样条会引起一定的数学问题， 因此在实际应用中都用 B 样条—一种特殊的样条函数。 B 样条是感兴趣区间以外均为零的其它样条的和，因此可按简单的方法用低次多项式进行局部拟合。

B 样条经常用于数字化的线划在显示之前进行平滑处理，例如土壤、地质图上的各种边界，传统的制图总希望绘出较平滑的曲线。但是用 B 样条做多边形边界平滑也存在一些问题，特别是多边形面积和周长的计算，结果会与平滑前的不同。

综上所述，样条函数是分段函数，每次只用少量数据点，故插值速度快。 样条函数与趋势面分析和移动平均方法相比，它保留了局部的变化特征。 线性和曲面样条函数都在视觉上上得到了令人满意的结果。 样条函数的一些缺点是： 样条内插的误差不能直接估算， 同时在实践中要解决的问题是样条块的定义以及如何在三维空间中将这些“块”拼成复杂曲面， 又不引入原始曲面中所没有的异常现象等问题。

4）空间自协方差最佳插值方法：克里金插值

前面介绍的几个插值方法对影响插值效果的一些敏感性问题仍没有得到很好的解决， 例如趋势面分析的控制参数和距离倒数插值方法的权重对结果影响很大，这些问题包括：

• 需要计算平均值数据点的数目；

• 搜索数据点的邻域大小、方向和形状如何确定；

• 有没有比计算简单距离函数更好的估计权重系数的方法；

• 与插值有关的误差问题。

为解决这些问题， 法国地理数学学家GeorgesMatheron和南非矿山工程师D.G.Krige研究了一种优化插值方法， 用于矿山勘探。 这个方法被广泛地应用于地下水模拟、土壤制图等领域，成为GIS软件地理统计插值的重要组成部分。 这种方法充分吸收了地理统计的思想，认为任何在空间连续性变化的属性是非常不规则的， 不能用简单的平滑数学函数进行模拟，可以用随机表面给予较恰当的描述。 这种连续性变化的空间属性称为“区域性变量”，可以描述象气压、高程及其它连续性变化的描述指标变量。 这种应用地理统计方法进行空间插值的方法，被称为克里金（Kriging）插值。 地理统计方法为空间插值提供了一种优化策略，即在插值过程中根据某种优化准则函数动态的决定变量的数值。 Matheron，Krige等人研究的插值方法着重于权重系数的确定，从而使内插函数处于最佳状态，即对给定点上的变量值提供最好的线性无偏估计。

克里金插值方法的区域性变量理论假设任何变量的空间变化都可以表示为下述三个主要成分的和（图8-17）：

1. 与恒定均值或趋势有关的结构性成分；

2. 与空间变化有关的随机变量，即区域性变量；

3. 与空间无关的随机噪声项或剩余误差项。

1. 地形的平均特性；

2. 空间相关的不规则变化；

3. 随机的、局部的变化

   另 x 为一维、二维或三维空间中的某一个位置，变量 z 在 x 处的值可又下式计算：

式中，m(x)是描述z(x)的结构性成分的确定性函数；是与空间变化有关的随机变化项，即区域性变量； 是剩余误差项，空间上具有零平均值、与空间无关的高斯噪声项。

克里金方法的第一步是确定适当的 m(x) 函数， 最简单的情况是 m(x) 等于采样区的平均值，距离矢量 h 分离的两点 x ， x+h 之间的数学期望等于零：

式中 z(x) , z(x+h) 是随机变量 z 在 x,x+h 处的值， 同时还假设两点之间的方差只与距离 h 有关，于是：

式中，是一种函数，称为半方差函数。

区域性变量理论的两个内在假设条件是差异的稳定性和可变性，一旦结构性成分确定后， 剩余的差异变化属于同质变化，不同位置之间的差异仅是距离的函数。 这样，区域性变量计算公式可以写成下式的形式：

半方差的估算公式如下：

式中， n 为距离为 h 的采样点对的数目（ n 对点）， 采样间隔 h 也叫延迟。 对应于 h 的的图被称为“半方差图”。 图8-18表示一个典型的半方差图。

半方差是定量描述区域性变化的第一步，它为空间插值、优化采样方案提供了有益的信息。 为了求得半方差图，必须先得到拟合半方差的理论模型，在半方差理论模型中：

1. 延迟 h 的值较大的部分曲线呈水平方向。曲线的水平部分成为“梁（Sill）”。说明在延迟的这个范围内数据点没有空间相关性，因为所有的方差不随距离增减而变化。

2. 曲线从的低值升到梁为止的延迟范围，称为“变程（Range）”。变程是半方差图最重要的部分，因为它描述了与空间有关的差异怎样随距离变化的。在变程范围内距离越近的点具有更相近的特征。变程给移动加权平均方法提供了一个确定窗口大小的方法。很显然，数据点和未知点之间的距离大于变程范围，表明该数据点与未知点距离太远，对插值没有作用。

3. 图中的拟合模型没有通过原点，而是在的正方向与坐标轴相截。

按半方差计算公式，当h=0时，必须为零。 模型中出现的正值是剩余误差的估计值，它是与空间无关的噪声。 称为“核(Nugget)”方差，是观测误差的和距离间隔很小的情况下的空间变化的组合。

当存在明显的变程和梁，同时核方差也很重要但数值不太大的情况下， 可用球面模型进行半方差拟合（图8-19-1）。 公式是：

0 < h < a

h >= a

h=0

式中，a是变程，h是延迟， c_0+c_1 为梁。 一般情况下用球面模型拟合效果比较理想。

如果存在明显的核方差和梁，而没有渐变的变程，则可用指数模型（图8-19-2）进行拟合：

如果核方差相对于与空间变化有关的随机变化很小的情况下，最好使用比较弯曲的曲线，如高斯曲线： （图8-19-4）

如果空间变化随变程渐变，但没有梁，则可用线性模型（图8-19-3）进行拟合：

式中，b为线的斜率。当变程的大小远超过人们希望的插值范围时，也用线性模型。

前面的讨论都假设地表特征在各个方向都是相同的，然而许多情况下空间变化中的都具有明显的方向性， 这是就要用不同参数的模型来拟合半方差图。

拟合后的半方差图重要的用途是确定局部内插需要的权重因子。 确定的过程与加权移动插值方法类似，但不是按一种固定的函数计算， 而是按采样点数据的半方差图的统计分析原理计算。 即：

权重的选择应使是无偏估计，且估计的方差小于观测值的其它线性组合产生的方差。

的最小方差为：

只有下式成立时，才可获得的最小方差：

式中，是z在采样点 x_i ,

x_j 之间的半方差；是采样点 x_j 和未知点 x_0 之间的半方差，这两个量均可从已拟合模型的半方差图上得到。量计算最小方差需要的拉格朗日算子。

这个方法叫常规克里金插值。 它是一个精确插值模型，内插值或最佳局部均值与数据点上的值一致。 制图中常用比采样间隔更细的规则格网进行插值，内插值又可用前边提到的方法转换成等值线图。 与此类似，估计的误差，又叫克里金方差，可以用来制图以反映在整个研究区内插值结果的可靠性。

一个克里金插值的示例：（图8-20）

计算图中(x=5,y = 5)处未知点I0的克里金方法的权重，使用球面模型对半方差进行拟和， 参数 c_0= 2.5 ， c_1 =7.5 ； 变程 a= 10.0。 这5个采样点的数据分别是：

解法过程用矩阵的表示如下：

其中，A为数据点之间的半方差矩阵，b是每个数据点与未知点之间的半方差向量， λ 为要计算的权重系数向量， Φ 为解方程的拉格郎日算子。

首先计算这5个数据点之间的距离矩阵：

*I 1 2 3 4 5*
*I1 0.0 5.099 9.899 5.000 3.162*
*I2 5.099 0.0 6.325 3.606 4.472*
*I3 9.899 6.325 0.0 5.0 7.211*
*I4 5.0 3.606 5.0 0.0 2.236*
*I5 3.162 4.472 7.211 2.236 0.0 *

和它们与未知点之间的距离向量：

*I I0*
*I1 4.243*
*I2 2.808 *
*I3 5.657 *
*I4 1.0*
*I5 2.0*

将上述数值带入球面模型，得到相应的半方差（矩阵A, b）：

*A= I 1 2 3 4 5 6*
*1 2.500 7.739 9.999 7.656 5.939 1.000*
*2 7.739 2.500 8.677 6.381 7.196 1.000*
*3 9.999 8.677 2.500 7.656 9.206 1.000*
*4 7.656 6.381 7.656 2.500 4.936 1.000*
*5 5.939 7.196 9.206 4.936 2.500 1.000*
*6 1.000 1.000 1.000 1.000 1.000 0.000*
*b = I 0 *
*1 7.151*
*2 5.597*
*3 8.815*
*4 3.621*
*5 4.720*
*6 1.000*

注意额外的第6行和第6列，是为了保证权重之和为1。

计算A的逆矩阵，得：

A:sup:`-1`= I 1 2 3 4 5 6

1 -0.172 0.050 0.022 -0.026 0.126 0.273

2 0.050 -0.167 0.032 0.077 0.007 0.207

3 0.022 0.032 -0.111 0.066 -0.010 0.357

4 -0.026 0.077 0.066 -0.307 0.190 0.030

5 0.126 0.007 -0.010 0.190 -0.313 0.134

6 0.273 0.207 0.357 0.003 0.134 -6.873

于是，权重 λ 为：

I 权重系数 距离

1 0.0175 4.423

2 0.2281 2.828

3 -0.0891 ∑ = 1 5.657

4 0.6437 1.000

5 0.1998 2.000

6 0.1182 φ

得未知点插值后得数值为：

z(I0) =0.0175*3 + 0.2281*4 - 0.0891*2 + 0.6437*4 + 0.1998*6

= 4.566

估计方差为：

σ_e_2 = [0.0175* 7.151+ 0.2281* 5.597 - 0.0891*8.815

+ 0.6437 *3.621 + .1998 *4.720 ] +φ

= 3.890 + 0.1182

= 4.008

克里金插值方法的目的是提供确定权重系数最优的方法和并能描述误差信息。 由于克里金点模型（常规克里金模型）的内插值与原始样本的容量有关，当样本少的情况下， 采用简单的点常规克里金插值的内插结果图会出现明显的凹凸现象。 可以通过修改克里金方程以估计子块 B 内的平均值来克服这一缺点。 该方法叫块克里金插值，该方法对估算给定面积试验小区的平均值或对给定格网大小的规则格网进行插值比较适用。

子块B内z的均值为：

式中， S_B 为子块B的面积。z(x)仍用相同的估计公式：

最小方差仍为：

但成立条件则变为：

块克里金插值估算的方差结果常常小于点克里金插值，所以生成的平滑插值表面不会发生点模型的凹凸现象。