高山復(fù)雜地形區(qū)域MODIS積雪產(chǎn)品的云移除算法

馬強(qiáng) 謝紹璞 張爽

摘? 要：免費(fèi)可用的MODIS積雪覆蓋區(qū)域圖（SCA）因其較高的時(shí)空分辨率以及積雪/云的辨別能力一直被廣泛用于積雪覆蓋變化的研究。然而，云遮擋阻礙了MODIS積雪產(chǎn)品的廣泛應(yīng)用，尤其是在高山區(qū)域。在高山復(fù)雜地形區(qū)域，空間局部異質(zhì)性相對(duì)較高并且在原始的MODIS SCA圖中，云覆蓋通常是連續(xù)的和持續(xù)的，尤其是在冬季?；谝陨系膯栴}，文章通過使用Terra和Aqua日結(jié)合、時(shí)間移動(dòng)平均（TMA）推測(cè)、極值雪線法以及K最近鄰 （KNN）算法對(duì)高山復(fù)雜地形區(qū)域的云覆蓋進(jìn)行移除，從而產(chǎn)生無云的日MODIS積雪覆蓋圖。在相對(duì)驗(yàn)證中，文章所提出的云移除算法的總體精度達(dá)到83.12%。因此，文章中提出的除云算法能夠有效的去除高山復(fù)雜地形區(qū)域中原始MODIS積雪產(chǎn)品的云覆蓋。

關(guān)鍵詞：MODIS積雪數(shù)據(jù);云移除;高山區(qū)域;復(fù)雜地形;時(shí)間移動(dòng)平均推測(cè);瞬時(shí)區(qū)域

中圖分類號(hào)：P426.635? ? ? 文獻(xiàn)標(biāo)志碼：A? ? ? ? ?文章編號(hào)：2095-2945（2020）29-0001-07

Abstract： The freely available MODIS Snow Cover Area （SCA） maps have been widely used to investigate snow coverage changes because of their higher spatio-temporal resolution and snow/cloud discrimination capabilities. However， cloud obscuration hinders their extensive application， especially in mountain regions. In alpine areas with complex to pography， spatial local heterogeneity is relatively higher and cloud coverage in original MODIS SCA maps is frequently continuous and persistent， especially in winter. Aimed at problems mentioned above， Terra and Aqua combination， Temporal Moving-Average（TMA） deduction， and extreme snow-line method and K-Nearest Neighbor （KNN） algorithm were used to remove clouds from MODIS SCA data in Alpine areas with complex to pography and yield cloud-free MODIS SCA data. In relative validation， overall accuracy of the cloud removal methodology proposed in this paper was 83.12%. Therefore， the solution proposed can effectively eliminate clouds from original MODIS snow cover data in alpine areas with complex to pography.

Keywords： MODIS snow data; cloud removal; alpine area; complex to pography; Temporal Moving-Average （TMA） deduction; transient zones

1 概述

由于實(shí)地測(cè)量的局限性，遙感被廣泛用于偏遠(yuǎn)地區(qū)或者高山區(qū)域從而探測(cè)積雪變化;同時(shí)遙感提供了多尺度的積雪觀測(cè)，這對(duì)于監(jiān)測(cè)積雪場(chǎng)的演化是非常有用的信息。積雪覆蓋是積雪場(chǎng)的一個(gè)重要因子。MODIS二值積雪覆蓋圖因其高的時(shí)間分辨率（一天）、中尺度的空間分辨率（500m）以及積雪/云的辨別能力，從而在眾多衛(wèi)星衍生的積雪覆蓋產(chǎn)品中具有明顯的優(yōu)勢(shì)[2-3]。然而，云遮擋阻礙了其廣泛的應(yīng)用，尤其是在高山區(qū)域。

在高山復(fù)雜地形區(qū)域，積雪累積和消融不僅受到海拔的影響，而且因?yàn)椴煌奶栞椛?、能量平衡以及迎風(fēng)/背風(fēng)等因素使得坡度和坡向也會(huì)對(duì)積雪積累和消融產(chǎn)生影響。這些因素會(huì)導(dǎo)致積雪在高山區(qū)域是非連續(xù)性分布的并且造成較大的空間局部異質(zhì)性。例如，在冬季坡度較陡的高海拔區(qū)域其地表覆蓋是陸地，但是毗鄰于該陸地區(qū)域的低洼平坦地區(qū)其地表覆蓋卻是積雪。即使將坡度、坡向和海拔這些地形因素融入到基于臨近像素的空間濾波中，該方法的除云精度和效果依然比較差。同時(shí)，高山區(qū)域的云覆蓋通常是持續(xù)的，從而使得時(shí)間濾波的精度或者除云效果將會(huì)降低。此外，經(jīng)過時(shí)空結(jié)合方法處理后MODIS 積雪覆蓋圖剩余的云主要集中在瞬時(shí)區(qū)域。然而，瞬時(shí)區(qū)域的積雪場(chǎng)是極其不穩(wěn)定的，尤其是在冬季以外的其他季節(jié)。因此，基于長時(shí)間序列分析的積雪循環(huán)方法不能夠有效的解決這個(gè)問題。

針對(duì)上述問題，本文的目的是提出一個(gè)有效的方法對(duì)高山復(fù)雜地形區(qū)域中MODIS日積雪產(chǎn)品的云覆蓋進(jìn)行去除并且產(chǎn)生無云的MODIS日積雪覆蓋數(shù)據(jù)。

2研究區(qū)和數(shù)據(jù)源

2.1 研究區(qū)

臺(tái)蘭河流域位于中國新疆烏魯木齊自治區(qū)溫宿縣境內(nèi)（圖1）。以臺(tái)蘭河水文站控制的流域面積總共是1324 km2，其空間范圍是從東經(jīng)80°04′E -80°35′E和北緯41°35′N-42°05′N。在30m DEM數(shù)據(jù)中所提取的流域面積是1361km2，因此，本研究中所認(rèn)定的流域面積也是1361km2。流域海拔的變化范圍是從1577到7421m asl，其平均海拔是3739m asl。在地形上，大約96%的區(qū)域高于2000m asl，73%的區(qū)域高于2000m asl以及46%的區(qū)域高于4000m asl。在山麓平原區(qū)，植被類型主要以荒漠草場(chǎng)和灌叢為主;在中海拔區(qū)域分布著大片的云杉和亞高山草甸;在高海拔區(qū)域普遍發(fā)育著永久性積雪覆蓋和冰川冰。

流域受到地形抬升、西風(fēng)環(huán)流和北冰洋水汽的影響，造成大氣水汽傳送從西到東、從北到南逐漸減少。因此該研究區(qū)屬于典型的干旱性大陸氣候，并且春夏兩季高山區(qū)域的冰雪融水是臺(tái)蘭河的主要補(bǔ)給源。冰雪融水約占到年徑流量的60%以上，雨水的貢獻(xiàn)所占比例相對(duì)較小并且主要降雨發(fā)生在暴雨時(shí)期。

2.2 MODIS數(shù)據(jù)

本研究的除云算法是對(duì)MODIS日積雪覆蓋數(shù)據(jù)MOD10A1（V005）和MYD10A1（V005）進(jìn)行處理。通過訪問美國國家冰雪信息中心（National Snow and Ice Data Center， NSIDC）網(wǎng)址，可以獲取免費(fèi)的MODIS日積雪覆蓋數(shù)據(jù)，該數(shù)據(jù)采用的正弦投影以HDF的格式發(fā)布。整個(gè)研究區(qū)需要兩個(gè)MODIS塊（H23V04和H24V04）進(jìn)行拼接才能對(duì)其完全覆蓋。我們使用的MODIS積雪覆蓋數(shù)據(jù)其日期范圍是從2007年1月1日開始到2011年12月31日截止。1820幅MOD10A1影像和1825幅MYD10A1影像將會(huì)用于該研究（表1）。缺失的影像將會(huì)用完全被云覆蓋的影像補(bǔ)充。MODIS重投影工具（MODIS Reprojection Tool， MRT） 將用于數(shù)字拼接并且將HDF格式轉(zhuǎn)化成GeoTIF格式;轉(zhuǎn)化后生成的影像將會(huì)被重投影到WGS1984 UTM 43N帶，并且采用最近鄰插值將重投影的影像重采樣成500m。最后，重采樣影像的11種像素類別會(huì)被重新分為3類，“湖冰”類別將會(huì)被納入到“積雪”類別;“湖泊”和“海洋”類別將會(huì)劃歸到“非積雪”（陸地）類別;“云”和其他種類將會(huì)劃分到“云”類別。

2.3 Landsat數(shù)據(jù)

Landsat TM/ETM+有六個(gè)空間分辨率為30m的波段影像，通過訪問USGS數(shù)據(jù)服務(wù)器（http：//earthexplorer.usgs.gov/）可以進(jìn)行免費(fèi)下載。我們選擇了13幅Landsat TM/ETM+影像與同日Terra衛(wèi)星所提供的MOD10A1積雪產(chǎn)品進(jìn)行多傳感器結(jié)合，從而生成非填充驗(yàn)證數(shù)據(jù)對(duì)本文中所描述的除云算法進(jìn)行精度評(píng)估（表2）。Landsat TM/ETM+影像的選取標(biāo)準(zhǔn)如下：（1）Landsat TM/ETM+影像是無云的或者云覆蓋極少的;（2）同日的原始MODIS積雪產(chǎn)品（MOD10A1和MYD10A1）的云覆蓋高于30%。

因?yàn)長andsat TM和ETM+傳感器具有相似的光譜波段，積雪制圖算法生成二值積雪覆蓋圖的時(shí)候會(huì)把它們當(dāng)成相同的傳感器。因此，在本文中，Landsat TM和ETM+將會(huì)被簡(jiǎn)稱為“Landsats”。Landsats影像預(yù)處理的過程如下：（1）首先將Landsats影像重投影到WGS1984 UTM 43N帶;（2）為了能夠適應(yīng)積雪制圖算法，將Landsats波段影像初始的DN值轉(zhuǎn)化為大氣頂部（Top of the Atmosphere ，TOA）反射率。積雪在可見光波段（Landsats band2）有高的反射率并且在近紅外波段（Landsats band5）是接近于0的反射率，積雪制圖算法利用這些光譜反射特性探測(cè)積雪。歸一化積雪指數(shù)（Normalized Difference Snow Index， NDSI）的計(jì)算公式如下：

一個(gè)二值的積雪覆蓋圖通過設(shè)定NDSI≥0.40、近紅外波段 （Landsat Band 4/MODIS Band 2）的反射率大于0.11和可見光波段（Landsats Band 3） 的反射率大于0.10而得到。然而，積雪制圖算法通常會(huì)將陰影區(qū)域的積雪錯(cuò)分為陸地即產(chǎn)生低估誤差。Macander對(duì)誤差的原因進(jìn)行了詳細(xì)的分析。為了避免這種低估誤差，我們對(duì)Landsats影像中的陰影進(jìn)行消除，然后再使用積雪制圖算法獲取Landsats二值積雪覆蓋圖。最后，將獲取的積雪覆蓋圖通過判定一個(gè)500m柵格內(nèi)其對(duì)應(yīng)的30m積雪像素所占比例，將其重采樣成500m。判定條件如下：如果一個(gè)500m柵格內(nèi)其對(duì)應(yīng)的30m積雪像素所占比例超過50%，則將該柵格像素劃分為積雪;如果所占比例小于50%，則將該像素劃分為陸地;否則，該500m像素將會(huì)被指定為云像素。

2.4 數(shù)字高程模型 （Digital Elevation Model， DEM）

ASTER GDEM 版本2數(shù)據(jù)（GDEM 2）由ASTER GDEM團(tuán)隊(duì)研發(fā)并且自2011年11月17日公布使用。數(shù)據(jù)來源于中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(tái)（http：//www. gscloud.cn）。整個(gè)研究區(qū)域需要兩景GDEM（N41E080和N42E080）才能夠被完全覆蓋。首先對(duì)兩景GDEM進(jìn)行拼接，然后將其重投影到WGS1984 UTM 43N帶，最后重投影的數(shù)據(jù)將會(huì)用于臺(tái)蘭河流域的提取。隨后使用雙線性插值方法將流域DEM圖由30m重采樣成500m，同時(shí)利用重采樣的DEM圖獲取一個(gè)坡向圖。

3 方法

3.1 除云算法

如概述所說，本文的目的是提出一個(gè)有效的除云算法用于去除高山復(fù)雜地形區(qū)域中MODIS積雪產(chǎn)品的云覆蓋。該方法由四個(gè)步驟組成：Terra和Aqua的日結(jié)合、時(shí)間移動(dòng)平均（Temporal Moving-Average， TMA）推測(cè)、極值雪線法和K最近臨（K-Nearest Neighbor， KNN）算法。方法按順序依次執(zhí)行每個(gè)步驟，每個(gè)步驟的輸出是后續(xù)步驟的輸入。整個(gè)方法通過Python語言實(shí)現(xiàn)。

3.1.1 日結(jié)合

步驟1是將MOD10A1（V005版本）和MYD10A1（V005版本）提供的MODIS / Terra和MODIS / Aqua日積雪覆蓋產(chǎn)品進(jìn)行日結(jié)合。在本研究中，我們使用MOD10A1優(yōu)先的合成規(guī)則，即MOD10A1中的云覆蓋像素將會(huì)被MYD10A1中的無云像素替代。

3.1.2 時(shí)間移動(dòng)平均推測(cè)

步驟2是基于積雪覆蓋的周期性變化對(duì)云遮擋下真實(shí)的地表覆蓋進(jìn)行恢復(fù)。一個(gè)區(qū)域的積雪覆蓋通常從冬季至夏季逐漸減少，從夏季至冬季逐漸增加，呈現(xiàn)出區(qū)間變化。在一個(gè)水文年內(nèi)，積雪覆蓋面積呈現(xiàn)出四個(gè)區(qū)間變化：?jiǎn)握{(diào)遞增區(qū)間、單調(diào)遞減區(qū)間、凸區(qū)間和凹區(qū)間。在單調(diào)遞減區(qū)間，當(dāng)日降雪的概率大于右鄰域降雪的概率;在單調(diào)遞增區(qū)間，當(dāng)日降雪的概率大于左鄰域降雪的概率;在凸區(qū)間，當(dāng)天降雪的概率大于左鄰域和右鄰域降雪的概率;然而，在凹區(qū)間，當(dāng)日降雪的概率小于左鄰域和右鄰域降雪的概率。通過以上總結(jié)，我們發(fā)現(xiàn)當(dāng)日降雪的概率大于增區(qū)間、減區(qū)間和凸區(qū)間左鄰域或右鄰域的降雪概率，從而可以推測(cè)出當(dāng)日云遮擋下真實(shí)的地表覆蓋。我們使用上述的積雪區(qū)間變化特征，可以檢測(cè)出位于凹區(qū)間的云覆蓋像素并且對(duì)這些云像素不進(jìn)行處理。

如圖2所示，深青色的線為積雪消耗曲線（A），可以發(fā)現(xiàn)日積雪覆蓋變化較大; 紅線為積雪消耗月移動(dòng)均值線（B），月移動(dòng)平均積雪覆蓋變化趨勢(shì)比較穩(wěn)定，從1月至8月其月移動(dòng)平均積雪覆蓋逐漸減少，8月至12月其月移動(dòng)平均積雪覆蓋逐漸增加，在8月初至8月中旬其平均積雪覆蓋達(dá)到最低并且較為穩(wěn)定;藍(lán)線是決策線（C）并且通常低于紅線，但是藍(lán)線在凹區(qū)間要高于紅線。決策線是前一個(gè)月和后一個(gè)月（31天）的月平均積雪覆蓋的最小值確定得到的。

基于以上的理論研究，我們對(duì)云像素真實(shí)的地表覆蓋進(jìn)行推測(cè)。首先，對(duì)當(dāng)日云像素前31天和后31天的降雪概率分別進(jìn)行計(jì)算。降雪概率計(jì)算的公式如下：

NS是31天內(nèi)積雪出現(xiàn)的個(gè)數(shù); NL是31天內(nèi)陸地出現(xiàn)的個(gè)數(shù)。當(dāng)日云像素前31天和后31天積雪出現(xiàn)的概率分別標(biāo)記為P-31→0和P0→+31。其次，選擇P-31→0和 P0→+31之間的最小值，并且標(biāo)記為Pmin。第三，計(jì)算當(dāng)日云像素前15日至后15日積雪出現(xiàn)的概率，作為當(dāng)日積雪出現(xiàn)概率并且標(biāo)記為P-15→+15。最后，如果P-15→+15Pmin>0.60，則將該云像素指定為積雪像素。Pmin>0.60主要是為了避開瞬時(shí)區(qū)域，因?yàn)樵搮^(qū)域的積雪場(chǎng)不穩(wěn)定并且存在多個(gè)積雪/陸地循環(huán)。相同的原理將用于把云像素重新分類為陸地像素。

3.1.3 極值雪線法

在本研究中，我們采用了Gafurov和Bárdossy提出的積雪過渡高程法并且對(duì)其進(jìn)行了改進(jìn)[17，21]。陸地線（Land Line）仍是積雪像素中最低的海拔作為閾值，然而雪線（Snow Line）是陸地像素中最高的海拔作為閾值。通常，雪線高于陸地線;如果陸地線高于雪線，則表明雪線和陸地線之間的像素全部被云遮擋。因此，我們利用如下規(guī)則對(duì)位于極高海拔和極低海拔位置的云像素進(jìn)行分類：

（1）如果陸地線低于雪線，則雪線以上的云像素變成積雪，陸地線以下的云像素變?yōu)殛懙?

（2）如果陸地線高于雪線，則陸地線變?yōu)檠┚€，而雪線變?yōu)殛懙鼐€，將雪線以上的云像素指定為積雪像素，陸地線以下的云像素指定為陸地像素;

（3）如果只有雪線而沒有陸地線，則將雪線變?yōu)殛懙鼐€，并將陸地線以下的云像素重新分類積雪像素;

（4）如果只有陸地線而沒有雪線，則陸地線變?yōu)檠┚€，并且雪線以上的云像素重新劃分為積雪像素。

3.1.4 最近鄰（K-Nearest Neighbor， KNN）算法

在經(jīng)過極值雪線 （步驟3） 處理后， 云覆蓋主要集中在瞬時(shí)區(qū)域，在這個(gè)區(qū)域的積雪場(chǎng)是極其不穩(wěn)定的。因此，使用極少的云像素或者一個(gè)積雪循環(huán)去估算瞬時(shí)區(qū)域的云像素是不適合的。為了解決這個(gè)問題，我們引入了機(jī)器學(xué)習(xí)中的K最近鄰（KNN）算法。KNN算法是由Cover和Hart提出的，該算法的基本原理如下：在樣本空間中，如果一個(gè)未知像素周圍最鄰近的K個(gè)已知像素主體是屬于一個(gè)特定的類別，則這個(gè)未知樣本將會(huì)被指定為這個(gè)特定類別。一個(gè)已知樣本通常是已知樣本和未知樣本之間的歐式距離確定的。在三維空間中兩個(gè)點(diǎn)的歐式距離計(jì)算公式如下：

其中y是像素行號(hào)索引;x是像素列號(hào)索引;t是像素日期索引。KNN算法的樣本空間大小取決于樣本的個(gè)數(shù)。為了確保參與估算云像素的無云像素個(gè)數(shù)充足，我們?cè)O(shè)定K=27。在高山區(qū)域，我們假定時(shí)間自相關(guān)大于空間自相關(guān)。這是因?yàn)樵诟呱絽^(qū)域由于較大的地形起伏造成空間局部異質(zhì)性相對(duì)較高，然而，積雪覆蓋變化相對(duì)較低因?yàn)橄鄬?duì)積雪覆蓋位于相對(duì)較高的海拔。因此，必須確保參與估算云覆蓋像素的無云像素具有較小空間距離。我們?cè)O(shè)置時(shí)間權(quán)重為2，然后將原始的時(shí)間軸縮小到一半，空間軸保持不變。在確定27個(gè)無云像素之后，將會(huì)對(duì)云像素的真實(shí)地表覆蓋進(jìn)行推測(cè)。如果積雪像素的個(gè)數(shù)大于13，則云覆蓋像素將會(huì)被重新劃分為積雪像素;否則，該云像素將會(huì)被指定為陸地像素。

3.2 驗(yàn)證和評(píng)估

在本研究中，我們引入一種新的驗(yàn)證方法，將其稱作非填充驗(yàn)證。該方法使用Landsats影像對(duì)除云算法的精度進(jìn)行評(píng)估。使用Landsats影像用作驗(yàn)證是基于以下的考慮：Landsat衛(wèi)星的過境時(shí)間要早于MODIS衛(wèi)星，并且在白天隨著時(shí)間的推移云覆蓋逐漸增加。這很有可能在同一天MODIS影像的含云量較高而Landsat衛(wèi)星獲取的波段影像是無云的或者是含云量極少的，這將為非填充驗(yàn)證提供了可能。非填充驗(yàn)證的“地表真實(shí)”是通過多傳感器結(jié)合生成的。MOD10A1中的云像素將會(huì)被Landsats二值積雪覆蓋圖中相同位置的非云像素進(jìn)行替代，從而生成非填充驗(yàn)證的“地表真實(shí)”。在非填充驗(yàn)證中，我們采用一致度 （Degree Of Agreement， DA）這個(gè)指標(biāo)對(duì)本研究中所提出的算法精度進(jìn)行評(píng)估。一致度（Degree Of Agreement， DA）是指經(jīng)過除云算法處理正確分類的像素占總的分類像素的比例，該指標(biāo)均采用百分比形式，具體計(jì)算公式如下：

式中，a和d分別是指正確分類為積雪和陸地的像素個(gè)數(shù); c是指將積雪錯(cuò)分為陸地的像素個(gè)數(shù); b是將陸地錯(cuò)分為積雪的像素個(gè)數(shù)?？傮w一致性（Overall Agreement， OA）是通過各個(gè)測(cè)試日不同表現(xiàn)性能的平均值得出的，對(duì)每個(gè)測(cè)試日對(duì)所檢查像素正確率的百分比進(jìn)行加權(quán)得到的，其計(jì)算公式為：

其中i代表驗(yàn)證日;RNi代表在原始的Terra影像（MOD10A1）中是云像素而在對(duì)應(yīng)的“地表真實(shí)”中是非云像素的個(gè)數(shù)。

4 結(jié)果

本文中描述的四個(gè)除云步驟能夠完全消除原始的MODIS積雪覆蓋數(shù)據(jù)中的云覆蓋。圖3a顯示了MOD10A1、MYD10A1和使用步驟1-4之后月平均云覆蓋。MODIS日積雪覆蓋產(chǎn)品（MOD10A1和MYD10A1）以及經(jīng)過各個(gè)步驟后的月平均覆蓋如圖3b所示。如圖3a所示，原始的MODIS積雪產(chǎn)品 （MDO10A1和MYD10A1） 中的月平均云覆蓋與季節(jié)有關(guān)，在冬季的云覆蓋達(dá)到最高并且在其他季節(jié)逐漸減少。MYD10A1的云覆蓋高于MOD10A1中的云覆蓋是因?yàn)锳qua衛(wèi)星的過境時(shí)間晚于Terra衛(wèi)星并且云覆蓋隨著白天時(shí)間的推移在逐漸增加。與原始的MOD10A1產(chǎn)品相比，經(jīng)過日結(jié)合處理（步驟1）后，月平均云覆蓋減少了7.06%-11.29%并且月積雪覆蓋略有增加。當(dāng)MOD10A1和MYD10A1的云結(jié)構(gòu)不同時(shí)，日結(jié)合的除云效果較好。時(shí)間移動(dòng)平均推測(cè)（步驟2）是最有效的步驟并且云覆蓋又減少了26.58%-46.06%;同時(shí)，積雪覆蓋有了明顯增加，尤其是在雪季 （10月至3月）。極值雪線法 （步驟3） 能夠有效的消除位于極低海拔和極高海拔持續(xù)性的云覆蓋，并且不需要對(duì)經(jīng)過步驟2處理后的影像設(shè)定特定的無云比例進(jìn)行除云。在經(jīng)過第3步處理后，云覆蓋率降低了2.02%-7.41%。與除云效果相比，極值雪線法對(duì)積雪覆蓋的恢復(fù)在大多數(shù)月份都是較少的。這主要是因?yàn)闃O值雪線法去除的云的位置相對(duì)固定，并且低海拔區(qū)域出現(xiàn)積雪的概率比較低。K最近鄰算法（步驟4）能夠完全去除影像中剩余的云覆蓋，并且增加了更多的積雪。這表明該方法能夠有效的解決經(jīng)過先前三個(gè)步驟處理后位于瞬時(shí)區(qū)域剩余的、持續(xù)時(shí)間較長的云覆蓋。圖4是2009年3月6日原始的MODIS積雪產(chǎn)品以及經(jīng)過各個(gè)除云步驟處理后剩余的云覆蓋圖。

5 驗(yàn)證

在非填充驗(yàn)證中，驗(yàn)證樣本由13幅多傳感器結(jié)合影像組成用于評(píng)估本文中所提出的除云算法各個(gè)步驟的精度以及除云算法的總體精度，這些影像幾乎全部來自冬季（表3）?？梢钥闯鋈战Y(jié)合（step 1）的總體一致性低于時(shí)間移動(dòng)平均推測(cè)（step 2）和極值雪線法（step 3），其原因在3.1.1節(jié)進(jìn)行了詳細(xì)的分析。雖然移動(dòng)時(shí)間平均推測(cè)的輸入是精度相對(duì)較低的日結(jié)合影像，但是當(dāng)日云像素的推測(cè)是基于積雪＼陸地出現(xiàn)的概率從而確保經(jīng)過移動(dòng)時(shí)間平均推測(cè)后日結(jié)合方法的不確定性誤差能夠被最大程度的消除并且獲取可靠的結(jié)果。移動(dòng)時(shí)間平均推測(cè)的云移除效果最好并且其精度最高。同時(shí)，該步驟的方差僅僅是0.18%，表明該方法是極其穩(wěn)定的。極值雪線法的總體精度要比日結(jié)合的精度高，這是因?yàn)樵摲椒ㄏ脑聘采w位置固定的并且類型是確定的，而日結(jié)合去除的云覆蓋較為離散且有一定的不確定性。然而，極值雪線法的方差是四個(gè)步驟中最高的，這可能是因?yàn)樘綔y(cè)器飽和、云或云覆蓋的影響造成生成的Landsat 積雪覆蓋圖在較高的海拔有較大的低估誤差，尤其是在冬季。在經(jīng)過上述三步處理后，剩余的云覆蓋主要位于瞬時(shí)區(qū)域，這個(gè)區(qū)域的積雪場(chǎng)是極其不穩(wěn)定。因此，該方法的精度較低是可預(yù)見的并且能夠接受的。K最近鄰算法（step 4）的精度最低，但是其方差要低于日結(jié)合。這表明更多的無云像素參與瞬時(shí)區(qū)域云像素的估算，能夠確保獲取穩(wěn)定的結(jié)果。通過對(duì)本文所提出的除云算法進(jìn)行整體評(píng)估，其總體精度達(dá)到83.12%，同時(shí)該算法的方差僅有0.21%并且原始的MODIS SCA圖像中的云覆蓋都被消除。因此，本文所提出的算法是適用于去除高山復(fù)雜地形區(qū)域MODIS積雪產(chǎn)品中的云覆蓋。

TFC 和AFC分別代表Terra和Aqua影像中的云比例;RF代表經(jīng)過各個(gè)步驟處理后Terra影像中的剩余的云量;step 1-4是指各云移除步驟， CRM是指由各個(gè)步驟組成的整個(gè)的云移除算法（Cloud Removal Methodology）。

6 討論與結(jié)論

本文的研究目的是為了提出一個(gè)適用于高山復(fù)雜地形區(qū)域的除云算法。為了克服高山復(fù)雜區(qū)域連續(xù)性和持續(xù)性云覆蓋以及地物分布空間局部異質(zhì)性較大的問題，我們提出一個(gè)由四個(gè)步驟組成的除云算法并且原始的MODIS積雪產(chǎn)品（MOD10A1和MYD10A1） 的云覆蓋在經(jīng)過我們的方法處理后能夠被完全的消除。第一個(gè)步驟是將Terra和Aqua衛(wèi)星提供的MODIS日積雪覆蓋圖（MOD10A1和MYD10A1）進(jìn)行日結(jié)合，該方法是利用了MODIS衛(wèi)星一日兩次獲取數(shù)據(jù)的獨(dú)特優(yōu)勢(shì)。第二個(gè)步驟是時(shí)間移動(dòng)平均推測(cè)，該方法使用一個(gè)水文年內(nèi)積雪覆蓋變化去推測(cè)非瞬時(shí)區(qū)域的云覆蓋并且在四個(gè)云移除步驟中精度最高、效果最好。同時(shí)，該方法具有最小的方差、線性時(shí)間復(fù)雜度并且易于實(shí)施。第三個(gè)步驟是極值雪線法，該方法使用陸地像素的最高海拔（雪線）和積雪像素的最低海拔（陸地線，Land Line）對(duì)位于極高海拔和極地海拔的云覆蓋進(jìn)行去除。這個(gè)方法能夠有效的避免了Parajka 提出的區(qū)域雪線法將低洼平原區(qū)的積雪覆蓋錯(cuò)分為陸地的低估誤差，并且其除云效果相對(duì)較好。位于瞬時(shí)區(qū)域云覆蓋的去除問題一直是困擾除云算法的一個(gè)難題，為了解決這個(gè)難題我們?cè)谧詈笠徊揭隟最近鄰（K-Nearest Neighbor， KNN）算法。KNN算法是一種相對(duì)簡(jiǎn)單的機(jī)器學(xué)習(xí)方法，它使用K個(gè)已知的最近鄰樣本來推斷未知樣本的特定類別。雖然KNN算法的精度相對(duì)較低，但是剩余的云覆蓋能夠被完全去除。

通過與其他除云算法相比，我們提出的解決方案具有最高的除云精度和最佳的除云效果。然而，使用本研究所提出的KNN算法對(duì)瞬時(shí)區(qū)域的云覆蓋進(jìn)行處理后所得結(jié)果的精度仍然較低，但是我們認(rèn)為可能會(huì)有更好的機(jī)器學(xué)習(xí)算法從而達(dá)到更高的精度。

參考文獻(xiàn)：

[1]ARSENAULT K R， HOUSER P R， DE LANNOY G J M. Evaluation of the MODIS snow cover fraction product[J]. Hydrological Processes， 2014， 28（3）： 980-998.

[2]HALL D K， RIGGS G A， SALOMONSON V V， et al. MODIS snow-cover products[J]. Remote Sensing of Environment， 2002， 83： 181-194.

[3]PARAJKA J， BLSCHL G. The value of MODIS snow cover data in validating and calibrating conceptual hydrologic models [J]. Journal of Hydrology， 2008， 358（3-4）： 240-258.

[4]AULT T W， CZAJKOWSKI K P， BENKO T， et al. Validation of the MODIS snow product and cloud mask using student and NWS cooperative station observations in the Lower Great Lakes Region[J]. Remote Sensing of Environment， 2006， 105（4）： 341-353.

[5]PARAJKA J， BLSCHL G. Validation of MODIS snow cover images over Austria[J]. Hydrology and Earth System Sciences Discussions， 2006， 10（5）： 1569-1601.

[6]KRAJí P， HOLKO L， PERDIGO R A P， et al. Estimation of regional snowline elevation （RSLE） from MODIS images for seasonally snow covered mountain basins[J]. Journal of Hydrology， 2014， 519：1769-1778.

[7]ZHOU X， XIE H， HENDRICKX J M H. Statistical evaluation of remotely sensed snow-cover products with constraints from streamflow and SNOTEL measurements[J]. Remote Sensing of Environment， 2005， 94（2）： 214-231.

[8]TONG J， DéRY S J， JACKSON P L. Topographic control of snow distribution in an alpine watershed of western Canada inferred from spatially-filtered MODIS snow products [J]. Hydrology and Earth System Sciences Discussions， 2008， 5（4）： 2347-2371.

[9]RIGGS G A， HALL D K， SALOMONSON V V. MODIS snow products user guide to collection 5[J/OL]. Boulder， Colorado USA： NASA National Snow and Ice Data Center Distributed Active Archive Center， 2006， https：//modis-snow-ice.gsfc.nasa.gov/？c=userguides.

[10]MACANDER M J， SWINGLEY C S， JOLY K， et al. Landsat-based snow persistence map for northwest Alaska [J]. Remote Sensing of Environment， 2015， 163： 23-31.

[11]Salomonson V V， Appel I. Development of the Aqua MODIS NDSI fractional snow cover algorithm and validation results [J]. IEEE Transactions on Geoscience and Remote Sensing， 2006，44 （7）： 1747-1756.

[12]Hall D K， Riggs G A. Accuracy assessment of the MODIS snow products[J]. Hydrological Processes，2007，21（12）：1534-1547.

[13]GAFUROV A， BáRDOSSY A. Cloud removal methodology from MODIS snow cover product [J]. Hydrology and Earth System Sciences， 2009， 13（7）： 1361-1373.

[14]COVER T M， HART P E. Nearest neighbor pattern classification [J]. IEEE Transactions on Information Theory， 1967， 13（1）： 21-27.

[15]KRAJí P， HOLKO L， PERDIGO R A P， et al. Estimation of regional snowline elevation （RSLE） from MODIS images for seasonally snow covered mountain basins[J]. Journal of Hydrology，2014，519：1769-1778.

[16]沈永平，劉時(shí)銀，丁永建，等.天山南坡臺(tái)蘭河流域冰川物質(zhì)平衡變化及其對(duì)徑流的影響[J].冰川凍土，2003，25（2）：124-129.

[17]蔣艷，周成虎，程維明.阿克蘇河流域徑流補(bǔ)給及徑流變化特征分析[J].自然資源學(xué)報(bào)，2005，20（1）：29-36.

[18]田龍.臺(tái)蘭河流域積雪覆蓋面積變化特征分析[J].西北水電，2015（02）：13-16.

[19]謝昌衛(wèi)，丁永建，劉時(shí)銀，等.托木爾峰南坡冰川水文特征及其對(duì)徑流的影響分析[J].干旱區(qū)地理，2004，27（4）：570-575.

[20]沈志寶，陳有虞.青藏高原平均云量的估算及其分布特征[J].高原氣象，1987，6（4）：49-58.

[21]孫美萍，姚曉軍，李忠勤，等.21世紀(jì)天山南坡臺(tái)蘭河流域徑流變化情景預(yù)估[J].氣候變化研究進(jìn)展，2012，8（5）：342-349.