基于模糊C均值聚类算法的去云处理

许 晨,徐 洋,康 雪,吕达仁

(1.成都市气象局，四川 成都 610000；2.四川省气象局，四川 成都 610000;3.中国科学院大气物理研究所，北京 100029)

在遥感数据中，薄云覆盖不仅会影响遥感数据的解译精度，严重时还会导致遥感数据不能得到有效利用。因此，去云处理是遥感图像改善技术的一个重要方面。本文重点探讨云区数据聚类的有效性，通过对同一轨不同通道的云区进行聚类处理，将云区与植被、土壤、水域等进行分离，生成云区和非云区数据，再运用中值滤波和云区替换的方法实现去云。

1 资料来源与方法

1.1 资料来源

选取NOAA-17卫星的1、4、5通道，分别对应可见光和两个远红外通道，如表1所示。其中，通道1获得白天云图及地表数据；通道4、5获得昼夜云图、海表、地表温度和土壤湿度[1]。

表1 NOAA-17卫星1、4、5通道的波长范围

1.2 模糊C均值聚类

1.2.1 算法介绍

模糊C均值(FCM，Fuzzy C-Means)算法是基于目标函数的聚类法，每个数据点属于某个聚类的程度用隶属度来确定，利用均方逼近理论构造带约束的非线性规划函数，该算法设计简单、容易实现，能借助非线性理论解决优化算法，应用面广[2-3]。

(1)隶属度函数

一个对象X隶属于集合A的程度用隶属度函数uA(X)表示，0≤uA(X)≤1，其自变量范围是所有可能属于集合的对象(即集合所在空间中的所有点)[4]，uA(X)=1代表X∈A，说明向量完全属于集合A；uA(X)=0代表XA，说明X完全不属于集合A。

(2)模糊集合

设X={x}是一个集合，uA(X)定义在X上，并且uA(X)∈[0,1]，则X上的一个模糊集合wA是指：对于任意x∈X，都确定了一个数uA(X)，称uA(X)为X对wA的隶属度。模糊集合完全由其隶属函数刻画，使一个元素隶属的集合类别非硬性[5]。

(3)模糊集合的C划分

把n个模式中的观测样本集X={x1,x2,…,xn}分为C个模糊组，求每个模糊组的聚类中心，使非相似性的价值函数最小，根据每个数据点的隶属函数得到样本属于各个类别的不确定性程度[6]。

根据隶属度定义可知，一个数据集的隶属度之和为1，即

(1)

因此，FCM的目标函数(价值函数)表示为：

(2)

其中uij∈[0,1]，ci表示模糊组的第i个聚类中心，dij=‖ci-xj‖是ci与第j个数据点间的欧几里德距离，m是加权指数,m∈[1,∞)[7]。

接下来求使式(2)达到最小值的必要条件，构造如下新目标函数：

(3)

其中，式(1)的n个约束式的拉格朗日乘子表示为λj(j=1,2,…,n)。对所有输入参量求导后得到使式(3)达到最小值的必要条件为[8]：

(4)

(4)加权指数m

在FCM算法中，加权指数m非常重要，它控制隶属度的分配和聚类模糊程度，m∈[1,∞)。m的选取以合理聚类为最终目标，要求FCM算法同时满足两个条件：实现模糊聚类及进一步将数据集划分清晰，这样才能准确分辨样本集中各对象的类属关系。关于参数m的优选方法可表示为：

(5)

Jm是聚类目标函数，U、P是不同m取值下FCM算法得到的最优划分，将最佳加权指数m的确定转化为一个带约束的非线性规划问题[9]。

1.2.2 FCM实现

(1)初始化[10]：设置聚类类别数c，2≤c≤n，n表示数据个数。设置迭代停止阈值ε，将聚类原型模式P(0)初始化，将迭代计数器设为b=0。

(2)用值在[0,1]的随机数初始化划分矩阵U，使其满足式(1)的约束条件。

(3)用式(4)更新划分矩阵U(b)。

(4)根据式(4)更新聚类原型模式矩阵P(0)。

(5)重复以上过程进行迭代，若‖c(b)-c(b+1)‖<ε，算法停止并输出聚类原型P和划分矩阵U。否则令b=b+1，转向步骤(3)。其中,‖·‖为某种合适的矩阵范数[11]。

上述算法也可先初始化模糊划分矩阵，然后再进行迭代。由于不能确保FCM收敛于一个最优解，因此算法的性能依赖于初始聚类中心[12]。

综上，FCM算法的聚类结果受参数控制，其主要影响参数分别是：聚类数目c，加权系数m和迭代阈值ε。经过FCM算法输出两个量：一个c*n的模糊划分矩阵，以及c个聚类中心点向量。前者表示各样本点属于每个类的隶属度，可根据模糊集合中的最大隶属原则确定各样本点的归类。聚类中心则是这个类的代表点，表示每个类的平均特征[13]。

1.3 去云算法实现

鉴于模糊C均值算法的聚类有效性，本文将其与去云技术有机地联系起来对遥感数据进行聚类处理，并在此基础上进行类属提取，将云与植被、土壤、水域等进行分离，生成云区和非云区数据，再运用中值滤波和云区替换的方法实现去云，如图1所示。

图1 算法整体实现框图

1.3.1 聚类初始化

读取HDF格式的NOAA-17遥感数据，显示通道1、4、5。据相关资料分析[14-15]，可找出NOAA卫星在可见光通道绿色植被的最小反射率和云区最高反射率之间的相对关系和特征，以此来确定云区反射率阈值，如表2所示。

表2 NOAA卫星在可见光波段的反射率判云阈值表(%)

待聚类的有云数据是灰度数据，将控制聚类结果的加权系数m转化为灰度值，并利用云、植被、水域、土壤等在三通道的灰度值的固有特性，将其聚成不同的类别，如表3所示。

表3 卫星数据反照率值

先将初始聚类中心灰度值设定在反射率灰度阈值和亮温阈值范围内，此时要针对不同季节和不同地区的地形特征进行具体分析，这样能使聚类效果更加准确，聚类速度更快。并将设定的初始聚类中心灰度值作为初始聚类的起始位置，计算聚类迭代次数，生成N个类别。

读入数据流程如图2所示，单通道归一化主要对图像像素做处理，以保证选择的待聚类数据有云且云区较清楚。

图2 读数据流程图

1.3.2 云区聚类

先读取待聚类数据，设置参数进行云区聚类。并分析聚类后数据，若不佳则改变参数，重新迭代聚类。然后提取类属，为去云做准备(图3、图4)。

图3 FCM算法初始化

图4 聚类提取云区

1.3.3 云区处理

鉴于云区的覆盖面积小，且它和椒盐噪声较类似，故选取对椒盐噪声有很好消除作用的中值滤波法处理云区。中值滤波的优点在于既能过滤噪声，又能保证细节不模糊、很好的保护边缘轮廓信息。它消除噪声的效果与模板尺寸和参与运算的像素数有关。其实现方法为[16]：首先用一个窗口在图像上扫描，并将窗口中心与待处理像素位置相重合，读取窗口内包含的图像像素的灰度值，将其按升序排列起来，取灰度值居中的像素灰度为窗口中心像素的灰度。

待中值滤波处理云区后，将有云区域进行替换，得到去云数据。需注意的是，在中值滤波处理云区时，使用的模板越大，去云会呈现一定的模糊。流程如图5所示。

图5 云区处理流程图

2 结果与分析

2.1 读取数据

通过ENVI软件读取HDF格式的NOAA-17遥感数据，选取7月25日在四川盆地的部分区域数据，并合成1、4、5通道数据，如图6所示。

图6 读取NOAA-17数据

分析图6，四幅图中的亮色区域都较小，故认为此亮色区域全为云区，忽略雪山等非云因素的影响。其中，(a)、(b)、(c)单通道数据很好地体现了云的特征：在可见光通道1中，白色区域为云区；在红外通道4、5中，亮色的小突起部分为云区；在合成通道数据中，亮色区域为云区。将上述四组数据作为待聚类数据。

2.2 数据聚类

采用FCM算法对图6的数据进行聚类。初始化聚类中心时，依据云在可见光和红外通道表现出的反射率差异大小，即数据中的灰度值差异，对数据进行聚类，效果如图7。

图7 FCM算法进行数据聚类

对比图7分析发现：对单通道的聚类结果(a)、(b)、(c)能真实反映聚类信息，而对合成通道的聚类结果(d)则在某些位置会产生信息重叠，造成信息量的部分丢失。

2.3 云区检测与提取

对聚类后数据进行云区检测和提取，云区检测的目的是保证云区信息量不丢失，保证其完整性。如图8所示。

图8 云区检测与提取

图8对比分析，可得三点结论：

(1)在单通道云区数据(a)、(b)、(c)中，(b)中通道4在图像左下方出现与(a)、(c)不相似的情况，这时参考云在三个通道的综合特征，忽略通道4左下部分的异常情况。

(2)合成通道云区数据图8(d)是从图7(d)中直接提取的云区，对比图8(e)发现，它丢失了少量云区信息。

(3)三通道云区交集数据(e)是三通道各自云区信息取其交集而得，真实全面地反映了云在合成通道的特征。

2.4 云区去除

将通道1作为辅助数据，通道4、5作为基准数据。然后采用中值滤波和替换融合处理通道4、5的云区数据，如图9所示。

图9 去云结果

如图9所示，(a)、(b)是对通道4、5的云区数据进行中值滤波后的结果，滤波模板为3×3大小。对比图6(b)、(c)分析，云区得到一定程度上的去除，但呈现一定模糊，丢失了细节信息；(c)、(d)是先对三通道进行中值滤波、再进行通道合成后的数据，分别采用3×3大小和7×7大小的滤波模板。可看出，在滤波过程中，采用的模板越大，图像越模糊；(e)、(f)是先对三通道进行中值滤波、然后参考三通道交集云区数据图8(e)，将通道1作为辅助数据、采用替换融合后的去云结果，分别对应3×3大小和7×7大小的滤波模板。可以看到，通过此方法，有效地对原数据进行了去云处理。

3 结论

提出了一种简单有效的多时相去云方法，通过模糊C均值聚类算法对云区数据聚类，再采用中值滤波、云区替换达到去云效果。(1)在对数据进行聚类时，对单通道的聚类结果能真实反映实际信息，而对合成通道的聚类结果则在某些位置会产生信息重叠，造成信息量的部分丢失。(2)进行云区检测和提取时，单通道云区信息可能出现不相似情况，这时参考云在三通道的综合特征，忽略异常情况。三通道云区交集数据较合成通道数据更能真实全面反映云在合成通道的特征。(3)对去云后的单通道数据直接采用中值滤波，云区虽得到较好去除，但呈现一定模糊，丢失了细节信息；在滤波过程中，采用的模板越大，图像越模糊；将中值滤波后的结果进行云区替换，使细节更加清晰。此外，本文的算法只针对薄云区域。