两种空间插值方法的比较研究

两种空间插值方法的比较研究摘要：距离倒数加权法算法简单，容易实现，适合分布较均匀的采样点集，但容易出现“牛眼”现象；克里金法是一种无偏最优估计法，精度较高，适合空间自相关程度高的数据，但其算法复杂，实现较难。

这两种方法各有其适用情形，本文比较了这两种方法的优劣并提出算法优化的思路。

关键字：距离倒数加权，克里金，优化1引言空间插值是根据一组已知的离散数据或分区数据，按照某种假设推求出其他未知点或未知区域的数据的过程，简单的说就是由已知空间特性推求未知空间特性。

它是地学研究中的基本问题，也是GIS 数据处理的重要内容。

在利用GIS 处理空间数据的过程中，需要进行空间插值的场合很多，如采样密度不够、采样分布不合理、采样存在空白区、等值线的自动绘制、数字高程模型的建立、区域边界分析、曲线光滑处理、空间趋势预测、采样结果的2.5维可视化等[1]。

通过归纳，空间插值可以简化为以下三种情形：（1）现有离散曲面的分辨率、像元大小或方向与所要求的不符，需要重新插值。

例如将一个扫描影像（航空像片、遥感影像）从一种分辨率或方向转换为另一种分辨率或方向的影像。

（2）现有连续曲面的数据模型与所需的数据模型不符，需要重新插值。

如将一个连续曲面从一种空间切分方式变为另一种空间切分方式，从TIN 到栅格、栅格到TIN 或矢量多边形到栅格。

（3）现有数据不能完全覆盖所要求的区域范围，需要插值。

如将离散的采样点数据内插为连续的数据表面[2]。

现有的空间插值方法多种多样，但每一种方法都有其适用情形和无法避免的缺陷，本文分析了距离倒数加权法和克里金法的插值结果，并提出改进的思路。

2方法距离倒数加权法和克里金法都是建立在地理学第一定律之上的，即：空间距离越近，地理事物的相似性越大[3]。

它们都是通过确定待插点周围采样点的权重来求取待插点的估计值，可统一表示。

设n x x ,,1 为区域上的一系列观测点，)(,),(1n x Z x Z 为相应的观测值。

待插点0x 处的值)(0x Z 可采用一个线性组合来估计：∑==ni i i x Z x Z 10)()(λ （1）但距离倒数加权法只考虑采样点与待插点之间的距离，而克里金法不仅考虑距离，还要考虑采样点的空间分布及其与待插点的空间方位关系[4]。

2.1距离倒数加权法距离倒数加权权重i λ的赋值表达式为[][]0;,,2,1),(),(100>==∑=--αλααm i x x d x x d m i i i i （2）式中，幂指数α越小，权重越趋向取平均值；α越大，越近的点权重越大，越远的点权重越小[5]。

当α为零时，就是等权模型，即m i /1=λ，等权虽然简单易操作，但忽略了地理学第一定律。

有的文献[6]采用下列方案确定权值i λ[]⎩⎨⎧==else x x d x x d x x d x x d n i i ,0),(),...,,(),,(min ),(,121λ （3）它相当于待插点取最邻近点的值，即泰森多边形法（最近邻点插值法）。

2.2克里金法地质统计学是以区域化变量为基础，借助变异函数，研究既具有随机性又具有结构性，或空间相关性和依赖性的自然现象的一门科学[4]。

克里金插值法是地质统计学的重要组成部分，也是地质统计学的核心。

2.2.1区域化变量能用空间分布来表征一个自然现象的变量称为区域化变量，它反映了区域内的某种特征或现象。

区域化变量根据区域内位置的不同而取不同的值，可以说，它是与位置有关的随机变量。

在进行采样观测以后，可以将其表示为一个空间点函数),,()(w v u x x x Z x Z = （4）式中， w v u x x x ,, 为三维直角坐标系中的三轴。

区域化变量具有以下几个特性：随机性：区域化变量是一个随机变量，它具有局部的、随机的、异常的特征；结构性：区域化变量在一定范围内具有某种程度的相似性，即自相关性，当超出这一范围时，自相关性消失；空间局限性：即这种结构性的表现被限定在一定的空间内；连续性：不同的区域化变量具有不同程度的连续性，其连续性是有变异函数来表示的；异向性：区域化变量可能表现为各向同性，也可能表现为各向异性[4,7]。

2.2.2平稳假设克里金插值法是一种无偏最优估计法[8]，其中，无偏是指偏差为0，即要服从二阶平稳或本证假设。

1.二阶平稳当区域化变量)(x Z 满足下列两个条件时，称其为二阶平稳或弱平稳，在整个研究区内有)(x Z 的数学期望存在，且等于常数，即：m h x Z E x Z E =+=)]([)]([ （5）在整个研究区内，)(x Z 的协方差函数存在且平稳，即只依赖于滞后h ，而与x 无关：)()]()([)]([)]([)]()([)}(),({2h C m h x Z x Z E h x Z E x Z E h x Z x Z E h x Z x Z Cov =-+=+-+=+ （6）特殊的，当0=h 时，上式变为)0()]([C x Z Var =，即方差存在且为常数。

2.本征假设是比二阶平稳更弱的平稳假设，当区域化变量)(x Z 的增量)()(h x Z x Z +-满足下列两条件时，称其为满足本征假设或内蕴假设。

在整个研究区内有0)]()([=+-h x Z x Z E （7）增量)()(h x Z x Z +-的方差函数存在且平稳（即不依赖于x ）：)(2),(2)]()([)]}()([{)]()([)]()([222h h x h x Z x Z E h x Z x Z E h x Z x Z E h x Z x Z Var γγ==+-=+--+-=+- （8）2.2.3变异函数变异函数是地统计学特有的基本工具，它既能描述区域化变量的空间结构性变化,又能描述其随机性变化。

区域化变量)(x Z 在点x 和h x +处的值)(x Z 与)(h x Z +差的方差的一半称为区域化变量)(x Z 的变异函数，记为),(h x γ。

在二阶平稳假设或本征假设的条件下，有222)]()([21)]}([)([{)]()([21)]()([21),(h x Z x Z E h x Z E x Z E h x Z x Z E h x Z x Z Var h x +-=+--+-=+-=γ （9） 由上式可知，变异函数依赖于自变量x 和h ，当变异函数),(h x γ仅仅依赖于距离h 而与位置x 无关时，),(h x γ可改写为)(h γ，即2)]()([21)(h x Z x Z E h +-=γ （10） 具体表示为∑=+-=)(12)]()([)(21)(h N i i i h x Z x Z h N h γ （11） 在变异函数中有四个参数，它们的定义如下。

变程（range ）：指区域化变量在空间上具有相关性的范围。

在变程范围之内，数据具有相关性；而在变程之外，数据之间互不相关，即在变程以外的观测值不对估计结果产生影响。

块金值(nugget)：变异函数如果在原点间断，在地质统计学中称为“块金效应”，表现为在很短的距离内有较大的空间变异性，无论h 多小，两个随机变量都不相关。

测量误差与自然现象的微观变异性任意一方或两者共同作用产生了块金值。

基台值(sill)：代表变量在空间上的总变异性大小。

即为变差函数在h 大于变程时的值，为块金值和拱高之和。

拱高(partial sill)：基台值与块金值之差。

当块金值等于0时，基台值即为拱高。

当变异函数确定之后，执行克里金系统就只是一个简单的计算过程。

变异函数的理论模型主要有以下三个[9,10]。

球状模型()⎪⎪⎩⎪⎪⎨⎧≥≤⎥⎥⎦⎤⎢⎢⎣⎡⎪⎭⎫ ⎝⎛-⋅==⎪⎭⎫ ⎝⎛⋅=ah c a h a h a h c h a h Sph c h ,,2123003γ （12） 指数模型 ()⎥⎦⎤⎢⎣⎡⎪⎭⎫ ⎝⎛--⋅=⎪⎭⎫ ⎝⎛⋅=a h c a h Exp c h 3exp 1γ （13） 高斯模型 ()()⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛--⋅=223ex p 1a h c h γ （14） 式中，c 为基台值，a 为变程，h 为滞后距变异函数的形式是内插质量的关键，需要注意的是，由于不同的区域有不同的空间模式，因而也就有不同的变异函数，而空间内插都有一个隐含的假定，即空间是连续的，因此，在选择变异函数模型之前，检查数据以确定空间连续性是十分必要的[11]。

2.2.4普通克里金依据本征假设和无偏性要求，有∑==n i i 11λ（15）在以上的条件下要达到估计方差最小[10]，即()()()()()(){}[]()()()[]min 200*200*00*2=-=---=x Z x Z E x Z x Z E x Z x Z E k σ （16）再应用拉格朗日乘数法求极值，进一步推导可得()()()⎪⎪⎩⎪⎪⎨⎧==-=--∑∑==n i i j ni i j i n j x x C x x C 1011,,1λμλ （17）最小的估计方差,即克里金方差可用以下公式求解 ()()∑=--+-=-=n i i i k x x C x x C x Z x Z Var 100000*2)]()([λμσ （18）或()()00102x x x x ni i i k --+-=∑=γμγλσ （19） 3实验3.1实验数据福建地处我国东南部，山地、丘陵占陆域的80%，属温暖湿润的亚热带海洋性季风气候[12]，其地形复杂，降雨丰富，本文选取福建省为空间插值的实验区，收集该省67个气象站点的历史数据，以2005年12月月平均降水为试验指标，进行内插分析比较研究。

3.2距离倒数加权法应用ArcGIS9.2 Geostatistcal Analyst 功能模块进行插值，搜索半径取1.2，幂指数取2， 从CrossValidation 的统计来看，最大误差高达-200.667774109（武夷山），最小误差只有-1.14440513614（龙海），出现这一结果的原因是在武夷山站周围气象观测站点分布较稀疏，可选点较少；而在龙海站周围的分布较均匀，可选点较多。

在结果图中“牛眼”现象比较明显的是九仙山站（269）（见图1），与它相邻的站点分别是永春站（120）、大田站（182）、尤溪站（182）、永泰站（124），产生这一现象的原因是距离倒数加权法在站点数据较少时容易受采样极值点的影响，会在网络区域产生围绕采样点位置的“靶心”，即所谓的“牛眼”现象[13]。

图1 距离倒数加权法，九仙山站及其周围四县市的插值结果图3.3克里金法在克里金精度评价系统中，符合下列标准的结果是最优的：标准平均值（Mean Standardized）最接近于0，均方根预测误差（Root-Mean-Square）最小，平均标准误差（Average Mean Error）最接近于均方根预测误差（Root-Mean-Square），标准均方根预测误差（Root-Mean-Square Standardized）最接近于1[4]。