克里金插值-Kriging插值-空间统计-空间分析

克里金插值方法-Kriging 插值-空间统计-空间分析
1.1 Kriging 插值
克里金插值（Kriging 插值）又称为地统计学，是以空间自相关为前提，以区域化变量理论为基础，以变异函数为主要工具的一种空间插值方法。

克里金插值的实质是利用区域化变量的原始数据和变异函数的结构特点，对未采样点的区域化变量的取值进行线性无偏、最优估计。

克里金插值包括普通克里金插值、泛克里金插值、指示克里金插值、简单克里金插值、协同克里金插值等，其中普通克里金插值是最为常用的克里金插值方法。

以下介绍普通克里金插值的原理。

包括普通克里金方法在内的各种克里金插值方法的使用前提是空间数据存在着显著的空间相关性。

判断数据空间相关性是否显著的工具是半变异函数（semi-variogram ），该函数以任意两个样本点之间的距离h 为自变量，在h 给定的条件下，其函数值估计方法如下：
2||||1()[()()]2()i j i j s s h h z s z s N h γ-==-∑
其中()N h 是距离为h 的样本点对的个数。

()h γ最大值与最小值的差m a x m i n γγ-可以度量空间相关性的强度。

max min γγ-越大，空间相关性越强。

如果()h γ是常数，即max min 0γγ-=，则说明无论样本点之间的距离是多少，样本点之间的差异不变，也就是说样本点上的值与其周围样本点的值无关。

在实际操作中，会取一些离散的h 值，当
||s s ||i j -接近某个h 时，即视为||||i j s s h -=。

然后会通过这些离散点拟合成连续的半变异函数。

拟合函数的形式有球状、指数、高斯等。

在数据存在显著的空间相关性的前提下，可以采用普通克里金方法估计未知点上的值。

普通克里金方法的基本公式如下：
01ˆ()()()n i i
i Z s w s Z s ==∑
普通克里金方法的基本思想是：通过调整i s 的权重()i w s ，使未知点的估计
值0ˆ()Z s 满足两个要求：1.
0ˆ()Z s 是无偏估计，即估计误差的期望值为0，2.估计误差的方差达到最小。

上述问题实际是一个条件极值问题，可以用拉格朗日乘子法求解。

构造出拉格朗日方程组如下：
11()(||s s ||)(||s s ||)()1n
i i j i j j n i i w s w s γμγ==⎧-+=-⎪⎪⎨⎪=⎪⎩∑∑
其中，1,,i n =，μ和()i w s 为方程的未知数。

上述方程组共有1n +个方程
和1n +个未知数，通常存在唯一解。

通过方程组的解即可计算0ˆ()Z s 。

【案例6-3】
模拟已知2008年北京某些监测点PM10的浓度，估计未知点的PM10浓度，并将北京的PM10浓度的空间分布进行可视化。

本例的输入数据包括已知采样点的经纬度坐标和PM10的浓度（C:\Example\Data\6.2Kriging\BeijingPM10.xls ）、未采样点的经纬度坐标（C:\Example\Data\6.2Kriging\forecastpoint.xls ）；为了后期进行可视化，另需北京各区县矢量地图（C:\Example\Data\Geodata \BeijingCounty.shp ）剪切而出。

现将部分数据展示如下：
表 错误！文档中没有指定样式的文字。

-1已知采样点部分数据
X
Y PM10_2008116.4141
39.9885131115.9568
40.351095116.2131
40.2007113116.1854
40.278296116.4142
39.9313132116.2842
39.8574138116.1708
39.9142135116.3547
39.9271118116.1970
40.0780130116.642740.3813100
表 错误！文档中没有指定样式的文字。

-2未采样点数据
X Y
116.629539.93803
116.224139.76457
116.455539.70717
116.326439.99417
116.459140.26683
步骤1：展示北京各区县地图
将北京市各区县地图（C:\Example\Data\Geodata \BeijingCounty.shp）从目录中拖至地图编辑窗口，结果如图错误！文档中没有指定样式的文字。

-1。

图错误！文档中没有指定样式的文字。

-1北京市各区县矢量地图
步骤2：生成已知采样点的矢量地图
点击【文件】-【添加数据】-【添加XY数据】，在对话框内容选择如图错误！文档中没有指定样式的文字。

-2，选择C:\Example\Data\6.2Kriging\BeijingPM10.xls 中的BeijingPM10$，此为已知采样点的数据表；X字段选X（经度），Y字段选Y（纬度）；接着设置地理坐标系，点击【编辑】-【地理坐标系】-【world】-【WGS1984】。

图错误！文档中没有指定样式的文字。

-2添加XY数据的对话框此时地图窗口会呈现新添加的点，但到这步还未生成矢量地图。

接下来右击图层中的【BeijingPM10$ 个事件】-【数据】-【导出数据】，输出要素类为C:\Example\Train\6.2Kriging\BeijingPM10.shp，点击【确定】。

在新弹出的对话框会询问是否添加至图层，点击【是】，新生成的BeijingPM10.shp即在地图窗口展现。

现可右击图层中的【BeijingPM10$个事件】将其【移除】，结果如图错误！文档中没有指定样式的文字。

-3。

图错误！文档中没有指定样式的文字。

-3北京市各区县已知采样点的矢量地图步骤3：克里金插值（图错误！文档中没有指定样式的文字。

-4）
打开【地统计向导】按钮，对话框如图错误！文档中没有指定样式的文字。

-4a，点击【克里金法/协同克里金法】，数据集的【源数据集】选择BeijingPM10，【数据字段】选择PM10_2008，【下一步】。

选择【普通克里金法】，输出表面类型选【预测】-【下一步】（图错误！文档中没有指定样式的文字。

-4b）。

图错误！文档中没有指定样式的文字。

-4c显示有半变异函数图，若图中的曲线越倾向于水平，说明空间相关性越弱，点击【下一步】。

图错误！文档中没有指定样式的文字。

-4d，可在X和Y中输入经纬度坐标，即可显示预测值，或者在左图中点击也可显示。

此时内容列表生成名为【克里金法】的图层，如图错误！文档中没有指定样式的文字。

-4f。

a
b
c
d
e
f
图错误！文档中没有指定样式的文字。

-4克里金插值
步骤4：未采样点的PM10浓度预测
生成未采样点的矢量地图：导入未采样点数据（C:\Example\Data\6.2Kriging\forecastpoint.xls）,生成预测点的矢量地图forecastpoints.shp，方法同步骤2。

预测（图错误！文档中没有指定样式的文字。

-5）：右击图层中的【克里金法】-【验证\预测】，在弹出的对话框中，输入观测点选择刚生成的forecastpoints，点击【确定】，输出点位置处的统计数据可选择默认。

点击【确定】，此时会新生成图层GALayerToPoints，右击该图层，点击【打开属性表】，已预测字段即为未采样点的预测值。

a
b
c
图错误！文档中没有指定样式的文字。

-5未采样点的PM10浓度预测
步骤5：插值结果的地图可视化
设置覆盖范围（图错误！文档中没有指定样式的文字。

-6）：右击图层中【克里金法】，点击【属性】，【范围】选择【矩形范围BeijingCounty】，使其覆盖北京所有区县，结果如图错误！文档中没有指定样式的文字。

-6c。

a
b
c
图错误！文档中没有指定样式的文字。

-6设置覆盖范围导出至删格（图错误！文档中没有指定样式的文字。

-7）：右击图层中【克里金法】，点击【数据】-【导出至栅格】，弹出GA图层至网格对话框，输出表面删格可设为C:\Example\Train\6.2Kriging\GALayerToGri。

点击【确定】后，结果如图错误！文档中没有指定样式的文字。

-7c，
a
b
c
图错误！文档中没有指定样式的文字。

-7导出至删格按掩膜提取（图错误！文档中没有指定样式的文字。

-8）：点击【地理处理】-【ArcToolbox】，打开ArcToolbox 工具箱，依次选择【Spatial Analysis Tools】-【提取分析】-【按掩膜提取】，输入删格选择GALayerToGri，输入删格数据或要素掩膜数据选择BeijingCounty，输出设置可为C:\Example\Train\6.2Kriging\Extract_GAla。

点击【确定】，图层中又出现名为Extract_GAla的图层。

内容列表的图层中，只打勾已知采样点BeijingPM10、未采样点forecastpoints、Extract_GAla、BeijingCounty的图层，且Extract_GAla拖放至BeijingCounty上面，此时PM10的浓度显示范围和北京市一样（图错误！文档中没有指定样式的文字。

-8d）。

a
b
c
d
图错误！文档中没有指定样式的文字。

-8按掩膜提取修改属性（图错误！文档中没有指定样式的文字。

-9）：这一步是修改图层Extract_GAla的属性。

右击图层中的【Extract_GAla】-【属性】；在图层属性对话框的【符号系统】卡中调节色带，本例选择由绿到红的色带，表示浓度由低到高的，高浓度的用红色警示；在【显示】卡中，将透明度调至50%，点击【确定】，
效果如图错误！文档中没有指定样式的文字。

-9d。

a
b
c
d
图错误！文档中没有指定样式的文字。

-9修改属性最后可参照3.1专题地图的内容，添加地图要素并导出，结果如图错误！文档中没有指定样式的文字。

-10所示。

根据克里金的插值结果，大体可知2008年北京南部的PM10浓度高于北部。

图 错误！文档中没有指定样式的文字。

-10 PM10浓度克里金插值估计的整体分
布图
!!!
!
!
&&
&
&&&&&&
&
&
&
&
&
&
&&
&
&
&
&&&&&&&±
010203040
5千米
Extract_GALa
Res: 1 : 78845.602像元大小: 510.110
值
高 : 140.977低 : 68.2002。