子空间与KL信息结合的FCM多光谱遥感图像分割

Fuzzy C-mean Multi-spectral Remote Sensing Image Segmentation with Combined Subspace and KL Information

摘要

针对传统模糊C-均值聚类（FCM）算法用于含噪声多光谱遥感图像分割时存在的精度不足问题，提出一种自适应模糊子空间与增强KL信息相结合的FCM多光谱遥感图像分割算法.首先，使用局部模糊因子，在不依赖参数的情况下，通过相似性度量和自适应约束参数自动消除噪声干扰，并提取图像的局部空间信息.其次，将原始图像信息和模糊因子处理过的局部空间信息统一整合到模糊子空间聚类中，对图像的多个通道进行自适应加权处理，以提高分割精度.最后，将KL信息以正则项的形式引入FCM目标函数中进行聚类计算，并通过ESD（Extreme Studentized Deviate）检测模型剔除隶属度矩阵中的离群值，以增强KL先验信息，降低隶属度模糊性.AID数据库和真实环境下的多光谱遥感图像分割实验表明，在模拟噪声环境中，所提出算法不仅可以抑制噪声，而且能得到较高的分割精度.此外，本文算法在分割精度、模糊系数和峰值信噪比等评价指标方面也均优于其他几种变体式FCM算法.

Abstract

For the problem of insufficient accuracy of traditional fuzzy C-means clustering （FCM） algorithm for noise-containing multi-spectral remote sensing image segmentation， an FCM multi-spectral remote sensing image segmentation algorithm combining adaptive local fuzzy subspace and enhanced KL is proposed. Firstly， the local fuzzy factor is used to automatically eliminate the noise interference and extract the local spatial information of the image by similarity metric and adaptive constraint parameters without relying on parameters. Secondly， the original image information and the local spatial information processed by the fuzzy factor are unified and integrated into the fuzzy subspace clustering， and the multiple channels of the image are adaptively weighted to enhance the segmentation accuracy. Finally， the KL information is introduced into the FCM objective function in the form of regular terms for clustering calculation， and the outliers in the membership matrix are removed by ESD （Extreme Studentized Deviate） detection model to enhance the KL prior information and reduce the ambiguity of the membership. The experiments of real multi-spectral remote sensing image segmentation show that in the simulation of noise environments， the algorithm in this paper can suppress the noise and can guarantee the segmentation accuracy better at the same time. In addition， the algorithm in this paper outperforms several other variant FCM algorithms in terms of evaluation indexes such as segmentation accuracy， fuzzy coefficient， and peak signal-to-noise ratio.

Graphical abstract

关键词

模糊聚类 / 子空间 / 模糊因子 / KL信息 / 图像分割

Key words

fuzzy clustering / subspace / fuzzy factor / KL information / image segmentation

随着人工智能和卫星传感器的发展，遥感技术突破了众多挑战，进入一个能够准确、高效地提供多种地物信息的新阶段^［1］.但在很多虚拟场景和实践当中^［2-4］，因采集工作和传输过程的影响，接收到的多光谱遥感图像通常会受到不同强度噪声干扰，从而影响图像分割精度.在过去几十年中，已产生了多种含噪声的多光谱遥感图像分割算法.在这些算法中，模糊C-均值聚类（FCM）^［5］算法结构简单，不需要训练样本，仍具有很大的研究价值.但如何进一步提升FCM算法的分割精度和噪声的抑制能力，目前还有许多可研究的空间.

为抑制FCM算法对噪声的敏感性，同时提高其在多光谱遥感图像应用中的分割精度和实用普适性，近几年来，学者们进行了大量研究和探索，在利用FCM对多光谱遥感图像分割方面取得了突破性进展.文献［6］提出了一种基于局部空间信息的算法，该算法利用各像素邻域集合的统计信息与聚类中心的相似度修正原始图像中像素的隶属度，在去除噪声的同时进行聚类分割.但随着噪声密度的增加，多光谱遥感图像像素的局部邻域信息可能被过多噪声掩盖，因此仅考虑中心像素的邻域信息是不充分的.文献［7］提出一种基于图像空间信息自适应变化窗口参数的算法，该算法针对搜索窗口和邻域窗口固定的问题进行合理的改进，并利用多光谱遥感图像的空间信息自适应地计算出搜索窗口和邻域窗口的取值，对于不同的图像，算法可以灵活计算其尺寸，在一定程度上提高了算法的适用场景，但该算法计算复杂度较高.为进一步挖掘空间信息在多光谱遥感图像分割中的作用，文献［8］通过寻找图像中相似邻域配置的特征来获取非局部空间信息，这样非局部空间信息中就包含了与当前像素邻域具有一定相似性的其他像素邻域，进而可通过该信息对噪声进行抑制，但是出现噪声破坏相似区域的小概率事件时，仍然会导致分割效果较差的问题.文献［9］提出一种结合形状信息的静态MST区域划分和隐马尔可夫FCM算法的高分辨率遥感图像分割算法，解决了算法抗噪声性差以及对地物复杂边界分割精度低的问题，但当其应用于较多场景的遥感图像时，该算法仍有不足之处.文献［10］提出一种基于多特征融合的FCM算法，该类算法首先提取图像的颜色通道或其他简单特征，再将其相加或直接拼接形成一个大的特征向量组，最后将该特征向量组输入聚类模型参与后续分割.但多特征融合方式主要采用平均加权，特征对聚类的主次作用均化，导致自适应能力不足，出现算法对部分图像欠分割的现象.

虽然上述算法对噪声具有一定的鲁棒性，但仍存在以下问题：1）在构造目标函数时，原始图像和修正图像（空间修正运算消除噪声）之间的约束参数为人为设定，就遥感图像本身而言，其地物种类和场景众多，对每一类场景设置不同的约束参数，会削弱算法的普适性和对噪声的抑制能力；2）多数FCM变体算法，都将通道特征对聚类结果的贡献做平均化处理，然而，在多光谱遥感图像分割中，特征对聚类结果的贡献程度可能是不同的，平均思想可能会降低算法的聚类效率，影响分割精度；3）在高密度噪声环境下，隶属度模糊化现象也会随之加重，带有隶属度均值先验信息的KL正则项不能更好地消除隶属度模糊现象.

针对上述问题，提出一种自适应模糊子空间与增强KL信息相结合的FCM多光谱遥感图像分割算法.首先，使用局部模糊因子，通过参数自适应的方法，进行相似性度量和约束参数自适应计算自动消除噪声干扰，并提取图像的局部空间信息.其次，将经过模糊聚类计算的原始图像信息和模糊因子处理过的局部空间信息统一整合到模糊子空间中，对图像的多个通道进行自适应加权处理，以提高分割精度.最后，将ESD模型增强的KL信息作为正则项融入目标函数中，剔除隶属度离群值，减少模糊现象.

1 相关工作

1.1 模糊聚类算法

FCM算法^［5］是图像分割领域中一种常用的算法，相比其他聚类算法，它可以保留更多的图像信息.该算法通过每个像素点到聚类中心的距离来考虑隶属度值，并通过目标函数的收敛将像素分配到合适的集群中，其目标函数为：

J^{(a)}={\displaystyle \sum _{i=\mathrm{1}}^K}{\displaystyle \sum _{j=\mathrm{1}}^N}(u_{ij}^{(a)}{)}^m{‖y_j-c_i^{(a)}‖}^{\mathrm{2}},s.t.{\displaystyle \sum _{i=\mathrm{1}}^K}{u}_{ij}^{(a)}=\mathrm{1}（1）

式中：K表示分割区域（聚类）的数量，y_j是图像I的第j个像素的灰度值，c_i^{(a)}是第a次迭代中的第i个聚类中心，u_{ij}^{(a)}是第a次迭代中第j个像素对第i个聚类簇的隶属度，m是隶属度加权指数，通常情况下为2.应用拉格朗日乘子法，可推导出隶属度u_{ij}^{(a)}和聚类中心c_i^{(a)}如下：

u_{ij}^{(a)}={{\displaystyle \sum _{r=\mathrm{1}}^K}\left(\frac{‖y_j-c_i^{(a)}‖}{‖y_j-c_r^{(a)}‖}\right)}^{-\frac{\mathrm{2}}{m-\mathrm{1}}}（2）

c_i^{(a)}=\frac{{\displaystyle \sum _{j=\mathrm{1}}^K}{\left(u_{ij}^{(a)}\right)}^m{y}_j}{{\displaystyle \sum _{j=\mathrm{1}}^K}{\left(u_{ij}^{(a)}\right)}^m}（3）

1.2 子空间聚类算法

传统的FCM算法只能对彩色图像的所有通道进行平均处理，但在实际聚类过程中，图像的不同通道信息对聚类中心的贡献是不相同的.为进一步优化FCM算法对彩色图像的分割性能，文献［11］提出了一种基于模糊子空间的聚类（Fuzzy Subspace Clustering，FSC）算法，该算法通过对每个通道的像素分配适当的权重，提升FCM算法的彩色图像分割能力.FSC算法的目标函数为：

\begin{array}{l}{J}_{\mathrm{F}\mathrm{S}\mathrm{C}}={\displaystyle \sum _{i=\mathrm{1}}^k}{\displaystyle \sum _{j=\mathrm{1}}^N}{u}_{ij}{\displaystyle \sum _{z=\mathrm{1}}^D}{w}_{iz}^{\tau }{\left(y_{jz}-c_{iz}\right)}^{\mathrm{2}}+\gamma {\displaystyle \sum _{i=\mathrm{1}}^k}{\displaystyle \sum _{z=\mathrm{1}}^D}{w}_{iz}^{\tau },\\ {\displaystyle \sum _{i=\mathrm{1}}^k}{u}_{ij}=\mathrm{1},{\displaystyle \sum _{z=\mathrm{1}}^D}{w}_{iz}=\mathrm{1}\end{array}（4）

式中：D为数据维度（图像通道），w_{iz}^{\tau }为加权因子，\gamma为子空间约束参数，\tau为模糊加权指数，通常\tau >\mathrm{1}， 式（4）中i=\mathrm{1,2},\cdots ,k；j=\mathrm{1,2},\cdots ,N.通过拉格朗日乘子法，可以得到u_{ij}、w_{iz}和c_{iz}的迭代公式如下：

u_{ij}=\left\{\begin{array}{l}\mathrm{1},{\displaystyle \sum _{z=\mathrm{1}}^D}{w}_{iz}^{\tau }{\left(y_{jz}-c_{iz}\right)}^{\mathrm{2}}\le {\displaystyle \sum _{z=\mathrm{1}}^D}{w}_{iz}^{\tau }{\left(y_{jz}-c_{pz}\right)}^{\mathrm{2}}\\ \mathrm{0}, \mathrm{其他}\end{array}\right.（5）

w_{iz}=\frac{{\left({\displaystyle \sum _{j=\mathrm{1}}^N}{u}_{ij}{\left(y_{jz}-c_{iz}\right)}^{\mathrm{2}}+\gamma \right)}^{\frac{-\mathrm{1}}{\tau -\mathrm{1}}}}{{\displaystyle \sum _{z=\mathrm{1}}^D}{\left({\displaystyle \sum _{j=\mathrm{1}}^N}{u}_{ij}{\left(y_{jz}-c_{iz}\right)}^{\mathrm{2}}+\gamma \right)}^{\frac{-\mathrm{1}}{\tau -\mathrm{1}}}}（6）

c_{iz}={\displaystyle \sum _{j=\mathrm{1}}^N}{u}_{ij}{y}_j/{\displaystyle \sum _{j=\mathrm{1}}^N}{u}_{ij}（7）

2 本文算法

2.1 整合局部模糊因子的子空间

为提升FCM算法对噪声的鲁棒性，文献［12］提出一种基于局部模糊因子的聚类算法.该算法通过局部运算修正FCM的隶属度矩阵，用于度量像素之间的局部相似性，在去除噪声的同时保持图像的细节信息.局部模糊因子G_{ir}为：

G_{ir}=\sum _{r\in N_j^s}\frac{\mathrm{1}}{d_{rj}+\mathrm{1}}{\left(\mathrm{1}-u_{ir}^{(a)}\right)}^m{‖y_r-c_i^{(a)}‖}^{\mathrm{2}}（8）

式中：N_j^s是以像素y_r为中心半径为s\times s邻域，y_r为该领域中的其他像素，s一般取值为3，中心像素与邻域像素之间的欧氏距离为d_{rj}.局部模糊因子通过定义每个像素的模糊性自动控制像素灰度级和空间结构性，达到抑制噪声的目的.本文将G_{ir}用于多光谱遥感图像分割中，通过其相性度量特性自动消除噪声干扰，经过G_{ir}处理后生成一幅新的图像，称为图像的修正信息.多数改进FCM算法在构建目标函数时，原始图像与修正图像之间的约束参数多为人为设定，如式（12）所示.对于不同的遥感图像，因其地物目标复杂多样，像素之间差异性较大，应该对不同的场景设置不同的约束参数，而这些参数通常是根据经验预先设定，需要多次重复实验修改和估计，人为设定不当，会影响对噪声的抑制能力和分割精度.因此，本文通过定义平滑度（局部方差）代替单个像素，实现参数\alpha和\beta的自适应计算，而平滑度又可较好地描述局部邻域像素的偏离程度，平滑度越小，表明图像被噪声污染越轻，反之亦然.将\alpha和\beta设为互为倒数的关系，便可以进一步增强自适应参数与图像自身信息之间的相关性.本文使用原始图像信息与局部空间信息通过式（10）进行了平滑度计算，这样做是因为单个像素难以反映图像局部的噪声污染情况，如果使用平滑度来代替单个像素可较好地描述局部邻域像素混乱程度，即噪声污染程度.其中，W_x^l是以像素y_j为中心的l\times l邻域，N_h为该邻域的基数.通过多次实验表明，l和h取3能得到较好的分割结果.

而对于约束参数\alpha和\beta的讨论可以根据图像受到噪声或干扰的强弱分为以下两种情况： 1）当图像被噪声等干扰因素破坏严重时，图像像素值和图像的空间结构改变较大，需要通过空间运算修正这种改变，而此时原始图像和修正图像之间的差值的绝对值也相应较大.由于被噪声等干扰破坏的原始图像无法有效表征图像的原始信息，需要更多地利用修正图像信息，因此在构造目标函数时，修正图像的约束参数值应较大，而原始图像的约束参数值应较小.本文先将通过平滑运算的原始图像信息\eta \left(y_j\right)和修正图像信息\eta \left(G\right)（通过局部模糊因子计算的图像局部空间信息）作差，然后对它们的差值进行了指数函数映射，用差值和其相应的倒数自适应约束修正图像和原始图像.为防止上述两种信息在计算中数值过分接近，使得\mathrm{\Delta }\eta出现极值的情况，在式（11）中加了一个极小值\mathrm{e}\mathrm{p}\mathrm{s}（可取两位小数点后的极小数值），式（11）中，k为一个常数，设置为40，这样便可实现约束参数的自适应计算. 2）当图像被噪声等各种干扰因素破坏较轻时，图像像素值和图像的空间结构改变较小，此时原始图像和修正图像之间的差值的绝对值也较小，由于原始图像破坏程度较轻，在构造目标函数时，修正图像的约束参数应较小，而原始图像的约束参数值应较大.构造的约束关系的目标函数，以及平滑度和自适应约束参数定义如下：

\alpha =\mathrm{\Delta }\eta ,\beta =\frac{\mathrm{1}}{\mathrm{\Delta }\eta }（9）

\eta =\frac{\mathrm{1}}{N_h}\sum _{i\in W_x^l}{‖y_j-y_i‖}^{\mathrm{2}},x=\mathrm{1,2},\dots ,N（10）

\mathrm{\Delta }\eta =k\mathrm{e}\mathrm{x}\mathrm{p}\left[‖\eta \left(y_j\right)-\eta \left(G\right)‖+\mathrm{e}\mathrm{p}\mathrm{s}\right]（11）

\begin{array}{l}{J}_{\alpha \beta }^{(a)}={\displaystyle \sum _{i=\mathrm{1}}^K}{\displaystyle \sum _{j=\mathrm{1}}^N}\alpha {\left(u_{ij}^{(a)}\right)}^m{‖y_j-c_i^{(a)}‖}^{\mathrm{2}}+\\ \beta \sum _{r\in N_j}\frac{\mathrm{1}}{d_{rj}+\mathrm{1}}{\left(\mathrm{1}-u_{ir}^{(a)}\right)}^m{‖y_r-c_i^{(a)}‖}^{\mathrm{2}}\end{array}（12）

为说明自适应局部模糊因子的噪声抑制能力，使用Gibbs采样器产生两幅不同的合成图像进行对比实验，如图1和图2所示.其中，（a）是对合成图像使用5%的混合噪声污染后的图像，（b）为（a）中用红色框标记的像素值，从（a）中可以明显看出，被红色框标记的像素是一类像素，（c）为原始图像.（b）、（d）和（e）中灰色背景标记的像素为被噪声污染的像素，（d）和（e）分别对应本文算法和对比算法^［12］处理后的3个聚类簇隶属度值，分两种情况讨论两种算法的性能，情况1对应图1，情况2对应图2.

情况1：所标注图像像素区域的中心像素没有噪声，但它周围的邻域像素部分被噪声污染.图中的噪声像素的灰度值分别为133、194、101，这与中心像素值（226）相差较大.通过本文算法和文献［12］算法处理后，理论上迭代后属于某一类的隶属度值应大于0.5，但文献［12］算法在迭代后，对被噪声污染的像素没有正确分类，部分像素的隶属度值都小于0.5，如图1（e）中红色标记框所示，说明这种算法难以分割被噪声污染严重的图像.而本文算法在迭代完成之后，隶属度值均大于0.91，如图1（d）所示，说明本文算法可以较好地分割受到噪声污染的图像.

情况2：所标注图像像素区域的中心像素被噪声污染，而它的邻域没有噪声.这种情况中心像素值与邻域像素差别较大.通过本文算法和文献［12］处理后，文献［12］算法在迭代后，部分像素被误分割，中心噪声像素也未能被正确分割，如图2（e）中红色标记框所示，只有部分的像素被正确分割.而本文算法的隶属度值多数都大于0.94，如图2（d）所示，说明在此情况下，本文算法仍然具有同样的噪声抑制能力.

在使用局部模糊因子和自适应约束参数消除噪声的同时，本文将模糊聚类计算的原始图像信息和模糊因子处理过的空间信息统一整合到模糊子空间中，并对图像的多个通道进行自适应加权处理，数学表达式如下：

J_{S\alpha \beta }^{(a)}={\displaystyle \sum _{i=\mathrm{1}}^K}{\displaystyle \sum _{j=\mathrm{1}}^N}\left[{\displaystyle \sum _{z=\mathrm{1}}^D}{\left(w_{iz}^{(a)}\right)}^{\tau }\left(\alpha u_{ij}^{(a)}{‖y_{jz}-c_{iz}^{(a)}‖}^{\mathrm{2}}+\beta G_{ir}\right)\right]+

\gamma {\displaystyle \sum _{i=\mathrm{1}}^k}{\displaystyle \sum _{z=\mathrm{1}}^D}{\left(w_{iz}^{(a)}\right)}^{\tau }（13）

在传统的FCM算法中，一般是将每个通道的特征对聚类的贡献度做平均化处理.但在实际图像分割场景中，图像的每个通道对聚类过程的贡献程度可能是不相同的，存在某个通道比其他通道特征更加重要.如图3所示，其中，图3（a）为彩色图像，图3（b）~（d）为图3（a）的三通道灰度图，从中可以看出，图3（b）的目标与背景可以清楚区分，相比之下，图3（c）的目标与背景只有部分重叠，但与前两个通道比较，图3（d）的目标与背景大部分被混叠在一起难以区分.上述举例充分证明在实际分割中，应该对不同的通道特征赋予不同的权重，这样可以进一步提升图像的分割效率.

2.2 基于ESD模型增强的KL信息

统计学中KL信息通常被用于度量变量在近似估计中的信息损失.而FCM算法在分割图像时，因其随机初始化隶属度，导致隶属度模糊，进而影响后续的图像分割结果.为降低隶属度的模糊性，增加隶属度之间的差异，文献［13］将KL信息作为一种约束关系引入目标函数中，其数学表达式如下：

\begin{array}{l}{J}_{\mathrm{K}\mathrm{L}\mathrm{F}\mathrm{C}\mathrm{M}}={\displaystyle \sum _{i=\mathrm{1}}^K}{\displaystyle \sum _{j=\mathrm{1}}^N}{u}_{ij}^{(a)}\left({‖y_j-c_i^{(a)}‖}^{\mathrm{2}}+\phi \mathrm{l}\mathrm{n}\frac{u_{ij}^{(a)}}{{\pi }_i^{(a)}}\right)\\ \mathrm{s}.\mathrm{t}.{\displaystyle \sum _{i=\mathrm{1}}^K}{u}_{ij}^{(a)}=\mathrm{1},{\displaystyle \sum _{i=\mathrm{1}}^K}{\pi }_i^{(a)}=\mathrm{1}\end{array}（14）

{\pi }_i^{(a)}=\frac{\mathrm{1}}{N}{\displaystyle \sum _{j=\mathrm{1}}^N}{u}_{ij}^{(a-\mathrm{1})}（15）

式中：\phi为KL信息的约束因子，{\pi }_i^{(a)}表示第i个聚类簇的隶属度先验信息的平均值，计算公式为式（15），以上约束因子与对数函数的映射关系可以作为一种降低隶属度的模糊性的修正项.从式（14）中可以看出，修正项的主要因子是{\pi }_i，如果{\pi }_i与u_{ij}近似相等，就可以很好地估计u_{ij}，以降低隶属度模糊性.文 献［14］和［15］均通过计算隶属度矩阵的均值来获取先验信息，这些方法虽然简单高效，但当图像被噪声污染严重时，隶属度矩阵中一定会有异常值出现.而通过均值方法计算先验信息会导致{\pi }_i出现偏大或者偏小的情况，这是因为如果在一个二维矩阵中存在一些较大或较小的数值，计算其均值就必须将矩阵中的所有数值相加求和，得到的结果会因其中大数或小数（离群值）的存在而偏大或偏小，这样便会导致其平均的结果也变大或变小，从而最终影响先验信息的准确性.为增强KL信息的准确性，本文使用ESD（Extreme Studentized Deviate）^［16］数据异常值检测模型剔除先验信息中的离群值，从而改善隶属度先验信息.通过数据直方图分析法，得出隶属度数值服从正态分布，满足使用ESD模型的先决条件.

基于ESD模型的增强KL信息方法，首先计算隶属度矩阵的平均值\overline{U}和方差{\sigma }_U，然后通过式（16）计算出临界值R_x，再次通过式（17）检验统计量{\lambda }_x，p=\mathrm{1}-a/\left[\mathrm{2}(n-x+\mathrm{1})\right].最后比较检验统计量{\lambda }_x与临界值R_x的大小，如果R_x>{\lambda }_x，则可以找出隶属度矩阵中的离群值x，这时只需将异常值剔除并得到新的隶属度矩阵\widetilde{u_{ij}}.由式（15）可以看出，\widetilde{u_{ij}}是将上一次迭代运算的隶属度矩阵u_{ij}^{(a-\mathrm{1})}进行纠正得到的，将N-x作为样本总数，并通过式（18）计算得到修正的隶属度先验信息，最后将修正结果h({\pi }_i^{(a)})作为隶属度先验信息加入KL散度中.上述过程用公式表示为：

R_x=\frac{\mathrm{m}\mathrm{a}{\mathrm{x}}_i\left|U_i-\overline{U}\right|}{{\sigma }_U}（16）

{\lambda }_x=\frac{(n-x)t_{p,n-x-\mathrm{1}}}{\sqrt[]{(n-x-\mathrm{1}+t_{p_{,n-x-\mathrm{1}}}^{\mathrm{2}})(n-x+\mathrm{1})}}（17）

h({\pi }_i^{(a)})=\frac{\mathrm{1}}{N-x}{\displaystyle \sum _{j=\mathrm{1}}^k}\widetilde{u_{ij}}（18）

图4表示算法在收敛之后隶属度值的分布结果，其中，图4（a）左侧为噪声图像、右侧为它红线框中对应3×3区域内的灰度值；图4（b）为原始FCM算法的隶属度值；图4（c）为文献［14］提出引入KL信息项后的隶属度值；图4（d）为本文增强KL信息项的隶属度值.FCM算法的隶属度值在［0，1］之间，理论上，正确聚类簇后，图4（b）、（c）、（d）中的第一列（3×3矩阵）隶属度值至少应该大于0.5或最好大于0.9，而原始FCM算法因很难将受到噪声污染的图像正确分割，其隶属度值均低于0.5.文献［14］提出的算法虽对噪声有一定的抑制能力，但隶属度中仍有异常值的出现（图中红色线框标记），这说明隶属度先验信息也会受到噪声的影响.本文算法的隶属度值均大于0.96，且没有离群值的出现，这表明基于ESD模型增强的KL信息可以有效消除隶属度先验信息的离群值，减少隶属度模糊现象.

引入增强KL信息的目标函数如下：

J_{\mathrm{E}\mathrm{S}\mathrm{D}}^{(a)}={\displaystyle \sum _{i=\mathrm{1}}^K}{\displaystyle \sum _{j=\mathrm{1}}^N}\left(u_{ij}^{(a)}{‖y_j-c_i^{(a)}‖}^{\mathrm{2}}\right)+

\phi {\displaystyle \sum _{i=\mathrm{1}}^K}{\displaystyle \sum _{j=\mathrm{1}}^N}{u}_{ij}^{(a)}\mathrm{l}\mathrm{n}\frac{u_{ij}^{(a)}}{h({\pi }_{ij}^{\left(a\right)})},{\displaystyle \sum _{i=\mathrm{1}}^K}{u}_{ij}=\mathrm{1}（19）

2.3 本文算法目标函数

为提升FCM算法用于含噪声多光谱遥感图像分割的精度，提出一种自适应模糊子空间与增强KL信息相结合的FCM多光谱遥感图像分割算法，目标函数如下：

\begin{array}{l}{J}_{}^{(a)}={\displaystyle \sum _{i=\mathrm{1}}^K}{\displaystyle \sum _{j=\mathrm{1}}^N}\left[{\displaystyle \sum _{z=\mathrm{1}}^D}{\left(w_{iz}^{(a)}\right)}^{\tau }\left(\alpha u_{ij}^{(a)}{‖y_{jz}-c_{iz}^{(a)}‖}^{\mathrm{2}}+\beta G_{ij}\right)\right]+\\ \gamma {\displaystyle \sum _{i=\mathrm{1}}^k}{\displaystyle \sum _{z=\mathrm{1}}^D}{\left(w_{iz}^{(a)}\right)}^{\tau }+\phi {\displaystyle \sum _{i=\mathrm{1}}^K}{\displaystyle \sum _{j=\mathrm{1}}^N}{u}_{ij}^{(a)}\mathrm{l}\mathrm{n}\frac{u_{ij}^{(a)}}{h({\pi }_{ij}^{\left(a\right)})},\\ {\displaystyle \sum _{i=\mathrm{1}}^K}{u}_{ij}=\mathrm{1},{\displaystyle \sum _{z=\mathrm{1}}^D}{w}_{iz}=\mathrm{1}\end{array}（20）

对式（20）使用拉格朗日乘子法得到含约束条件的拉格朗日函数如下：

\begin{array}{l}F={\displaystyle \sum _{i=\mathrm{1}}^K}{\displaystyle \sum _{j=\mathrm{1}}^N}\left[{\displaystyle \sum _{z=\mathrm{1}}^D}{\left(w_{iz}^{(a)}\right)}^{\tau }\left(\alpha u_{ij}^{(a)}{‖y_{jz}-c_{iz}^{(a)}‖}^{\mathrm{2}}+\beta G_{ij}\right)\right]+\\ \gamma {\displaystyle \sum _{i=\mathrm{1}}^k}{\displaystyle \sum _{z=\mathrm{1}}^D}{\left(w_{iz}^{(a)}\right)}^{\tau }+\phi {\displaystyle \sum _{i=\mathrm{1}}^K}{\displaystyle \sum _{j=\mathrm{1}}^N}{u}_{ij}^{(a)}\mathrm{l}\mathrm{n}\frac{u_{ij}^{(a)}}{h({\mathrm{\pi }}_{ij}^{\left(a\right)})}+\\ {\displaystyle \sum _{j=\mathrm{1}}^N}{\lambda }_{j\mathrm{1}}\left(\mathrm{1}-{\displaystyle \sum _{i=\mathrm{1}}^k}{u}_{ij}^{(a)}\right)+{\displaystyle \sum _{j=\mathrm{1}}^N}{\lambda }_{j\mathrm{2}}\left(\mathrm{1}-{\displaystyle \sum _{z=\mathrm{1}}^D}{w}_{iz}^{(a)}\right)\end{array}（21）

式中：{\lambda }_{j\mathrm{1}}和{\lambda }_{j\mathrm{2}}为拉格朗日乘子在多通道情况下的约束参数.式（21）中的G_{ij}应表示为：\sum _{r\in N_j}\frac{\mathrm{1}}{d_{rj}+\mathrm{1}}{\left(\mathrm{1}-u_{ir}^{(a)}\right)}^m{‖y_{rz}-c_{iz}^{(a)}‖}^{\mathrm{2}}.接下来，将函数F分别对拉格朗日乘子{\lambda }_{j\mathrm{1}}和{\lambda }_{j\mathrm{2}}求偏导，并令其等于0，得到下面两个约束条件：

\frac{\partial F}{\partial {\lambda }_{j\mathrm{1}}}=\mathrm{0}\Rightarrow {\displaystyle \sum _{i=\mathrm{1}}^K}{u}_{ij}^{(a)}=\mathrm{1};\frac{\partial F}{\partial {\lambda }_{j\mathrm{2}}}=\mathrm{0}\Rightarrow {\displaystyle \sum _{i=\mathrm{1}}^K}{w}_{ij}^{(a)}=\mathrm{1}（22）

然后，对式（21）关于u_{ij}^{(a)}求偏导，得到：

\frac{\partial F}{\partial u_{ij}^{(a)}}=Q_{ij}+\phi \mathrm{l}\mathrm{n}\frac{u_{ij}^{\left(a\right)}}{{\pi }_{ij}^{\left(a\right)}}+\phi -{\lambda }_{j\mathrm{1}}（23）

为了让公式简洁直观，将重复的多项式用字母Q_{ij}表示，其中，V_{rj}=\frac{\mathrm{1}}{d_{rj}+\mathrm{1}}(\mathrm{1}-u_{ir}).

Q_{ij}={\displaystyle \sum _{z=\mathrm{1}}^D}{w}_{iz}^{\tau }\left[\alpha {‖y_{jz}-c_{iz}^{(a)}‖}^{\mathrm{2}}-\beta \sum _{r\in N_j}{V}_{rj}{‖y_{rz}-c_{iz}^{(a)}‖}^{\mathrm{2}}\right]（24）

\widetilde{Q_{ij}}={\displaystyle \sum _{z=\mathrm{1}}^D}{w}_{iz}^{\tau }\left[\alpha {‖y_{jz}-c_{sz}^{(a)}‖}^{\mathrm{2}}-\beta \sum _{r\in N_j}{V}_{rj}{‖y_{rz}-c_{sz}^{(a)}‖}^{\mathrm{2}}\right]（25）

令式（23）等于0，得到：

u_{ij}^{(a)}={\pi }_{ij}^{(a)}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{{\lambda }_{j\mathrm{1}}}{\phi }-\mathrm{1}-\frac{Q_{ij}}{\phi }\right)（26）

最后，将第一约束条件代入式（26）得到：

u_{ij}^{(a)}=\frac{{\pi }_{ij}^{(a)}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{1}-\frac{Q_{ij}}{\phi }\right)}{{\displaystyle \sum _{s=\mathrm{1}}^K}{\pi }_{is}^{(a)}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{1}-\frac{\widetilde{Q_{ij}}}{\phi }\right)}（27）

接下来，让函数F对w_{iz}^{(a)}求偏导，并令其等于0，再将第二约束条件代入，得到式（28）.因该求解过程与u_{ij}^{(a)}类似，因此不再赘述推导过程.

w_{iz}^{(a)}={\displaystyle \sum _{s=\mathrm{1}}^D}\left[\frac{{\displaystyle \sum _{j=\mathrm{1}}^N}\alpha {u_{ij}}^{(a)}{‖y_{jz}-c_{iz}^{(a)}‖}^{\mathrm{2}}+\beta G_{ij}+\gamma }{{\displaystyle \sum _{j=\mathrm{1}}^N}\alpha {u_{ij}}^{(a)}{‖y_{js}-c_{is}^{(a)}‖}^{\mathrm{2}}+\beta \widetilde{G_{ij}}+\gamma }\right]（28）

式中：G_{ij}=\sum _{r\in N_j}{V}_{rj}{\left(\mathrm{1}-u_{ir}^{(a)}\right)}^m{‖y_{rz}-c_{iz}^{(a)}‖}^{\mathrm{2}}，\widetilde{G_{ij}}=\sum _{r\in N_j}{V}_{rj}{\left(\mathrm{1}-u_{ir}^{(a)}\right)}^m{‖y_{rs}-c_{is}^{(a)}‖}^{\mathrm{2}}.最终，让函数F对c_{iz}^{(a)}求偏导，并令其等于0，得到c_{iz}^{(a)}的递推式如下：

c_{iz}^{(a)}=\frac{{\displaystyle \sum _{j=\mathrm{1}}^N}{\displaystyle \sum _{z=\mathrm{1}}^D}{\left({w_{iz}}^{(a)}\right)}^{\tau }\left[\alpha {u_{ij}}^{(a)}{y}_{jz}+\beta \sum _{r\in N_j}{V}_{rj}{y}_{rz}\right]}{{\displaystyle \sum _{j=\mathrm{1}}^N}{\displaystyle \sum _{z=\mathrm{1}}^D}{\left({w_{iz}}^{(a)}\right)}^{\tau }\left[\alpha {u_{ij}}^{(a)}+\beta \sum _{r\in N_j}{V}_{rj}\right]}（29）

通过以上u_{ij}^{(a)}、{w_{iz}}^{(a)}、c_{iz}^{(a)}的递推公式，就可以实现目标函数最小化，并依据每个像素的隶属度，实现样本的同质分区（簇），最终达到图像分割的目的.

本文算法详细步骤如下.

步骤1：参数设置，包括聚类数目K，最小误差e，最大迭代次数T，模糊隶属度指数m，子空间约束参数\gamma，正则项控制参数\phi；

步骤2：初始化隶属度矩阵u_{ij}^{(a)}，并令a=\mathrm{1}；

步骤3：利用式（29）计算聚类中心c_{iz}^{(a)}；

步骤4：利用式（8）计算模糊因子G_{ir}；

步骤5：利用式（9）~（11）计算\alpha和\beta的值；

步骤6：利用式（28）计算子空间加权因子w_{iz}^{(a)}；

步骤7：利用式（16）~（18）计算隶属度先验信息h({\pi }_i^{(a)})；

步骤8：利用式（27）更新隶属度矩阵u_{ij}^{(a)}；

步骤9：利用式（20）计算目标函数J_{}^{(a)}；

步骤10：判断是否成立，如成立则退出循环，否则返回步骤3；

步骤11：将每个像素划分到不同的聚类簇中完成图像分割.

3 实验结果与分析

为验证算法的实用性和有效性，采用真实多光谱遥感图像与AID数据库^［17］中的遥感图像进行实验，并使用FCM^［5］、ADFLICM^［18］、CGFFCM^［19］、FCM_SCIM^［15］、FRFCM^［20］和FSC_LNML^［21］算法进行对比实验.本实验在MATLAB上执行，台式MSI的配置为Intel（R） Core （TM） CPU、i9-9900、4.7GHz和16GB RAM.

为测试算法抑制噪声的能力，将混合噪声添加到全部的测试图像中，然后通过分割对比实验验证每种算法的抗噪声能力.模拟噪声设计如下.

1） 输入原始图像I.

2） 将输入的原始图像转化为“double”类型的数据，将所有像素归一化在［0，1］区间.

3） 设置噪声密度\sigma.

4） 采用平均值为0，方差为\sigma的高斯噪声与图像I混合.

5） 采用密度为\sigma的椒盐噪声与图像I混合.

6） 采用密度为\sigma的均匀乘性噪声与图像I混合.

7） 将图像I的像素规范化为［0，255］之间，并将图像的数据类型转换为“uint8”.

8） 在接下来的实验中，所有的被噪声污染的图像描述为：混合\sigma \times \mathrm{100}\mathrm{\%}的噪声.

3.1 分割评价指标

分割评价指标采用模糊分割系数（V_PC）、模糊分割熵（V_PE）、分割精度（SA）、平均交互比（mIoU）、峰值信噪比（PSNR）.V_PC是衡量隶属度紧密程度的一个指标，其数值越大，表示每个像素的属性越独特，被聚类的程度越好；模糊分割熵V_PE表示分割的模糊性，其数值越小，表示被聚类划分像素的模糊性越小；SA和mIoU是图像分割中的常用指标，SA和mIoU的值越大，也表示分割结果越好；PSNR可用来衡量分割结果与分割参考图像之间的差异，PSNR数值越大，表示分割结果与参考结果越接近.

3.2 多光谱遥感图像分割实验

为测试算法的性能，分别选取AID数据库中的4幅，真实场景中的2幅，大小都为600×600的多光谱遥感图像进行实验，如图5所示.图5（a）~（d）为AID数据库图像，图5（e）和（f）为真实场景图像，之后用混合5%的噪声破坏它们，如图6所示.分割后的视觉效果如图7和图8所示.它们的噪声变化测试如表1所示.

视觉效果上，FCM、ADFLICM和CGFFCM的分割结果中充满噪声，FRFCM的分割结果多数被噪声淹没.这是因为FCM和CGFFCM只考虑图像的自身信息，缺少空间信息的引入，导致它们对噪声的抑制能力不足，使得分割结果受到较大的影响.ADFLICM和FRFCM虽然都考虑了图像的空间信息，但它们缺少空间信息和通道特征的参数自适应环节，因此，它们在处理多目标和地物特征相对复杂的多光谱遥感图像时，其抑制噪声的能力和分割精度都不能达到令人满意的结果.而本文算法、FCM_SCIM和FSC_LNML的分割结果都相对清晰，这是因为这类算法不光结合图像空间信息，而且增加了自适应参数，所以它们对噪声的抑制效果较好.图7中，本文算法在抑制绿色森林区域的噪声方面效果更好，并将林地与道路很好地区分，在处理道路、桥梁和广场时，可将桥梁、普通道路和圆形道路完整地分割.而FCM_SCIM和FSC_LNML算法，它们都不能对周围的部分小目标进行细节分割，部分图像小的细节被误分割.从真实场景的分割结果看（图8），FCM_SCIM的分割结果十分模糊，FSC_LNML对噪声的抑制能力明显下降，并存在边缘模糊的情况；而本文算法不仅可以很好地抑制噪声，而且将图中的房屋和道路正确分割，而它们的边缘细节也相对完整.从表1反映的数据看，当噪声的密度增加时，其他算法的SA和PSNR指标不断下降，但本文算法仍然可以保持较好的分割精度和抗噪声性能，再次验证本文算法对噪声多光谱遥感图像有更好的分割性能.

为进一步验证本文算法的性能，使用以上各类算法对AID数据库中随机选取的52幅图像进行对比实验，分别混入5%，7.5%，10%，15%的噪声，并计算各类指标的平均值，实验结果如表2和表3所示（时间单位：s）.从表中可以看出，随着噪声密度增加，各类算法的分割精度都会下降，运算时间会变长.本文算法与FCM_SCIM和FCM_LNML运算时间近似，但抑制噪声效果更好.在混有10%噪声的情况下，本文算法的分割精度仍然可以达到92%以上，mIoU数值也维持在较高水平，而FCM_SCIM和FCM_LNML的分割精度只能达到90%或低于90%.这表明本文算法所使用的自适应模糊因子的噪声去除能力更为突出.此外，在两个模糊系数（V_PC和V_PE）指标上，本文算法仍然具有一定优势，这再次验证本文所提的增强KL信息方法可以减少隶属度的模糊性，进而优化算法的分割能力.

4 结 论

本文提出了一种自适应模糊子空间与增强KL信息相结合的FCM多光谱遥感图像分割算法，针对传统FCM算法难以分割含噪声多光谱遥感图像的问题，通过使用自适应局部模糊因子达到对噪声的抑制目的.并使用模糊子空间理论将原始图像信息和空间修正信息整合到聚类中，进而实现了图像通道权重自适应计算.最后，在目标函数中引入了增强KL信息正则项以消除隶属度中的离群值，减少隶属度模糊现象.模拟混合噪声的多光谱遥感图像实验证明，本文算法在抑制噪声能力方面更强.此外，相比FCM_SCIM和FSC_LNML，本文算法在满足去除噪声的同时，仍有着较高的分割精度.本文算法虽然在抗噪声性能和分割精度方面取得了一定表现，但算法的计算复杂度较高，运算时间较长，使用的噪声仅为模拟噪声，下一步研究目标是如何降低算法计算复杂度和缩短运算时间以及如何在真实噪声环境中让算法发挥更好的性能.

参考文献

[1] ZHAO X M， WANG H J， WU J，et al ．Remote sensing image segmentation using geodesic-kernel functions and multi-feature spaces［J］．Pattern Recognition，2020，104：107333．

[2] 赵小强，李希尧 ．基于信息蒸馏级联伸缩网络的图像超分辨率重建［J］．湖南大学学报（自然科学版），2023，50（8）：52-61．

[3] ZHAO X Q， LI X Y ．Image super-resolution reconstruction based on information distillation cascade expansion and compression network［J］．Journal of Hunan University （Natural Sciences），2023，50（8）：52-61．（in Chinese）

[4] 杨宏宇，杨帆 ．基于图像去噪和图像生成的对抗样本检测方法［J］．湖南大学学报（自然科学版），2023，50（8）：72-81．

[5] YANG H Y， YANG F ．Adversarial example detection method based on image denoising and image generation［J］．Journal of Hunan University （Natural Sciences），2023，50（8）：72-81．（in Chinese）

[6] HOSSAIN M D， CHEN D M ．Segmentation for Object-Based Image Analysis （OBIA）：a review of algorithms and challenges from remote sensing perspective［J］. ISPRS Journal of Photogrammetry and Remote Sensing，2019，150：115-134．

[7] DUNN J C ．A fuzzy relative of the ISODATA process and its use in detecting compact well-separated clusters［J］．Journal of Cybernetics，1973，3（3）：32-57．

[8] 杨蕴，李玉，赵泉华 ．基于局部空间信息的可变类模糊阈值光学遥感图像分割［J］．自动化学报，2022，48（2）：582-593．

[9] YANG Y， LI Y， ZHAO Q H ．Fuzzy threshold optical remote sensing image segmentation with variable class number based on local spatial information［J］．Acta Automatica Sinica，2022， 48（2）：582-593．（in Chinese）

[10] 王小鹏，魏统艺，房超，等 ．自适应非局部空间约束与K-L信息的模糊C-均值噪声图像分割算法［J］．控制理论与应用，2022，39（7）：1261-1271．

[11] WANG X P， WEI T Y， FANG C，et al ．Fuzzy C-means algorithm with adaptive non-local spatial constraints and K-L information for noisy image segmentation［J］．Control Theory & Applications，2022，39（7）：1261-1271．（in Chinese）

[12] ZHANG H， BRUZZONE L， SHI W Z，et al ．Enhanced spatially constrained remotely sensed imagery classification using a fuzzy local double neighborhood information C-means clustering algorithm［J］．IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2018，11（8）：2896-2910．

[13] 林文杰，李玉，赵泉华 ．结合MST划分和RHMRF-FCM算法的高分辨率遥感图像分割［J］．测绘学报，2019，48（1）：64-74．

[14] LIN W J， LI Y， ZHAO Q H ．High-resolution remote sensing image segmentation using minimum spanning tree tessellation and RHMRF-FCM algorithm［J］．Acta Geodaetica et Cartographica Sinica，2019，48（1）：64-74．（in Chinese）

[15] XU J D， FENG G Z， ZHAO T Y，et al ．Remote sensing image classification based on semi-supervised adaptive interval type-2 fuzzy c-means algorithm［J］．Computers and Geosciences，2019，131：132-143．

[16] WANG J， WANG S T， CHUNG F，et al ．Fuzzy partition based soft subspace clustering and its applications in high dimensional data［J］．Information Sciences，2013，246：133-154．

[17] KRINIDIS S， CHATZIS V ．A robust fuzzy local information C-means clustering algorithm［J］．IEEE Transactions on Image Processing： a Publication of the IEEE Signal Processing Society，2010，19（5）：1328-1337．

[18] ICHIHASHI H， MIYAGISHI K， HONDA K ．Fuzzy C-means clustering with regularization by K-L information［C］//10th IEEE International Conference on Fuzzy Systems．Melbourne，VIC，Australia． IEEE，2001：924-927．

[19] GHARIEB R R， GENDY G， ABDELFATTAH A，et al ．Adaptive local data and membership based KL divergence incorporating C-means algorithm for fuzzy image segmentation［J］．Applied Soft Computing，2017，59（C）：143-152．

[20] WANG Q S， WANG X P， FANG C，et al ．Fuzzy image clustering incorporating local and region-level information with Median memberships［J］．Applied Soft Computing，2021，105：107245．

[21] JAIN R B ．A recursive version of Grubbs’ test for detecting multiple outliers in environmental and chemical data［J］．Clinical Biochemistry，2010，43（12）：1030-1033．

[22] XIA G S， HU J W， HU F，et al ．AID：a benchmark data set for performance evaluation of aerial scene classification［J］．IEEE Transactions on Geoscience and Remote Sensing，2017，55（7）：3965-3981．

[23] ZHANG H， WANG Q M， SHI W Z，et al ．A novel adaptive fuzzy local information C-means clustering algorithm for remotely sensed imagery classification［J］．IEEE Transactions on Geoscience and Remote Sensing，2017，55（9）：5057-5068．

[24] GOLZARI OSKOUEI A， HASHEMZADEH M， ASHEGHI B，et al ．CGFFCM：cluster-weight and Group-local Feature-weight learning in Fuzzy C-means clustering algorithm for color image segmentation［J］．Applied Soft Computing，2021，113：108005．

[25] LEI T， JIA X H， ZHANG Y N，et al ．Significantly fast and robust fuzzy C-means clustering algorithm based on morphological reconstruction and membership filtering［J］．IEEE Transactions on Fuzzy Systems，2018，26（5）：3027-3041．

[26] WEI T Y， WANG X P， LI X N，et al ．Fuzzy subspace clustering noisy image segmentation algorithm with adaptive local variance & non-local information and mean membership linking［J］．Engineering Applications of Artificial Intelligence，2022，110：104672．

基金资助

国家自然科学基金资助项目(61761027)

National Natural ScienceFoundation of China(61761027)

兰州市科技计划资助项目