一种平衡面积、拓扑和形状特征的矢量多边形栅格化方法

本发明涉及一种平衡面积、拓扑和形状特征的矢量多边形栅格化方法，矢量数据向栅格数据转换方法技术领域。

技术背景

矢量和栅格是地理信息系统(gis)中空间数据的基本格式(maguireetal.,1991)。这两种数据结构存在着根本差异(mattikalli,1995)，导致许多数据处理和分析工具不能通用，因此常面临矢量和栅格相互转换的问题。栅格数据由于结构简单、数据采集和分析方法简便等优点，通常比矢量数据更适合科学研究(peuquet,1984；goodchild,2011)。近年来，遥感技术极大拓展了栅格数据的应用场景。尤其基于gis的空间分析和空间建模等应用(demers,2001；kuhnertetal.,2005)，如资源与环境监管(hongetal.,2016)、土地利用/土地覆被变化模拟(mustafaetal.,2017)，促使矢量-栅格转换的研究更加深入。多边形栅格化时内部像元取值一般来自该多边形，而关键问题在于边界像元的确定和赋值。已经开发了边界代数法(feito,1995)、边界追踪法(wu,1998)等算法寻找边界像元。由于每个栅格像元只能被赋予一个属性值，当边界像元与多个多边形重叠时，必须制定规则对多个属性值进行取舍(wehdeetal.,1980)。此类像元常规的取舍规则包括中心点法和最大面积法，前者依据包含像元中心点的多边形，而后者根据与像元重叠面积最大的多边形来确定该像元取值。

栅格数据的离散像元无法表达多边形的平滑边界，导致矢量-栅格转换是一个有损过程(vanderknaap,1992)。栅格化误差主要表现为多边形的周长、面积、形状、位置和拓扑等发生变化(carver&brunsdon,1994；galton,2003；shortridge,2004)。尽管减小像元尺寸可以有效降低各种栅格化误差，但会导致数据量和处理时间显著增加(bettingeretal.,1996；baietal.,2011)，因此在实际应用中栅格像元并非越小越好。分析与评价栅格化误差是gis领域的一项重要工作，其中以面积误差为主(wadeetal.,2003)。bregt等(1991)提出基于双重转换的面积误差分析方法，发现该误差与边界指数(单位面积内多边形的边界长度)呈线性正相关关系。kam等(2000)发现面积误差因各种gis软件对栅格像元的取舍规则不同而存在差异，这会对依赖于像元面积和数量的空间分析带来问题。zhou等(2006)根据多边形与其边界像元的空间关系，采用领域补偿原则最小化面积误差。liao等(2012)通过分析栅格化误差的影响因素，建立了面积误差与弧长密度、多边形密度及栅格像元大小之间的函数关系。wei等(2019)基于网格纯度指数开展多尺度栅格化试验，面积误差分析结果表明该方法明显优于最大面积法。

矢量-栅格转换也会显著影响多边形的形状和拓扑特征。congalton(1988)研究了耕地、草地和林地3种景观的空间复杂性和误差模式，发现周长-面积比较大的多边形在栅格表达中容易丢失或严重变形。bettinger等(1996)指出，当栅格像元大小相对多边形过于粗糙时，后者的形状和拓扑会发生很大变化，例如某个多边形因狭窄通道丢失而破碎、两个相邻的多边形分离。congalton(1997)探讨了多尺度栅格化后多边形的面积和形状变化，指出保留多边形的形状相比仅仅维持多边形存在需要更小的栅格像元。galton(2003)详细分析了栅格化对多边形面积、形状、拓扑、位置与方向等特征的影响，其中根据对称性和凸性讨论形状损失，根据区域连通性判断拓扑差异。zhou等(2017)比较了不同类型的拓扑变化及其形成原因，然后提出一种保留拓扑的栅格化算法。

基于简单的中心点法或最大面积法的常规栅格化算法仍然应用广泛(longleyetal.,2005)，大量栅格化误差研究表明多边形的空间特征在转换过程中损失严重。少数学者致力于保留多边形的面积或拓扑，忽视了该方法会给其他特征造成更大破坏。面积、形状及拓扑是多边形重要的空间特征，本发明旨在开发一种矢量-栅格转换中平衡保留这三种特征的算法，以达到综合降低栅格化误差的目标。评价了多种栅格像元大小下算法的精度和性能，并与修复单一误差的情景进行比较。分析了常规栅格化结果及其修复结果中各类型土地的面积变化。最后讨论不同的像元大小和矢量数据对算法的影响，算法的一些应用场景以及今后的工作。

技术实现要素：

本发明要解决的技术问题是：克服现有的技术缺点，提出一种平衡面积、拓扑和形状特征的矢量多边形栅格化方法。

为了解决上述技术问题，本发明提出的技术方案是：平衡面积、拓扑和形状特征的矢量多边形栅格化方法，将多边形p1和p2分别转换为区域集r＝{r1,r2,...,rn}和区域r2，执行以下步骤修复栅格化误差：

步骤1、获取多边形的边界像元，作为待分类的候选像元；

步骤2、若或n>1，则进行形状修复，优先保障多边形p1的栅格区域存在且连通，得到区域r1，具体执行步骤如下：

获取多边形p1损失的边界像元，作为候选集po；

若即不存在栅格区域，则从候选集po中取出一个像元分配给p1；

若n>1，即栅格区域不连通，则搜索多边形p1断开部位的像元，将它们从po中取出后分配给p1；

识别多边形p1丢失的凸起部位，从po中取出该部位的像元分配给p1；

至此，获得多边形p1所对应的区域r1；

步骤3、若区域r1与区域r2分离而多边形p1与多边形p2邻接，或者区域r1与区域r2邻接而多边形p1与多边形p2分离，则修复区域的拓扑关系，具体步骤如下：

若区域r1与区域r2分离而多边形p1与多边形p2邻接，获取多边形公共边界处损失的像元，作为候选集po；分别以多边形p1、多边形p2占有的像元作为起止点，基于最短路径思想在po中搜索像元并分配给与之重叠的多边形p1或多边形p2；

若区域r1与区域r2邻接而多边形p1与多边形p2分离，查找造成区域邻接的像元对，作为候选集pc＝{(p11,p12),(p21,p22),...,(pn1,pn2)}，将pc中每对像元中的一个分配给多边形p1、多边形p2之外的其他多边形；

步骤4、若多边形狭长的凸起部位丢失，将该部位的边界像元重新分配给该多边形；

步骤5、结合多边形面积与栅格像元大小，利用其余边界像元协调区域r1、r2的面积，使得多边形面积与其栅格区域的面积差异不大于像元面积的一半。

本发明根据矢量多边形的真实拓扑修复栅格区域的拓扑变化。保留多边形狭窄通道、凸起与凹陷部位的栅格像元，可以减少形状特征的损失。比较相邻栅格区域与原始多边形的面积差异，重新协调栅格化面积。实验结果表明，本发明方法相比常规方法能够有效保留多边形的面积、形状和拓扑特征，总体栅格化误差明显降低。

附图说明

下面结合附图对本发明作进一步的说明。

图1是实验区数据图(s＝10m)，其中，(a)为矢量数据图，(b)是使用平衡策略得到的栅格数据图。

图2是栅格化结果的比较图(s＝10m)，其中，(a)为矢量多边形；(b)为初始栅格化结果；(c)-(e)为三种常规策略处理结果图；(f)平衡策略的处理结果图。

具体实施方式

下面根据附图详细阐述本发明，使本发明的技术路线和操作步骤更加清晰。

本发明平衡面积、拓扑和形状特征的矢量多边形栅格化方法，将多边形p1和p2分别转换为区域集r＝{r1,r2,...,rn}和区域r2，执行以下步骤修复栅格化误差：

步骤1、获取多边形的边界像元，作为待分类的候选像元；

步骤2、若或n>1，则进行形状修复，优先保障多边形p1的栅格区域存在且连通，得到区域r1，具体执行步骤如下：

获取多边形p1损失的边界像元，作为候选集po；

若即不存在栅格区域，则从候选集po中取出一个像元分配给p1；

若n>1，即栅格区域不连通，则搜索多边形p1断开部位的像元，将它们从po中取出后分配给p1；

识别多边形p1丢失的凸起部位，从po中取出该部位的像元分配给p1；

至此，获得多边形p1所对应的区域r1；

步骤3、若区域r1与区域r2分离而多边形p1与多边形p2邻接，或者区域r1与区域r2邻接而多边形p1与多边形p2分离，则修复区域的拓扑关系，具体步骤如下：

若区域r1与区域r2分离而多边形p1与多边形p2邻接，获取多边形公共边界处损失的像元，作为候选集po；分别以多边形p1、多边形p2占有的像元作为起止点，基于最短路径思想在po中搜索像元并分配给与之重叠的多边形p1或多边形p2；

若区域r1与区域r2邻接而多边形p1与多边形p2分离，查找造成区域邻接的像元对，作为候选集pc＝{(p11,p12),(p21,p22),...,(pn1,pn2)}，将pc中每对像元中的一个分配给多边形p1、多边形p2之外的其他多边形；

步骤4、若多边形狭长的凸起部位丢失，将该部位的边界像元重新分配给该多边形；

步骤5、结合多边形面积与栅格像元大小，利用其余边界像元协调区域r1、r2的面积，使得多边形面积与其栅格区域的面积差异不大于像元面积的一半。

步骤5具体执行过程如下：

a1、计算多边形p与栅格化后得到的区域r之间的面积差，并据此确定需要重分类的像元数量n，属性值相同的邻接像元组成一个区域r；

式中，sp为多边形p的面积，sr为区域r的面积，s为栅格像元的面积；

a2、获取多边形p的边界像元，将其中属于区域r的像元作为候选集pa，其他像元作为候选集po；

a3、若n>0，从候选集pa中取出一个像元分配给其他多边形；若n<0，从候选集po中取出一个像元分配给多边形p；

a4、重复步骤a3，直至区域r及其邻接区域的平均面积误差不再降低或候选集po为空。

栅格化完成后进行栅格化误差评价，分别利用下式计算面积误差ae、形状误差se和拓扑误差te，

其中，av、ar为多边形栅格化前后的面积，max(av,ar)表示取av与ar两者中的较大值；cv、cr分别为多边形栅格化前后的紧凑度，若多边形转换后消失或破碎为多个小区域，则|cr-cv|＝1.0；

nr为真实发生的拓扑变化的多边形对或多边形个数，nv为可能发生的拓扑变化的多边形对或多边形个数；

多边形的紧凑度计算公式如下：

式中，a和p分别为多边形的面积和周长，c取值范围为(0,1]。

可能发生拓扑变化的多边形对或多边形个数总数根据下式求得：

其中i和j表示两个多边形，d(i,j)表示多边形i与多边形j边界的最短距离，s为栅格像元的面积。

若初始栅格化结果的误差为ae0、se0和te0，修复结果的误差为ae1、se1和te1，利用下式评价算法的总体修复能力ofc：

试验数据为中国常州市魏村街道的土地利用数据，来源于2018年常州市土地利用变更调查成果。该矢量数据的格式为esrishapefile，地理位置如图1所示。图层中有14,000个多边形，总面积为7,462公顷。包括6种土地利用类型：耕地、园地、水域、建设用地、交通用地和其他用地。表1列出了不同类型土地的多边形数量、平均面积以及平均紧凑度。

表1不同类型土地的多边形数量、平均面积以及平均紧凑度

算法使用python编程语言实现。地理空间数据抽象库(geospatialdataabstractionlibrary，gdal)提供了一系列空间类和空间操作函数，用于矢量和栅格数据的读取、写入、转换与处理。调用已有的矢量-栅格转换算法获取初始栅格化结果。矢量数据中每个多边形拥有唯一的id值，栅格化过程中将该值赋予栅格像元。根据多边形及其初始栅格区域的面积，计算最小化其面积误差所需重分类像元的最佳数量。为了高效地查询邻接多边形，针对试验数据设计并应用了四叉树空间索引(manolopoulosetal.,2018)。基于bresemhan算法的基本思想(gaol,2013)，计算多边形边界经过的所有像元。通过查找边界像元及其4-邻居获取栅格区域的邻接状态，用于识别已发生的拓扑变化。此外，重分类栅格像元时避免引入新的拓扑变化。

算法运行在windows10平台上，配有intel(r)core(tm)i7-6700处理器(时钟频率3.4ghz)、8gb内存和1tb硬盘。

在多种栅格化尺度下开展对照试验。基于中心点规则获取初始栅格化结果。分别利用三种常规策略(下文简称为面积策略、形状策略和拓扑策略)和平衡策略修复对应类型的栅格化误差。详细比较并分析了这些修复结果。同时评价了算法的精度和性能，以及修复误差后不同类型土地的面积变化。

以10m的栅格像元为例，对图1(a)的矢量多边形进行栅格化。三种常规策略以及平衡策略得到四种修复后的栅格化结果。图1(b)为平衡策略的结果。

如图2所示，在各栅格化结果中选择同一处位置，用于比较不同策略对多边形及其邻居的空间特征的作用。图2(a)为矢量图层，多边形p2的形状近似三角形。图2(b)为初始栅格化结果，多边形p1转换为两个不连通的区域，面积误差为2.72％；区域r2的形状近似矩形。图2(c)中减少了面积变化，多边形p1的面积误差降至0.02％，但导致区域r1与r4邻接、r5与r6分离。在图2(d)中，区域r2保留了多边形p2的一个凸角，但它与r3邻接；p1、p2的面积误差增加。图2(e)中修复了拓扑变化，但多边形p1的面积误差仍有2.61％，并且倾斜的狭窄通道变得竖直。图2(f)为平衡策略的结果，多边形p1的面积误差降至0.02％；p1、p2的狭窄通道和凸角损失减少；区域的拓扑与多边形的一致。以上比较表明，平衡策略对保留面积、形状和拓扑的综合效果较好。

使用5、10、15、20、25、30m六种像元尺寸开展栅格化试验。表2记录了常规策略和平衡策略得到的栅格化结果的面积误差ae、形状误差se及拓扑误差te。计算了反映各策略对栅格化误差的修复能力的指标ofc。

结果表明，栅格化误差对像元大小敏感。随着像元增大，ae、se和te均显著增大。常规策略对单一类型栅格化误差的修复效果较好。但面积策略导致形状误差和拓扑误差增大，像元尺寸较小时尤其明显。形状策略或拓扑策略有助于减小形状误差和拓扑误差，但都引入了新的面积误差。平衡策略可以同时降低这三种栅格化误差。ofc指标表明，随着像元增大，面积策略和形状策略的修复能力有所增长，拓扑策略和平衡策略的修复能力逐渐降低。主要原因是过于粗糙的栅格像元使多边形的空间特征损失严重。许多拓扑变化难以修复，并且不同误差之间的相互影响加强。综合来看，尽管栅格化误差始终存在，但平衡策略的结果明显优于其他策略。当像元尺寸小于10m时，一半以上的误差可被修复。

表2不同栅格像元大小下常规策略和平衡策略的栅格化误差

矢量-栅格转换会破坏多边形的面积、形状和拓扑等空间特征，从而引入多种类型的栅格化误差。本发明设计了三种单独修复面积误差、形状误差和拓扑误差的常规策略，并在此基础上开发一种可以同时修复这三种误差的平衡策略。在不同的栅格像元大小下，平衡策略能够更加有效地修复多种类型的栅格化误差。常规策略在修复某一种误差的同时会增大其他类型的误差。当像元尺寸小于10m时，平衡策略的结果中超过50％的栅格化误差被修复，明显优于常规策略的结果，但前者需要较长的执行时间。

除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。