CN103034858A - 一种卫星云图的二次聚类分割方法 - Google Patents

Abstract

本发明公开了一种卫星云图的二次聚类分割方法，该方法首先将大范围的整体卫星云图进行分块化处理；其次依次对个子区域云图每个样本点提取对应的多通道光谱特征及TPLBP纹理特征进行精细的初次核聚类分割，从而得到若干子区域云图分割结果；最后在初次聚类分割结果的基础上，利用其作为先验知识提取初次聚类结果各类对应于原子区域云图范围内的灰度均值特征和DI特征进行全局云图的二次核聚类分割，以保证云分类的完整性；优点是精确度高、鲁棒性强，可以针对大范围气象卫星云图所包含的庞大繁杂的地理、海洋和大气等信息进行精细的云类识别。

Description

一种卫星云图的二次聚类分割方法

技术领域

本发明涉及数气象卫星云图的处理技术，尤其是涉及一种卫星云图的二次聚类分割方法。

背景技术

天气预报是气象学中的一个重要应用，它与每个人的日常生活息息相关，甚至关乎到人们的生命财产安全。精确的云分类是天气预报的基础，因此，对卫星云图进行分割而实现正确的云分类就显得至关重要。最早的云图分类是由气象专家通过丰富的气象知识进行肉眼分析判别。然而随着气象卫星技术的发展，大型气象卫星每天向地面发送GB数量级以上的云图数据，以我国风云二号D星(FY2-D)为例，2006年12月份它的发射成功实现了与风云二号C星(FY2-C)的“双星观测”。每颗卫星每半个小时观测一次地球，把双星观测时间错开15分钟，这样重叠区域就可以每隔15分钟接收云图数据并且全天候工作，如此海量的云图数据如果仅靠肉眼、手工来实现云图分割不仅耗时乏味，而且分割结果易受主观因素的影响。因此，人们渐渐把目光聚焦到利用计算机实现自动云分类上。

云图的自动分割主要包括三个步骤：云图预处理、特征提取和分类。其中云图特征的提取与选择及分类方法的选择尤显重要，已成为中外学者关注研究的重点。判别卫星云图中云类的主要依据有：云的结构、云的范围大小、云的边界状况、云的色调、云的纹理及云顶亮温等。同样，这些特征也成为计算机智能云分类的依据，但由于云的结构边界等视觉特征难于描述，因此云的形态特征很少用于云图的分类，所以一般用于计算机自动云分类的云图特征主要为光谱特征及纹理结构特征两大类。在光谱特征上，HongG等人利用红外云图反演云顶亮温用于强对流云团的识别。但此种单一利用红外通道云顶亮温的阈值分割法易受时空变化的影响,很难实现将强对流云从卷云层中分离出来。因此，MecikalskiJR等人提出亮温差的概念用于强对流云团的识别，利用红外一通道与水汽通道亮温差值大于0K，以及红外一通道与红外二通道亮温差值在-25K—5K范围作为阈值对强对流云进行判别，取得较好的效果。在纹理结构特征上，R.M.Welch等人利用灰度共生矩阵(GLCM)提取云图纹理特征用于层积云、积雨云和卷云的分类。K.S Kuo等人利用灰度差分矢量(GLDV)法提取特征用于区分普通云和冰云。

近年来，多特征融合思想逐渐被重视，人们更多的是把光谱及纹理结构特征相结合用于云分类。Tian B与Rashpal K等人利用奇异值分解(SVD)同时提取云图光谱及纹理特征取得较好效果，后期Tian B另辟蹊径采用二维Harr小波包对卫星云图进行3层分解分别提取特征进行云分类。国内周毓荃、陈英英等人利用云顶亮温、云顶高度、云光学厚度等云微观物理特征对云分类进行过深入研究，并取得一系列成果。目前，卫星云图的自动分割及云分类研究处于高速发展时期。云分类中的特征提取，已然从单一阈值过渡到多阈值相结合；从起初的光谱特征发展到光谱特征、纹理特征及云微物理特征相结合。然而，现有的大多数方法仅仅适用于小范围、类别少的云图进行分类，并且大多数方法在抗噪性、分割精确性上不够理想，这将不利于其在气象领域上的实际运用。

发明内容

本发明所要解决的技术问题是提供一种精确度高、鲁棒性强的卫星云图的二次聚类分割方法，该方法可以针对大范围气象卫星云图包含庞大繁杂的地理、海洋和大气等信息的特点，实现对其进行精细的云类识别。

本发明解决上述技术问题所采用的技术方案为：一种卫星云图的二次聚类分割方法，包括以下步骤：

①对卫星云图进行分块处理，将卫星云图分为L×L共L²块大小相等的子区域云图，其中2≤L≤10；

②对各个子区域云图进行初次聚类分割：

②-1将子区域云图中每个像素点作为输入样本点，提取每个子区域云图中所有样本点在多个通道云图中相应区域内对应点的邻域灰度均值，该邻域灰度均值作为相应样本点的光谱特征，其中多个通道云图的分辨率一致；

②-2将每个子区域云图中的某个像素点看做中心生成滑动窗口，将此滑动窗口看做一张子图，计算该滑动窗口内所含的所有像素点的三模块局部二元模式码值，并基于三模块局部二元模式码值构建得到该滑动窗口的三模块局部二元模式码图，统计该滑动窗口的三模块局部二元模式码图的一阶直方图并计算其熵、均值、一致性和三阶距，得到子区域云图中该像素点的二次统计量纹理特征；以上述方法得到每个子区域云图中所有的像素点的二次统计量纹理特征；

②-3采用最大最小距离法自适应确定初始聚类中心及其个数；

②-4将步骤②-1和步骤②-2中每块子区域云图中提取的每个像素点的两个纹理特征和一个灰度均值特征构成的特征矩阵作为输入样本和初始聚类中心通过非线性映射函数映射到高维特征空间，在高维空间中对卫星云图中所有样本点进行模糊核聚类，得到每块子区域云图的初次聚类分割结果；

③对各个子区域云图进行二次聚类分割：

③-1将初次聚类结果中每个类作为输入样本点，提取该输入样本点的类内灰度均值特征；

③-2提取初次聚类分割结果的类内多样性指数特征；

③-3将初次聚类分割结果的类内灰度特征和多样性指数特征作为样本点的二维特征，构成特征矩阵作为二次聚类分割的输入样本，采用最大最小距离法自适应确定最终聚类中心个数；

③-4将二次聚类分割算法的样本点通过非线性映射函数映射到高维特征空间，在高维空间中对所有样本点进行模糊核聚类，最终得到整张卫星云图的二次聚类分割结果。

步骤②-2中得到每个子区域云图的二次统计量纹理特征的具体步骤为：

a、选取子区域云图中的任一像素点O，以该像素点为中心生成(2M+1)×(2M+1)的滑动窗口，此时，将此滑动窗口看做一张子图，在子图中任选一个像素点p，以该像素点p为中心生成大小为3×3的邻域模块，记为C_p，同时在以该模块C_p为中心，半径为r圆周上额外生成s个与该模块C_p大小相同且均匀分布的模块，s个模块分别记为C₀，C₁，……C_j……C_s-1，其中4＜r＜M，r＜M＜(F/L)÷10，s≥2，j=0，1，……，s-1，F表示卫星云图的分辨率；

b、在s个模块中选取C_j和C_j′两个模块，其中j′=（j+a）mod（s），a为整数且2≤a≤s-1，mod表示取两者相除后的余数；

c、计算C_p、C_j和C_j‘的局部二元模式码值，分别记为LBP(C_p)、LBP(C_j)和LBP(C_j′)；

d、利用公式φ(h,q)=LBP(h)-LBP(q)计算LBP(C_j)和LBP(C_j′)与LBP(C_p)的差值， ${\mathrm{\φ}}_{C_j}(C_j,C_p)=\mathrm{LBP}\left(C_j\right)-\mathrm{LBP}\left(C_p\right),$ ${\mathrm{\φ}}_{C_{j^{\′}}}(C_{j^{\′}},C_p)=\mathrm{LBP}\left(C_{j^{\′}}\right)-\mathrm{LBP}\left(C_p\right);$

e、将

和代入公式 $I_p={\mathrm{\&phi;}}_{C_j}(C_j,C_p)-{\mathrm{\&phi;}}_{C_{j^{\&prime;}}}(C_{j^{\&prime;}},C_p)$ 计算得到一组I_p，利用公式 $T\left(I_p\right)=\left\{\begin{array}{c}1,I_p\&GreaterEqual;0\\ 0,I_p=<0\end{array}\right.$ 对I_p进行二值化，得到一组二进制数串T(I_p)；

f、将T(I_p)代入公式 ${\mathrm{TPLBP}}_{r,s,a}\left(p\right)=\underset{j=0}{\overset{s-1}{\mathrm{\&Sigma;}}}T({\mathrm{\&phi;}}_{C_j}(C_j,C_p)-{\mathrm{\&phi;}}_{C_{j^{\&prime;}}}(C_{j^{\&prime;}},C_p))\&times;2^j$ 进行加权计算得到TPLBP_r，s，a(p)，TPLBP_r，s，a(p)即为像素点p的十进制三模块局部二元模式码值；

g、同理计算(2M+1)×(2M+1)大小滑动窗口内所含的所有像素的十进制三模块局部二元模式码值并生成窗口内三模块局部二元模式码图，统计窗口三模块局部二元模式码图的一阶直方图并计算其熵、均值、一致性和三阶距，得到子区域云图中像素点O的二次统计量纹理特征，依次利用滑动窗口对整个子区域云图进行操作，得到子区域云图中每个像素点的三模块局部二元模式二次统计量纹理特征，至此每个像素点可以提取到2个纹理特征，作为初次聚类分割输入样本点的二维特征；

h、按照步骤a~g利用滑动窗口可对每个子区域云图中进行纹理特征的提取。

步骤②-3中采用最大最小距离法自适应确定聚类中心个数的具体步骤为：

a、将卫星云图表示为图像的灰度直方图{H_l}输入，l为图像的像素灰度级，H_l为图像中像素灰度级为l的像素总个数，l0，1，…，255；

b、选择输入样本中任意一样本点的像素灰度级作为第一个已知的聚类中心，记为l₁；

c、将已知聚类中心l₁代入公式d_ll'=|l-l'|计算灰度级为l的样本点与已知聚类中心的距离，其中l'表示已知聚类中心；

d、在{H_l}中找出离l₁距离最远的样本点的像素灰度级的作为第二个已知的聚类中心，记为l₂；

e、分别计算{H_l}中除已知聚类中心对应的样本点以外的剩余样本点到所有的已知聚类中心的距离；

f、比较同一样本点到各已知聚类中心的距离，选择其中的最小值；

g、比较步骤f中所有样本点到各已知聚类中心的距离最小值，选择所有距离最小值中的最大值d_max，若d_max＞E×d_ll'，其中0.5<E<1，则d_max。对应的样本点的像素灰度级作为新的已知的聚类中心，并返回步骤e；若d_max≤E×d_ll'，则统计所有已知的聚类中心，得到初始聚类中心v_g及其数量z，g=1，2，……，z。

步骤②-4中获得卫星云图初次聚类结果的具体步骤为：

a、将图像样本点和初始聚类中心v_g通过非线性映射函数Φ(x)从输入空间映射到高维特征空间，图像中各样本点用x_k表示，k=1，2，……，n，n表示图像包含的样本点的总数，x_k表示图像中第k个样本点；

b、将x_k和v_g代入公式D=||Φ(x_k)-Φ(v_g)||中，计算得到样本点x_k与初始聚类中心v_g之间的欧式距离D；

c、将||Φ(x_k)-Φ(v_g)||²=Φ(x_k)^TΦ(x_k)-2Φ(x_k)^TΦ(v_g)+Φ(v_g)^TΦ(v_g)代入特征空间模糊核聚类算法目标函数 $J_m(U,V)=\underset{g=1}{\overset{z}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m{\left|\right|\mathrm{\&Phi;}\left(x_k\right)-\mathrm{\&Phi;}\left(v_g\right)\left|\right|}^2$ 中，其中U＝{u_gk}，V＝(Φ(v₁)，Φ(v₂)，...,Φ(v_g))，u_gk表示第k个样本点对应于第g类的隶属度，常数m＞1为隶属度矩阵指数，令核函数K(x,y)=Φ(x)^TΦ(y)，则 $J_m(U,V)=\underset{g=1}{\overset{z}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m(K_{\mathrm{kk}}-{2K}_{\mathrm{kg}}+K_{\mathrm{gg}}),$ 其中K_kk=K(x_k，x_k)，K_kg=K(x_k，v_g)，K_gg=K(v_gv_g)，

d、在约束条件

下，计算得到特征空间模糊核聚类隶属度u的函数式为 $u_{\mathrm{gk}}=\frac{{(1/(K_{\mathrm{kk}}-{2K}_{\mathrm{kg}}+K_{\mathrm{gg}}))}^{1/(m-1)}}{\underset{g=1}{\overset{z}{\mathrm{\&Sigma;}}}{(1/(K_{\mathrm{kk}}-{2K}_{\mathrm{kg}}+K_{\mathrm{gg}}))}^{1/(m-1)}}$ 及聚类中心v的函数式为 $v_g=\frac{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m{K}_{\mathrm{kg}}{x}_k}{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m{K}_{\mathrm{kg}}};$

e、初始化参数m、U_gk和核函数K(x,y)，计算得到核函数矩阵K，再对特征空间模糊核聚类隶属度u及聚类中心v进行反复迭代更新，求得最终的U_gk及V_g，根据U_gk及V_g判别得到图像中各样本点的隶属类别。

所述的核函数为高斯核函数、多项式核函数或者Sigmoid核函数。

步骤③-2中提取多样性指数特征的具体步骤为：

a、计算每个子区域云图中感兴趣区域块内灰度的最小值Gray_min及最大值Gray_max，将区间[Gray_min，Gray_max]等间隔的划分成10个区段，区段号分别为N₁，N₂，……，N₁₀，其中，子区域云图初次聚类结果图中一个类范围块所对应于聚类前原子区域云图相同范围内的云图块为感兴趣区域块；

b、计算每个子区域云图中所有的感兴趣区域内的灰度值Gray，将得到的灰度值Gray分别映射到相应的区段，并用对应的区段号代替原灰度值Gray；

c、统计每个区段中各区段号N_t的数量I_t，t=1，2，……，10，设每个子区域云图中所有的感兴趣区域的总像素个数为N'，

将P_t代入公式

得到各子区域云图的多样性指数特征。

与现有技术相比，本发明的优点在于：

1）在整体聚类分割流程上，相比于传统基于整张云图直接聚类分割技术，本发明提出大范围云图的分块二次聚类法首先将大范围的整体云图进行分块化处理；其次，依次对小块云图进行精细的初次核聚类分割，从而得到若干小块云图分割结果；最后，在初次聚类分割的基础上，利用其结果作为先验知识进行全局云图的二次核聚类分割以保证云分类的完整性，这样很大程度上避免了直接利用原图进行分割处理带来的计算量大、分割精确性降低、错分误分等缺陷；

2）在特征的分析选取上，相比于现有单通道光谱特征及共生矩阵提取纹理特征法，本发明在子区域云图初次聚类时，采用样本点对应于此点在多个通道云图相应区域内该点邻域灰度均值作为其光谱特征，即多通道光谱特征；同时，提出一种三模块局部二元模式（TPLBP）纹理特征参与初次聚类分割。TPLBP纹理特征具有良好的旋转不变性，并且能克服光照、位移变化的影响，与多通道光谱特征一样更具鲁棒性，大大提升云分类的正确率；在全局二次聚类中，本发明除提取子区域云图分割结果中各类范围内灰度均值特征外，更引入生态系统多样性研究中的多样性指数(Diversity Index)特征用于二次聚类分割中；DI特征反映的是感兴趣区域块内云的复杂度，能更好的保留云的细节信息；

3）在聚类算法上，本发明结合最大最小距离法、核方法及FCM构造聚类中心个数自适应的特征空间核聚类法，实现了聚类中心个数的初始自适应确定，同时有效的降低了对噪声的敏感性，提高了算法的鲁棒性；

4）本发明二次聚类算法理论上可以实现n(n>0)次聚类扩展，也即在分块化的子云图区域再次进行二次聚类分割实现三次聚类分割，以及扩展到n(n>0)次聚类分割，这样方便了对大范围气象卫星云图的云分类处理。

附图说明

图1为二次聚类法流程框图；

图2a为FY2-D IR1通道云图；

图2b为FY2-D IR1的分块处理图；

图3a为Brodatz纹理图像；

图3b为加入椒盐噪声的图4a Brodatz纹理图像；

图3c为图4a TPLBP码图；

图3d为图4b TPLBP码图；

图4a为晴空陆地样本；

图4b为晴空海面样本；

图4c为密层云样本；

图4d为卷层云样本；

图4e为积雨云样本；

图4f为高层或雨层云样本；

图4g为层积或高积云样本；

图5为TPLBP纹理特征曲线；

图6为DI特征曲线；

图7为3×3块子云图初次聚类分割结果；

图8为本发明算法二次聚类全局云分类结果；

图9为利用整幅云图直接聚类分割结果；

图10为国家卫星气象中心提供专家分类结果图。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

如图1所示，本发明提出了一种卫星云图的二次聚类分割方法，包括以下步骤：

①对卫星云图进行分块处理，将卫星云图分为L×L共L²块大小相等的子区域云图，其中2≤L≤10；

②对各个子区域云图进行初次聚类分割：

②-1将子区域云图中每个像素点作为输入样本点，提取每个子区域云图中所有样本点在多个通道云图中相应区域内对应点的邻域灰度均值，该邻域灰度均值作为相应样本点的光谱特征，其中多个通道云图的分辨率一致；

②-2将每个子区域云图中的某个像素点看做中心生成滑动窗口，将此滑动窗口看做一张子图，计算该滑动窗口内所含的所有像素点的三模块局部二元模式码值，并基于三模块局部二元模式码值构建得到该滑动窗口的三模块局部二元模式码图，统计该滑动窗口的三模块局部二元模式码图的一阶直方图并计算其熵、均值、一致性和三阶距，得到子区域云图中该像素点的二次统计量纹理特征；以上述方法得到每个子区域云图中所有的像素点的二次统计量纹理特征；

②-3采用最大最小距离法自适应确定初始聚类中心及其个数；

②-4将步骤②-1和步骤②-2中每块子区域云图中提取的每个像素点的两个纹理特征和一个灰度均值特征构成的特征矩阵作为输入样本和初始聚类中心通过非线性映射函数映射到高维特征空间，在高维空间中对卫星云图中所有样本点进行模糊核聚类，得到每块子区域云图的初次聚类分割结果；

③对各个子区域云图进行二次聚类分割：

③-1将初次聚类结果中每个类作为输入样本点，提取该输入样本点的类内灰度均值特征；

③-2提取初次聚类分割结果的类内多样性指数特征；

③-3将初次聚类分割结果的类内灰度特征和多样性指数特征作为样本点的二维特征，构成特征矩阵作为二次聚类分割的输入样本，采用最大最小距离法自适应确定最终聚类中心个数；

③-4将二次聚类分割算法的样本点通过非线性映射函数映射到高维特征空间，在高维空间中对所有样本点进行模糊核聚类，最终得到整张卫星云图的二次聚类分割结果。

实施例：本实施例中利用中国国家卫星气象中心提供的FY2-D卫星采集于2012年5月29日凌晨5点30分IR1（红外一通道）通道分区图，辅之以IR2（红外二通道）、IR3（红外三通道）、水汽和VIS四通道分区图进行实验。根据本发明算法对IR1通道卫星云图进行二次聚类分割。

一种卫星云图的二次聚类分割方法，包括以下步骤：

①对如图2a所示的卫星云图进行分块处理，将卫星云图分为L×L共L²块大小相等的子区域云图；如图2b所示，本实施例中将大范围的卫星云图分为3×3共9个大小相等子区域云图供后续初次精细聚类分割；需要注意的是，分块处理不是分的个数越多、范围越小越好；过小的分块处理可能导致子区域内云图类别个数稀少，甚至只有一类，这样聚类分割就失去意义，分块越多更增加二次聚类的复杂度；综上所述，本发明经过多次试验对FY2-D卫星云图IR1通道进行3×3分块处理最佳，即L=3；

②对各个子区域云图进行初次聚类分割：

②-1将子区域云图中每个像素点作为输入样本点，提取每个子区域云图中所有样本点在多个通道云图中相应区域内对应点的邻域灰度均值，该邻域灰度均值作为相应样本点的光谱特征，其中多个通道云图的分辨率一致；需要说明的是：大型气象卫星扫描辐射计一般拥有多个通道，本实施例中，我国FY2-D星扫描辐射计拥有5个通道，分别为3个红外通道、1个水汽通道及1个可见光通道。传统云图分割法利用单个通道光谱特征进行分割，结果精确度小，鲁棒性弱，因此，在对单个通道子区域云图进行分割前，本发明提取的光谱特征为子区域云图每个样本点的多通道光谱特征，即此点在5个通道云图相同区域内该点的邻域灰度均值特征。针对多通道云图分辨率不同的问题，以IR1通道为基准采用双线性插值法使5通道云图分辨率一致。如此，每个样本点就提取到对应于5个通道子区域的5个光谱特征供后续分割；

②-2将每个子区域云图中的某个像素点看做中心生成滑动窗口，将此滑动窗口看做一张子图，计算该滑动窗口内所含的所有像素点的三模块局部二元模式码值，并基于三模块局部二元模式码值构建得到该滑动窗口的三模块局部二元模式码图，统计该滑动窗口的三模块局部二元模式码图的一阶直方图并计算其熵、均值、一致性和三阶距，得到子区域云图中该像素点的二次统计量纹理特征；以上述方法得到每个子区域云图中所有的像素点的二次统计量纹理特征；具体步骤为：

a、选取子区域云图中的任一像素点O，以该像素点为中心生成(2M+1)×(2M+1)的滑动窗口，此时，将此滑动窗口看做一张子图，在子图中任选一个像素点p，以该像素点p为中心生成大小为3×3的邻域模块，记为C_p，同时在以该模块C_p为中心，半径为r圆周上额外生成s个与该模块C_p大小相同且均匀分布的模块，s个模块分别记为C₀，C₁，……C_j……C_s-1，其中4＜r＜M，r＜M＜(F/L)÷10，s≥2，j=0，1，……，s-1，F表示卫星云图的分辨率；

b、在s个模块中选取C_j和C_j′两个模块，其中j′=（j+a）mod（s），a为整数且2≤a≤s-1，mod表示取两者相除后的余数；

c、计算C_p、C_j和C_j‘的局部二元模式码值，分别记为LBP(C_p)、LBP(C_j)和LBP(C_j′)；

d、利用公式φ(h,q)=LBP(h)-LBP(q)计算LBP(C_j)和LBP(C_j′)与LBP(C_p)的差值， ${\mathrm{\&phi;}}_{C_j}(C_j,C_p)=\mathrm{LBP}\left(C_j\right)-\mathrm{LBP}\left(C_p\right),$ ${\mathrm{\&phi;}}_{C_{j^{\&prime;}}}(C_{j^{\&prime;}},C_p)=\mathrm{LBP}\left(C_{j^{\&prime;}}\right)-\mathrm{LBP}\left(C_p\right);$

e、将

和

代入公式 $I_p={\mathrm{\&phi;}}_{C_j}(C_j,C_p)-{\mathrm{\&phi;}}_{C_{j^{\&prime;}}}(C_{j^{\&prime;}},C_p)$ 计算得到一组I_p，利用公式 $T\left(I_p\right)=\left\{\begin{array}{c}1,I_p\&GreaterEqual;0\\ 0,I_p=<0\end{array}\right.$ 对I_p进行二值化，得到一组二进制数串T(I_p)；

f、将T(I_p)代入公式 ${\mathrm{TPLBP}}_{r,s,a}\left(p\right)=\underset{j=0}{\overset{s-1}{\mathrm{\&Sigma;}}}T({\mathrm{\&phi;}}_{C_j}(C_j,C_p)-{\mathrm{\&phi;}}_{C_{j^{\&prime;}}}(C_{j^{\&prime;}},C_p))\&times;2^j$ 进行加权计算得到TPLBP_r，s，a(p)，TPLBP_r，s，a(p)即为像素点p的十进制三模块局部二元模式码值；

g、同理计算(2M+1)×(2M+1)大小滑动窗口内所含的所有像素的十进制三模块局部二元模式码值并生成窗口内三模块局部二元模式码图，统计窗口三模块局部二元模式码图的一阶直方图并计算其熵、均值、一致性和三阶距，得到子区域云图中像素点O的二次统计量纹理特征，依次利用滑动窗口对整个子区域云图进行操作，得到子区域云图中每个像素点的三模块局部二元模式二次统计量纹理特征，至此每个像素点可以提取到2个纹理特征，作为初次聚类分割输入样本点的二维特征；

h、按照步骤a~g利用滑动窗口可对每个子区域云图中进行纹理特征的提取。

②-3采用最大最小距离法自适应确定初始聚类中心及其个数，具体步骤为：

a、将卫星云图表示为图像的灰度直方图{H_l}输入，l为图像的像素灰度级，H_l为图像中像素灰度级为l的像素总个数，l0，1，…，255；

b、选择输入样本中任意一样本点的像素灰度级作为第一个已知的聚类中心，记为l₁；

c、将已知聚类中心l₁代入公式d_ll′=|l-l'|计算灰度级为l的样本点与已知聚类中心的距离，其中l'表示已知聚类中心；

d、在{H_l}中找出离l₁距离最远的样本点的像素灰度级的作为第二个已知的聚类中心，记为l₂；

e、分别计算{H_l}中除已知聚类中心对应的样本点以外的剩余样本点到所有的已知聚类中心的距离；

f、比较同一样本点到各已知聚类中心的距离，选择其中的最小值；

g、比较步骤f中所有样本点到各已知聚类中心的距离最小值，选择所有距离最小值中的最大值d_max，若d_max＞E×d_ll'，其中0.5<E<1，则d_max。对应的样本点的像素灰度级作为新的已知的聚类中心，并返回步骤e；若d_max≤E×d_ll'，则统计所有已知的聚类中心，得到初始聚类中心v_g及其数量z，g=1，2，……，z；

②-4将步骤②-1和步骤②-2中每块子区域云图中提取的每个像素点的两个纹理特征和一个灰度均值特征构成的特征矩阵作为输入样本和初始聚类中心通过非线性映射函数映射到高维特征空间，在高维空间中对卫星云图中所有样本点进行模糊核聚类，得到每块子区域云图的初次聚类分割结果；具体步骤为：

a、将图像样本点和初始聚类中心v_g通过非线性映射函数Φ(x)从输入空间映射到高维特征空间，图像中各样本点用x_k表示，k=1，2，……，n，n表示图像包含的样本点的总数，x_k表示图像中第k个样本点；

b、将x_k和v_g代入公式D=||Φ(x_k)-Φ(v_g)||中，计算得到样本点x_k与初始聚类中心v_g之间的欧式距离D；

c、将||Φ(x_k)-Φ(v_g)||²=Φ(x_k)^TΦ(x_k)-2Φ(x_k)^TΦ(v_g)+Φ(v_g)^TΦ(v_g)代入特征空间模糊核聚类算法目标函数 $J_m(U,V)=\underset{g=1}{\overset{z}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m{\left|\right|\mathrm{\&Phi;}\left(x_k\right)-\mathrm{\&Phi;}\left(v_g\right)\left|\right|}^2$ 中，其中U＝{u_gk}，Ｖ＝(Φ(v₁)，Φ(v₂)，...,Φ(v_g))，u_gk表示第k个样本点对应于第g类的隶属度，常数m＞1为隶属度矩阵指数，令核函数K(x,y)=Φ(x)^TΦ(y)，则 $J_m(U,V)=\underset{g=1}{\overset{z}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m(K_{\mathrm{kk}}-{2K}_{\mathrm{kg}}+K_{\mathrm{gg}}),$ 其中K_kk=K(x_k，x_k)，K_kg=K(x_k，v_g)，K_gg=K(v_g，v_g)，

目前常用的核函数K(x，y)有高斯核函数、多项式核函数或者Sigmoid核函数等，本实施例中选择高斯核函数，即其中e为常数，即自然对数函数的底数，σ为高斯函数尺度，即宽度；

d、在约束条件

下，计算得到特征空间模糊核聚类隶属度u的函数式为 $u_{\mathrm{gk}}=\frac{{(1/(K_{\mathrm{kk}}-{2K}_{\mathrm{kg}}+K_{\mathrm{gg}}))}^{1/(m-1)}}{\underset{g=1}{\overset{z}{\mathrm{\&Sigma;}}}{(1/(K_{\mathrm{kk}}-{2K}_{\mathrm{kg}}+K_{\mathrm{gg}}))}^{1/(m-1)}}$ 及聚类中心v的函数式为 $v_g=\frac{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m{K}_{\mathrm{kg}}{x}_k}{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m{K}_{\mathrm{kg}}};$

e、初始化参数m、U_gk和核函数K(x,y)的尺度σ，计算得到核函数矩阵K，再对特征空间模糊核聚类隶属度u及聚类中心v进行反复迭代更新，求得最终的U_gk及V_g，根据U_gk及V_g判别得到图像中各样本点的隶属类别；本实施例中，结合样本点多通道光谱特征及TPLBP纹理特征，即每个样本点共7维的特征矩阵作为输入，依次对3×3共9个子区域云图进行②-3、②-4步骤的核聚类分割，从而得到9块子区域云图的初次聚类分割结果；

③对各个子区域云图进行二次聚类分割：

③-1将初次聚类结果中每个类作为输入样本点，提取该输入样本点的类内灰度均值特征；在得到初次聚类结果后，以3×3共9块子区域云图分割结果中各类范围内样本点的灰度均值作为二次聚类的特征，这里所说的灰度均值是指初次聚类分割结果中各类对应于原子区域云图相同范围内的样本点灰度均值；

③-2提取初次聚类分割结果的类内多样性指数特征，具体步骤为：

a、计算每个子区域云图中感兴趣区域块内灰度的最小值Gray_min及最大值Gray_max，将区间[Gray_min，Gray_max]等间隔的划分成10个区段，区段号分别为N₁，N₂，……，N₁₀；其中，子区域云图初次聚类结果图中一个类范围块所对应于聚类前原子区域云图相同范围内的云图块为感兴趣区域块；

b、计算每个子区域云图中所有的感兴趣区域内的灰度值Gray，将得到的灰度值Gray分别映射到相应的区段，并用对应的区段号代替原灰度值Gray；

c、统计每个区段中各区段号N_t的数量I_t，t＝1，2，……，10，设每个子区域云图中所有的感兴趣区域的总像素个数为N'，

将P_t代入香农公式

得到各子区域云图的多样性指数特征；

③-3将初次聚类分割结果的类内灰度特征和多样性指数特征作为样本点的二维特征，构成特征矩阵作为二次聚类分割的输入样本，采用最大最小距离法自适应确定最终聚类中心个数，其原理与步骤②-3相同；

③-4将二次聚类分割算法的样本点通过非线性映射函数映射到高维特征空间，在高维空间中对所有样本点进行模糊核聚类，最终得到整张卫星云图的二次聚类分割结果，其原理与步骤②-4相同。

对本发明的方法进行实验仿真，实验仿真在Matlab平台上进行。

1）TPLBP纹理特征特性分析

TPLBP算子是对基本局部二元模式(Local Binary Patterns；LBP)算子的改进，能够很好的描述图像的局部微纹理特征，同时它克服了LBP算子无法进行大尺度范围纹理特征计算的缺陷，并且由于TPLBP考虑到半径为r的圆形邻域模块LBP关系，有效的减少了噪声点的影响，大大提高其鲁棒性。以Brodatz自然纹理图像为例(如图3a)，生成其TPLBP码图(如图3c)，并对原图加入密度为0.02的椒盐噪声(如图3b)，再次生成其TPLBP码图(如图3d)。由图3c可以看到，TPLBP码图能很好的反映图像的纹理信息，并且在加入椒盐噪声后其TPLBP码图图3d相比于图3c变化甚微，这就说明了TPLBP算子抗噪声强的优势。

回归到云图，仿真以晴空陆地、晴空海面、密层云、卷层云、积雨云、高层或雨层云、层积或高积云7种云图类别建立测试样本(如图4a~4g)，利用TPLBP算子对7类样本生成TPLBP码图，并统计其一阶直方图，进而计算其熵、均值、一致性、三阶距四个二次统计量，最后绘制其纹理特征曲线如图5所示。由图5四条特征曲线我们可以看到，熵和一致性曲线波动小，一致性甚至在前6类几乎呈直线，说明这两类特征对于7类样本区分度低。相反，均值与三阶距特性曲线波动大，各类之间特征值差异大，能够很好的区分7类云图样本。因此本实施例在子区域云图进行初次聚类时，每个样本点只提取均值、三阶距两个TPLBP纹理特征。

2）DI特征特性分析

DI特征相对于其他云图特征能更好的保留云的细节信息。常用特征计算法或多或少存在均值化处理而弱化目标图像的细节信息，而DI特征基于香农公式反映的是地区生物种类的复杂度，也即感兴趣区域块内云类的复杂度，因此避免了这一问题。同样以图4a~4g七类云图样本作为测试数据提取其各自DI特征，并绘制DI特征曲线如图6所示。图6中可以看到，7类云图样本各自DI特征值差异明显，也即能很好的区分7类样本，但是4与6类(卷层云与高层或雨层云)DI特征值差异小，这与第6类样本在手工选取时人为去除右下角干扰云而导致像素丢失有关。DI特征曲线总体上表现出良好的区分7类云图样本的特性。

3）FY2-D云图云分类结果

利用中国国家卫星气象中心提供的FY2-D卫星采集于2012年5月29日凌晨5点30分IR1通道1900×1300×8bit分区图，辅之以IR2、IR3、水汽、VIS四通道分区图进行实验。根据上述算法对IR1通道云图进行二次聚类分割。

实验初次聚类结果如图7所示，二次聚类分割最终云分类结果如图8所示。图9展示了利用整幅IR1通道云图直接进行聚类分割处理结果，利用中国国家卫星气象中心提供的同一时刻云图专家分类结果图(如图10)作为基准进行对比，验证本发明提出的二次聚类分割算法在FY2-D云图分类的准确性、优越性。由图对比可以看到，直接利用整幅图进行聚类分割结果左下方、左上方、中下方及右上方没能将晴空陆地与海洋正确区分开，左下方和左下方的层积云、高层云也没能识别出,同时在整体上云识别出来的范围普遍偏小以至于占比例小的卷层云几乎没有被识别出来。而本发明算法结果在此些方面就显示出其优越性，图中结果克服了基于整张云图直接聚类分割误分、错分的缺点，它能正确的区分晴空陆地与海面，各类云之间也能正确精细的区分出，在整体云分类结果上与卫星气象中心提供专家分类结果基本吻合。但是，在某些细节上分割还有待完善，比如本文算法结果图8相比于专家分割结果，在右上方一小块晴空陆地与海面被误分，同时在中国版图内卷层云、积雨云没能正确的识别出来，当然这与卷层云本身所占比例小而难于区分有关。

表1对比了采用本发明方法的分割结果与基于整张云图直接聚类分割结果的错分率(错分率定义为结果中每类样本错分像素点个数占图像总像素个数的百分比)。从错分率指标上来看，本发明方法云分类结果正确率相比于直接聚类分割云分类而言，优势相当明显。这主要是因为本发明方法在初次聚类分割时仅仅是针对小范围云图块进行精细分割，一定程度上避免了大范围云图分割带来的错分误分。再者，全局范围的第二次聚类分割在保证分割的完整性的同时又实现了对初次聚类分割结果的适当修正，使同种云最大程度归为同一类。

表1云分类错误率计算(%)

Table1(%)

Claims (6)

1.一种卫星云图的二次聚类分割方法，其特征在于包括以下步骤：

①对卫星云图进行分块处理，将卫星云图分为L×L共L²块大小相等的子区域云图，其中2≤L≤10；

②对各个子区域云图进行初次聚类分割：

②-1将子区域云图中每个像素点作为输入样本点，提取每个子区域云图中所有样本点在多个通道云图中相应区域内对应点的邻域灰度均值，该邻域灰度均值作为相应样本点的光谱特征，其中多个通道云图的分辨率一致；

②-2将每个子区域云图中的某个像素点看做中心生成滑动窗口，将此滑动窗口看做一张子图，计算该滑动窗口内所含的所有像素点的三模块局部二元模式码值，并基于三模块局部二元模式码值构建得到该滑动窗口的三模块局部二元模式码图，统计该滑动窗口的三模块局部二元模式码图的一阶直方图并计算其熵、均值、一致性和三阶距，得到子区域云图中该像素点的二次统计量纹理特征；以上述方法得到每个子区域云图中所有的像素点的二次统计量纹理特征；

②-3采用最大最小距离法自适应确定初始聚类中心及其个数；

②-4将步骤②-1和步骤②-2中每块子区域云图中提取的每个像素点的两个纹理特征和一个灰度均值特征构成的特征矩阵作为输入样本和初始聚类中心通过非线性映射函数映射到高维特征空间，在高维空间中对卫星云图中所有样本点进行模糊核聚类，得到每块子区域云图的初次聚类分割结果；

③对各个子区域云图进行二次聚类分割：

③-1将初次聚类结果中每个类作为输入样本点，提取该输入样本点的类内灰度均值特征；

③-2提取初次聚类分割结果的类内多样性指数特征；

③-3将初次聚类分割结果的类内灰度特征和多样性指数特征作为样本点的二维特征，构成特征矩阵作为二次聚类分割的输入样本，采用最大最小距离法自适应确定最终聚类中心个数；

③-4将二次聚类分割算法的样本点通过非线性映射函数映射到高维特征空间，在高维空间中对所有样本点进行模糊核聚类，最终得到整张卫星云图的二次聚类分割结果。

2.根据权利要求1所述的一种卫星云图的二次聚类分割方法，其特征在于步骤②-2中得到每个子区域云图的二次统计量纹理特征的具体步骤为：

a、选取子区域云图中的任一像素点O，以该像素点为中心生成(2M+1)×(2M+1)的滑动窗口，此时，将此滑动窗口看做一张子图，在子图中任选一个像素点p，以该像素点p为中心生成大小为3×3的邻域模块，记为C_p，同时在以该模块C_p为中心，半径为r圆周上额外生成s个与该模块C_p大小相同且均匀分布的模块，s个模块分别记为C₀，C₁，……C_j……C_s-1，其中4＜r＜M，r＜M＜(F/L)÷10，s≥2，j=0，1，……，s-1，F表示卫星云图的分辨率；

b、在s个模块中选取C_j和C_j′两个模块，其中j′=（j+a）mod（s），a为整数且2≤a≤s-1，mod表示取两者相除后的余数；

c、计算C_p、C_j和C_j‘的局部二元模式码值，分别记为LBP(C_p)、LBP(C_j)和LBP(C_j′)；

d、利用公式φ(h,q)=LBP(h)-LBP(q)计算LBP(C_j)和LBP(C_j′)与LBP(C_p)的差值， ${\mathrm{\&phi;}}_{C_j}(C_j,C_p)=\mathrm{LBP}\left(C_j\right)-\mathrm{LBP}\left(C_p\right),$ ${\mathrm{\&phi;}}_{C_{j^{\&prime;}}}(C_{j^{\&prime;}},C_p)=\mathrm{LBP}\left(C_{j^{\&prime;}}\right)-\mathrm{LBP}\left(C_p\right);$

e、将

和

代入公式 $I_p={\mathrm{\&phi;}}_{C_j}(C_j,C_p)-{\mathrm{\&phi;}}_{C_{j^{\&prime;}}}(C_{j^{\&prime;}},C_p)$ 计算得到一组I_p，利用公式 $T\left(I_p\right)=\left\{\begin{array}{c}1,I_p\&GreaterEqual;0\\ 0,I_p=<0\end{array}\right.$ 对I_p进行二值化，得到一组二进制数串T(I_p)；

f、将T(I_p)代入公式 ${\mathrm{TPLBP}}_{r,s,a}\left(p\right)=\underset{j=0}{\overset{s-1}{\mathrm{\&Sigma;}}}T({\mathrm{\&phi;}}_{C_j}(C_j,C_p)-{\mathrm{\&phi;}}_{C_{j^{\&prime;}}}(C_{j^{\&prime;}},C_p))\&times;2^j$ 进行加权计算得到TPLBP_r，s，a(p)，TPLBP_r，s，a(p)即为像素点p的十进制三模块局部二元模式码值；

g、同理计算(2M+1)×(2M+1)大小滑动窗口内所含的所有像素的十进制三模块局部二元模式码值并生成窗口内三模块局部二元模式码图，统计窗口三模块局部二元模式码图的一阶直方图并计算其熵、均值、一致性和三阶距，得到子区域云图中像素点O的二次统计量纹理特征，依次利用滑动窗口对整个子区域云图进行操作，得到子区域云图中每个像素点的三模块局部二元模式二次统计量纹理特征，至此每个像素点可以提取到2个纹理特征，作为初次聚类分割输入样本点的二维特征；

h、按照步骤a~g利用滑动窗口可对每个子区域云图中进行纹理特征的提取。

3.根据权利要求2所述的一种卫星云图的二次聚类分割方法，其特征在于步骤②-3中采用最大最小距离法自适应确定聚类中心个数的具体步骤为：

a、将卫星云图表示为图像的灰度直方图{H_l}输入，l为图像的像素灰度级，H_l为图像中像素灰度级为l的像素总个数，l0，1，…，255；

b、选择输入样本中任意一样本点的像素灰度级作为第一个已知的聚类中心，记为l₁；

c、将已知聚类中心l₁代入公式d_ll'=|l-l'|计算灰度级为l的样本点与已知聚类中心的距离，其中l'表示已知聚类中心；

d、在{H_l}中找出离l₁距离最远的样本点的像素灰度级的作为第二个已知的聚类中心，记为l₂；

e、分别计算{H_l}中除已知聚类中心对应的样本点以外的剩余样本点到所有的已知聚类中心的距离；

f、比较同一样本点到各已知聚类中心的距离，选择其中的最小值；

g、比较步骤f中所有样本点到各已知聚类中心的距离最小值，选择所有距离最小值中的最大值d_max，若d_max。＞E×d_ll'，其中0.5<E<1，则d_max。对应的样本点的像素灰度级作为新的已知的聚类中心，并返回步骤e；若d_max≤E×d_ll′，则统计所有已知的聚类中心，得到初始聚类中心v_g及其数量z，g=1，2，……，z。

4.根据权利要求3所述的一种卫星云图的二次聚类分割方法，其特征在于步骤②-4中获得卫星云图初次聚类结果的具体步骤为：

a、将图像样本点和初始聚类中心v_g通过非线性映射函数Φ(x)从输入空间映射到高维特征空间，图像中各样本点用x_k表示，k=1，2，……，n，n表示图像包含的样本点的总数，x_k表示图像中第k个样本点；

b、将x_k和v_g代入公式D=||Φ(x_k)-Φ(v_g)||中，计算得到样本点x_k与初始聚类中心v_g之间的欧式距离D；

c、将||Φ(x_k)-Φ(v_g)||²=Φ(x_k)^TΦ(x_k)-2Φ(x_k)^TΦ(v_g)+Φ(v_g)^TΦ(v_g)代入特征空间模糊核聚类算法目标函数 $J_m(U,V)=\underset{g=1}{\overset{z}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m{\left|\right|\mathrm{\&Phi;}\left(x_k\right)-\mathrm{\&Phi;}\left(v_g\right)\left|\right|}^2$ 中，其中U＝{u_gk}，V＝(Φ(v₁)，Φ(v₂)，...,Φ(v_g))，u_gk表示第k个样本点对应于第g类的隶属度，常数m＞1为隶属度矩阵指数，令核函数K(x,y)=Φ(x)^TΦ(y)，则 $J_m(U,V)=\underset{g=1}{\overset{z}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m(K_{\mathrm{kk}}-{2K}_{\mathrm{kg}}+K_{\mathrm{gg}}),$ 其中K_kk=K(x_k，x_k)，K_kg=K(x_k，v_g)，K_gg=K(v_g，v_g)，

d、在约束条件

下，计算得到特征空间模糊核聚类隶属度u的函数式为 $u_{\mathrm{gk}}=\frac{{(1/(K_{\mathrm{kk}}-{2K}_{\mathrm{kg}}+K_{\mathrm{gg}}))}^{1/(m-1)}}{\underset{g=1}{\overset{z}{\mathrm{\&Sigma;}}}{(1/(K_{\mathrm{kk}}-{2K}_{\mathrm{kg}}+K_{\mathrm{gg}}))}^{1/(m-1)}}$ 及聚类中心v的函数式为 $v_g=\frac{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m{K}_{\mathrm{kg}}{x}_k}{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{gk}}^m{K}_{\mathrm{kg}}};$

e、初始化参数m、U_gk和核函数K(x,y)，计算得到核函数矩阵K，再对特征空间模糊核聚类隶属度u及聚类中心v进行反复迭代更新，求得最终的U_gk及V_g，根据U_gk及V_g判别得到图像中各样本点的隶属类别。

5.根据权利要求4所述的一种卫星云图的二次聚类分割方法，其特征在于所述的核函数为高斯核函数、多项式核函数或者Sigmoid核函数。

6.根据权利要求4或5所述的一种卫星云图的二次聚类分割方法，其特征在于步骤③-2中提取多样性指数特征的具体步骤为：

a、计算每个子区域云图中感兴趣区域块内灰度的最小值Gray_min及最大值Gray_max，将区间[Gray_min，Gray_max。]等间隔的划分成10个区段，区段号分别为N₁，N₂，……，N₁₀，其中，子区域云图初次聚类结果图中一个类范围块所对应于聚类前原子区域云图相同范围内的云图块为感兴趣区域块；

b、计算每个子区域云图中所有的感兴趣区域内的灰度值Gray，将得到的灰度值Gray分别映射到相应的区段，并用对应的区段号代替原灰度值Gray；

c、统计每个区段中各区段号N_t的数量I_t，t＝1，2，……，10，设每个子区域云图中所有的感兴趣区域的总像素个数为N'，

将P_t代入公式

得到各子区域云图的多样性指数特征。

Images