第六章 GIS空间分析原理与方法
第一节 GIS空间分析模型
地学模型概述
地理信息系统以数字世界表示自然世界,具有完备的空间特性,可以存储和处理不同地理发展时期的大量地理数据、并具有极强的空间系统综合分析能力,是地理分析的有力工具。因此,地理信息系统不仅要完成管理大量复杂的地理数据的任务,更为重要的是要完成地理分析、评价、预测和辅助决策的任务,必须发展广泛的适用于地理信息系统的地理分析模型,这是地理信息系统走向实用的关键。
所谓模型,就是将系统的各个要素,通过适当的筛选,用一定的表现规则所描写出来的简明映像。模型通常表达了某个系统的发展过程或发展结果。地学模型是用来描述地理系统各地学要素之间的相互关系和客观规律信息的语言的或数学的或其他表达形式,通常反映了地学过程及其发展趋势或结果。地学模型也称为专题分析模型。对于地理信息系统来说,专题分析模型是根据关于目标的知识将系统数据重新组织,得出与目标有关的更为有序的新的数据集合的有关规则和公式。这是应用地理信息系统进行生产和科研的重要手段。模型化是将主观性的思考,以模型的形式反映出来,不同的理论观点,不同的体系可以产生不同的结果。
地学分析模型主要包含以下几种形式:1、逻辑模型:由地理名词和逻辑运算符组成的逻辑表达式表示;2、物理模型:由物理模拟过程表达;3、数学模型:由常数、参数、变量和函数关系等组成的数学表达式表示;4、图像模型:由某种图像或图像运算的集合表达,如各种专题地图。
专题分析模型在地理信息系统中的作用表现在以下几个方面:
地理信息系统的设计
任何地理信息系统都是为一定的应用目的而建立的,必须根据具体需要采用适用的分析模型指导地理信息系统总体设计。主要包括:①数据项的选择,数据的范围、精度、量测方法等,如果毫无选择地录入数据,只会使系统增加负担,降低效率,无法突出主要因素,甚至因为数据采集周期过长而失去意义;数据结构应以最好地表示地理现象和易于模型实现为标准;②硬件环境的选择,根据模型的输入、输出和运算方法选择经济实用的硬件支持;③软件功能的选择,根据模型的管理和运行设计适用的软件功能。
地理信息系统的应用
目前地理信息系统技术的推广应用遇到了三个方面的困难:①硬件环境特殊,不易配备;②地理信息系统知识没有为许多用户掌握;③缺乏足够的专题分析模型。而最重要的因素在于地理信息系统是否具有实用价值,实用性则必须依靠正确地应用专题分析模型。
利于信息交流
模型是表达思维对自然界认识的工具,因此地理信息系统的各种分析模型则有利于完整准确地表达使用者对问题的认识和处理方法,既利于使用者与系统设计者之间的交流以发展系统功能,又利于使用者之间交流以增强系统的共享性。
二、地理信息系统模型化的一般方法
模型的建立过程可由下式表示:
XOY=M
其中X表示某个体系,可以看作是地理系统中被主观选取的一个局部;Y表示某种介体,具体讲就是某种模型化方法;O代表Y对X产生的作用;M是体系X通过介体Y产生的作用O所建立的模型。
分析模型按建立的方法主要有以下三种类型:
1、概念模型:又称为逻辑模型,主要指通过观察、总结、提炼而得到的文字描述或逻辑表达式,常由此构成专家系统的知识库;
2、数学模型:又称为理论模型,是应用数学分析方法建立的数学表达式,反映地理过程本质的物理规律,如区位模型(Location Models)就是解决地理空间问题的很有价值的理论模型;
3、统计模型:包括经验模型,是通过数理统计方法和大量观察实验得到的定量模型,具有简单实用的优点。
通常需综合各种方法:概念模型比较灵活,可以引入许多模糊概念,适用范围很广,易于为多数人接受,但难以进行精确定量分析;数学模型因果关系清楚,可以精确地反映系统内各要素之间的定量关系,易于用来对自然过程施加控制,但通常难以包括太多的要素,而常常是大大简化的理想情形,削弱了其实用性;统计模型可以通过大量的实践建立,具有简单实用、适用性广、可以处理大量相关因素的特点,缺点是过程不清,一般是采用“黑箱”或“灰箱”方法建立的。
作为一般规则,首先应在实践中不断观察总结,形成愈来愈丰富的概念模型,在积累经验的基础上采用数理统计方法摸索统计规律;最后上升到理论模型;再采用综合方法建立实用的分析模型。
运用综合方法建立地理信息系统分析模型可采用以下步骤: ①系统描述与数据分析:对模型所要分析的系统,选择可以描述系统的状态、与外部关系、随时间变化等方面的数据,构造该系统的数据体系。②理论推导:根据地理规律和系统的特点,进行理论推导,确定上面的数据体系中多因子之间的量纲关系,作为分析模型的基本框架;③简化表达:根据理论分析和具体应用要求,筛选去除相对影响较小和不重要的要素,或采用主成分分析法等数学方法简化表达形式,使模型接近实用;④参数确定:模型参数的确定可采用参数试验方法,或采用层次分析法(AHP)、专家打分法、确定模糊隶属度等方法。形式和参数确定后,分析模型可在应用中完善。
由于理论和实践方面的原因,有时可采用递归模型。递归模型便于导出地理系统在任一演变时期的状态和演变过程,在较短的间隔周期内可以作为线性问题处理,并且可以参照假设条件的变化随时间调整模型参数。
第二节 栅格数据分析的基本模式
栅格数据由于其自身数据结构的特点,在数据处理与分析中通常使用线性代数的二维数字矩阵分析法作为数据分析的数学基础。因此,具有自动分析处理较为简单,而且分析处理模式化很强的特征。一般来说,栅格数据的分析处理方法可以概括为聚类聚合分析、多层面复合叠置分析、窗口分析及追踪分析等几种基本的分析模型类型。以下分别进行描述与讨论。
一、栅格数据的聚类、聚合分析
栅格数据的聚类、聚合分析均是指将一个单一层面的栅格数据系统经某种变换而得到一个具有新含义的栅格数据系统的数据处理过程。也有人将这种分析方法称之为栅格数据的单层面派生处理法。
1、聚类分析
栅格数据的聚类是根据设定的聚类条件对原有数据系统进行有选择的信息提取而建立新的栅格数据系统的方法。
图6—1(a)为一个栅格数据系统样图,1、2、3、4为其中的四种类型要素, 图6—1(b)为提取其中要素“2”的聚类结果。
2、聚合分析
栅格数据的聚合分析是指根据空间分辨力和分类表,进行数据类型的合并或转换以实现空间地域的兼并。
空间聚合的结果往往将较复杂的类别转换为较简单的类别,并且常以较小比例尺的图形输出。当从地点、地区到大区域的制图综合变换时常需要使用这种分析处理方法。对于图6—1(a),如给定聚合的标准为1、2类合并为b,3、4类合并为a,则聚合后形成的栅格数据系统如图6—2(a)所示, 如给定聚合的标准为2、3类合并为c,1、4类合并为d,则聚合后形成的栅格数据系统如图6—2(b)所示。
栅格数据的聚类聚合分析处理法在数字地形模型及遥感图象处理中的应用是十分普遍的。例如,由数字高程模型转换为数字高程分级模型便是空间数据的聚合,而从遥感数字图象信息中提取其一地物的方法则是栅格数据的聚类。
二、栅格数据的信息复合分析
能够极为便利地进行同地区多层面空间信息的自动复合叠置分析,是栅格数据一个最为突出的优点。正因为如此,栅格数据常被用来进行区域适应性评价、资源开发利用、规划等多因素分析研究工作。在数字遥感图象处理工作中,利用该方法可以实现不同波段遥感信息的自动合成处理;还可以利用不同时间的数据信息进行某类现象动态变化的分析和预测。因此该方法在计算机地学制图与分析中具有重要的意义。信息复合模型(overlay)包括两类,即简单的视觉信息复合和较为复杂的叠加分类模型。
1、视觉信息复合
视觉信息复合是将不同专题的内容叠加显示在结果图件上,以便系统使用者判断不同专题地理实体的相互空间关系,获得更为丰富的信息。地理信息系统中视觉信息复合包括以下几类:
(1) 面状图、线状图和点状图之间的复合;
(2) 面状图区域边界之间或一个面状图与其他专题区域边界之间的复合;
(3) 遥感影像与专题地图的复合;
(4) 专题地图与数字高程模型复合显示立体专题图;
(5) 遥感影像与DEM复合生成真三维地物景观。
2、叠加分类模型
简单视觉信息复合之后,参加复合的平面之间没发生任何逻辑关系,仍保留原来的数据结构;叠加分类模型则根据参加复合的数据平面各类别的空间关系重新划分空间区域,使每个空间区域内各空间点的属性组合一致。叠加结果生成新的数据平面,该平面图形数据记录了重新划分的区域,而属性数据库结构中则包含了原来的几个参加复合的数据平面的属性数据库中所有的数据项。叠加分类模型用于多要素综合分类以划分最小地理景观单元,进一步可进行综合评价以确定各景观单元的等级序列。
以下按复合运算方法的不同进行分类讨论。
(1)逻辑判断复合法
设有A、B、C三个层面的栅格数据系统,一般可以用布尔逻辑算子以及运算结果的文氏图(见图6-3)表示其一般的运算思路和关系。
(2)数学运算复合法
是指不同层面的栅格数据逐网格按一定的数学法则进行运算,从而得到新的栅格数据系统的方法。其主要类型有以下几种:
1)算术运算
指两层以上的对应网格值经加、减运算,而得到新的栅格数据系统的方法。这种复合分析法具有很大的应用范围。图6-4给出了该方法在栅格数据编辑中的应用例证。
2)函数运算
指两个以上层面的栅格数据系统以某种函数关系作为复合分析的依据进行逐网格运算,从而得到新的栅格数据系统的过程。
这种复合叠置分析方法被广泛地应用到地学综合分析、环境质量评价、遥感数字图像处理等领域中。
例如利用土壤侵蚀通用方程式计算土壤侵蚀量时,就可利用多层面栅格数据的函数运算复合分析法进行自动处理。一个地区土壤侵蚀量的大小是降雨(R)、植被覆度(C)、坡度(S)、坡长(L)、土壤抗蚀性(SR)等因素的函数。可写成
E=F(R,C,S,L,SR…)
逐网格的复合分析运算如图6-5所示。
类似这种分析方法在地学综合分析中具有十分广泛的应用前景。只要得到对于某项事物关系及发展变化的函数关系式,便可运用以上方法完成各种人工难以完成的极其复杂的分析运算。这也是目前信息自动复合叠置分析法受到广泛应用的原因。
值得注意是,信息的复合法只是处理地学信息的一种手段,而其中各层面信息关系模式的建立对分析工作的完成及分析质量的优劣具有决定性作用。这往往需要经过大量的试验研究,而计算机自动复合分析法的出现也为获得这种关系模式创造了有利的条件。
三、栅格数据的追踪分析
所谓栅格数据的追踪分析是指,对于特定的栅格数据系统,由某一个或多个起点,按照一定的追踪线索进行追踪目标或者追踪轨迹信息提取的空间分析方法。如图6-6所示例,栅格所记录的是地面点的海拔高程值,根据地面水流必然向最大坡度方向流动的基本追踪线索,可以得出在以上两个点位地面水流的基本轨迹。此外,追踪分析法在扫描图件的矢量化、利用数字高程模型自动提取等高线、污染源的追踪分析等方面都发挥着十分重要的作用。
四、栅格数据的窗口分析
地学信息除了在不同层面的因素之间存在着一定的制约关系之外,还表现在空间上存在着一定的关联性。对于栅格数据所描述的某项地学要素,其中的(I,J)栅格往往会影响其周围栅格的属性特征。准确而有效地反映这种事物空间上联系的特点,也必然是计算机地学分析的重要任务。窗口分析是指对于栅格数据系统中的一个、多个栅格点或全部数据,开辟一个有固定分析半径的分析窗口,并在该窗口内进行诸如极值、均值等一系列统计计算,或与其它层面的信息进行必要的复合分析,从而实现栅格数据有效的水平方向扩展分析。
分析窗口的类型
按照分析窗口的形状,可以将分析窗口划分为以下类型:
(1)矩形窗口:是以目标栅格为中心,分别向周围八个方向扩展一层或多层栅格,从而形成如图6-7所示的矩形分析区域。
(2)圆型窗口:是以目标栅格为中心,向周围作一等距离搜索区,构成一圆型分析窗口,见图6-7。
(3)环型窗口:是以目标栅格为中心,按指定的内外半径构成环型分析窗口,见图6-7。
(4)扇型窗口:是以目标栅格为起点,按指定的起始与终止角度构成扇型分析窗口,见图6-7。
窗口内统计分析的类型
栅格分析窗口内的空间数据的统计分析类型一般有以下几种类型:
(1)Mean;(2)Maximum; (3)Minimum;(4)Median;(5)Sum;
(6)Range;(7)Majority;(8)Minority;(9)Variety
3×3窗口 5×5窗口 7×7窗口
圆形窗口 环形窗口 扇形窗口
图6-7 分析窗口的类型
在实际工作中,为解决某一个具体的应用命题,以上4种栅格数据的分析模式往往综合使用。
第三节 矢量数据分析的基本方法
与栅格数据分析处理方法相比, 矢量数据一般不存在模式化的分析处理方法, 而表现为处理方法的多样性与复杂性。本节选择几种最为常见的几何分析法, 并以其为例, 说明矢量数据分析处理的基本原理与方法。
包含分析
确定要素之间是否存在着直接的联系,即矢量点、线、面之间是否存在在空间位置上的联系,这是地理信息分析处理中常要提出的问题,也是在地理信息系统中实现图形—属性对位检索的前提条件与基本的分析方法。例如,若在计算机屏幕上利用鼠标点击对应的点状、线状或面状图形,查询其对应的属性信息;或需要确定点状居民地与线状河流或面状地类之间的空间关系(如是否相邻或包含),都需要利用矢量数据的包含分析与数据处理方法。例如,要确定某个井位属于哪个行政区;要测定某条断裂线经过哪些城市建筑,都需要通过GIS信息分析方法中对已有矢量数据的包含分析来实现以上目标。
在包含分析的具体算法中,点与点、点与线的包含分析一般均可以分别通过先计算点到点,点到线之间的距离,然后,利用最小距离阈值判断包含的结果。点与面之间的包含分析,或称为Point-Polygon分析,具有较为典型的意义。可以通过著名的铅垂线算法来解决,如图6-9所示,由Pt点作一条铅垂线。现在要测试Pt是在该多边形之内或之外。 其基本算法的思路是,如果该铅垂线与某一图斑有奇数交点,则该Pt点必位于该图斑内(某些特殊条件除外)。
利用这种包含分析方法,还可以解决地图的自动分色,地图内容从面向点的制图综合,面状数据从矢量向栅格格式的转换,以及区域内容的自动计数(例如某个设定的森林砍伐区内,某一树种的颗数) 等等。例如,确定某区域内矿井的个数,这是点与面之间的包含分析,确定某一县境内公路的类型以及不同级别道路的里程,是线与面之间的包含分析。分析的方法是:首先对这些矿井、公路要点、线要素数字化,经处理后形成具有拓扑关系的相应图层,然后和已经存放在系统中的多边形进行点与面、线与面的叠加;最后对这个多边形或区域进行这些点或线段的自动计数或归属判断。
二、矢量数据的缓冲区分析
缓冲区分析是研究根据数据库的点、线、面实体,自动建立其周围一定宽度范围内的缓冲区多边形实体,从而实现空间数据在水平方向得以扩展的信息分析方法。它是地理信息系统重要的和基本的空间操作功能之一。例如,城市的噪音污染源所影响的一定空间范围、交通线两侧所划定的绿化带,即可分别描述为点的缓冲区与线的缓冲带。而多边形面域的缓冲带有正缓冲区与负缓冲区之分(如图6-10所示)。
三、多边形叠置分析
多边形叠置分析也称为Polygon-on-polygon 叠置,它是指同一地区、同一比例尺的两组或两组以上的多边形要素的数据文件进行叠置。参加叠置分析的两个图层应都是矢量数据结构。若需进行多层叠置, 也是两两叠置后再与第三层叠置,依次类推。其中被叠置的多边形为本底多边形,用来叠置的多边形为上覆多边形,叠置后产生具有多重属性的新多边形。
其基本的处理方法是,根据两组多边形边界的交点来建立具有多重属性的多边形或进行多边形范围内的属性特性的统计分析。其中,前者叫做地图内容的合成叠置(图6-11),后者称为地图内容的统计叠置(图6-12)。
合成叠置的目的,是通过区域多重属性的模拟,寻找和确定同时具有几种地理属性的分布区域。或者按照确定的地理指标,对叠置后产生的具有不同属性的多边形进行重新分类或分级,因此叠置的结果为新的多边形数据文件。统计叠置的目的,是准确地计算一种要素(如土地利用)在另一种要素(如行政区域)的某个区域多边形范围内的分布状况和数量特征(包括拥有的类型数、各类型的面积及所占总面积的百分比等等),或提取某个区域范围内某种专题内容的数据。
多边形叠置方法在国内外已有发展并得到较为广泛的应用,如ARC/INFO地理信息系统中,多边形叠置是该系统的关键性软件,英国运用多边形叠置技术进行了土地适宜性评价,此外,国际上已建立起来的地理信息系统中,有许多具备了多边形叠置分析的功能。
四、矢量数据的网络分析
1、基本概念
网络分析的主要用途是:选择最佳路径;选择最佳布局中心的位置。所谓最佳路径是 指从始点到终点的最短距离或花费最少的路线(图6-13);最佳布局中心位置是指各中心 所覆盖范围内任一点到中心的距离最近或花费最小;网流量是指网络上从起点到终点的某个函数,如运输价格,运输时间等。网络上任意点都可以是起点或终点。其基本思想则在于人类活动总是趋向于按一定目标选择达到最佳效果的空间位置。这类问题在生产、社会、经济活动中不胜枚举,因此研究此类问题具有重大意义。
网络中的基本组成部分和属性如下:
(1)链(Link),网络中流动的管线,如街道,河流,水管等,其状态属性包括阻力和需求。
(2)障碍,禁止网络中链上流动的点。
(3)拐角点,出现在网络链中所有的分割结点上状态属性的阻力,如拐弯的时间和限制(如不允许左拐)。
(4)中心,是接受或分配资源的位置,如水库、商业中心、电站等。其状态属性包括资源容量,如总的资源量;阻力限额,如中心与链之间的最大距离或时间限制。
(5)站点,在路径选择中资源增减的站点,如库房、汽车站等,其状态属性有要被运输的资源需求,如产品数。
网络中的状态属性有阻力和需求两项,实际的状态属性可通过空间属性和状态属性的转换,根据实际情况赋到网络属性表中。
2、网络分析的基本方法
(1) 路径分析
a. 静态求最佳路径:由用户确定权值关系后,即给定每条弧段的属性,当需求最佳路径时,读出路径的相关属性,求最佳路径。
b. 动态分段技术:给定一条路径由多段联系组成,要求标注出这条路上的公里点或要求定位某一公路上的某一点,标注出某条路上从某一公里数到另一公里数的路段。
c. N条最佳路径分析:确定起点、终点,求代价较小的几条路径,因为在实践中往往仅求出最佳路径并不能满足要求,可能因为某种因素不走最佳路径,而走近似最佳路径。
d. 最短路径:确定起点、终点和所要经过的中间点、中间连线,求最短路径。
e. 动态最佳路径分析:实际网络分析中权值是随着权值关系式变化的,而且可能会临时出现一些障碍点,所以往往需要动态地计算最佳路径。
(2) 地址匹配
地址匹配实质是对地理位置的查询,它涉及到地址的编码。地址匹配与其它网络分析功能结合起来,可以满足实际工作中非常复杂的分析要求。所需输入的数据,包括地址表和含地址范围的街道网络及待查询地址的属性值。
(3)资源分配
资源分配网络模型由中心点(分配中心)及其状态属性和网络组成。分配有两种方式,一种是由分配中心向四周输出,另一种是由四周向中心集中。这种分配功能可以解决资源的有效流动和合理分配。其在地理网络中的应用与区位论中的中心地理论类似。在资源分配模型中,研究区可以是机能区,根据网络流的阻力等来研究中心的吸引区,为网络中的每一连接寻找最近的中心,以实现最佳的服务。还可以用来指定可能的区域。
资源分配模型可用来计算中心地的等时区、等交通距离区、等费用距离区等。可用来进行城镇中心、商业中心或港口等地的吸引范围分析,以用来寻找区域中最近的商业中心,进行各种区划和港口腹地的模拟等。
第三节 空间数据的其它分析方法
一、空间数据的量算
空间信息的自动化量算是地理信息系统所具有的重要功能,也是进行空间分析的定量化基础。其中的主要量算有:
质心量算
描述地理目标空间分布的最有用的单一量算量是目标的质心位置。地理目标的质心是目标的半径位置。它的目标保持均匀分布的平衡点,它可通过对目标坐标值加权平均求得。
其中,i为离散目标物,Wi为该目标权重,Xi、Yi 为目标。
质心的量算,可以跟踪某些地理分布的变化,例如人口的变迁、土地类型的变化,也可以简化某些复杂目标, 在某些情况下,可以方便的导出某些预测模型。
2、几何量算
几何量算对点、线、面、体4类目标物而言,其含义是不同的:
·点状目标:坐标;
·线状目标:长度、曲率、方向;
·面状目标:面积、周长等;
·体状目标:表面积、体积等。
线由点组成,而线长度可由两点间直线距离相加得到。
面积和周长的计算。在平面直角坐标系中,计算面积时,计算y值以下面积.按矢量方向,分别求出向右向左两个方向各自的面积,它们的绝对值之差,便是多边形面积值,周长则是线段之和。
3、形状量算
目标物的外观是多变的,很难找到一个准确的量对其进行描述。因此,对目标属紧凑型的或膨胀型的判断极其模糊。
如果认为一个标准的圆目标既非紧凑型也非膨胀型,则可定义其形状系数为:
其中,P为目标物周长,A为目标物面积。
如果 r<1,目标物为紧凑型;
r=1,目标物为一标准圆;
r>1,目标物为膨胀型。
二、空间数据的内插
空间数据往往是根据自己要求所获取的采样观测值,诸如土地类型、地面高程等。这些点的分布往往是不规则的,在用户感兴趣或模型复杂区域可能采样点多,反之则少。由此而导致所形成的多边形的内部变化不可能表达得更精确、更具体,而只能达到一般的平均水平或“象征水平”。但用户在某些时候却欲获知未观测点的某种感兴趣特征的更精确值,这就导致了空间内插技术的诞生。一般来讲,在已存在观测点的区域范围之内估计未观测点的特征值的过程称内插;在已存在观测点的区域范围之外估计未观测点的特征值的过程称推估。
前面已经提到,现实空间可以分为具有渐变特征的连续空间和具有跳跃特征的离散空间。举例来讲,土地类型分布属离散空间,而地形表面分布则是连续空间,见图6-14。
对于离散空间,假定任何重要变化发生在边界上,如bc段上方为土地类型B,下方则为类型C,其边界内的变化则是均匀的,同质的,即在各个方面都是相同的。对于这种空间的最佳内插方法是邻近元法,即以最邻近图元的特征值表征未知图元的特征值。这种方法在边界会产生一定的误差,但在处理大面积多边形时,则十分方便。但是,对于连续空间表面,上述处理方法则不合适。连续表面的内插技术必须采用连续的空间渐变模型实现这些连续变化,可用一种平滑的数学表面加以描述。这类技术可分为整体拟合和局部拟合技术两大类。整体拟合技术即拟合模型是由研究区域内所有采样点上的全部特征观测值建立的。通常采用的技术是整体趋势面拟合。这种内插技术的特点是不能提供内插区域的局部特性,因此,该模型一般用于模拟大范围内的变化。而局部拟合技术则是仅仅用邻近的数据点来估计未知点的值,因此可以提供局部区域的内插值,而不致受局部范围外其它点的影响。这类技术包括:双线性多项式内插、移动拟合法、最小二乘配置法等等。
1、趋势面拟合技术
多项式回归分析是描述大范围空间渐变特征的最简单方法。多项式回归的基本思想是用多项式表示线(数据是一维时)或面(数据是二维时),并按最小二乘法原理对数据点进行的拟合,拟合时假定数据点的空间坐标X、Y为独立变量,而表征特征值的Z坐标为因变量。当数据为一维(X)时,这种变化可用回归线近似表示为:
Z=a0十a1X
式中,a0、a1为多项式系数。当n个采样点上的方差和为最小时,则认为线性回归方程与被拟合曲线达到了最佳配准(见图6-15)。即:
但实际空间中,数据往往是二维的,而且以更为复杂的方式变化,如图6-16所示,在这种情况下需用二次或高次多项式:
趋势面的次数并非越高越好,超过3次的复项多项式往往会导致解的奇异,因此,一般控制在二次变化曲面。
趋势面是一种平滑函数,很难正好通过原始数据点,这就是说在多重回归中的残差属正常分布的独立误差,而且趋势面拟合产生的偏差几乎都具有一定程度的空间非相关性。
整体趋势面拟合除应用整体空间的独立点内插外,另一个最有成效的应用之一是揭示区域中不同于总趋势的最大偏离部分。因此,在利用某种局部内插方法以前,可以利用整体趋势面拟合技术从数据中去掉一些宏观特征(例如最小二乘配置法)。
2、局部拟合技术
实际连续空间表面很难用一种数学多项式表达,因此,往往采用局部拟合技术利用局部范围内的已知采样点拟合内插值。这在表达地形变化特征的数字高程模型(DEM)内插中应用尤为广泛。
(1)双线性多项式内插
根据最近邻的四个数据点,确定一个双线性多项式:
当四个数据为正方形排列时,设边长为1,内插点相对于A点的坐标为X、Y,则有
当推广到双三次多项式时,采用分块方式,每一分块可以定义出一个不同的多项式曲面,当n次多项式与其相邻分块的边界上所有n-1次导数都连续时,称之为样条函数。
在数据点为方格网的情况下,采用三次曲面来描述格网内的内插值时,待定点内插值Zp为:
样条函数可用于精确的局部内插(即通过所有的已知采样点)。由于采用分块技术,每次只采用少量已知数据点,故内插运算速度很快,此外由于保留了局部微特征,在视觉上也有令人满意的效果。
(2) 移动拟合法
这种方法以待定点为中心进行内插。其原理是:定义一个合适的局部函数去拟合周围的数据点,通过解求拟合函数,解求出待定点的内插值。这种方法一般采取多余观测,利用最小二乘原理求解。通常做法是将坐标原点放置在待定点上,而采用的数据点应落在半径为R的圆内(如图6-17)。
图6-17移动拟合法原理示意图
局部函数可以为一次多项式(平面),但通常考虑到数据表面的光滑性,可采用二次多项式(曲面):
所采用的数据点坐标(X,Y)应满足:
R为数据圆半径。
数据圆的大小应以保证落入圆内的己知数据点数目不少于6个为基本准则,同时还应考虑数据表面的连续变化特征,在出现跃变的数据范围,应另选局部函数。
落入数据圆中的数据点,由于其对中心原点处的作用大小各不相同,因此,规定不同大小的权。根据通常数据内插的原则,距内插点越近其值越相近,越远则相差越大,可以依据内插点与已知点之间的空间距离分配权的大小。具体地说,当内插点无限接近于某个数据点时,则该数据点的权应无限大。
通常可采用的权的形式有:
究竟采取何种加权方式,应视具体的情况而定。
除上述介绍的两种主要局部内插法外,还有其它的一些内插方法,如最小二乘配置法,有限元内插法等。但是这些方法一般用于数据表面复杂,待求点众多的地形表面,用于生成规则的格网数字地面模型(DTM),有关这方面的介绍有专门的章节。
三、空间信息分类
空间信息分类方法是地理信息系统功能组成的重要组成部分。与地图相比较,地图上 所载负的数据是经过专门分类和处理过的,而地理信息系统存储的数据则具有原始数据的性质,这样用户就可以根据不同的使用目的对数据进行任意提取和分析。对于数据分析来说,随着采用的分类方法和内插方法的不同,得到的结果会有很大的差异,因此,在大多数情况下,首先是将大量未经分类的数据输入地理信息系统的数据库,然后根据用户建立的具体分类算法来获得所需要的信息。以下介绍空间信息分类中常用的几种数学方法:
1、主成分分析法
地理问题往往涉及大量相互关联的自然和社会要素,众多的要素常常给分析带来很大 困难,同时也增加了运算的复杂性。主成分分析法通过数理统计分析,将众多要素的信息 压缩表达为若干具有代表性的合成变量,这就克服了变量选择时的冗余和相关,然后选择 信息最丰富的少数因子进行各种聚类分析。
设有m个样本,n个变量,构造矩阵
Z=(Xij)nxm
其斜方差方阵R为实对称矩阵
R = Z·ZT/n = (rij)nxm
用Jacobi方法找出线性变换
使得y1, y2,…,yn互不相关,R矩阵的特征值越大,该主成分的贡献越大,因而可以选择累计贡献百分比在一定阈值以内的若干因子作为主因子参加分析运算。
2、层次分析法(AHP)
在分析涉及大量相互关联、相互制约的复杂因素时,各因素对问题的分析有着不同程度的重要性,决定它们对目标的重要性序列对问题的分析十分重要。AHP方法把相互关联的要素按隶属关系划分为若干层次,请有经验的专家们对各层次各因素的相对重要性给出定量指标,利用数学方法,综合众人意见给出各层次各要素的相对重要性权值,作为综合分析的基础。
3、系统聚类分析
系统聚类是根据多种地学要素对地理实体划分类别的方法。对不同的要素划分类别往往反映不同目标的等级序列,如土地分等定级、水土流失强度分级等。
系统聚类根据实体间的相似程度,逐步合并为若干类别,其相似程度由距离或相似系数定义,主要有绝对值距离、欧氏距离、切比雪夫距离、马氏距离等。
4、判别分析
判别分析与聚类分析同属分类问题,所不同的是,判别分析是根据理论与实践,预先确定出等级序列的因子标准,再将分析的地理实体安排到序列的合理位置上。对于诸如水土流失评价、土地适宜性评价等有一定理论根据的分类系统的定级问题比较适用。常规的判别分析主要有距离判别法和Bayes最小风险判别法等。
四、空间统计分析
1、常规统计分析
常规统计分析主要完成对数据集合的均值、总和、方差、频数、峰度系数等参数的统计分析。
2、空间自相关分析
空间自相关分析是认识空间分布特征、选择适宜的空间尺度来完成空间分析的最常用的方法。目前,普遍使用空间自相关系数—Moran I指数,其计算公式如下:
其中,N表示空间实体数目;xi表示空间实体的属性值; 是x的平均值;Wij=1表示空间实体i与j相邻,Wij=0表示空间实体i与j不相邻。
I的值介于-1与1之间,I=1表示空间自正相关,空间实体呈聚合分布;I=-1表示空间自负相关,空间实体呈离散分布;I=0则表示空间实体是随机分布的。Wij表示实体i与j的空间关系,它通过拓扑关系获得。
3、回归分析
回归分析用于分析两组或多组变量之间的相关关系,常见回归分析方程有线性回归、指数回归、对数回归、多元回归等。
4、趋势分析
通过数学模型模拟地理特征的空间分布与时间过程,把地理要素时空分布的实测数据点之间的不足部分内插或预测出来。
5、专家打分模型
专家打分模型将相关的影响因素按其相对重要性排队,给出各因素所占的权重值;对 每一要素内部进行进一步分析,按其内部的分类进行排队,按各类对结果的影响给分,从 而得到该要素内各类别对结果的影响量,最后系统进行复合,得出排序结果,以表示对结 果影响的优劣程度,作为决策的依据。其数学表达式为:
Gp=WiCip
式中,Gp表示声点的最终复合结果值,Wi表示第i个要素的权重,Cip表示第i个要素在p 点的类别的专家打分分值。
专家打分模型可分二步实现。第一步,打分:用户首先在每个feature的属性表里增加一个数据项,填入专家赋给的相应的分值;第二步,复合:调用加权复合程序,根据用户对各个feature所给定的权重值进行叠加,得到最后的结果。
第四节 数字地面模型及其应用
数字地面模型是地理信息系统地理数据库中最为重要的空间信息资料和赖以进行地形分析的核心数据系统。数字地面模型已经在测绘、资源与环境、灾害防治、国防等与地形分析有关的科研及国民经济各领域发挥着越来越巨大的作用。这里特别需要强调的是,数字地面模型的基本理论与数据处理方法,相当全面反映了地理信息系统空间信息分析的基本方法,本节试图通过较为细致、全面地对数字地面模型的概念、数据分析处理原理以及应用方法的叙述,使读者对数字地面模型有较为全面的认识,并由此强化对地理信息系统空间信息分析方法的理解与认识。
一、DTM与DEM的概念
数字地面模型是利用一个任意坐标场中大量选择的已知x、y、z的坐标点对连续地面的一个简单的统计表示,或者说,DTM就是地形表面简单的数字表示。
自从提出DTM的概念以后,相继又出现了许多其他相近的术语。如在德国使用的DHM (Digital Height Model)、英国使用的DGM(Digital Ground Model)、美国地质测量局USGS使用的DTEM (Digital Terrain Elevation Model)、DEM (Digital Elevation Model)等。这些术语在使用上可能有些限制,但实质上差别很小。比如height和elevation本身就是同义词。当然,DTM趋向于表达比DEM和DHM更广意义上的内容,如河流、山脊线、断裂线等也可以包括在内。
数字地面模型更通用的定义是描述地球表面形态多种信息空间分布的有序数值阵列,从数学的角度,可以用下述二维函数系列取值的有序集合来概括地表示数字地面模型的丰富内容和多样形式:
Kp=fk(up,vp) (k=l,2,3,…,m; p=1,2,3,…,n)
式中,Kp为第p号地面点(可以是单一的点,但一般是某点及其微小邻域所划定的一个地表面元)上的第k类地面特性信息的取值;up,vp为第p号地面点的二维坐标,可以是采用任一地图投影的平面坐标,或者是经纬度和矩阵的行列号等;m(m大于等于1)为地面特性信息类型的数目;n为地面点的个数。当上述函数的定义域为二维地理空间上的面域、线段或网络时,n趋于正无穷大;当定义域为离散点集时,n一般为有限正整数。例如,假定将土壤类型编作第i类地面特性信息,则数字地面模型的第i个组成部分为:
Ip=fi(up,vp) (p=l,2,3,…,n) (6-1)
地理空间实质是三维的,但人们往往在二维地理空间上描述并分析地面特性的空间分布,如专题图大多是平面地图。数字地面模型是对某一种或多种地面特性空间分布的数字描述,是叠加在二维地理空间上的一维或多维地面特性向量空间,是地理信息系统(GIS)空间数据库的某类实体或所有这些实体的总和。数字地面模型的本质共性是二维地理空间定位和数字描述。
在式6-1中,当m=l且fi为对地面高程的映射,(up,vp)为矩阵行列号时,式(6-1)表达的数字地面模型即所谓的数字高程模型(Digital Elevation Model,简称DEM)。显然,DEM是DTM的一个子集。实际上,DEM是DTM中最基本的部分,它是对地球表面地形地貌的一种离散的数字表达。
总之,数字高程模型DEM是表示区域D上的三维向量有限序列,用函数的形式描述为:
Vi=(xi,yi,zi)(i=l,2,3,…,n) (6-2)
式(6-2)中,xi、yi是平面坐标,zi是(xi,yi)对应的高程。当该序列中各平面向量的平面位置呈规则格网排列时,其平面坐标可省略,此时Vi就简化为一维向量序列zi, i=1,2,3,…n。
数字高程模型既然是地理空间定位的数字数据集合,因此凡牵涉到地理空间定位,在研究过程中又依靠计算机系统支持的课题,一般都要建立数字高程模型。从这个角度看,建立数字高程模型是对地面特性进行空间描述的一种数字方法途径,数字高程模型的应用可遍及整个地学领域。在测绘中可用于绘制等高线、坡度图、坡向图、立体透视图、立体景观图,并应用于制作正射影像图、立体匹配片、立体地形模型及地图的修测。在各种工程中可用于体积和面积的计算、各种剖面图的绘制及线路的设计。军事上可用于导航(包括导弹及飞机的导航)、通讯、作战任务的计划等。在遥感中可作为分类的辅助数据。在环境与规划中可用于土地现状的分析、各种规划及洪水险情预报等。
一般而言,可将DEM的主要应用归纳为:
(1)作为国家地理信息的基础数据:我国现在强调4D产品的建设。即:数字线化图(Digital Linear Graphs,简称DLG);数字高程模型(Digital Elevation Models,简称DEM);数字正射影像(Digital Orthophoto Quadrangles,简称DOQ);数字栅格图(Digital Raster Graphs,简称DRG)。并以前3D作为国家空间数据基础设施(NSDl)的框架数据。
(2)土木工程、景观建筑与矿山工程的规划与设计
(3)为军事目的(军事模拟等)而进行的地表三维显示
(4)景观设计与城市规划
(5)流水线分析、可视性分析
(6)交通路线的规划与大坝的选址
(7)不同地表的统计分析与比较
(8)生成坡度图、坡向图、剖面图,辅助地貌分析,估计侵蚀和径流等。
(9)作为背景叠加各种专题信息如土壤、土地利用及植被覆盖数据等,以进行显示与分析等等
二、DTM的数据采集与表示
1、DTM的数据源与采集方法
(1)以航空或航天遥感图像为数据源
这种方法是由航空或航天遥感立体像对,用摄影测量的方法建立空间地形立体模型,量取密集数字高程数据,建立DTM(见图6-18)。采集数据的摄影测量仪器包括各种解析的和数字的摄影测量与遥感仪器。摄影测量采样法还可以进一步分成:
a、选择采样 在采样之前或采样过程中选择所需采集高程数据的样点(地形特征点:如断崖、沟谷、脊等)。
b、适应性采样 采样过程中发现某些地面没有包含必要信息时,取消某些样点,以减少冗余数据(如平坦地面)。
c、先进采样法 采样和分析同时进行,数据分析支配采样过程。先进采样在产生高程矩阵时能按地表起伏变化的复杂性进行客观、自动地采样。实际上它是连续的不同密度的采样过程,首先按粗略格网采样,然后在变化较复杂的地区进行精细格网(采样密度增加一倍)采样。由计算机对前两次采样获得的数据点进行分析后,再决定是否需要继续作高一级密度的采样。
计算机的分析过程是,在前一次采样数据中选择相邻的9个点作窗口,计算沿行或列方向邻接点之间的一阶和二阶差分。由于差分中包含了地面曲率信息,因此可按曲率信息选取阈值。如果曲率超过阈值时,则必需进行另一级格网密度的采样。
(2)以地形图为数据源
主要以比例尺不大于1:1万的国家近期地形图为数据源,从中量取中等密度地面点集的高程数据,建立DTM。其方法有下列几种:
a、手工方法采用方格膜片、网点板或带刻划的平移角尺叠置在地形图上,并使地形图的格网与网点板或膜片的格网线逐格匹配定位,自上而下,逐行从左到右量取高程。当格网交点落在相邻等高线之间时,用目视线性内插方法估计高程值。它的优点是几乎不需要购置仪器设备,而且操作简便。
b、手扶跟踪数字化仪采集采集方式有:沿主要等高线采集平面曲率极值点,并选采高程注记点和线性加密点作补充;逐条等高线的线方式连续采集样点,并采集所有高程注记点作补充,这种方式适用于等高线较稀疏的平坦地区;沿曲线和坡折线采集曲率极值点,并补采峰—鞍线和水边线的支撑点,分别以等高线,峰—鞍链和边界链格式存储。
c、扫描数字化仪采集这种方式采集速度最快,但目前仅能以扫描分版等高线图方式采集高程。随着研究的不断深入,一些难点和瓶颈问题被解决,从地图扫描数据中自动地建立DTM技术必将达到实用水平。
(3)以地面实测记录为数据源
用电子速测仪(全站仪)和电子手簿或测距经纬仪配合PC1500等袖珍计算机,在已知点位的测站上,观测到目标点的方向、距离和高差三个要素。计算出目标点的x、y、z三维坐标,存储于电子手簿或袖珍计算机中,成为建立DEM的原始数据。这种方法一般用于建立小范围大比例尺(比例尺大于1:5000)区域的DEM,对高程的精度要求较高。另外气压测高法获取地面稀疏点集的高程数据,也可用来建立对高程精度要求不高的DTM。
(4)其它数据源
采用近景摄影测量在地面摄取立体像对,构造解析模型,可获得小区域的DTM。此时,数据的采集方法与航空摄影测量基本相同。这种方法在山区峡谷、线路工程和露天矿山中有较大的应用价值。
另外,航空测高仪可获得精度要求不太高的高程数据,也可以依此来构造DTM。
2、DTM的表示方法
(1)数学分块曲面表示法
这种方法把地面分成若干个块,每块用一种数学函数,如傅立叶级数高次多项式、随机布朗运动函数等,以连续的三维函数高平滑度地表示复杂曲面,并使函数曲面通过离散采样点。这种近似数学函数表示的DTM不太适合于制图,但广泛用于复杂表面模拟的机助设计系统。
(2)规则格网表示法
规则格网表示方法是把DTM表示成高程矩阵
DTM={Hij}, i=1,2,… m-l,m; j=1,2,… n-1,n
图6-20 显示了几种规则格网的DTM。
此时,DTM来源于直接规则矩形格网采样点或由规则或不规则离散数据点内插产生。由于计算机对矩阵的处理比较方便,特别是以栅格为基础的GIS系统中高程矩阵已成为DTM最通用的形式。高程矩阵特别有利于各种应用,但规则的格网系统也有下列缺点:
a.地形简单的地区存在大量冗余数据;
b.如不改变格网大小,则无法适用于起伏程度不同的地区;
c.对于某些特殊计算如视线计算时,格网的轴线方向被夸大;
d.由于栅格过于粗略,不能精确表示地形的关键特征,如山峰、洼坑、山脊、山谷等。为了压缩栅格DTM的冗余数据,可采用游程编码或四叉树编码方法。
(3) 不规则三角网(TIN)表示法
不规则三角网(Triangulated Irregular Network,TIN)是专为产生DTM数据而设计的一种采样表示系统。它克服了高程矩阵中冗余数据的问题,而且能更加有效地用于各类以DTM为基础的计算。因为TIN可根据地形的复杂程度来确定采样点的密度和位置,能充分表示地形特征点和线,从而减少了地形较平坦地区的数据冗余。TIN表示法利用所有采样点取得的离散数据,按照优化组合的原则,把这些离散点(各三角形的顶点)连接成相互连续的三角面(在连接时,尽可能地确保每个三角形都是锐角三角形或是三边的长度近似相等),如图6—21所示。
在概念上,TIN模型类似于多边形网络中的矢量拓扑结构,只是TIN中不必要规定“岛屿”或“洞”的拓扑关系。TIN把结点看作数据库中的基本实体。拓扑关系的描述,是在数据库中建立指针系统来表示每个结点到邻近结点的关系、结点和三角形的邻里关系,列表是从每个结点的北方向开始按顺时针方向分类排列的。TIN模型区域之外的部分由拓扑反向的虚结点表示,虚结点说明该结点为TIN的边界结点,使边界结点的处理更为简单。
在图6—22中,给出了TIN网络模型数据结构中的一部分,包括4个结点和2个三角形,数据则由结点列表、指针列表和三角形列表组成。区域中包括了边界结点,故设置虚指针。
由于结点列表和指针列表包含了各种必要的信息和连接关系,因而可以完成DTM数
值的有效提取和计算。
三、DTM的空间内插方法
DTM空间内插的概念十分简单,即在一个由x、y坐标平面构成的二维空间中,由已知若干离散点Pi的高程,估算待内插点的高程值。由于DTM采样的数据点呈离散分布形式,或是数据点虽按格网排列,但格网的密度不能满足使用的要求,这就需要以数据点为基础进行插值运算。
DTM内插按插点分布范围,可分为分块内插、剖分内插和单点移面内插三类,见图6—23。
分块内插,是把需要建立DTM的地区,切割成一定大小的规则方块,形状通常为正方形。它的尺寸应根据地形复杂程度和数据源的比例尺确定。在每一个分块上展铺一张数学面,相邻分块之间有适当宽度的重叠带,以使重叠带内全部数据点成为相邻块展铺数学面时的共用数据,保证一张数学面能够较平滑地与相邻分块的数学面拼接。这种内插方法的优点是可以得到光滑连续的空间曲面。
剖分内插是把需要建立DTM的地区切割成大小和形状不同的子区(剖分),子区间拥有公共边但不重叠,在该区内展铺一个数学面,内插剖分区内任意点的高程。该法只在剖分间边界端点处重合,通常没有严格重合的边界,所以既不连续,也不光滑。剖分多边形的顶点都是数据点,最常见的数据点个数为3,与TIN结构相同。
单点移面内插是以待插点为中心,以适当半径或边长的圆或正方形作为移动面去捕捉适当数目的数据点,并以此展铺一张数学面,内插该中心的高程。图6—24(a)分块内插区域;图6—24(b)剖分内插区域;图6—24(c)单点移面插值区域.
下面介绍几种常用的内插方法。
1、二元多项式拟合内插
用二元多项式来拟合地形表面,则有
在满项时,多项式的项数N=(m+1)2, 缺项时N=(m+1)2-N’ (N’>=1)。如果在内插区域内,数据点个数n大于多项式的项数k,这时可利用n个数据点的三维坐标(xk,yk,zk)以及逼近面和实际地面在数据点处的高程较差VK(k=l,2 … n)。列出n个误差方程,按最小二乘法解算出多项式的N个系数,即
按[VV]=min的最小二乘法原理,解得
C=[AT]-1[ATZ]
将所求参数代入二元高次多项式,可得二元高次多项式曲面拟合方程。若要求内插区任一点P的高程,只需把(xp,yp)代入该方程,就可求出Zp。
这种内插的特点是拟合曲面不通过所有数据点,而是取得最靠近数据点的光滑曲面,以保证邻块问的光滑连续拼接。
2、二元样条函数内插
二元样条函数是在分块范围内,按一定规则,用相邻数据点连线将块分割成若干个多边形分片(当数据点组呈正方形格网结点分布时,各分片是大小相等的正方形),通过每—分片上的全部数据点,展铺一张光滑的数学曲面,并使相邻分片间保持连续光滑的拼接。
对于DTM的内插,一般采用二元三次样条函数
写成矩阵形式为:
设分块范围内的数据点按单位边长正方形格网结点排列,一个单位边长的正方形为一个分片(图6—25)。取分片的左下角点为该分片平面直角坐标系统的原点,分片内任一点P的平面直角坐标为0<=xp<= 1,0<=yp<=1为了保证展铺的曲面在相邻分片上连续且光滑,必须满足弹性材料的力学条件,即
(1)相邻分片拼接处在x和y轴方向的斜率都应保持连续;
(2)相邻分片拼接处的扭矩连续。
拼接后整个分块的逼近面,就是二元三次样条函数曲面。
由于每个分片仅有4个格网结点信息(x,y,z),只能列出4个方程式,而函数的待定系数为16个,因此其余12个方程只能根据上述力学条件建立。据此建立的12个线性方程中,要用到沿x轴方向的斜率R,沿y轴方向的斜率S以及扭矩T,它们可由下式求得:
对于图6—25中画有阴影线的分片,其4个数据点的三维坐标分别是:0(0,0,z0),1(1,O,z1),2(1,1,z2),3(0,1,z3)。以O点为例,所建立的4个方程为:
Z0 = C00
按同样的方法可建立分别以1、2、3为顶点的共12个上述类型的方程。
用分片四个角点的高程,以及由各相关数据点高程计算得到的两个方向的斜率和扭矩数值组成一个4×4的常数矩阵A
按照对斜率R、S和扭矩T的定义,以及二元三次样条函数定义,得
把分片的4个角点的平面直角坐标代入该点,得
写成紧凑矩阵形式
A=XCYT
解此方程,有
C=X-1A(Y-1)T
把解得的系数阵代入式(6—3),则建立了二元三次样条函数式。对于分片中任意点P,把它的平面直角坐标(xp,yp)代入,就可求出其高程Zp。
3.移动拟合法内插,
移动拟合法是典型的单点移面内插方法。对每一个待定点取一个多项式曲面来拟合该 点附近的地形表面。此时,取待定点作为平面坐标的原点,并用待定点为圆心,以r为半径的圆内诸数据点来定义函数的待定系数(图6—24(c))。
设取二次多项式来拟合数据,则待求点的高程可写成下列的一般式:
zp=Az2 + Bxy + Cy2 + Dx + Ey + F
把坐标原点平移到待定点处,则 ,并代入上式,得移动拟合法二次多项式公式为:
为了求取式(6—4)中的待定系数,应以xp为圆心,以r为半径作圆(r的大小根据落入圆内的数据点个数n来确定,6<n<20),圆内的数据点均被采用,则可建立误差方程式
式中 为坐标原点平移至P点后,各数据点的三维坐标。
在求解系数时,可根据数据点至待定点的距离给予适当的权。权的给定按下式进行:
式中 。
根据最小二乘法原理,按vwvT=min的方法,建立法方程式,求解各待定系数。把系数代入插值公式,就能求得待定点高程zp。
4.加权平均内插法
加权平均法是最简单的单点移面内插方法,它是搜索区域内的高程数据点,并求得加 权平均值作为待定点的高程值。一般情况下,所考虑的权仅是距离的单调下降函数。为了提高插值精度,还应考虑数据点的分布方向,权衡搜索圆内所有点的方向和距离的分布情况赋权。搜索圆内数据点一般为4个,最多10个,这由调节搜索圆半径r的方法确定:
式中 xp,yp —— 待定点的平面直角坐标
zp —— 待定点高程
(xi,yi,zi)—— 搜索圆内数据点i的三维坐标,i=1,2 ... n, 4<n<10;
wTi —— 数据点的权。
四、DTM在地图制图与地学分析中的应用
DTM在科学研究与生产建设中的应用是多方面的,这里不可能将其所有的应用方面进行全面、系统的探讨,而仅以DTM在地学分析与地图制图中有典型意义的几个应用为例证说明其应用的基本思路和方法。它也将向我们展示栅格数据系统分析和应用的基本要点,对于帮助我们增强对栅格数据在地学信息自动处理中的作用和意义的理解有着十分重要的意义。
1、利用DEM绘制等高线图
如图6-26所示,利用DEM绘制等高线图,是以格网点高程数据或者将离散的高程数据由栅格追踪法原理转换为矢量等值线所产生的。该方法可以适用于所有的利用格网数据方法绘制等值线图。
2、利用DEM绘制地面晕渲图
晕渲图是以通过模拟实际地面本影与落影的方法有效反映地形起伏的重要的地图制图学方法。在各种小比例尺地形图、地理图,以及各类有关专题地图上得到非常广泛的应用。但是,传统的人工描绘晕渲图的方法不但费工、费时,而且带有很大的主观因素。而利用DEM数据作为信息源,在地面光照通量数学函数为自变量,计算该栅格应选用输出的灰度值。由此产生的晕渲图具有相当逼真的立体效果(如图6-27所示)。
3、透视立体图的绘制
立体图是表现物体三维模型最直观形象的图形,它可以生动逼真地描述制图对象在平 面和空间上分布的形态特征和构造关系。通过分析立体图,我们可以了解地理模型表面的平缓起伏,而且可以看出其各个断面的状况,这对研究区域的轮廓形态、变化规律以及内部结构是非常有益的。然而长期以来,人们为了在地图上形象地表示立体效果,制作了鸟瞰图、透视剖面图、写景图等。这些图解在较高艺术技巧的条件下,是可以得到好的效果。但表现它们要花费许多时间和人力,要有较高的艺术修养,因而难以普遍推广应用。机助制图为解决这方面的问题。提供了新的途径,而且在几何精度和实际艺术效果上,都能得到较好的保证。
计算机自动绘制透视立体图的理论基础是透视原理(见图6-29),而DEM是其绘制的数据基础。图6-28为制作透视立体图的基本流程。
4、DTM的地形分析
尽管DTM的应用十分广泛,但地形分析是其基本应用,其它应用都可由此推演、扩展。地形分析的内容有地形因子提取、地表类型分类以及剖面图的绘制等。现以栅格结构的DTM为例讨论地形分析。
(1) 坡度和坡向分析
坡度定义为水平面和地形表面之间夹角的正切值;坡向为坡面法线在水平面上的投影与正北方向的夹角(图6—32)。
坡度和坡向的计算通常在3×3个DTM格网窗口中进行。窗口在DTM数据矩阵中连续移动后完成整幅图的计算工作。
坡度slope=tgP=[ ]
坡向Dir=
式中的 一般采用2阶差分方法计算。对图6—33所示的格网,有
对于(i,j)点
式中 , 为格网结点x、y方向的间隔。
(2)地表粗糙度的计算
地表粗糙度是反映地表的起伏变化和侵蚀程度的指标,一般定义为地表单元的曲面面积与其在水平面上的投影面积之比。但根据这种定义,对光滑而倾角不同的斜面所求出的粗糙度,显然不妥当。实际应用中,以格网顶点空间对角线L1和L2的中点距离D来表示地表粗糙度(图6—34),D值愈大,说明4个顶点的起伏变化也愈大。其计算公式为
(3)地表曲率的计算
地面剖面曲率计算
地面的剖面曲率(profile curvature)其实质是指地面坡度的变化率,可以通过计算地面坡度的坡度而求得。
b、地面平面曲率计算
地面的平面曲率(plan curvature)是指地面坡度的变化率,可以通过计算地面坡向的坡度而求得。
5、谷脊特征分析
谷和脊是地表形态结构中的重要部分。谷是地势相对最低点的集合,脊是地势相对最 高点的集合。在栅格DEM中,可按照下列判别式直接判定谷点和脊点:
(1) 当(zi,(j-1) –zi,j)(zi,(j+1) –zi,j)>0 时
若 zi,(j+1) > zi,j 则VR(i, j) = -1
若 zi,(j+1) < zi,j 则VR(i, j) = 1
(2) 当(z(i-1),j – zi,j)(z(i+1),j –zi,j)>0 时
若 z(i+1),j > zi,j 则VR(i, j) = -1
若 zi,(j+1) < zi,j 则VR(i, j) = 1
在其它情况下,VR(i, j) = 0
其中 VR(i, j) = -1表示谷点
VR(i, j) = 1表示脊点
VR(i, j) = 0表示其它点
这种判定只能提供概略的结果。当需对谷脊特征作较精确分析时,应由曲面拟合方程建立地表单元的曲面方程,然后,通过确定曲面上各种插点的极小值和极大值,以及当插值点在两个相互垂直的方向上分别为极大值或极小值时,则可确定出谷点或脊点。
6、DEM水文分析
从DEM生成的集水流域和水流网络数据,是大多数地表水文分析模型的主要输入数据。表面水文分析模型用于研究与地表水流有关的各种自然现象如洪水水位及泛滥情况,或者划定受污染源影响的地区,以及预测当改变某一地区的地貌时对整个地区将造成的后果等。在城市和区域规划、农业及森林等许多领域,对地球表面形状的理解具有十分重要的意义。这些领域需要知道水流怎样流经某一地区,以及这个地区地貌的改变会以什么样的方式影响水流的流动。
地表的物理特性决定了流经其上的水流的特性,同时水流的流动将反过来影响地表的特性。对地表影响最大的水流特性为水流的方向和速度。水流方向由地表上每一点的方位决定。水流能量由地表坡度决定,坡度越大,水流能量也越大。当水流能量增加时,其携带更多和更大泥沙颗粒的能力也相应增加,因此更陡的坡度意味着对地表更大的侵蚀能力。另外由不同地表曲率决定的凸形或凹形地表也会对水流的流动产生影响,在凸形地表区域,水流加速,能量增大,其携带泥沙的能力增加,因而凸形剖面的区域为水流侵蚀地区。与此相反,在凹形剖面处水流流速降低,能量减少,导致泥沙的沉积。因此对水文分析来说,关键在于确定地表的物理特征,然后在此特征之上再现水流的流动过程,最终完成水文分析的过程。
从数字高程模型中可提取大量的陆地表面形态信息,这些形态信息包括坡度、方位以及阴影等。在大多数栅格处理系统中,使用传统的邻域操作便可以提取这些信息。集水流域和陆地水流路径与坡度、方位之类的信息密切相关,但同时也需要一些非邻域的操作计算,比如确定大的平坦地区范围内的水流方向等,因此简单的邻域操作对这些计算是不够的。为克服这些限制,达到提取地形形态的目的,一些研究者提出了既使用邻域技术又使用可称之为区域生长过程的空间迭代技术的算法,这些算法提供了从DEM中提取集水流域、地表水流路径以及排水网络等形态特征的能力。
上述算法的发展大体上经历了两个阶段,前一个阶段的算法一般基于格网点与空间相邻的8个格网之间的邻域操作,但不能很好地处理洼地;后一阶段的算法与此类似,但能完整地处理洼地与平坦地区。
以前的研究普遍认为被高程较高的区域围绕的洼地是进行水文分析的一大障碍,因为在决定水流方向以前,必须先将洼地填充。有些洼地是在DEM生成过程中带来的数据错误,但另外一些却表示了真实的地形如采石场或岩洞等。一些研究者曾试图通过平滑处理来消除洼地,但平滑方法只能处理较浅的洼地,更深的洼地仍然得以保留。处理洼地的另一种方法是通过将洼地中的每一格网赋以洼地边缘的最小高程值,从而达到消除洼地的目的。
下面介绍的算法以第二种方法为基础。通过将洼地填充,这些算法使洼地成为水流能通过的平坦地区。整个水文因子的计算由三个主要步骤组成,即无洼地DEM的生成、水流方向矩阵的计算和水流累积矩阵的计算,下面将对此分别进行介绍。需要指出的一点是,在整个DEM水文分析基础数据的计算过程中,虽然无洼地的DEM数据应首先生成,但在确定DEM洼地的过程中,使用了每一格网的方向数据,因此DEM水流方向矩阵的计算应最先进行,作为洼地填平算法的输人数据,在无洼地DEM的计算完成之后,重新计算经填平处理的格网的水流方向,生成最终的水流方向矩阵。
三个数字矩阵的获取的具体方法如下:
(1)无洼地DEM的生成
地形洼地是区域地形的集水区域,洼地底点(谷底点)的高程通常小于其相邻近点(至少八邻域点)的高程。对原始DEM先进行水流方向矩阵的计算,将结果矩阵中方向值满足下列条件的格网点作为洼地底点:①格网点的方向值为负值(方向值的具体意义在下面介绍);②八邻域格网点对的水流方向互相指向对方。对于自然地形进行分析不难知道,地形洼地一般有三种,分别是单点洼地、独立洼地区域和复合洼地区域。对于这三种洼地区域分别采用以下三种方法进行填平:
a、单格网洼地的填平的方法
数字地面高程模型中的单格网洼地是指数字地面高程模型中的某一点的八邻域点的高程都大于该点的高程,并且该点的八邻域点至少有一个点是该洼地的边缘点(即洼地区域集水流水的出口),对于这样的单格网洼地可直接赋以其邻域格网中的最小高程值或邻域格网高程的平均值。
b、独立洼地区域的填平方法
独立洼地区域是指洼地区域内只有一个谷底点,并且该点的八邻域点中没有一个是该洼地区域的边缘点。对独立洼地区域的填平可采用以下方法:首先以谷底点为起点,按流水的反方向采用区域增长算法,找出独立洼地区域的边界线,即水流流向该谷底点的区域边界线。在该独立洼地区域边缘上找出其高程最小的点,即该独立洼地区域的集水流出点,将独立洼地区域内的高程值低于该点高程值的所有点的高程用该点的高程代替,这样就实现了独立洼地区域的填平。
c、复合洼地的区域的填平方法
复合洼地区域是指洼地区域中有多个谷底点,并且各个谷底点所构成的洼地区域相互
邻接。复合洼地区域是地形洼地区域的一种主要表现形式。对于复合洼地的填平可采用下述方法:
首先以复合洼地区域的各个谷底点为起点,按水流的反方向应用区域增长算法,找出各个谷底点所在的洼地的边缘和它们之间的相互关联关系以及各个谷底点所在洼地的集水出水口所在的点位。出水口点的位置有两种,即在与“0”区域(非洼地区域)关联的边上或在与非“0”区域(洼地区域)相关联的边上。对于出水口位于与“0”区域相关联的边上的洼地区域,找出其出水口的高程最小的洼地区域,并将该区域内高程值低于该点的那些点的高程用该出水口的高程值代替。与该洼地区域相邻的洼地区域的集水出水口位于其所在洼地区域与该区域相邻的边缘,且其高程值低于该洼地区域集水出水口时,将这个洼地区域斤十时集水出水口点的高程值用该洼地区域集水出水口点的高程值代替。这样就将“0”区域复合洼地区域中的一个谷底点所构成的洼地区域填平,将所剩复合洼地区域用同样的办法依次对各个谷底点所构成的洼地区域进行填平,最后可将整个复合洼地区域填平。
用上述方法对数字高程模型区域中存在的洼地及洼地区域进行填平,可以得到一个与原数字高程模型相对应的无洼地区域的数字高程模型。在这个数字高程模型中由于无洼地区域存在,自然流水可以畅通无阻地流至区域地形的边缘。因此,我们可借助这个无洼地的数字高程模型对原数字模型区域进行自然流水模拟分析。
(2)水流方向矩阵的计算
水文因子计算的第二步是生成水流方向数据。对每一格网,水流方向指水流离开此格网时的指向。通过将格网x的8个邻域格网编码,水流方向便可以其中一值来确定,格网方向编码为:
64 128 1
32 x 2
16 8 4
例如,如果格网2的水流流向左边,则其水流方向被赋值32。方向值以2的幂值指定是因为存在格网水流方向不能确定的情况,需将数个方向值相加,这样在后续处理中从相加结果便可以确定相加时中心格网的邻域格网状况。另外需要说明的是出现在下面步骤中的距离权落差概念,距离权落差通过中心格网与邻域格网的高程差值除以两格网间的距离决定,而格网间的距离与方向有关,如果邻域格网对中心格网的方向值为l、4、16、64,则格网间的距离为√2,否则距离为1。确定水流方向的具体步骤是:
a、对所有DEM边缘的格网,赋以指向边缘的方向值。这里假定计算区域是另一更大数据区域的一部分;
b、对所有在第一步中未赋方向值的格网,计算其对8个邻域格网的距离权落差值;
c、确定具有最大落差值的格网,执行以下步骤:
①如果最大落差值小于0,则赋以负值以表明此格网方向未定(这种情况在经洼地填充处理的DEM中不会出现)。
②如果最大落差值大于或等于0,且最大值只有一个,则将对应此最大值的方向值作为中心格网处的方向值。
③如果最大落差值大于0,且有一个以上的最大值,则在逻辑上以查表方式确定水流方向。也就是说,如果中心格网在一条边上的三个邻域点有相同的落差,则中间的格网方向被作为中心格网的水流方向,又如果中心格网的相对边上有两个邻域格网落差相同,则任选一格网方向作为水流方向。
④如果最大落差等于0,且有一个以上的0值,则以这些0值所对应的方向值相加。在极端情况下,如果8个邻域高程值都与中心格网高程值相同,则中心格网方向值赋以255。
d、对没有赋以负值,0,l,2,4,…,128的每一格网,检查对中心格网有最大落差值的邻域格网。如果邻域格网的水流方向值为l,2,4,…,128,且此方向没有指向中心格网,则以此格网的方向值作为中心格网的方向值;
e、重复第4步,直至没有任何格网能被赋以方向值;对方向值不为1,2,4,…,128的格网赋以负值(这种情况在经洼地填充处理的DEM中不会出现)。
(3)水流累积矩阵的计算
区域流水量累积数值矩阵表示区域地形每点的流水累积量,它可以用区域地形曲面的流水模拟方法获得。流水模拟可以用区域的数字地面高程模型区域的流水方向数值矩阵来进行。其基本思想是,它认为以规则格网表示的数字地面高程模型每点处有一个单位的水量,按照自然水流从高处流往低处的自然规律,根据区域地形的水流方向数字矩阵计算每点处所流过的水量数值,便可得到该区域水流累积数字矩阵。在此过程中实际上使用了权值全为1的权矩阵,如果考虑特殊情况如降水并不均匀的因素,则可以使用特定的权矩阵,以更精确地计算水流累积值。图6-41(a)、(b)、(c)分别给出了一个简单的原始DEM矩阵以及计算出来的水流方向矩阵和水流累积矩阵。
6、基于DTM的可视性分析
可视性分析的基本内容如下:
(1)可视性分析的基本因子
可视性分析也称道视分析,它实质属于对地形进行最优化处理的范畴,比如设置雷达站、 电视台的发射站、道路选择、航海导航等,在军事上如布设阵地(如炮兵阵地、电子对抗阵地)、设置观察哨所、铺架通信线路等。
可视性分析的基本因子有两个,一个是两点之间的通视性(Intervisibility),另一个是可视域(ViewShed),即对于给定的观察点所覆盖的区域。
a、判断两点之间的可视性的算法
比较常见的一种算法基本思路如下:
①确定过观察点和目标点所在的线段与XY平面垂直的平面S;
②求出地形模型中与S相交的所有边;
③判断相交的边是否位于观察点和目标点所在的线段之上,如果有一条边在其上,则观察点和目标点不可视。
另一种算法是所谓的“射线追踪法”。这种算法的基本思想是对于给定的观察点V和某个观察方向,从观察点V开始沿着观察方向计算地形模型中与射线相交的第一个面元,如果这个面元存在,则不再计算。显然这种方法既可用于判断两点相互间是否可视,又可以用于限定区域的水平可视计算。
需要指出的是,以上两种算法对于基于规则格网地形模型和基于TIN模型的可视分析都适用。对于基于等高线的可视分析,适宜使用前一种方法。
对于线状目标和面状目标,则需要确定通视部分和不通视部分的边界。
b、计算可视域的算法
计算可视域的算法对于规则格网DEM和基于TIN的地形模型则有所区别。基于规则格 网DEM的可视域算法在GlS分析中应用较广。在规则格网DEM中,可视域经常是以离散的形式表示,即将每个格网点表示为可视或不可视,这就是所谓的“可视矩阵”。
计算基于规则格网DEM的可视域,一种简单的方法就是沿着视线的方向,从视点开始到目标格网点,计算与视线相交的格网单元(边或面),判断相交的格网单元是否可视,从而确定视点与目标视点之间是否可视。显然这种方法存在大量的冗余计算。Van和Kreveld提出了一种基于“线扫描”的算法,对于n个视点,算法的时间复杂度为O(nlgn)。总的来说,由于规则格网DEM的格网点一般都比较多,相应的时间消耗比较大。针对规则格网DEM的特点,比较好的处理方法是采用并行处理。
基于TIN地形模型的可视域计算一般通过计算地形中单个的三角形面元可视的部分来 实现。Lee讨论了离散的可视域的计算方法,实际上基于TIN地形模型的可视域计算与三维场景中的隐藏面消去问题相似,可以将隐藏面消去算法加以改进,用于基于TIN地形模型的可视域计算。这种方法在最复杂的情形下,时间复杂度为O(n2)。各种改进的算法基本都是围绕提高可视分析的速度展开的。
(2)可视性分析的基本用途
可视性分析最基本的用途可以分为三种:
a、可视查询。可视查询主要是指对于给定的地形环境中的目标对象(或区域),确定从某个观察点观察,该目标对象是可视还是某一部分是可视。可视查询中,与某个目标点相关的可视只需要确定该点是否可视即可。对于非点的目标对象,如线状、面状对象,则需要确定对象的某一部分可视或不可视。由此,也可以将可视查询分为点状目标可视查询、线状目标可视查询和面状目标可视查询等。
b、地形可视结构计算(即可视域的计算)。地形可视结构计算主要是针对环境自身而言,计算对于给定的观察点,地形环境中通视的区域及不通视的区域。地形环境中基本的可视结构就是可视域,它是构成地形模型的点中相对于某个观察点所有通视的点的集合。利用这些可视点,即可以将地形表面可视的区域表示出来,从而为可视查询提供丰富的信息。
c、水平可视计算。水平可视计算是指对于地形环境给定的边界范围,确定围绕观察点所有射线方向上距离观察点最远的可视点。水平可视计算是地形可视结构计算的一种特殊形式,但它在一些特殊领域中有着广泛的应用,而且需要的存储空间很小。