CN103824283B - 基于随机概率采样与多级融合的大幅面地形图分割方法 - Google Patents

Abstract

本发明公开了一种基于随机概率采样与多级融合的大幅面地形图分割方法：步骤1，输入原始地形图图像；步骤2，对原始地形图图像进行随机概率采样；步骤3，确定地形图分割种类数；步骤4，计算聚类中心，得到聚类中心矩阵；步骤5，计算原始地形图的隶属矩阵，得到隶属矩阵；步骤6，利用隶属矩阵对原始地形图进行模糊分类，得到分割后的图像；步骤7，对分割后的图像进行多级融合，得到的融合后图像；步骤8，输出分割后的分版图像。该方法将随机概率采样与图像多级融合相结合以进行地形图分割，最终既能够更准确地获得地形图的分割版图，又能大幅度提高分割效率。尤其适用于大幅面的地形图分割。

Description

基于随机概率采样与多级融合的大幅面地形图分割方法

技术领域

本发明属于地理信息处理技术领域，具体涉及一种基于随机概率采样与多级融合的大幅面地形图分割方法。本发明应用于大幅面地形图的图像分割。

背景技术

在地形图处理领域中，为了快速准确地对地形图进行分割得到分割版图，并为地形图中地理要素的进一步提取提供有利条件，通常需要对地形图分割。目前，地形图分割方法主要采用基于颜色空间的聚类和基于样本学习的分类法。

郑华利等人在文献“郑华利，周献中，王建宇，‘基于色彩空间转换及模糊约束聚类的地形图自动分色’[J]，测绘学报，32(2)，183-187(2003)”中提出了一种模糊约束聚类算法。该方法首先对地形图的色彩空间进行转换，使其直方图的峰谷特性更加明显，然后对模糊C均值聚类算法进行模糊约束，从而较好地解决了地图的自动分色问题。由于经过色彩空间的转换后直方图的峰值比较明显，所以该算法对于初始聚类中心选择的依赖性较小。此外，采用了模糊约束的颜色聚类，对于那些在地形图直方图中占有较小比例的颜色信息也能正确地进行分层。该方法虽然能够得到较好的地形图分割效果，但是仍然存在的缺点是，模糊C均值聚类算法对地形图进行4分类，没有考虑地形图中存在的过渡色、渐变色和覆盖色等复杂问题，导致分割结果不准确。此外，该方法对地形图中所有数据进行一次性处理，没有考虑到大幅面地形图数据量大的问题，因此分割效率不高。

Alireza Khotanzad等人在文献“Alireza Khotanzad and Edmund Zink，‘Contour Line and Geographic Feature Extraction from USGS Color TopographicalPaper Maps’[J]，IEEE Transactions on Pattern Analysis and MachineIntelligence，25(1)，18-31(2003)”中提出了一种基于样本学习的分类方法对地形图进行分割。该方法首先根据地形图中不同地理要素之间的重叠覆盖情况，制作相应的模板，然后根据模板的颜色特征与地形图中像素的颜色特征之间的差异对地形图进行分割。此方法虽然能够较好的处理面状地理要素信息以及地形图中线划要素信息。但是仍然存在的缺点是，地形图中不同地理要素的模板样本难以制作齐全，此外，图像分割时仅仅使用颜色距离，难以很好的处理地形图中存在的过渡色、渐变色和覆盖色等复杂问题，这些因素都是导致最终分割结果不准确的原因，从而不能得到较好的地形图分割效果。

综上所述，基于颜色空间的聚类和基于样本学习的分类方法在地形图分割方面虽然获得了较好的效果，但这些方法都难以快速准确地进行大幅面地形图图像的分割，此外，在地形图分割种类数的确定方面没有一个准确的自适应的确定方法。所以在利用传统的基于颜色空间的聚类和基于样本学习的分类方法得到的分割版图中，存在分割不准确以及效率低下的问题。

发明内容

针对上述现有技术中基于颜色空间的聚类和基于样本学习的分类方法导致的地形图分割效果不理想和效率低下的缺点，本发明的目的在于，提出一种基于随机概率采样与多级融合的大幅面地形图分割方法，该方法中对地形图的分割主要包括获得聚类中心，模糊分类和多级融合三部分。本发明充分利用随机概率采样能够较大程度减少获得聚类中心所需要的数据量的优点，以及多级融合能够更准确地对分割后图像进行融合的特点，将随机概率采样与图像多级融合相结合以进行地形图分割，最终既能够更准确地获得地形图的分割版图，又能大幅度提高分割效率。尤其适用于大幅面的地形图分割。

为了达到上述目的，本发明采用如下技术方案予以解决：

一种基于随机概率采样与多级融合的大幅面地形图分割方法，具体包括如下步骤：

步骤1，输入原始地形图图像；

步骤2，对原始地形图图像进行随机概率采样，得到采样数据矩阵；

步骤3，确定地形图分割种类数；

步骤4，计算聚类中心，得到聚类中心矩阵；

步骤5，计算原始地形图的隶属矩阵，得到隶属矩阵；

步骤6，利用隶属矩阵对原始地形图进行模糊分类，得到分割后的图像；

步骤7，对分割后的图像进行多级融合，得到的融合后图像；

步骤8，输出分割后的分版图像。

进一步的，所述步骤2的具体步骤如下：

2a)生成一个符合正态分布的，且和原始地形图A等大小的随机矩阵R_M×N作为随机概率矩阵,该矩阵中所有元素的值为0和1之间的任意随机数；其中，M为原始地形图的宽度，N为原始地形图的长度；

2b)计算概率阈值P，并利用步骤2a)生成的随机概率矩阵R_M×N对原始地形图进行数据采样；概率阈值P的计算如下：

$P=\frac{{10}^5}{M\×N}$

其中，M为原始地形图的宽度，N为原始地形图的长度。

数据采样过程如下：遍历随机概率矩阵R_M×N，如果矩阵中某一元素的值大于等于概率阈值P，则原始地形图A中相同位置像素点的颜色信息被选择为采样信息；否则原始地图A中相同位置像素点的颜色信息为0；遍历结束后得到采样数据矩阵A'；

进一步的，所述步骤3中，利用下式确定地形图的分割种类数：

Cluster_n=C_l ¹+C_l ²+C_l ³+…C_l ^l+1. (5)

其中，Cluster_n为地形图的分割种类数，C_l ¹为排列组合表达式，表示在由l种颜色绘制的地形图中的某区域有1种颜色相互覆盖所得到的颜色种类数，l≥1。

进一步的，所述步骤4的具体步骤如下：

4a)初始化隶属矩阵：初始化一个与原始地形图A相同大小的随机矩阵作为初始隶属矩阵U₀；并使其满足下式：

$\underset{k=1}{\overset{\mathrm{Cluster}\_n}{\mathrm{\Σ}}}{u}_{0,k(i,j)}=1$

其中，u_0,k（i,j）为初始隶属矩阵U₀中坐标（i,j）位置第k类元素的值；

4b)利用下式计算得到矩阵大小为Cluster_n×3的聚类中心矩阵v_I，其中，I为迭代次数，当迭代次数I=1时，隶属度矩阵U=U₀，u_k（i,j）=u_0,k（i,j）；

$v_{I,k}=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left(u_{k(i,j)}\right)}^m{x}_{(i,j)}}{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left(u_{k(i,j)}\right)}^m}$

其中，v_I,k为第I次迭代中聚类中心矩阵v_I的第k个聚类中心，k=1,…,Cluster_n；m=2；x_(i,j)为原始地形图中坐标(i,j)位置的像素值；（u_k（i,j））^m为隶属矩阵U中坐标（i,j）位置第k类元素值的m次方；

4c)根据下式计算第I次迭代的总体价值函数J_I：

$J_I=\underset{k=1}{\overset{\mathrm{Cluster}\_n}{\mathrm{\Σ}}}{J}_{I,k}=\underset{k=1}{\overset{\mathrm{Cluster}\_n}{\mathrm{\Σ}}}\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left(u_{k(i,j)}\right)}^m{\left|\right|v_{I,k}-x_{(i,j)}\left|\right|}^2\right)$

其中，J_I,k为第I次迭代的第k类价值函数；

如果总体价值函数J_I小于阀值J_min=10^-5，或总体价值函数J_I相对上次得到的价值函数值J_I-1的改变量小于阀值J_{dif_min}=10^-5，则步骤4的计算聚类中心结束，最终得到聚类中心矩阵v＝v_I，然后执行步骤5；否则，进行步骤4d)；

4d)使用下式计算新的隶属度矩阵U'；

${u^{\′}}_{k(i,j)}=\frac{1}{\underset{p=1}{\overset{\mathrm{Cluster}\_n}{\mathrm{\Σ}}}{\left(\frac{\left|\right|v_{I,k}-x_{(i,j)}\left|\right|}{\left|\right|v_{I,p}-x_{(i,j)}\left|\right|}\right)}^{2/(m-1)}}$

其中，u'_k(i,j)为新的隶属矩阵U'中坐标（i,j）位置第k类元素的值；v_I,k和v_I,p分别表示聚类中心矩阵v_I中第k个和第p个聚类中心；

令U=U'，u_k（i,j）=u'_k（i,j），I=I+1，返回步骤4b)。

进一步的，所述步骤5的具体步骤如下：

5a)读取步骤(4)中得到的聚类中心矩阵v；

5b)利用模糊C均值聚类算法中的隶属度计算公式，计算原始地形图A中各像素的隶属度，得到所有像素的隶属度组成隶属矩阵U_F；

$u_{F,k(i,j)}=\frac{{(1/{\left|\right|v_k-x_{(i,j)}\left|\right|}^2)}^{1/(m-1)}}{\underset{p=1}{\overset{\mathrm{Cluster}\_n}{\mathrm{\Σ}}}{(1/{\left|\right|v_p-x_{(i,j)}\left|\right|}^2)}^{1/(m-1)}}$

其中，u_F,k（i,j）为隶属矩阵U_F中位置为（i,j）处第k类元素的值，对应于原始地形图中坐标位置的像素属于k类的隶属度；v_k和v_p分别为聚类中心矩阵v的第k个和第p个聚类中心；x_(i,j)为原始地形图中坐标(i,j)位置的像素值；m=2。

进一步的，所述步骤6的具体步骤如下：

利用隶属矩阵U_F，按以下公式所示的模糊分类方法对原始地形图A进行图像分割得到分割后的图像：

x_(i,j)∈image_k如果u_F,k（i,j）≥u_F,p（i,j）(p＝1,2,…Cluster_n)

其中，image_k为原始地形图中像素x_(i,j)被划入的第k个分割后的图像。

进一步的，所述步骤7的具体步骤如下：

7a)读取聚类中心矩阵；

7b)计算每个聚类中心和其他聚类中心之间的距离；

7c)利用聚类中心之间的距离以及多级融合方法对分割后的图像进行融合，直至融合后得到的分版图像数目和地形图的分割种类数Cluster_n相等时停止，得到融合后的分版图像；

进一步的，所述步骤7b）中的距离采用马氏距离：

${\mathrm{Dist}}_{k1,k2}=\sqrt{{(v_{k1}-v_{k2})}^T{\mathrm{\Σ}}^{-1}(v_{k1}-v_{k2})}$

其中，Dist_k1,k2为聚类中心矩阵v中第k₁个聚类中心v_k1和第k₂个聚类中心v_k2之间的马氏距离，()^T为矩阵的转置。

进一步的，所述步骤7c）中的多级融合方法的规则：当且仅当与第i类的最小距离是第j类，同时与第j类的最小距离是第i类时，将第i类和第j类融合为一新类。

与现有技术相比，本发明具有以下优点：

第一，在获得聚类中心时，本发明应用随机概率采样的方法对原始地形图进行数据采样，能够较大程度上减少计算聚类中心所需要的数据量，因而能够快速获得聚类中心，以提高分割效率。克服了现有技术中对大幅面地形图分割效率低下的缺点,使得使用本发明中的方法能够快速获得大幅面地形图的分割结果。

第二，本发明在确定分割种类数时，采用基于地形图颜色覆盖的自适应确定方法，能够根据地形图中颜色的种类数及其覆盖情况自适应确定图像分割种类数，因而能够得到更准确的分割种类数。

第三，本发明在对分割后图像进行融合处理时，应用了多级融合的方法。由于多级融合采用了基于马氏距离的相互融合规则，克服了现有技术中融合结果不准确的问题，使得使用本发明中的方法能够获得更准确的地形图分版图像。

附图说明

图1为本发明的流程图。

图2为本发明用于对小幅面地形图进行分割的仿真图。其中，(a)是原始小幅面地形图图像；(b₁)是本发明方法得到的绿版图；(b₂)是基于地形图空间信息的图像分割方法得到的绿版图；(b₃)是基于遗传模糊C均值聚类方法得到的绿版图；(b₄)是基于改进的模糊C均值聚类方法得到的绿版图；(c₁)是本发明方法得到的棕版图；(c₂)是基于地形图空间信息的图像分割方法得到的棕版图；(c₁’)是c₁局部放大图；(c₂’)是c₂局部放大图；(c₃’)是c₃局部放大图；(c₄’)是c₄局部放大图；(c₃)是基于遗传模糊C均值聚类方法得到的棕版图；(c₄)是基于改进的模糊C均值聚类方法得到的棕版图；(d₁)是本发明方法得到的黑版图；(d₂)是基于地形图空间信息的图像分割方法得到的黑版图；(d₁’)是d₁局部放大图；(d₂’)是d₂局部放大图，(d₃’)是d₃局部放大图；(d₄’)是d₄局部放大图；(d₃)是基于遗传模糊C均值聚类方法得到的黑版图；(d₄)是基于改进的模糊C均值聚类方法得到的黑版图；(e₁)是本发明方法得到的蓝版图；(e₂)是基于地形图空间信息的图像分割方法得到的蓝版图；(e₁’)是e₁局部放大图；(e₂’)是e₂局部放大图；(e₃’)是e₃局部放大图；(e₄’)是e₄局部放大图；(e₃)是基于遗传模糊C均值聚类方法得到的蓝版图；(e₄)是基于改进的模糊C均值聚类方法得到的蓝版图。

图3为本发明用于对大幅面地形图进行分割的仿真图。其中，(a)是原始地形图图像；（b）是本发明的方法得到的黑版图；（c）是本发明方法得到的蓝版图；（d）是本发明的方法得到的棕版图；（e）是本发明的方法得到的红版图；（f）是本发明的方法得到的绿版图。

以下结合附图和具体实施方式对本发明进一步解释说明。

具体实施方式

参照图1，本发明的基于随机概率采样与多级融合的大幅面地形图分割方法，具体包括如下步骤：

步骤1，输入原始地形图图像A

在计算机中应用matlab软件读取存储在计算机硬盘空间中的原始地形图A。

步骤2，随机概率采样

为了减少聚类算法在计算聚类中心的数据量，先对原始地形图进行随机概率采样，具体包括如下步骤：

2a)利用matlab软件生成一个符合正态分布的，且和原始地形图A等大小的随机矩阵R_M×N作为随机概率矩阵,该矩阵中所有元素的值为0和1之间的任意随机数。其中，M为原始地形图的宽度，N为原始地形图的长度。

2b)计算概率阈值P，并利用步骤2a)生成的随机概率矩阵R_M×N对原始地形图进行数据采样。概率阈值P的计算如下：

$P=\frac{{10}^5}{M\×N}$

其中，M为原始地形图的宽度，N为原始地形图的长度。

数据采样过程如下：遍历随机概率矩阵R_M×N，如果矩阵中某一元素的值大于等于概率阈值P，则原始地形图A中相同位置像素点的颜色信息（即该像素点的R,G,B值）被选择为采样信息；否则原始地图A中相同位置像素点的颜色信息为0；遍历结束后得到采样数据矩阵A'。即：

其中，A'(i,j)为采样后图像中坐标(i,j)位置的像素值，A(i,j)为原始地形图中坐标(i,j)位置的像素值，R(i,j)为随机概率矩阵R_M×N中坐标(i,j)位置的元素值，1≤i≤M,1≤j≤N。

步骤3，地形图分割种类数的自适应确定

利用地形图分割种类数的自适应方法确定地形图的分割种类数。

Cluster_n=C_l ¹+C_l ²+C_l ³+…C_l ^l+1. (5)

其中，Cluster_n为地形图的分割种类数，C_l ¹为排列组合表达式，表示在由l种颜色绘制的地形图中的某区域有1种颜色相互覆盖所得到的颜色种类数，l≥1，C_l ²、C_l ³、…C_l ^l同理依次类推。地形图中黑色和其他所有颜色混合都为黑色，因此最后加1。

步骤4，计算聚类中心，得到聚类中心矩阵

利用模糊C均值聚类算法对步骤2b)得到采样数据矩阵A'进行聚类处理，得到和分割种类数Cluster_n相同个数的聚类中心矩阵v，并存入计算机。v是一个大小为Cluster_n×3的矩阵。具体包括以下步骤：

4a)初始化隶属矩阵：利用matlab软件初始化一个与原始地形图A相同大小的随机矩阵作为初始隶属矩阵U₀，初始隶属矩阵U₀中所有元素的值为0和1之间的任意随机数。并使其满足下式：

$\underset{k=1}{\overset{\mathrm{Cluster}\_n}{\mathrm{\Σ}}}{u}_{0,k(i,j)}=1$

其中，u_0,k（i,j）为初始隶属矩阵U₀中坐标（i,j）位置第k类元素的值。

4b)利用下式计算得到矩阵大小为Cluster_n×3的聚类中心矩阵v_I，其中，I为迭代次数，当迭代次数I=1时，隶属度矩阵U=U₀，u_k（i,j）=u_0,k（i,j）；

$v_{I,k}=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left(u_{k(i,j)}\right)}^m{x}_{(i,j)}}{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left(u_{k(i,j)}\right)}^m}$

其中，v_I,k为第I次迭代中聚类中心矩阵v_I的第k个聚类中心，k=1,…,Cluster_n；m为权值参数，本发明中m=2；x_(i,j)为原始地形图中坐标(i,j)位置的像素值；（u_k（i,j））^m为隶属矩阵U中坐标（i,j）位置第k类元素值的m次方。

4c)根据下式计算第I次迭代的总体价值函数J_I：

$J_I=\underset{k=1}{\overset{\mathrm{Cluster}\_n}{\mathrm{\Σ}}}{J}_{I,k}=\underset{k=1}{\overset{\mathrm{Cluster}\_n}{\mathrm{\Σ}}}\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left(u_{k(i,j)}\right)}^m{\left|\right|v_{I,k}-x_{(i,j)}\left|\right|}^2\right)$

其中，J_I,k为第I次迭代的第k类价值函数。

如果总体价值函数J_I小于阀值J_min=10^-5，或总体价值函数J_I相对上次得到的价值函数值J_I-1的改变量小于阀值J_{dif_min}=10^-5，则步骤4的计算聚类中心结束，最终得到聚类中心矩阵v＝v_I，然后执行步骤5。否则，进行步骤4d)。

4d)使用下式计算新的隶属度矩阵U'；然后令U=U'，u_k（i,j）=u'_k（i,j），I=I+1，返回步骤4b)。

${u^{\′}}_{k(i,j)}=\frac{1}{\underset{p=1}{\overset{\mathrm{Cluster}\_n}{\mathrm{\Σ}}}{\left(\frac{\left|\right|v_{I,k}-x_{(i,j)}\left|\right|}{\left|\right|v_{I,p}-x_{(i,j)}\left|\right|}\right)}^{2/(m-1)}}$

其中，u'_k(i,j)为新的隶属矩阵U'中坐标（i,j）位置第k类元素的值；v_I,k和v_I,p分别表示聚类中心矩阵v_I中第k个和第p个聚类中心。

步骤5，计算原始地形图的隶属矩阵

5a)读取步骤(4)中得到的聚类中心矩阵v；

5b)利用模糊C均值聚类算法中的隶属度计算公式，计算原始地形图A中各像素的隶属度，得到所有像素的隶属度组成隶属矩阵U_F。

$u_{F,k(i,j)}=\frac{{(1/{\left|\right|v_k-x_{(i,j)}\left|\right|}^2)}^{1/(m-1)}}{\underset{p=1}{\overset{\mathrm{Cluster}\_n}{\mathrm{\Σ}}}{(1/{\left|\right|v_p-x_{(i,j)}\left|\right|}^2)}^{1/(m-1)}}$

其中，u_F,k（i,j）为隶属矩阵U_F中位置为（i,j）处第k类元素的值，对应于原始地形图中坐标位置的像素属于k类的隶属度；v_k和v_p分别为聚类中心矩阵v的第k个和第p个聚类中心；x_(i,j)为原始地形图中坐标(i,j)位置的像素值；m为权值参数，本发明中m=2；

步骤6，模糊分类，得到分割后的图像。具体操作如下：

利用步骤5b）中得到的隶属矩阵U_F，按以下公式所示的模糊分类方法对原始地形图A进行图像分割得到分割后的图像。

x_(i,j)∈image_k如果u_F,k（i,j）≥u_F,p（i,j）(p＝1,2,…Cluster_n)

其中，image_k为原始地形图中像素x_(i,j)被划入的第k个分割后的图像。

步骤7，多级融合，得到的融合后图像。具体操作如下：

7a)读取步骤(4)中得到的聚类中心矩阵v；

7b)利用下式计算每个聚类中心和其他聚类中心之间的距离，本发明中采用马氏距离公式：

${\mathrm{Dist}}_{k1,k2}=\sqrt{{(v_{k1}-v_{k2})}^T{\mathrm{\Σ}}^{-1}(v_{k1}-v_{k2})}$

其中，Dist_k1,k2为聚类中心矩阵v中第k₁个聚类中心v_k1和第k₂个聚类中心v_k2之间的马氏距离，()^T为矩阵的转置；

7c)利用步骤7b)中得到的聚类中心之间的马氏距离以及多级融合方法对步骤6得到的分割后的图像进行融合，直至融合后得到的分版图像数目和地形图的分割种类数Cluster_n相等时停止，得到融合后的分版图像；

定义：多级融合方法的规则是：当且仅当与第i类的最小距离是第j类，同时与第j类的最小距离是第i类时，将第i类和第j类融合为一新类。即如果同时满足以下两个条件：

Dist_i,j=min({Dist_i,C|C＝1,2，...，Cluster_n})

Dist_i,j=min({Dist_C,j|C＝1,2，...，Cluster_n})

则

map_c(i,j)=image_k1(i,j)+image_k2(i,j)

其中，map_c(i,j)为第c个融合后的图像中坐标(i,j)位置像素值，image_k1(i,j)和image_k2(i,j)分别为第k₁和第k₂个分割后图像中坐标(i,j)位置像素值。

步骤8，输出分割后的分版图像。

本发明的有益效果可以通过以下仿真进一步说明。

仿真试验1，该试验对本发明中小幅面地形图进行分割的仿真。

仿真条件：在MATLAB7.0软件下进行。

参照图2，分别使用本发明方法与基于地形图空间信息的图像分割方法、基于遗传模糊C均值聚类方法、基于改进的模糊C均值聚类方法等三种传统的地形图分割方法对大小为795×680像素的小幅面地形图进行仿真实验。图2(a)是原始小幅面地形图图像，主要使用四种颜色表示不同的地理要素信息，但是同时也存在多种渐变色，过渡色和覆盖色。

从图2(d₂,e₂，e₄)可以看出,传统图像分割算法得到的结果图像中存在大量的本不应该属于此版图的杂色。特别是在图2(e₄)中，蓝色分版图像中含有大量的黑色像素，导致在黑色分版图像中的线划要素出现较为严重的断裂的情况，而与之相对应的由本发明方法得到的分割结果则具有较好的效果。在图2(c₁’,c₂’,c₃’,c₄’),图2(d₁’,d₂’,d₃’,d₄’)和图2(e₁’,e₂’,e₃’,e₄’)中，这些图像时分割结果的细节图像，从这些图像中可以看出使用本发明方法得到的分割结果与其它方法得到的分割结果相比，结果图像中残留的噪声更少，线划要素更完整，更连续。

仿真试验2，该试验对本发明中大幅面地形图进行分割的仿真。

仿真条件：在MATLAB7.0软件下进行。

参照图3，对大小为5000×5000像素大幅面地形图进行仿真实验。从图3（d）中可以看出，本发明方法分割得到的棕色版图中基本不含其他杂色，且等高线信息完整，线划连续；从图3（f）中可以看出，本发明方法分割得到的绿色版图中面状区域平滑连续，边界清晰。此外，本发明方法对大小为5000×5000的地形图进行分割需要的时间为325.7302(s)，是传统方法所需时间1513.2350(s)的1/5。因此，本发明方法能够对大幅面地形图进行快速准确的分割。

分别对传统的图像分割方法和本发明方法计算分割效率评价参数，最终数据如表1所示。

表1不同大小地形图分割方法使用时间表

地形图大小	传统的图像分割方法	本发明方法
			500×500	32.1287(s)	13.8662(s)
600×600	44.1108(s)	19.6646(s)
			700×700	59.6015(s)	26.5989(s)
800×800	75.8195(s)	20.1525(s)
			900×900	96.3398(s)	26.0943(s)
1000×1000	122.2289(s)	31.3653(s)
			1100×1100	152.0572(s)	26.3407(s)
1200×1200	178.4432(s)	31.2745(s)
			1300×1300	205.7123(s)	36.5182(s)
1400×1400	240.4408(s)	42.1408(s)
			1500×1500	286.1583(s)	49.0366(s)
2000×2000	492.9253(s)	73.1945(s)
			2500×2500	805.2184(s)	112.6271(s)

从表1的客观评价测量值可以看出，本发明在地形图分割效率上优于传统的地形图分割方法，能够对大幅面地形图进行快速分割。且随着图像大小的增加，本发明的高效性更明显。

Claims (8)

1.一种基于随机概率采样与多级融合的大幅面地形图分割方法，其特征在于，具体包括如下步骤：

步骤1，输入原始地形图图像；

步骤2，对原始地形图图像进行随机概率采样，得到采样数据矩阵；

步骤3，确定地形图分割种类数；

步骤4，计算聚类中心，得到聚类中心矩阵；

步骤5，计算原始地形图的隶属矩阵，得到隶属矩阵；

步骤6，利用隶属矩阵对原始地形图进行模糊分类，得到分割后的图像；

步骤7，对分割后的图像进行多级融合，得到的融合后图像；

步骤8，输出分割后的分版图像；

所述步骤2的具体步骤如下：

2a)生成一个符合正态分布的，且和原始地形图A等大小的随机矩阵R_M×N作为随机概率矩阵,该矩阵中所有元素的值为0和1之间的任意随机数；其中，M为原始地形图的宽度，N为原始地形图的长度；

2b)计算概率阈值P，并利用步骤2a)生成的随机概率矩阵R_M×N对原始地形图进行数据采样；概率阈值P的计算如下：

$P=\frac{{10}^5}{M\×N}$

其中，M为原始地形图的宽度，N为原始地形图的长度；

数据采样过程如下：遍历随机概率矩阵R_M×N，如果矩阵中某一元素的值大于等于概率阈值P，则原始地形图A中相同位置像素点的颜色信息被选择为采样信息；否则原始地图A中相同位置像素点的颜色信息为0；遍历结束后得到采样数据矩阵A'。

2.如权利要求1所述的方法，其特征在于，所述步骤3中，利用下式确定地形图的分割种类数：

Cluster_n＝C_l ¹+C_l ²+C_l ³+…C_l ^l+1

其中，Cluster_n为地形图的分割种类数，C_l ¹为排列组合表达式，表示在由l种颜色绘制的地形图中的某区域有1种颜色相互覆盖所得到的颜色种类数，l≥1。

3.如权利要求1所述的方法，其特征在于，所述步骤4的具体步骤如下：

4a)初始化隶属矩阵：初始化一个与原始地形图A相同大小的随机矩阵作为初始隶属矩阵U₀；并使其满足下式：

$\underset{k=1}{\overset{Cluster\_n}{\Σ}}{u}_{0,k(i,j)}=1$

其中，u_0,k(i,j)为初始隶属矩阵U₀中坐标(i,j)位置第k类元素的值；

4b)利用下式计算得到矩阵大小为Cluster_n×3的聚类中心矩阵v_I，其中，I为迭代次数，当迭代次数I＝1时，隶属度矩阵U＝U₀，u_k(i,j)＝u_0,k(i,j)；

$v_{I,k}=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left(u_{k\left(i,j\right)}\right)}^m{x}_{\left(i,j\right)}}{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left(u_{k\left(i,j\right)}\right)}^m}$

其中，v_I,k为第I次迭代中聚类中心矩阵v_I的第k个聚类中心，k＝1,…,Cluster_n；m＝2；x_(i,j)为原始地形图中坐标(i,j)位置的像素值；(u_k(i,j))^m为隶属矩阵U中坐标(i,j)位置第k类元素值的m次方；

4c)根据下式计算第I次迭代的总体价值函数J_I：

$J_I=\underset{k=1}{\overset{Cluster\_n}{\mathrm{\Σ}}}{J}_{I,k}=\underset{k=1}{\overset{Cluster\_n}{\Σ}}\left(\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}{\left(u_{k(i,j)}\right)}^m\right||v_{I,k}-x_{(i,j)}|{|}^2)$

其中，J_I,k为第I次迭代的第k类价值函数；

如果总体价值函数J_I小于阀值J_min＝10^-5，或总体价值函数J_I相对上次得到的价值函数值J_I-1的改变量小于阀值J_{dif_min}＝10^-5，则步骤4的计算聚类中心结束，最终得到聚类中心矩阵v＝v_I，然后执行步骤5；否则，进行步骤4d)；

4d)使用下式计算新的隶属度矩阵U'；

${u^{\′}}_{k(i,j)}=\frac{1}{\underset{p=1}{\overset{Cluster\_n}{\Σ}}{\left(\frac{\left|\right|v_{I,k}-x_{(i,j)}\left|\right|}{\left|\right|v_{I,p}-x_{(i,j)}\left|\right|}\right)}^{2/(m-1)}}$

其中，u'_k(i,j)为新的隶属矩阵U'中坐标(i,j)位置第k类元素的值；v_I,k和v_I,p分别表示聚类中心矩阵v_I中第k个和第p个聚类中心；

令U＝U'，u_k(i,j)＝u'_k(i,j)，I＝I+1，返回步骤4b)。

4.如权利要求1所述的方法，其特征在于，所述步骤5的具体步骤如下：

5a)读取步骤(4)中得到的聚类中心矩阵v；

5b)利用模糊C均值聚类算法中的隶属度计算公式，计算原始地形图A中各像素的隶属度，得到所有像素的隶属度组成隶属矩阵U_F；

$u_{F,k(i,j)}=\frac{{(1/||v_k-x_{(i,j)}|{|}^2)}^{1/(m-1)}}{\underset{p=1}{\overset{Cluster\_n}{\Σ}}{(1/||v_p-x_{(i,j)}|{|}^2)}^{1/(m-1)}}$

其中，u_F,k(i,j)为隶属矩阵U_F中位置为(i,j)处第k类元素的值，对应于原始地形图中坐标位置的像素属于k类的隶属度；v_k和v_p分别为聚类中心矩阵v的第k个和第p个聚类中心；x_(i,j)为原始地形图中坐标(i,j)位置的像素值；m＝2。

5.如权利要求1所述的方法，其特征在于，所述步骤6的具体步骤如下：

利用隶属矩阵U_F，按以下公式所示的模糊分类方法对原始地形图A进行图像分割得到分割后的图像：

x_(i,j)∈image_k如果u_F,k(i,j)≥u_F,p(i,j) (p＝1,2,…Cluster_n)

其中，image_k为原始地形图中像素x_(i,j)被划入的第k个分割后的图像。

6.如权利要求1所述的方法，其特征在于，所述步骤7的具体步骤如下：

7a)读取聚类中心矩阵；

7b)计算每个聚类中心和其他聚类中心之间的距离；

7c)利用聚类中心之间的距离以及多级融合方法对分割后的图像进行融合，直至融合后得到的分版图像数目和地形图的分割种类数Cluster_n相等时停止，得到融合后的分版图像。

7.如权利要求6所述的方法，其特征在于，所述步骤7b)中的距离采用马氏距离：

${\mathrm{Dist}}_{k1,k2}=\sqrt{{(v_{k1}-v_{k2})}^T{\mathrm{\Σ}}^{-1}(v_{k1}-v_{k2})}$

其中，Dist_k1,k2为聚类中心矩阵v中第k₁个聚类中心v_k1和第k₂个聚类中心v_k2之间的马氏距离，()^T为矩阵的转置。

8.如权利要求6所述的方法，其特征在于，所述步骤7c)中的多级融合方法的规则：当且仅当与第i类的最小距离是第j类，同时与第j类的最小距离是第i类时，将第i类和第j类融合为一新类。