一种专题地图地理面要素自动提取和矢量化方法

本发明涉及一种专题地图地理面要素自动提取和矢量化方法，属于计算机技术应用领域。

背景技术：

专题地图作为一种流行的数据可视化形式被广泛使用，随着互联网技术的不断发展和地理数据的不断丰富，互联网平台上发布的大量专题地图涵盖了极其丰富的地理信息和地理内容，快速高效地提取专题地图的面状地理要素及其矢量化工作可以为地理空间智能任务构建相关的数据集，为地理知识图谱构建、时空数据挖掘等技术的应用提供更加广泛的数据来源，具有巨大的数据潜力和发展前景。

目前互联网中专题地图面状地理要素的提取和矢量化过程主要存在的问题有：

1)专题地图通常以位图发布，无法直接获取相关地理要素的矢量化信息；

2)图像质量参差不齐，只能访问像素级信息，通常存在图像失真、边缘模糊等问题；

3)互联网中专题地图的场景内容、尺度范围、图例标注不统一；

4)专题地图中地理要素类型通常不止一类，面要素的提取会受线要素或点要素的影响；

5)专题地图中包含的文字标注、图形标注等地图要素，会对面要素的提取造成干扰；

现有的专题地图矢量化方法主要分为两类，一类是利用专业地理软件进行矢量化，但是对图像质量要求很高，且需要大量的人工操作，不同尺度和内容的图像需要设置不同的参数，矢量化过程复杂；第二类是利用图像处理方法实现专题地图的矢量化，但此类方法仅对某一类图像有效，无法同时实现多种专题地图的地理要素提取和矢量化。

以上两类矢量化方法在对互联网中多尺度多内容专题地图进行矢量化时的不足之处在于：

1)普适性差，方法应用场景单一，仅适用于同一尺度和内容的专题地图；

2)矢量化效率低，对于多尺度多内容专题地图的处理过程自动化程度低，并且矢量化精度和处理时长通常不能同时兼顾；

3)抗干扰性差，对于存在线要素、点要素、地图标注、行政区划线等的地图，在对面要素的抽取和矢量化过程中通常无法排除以上要素的干扰；

4)无法满足低质量图像的矢量化要求，仅对清晰的高质量图像处理效果好，不能正确获取互联网平台上大量的低质量专题地图的矢量化信息。

随着互联网专题地图的日益丰富以及各行业对于地理知识集成、挖掘、智能化等需求的不断增长，如何快速高效的实现互联网中多尺度多内容专题地图中地理面要素的自动获取和矢量化问题仍然亟待解决。

技术实现要素：

针对现有技术的不足，本发明的目的在于提供一种多尺度多内容的专题地图地理要素的自动提取和矢量化方法，具有效率高、普适性强、自动化程度高、抗干扰性强等优点，同时能够有效改进低质量专题地图的内容表达，增强其利用价值和地理可读性，能够快速有效地提取互联网中各类专题地图面状地理要素，并获取其矢量化数据。

本发明解决其技术问题所采用的技术方案是：

一种专题地图地理要素的自动提取和矢量化方法，包括以下步骤：

第一步，获取专题地图并对图像预处理，判断每一专题地图中是否包含地理投影和空间参考信息，若没有则对该专题地图进行地理配准；获取的专题地图包括各种不同尺度、内容或质量的地图。

第二步，自适应图像聚类，首先是进行图像类别预判，根据图像大小和图像的灰度统计特征分别计算步长s和峰值，其中，图像的灰度统计特征是图像灰度分布图的特征，图像灰度分布图即为图像的灰度统计，灰度分布图中峰值的个数即为预判的类别数k；然后将专题地图rgb三通道的像素值集合以及类别数k作为输入，利用k-means算法对图像进行聚类，并对聚类后图像(即聚类图像)进行中值滤波以去除噪声点；自适应图像聚类，能够改善低质量专题地图，一定程度上解决了图像模糊、失真的问题，且不受专题地图的尺度、内容、图例标注等的影响。

第三步，干扰要素自动检测和去除，首先利用canny边缘检测方法获取聚类图像中所有要素的边界(包含点要素、线要素和面要素以及地图标注)；然后利用图像膨胀扩展边界像素点的范围，获取当前图像所有边缘像素的位置索引；利用滑动窗口和位置索引，在聚类图像的基础上，遍历所有的边界像素点，并对其周围的s×s像素进行灰度值统计，将出现频率最高的灰度值赋值给新图像同位置的中心像素点；遍历结束后，对于新图像未赋值的像素点，其像素值与聚类图像同位置的像素相同；最后对新图像进行中值滤波。干扰要素自动检测和去除，可以排除专题地图中点要素、线要素以及地图标注等的影响，正确抽取地图的面要素，抗干扰性强。

第四步，图像矢量化，首先根据新图像的像素灰度值生成各类的二值图，然后根据每一类二值图的像素点的灰度值和空间位置的连续性生成表示面要素的多边形边界，并根据第一步的地理投影和空间参考信息为每个多边形对象创建相应的投影和空间参考，最后输出专题地图各类面要素的矢量文件。根据各类面要素多边形边界严格生成相应的矢量文件，拓扑关系正确，无需进行拓扑关系检查和修正，矢量化效率高。

所述第二步中，自适应图像聚类对一个图像内的信息进行聚类，其方法流程如下：

(1)输入待处理图像img；

(2)对img进行灰度转换，获得其灰度图像imgg；

(3)统计imgg的灰度特征，得到对应的灰度分布图，横轴代表每个灰度值dn的分布(dn＝0,1,…,255)，纵轴表示对应灰度值的像素点的总数；

(4)计算图像imgg的步长s:

其中，符号代表向下取整，n代表图像imgg中像素点的总数，pi表示图像imgg中第i个像素的灰度值，μ代表imgg中所有像素灰度值的均值。

(5)计算imgg图像类别的个数(灰度值集中且相近的可分成一类)：

1)k＝0；

2)根据灰度值及其对应的像素点个数，对每个灰度值(dn＝0,1,...,255)进行峰值判断：

当peak(dn)＝1时，k＝k+1,

其中，tdn代表灰度图像imgg中灰度值为dn的像素点的个数，s是步骤(4)计算的步长。在遍历所有的dn之后，最终得到的k值即为图像imgg的类别数。

(6)对图像img进行聚类：

1)将图像img各像素点的rgb值作为输入集合：{p0,p1,…,pn}，即{(r0,g0,b0)(r1,g1,b1),…,(rn,gn,bn)}，其中n代表图像img像素点的总数；

2)从输入集合{p0,p1,…,pn}中随机选取k个点作为初始聚类中心(c1,c2,…,ck)，即{(r0,g0,b0)(r1,g1,b1),…,(rk,gk,bk)}，分别对应着k个类别(c1,c2,…,ck)；

3)计算各点到中心点的距离,并将其归类至距离最短的类别内：

当dist(pi,cj)＝min(dist(pi,c1),dist(pi,c2),...,dist(pi,ck))时，cj＝cj∪pi，

其中，0≤i≤n,0≤j≤k。

4)更新中心点(c1,c2,…,ck)：

5)重复步骤3)和4)，直至中心点不再变化。

(7)对每一类像素重新赋值，不同类别的像素赋予不同的灰度值，形成聚类后图像imgk；在聚类后图像中，共有k个灰度值，不同的灰度值代表不同的类别；

(8)对聚类后图像imgk进行中值滤波，去除噪声点，输出最终灰度图像imgc。

所述第三步中，干扰要素的自动检测和去除方法流程如下：

(1)对imgc进行canny边缘检测，生成二值图像imgcanny，其灰度值取值仅有两个，0和1，灰度值为1的像素代表边缘像素，灰度值为0的像素代表非边缘像素；

(2)对imgcanny进行图像膨胀操作，扩展边缘像素的范围，得到新的二值图像imgdilate：

其中，“⊕”代表膨胀运算，d代表膨胀操作的卷积核，p代表imgcanny中的像素点，

将imgcanny与d进行卷积运算，即可扩展边缘像素的范围。

(3)根据步骤(2)中获取的边缘像素，利用滑窗统计对imgc中同位置的像素重新赋值，从而消除图像imgc中的干扰要素(线状要素和标注要素)，最后进行中值滤波，获得图像imgl：

1)获取imgdilate中表示边缘要素的像素索引，并建立imgc中同索引的像素点集合：

edge＝{p(m,n)|dn(pi(m,n))＝1,pi(m,n)∈imgdilate,p(m,n)∈imgc}

其中，(m,n)表示像素点的索引，pi(m,n)和p(m,n)分别为二值图imgdilate和灰度图imgc中索引为(m,n)的像素点，dn(pi(m,n))表示图像imgdilate中相应像素点的灰度值。

2)创建灰度图imgn：

imgn＝imgc

3)对灰度图imgc进行滑窗统计，并将统计结果赋值给imgn的同位置中心像素点：

dn(pn(i,j))＝dn(i,j)

其中，s是第二步计算得到的步长，maxdn()表示该窗口内统计次数最多的灰度值，dn(pn(i,j))表示imgn中索引为(i,j)的像素点的灰度值。

4)对步骤3)的结果进行中值滤波去除噪声点，获得灰度图像imgl。

与现有技术方案相比，本发明的有益效果是：

(1)方法普适性强，相比于仅针对同一尺度或同一类型的传统专题地图的矢量化方法，本发明的应用场景更加广泛，基本上不受地理范围、主题内容、地图形式、其他地理要素等因素的限制；

(2)方法的抗干扰性强，能够排除地图标注(文字、数字等)、要素边界线、线要素(河流、道路等)和点要素等的干扰，准确识别并提取面要素；

(3)矢量化效率高，自动化程度高，避免了传统方法中参数设置、拓扑修正等操作，能够快速准确的提取面状地理要素并生成相应的矢量化数据；

(4)能够保证低质量专题地图的处理效果和精度，实现互联网平台中低质量专题地图中面状地理要素的自动提取和矢量化。

附图说明

图1为本发明方法和传统矢量化方法的流程示意图。

图2为本方法的总体流程图。

图3为自适应图像聚类方法流程图。

图4为干扰要素自动检测和去除方法流程图。

图5为图像矢量化方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清晰明了，以下结合附图说明，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明的基本思路在于，根据专题地图的分布特征和颜色特征，利用图像处理技术和计算机技术，消除地图内的干扰要素，准确识别和提取同类地理面要素，生成拓扑关系正确的矢量文件，实现专题地图面要素的自动化或半自动化提取和快速矢量化。如图1所示，相比于传统矢量化方法对于图像质量、方法参数或者人工操作的强依赖，本发明方法能够通过图像处理技术增强低质量专题地图的表达，降低其模糊度，正确提取其地理要素信息；根据图像特征自动调整步长、聚类数等参数，无需针对地图的尺度或内容反复调整相关参数；严格根据处理后图像的像素边界进行矢量化，拓扑关系正确，免除后续的手动调整过程。

本发明方法的总体流程图如图2示，首先判断文件是否包含地理信息，若没有则进行手动地理配准，然后提取文件的地理信息和图像信息，将所有像素点整理成集合{p0,p1,…,pn}，集合元素为像素点在rgb各波段的值，即第i个像素点pi＝(ri,gi,bi)；通过图像的灰度转换，将彩色图转为灰度图，根据其灰度分布特征和像素个数计算步长s；再由灰度分布特征和步长s判断图像峰值，峰值个数即为聚类数k；将像素点集合{p0,p1,…,pn}和聚类数k作为输入参数，对图像进行k-means聚类，并对结果进行中值滤波，各类要素即可根据灰度值进行区分；利用canny边缘检测算法识别聚类图像的线要素，并通过图像膨胀算法扩大线要素的标注范围，之后对聚类图像进行滑窗统计，即以有标注的像素为中心，以聚类图像为基础，对其周围的像素点进行统计，将出现次数最多的灰度值赋给中心像素，所有标注像素点完成赋值后，再对其结果图像进行中值滤波，以上步骤可检测并去除图中的线要素，避免干扰面要素的识别和提取；将上一步获取的最终图像转换为各类的二值图，根据文件的地理信息为各图创建地理空间参考信息，并严格按照像素值进行矢量化，创建相应的矢量文件，最终实现专题地图中各类面要素的矢量化。

本发明方法具体步骤如下所示：

第一步，图像预处理

判断文件是否包含地理投影信息和空间参考信息，如果有，则提取其地理信息和图像波段信息；如果没有，则需通过专业软件进行地理配准，为文件创造地理投影信息和空间参考，之后再提取其地理信息和图像信息。

第二步，自适应图像聚类

图像聚类算法需要两个输入参数，分别是聚类数和待分类数据集。如图3所示，根据灰度转换后图像大小和灰度统计特征计算步长s，再由步长s判断图像灰度直方图的峰值，峰值个数即为聚类数；把所有彩色图像的像素点输出为集合{p0,p1,…,pn}，集合元素为像素点在rgb各波段的值，即pi＝(ri,gi,bi)；最后将像素点集合{p0,p1,…,pn}和聚类数k作为输入参数，对图像进行k-means聚类，并进行中值滤波。通过自适应图像聚类，获取专题地图各像素点的聚类结果，能够降低图像模糊度，可以改善低质量专题地图，提升其面要素识别和的准确性。

自适应图像聚类的具体步骤如下：

(1)输入待处理图像img；

(2)对img进行灰度转换，获得其灰度图像imgg；

(3)统计imgg的灰度特征，得到对应的灰度分布图，横轴代表每个灰度值dn的分布(dn＝0,1,…,255)，纵轴表示对应灰度值的像素点的总数；

(4)计算图像的步长s:

其中，符号代表向下取整，n代表图像像素点的总数，pi表示图像第i个像素的灰度值，μ代表所有像素灰度值的均值。

(5)计算imgg图像类别的个数：

1)k＝0；

2)根据灰度值及其对应的像素点个数，对每个灰度值(dn＝0,1,...,255)进行峰值判断：

当peak(dn)＝1时，k＝k+1,

其中，tdn代表灰度图像中灰度值为dn的像素点的个数，s是步骤(4)计算的步长。

在遍历所有的dn之后，最终得到的k值即为图像的类别数。

(6)对图像img进行聚类：

1)将图像每个像素点的rgb值作为输入集合：{p0,p1,…,pn}，即{(r0,g0,b0)(r1,g1,b1),…,(rn,gn,bn)}，其中n代表图像像素点的总数；

2)从输入集合{p0,p1,…,pn}中随机选取k个点为初始聚类中心(c1,c2,…,ck)，即{(r0,g0,b0)(r1,g1,b1),…,(rk,gk,bk)}，分别对应着k个类别(c1,c2,…,ck)；

3)计算各点到中心点的距离,并将其归类至距离最短的类别内：

当dist(pi,cj)＝min(dist(pi,c1),dist(pi,c2),...,dist(pi,ck))时，cj＝cj∪pi，

其中，0≤i≤n,0≤j≤k。

4)更新中心点(c1,c2,…,ck)：

5)重复步骤3)和4)，直至中心点不再变化

(7)对每一类像素重新赋值，不同类别的像素赋予不同的灰度值，形成聚类后图像imgk；

(8)对聚类后图像imgk进行中值滤波，去除噪声点，输出最终灰度图像imgc。

第三步，干扰要素自动检测和去除

由于专题地图的内容和形式多样化，针对专题地图地理面要素的提取通常受到线要素或点要素的干扰，例如水系、河流、道路、行政边界、地图标注等。为保证面要素提取结果的正确性，需要对以上要素进行自动检测和去除。如图4所示，对由第二步获取的聚类图像进行干扰要素的检测和去除，首先利用canny边缘检测方法获取专题地图中所有要素的边界；然后利用图像膨胀扩展边界像素点的范围，获取所有边缘像素的位置索引；利用滑动窗口和位置索引，在聚类图像的基础上，遍历所有的边界像素点，并对其周围的s×s像素进行灰度值统计，将出现频率最高的灰度值赋值给新图像imgn同位置的中心像素点；遍历结束后，对于新图像未赋值的像素点，其像素值与聚类图像同位置的像素相同；最后进行中值滤波。

干扰要素的自动检测和去除方法流程如下：

(1)对imgc进行canny边缘检测，生成二值图像imgcanny，其灰度值取值仅有两个，0和1，灰度值为1的像素代表检测出的边缘像素，灰度值为0的像素代表非边缘像素；

(2)对imgcanny进行图像膨胀操作，扩展边缘像素的范围，得到新的二值图像imgdilate：

其中，“⊕”代表膨胀运算，d代表膨胀操作的卷积核，p代表imgcanny中的像素点，将imgcanny与d进行卷积运算，即可扩展边缘像素的范围。

(3)根据步骤(2)中获取的边缘像素的索引，利用滑窗统计对imgc中同索引的像素重新赋值，从而消除图像干扰要素(线状要素和标注要素)，最后进行中值滤波，获得图像imgl：

1)获取imgdilate中表示边缘要素的像素索引，并建立imgc中同索引的像素点集合：

edge＝{p(m,n)|dn(pi(m,n))＝1,pi(m,n)∈imgdilate,p(m,n)∈imgc}

其中，(m,n)表示像素点的索引，pi(m,n)和p(m,n)分别为二值图imgdilate和灰度图imgc中索引为(m,n)的像素点，dn(pi(m,n))表示图像imgdilate中相应像素点的灰度值。

2)创建灰度图imgn：

imgn＝imgc

3)对灰度图imgc进行滑窗统计，并将统计结果赋值给imgn的同位置中心像素点：

dn(pm(i,j))＝dn(i,j)

其中，s是第二步计算得到的步长，maxdn()表示该窗口内统计次数最多的灰度值，dn(pn(i,j))表示imgn中索引为(i,j)的像素点的灰度值

4)对步骤3)的结果进行中值滤波去除噪声点，获得灰度图像imgl。在imgl中，不同的灰度值代表不同的类别，根据灰度值可获取各个类别的二值图。

第四步，图像矢量化

如图5所示，经自适应聚类和干扰要素去除等处理后，可获取各类面要素的灰度图像，首先从灰度图像imgl中根据灰度值分别提取出各类面要素的二值图，然后根据像素点的灰度值和空间位置的连续性生成面要素的多边形边界，并根据第一步获取的地理投影和空间参考信息为每个多边形对象创建相应的投影和空间参考，最后输出专题地图各类面要素的矢量文件。

为证明本发明方法的有效性，选取多个不同尺度和内容的专题地图作为实验数据，其中分布包含文字标注、线状要素(行政边界、水系、道路等)以及图例样式(离散、连续)等要素。利用本发明方法对实验数据进行处理，同时采用专业地理软件(arcmap转换工具)进行对比实验，实验结果如表1所示。

表1

以上实施例仅用以说明本发明的技术方案而非对其进行限制，本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明的原理和范围，本发明的保护范围应以权利要求书所述为准。