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

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

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

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

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

引 言

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

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

為提高 MODIS積雪產(chǎn)品的監(jiān)測精度，大量學(xué)者對 MODIS積雪產(chǎn)品中的云覆蓋進(jì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)品，因?yàn)槔肗DSI值可以很方便地推算出二值積雪覆蓋和FSC值[13]。目前，已有研究基于相似像元/塊的思想，對MODIS NDSI積雪產(chǎn)品中的云覆蓋進(jìn)行了重建[28-30]，但受青藏高原云系旺盛、積雪斑塊化嚴(yán)重以及積雪消融較快等因素影響，在該區(qū)域搜索到合適的相似像元/塊并非易事，有時甚至根本就不存在合理的相似像元。XING等[31]利用深度學(xué)習(xí)技術(shù)深入挖掘時空特征，充分利用MODIS NDSI積雪產(chǎn)品中可用的時空信息，構(gòu)建了基于部分卷積神經(jīng)網(wǎng)絡(luò)（Partial Convolution Neural Network，PCNN）的MODIS NDSI積雪產(chǎn)品云像元重建模型，并在黃河源區(qū)進(jìn)行了驗(yàn)證，結(jié)果表明PCNN模型重建后的MODIS NDSI產(chǎn)品達(dá)到了可靠的較高精度。

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

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

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

黃河源區(qū)MODIS FSC數(shù)據(jù)集的制作流程主要分為5個步驟：獲取MODIS積雪產(chǎn)品、數(shù)據(jù)預(yù)處理、MODIS NDSI云像元重建（包括：MODIS/Terra和MODIS/Aqua數(shù)據(jù)合成、臨近三天時間濾波和基于PCNN的云像元時空重建）、FSC估計和精度驗(yàn)證，具體流程如圖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/）免費(fèi)下載了研究區(qū)2000-2021年的數(shù)據(jù)。覆蓋整個黃河源區(qū)需要二個瓦片（即：h25v05和h26v06），影像數(shù)據(jù)均為HDF文件存儲格式、Sinusoidal投影。

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

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

利用MODIS MRT批處理工具對兩個瓦片進(jìn)行拼接、重投影、坐標(biāo)轉(zhuǎn)換、重采樣等相關(guān)預(yù)處理，將原始數(shù)據(jù)轉(zhuǎn)換為地理坐標(biāo)系，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ù)進(jìn)行云掩膜，并將Raw NDSI實(shí)際值小于0的像元的NDSI值設(shè)置為0，將處理后的Raw NDSI數(shù)據(jù)作為基礎(chǔ)數(shù)據(jù)層。

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

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

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

1.4.2 鄰近時間濾波

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

1.4.3 基于PCNN的云像元重建

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

1.5 MODIS FSC估計

使用NASA第五版MODIS FSC的標(biāo)準(zhǔn)算法[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值，得到對應(yīng)的FSC產(chǎn)品，分別簡稱為MOD、MYD、MOYD、ATF和PCNN產(chǎn)品。

1.6 精度驗(yàn)證

1.6.1 基于站點(diǎn)雪深觀測值的驗(yàn)證

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

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

3個驗(yàn)證指標(biāo)總體精度（OA）、高估積雪事件（MO）和低估積雪事件（MU）分別定義如下：

1.6.2 基于云假設(shè)的驗(yàn)證

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

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

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

黃河源區(qū)2000-2021年無云積雪覆蓋比例數(shù)據(jù)集包含從2000-2001積雪季開始（即2000年11月1日）至2020-2021積雪季結(jié)束（即2021年4月30日）的7485個逐日的柵格數(shù)據(jù)文件。為便于數(shù)據(jù)在通用遙感軟件下處理與展示，將數(shù)據(jù)存儲為tif格式，投影方式為WGS84地理坐標(biāo)系，文件命名規(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ū)上的實(shí)際情況相符合。積雪季上不同時間的 FSC空間分布很好地展示出了積雪的積累-消融變化過程，一些積雪事件（如2018年11月15日以及2019年2月底的新降雪），MODIS FSC產(chǎn)品都有很好的及時響應(yīng)。

圖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 基于站點(diǎn)雪深觀測值的評估

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

表2 基于站點(diǎn)觀測雪深值的驗(yàn)證結(jié)果Table 2 Performance evaluation results based on snow depth observations from stations

3.2 基于云假設(shè)的評估

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

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

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

表3 5個晴空日上基于云假設(shè)的驗(yàn)證結(jié)果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)品為基礎(chǔ)，對 MODIS NDSI產(chǎn)品中的云像元重建（包括：MODIS/Terra和MODIS/Aqua數(shù)據(jù)合成、臨近三天時間濾波和基于PCNN的云像元時空重建）的基礎(chǔ)上，進(jìn)一步計算對應(yīng)的FSC，制備了2000-2021年黃河源區(qū)積雪覆蓋比例數(shù)據(jù)集，可以為黃河源區(qū)積雪狀況評估、積雪演變特征分析、水資源綜合管理、雪災(zāi)風(fēng)險評估及預(yù)測等研究工作提供數(shù)據(jù)支撐。