空间分析.doc

第五章 空间分析 空间分析是对分析空间数据有关技术的统称。 根据作用的数据性质不同，可以分为 （1）基于空间图形数据的分析运算； （2）基于非空间属性的数据运算； （3）空间和非空间数据的联合运算。 一、 空间查询与量算简介 查询和定位空间对象，并对空间对象进行量算是地理信息系统的基本功能之一，它是地理信息系统进行高层次分析的基础。 例如在地理信息系统中，为进行高层次分析，往往需要查询定位空间对象，并用一些简单的量测值对地理分布或现象进行描述，如长度，面积，距离，形状等。实际上，空间分析首先始于空间查询和量算，它是空间分析的定量基础。 1． 空间查询 主要有两类 第一类是按属性信息的要求来查询定位空间位置，称为“属性查图形”。 例如在中国行政区划图上查询人口大于4000万且城市人口大于1000万的省有哪些 这和一般非空间的关系数据库的SQL查询没有区别，查询到结果后，再利用图形和属性的对应关系，进一步在图上用指定的显示方式将结果定位绘出。 第二类是根据对象的空间位置查询有关属性信息，称为“图形查属性”。 例如一般地理信息系统软件都提供一个“INFO”工具，让用户利用光标，用点选、画线、矩形、圆、不规则多边形等工具选中地物，并显示出所查询对象的属性列表，可进行有关统计分析。 该查询通常分为两步，首先借助空间索引，在地理信息系统数据库中快速检索出被选空间实体，然后根据空间实体与属性的连接关系即可得到所查询空间实体的属性列表 。 基于空间关系查询 空间实体间存在着多种空间关系，包括拓扑、顺序、距离、方位等关系。通过空间关系查询和定位空间实体是地理信息系统不同于一般数据库系统的功能之一。 例如查询满足下列条件的城市 – 在京沪线的东部 – 距离京沪线不超过50公里 – 城市人口大于100万 – 城市选择区域是特定的多边形； 整个查询计算涉及了空间顺序方位关系（京沪线东部），空间距离关系（距离京沪线不超过50公里），空间拓扑关系（使选择区域是特定的多边形），甚至还有属性信息查询（城市人口大于100万）。 简单的面、线、点相互关系的查询包括 面面查询，如与某个多边形相邻的多边形有哪些。 面线查询，如某个多边形的边界有哪些线。 面点查询，如某个多边形内有哪些点状地物。 线面查询，如某条线经过（穿过）的多边形有哪些，某条链的左、右多边形是哪些。 线线查询，如与某条河流相连的支流有哪些，某条道路跨过哪些河流。 线点查询，如某条道路上有哪些桥梁，某条输电线上有哪些变电站。 点面查询，如某个点落在哪个多边形内。 点线查询，如某个结点由哪些线相交而成。 2．空间量算 2．1几何量算 几何量算对不同的点、线、面地物有不同的含义 点状地物（0维）坐标； 线状地物（1维）长度，曲率，方向； 面状地物（2维）面积，周长，形状，曲率等； 体状地物（3维）体积，表面积等。 一般的GIS软件都具有对点、线、面状地物的几何量算功能，或者是针对矢量数据结构，或者是针对栅格数据结构的空间数据 。 1）线的长度计算 a 在矢量数据结构下 线表示为点对坐标（X，Y）或（X，Y，Z）的序列，在不考虑比例尺情况下，线长度的计算公式为 b 在栅格数据结构下 线状地物的长度就是累加地物骨架线通过的格网数目，骨架线通常采用8方向连接，当连接方向为对角线方向时，还要乘上 2）面状地物的面积 a 在矢量结构下面状地物以其轮廓边界弧段构成的多边形表示的。对于没有空洞的简单多边形，假设有N个顶点，其面积计算公式为 注意所采用的是几何交叉处理方法，即沿多边形的每个顶点作垂直与X轴的垂线，然后计算每条边、它的两条垂线及这两条垂线所截得X轴部分所包围的面积，所求出的面积的代数和，即为多边形面积 。对于有孔或内岛的多边形，可分别计算外多边形与内岛面积，其差值为原多边形面积。此方法亦适合于体积的计算。 b 栅格结构，多边形面积计算就是统计具有相同属性值的格网数目。 形状量算 线状地物形状量算 测量弯曲度 波状起伏的河流（23.5km） 直线距离18.5km 直线距离 起伏距离 ＝23.5/18.5＝1.27 面状地物形状（多边形的形状）量测 两个基本考虑 1）空间一致性问题，即有孔多边形和破碎多边形的处理； 2）多边形边界特征描述问题。 1.度量空间一致性最常用的指标欧拉函数 欧拉函数用来计算多边形的破碎程度和孔的数目。欧拉函数的结果是一个数，称为欧拉数。欧拉函数的计算公式为 欧拉数（孔数）-（碎片数-1） 对于图a，欧拉数4-1-14或欧拉数4-04；对于图b欧拉数4-2-13或欧拉数4-13；图c欧拉数5-3-13。 2.描述多边形边界特征的指标 关于多边形边界描述的问题，由于面状地物的外观是复杂多变的，很难找到一个准确的指标进行描述。 最常用的指标包括 a多边形长、短轴之比； b周长面积比； c面积长度比等。其中绝大多数指标是基于面积和周长的； d.通常认为圆形地物既非紧凑型也非膨胀型，则可定义其形状系数r为，其中P为地物周长，A为面积。如果r1为膨胀型。 质心量算 质心定义 质心是描述地理对象空间分布的一个重要指标。质心通常定义为一个多边形或面的几何中心。 例如要得到一个全国的人口分布等值线图，而人口数据只能到县级，所以必须在每个县域里定义一个点作为质心，代表该县的数值，然后进行插值计算全国人口等值线。 其中，Wi为第i个离散目标物权重，Xi，Yi为第i个离散目标物的坐标1 5 0.5 加权平均中心 平均中心或重心 几何中心 2 1.1 1.2 1 1.2 质心量测经常用于宏观经济分析和市场区位选择，还可以跟踪某些地理分布的变化，如人口变迁，土地类型变化等。 距离量算 “距离”是人们日常生活中经常涉及到的概念，它描述了两个事物或实体之间的远近程度。 1. 欧氏距离最常用的距离概念是欧氏距离，无论是矢量结构，还是栅格结构都很容易实现。 2.耗费距离 考虑到阻力影响，计算的距离称为耗费距离。物质在空间中移动总要花费一些代价，如资金、时间等。阻力越大耗费也越大。 3.距离表面在GIS中，距离通常是两个地点之间的计算，但有时人们想知道一个地点到所有其它地点的距离，这时得到的距离是一个距离表面 a.各向同性区域如果一区域中所有的性质与方向无关，则称为各向同性区域。以旅行时间为例，如果从某一点出发，到另一点的所耗费的时间只与两点之间的欧氏距离成正比，则从一固定点出发，旅行特定时间后所能达到的点必然组成一个等时圆。 b.各向异性距离表面而现实生活中，旅行所耗费的时间不只与欧氏距离成正比，还与路况、运输工具性能等有关，从固定点出发，旅行特定时间后所能到达的点则在各个方向上是不同距离的，形成各向异性距离表面。 通过耗费距离得到的距离表面称为阻力表面或耗费表面，其属性值代表一耗费或阻力大小。可以根据阻力表面计算最小耗费距离 障碍物 相对障碍物 绝对障碍物 栅格结构对耗费距离的表示 相对障碍物 绝对障碍物 矢量结构对非欧氏距离的表示 对于描述点、线、面坐标的矢量结构，也有一系列的不同于欧氏距离的概念。 欧氏距离通常用于计算两点的直线距离 当有障碍或阻力存在时，两点之间的距离就不能用直线距离，计算非标准欧氏距离的一般公式为 当k2时，就是欧氏距离计算公式。 当k1时，得到的距离称为曼哈顿距离。 欧氏距离、曼哈顿距离和非欧氏距离的计算如图所示 二、缓冲区分析 缓冲区分析是解决邻近度问题的空间分析工具之一 。 邻近度（Proximity）描述了地理空间中两个地物距离相近的程度，其确定是空间分析的一个重要手段。 例如交通沿线或河流沿线的地物有其独特的重要性，公共设施（商场，邮局，银行，医院，车站，学校等）的服务半径，大型水库建设引起的搬迁，铁路，公路以及航运河道对其所穿过区域经济发展的重要性等，均是一个邻近度问题。 1. 缓冲区的定义 所谓缓冲区就是地理空间目标的一种影响范围或服务范围。从数学的角度看，缓冲区分析的基本思想是给定一个空间对象或集合，确定它们的邻域，邻域的大小由邻域半径R决定。因此对象Oi的缓冲区定义为 即对象 Oi 的半径为R的缓冲区Bi为距Oi的距离d小于R的全部点的集合。d一般是最小欧氏距离，但也可是其它定义的距离。 对于对象集合 其半径为R的缓冲区是各个对象缓冲区的并集，即 特殊形态的缓冲区 点对象有三角形，矩形和圈形等； 对于线对象有双侧对称，双侧不对称或单侧缓冲区； 对于面对象有内侧和外侧缓冲区。 这些适合不同应用要求的缓冲区，尽管形态特殊，但基本原理是一致的。 2. 基于矢量数据的缓冲区分析 2.1. 双线问题 缓冲区计算的基本问题是双线问题。 双线问题有很多另外的名称，如图形加粗，加宽线，中心线扩张等，它们指的都是相同的操作。 角分线法 双线问题最简单的方法是角分线法（简单平行线法）。算法是在轴线首尾点处，作轴线的垂线并按缓冲区半径R截出左右边线的起止点；在轴线的其它转折点上，用与该线所关联的前后两邻边距轴线的距离为R的两平行线的交点来生成缓冲区对应顶点。 角分线法的缺点 当缓冲区半径不变时，d随张角B的减小而增大，结果在尖角处双线之间的宽度遭到破坏。 凸角圆弧法 在轴线首尾点处，作轴线的垂线并按双线和缓冲区半径截出左右边线起止点； 在轴线其它转折点处，首先判断该点的凸凹性，在凸侧用圆弧弥合，在凹侧则用前后两邻边平行线的交点生成对应顶点。 这样外角以圆弧连接，内角直接连接，线段端点以半圆封闭。 折点凸凹性的自动判断 该算法非常重要的一环是折点凸凹性的自动判断。此问题可转化为两个矢量的叉积把相邻两个线段看成两个矢量，其方向取坐标点序方向。若前一个矢量以最小角度扫向第二个矢量时呈逆时针方向，则为凸顶点，反之为凹顶点。 具体算法过程如下 由矢量代数可知，矢量AB，BC可用其端点坐标差表示 矢量代数叉积遵循右手法则，如果拐点呈逆时针方向时 若S0，则ABC呈逆时针，顶点为凸； 若S1000mm ，土厚50cm 的地区。则表达式为 E （降雨量1000∩土厚50cm） （2）根据专家命题模型来确定条件 。例如根据农业专家经验知道适宜于种植水稻的条件为 积温3200 度 降雨800mm 坡度200 天 E 积温3200 ∩ 降雨800mm ∩ 坡度200 天 叠置方法 下面以游程编码为例说明叠置分析的实现 已知某地区同一比例尺的降雨量图和土层厚度图（以第K 行为例） 若对第K行按如下关系表达式进行条件叠置，其条件表达式为 E（降雨量＝1000）∩（土厚＝50） 游程号 游程属性 游程最右列 1 0 400 2 1000，50 680 3 0 800 网络分析 网络分析是运筹学模型中的一个基本模型，它的根本目的是研究、筹划一项网络工程如何安排，并使其运行效果最好。其基本思想则在于人类活动总是趋于按一定目标选择达到最佳效果的空间位置。 网络结构的评估指标 γ指数一个网络实际连接数目和最大可能连接数目的比值。 α指数网络中实际的回路数与最大回路数的比值。 网络直径一个连通的网络中任意结点到其他结点最短距离所需的最大步测数。 网络的连通性 一个网络的连通性可通过建立一个被称为C矩阵的矩阵集合来检验。 C1 表示第一阶连通矩阵，是结点之间的直接连接。简单说，如果两个结点直接连接则相应单元为1，否则为0。 C2 矩阵中的每一个元素表示，从一个结点到另一个结点跨两步的路线数目。 T C1 C2 C3 主要网络分析功能 路径分析 1）静态求最佳路径在给定每条链上的属性后，求最佳路径。 2）N条最佳路径分析确定起点或终点，求代价最小的N条路径，因为在实践中最佳路径的选择只是理想情况，由于种种因素而要选择近似最优路径。 3）最短路径或最低耗费路径确定起点、终点和要经过的中间点、中间连线，求最短路径或最小耗费路径。 4）动态最佳路径分析实际网络中权值是随权值关系式变化的，可能还会临时出现一些障碍点，需要动态的计算最佳路径 计算最短路径的Dijkstra算法 下面是该算法的描述。 1）用带权的邻接矩阵Cost来表示带权的n个节点的有向图，Cost[i,j]表示弧的权值，如果从vi到vj不连通，则Cost[i,j]∞。下图表示了一个带权有向图以及其邻接矩阵。 然后，引进一个辅助向量Dist，每个分量Dist[i]表示从起始点到每个终点vi的最短路径长度。假定起始点在有向图中的序号为i0，并设定该向量的初始值为 Dist[i]Cost[i0,i] vi∈V。 令S为已经找到的从起点出发的最短路径的终点的集合。 2）选择Vj，使得 Dist[j]Min{ Dist[i]|Vi∈V-S} vi∈V vj就是当前求得的一条从vi0出发的最短路径的终点，令 SS∪{vj} 3）修改从vi0出发到集合V-S中任意一顶点vk的最短路径长度。如果 Dist[j]Cost[j,k]CkfB或fACkfB.换言之，当fA-CK*fB-Ck0表示等值线和棱边AB相交 n b. 求出交点 n 2）起始点的搜索 n 3）等值线的跟踪 a. 跟踪状态分析 b. 方向的选择 n 4）终点的确定 4．不规则三角网DEM 曲面数据结构 n 曲面是指连续分布现象的覆盖表面，具有这种覆盖表面的要素有地形、降水量、温度、磁场等。表示和存储这些要素的基本要求是必须便于连续现象在任一点的内插计算，因此经常采用不规则三角网来拟合连续分布现象的覆盖表面，称为TIN（Triangulated Irregular Network）数据结构. 不规则三角网（TIN） 不规则三角网数字高程特点 n 不规则三角网数字高程由连续的三角面组成，三角面的形状和大小取决于不规则分布的测点，或节点的位置和密度。 不规则三角网DEM与格网DEM的区别 n 不规则网DEM可以随地形起伏变化的复杂性而改变采样点的密度和决定采样点的位置，因而它能够避免地形平坦时的数据冗余，又能按地形特征点如山脊、山谷线、地形变化线等表示数字高程特征。 不规则三角网DEM的建立步骤 算法要点1.基本依据三角形余弦定理 狄洛尼（Delaunay）三角网 n 三角网中每个三角形要求尽量接近等边形状，并保证由最邻近的点构成的三角形，即三角形的边长之和最小。 n 在所有可能的三角网中，狄洛尼（Delaunay）三角网在地形拟合方面表现最为出色。 狄洛尼（Delaunay）三角网定义 n 狄洛尼（Delaunay）三角网为相互邻接且互相不重叠的三角形的集合，每一个三角形的外接圆内不不含其他的点。 n 狄洛尼三角形外接圆不包含其他点的特性被用作从一系列不重合的平面点建立狄洛尼三角网的基本法则，可以称为狄洛尼法则 对于给定的初始点集P，有多种三角网剖分方式，而Delaunay三角网有以下特性 n 1）其Delaunay三角网是唯一的； n 2）三角网的外边界构成了点集P的凸多边形“外壳”； n 3）没有任何点在三角形的外接圆内部，反之，如果一个三角网满足此条件，那么它就是Delaunay三角网。 n 4）如果将三角网中的每个三角形的最小角进行升序排列，则Delaunay三角网的排列得到的数值最大，从这个意义上讲，Delaunay三角网是“最接近于规则化”的三角网。 Delaunay三角形产生的基本准则 n 1.Delaunay三角形产生准则的最简明的形式是任何一个Delaunay三角形的外接圆的内部不能包含其它任何点[Delaunay 1934]。 n 2。最大化最小角原则 Lawson[1972]提出了最大化最小角原则每两个相邻的三角形构成的凸四边形的对角线，在相互交换后，六个内角的最小角不再增大。 局部优化过程LOP（Local Optimization Procedure）（ Lawson [1977]） n 如图所示。先求出包含新插入点p的外接圆的三角形，这种三角形称为影响三角形（Influence Triangulation）。删除影响三角形的公共边（图b中粗线），将p与全部影响三角形的顶点连接，完成p点在原Delaunay三角形中的插入。 格网DEM转成TIN 格网DEM转成TIN可以看作是一种规则分布的采样点生成TIN的特例，其目的是尽量减少TIN的顶点数目，同时尽可能多地保留地形信息，如山峰、山脊、谷底和坡度突变处。规则格网DEM可以简单地生成一个精细的规则三角网，针对它有许多算法，绝大多数算法都有两个重要的特征 1）筛选要保留或丢弃的格网点； 2）判断停止筛选的条件。 其中两个代表性的方法算法是保留重要点法和启发丢弃法。 保留重要点法 [Chen、Gauvara（1987）]。 n 思想 保留规则格网DEM中的重要点，然后用其构造TIN。 重要点（VIP，Very Important Point）确定 n 如图所示，由3*3的模板得到中心点P和8邻点的高程值，计算中心点P到直线AE，CG，BF，DH的距离，再计算4个距离的平均值。如果平均值超过阈值，P点为重要点，则保留，否则去除P点。 等高线转成格网DEM n 方法使用局部插值算法 n 问题 “阶梯”地形。 数字化等高线时越小心，采样点越多，问题越严重。以带“阶梯”地形的DEM为基础，计算坡度往往会出现不自然的条斑状分布模式。 图等值线插值造成“阶梯地形”的原因 最好把等高线数据点减少到最少，增加标识山峰、山脊、谷底和坡度突变的数据点，同时使用一个较大的搜索窗口。 利用格网DEM提取等高线 n 1.根据格网DEM中相邻四个点组成四边形进行等高线跟踪。 n 2.将每个矩形分割成为两个三角形，并应用TIN提取等高线算法。 但是由于矩形有两种划分三角形的方法，在某些情况下，会生成不同的等高线（图），这时需要根据周围的情况进行判断并决定取舍。 产生问题1 n 使用四边形跟踪等高线，也会出现等高线跟踪的二义性。 n 解决方法计算距离，距离近的连线方式优于距离远的连线方式。在图中，就要采用（b）图所示的跟踪方式。 产生问题2 n 如果网格点的数值恰好等于要提取的等高线的数值，会使判断过程变得复杂，并且会生成不闭合的等高线。 n 解决办法将这些网格点的数值增加一个小的偏移量。 TIN转成格网DEM n TIN转成格网DEM可以看作普通的不规则点生成格网DEM的过程。方法是按要求的分辨率大小和方向生成规则格网，对每一个格网搜索最近的TIN数据点，按线性或非线性插值函数计算格网点高程。 DEM应用 地形曲面拟合 n DEM最基础的应用是求DEM范围内任意点的高程，在此基础上进行地形属性分析 。 由于已知有限个格网点的高程，可以利用这些格网点高程拟合一个地形曲面，推求区域内任意点的高程 立体透视图 n 绘制透视立体图是DEM的一个极其重要的应用。透视立体图能更好地反映地形的立体形态，非常直观。 n 从一个空间三维的立体的数字高程模型到一个平面的二维透视图，其本质就是一个透视变换。将“视点”看作为“摄影中心”，可以直接应用共线方程从物点（X，Y，Z）计算“像点”坐标（X,Y）。透视图中的另一个问题是“消隐”的问题，即处理前景挡后景的问题。 n 调整视点、视角等各个参数值，就可从不同方位、不同距离绘制形态各不相同的透视图制作动画。计算机速度充分高时，就可实时地产生动画DTM透视图。 通视分析 通视问题可以分为五类 n 1）已知一个或一组观察点，找出某一地形的可见区域。 n 2）欲观察到某一区域的全部地形表面，计算最少观察点数量。 n 3）在观察点数量一定的前提下，计算能获得的最大观察区域。 n 4）以最小代价建造观察塔，要求全部区域可见。 n 5）在给定建造代价的前提下，求最大可见区。 通视可分为点的通视，线的通视和面的通视。 1.点的通视指计算视点与待判定点之间的可见性问题； 2.线的通视指已知视点，计算视点的视野问题； 3.区域的通视指已知视点，计算视点能可视的地形表面区域集合的问题。 剖面分析 剖面积 n 根据工程设计的线路，可计算其与DEM各格网边交点Pi（Xi，Yi，Zi），则线路剖面积为 其中n为交点数；Di，i1为Pi与P i1之距离。同理可计算任意横断面及其面积。 体积 n 其中S3与S4分别是三棱柱与四棱柱的底面积。 n 根据两个DEM可计算工程中的挖方、填方及土壤流失量。 表面积 n 对于含有特征的格网，将其分解成三角形，对于无特征的格网，可由4个角点的高程取平均即中心点高程，然后将格网分成4个三角形。由每一三角形的三个角点坐标（xi，yi，zi）计算出通过该三个顶点的斜面内三角形的面积，最后累加就得到了实地的表面积。 坡度和坡向的计算 图 提取坡度 图 提取坡向 其他的地形分析 n 表面曲率 n 流域分析 空间插值 n 点插值（全局方法和局部方法） 一、全局方法 n 趋势面分析 1.线性插值 2.双线性多项式插值 3.双三次多项式（样条函数）插值 n 回归模型 二、局部方法（泰森多边形、密度估算、反距离权重插值、薄板样条函数法、克里金法） n 区域插值 一.叠合法 二、比重法