基于矢量数据的GIS分析（2）

四、矢量网络分析

在现实世界中，若干线状要素(如某种资源、物质或信息在地理空间上的流动) 相互连接构成网络，组成网络。GIS网络分析就是通过模拟、分析网络的状态以及资源在网络上的流动和分配等，研究网络结构、流动效率及网络资源等的优化问题的一种方法。对地理网络（交通网）、城市基础设施网络（电力线、电话线、给排水管线等）进行地理分析和模块化，是GIS网络分析的主要目的。第3章介绍了在GIS中表示和存储网络的数据模型与数据结构，这里，重点讨论矢量数据的网络分析。在网络分析中，路径分析是GIS网络分析中最基本的功能，其核心是对最短路径、最佳路径以及服务范围进行分析。在实际工作中，网络数据随着应用的目的、范围、对象而变化，相互间的拓扑关系更是动态、复杂而且数量庞大。对于基于矢量数据的网络分析系统而言，还需事先做好显式数据，组织各顶点数据，然后取得所有两结点之间的路径长度，再通过Dijkstra或其他算法取得全部顶点两两之间的路径长度，一旦一个数据变化，全部组织需要重新进行。

1.最短路径分析

最短路径分析是计算和寻找网络中任何两点之间距离最短的路径。然而，最短路径不一定以距离衡量，也可以时间、费用等衡量，衡量标准依赖于分析的目的。例如，一辆救护车出发到一个出事地点，所关心的是最快的路径，最快路径往往需要考虑道路的交通拥挤情况等限制因素。

最短路径分析需要计算网络中从起点到终点所有可能的路径，从中选择距离最短的一条。用于最短路径分析的算法很多，其中最著名的是Dijkstra算法（Dijkstra,1959）,该算法可描述如下。

设一个网络由k个结点组成，以N＝{n_i│i=1,2,…，k}表示结点集，其中一个结点为起始结点，设其为n。Dijkstra将N划分为两个子集，一个子集包含那些到起始点的最短距离已确定的结点，称这些结点为已确定结点，以S表示这一子集；另一个子集包含未确定结点，即它们到起点结点的最短距离尚未确定，以Q表示这一子集。又设d为一个距离矩阵（Array），存放每个结点到n_s的最短距离，d（i）表示结点n_i到n_s的最短距离；p为一前置结点矩阵，存放由n_s到其他结点的最短路径上每个结点的前一个结点，p（i）表示结点n_i在最短路径上的前一个结点。已知每两个相连结点之间的距离（或它们之间路径的长度），Dijkstra算法按如下几个步骤运行。

（1）将d和p初始化，使d的每个元素值为无穷大，p的每个元素值为空值，并设S和Q为空集；

（2）将n_s加入Q，令d（s）＝0。

（3）从Q中找出到n_s最短距离为最小的结点，设该结点为n_u。

（4）将n_u加入S，并将它从Q中删除。

（5）找出与n_u相连的所有结点，从这些结点中取出一个，令其为n_v，

①如nv已存在于S中，则执行下面第②步，否则，作如下判断：

如果d(v)< d (u)+n_u和n_v之间的距离，执行第（2）步；

否则令d（v）＝d（u）＋n_u和n_v之间的距离；P（v）=n_u；将n_v加进Q。

②如果与nu相连的所有结点都已做过上述判断，继续执行第（6）步；否则，取下一个为判断结点，令其为nv，执行上面的第（1）步；

（6）判断Q是否为空集，若不是，或到第（3）步；否则，停止运算。

以图5.15为例说明这一算法。

图5.15最短路径分析

本例使用的网络由四个结点组成，这四个结点为n₁、n₂、n₃和n₄,如图5.15（a）所示。设n₁为起始结点，d（1）＝0。与n₁相连的结点包括n₂、n₃和n₄，先计算n₁到n₂的最短距离：

d（1）＋n₁到n₂的距离=0+1=1

由于d（2）的初值无穷大，因此，令

d（2）＝1

p（2）＝n₁

将n₂加入Q。 类似地，得出d（3）＝4，p（3）＝n₁，d（4）＝5，p（4）＝n₁。n₃和n₄也先后被加入Q。到此，Q中有三个结点n₂、n₃和n₄，因此，进入下一个环节。

在Q包含的三个结点中，n₂到n₁的最短距离d（2）＝1，是Q中当前到n₁最短距离的最小值，因此，将n₂加入S，并从Q中将其删除。寻找与n₂相连的结点，这些包括n₁，n₃和n_4。n₁已存在于S中，故不予考虑。首先计算n₃到n₁的最短距离：

d（2）＋n₂到n₃的距离＝1＋2＝3

但是，在上个循环中获得的d（3）＝4 >3，可见在n₁和n₃之间还有一条更近的路径，为此，给d（3）和p（3）重新赋值，即令

d（3）＝3

p（3）＝n₂

类似地，重新计算d（4）和p（4），令

    d（4）＝4

p（4）＝n₂

这时，Q中有两个结点n₃和n₄，因此，进入第三个循环。

在当前Q包含的两个结点中，n₃到n₁的距离d（3）＝3，小于d（4），因此，将n₃加入S，然后，从Q中将它删除。找出与n₃相连的结点，即n₁、n₂、和n₄，由于n₁和n₂都已存在于S中，故只需计算n₄到n₁的距离：

d（3）＋n₃到n₄的距离＝3＋2＝5

因为d（4）＝4<5，因此，d（4）和p（4）的值都不变。现在Q中只剩下n₄，与n₄相连接的所有结点n₁，n₂，n₃都已在S中，将它从Q中删除，算法到此结束。

上述过程产生的矩形d＝（0，1，3，4）包含了每个结点到起始结点n₁的最短距离，而矩阵p＝（n₁、n₂、n₃）则包含了每个结点到n₁的最短路径上它的前一个结点，由p可以根据前置点跟踪出每个结点到n₁的最短路径，图5.15（b）以箭头表示了本例中推算出来的到n₁的最短路径。

在Dijkstra算法中，若将两个结点的距离改成通过这两个结点之间路径所需的时间或费用，则获得的结果将是最快路径或最少费用路径。Dijkstra算法适用于起点和终点都是网络结点的情况，如果起点和终点位于网络路径上，但它们不是结点，就需要将它们所在的路径在起点和终点处分割开来，使起点和终点成为网络结点，然后重建拓扑结构，这样才能应用Dijkstra算法。在ArcView 3.3 的网络分析（Network Analyst）模块中能执行这一算法，可用于描述从一地到另一地（如某一居民住址到一所学校沿街的）最短路径或一种资源的最佳分配等。

现代GIS的网络分析功能还可输出从起点到终点沿最短路径前进的方向、经过的街道名称和距离等方面的信息。

除了寻找起点到终点的最短路径以外，GIS一般还提供另外两种类型的最短路径分析功能。一种是寻找最近设施（Nearest Facility），即从分布于网络中的提供某种服务的一系列设施或服务中心中，寻找最靠近一个给定地点的服务设施或中心。例如，寻找距离某一居住区的最近的医疗诊所等。另一种是路径规划（Routing Planning），指给定全程所需访问或停靠的地点，制定访问这些地点的顺序和最短路径）。

例如，以福州市区为例，需要为游客提供从乌山历史风貌区到闽江公园的最佳旅游路线，利用Dijkstra算法构建最佳路径模型，可假设：

（1）公交车已经覆盖福州五区的任何一条道路，其和自行车在任何一条道路上都可以通行；

（2）不考虑道路拥堵的情况，不考虑路口存在红绿灯；

（3）不考虑道路存在高架桥，道路可以双向通行。

再分别以时间为权重，或以道路长度为权重，或以时间和道路长度两者为权重算出三条最佳路径（如图5.16）。最后考虑计算以不同的旅游路线从乌山历史风貌区到闽江公园所需要的成本。给游客方便出行提供参考。

图5.16  某区域基于模型的三条最佳路径显示

2.服务范围分析

服务范围分析功能是用于寻找或确定围绕任何一个位于网络中的服务设施的服务范围（Service Area）。这里的服务设施可以是商店、学校、派出所、救护站、消防站等，服务范围是指从一个服务设施在一定的时间或距离内所能到达的区域。服务范围分析的结果可以用来衡量某类服务实施在网络中的分布的有效性和合理性，如确定某个城市现有的消防站在两分钟之内所到达的区域。若用ArcView 或ArcGIS中的网络分析（Network Analyst）模块可计算和输出围绕三个购物中心的1km服务范围，若将这些服务范围与居民住宅区及其居住人口数据叠置起来，可以统计出这三个购物中心在1km范围内可服务的人口总数。

服务范围分析的算法相对比较简单，其中一种算法称为宽度优先搜索法（Breadth First Search）。给定一个服务设施所在的网络结点以及一定的距离或时间界限（即搜索范围），这种算法以该服务设施为起始结点，首先找出与它相连接的所有结点，对于这些相邻结点中每一个位于搜索范围内的结点，寻找与它们相连接的结点，判断哪些结点是否位于搜索范围内。如此下去，直到不再存在任何一个结点，其相连结点位于搜索范围内。根据此算法搜索到的所有位于搜索范围内的结点，由它们之间的路径构成的网格所覆盖的区域即为服务范围。图5.17是设置设置响应时间和服务范围的网络分析图示（具体基于软件实现操作见《地理信息系统导论实验指导》一书）

图5.17 设置响应时间和服务范围的网络分析图

上述算法如考虑服务设施的供应量（如一所小学所能招收的一年级学生的数量）和每节路径的需求量（如每个街道现有的小学适学儿童数量），则可以用于划分服务设施的资源分配范围。在这种情况下，搜索范围根据供应量定义，算法将在需求量快要超过供应量是停止结点的搜索。现有的GIS软件所提供的服务范围分析功能虽在网络数据模型、功能大小、运算效率、应用范围上都有较大的不同，但一般都可利用网络数据结构中所定义的网络特征数据进行计算。例如，对于道路网，网络的特征数据可以包括路段的长度、限速、规定的行驶方向、转弯方向、转弯阻强等，这些特征数据对于最短路径的计算、服务范围的划分都有很大的影响。通过运用网络的特征数据进行网络分析，可以更加真实地描述和模拟网络中资源、物质、能量或信息的流动，为规划和决策提供可靠的数据。

值得注意的是，在矢量数据的网络分析中，一方面，由于矢量数据本身结构和矢量网络分析所基于的图论知识的限制。因图论概念众多、结构复杂、方法多、组织形式多、变化多、地学数据量大、精度要求相对苛刻，数据组织和输入的困难极大地影响了它的广泛应用；另一方面，在基于矢量数据进行大型网络分析时，由于矢量数据基于的是点、线、面，分析计算复杂，其算法效率问题更显著。虽然各种局部的改良方法很多，但这些方法的共同弱点是连接的拓扑关系的数据以及相应的几何数据（距离）是运算的必要数据，因而在庞大网络经常性地动态变化时，例如故障和维修引起的中断，线路和权重改变等，维护和更新这些拓扑数据、几何数据是十分困难的，并且这些变化将引起结构的整体变化。因此算法本身效率不高以及数据维护、更新困难等问题使得矢量数据在进行网络分析应用和推广上碰到了巨大困难。目前，人们对于如何有效进行GIS网络分析还在不断探索。

五、矢量地形分析

地形分析是对地形环境认知的一种重要手段，GIS地形分析(也称“GIS表面分析”)与数据结构关系密切。不同的数据结构对应不同的模型。不同的模型（如等高线、TIN、DEM等）对地面真实地形的模拟状况是不同的。在矢量GIS中，面对复杂地形分析多以不规则三角网数据结构（TIN）为基础（见第3章），提取基本地形要素，包括坡度、坡向、地面点高程、地势剖面、通视情况等。

1.基于矢量GIS的坡度和坡向的计算

坡度是地表单元陡缓的程度，由地表单元上两点之间的垂直距离（vd）和水平距离(hd)计算，可以度数或百分比表示，如图5.18（a）所示。坡向为地表单元的朝向，定义为地表单元法线（与地表单元表面相垂直的线）在平面上的投影与正北方向的夹角，自正北（0°）起按顺时针方向计算，如图5.18（b）所示。在GIS分析中，通常将坡向划分为四个方向：东（E）、西（W）、南（S）、北（N），或八个方向：东、东南（SE）、东北（NE）、西、西南（SW）、西北（NW）、南、北。

图5.18坡度和坡向

TIN是由不规则的、互不重叠的三角形构成的连续表面，每个三角形代表一个地形单元，三角形的顶点为高程数据点，每个点具有一个三维坐标（x，y，z），（x，y）为点在某一坐标系统中的平面坐标，z为其高程。设一个三角形的三个顶点分别为A（x₁，y₁，z₁），B(x₂，y₂，z₂)和C（x₃，y₃，z₃），则该三角形的坡度S（%）可由（5.3）式计算。

         (5.3)

则地表坡度D（°）可由（5.4）式计算.

  

根据推算坡向的规则可计算出该三角形的坡向。

2.等高线和地形剖面的绘制

等高线（Contour）是地面上高程相等的点连接起来所形成的闭合曲线，是地形图上表示地形特征的主要手段。一条等高线反映的是某一高度的地形平面轮廓，而一组等高线则以其疏密变化反映地形坡度的变化。等高线密度越密集的地方，其地面坡度愈陡；越稀疏的地方，地面坡度愈缓。相邻两等高线的高差称为等高距。传统的地形分析是通过等高线的阅读和手工测量，在当今的GIS时代，运用等高线手工测量、计算地形要素已不多见，但在许多情况下，等高线仍然是一种有效手段。据根TIN（或 DEM），GIS可以自动跟踪和描绘等高线。 

在TIN中，假设沿着每个三角形的边高程是连续变化的，那么运用TIN绘制等高线的方法可简单地描述如下：

给定某一高程h，检查每个三角形的边，根据三角形边两个端点的高程判断高程h的等高线是否穿过它们。如果等高线穿过三角形的边，那么就使用线性插值的方法确定这条等高线在该三角形边上的位置。设一个三角形边的两端顶点坐标为（x₁，y₁，z₁）和(x₂，y₂，z₂)，穿过该边的高程h的等高线在其上的位置为（x，y，h），则：

在确定了该等高线在所有其穿过的三角形边上的位置以后，将这些位置连接起来即构成了这根等高线，如图5.19所示。起初的等高线是以直线段相连的，运用GIS线段光滑技术可将它转变成光滑的曲线。

图5.19  根据TIN内插500米等高线

就像传统的利用等高线绘制地势剖面图一样，GIS可以根据TIN和DEM绘制地势剖面图。地势剖面图以曲线表示出某一方向线上地势的起伏和坡度的陡缓。GIS根据TIN绘制的地势剖面图的方法如下。

(1)在表示地势的TIN上作一条方向线。

(2)计算方向线与每个三角形边的交点，记录每个交点在方向线上的平面位置，然后利用线性内插的方法，根据相交三角形边两端顶点的坐标和高程计算出交点的高程。

(3)选择合适的垂直比例尺以决定剖面的高度。

(4)作与方向线等长的水平线，以此为x轴，以高程为y轴建立平面直角坐标系。

以第（2）步计算出的交点在方向线上的位置为x，高程乘以垂直比例尺的积为y，绘出每个交点，再将它们连接成一条曲线，完成由TIN绘制而成的地势剖面图。

由TIN可生成DEM。图5.20为闽西某流域地势剖面示意图。

图5.20 闽西某流域地势剖面示意图

有关坡度、坡向及地形剖面提取具体操作参见《地理信息系统导论实验指导》一书

3.通视情况分析

根据TIN模型进行视线运算和视域分析，可视区表述一般通过计算地形中单个的三角形面元可视的部分来实现。Clarke（1990）介绍了一种用于通视情况分析的“光线追溯法”(Ray Tracing)，这种方法首先以方向线或视线将观察点与目标点相连，从目标点向观察点回溯，或从观察点向目标点搜索，只要发现沿方向线有一地面点高出同一点位的方向线，就说明这一点的地形挡住了目标点，即观察点与目标点是不相通视的。

图5.21 光线追溯法

如图5.21所示，设观察点与目标之间有一地面点p，p的高程为H_p，然而在同一点方向线的高度为H_r，H_p>   H_r，p高出了方向线，因此，从观察点看不到目标点。若已知目标点到观察点的水平距离为D，p到观察点的距离为D_p，则H_r可通过计算如下：

在以TIN判断两点间通视情况时，在两点之间作一视线，计算视线与每个三角形边的交点，并内插出交点的高程，如果一个交点的高程大于H_r，那么这两点是不相通视的。在使用TIN进行视域分析中，为了计算上的简便，往往是将研究区域划分成具有一定网格大小的格网，以每个格网作为目标，根据上述方法从给定的观察点出发，判断观察点到每个网格的通视情况，并将网格分成两类：可见区和不可见区，最后以栅格数据形式输出视域分析结果。

六、矢量邻域分析

从给定位置的某现象的“值”去推算给定邻域相关变量的值，这在GIS中通常归结为“空间插值（或逼近）”，属于邻域分析或趋势分析。在实际中应用较广。例如，根据气象站观测点雨量、气温、湿度等气象数据，估算区域内其他各个地点的雨量、气温或湿度等值；根据从某一田块不同地点获取的土壤样本测取的土壤pH值，估算田块其他地点的土壤pH值，等等。

1.趋势面分析

趋势面分析（Trend Surface Analysis）是以数学模型来拟合观测点数据、建立光滑数学曲面的方法。这种数学曲面称为趋势面，根据趋势面可以估算出未知点的值。最常用的建立趋势面的数学方法是使用多项式函数（Polynomial Function）。根据多项式的次数可以将趋势面划分为一次趋势面、二次趋势面、三次趋势面等。

一次趋势面以一个在空间上倾斜的平面来拟合观测点数据，其数学模型可表达为：

        (5.7)

(5.7)式中（x，y）为某一点的平面直角坐标，z为某一类地理实体在该点的属性值，b₀，b₁，b₂为多项式系数。当某一类地理实体的数值由一个方向向另一个方向递增或递减时，一次趋势面可以很好地逼近观测点数据。然而，在现实世界中，很少的地理变量值是随着地点呈线性变化的，常常需要使用高于一次的多项式拟合。二次趋势面以二次多项式拟合，适合于当某一类地理实体的数值在空间上呈抛物面时使用，其数学公式为：

通过增加趋势面次数，可以用更加复杂的曲面来拟合观测点数据，一个p次趋势面可以写成通式：

该多项式模型有（p+1）(p+2)/2个系数，{b_rs | r=0,1,…,p；s=0,1, …,p；r+s≤p}。求解b_rs是多元回归中的一个标准问题，常采用最小二乘法原理。假设有n个观测点数据{z_i| i = 1，2，…，n}，在点（x_i，y_i）处的数值为理论值，记为 ^z_i，令：

为了使趋势面能很好地拟合观测点数据，必须使Q趋于最小，然后建立多元线性方程，对每个b_rs求偏导，令这些偏导数等于0，得到正规方程组，通过解正规方程组即可求得每个b_rs，再将它们代入趋势面方程公式z即可。用它求解区域内任意一点的数值，从而达到空间插值的目的。

在GIS中，趋势面通常以栅格数据输出，每个网格值都由趋势面方程计算。根据输出的栅格数据可以追踪出等值线，因此，趋势面可以栅格图层或等值线表示。

在图5.22中，以降水量值的空间插值为例，给出趋势面分析。趋势面的拟合程度可用统计分析中的F分布进行检验。

图5.22  趋势面分析等降雨量线

在GIS中，趋势面通常以栅格数据输出，每个网格值都由趋势面方程计算。根据输出的栅格数据可以追踪出等值线，因此，趋势面可以栅格图层或等值线表示。

在图5.22中，以降水量值的空间插值为例，给出趋势面分析。趋势面的拟合程度可用统计分析中的F分布进行检验。

2.按距离加权法（距离倒数权重法）

按距离加权法（Inverse Distance Weighting，简称IDW）假设某一未知点的数值是其一定大小的邻域内所有观测点数据的按距离加权平均。设待计算点的值为z₀，s为与其最邻近的或在其邻域中的，将用于插值的观测点数目，这些观测点的值分别为z₁，z₂，…，z_s，它们到待计算点的距离分别d₁，d₂，…，d_s，以（5.11）式表示。

（5.11）式中r为幂方。当r=1时，上式为线性插值（LinearInterpolation），随着r的增大，距离对z₀的影响也相应地增大。按距离加权插值法通常也以栅格数据输出，并可由此产生等值线用于结果的表示。图5.23是利用有限站点的数据，使用距离加权法内插生成的福建平均气温趋势图。

图5.23 基于距离加权法生成的福建省平均气温趋势图

由于“按距离加权插值法”计算较简单，在估值方面应用较多。但要注意以下几个问题。

（1）插值结果受r值的影响很大。若根据不同r值估算的同一未知点的值会有很大的差别。

（2）由于当任何一个d_i=0（1≤i≤s）时，公式无意义，因此，在内插一点时，如果该点与某一已知观测点同点，必须以同点的观测点数作为该点的输出值，从而可能导致输出数据在这一点上的不连续。

（3）按距离加权插值的结果一般都会产生许多围绕某些观测点的封闭等值线，像坐落在丘陵地带的一个个圆形小山包或圆形洼地。由于公式使用的距离都为正值，插值结果都必然在观测点数据极值范围内，既不可能大于观测点数据的最大值，也不可能小于观测点数据的最小值。因此，如果已知的观测点数据中不包含有某个局部地区的峰值时，在该出现峰值的地方获得的内插值将会低于附近周围其他点的值，从而在输出的结果中，应该表示为山峰的地方却被表示为一块低洼盆地，如见图5.22所示。因此，用户应认真评价使用该方法产生的插值结果，以避免作出不合理的判断。目前尚无数学方法用于检验这种插值方法的效果。

3.样条函数法

样条函数（Spline）是模仿手工样条经过一系列数据点绘制光滑曲线的数学方法，也可用于根据一系列观测点内插出一个光滑曲面表示连续分布的面状实体。在GIS软件中，常用的用于曲面内插的样条函数为薄板样条（Thin_Plate Spline）。薄板样条以一个最小曲率的光滑曲面拟合观测点数据，即曲面各处的梯度变化达到最小。薄板样条函数的数学表达（Franke,1982）见（5.12）式。

（5.12）式中z₀为待计算点的值，（x，y）为待计算点的平面坐标，n为用于计算z₀的一个局部区域内观测点数目，d_i为待计算点到第i个观测点之间的距离，A_i，a，b和c为待定系数，它们可以通过求解下面n+3个线性方程组获得：

式中d_j,i表示第j个观测点到第i个观测点的距离，z₁，z₂，…z_n为n个观测点的已知值。

薄板样条函数的主要问题是在数据点或观测点稀少的地区会产生很大的梯度变化。为了纠正这一问题，人们已提出了几种改进的样条函数，包括规则化样条函数（Regularized Spline）(Mitas and Mitasova,1988)和张力薄板样条函数（Thinplate Spline with Tension）(Franke,1985)。图5.24显示的是由规则化样条函数内插出来的气温。由样条函数产生的等值线显得非常平滑。样条函数比较适合于呈连续、平滑曲面状分布的面状实体，如地势、雨量等。在某些情况下，则可能会产生不切合实际的平滑效果。

图5.24使用规则化样条函数内插的气温值

4.克里（金）格法

克里（金）格法（Kriging）是根据该方法的创始人，南非采矿工程师Danie Krige命名的，是一种用于空间插值的地学统计方法（GeostatisticalMethod）。克里格法使用区域化变量（Regionalized Variable）的概念，这种变量的值随着地点而变化，具有一定的连续性，但不能以一个单一的平滑数学方程来模拟，许多地形表面、土壤质量的区域变化等都具有这种特性。这类现象的空间变化既不是随机的，也不是确定的，而是由三个部分组成：一是变化的总体趋势（Drift），即现象成连续面状分布的总体结构。二是空间自相关（Spatial Autocorrelation）变化，即偏离总体变化趋势的一些小的变化，这些变化虽然是随机的，但在空间上是相互关联的。例如，在地形表面出现的一些小峰谷，它们的高程与周围地点的高程是呈连续变化的，即相关的。三是随机误差或随机噪音（Random Noise），它们既与总体变化趋势无关，又没有空间上的相互关系。例如地面上突然出现的巨大砾石。克里格法对区域化变量的这三个成分分别采用不同的方法进行模拟和估算。总体趋势以模拟变量区域性变化的数学方程表达和计算，如趋势面方程；空间自相关变化和随机噪音通过使用称为半方差图（Semivariogram）的统计方法估算。

半方差（Semivariance）是衡量观测点或样本点数据之间空间依赖性的一个统计量，它的大小取决于观测点之间的距离。观测点数据在距离为h的半方差r（h）定义为所有相距h距离的观测点数据的方差除以一半，可用（5.13）数学公式表达：

(5.13)式中:h为观测点之间的距离，n为相距h距离的观测点的个数，z（x_i）为观测点x_i的数据值，z（x_i+h）为与观测点x_i相距h距离的另一观测点的数据值。半方差应随着距离的增大而增大。利用此曲线可以估算任何距离的半方差。利用这条拟合曲线可以模拟了区域变量的自相关变化。

半方差图中的拟合曲线通常采用以下五种数学曲线模型之一：球形曲线模型（Spherical）、圆环曲线模型（Circular）、指数模型（Exponential）、高斯模型（Gaussian）和线性模型（Linear）。

当h很小的时候，半方差很小，这说明两个距离很近的点，它们的值很相近，因此它们在空间上高度相关。随着h的增大，半方差迅速增大，即数据点之间值的差异快速增大，相应的，它们之间的空间相关性迅速减小。当h增大到一个临界值时，半方差开始趋于平稳，也即当h超过这个临界值时，随着h的继续增大，半方差变化很小，拟合曲线几乎近乎水平，这表面两点之间的距离超过h临界值以后它们互不相关，它们之间的值没有任何联系。h的这一临界值称为相关范围（Range），以a表示，在空间插值中，它可以用于确定待计算点邻域的大小，以选择邻域内所有的与其相关的观测点来估算该点的值。

此外，当h=0时，半方差应当为0，从半方差图可以注意到，拟合曲线与y轴有一个交点，即当h=0时，根据拟合曲线估算的半方差是一个正值，并非0。这个半方差值为残差，是由区域化变量中的随机噪音引起的不相关方差，称为矿块方差（Nugget Variance）以c₀表示。h等于相关范围时的半方差与矿块方差之差，以c表示。

克里格法使用半方差图来估算插值过程中所需要的已知点权重值，以使插值结果的方差达到最小，并运用半方差图估算出插值结果的方差值。根据对区域化变量分布特性的假设不同，克里格法可划分为好几类，GIS中常用的有两类：普通克里格法和广义克里格法。

除了上述两种常见的克里格插值法以外，还有块克里格法（BlockKriging）、普通克里格法（Co_Kriging）等。块克里格法用于估算一个区域化变量在一个小范围区域内（如一个田块）的平均值，而非某一点的值。普通克里格是在内插某一个变量值的过程中考虑另一个相关变量的值。往往在一个观测点可以获取多种变量值，如果两个变量在空间上是相关的话，有关其中一个变量空间变化方面的信息可以用来帮助内插另一个变量在空间各点的值，或改善对另一个变量的内查结果。例如，使用补充克里格法内插雨量数据时，可以将高程数据考虑进去。

图5.25 基于普通克里格生成的福建省平均气温趋势图

图5.26广义克里格法一次拟合和二次拟合的温度值内插图

5.密度估算

密度估算（Density Estimation）是根据点状实体的分布估算它们在一个区域内各地的分布密度。上述的四种空间插值方法用于估算连续性面状实体在每一点的值，而这里的密度估算则针对点状实体。例如:根据人口的点状分布，估算人口密度；根据动物的点状分布，估算动物的密度。密度估算产生一个面，可以栅格数据或等值线的形式输出。

最简单的密度估算方法是将一个格网覆盖在一个点状实体（如人口）分布图上，计算位于每个网格的点数，根据每个点所代表的点状实体的个数（如一个点代表10个人）计算每个网格所包含的点状实体的总数，将这个总数除以网格的面积即得一个网格的密度。

密度估算的另一种方法是采用核函数（Kernel Function）计算。核函数k（）在一个三维空间的形状类似于一个隆起的圆形沙丘，其底部圆形半径称为带宽。在密度估算过程中，可以假想将核的中心（底部圆心）移到待计算点上，以带宽为搜索半径，寻找核所包含的所有观测点，根据各观测点到待计算点的距离由核函数估算它们对待计算点影响的权重，然后用(5.17)式估算待计算点的密度：

使用上述函数可估算的空间物体分布的密度，例如课程实验项目3中的鹿密度估算。

在GIS中，通过空间插值获取的连续分布数据可以是TIN、栅格数据模型(Grid)或等高(值)线（Isoline）表示，空间插值技术也可用于TIN和栅格数据模型之间的相互转换以及栅格数据精度的变换。

值得强调的是：所有的空间插值技术都有一个基本的假设，即一个点的值受其附近已知点值的影响大于远距离点值的影响。例如，某一地点的雨量很可能接近于离它最近的气象站观测的雨量，与离它较远的气象站测得的数据相差较大。空间插值技术有很多种，在GIS中常用的五种是：趋势面分析、按距离加权、样条函数、克里格法和密度估算。但不管使用哪一种方法，通过插值获取的数据都是估算的近似数据，它们与实际值会有一定的差距，任何根据内插数据进行分析的结果都具有某种程度上的不确定性，因此，在使用内插数据时必须认识到数据的局限性。