遥感影像数据中经常会有云覆盖, 云会遮挡地表真实反射, 干扰和影响数据反演的精度, 降低数据利用率；因此, 影像预处理中的一个关键环节是完成对云的识别[1].
目前有许多对云识别方法的研究[2-4], 物理阈值法因其简单高效, 应用最广泛.如Irish等[5-6]提出自动云覆盖评估(automatic cloud cover assessment, ACCA)方法, 该方法通过两组滤波器, 利用地物的先验特征进行影像的云识别.Zhu等[7-9]提出Fmask方法.该方法利用地类的光谱特性进行一系列的光谱测试, 选取NDVI和NDSI等参数的最优阈值以识别较明显的云, 从而得到潜在云层；然后根据统计概率学中的概率因子计算云概率参数, 得到最终云层.美国地质勘探局提供Landsat 8云掩膜产品, 该产品包含于影像的质量评估波段(quality assessment band, QAB)中[10], 由ACCA、See-5、Cirrus和AT-ACCA算法计算得到[11], Shen等[12]已在实际研究使用该产品.
针对下垫面和云混淆的问题, 目前主要有以下2种解决方法：1)针对某一特定下垫面的云识别方法[13]；2)将下垫面分类后, 再针对某一具体类型提出相应的方法[14].这些方法虽然能够在一定程度上解决问题, 但它们或是只针对某一特定下垫面有效, 不具有普适性, 或是操作复杂, 不能满足快速云识别的要求.散点图具有直观、便于定量化研究等优点, 于敏等[15-18]已将散点图应用于多种领域.
阈值的选取将影响最终的云识别效果, 常规的做法是利用经验或统计数据设置一固定阈值[19-21].固定阈值法能够识别包括薄云在内的云像元, 但在数据源较少、时间要求较高的情况下, 薄云区影像具有相当大的应用价值.阈值会随着影像的季节、地区气候等因素的不同而发生变化, 固定阈值法的普适性较差.
针对上述问题, 本文提出基于Landsat 8影像的云识别方法——SARM(spectral area ratio method).首先利用云和不同地物在可见光和热红外波段辐射率的差异, 以波谱面积表征像元总辐射强度, 计算得到基于像元的波谱面积比值(spectral area ratio, Sratio).结合波谱面积比值和归一化植被指数(normalized difference vegetation index, NDVI), 构建影像的散点图.针对固定阈值法影像利用率低、普适性差的不足, 本文对每景影像单独设置阈值, 采用高、中、低3种云识别置信区间.
1 方法简介 1.1 波谱面积比值的计算Landsat 8影像中不同地物类型在光谱空间中的分布特征存在明显差异[22], 分别对影像(p118_r39, 2014-10-3)中云、植被、裸地和水体4种类型进行采样, 选择均匀分布的像元样本, 每类样本数量均超过25个.分别对采样对象在光谱空间中的分布特征进行统计分析, 得到不同覆盖类型的辐射光谱曲线, 如图 1所示.图中, R为辐射率，λ为波长，实线表示本文方法使用的波谱曲线, 虚线表示未使用的波谱曲线.从图 1可以看出, 在可见光和近红外波段(波段1~9), 云的辐射率远高于其他地物, 但是随着波长的增加, 云的辐射率会逐渐降低, 到热红外波段(波段10、11), 云的辐射率低于其他地物；裸地的辐射率会随着波长的增加而缓慢降低, 但在热红外波段, 裸地的辐射率高于其他地物；由于受叶绿素和含水量的影响, 植被在红波段(0.64~0.67 μm)和卷云波段(1.36~1.38 μm)附近有明显的吸收谷, 中间是一个反射峰；在可见光和热红外波段, 水体和植被的辐射率近似.图 1中, 每类地物在波谱空间中都有连续的特征波谱曲线, 这是利用光谱区分地物的物理基础.传统云识别方法[5-10]大都利用几个波段的差值或比值来强化云和下垫面的差异, 但当下垫面的光谱和云类似时, 离散的波段组合很难将两者区分.
|
图 1 Landsat 8影像中云、薄云、裸地、植被和水体的辐射波谱特征 Fig. 1 Radiation spectral characteristics of clouds, thin clouds, bare land, vegetation and water bodies in Landsat 8 images |
针对这一问题, 引入波谱面积比值的方法：地物在可见光、近红外波段处的波谱面积与热红外波段处波谱面积的比值(此处的面积指的是在笛卡尔坐标系上特征波谱曲线与横轴所围限的范围).该方法的优势如下：波谱面积能够更好地贴近地物的实际波谱曲线, 面积比值能够更好地强化云和具有类似光谱特征地物的差异, 达到较好的区分效果.具体操作如下：Landsat 8拥有11个波段, 使用可见光到近红外(波段1~5)和热红外(波段10、11)这7个波段.将这7个波段划分为两部分：A部分包括5个波段(0.43~0.88 μm)；B部分包括2个波段(10.60~12.51 μm).A、B两部分波谱面积及波谱面积比值如下：
| $ {S_{\rm{A}}} = \sum\limits_{i = 1}^4 {\frac{1}{2}} \left( {{R_i} + {R_{i + 1}}} \right)\left( {{\lambda _{i + 1}} - {\lambda _i}} \right), $ | (1) |
| $ {S_{\rm{B}}} = \frac{1}{2}({R_{10}} + {R_{11}})({\lambda _{11}} - {\lambda _{10}}), $ | (2) |
| $ {S_{{\rm{ratio}}}} = \ln \frac{{{S_{\rm{A}}}}}{{{S_{\rm{B}}}}}. $ | (3) |
式中：Ri为第i波段的辐射率, λi为第i波段的中心波长, SA为前5个波段的波谱面积, SB为后2个波段的波谱面积, Sratio为波谱面积比值.
1.2 NDVI的选取及散点图的构建当下垫面类型单一时(只含陆地或水体), 可以构建波谱面积比值频率的分布曲线, 选择曲线的最低点作为区分云和下垫面的阈值[23].当下垫面条件复杂时, 由于水、陆和云的光谱特征差异较大, 波谱面积比值频率分布曲线上会同时存在多个峰(水、陆和云), 不同下垫面可能具有相似的Sratio而混合在一起, 频率分布曲线图存在无法选取阈值的问题.因此, 需引入另一变量来构建散点图, 将各地物区分开.归一化植被指数(NDVI)能够有效地区分各类地物, 选择Sratio为横轴变量, NDVI为纵轴变量, 构建散点图.NDVI的定义式为
| $ {\rm{NDVI}} = \frac{{{\rm{NIR}} - R}}{{{\rm{NIR + }}R}}. $ | (4) |
式中：NIR和R分别代表近红外和红波段, 负值表示地面覆盖为水、雪等；0表示有岩石、裸土和云等；正值表示有植被覆盖.影像(见图 2(a))中每个像元都能够计算得到相应的Sratio和NDVI, 将各像元计算结果投影到散点图中, 得到图 2(b)所示的效果.图中, P为像元数.
|
图 2 不同地物类型在散点图中的分布(p118_r39, 2014-10-3) Fig. 2 Distribution of different land cover types in scatter plot(p118_r39, 2014-10-3) |
图 2(b)中, 陆-云过渡带(陆上薄云、云边缘等)不仅具有云的光谱特征, 而且混合了陆地的辐射特征, 因此呈条带状分布于陆地和云之间；水-云过渡带(水上薄云、云边缘等)混合了水体和云的信息, 分布于水体和云之间.
1.3 阈值的确定由于水体和陆地的光谱特征不同, 两者的Sratio和NDVI具有明显的差异.水-云和陆-云过渡带混合了水陆的辐射信息, 两者在散点图中的分布不一致, 因此需要针对这两种不同的过渡带, 分别设置判别阈值.阈值的确定以下垫面是陆地时为例, 具体过程如下.
1) 在NDVI-Sratio散点图(见图 3(a))上沿陆地向云方向作直线.
|
图 3 阈值点的确定 Fig. 3 Determination of thresholds |
2) 对直线进行采样计算, 得到陆-云像元频率分布曲线(见图 3(b)).像元频率分布曲线为典型的双峰分布, 左侧峰为陆地像元投影, 右侧峰为云像元投影, 中间谷为含陆地信息的云像元投影(陆上薄云、云边缘等), 即陆-云过渡带.
3) 在像元频率曲线(见图 3(b))上存在一最低点BL, 左侧峰与谷之间存在一分界点AL(可以用两切线求交法确定分界点位置), 右侧峰与谷之间存在一分界点CL.AL、BL、CL 3个阈值点分别对应Landsat 8 QA波段中云低、中、高3种置信区间.
同理可得水-云像元频率分布曲线(见图 3(c))以及AW、BW、CW 3个阈值点.设置阈值点的目的如下.a)低置信区间阈值：当云对研究精度的影响较大时(如反演地面温度)时, 可以选择AL、AW两阈值点, 尽可能去除薄云对反演结果的影响, 保证研究精度.b)高置信区间阈值：当云对研究精度的影响较小(如土地分类)时, 可以选择CL、CW两阈值点, 尽可能保留更多的有效像元, 提高影像利用率.c)中置信区间阈值：通常情况下, 可以选择BL、BW两阈值点.确定阈值点后, 过阈值点可以刻画两条判别函数线, 函数线的右侧为剔除的云像元.
如图 4所示为采用高、低两种不同置信区间的云识别效果.图 4(d)中的判别函数线参考CL、CW点, 采用高置信区间, 去云影像(见图 4(b))中只包括厚云；图 4(g)中的判别函数线参考AL、AW点, 采用低置信区间, 去云影像(见图 4(e))中包括厚云、水/陆-云过渡带上的薄云和云边缘.通过实验可以看出, 利用提出的SARM方法可以动态调整阈值, 满足不同研究目的对云识别的需求.
|
图 4 SARM方法高、低置信区间云识别效果图(p118_r39, 2014-10-3) Fig. 4 High and low confidence intervals of cloud detection effect diagrams (p118_r39, 2014-10-3) |
选择USGS网站(https://www.usgs.gov/)上3景不同地区的Landsat 8影像作为实验数据：热带季风气候区(116/50)、温带大陆性气候区(13/32)和高原山地气候区(138/39), 具体信息如表 1所示.
|
表 1 Landsat 8影像基本信息 Table 1 Basic information of Landsat 8 images |
为了使结果更客观, SARM方法在3景影像中统一采用中置信区间阈值.图 5(a)中下垫面以植被、水体、熔岩和城市为主, 由于熔岩拥有较高的反射率, 通常容易被误判为云.在散点图中, 熔岩呈团块状分布于中央, 云呈带状分布于右侧；图 5(b)中城市不仅拥有高反射率, 而且表现为低频率、弱纹理特征, 通常容易和云发生混淆.在散点图中, 城市和云被有效地区分开；如图 5(c)所示为念青唐古拉山脉区域, 山脉上覆盖大量冰川.由于冰川和云在可见光波段拥有相似的光谱特性, 很容易被误判为云.在散点图中, 虽然云和冰川的Sratio相似, 但因云的NDVI更高, 两者得到了有效的区分.通过实验可以看出：在面对异物同谱的情况时, SARM方法能够取得令人满意的云识别效果.
|
图 5 SARM方法云识别的效果 Fig. 5 Cloud detection effects of SARM method |
选取Landsat 8影像QA波段的云掩膜图像[11](以下简称QA波段法)、ACCA方法[5-6]和Cirrus波段相结合的ACCA+Cirrus方法以及Fmask云识别方法[9]作为参照, 选取3景Landsat 8影像进行方法对比(影像信息如表 1所示), 从目视效果和识别精度上对4种方法进行比较和评价.
3.1 目视效果图 6(a)中, QA波段法和ACCA+Cirrus方法将熔岩误判为云, Fmask方法只识别了厚云, 不能识别影像中大量的积云.图 6(b)中, QA波段法将高反射率的城市道路和海岸线误判为云, Fmask和ACCA+Cirrus方法存在对云边缘的误判.图 6(c)中, QA波段法、Fmask及ACCA+Cirrus方法都把冰川误判为云.与其他方法相比, SARM方法的目视效果最好, 在3景实验影像中都表现了较强的抗噪声能力.
|
图 6 4种云识别方法的视觉对比 Fig. 6 Comparison of vision for four cloud detection methods |
为了定量比较各方法的识别精度, 在3景实验影像中选取具有代表性的9个区域对各方法进行精度评价, 每个评价区域为10 000个像元.由于无法获取完全真实的云图, 依据经验[24]可知, 真实云图由手工勾画获取.使用混淆矩阵的精度评价方法[8].具体计算公式为
| $ {A_{\rm{P}}} = \frac{{{T_{\rm{C}}}}}{{{T_{\rm{C}}} + {T_{\rm{F}}}}}, $ | (5) |
| $ {A_{\rm{U}}} = \frac{{{T_{\rm{C}}}}}{{{T_{\rm{C}}} + {F_{\rm{T}}}}}, $ | (6) |
| $ {A_{\rm{O}}} = \frac{{{T_{\rm{C}}} + {T_{\rm{G}}}}}{{{T_{\rm{P}}}}}, $ | (7) |
式中：AP为生产者精度, TC为被正确分类的云像元数, TF为云判定为非云的像元个数, AU为用户精度, FT为非云判定为云的像元个数, AO为总体精度, TG为被正确分类的地表像元数, TP为总像元数.各方法的生产者精度、用户精度和总体精度如表 2所示.
|
|
表 2 4种云识别方法的精度对比 Table 2 Comparison of accuracy for four cloud detection methods |
结合表 2分析可知, QA波段法、Fmask及ACCA+Cirrus方法难以兼顾生产者精度和用户精度, 漏判和误判较严重.QA波段法在9个评价区的生产者精度均高于80%, 但误判率较高, 用户精度偏低.Fmask法在4个评价区中的生产者精度达到100%, 但由于对冰川和人造地物存在许多误判, 用户精度最低仅为9.7%.在评价区a, Fmask方法在该区域存在大量漏判, 生产者精度仅为20.5%.同样, ACCA+Cirrus在评价区g和h的生产者精度达到100%, 但因误判较多, 用户精度分别仅为28%和12.2%.如图 7所示为4种方法的平均精度Aa.
|
图 7 识别方法的平均精度 Fig. 7 Average accuracy of four cloud detection methods |
(1) Landsat 8影像可见光到近红外波段(波段1~5)中, 云的辐射率远高于其他地物；在热红外波段(波段10、11), 云的辐射率比其他地物低.根据这一特征, 提出基于像元的波谱面积比值Sratio.该比值能够准确刻画不同地物在波谱空间中的整体特征, 有效地排除高反射率城镇、熔岩和海岸线等干扰.
(2) 构建基于归一化植被指数(NDVI)和波谱面积比值的散点图.该散点图能够适用于多种下垫面条件下的云识别, 并且操作简便, 效率较高.
(3) 采用高、中、低3种云识别置信区间, 可以通过调整阈值来满足不同研究目的对云识别的需求, 提高了影像利用率, 弥补了传统方法阈值固定、普适性差的不足.
(4) 在3景不同地区验证了SARM方法, 并和其他3种云识别方法在9个不同实验区进行精度评价分析.结果表明：SARM方法在目视效果和识别精度上都表现了较强的优越性, 总体精度提高10%左右.
| [1] |
HAGOLLE O, HUC M, PASCUAL D V, et al. A multi-temporal method for cloud detection, applied to FORMOSAT-2, VENμS, LANDSAT and SENTINEL-2 images[J]. Remote Sensing of Environment, 2010, 114(8): 1747-1755. DOI:10.1016/j.rse.2010.03.002 |
| [2] |
GOODWIN N R, COLLETT L J, DENHAM R J, et al. Cloud and cloud shadow screening across Queensland, Australia:an automated method for Landsat TM/ETM plus time series[J]. Remote Sensing of Environment, 2013, 134: 50-65. DOI:10.1016/j.rse.2013.02.019 |
| [3] |
JIN S, HOMER C, YANG L, et al. Automated cloud and shadow detection and filling using two-date Landsa-t imagery in the USA[J]. International Journal of Remote Sensing, 2013, 34(5): 1540-1560. DOI:10.1080/01431161.2012.720045 |
| [4] |
ZHENG L J, WU Y, YU T, et al. Object-based cloud detection of multitemporal high-resolution stationary satellite images[J]. Optical Engineering, 2017, 56(7): 73103. DOI:10.1117/1.OE.56.7.073103 |
| [5] |
IRISH R R. Landsat 7 automatic cloud cover assessment[C]//Proceedings of SPIE. Orlando: SPIE, 2000: 348-355. http://www.researchgate.net/publication/228946885_Landsat_7_automatic_cloud_cover_assessment
|
| [6] |
IRISH R R, BARKER J L, GOWARD S N, et al. Characterization of the Landsat-7 ETM+ automated cloud cover assessment (ACCA) algorithm[J]. Photogrammetric Engineering and Remote Sensing, 2006, 72(10): 1179-1188. DOI:10.14358/PERS.72.10.1179 |
| [7] |
ZHU Z, WOODCOCK C E. Object-based cloud and cloud shadow detection in Landsat imagery[J]. Remote Sensing of Environment, 2012, 118: 83-94. DOI:10.1016/j.rse.2011.10.028 |
| [8] |
ZHU Z, WANG S X, WOODCOCK C E. Improvement and expansion of the Fmask algorithm:cloud, cloud shadow, and snow detection for Landsats 4-7, 8, and Sentinel 2 images[J]. Remote Sensing of Environment, 2015, 159: 269-277. |
| [9] |
FOGA S, SCARAMUZZA P L, GUO S, et al. Cloud detection algorithm comparison and validation for operational Landsat data products[J]. Remote Sensing of Environment, 2017, 194: 379-390. DOI:10.1016/j.rse.2017.03.026 |
| [10] |
ZANTER K. Landsat 8(L8) data users handbook[M]. 2th ed. Sioux Falls, South Dakota: EROS, 2016, 47-54.
|
| [11] |
SCARAMUZZA P L, BOUCHARD M A, DWYER J L. Development of the Landsat data continuity mission cloud-cover assessment algorithms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(4): 1140-1154. DOI:10.1109/TGRS.2011.2164087 |
| [12] |
SHEN Y, WANG Y, LV H, et al. Removal of thin clouds using cirrus and QA bands of Landsat-8[J]. Photogrammetric Engineering and Remote Sensing, 2015, 81(9): 721-731. DOI:10.14358/PERS.81.9.721 |
| [13] |
蒋嫚嫚, 邵振峰. 采用主成分分析的改进云检测算法[J]. 测绘科学, 2015, 40(2): 150-154. JIANG Man-man, SHAO Zhen-feng. Advanced algorithm of PCA-based Fmask cloud detection[J]. Science of Surveying and Mapping, 2015, 40(2): 150-154. |
| [14] |
夏浪, 毛克彪, 孙知文, 等. 针对NPP VIIRS数据的云检测方法研究[J]. 中国环境科学, 2014, 34(3): 574-580. XIA Lang, MAO Ke-biao, SUN Zhi-wen, et al. Cloud detection application on NPP VIIRS[J]. China Environmental Science, 2014, 34(3): 574-580. |
| [15] |
于敏, 程明虎, 刘辉. 地表温度-归一化植被指数特征空间干旱监测方法的改进及应用研究[J]. 气象学报, 2011, 5(69): 922-931. YU Min, CHENG Ming-hu, LIU Hui. An improvement of the land surface temperature-NDVI space drought monitoring method and its applications[J]. Acta Meteorologica Sinica, 2011, 5(69): 922-931. |
| [16] |
王娇, 丁建丽, 袁泽, 等. 基于Ts-NDVI特征空间的绿洲土壤水分监测算法改进[J]. 中国沙漠, 2016, 6(36): 1606-1612. WANG Jiao, DING Jian-li, YUAN Ze, et al. Improvement and comparison of soil moisture monitoring algorithm in oasis based on Ts-NDVI feature space[J]. Journal of Desert Research, 2016, 6(36): 1606-1612. |
| [17] |
PAN P P, CHEN G Y, SARUTA K, et al. Snow cover detection based on two-dimensional scatter plots from MODIS imagery data[J]. Journal of Applied Remote Sensing, 2015, 9(1): 096083. DOI:10.1117/1.JRS.9.096083 |
| [18] |
CAO X M, FENG Y M, WANG J L. Remote sensing monitoring the spatio-temporal changes of aridification in the Mongolian Plateau based on the general Ts-N-DVI space, 1981-2012[J]. Journal of Earth System Science, 2017, 126(4): 58. DOI:10.1007/s12040-017-0835-x |
| [19] |
ACKERMAN S A, STRABALA K I, MENZEL W P, et al. Discriminating clear sky from clouds with MODI-S[J]. Journal of Geophysical Research:Atmospheres, 1998, 103(24): 32141-32157. |
| [20] |
CHEN P Y, SRINIVASAN R, FEDOSEJEVS G, et al. An automated cloud detection method for daily NOAA-14 AVHRR data for Texas, USA[J]. International Journal of Remote Sensing, 2002, 23(15): 2939-2950. DOI:10.1080/01431160110075631 |
| [21] |
LIU Y H, KEY J R, FREY R A, et al. Nighttime polar cloud detection with MODIS[J]. Remote Sensing of Environment, 2004, 92(2): 181-194. DOI:10.1016/j.rse.2004.06.004 |
| [22] |
宋瑞祥, 张庆国, 于海敬, 等. 遥感数据的城市不透水面估算及增温效应[J]. 浙江大学学报:工学版, 2017, 51(5): 1051-1056. SONG Rui-xiang, ZHANG Qing-guo, YU Hai-jing, et al. Estimations to impervious surface and their effects of warming for city using remote sensing data[J]. Journal of Zhejiang University:Engineering Science, 2017, 51(5): 1051-1056. |
| [23] |
GUO F, SHEN X H, ZOU L J, et al. Cloud detection method based on spectral area ratios in MODIS data[J]. Canadian Journal of Remote Sensing, 2015, 41(6): 561-576. |
| [24] |
SUN L, MI X T, WEI J, et al. A cloud detection algorithm-generating method for remote sensing data at visible to short-wave infrared wavelengths[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2017, 124: 70-88. DOI:10.1016/j.isprsjprs.2016.12.005 |