国产高分辨率卫星影像云检测

刘燕，张力，王庆栋，王春青，韩晓霞

(1.兰州交通大学 测绘与地理信息学院，兰州 730070；2.中国测绘科学研究院，北京 100036；3.青海省基础测绘院，西宁 810001；4.地理国情监测技术应用国家地方联合工程研究中心，兰州 730070；5.甘肃省地理国情监测工程实验室，兰州 730070)

0 引言

随着计算机技术、空间技术和信息技术的发展，遥感作为探索世界的重要手段，为人类的生存、繁荣和可持续发展提供服务。遥感数据作为地物目标识别和分类的主要数据源，其处理趋向于计算机自动化、智能化[1]。云作为正常的自然现象，覆盖了地球表面百分之六十以上的区域，在获取的遥感影像中存在大量区域被云覆盖的现象，导致部分成像信息丢失[2]，严重影响地物的识别和分类、影像匹配和镶嵌等后期处理，也对影像的配准、融合等处理造成不便。因此，云的自动识别检测是遥感影像数据处理中十分重要的过程。

云在遥感影像具有特殊的光谱特性，根据这种特性进行光谱测试，制定出最优阈值，从而检测云区域，如在Landsat数据中广泛使用的Fmask算法[3]和基于云的波谱特性，尤其在红外波段的特性对MODIS数据进行云检测等方法[4]。这些方法需要影像具有较多波段，对于单波段或红绿蓝三波段的影像并不适用。故研究者提出基于较少波段(四个波段以下)的遥感影像云检测方法，这些方法可分为物理方法、基于纹理和空间特性的方法及模式识别三种方法[5]。物理方法是对遥感影像的灰度值进行灰度均值、方差等的统计和计算，结合动态阈值、自适应阈值等方法检测云区域[6-8]，如周雪珺等[9]将遥感影像切分成小块子图，统计子图的灰度均值和方差，根据阈值将云进行分类，快速地进行云检测；陈振炜等[10]将影像分成小块，根据影像的二阶矩值和一阶差分退化率的阈值将影像分为有云影像和无云影像。这种方法可以快速地进行云检测，但单纯依赖于遥感影像的灰度特征，对于具有复杂地物的遥感影像有可能会检测失误。基于遥感影像纹理和空间特性的方法是通过图像灰度共生矩阵、LBP等算法揭示云的纹理特征，并结合阈值算法得到云检测结果[11]，如黄宇[12]提出一种改进的 Itti 显著性模型，并研究云层和下垫面样本的各类纹理特征，利用较少波段信息快速有效地实现遥感影像云检测；董志鹏等[13]提出一种基于对象光谱与纹理的高分辨率遥感影像云检测方法，通过均衡化影像直方图、简单线性迭代聚类算法、纹理均值及角二阶矩等方法获得良好的影像云检测结果，由于云的种类较复杂，因此难以准确地区分与云纹理相近的地物。模式识别中的支持向量机、监督分类、深度学习等方法可以获得较好的云检测结果[14-15]，如赫英明等[16]提出了基于支持向量机的遥感影像分类方法，分析了云检测过程中的特征提取和选择，建立了基于支持向量机的遥感影像分类模型，可将云与陆地、水体、积雪准确地区分；高军等[17]以全卷积神经网络为基础进行云检测，利用其自动提取深层隐含特征等特性保留特征信息，最后结合全连接条件随机场模型进行云系边缘优化获得效果好的云检测结果。但这种方法需要大量的云检测样本且受限于对训练样本的依赖，对于与样本反差较大的遥感影像效果不太理想。

基于上述原因，本文通过分析遥感影像中云区域光谱和纹理特征，提出一种快速、普适性强的多源高分辨率遥感影像云检测方法。综合多项式阈值法和大津法的优势，提出了一种联合阈值分割算法，不仅能够检测多云影像，而且解决了少云影像无法检测的问题，从而提高云检测的正确率。同时根据类云地物边缘梯度值变化大而云边缘区域变化平滑的特点实现了类云地物的去除。

1 理论基础

1.1 影像云判别

影像直方图可以表示数字图像中亮度分布，标绘出图像中每个亮度值的像素数。将彩色影像转为灰度影像，并做直方图统计，为了更直观地显示，将纵坐标(像素数量)归一化到0～255。从图1中可以看出，当影像中含有云区域时，直方图大致呈“大U”形，且分布主要集中在横坐标[200，255]区间内，纵坐标不为零。从图2中可以看出，当影像不含云时，

图1 有云遥感影像和直方图

图2 无云遥感影像和直方图

直方图分布主要聚集在横坐标[200，255]区间内，纵坐标为零。由此，本文根据直方图横坐标在大于直方图阈值(TH=230)时，纵坐标是否有值，将影像区分为有云影像和无云影像。

1.2 梯度增强

如图3所示，遥感影像云区域边缘像素值较中心区域小，但与非云区域的像素值相差较大，因此为了提高云范围检测准确率，对影像进行梯度增强，增强云边缘区域强度。将灰度影像和梯度影像按照0.7和0.3的权重进行叠加并拉伸到0～255，得到梯度增强影像，如图3所示，梯度增强后的影像云过渡区域和云层较薄区域得到明显增强。

图3 梯度增强对比图

1.3 联合阈值分割

大津法又称为最大类间方差法，是经典的阈值分割算法，经常用于遥感影像中云、阴影和水体等提取，它是基于某一阈值将图像直方图的灰度像素分为两类，计算两类的类间方差，通过不断迭代，使得类间方差达到极小值，从而得到阈值。类间方差法对噪声和目标大小十分敏感。当目标与背景的大小比例悬殊，比如云量较少，或者同时有水体和阴影等多种地物的遥感影像时，阈值难以确定，如图4(b)和图5(b)所示，大津法阈值分割将部分地物误检测为云。

多项式拟合直方图阈值分割法是用五次多项式对直方图曲线通过最小二乘的计算方法进行曲线拟合，倒序排列求拟合曲线的第一个最小值作为阈值，进行影像分割。这种方法对于检测云量较少的遥感影像效果较好。对于有大量水域的影像，大津法检测效果较差，将大部分陆地误检测为云区域，而多项式阈值分割方法可以较准确地分割出云区域，如图5(c)所示。但是多项式阈值分割对于云过渡区域检测效果不理想，如图6(c)所示，多项式阈值分割法只能检测出云的主体区域，无法检测云边缘过渡区和云量较薄区域，相反，如图6(b)所示的大津法能够对以上情况进行检测。因此，当云区域较大时，大津法检测效果要优于多项式阈值分割法，当云区域较小时，大津法效果要劣于多项式阈值分割法。

图4 少云影像阈值分割

图5 有水域影像分割

图6 多云图像阈值分割

统计原始影像中云区域的平均灰度值，并生成云样本平均灰度离散图(图7)，发现99%的云区域平均灰度值都大于阈值(TA2=150)，根据这个特性可验证联合分割阈值云检测结果的正确性。85%云区域的平均灰度值大于阈值(TA1=180)，故用TA1验证Otsu阈值分割的正确性。本文首先选用大津法阈值分割，通过计算原始影像中阈值分割云结果区域的灰度均值，作为验证检测结果是否正确的条件。如果检测出来的云区域灰度均值小于TA1，表明大津法阈值分割错误，重新选用多项式拟合直方图阈值分割。最终，如果云区域灰度均值小于TA2，则进一步进行分块迭代云检测。

图7 云样本平均灰度值散点图

1.4 波段光强差误检测区域去除

云作为一种特殊的自然形态，具有较强的反射率，在遥感影像上为白色，对于具有三波段以上的影像，在红、绿、蓝三个波段上的强度差异较小，根据样本计算得出云区域的像素色彩强度差值都小于阈值(TRGB=25)，即(|红-蓝|<25)且(|红-绿|<25)且(|绿-蓝|<25))。根据这个特性，可进一步将云和类云地物区分，如图8所示。

图8 根据红绿蓝波段色彩差值去除类云地物

1.5 边缘梯度类云地物去除

遥感影像云检测中云与类云地物(包括白色的建筑物、岩石、冰、雪等)的区分一直是一个技术难题。云由于种类不同，具有不同的形态和灰度特征，类云地物也由于气候、地形的复杂具有不同的特征。如图9(a)所示，云区域从中心到边界呈现慢慢过渡，从厚变薄，颜色从白色过渡到地面颜色；如图9(b)所示，边界清晰的白色类云地物边界颜色过渡较快，在五个像素内就从白色过渡为地面颜色。因此，根据上述特征，针对其中边界清晰的冰雪和建筑物可以制定算法区分云和类云地物。

图9 云与类云地物影像

像素值变化的快慢，即图像边缘梯度可作为判断准则。如图10所示，分别为云区域和类云地物从中心到边缘像素值的变化图。可以看出，边界清晰的白色地物灰度值变化剧烈，故梯度值大；而云边界变化缓慢，故梯度值较小。

图10 过渡区域灰度值变化图

图11 单独云块过渡区梯度均值统计

根据云与类云地物过渡区域，计算过渡区域原始影像的梯度值，同时统计出每个单独云块和类云块，即单独连通域边缘的梯度均值。如图11(a)所示，横坐标为单独过渡区域的样本号，竖坐标为单独过渡区域的梯度均值，从图中可以看出梯度均值分为上下两大类，并且中间有明显的间隔，故推断出云过渡区域的梯度均值为下部分，类云地物过渡区域的梯度均值为上部分。对梯度均值的数量进行统计，如图11(b)所示，画出梯度均值数量的曲线图，曲线有两个波峰，梯度均值小的波峰为云区域，梯度均值大的波峰为类云区域，两个波峰之间的波谷最低点即为区分云和类云地物的阈值TG。

2 算法描述

根据上述分析，本文提出了一种针对多类型的国产高分辨率卫星影像云检测算法，主要分为四个部分：影像云判别、大津法和多项式拟合直方图联合云分割、波段光强差误检测区域去除和边缘梯度类云地物去除，如图12所示。

图12 算法流程图

具体步骤如下所示。

第一步影像预处理。首先，将影像降位到8 bit并四倍降采样，以减少数据量；其次，对影像进行中值滤波去除噪声，生成彩色影像A；然后，将影像A转化为灰度影像，并计算灰度影像的梯度影像；最后，将灰度影像和梯度影像按照0.7和0.3的权重相加，得到梯度增强后的影像B并生成直方图。如果当直方图横坐标大于TH时，纵坐标不为零，表明影像存在像素值大于TH的区域，则该影像为有云影像，否则为无云影像。

第二步大津法和多项式拟合直方图联合阈值分割。利用大津法对影像B进行阈值分割，得到初步云检测结果。对于单波段影像，此步骤跳过，对于三波段以上影像，计算彩色影像A红、绿、蓝波段之间像素值之差，如果云检测区域内像素值之差都大于TRGB，则判断为非云区域并将其去除。通过计算云区域部分影像B的灰度均值来判断云检测结果是否正确，如果灰度均值大于TA1，则保留云初步检测结果；如果灰度均值小于TA1，则表明大津法不适合该影像云检测，则对影像B进行多项式阈值分割并根据阈值TRGB去除非云区域，剩下区域作为云初步检测结果。

第三步判断云检测的正确性。计算云初步检测区域内影像B的灰度均值，若灰度均值

第四步类云地物去除。对云检测得到的二值图像进行腐蚀运算，并与腐蚀前的二值图像进行逻辑运算得到云过渡区域。计算影像B中每块云过渡区域的梯度均值，如果梯度均值>TG，则表明该区域为类云地物，并将该区域去除；如果梯度均值≤TG，则该区域为云区域。

3 实验与分析

本文选用高分一号、高分二号、资源三号和天绘一号高分辨率影像数据进行实验，如表1所示，所选择的影像涉及多种地表类型，包括平坦的裸地、大型水域、雪及山区等，具有不同的纹理特征，包括少云影像和多云影像，因此有很大的光照差异和场景变化。选取经典的大津法阈值算法[18]、文献[10]的树状结构云检测算法，以及人工经验阈值提取结果和本文算法结果进行对比。如图13所示，对于GF-1-1少云影像，大津法无法有效检测云区域，树状结构算法将部分地物检测为云区域，本文算法正确地检测到了云区域，但是云区域不够完整。对于具有大面积水域的影像，如GF-1-2，大津法和树状结构法将大量陆地区域误检测为云区域，而本文算法能够准确地检测云区域，并且和人工提取结果相似。对于具有类云地物(如雪、白色建筑物等，并且地物边界清晰)的影像，如GF-1-3，大津法和树状结构法都将类云地物误检测为云区域，而且树状结构法漏检测大量的云区域，而本文算法不仅能有效地去除类云地物，而且云检测效果较为理想。对于多云影像，如GF-1-4，大津法和本文算法准确地检测出云区域，而树状结构法将部分地物误检测为云区域。对于具有单波段影像，如ZY-3，树状结构法将部分地物误检测为云区域，本文算法和大津法能准确地检测出云区域。对于具有高亮度建筑物的影像，如TH-1，大津法将高亮建筑物误检测为云区域，树状结构法检测的云区域比实际云区域小，本文算法虽然将少部分建筑物误检测为云，但总体效果要优于大津法和树状结构法。

表1 实验数据介绍

为了客观地分析云检测结果，本文选用目视判别与定量评价相结合的方式验证云检测结果的有效性[19]，包括正确率(accuracy rate，AR)、总体正确率(overall accuracy，OA)、漏检率(miss rate，MR)和错误率(false rate，FR)来定量分析检测结果，计算如式(1)至式(4)所示。

(1)

(2)

(3)

(4)

式中：MFN代表云区域被误检测为非云区域的像素数；MFP代表非云区域被误检测为云区域的像素数；MTP代表检测正确的像素数。

根据上述评价指标，本文将大津算法、树状结构算法和本文算法检测结果进行统计，如表2所示，发现：对于少云影像，具有特殊地形(有多面积水域、类云地物)的影像，本文算法的正确率、总体正确率都高于大津法和树状结构法，漏检率也低于大津法和树状结构法；对于云量较多，并且无类云地物的影像，本文算法的正确率和总体正确率与大津法相差不大，高于树状结构法；对于具有高亮建筑物的有云影像，本文算法的正确率和总体正确率高于大津法和树状结构法，漏检率和错误率低于大津法和树状结构法，具有明显的优势。

图13 云检测算法视觉对比

表2 云检测结果精度评价表 %

上述实验表明，本文云检测算法适用性强，能够适用于多种类型的高分辨率遥感影像，云检测正确率高，误检率低，能够去除多种类云地物，对于复杂地物的有云影像具有较好的云检测结果。

4 结束语

本文在详细分析了多云影像和少云影像的光谱特征以及云与类云地物边缘梯度差异的基础上，提出了一种多源高分辨率遥感影像云检测算法。该算法根据影像中云区域的光谱特性并联合大津阈值分割法和多项式拟合阈值分割法的优点对云进行初步检测，然后利用云和类云过渡区域纹理特征差异去除类云地物，并根据光谱之间的差异特征对云区域进一步检查，提高云检测结果的准确性。本算法有以下优点。

1)能够自动实现多源影像云检测，不仅适用于具有红、绿、蓝三个波段的影像，而且适用于单波段影像，具有较强的普适性。

2)该算法不仅能够准确地检测多云影像，而且对于少云影像也具有理想的检测效果。

3)该算法能够去除边缘清晰的类云地物。

4)算法复杂度小，耗时短，可实时对影像进行云检测，并且有较高正确率和较低的误检率。

然而，本文提出的算法依然存在一定的缺陷，需要在今后的工作中继续改进。一方面，对于较薄云层会出现漏检；另一方面，无法去除过渡区不是很明显的类云地物。因此，后期会详细分析类云地物和云区域的特征差异，继续对算法进行改进。