矢量地图叠加分析算法研究朱效民Word格式文档下载.docx
《矢量地图叠加分析算法研究朱效民Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《矢量地图叠加分析算法研究朱效民Word格式文档下载.docx（22页珍藏版）》请在冰豆网上搜索。
上述实现已经应用于实际的GIS系统中，取得了良好的效果。
关键词:
地理信息系统;
空间分析;
叠加分析;
矢量地图叠加;
点的包含性测试;
线裁剪
Researchonvectormapoverlay
ZhuXiaomin1），2），ZhaoHongchao1），LiuYan1），2），LiaoHaojun1），2），FangJinyun1），HanChengde1）
1）（InstituteofComputingTechnology，ChineseAcademyofSciences，Beijing100190）
2）（GraduateUniversityofChineseAcademyofSciences，Beijing100049）
Abstract:
Inthispaper，weintroduceaseriesofalgorithmsforvectormapoverlay:
Weproposedanalgorithmforlargeamountsofpointsinclusiontest，whichfirstpreprocessesthepolygonandthenadoptstheray-crossingideaforeachpointinclusiontest．Weproposedanalgorithmforlargeamountsoflinesclipping，whichbuildsindexesonthewholelineandthesegment，andavoidsunnecessaryintersection-computing．Formapoverlayonpolygons，weusetheimprovedplanesweepalgorithmtogetalltheintersectionpoints，andthecompletecategoriesondistributionofsegmentspassingthesamepointcoverallthespecialcasesofintersectionpoints．WerecordtheIDintheprocessofconstructingtheresultrings，whichsimplifiesthetwoprocesses-findingouterringforholesandattributepropagation．Allthegeometricalgorithmsandthecorrespondingoverlayfunctionsareimplemented，andthecomparisonswiththestate-of-the-artalgorithmsandArcGIS’sfunctiondemonstrateitscorrectness，robustness，efficiencyandusability．TheabovefunctionshavebeenusedinGISplatformandworkwell．
Keywords:
GIS;
spatialanalysis;
overlayanalysis;
vectormapoverlay;
pointinclusiontest;
lineclipping
0引言
GIS自从20世纪60年代在加拿大和美国诞生
以来，取得了飞速的发展，随着空间分析技术的不断发展，GIS也将从一般的空间事务处理向分析型空间决策支持方向迈进［1］，GIS已从数据库型GIS进入分析型GIS的阶段［2］。
矢量地图叠加分析（也叫
叠置分析，即vectormapoverlay）作为空间分析的一种，在GIS中扮演了重要的角色。
它是最重要的分析操作之一，是GIS中的核心操作之一，是一项组织来自不同数据源的数据的基本技术。
而且矢量地图叠加分析算法也是GIS中的难点，叠加过程被认为是GIS发展过程中的最具挑战性的问题之一;
退化情形以及浮点计算的有限精度使得矢量数据叠加计算非常困难。
因此，从学科的角度看，地图叠加分析算法是地理信息科学的重点与难点。
一般而言，叠加分析常常是对输入图层和叠加图层（一般为面图层）进行分析。
根据输入图层的类型不同，叠加分析一般有3种:
点面叠加、线面叠加和面面叠加。
点面叠加即得到在面内的所有点，并得到点在哪个面内。
线面叠加即区分在面内的线和在面外的线;
对于跨越面边界的线，将线在边界处打断;
对于在面内的线，需要得到线在哪个面内。
面面叠加即将两个面图层覆盖的区域区分为3部分，即公共部分（交）、只属于第1个图层的部分（差）以及只属于第2个图层的部分（差）。
对于交，需要得到是由哪两个区域相交产生的;
对于差的部分，需要计算此区域是由哪一个源区域产生的，至于是被哪些区域裁切产生的，常常不需要知道，因为一个区域产生的差区域常常是此区域被多个区域裁切得到的。
基于叠加分析对于系统及科学的重要性，许多专家学者将矢量地图叠加分析作为学位论文的研究工作做了深入系统的研究［3-12］。
大部分的文献只是简单介绍了叠加分析，并没有提出太多新的计算方法。
点面叠加，大部分采用的是经典的射线法［4，8-9］;
文献［10］介绍了利用堆栈结构和扫描线方法相结合，扫描结束即可得到所有点与所有多边形的关系;
文献［12］提出了基于单调链的点与多边形关系判断的方法。
线面叠加方面，阐述的较少。
面面叠加，大部分采用的是以单个多边形裁剪为基础［3，8］（常采用出点入点的思路［7］），一一循环计算的思路;
部分采用了一些空间索引策略（如R+树［4］）。
程振林的Han-E方法［5］采用了一种整体计算的思路，但在两组多边形都为复杂多边形集合时，内外环的匹配、属性继承、根据特定操作提取对应的多边形结果等存在一定的难度。
以上是国内近些年对叠加分析方面进行系统研究的相关文献。
此之外，针对单个叠加分析算法（点面、线面及面面叠加），也有许多相关研究，这些
相关工作将在下面对应主题中进行介绍。
1点面叠加
点的包含性测试即判断一个或者一系列点是否位于多边形的内部。
随着GIS的发展，大量点与多边形的包含性测试成为空间数据处理的一个核心问题。
诸多专家学者提出了一系列的算法来解决此问题。
这些算法有的不对多边形进行预处理，如角度之和法和射线法［13］等，有的则预处理多边形，如将多边形分解成许多简单单元（如凸多边形［14］）。
根据现有文献对算法的分析，最有效的方法是基于Cell的方法［15］和基于多边形凸剖分的方法［14］。
基于Cell的方法将多边形栅格化，然后对不含多边形边的栅格计算其在多边形内部还是外部。
对点进行包含性测试时，首先计算点在哪个栅格内;
若栅格不含任何边，则可以在O
（1）时间内得到点与多边形的关系;
否则，需要计算点与在当前栅格内边的关系。
此方法对部分点的包含性测试时间复杂度是O
（1），因此对于一些简单的多边形，它是非常高效的。
基于多边形凸剖分的方法是目前已知最高效的方法。
该方法将多边形进行凸剖分，并用BSP（binaryspacepartition）树管理这些凸剖分。
对任何一个点，查询BSP树，计算可能包含其的凸剖分，然后通过点与凸剖分关系得到点与多边形的关系。
其包含性测试较迅速，而其预处理时间相对较多，尤其对于边数较多的复杂多边形。
1．1本文方法
提出了一种新的大量点的包含性测试方法。
方法包含两个步骤:
预处理和包含性测试。
预处理即将多边形的MBR（minimumboundingrectangle）分成若干相等的竖直条带，并且计算每条边的位置。
如图1所示，多边形有16条边，标记为0～15;
将MBR分成16个条带，每个条带包含上下两部分，然后计算每条边经过的条带。
经过预处理后，可以得到每个条带关联的边。
如图1所示，条带0的上部、下部分别包含的边为0、15和12、13。
经过预处理后，算法对每个点进行包含性测试。
首先计算点所在的条带，然后构造辅助射线，并与所有当前条带内的边进行相交测试，根据相交次数的奇偶性得到点与多边形的关系。
对于位于MBR上半部分的点，辅助射线竖直向上;
否则竖直向下。
这样的策略以一个判断浮点数大小的Boolean操作为
代价，避免了随意构造射线可能造成的更多相交
情形。
图1多边形的预处理
Fig．1Preprocessingonapolygon
以射线竖直向上的情形为例，介绍相交判断方法。
如图2所示，点所在的水平直线与构造的射线将整个平面划分为3部分:
Ⅰ、Ⅱ和Ⅲ。
对于每条边，首先计算其两个端点的位置。
如果两个端点在同一区域，那么边与射线不相交;
否则，边的两个端点位置可能的组合为:
Ⅰ和Ⅱ，Ⅰ和Ⅲ，Ⅱ和Ⅲ。
两个端点在区域Ⅰ和Ⅱ，边必然与射线相交。
两个端点在Ⅰ和Ⅲ，如果在Ⅲ的端点位于Ⅲ的左半部分，边不与射线相交;
否则，假设当前边方向是右下往左上，如箭头所示，那么顶端点T一定位于当前边的右侧。
因此，只需要计算底部端点（当前点）与当前边的位置关系（左侧或右侧）。
如果点位于边的右侧，则边与射线不相交，否则即相交。
第3种情形与第2种情形类似。
上述步骤中左右位置关系计算需要乘法和减法操作，其他计算只需要浮点数的大小比较操作，因此具有较高的效率。
射线竖直向下时，同理可得。
对于一些特殊情形，如:
相交在端点或部分重叠等情形，采用微小平移的方法，将构造的射线
图2判断边是否与射线相交
Fig．2Judgingwhetheranedgeintersectwiththeray
左移微小的距离。
只要射线平移的方向是一致的，这些特殊的相交情形都能得到正确的判断结果。
1．2算法分析
下面从最坏和平均两种情形分析算法的时间和空间复杂度。
假设多边形的边数为N，使用的条带数目是N/2。
预处理过程的基础操作即对每条边计算一次，然后将其附加到对应的条带中去。
其时间空间复杂度为O（N×
M），其中M表示每个条带包含的边数。
因此最坏情形下其复杂度为O（N×
N）。
真实世界数据中，每个条带包含的边是一定的，其实际空间复杂度是O（LM×
N），LM为有限数值，一般与N无关，是常数或者接近常数。
基于真实地理数据的分析如图3所示。
3组数据分别有2525、2501和
122552个多边形。
图中“aver”是指平均每个条带内边的数目，即包含性测试的时间复杂度;
“times”是所有边的存储次数，即预处理阶段的空间复杂度。
如图所示，times接近线性于N，即平均情形下，其时间空间复杂度约为O（N）。
包含性测试的空间复杂度为O
（1）。
时间复杂度的最差情形是O（N）。
然而，基于预处理阶段的空间复杂度分析，我们认为其时间复杂度接近一个有限值，这个值一般远小于N。
图4所示是对前面3组数据的时间复杂度分析，也就是aver随多边形边数的变化趋势。
X轴是边数N，对应的aver是平均每个条带内的边的数目，log是指logN，用来衡量aver的。
如4（a）（b）所示，aver值是阶段性的，其值与N有一定关系，但总体上是一个接近常量的值;
而且，aver总体上小于Log，因此我们可以认为其时间复杂度接近O
（1）。
图4（c）中，可以看到aver要大于Log，但是小于log2，近似为O（（logN）1．5）。
我们认为O（logN～（（logN）1．5））为其时间复杂度。
这与数据集1和2有所不同，究其原因，是因为数据集1和2都是较为简单的多边形，而数据集3的多边形较为复杂，且常带洞。
综上，对于点与有N条边的多边形进行包含性测试时，平均情形下，预处理的时间和空间复杂度都约为O（N），包含性测试的空间复杂度为O
（1），时间复杂度对于简单多边形接近O
（1），对于复杂多边形约为O（（logN）1．5）。
1．3算法对比测试
测试的输入是由多边形MBR过滤得到的点。
本文算法（Our）与Cell［15］算法（其源码为地址http:
//www．iamg．org/index．php/publisher/articleview/
frmArticleID/104）进行了对比。
对比测试在同一台IBMX31笔记本（IntelPM处理器1．7GHZ，512M内存）上进行，并且两个算法在VC++．NET中实现，使用相同的编译、链接、执行环境。
对比的结果如
表1所示。
3个多边形是上述源码带的数据，点数据是一个均匀分布于多边形的MBR中的1000×
1000的点集。
从表1对比可得，本文方法预处理及点包含性测试的时间都少于Cell方法。
此外，本文方法得到的结果是精确的，Cell的计算结果正确率是98%［15］。
基于凸剖分的方法较为复杂，不同的实现有不同的效率，因此我们没有实现其算法。
凸剖分的方法也用这3个多边形做了测试，其测试结果在文献［14］中。
文献［14］中的Grid方法即Cell方法，而且对比也是基于文献［15］带的源码及数据。
以Grid法为基准，计算可得对应3个数据的预处理以及包含性测试的相对值如表2所示，可得本文算法（Our）具有较少的初始化时间，对于复杂多边形更为明显。
点的包含性测试而言，本文算法也稍优于方法基于凸剖分（ConvexA，ConvexB）的方法。
表2以Grid法为基准的相对对比结果
Tab．2TherelativecomparisonresultbasedontheGridmethod
端点预处理时间相对值包含性测试相对值/%
个数GridConvexAConvexBOurGridCovexAConvexBOur
理和包含性测试方面都要优于Cell、ConvexA、ConvexB方法。
以此几何计算算法为基础，辅以数据读写模块，实现了基于ShapeFile文件的矢量点面叠加分析功能。
与ArcGIS的点面叠加分析（点面求交功能）功能对比如表3所示。
正确性方面，本文方法（VegaGIS）与ArcGIS得到的结果要素数是一致的;
效率方面，本文的叠加分析功能整体计算时间为ArcGIS同样功能的20%～25%。
以MBR内的点数
表3点面叠加分析功能对比测试表
Tab．3Comparisontableofpoint-polygon-overlay
10010．170．060129．535．923．0
MBR中
多边形
时间/s结果要素数裁剪
129412．440．8970154．492．639．7
点的个数的个数
VegaGISArcGISVegaGISArcGIS区域
28105111．897．990．22119．128．514．9
结合以上两表的绝对对比与相对对比数据，可以得到一个凸剖分方法的近似计算值。
以相对值为基础，可得其预处理、包含性测试的时间变化曲线分别如图5、图6所示。
从图示可得，本文方法在预处
2620912．29．01971019710省域
417726
1
21．9
96
272536
区域2
597531
2
32．9
144
407147
工作范围
649567
2525
52．9
223
639684
行政区划
264722112．051．0134611134611区域1
为基础，可得其随着点数变化计算时间的变化曲线
如图7所示。
图7点面叠加分析时间随点数变化曲线图Fig．7Transformationcurveofpoint-polygonoverlaytimecostonpointnumber
2线面叠加
GIS中的线裁剪的对象往往是大量的线，其核心问题是如何避免不必要的求交运算、快速求取所有交点。
如果采用传统的线裁剪算法，则需要一一计算，即使采用了高效的单条线裁剪算法，如采用错切变换法［16］求交点，或者将多边形进行凸剖分等［17］，仍然不够高效。
有一系列高效的算法解决了大量线段的求交问题，如扫描线算法［18］、基于空间索引的方法等。
文献［19］分类总结了这些高效的求交点算法。
真实地理数据中，有很多线不与多边形的任何边相交，因此这些线无需参与求交运算。
传统的扫描线算法、空间索引算法等，对所有输入线段不加区别的进行处理，这些线同样参与了运算—在扫描线算法中，它们被插入、删除;
在空间索引方法中，它们需要参与建立索引以及查询的时间—这些都导致了不必要的运算。
2．1本文方法
基于以上分析，采用双层的索引机制，基于线、线段建立索引。
本文方法依赖于一个高效的R-tree实现。
R-tree提供索引建立以及索引查询的功能。
第1层索引是基于整条线MBR的。
索引建立之后，基于多边形的每条边的BR（boundingrectangle）进行查询，得到了所有可能与多边形相交
的线。
对于这些可能与多边形相交的线，其线段的
数量仍然是相当多的，因此基于线段的第2层索引
将被建立。
然后针对多边形的每条边进行查询，得到可能与当前边相交的所有线段，并计算交点。
之后对这些交点按照所属线的不同以及在线上出现的顺序进行排序，计算交点的出入属性，并根据奇偶性设置所有交点的属性，最后根据不同的操作，截取在多边形内部或者外部的部分。
对于那些没有产生任何交点的线，为每条线选取一个点作为代表点，然后执行一次点与多边形的包含性测试，即可得到点代表的线与多边形的关系，然后将其完整地插入到内部或者外部结果集合中去。
2．2算法分析
就复杂度而言，线裁剪算法中复杂度最高的即为求交点部分。
对于M条线段，假设平均具有N个点，这些线段被一个具有P条边的多边形进行裁剪，其求交部分的复杂度为O（M×
N×
P）。
对最坏的情形（任何一条线的MBR都与多边形相交，线的任意一条线段都与多边形相交），本文方法具有相同的时间复杂度。
然而在真实世界的地理数据中，最普遍的情况即用一个多边形对大量线进行裁剪，相对于多边形的MBR，线的MBR一般是很小的;
只有一部分线与多边形的边相交;
对于多边形相交的线，只有部分线段与多边形的边相交。
从理论上分析，面积的比值在一定程度上代表了过滤率。
将MBR的面积（多边形所有边的MBR的面积与整个多边形MBR的面积的比值）作为一个衡量因素，来分析过滤率。
分别选取两组地理数据，一组是描述行政区域的多边形数据，边数较少，较为简单;
另一组是描述土地利用的数据，边数较多且多含有洞，较为复杂。
图8中两图分别描述了简单、复杂多边形的过滤率随着多边形边数的变化而变化的趋势。
可得大部分多边形的面积比值都在0．1左右，从而可以利用双重索引策略，达到一个较为理想的过滤率，减少不必要的运算。
2．3线裁剪功能测试
我们没有从算法角度与线面一一裁剪的方法进行对比，而是辅以必要的读写过程，实现了整个线面叠加分析方法，并与ArcGIS的对应功能进行了对比，如表4所示。
计算时间随要素个数的变化趋势图如图9
升级会员 