一种改进的遥感影像云检测方法

张昊,焦瑞莉,乔聪聪,霍娟,宗雪梅

(1.北京信息科技大学 信息与通信工程学院,北京 100101;2.中国科学院大气物理研究所,北京 100029)

0 引言

利用云和云阴影区域与其他区域的差别实现气溶胶参数的反演是一种可行的方法[1],并且由于云阴影检测通常由云检测实现,因此在遥感影像中准确的云检测是后续气溶胶研究的重要前提。

近年来,国内外在遥感影像的云检测方面进行了大量的研究,并提出多种不同的方法。云检测的方法可分为两大类:光谱阈值方法和深度学习方法。

光谱阈值方法利用卫星遥感影像的若干个光学波段的反射率进行阈值判别,从而进行遥感影像的云检测。光谱阈值方法可分为单时相方法和多时相方法。单时相方法通过对单幅遥感影像进行光谱反射率的阈值判别得到云检测结果。Luo等[2]提出使用阈值决策矩阵来识别MODIS卫星影像中的云。Zhu等[3]提出Fmask方法,对Landsat卫星影像采用光谱阈值和概率阈值相结合的两步阈值法进行云检测,并结合卷云波段和多种辅助数据提升检测准确率[4-6]。多时相方法是利用多幅不同时间、相同地点的卫星影像组成时间序列,并结合光谱反射率阈值进行判别得到云检测结果。文献[7-8]分别在Fmask方法的云检测结果基础上,利用时间序列进一步检测。然而,单时相光谱阈值方法主要通过光谱特征进行检测,因此针对某些和云具有相似光谱特征的区域效果较差;多时相光谱阈值方法虽然通过时间序列在一定程度上提高了云检测的准确率,但是其对数据量需求较大,运算效率较低。

深度学习方法利用庞大的数据集对云检测模型进行训练,然后通过训练完成的模型对被测影像的每一个像素进行预测分类。Wieland等[9]使用U-Net模型进行遥感影像中的云检测。Li等[10]提出了一种以多个残差单元为编码器,以相同结构的多个残差单元和反卷积为解码器的多尺度特征融合云检测方法。Zhang等[11]提出了一种以多个残差单元和小波变换模块为编码器,以ASPP金字塔和空间注意力机制为解码器的多尺度特征融合模型。Mehdi等[12]提出了一种以不同尺寸的卷积核组成的并联结构为编码器,以反卷积为解码器的多尺度特征融合模型。Qu等[13]提出了一种以条状残差网络为编码器,以双线性插值和双注意力机制为解码器的多尺度特征融合模型。然而,深度学习方法在具有良好检测结果的同时,需要庞大的数据量和较高的硬件需求,且对训练参数的初始化较为敏感。

多特征联合(multi-feature combined,MFC)是一种针对单幅遥感影像,并结合了光谱特征、几何特征和纹理特征,应用在高分一号的云检测方法[14]。由于MFC方法在高分一号数据上表现出良好的性能,因此本文采用波段范围更加丰富的Landsat 8数据进行MFC方法的性能评价,并提出两点改进。一是针对几何特征与云相似的非云区域,利用灰度共生矩阵提取纹理特征,提高云识别正确率。二是针对冰雪检测的局限性,引入一种光谱特征对冰雪区域进行检测和滤除。

1 理论基础

1.1 光谱特征和引导滤波生成初步云检测结果

MFC方法首先通过3个光谱特征的阈值判别生成一个包含大部分厚云的云检测结果,在最大限度上减小非云像素引起的错分误差。

第1个光谱特征是雾霾最优变换(haze optimized transformation,HOT)。其原理是在蓝色波段和红色波段中,云和大部分地表以及薄雾有着不同的光谱特性,前者的HOT值通常大于后者的HOT值。

第2个光谱特征是可见光波段比值(visible band ratio,VBR)。VBR能够排除具有显著颜色特征的非云像素。云在可见光波段通常为灰色或白色,其VBR值接近1。

第3个光谱特征是红色波段反射率,其原理是云在红色波段的反射率较高。

MFC方法分别将3个光谱特征的阈值设置为0.13、0.7和0.07。如果像素同时满足上述3个阈值条件,则将其检测为云。

在得到包含大部分厚云的云检测结果之后,MFC方法认为薄云通常分布在厚云周围,即在云边界周围存在一个从厚云到薄云的过渡区域。因此在光谱特征检测的基础上,采用引导滤波器检测每个云边界遗漏的薄云。

引导滤波器包括引导图像、输入图像和输出图像。在MFC方法中,引导图像是蓝色波段、绿色波段和红色波段的合成影像,输入图像是光谱特征阈值判别得到的二值云检测图像,输出图像是引导滤波后的二值云检测图像。

引导滤波器描述了引导图像和输出图像之间的局部线性关系,即输出图像是引导图像以某个像素为中心的窗口的线性变换。对于以某个像素为中心的窗口,线性变换的系数包括ak和bk。其中ak与引导图像在此窗口内的所有像素的总数、均值、方差,以及正则化参数有关;bk与引导图像和输入图像在此窗口内的所有像素的均值,以及ak有关。

对于任意像素i,其包含在多个不同的窗口,因此像素i的输出应考虑到包含其本身的所有窗口,即系数a和b分别为包含其本身的所有窗口的ak和bk的平均值。

MFC方法分别针对水区域和陆地区域设置了不同的云边界检测方法。首先采用归一化植被指数(normalized difference vegetation index,NDVI)和近红外波段反射率将水区域和陆地区域进行划分,原理在于后者的NDVI值通常高于前者的NDVI值,且前者在近红外波段的反射率较低。MFC方法使用了两组不同的阈值判别进行水区域的检测。第1组将NDVI阈值设置为0.15,近红外波段阈值设置为0.2;第2组将NDVI阈值设置为0.2,近红外波段阈值设置为0.15。如果满足其中任意一组阈值条件,则检测为水区域,否则为陆地区域。

对于水区域,MFC方法采用大津阈值法对引导滤波器输出的灰度图像进行阈值分割生成二值云检测图像,阈值设置为0.12;对于陆地区域,在对引导滤波器输出图像进行阈值分割的基础上,再次使用HOT进行阈值分割,以防止云周围非云像素的干扰,阈值设置为0.08。

然而,初步云检测结果可能包含一些非云区域,如冰雪、建筑物等,它们与云区域具有相似的光谱特征,易与云区域混淆。因此在初步云检测结果中对这些非云区域进行滤除是十分必要的。

1.2 滤除非云区域

在滤除非云区域之前,MFC方法首先在初步云检测二值图像中进行8邻域连通,形成若干个连通域,以连通域为单位提取多种特征并逐一进行判别。

MFC方法通过利用连通域提取分形维数指数(fractal dimension index,FRAC)和长宽比(length width ratio,LWR)这两种几何特征对海岸线、道路、建筑物等具有几何形状的区域进行滤除。FRAC由连通域的周长和面积确定,反映了连通域形状的复杂性。LWR由连通域最小外接矩形的长和宽确定,反映了连通域最小外接矩形的长与宽之间的关系。相比于云区域,这些非云区域的FRAC值和LWR值较大。

然而,其中一些非云区域的FRAC和LWR和云区域较为相似,仅通过几何特征无法将它们与云进行区分。因此,本文在几何特征的基础上采用灰度共生矩阵进一步提取连通域的纹理特征进行判别。

灰度共生矩阵反映了图像灰度关于方向、距离、变化幅度的综合信息,它是分析图像局部模式结构及其排列规则的基础,其统计了图像在某一方向上相距一定距离的两个灰度值同时出现的概率。

作为纹理分析的特征量,往往不是直接应用灰度共生矩阵,而是在灰度共生矩阵的基础上再提取纹理特征量。均值和[15]是纹理特征量之一,是区域内像素平均灰度值的度量,反映了区域的明暗深浅。对于具有较高灰度值的区域,均值和较大;对于具有较低灰度值的区域,均值和较小。

在可见光波段,云具有较高反射率,即较高灰度值,均值和较大;大部分地表具有较低反射率,即较低灰度值,均值和较小。因此选取蓝色波段用于计算灰度共生矩阵的均值和,提取纹理特征。

考虑到方法的效率以及灰度共生矩阵的计算量,本文从4个方向中选取45°方向,相对距离设置为1,灰度级数设置为16,窗口大小设置为11。为了避免个别异常值的影响,采用以连通域为单位的判别方式,即统计每个连通域中灰度共生矩阵均值和大于阈值T的像素总数,并计算其与所在连通域的面积的比值,如果比值小于0.5,则判别为非云区域并进行滤除,否则判别为云区域。对于阈值T的选择,本文选取了239个误判为云的非云连通域和276个云连通域,并计算在不同阈值T的选取下分类的正确率。首先在区间[2,3]之间以0.1为间隔对T进行取值,计算分类正确率,结果如图1所示。可以看出,阈值在区间[2.5,2.6]时,云分类正确率所受影响较小,且非云分类正确率较高。为了得到更加精确的阈值,在区间[2.5,2.6]内以0.01为间隔对阈值T进行取值,计算分类正确率,结果如图2所示。可以看出,若要同时兼顾二者的分类正确率,则将阈值T设置为2.56时最佳。

图1 阈值区间为[2,3]的分类正确率曲线

图2 阈值区间为[2.5,2.6]的分类正确率曲线

对于冰雪区域,MFC方法选取了图像尺寸为100像素×100像素的84个样本,计算它们的局部二值模式(local binary patterns,LBP)纹理直方图,并依据卡方距离进行判断。分析认为此方法具有一定的局限性。一是由于MFC方法针对面积小于40 000个像素点的连通域进行判别,而遥感影像中有可能存在大面积的冰雪区域,因此可能会造成漏判。二是检测结果取决于样本的选择,如果样本选择不佳,则可能影响最终的检测结果。因此,本文以像素点为单位,引入一种光谱特征对冰雪区域进行判别并滤除。

图3是云和雪从可见光波段到短波红外波段的反射率折线图[16],其中蓝色波段的波长范围为0.45～0.52 μm,短波红外波段的波长范围为1.56～1.66 μm。由图3可知,在蓝色波段,云和雪的反射率较高,通常大于0.6;在短波红外波段,云和雪的反射率相对于可见光波段有所降低。不同的是,云的反射率在0.6浮动,下降量较小,而雪的反射率小于0.2,下降量较大。因此,通过蓝色波段和短波红外波段进行光谱特征提取,特征表达式见式(1)。

图3 云和冰雪的反射率曲线

SI=2×B6-B2

(1)

式中:B6表示短波红外波段反射率;B2表示蓝色波段反射率。如果像素的SI值小于阈值0,则将其检测为冰雪,否则将其检测为云。

2 实验与分析

2.1 实验方案

实验所采用的数据是Landsat 8 OLI卫星遥感数据。为了验证改进MFC方法在不同地表类型下的性能,本文选取了涵盖森林、冰雪、水、城市、瘠地等不同地表类型的遥感影像进行性能评价。

用于计算性能参数的参照结果采用LandsatLook 8-bit质量评估图像。LandsatLook 8-bit质量评估图像采用不同颜色标注了遥感影像中每个像素的类别。因此,通过对LandsatLook 8-bit质量评估图像中的云像素进行标记得到云参照结果,将其作为真实标记,并分别与原始MFC方法和改进MFC方法的云检测结果进行性能参数的计算。

实验通过计算总体精度(overall accuracy,OA)、用户精度(user accuracy,UA)、生产精度(producer accuracy,PA)这3个性能参数对改进MFC方法进行性能评价,并与原始MFC方法进行对比。将所有像素划分为云类或非云类,通过衡量云检测结果和云参照结果中像素类别的一致性和差异性进行精度计算,计算方法如式(2)、式(3)、式(4)所示。

(2)

(3)

(4)

式中:MP和MN分别表示参照结果中云的像素数量和参照结果中非云的像素数量;MTP表示参照结果中标记为云且检测结果标记为云的像素数量;MTN表示参照结果标记为非云且检测结果标记为非云的像素数量;MFP表示参照结果标记为非云且检测结果标记为云的像素数量;MFN表示参照结果标记为云且检测结果标记为非云的像素数量。

总体精度OA衡量了总体分类的正确性;生产精度PA衡量了遗漏分类的程度;用户精度UA衡量了错误分类的程度。

图4是实验流程图,总共分为4个步骤。第1步是数据预处理,将遥感影像的灰度值转换为反射率;第2步是结合光谱阈值分割和引导滤波器生成初始云检测结果;第3步是利用几何特征、灰度共生矩阵的纹理特征,以及本文引入的光谱特征滤除初始云检测结果中的建筑物、冰雪等非云区域;第4步是进行实验结果分析与对比。

图4 实验流程图

2.2 实验结果与分析

图5、图6、图7列出了3组同时含有云和建筑物、道路等非云区域的不同地点的结果,从左至右依次为遥感影像、原始MFC方法的云检测结果以及改进MFC方法的云检测结果,其中白色表示云区域,黑色表示非云区域。每组图像中,图(a)是区域尺寸为1 024像素×1 024像素的影像,其中红色矩形标记了区域尺寸为256像素×256像素的部分区域。图(b)是图(a)中红色矩形区域的放大图像,其中黄色矩形标记了其中部分非云区域。

图5 含有建筑物等非云区域的遥感影像云检测结果一

图6 含有建筑物等非云区域的遥感影像云检测结果二

图7 含有建筑物等非云区域的遥感影像云检测结果三

从上述图中可以看出,原始MFC方法的云检测结果中仍然存在着一些非云区域,说明仅通过几何特征并不能有效地从云检测结果中滤除这些非云区域。而改进MFC方法的云检测结果中存在着更少的非云区域,说明进一步通过灰度共生矩阵提取纹理特征能够对这些非云区域进行更加有效的滤除。

本文针对冰雪区域的检测和滤除所引入的光谱特征也表现出良好的效果。图8列出了3组不同地点的遥感影像,每组从左至右依次为尺寸1 024像素×1 024像素的遥感影像及对应的改进MFC方法的云检测结果,其中白色表示云区域,黑色表示非云区域。图中红色矩形标记了部分冰雪区域。图8(a)是一幅含有大面积冰雪区域的遥感影像,可以看出,改进MFC方法滤除了大面积冰雪区域,防止其与云区域混淆。图8(b)和图8(c)是两幅同时含有云和冰雪的遥感影像,可以看出,改进MFC方法有效地区分了云和冰雪,并将冰雪区域进行了滤除。

图8 含有冰雪的遥感影像云检测结果

表1和表2分别针对云量少于50%和云量高于50%的遥感影像,列出了原始MFC方法和改进MFC方法的总体精度、用户精度和生产精度。表3列出了不同云量条件下的改进MFC方法性能参数的增长百分比。

表1 云量>50%的遥感影像云检测结果的性能参数 %

表2 云量<50%的遥感影像云检测结果的性能参数 %

表3 改进MFC方法云检测结果性能参数的增长百分比 %

根据表1和表2可得出,相比于原始方法,改进方法在总体精度、用户精度和生产精度上都有所提升;且相比于总体精度和生产精度,用户精度的增长百分比较大,原因在于改进方法能够滤除更多的建筑物、道路等非云区域,使得云检测结果中包含更少的非云像素。

根据表3可得出,对于云量大于50%的遥感影像,改进MFC方法的总体精度、用户精度和生产精度的增长百分比分别为0.46%、0.59%和0.27%,说明在云量较多的情况下改进方法的提升效果不显著;对于云量小于50%的遥感影像,总体精度、用户精度和生产精度的增长百分比分别为1.21%、4.37%和1.29%,说明在云量较少的情况下改进方法的提升效果较为显著。

分析获得上述结论的原因在于,通常相比于云区域,遥感影像中建筑物、道路等非云区域的面积较小。对于云量较多的遥感影像,这些非云区域所占的像素数量比例较小,其误检对云检测结果的影响较小,即对精度的影响较小,从而改进方法的提升效果不明显;反之,对于云量较少的遥感影像,这些非云区域所占的像素数量比例较大,其误检对云检测结果的影响较大,即对精度的影响较大,从而改进方法的提升效果较为明显。

3 结束语

本文验证了MFC方法在Landsat 8遥感影像数据上的性能,并针对该方法进行了两点改进。一是采用灰度共生矩阵提取纹理特征,并结合原始方法采用的几何特征对初始云检测结果中检测为云的非云区域进行检测和滤除。二是引入一种光谱特征,对初始云检测结果中可能存在的冰雪区域进行检测和滤除。通过实验数据的收集和分析,验证了MFC方法能够有效检测Landsat 8遥感影像中的云,并且验证了改进方法的有效性,在一定程度上弥补了原始方法的缺点。但是改进方法仍存在局限性。一是虽然使用灰度共生矩阵提取纹理特征滤除了更多的非云区域,但是同时也会滤除一些极薄的薄云,原因在于一些薄云区域在蓝色波段的反射率较低,与大部分非云区域相似。二是虽然引入的光谱特征能够有效区分云区域和冰雪区域,但是其无法使用在一些没有短波红外波段的卫星遥感影像上。三是对于云量较多或者含有较少易与云区域混淆的非云区域的遥感影像,改进方法对云检测结果的提升效果有限。在后续的研究中,可以结合时间序列或其他方法,以进一步提高云检测精度。