基于改进CDAG云检测算法的Landsat8 OLI云雪识别

胡学谦，孙 林，隋淞蔓，顾晓鹤，孙 宇,何玉龙

(1.山东科技大学 测绘与空间信息学院，山东 青岛 266590；2.国家农业信息化工程技术研究中心，北京 100097)

0 引言

云层的存在，阻挡了太阳辐射能量的传输，严重影响了遥感图像的利用率，给地物分类、地表参数定量反演等应用带来较大的困难，但在遥感成像中，云是难以避免的，为了减少云的影响，需要对云进行检测。由于雪和云在可见光和近红外波段具有相似的反射率特征，当遥感影像中存在积雪时，易导致云检测精度下降[1]。

目前主要的云检测方法有阈值法、人工神经网络法和统计学方法。阈值法因为具有算法简单、速度快和检测精度高的优点而被广泛使用[2]。阈值法云检测的基本思想是基于云和晴空像元的反射率差异和亮度温度差异进行波谱分析和特征通道选择，使用单通道、多通道组合的反射率、亮温等设置一系列云检测阈值。常用的阈值法有ISCCP方法[3]、APOLLO法[4]、CLAVR[5]法和CO2薄片法[6]等。传统阈值法使用典型地物和云的光谱分析或者基于图像统计得出阈值，没有充分考虑下垫面和云形态的复杂性，云检测精度较低。国内外对积雪识别也进行了广泛研究。归一化差值积雪指数(NDSI)的方法由于其简单且精度较高而被广泛使用[7]，并在此基础上发展出了增强型积雪指数(ENDSI)法[8]。1995年，Hall等人在NDSI方法的基础上提出了SNOWMAP算法[9]，用于使用MODIS数据绘制全球积雪地图，是目前最常用的积雪识别方法之一。以上这些方法都没有与云检测有效结合起来。

先验数据库支持的云检测阈值自动生成(CDAG)算法是一种针对不同卫星传感器的云检测阈值生成方法[10]。算法通过在高光谱分辨率的AVIRIS数据上人工标注的云像元和晴空像元，模拟不同传感器的云和晴空像元库，自动、充分分析云和地表不同组分的光谱特征差异，找到最合适的云检测阈值，具有检测效率高、适应性强和检测效果好的优点。由于CDAG算法是基于AVIRIS高光谱数据建立的像元库运行的，受飞行平台的限制，像元库中雪像元较少，并且在构建像元库时，忽略了云和雪反射率特征的相似性，导致对存在雪覆盖的影像云检测精度下降。

为改进原算法在云检测中对雪像元识别效果差，易误判的问题，本文提出了一种增强雪识别的改进CDAG算法。采用Landsat8陆地成像仪(Operational Land Imager,OLI)数据进行验证，结合多种云检测算法，利用改进后的像元库进行迭代计算，自动生成云检测阈值，实现增强雪识别的云检测。经验证明，算法能够准确识别出影像中的云、雪像元。

1 研究区与数据准备

1.1 研究区概况

为了验证增强雪识别的改进CDAG算法的有效性，需要研究区域同时存在积雪和云层覆盖，而青藏高原地区平均海拔4 000 m以上，部分地区终年积雪，地貌类型简单，以岩石、冻土、冰川和裸地为主，并且太阳辐射强，大气条件纯净，适合进行增强云雪识别云检测算法的验证，是本研究的良好研究区域[11]。选择青藏高原地区78°E～85°E，29°N～36°N区域作为研究区，所选区域地处青藏高原西部，覆盖喀喇昆仑山、帕米尔高原、昆仑山西麓和喜马拉雅山西麓，年积雪覆盖日数(SCD)60～120天，平均雪深10 cm左右，冬春季节(11月—次年5月)云量45%～70%，盛行高云[12]。研究区内存在不同类型的云层和积雪场景，存在大量混合像元，增加了云雪识别的难度,是检验算法的理想区域。

1.2 数据介绍

CDAG云检测算法的基础在于高光谱像元库的构建，像元库是利用AVIRIS数据构建的。如表1所示，AVIRIS是第一个测量370～2 507 nm太阳辐射光谱的成像光谱仪，具有光谱分辨率、信噪比和空间分辨率高的特点，并通过224个连续的光谱通道以10 nm为间隔测量上行辐射[13]，能够探测到地物在光谱特性上更微小的差异[14]。

表1 AVIRIS传感器特征

本研究选择Landsat8 OLI数据作为实验数据。Landsat8卫星于2013年2月在美国加州发射，其搭载的OLI共有11个波段，刈幅宽度为185 km，重访周期为16 d，全色波段的空间分辨率为15 m，其余波段为30 m[15]。云和雪的反射波谱特征差异在波长0.38～2.5 μm的可见光波段和近红外波段最为明显，也是云检测的常用波段[16-17]，实验数据选择的Landsat8 OLI 1～7波段和第9波段即位于这一波段范围内。表2为选择的波段及其波长范围。

表2 选用的Landsat8 OLI波段及波长范围

选择研究区域冬春季节的6景Landsat8 OLI数据进行验证，所选数据均存在不同类型积雪及云覆盖。数据如表3所示。

表3 验证数据信息

本文选用的数据为L1T级数据产品，已经过预处理，但为了保证数据的一致性，进行实验验证前，对数据做统一的辐射定标[18]。使用影像元数据中的辐射校正系数，将DN值转化为大气上层反射率。

定标公式为：

(1)

式中，ρλ为大气上层表观反射率；Mρ,Aρ为从Landsat8 OLI元数据中获得的定标系数；Qcal为Landsat8 OLI原始影像DN值；θsz为太阳天顶角。

2 研究方法

2.1 CDAG的云检测算法介绍

CDAG算法是一种自动阈值的云检测算法。云检测的基本原理是根据云在可见光、近红外和短波红外波段与其他下垫面(植被、土壤、水体、岩石和建筑等)的光谱特征差异对云进行检测[19]。CDAG算法使用预先标识出云和典型地物的高光谱像元库，利用数据模拟的方式来实现不同卫星的云检测阈值确定，其实现主要包括以下几个部分：

(1) 高光谱像元库的构建。在选择AVIRIS高光谱数据建立云像元库和晴空像元库时，为了让云检测算法更加完善，构建的阈值有广泛适应性，选择的云和地表类型的种类要全面。云像元库要同时包含厚云、薄云和碎云等像元；晴空像元库要同时包括裸地、岩石、建筑、水体和林地等像元，并且要有足够的像元数量。

(2) 利用高光谱像元库模拟不同的多光谱传感器的像元。CDAG算法的模拟方法采用了Lian等人[20]提出的方法。根据波谱响应函数将窄波段的AVIRIS像元库数据进行模拟得到宽波段的待测数据像元库，模拟公式如下[21]：

(2)

(3) 基于模拟得到的多光谱传感器的云和晴空像元库，进行统计分析，采用机器学习，自动计算不同波段和波段组合在不同阈值下的云像元识别正确率和晴空像元误判率。阈值确定原理如图1所示。

图1 阈值确定原理Fig.1 Principle of threshold determination

为了确定每种算法的最佳阈值，将云像元库和晴空像元库中的像元数据作为训练样本，统计反射率的最大值和最小值，使阈值在2个最值之间以步长0.01变化，云像元正确率和地表误判率也随之变化，当云像元的识别正确率达到最高，且晴空像元误判率开始随阈值变化明显增加时,对应阈值为此种算法的最佳阈值。

2.2 增强雪识别的改进CDAG算法

CDAG算法在达到较高云检测效率和精度的同时，也存在没有充分考虑云和雪反射率特征的相似性，容易对雪像元产生误判的问题，导致影像存在积雪时云检测精度降低。

针对CDAG算法存在的问题，本文提出了增强雪识别的改进CDAG算法。云与雪的识别主要是利用云雪在可见光和近红外波段存在的反射率差异。表4为利用AVIRIS高光谱数据绘制的400～2 500 nm波段范围的云、雪、岩石和裸地等研究区常见目标的光谱曲线。

表4 研究区典型目标像元

由表4可以看出，云在950，1 150 nm处有2个明显的吸收谷，在1 700，2 200 nm处有2个明显的反射峰，与雪在对应波段的反射率存在明显的差异，可以利用这些波段对云和雪进行识别，而岩石和裸土与前二者反射率差异明显，不易产生混淆算法的改进包括以下几部分：

(1) 利用AVIRIS数据构建的像元库，模拟待检测的Landsat8 OLI数据的云像元库和晴空像元库[22-24]。与高光谱的AVIRIS传感器相比，Landsat8 OLI传感器光谱分辨率较低且每一个波段都覆盖了AVIRIS的多个波段。通过数据模拟将对应的若干个连续的窄波段高光谱数据合成对应的Landsat8 OLI传感器的宽波段多光谱数据，得到Landsat8 OLI数据的云像元库和晴空像元库。

(2) Landsat8 OLI像元库的优化。由于AVIRIS数据构建的高光谱像元库并未充分考虑雪对云检测的影响，为提高云检测中云雪识别的效果，选择研究区内经过预处理的Landsat8 OLI数据，在影像上选取总数大于10万个的不同地表覆盖的积雪、碎雪等具有代表性的雪像元，提取其在不同波段下的表观反射率，并将获取的反射率数据加入模拟得到的Landsat8 OLI晴空像元库，在不同波段组合的云检测算法下以机器学习的方式自动生成改进后的最佳云检测阈值。

(3) 基于Landsat8 OLI云像元库和晴空像元库的反射率差异，对其进行波段运算、统计分析、自动生成阈值以及重合度判断等处理，自动生成云检测算法。本文主要选取了单波段、双波段组合、波段比值和波段差值作为云检测算法的公式，最终的云检测算法如式(3)～式(6)所示：

ρi>Ti，

(3)

ρi>Ti&ρj>Tj，

(4)

Ta<ρi/ρj

(5)

Tc<ρi-ρj

(6)

式中，i,j为1～7,9；ρi,j代表第i,j波段的反射率；Ti,j代表第i,j波段的最佳阈值;Ta,b代表不同波段比值算法对应的最佳阈值；Tc,d代表不同波段差值算法对应的最佳阈值。不同波段组合下生成的最佳阈值如表5所示。

表5 不同算法的最佳阈值

(4) 云检测结果的生成。利用生成的算法对待测数据进行实验，每一种算法都会生成相应的云检测结果，利用式(7)对不同算法得到的云检测结果加权合成。云概率计算公式为：

(7)

式中，Fi为云检测结果；Qi为对应云检测结果的权重；N为云检测算法总数。通过加权合成，得到最终的云检测结果，并对检测结果进行精度验证和分析。

2.3 精度验证

为了评价增强雪识别的改进CDAG算法对于Landsat8 OLI数据的云检测效果，本文对实验数据的云检测精度进行定量分析。采用目视解译的方式，在测试数据的标准假彩色影像上区分出云覆盖区域和雪覆盖区域并将其矢量化，与云检测结果进行比较，分别统计识别正确和识别错误的云像元数量和晴空像元数量，将像元分为表6中的4类。

表6 像元分类

采用以下几个指标对云检测精度进行评价：云像元识别正确率(CRA)、云像元漏判率(CRM)、晴空像元识别正确率(SRA)和晴空像元漏判率(SRM)，4个指标的计算如下：

(8)

(9)

(10)

(11)

在以上4个指标中，CRA和SRA可以反映增强雪识别的改进CDAG算法的可靠性；CRM指漏判的云像元在总的云像元中的占比，反映了算法低估云的程度；SRM指漏判的晴空像元在总晴空像元中的占比，反映了算法高估云的程度。

3 结果与分析

图2为增强雪识别的改进CDAG算法的云检测结果示例影像，其中左图为Landsat8 OLI数据5,4,3波段合成的假彩色影像，右图为对应的云雪识别结果，白色为云，灰色为雪，黑色为下垫面。通过对比左图的假彩色影像和右图的检测结果，可以看出原数据中的云和雪都被几乎完整地提取出来，并赋予不同颜色加以区分，体现出算法具有较好的识别效果。

图2 云雪识别结果示意Fig.2 Result of cloud and snow identification

根据云的几何形态和密度，可以将云划分为厚云、薄云和碎云[25]。当云层呈大面积的块状，分布均匀，且无法透过下垫面信息，则认为该云层为厚云。厚云主要包括低云和中云中较厚的云层，一般呈亮白色或灰白色，形状大小没有规律。薄云与厚云的区别在于能否从影像中透过云层观察到下垫面信息。由于薄云轻薄和半透明的特性，其光谱信息往往与地表混合。碎云指面积较小，厚度不均匀，在影像中呈现离散状分布的小云块，由于体积较小，在云检测时常常会被遗漏。不同类型云的反射率特征存在差异，本研究对不同云层类型覆盖条件对改进CDAG算法的适应性进行分析。

3.1 厚云的检测

图3为改进CDAG算法对不同积雪覆盖区域上空厚云的检测结果，左图为5,4,3波段合成的标准假彩色影像；中图为对云和雪的识别效果图，白色为云，灰色为雪，黑色为除雪像元外其他晴空像元；右图为云检测结果图，每种云层条件分别选取4个不同区域和尺度进行分析。

(a) 小尺度区域

厚云是最常见最典型的云层，与其他地物反射率差异大，与积雪的区别也相对明显，较容易识别，识别正确率较高。图3(a)、图3(b)为较小尺度下存在厚云的区域，可以看到，积雪上空的云层可以被完整地识别，积雪也可以较为准确地识别，云和雪具有清晰的边界，极少出现误判。图3(c)、图3(d)为较大尺度区域，可以看到，检测结果中云被准确地识别，图3(d)中云层空隙漏出的积雪地表也被准确识别，并在云检测结果中剔除，体现出方法的有效性。

3.2 薄云的检测

图4为改进CDAG算法对Landsat8 OLI数据不同积雪覆盖区域上空薄云的检测结果。薄云具有可以部分显示下垫面的特点，在常用云检测波段的反射率较厚云更低，在存在积雪的区域更容易产生误判。

(a) 小尺度区域

图4(a)、图4(b)为存在薄云的较小尺度积雪覆盖区域，可以看出，云和雪被准确地识别，极少出现误判。图4(c)、图4(d)为存在积雪和其他地物类型的大尺度区域，从检测结果可以看到，绝大部分积雪和薄云均被准确识别，但也存在少量小块积雪被误判为云。在图像中薄云和积雪地表的过渡部分，存在一定程度的误判，这是由于薄云可以部分透过下垫面信息的特点，无法与下垫面划定明确的边界，对最终的云检测精度产生负面影响。

3.3 碎云的检测

碎云在影像上一般表现为小面积图斑，相互独立且形状不规则，同时具有厚云和薄云的特点，反射率变化较大且边缘较复杂。图5为碎云覆盖区域的云检测结果。

(a) 小尺度区域

图5(a)、图5(b)为存在碎云且有大面积积雪覆盖的较小尺度影像，可以从中看出，碎云和雪均被较为准确地提取，较少出现误判。图5(c)、图5(d)为存在碎云和积雪且尺度相对较大的影像，在这2组图上可以较为清晰地看到碎云和雪的边界，碎云被较好地识别，受到雪的影响较小，误判较少。图5(b)中碎云相对图5(c)较薄，对比2幅图像可以直观看出，图5(b)的云检测结果中出现较多被误判为云的白色亮点，而图5(c)检测结果更为准确，这是由于当碎云较薄时，体现出薄云的特点，更易产生混合像元，增加了误判，降低了云检测的精度，说明对碎云的检测更为困难。

3.4 精度评定

为了定量评价改进云检测算法的准确性，对其进行精度验证。计算不同云雪覆盖条件下云检测结果的4个指标，结果如表7所示。

表7 云检测结果定量评价

表7是图3、图4和图5对应的存在积雪覆盖的厚云、薄云和碎云3类影像云检测结果的定量评价和总体的精度评定。从表7可以看到，改进的云检测算法对于Landsat8 OLI数据的云检测达到了较高的精度，总体的云像元识别正确率和晴空像元识别正确率均达到了90%以上。参考图3、图4和图5的云检测结果影像可以看到，改进的算法较少出现云像元和雪像元的误判，并且对其他晴空像元的误判率也较低。分别分析厚云、薄云和碎云的定量评价结果可以看到，不同类型的云对检测结果存在影响，厚云的整体检测精度较高且检测结果较稳定，云识别的平均正确率达到了95.2%，晴空像元的平均漏判率为4.6%；而薄云和碎云的检测精度相对低，云识别的平均正确率分别为91.9%和92.2%，晴空像元的平均漏判率分别为3.8%和4.6%。总体的云识别平均正确率为93.1%，总体的晴空像元平均漏判率为4.3%。通过精度验证可以看出，本文提出的增强雪识别的改进云检测算法对Landsat8 OLI数据的云检测结果具有很高的精度，并且能够正确区分云雪像元。

4 讨论

云的存在阻碍了太阳辐射能量的传输，增加了地物识别和地表参数反演的难度，导致遥感图像利用率下降，为了降低云的影响，需要进行云检测[26]。而在高海拔地区和中高纬度地区的寒冷季节，地表常会存在积雪，雪和云的反射波谱特征具有相似性，增加了此类区域的云检测难度，对遥感影像的判读和使用产生消极影响[27]。为了提高对存在积雪覆盖区域的云检测精度，本文提出了一种增强雪识别的改进CDAG算法，并对这种方法的有效性进行了验证。本文利用青藏高原地区的多景Landsat8 OLI影像数据对改进CDAG算法进行了实验，并对厚云、薄云和碎云3种典型云覆盖条件下的积雪覆盖区域分别进行了云检测。

从图3、图4和图5的云检测结果可以看出，本文方法针对不同云层类型均能得到准确的云检测结果，能够正确地识别出云和雪，在云检测结果中云和雪及其他地貌类型具有清晰的边界，对于云的识别完整、清晰自然，达到了很好的效果。从表7的精度评价结果可见，总体的云像元识别正确率达到了93.11%，晴空像元的识别正确率达到了95.65%，较传统的云检测方法有了很大的提升，显示了本文方法的有效性。但不同类型云覆盖下的检测精度仍存在差异，其中厚云的云像元识别正确率达到了95.2%，明显高于薄云和碎云的92%左右，对其进行分析，有以下原因：① 厚云的形态较为稳定，云层较厚且与地面有明确的边界，不易产生混合像元；② 薄云能够部分反映地面信息，碎云形态破碎边界杂乱，这些特性导致了混合像元的产生,从图中也可以看到这些特征，混合像元的存在模糊了云与雪的波谱特征差异，导致薄云和碎云的检测精度较厚云偏低。

5 结束语

本文在CDAG算法的基础上改进得到了增强雪识别的改进CDAG算法，主要改进包括以下2个方面：① 针对雪覆盖地表，扩充了像元库，增加雪像元在晴空像元中的占比，通过机器学习提高了云检测阈值的准确性，提高了对雪像元的识别正确率；② 采用了单波段、多波段组合、波段比值和波段差值多种算法，提高了云检测的精度。改进后的算法在保留了CDAG算法识别正确率高、自动生成阈值等优点的基础上，降低了对雪像元的误判率，在实验中取得了良好的检测效果，在积雪覆盖的区域平均云像元识别正确率达到93.11%。目前算法存在的在薄云和碎云条件下，云检测精度较厚云条件下降以及像元库数据由人工选择组成，主观性较强，可能导致数据库的构建存在偏差的问题，是下一步需要研究的内容。