ABSTRACT

Spatial data interpolation is a kind of scientific algorithm which is widely used in meteorology, agriculture, geological exploration and some other fields. People can get a lot of data which can meet the needs of research and production with a few investments by using the spatial data interpolation algorithm to estimate the attribute value of unknown point. Therefore, spatial data interpolation can meet the needs of human to a great extent, which has great practical significance. In this paper, we will introduce two classical spatial data interpolation algorithm—the Inverse Distance Weighted interpolation algorithm (IDW) and the Ordinary Kriging interpolation algorithm (OK). Furthermore, on the grounds of researching in the Inverse Distance Weighted interpolation algorithm, an improved Inverse Distance Weighted interpolation algorithm is proposed for the distance power exponent k, which is an optimizable parameter. By applying the spatial data interpolation algorithms which are introduced in this paper into air quality analysis of Jinan city, we can derive the following conclusions: (1) the Ordinary Kriging interpolation algorithm is better than the Inverse Distance Weighted interpolation algorithm; (2) the exponential model has the best interpolation effect in the Ordinary Kriging interpolation algorithm; (3) the improved Inverse Distance Weighted interpolation algorithm is better than the traditional Inverse Distance Weighted interpolation algorithm.

Keywords:Spatial Interpolation, Inverse Distance Weighted Interpolation, Kriging Interpolation, Variogram

空间数据插值方法研究

张婕

天津师范大学数学科学学院，天津

收稿日期：2019年11月6日；录用日期：2019年11月21日；发布日期：2019年11月28日

摘 要

空间数据插值算法是目前普遍应用的一种科学算法，其广泛应用于气象、农业及地质勘探等领域。通过空间数据插值，对缺省或非有效空间数据进行估计和推测，在投资较少的情况下得到大量精度能够满足研究及生产需要的数据，可以在很大程度上满足人类的需要，具有很强的现实意义。本文详细介绍了两种常用的空间数据插值算法——反距离权重插值算法(IDW)和普通克里金插值算法(OK)的原理，并且基于对反距离权重插值算法的研究，针对其可优化参数——距离幂指数k，提出了一种改进的反距离权重插值法。通过数值模拟对济南市的空气质量状况进行空间插值分析，可以得出结论：(1) 普通克里金插值算法优于反距离权重插值算法；(2) 在普通克里金插值算法中，指数模型的插值效果较好；(3) 本文提出的改进的反距离权重插值算法较原算法的插值效果更佳。

关键词 :空间插值，反距离权重插值，克里金插值，变异函数

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

随着空间观测技术的快速发展，传递给人类的信息越来越多，人们对于空间地质信息预测的精度要求不断提高，也越来越意识到空间数据的采集、处理与分析对于日常生产生活所起到的重大作用。由于多种原因，有时无法通过观测获取完整有效的空间数据。例如在获取遥感影像数据时受到云层或阴影的影响，在地质资源探测时受到特殊地形地势的影响等都可能导致所获得的数据不够完整。如何利用非有效空间数据成为一个必须要解决的难题。如果直接将无效数据过滤掉，势必会影响数据的完整性，因此应尽可能地使用临近空间的已知采样点数据对目标空间的相关数据进行推测。这是解决缺省或者无效空间数据的非常有效的手段之一，可以在很大程度上满足人们的需要，即空间数据插值。通过空间数据插值方法，可以构造出更加贴近实际的模型，从而在投资较少的情况下得到大量精度能够满足研究及生产需要的数据，具有很强的现实意义。

空间数据插值算法有很多的种类，也具有相对成熟的理论基础，对很多空间数据都具有很好的模拟构造功能。其主要目标是对缺省或非有效的数据的估计和推测，从而实现数据的网格化，并显示数据的空间分布特征 [1]。虽然数据插值一般要以数据质量的损失为代价，但是国内外学者不断地改进各种方法，完成了很多成效甚佳的工作，这也表明空间数据插值算法切实可行。胡庆芳、胡艳等 [2] 将普通克里金插值应用于降水研究，开发了基于MATLAB平台的面向大范围降水空间插值的普通克里金模型。曾晖 [3] 利用泛克里金插值法对南京和江门两市的住宅小区均价进行插值分析，揭示了两市住宅小区均价的时空分布规律。刘敏 [4] 将空间数据插值算法应用于和田绿洲地下水时空分布规律的研究，为和田地区农业灌溉和地下水的开发利用提供了科学依据。吴健生，王仰麟等 [5] 利用空间数据插值算法对矿体内部的特性数据进行插值估算，从中得到矿体的内部特征规律，指导矿山的生产。虽然众多学者在空间数据插值算法方面做了大量的统计分析和改进工作，但实际应用往往不能与理论模型完全相符，目前的技术方法也并非能完全满足人们的实际应用。因此在选择插值方法时，需要联系实际需要对算法进行分析、研究、改进，然后应用 [6]。

本文主要介绍了两种常用的空间数据插值算法——反距离权重插值法(IDW)和普通克里金(OK)插值法，将上述两种算法应用于济南市的空气质量状况研究，并且提出了一种改进的反距离权重插值算法。本文的第二部分将对IDW插值法及OK插值法进行介绍；在第三部分，笔者基于对IDW插值法的研究，提出一种改进的IDW插值法；第四部分将对数值试验的结果进行展示；第五部分对本文的主要内容进行总结。

2. 反距离权重插值法和普通克里金插值法

本章主要介绍了反距离权重(IDW)插值算法和普通克里金(OK)插值法，作为经典的两种插值方法，它们被广泛应用于遥感影像、地质勘探等众多领域。

2.1. 反距离权重插值法

反距离权重(IDW)插值算法是空间数据插值领域应用最广的一种算法。它属于距离权重系数方法的一种，是以相近相似原理为理论基础，通过控制权重的大小来决定已知点的值对待估点的值的影响大小，给予距离远的点的权重小于距离近的点的权重 [7] [8]，即待估点的值是根据待估点周围的样本点属性值随着其到样本点距离的变化而发生反相关变化。随着待估点与某样本点的接近，该样本点所占权重趋近于一，故待估点的属性值和样本点的属性值随着距离的接近而接近。

表示如下：

$z^{*}\left(x_0\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\omega }_{i}z\left(x_i\right)}$ (2.1.1)

其中，x₀为待估值点，z^*(x₀)为待估值，x_i为区域内已知样本点， ${\omega }_i$ 为该点对应的权重系数，z(x_i)为x_i点的观测值。

权重系数 ${\omega }_i$ 一般由下式给出：

${\omega }_i=f\left(D_i\right)/{\displaystyle \underset{i=1}{\overset{n}{\sum }}f\left(D_i\right)}$ (2.1.2)

其中，n为参与插值的样本点个数， $f\left(D_i\right)$ 为待估值点与已知点之间的距离D_i的权重系数，最常用的 $f\left(D_i\right)$ 的形式为：

$f\left(D_i\right)=\frac{1}{D_i^k}$ (2.1.3)

其中，D_i是待估值点与已知点间的距离，一般取欧氏距离，即

$D_i=\sqrt{{\left(x_0-x_i\right)}^2+{\left(y_0-y_i\right)}^2}$ (2.1.4)

k > 0，为距离的幂指数，显著地影响插值的结果。IDW插值法存在一个特例，当IDW插值法的幂指数k取值比较大时(如k > 10)，其插值的结果与最近距离插值法的估值结果更接近 [9]。实际应用中，一般取k = 1或k = 2进行插值，反距离平方加权法就是k的取值为2的IDW插值法。

结合式(2.1.1)~(2.1.3)，IDW插值法可由下式表示：

$z^{*}\left(x_0\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{1}{{\left(D_i\right)}^k}z\left(x_i\right)}/{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{1}{{\left(D_i\right)}^k}}$ (2.1.5)

可见，运用IDW插值法进行插值估计的难点在于对所研究的空间范围(或空间半径)以及距离的幂指数k等多个参数进行合理选择，这也是在利用该方法解决实际问题时的关键。简单、易实现、效率高、运行时所占存储空间小是IDW插值法的优点。但该法的缺点是没有考虑点的构性，而只考虑了距离因素，所得插值结果的误差可能较大。

2.2. 克里金插值原理

克里金插值算法是由南非矿产学家RG. Krige首先提出的一种插值算法，并且以此而命名。其理论基础为变异函数理论分析，它是对有限区域内的区域化变量的属性值进行无偏最优估计的一种方法。概率论与数理统计中所涉及的无偏、最小方差条件是克里金插值法的基本原理。此外，它还是一种线性插值算法，是改进距离权重插值算法的一种方法。该法考虑了空间数据插值中的多种因素：(1) 已知点的空间分布结构特点；(2) 邻近的已知点与未知点间的空间位置关系；(3) 各邻近点之间的位置关系，而且对空间结构进行了研究，从而使插值算法所得结果更精确，通过插值所建立的模型更加符合实际。

普通克里金(OK)插值法是使用最广泛的一种最优无偏估计方法，其适用的条件是：所研究的区域化变量的数学期望未知，但是一定为一常数m。它的计算公式与反距离权重插值类似，具体为：

$Z^{*}\left(x_0\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_{i}Z\left(x_i\right)}$ (2.2.1)

其中， $Z\left(x_i\right)\left(i=1,2,\cdots ,n\right)$ 为数据点x_i的值，x₀为所求的未知点， ${\lambda }_i$ 为权重系数。此处的 ${\lambda }_i$ 不是单一地由距离来决定，它是在无偏最小方差的条件下，依赖于变异函数 [10] 的计算而确定的，这也是它与IDW插值的不同之处，也是它的优势的体现。

为了使Z^*(x₀)为Z(x₀)的无偏估计值，需要

$E\left[Z^{*}\left(x_0\right)-Z\left(x_0\right)\right]=0$，

而 $E\left[Z\left(x_0\right)\right]=m$，

又 $E\left[Z^{*}\left(x_0\right)\right]=E\left[{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_{i}Z\left(x_i\right)}\right]={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_{i}E\left[Z\left(x_i\right)\right]}=m{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_i}$，

故： ${\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_i}=1$。 (2.2.2)

其估计方差的计算公式为：

$E\left\{{\left[Z\left(x_0\right)-Z^{*}\left(x_0\right)\right]}^2\right\}={\displaystyle {\sum }_{i=0}^n{\displaystyle {\sum }_{j=0}^n{\alpha }_i{\alpha }_{j}Cov\left(x_i,x_j\right)}}$，

如果令 ${\alpha }_0=1$，${\alpha }_i=-{\lambda }_i\left(i=1,2,\cdots ,n\right)$，则

$\begin{array}{l}E\left\{{\left[Z\left(x_0\right)-Z^{*}\left(x_0\right)\right]}^2\right\}={\displaystyle \underset{i=0}{\overset{n}{\sum }}{\displaystyle \underset{j=0}{\overset{n}{\sum }}{\alpha }_i{\alpha }_{j}Cov\left(x_i,x_j\right)}}\\ =Cov\left(x_0,x_0\right)-2{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\lambda }_{j}Cov\left(x_0,x_j\right)}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\lambda }_i{\lambda }_{j}Cov\left(x_i,x_j\right)}}\end{array}$ (2.2.3)

则需要在(2.2.2)式的约束下，求出使得估计方差为最小的 ${\lambda }_i$，由(2.2.3)式，利用拉格朗日乘数法构造函数：

$F=E\left\{{\left[Z\left(x_0\right)-Z^{*}\left(x_0\right)\right]}^2\right\}-2\mu \left({\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_i}-1\right)$ (2.2.4)

其中， $\mu >0$ 为拉格朗日乘子。如果有一组 ${\lambda }_i\left(i=1,2,\cdots ,n\right)$ 满足上述要求，则这组数值可以使得方程组(2.2.4)满足：

$\{\begin{array}{l}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_i}=1\hfill \\ {\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_{i}Cov\left(x_i,x_j\right)}-\mu =Cov\left(x_0,x_j\right)\left(j=1,\cdots ,n\right)\hfill \end{array}$ (2.2.5)

其中

这n + 1个方程组就是普通克里金方程组，由此可得其矩阵形式为：

$\left[K\right]\left[\lambda \right]=\left[M\right]$ (2.2.6)

其中，

$\left[\lambda \right]=\left[\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ ⋮\\ {\lambda }_n\\ -\mu \end{array}\right]$，$\left[M\right]=\left[\begin{array}{c}{C}_{01}\\ C_{02}\\ ⋮\\ C_{0n}\\ 1\end{array}\right]$，$\left[K\right]=\left[\begin{array}{ccccc}{C}_{11}& C_{12}& \cdots & C_{1n}& 1\\ C_{21}& C_{22}& \cdots & C_{2n}& 1\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ C_{n1}& C_{n2}& \cdots & C_{nn}& 1\\ 1& 1& \cdots & 1& 0\end{array}\right]$，

$C_{ij}$ 表示 $Cov\left(x_i,x_j\right)$。

进行普通克里金法估计时所得的估计方差被称为普通克里金方差，其表达式为：

$\begin{array}{c}{\sigma }_{OK}^2=Var{\left[Z^{*}\left(x_0\right)-Z\left(x_0\right)\right]}^2\\ =Cov\left(x_0,x_0\right)-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_{i}Cov\left(x_0,x_i\right)}+\mu \end{array}$ (3.2.7)

变异函数 $\gamma \left(x_i,x_j\right)$ 定义为：

$\gamma \left(x_i,x_j\right)=\gamma \left(x_i-x_j\right)=\frac{1}{2}E{\left[Z\left(x_i\right)-Z\left(x_j\right)\right]}^2$ (3.2.8)

若用变异函数 $\gamma \left(x_i,x_j\right)$ 来表示协方差 $Cov\left(x_i,x_j\right)$，则由方程组(3.2.5)得：

$\{\begin{array}{l}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_i}=1\hfill \\ {\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_i\gamma \left(x_i,x_j\right)}-\mu =\gamma \left(x_0,x_j\right)\left(j=1,\cdots ,n\right)\hfill \end{array}$ (3.2.9)

它对应的普通克里金方差为：

${\sigma }_{OK}^2={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_i\gamma \left(x_0,x_i\right)}-\gamma \left(x_0,x_0\right)+\mu$ (3.2.10)

在实现OK插值算法时，首先对获取的样本数据进行预处理，然后根据处理后的数据计算实验的变异函数，选择最佳的变异函数理论模型来进行参数拟合，确定变异函数中的相关参数，之后使用普通克里金方程进行计算，最后验证插值结果是否合理，若结果合理，则显示插值结果。

3. 改进的反距离权重插值法

反距离权重插值法的公式为：

$z^{*}\left(x_0\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{1}{{\left(D_i\right)}^k}z\left(x_i\right)}/{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{1}{{\left(D_i\right)}^k}}$ (3.1.1)

其中，距离幂指数k为可优化参数，针对不同的地点，距离幂指数k会有所不同，故改进的反距离权重插值法以研究区域之前已知的数据建立数据库，通过遍历指定区间内的k值，利用统计学的方法确定该区域上的最佳k值，应用于插值计算，从而使插值结果更加精确 [11]。

3.1. 统计分析方法

采用交叉验证法(Cross Validation)来对不同的k值对应的插值结果进行验证比较。即首先假定所研究区域内部分点的属性值未知，使用周围点的属性值来进行估算，然后计算实际观测值与估测值的误差，以此来找出该研究区域所对应的最佳k值。采用平均误差(ME)、平均绝对误差(MAE)、插值平均误差平方和的平方根(RMSIE)作为分析比较的标准，从估计误差的大小、估计值可能的误差范围、利用样本点数据的估值灵敏度和极值效应等方面对不同k值所对应的插值结果进行验证比较 [12]，从而找出最佳k值。

3.2. 实现流程

具体实现流程如下：

4. 数值实验

空间数据插值技术在地质、医学、气象、空气质量监测等不同领域均有很大应用，本章主要探索空间数据插值在空气质量监测领域的应用。在本次数值实验中，我们选取的研究区域为山东省的济南市、泰安市、德州市等17个地市，利用泰安市、德州市等16个地市自2016年每月1号的空气质量指数对济南市2016年每月1号的空气质量进行空间插值分析。本文所研究数据的来源为中华人民共和国环保部数据中心 [13]。

4.1. 反距离权重插值法和普通克里金插值法

4.1.1. 反距离权重插值法的数值实验结果

用反距离权重插值法对济南市2016年每月1号的空气质量进行插值，为使计算简便，取反距离权重插值法中的幂指数k为2进行计算，计算结果见表1。

4.1.2. 普通克里金插值法的数值实验结果

本文分别用普通克里金插值法的球状模型、高斯模型和指数模型为变异函数模型对济南市2016年每月1号的空气质量进行插值 [14] [15]，计算结果见表2、表3、表4。

4.1.3. 反距离权重插值与普通克里金插值的比较

通过对表1、表2、表3、表4的比较，可以看出，在本次实例研究中，除了2016年1月1日外，普通克里金插值均优于反距离权重插值，通过进一步的分析比较，可以看出，在普通克里金插值中，指数模型的插值准确率最高，效果是最优的。而对于高斯模型和球状模型，球状模型插值的绝对误差大大低于高斯模型，故球状模型的插值效果优于高斯模型。为了更加直观地比较插值的效果，图1展示了反距离权重插值与三种普通克里金插值结果的比较图。

表1. 反距离权重插值法(IDW)

表2. 普通克里金插值(球状模型)

表3. 普通克里金插值(高斯模型)

表4. 普通克里金插值(指数模型)

图1. 插值结果比较图

4.2. 改进的反距离权重插值法

将改进的反距离权重插值法应用于4.1所研究的实例，以所研究区域2013年至2015年每一天的AQI值为原始数据库，设定k的取值范围为 [1] [4]，步长为0.1，并将该方法的结果与反距离权重插值方法结果进行比较，比较结果见表4和表5。

表5. 改进的反距离权重插值

将反距离权重插值的结果与上述结果进行比较，比较结果见图2。

由此，我们可以看出经过改进的反距离权重插值法的插值结果优于反距离权重插值法。

图2. 反距离权重插值与改进的反距离权重插值比较

5. 总结

本文详细介绍了反距离权重插值法(IDW)与普通克里金插值算法(OK)的原理和实现，并将这两种方法应用于济南南市2016年每月1日的空气质量状况研究。经过实验，我们得出：普通克里金插值的插值效果优于反距离权重插值；在普通克里金插值中，指数模型的插值准确率最高，效果最优，球状模型插值的绝对误差大大低于高斯模型，故其插值效果优于高斯模型。在对反距离权重插值算法研究的基础上，针对其可优化参数——距离幂指数k，提出了一种改进的反距离权重插值法。通过找出所研究区域上的最佳k值，将其应用于插值计算，以提高插值效果。利用该插值算法对济南市的空气质量状况进行研究，并与反距离插值分析的结果进行比较，可以看出，经过改进的反距离权重插值算法的插值效果明显优于反距离权重插值算法，具有很强的可行性。

文章引用

张 婕. 空间数据插值方法研究
Research on Spatial Data Interpolation[J]. 应用数学进展, 2019, 08(11): 1859-1869. https://doi.org/10.12677/AAM.2019.811216

参考文献