一种基于栅格方法的空间谓词判断方法及其系统

一种基于栅格方法的空间谓词判断方法及其系统董慧***，程振林*，方金云*，赵红超*（*中国科学院计算技术研究所北京市海淀区科学院南路6号100190）（**中国科学院研究生院北京市石景山区玉泉路19号（甲）100049）摘要：地理信息获得了越来越广泛与深入的应用。空间分析是地理信息系统平台最核心的计算之一。矢量方法是面向物体的描述，物体间的几何关系隐含，关系判断需要基于计算几何算法定位、分析和检索，时间与空间复杂度较高。为避免其缺点，本文提出了基于栅格的空间谓词算法实现框架。此算法基于亚像素精度，可以较准确记录边界栅格的覆盖面积，通过判断两个图层相应栅格的覆盖面积，即可得出空间关系。此方法的正确率大大高于四色栅格签名（4CRS）。同时，栅格索引中保存了要素属性信息等，可以为空间谓词判断提供更有用的结果信息。关键词：空间谓词，栅格化方法，4CRS（fourcolorrastersignature），亚像素精度一．介绍随着GIS自身的发展和经济与社会的信息化，GIS开始融入信息技术的主流。由于GIS技术能较好地解决基于时空框架的数据建模问题，填补了传统信息技术在这方面的空白，逐步成为信息技术的核心支撑技术[1]。基于Web的地图应用使得GIS的用户从专业人士迅速扩大到公众。以Web编程接口的形式提供空间信息服务成为GIS与其他的业务信息系统进行应用整合的重要途径，这为GIS应用开辟了更广阔的应用范围和场景。但是，基于Internet的GIS的体系结构决定了大量的业务逻辑集中在服务器端。满足众多用户（包括Web服务客户端）的访问并保证服务质量，给后端服务器的性能、可扩展性提出了更高的要求。空间谓词是比较两个空间对象并返回一个布尔变量值作为结果，它表明了存在于两个空间对象之间特殊的关系。如：是否相交、是否相互包含等。空间谓词是GIS的核心，是GIS区别于一般的信息系统、CAD（计算设计辅助设计）或者电子地图系统的主要标志之一。OGC(OpenGeospatialConsortium，开放地理信息联盟)的web要素服务（WebFeatureService）规范中的空间过滤器是通过空间谓词方式获取要素数据的有力方式，规范中提出了Disjoint/Intersect、Equals、Within/Contains、Overlaps、BBOX等多种空间谓词过滤器。空间谓词，配合空间数据的属性信息，能提供强大、丰富的空间数据查询功能。空间谓词以及空间分析功能具有算法复杂、计算密集等特点，如何在Web上提供空间谓词的功能，在学术研究与工程实践上，都具有重要的意义，能促进GIS应用的Web迁移并促进GIS应用与其他应用的融合。随着网络地图服务的流行，如何在网络地图服务器上提供空间谓词功能成为需要解决的问题。一种常见的空间谓词判断方法是利用计算几何来实现。作为计算机科学的一个分支，计算几何主要研究解决几何问题的算法。常见的做法是针对两个多边形进行，在大量的多边形计算面前无能为力。如果采用“暴力”算法，通过反复调用两个多边形空间谓词的算法来完成，则算法实现计算复杂度高，实用性差。如基于出入点判别的空间谓词实现方法，如何确定交点的进点、出点属性在实际的图形中会遇到众多的特殊情况。特别是在发生了线段与线段交在端点、线段与线段重叠的情况下，如何区分交点的出点、入点情况非常复杂导致效率降低。这类做法中采用的线段求交算法一般是采用平面扫描算法，优点是结果比较精确，缺点是由于要进行频繁的坐标排序、角度计算等操作，计算量大[2]。为解决上述问题，本文提供了一种基于栅格方法的空间谓词判断方法及其系统，能够减少计算量，提高计算效率。四色栅格签名（fourcolorrastersignature--4CRS)[3][4]可以实现空间关系的判断，但是4CRS结果不够精确，只能以四种栅格属性区分栅格面积--空、满、强（一半以上）、弱（一半以下），只能判断出强*强是确定相交的，而强*弱、弱*强、弱*弱则都是不确定的结果，还需要精确计算。在计算近似面积、置信度区间等都是用数学期望和概率等公式来估计，而真实数据可能并不符合这种规律。本文使用的基于亚像素精度的方法能准确记录边界格子的面积，因此可以根据两个对象占用的单元格的面积来确定是否相交。如同一个位置的栅格，第一个图层占有此栅格为49%，第二个图层为52%，相加100%,所以确定相交，而4CRS则不能确定。而且，4CRS只能实现两个多边形的判断，而不能处理多个多边形的图层。而本文采用的方法不但能精确记录栅格面积百分比，提高结果的准确性，并能把百分比和图层属性信息等都保存下来，给用户返回更多有效的结果信息（结果多边形的ID号等）。二．算法框架整个算法分为两步：栅格索引生成和基于索引的栅格空间谓词实现。整个栅格索引生成过程算法框架如图1所示。首先，对地理要素的矢量点进行坐标转换，便于利用亚像素精度进行后续的计算。然后，进行轮廓扫描和绘制控制器计算填充单元跨段，在此过程中计算出图形轮廓处的该图形占栅格单元的面积的百分比，该百分比称为栅格单元对应于该图形的实际占用面积。最后，将要素ID信息和实际占用面积信息存入以像素为索引的结构里。坐标转换通道轮廓扫描绘制控制器计算填充单元及跨段水平扫描线填充单元及跨段要素信息（栅格ID号、面积百分比等）信息按扫描线要求存入缓存输入矢量要素生成底图索引图1.步骤一的算法框架经过底图索引生成这个步骤之后，即完成了栅格化过程，以栅格矩阵形式存储，其中坐标值为索引，所属的多边形ID号和所占的栅格面积百分比等为图层信息，以文件形式保存下来。这样，两个图层的空间谓词判断的时候，只需要找到同一栅格，对其面积等进行判断，从而获得符合条件的多边形ID号返回结果。示意图如图2。图2步骤二的流程图三．基于栅格方法的空间谓词实现栅格索引生成栅格索引生成部分，用于对于输入的图层生成对应的栅格底图，栅格底图的栅格单元以压盖所述栅格单元的图层的图形的要素ID为要素索引，栅格单元以坐标值为位置索引，每个栅格单元具有对其压盖的图形在所述栅格单元的实际占用面积的比值。图3(a)地理要素的几何表示(b)栅格近似实现图片像素与矢量点关联，具体指在绘制地理要素的图形同时，生成一张与该图形匹配的栅格底图，该栅格底图就是与该图形对应的索引图，如图3所示。栅格底图的每个栅格单元都对应了该地理要素的属性信息内容，例如地理要素的要素ID号等，以要素ID号为栅格单元的要素索引。因此完全不存在搜索存在重叠区域这种情况。可作为索引快速提取要素信息，也是实现下述空间谓词的基础。生成栅格索引需要以下步骤：（1）输入矢量方式表示的图层中地理要素的图形的矢量点，按显示屏幕的分辨率对地理要素的矢量点坐标进行坐标转换，按显示屏幕的像素点进行栅格划分，栅格单元以坐标值为位置索引，栅格单元以对其压盖的图形的要素ID为要素索引。如果没有图形对栅格单元进行压盖，则栅格单元的要素索引和实际占用面积的信息都为初始化时的缺省值。坐标值是根据实际地理坐标投影到屏幕坐标的栅格行列号。在本文中，将double类型的矢量点坐标都乘以256，相当于将该坐标的二进制表示左移8位。这种坐标转换的优点是考虑了小数部分对像素的栅格单元权值(cover)的影响，便于利用亚像素精度进行后续的计算。（2）对图层中的图形进行轮廓扫描，对于每个图形，按如下公式计算图形的轮廓线经过的每一栅格单元的权值和覆盖面积[5]，cov21erfyfy（1）(21)covareafxfxer（2）其中，（fx1，fy1）为经过该栅格单元的轮廓的线段起始点的小数坐标部分，（fx2，fy2）为经过该栅格单元的轮廓的线段终止点的小数坐标部分，cover为权值，area为面积。运用了亚像素精度（subpixelaccuracy）的布兰森汉姆（Bresenham）[6]生成直线算法进行轮廓所描，Bresenham生成直线算法是一种基于误差判别式来生成直线的方法。与传统Bresenham所不同的是，该算法利用误差判别选择像素的过程是基于亚像素的，把一个像素分成N×N个小像素，例如，将单位栅格平均分成了256×256个子像素。（3）对于每个图形，绘制控制器遍历图形的轮廓经过的栅格单元，依据栅格单元的覆盖面积判断所述栅格单元是否被图形完全填充，对完全填充的栅格单元和未完全填充的栅格单元分别进行标记，将轮廓内的栅格单元进行跨度填充，将这些栅格单元标记为完全填充。首先，绘制控制器将轮廓经过的所有栅格单元，例如上述计算的A和B等整数坐标属于此处的单元格，进行排序，将X坐标相同的点按照Y坐标由小到大排序。然后，绘制控制器按照行扫描顺序，扫描顺序指按照X坐标由小到大扫描每行，从最小行到最大行逐行扫描，利用图形中每个扫描栅格单元的area进行判断栅格单元是否完全被图形填充，对于未完全填充的栅格单元，进行标记，例如利用add_cell（）函数进行标记，对于完全填充的栅格进行标记，例如利用add_span（）函数进行标记。未完全填充的栅格指图中所计算的area没有完全覆盖当前栅格，如图4所示。图4area和cover示意图由于cover有正负，所以在扫描每一行的时候都将所有轮廓的单元格的cover相加，因为图形的轮廓是闭合的，所以当该跨段标记完以后遍历的所有栅格单元的cover加和会为0，通过这样的判断，能够自动找到哪部分需要填充，哪部分不需要填充。（4）将完全填充的栅格单元的图形在栅格单元的实际占用面积的比值设置为1；对于未完全填充的栅格单元，根据所述栅格单元的权值和覆盖面积计算所述栅格单元的图形在栅格单元的实际占用面积的比值；保存栅格单元的位置索引、要素索引和实际占用面积的比值，进而生成所述图形对应的栅格底图。轮廓扫描的时候保存的栅格单元的area并非是有效面积，因为不知道图形的内部还是外部。内部还是外部同轮廓的走向相关，由图4可以看出，area是由轮廓跟栅格单元的左侧围成的面积，而当轮廓的走向是顺时针时，真正的面积应该是栅格单元的面积减去area得到的，剩余一部分的面积，所以真正的轮廓的栅格单元实际被图形包含的面积需要在绘制的时候才能确定。计算每个未完全填充的栅格单元的实际面积。确定出一个图形轮廓所经过的每个栅格单元,然后再计算每个栅格单元中有多少面积是落在图形内,把整数坐标相同的栅格单元的area和cover累加，然后通过cover和area共同作用计算出轮廓经过的栅格单元实际被图形占用面积，用alpha表示。轮廓走向由cover记录，cover正负值表示不同的轮廓走向，顺时针或逆时针，所以真正需要计算的属于图形内部的面积是需要这两个变量共同计算得出，通过area和cover共同计算出alpha。alpha=cover-256256area=cover（1-2562562fx1fx）（3）这个值越大，则代表图形压盖的格子百分比越多。经过以上步骤，能更好的利用已有的矢量数据栅格化近似，提高处理速度，也就是完成了栅格化过程，以栅格矩阵形式存储，其中坐标值为位置索引，属性信息等为图层信息，以文件形式保存下来。基于索引的栅格空间谓词实现对于两个待比较的地理要素的图层，两个图层的栅格底图中坐标相同的栅格单元相互对应，将相对应的两个栅格单元的实际占用面积的信息进行比较，得出所述两个图层的空间谓词判断结果。对于遍历的栅格单元，计算所述栅格单元对应于两个图层的实际占用面积的比值的和，判断所述和是否大于100%，如果是，停止遍历，空间谓词判断结果为两个图层相交。并且在得出空间谓词判断结果后还用于返回所述栅格单元在两个栅格底图中对应的要素ID。首先，分别确定两个待比较的图层的外包，以两个图层的外包的相交区域为遍历区域；对遍历区域中的栅格单元进行遍历。在遍历的所有栅格单元的所述实际占用面积的比值的和都不大于1时，按如下公式计算平均期望值，countcellareacellarea21（4）其中，1cellarea为遍历的栅格单元对应于一个图层的实际占用面积的比值，2cellarea为所述栅格单元对应于另一个图层的实际占用面积的比值