高分辨率遥感图像厚云检测算法

郭 秋 ， 李大成 ， 刘贺春

(1.晋中职业技术学院，山西 晋中 030600； 2.太原理工大学矿业工程学院，太原 030024； 3．山西华晋岩土工程勘察有限公司，太原 030021)

高分辨率遥感图像厚云检测算法

郭 秋1， 李大成2， 刘贺春3

(1.晋中职业技术学院，山西 晋中 030600； 2.太原理工大学矿业工程学院，太原 030024； 3．山西华晋岩土工程勘察有限公司，太原 030021)

针对高空间分辨率多光谱图像，设计了一种基于缨帽变换、数据正则化以及多尺度光谱分析技术的多尺度循环检测算法。实验结果表明：对于Landsat-7 ETM+图像，其检测精度与基于热波段的自动化云量估计算法相近，但优于传统的直接分类方法(最大似然分类与ISODATA分类)。由于该算法无需借助热波段及其他辅助数据，因而在实现高分辨率图像厚云掩膜的高效与高精度提取方面具有很高的应用潜力。

厚云检测；缨帽变换；多尺度光谱分析；自动化云量估计

0 引言

在基于遥感技术的对地观测研究中，厚云及其阴影不仅限制了图像中的可用信息量,还会严重影响图像分类、大气校正、气溶胶反演以及生物物理参数的估算等过程，而去除厚云污染的首要任务则是高效、精确地检测出图像中的厚云像元。

鉴于厚云的高反射和低温特性，常用的检测策略是依据基于多光谱特征的经验阈值来过滤出云像元，如针对MODIS数据的云掩膜算法需利用基于辐射率或亮度温度的一系列经验阈值[1-2]；R. R. Irish在2000年针对Landsat-7 ETM+数据发展了一种基于热波段数据的自动化云量估计算法(automated cloud coverage assessment，ACCA)[3]，该算法因检测精度高已被美国地质调查研究局用来估算Landsat TM/ETM+图像的云量。H.Choi等[4]和G.R.Watmough等[5]拓展了ACCA算法，使其对云边界、卷云以及海冰上的厚云检测更为有效。

通过分析不同时相下同一传感器图像的厚云分布特征[6-7]或相同时相下不同传感器数据的一致性特征[8]，多图像分析方法可以获得较佳的检测效果，如最大值合成法[9-10]、基于时间序列分析的厚云检测算法[11]等。但此类方法一方面要求待处理图像与辅助数据间应进行精确的几何配准与辐射定标(此情形下，辐射正则化方法可以提供较为有效的相对辐射校正处理[12])；另一方面，由于需要借助辅助数据，因而无法处理研究区存在持续云覆盖的情形。

上述厚云检测方法要求待处理的传感器数据应包含热波段图像或者需要覆盖同一研究区且分辨率相近的辅助遥感数据，但高空间分辨率影像由于成像系统的限制一般都缺少热红外谱段，并且因数据采集周期的影响而难以搜集到可用的辅助参考影像。当前针对此问题的主要解决途径是借助基于统计学思想的厚云检测 策略，其中以分类方法应用最为广泛，如基于神经网络分类器的方法[13]、基于模糊逻辑学的方法[14]、最大似然估计法[15]以及基于支持向量机的分类器方法[16]等。虽然监督分类方法比非监督分类方法的精度更高，但它需人工参与且待处理数据应满足特定的数学模型，如贝叶斯估计和最大似然估计均要求实际数据满足正态分布；而非监督分类法虽然具有自动化程度高的特点，但也因缺少专家知识而造成检测的精度水平较低，如基于迭代自组织数据分析算法的分类方法(iterative self-organizing data analysis techniques algorithm，ISODATA)[17]。

综上所述，对于高空间分辨率遥感影像，基于分类策略的厚云检测方法在精度与效率上无法充分兼顾。本研究通过综合缨帽(tasseled cap，TC)变换[18]、数据正则化以及多尺度光谱分析技术设计一种通用性较强的多尺度循环检测算法，以期在保持高检测效率的同时提升厚云的检测精度，并评估该方法的精度。

1 研究区概况

研究区位于广东省西南部的阳春市附近，中心位置经纬度为东经111°40′54″和北纬21°40′02″。该区地形以山地丘陵为主，地表覆盖类型包括植被、水体、居民地、农田、山地等，有利于验证厚云检测方法在复杂地形条件下的处理效果。为便于对比分析，采用Landsat-7 ETM+图像作为实验数据，轨道编号123/45(列/行)，采集时间为2001年8月7日，在此时相下影像中的云类别以中等尺寸的厚云以及少量的碎云/薄云为主。

根据不同厚云检测方法对数据预处理的要求，首先对覆盖研究区的ETM+ 1～7个波段的数据进行辐射校正，将原始的DN(digital number)值转换为TOA(top of atmosphere)辐射率。然后，非热波段(波段1～5，7)数据转换为TOA反射率，热波段(波段6)数据转换为星上温度[3]。

2 基于多尺度循环的厚云检测算法原理

所设计的多尺度循环检测算法的整体处理流程如图1所示。

图1 基于多尺度循环的厚云检测流程

2.1 缨帽变换

对待处理的高分辨率图像进行缨帽变换，缨帽变换是根据土壤、植被等光谱信息在多维光谱空间中的信息分布结构对图像所做的经验性线性正交变换，其变换模型如下：

(1)

式中：x和u分别是变换前后多光谱图像的信号矢量；R是由多个相互正交的单位矢量组成的单位矩阵；r是为避免出现负值而设置的补偿值。变换后所得的信号矢量u中，前3个分量图像分别称为亮度、绿度和湿度，且它们集中了原始影像中绝大部分的信息量，因而缨帽变换后得到的光谱空间能够为地表特征提供更好的区分性与解析性。此外，依赖于传感器、季相条件以及地理生态分区的缨帽转换系数是进行自动化与高精度缨帽变换的必要条件，本研究针对待处理的研究区数据选取了多幅符合上述要求的ETM+数据，并进一步提取了基于该研究区的ETM+图像在相同季相条件下的缨帽转换系数，以保证参与计算的转换系数对于研究区的ETM+数据具有代表性。

由于厚云在来自不同传感器、不同地理区域或不同时相的图像上所表现出的光谱分布特征存在差异，直接分析厚云在传统的缨帽空间中的光谱分布特征一般不具有代表性，如亮度分量中就存在较高的噪声干扰。为强化缨帽空间对数据的光谱解析能力并且提高算法的通用性，这里采用一种数据正则化方法将原始的缨帽三分量转换为光谱结构特征更强且噪声干扰较少的正则化分量，并利用转换后的正则化分量图像进行后续的云检测过程，其转换过程可表达为：

(2)

(3)

(4)

式中：BI，GI,WI分别是经缨帽变换得到的原始三分量(亮度、绿度与湿度)数据；Bμ，Gμ,Wμ分别为原始三分量的图像均值；Bσ，Gσ,Wσ分别为原始三分量的图像标准差；BR，GR,WR是经正则化后的亮度、绿度和湿度分量。

2.2 ISODATA非监督分类

利用ISODATA非监督分类[17]算法处理经缨帽变换与数据正则化后的前3个缨帽分量数据，并以此作为多尺度循环检测的起始步骤。通过结合最小距离技术与迭代策略，ISODATA分类方法可以在未预先指定类别数的情况下得到满足分类条件的最优类别数。这里利用基于类别尺度的光谱分布特征将所有得到的N个类别K1，K2，…，KN初步划分为3组，即“云类别”、“待定类别”以及“非云类别”。鉴于云的高反射特性，这里采用缨帽变换后的亮度分量BR图像作为基于类别尺度的光谱分析的依据，并使用两个基于目标传感器地物光谱分布特征的经验阈值Tb1和Tb2来区分出上述3个组别。

(5)

(6)

(7)

式中：[BR(Ki)]μ是类别Ki(i∈[1,N])在亮度分量BR图像上的光谱均值。

若类别Ki能够满足式(5)、(6)或(7)，则分别划分为“云类别”、“待定类别”或“非云类别”，其中非云类别中的所有像元将被直接剔除而不参与后续检测流程。此外，为剔除“云类别”中那些与云块相似的高亮地表特征(如沙滩、河床以及裸地等)，这里利用基于像元尺度的过滤器来进一步过滤出“云类别”中的非云像元。具体而言，本算法通过分析“云类别”中的目标像元在由亮度分量BR、绿度分量GR以及湿度分量WR所构成的光谱空间中的分布特征来设置作为过滤参数的3个约束阈值TBR，TGR与TWR。

(8)

式中：BR(Px,y)，GR(Px,y)以及WR(Px,y)分别是图像位置(x,y)上的待定像元P在数据正则化后的亮度、绿度以及湿度分量上的光谱值。若“云类别”中的目标像元满足式(8)，则该像元被标记为最终的云像元，否则将和“待定类别”一起被过滤到“候选像元集”并加入后续的多尺度循环检测过程进行进一步识别。

2.3 多尺度循环检测

多尺度循环检测的流程与上述过程基本一致，但其中基于类别和基于像元尺度的检测参数将进行自适应调整，调整方式既可以在原始参数的基础上按一定比例调节，也可以通过分析“候选像元集”中所有像元所构成的光谱直方图来设定新的检测参数。基于类别和基于像元的多尺度循环检测策略能通过上述方式来不断优化检测效果，同时可以实现云边界的拓展以及降低漏检率，最终实现厚云像元自动化的精确识别。而对于ISODATA分类结果中的“碎类”以及检测结果中存在的“碎云”，则可以采用形态学方法来消除(如开操作处理)。

3 实验结果与分析

ACCA算法是当前针对ETM+数据精度最高的厚云检测方法，这里以其作为精度验证的比较依据。图2是选取的3个研究区子图及4种厚云检测方法(ACCA算法、最大似然分类、ISODATA分类以及多尺度循环检测算法)的处理结果。其中，多尺度循环检测算法的运行参数为：Tb1=2.0，Tb2=4.0，TBR=1.5，TGR=1.0，TWR=1.5，循环次数为1。

图2 研究区子图及4种厚云检测方法的云掩膜结果

图3是以ACCA算法为比较依据所展示的最大似然分类、ISODATA分类以及多尺度循环检测算法所得到的不同尺寸水平下的云块数量及其一致性对比图，其中横坐标代表云块的尺寸水平，左侧纵坐标代表所检测出的云块数量，右侧纵坐标为该方法与作为对比依据的ACCA算法的检测结果的总体一致性。图中的研究区子图1、子图2、子图3分别是与图2对应的3个研究区子图。

从图2可知：多尺度循环检测算法的云掩膜结果在整体上具有比ACCA算法和传统的直接分类方法(最大似然与ISODATA分类)更为拓展的云边界，但更容易产生一些碎小的误检点，这表明该算法对于厚云(尤其是碎云)的敏感性较高，并且可通过形态学操作得到有效改善。从图3中3个研究区子图的对比结果可以发现：多尺度循环检测算法在不同尺寸水平下所检测到的云块数量更为接近ACCA算法，但在检测结果的总体一致性上要低于传统的直接分类方法(最大似然分类与ISODATA分类)。

表1统计了3个研究区子图下最大似然分类、ISODATA分类、多尺度循环检测算法与作为对比依据的ACCA算法的检测结果的一致性详细信息。从中可知：一方面，传统的直接分类方法在误检率与总体一致性上要稍优于本研究的算法，这与图2、图3对比分析后的结果一致；但另一方面，多尺度循环检测算法在漏检率与Kappa系数这两项对比数据上占有较大优势。

图3 最大似然分类、ISODATA分类、多尺度循环检测算法与参考方法在不同云块尺寸下的一致性

综上，多尺度循环检测算法与ACCA算法对于不同尺寸云块的检测结果差异(对应于误检率)主要来自于小块云团或碎云(小于50个像元，见图3)，两者在有效厚云像元的检测上具有很高的一致性，这主要体现在表1中多尺度循环检测算法的低漏检率(小于4%)与高Kappa系数(大于92%)，即该算法与ACCA算法的实际检测结果更为接近；与ACCA算法相比，传统的直接分类方法所检测到的有效云像元数量较少，因而具有更高的总体一致性，但它们的高漏检率说明传统的直接分类方法对于有效厚云像元的检测效果较差(表1)。

表1 不同厚云检测算法与ACCA算法的一致性统计信息表

4 结论

1)对于Landsat ETM+图像，所提出的多尺度循环检测算法的精度水平与ACCA算法相近，但优于传统的直接分类方法(以最大似然和ISODATA分类为例)；此外，该算法还可通过设置循环参数不断优化检测结果。

2)多尺度循环检测算法在一定程度上受缨帽变换对光谱空间重构能力的影响，而转换系数的提取又要依赖于研究区的地理生态和季相特征，因此，在实际应用中应针对研究区所处的区域生态和季相因子来统计出可靠的缨帽转换系数。

3)多尺度循环检测算法既不需要热波段数据，也不需要其他的辅助参考数据。当研究区具有可靠的缨帽转换系数时，该算法能够实现高度的自动化处理，因而可以对常用的高分辨率遥感数据进行快速而精确的厚云检测处理(如SPOT[18]、QuickBird等卫星)。

[1] Frey R A,Ackerman S A,Liu Y H,etal.Cloud Detection with MODIS.Part I:Improvements in the MODIS Cloud Mask for Collection 5[J].Journal of Atmospheric and Oceanic Technology,2008,25(7):1057-1072.

[2] Liu R,Liu Y.Generation of New Cloud Masks from MODIS Land Surface Reflectance Products[J].Remote Sensing of Environment,2013,133(12):21-37.

[3] Irish R R,Barker J L,Samuel N G,etal.Characterization of the Landsat-7 ETM+ Automated Cloud-cover Assessment (ACCA) Algorithm[J].Photogrammetric Engineering & Remote Sensing,2006,72(10):1179-1188.

[4] Choi H,Bindschadler R.Cloud Detection in Landsat Imagery of Ice Sheets Using Shadow Matching Technique and Automatic Normalized Difference Snow Index Threshold Value Decision[J].Remote Sensing of Environment,2004,91(2):237-242.

[5] Watmough G R,Atkinson P M,Hutton C W.A Combined Spectral and Object-based Approach to Transparent Cloud Removal in An Operational Setting for Landsat ETM+ [J].International Journal of Applied Earth Observation and Geoinformation,2011,13(2):220-227.

[6] 葛艳琴,郭丽丽.Landsat 7影像厚云及其阴影的自动检测方法[J].中国科技博览,2012(27):610-610.

[7] Derrien M,Gleau H L.Improvement of Cloud Detection near Sunrise and Sunset by Temporal-Differencing and Region-Growing Techniques with Real-Time SEVIRI[J].International Journal of Remote Sensing,2010,31(7):1765-1780.

[8] Fernando S,Pieter K,Peter S,etal.A Cloud Mask Methodology for High Resolution Remote Sensing Data Combining Information from High and Medium Resolution Optical Sensors[J].ISPRS Journal of Photogrammetry and Remote Sensing,2011,66(5):588-596.

[9] Holben B N.Characteristics of Maximum-Value Composite Images from Temporal AVHRR Data[J].International Journal of Remote Sensing,1986,7(11):1417-1434.

[10] Stowe L L,Mcclain E P,Carey R,etal.Global Distribution of Cloud Cover Derived from NOAA/AVHRR Operational Satellite Data[J].Advances in Space Research,1991,11(3):51-54.

[11] Hagolle O,Huc M,Pascual D V,etal.A Multi-temporal Method for Cloud Detection,Applied to Formosat-2,VENUS,Landsat and Sentinel-2 Images[J].Remote Sensing of Environment,2010,114(8):1747-1755.

[12] Caselles V.An Alternative Simple Approach to Estimate Atmospheric Correction in Multitemporal Studies[J].International Journal of Remote Sensing,1989,10(6):1127-1134.

[13] Turaga S C,Murray J F,Jain V,etal.Convolutional Networks Can Learn to Generate Affinity Graphs for Image Segmentation[J].Neural Computation,2010,22(2):511-538.

[14] Fernández A,Garcia S,Del J,etal.A Study of the Behaviour of Linguistic Fuzzy Rule Based Classification Systems in the Framework of Imbalanced Data-sets[J].Fuzzy Sets and Systems,2008,159(18):2378-2398.

[15] Strahler A H.The Use of Prior Probabilities in Maximum Likelihood Classification of Remotely Sensed Data[J].Remote Sensing of Environment,1980,10(2):135-163.

[16] Mountrakis G,Im J,Ogole C.Support Vector Machines in Remote Sensing:A Review[J].ISPRS Journal of Photogrammetry and Remote Sensing,2011,66(3):247-259.

[17] Melesse A M,Jordan J D.A Comparison of Fuzzy vs. Augmented-ISODATA Classification Algorithms for Cloud-Shadow Discrimination from Landsat Images[J].Photogrammetric Engineering & Remote Sensing,2002,68(9):905-911.

[18] Ivits E,Lamb A,Langar F,etal.Orthogonal Transformation of Segmented SPOT5 Images Seasonal and Geographical Dependence of the Tasselled Cap Parameters[J].Photogrammetric Engineering & Remote Sensing,2008,74(11):1351-1364.

Cloud Masking for Multi-spectral Remotely Sensed Imagery with High Resolution

Guo Qiu1, Li Dacheng2, Liu Hechun3

(1.JinzhongVocationalTechnicalInstitute,Jinzhong030600,China; 2.MiningEngineeringInstitute,TaiyuanUniversityofTechnology,Taiyuan030024,China; 3.ShanxiHuajinEngineeringReconnaissanceLtd.,Taiyuan030021,China)

As to multi-spectral images with high spatial resolution, an automated cloud-masking algorithm was designed based on the Tasseled Cap transformation, data normalization and multi-scale spectrum analysis technique. The preliminary experimental result showed that, for the Landsat-7 ETM plus imagery, the proposed algorithm had a similar precision level with the automated cloud coverage assessment algorithm which relied on thermal band, but higher precision could be expected than direct classification methods (maximum-likelihood and ISODATA). Since this algorithm needs no thermal band or other referenced data, it has high potential in retrieving cloud masks from high spatial resolution imagery effectively and accurately.

thick-cloud detection; tasseled cap transformation; multi-scale spectrum analyses; automated cloud coverage assessment

2015-08-27；

2016-06-28

国家自然科学基金青年项目(41501372)

郭秋(1980-)，女，山西太原市人，讲师，硕士，主要从事3S技术应用研究，(E-mail)guoqiu1980@126.com。

李大成(1983-)，男，山西太原市人，讲师，博士，主要从事遥感信息处理及地学应用研究，(E-mail)lidacheng@tyut.edu.cn 。

TP751

A

1003-2363(2016)04-0172-05