基于改進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%。目前算法存在的在薄云和碎云條件下，云檢測精度較厚云條件下降以及像元庫數據由人工選擇組成，主觀性較強，可能導致數據庫的構建存在偏差的問題，是下一步需要研究的內容。