Chapter03第三章 空间平滑和空间插值.docx

1、Chapter03第三章 空间平滑和空间插值35第三章 空间平滑和空间插值本章介绍基于GIS的空间分析中两个常用操作：空间平滑和空间插值。空间平滑和空间插值关系密切，它们都可以用于显示空间分布态式及空间分布趋势，二者还共享某些算法（如核密度估计法Find/Replace All）。空间平滑和空间插值的方法有很多种，本章只介绍其中最常用的几种。空间平滑与移动平均在概念上类似（移动平均是求一个时间段内的均值），而空间平滑术是一个空间窗口内计算平均值。第3.1节介绍空间平滑的概念和方法，第3.2节是案例分析3A，用空间平滑法研究中国南方/泰语地名(Find/Replace all)分布。空间插值是用

2、某些点的已知数值来估算其他点的未知数值。第3.3节介绍了基于点的空间插值，第3.4节为案例3B，演示了一些常用的点插值法。案例3B所用数据与3A相同，是案例3A工作的延伸。第3.5节介绍基于面的空间插值，用一套面域数值（一般面单元较小）来估算另一个面域的数值（范围较大）。面插值可用于数据融合以及不同面域单元的数据整合。第3.6节为案例3C，介绍两种简单的面插值法。第3.7节为小结。3.1空间平滑与移动平均法计算一个时间段的平均值（例如：五日平均温度）相似，空间平滑是将某点周围地区（定义为一个空间窗口）的平均值作为该点的平滑值，以此减少空间变异。空间平滑适用面很广。其中一种应用是处理小样本问题，

3、我们在第八章会详细讨论。对于那些人口较少的地区，由于小样本事件中随机误差的影响，癌症或谋杀等稀有事件发生率的估算不够可靠。对于某些地区，这样的事情发生一次就可导致一个高发生率，而对于另外许多地区，没有发生这种事情的结果是零发生率。另外一种应用是将离散的点数据转化为连续的密度图，从而考察点数据的空间分布模式，可参见下面的第3.2节。本节介绍两种空间平滑方法（移动搜索法及核密度估计法），附录3介绍经验贝叶斯估计。363.1.1移动搜索法移动搜索法（FCA）是以某点为中心画一个圆或正方形作为滤波窗口，用窗口内的平均值（或数值密度）作为该点的值。将窗口在研究区内移动，直到得到所有位置的平均值。平均值的

4、变动性较小，从而实现空间上的平滑效果。FCA也可以用于可达性测量（见第五章第5.2节）等其他研究中。图3.1是由72个网格单元组成的研究区。以网格53为中心的圆定义了一个包含33个网格的窗口（如果网格的中心在圆内，则它属于这个圆），从而这33个网格的平均值为网格33的空间平滑值。将圆心在研究区内不同网格中心之间移动，就得到所有网格的空间平滑值。例如，以网格56为中心的等半径的圆包含的33个网格定义了网格56的滤波窗口。值得注意的是，研究区边界附近的滤波窗口包含网格较少，从而平滑度较低。这种效果称为边缘效应。窗口的大小很重要，需要仔细确定。较大的窗口得到较强的空间平滑效果，从而更好地反映了区域全

5、局分布态势，而不是局部差异；较小窗口则得到相反的结果。我们可以通过试验不同大小的窗口来寻求一个合适的窗口。37案例3A详细介绍了FCA法在ArcGIS中的应用。我们先计算所有点之间的距离（例如欧式距离），然后提取那些小于或等于阈值的距离。在计算距离时，我们可以选用阈值距离作为搜索半径，从而直接得到阈值范围内的距离。这里我们先计算所有点之间的距离，然后通过属性值(距离)大小提取不同窗口内的观察值, 比较灵活些。在ArcGIS中，我们可以基于提取的距离表通过汇总起始点计算属性值的均值来计算在起始总周边范围内观察点的平均值。因为距离表只包含那些阈值范围内的距离，只有在窗口之内的观察值才参与了汇总操作

6、，实现了搜索的效果。这样，我们可以省掉逐个画圆并查询圆内点的反复计算过程。3.1.2核估计核估计与FCA的方法类似。两种方法都要用一个滤波窗口来定义近邻对象。所不同的是，在FCA法中，所有对象的权重相同，而在核估计法中，距离较近的对象，权重较大。这种方法在分析和显示点数据时尤其有用。离散点的数据直接用图表示, 空间趋势往往不明显。核估计可以得到研究对象的一个连续密度图示，即以“波峰”和“波谷”的方式强化的空间分布模式。这种方法也可以用于空间插值。核密度方程的几何意义为：密度分布在每个xi点中心处最高，向外不断降低，当距离中心达到一定阈值范围（窗口的边缘）处密度为0，参见图3.2。网格中心x处的

7、核密度为窗口范围内的密度和：这里K( )为核密度方程，h为阈值，n为阈值范围内的点数，d为数据的维数。参见相关文献中 (Silverman, 1986: 43)常用的核密度方程。例如，当d = 2时，一个常用的核密度方程可以定义为 38这里，为点(xi, yi)和(x, y)之间的离差。与FCA法中窗口的作用类似，这里较大的阈值揭示一种区域分布态势，而较小的阈值则强调局部分布差异(Fotheringham et al., 2000: 46)。ArcGIS内置有核估计工具。在调用之前，我们先要打开空间分析扩展模块，可以通过主工具条的扩展键来实现，即点击空间分析下拉箭头，点击Density，对于对

8、话框中的Density Type，选择Kernel。 3.2案例3A：用空间平滑法分析华南地区的傣族姓氏分布本例研究华南地区傣族姓氏的分布模式，它是一个研究华南地区傣族历史起源研究项目之一部分。本项目的合作者包括北伊利诺伊大学的约翰哈特曼（John Hartmann）和罗卫。在中国，少数民族（例如傣族）的汉化是一个长期持续的过程。历史变迁的痕迹可以在地名上反映出来。我们发现，一些早期的傣族地名常常以地理或自然界的事物而命名，如“稻田”、“乡村”、“河口”、“山”等。另外的一些傣族地名则在汉化过程中逐渐湮没或改变。我们的研究项目主要是为了重构早期的傣族地名，以便揭示汉族南迁之前华南地区原始傣族居

9、民点的分布范围。本案例用来演示GIS技术在历史语言文化研究中的应用，这是一个学者较少涉猎的领域。我们的研究区为广西钦州市（参见图3.3）。地图是研究空间分布模式的一种重要方法，但是直接标绘傣族地名能够读取的信息不多。图3.3是傣族与非傣族地名的分布，由此可以粗略地看到傣族地名在空间分布上疏密有别。空间平滑技术可以形象地显示傣族地名的空间分布模式。本例需要用到光盘中的下述数据：1.钦州市乡镇地名的点图层qztai，属性TAI为地名的傣族（=1）或非傣族（=0）标记。2.qzcnty为研究区内6个县的边界图层。3.2.1第一部分：基于移动搜索法的空间平滑首先应用移动搜索法进行研究。我们要试验不同的

10、窗口大小，寻找一个适中窗口，在这个过程,我们希望既有一定的平滑程度以便显示总体分布态势，又能揭示地区差异。在围绕某点的窗口内，傣族地名在所有地名中的比率代表该点周围傣族地名的集中度。在实际应用时，关键的一步是利用任意两点的距离矩阵来提取某点周围一定半径范围内的地名点。1.39计算各点之间的距离矩阵：参照第2.3.1节的办法测算欧式距离。在ArcToolbox中，依次选择Analysis Tools Proximity Point Distance。在“Input Features”和 “Near Features”栏都输入qztai（point），将输出表命名为Dist_50km.dbf。“s

11、earch radius”取50km，这样我们可以利用距离表处理50km以内的不同窗口。在距离表Dist_50km.dbf中，列数据INPUT_FID为起点，而NEAR_FID为终点。2.将傣族地名连接到距离矩阵：以qztai中的FID和Dist_50km.dbf 中的NEAR_FID为连接指针，将属性数据表qztai连接(join)到距离表Dist_50km.dbf。这样，每个终点可以通过属性数据point:Tai来判断是否为傣族地名。3.40提取窗口内的距离矩阵：例如，定义一个10km半径的窗口，打开表Dist_50km.dbf，进行下述操作：单击右下边的“option”，选择Select

12、 By Attributes ，输入条件Dist_50km.DISTANCE Export，输出新表，命名为Dist_10km.dbf，里面为所有距离小于10km的数据。那些距离值为0（distance = 0）的点（即起点和终点相同）为圆心。4.计算窗口内傣族地名的比重：打开表Dist_10km.dbf，右键单击列INPUT_FID选择Summarize，在弹出的对话框中，第一栏（field to summarize）里为INPUT_FID，在第二栏（summary statistics）中选择TAI(Sum)，并将输出表命名为Sum_10km.dbf，所得表中，列Sum_TAI为10km距

13、离内的傣族地名数，而列Count_INPUT_FID为10km内的总地名数。在表Sum_10km.dbf中添加新的一列Tairatio，按照公式Tairatio = Sum_TAI / Cnt_INPUT_计算数值。这里，Cnt_INPUT_为列名Count_INPUT_FID的简写。所得比值为窗口内傣族地名数占所有地名数的比重。5.将傣族地名比重值连接到点图层：以Sum_10km.dbf中的INPUT_FID及qztai中的FID为连接指针，将Sum_10km.dbf连接到qztai的属性数据表。6.绘制傣族地名比重图：用“proportional point symbols”方式绘制傣族地

14、名比重图（比重值代表每点周围10km范围内傣族地名的比重)，见图3.4。上面演示的即为FCA空间平滑法，它将二值变量TAI转化为比值变量Tairatio。7.敏感性分析：用其他窗口值如5 km或15 km，重复上述3-6步的操作，将所得结果与图3.4对比，以考察窗口大小的影响。表3.1是所得数据的一些统计描述值。由此可知，当窗口值增加时，傣族地名比重值的标准离差降低，表明空间平滑性增加。413.2.2第二部分：基于核估计的空间平滑1.42执行核估计：打开ArcMap的空间分析（Spatial Analyst）扩展模块：单击Tools菜单 选择Extensions 选中Spatial Analy

15、st，单击View菜单 选择Toolbars 选中Spatial Analyst。单击Spatial Analyst下拉箭头 选择Density，弹出新的对话框。在对话框中， Input data栏选择qztai(point)，在Population field栏选择TAI，选择kernel作为Density type，设置Search radius为10000 (meters)，Area units为square kilometers，Output cell size为1000 (meters)，将输出栅格数据命名为kernel_10k。2.绘制核密度图：默认状况下，核密度图以9级分色显示，

16、图3.5是按5级显示的结果，背景为县域边界。核密度图上傣族地名的分布为一个连续的面，显示了波峰与波谷的分布态势。但是，图上的密度值只表示相对的集中度，并不象FCA法得到的Tairatio（泰语地名在一定范围内的百分比）有实际的意义。3.3基于点的空间插值基于点的空间插值包括整体和局部两种方法。整体插值借助所有已知点（控制点）的数据来估计未知值。局部插值借助未知点周边的样本来估计未知值。正如托布罗第一地理定律（Tobler, 1970）所述，“所有事物彼此相关，距离越近关系越强”。是用整体插值还是局部插值，取决于远处的控制点是否对待估未知数据点有作用。究竟选取哪一种方法并没有明确的规律可循。可以

17、认为，从整体到局部的尺度是连续的。如果数值主要受邻近的控制点影响，可以用局部插值法。局部插值法的计算量要比整体插值法小得多（Chang, 2004: 277）。我们也可以用检验技术来比较不同的方法。例如，控制点可以分成两个样本：一个用于构建模型，另一个用于检验已经构建的模型的准确性。本节在简单介绍2个整体插值法之后，重点讲解3种局部插值法。3.3.1整体插值法整体插值法包括趋势面分析和回归模型分析两种。趋势面分析是用多项式模型拟合已知数据点 这里，属性数值z被认为是坐标x和y的函数（Bailey and Gatrell, 1995）。例如，一个三次趋势面模型可以表作 上述方程通常用最小二乘法进

18、行估计，然后将拟合所得方程用于估算其他点的值。43一般来说，高阶模型可用于描述复杂表面，从而得到较高的拟合优度R2或较低的RMS。其中RMS（root mean square均方根）的计算方法为：。但是，对控制点拟合较好的模型并不一定是估计未知数值的好模型。有必要对不同模型进行比较检验。如果因变量（即待估属性）是二值变量（即0和1），则模型为一逻辑斯蒂趋势面模型，表征一个概率曲面。局部趋势面分析是用一个未知点周边的控制样本来估计该点的未知数值，通常称为局部多项式插值。ArcGIS提供了最高12阶的趋势面模型。为了实现这种方法，需要打开Geostatistical Analyst扩展模块。在Ar

19、cMap中，操作过程为：单击Geostatistical Analyst下拉箭头 Explore Data Trend Analysis。回归模型是用线形回归法来得到因变量与自变量之间的方程，然后用来估计未知点的数值（Flowerdew and Green, 1992）。回归模型既可以用空间变量（不一定是上述趋势面模型中用到的x-y坐标），也可以用属性变量，而趋势面分析只能用x-y坐标进行预测。3.3.2局部插值法下面讨论三种局部插值法：反距离加权法、薄片样条插值法、克里金法。反距离加权法（IDW）用周边点的加权平均值作为未知点的估计值，这里的权重按距离的幂次衰减（Chang, 2004: 2

20、82）。因此，IDW法是托布罗第一地理定律的例证。IDW模型可以表作 这里zu为待估u点的未知值，zi为控制点i的属性值，diu是点i与u之间的距离，s为所用控制点的数目，k为幂次。幂次越高，距离衰减作用越强（越快）（即近邻点的权重比远处点的权重高得多）。换言之，距离的幂次越高，局部作用越强。薄片样条插值是通过拟合得到一个曲面，对所有控制点的预测值完全拟合，并在所有点的变化率最小（Franke, 1982）。其模型可表作44这里，x和y是未知数据点的坐标，是到控制点（xi, yi）的距离，Ai, a, b和 c待估的n+3个参数。这些参数可以通过解一个由n+3个线形方程组来得到（参见第十一章）

21、，则有 ; ; and .需要注意的是，上面第一个方程实际上代表了i = 1, 2, , n 取值时的n个方程，zi为已知i点的属性值。薄片样条插值法在数据稀少的区域将产生陡峭的梯度值，可以用张力薄片样条插值、规则样条插值、紧缩规则样条插值来减轻这个问题（参见Chang, 2004: 285）。这些高级插值法都可归为径向基函数。克里金法（Krige, 1966）认为空间变异包含三个部分：空间相关组分，代表区域化变量；“漂移”或结构，代表趋势；随机误差。克里金法借助半方差函数（1/2方差 ）来检验自相关： 这里n为相距（或称空间滞后）h的控制点对的数目，z为属性值。由于空间依赖关系，(h) 随h

22、增加而增大，即近邻物体之间的相似性大于远距离物体。可以用半方差图来显示(h) 随h变化的情况。克里金法通过拟合半方差图得到一个数学模型，以此来估计任意给定距离的半方差函数，从而用之计算空间权重。这里所用空间权重的效果与IDW法相似，即近邻控制点的权重比远点的权重高。例如，对于某个未知点s（需要插值的点），控制点i的权重为Wis，则s点的插值为：这里，ns为s周围样本控制点数，zs和zi分别为s和i点的属性值。与核估计相似，克里金法可以基于点数据得到一个连续的面。45在ArcGIS中，三种局部插值都可以在Geostatistical Analyst扩展模块中实现。在ArcMap里，单击Geost

23、atistical Analyst下拉箭头 Geostatistical Wizard 选择 Inverse Distance Weighting、Radial Basis Functions或Kriging来分别调用IDW法、各种薄片样条插值法、或克里金法。这三种局部插值法也可以用Spatial Analyst或3D Analyst来实现。这里推荐Geostatistical Analyst法，因为它提供更多信息和更好的交互界面（Chang, 2004: 298）。3.4案例3B：表面建模及华南傣族地名图的绘制这里延续案例3A的工作，用各种表面建模技术绘制钦州市傣族地名的空间集聚图。所用数据

24、不变。同时还会用到前面案例3A第一部分所生成的数据，尤其是用FCA法计算得到的傣族地名比重的数据。3.4.1第一部分：用趋势面分析法制图1.激活Geostatistical Wizard对话框：在案例3A中，如果退出ArcMap时没有保存工程，则需要重复案例3A第一部分第5步的工作：将表Sum_10km.dbf连接到属性表qztai。在ArcMap中，打开Geostatistical Analyst和Spatial Analyst扩展模块。单击Geostatistical Analyst下拉箭头选择Geostatistical Wizard，弹出对话框。2.用趋势面分析生成趋势面：在第一步弹出

25、的对话框中，在输入数据框（Input Data）中选择qztai，属性框（Attribute）中选择Sum_10km.Tairatio，方法框（Methods）中选择Global Polynomial Interpolation。在下一个对话框中，用不同的幂次（例如1，3，4，8，10）分别试验，观察所得趋势面及对应的RMS值。随着幂次的增加，趋势面包含的局部信息越多，得到的RMS值越小。这里，我们取幂次为10，得到的RMS=0.1124，生成的趋势面图层Global Polynomial Interpolation Prediction Map将自动添加到图层中。3.绘制研究区的趋势面图：右

26、键单击趋势面图层，选择Data Export to Raster。将输出栅格图命名为trend10。单击Spatial Analyst下拉箭头选择Options 设置qzcnty为Analysis mask。再次单击Spatial Analyst下拉箭头 选择Raster Calculator 双击图层一栏中的trend10，将其添加到计算框中，再单击Evaluate。计算得到的栅格图Calculation即为研究区内的趋势面图。右键单击Calculation图层，选择属性（Properties）以改善显示效果（例如，在Display选项卡中，定义透明度为30%；在Symbology选项卡中，

27、修改图例等）。图3.6为傣族地名比重的趋势面图，它跟用核估计得到的图3.5的分布态势相似，但显示了更多整体趋势。它清楚地表明，傣族地名高度集中在西南地区，并向东北及其他方向延伸。值得注意的是，趋势面分析的有些插值如负值和大于1的值在现实中并不存在。4.46可选操作：逻辑斯蒂趋势面分析：ArcGIS也可以直接用原始的二值变量Tai来生成趋势面。在第二步的对话框中，属性（Attribute）一栏选择point:Tai，其他选项不变，重复上述操作，即可得到基于逻辑斯蒂趋势面分析的概率面（即某点为傣族地名的概率）。3.4.2第二部分：用局部插值法绘制分布图1.用IDW法绘制趋势面：与上面第一部分第1步

28、类似，打开Geostatistical Wizard对话框。选择qztai为输入数据（Input），Sum_10km.Tairatio为属性（Attribute），选择方法框里面的Inverse Distance Weighting。使用默认的幂次为2，设置邻近点数（neighbor to选项）为15，并用圆域来选择控制点。计算得到的RMS = 0.0844。将所得趋势面输出为栅格图idw2。与上述第一部分第3步类似，我们将得到研究区内的趋势面，如图3.7所示。需要注意的是，所有插值与原始数据的范围相同，即在01之间。与图3.6相比，图3.7的局部分布态势更明显。 2.47用薄片样条线法绘制趋

29、势面：类似的，在Geostatistical Wizard对话框中，选择Radial Basis Functions作为插值方法，其他选项不变。单击next，设置插值的参数，使用默认的Completely Regularized Spline法，其他项使用默认设置即可，于是得到一个新的趋势面，将所得栅格图命名为regspline。所得图与3.7略有不同，在此从略。 3.用克里金法绘制趋势面：类似的，在Geostatistical Wizard对话框中，选择Kriging作为插值方法，其他选项不变。单击next，设置插值的参数，使用默认的Ordinary Kriging Prediction M

30、ap法，其他项使用默认设置即可，于是得到一个新的趋势面，将所得栅格图命名为ordkrig。所得图形与图3.7类似，在此从略。 3.5基于面域的空间插值基于面域的插值（面域插值）也称为交叉面域聚合，它是将数据从一种面域单元系统（源区域）转换到另一种面域系统（目标区域）。点插值中用到的一些方法如克里金法或多阶趋势面分析也可用于栅格表面的插值，原始栅格单元的值转化为目标区域的值。换言之，假设面状单元可以用它的质心代替，基于点的插值方法可以近似地用于面状单元的属性插值。48面域插值的方法有很多种（Goodchild et al., 1993），其中最简单也最常的是面积权重插值（Goodchild an

31、d Lam, 1980）。这种方法将源区域的属性值按面积比例分配到目标区域。它假设属性值在每个源区域内均匀分布。如果研究区的信息较多，我们还可以用一些更高级的方法来改进插值。下面介绍一种对美国人口普查数据进行插值时尤其有用的方法。利用美国统计局TIGER数据中的路网数据，谢一春（Xie, 1995; 也可参见Batty and Xie, 1994a, 1994b）发展了一些网络覆盖算法将人口或其他基于居民的属性数据从一个面域单元投影到另一个面域单元。居民住房通常沿街道或道路分布，从而人口分布与街道网络紧密相关。在“网络长度”、“网络等级权重”和“网络住房载荷”三种算法中，网络等级权重（NHW）

32、法得到的结果最为理想。NHW法的关键部分由一系列GIS叠加操作组成。下面举例说明。在研究城市问题时，我们常用交通规划普查数据库（CTPP）来分析土地利用和交通问题。1990年CTPP数据：http:/www.bts.gov/publications/census_transportation_planning_package_1990。2000年CTPP数据：http:/www.fhwa.dot.gov/ctpp/在1990年的CTPP城市要素数据中，大部分区域的数据是按交通分析小区（TAZ）这一层次进行汇总的。例如，在1990年303个CTPP城市要素区域中，265个区域的数据是按TAZ汇总的，13个是按普查小区进行的，25个按街区进行（作者根据美国交通统计局发布的文件Regncode.asc进行汇总）。因为种种原因（例如