2000-2021年黃河源區(qū)無云MODIS積雪覆蓋面積比例數(shù)據(jù)集

楊映，唐忠喜，邢德，侯金亮

1.甘肅洮河國家級自然保護區(qū)管護中心，甘肅甘南藏族自治州 747600

2.甘肅省地質(zhì)礦產(chǎn)勘查開發(fā)局第三地質(zhì)礦產(chǎn)勘查院，蘭州 730050

3.國防科技大學海洋氣象學院，長沙 050024

4.中國科學院西北生態(tài)環(huán)境資源研究院甘肅省遙感重點實驗室，蘭州 730000

引 言

2019年9月18日，習近平總書記在主持召開黃河流域生態(tài)保護和高質(zhì)量發(fā)展座談會時，將黃河流域生態(tài)保護和高質(zhì)量發(fā)展提升為重大國家戰(zhàn)略，特別強調(diào)了黃河流域的水源涵養(yǎng)問題。黃河源區(qū)位于青藏高原東北部，平均海拔為4026 m，海拔3000 m以上的區(qū)域占整個研究區(qū)的95.80%，年平均徑流量約為 204.0億立方米，占黃河流域年徑流量的三分之一，是黃河流域的主要產(chǎn)水區(qū)和水源涵養(yǎng)區(qū)，被譽為“黃河水塔”[1-2]。而青藏高原是北半球中低緯度海拔最高、積雪覆蓋最大的地區(qū)，既是氣候變化的敏感區(qū)，又對水資源系統(tǒng)產(chǎn)生重要影響[3-5]。已有研究表明，青藏高原東部的積雪變化最顯著，其積雪主導了整個青藏高原積雪的變化[6-8]。因此，獲取黃河源區(qū)準確的積雪覆蓋變化信息，具有十分重要的意義。

衛(wèi)星遙感是持續(xù)監(jiān)測積雪時空變化的有力工具[9]。光學傳感器利用積雪與云以及其他地物在可見光（VIS, 0.4-0.7 μm）和近紅外（NIR, 1-4 μm）波段的較大反射率差異特性，在全球、區(qū)域尺度上的積雪監(jiān)測中發(fā)揮了重要作用[10]。目前，國內(nèi)外學者利用甚高分辨率輻射計（Advanced Very High Resolution Radiometer, AVHRR）、中分辨率成像光譜儀（Moderate Resolution Imaging Spectroradiometer,MODIS）、可見光紅外成像輻射儀（Visible infrared Imaging Radiometer Suite, VIIRS）和風云三號氣象衛(wèi)星（FY-3）上搭載的中分辨率光譜成像儀（Medium Resolution Spectral Imager, MERSI）等光學傳感器研制了不同時間分辨率和空間分辨率的多種全球和區(qū)域尺度的積雪覆蓋產(chǎn)品[11-14]。這些積雪產(chǎn)品不僅可以用于監(jiān)測積雪覆蓋的時空變化，還可以作為各種大氣模式和陸面/水文過程模型的輸入?yún)?shù)。其中，AVHRR和FY-3積雪面積產(chǎn)品的空間分辨率較粗，在積雪斑塊化分布、消融變化劇烈以及地形復雜山區(qū)等情況下，難以精準反映積雪的空間分布細節(jié)特征。VIIRS和MODIS積雪產(chǎn)品非常相似，但VIIRS時間序列較短（2012年起）。簡而言之，已有的各種積雪覆蓋產(chǎn)品中，MODIS系列產(chǎn)品因其較優(yōu)的質(zhì)量、較高的時空分辨率（可達到500m，逐日）、全球覆蓋、長時間序列（2000年起）、免費使用等優(yōu)勢[9,15-16]，在積雪研究領域應用最廣。全球不同區(qū)域上的大量驗證結果表明，MODIS積雪產(chǎn)品具有很好的精度，晴空條件下的總體精度一般在 90%-96%之間[17-19]。然而，由于云覆蓋的影響，導致MODIS積雪產(chǎn)品中存在大量的數(shù)據(jù)缺失。已有研究表明，逐日MODIS積雪產(chǎn)品中由于云覆蓋造成的數(shù)據(jù)缺失在 39%-70%之間[20-21]。在所有天氣條件下，MODIS積雪產(chǎn)品的總體平均精度僅為34%-50%[22]。顯然，MODIS積雪產(chǎn)品中的數(shù)據(jù)缺失嚴重阻礙了其在水文模擬、氣候變化以及雪災預警預報等方面的應用。

為提高 MODIS積雪產(chǎn)品的監(jiān)測精度，大量學者對 MODIS積雪產(chǎn)品中的云覆蓋進行了重建研究，如黃曉東等[23]、YU等[24]和邱玉寶等[25]均研制了青藏高原區(qū)域的逐日無云二值積雪產(chǎn)品，HAO等[26]基于MODIS發(fā)展了中國逐日無云二值積雪產(chǎn)品。唐志光等[27]發(fā)展了青藏高原區(qū)域的MODIS積雪面積比例（FSC）產(chǎn)品。自2016年起，新版本的MODIS積雪產(chǎn)品不再直接提供二值積雪產(chǎn)品和FSC產(chǎn)品，而是提供歸一化積雪指數(shù)（NDSI）產(chǎn)品，因為利用NDSI值可以很方便地推算出二值積雪覆蓋和FSC值[13]。目前，已有研究基于相似像元/塊的思想，對MODIS NDSI積雪產(chǎn)品中的云覆蓋進行了重建[28-30]，但受青藏高原云系旺盛、積雪斑塊化嚴重以及積雪消融較快等因素影響，在該區(qū)域搜索到合適的相似像元/塊并非易事，有時甚至根本就不存在合理的相似像元。XING等[31]利用深度學習技術深入挖掘時空特征，充分利用MODIS NDSI積雪產(chǎn)品中可用的時空信息，構建了基于部分卷積神經(jīng)網(wǎng)絡（Partial Convolution Neural Network，PCNN）的MODIS NDSI積雪產(chǎn)品云像元重建模型，并在黃河源區(qū)進行了驗證，結果表明PCNN模型重建后的MODIS NDSI產(chǎn)品達到了可靠的較高精度。

本數(shù)據(jù)集利用2000-2021年逐日的MODIS NDSI產(chǎn)品，使用上述基于PCNN的MODIS NDSI云像元重建模型，生成了時空連續(xù)的MODIS NDSI數(shù)據(jù)產(chǎn)品；其次，再利用NASA FSC產(chǎn)品的標準算法，制備了黃河源區(qū)2000-2021年無云逐日MODIS FSC遙感監(jiān)測數(shù)據(jù)集，將為黃河源區(qū)的積雪分布、雪水儲量估計、積雪變化分析和雪災風險評估等研究工作提供數(shù)據(jù)支撐。

1 數(shù)據(jù)采集和處理方法

1.1 數(shù)據(jù)生產(chǎn)流程

黃河源區(qū)MODIS FSC數(shù)據(jù)集的制作流程主要分為5個步驟：獲取MODIS積雪產(chǎn)品、數(shù)據(jù)預處理、MODIS NDSI云像元重建（包括：MODIS/Terra和MODIS/Aqua數(shù)據(jù)合成、臨近三天時間濾波和基于PCNN的云像元時空重建）、FSC估計和精度驗證，具體流程如圖1所示。

圖1 數(shù)據(jù)處理流程Figure 1 Data processing flow

1.2 數(shù)據(jù)獲取

本研究使用的是MODIS/Terra和MODIS/Aqua逐日500 m積雪覆蓋L3級產(chǎn)品（MOD10A1和MYD10A1, V6）。從NASA的雪冰數(shù)據(jù)中心（https://nsidc.org/）免費下載了研究區(qū)2000-2021年的數(shù)據(jù)。覆蓋整個黃河源區(qū)需要二個瓦片（即：h25v05和h26v06），影像數(shù)據(jù)均為HDF文件存儲格式、Sinusoidal投影。

雪深驗證數(shù)據(jù)來自中國氣象局（http://data.cma.cn）地面臺站觀測的2000-2020年每日地面氣候積雪資料數(shù)據(jù)集，雪深觀測以厘米（cm）為單位，取整數(shù)，雪深小于1 cm的記錄為無雪。黃河源區(qū)上6個氣象站記錄的雪深資料被用來對積雪覆蓋比例數(shù)據(jù)集進行基于離散值的精度評估。

1.3 數(shù)據(jù)預處理

利用MODIS MRT批處理工具對兩個瓦片進行拼接、重投影、坐標轉(zhuǎn)換、重采樣等相關預處理，將原始數(shù)據(jù)轉(zhuǎn)換為地理坐標系，0.005°空間分辨率。本研究選取Raw NDSI和NDSI snow cover作為研究數(shù)據(jù)。在Raw NDSI的值在-10000-10000，轉(zhuǎn)換系數(shù)為0.0001。NDSI snow cover數(shù)據(jù)的編碼為整數(shù)（0，10-100），其他類型包括云、內(nèi)陸水體、海洋、缺失值、無決策、夜間和傳感器飽和數(shù)據(jù)。由于通過Raw NDSI的值無法判斷出一個像元是否為有效像元，本研究利用NDSI snow cover數(shù)據(jù)中的分類，將NDSI snow cover數(shù)據(jù)中所有其他類型都重分類為云，再使用重分類后的NDSI snow cover數(shù)據(jù)對Raw NDSI數(shù)據(jù)進行云掩膜，并將Raw NDSI實際值小于0的像元的NDSI值設置為0，將處理后的Raw NDSI數(shù)據(jù)作為基礎數(shù)據(jù)層。

1.4 MODIS NDSI積雪產(chǎn)品云像元重建

1.4.1 MOD10A1和MYD10A1數(shù)據(jù)的逐日合成

根據(jù)云會移動的特點，對上午星的Terra/MODIS和下午星的Aqua/MODIS的NDSI數(shù)據(jù)進行逐日合成處理，采用“最大化”合成準則，具體是將一日內(nèi)兩幅NDSI影像進行對比，每個像元NDSI值取兩者中的較大值，合成后的數(shù)據(jù)簡稱為MOYD產(chǎn)品。

1.4.2 鄰近時間濾波

根據(jù)積雪降落之后，會在地表保留一段時間的特點，對逐日合成后的MOYD積雪產(chǎn)品進行時間上的濾波分析?？紤]到黃河源區(qū)存在大量斑塊狀積雪、消融變化快的特點，選擇最嚴格的時間窗口限制，即鄰近三天時間濾波。具體方法是：當一個像元被云覆蓋時，若該像元前一天和后一天均未被云覆蓋，取前、后兩天NDSI的平均值，作為該云像元的NDSI值。

1.4.3 基于PCNN的云像元重建

使用XING等[31]提出的基于端到端深度學習的三維PCNN的NDSI預測模型，利用時空信息估計云像元的地表覆蓋信息。首先，選擇塊的大小為 64×64，在研究區(qū)上遍歷，尋找不包含云的完整NDSI塊，研究區(qū)上2000-2021年間21個積雪季（即前一年11月1日至后一年4月30日）上共搜索到103630個完整的塊。其次，隨機選擇103630個有云的塊（云覆蓋率0%-10%，……，90%-100%各占1/10），以七天為時間窗口（即前、后各三天），提取完整塊七天的數(shù)據(jù)，形成時空信息塊組，并用對應有云塊組對完整塊組進行云掩膜。再次，將103 630個塊組，按7:3劃分為訓練集和測試集（即72 541個訓練塊組，31 089個測試塊組），在訓練集上訓練PCNN模型，獨立測試集上的數(shù)據(jù)用于驗證模型的性能。最后，將所構建的模型應用到整個黃河源區(qū)上，重建所有缺失像元的NDSI值，得到黃河源區(qū)2000-2021年時空連續(xù)的NDSI數(shù)據(jù)集。

1.5 MODIS FSC估計

使用NASA第五版MODIS FSC的標準算法[15]，根據(jù)像元的NDSI值推算其FSC值，公式如下：

根據(jù)公式（1）分別利用Terra/MODIS原始NDSI產(chǎn)品、Aqua/MODIS原始NDSI產(chǎn)品、雙星合成后的NDSI數(shù)據(jù)、臨近三天時間濾波的NDSI數(shù)據(jù)以及PCNN模型重建的時空連續(xù)NDSI數(shù)據(jù)計算FSC值，得到對應的FSC產(chǎn)品，分別簡稱為MOD、MYD、MOYD、ATF和PCNN產(chǎn)品。

1.6 精度驗證

1.6.1 基于站點雪深觀測值的驗證

利用黃河源區(qū)上6個氣象觀測站點（瑪多、達日、河南、久治、若爾蓋、紅原）2000年1月1日至2020年3月31日觀測的雪深值為“真值”，對計算出的MODIS FSC進行評估。提取6個站點對應像元的MODIS FSC值，與站點實測的SD值進行比較，構建混淆矩陣，如表1。

表1 MODIS FSC和站點觀測雪深的對比混淆矩陣Table 1 Comparative confusion matrix between MODIS FSC value and snow depth observations from stations

3個驗證指標總體精度（OA）、高估積雪事件（MO）和低估積雪事件（MU）分別定義如下：

1.6.2 基于云假設的驗證

本研究進一步基于云假設，即用“假想”的云覆蓋，對原始晴空 NDSI影像進行掩膜，再利用PCNN方法對云像元重建后，對FSC的估計結果進行精確的定量評價，定義兩個定量評價指標平均絕對偏差（MAE）和相關系數(shù)（R）如下：

其中，N是總的云像元數(shù)；xis是第i個假設“云”像元FSC的估計值；xiT是對應的參考“真實”FSC值；?和?是所有xis和xiT的平均值。

2 數(shù)據(jù)樣本描述

黃河源區(qū)2000-2021年無云積雪覆蓋比例數(shù)據(jù)集包含從2000-2001積雪季開始（即2000年11月1日）至2020-2021積雪季結束（即2021年4月30日）的7485個逐日的柵格數(shù)據(jù)文件。為便于數(shù)據(jù)在通用遙感軟件下處理與展示，將數(shù)據(jù)存儲為tif格式，投影方式為WGS84地理坐標系，文件命名規(guī)則為：YYYYDDD.tif，其中 YYYY代表年，DDD代表儒略日（001-365），空間分辨率為0.005°，屬性值FSC為像元上積雪覆蓋的百分比，覆蓋范圍為整個黃河源區(qū)。積雪覆蓋比例產(chǎn)品值的取值范圍均為0-100，填充值為200（即黃河源區(qū)范圍之外像元上的值），單位：無。注意，由于MODIS 原始數(shù)據(jù)的缺失（2000301、2001166-2001184、2001221、2001310、2001339、2002079-2002087以及2002105共33天無數(shù)據(jù)），故2000年11月1日至2021年4月30日的7518天中，共有7485天有數(shù)據(jù)。

圖2展示了2018-2019積雪季上6個樣例日MODIS FSC的空間分布狀況?？梢钥闯?，本無云MODIS FSC數(shù)據(jù)集從視覺上看，積雪分布變化平滑，黃河源區(qū)的積雪主要分布在中部的阿尼瑪青山和南部的巴彥喀拉爾山主峰；同時，黃河源區(qū)上分布著大量斑塊狀的較淺積雪，這與黃河源區(qū)上的實際情況相符合。積雪季上不同時間的 FSC空間分布很好地展示出了積雪的積累-消融變化過程，一些積雪事件（如2018年11月15日以及2019年2月底的新降雪），MODIS FSC產(chǎn)品都有很好的及時響應。

圖2 MODIS FSC分布圖。(a) 2018年10月15日；(b) 2018年11月15日；(c) 2018年12月15日；(d) 2019年1月15日；(e) 2019年2月15日；(f) 2019年4月15日Figure 2 MODIS FSC distribution map.(a) October 15, 2018; (b) November 15, 2018; (c) December 15, 2018; (d)January 15, 2019; (e) January 15, 2019; (f) April 15, 2019

3 數(shù)據(jù)質(zhì)量控制和評估

3.1 基于站點雪深觀測值的評估

基于站點雪深觀測的驗證，是最直接的驗證方式。利用站點實測雪深值的驗證結果如表2所示。表中平均云覆蓋日數(shù)百分比表示6個站點平均云覆蓋日數(shù)占總?cè)諗?shù)的百分比?？梢钥闯?，由于云覆蓋造成的嚴重數(shù)據(jù)缺失，導致MODIS原始產(chǎn)品超過一半的時間被云覆蓋。雙星合成和臨近時間窗口濾波，均可以將 MODIS FSC產(chǎn)品的云覆蓋降低 10%左右，PCNN重建所有云像元后，時空連續(xù)MODIS FSC產(chǎn)品的總體精度顯著提高，可以達到94%，并且跟原始產(chǎn)品晴空像元條件下比較，并沒有顯著增加高估積雪事件和低估積雪事件發(fā)生的概率（高估和低估均增加1%）。

表2 基于站點觀測雪深值的驗證結果Table 2 Performance evaluation results based on snow depth observations from stations

3.2 基于云假設的評估

受黃河源區(qū)地形復雜、氣象站點稀少且站點主要分布在低海拔平地（只有瑪多站位于海拔4000 m以上的區(qū)域）等因素影響，導致站點雪深對積雪信息描述存在較大不確定性?；谡军c雪深值的驗證，可能并不能準確代表整個研究區(qū)上的情況?；谠蒲谀さ尿炞C，是對基于 PCNN的 MODIS FSC估計模型的精確定量評價，結果如圖3所示。可以看出，PCNN模型重建的MODIS FSC值的平均MAE為10.43%，平均R約為0.93，模型在FSC值40%-70%時，性能稍差一些，主要是因為此時雪層較淺，積雪分布不穩(wěn)定。

此外，我們還進一步選擇了研究區(qū)上5天無云的影像（2011314、2013320、2014084、2016306、2017348）作為真值，再假設每景都被20%、50%、80%的云覆蓋，使用本研究提出的基于PCNN的方法去云后，將去云的影像與原始真值進行了比較，性能定量評價結果如表3。以2017328這一天為例，原始影像、假設的云覆蓋情況以及去云后的影像如圖4所示。結果表明：PCNN方法的魯棒性較好，即使云覆蓋率達到80%，基于PNCC的云像元重建方法也可以有效利用時空信息，將云像元的積雪覆蓋信息高質(zhì)量地重建。

圖3 基于云假設的性能定量評價結果。(a) MAE；(b) RFigure 3 Performance quantitative evaluation results based on cloud assumption.(a) MAE; (b) R

表3 5個晴空日上基于云假設的驗證結果Table 3 Performance evaluation results based on cloud assumption for five selected cloudless days

圖4 2017年12月14日（即2017328）的FSC分布圖。(a1)真值；(b1)-(b3)分別指20%、50%、80%的云覆蓋（白色區(qū)域）掩膜真值影像；(c1)-(c3) 分別指PCNN方法去云后的影像Figure 4 FSC distribution map on 14 December 2017 (i.e.2017328).(a1) Truth FSC image on December 14,2017; (b1)-(b3) indicate that the 20%, 50% and 80% of the original truth image are covered by artificial clouds (white area); (c1)~(c3) indicate the reconstructed FSC distribution map of (b1)-(b3) by PCNN method.

4 數(shù)據(jù)價值

本數(shù)據(jù)以 MODIS V6版積雪產(chǎn)品為基礎，對 MODIS NDSI產(chǎn)品中的云像元重建（包括：MODIS/Terra和MODIS/Aqua數(shù)據(jù)合成、臨近三天時間濾波和基于PCNN的云像元時空重建）的基礎上，進一步計算對應的FSC，制備了2000-2021年黃河源區(qū)積雪覆蓋比例數(shù)據(jù)集，可以為黃河源區(qū)積雪狀況評估、積雪演變特征分析、水資源綜合管理、雪災風險評估及預測等研究工作提供數(shù)據(jù)支撐。