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

35第三章 空间平滑和空间插值本章介绍基于GIS的空间分析中两个常用操作：空间平滑和空间插值。空间平滑和空间插值关系密切，它们都可以用于显示空间分布态式及空间分布趋势，二者还共享某些算法（如核密度估计法Find/Replace All）。空间平滑和空间插值的方法有很多种，本章只介绍其中最常用的几种。空间平滑与移动平均在概念上类似（移动平均是求一个时间段内的均值），而空间平滑术是一个空间窗口内计算平均值。第3.1节介绍空间平滑的概念和方法，第3.2节是案例分析3A，用空间平滑法研究中国南方/泰语地名(Find/Replace all)分布。空间插值是用某些点的已知数值来估算其他点的未知数值。第3.3节介绍了基于点的空间插值，第3.4节为案例3B，演示了一些常用的点插值法。案例3B所用数据与3A相同，是案例3A工作的延伸。第3.5节介绍基于面的空间插值，用一套面域数值（一般面单元较小）来估算另一个面域的数值（范围较大）。面插值可用于数据融合以及不同面域单元的数据整合。第3.6节为案例3C，介绍两种简单的面插值法。第3.7节为小结。3.1空间平滑与移动平均法计算一个时间段的平均值（例如：五日平均温度）相似，空间平滑是将某点周围地区（定义为一个空间窗口）的平均值作为该点的平滑值，以此减少空间变异。空间平滑适用面很广。其中一种应用是处理小样本问题，我们在第八章会详细讨论。对于那些人口较少的地区，由于小样本事件中随机误差的影响，癌症或谋杀等稀有事件发生率的估算不够可靠。对于某些地区，这样的事情发生一次就可导致一个高发生率，而对于另外许多地区，没有发生这种事情的结果是零发生率。另外一种应用是将离散的点数据转化为连续的密度图，从而考察点数据的空间分布模式，可参见下面的第3.2节。本节介绍两种空间平滑方法（移动搜索法及核密度估计法），附录3介绍经验贝叶斯估计。36移动搜索法移动搜索法（FCA）是以某点为中心画一个圆或正方形作为滤波窗口，用窗口内的平均值（或数值密度）作为该点的值。将窗口在研究区内移动，直到得到所有位置的平均值。平均值的变动性较小，从而实现空间上的平滑效果。FCA也可以用于可达性测量（见第五章第5.2节）等其他研究中。图3.1是由72个网格单元组成的研究区。以网格53为中心的圆定义了一个包含33个网格的窗口（如果网格的中心在圆内，则它属于这个圆），从而这33个网格的平均值为网格33的空间平滑值。将圆心在研究区内不同网格中心之间移动，就得到所有网格的空间平滑值。例如，以网格56为中心的等半径的圆包含的33个网格定义了网格56的滤波窗口。值得注意的是，研究区边界附近的滤波窗口包含网格较少，从而平滑度较低。这种效果称为边缘效应。窗口的大小很重要，需要仔细确定。较大的窗口得到较强的空间平滑效果，从而更好地反映了区域全局分布态势，而不是局部差异；较小窗口则得到相反的结果。我们可以通过试验不同大小的窗口来寻求一个合适的窗口。37案例3A详细介绍了FCA法在ArcGIS中的应用。我们先计算所有点之间的距离（例如欧式距离），然后提取那些小于或等于阈值的距离。在计算距离时，我们可以选用阈值距离作为搜索半径，从而直接得到阈值范围内的距离。这里我们先计算所有点之间的距离，然后通过属性值(距离)大小提取不同窗口内的观察值, 比较灵活些。在ArcGIS中，我们可以基于提取的距离表通过汇总起始点计算属性值的均值来计算在起始总周边范围内观察点的平均值。因为距离表只包含那些阈值范围内的距离，只有在窗口之内的观察值才参与了汇总操作，实现了搜索的效果。这样，我们可以省掉逐个画圆并查询圆内点的反复计算过程。核估计核估计与FCA的方法类似。两种方法都要用一个滤波窗口来定义近邻对象。所不同的是，在FCA法中，所有对象的权重相同，而在核估计法中，距离较近的对象，权重较大。这种方法在分析和显示点数据时尤其有用。离散点的数据直接用图表示, 空间趋势往往不明显。核估计可以得到研究对象的一个连续密度图示，即以“波峰”和“波谷”的方式强化的空间分布模式。这种方法也可以用于空间插值。核密度方程的几何意义为：密度分布在每个xi点中心处最高，向外不断降低，当距离中心达到一定阈值范围（窗口的边缘）处密度为0，参见图3.2。网格中心x处的核密度为窗口范围内的密度和：这里K( )为核密度方程，h为阈值，n为阈值范围内的点数，d为数据的维数。参见相关文献中 (Silverman, 1986: 43)常用的核密度方程。例如，当d = 2时，一个常用的核密度方程可以定义为38这里，为点(xi, yi)和(x, y)之间的离差。与FCA法中窗口的作用类似，这里较大的阈值揭示一种区域分布态势，而较小的阈值则强调局部分布差异(Fotheringham et al., 2000: 46)。ArcGIS内置有核估计工具。在调用之前，我们先要打开空间分析扩展模块，可以通过主工具条的扩展键来实现，即点击空间分析下拉箭头，点击Density，对于对话框中的Density Type，选择Kernel。 3.2案例3A：用空间平滑法分析华南地区的傣族姓氏分布本例研究华南地区傣族姓氏的分布模式，它是一个研究华南地区傣族历史起源研究项目之一部分。本项目的合作者包括北伊利诺伊大学的约翰哈特曼（John Hartmann）和罗卫。在中国，少数民族（例如傣族）的汉化是一个长期持续的过程。历史变迁的痕迹可以在地名上反映出来。我们发现，一些早期的傣族地名常常以地理或自然界的事物而命名，如“稻田”、“乡村”、“河口”、“山”等。另外的一些傣族地名则在汉化过程中逐渐湮没或改变。我们的研究项目主要是为了重构早期的傣族地名，以便揭示汉族南迁之前华南地区原始傣族居民点的分布范围。本案例用来演示GIS技术在历史语言文化研究中的应用，这是一个学者较少涉猎的领域。我们的研究区为广西钦州市（参见图3.3）。地图是研究空间分布模式的一种重要方法，但是直接标绘傣族地名能够读取的信息不多。图3.3是傣族与非傣族地名的分布，由此可以粗略地看到傣族地名在空间分布上疏密有别。空间平滑技术可以形象地显示傣族地名的空间分布模式。本例需要用到光盘中的下述数据：1. 钦州市乡镇地名的点图层qztai，属性TAI为地名的傣族（=1）或非傣族（=0）标记。2. qzcnty为研究区内6个县的边界图层。第一部分：基于移动搜索法的空间平滑首先应用移动搜索法进行研究。我们要试验不同的窗口大小，寻找一个适中窗口，在这个过程,我们希望既有一定的平滑程度以便显示总体分布态势，又能揭示地区差异。在围绕某点的窗口内，傣族地名在所有地名中的比率代表该点周围傣族地名的集中度。在实际应用时，关键的一步是利用任意两点的距离矩阵来提取某点周围一定半径范围内的地名点。391. 计算各点之间的距离矩阵：参照第节的办法测算欧式距离。在ArcToolbox中，依次选择Analysis Tools Proximity Point Distance。在“Input Features”和 “Near Features”栏都输入qztai（point），将输出表命名为Dist_50km.dbf。“search 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来判断是否为傣族地名。403. 提取窗口内的距离矩阵：例如，定义一个10km半径的窗口，打开表Dist_50km.dbf，进行下述操作：单击右下边的“option”，选择Select 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距离内的傣族地名数，而列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”方式绘制傣族地名比重图（比重值代表每点周围10km范围内傣族地名的比重)，见图3.4。上面演示的即为FCA空间平滑法，它将二值变量TAI转化为比值变量Tairatio。7. 敏感性分析：用其他窗口值如5 km或15 km，重复上述3-6步的操作，将所得结果与图3.4对比，以考察窗口大小的影响。表3.1是所得数据的一些统计描述值。由此可知，当窗口值增加时，傣族地名比重值的标准离差降低，表明空间平滑性增加。41第二部分：基于核估计的空间平滑421. 执行核估计：打开ArcMap的空间分析（Spatial Analyst）扩展模块：单击Tools菜单 选择Extensions 选中Spatial Analyst，单击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级分色显示，图3.5是按5级显示的结果，背景为县域边界。核密度图上傣族地名的分布为一个连续的面，显示了波峰与波谷的分布态势。但是，图上的密度值只表示相对的集中度，并不象FCA法得到的Tairatio（泰语地名在一定范围内的百分比）有实际的意义。3.3基于点的空间插值基于点的空间插值包括整体和局部两种方法。整体插值借助所有已知点（控制点）的数据来估计未知值。局部插值借助未知点周边的样本来估计未知值。正如托布罗第一地理定律（Tobler, 1970）所述，“所有事物彼此相关，距离越近关系越强”。是用整体插值还是局部插值，取决于远处的控制点是否对待估未知数据点有作用。究竟选取哪一种方法并没有明确的规律可循。可以认为，从整体到局部的尺度是连续的。如果数值主要受邻近的控制点影响，可以用局部插值法。局部插值法的计算量要比整体插值法小得多（Chang, 2004: 277）。我们也可以用检验技术来比较不同的方法。例如，控制点可以分成两个样本：一个用于构建模型，另一个用于检验已经构建的模型的准确性。本节在简单介绍2个整体插值法之后，重点讲解3种局部插值法。整体插值法整体插值法包括趋势面分析和回归模型分析两种。趋势面分析是用多项式模型拟合已知数据点这里，属性数值z被认为是坐标x和y的函数（Bailey and Gatrell, 1995）。例如，一个三次趋势面模型可以表作上述方程通常用最小二乘法进行估计，然后将拟合所得方程用于估算其他点的值。43一般来说，高阶模型可用于描述复杂表面，从而得到较高的拟合优度R2或较低的RMS。其中RMS（root mean square均方根）的计算方法为：。但是，对控制点拟合较好的模型并不一定是估计未知数值的好模型。有必要对不同模型进行比较检验。如果因变量（即待估属性）是二值变量（即0和1），则模型为一逻辑斯蒂趋势面模型，表征一个概率曲面。局部趋势面分析是用一个未知点周边的控制样本来估计该点的未知数值，通常称为局部多项式插值。ArcGIS提供了最高12阶的趋势面模型。为了实现这种方法，需要打开Geostatistical Analyst扩展模块。在ArcMap中，操作过程为：单击Geostatistical Analyst下拉箭头 Explore Data Trend Analysis。回归模型是用线形回归法来得到因变量与自变量之间的方程，然后用来估计未知点的数值（Flowerdew and Green, 1992）。回归模型既可以用空间变量（不一定是上述趋势面模型中用到的x-y坐标），也可以用属性变量，而趋势面分析只能用x-y坐标进行预测。局部插值法下面讨论三种局部插值法：反距离加权法、薄片样条插值法、克里金法。反距离加权法（IDW）用周边点的加权平均值作为未知点的估计值，这里的权重按距离的幂次衰减（Chang, 2004: 282）。因此，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个线形方程组来得到（参见第十一章），则有; ; and .需要注意的是，上面第一个方程实际上代表了i = 1, 2, , n 取值时的n个方程，zi为已知i点的属性值。薄片样条插值法在数据稀少的区域将产生陡峭的梯度值，可以用张力薄片样条插值、规则样条插值、紧缩规则样条插值来减轻这个问题（参见Chang, 2004: 285）。这些高级插值法都可归为径向基函数。克里金法（Krige, 1966）认为空间变异包含三个部分：空间相关组分，代表区域化变量；“漂移”或结构，代表趋势；随机误差。克里金法借助半方差函数（1/2方差 ）来检验自相关： 这里n为相距（或称空间滞后）h的控制点对的数目，z为属性值。由于空间依赖关系，(h) 随h增加而增大，即近邻物体之间的相似性大于远距离物体。可以用半方差图来显示(h) 随h变化的情况。克里金法通过拟合半方差图得到一个数学模型，以此来估计任意给定距离的半方差函数，从而用之计算空间权重。这里所用空间权重的效果与IDW法相似，即近邻控制点的权重比远点的权重高。例如，对于某个未知点s（需要插值的点），控制点i的权重为Wis，则s点的插值为：这里，ns为s周围样本控制点数，zs和zi分别为s和i点的属性值。与核估计相似，克里金法可以基于点数据得到一个连续的面。45在ArcGIS中，三种局部插值都可以在Geostatistical Analyst扩展模块中实现。在ArcMap里，单击Geostatistical Analyst下拉箭头 Geostatistical Wizard 选择 Inverse Distance Weighting、Radial Basis Functions或Kriging来分别调用IDW法、各种薄片样条插值法、或克里金法。这三种局部插值法也可以用Spatial Analyst或3D Analyst来实现。这里推荐Geostatistical Analyst法，因为它提供更多信息和更好的交互界面（Chang, 2004: 298）。3.4案例3B：表面建模及华南傣族地名图的绘制这里延续案例3A的工作，用各种表面建模技术绘制钦州市傣族地名的空间集聚图。所用数据不变。同时还会用到前面案例3A第一部分所生成的数据，尤其是用FCA法计算得到的傣族地名比重的数据。第一部分：用趋势面分析法制图1. 激活Geostatistical Wizard对话框：在案例3A中，如果退出ArcMap时没有保存工程，则需要重复案例3A第一部分第5步的工作：将表Sum_10km.dbf连接到属性表qztai。在ArcMap中，打开Geostatistical Analyst和Spatial Analyst扩展模块。单击Geostatistical Analyst下拉箭头选择Geostatistical Wizard，弹出对话框。2. 用趋势面分析生成趋势面：在第一步弹出的对话框中，在输入数据框（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. 绘制研究区的趋势面图：右键单击趋势面图层，选择Data Export to Raster。将输出栅格图命名为trend10。单击Spatial Analyst下拉箭头选择Options 设置qzcnty为Analysis mask。再次单击Spatial Analyst下拉箭头 选择Raster Calculator 双击图层一栏中的trend10，将其添加到计算框中，再单击Evaluate。计算得到的栅格图Calculation即为研究区内的趋势面图。右键单击Calculation图层，选择属性（Properties）以改善显示效果（例如，在Display选项卡中，定义透明度为30%；在Symbology选项卡中，修改图例等）。图3.6为傣族地名比重的趋势面图，它跟用核估计得到的图3.5的分布态势相似，但显示了更多整体趋势。它清楚地表明，傣族地名高度集中在西南地区，并向东北及其他方向延伸。值得注意的是，趋势面分析的有些插值如负值和大于1的值在现实中并不存在。464. 可选操作：逻辑斯蒂趋势面分析：ArcGIS也可以直接用原始的二值变量Tai来生成趋势面。在第二步的对话框中，属性（Attribute）一栏选择point:Tai，其他选项不变，重复上述操作，即可得到基于逻辑斯蒂趋势面分析的概率面（即某点为傣族地名的概率）。第二部分：用局部插值法绘制分布图1. 用IDW法绘制趋势面：与上面第一部分第1步类似，打开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的局部分布态势更明显。 472. 用薄片样条线法绘制趋势面：类似的，在Geostatistical Wizard对话框中，选择Radial Basis Functions作为插值方法，其他选项不变。单击next，设置插值的参数，使用默认的Completely Regularized Spline法，其他项使用默认设置即可，于是得到一个新的趋势面，将所得栅格图命名为regspline。所得图与3.7略有不同，在此从略。 3. 用克里金法绘制趋势面：类似的，在Geostatistical Wizard对话框中，选择Kriging作为插值方法，其他选项不变。单击next，设置插值的参数，使用默认的Ordinary Kriging Prediction Map法，其他项使用默认设置即可，于是得到一个新的趋势面，将所得栅格图命名为ordkrig。所得图形与图3.7类似，在此从略。 3.5基于面域的空间插值基于面域的插值（面域插值）也称为交叉面域聚合，它是将数据从一种面域单元系统（源区域）转换到另一种面域系统（目标区域）。点插值中用到的一些方法如克里金法或多阶趋势面分析也可用于栅格表面的插值，原始栅格单元的值转化为目标区域的值。换言之，假设面状单元可以用它的质心代替，基于点的插值方法可以近似地用于面状单元的属性插值。48面域插值的方法有很多种（Goodchild et al., 1993），其中最简单也最常的是面积权重插值（Goodchild and Lam, 1980）。这种方法将源区域的属性值按面积比例分配到目标区域。它假设属性值在每个源区域内均匀分布。如果研究区的信息较多，我们还可以用一些更高级的方法来改进插值。下面介绍一种对美国人口普查数据进行插值时尤其有用的方法。利用美国统计局TIGER数据中的路网数据，谢一春（Xie, 1995; 也可参见Batty and Xie, 1994a, 1994b）发展了一些网络覆盖算法将人口或其他基于居民的属性数据从一个面域单元投影到另一个面域单元。居民住房通常沿街道或道路分布，从而人口分布与街道网络紧密相关。在“网络长度”、“网络等级权重”和“网络住房载荷”三种算法中，网络等级权重（NHW）法得到的结果最为理想。NHW法的关键部分由一系列GIS叠加操作组成。下面举例说明。在研究城市问题时，我们常用交通规划普查数据库（CTPP）来分析土地利用和交通问题。1990年CTPP数据：。2000年CTPP数据：在1990年的CTPP城市要素数据中，大部分区域的数据是按交通分析小区（TAZ）这一层次进行汇总的。例如，在1990年303个CTPP城市要素区域中，265个区域的数据是按TAZ汇总的，13个是按普查小区进行的，25个按街区进行（作者根据美国交通统计局发布的文件Regncode.asc进行汇总）。因为种种原因（例如，一些普查小区的数据合并在一起），我们需要将基于TAZ的CTPP数据转普到普查小区上。在这里，TAZ为源区域，普查小区为目标区域。这一操作包括下述5步：1. 把TAZ图层与普查小区图层进行叠加以得到新的图层TAZ-tract（多边形），将TAZ-tract图层与道路网络图层叠加得到一个控制网（线）图层。2. 构造公路权重矩阵，不同类型公路沿线分布的人口和商业密度不同，从而赋予不同权重。3. 将TAZ图层和网络图层叠加，计算不同公路的长度以及每个TAZ内的加权长度，将人口及其他一些属性分配给网络。4. 将第3步得到的结果（人口及其他属性）连接到控制网图层，然后基于TAZ-tract图层计算每个多边形内属性值之和。5. 按普查小区计算汇总属性值，从而得到每个普查小区的属性插值。3.6案例3C：将俄亥俄克利夫兰市普查数据从普查小区转化到邻里单元和入学片区本例介绍两种常用的基于面域的数据整合方法。第一种本质上并不是插值法，它只是在每个目标区域包含多个源区域时（或者假设近似如此），简单地将数据从一种面域单元转换到另一种面域单元。第二种方法是面积权重插值。第一种方法将在下面第一部分介绍，我们将把克利夫兰市的普查数据从普查小区转换到邻里单元，每个邻里单元包括多个完整的普查小区。具体操作时以人口属性数据为例来演示这种插值方法。49本书光盘附有本例所需下述数据： 1. 克利夫兰市36个邻里单元（或称统计规划区，SPA）的shp文件clevspa2k。2. 库伊合各县统一入学片区的shp文件tgr39035uni。3. 库伊合各县的普查小区shp文件cuyautm及普查小区的质心文件cuya_pt。3.6.1克利夫兰市普查小区数据向邻里单元数据的简单融合这种简单融合通过空间连接的办法将数据从源区域转换到目标区域，直接操作即可得到结果。在ArcMap中，打开cuyautm的属性表，它包括人口密度数据，但没有人口数据。因为当时的属性连接是暂时性的，所以在前面案例1A中添加的数据列area和popuden都还在，但那些连接的列在退出工程时就没有了。我们可以用两种方法来恢复POP2000数据，一是重新连接表tgr39035trt00.dbf，二是添加一列POP2000，然后按照POP2000 = area * popuden / 1000000的公式计算得到。对质心shp文件cuya_pt进行同样的操作。右键单击目标图层clevspa2k，选择Joins and Relates Join。在弹出的对话框中，选择“Join data from another layer based on spatial location”，选择cuya_pt作为源图层，在数据汇总选项（其说明文字为“Each polygon will be given a summary of the numerical attributes of the points that fall inside it ”）中，选中“Sum”前面的复选框，将输出结果命名为spa_pop。如果选择cuyautm为源图层，则可执行多边形到多边形的连接。比较保险的做法是上述点到多边形连接。spa_pop属性表中的sum_pop2000数据列即为汇总的邻里单元人口数（其他列可删去）。 库伊合各县从普查小区到入学片区（Find/Relaceby“校区”）的面积权重数据融合501. 准备入学片区图层：在ArcToolbox中，调用Data Management Tools Projections and Transformations Feature Project，将地理坐标系的tgr39035uni转换到UTM坐标系的cuyauni文件，这里的UTM从cuyautm文件导入（import）。在cuyauni 的属性表中添加一个新的变量area列，并更新其数据（参照第1.2节第3步）。512. 普查小区与入学片区图层叠加：在ArcToolbox中，调用Analysis Tools Overlay Intersect 先选择cuyautm作为输入要素，然后再添加cuyauni作为输入要素，将输出结果命名为tmp_int。在tmp_int的属性表中，数据列area是从cuyautm得到的普查小区面积，数据列area_1是从cuyauni得到的每个入学片区的面积。图3.8以研究区右下角的一块区域为例演示了叠加过程。需要注意的是，普查小区184104被分成了两部分：多边形A属于入学片区10016，多边形B属于入学片区04660。3. 根据面积大小分配属性数据：在tmp_int的属性数据中添加一列area_2，使其面积等于新叠加生成的shp文件的面积。在tmp_int的属性数据中再新增一列popu_est，并按照公式popu_est = pop2000*area_2/area求值。这样，我们就按面积权重法得到了叠加操作后各多边形的人口插值。例如，多边形A的人口数等于普查小区184104的人口乘以多边形A的面积再除以普查小区184104的面积，即1468*1297600/16810900=113。为了检验所得数据，读者可以新添加一列popu_valid，并按照popu_valid = popuden*area_2 /1000000公式进行计算，所得结果应该与popu_est相同。面积权重法假设每个普查小区内的人口呈均匀分布，从而叠加后的多边形（普查小区的一部分）具有普查小区一样的人口密度。4. 将数据转换到入学片区：在tmp_int的属性表中，右键单击列unified（入学片区编码），在第二个对话框中选择Summarize 选择popu_est（sum），并将输出结果命名为uni_pop.dbf。uni_pop.dbf 表的数据列sum_popu_e即为入学片区的人口插值，它可以连接到cuyauni图层以便进行绘图等其他操作。3.7小结本章介绍的技巧包括下面几种：1. 用FCA法进行空间平滑；2. 用核估计绘制点状要素数据；3. 趋势面分析（包括逻辑斯蒂趋势面分析）；4. 局部插值法，如反距离权重法、薄片样条线法、克里金法；5. 简单数据融合，要求多个源数据区域全部包含在目标区域内；6. 面积权重数据融合，当源数据区域与目标数据区域边界不同时适用。空间平滑和空间插值经常用于绘制空间分布模式图，如案例3A和3B那样绘制华南地区傣族地名分布图。这两种技术可应用于许多基于点状数据的工作。例如，在案例4A中定义专业球队的影响区时，我们使用了一种简单空间插值法来得到居民选择不同球队的概率趋势面（参见图4.4）。但是，趋势面仅仅是描述性的。对空间集聚或分散的判断常常很武断。哪些地方的集聚具有统计上的显著性而不是偶然情况呢？要回答这个问题需要严格的统计分析，例如空间聚类分析本书第九章将讨论这个问题（案例9A基于同一套数据来判断傣族地名的空间聚类）。52基于面域的空间插值常常用于数据整合分析，即将数据从不同源区域转换到一个面域单元。它也可用于研究可变地域单元问题（MAUP）时将数据从一个高精度面域到一个低精度面域。例如，在案例6关于城市密度方程的问题中，将数据从普查小区转换到乡镇单元，从而使基于不同面域单元的方程具有可比性。附录3：空间平滑的经验贝叶斯估计经验贝叶斯估计是另外一种常用于调整或平滑面域变量（尤其是比率）的方法（例见Clayton and Kaldor, 1987; Cressie, 1992）。因为两个事件的联合概率等于第一个事件的概率与第二个事件条件概率（基于第一个事件）之积，在估计数据的概率分布时，贝叶斯推断可看作关于数据集内生的先验信息或推断（Langford, 1994: 143），即：似然函数先验概率=后验概率以疾病风险为例，研究区的观测数据可用泊松分布的似然函数表示。先验信息是基于研究区观测数据的相对风险（率）分布：例如，人口较多地区的相对风险估计比人口较少地区的估计要可靠得多。总而言之：（1）研究区内的平均风险率是可靠和无偏的；（2）跟对小规模人口的估计相比，对大规模人口的疾病风险率估计更准确；（3）疾病风险率服从一种已知的概率分布。假定用一种概率分布如分布来描述先验的风险率分布。分布有两个参数，即形态参数和尺度参数，均值为/，方差为/2。参数和可以综合用马歇尔（Marshall, 1991）提出的最大似然法和力矩法进行估算。对于某个人口数为Pi，病例数为ki的 i地区，发病的原始概率为ki /Pi。由此可以得到，后验期望率或经验贝叶斯估计为如果地区i的人口数较少，与和相比，ki和Pi都较小，从而经验贝叶斯估计Ei 接近于/（全研究区的比率）。相反地，如果地区i的人口较多，ki和Pi相对较大，从而经验贝叶斯估计Ei接近于原始概率ki /Pi。经验贝叶斯估计Ei是原始概率ki /Pi经和两个内生参数平滑的结果。53当经验贝叶斯估计（EB）用于整个研究区时，所有地区的比率都基于整个研究区的比率进行了平滑，这即是全局经验贝叶斯平滑。当它用于局部地区时，基于每个地区定义一个邻域，将比率按邻域地区的比率进行平滑，此即为局部经验贝叶斯平滑。一个地区的邻域定义为它的邻近地区与它本身之和。邻域可以按“边邻域”或“任意邻域”（参见第1.4节）、一阶或二阶邻域等方式定义。鲁科安索林（Luc Anselin）和同事合作开发了一种免费软件包GeoDa（），它可以用来做EB估计的空间平滑。在GeoDa-i中，依次选择Map Smooth Empirical Bayes（或Spatial Empirical Bayes）即可。这里的Empirical Bayes命令将比率按整个研究区均值进行平滑，因而是一种全域经验贝叶斯平滑。Spatial Empirical Bayes命令是将比率按一个地区空间“窗口”（基于一个空间权重文件，按某个地区及其邻域地区进行定义）进行平滑，因而是局部经验贝叶斯平滑。 表3.1 基于不同窗口大小的FCA空间平滑窗口大小（半径）傣族地名的比重最小值最大值均值标准差5 km01.00.18680.300510 km01.00.18860.198615 km00.83330.18780.1642