侯小剛， 鄭照軍， 李 帥， 陳雪華， 崔 宇

(1.北京郵電大學(xué)網(wǎng)絡(luò)技術(shù)研究院，北京 100876； 2.中國(guó)氣象局國(guó)家衛(wèi)星氣象中心，北京 100081； 3.中國(guó)氣象局烏魯木齊沙漠氣象研究所，烏魯木齊 830002；4.新疆測(cè)繪科學(xué)研究院,烏魯木齊 830002)

0 引言

新疆是我國(guó)三大穩(wěn)定積雪區(qū)之一。積雪在新疆的氣候變化、水資源供給和農(nóng)牧業(yè)發(fā)展方面具有不可估量的作用； 但新疆也是雪災(zāi)多發(fā)區(qū)，由積雪引發(fā)的雪災(zāi)、暴風(fēng)雪、雪崩和融雪性洪水等自然災(zāi)害頻發(fā)，對(duì)當(dāng)?shù)卦斐闪藝?yán)重的危害[1-2]。新疆地域遼闊，但積雪監(jiān)測(cè)地面站點(diǎn)不足，大范圍積雪都分布在高山等地面臺(tái)站無法監(jiān)測(cè)的區(qū)域，因此，利用遙感技術(shù)對(duì)新疆積雪進(jìn)行實(shí)時(shí)動(dòng)態(tài)監(jiān)測(cè)和大范圍積雪資源調(diào)查是未來的發(fā)展方向[3-4]。利用遙感技術(shù)進(jìn)行積雪觀察迄今已經(jīng)超過半個(gè)世紀(jì)，形成了一系列積雪制圖算法及監(jiān)測(cè)產(chǎn)品[5]。在國(guó)外有SMMR，SSM/I和AMSR-E等微波產(chǎn)品，AVHRR，TM，ETM+，SPOT和MODIS等光學(xué)產(chǎn)品以及交互式多傳感器雪冰制圖系統(tǒng)(interactive multi-sensor snow ice mapping system，IMS)等多源數(shù)據(jù)合成產(chǎn)品； 在國(guó)內(nèi)則有FY-2C/2D/2E/2F/2G積雪范圍產(chǎn)品，F(xiàn)Y-3/MERSI和FY-3/VIRR等光學(xué)產(chǎn)品以及FY-3/MWRI等微波產(chǎn)品。MODIS因具有較高的時(shí)間和空間分辨率而在積雪監(jiān)測(cè)中具有絕對(duì)優(yōu)勢(shì)[6-9]。

目前，MODIS數(shù)據(jù)是積雪研究與實(shí)時(shí)監(jiān)測(cè)中獲取最便捷、應(yīng)用最廣泛、研究最深入的數(shù)據(jù)，已開發(fā)出了一系列基于MODIS的積雪產(chǎn)品，包括MOD10A2和MYD10A2的8 d積雪最大值合成產(chǎn)品(maximum snow extent，MSE)、MOD10A1和MYD10A1日積雪面積比例產(chǎn)品和日積雪面積二值產(chǎn)品[10-13]。但由于光學(xué)遙感產(chǎn)品會(huì)受到云的影響，很大程度上限制了其在積雪監(jiān)測(cè)中的應(yīng)用。據(jù)統(tǒng)計(jì)，在新疆地區(qū)積雪期(11月初—翌年3月底)云覆蓋占整個(gè)積雪期的80%，因而大大限制了MODIS數(shù)據(jù)10A1產(chǎn)品在新疆地區(qū)的應(yīng)用，即使部分MSE產(chǎn)品也有云像元的存在，因而如何獲取云覆蓋下的積雪信息(去云處理)，會(huì)直接影響MODIS產(chǎn)品數(shù)據(jù)在新疆地區(qū)積雪監(jiān)測(cè)中的精準(zhǔn)度。為了降低MODIS積雪產(chǎn)品中的云覆蓋，已有很多學(xué)者進(jìn)行了有效的去云算法研究。目前，去云算法大致可分為5類： ①基于積雪時(shí)間連續(xù)性法[6，9，14-20]，具體指上下午星合成、臨近日分析、靈活多日結(jié)合和基于時(shí)間序列分析等； ②基于積雪空間連續(xù)性法[6，9，14-20]，具體有臨近日4像元法和臨近日8像元法等； ③基于數(shù)字高程模型法[6，9，14-23]，具體有Snowl法、高程濾波法和高程掩模法等； ④多傳感器融合法[6，14-24]，大多采用的是與被動(dòng)微波融合去云法； ⑤其他方法[24-25]，包括計(jì)算積雪可能性和函數(shù)擬合積雪邊界等方法。通過對(duì)相關(guān)MODIS數(shù)據(jù)去云文獻(xiàn)[6，9，14-25]的總結(jié)，發(fā)現(xiàn)多數(shù)文獻(xiàn)對(duì)去云方法的結(jié)果驗(yàn)證主要集中在青藏高原等地區(qū)，而對(duì)新疆區(qū)域的積雪去云方法結(jié)果數(shù)據(jù)的精度驗(yàn)證較少。在多源遙感數(shù)據(jù)融合去云方面，主要采用微波數(shù)據(jù)具有的穿透性進(jìn)行云下信息的判識(shí)[26]，但微波數(shù)據(jù)低空間分辨率(25 km)的缺點(diǎn)是造成這類方法精度偏低的最主要原因，而很少有方法采用具有較高空間分辨率(4 km)且不受云影響的IMS雪冰產(chǎn)品數(shù)據(jù)進(jìn)行新疆區(qū)域內(nèi)逐日去云積雪產(chǎn)品數(shù)據(jù)的生成。針對(duì)上述情況，本文以新疆地區(qū)為研究區(qū)，綜合利用MODIS逐日積雪覆蓋產(chǎn)品MOD10A1和MYD10A1數(shù)據(jù)、IMS數(shù)據(jù)、新疆地區(qū)數(shù)字高程數(shù)據(jù)(digital elevation model,DEM)、新疆109個(gè)常規(guī)氣象站數(shù)據(jù)和野外實(shí)驗(yàn)人工實(shí)測(cè)數(shù)據(jù)，結(jié)合已有研究成果，擬發(fā)展一套適合于新疆地區(qū)的基于多源數(shù)據(jù)的積雪遙感數(shù)據(jù)去云算法，生成研究區(qū)2002—2016年間逐日無云積雪遙感產(chǎn)品數(shù)據(jù)集； 并對(duì)生成的結(jié)果數(shù)據(jù)進(jìn)行精度驗(yàn)證，給出精度評(píng)價(jià)驗(yàn)證結(jié)果。

1 數(shù)據(jù)源及其預(yù)處理

1.1 光學(xué)積雪遙感產(chǎn)品數(shù)據(jù)

本文采用的光學(xué)積雪遙感產(chǎn)品數(shù)據(jù)來自美國(guó)國(guó)家雪冰數(shù)據(jù)中心(National Snow and Ice Data Center,NSIDC)網(wǎng)站(https: //ladsweb.nascom.nasa.gov/)提供的MODIS/Terra和MODIS/Aqua衛(wèi)星日積雪產(chǎn)品數(shù)據(jù)(MOD10A1和MYD10A1)。此類積雪產(chǎn)品是基于歸一化積雪指數(shù)(normalized difference snow index，NDSI)閾值法生成的全球每日積雪覆蓋產(chǎn)品[27]，是目前積雪研究中應(yīng)用最為廣泛的長(zhǎng)時(shí)間序列數(shù)據(jù)。已有學(xué)者對(duì)MODIS積雪制圖算法的精度進(jìn)行了驗(yàn)證，結(jié)果表明在晴空狀態(tài)下上述2種產(chǎn)品的積雪識(shí)別率達(dá)到90%[14,26]。目前全球發(fā)布的MODIS積雪制圖算法的版本有第四版本(V004)和第五版本(V005)。本文采用V005版本算法生成的每日積雪覆蓋產(chǎn)品，空間分辨率為500 m，時(shí)間分辨率為1 d，數(shù)據(jù)文件為hdf4格式，投影方式為Inergrized Sinosoida GRID。軌道編號(hào)分別為h23v04,h23v05,h24v04,h24v05,h25v04和h25v05 的6塊數(shù)據(jù)覆蓋了整個(gè)研究區(qū)。本文下載了2002年8月1日—2016年8月1日共15 a間的所有可下載數(shù)據(jù)，其中，MOD10A1共30 602景(缺失82景)，MYD10A1共30 648景(缺失36景)； 對(duì)于同一天的數(shù)據(jù)，如果有數(shù)據(jù)缺失，就只采用1個(gè)傳感器的數(shù)據(jù)參與去云。通過編程，對(duì)所有數(shù)據(jù)自動(dòng)進(jìn)行拼接與投影轉(zhuǎn)換(轉(zhuǎn)換為經(jīng)緯度投影)，得到覆蓋研究區(qū)(E73°～95°，N34°～50°)的hdf4格式積雪產(chǎn)品數(shù)據(jù)。由于MODIS數(shù)據(jù)有時(shí)存在數(shù)據(jù)缺失及無法判識(shí)現(xiàn)象，因此把這一類數(shù)據(jù)也定義為疑似云像元，進(jìn)行后續(xù)去云處理； 并結(jié)合IMS數(shù)據(jù)的編碼意義，對(duì)MODIS數(shù)據(jù)進(jìn)行重新分類，結(jié)果如表1所示。

表1 MODIS標(biāo)準(zhǔn)積雪產(chǎn)品類別及重新劃分標(biāo)準(zhǔn)Tab.1 MODIS standardized snow cover productclasses and reclassification rules

1.2 IMS雪冰產(chǎn)品數(shù)據(jù)

IMS雪冰產(chǎn)品數(shù)據(jù)由美國(guó)國(guó)家海洋和大氣管理局(Nation Oceanic and Atmospheric Administration，NOAA)制作，通過NSIDC網(wǎng)站(http: //nsidc.org/data/docs/noaa/g02156_ims_snow_ ice_analysis/)可以免費(fèi)獲取。該產(chǎn)品是由NOAA極地軌道衛(wèi)星、地球靜止軌道衛(wèi)星、日本GMS同步氣象衛(wèi)星、日本MTSAT多功能衛(wèi)星、歐洲同步氣象衛(wèi)星、美國(guó)國(guó)防部極地衛(wèi)星、美國(guó)國(guó)防部氣象衛(wèi)星計(jì)劃、NIC每周海冰分析產(chǎn)品、美國(guó)空軍雪深/雪冰分析產(chǎn)品和微波傳感器數(shù)據(jù)等衛(wèi)星數(shù)據(jù)及產(chǎn)品融合而成，并通過分析師判識(shí)確定的多源數(shù)據(jù)雪冰覆蓋產(chǎn)品[27-28]，因此IMS數(shù)據(jù)不受云的影響。文獻(xiàn)[28]得出IMS數(shù)據(jù)在新疆積雪的年平均總體準(zhǔn)確率超過了92%，積雪季平均總體準(zhǔn)確率均超過了88%。IMS雪冰產(chǎn)品采用以北極點(diǎn)為中心的極地方位投影，其產(chǎn)品包括ASCII格式和TIFF格式，時(shí)間分辨率為1 d，空間分辨率有24 km(自1997年至今)、4 km(自2004年至今)和1 km(自2014年至今)3種。不同的像元值代表不同的地物： 1為海洋，2為陸地，3為海冰與湖冰，4為積雪。本文下載了2002年8月1日—2014年2月23日的24 km空間分辨率和2014年2月24日—2016年8月1日的1 km空間分辨率的數(shù)據(jù)，并通過編程對(duì)研究區(qū)范圍內(nèi)的數(shù)據(jù)進(jìn)行提取與投影轉(zhuǎn)換。轉(zhuǎn)換后的投影方式為經(jīng)緯度投影，數(shù)據(jù)格式為hdf4。

1.3 DEM數(shù)據(jù)

DEM選取的是美國(guó)地質(zhì)調(diào)查局(United States Geological Survey，USGS)地球資源觀測(cè)系統(tǒng)(the earth resources observation systems，EROS)數(shù)據(jù)中心(ftp: //eosrp01u. ecs.nasa.gov/)發(fā)布的空間分辨率為90 m的SRTM DEM數(shù)據(jù)。為了與積雪產(chǎn)品匹配，將其空間分辨率重采樣為500 m。

1.4 實(shí)測(cè)雪深數(shù)據(jù)

實(shí)測(cè)雪深數(shù)據(jù)采用新疆氣象局提供的新疆109個(gè)常規(guī)氣象站2002—2016年間觀測(cè)的日值數(shù)據(jù)。氣象站位置如圖1所示。

圖1 研究區(qū)及氣象臺(tái)站位置

1.5 氣象數(shù)據(jù)

氣象數(shù)據(jù)包括測(cè)站基本信息、氣壓、溫度、降水、風(fēng)、地溫、人工觀測(cè)數(shù)據(jù)和其他重要天氣數(shù)據(jù)，其中人工觀測(cè)數(shù)據(jù)包括能見度、云等數(shù)據(jù)，其他重要天氣數(shù)據(jù)包括積雪深度(0.1 cm)、雪壓(0.1 g/cm2)、凍土、龍卷風(fēng)、電線積冰和冰雹直徑等信息。每個(gè)氣象站日值為一個(gè)數(shù)據(jù)樣本，若雪深大于0 cm，表示有積雪覆蓋； 若雪深為0 cm，表示無雪樣本。

1.6 野外實(shí)測(cè)數(shù)據(jù)

野外實(shí)測(cè)數(shù)據(jù)是指近5 a來積雪期在研究區(qū)野外積雪觀察實(shí)驗(yàn)中所測(cè)積雪信息數(shù)據(jù)，包括觀察時(shí)間、觀察地點(diǎn)經(jīng)緯度、天氣、SnowFork數(shù)據(jù)、積雪光譜數(shù)據(jù)、雪深(cm)、雪壓(g/cm2)、雪層溫度(℃)、雪粒徑大小、雪層情況、雪層底地表特性、積雪覆蓋度、地形和環(huán)境照片等信息。野外實(shí)測(cè)數(shù)據(jù)作為隨機(jī)數(shù)據(jù)，與氣象站數(shù)據(jù)一同作為精度驗(yàn)證數(shù)據(jù)。

2 研究方法

本文綜合已有去云方法的優(yōu)點(diǎn)，引入IMS數(shù)據(jù)和DEM數(shù)據(jù)，分析不同影像的空間和時(shí)間信息，設(shè)計(jì)了一套適合于研究區(qū)的去云算法，算法流程如圖2所示。

圖2 去云算法流程

2.1 上下午星結(jié)合去云

同一天的MODIS/Terra衛(wèi)星與MODIS/Aqua衛(wèi)星在新疆區(qū)域成像時(shí)間相差4 h左右，這段時(shí)間內(nèi)云會(huì)產(chǎn)生變化和發(fā)生移動(dòng)，而云下積雪相對(duì)變化較小。利用2顆衛(wèi)星的這一時(shí)間差，有可能將原來被云遮擋的積雪信息識(shí)別出來。在長(zhǎng)期的業(yè)務(wù)運(yùn)行中發(fā)現(xiàn)以及文獻(xiàn)[24]表明，在新疆地區(qū)MODIS/Terra衛(wèi)星資料的質(zhì)量要明顯高于MODIS/Aqua，故本文以MODIS/Terra衛(wèi)星的積雪產(chǎn)品MOD10A1為初步去云數(shù)據(jù)。當(dāng)MOD10A1重新編碼后像元值為50時(shí)，進(jìn)行如下操作，否則保持不變： 如果相同位置MYD10A1產(chǎn)品數(shù)據(jù)重新編碼后像元值也為50，則保持MOD10A1值不變； 如果相同位置MYD10A1產(chǎn)品數(shù)據(jù)的值為非50，則用MYD10A1值替代相同位置MOD10A1的值，上下午星結(jié)合去云最終結(jié)果的產(chǎn)品記作MOYD10A1。

2.2 前后日合成去云

前后日合成去云通過考慮云體在不同時(shí)相的位置差異和在這一時(shí)段的降融雪過程，獲取當(dāng)日積雪信息。此方法最大的缺點(diǎn)是隨著合成日數(shù)的增加，對(duì)持續(xù)較短的積雪時(shí)間探測(cè)能力有所下降，而新疆積雪除了降融雪時(shí)期外，整個(gè)積雪期的積雪都相對(duì)穩(wěn)定，這為此方法的應(yīng)用提供了可能?？紤]到結(jié)果數(shù)據(jù)應(yīng)用的時(shí)效性及研究區(qū)積雪存在的機(jī)理，本文采用探測(cè)前一日(d-1)與后一日(d+1)和當(dāng)日(d)MOYD10A1產(chǎn)品中被云遮擋的積雪信息，具體的積雪判識(shí)方法如表2所示。

表2 前后日合成去云判識(shí)Tab.2 Cloud removal judgement by combination of beforeand after day

判識(shí)過程為： 類型1表示d±1日晴空有雪，則將有云d日判識(shí)為有雪； 類型5表示d±1日晴空無雪，則有云d日判識(shí)為無雪； 類型4表示d-1日晴空無雪，d+1日晴空有雪，則有云d日表示有降雪過程，判識(shí)為有雪； 類型7表示d-1日有云，d+1日晴空有雪，即連續(xù)2日有云后的晴天監(jiān)測(cè)到有雪，說明晴天之前有積雪存在，表示d日有雪； 其他類型尚無法判識(shí)，待后續(xù)進(jìn)行判識(shí)。前后日合成去云最終結(jié)果的產(chǎn)品記做MOYD10A1_T。

2.3 空間濾波法去云

因?yàn)榉e雪具有一定的地帶分布規(guī)律，在同一區(qū)域的相同地理環(huán)境下，積雪的分布具有一定的相似性(即相近相似原理)。空間濾波法去云即依據(jù)積雪在局部區(qū)域的空間分布相似性原理，消除由于IMS數(shù)據(jù)空間分辨率比MOYD10A1_T低而引起的對(duì)大面積云體邊緣及零星云體的云下信息的誤判。具體判識(shí)過程為： 選擇窗口大小為9像元×9像元空間濾波器，能夠判識(shí)的大面積云體邊緣及零星云體像元類別，如圖3所示(C為云像元，S為積雪像元，L為陸地和水體等其他晴空非雪像元)。

(a) 類別1 (b) 類別2(c) 類別3

圖3濾波窗口像元類別

Fig.3Pixeltypesinfilteringwindow

圖3中，類別1濾波窗口內(nèi)僅存在L和C像元，對(duì)比窗口內(nèi)每個(gè)像元的DEM，若某C像元的DEM小于所有L像元的平均DEM，則可判識(shí)這個(gè)C像元為最近的L像元； 否則表示無法判識(shí)。類別2濾波窗口內(nèi)僅存在S和C像元，對(duì)比窗口內(nèi)每個(gè)像元的DEM，若某C像元的DEM大于所有S像元的平均DEM，則可判識(shí)這個(gè)C像元為最近的S像元； 否則表示無法判識(shí)。類別3濾波窗口內(nèi)同時(shí)存在S，L與C像元，對(duì)比窗口內(nèi)每個(gè)像元的DEM，若某C像元的DEM大于所有S像元的平均DEM，則可判識(shí)這個(gè)C像元為最近的S像元； 若某C像元的DEM小于所有L像元的平均DEM，則可以判識(shí)這個(gè)C像元為最近的L像元； 否則表示無法判識(shí)?？臻g濾波法去云最終結(jié)果的產(chǎn)品記作MOYD10A1_T_SF。

2.4 多源數(shù)據(jù)合成去云

雖然通過上述步驟可以消除一部分云量，但是對(duì)時(shí)間濾波窗口期內(nèi)沒有明顯變化的大面積云量還是無法去除； 所以本文引入不受云影響的IMS數(shù)據(jù)，對(duì)MOYD10A1_T_SF數(shù)據(jù)進(jìn)行最終去云處理。對(duì)重新編碼后MOYD10A1_T_SF數(shù)據(jù)中值為50的像元與IMS數(shù)據(jù)進(jìn)行融合，具體規(guī)則為： 若IMS=1，則表示云下信息為海洋，編碼為39； 若IMS=2，則表示云下信息為無雪覆蓋陸地，編碼為25； 若IMS=3，則表示云下信息為海冰與湖冰，編碼為37； 若IMS=4，則表示云下信息為積雪覆蓋，編碼為200。通過多源數(shù)據(jù)合成法去除以上方法無法去除的云體，最終得到無云的積雪產(chǎn)品MOYD10A1_T_SF_IMS。

3 結(jié)果驗(yàn)證

3.1 基于統(tǒng)計(jì)分析的去云效果評(píng)價(jià)

云覆蓋度Pc和雪覆蓋度Ps可以反映不同融合方法的去云效果，其表達(dá)式為

(1)

(2)

式中：Nc為影像中的云像元數(shù)；Ns為影像中的雪像元數(shù)；N為像元總數(shù)。

圖4以2013年3月9日MODIS積雪產(chǎn)品數(shù)據(jù)為例，給出了通過本文算法進(jìn)行逐步去云后所得產(chǎn)品的空間影像圖。圖5給出了各去云步驟所得產(chǎn)品數(shù)據(jù)近15 a間逐日平均Pc。

(a) MYD10A1 (b) MOD10A1 (c) MOYD10A1

(d) MOYD10A1_T(e) MOYD10A1_T_SF

圖4研究區(qū)MODIS積雪產(chǎn)品去云過程影像

Fig.4CloudremovalprocessimagesofMODISsnowcoverproductofstudyarea

圖5 研究區(qū)近15 a不同去云方法逐日平均云覆蓋度(2002—2016年)

對(duì)研究區(qū)近15 a的MOD10A1和MYD10A1數(shù)據(jù)進(jìn)行逐日云像元統(tǒng)計(jì)分析，發(fā)現(xiàn)研究區(qū)內(nèi)全年大部分時(shí)間的Pc占比在50%左右，尤其在積雪季(11月─翌年3月)，大部分時(shí)間逐日云覆蓋占比高達(dá)60%左右，如圖5中MOD10A1曲線與MYD10A1曲線所示。以2013年3月9日數(shù)據(jù)為例， MYD10A1的Pc可達(dá)68.35%，Ps僅為3.96%(圖4(a))； 而MOD10A1的Pc達(dá)65.15%，Ps為4.03%(圖4(b))。顯然，云覆蓋已成為影響MODIS數(shù)據(jù)在研究區(qū)進(jìn)行積雪監(jiān)測(cè)的最大制約因素，通過本文去云算法，最終解決了研究區(qū)MODIS積雪產(chǎn)品中的云覆蓋問題。

從逐日平均Pc統(tǒng)計(jì)結(jié)果可知，MOD10A1和MYD10A1數(shù)據(jù)近15 a間平均Pc分別為46.22%和48.12%(如圖5中MOD10A1曲線與MYD10A1曲線所示)。通過上下午星結(jié)合去云法得到的產(chǎn)品MOYD10A1，其近15 a間逐日平均Pc為34.26%(如圖5中MOYD10A1曲線所示)，與MOD10A1數(shù)據(jù)相比，Pc降低了11.96%。從單日統(tǒng)計(jì)結(jié)果可以看出，2013年3月9日MOYD10A1的Pc為56.02%，與同日MOD10A1數(shù)據(jù)相比，Pc降低了9.13%(圖4(c))。從圖4可以看出，上下午星結(jié)合去云主要是去掉了同一天內(nèi)上下午星過境時(shí)間差內(nèi)變化的云和移動(dòng)的云。前后日合成去云在幾種方法中去云量最多，去除云量的多少主要取決于d±1日內(nèi)云體的變化及其移動(dòng)的范圍，其結(jié)果產(chǎn)品MOYD10A1_T的近15 a間逐日年平均Pc為17.86%(如圖5中MOYD10A1_T曲線所示)，與MOYD10A1相比，Pc降低了16.40%。2013年3月9日單日Pc為40.22%(圖4(d))，減少了15.80%。從圖4(e)中的空間濾波結(jié)果產(chǎn)品空間影像可以看出，空間濾波法去云主要去除了空間零星云體，對(duì)云體—積雪、云體—晴空陸地等云體與其他非云體相接邊緣處的云體也有較好的去除效果，這為后續(xù)去云方法提供了較大幫助。結(jié)果產(chǎn)品MOYD10A1_T_SF近15 a間逐日年平均Pc為13.29%(如圖5中MOYD10A1_T_SF曲線所示)，與MOYD10A1_T相比，Pc減少了4.57%。2013年3月9日單日Pc為31.53%(圖4(e))，減少了8.69%。多源數(shù)據(jù)融合去云法結(jié)合IMS數(shù)據(jù)去除了以上步驟無法去除的剩余云體，最終得到研究區(qū)內(nèi)晴空積雪產(chǎn)品數(shù)據(jù)MOYD10A1_T_SF_IMS。與無云遮擋的IMS數(shù)據(jù)和原始MODIS積雪產(chǎn)品數(shù)據(jù)相比，該產(chǎn)品數(shù)據(jù)在保持了原有數(shù)據(jù)時(shí)空分辨率的同時(shí)，得到了研究區(qū)內(nèi)無云遮擋的晴空積雪覆蓋數(shù)據(jù)，圖6給出了研究區(qū)2013年3月9日IMS數(shù)據(jù)與MOYD10A1_T_SF_IMS數(shù)據(jù)積雪空間分布。

(a) IMS (b) MOYD10A1_T_SF_IMS

圖6研究區(qū)晴空積雪覆蓋空間分布

Fig.6Spacedistributionofsnowcoverinclearskyinstudyarea

對(duì)研究區(qū)內(nèi)近15 a間去云前MOD10A1逐日數(shù)據(jù)和去云后MOYD10A1_T_SF_IMS逐日無云數(shù)據(jù)進(jìn)行了Ps統(tǒng)計(jì)。去云前因受云的影響，研究區(qū)內(nèi)逐日Ps呈現(xiàn)出比較劇烈的波動(dòng)，如圖7(a)所示研究區(qū)代表性年份逐日Ps曲線； 去云后變化趨于平穩(wěn)，恢復(fù)了自然狀態(tài)積雪變化的波動(dòng)，如圖7(b)所示研究區(qū)代表性年份逐日Ps曲線。統(tǒng)計(jì)近15 a間各去云方法的逐日平均Ps，結(jié)果如圖7(c)所示，去云前，MYD10A1和MOD10A1的逐日年平均Ps分別為6.39%和6.99%； 通過各去云方法依次去云后，得到的研究區(qū)各產(chǎn)品數(shù)據(jù)MOYD10A1，MOYD10A1_T，MOYD10A1_T_SF和MOYD10A1_T_SF_IMS的逐日年平均Ps分別為8.97%，12.54%，13.94%和19.31%，隨著云量的減少，積雪覆蓋逐漸增加，直至接近真實(shí)情況。

(a) 去云前

(b) 去云后

(c) 不同去云方法的逐日平均Ps

圖7研究區(qū)近15a各去云方法逐日積雪覆蓋度(2002—2016年)

Fig.7Dailyaveragesnowcoverproportionbydifferentcloudremovalmethodsfrom2002to2016instudyarea

3.2 基于實(shí)測(cè)數(shù)據(jù)的去云精度評(píng)價(jià)

MODIS積雪產(chǎn)品精度是指在晴空狀態(tài)下進(jìn)行的精度驗(yàn)證，對(duì)于云覆蓋的區(qū)域，應(yīng)用本文算法去云后，與去云前晴空狀態(tài)數(shù)據(jù)相比，原來為云覆蓋的區(qū)域在MOYD10A1_T_SF_IMS數(shù)據(jù)中表現(xiàn)為新增晴空狀態(tài)數(shù)據(jù)。選擇研究區(qū)內(nèi)109個(gè)常規(guī)氣象站近15 a間的日值數(shù)據(jù)，分別對(duì)近15 a間去云前MODIS晴空狀態(tài)下積雪產(chǎn)品數(shù)據(jù)(表3)和去云后新增晴空狀態(tài)下積雪產(chǎn)品數(shù)據(jù)(表4)進(jìn)行精度驗(yàn)證。

表3 去云前MODIS數(shù)據(jù)與氣象站點(diǎn)積雪觀察誤差矩陣Tab.3 Error matrix of MODIS data before cloud removaland in-situ observation data at meteorological station

表4 去云后MODIS數(shù)據(jù)與氣象站點(diǎn)積雪觀察誤差矩陣Tab.4 Error matrix of MODIS data after cloud removaland in-situ observation data at meteorological station

去云前晴空狀態(tài)下MODIS積雪產(chǎn)品數(shù)據(jù)能與氣象站觀察數(shù)據(jù)匹配的數(shù)據(jù)共有377 472條記錄，實(shí)測(cè)數(shù)據(jù)共有876條記錄，共288 348對(duì)數(shù)據(jù)，精度結(jié)果驗(yàn)證(表3)表明，積雪一致性為91.5%，陸地一致性為95.1%，總體精度為93.3%，說明在晴空狀態(tài)下，MODIS積雪產(chǎn)品數(shù)據(jù)在研究區(qū)有較高的精度。從分析結(jié)果也可以看出MODIS積雪產(chǎn)品對(duì)研究區(qū)積雪覆蓋存在11.1%高估現(xiàn)象。

去云后對(duì)MOYD10A1_T_SF_IMS產(chǎn)品中新增晴空狀態(tài)數(shù)據(jù)精度驗(yàn)證的實(shí)測(cè)數(shù)據(jù)應(yīng)同時(shí)滿足如下要求： ①去云前MOD10A1數(shù)據(jù)云覆蓋下的氣象臺(tái)站； ②氣象臺(tái)站云總量記錄值超過60%； ③氣象臺(tái)站積雪深度記錄值超過1 cm。符合條件的數(shù)據(jù)共有234 137條記錄，符合條件的野外實(shí)測(cè)數(shù)據(jù)有974條記錄，共235 111對(duì)數(shù)據(jù)，精度驗(yàn)證(表4)表明，積雪一致性為88.12%，陸地一致性為93.11%，總體精度為90.61%。從分析結(jié)果也可以看出，去云后的積雪產(chǎn)品對(duì)研究區(qū)積雪覆蓋存在著15.62%高估現(xiàn)象，說明原來被云覆蓋的區(qū)域經(jīng)過本文算法去云后，生成的研究區(qū)逐日無云積雪產(chǎn)品MOYD10A1_T_SF_IMS的精度基本接近MODIS晴空積雪產(chǎn)品數(shù)據(jù)在研究區(qū)的積雪監(jiān)測(cè)精度，驗(yàn)證了本文算法對(duì)研究區(qū)MODIS數(shù)據(jù)積雪產(chǎn)品去云的有效性。

4 結(jié)論

1)MODIS每日積雪產(chǎn)品一般存在著嚴(yán)重的云覆蓋，大大影響了對(duì)積雪覆蓋的實(shí)時(shí)監(jiān)測(cè)及應(yīng)用； 通過本文去云算法處理之后，在保證產(chǎn)品時(shí)空分辨率不變的前提下，解決了MODIS積雪產(chǎn)品在研究區(qū)的云覆蓋問題。通過對(duì)去云后產(chǎn)品的精度驗(yàn)證，得出去云后產(chǎn)品總體精度為90.61%，接近MODIS原始數(shù)據(jù)在研究區(qū)的晴空積雪總體精度(93.3%)，說明經(jīng)本文算法去云后的逐日無云積雪產(chǎn)品可以用于研究區(qū)積雪監(jiān)測(cè)。通過建立2002—2016年近15 a間研究區(qū)逐日無云數(shù)據(jù)集，為后續(xù)新疆積雪研究提供了較高空間分辨率的逐日積雪覆蓋產(chǎn)品數(shù)據(jù)。

2)不同的去云方法依據(jù)不同原理進(jìn)行去云，因而造成的誤差也不同。在本文的4步去云法中，第2步前后日合成去云法會(huì)在積雪不穩(wěn)定區(qū)、降雪期或融雪期相對(duì)帶來更大誤差； 而新疆積雪除了降、融雪期，其他時(shí)期比較穩(wěn)定，這為本算法在新疆地區(qū)的應(yīng)用創(chuàng)造了條件。第4步多源數(shù)據(jù)合成去云法由于不同數(shù)據(jù)空間分辨率的差異，在去云過程中造成了因空間分辨率不同而產(chǎn)生的誤差。

3)云覆蓋是影響高時(shí)空分辨率光學(xué)積雪遙感產(chǎn)品應(yīng)用的瓶頸。利用高時(shí)間分辨率和積雪存在特點(diǎn)，并結(jié)合多源積雪產(chǎn)品數(shù)據(jù)，可以達(dá)到去云效果。通過改進(jìn)本文算法，也可適應(yīng)于其他高時(shí)間分辨率光學(xué)積雪遙感產(chǎn)品的去云。但去云算法只是一種手段，如何采用高時(shí)空分辨率且不受云影響的遙感資料直接進(jìn)行積雪監(jiān)測(cè)是未來發(fā)展的趨勢(shì)。

參考文獻(xiàn)(References)：

[1] 劉 艷,汪 宏,張 璞,等.MODIS大氣校正精度評(píng)價(jià)及其對(duì)表層雪密度提取影響[J].國(guó)土資源遙感,2011,23(1):128-132.doi:10.6046/gtzyyg.2011.01.26.

Liu Y,Wang H,Zhang P,et al.Accuracy evaluation of MODIS atmospheric correction and its effects on surface-snow density extraction[J].Remote Sensing for Land and Resources,2011,23(1):128-132.doi:10.6046/gtzyyg.2011.01.26.

[2] 馬麗云,李建剛,李 帥.基于FY-3/MERSI數(shù)據(jù)的新疆融雪性洪水災(zāi)害監(jiān)測(cè)[J].國(guó)土資源遙感,2015,27(4):73-78.doi:10.6046/gtzyyg.2015.04.12.

Ma L Y,Li J G,Li S.Snowmelt flood disaster monitoring based on FY-3/MERSI in Xinjiang[J].Remote Sensing for Land and Resources,2015,27(4):73-78.doi:10.6046/gtzyyg.2015.04.12.

[3] 胡列群,武鵬飛,梁鳳超,等.新疆冬春季積雪及溫度對(duì)凍土深度的影響分析[J].冰川凍土,2014,36(1):48-54.

Hu L Q,Wu P F,Liang F C,et al.Analyzing the effect of snow cover in spring and winter and air temperature on frozen ground depth in Xinjiang[J]Journal of Glaciology and Geocryology,2014,36(1):48-54.

[4] 白淑英,吳 奇,史建橋,等.青藏高原積雪深度時(shí)空分布與地形的關(guān)系[J].國(guó)土資源遙感,2015,27(4):171-178.doi:10.6046/gtzyyg.2015.04.26.

Bai S Y,Wu Q,Shi J Q,et al.Relationship between the spatial and temporal distribution of snow depth and the terrain over the Tibetan Plateau[J].Remote Sensing for Land and Resources,2015,27(4):171-178.doi:10.6046/gtzyyg.2015.04.26.

[5] 李三妹,閆 華,劉 誠(chéng).FY-2C積雪判識(shí)方法研究[J].遙感學(xué)報(bào),2007,11(3):406-413.

Li S M,Yan H,Liu C.Study of snow detection using FY-2C satellite data[J].Journal of Remote Sensing,2007,11(3):406-413.

[6] 梁守真,施 平,邢前國(guó).MODIS NDVI時(shí)間序列數(shù)據(jù)的去云算法比較[J].國(guó)土資源遙感,2011,23(1):33-36.doi:10.6046/gtzyyg.2011.01.06.

Liang S Z,Shi P,Xing Q G.A comparison between the algorithms for removing cloud pixel from MODIS NDVI time series data[J].Remote Sensing for Land and Resources,2011,23(1):33-36.doi:10.6046/gtzyyg.2011.01.06.

[7] 孫知文,于鵬珊,夏 浪,等.被動(dòng)微波遙感積雪參數(shù)反演方法進(jìn)展[J].國(guó)土資源遙感,2015,27(1):9-15.doi:10.6046/gtzyyg.2015.01.02.

Sun Z W,Yu P S,Xia L,et al.Progress in study of snow parameter inversion by passive microwave remote sensing[J].Remote Sensing for Land and Resources,2015,27(1):9-15.doi:10.6046/gtzyyg.2015.01.02.

[8] Xie H J,Wang X W,Liang T G.Development and assessment of combined Terra and Aqua snow cover products in Colorado Plateau,USA and northern Xinjiang,China[J].Journal of Applied Remote Sensing,2009,3(1):033559.

[9] 王增艷，車 濤.中國(guó)干旱區(qū)積雪面積產(chǎn)品去云處理方法驗(yàn)證與評(píng)估[J].干旱區(qū)研究，2012,29(2)：312-319.

Wang Z Y,Che T.Validation and assessment of cloud obscuration reduction of snow cover products in arid areas in China[J].Arid Zone Research,2012,29(2):312-319

[10] 劉 艷,李 楊,張 璞.基于實(shí)測(cè)光譜和SRF的稀疏植被區(qū)MODIS積雪信息提取[J].國(guó)土資源遙感,2013,25(1):26-32.doi:10.6046/gtzyyg.2013.01.05.

Liu Y,Li Y,Zhang P.Extraction of snow cover information in sparse vegetation area based on spectral measurement and SRF by using MODIS data[J].Remote Sensing for Land and Resources,2013,25(1):26-32.doi:10.6046/gtzyyg.2013.01.05.

[11] 許劍輝,舒 紅,李 楊.基于MODIS雪蓋數(shù)據(jù)的北疆雪深多元非線性回歸克里金插值[J].國(guó)土資源遙感,2015,27(3):84-91.doi:10.6046/gtzyyg.2015.03.15.

Xu J H,Shu H,Li Y.Mapping of monthly mean snow depth in Northern Xinjiang using a multivariate nonlinear regression Kriging model based on MODIS snow cover data[J].Remote Sensing for Land and Resources,2015,27(3):84-91.doi:10.6046/gtzyyg.2015.03.15.

[12] 張瑞江,趙福岳,方洪賓,等.青藏高原近30年現(xiàn)代雪線遙感調(diào)查[J].國(guó)土資源遙感,2010,22(s1):59-63.doi:10.6046/gtzyyg.2010.s1.14.

Zhang R J,Zhao F Y,Fang H B,et al.Remote sensing survey of existing snowlines in the past 30 years in Qinghai-Tibet Plateau[J].Remote Sensing for Land and Resources,2010,22(s1):59-63.doi:10.6046/gtzyyg.2010.s1.14.

[13] 李光輝,王 成,習(xí)曉環(huán),等.機(jī)載LiDAR和高光譜數(shù)據(jù)融合提取冰川雪線[J].國(guó)土資源遙感,2013,25(3):79-84.doi:10.6046/gtzyyg.2013.03.14.

Li G H,Wang C,Xi X H,et al.Extraction of glacier snowline based on airborne LiDAR and hyperspectral data fusion[J].Remote Sensing for Land and Resources,2013,25(3):79-84.doi:10.6046/gtzyyg.2013.03.14.

[14] 黃曉東,郝曉華,王 瑋,等.MODIS逐日積雪產(chǎn)品去云算法研究[J].冰川凍土,2012,34(5):1118-1126.

Huang X D,Hao X H,Wang W,et al.Algorithms for cloud removal in MODIS daily snow products[J].Journal of Glaciology and Geocryology,2012,34(5):1118-1126.

[15] 唐志光,王 建,李弘毅,等.青藏高原MODIS積雪面積比例產(chǎn)品的精度驗(yàn)證與去云研究[J].遙感技術(shù)與應(yīng)用,2013,28(3):423-430.

Tang Z G,Wang J,Li H Y,et al.Accuracy validation and cloud obscuration removal of MODIS fractional snow cover products over Tibetan Plateau[J].Remote Sensing Technology and Application,2013,28(3):423-430.

[16] 張 歡,邱玉寶,鄭照軍,等.基于MODIS的青藏高原季節(jié)性積雪去云方法可行性比較研究[J].冰川凍土,2016,38(3):714-724.

Zhang H,Qiu Y B,Zheng Z J,et al.Comparative study of the feasibility of cloud removal methods based on MODIS seasonal snow cover data over the Tibetan Plateau[J].Journal of Glaciology and Geocryology,2016,38(3):714-724.

[17] 趙明洋,別 強(qiáng),何 磊,等.基于去云處理的祁連山積雪覆蓋遙感監(jiān)測(cè)研究[J].干旱區(qū)地理,2014,37(2):325-332.

Zhao M Y,Bie Q,He L,et al.Estimating snow cover accuracy from MODIS and AMSR-E with cloud removal methodology in Qilian Mountains[J].Arid Land Geography,2014,37(2):325-332.

[18] 張 穎,黃曉東,王 瑋,等.MODIS逐日積雪覆蓋率產(chǎn)品驗(yàn)證及算法重建[J].干旱區(qū)研究,2013,30(5):808-814.

Zhang Y,Huang X D,Wang W,et al.Validation and algorithm redevelopment of MODIS daily fractional snow cover products[J].Arid Zone Research,2013,30(5):808-814.

[19] 黃曉東,張學(xué)通,李 霞,等.北疆牧區(qū)MODIS積雪產(chǎn)品MOD10A1和MOD10A2的精度分析與評(píng)價(jià)[J].冰川凍土,2007,29(5):722-729.

Huang X D,Zhang X T,Li X,et al.Accuracy analysis for MODIS snow products of MOD10A1 and MOD10A2 in Northern Xinjiang Area[J].Journal of Glaciology and Geocryology,2007,29(5):722-729.

[20] Dankers R,Jong S M D.Monitoring snow-cover dynamics in Northern Fennoscandia with SPOT VEGETATION images[J].International Journal of Remote Sensing,2004,25(15):2933-2949.

[21] Hartman R K,Rost A A,Anderson D M.Operational processing of multi-source snow data[C]//Proceedings of the 63rd Annual Western Snow Conference.Sparks,Nevada:Western Snow Conference,1995:147-151.

[22] Hall D K,Riggs G A,Salomonson V V,et al.MODIS snow-cover products[J].Remote Sensing of Environment,2002,83(1/2):181-194.

[23] Parajka J,Pepe M,Rampini A,et al.A regional snow-line method for estimating snow cover from MODIS during cloud cover[J].Journal of Hydrology,2010,381(3/4):203-212.

[24] Gao Y,Xie H J,Lu N,et al.Toward advanced daily cloud-free snow cover and snow water equivalent products from Terra-Aqua MODIS and Aqua AMSR-E measurements[J].Journal of Hydrology,2010,385(1/4):23-35.

[25] 付文軒,沈煥鋒,李星華,等.時(shí)空自適應(yīng)加權(quán)的MODIS積雪產(chǎn)品去云方法[J].遙感信息,2016,31(2):36-43.

Fu W X,Shen H F,Li X H,et al.Adaptively spatio-temporal weighted method for removing cloud obscuration from MODIS daily snow cover products[J].Remote Sensing Information,2016,31(2):36-43.

[26] Hall D K,Riggs G A,Salomonson V V.Algorithm Theoretical Basis Document (ATBD) for the MODIS Snow and Sea Ice-Mapping Algorithms[R].Washington:NASA,2001.

[27] 侯小剛,張 璞,鄭照軍,等.基于多源數(shù)據(jù)的阿勒泰地區(qū)雪深反演研究[J].遙感技術(shù)與應(yīng)用,2015,30(1):178-185.

Hou X G,Zhang P,Zheng Z J,et al.Study of snow depth retrieval based on multi-source data about Aletai Area[J].Remote Sensing Technology and Application,2015,30(1):178-185.

[28] 劉 洵,金 鑫,柯長(zhǎng)青.中國(guó)穩(wěn)定積雪區(qū)IMS雪冰產(chǎn)品精度評(píng)價(jià)[J].冰川凍土,2014,36(3):500-507.

Liu X,Jin X,Ke C Q.Accuracy evaluation of the IMS snow and ice products in stable snow covers regions in China[J].Journal of Glaciology and Geocryology,2014,36(3):500-507.