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

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

本章介绍基于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介绍经验贝叶斯估计。

35

3.1.1移动搜索法

移动搜索法（FCA）是以某点为中心画一个圆或正方形作为滤波窗口，用窗口内的平均值（或数值密度）作为该点的值。将窗口在研究区内移动，直到得到所有位置的平均值。平均值的变动性较小，从而实现空间上的平滑效果。FCA也可以用于可达性测量（见第五章第5.2节）等其他研究中。

图3.1是由72个网格单元组成的研究区。以网格53为中心的圆定义了一个包含33个网格的窗口（如果网格的中心在圆内，则它属于这个圆），从而这33个网格的平均值为网格33的空间平滑值。将圆心在研究区内不同网格中心之间移动，就得到所有网格的空间平滑值。例如，以网格56为中心的等半径的圆包含的33个网格定义了网格56的滤波窗口。值得注意的是，研究区边界附近的滤波窗口包含网格较少，从而平滑度较低。这种效果称为边缘效应。

窗口的大小很重要，需要仔细确定。较大的窗口得到较强的空间平滑效果，从而更好地反映了区域全局分布态势，而不是局部差异；较小窗口则得到相反的结果。我们可以通过试验不同大小的窗口来寻求一个合适的窗口。

案例3A详细介绍了FCA法在ArcGIS中的应用。我们先计算所有点之间的距离（例如欧式距离），然后提取那些小于或等于阈值的距离。在计算距离时，我们可以选用阈值距离作为搜索半径，从而直接得到阈值范围内的距离。这里我们先计算所有点之间的距离，然后通过属性值(距离)大小提取不同窗口内的观察值, 比较灵活些。在ArcGIS中，我们可以基于提取的距离表通过汇总起始点计算属性值的均值来计算在起始总周边范围内观察点的平均值。因为距离表只包含那些阈值范围内的距离，只有在窗口之内的观察值才参与了汇总操作，实现了搜索的效果。这样，我们可以省掉逐个画圆并查询圆内点的反复计算过程。

3.1.2核估计

核估计与FCA的方法类似。两种方法都要用一个滤波窗口来定义近邻对象。所不同的是，在FCA法中，所有对象的权重相同，而在核估计法中，距离较近的对象，权重较大。这

36

37

种方法在分析和显示点数据时尤其有用。离散点的数据直接用图表示, 空间趋势往往不明显。核估计可以得到研究对象的一个连续密度图示，即以“波峰”和“波谷”的方式强化的空间分布模式。这种方法也可以用于空间插值。

核密度方程的几何意义为：密度分布在每个xi点中心处最高，向外不断降低，当距离中心达到一定阈值范围（窗口的边缘）处密度为0，参见图3.2。网格中心x处的核密度为窗口范围内的密度和：

1nf?(x)?x?xinhd?K(i?1h) 这里K( )为核密度方程，h为阈值，n为阈值范围内的点数，d为数据的维数。参见相关文献中 (Silverman, 1986: 43)常用的核密度方程。例如，当d = 2时，一个常用的核密度方程可以定义为

1n(x?x22

f?(x)?i)?(y?y)2nh2??[1?i2] i?1h这里，(x?x)2?(y?y)2ii为点(xi, yi)和(x, y)之间的离差。 与FCA法中窗口的作用类似，这里较大的阈值揭示一种区域分布态势，而较小的阈值则强调局部分布差异(Fotheringham et al., 2000: 46)。

ArcGIS内置有核估计工具。在调用之前，我们先要打开空间分析扩展模块，可以通过主工具条的扩展键来实现，即点击空间分析下拉箭头，点击Density，对于对话框中的Density Type，选择Kernel。

3.2案例3A：用空间平滑法分析华南地区的傣族姓氏分布

本例研究华南地区傣族姓氏的分布模式，它是一个研究华南地区傣族历史起源研究项目之一部分。本项目的合作者包括北伊利诺伊大学的约翰·哈特曼（John Hartmann）和罗卫。在中国，少数民族（例如傣族）的汉化是一个长期持续的过程。历史变迁的痕迹可以在地名上反映

38

出来。我们发现，一些早期的傣族地名常常以地理或自然界的事物而命名，如“稻田”、“乡村”、“河口”、“山”等。另外的一些傣族地名则在汉化过程中逐渐湮没或改变。我们的研究项目主要是为了重构早期的傣族地名，以便揭示汉族南迁之前华南地区原始傣族居民点的分布范围。本案例用来演示GIS技术在历史－语言－文化研究中的应用，这是一个学者较少涉猎的领域。

我们的研究区为广西钦州市（参见图3.3）。地图是研究空间分布模式的一种重要方法，但是直接标绘傣族地名能够读取的信息不多。图3.3是傣族与非傣族地名的分布，由此可以粗略地看到傣族地名在空间分布上疏密有别。空间平滑技术可以形象地显示傣族地名的空间分布模式。

本例需要用到光盘中的下述数据：

1. 钦州市乡镇地名的点图层qztai，属性TAI为地名的傣族（=1）或非傣族（=0）标

记。

2. qzcnty为研究区内6个县的边界图层。 3.2.1第一部分：基于移动搜索法的空间平滑

首先应用移动搜索法进行研究。我们要试验不同的窗口大小，寻找一个适中窗口，在这个过程,我们希望既有一定的平滑程度以便显示总体分布态势，又能揭示地区差异。在围绕某点的窗口内，傣族地名在所有地名中的比率代表该点周围傣族地名的集中度。在实际应用时，关键的一步是利用任意两点的距离矩阵来提取某点周围一定半径范围内的地名点。

1. 计算各点之间的距离矩阵：参照第2.3.1节的办法测算欧式距离。在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

39

为连接指针，将属性数据表qztai连接(join)到距离表Dist_50km.dbf。这样，每个终点可以通过属性数据point:Tai来判断是否为傣族地名。

3. 提取窗口内的距离矩阵：例如，定义一个10km半径的窗口，打开表

Dist_50km.dbf，进行下述操作：单击右下边的“option”，选择Select By Attributes ，输入条件Dist_50km.DISTANCE <=10000，执行操作后，对每个起点，所有10km距离内的终点将被选中。点击Options > 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是所得数据的一些统计描述值。由此可知，当窗口值增加时，傣族地名比重值的标准离差降低，表明空间平滑性增加。

40

3.2.2第二部分：基于核估计的空间平滑

1. 执行核估计：打开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种局部插值法。

41

42

3.3.1整体插值法

整体插值法包括趋势面分析和回归模型分析两种。趋势面分析是用多项式模型拟合已知数据点

z?f(x,y)

这里，属性数值z被认为是坐标x和y的函数（Bailey and Gatrell, 1995）。例如，一个三次趋势面模型可以表作

z(x,y)?b2232230?b1x?b2y?b3x?b4xy?b5y?b6x?b7xy?b8xy?b9y 上述方程通常用最小二乘法进行估计，然后将拟合所得方程用于估算其他点的值。

一般来说，高阶模型可用于描述复杂表面，从而得到较高的拟合优度R2

或较低的RMS。

n其中RMS（root mean square均方根）的计算方法为：RMS??(z?z)2i,obsi,est/n。但i?1是，对控制点拟合较好的模型并不一定是估计未知数值的好模型。有必要对不同模型进行比较检验。如果因变量（即待估属性）是二值变量（即0和1），则模型为一逻辑斯蒂趋势面模型，表征一个概率曲面。局部趋势面分析是用一个未知点周边的控制样本来估计该点的未知数值，通常称为局部多项式插值。

ArcGIS提供了最高12阶的趋势面模型。为了实现这种方法，需要打开Geostatistical Analyst扩展模块。在ArcMap中，操作过程为：单击Geostatistical Analyst下拉箭头> Explore Data > Trend Analysis。

回归模型是用线形回归法来得到因变量与自变量之间的方程，然后用来估计未知点的数值（Flowerdew and Green, 1992）。回归模型既可以用空间变量（不一定是上述趋势面模型中用到的x-y坐标），也可以用属性变量，而趋势面分析只能用x-y坐标进行预测。 3.3.2局部插值法

下面讨论三种局部插值法：反距离加权法、薄片样条插值法、克里金法。

43

反距离加权法（IDW）用周边点的加权平均值作为未知点的估计值，这里的权重按距离的幂次衰减（Chang, 2004: 282）。因此，IDW法是托布罗第一地理定律的例证。IDW模型可以表作

s?zidiuzu?i?1s?k

?kiu?di?1这里zu为待估u点的未知值，zi为控制点i的属性值，diu是点i与u之间的距离，s为所用控制点的数目，k为幂次。幂次越高，距离衰减作用越强（越快）（即近邻点的权重比远处点的权重高得多）。换言之，距离的幂次越高，局部作用越强。

薄片样条插值是通过拟合得到一个曲面，对所有控制点的预测值完全拟合，并在所有点的变化率最小（Franke, 1982）。其模型可表作

nz(x,y)?Adlnd?a?bx?cy ?iiii?12?(x?x)?(y?y)是到控制点（xi, yi）的距这里，x和y是未知数据点的坐标，diii22离，Ai, a, b和 c待估的n+3个参数。这些参数可以通过解一个由n+3个线形方程组来得到（参见第十一章），则有

n44

Adlnd?a?bx?cy?z;?

iiiii?1nninii2ii

?Ai?1?0;

?Axi?1?0; and

?Ayii?1i?0.

需要注意的是，上面第一个方程实际上代表了i = 1, 2, …, n 取值时的n个方程，zi为已知i点的属性值。

薄片样条插值法在数据稀少的区域将产生陡峭的梯度值，可以用张力薄片样条插值、规则样条插值、紧缩规则样条插值来减轻这个问题（参见Chang, 2004: 285）。这些高级插值法都可归为径向基函数。

克里金法（Krige, 1966）认为空间变异包含三个部分：空间相关组分，代表区域化变量；“漂移”或结构，代表趋势；随机误差。克里金法借助半方差函数（1/2方差 ）来检验自相关：

?n12(h)??[z(x)?z(x?h)] ii2ni?1这里n为相距（或称空间滞后）h的控制点对的数目，z为属性值。由于空间依赖关系，γ(h) 随h增加而增大，即近邻物体之间的相似性大于远距离物体。可以用半方差图来显示γ(h) 随h变化的情况。

克里金法通过拟合半方差图得到一个数学模型，以此来估计任意给定距离的半方差函数，从而用之计算空间权重。这里所用空间权重的效果与IDW法相似，即近邻控制点的权重比远点的权重高。例如，对于某个未知点s（需要插值的点），控制点i的权重为Wis，则s点的插值为：

nszs??Wzisi

i?1这里，ns为s周围样本控制点数，zs和zi分别为s和i点的属性值。与核估计相似，克里金法可以基于点数据得到一个连续的面。

在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）。

45

3.4案例3B：表面建模及华南傣族地名图的绘制

这里延续案例3A的工作，用各种表面建模技术绘制钦州市傣族地名的空间集聚图。所用数据不变。同时还会用到前面案例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. 用趋势面分析生成趋势面：在第一步弹出的对话框中，在输入数据框（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的分布态势相似，但显示了更多

业球队的影响区时，我们使用了一种简单空间插值法来得到居民选择不同球队的概率趋势面（参见图4.4）。但是，趋势面仅仅是描述性的。对空间集聚或分散的判断常常很武断。哪些地方的集聚具有统计上的显著性而不是偶然情况呢？要回答这个问题需要严格的统计分析，例如空间聚类分析——本书第九章将讨论这个问题（案例9A基于同一套数据来判断傣族地名的空间聚类）。

基于面域的空间插值常常用于数据整合分析，即将数据从不同源区域转换到一个面域单元。它也可用于研究可变地域单元问题（MAUP）时将数据从一个高精度面域到一个低精度面52 域。例如，在案例6关于城市密度方程的问题中，将数据从普查小区转换到乡镇单元，从而使基于不同面域单元的方程具有可比性。

附录3：空间平滑的经验贝叶斯估计

经验贝叶斯估计是另外一种常用于调整或平滑面域变量（尤其是比率）的方法（例见Clayton and Kaldor, 1987; Cressie, 1992）。因为两个事件的联合概率等于第一个事件的概率与第二个事件条件概率（基于第一个事件）之积，在估计数据的概率分布时，贝叶斯推断可看作关于数据集内生的先验信息或推断（Langford, 1994: 143），即：

似然函数×先验概率=后验概率

以疾病风险为例，研究区的观测数据可用泊松分布的似然函数表示。先验信息是基于研究区观测数据的相对风险（率）分布：例如，人口较多地区的相对风险估计比人口较少地区的估计要可靠得多。总而言之：（1）研究区内的平均风险率是可靠和无偏的；（2）跟对小规模人口的估计相比，对大规模人口的疾病风险率估计更准确；（3）疾病风险率服从一种已知的概率分布。

假定用一种概率分布如Γ分布来描述先验的风险率分布。Γ分布有两个参数，即形态参数α和尺度参数ν，均值为ν/α，方差为ν/α。参数α和ν可以综合用马歇尔（Marshall, 1991）提出的最大似然法和力矩法进行估算。对于某个人口数为Pi，病例数为ki的 i地区，发病的原始概率为ki /Pi。由此可以得到，后验期望率或经验贝叶斯估计为

2

ki??P??iEi? 如果地区i的人口数较少，与ν和α相比，ki和Pi都较小，从而经验贝叶斯估计Ei 接近于ν/α（全研究区的比率）。相反地，如果地区i的人口较多，ki和Pi相对较大，从而经验贝叶斯估计Ei接近于原始概率ki /Pi。经验贝叶斯估计Ei是原始概率ki /Pi经ν和α两个内生参数平滑的结果。

当经验贝叶斯估计（EB）用于整个研究区时，所有地区的比率都基于整个研究区的比率进行了平滑，这即是全局经验贝叶斯平滑。当它用于局部地区时，基于每个地区定义一个邻域，将比率按邻域地区的比率进行平滑，此即为局部经验贝叶斯平滑。一个地区的邻域定义为53 它的邻近地区与它本身之和。邻域可以按“边邻域”或“任意邻域”（参见第1.4节）、一阶或二阶邻域等方式定义。

鲁科·安索林（Luc Anselin）和同事合作开发了一种免费软件包GeoDa（http://sal.agecon.uiuc.edu/geoda_main.php），它可以用来做EB估计的空间平滑。在GeoDa0.9.5-i中，依次选择Map > Smooth > Empirical Bayes（或Spatial Empirical Bayes）即可。这里的Empirical Bayes命令将比率按整个研究区均值进行平滑，因而是一种全域经验贝叶斯平滑。Spatial Empirical Bayes命令是将比率按一个地区空间“窗口”（基于一个空间权重文件，按某个地区及其邻域地区进行定义）进行平滑，因而是局部经验贝叶斯平滑。

表3.1 基于不同窗口大小的FCA空间平滑

窗口大小 （半径） 5 km 10 km 15 km 最小值 0 0 0 最大值 1.0 1.0 0.8333 傣族地名的比重 均值 0.1868 0.1886 0.1878 标准差 0.3005 0.1986 0.1642

本文来源：https://www.bwwdw.com/article/q1c8.html