CN106408571A - 一种基于最优模糊因子选择的可变类遥感图像分割方法 - Google Patents

Abstract

本发明提供一种基于最优模糊因子选择的可变类遥感图像分割方法，包括：读取待分割的遥感图像；确定待分割遥感图像的最优同质区域类别数；通过反模糊化找到待分割遥感图像中每个像素光谱测度矢量的所属同质区域类别，得到待分割遥感图像的分割结果。本发明采用划分熵指数作为优选模糊因子的指标，当对待分割遥感图像的模糊因子小于最优模糊因子时，PE指数值较大；当模糊因子恰好等于最优模糊因子时，PE指数会跳变到一个较小的值，随着模糊因子的进一步增多，PE指数值逐渐趋于稳定；选择PE指数收敛时所对应的最小模糊因子为最优模糊因子，在最优模糊因子时对应的类别数为最优类别数，实现确定遥感图像中同质区域的类别数，获得较好的分割结果。

Description

一种基于最优模糊因子选择的可变类遥感图像分割方法

技术领域：

本发明涉及图像分割领域，尤其涉及一种基于最优模糊因子选择的可变类遥感图像分割方法。

背景技术：

图像分割是图像处理的重要步骤之一，其关键问题在于如何准确确定图像类别数及有效分割同质区域。大多数现有的图像分割算法需人为给出恰当的类别数，其中应用较为广泛的为FCM及其改进算算法。遥感图像由于覆盖范围广、地物类型复杂，人为确定图像类别数较为困难。因此，实现图像的可变类分割具有重要意义。目前，实现可变类图像分割的方法主要有统计方法和聚类方法，相对于统计方法，基于聚类原理的可变类图像分割方法众多，且应用较为广泛。其中经典的聚类方法有迭代自组织数据分析技术方法(ISODATA)，但是该方法抗噪性能较差。可见，现有的分割方法均不能实现图像的准确分割。

发明内容：

针对现有技术的缺陷，本发明提供一种基于最优模糊因子选择的可变类遥感图像分割方法。

本发明的技术方案是：

基于最优模糊因子选择的可变类遥感图像分割方法，包括以下步骤：

步骤1：读取待分割的遥感图像，测量待分割遥感图像中每个像素的光谱测度矢量，将待分割遥感图像表示为各个像素的光谱测度矢量的集合；

步骤2：确定待分割遥感图像的最优同质区域类别数；

步骤2.1：初始化参数：模糊因子初始值γ⁽⁰⁾＝1，迭代次数t＝0，模糊因子最大值为M，模糊因子的变化步长为Δγ；

步骤2.2：对待分割的遥感图像进行自组织分割；

步骤2.3：归一化第t+1次迭代后待分割遥感图像中各个像素的光谱测度矢量属于各聚类的隶属度，并计算第t+1次迭代后的划分熵指数；

步骤2.4：计算第t+1次迭代后的不同聚类中心的个数即待分割的遥感图像中同质区域类别数；

步骤2.5：迭代判定：当第t+1次迭代的模糊因子值大于模糊因子最大值M时，停止迭代，进入步骤2.6，所述第t+1次迭代的模糊因子值为第t次迭代产生的模糊因子γ^(t)与模糊因子变化步长Δγ之和；否则，令t＝t+1，，返回步骤2.2；

步骤2.6：以第t+1次迭代产生的划分熵指数为纵轴、以第t+1次迭代产生的模糊因子为横轴做划分熵指数随模糊因子变化的曲线，以第t+1次迭代产生的不同的聚类中心个数即同质区域类别数为纵轴，以第t+1次迭代产生的模糊因子为横轴做同质区域类别数随模糊因子变化的曲线；

步骤2.7：通过步骤2.6所得的曲线，选择划分熵指数达到稳态时对应的最小的模糊因子γ作为最优模糊因子γ^*，最优模糊因子γ^*对应的不同的聚类中心个数作为最优同质区域类别数c*，完成待分割图像的类别数确定；

步骤3：通过反模糊化c*个聚类中心找到待分割遥感图像中每个像素光谱测度矢量的所属同质区域类别，得到待分割遥感图像的分割结果。

所述步骤2.2包括：

步骤2.2.1：设置内循环迭代指示器l，并初始化为0；

步骤2.2.2：将待分割图像的所有光谱测度矢量作为初始聚类中心；

步骤2.2.3：计算第l+1次迭代后第i个像素点的光谱测度矢量属于第j个聚类的隶属度||x_i-z_j ^(l)||为待分割遥感图像中第i个像素的光谱测度矢量x_i与第l次迭代后第j个聚类的聚类中心之间的欧氏距离，待分割遥感图像中各像素的光谱测度矢量的方差其中为待分割遥感图像中各像素的光谱测度矢量的均值；

步骤2.2.4：计算第l+1次迭代后的聚类中心

步骤2.2.5：重复步骤2.2.1～2.2.4，直到当max||z_j ^(l+1)-z_j ^(l)||＜ε时停止，ε表示判定阈值，否则更新迭代次数l，l＝l+1，并进入步骤2.2.3。

由上述技术方案可知，本发明的提出的基于最优模糊因子选择的可变类遥感图像分割方法，采用划分熵(Partition Entropy，PE)指数作为优选模糊因子的指标，当对待分割遥感图像的模糊因子小于最优模糊因子时，PE指数值较大；当对待分割遥感图像的模糊因子恰好等于最优模糊因子时，PE指数会跳变到一个较小的值；随着模糊因子的进一步增多，PE指数值逐渐趋于稳定。因此，选择PE指数收敛时所对应的最小模糊因子为最优模糊因子，在最优模糊因子时对应的类别数为最优类别数。以最优化模糊因子来确定最优类别数，实现了自动确定遥感图像中同质区域的类别数，可获得较好的分割结果。通过实验本发明和ISODATA方法的合成图像分割结果进行定量精度评价本发明各个精度指标均大于ISODATA方法，且最小精度指标值为98.1％。

附图说明：

图1为本发明具体实施方式中基于最优模糊因子选择的可变类遥感图像分割方法流程图；

图2为本发明具体实施方式中步骤2具体流程图；

图3为本发明具体实施方式中步骤2.2具体流程图；

图4为本发明具体实施方式中以数据点为例的同质区域类别数和模糊因子的关系，其中(a)为实验中随机生成的300个数据点；(b)为模糊因子γ＝1的实验得到同质区域类别数c＝1的结果；(c)为模糊因子γ＝5的实验得到同质区域类别数c＝3的结果；(d)为模糊因子γ＝50的实验得到同质区域类别数c＝28的结果；

图5为本发明具体实施方式中对图4(a)中的数据点实验得到划分熵指数和同质区域类别数分别随模糊因子变化的曲线，其中，(a)为PE随γ变化的曲线，(b)为c随γ变化的曲线；

图6为本发明具体实施方式中利用最优模糊因子下最优同质区域类别数进行图像分割的实验结果；

图7为本发明具体实施例中的合成图像，其中，(a)为合成图像模板，(b)为合成图像；

图8为本发明具体实施例中对图7(b)实验得到的划分熵指数和同质区域类别数分别随模糊因子变化的曲线；

图9为本发明具体实施方式中利用本发明方法和ISODATA方法对图7(b)的分割结果；

图10(a)～(b)分别为本发明具体实施例中应用的两幅真实遥感图像；

图11为本发明具体实施方式中对真实遥感图像实验得到的划分熵指数和同质区域类别数分别随模糊因子变化的曲线，其中，(a)与(c)分别为本发明方法对图10(a)实验得到的划分熵指数和同质区域类别数分别随模糊因子变化的曲线，(b)与(d)分别为本发明方法对图10(b)实验得到的划分熵指数和同质区域类别数分别随模糊因子变化的曲线；

图12为本发明具体实施方式中应用本发明方法和ISODATA方法对真是遥感图像的分割结果，其中，(a)与(b)分别为本发明方法对图10(a)与图10(b)的分割结果，(c)与(d)分别为ISODATA方法对图10(a)与图10(b)的分割结果。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

基于最优模糊因子选择的可变类遥感图像分割方法，如图1所示，包括以下步骤：

步骤1：读取待分割的遥感图像，测量待分割遥感图像中每个像素的光谱测度矢量，将待分割遥感图像表示为各个像素的光谱测度矢量的集合；

本实施方式中，定义待分割的遥感图像X＝{x_i，i＝1，...，n}，其中，i为像素索引，n为像素数，x_i＝(x_i1，...，x_id)为像素i的光谱测度矢量，d为波段数。

步骤2：确定待分割遥感图像的最优同质区域类别数；

如图2所示，步骤2包括：

步骤2.1：初始化参数：模糊因子初始值γ⁽⁰⁾＝1，迭代次数t＝0，模糊因子最大值为M，模糊因子的变化步长为Δγ；

步骤2.2：对待分割的遥感图像进行自组织分割；

如图3所示，所述步骤2.2包括：

步骤2.2.1：设置内循环迭代指示器l，并初始化为0；

步骤2.2.2：将待分割遥感图像的所有光谱测度矢量作为初始聚类中心；本实施方式中，z＝{z_j，j＝1，...，c}为聚类中心集合，其中，j为聚类索引，c为同质区域类别数，z_j＝(z_j1，...，z_jd)为第j个聚类中心矢量。

步骤2.2.3：计算第l+1次迭代后第i个像素点的光谱测度矢量属于第j个聚类的隶属度||x_i-z_j ^(l)||为待分割遥感图像中第i个像素的光谱测度矢量x_i与第l次迭代后第j个聚类的聚类中心之间的欧氏距离，待分割遥感图像中各像素的光谱测度矢量的方差其中为待分割遥感图像中各像素的光谱测度矢量的均值；

步骤2.2.4：计算第l+1次迭代后的聚类中心

本实施方式中，U＝[μ_ij]_n×c表达X的模糊分割，为了得到最优模糊分割，定义如下目标函数，

其中，γ(＞0)是模糊因子，表征模糊程度。J(z)描述所有像素的光谱测度矢量对所有同质区域类别的隶属度指数之和。为实现最优分割，需最大化目标函数。

目标函数对z_j求偏导，

令其为零时求得z_j，

即

步骤2.2.5：重复步骤2.2.1～2.2.4，直到当max||z_j ^(l+1)-z_j ^(l)||＜ε时停止，ε表示判定阈值，否则更新迭代次数l，l＝l+1，并进入步骤2.2.3。

步骤2.3：归一化第t+1次迭代后待分割遥感图像中各个像素的光谱测度矢量属于各聚类的隶属度，并计算第t+1次迭代后的划分熵指数PE^(l+1)；

本实施方式中，首先将μ_ij归一化为μ_ij′，并使μ_ij′满足

PE指数定义为：

步骤2.4：计算第t+1次迭代后的不同聚类中心的个数即待分割的遥感图像中同质区域类别数c^(l+1)；

步骤2.5：迭代判定：当第t+1次迭代的模糊因子值大于模糊因子最大值M时，停止迭代，得到PE⁽¹⁾...PE^(t+1)和c⁽¹⁾...c^(t+1)，进入步骤2.6，所述第t+1次迭代的模糊因子值γ^(t+1)为第t次迭代产生的模糊因子γ^(t)与模糊因子变化步长Δγ之和；否则，令t＝t+1，，返回步骤2.2；

步骤2.6：以第t+1次迭代产生的划分熵指数PE⁽¹⁾……PE^(l+1)为纵轴、以第t+1次迭代产生的模糊因子γ⁽¹⁾……γ^(t+1)为横轴做划分熵指数随模糊因子变化的曲线，以第t+1次迭代产生的不同的聚类中心个数即同质区域类别数c⁽¹⁾……c^(l+1)为纵轴，以第t+1次迭代产生的模糊因子γ⁽¹⁾……γ^(t+1)为横轴做同质区域类别数随模糊因子变化的曲线；

步骤2.7：通过步骤2.6所得的曲线，选择划分熵指数达到稳态时对应的最小的模糊因子γ作为最优模糊因子γ^*，最优模糊因子γ^*对应的不同的聚类中心个数作为最优同质区域类别数c*，完成待分割图像的类别数确定；

在本实施方式中，PE的稳定性可由连续记录δ次γ对应的PE值{PE^(γ)，...，PE^(γ+δ)}判断，即，如果(max{PE^(γ)，...，PE^(γ+δ)}-min{PE^(γ)，...，PE^(γ+δ)})/max{PE^(γ)，...，PE^(γ+δ)}＜ε，则认为PE收敛，其中，δ和ε为预先指定常数。进一步，则有γ*＝arg min{PE^(γ)，...，PE^(γ+δ)}。

步骤3：通过反模糊化c*个聚类中心找到待分割遥感图像中每个像素光谱测度矢量的所属同质区域类别，得到待分割遥感图像的分割结果。

为验证本发明的有效性，本发明在CPU为Core(TM)i5-3470 3.20GHz、内存4GB、Windows 7旗舰版系统上使用MATLAB R2014b软件编程实现仿真，并通过以下实验检验本发明的分割精度。

图4为本发明实施方式中以数据点为例说明同质区域类别数和模糊因子的关系，其中(a)为实验中随机生成的300个数据点；(b)为设置模糊因子γ＝1的实验，得到同质区域类别数c＝1的结果；(c)为设置模糊因子γ＝5的实验，得到同质区域类别数c＝3的结果；(d)为设置模糊因子γ＝50的实验，得到同质区域类别数c＝28的结果，图4(a)、(b)、(c)、(d)中横坐标x与纵坐标y代表生成数据点的横轴坐标。由图4可知，当模糊因子γ越小时得到的同质区域类别数c越少，当模糊因子γ越大时得到的同质区域类别数c越多。

图5为本发明具体实施方式中对图4(a)中的数据点实验得到划分熵指数和同质区域类别数分别随模糊因子变化的曲线；其中，(a)为划分熵指数随模糊因子变化的曲线，(b)为同质区域类别数随模糊因子变化的曲线。由图5可知，本发明得到的最优模糊因子γ*＝5，γ*＝5对应的最优类别数c*＝3，符合实际同质区域类别数。本发明实现了数据点的自动确定类别数。

图6为本发明具体实施方式中利用最优模糊因子γ*＝5下最优类别数c*＝3进行图像分割的实验结果，由图6可知，本发明的分割结果较好。

图7为本发明具体实施方式中的合成图像，其中，(a)为合成图像模板，(b)为合成图像；其中，(a)中I～V代表不同的同质区域；(b)为以(a)为模板分别截取草地、林地、绿地、海水和裸地5类合成的图像。

图8为本发明具体实施方式中对图7(b)实验得到的划分熵指数和同质区域类别数分别随模糊因子变化的曲线；其中(a)和(b)分别为对图7(b)实验得到的划分熵指数和同质区域类别数分别随模糊因子变化的曲线，通过曲线可以看出，合成图像最优模糊因子γ*＝13，相应的最优同质区域类别数c*＝5，符合实际同质区域类别数。

图9为本发明具体实施方式中利用本发明方法和ISODATA方法对图7(b)的分割结果；其中，(a)和(b)分别为利用本发明方法和ISODATA方法对图7(b)实验的分割结果；由图可知，本发明方法和ISODATA方法均能正确确定同质区域类别数，但本发明分割结果较好，而ISODATA方法分割结果有很多被误分区域，且噪声较多。

为了验证本发明的有效性，对本发明和ISODATA方法的合成图像分割结果进行定量精度评价，以图7(a)合成图像模板中各同质区域为标准计算混淆矩阵，并根据混淆矩阵进一步计算产品精度、用户精度、总精度和kappa值。精度计算结果如表1所示，可以看出，本发明各个精度指标均大于ISODATA方法，且最小精度指标值为98.1％。

表1：ISODATA聚类方法和本发明方法分割结果进行定量评价

图10(a)～(b)分别为本发明具体实施例中对应用的两幅遥真实感图像；

图11为本发明具体实施例中对真实遥感图像实验得到的划分熵指数和类别数分别随模糊因子变化的曲线，其中，(a)与(c)分别为本发明方法对图10(a)实验得到的划分熵指数和类别数分别随模糊因子变化的曲线，(b)与(d)分别为本发明方法对图10(b)实验得到的划分熵指数和类别数分别随模糊因子变化的曲线；由图可以看出，两幅真实遥感图像最优模糊因子分别为γ*＝8、13，相应的最优类别数c*＝3、5。本发明对两幅真实遥感图像均实现了自动确定类别数。

图12为本发明具体实施例中应用本发明方法和ISODATA方法对真实遥感图像的分割结果，其中，(a)与(b)分别为本发明方法对图10(a)与图10(b)的分割结果，(c)与(d)分别为ISODATA方法对图10(a)与图10(b)的分割结果。由图可以看出，本发明在正确确定类别数的基础上，对各同质区域分割效果均较好，而IOSDATA方法不能克服噪声影响，分割结果较差。

本发明提供的基于最优模糊因子选择的可变类遥感图像分割方法，采用划分熵(Partition Entropy，PE)指数作为优选模糊因子的指标，当对待分割遥感图像的模糊因子小于最优模糊因子时，PE指数值较大；当对待分割遥感图像的模糊因子恰好等于最优模糊因子时，PE指数会跳变到一个较小的值；随着模糊因子的进一步增多，PE指数值逐渐趋于稳定。因此，选择PE指数收敛时所对应的最小模糊因子为最优模糊因子，在最优模糊因子时对应的类别数为最优类别数。以最优化模糊因子来确定最优类别数，实现了自动确定遥感图像的类别数，可获得较好的分割结果。通过与ISODATA方法进行对比进一步说明了本发明的可行性和有效性。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (2)

1.基于最优模糊因子选择的可变类遥感图像分割方法，其特征在于，包括以下步骤：

步骤1：读取待分割的遥感图像，测量待分割遥感图像中每个像素的光谱测度矢量，将待分割遥感图像表示为各个像素的光谱测度矢量的集合；

步骤2：确定待分割遥感图像的最优同质区域类别数；

步骤2.1：初始化参数：模糊因子初始值γ⁽⁰⁾＝1，迭代次数t＝0，模糊因子最大值为M，模糊因子的变化步长为Δγ；

步骤2.2：对待分割的遥感图像进行自组织分割；

步骤2.3：归一化第t+1次迭代后待分割遥感图像中各个像素的光谱测度矢量属于各聚类的隶属度，并计算第t+1次迭代后的划分熵指数；

步骤2.4：计算第t+1次迭代后的不同聚类中心的个数即待分割的遥感图像中同质区域类别数；

步骤2.5：迭代判定：当第t+1次迭代的模糊因子值大于模糊因子最大值M时，停止迭代，进入步骤2.6，所述第t+1次迭代的模糊因子值为第t次迭代产生的模糊因子γ^(t)与模糊因子变化步长Δγ之和；否则，令t＝t+1，，返回步骤2.2；

步骤2.6：以第t+1次迭代产生的划分熵指数为纵轴、以第t+1次迭代产生的模糊因子为横轴做划分熵指数随模糊因子变化的曲线，以第t+1次迭代产生的不同的聚类中心个数即同质区域类别数为纵轴，以第t+1次迭代产生的模糊因子为横轴做同质区域类别数随模糊因子变化的曲线；

步骤2.7：通过步骤2.6所得的曲线，选择划分熵指数达到稳态时对应的最小的模糊因子γ作为最优模糊因子γ^*，最优模糊因子γ^*对应的不同的聚类中心个数作为最优同质区域类别数c*，完成待分割图像的类别数确定；

步骤3：通过反模糊化c*个聚类中心找到待分割遥感图像中每个像素光谱测度矢量的所属同质区域类别，得到待分割遥感图像的分割结果。

2.根据权利要求1所述的方法，其特征在于，所述步骤2.2包括：

步骤2.2.1：设置内循环迭代指示器l，并初始化为0；

步骤2.2.2：将待分割遥感图像的所有光谱测度矢量作为初始聚类中心；

步骤2.2.3：计算第l+1次迭代后第i个像素点的光谱测度矢量属于第j个聚类的隶属度||x_i-z_j ^(l)||为待分割遥感图像中第i个像素的光谱测度矢量x_i与第l次迭代后第j个聚类的聚类中心之间的欧氏距离，待分割遥感图像中各像素的光谱测度矢量的方差其中为待分割遥感图像中各像素的光谱测度矢量的均值；

步骤2.2.4：计算第l+1次迭代后的聚类中心

步骤2.2.5：重复步骤2.2.1～2.2.4，直到当对停止，ε表示判定阈值，否则更新迭代次数l，l＝l+1，并进入步骤2.2.3。