二、局部内插法

1. 样条函数：计算机用于曲线与数据点拟合以前，制图员用曲线规逐段地拟合出平滑的曲线。这种灵活的曲线规绘成的分段曲线称为样条。与样条匹配的那些数据点称为桩点，绘制曲线时桩点控制曲线的位置。曲线规绘出的曲线近似于分段的三次多项式曲线，该曲线为连续而有一阶和二阶连续导数。

样条函数是灵活曲线规的数学等式，也是分段函数，一次拟合只有少数数据点配准，同时保证曲线段的连接处为平滑连续曲线。这就意味着样条函数可以修改曲线中的某一段而不必重新计算整条曲线，趋势面和傅里叶级数都作不到这一点（图 4—7）。

图中（a）表示 1 个点移动后，二次样条必须重新计算 4 个点；（b）表

示线性样条只需重新计算 2 个点。分段多项式 p（x）的定义为：

P（x）=Pi（x）xi＜x＜xi+1（i=1，2，⋯，k－1）（4—13）

p ) （xi ） = P ^(j) （x ）

i(j

i+1 i

（j=0，1，⋯，r－1；i=1，2，⋯，k-1）（4-14）

x1，⋯xk-1 将区间 x0，xk 分成 k 个子区间，这些分割点称为“断点”， 曲线上具有这些 x 值的点又可称为“节”。函数 pi（x）小于等于 m 次多项式。r 项用来表示样条函数约束条件。r=0 时，无约束；r=1 时，函数连续且对它的导数无任何约束；r=m-1 时，区间 x0，xk 可用一个多项式表示。所以r=m 时，约束条件最多。m=1，2，3 时的样条分别称为线性、二次、三次样条函数，其导数分别为 0，1，2 阶导数。于是二次样条函数的每个节点处必须有一阶连续导数；三次样条则应有二阶连续导数等。r=m 的简单样条只有 k+m 自由度，r=m=3 时，因为它是三次分段多项式函数，该函数首次被人们称为样条函数。术语“三次样条”用于三维情况，此时曲面代替了曲线进行内插。

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

数字化的线划在显示之前常用样条进行平滑。例如土壤、地质图上的各种边界，传统的制图总希望绘出较平滑的曲线，因此常用 B 样条来处理复杂形状边界问题。综上所述，样条既可用于精确的内插（通过所有的数据点） 也可用于平滑处理。样条函数是分段函数，每次只用少量数据点，故内插速度快。样条函数与趋势面相反，它保留了微地物特征。线性和曲面样条函数都在美学上得到令人满意的结果。但样条内插的误差不能直接估算，同时在实践中要解决的问题是样条块的定义以及如何在三维空间中将这些“块”拼成复杂曲面而又不至于引入原始曲面中所没有的异常现象等问题。

1. 移动平均法：在未知点 x 处内插变量 z 值时，最常用的方法之一是在局部邻域（或称窗口）中计算各数据点的平均值。在数据点沿某断面规则分布的情况下，可用移动平均法的最简单形式即对窗口中心点 x 的移动平均值，其计算式为：

^ 1 ∑n

z(x) = z(x i ) (4 - 15)

i=1

二维移动平均值可用相同的公式（4－15）计算，但位置 xi 应被坐标矢量 xi 代替。

窗口的大小对平滑的输出形式有决定性的影响。窗口增宽将增强长距离变化的特征而减少短距离变化特征。

前已述及，观测点的相互位置越是接近，其相似性越强；点距越远，即使落在同一个制图单元的轮廓线内，其相似性也不如接近的点。特定的采样点对未知点的平均内插值的确定，要按采样点到内插点的距离函数加权计算。加权移动平均值可按下式计算：

z(x _i )d ^-2

^

z(x j ) = i=1

−2 ij

(4 - 16)

i=1

式中x _j 是曲面上内插的点，通常这些点都处于规则格网上。d -2是距离的

平方倒数权。

加权平均内插的结果随使用的函数及其参数、数据采集点的定义域窗口的大小不同而变化，见图 4—8。

移动平均法计算值易受数据点聚群的影响，亦易为观测的二维趋势支配与否所影响，有人提出用加权距离最小二乘法。函数域（窗口）不仅影响平均值的估计，还关系到内插的时间。对窗口大小的要求是窗口内应包括数据的极大极小值，使计算效率与精度达到要求。计算的点数 n 可在 4—12 之间， 通常为 6—8 点。特别是原始数据为规则的栅格数据是常用 6—8 点。另一方面，人们在计算平均值时可能将数据点数固定一致，这样处理对规则的栅格数据不会有问题，而对不规则分布的数据来说，得到固定数量的数据点就要不断地改变窗口的大小、形状和方向。

1. 空间自协方差最佳内插法：局部加权移动平均内插法能在很多场合下使用，并能取得良好效果，但仍有一些重要问题留待解决，例如①函数域（窗口）应多大？②窗口的形状和方向应如何确定才能达到最佳内插效果？③有无比简单距离函数更好的方法来估计权的大小？④内插的误差有多大等。为解决这些问题，地理数学家们把发展内插方法着重于权的选择，从而使内插函数处于最佳状态，亦即对给定点上的变量值提供最好的无偏估计。

这种方法的依据是：由于任何地质、土壤、水文等区域特性变量过于杂乱，不能用平滑数学函数进行模拟，但能用随机表面给予较适当的描述。该方法的内插过程是：首先探查区域性变量的随机状况，然后模拟这些变量的随机状况，最后用前两步产生的信息估计内插的权因子λi。

克里金法作为成功的最佳内插法，其依据是有效地假设变量的统计特征，这种假设是区域变化理论的内容之一，也是克里金法的基础。

区域变化理论假设任何变量的空间变化都可以用下述三个主要成分的和来表示：①与均值即趋势有关的结构成分；②与局部变化有关成分；③随机噪声项即剩余误差项。令 x 是一维、二维或三维空间中的某一位置，变量 z 在 x 处的值由下式计算：

z（x）=m（x）+ε′（x）+ε″（4—17）

式中 m（x）是描述 x 处 z 的结构成分的确定性函数；ε′（x）是随机指示项，表示局部变化；ε″是剩余误差项，空间上具有零平均值、σ ² 方差的独立高斯噪声项（图 4－9）。

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

E[z（x）-z（x+h）]=0（4—18）

同时还假设差分的方差只与两位置之间的距离 h 有关，于是：

E[z（x）-z（x＋h）²]=E{[ε′（x）-ε′（x+h）²]}=2r（h）（4—19） 式中 r（h）是一种函数，称为半方差函数。

差分的可变性和稳定性两个条件确定了区域变化理论的内在要求，这就是说一旦结构影响确定以后，变量在变化范围内的剩余变化是同性变化，因此位置之间的差异仅是位置间距离的函数。这样就可用下式估算采样数据的半方差：

_^ 1 n 2

r(h) = ∑[z(x i ) - z(x i + h)] (4 - 20)

i=1

式中 n 是距离为 h 的采样点的对数（n 对点），采样间

^

隔h也叫延迟。对应于h的r （h）的图形称为半可变图。计算半方差是确

定最佳内插权因子的重要步骤。

图 4－10 表示典型半可变图。该图表示野外测量的表土层（0—20cm）中含粘土的百分比，所测数据是在 6×25 个正方形采样格网上观测的，格网长轴方向为东南—西北向，小区域的格网像元为 125m×125m，150 个像元中的136 个像元用钻孔法采样，其余 14 个像元因野外障碍未采样而遗漏数据。该地区在地貌上分成两类土地景观单元：排水良好的阶地和排水不良的洪积平原。以航空像片判读及 1∶50000 比例尺的地面测量为基础的土地调查勾绘出4 个制图单元。后来又进行了更详细的 1∶25000 比例尺的地面测量，将该地区分成 6 个类别单元（图 4－11）。该地区的排水河道沿东北—西南方向流动。

数据的半可变图形展现出一种特殊模式：短距离（h=125m）的半方差小； 但采样间隔增大到 4 到 5 倍像元大小时（h=500m 到 625m）半方差增大。此后半方差保持常数水平，即使采样间隔增加到 7—10 倍 h，平均值下降不多。

图 4—10 中的两条曲线从实验获取的数据点中通过时，两者已较好地匹配一起。这两条曲线都是数学模型，为了描述半方差随距离 h 而变化的情况， 该数学模型已与实验产生的半方差匹配。首先考虑用实线表示的曲线，它显示出一些重要特性：①延迟（距离）h 的值较大的部分，曲线呈水平走向。曲线的水平部分称为“梁”（sill），它说明在这些延迟范围内数据点没有空间依赖性，因为所有的差分方差不随距离增减而变化。②曲线从 r（h）的低值处开始上升，直升到“梁”的位置为止。达到“梁”时的λ值称为变程。这是半可变图的最重要部分，因为它描述了距离什么值的范围内其内部差分方差为空间非独立方差。变程范围通常用符号α表示，位置越靠近就越相似。加权移动平均内插法就是用变程范围来确定窗口的大小。很显然，数据点和未知点之间的距离大于变程范围时，该数据点在内插未知点的值时将毫无用处，因为它们离得太远。③图中的拟合模型没有通过原点，而是在 r（h）的正方向与 y 轴相截。按（4—20）式 h=0′时，r（h）必须是零。模型中出现的正值是剩余误差ε″的估计值，它是空间上无关噪声。ε″称为“核”方差，它是观测误差的剩余变差与空间变差的组合体。

球面模型常被认为是土壤数据的实验性可变图的最佳描绘形式。其它的模型还有指数模型、线性模型等。

球面模型为：

r（h）=c0+c1{3h/2α-（h/α）³/2} 0＜h＜α r（h）=c0+c1 h＞α（4-21）

r（0）=0 h=0

式中α是变程，h 是延迟，c0 为核方差，而 c0+c1 为梁，在这例子中，球面拟合模型的参数是α=4.598，c0=2.959，c0+c1=41.744。这些参数是由加权最小二乘法计算的。

指数模型为：

r（h）=c0+c1[1-exp（-h/α）]（4－22）

图 4－10 中指数模型的各参数为α=1.942，c0=0.023，c0+c1=42.957。线性模型为：

r（h）=c0+bh（4-23）

式中 b 为线的斜率。当半可变图不出现“梁”的情况下使用线性模型。另外， 当变程的大小远远超过人们希望的内插范围时也用线性模型，统一克里金技术也要用到线性模型。

前面的讨论都假设地表特征的变化在各个方向上都是相同的，然而许多情况下空间变化中的ε″（x）都具有明显的方向性和空间上的依存关系，这时就要用模型不同参数来描述半可变图。

拟合后的半可变图可以用于决定局部内插需要的权因子λ1。决定λ1 的过程与加权移动平均法类似，但不是按任何一种方便的固定空间函数计算λ 1，而是按采样半可变图的地理空间统计分析原理计算。即

^ n

z(x 0 ) = ∑λ i z(x i )

i=1

（4 - 24 ）

其中∑λi = 1。权λ i 的选取应使z（x 0 ）无偏，且估计的方差σ ² 小于

观测值的其它线性组合产生的方差。

下式成立时可获得 z（x0）的极小方差：

n

∑λ _j·r(x_i ，x _j ) + Ψ = r(x _i ，x₀ )

j=1

极小方差为：

i = 1，2， ，m

(4 − 25)

σ ² =



j=1

λj·r(x _j , x ₀ ) + Ψ （4－ 26）

式中 r（xi，xj）是 z 在采样点 xi 和 xj 之间的半方差；r（xi，x0）是样点 xi 和未知点 x0 之间的半方差，这两个量都从已拟合的经验模型产生的半可变图上获取。量Ψ是计算最小方差需要的拉格朗日乘法算子。

如果克里金模型使用以上公式，那么它就是最精确的内插模型，内插值或最佳局部均值与数据点上的值恰好一致。制图中常用规则格网内插未知点的值，因为内插精度比其它格网形式好。内插值可用前面提到过的技术，转

换成等值线地图。与此类似，估算的误差σ ²也可用于制图，误差制图的

e

结果能以数量信息反映整个研究区域内内插值的可靠度。

图 4－12 为考虑各向异性的线性模型对图 4—11 的土壤数据进行内插的结果。图 4－12a 是克里金点模型对 37.5×37.5m 像元进行运算后获得的结果。图 4—12b 是克里金方差的制图结果。

显然，克里金模型可以按更好的办法解决内插中权的估算问题，而且能够提供误差方面的信息。但内插值产生的地图不一定能满足需要，因为克里金点模型或简单模型（4—25 式和 4—26 式）的内插值都与原始样本的面积和容量有关，而土壤采样时，样本经常只有几厘米宽。对该方差起决定作用的土壤变化的大范围和短距离特点，使简单克里金模型产生的地图上，造成了数据点上的凹凸现象。克服这一现象的办法是修改克里金方程，使其能估算子块 B 内的 z 均值，或者只对相对粗略的栅格数组进行内插时，这种方法比较适用。

子块 B 内的 z 均值为：

z（xB）=B×z（x）dx/面积 B

z（xB）由式

^ n

z（x B） = ∑ λi z（xi ） （4 - 27）

i=1

n

来估算。其中∑λi = 1，与前述一样，最小方差变成：

i=1

n

σ² = ∑λ r（x ，x ） + Ψ - r（x ，x ） （4－28）

2 j j B j=1

当

n

B j B

∑λ j r(x i , x j ) + Ψ B = r(xi , xB ) i = 1,2, ，n

j=1

成立时，可得σ² 。

（4 - 29）

由克里金块模型获得的方差估计值，常小于它的点模型估计值，使用这些块模型公式生成的平滑表面就不会发生点模型的凹凸现象。

图 4－13 表示克里金块模型估算的方差结果，与前面的方法一样，使用各向异性线性半可变模型进行内插，得到的平滑曲面反映出明显的长距离渐变特征。

克里金内插模型是一种先进的内插技术，它强烈依赖于统计理论和大量的计算工作，特别是用于制图的克里金内插结果效果更好。克里金模型是理想的内插工具，它在计算过程中妥善选取表明空间相关联的形状和大小，使点或面的局部估计值得以改善。另一个好处是，它能同时产生与内插值有关的误差估计值，这是其它内插法不能做到的。

克里金模型的缺点是：从数据中清除不稳定性，使其达到与固有假设相同的状态是极其困难的研究领域，GIS 的新用户必然会向专家咨询，如何才能达到必要的数据转换。用于矿业、土壤调查等的地理信息系统已开始把最佳内插法综合起来使用。实际上，地矿勘探用的地理信息系统几乎全是地理统计软件包。