利用对象光谱与纹理实现高分辨率遥感影像云检测方法

董志鹏，王 密，2，李德仁，2，王艳丽，张致齐

1. 武汉大学测绘遥感信息工程国家重点实验室，湖北 武汉 430079； 2. 地球空间信息协同创新中心，湖北 武汉 430079

随着对地观测卫星技术的发展，高分辨率的遥感影像已经应用到城市规划、作物分类、灾害监测等众多领域[1-2]。由于云层覆盖了地表上空约66%的面积，造成高分辨率遥感影像中有大量云区域存在，降低了高分辨率遥感影像的使用价值，对基于高分辨率遥感影像的目标识别、图像分类等处理产生负面影响[3-5]。因此，云检测已成为高分辨率遥感影像处理中非常重要的内容[6-7]。针对遥感影像的云检测处理，国内外学者开展了大量的研究，由于高分辨率遥感影像光谱波段范围的限制，其云检测算法多采用影像多特征融合处理的思想实现高分辨率遥感影像中云层区域的检测[8-9]。如文献[9]首先基于光谱特征得到影像中云区域的粗提取结果，其次结合几何特征与纹理特征对提取结果进行提纯，最后对提纯结果进行连通性处理获得最终的影像云检测结果，但该方法的云检测结果易受影像中山脉、雪等类云地物影响。文献[10]根据影像中云层对象和下垫面对象的光谱、纹理等属性训练SVM分类器，使用训练好的SVM分类器对影像中图像块对象进行是否为云的判断，提取影像中的云区域，但该方法需要大量的云检测训练样本。文献[11]首先将影像分割为64×64大小的图块对象，根据每个对象的光谱均值与方差、纹理角二阶矩与一阶差分等对对象是否为云进行判断，实现影像中云区域的检测，但该方法中合适的影像云检测光谱阈值难以确定，云检测精度较低。文献[12]将影像从RGB转换到HIS颜色空间，利用光谱特征实现影像中云区域的粗提取，利用纹理属性对粗提取结果进行提纯处理，最后对云层提取结果膨胀处理，获得最终的云检测结果，但该方法中合适的影像云检测光谱阈值难以确定。综上所述，在基于影像多特征融合处理的高分辨率遥感影像云检测思想中，存在合适的影像云检测光谱阈值难以确定的问题，以及影像中雪、山脉、建筑物等类云地物影响影像云检测精度等问题。

针对以上问题，本文提出一种基于对象光谱与纹理的高分辨率遥感影像云检测方法。该方法首先对影像进行直方图均衡化处理，根据均衡化影像直方图获得合适的云检测光谱阈值。其次用简单线性迭代聚类(simple linear iterative clustering，SLIC)算法对影像进行分割，生成超像素对象。根据云检测光谱阈值，以超像素对象为处理单元获得初始影像云检测结果。然后求得影像旋转不变局部二值模式(local binary patterns，LBP)纹理图，根据超像素的LBP纹理均值与角二阶矩对初始影像云检测结果提纯，消除类云地物对云检测的影响。最后对提纯后的影像云区域进行区域增长及膨胀处理，获得最终的影像云检测结果。定性对比试验和定量评价验证了本文方法的有效性。

1 方 法

本文方法主要分为4个步骤：①影像直方图均衡化处理，根据均衡化影像直方图获得合适的云检测光谱阈值；②SLIC算法对影像分割生成超像素对象，根据云检测光谱阈值和超像素光谱属性获得初始影像云检测结果；③求得均衡化影像的LBP纹理图，根据超像素LBP纹理均值与角二阶矩对初始影像云检测结果提纯；④对提纯后云区域进行区域增长与膨胀处理，获得最终的云检测结果。总体流程如图1所示。

图1 影像云检测流程Fig.1 The flowchart of cloud detection

1.1 云区域粗提取

1.1.1 确定自适应云检测光谱阈值

云在可见光和近红外波段对于光线的反射率比大多数地物强，在影像上表现为云相对于地面目标有较高的灰度值[13-14]。因此利用该特征，采用基于光谱特征的阈值判断可以有效地实现影像中的云与地面目标的分类[6]。而合适的影像云检测光谱阈值多采用人为试错求得，难以准确自动获得。高分辨率遥感影像云区域的云边界到云中心存在着薄云到厚云的过渡带区域的特征，如图2 所示。通过求得云过渡带区域像素对应的光谱值，利用云过渡带区域像素光谱值求得云边界像素光谱值，将大于云边界像素光谱值的影像区域作为云区域，可以有效提取影像中的云区域。

图2 试验影像Fig.2 Experimental image

本文使用的试验数据为包含蓝(0.45～0.52 μm)、绿(0.52～0.59 μm)、红(0.63～0.69 μm)和近红外(0.77～0.89 μm)4个波段的10 bits多光谱遥感影像，影像的光谱值范围为0～1023。选取影像的近红外、红和绿波段作为试验中图像处理的R、G、B波段，并将影像的光谱属性除以4，将影像的光谱值范围压缩为0～255，便于后续的试验计算处理。根据式(1)求得图2影像的灰度直方图，图3(a)为该影像灰度直方图走势图。将影像的R、G、B 3个波段分别进行直方图均衡化处理，均衡化的目的是突出隐含有纹理细节的图像，由于地物包含较丰富的纹理信息，其细节清晰度比均衡化之前有较大提高[11]。根据式(1)求得影像均衡化后的灰度直方图[15]，图3(b)为该影像均衡化后灰度直方图走势图。

gray=R×0.299+G×0.587+B×0.114

(1)

图3 影像灰度直方图走势Fig.3 The image gray histogram trend diagram

在图3(a)中看似不存在明显的可获得云检测光谱阈值的变化规律。在图3(b)中当灰度为219时直方图剧烈下降，在灰度为224时直方图上升，灰度大小属于[219，224]的影像云掩膜结果如图4(a)所示。图4(b)为图4(a)影像云掩膜结果与原始影像叠加示意图，当灰度大小处于[219，224]时，对应的影像云掩膜结果处于影像云边界与云中心的过渡区域。试验结果表明，通过捕获影像均衡化直方图中的突变点，可以准确获得云过渡带区域像素对应的光谱值。利用云过渡带区域像素光谱值求得云边界像素光谱值，将大于云边界像素光谱值的区域作为云，实现影像中云区域的提取。统计157景资源三号02星、高分一号和高分二号多光谱影像通过均衡化直方图求影像

云检测光谱阈值试验结果，总结出以下求自适应云检测光谱阈值准则:

(1) 在[160，253]内TD1为均衡化直方图中的极值点，其对应的像素数为count[TD1]；TD2为大于TD1的第一个极值点，其对应的像素数为count[TD2]。当TD1、TD2为[160，253]内首次满足式(2)的灰度值时，[TD1,TD2]为该影像云过渡区域像素对应的光谱值范围。

(2)

式中，max(count[TD1]、count[TD2])为count[TD1]与count[TD2]中的较大值。

(2) 影像云边界像素的光谱值TD=TD1-2，则影像自适应云检测光谱阈值为TD。

图4 试验结果Fig.4 Experimental results

1.1.2 获取初始云检测结果

由于高分辨率遥感影像中存在“椒盐”噪声，基于像素的光谱阈值云检测会将光谱值高的“椒盐”噪声检测为云，而基于对象的高分辨率遥感影像处理方法可以有效地消除“椒盐”噪声对云检测的影响[16-17]。文献[18]提出了超像素这一概念，所谓超像素是指具有相似纹理、颜色、亮度等特征的相邻像素构成的图像块。文献[19]中提出SLIC算法生成超像素，将彩色图像转换为CIELAB颜色空间和X、Y坐标下的5维特征向量，然后对5维特征向量构造度量标准，对图像像素进行局部聚类生成超像素。相对于分水岭算法、区域增长算法和基于图的图像分割算法等传统的影像分割算法生成的超像素，SLIC算法生成的超像素具有更好的地物边界依附性、更加规则紧凑的形状，并且可以人为控制生成超像素的个数并具有更好的抗噪性[15,20]。使用SLIC算法对影像进行分割生成超像素对象，以超像素对象为处理单元，利用式(1)求得超像素的灰度值，将灰度值大于等于云检测光谱阈值TD的超像素作为云，可以有效地消除“椒盐”噪声对影像云检测结果的影响，并获得初始的影像云检测结果。

经试验统计分析得到当超像素包含的像素数为500时，影像的云检测效果最佳，则影像分割生成超像素个数为N/500，其中N为影像中像素数。统计每个超像素包含的像素的光谱均值作为超像素的光谱属性值，则超像素的光谱属性为Si(RiGiBi)，计算公式如下

(3)

式中，ni为第i个超像素包含的像素的个数；Rk、Gk、Bk为第i个超像素中包含的第k个像素的R、G、B属性值。

1.2 粗差剔除

沙漠、山脉等地物在可见光、近红外波段具有较高的反射率，在高分辨率遥感影像上表现为具有较高的灰度值，如图5(a)、(c)。基于对象光谱阈值的云检测方法会将光谱属性值高的沙漠、山脉等地物识别为云，则需要消除初始云检测结果中的沙漠、山脉等类云地物。文献[21]提出采用局部二值模式用于图像的纹理分析，即求得图像的LBP纹理。LBP通过局部模式分析纹理，打破了传统的在像素层上纹理分析的方式，且相比其他局部特征描述子(如SIFT、RIFF等)，LBP具有理论简单、容易计算和特征提出无须训练的优势[22]。为了使图像具有旋转不变性，文献[23]对原始LBP进行改进，提出旋转不变的LBP，通过旋转圆形邻域得到一系列的初始定义的LBP，取其最小值作为旋转不变的LBP的值。大量旋转不变和光照不变的纹理试验表明，旋转不变LBP具有良好的低维性和不变性，可以很好地用于图像分类[22-23]。本文使用旋转不变的LBP纹理用于沙漠、山脉等类云地物与云之间的差异区分。图5(b)、(d)分别为图5(a)、(c)影像直方图均衡化后的旋转不变LBP纹理图。在图5(a)、(c)中沙漠、山脉与云在影像上均表现出具有较高的光谱属性，但在图5(b)、(d)中沙漠、山脉等类云地物表现出灰度值较低，且灰度分布不均匀及纹理较细的特征，而云具有较高的灰度值，且灰度分布均匀及纹理较粗的特征。据此，通过LBP纹理图中的灰度值大小与分布的均匀程度可实现云与类云地物间的分类。

本文以超像素作为处理单元，求得各超像素包含像素的LBP值均值，用来描述超像素在LBP图像中灰度值的大小，计算公式如式(4)所示。角二阶矩能良好地反映图像灰度分布均匀程度和纹理粗细度[24]，则求得超像素包含像素LBP值的灰度共生矩阵的角二阶矩，用来描述超像素包含像素的LBP值分布的均匀性，计算公式如式(5)所示。其中云的纹理较粗，角二阶矩较大；类云地物的纹理较细，角二阶矩较小[25-26]。影像中所有超像素LBP值的均值计算公式如式(6)所示，所有超像素角二阶矩的均值计算公式如式(7)所示。通过大量试验得出，根据超像素的LBP值与角二阶矩实现云与类云地物间分类的判断准则如式(8)所示。在初始云检测结果中，当超像素的LBP值与角二阶矩同时满足式(8)时，该超像素为云，否则为非云地物，经过此判断消除初始云检测结果中的类云地物。

图5 试验影像Fig.5 Experimental images

(4)

(5)

(6)

(7)

(8)

式中，SLBPi为第i个超像素的LBP值；ni为第i个超像素包含的像素的个数；LBPk为第i个超像素中包含的第k个像素的LBP值；SASMi为第i个超像素的角二阶矩；P(l,m)为第i个超像素中像素的LBP值归一化灰度共生矩阵中位置(l,m)处值；aveLBP为所有超像素LBP值的均值；aveASM为所有超像素角二阶矩的均值；n为影像中超像素的个数。

1.3 云区域精检测

在影像旋转不变LBP纹理图像中，云区域边界像素的LBP值比中心像素低，云区域经过粗差剔除后，云区域薄云边界也被剔除，需要恢复云区域的薄云边界。以粗差剔除后云区域中包含的超像素为种子点进行区域增长，判断作为种子点的超像素的邻接超像素光谱属性值是否大于等于初始云检测光谱阈值TD，如果大于等于则将该邻接超像素加入云区域且作为云区域增长的种子点。循环区域增长过程，直到没有超像素加入云区域时停止，从而恢复云区域的薄云边界。

云区域增长结束后，云区域中存在少量间隙孔洞，则对云区域进行膨胀处理，消除其中的孔洞[11]。以云区域中的单个像素为种子点，判断其8邻域的像素是否均为云，将非云像素加入云区域并作为云区域膨胀处理的种子点。试验证明循环膨胀处理5次可以有效地消除云区域中的间隙孔洞。经区域增长与膨胀处理后得到最终的云区域检测结果。

1.4 云检测结果评价

云检测结果评价是影像云检测研究中必不可少的一步，目前多采用目视判别与定量评价相结合的方式验证云检测结果的有效性[3-9,11-12]。本文同样采用目视判别与定量评价相结合的方式对本文方法的有效性进行验证。

目视判别是一种最基本、常用的评价方法。通过目视判别可以直观地观察云检测结果中的漏检、错检等情况，而且只有云检测结果和目视评价的效果相吻合时，定量评价指标才能使人信服。因此，在目视判别的基础上采用召回率(recall)、虚警率(falsealarm)和准确率(accuracy)对本文方法进行定量评价。召回率的取值范围为[0，1]，召回率越大，说明算法识别为云的像素占影像中真云像素的比例越高。虚警率的取值范围为[0，1]，虚警率越小，说明在影像中算法识别为云的像素中非云像素的比例越低。准确率的取值范围为[0，1]，准确率越大，说明算法的云识别能力越高[5,7]，当准确率为1时，说明算法的云检测结果与实际云层分布完全一致。召回率、虚警率和准确率的计算公式如下

(9)

(10)

(11)

式中，TP为算法识别为云的像素中真实云像素的个数；FN为算法将影像中云像素识别为非云像素个数；FP为算法将影像中非云像素识别为云像素的个数；TN为算法识别为非云像素中真实非云像素的个数；N为影像中像素的个数。

为了更加全面地评价本文方法的有效性，与经典的多光谱阈值法[27]与文献[11]中的树状结构云检测方法进行对比评价。

2 试验结果与分析

通过两组试验对本文方法的有效性进行验证。对多光谱阈值法和树状结构云检测方法调整参数进行多次云检测试验，且采用目视判别的方法选择整体云检测效果最佳的结果与本文方法进行定性与定量对比评价。

2.1 试验1

试验1使用的数据为4景资源三号02星标准景多光谱影像，影像包含蓝(0.45～0.52 μm)、绿(0.52～0.59 μm)、红(0.63～0.69 μm)和近红外(0.77～0.89 μm)4个波段，各波段分辨率为5.8 m，大小为8816×8792像素，具体如图6所示。图6(a)影像中心地理坐标为80.3°E，36.2°N，影像中主要包括山脉、云和雪等。图6(b)影像中心地理坐标为82.7°E，44.1°N，影像中主要包括山脉、云和雪等。图6(c)影像中心地理坐标为58.1°E，42.5°N，影像中主要包括戈壁、湖泊和云等。图6(d)影像中心地理坐标为88.5°E，30.3°N，影像中主要包括山脉、雪、冰和云等。

2.1.1 检测结果的目视评价

图7为不同检测方法对图6中4景影像处理后得到的云检测结果。由遥感领域专家目视判别获得4景影像的云检测目视解译结果，如图7(a1)、(b1)、(c1)、(d1)所示。目视解译结果作为云检测算法效能评估的基准。图7(a2)、(b2)、(c2)、(d2)分别为多光谱阈值法对图6中4景影像的云检测试验结果。图7(a3)、(b3)、(c3)、(d3)分别为树状结构云检测方法对图6中4景影像的云检测试验结果。图7(a4)、(b4)、(c4)、(d4)分别为本文方法对图6中4景影像的云检测试验结果。图中的绿色区域为云检测结果。

在图7(a1)、(a2)、(a3)、(a4)云检测试验结果中，黄色矩形框圈定的山脉区域，由于山脉光谱与云光谱相似，多光谱阈值法将山脉区域识别为云，难以实现山脉与云间的区分。在树状结构云检测方法试验结果中，出现了大量的云区域错检现象，云检测识别精度较低。本文方法可以有效地识别该影像中的云区域，云检测结果与目视解译结果相似。

在图7(b1)、(b2)、(b3)、(b4)云检测试验结果中，图7(b3)中出现了云区域错检现象，树状结构云检测方法的精度较低。多光谱阈值法与本文方法均出现了云区域漏检现象，但本文方法的云检测结果与目视解译结果更相似。

在图7(c1)、(c2)、(c3)、(c4)云检测试验结果中，黄色圆形圈定的区域内，由于戈壁区域具有较高的反射率，多光谱阈值法将该区域误检为云，而本文方法可以有效地区分戈壁区域与云层间的差异。在图7(c3)中树状结构云检测方法出现了大量云错检与漏检现象。从整体上观察，本文方法的云检测结果与目视解译结果基本相似。

在图7(d1)、(d2)、(d3)、(d4)云检测试验结果中，黄色矩形框圈定的区域内，由于冰、雪具有较高的反射率，多光谱阈值法将冰、雪误检为云，而本文方法可有效地区分冰、雪与云间的差异。在多光谱阈值法的检测结果中，存在云误检现象。在树状结构方法的检测结果中，存在大量误检与漏检现象，检测精度较低。本文方法的云检测结果与目视解译结果基本一致，获得准确的云检测结果。

2.1.2 检测结果的定量评价

表1列出了3种不同算法对试验1中4景影像云检测结果的召回率、虚警率和准确率的定量评价结果。在表1中本文方法对于4景影像云检测结果的召回率均高于多光谱阈值法和树状结构云检测方法，说明本文方法比多光谱阈值法和树状结构云检测方法有更高的影像云层识别能力。在图6(a)、(c)、(d)影像云检测结果中，本文方法的虚警率均小于多光谱阈值法和树状结构云检测方法；在图6(b)影像的云检测结果中，本文方法与多光谱阈值法的虚警率基本一致，且均小于树状结构云检测方法，说明本文方法的云识别能力的错误率更低，精度更高。在4景影像的云检测结果中，本文方法的准确率均高于多光谱阈值法和树状结构云检测方法，说明本文方法的云检测结果与目视解译结果更相似。表1中本文方法召回率与准确率的平均值分别为0.827与0.969 4，均高于多光谱阈值法的0.736 2与0.945 9和树状结构云检测方法的0.485 9与0.826 5；本文方法虚警率的平均值为0.034 1，小于多光谱阈值法的0.121 8和树状结构法的0.497 3，表明本文方法对资源三号02星影像有更强的云检测能力和更高云识别精度。

2.2 试验2

为了进一步验证本文方法的有效性，采用两景高分一号标准景多光谱影像和两景高分二号标准景多光谱影像对本文方法进行试验验证。图8(a)、(b)为高分一号多光谱影像，影像包含蓝(0.45～0.52 μm)、绿(0.52～0.59 μm)、红(0.63～0.69 μm)和近红外(0.77～0.89 μm)4个波段，各波段分辨率为8 m，大小为4548×4296像素。图8(a)影像中心地理坐标为96.2°E，46.3°N，影像中主要包括戈壁、云、房屋和道路等。图8(b)影像中心地理坐标为108.8°E，47.2°N，影像中主要包括山脉、云、雪和戈壁等。图8(c)、(d)为高分二号多光谱影像，影像包含蓝(0.45～0.52 μm)、绿(0.52～0.59 μm)、红(0.63～0.69 μm)和近红外(0.77～0.89 μm)4个波段，各波段分辨率为4 m，大小为7300×6908像素。图8(c)影像中心地理坐标为102.4°E，36.1°N，影像中主要包括山、海洋、云、房屋和道路等。图8(d)影像中心地理坐标为115.4°E，41°N，影像中主要包括山脉、云、道路和房屋等。

表1 试验1云检测定量评价结果

2.2.1 检测结果的目视评价

图9为不同检测方法对图8中4景影像处理得到的云检测结果。图9(a1)、(b1)、(c1)、(d1)分别为图8中4景影像的云检测目视解译结果。图9(a2)、(b2)、(c2)、(d2)分别为多光谱阈值法对图8中4景影像的云检测试验结果。图9(a3)、(b3)、(c3)、(d3)分别为树状结构云检测方法对图8中4景影像的云检测试验结果。图9(a4)、(b4)、(c4)、(d4)分别为本文方法对图8中4景影像的云检测试验结果。在试验结果中，影像中的绿色区域为云检测结果。

在图9(a1)、(a2)、(a3)、(a4)云检测试验结果中，黄色圆形圈定的区域内，多光谱阈值法将山脉区域识别为云，出现错检现象。在图9(a3)中树状结构云检测方法出现大量云错检与漏检现象。本文方法的云检测结果与目视解译结果基本相似。

在图9(b1)、(b2)、(b3)、(b4)云检测试验结果中，红色圆形圈定的区域内，多光谱阈值法将雪识别为云，出现错检现象。在图9(b3)中树状结构云检测方法出现大量云错检与漏检现象。本文方法也有少量的漏检现象，但是基本与目视解译结果相似。

在图9(c1)、(c2)、(c3)、(c4)云检测试验结果中，黄色矩形圈定的区域内，多光谱阈值法将光谱值较高的房屋识别为云，出现错检现象。在图9(c3)中树状结构云检测方法出现大量云错检现象。本文方法也有少量的漏检现象，但整体上与目视解译结果更接近。

在图9(d1)、(d2)、(d3)、(d4)云检测试验结果中，树状结构云检测方法出现大量云错检与漏检现象，本文方法与多光谱阈值法均出现漏检现象，但是本文方法与目视解译结果更接近。

2.2.2 检测结果的定量评价

表2列出了3种不同算法对试验2中4景影像云检测结果的召回率、虚警率和准确率的定量评价结果。表2中本文方法的召回率与准确率均高于多光谱阈值法和树状结构云检测方法，说明本文方法有更强的影像云识别能力。本文方法的虚警率与多光谱阈值法基本一致，且小于树状结构云检测方法，说明本文方法有更高的云检测精度。表2中本文方法召回率与准确率的平均值分别为0.8和0.965 2，均高于多光谱阈值法和树状结果云检测方法；本文方法虚警率的平均值为0.052 2，小于多光谱阈值法和树状结果云检测方法，说明本文方法对高分一号、高分二号等高分系列影像有更准确的云检测结果。

3 结 论

针对高分辨率遥感影像云检测过程中合适的云检测光谱阈值难以确定，及影像中类云地物对云检测精度影响的问题，本文提出一种基于对象光谱与纹理的高分辨率遥感影像云检测方法。

图6 试验1影像数据Fig.6 Experimental images of test 1

图7 试验1的结果Fig.7 Experimental results of test 1

图8 试验2影像数据Fig.8 Experimental images of test 2

图9 试验2的结果Fig.9 Experimental results of test 2

表2 试验2云检测定量评价结果

该方法首先对影像进行直方图均衡化处理，根据均衡化影像直方图获得合适的影像云检测光谱阈值。其次使用SLIC算法对影像进行分割生成超像素对象，以超像素为处理单元，根据云检测光谱阈值和超像素光谱属性获得初始影像云检测结果。然后求得直方图均衡化影像的LBP纹理图，根据超像素LBP纹理均值及角二阶矩对初始云检测结果提纯，消除影像中类云地物对云检测的影响。最后对提纯云区域进行区域增长与膨胀处理获得最终的影像云检测结果。通过两组试验对本文方法与多光谱阈值法、树状结构云检测方法进行定性、定量对比评价，试验结果表明，本文方法可以消除山脉、房屋、小面积积雪等类云地物对影像云检测的影响，获得良好的影像云检测结果。下一步将对本文方法进行优化提高其运行效率。