高分辨率遙感圖像厚云檢測算法

郭 秋 ， 李大成 ， 劉賀春

(1.晉中職業(yè)技術學院，山西 晉中 030600； 2.太原理工大學礦業(yè)工程學院，太原 030024； 3．山西華晉巖土工程勘察有限公司，太原 030021)

高分辨率遙感圖像厚云檢測算法

郭 秋1， 李大成2， 劉賀春3

(1.晉中職業(yè)技術學院，山西 晉中 030600； 2.太原理工大學礦業(yè)工程學院，太原 030024； 3．山西華晉巖土工程勘察有限公司，太原 030021)

針對高空間分辨率多光譜圖像，設計了一種基于纓帽變換、數(shù)據(jù)正則化以及多尺度光譜分析技術的多尺度循環(huán)檢測算法。實驗結果表明：對于Landsat-7 ETM+圖像，其檢測精度與基于熱波段的自動化云量估計算法相近，但優(yōu)于傳統(tǒng)的直接分類方法(最大似然分類與ISODATA分類)。由于該算法無需借助熱波段及其他輔助數(shù)據(jù)，因而在實現(xiàn)高分辨率圖像厚云掩膜的高效與高精度提取方面具有很高的應用潛力。

厚云檢測；纓帽變換；多尺度光譜分析；自動化云量估計

0 引言

在基于遙感技術的對地觀測研究中，厚云及其陰影不僅限制了圖像中的可用信息量,還會嚴重影響圖像分類、大氣校正、氣溶膠反演以及生物物理參數(shù)的估算等過程，而去除厚云污染的首要任務則是高效、精確地檢測出圖像中的厚云像元。

鑒于厚云的高反射和低溫特性，常用的檢測策略是依據(jù)基于多光譜特征的經(jīng)驗閾值來過濾出云像元，如針對MODIS數(shù)據(jù)的云掩膜算法需利用基于輻射率或亮度溫度的一系列經(jīng)驗閾值[1-2]；R. R. Irish在2000年針對Landsat-7 ETM+數(shù)據(jù)發(fā)展了一種基于熱波段數(shù)據(jù)的自動化云量估計算法(automated cloud coverage assessment，ACCA)[3]，該算法因檢測精度高已被美國地質(zhì)調(diào)查研究局用來估算Landsat TM/ETM+圖像的云量。H.Choi等[4]和G.R.Watmough等[5]拓展了ACCA算法，使其對云邊界、卷云以及海冰上的厚云檢測更為有效。

通過分析不同時相下同一傳感器圖像的厚云分布特征[6-7]或相同時相下不同傳感器數(shù)據(jù)的一致性特征[8]，多圖像分析方法可以獲得較佳的檢測效果，如最大值合成法[9-10]、基于時間序列分析的厚云檢測算法[11]等。但此類方法一方面要求待處理圖像與輔助數(shù)據(jù)間應進行精確的幾何配準與輻射定標(此情形下，輻射正則化方法可以提供較為有效的相對輻射校正處理[12])；另一方面，由于需要借助輔助數(shù)據(jù)，因而無法處理研究區(qū)存在持續(xù)云覆蓋的情形。

上述厚云檢測方法要求待處理的傳感器數(shù)據(jù)應包含熱波段圖像或者需要覆蓋同一研究區(qū)且分辨率相近的輔助遙感數(shù)據(jù)，但高空間分辨率影像由于成像系統(tǒng)的限制一般都缺少熱紅外譜段，并且因數(shù)據(jù)采集周期的影響而難以搜集到可用的輔助參考影像。當前針對此問題的主要解決途徑是借助基于統(tǒng)計學思想的厚云檢測 策略，其中以分類方法應用最為廣泛，如基于神經(jīng)網(wǎng)絡分類器的方法[13]、基于模糊邏輯學的方法[14]、最大似然估計法[15]以及基于支持向量機的分類器方法[16]等。雖然監(jiān)督分類方法比非監(jiān)督分類方法的精度更高，但它需人工參與且待處理數(shù)據(jù)應滿足特定的數(shù)學模型，如貝葉斯估計和最大似然估計均要求實際數(shù)據(jù)滿足正態(tài)分布；而非監(jiān)督分類法雖然具有自動化程度高的特點，但也因缺少專家知識而造成檢測的精度水平較低，如基于迭代自組織數(shù)據(jù)分析算法的分類方法(iterative self-organizing data analysis techniques algorithm，ISODATA)[17]。

綜上所述，對于高空間分辨率遙感影像，基于分類策略的厚云檢測方法在精度與效率上無法充分兼顧。本研究通過綜合纓帽(tasseled cap，TC)變換[18]、數(shù)據(jù)正則化以及多尺度光譜分析技術設計一種通用性較強的多尺度循環(huán)檢測算法，以期在保持高檢測效率的同時提升厚云的檢測精度，并評估該方法的精度。

1 研究區(qū)概況

研究區(qū)位于廣東省西南部的陽春市附近，中心位置經(jīng)緯度為東經(jīng)111°40′54″和北緯21°40′02″。該區(qū)地形以山地丘陵為主，地表覆蓋類型包括植被、水體、居民地、農(nóng)田、山地等，有利于驗證厚云檢測方法在復雜地形條件下的處理效果。為便于對比分析，采用Landsat-7 ETM+圖像作為實驗數(shù)據(jù)，軌道編號123/45(列/行)，采集時間為2001年8月7日，在此時相下影像中的云類別以中等尺寸的厚云以及少量的碎云/薄云為主。

根據(jù)不同厚云檢測方法對數(shù)據(jù)預處理的要求，首先對覆蓋研究區(qū)的ETM+ 1～7個波段的數(shù)據(jù)進行輻射校正，將原始的DN(digital number)值轉(zhuǎn)換為TOA(top of atmosphere)輻射率。然后，非熱波段(波段1～5，7)數(shù)據(jù)轉(zhuǎn)換為TOA反射率，熱波段(波段6)數(shù)據(jù)轉(zhuǎn)換為星上溫度[3]。

2 基于多尺度循環(huán)的厚云檢測算法原理

所設計的多尺度循環(huán)檢測算法的整體處理流程如圖1所示。

圖1 基于多尺度循環(huán)的厚云檢測流程

2.1 纓帽變換

對待處理的高分辨率圖像進行纓帽變換，纓帽變換是根據(jù)土壤、植被等光譜信息在多維光譜空間中的信息分布結構對圖像所做的經(jīng)驗性線性正交變換，其變換模型如下：

(1)

式中：x和u分別是變換前后多光譜圖像的信號矢量；R是由多個相互正交的單位矢量組成的單位矩陣；r是為避免出現(xiàn)負值而設置的補償值。變換后所得的信號矢量u中，前3個分量圖像分別稱為亮度、綠度和濕度，且它們集中了原始影像中絕大部分的信息量，因而纓帽變換后得到的光譜空間能夠為地表特征提供更好的區(qū)分性與解析性。此外，依賴于傳感器、季相條件以及地理生態(tài)分區(qū)的纓帽轉(zhuǎn)換系數(shù)是進行自動化與高精度纓帽變換的必要條件，本研究針對待處理的研究區(qū)數(shù)據(jù)選取了多幅符合上述要求的ETM+數(shù)據(jù)，并進一步提取了基于該研究區(qū)的ETM+圖像在相同季相條件下的纓帽轉(zhuǎn)換系數(shù)，以保證參與計算的轉(zhuǎn)換系數(shù)對于研究區(qū)的ETM+數(shù)據(jù)具有代表性。

由于厚云在來自不同傳感器、不同地理區(qū)域或不同時相的圖像上所表現(xiàn)出的光譜分布特征存在差異，直接分析厚云在傳統(tǒng)的纓帽空間中的光譜分布特征一般不具有代表性，如亮度分量中就存在較高的噪聲干擾。為強化纓帽空間對數(shù)據(jù)的光譜解析能力并且提高算法的通用性，這里采用一種數(shù)據(jù)正則化方法將原始的纓帽三分量轉(zhuǎn)換為光譜結構特征更強且噪聲干擾較少的正則化分量，并利用轉(zhuǎn)換后的正則化分量圖像進行后續(xù)的云檢測過程，其轉(zhuǎn)換過程可表達為：

(2)

(3)

(4)

式中：BI，GI,WI分別是經(jīng)纓帽變換得到的原始三分量(亮度、綠度與濕度)數(shù)據(jù)；Bμ，Gμ,Wμ分別為原始三分量的圖像均值；Bσ，Gσ,Wσ分別為原始三分量的圖像標準差；BR，GR,WR是經(jīng)正則化后的亮度、綠度和濕度分量。

2.2 ISODATA非監(jiān)督分類

利用ISODATA非監(jiān)督分類[17]算法處理經(jīng)纓帽變換與數(shù)據(jù)正則化后的前3個纓帽分量數(shù)據(jù)，并以此作為多尺度循環(huán)檢測的起始步驟。通過結合最小距離技術與迭代策略，ISODATA分類方法可以在未預先指定類別數(shù)的情況下得到滿足分類條件的最優(yōu)類別數(shù)。這里利用基于類別尺度的光譜分布特征將所有得到的N個類別K1，K2，…，KN初步劃分為3組，即“云類別”、“待定類別”以及“非云類別”。鑒于云的高反射特性，這里采用纓帽變換后的亮度分量BR圖像作為基于類別尺度的光譜分析的依據(jù)，并使用兩個基于目標傳感器地物光譜分布特征的經(jīng)驗閾值Tb1和Tb2來區(qū)分出上述3個組別。

(5)

(6)

(7)

式中：[BR(Ki)]μ是類別Ki(i∈[1,N])在亮度分量BR圖像上的光譜均值。

若類別Ki能夠滿足式(5)、(6)或(7)，則分別劃分為“云類別”、“待定類別”或“非云類別”，其中非云類別中的所有像元將被直接剔除而不參與后續(xù)檢測流程。此外，為剔除“云類別”中那些與云塊相似的高亮地表特征(如沙灘、河床以及裸地等)，這里利用基于像元尺度的過濾器來進一步過濾出“云類別”中的非云像元。具體而言，本算法通過分析“云類別”中的目標像元在由亮度分量BR、綠度分量GR以及濕度分量WR所構成的光譜空間中的分布特征來設置作為過濾參數(shù)的3個約束閾值TBR，TGR與TWR。

(8)

式中：BR(Px,y)，GR(Px,y)以及WR(Px,y)分別是圖像位置(x,y)上的待定像元P在數(shù)據(jù)正則化后的亮度、綠度以及濕度分量上的光譜值。若“云類別”中的目標像元滿足式(8)，則該像元被標記為最終的云像元，否則將和“待定類別”一起被過濾到“候選像元集”并加入后續(xù)的多尺度循環(huán)檢測過程進行進一步識別。

2.3 多尺度循環(huán)檢測

多尺度循環(huán)檢測的流程與上述過程基本一致，但其中基于類別和基于像元尺度的檢測參數(shù)將進行自適應調(diào)整，調(diào)整方式既可以在原始參數(shù)的基礎上按一定比例調(diào)節(jié)，也可以通過分析“候選像元集”中所有像元所構成的光譜直方圖來設定新的檢測參數(shù)?；陬悇e和基于像元的多尺度循環(huán)檢測策略能通過上述方式來不斷優(yōu)化檢測效果，同時可以實現(xiàn)云邊界的拓展以及降低漏檢率，最終實現(xiàn)厚云像元自動化的精確識別。而對于ISODATA分類結果中的“碎類”以及檢測結果中存在的“碎云”，則可以采用形態(tài)學方法來消除(如開操作處理)。

3 實驗結果與分析

ACCA算法是當前針對ETM+數(shù)據(jù)精度最高的厚云檢測方法，這里以其作為精度驗證的比較依據(jù)。圖2是選取的3個研究區(qū)子圖及4種厚云檢測方法(ACCA算法、最大似然分類、ISODATA分類以及多尺度循環(huán)檢測算法)的處理結果。其中，多尺度循環(huán)檢測算法的運行參數(shù)為：Tb1=2.0，Tb2=4.0，TBR=1.5，TGR=1.0，TWR=1.5，循環(huán)次數(shù)為1。

圖2 研究區(qū)子圖及4種厚云檢測方法的云掩膜結果

圖3是以ACCA算法為比較依據(jù)所展示的最大似然分類、ISODATA分類以及多尺度循環(huán)檢測算法所得到的不同尺寸水平下的云塊數(shù)量及其一致性對比圖，其中橫坐標代表云塊的尺寸水平，左側縱坐標代表所檢測出的云塊數(shù)量，右側縱坐標為該方法與作為對比依據(jù)的ACCA算法的檢測結果的總體一致性。圖中的研究區(qū)子圖1、子圖2、子圖3分別是與圖2對應的3個研究區(qū)子圖。

從圖2可知：多尺度循環(huán)檢測算法的云掩膜結果在整體上具有比ACCA算法和傳統(tǒng)的直接分類方法(最大似然與ISODATA分類)更為拓展的云邊界，但更容易產(chǎn)生一些碎小的誤檢點，這表明該算法對于厚云(尤其是碎云)的敏感性較高，并且可通過形態(tài)學操作得到有效改善。從圖3中3個研究區(qū)子圖的對比結果可以發(fā)現(xiàn)：多尺度循環(huán)檢測算法在不同尺寸水平下所檢測到的云塊數(shù)量更為接近ACCA算法，但在檢測結果的總體一致性上要低于傳統(tǒng)的直接分類方法(最大似然分類與ISODATA分類)。

表1統(tǒng)計了3個研究區(qū)子圖下最大似然分類、ISODATA分類、多尺度循環(huán)檢測算法與作為對比依據(jù)的ACCA算法的檢測結果的一致性詳細信息。從中可知：一方面，傳統(tǒng)的直接分類方法在誤檢率與總體一致性上要稍優(yōu)于本研究的算法，這與圖2、圖3對比分析后的結果一致；但另一方面，多尺度循環(huán)檢測算法在漏檢率與Kappa系數(shù)這兩項對比數(shù)據(jù)上占有較大優(yōu)勢。

圖3 最大似然分類、ISODATA分類、多尺度循環(huán)檢測算法與參考方法在不同云塊尺寸下的一致性

綜上，多尺度循環(huán)檢測算法與ACCA算法對于不同尺寸云塊的檢測結果差異(對應于誤檢率)主要來自于小塊云團或碎云(小于50個像元，見圖3)，兩者在有效厚云像元的檢測上具有很高的一致性，這主要體現(xiàn)在表1中多尺度循環(huán)檢測算法的低漏檢率(小于4%)與高Kappa系數(shù)(大于92%)，即該算法與ACCA算法的實際檢測結果更為接近；與ACCA算法相比，傳統(tǒng)的直接分類方法所檢測到的有效云像元數(shù)量較少，因而具有更高的總體一致性，但它們的高漏檢率說明傳統(tǒng)的直接分類方法對于有效厚云像元的檢測效果較差(表1)。

表1 不同厚云檢測算法與ACCA算法的一致性統(tǒng)計信息表

4 結論

1)對于Landsat ETM+圖像，所提出的多尺度循環(huán)檢測算法的精度水平與ACCA算法相近，但優(yōu)于傳統(tǒng)的直接分類方法(以最大似然和ISODATA分類為例)；此外，該算法還可通過設置循環(huán)參數(shù)不斷優(yōu)化檢測結果。

2)多尺度循環(huán)檢測算法在一定程度上受纓帽變換對光譜空間重構能力的影響，而轉(zhuǎn)換系數(shù)的提取又要依賴于研究區(qū)的地理生態(tài)和季相特征，因此，在實際應用中應針對研究區(qū)所處的區(qū)域生態(tài)和季相因子來統(tǒng)計出可靠的纓帽轉(zhuǎn)換系數(shù)。

3)多尺度循環(huán)檢測算法既不需要熱波段數(shù)據(jù)，也不需要其他的輔助參考數(shù)據(jù)。當研究區(qū)具有可靠的纓帽轉(zhuǎn)換系數(shù)時，該算法能夠?qū)崿F(xiàn)高度的自動化處理，因而可以對常用的高分辨率遙感數(shù)據(jù)進行快速而精確的厚云檢測處理(如SPOT[18]、QuickBird等衛(wèi)星)。

[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