青藏高原MODIS逐日無云積雪范圍產品精度驗證

李 諾，韓其飛，馬 英，黃曉東

（1.南京信息工程大學 地理科學學院，江蘇 南京 210041；2.蘭州大學 草地農業生態系統國家重點實驗室 草地農業科技學院，甘肅 蘭州 730000）

0 引言

積雪作為冰凍圈的重要組成成分之一，對氣候變化十分敏感。根據政府間氣候變化專門委員會（Intergovernmental Panel on Climate Change，IPCC）于2021 年8 月發布的第六次評估報告，冰凍圈正呈現加速萎縮狀態，北半球積雪范圍也在加速減?。?］。積雪因為其特定的高反照率和低導熱率，對地球能量平衡及水循環有著深刻的影響［2-3］。青藏高原地處中緯度地區，作為我國三大穩定積雪區之一，該地區的積雪影響東亞地區的大氣環流及天氣系統，其冰雪融水更是農業灌溉及周邊眾多河流的重要來源［4］。因此準確獲取青藏高原地區的積雪信息對研究青藏高原的徑流變化及氣候演變有著重要意義。

遙感技術相較于傳統的氣象臺站觀測數據具有大尺度、獲取信息方便快捷、不受地理條件限制等特征，成為大范圍積雪監測的關鍵手段［5］。搭載在Terra 和Aqua 衛星的中分辨率成像光譜儀（Moderate Resolution Imaging Spectroradiometer，MODIS）因其在時間、空間和光譜分辨率的較大優勢，在區域及半球尺度積雪監測中應用廣泛。MODIS積雪產品第5 版本（V5）提供了二值積雪范圍產品和積雪覆蓋率產品［6］，兩種產品的算法均基于歸一化差值積雪指數（Normalized Difference Snow Index，NDSI）。二值積雪產品通過設置NDSI 閾值對每個像元進行分類，即一個MODIS 像元如果NDSI 值大于等于0.4 且波段2 反射率大于11%，則該像元被分類為積雪［7］，并引入歸一化植被指數（Normalized Difference Vegetation Index，NDVI）用于提高MODIS 在林區的積雪分類精度［8］。評估結果表明，MODIS 逐日二值積雪產品（V5）在晴空狀態下的總體分類精度介于85%～99%［9］。積雪覆蓋率（Fractional Snow Cover，FSC）產品通過建立Landsat 獲取的FSC 與NDSI 之間的線性關系生成［10］。2016年發布的MODIS 第6版本積雪產品（V6）不再提供積雪二值產品與積雪覆蓋率產品，僅提供網格的NDSI數據。V6的產品結合高程重新定義了地表溫度對高海拔地區積雪的分類閾值。另外，采用定量圖像恢復算法重建了Aqua/MODIS 短波紅外波段（波段6），替代了V5 版本使用波段7 計算NDSI的策略［11］。V6版本較V5版本的積雪產品，質量和精度均得到了有效的提升［12］。

云是影響MODIS 逐日積雪遙感產品準確獲取地表積雪覆蓋范圍的主要因素之一。雖然MODIS逐日積雪產品在晴空狀態下具有可靠的精度，但由于云的影響，導致該產品在區域積雪監測中受到極大的限制。因此針對MODIS 逐日積雪產品的云下信息恢復，一直是領域內學者關注的一個熱點問題。國內外許多學者基于V5 產品針對云下信息恢復開展了大量的研究。唐志光等［13］利用三次樣條函數插值法對MODIS 逐日積雪產品MOD10A1 進行去云處理，得到青藏高原MODIS 逐日無云的積雪覆蓋率產品，該產品獲取的積雪日數（Snow-Covered Days，SCD）與地面觀測值得到的SCD 平均一致性為87.0%。侯小剛等［14］通過融合MOD10A1 與MYD10A1 數據，發展了一套適合新疆地區的逐日無云積雪范圍數據集，去云后產品與氣象臺站積雪觀測的一致性為88.1%，總體精度達到90.6%。黃曉東等［15］利用MODIS 上下午星獲取的逐日積雪產品和被動微波數據AMSR-E 雪水當量產品，對MODIS 逐日積雪圖像進行上下午星合成和鄰近日像元分析，并結合數字高程數據利用雪線算法進行云下信息恢復，獲取了青藏高原時空連續的MODIS 積雪逐日產品，總體分類精度達到90.7%。Gafurov等［16］在阿富汗東北部的Kokcha盆地，綜合了雙星合成、時間濾波、空間濾波等五種方法組成連續六步去云算法，生成的MODIS 逐日積雪產品，其總體分類精度也達到90%。Huang 等［17］提出一種基于隱馬爾可夫隨機場（Hidden Markov Random Field，HMRF）的時空濾波模型，充分利用了立方體內的時空信息對MODIS 逐日積雪產品進行去云處理，總體分類精度達88.0%，云量減少到1%。Hall等［18］提出一種基于MOD10C1積雪覆蓋產品的時間濾波去云方法，生成MODIS Cloud-Gap-Filled（CGF）逐日云填補積雪產品。隨著V6 版本積雪產品的發布，陸續有學者基于V6 積雪產品開展了類似研究。邱玉寶等［19］于2021 年11 月在科學數據銀行（Science Data Bank，ScienceDB）更新了基于V6 版本的青藏高原MODIS 逐日無云積雪面積數據集，該數據集共使用八個連續步驟進行去云。Muhammad 等［20］于2021 年3 月在Data Publisher for Earth & Environmental Science（PANGAEA）發布了基于MODIS V6產品的亞洲高山區逐日無云二值積雪數據集。郝曉華［21］于2020 年11 月在國家冰川凍土沙漠科學數據中心發布中國逐日無云500 m 積雪面積產品數據集（2000—2020 年），該數據集基于MODIS V6 逐日地表反射率產品重建NDSI，通過訓練不同土地覆蓋類型下NDSI 識別閾值，然后利用隱馬爾可夫隨機場時空濾波模型進行了去云處理。上述產品自發布后，在青藏高原還未得到系統驗證。

青藏高原作為我國三大積雪分布地區之一，海拔高且輻射強，導致積雪消融迅速、積雪破碎化嚴重，使得目前不同積雪產品在該地區的積雪監測精度均不夠理想。因此，從應用角度出發，對MODIS逐日無云積雪產品在青藏高原進行系統評估，研究其在積雪監測中的精度及影響因素尤為關鍵。鑒于此，本文利用中高分辨率的Landsat影像數據作為參考值，對上述三套MODIS 逐日無云積雪產品進行驗證，并分析土地覆蓋類型和積雪覆蓋率對MODIS積雪產品精度的影響，為評估青藏高原積雪范圍及其時空動態變化提供參考依據。

1 研究區及數據

1.1 研究區概況

青藏高原位于我國西南部，平均海拔4 000 m左右，被稱為“世界屋脊”“第三極”［22］（圖1）。青藏高原西起帕米爾高原向東至橫斷山脈，自南喜馬拉雅山脈南緣延伸至北昆侖山脈北緣和祁連山脈，總面積約2.57×106km2［23］。青藏高原以冰川、積雪、凍土等形態儲存了巨大的水資源，我國的黃河、長江、怒江、瀾滄江、雅魯藏布江等均發源于青藏高原，有“中華水塔”和“亞洲水塔”之稱［24］。青藏高原積雪水儲量關系著所在區域及周邊區域的生活和灌溉用水，影響高原及周邊區域的植被生長，該地區的積雪具有重要的水文、生態及氣候意義［25］。

圖1 青藏高原土地覆蓋類型及Landsat-8影像分布Fig.1 Land cover types in Tibetan Plateau and Landsat-8 images spatial distribution

1.2 MODIS逐日無云積雪產品

（1）中國逐日無云500 m 積雪面積產品數據集（MODIS CGF SCE）：來源于國家冰川凍土沙漠科學數據中心（http：//www.ncdc.ac.cn/）。該產品基于MODIS反射率產品MOD09GA/MYD09GA，利用高分辨率Landsat TM 數據作為真值，結合MODIS土地覆蓋分類產品，確定林區和非林區下積雪判別的指標閾值，利用MODIS 積雪反演算法獲取初級產品，通過隱馬爾科夫隨機場時空濾波算法進行云下信息恢復。該產品以HDF5 文件格式存儲，每個HDF5文件包含18個數據要素，其中包括數據值、數據起始日期、經緯度等。本文采用無云積雪數據，為二值影像。其編碼有0、1、2、3、4和255，對應的編碼含義為陸地、影像識別積雪、去云插補積雪、雪深插補積雪、水體和填充值［21］。

（2）青藏高原MODIS 逐日無云積雪面積數據集（MODIS_Dysno_Cloudfree）：來源于ScienceDB（https：//www.scidb.cn/）。該產品以MODIS V6 積雪產品MOD10A1 和MYD10A1 為基礎，通過上下午星合成、三天合成、“長時間”積雪和陸地法、鄰近四像元法、高程濾波法、修改陰影區錯誤分類和最大積雪陸地范圍掩膜這7 個連續步驟，獲得MODIS少云積雪產品，而后將研究區根據坡度劃分為“印度平原”“高原北部”“高原腹地”“藏東南山區”“帕米爾高原”“天山山脈”以及“喜馬拉雅山脈”7 個區域，針對7 個區域采用擬合預期雪線方法去除全部云污染，獲得MODIS 逐日無云積雪產品。本次所用產品為該數據集中青藏高原MODIS 逐日無云積雪產品面積數據集第二次更新版（2002—2021 年MODIS積雪產品C6.1版，http：//www.csdata.org/p/15/）。該產品為二值產品，影像分類代碼與MOD10A1 一致，其編碼有25、37、40、100、150 和200，對應含義為無積雪覆蓋的陸地、湖泊、湖泊不確定、湖冰、湖冰不確定和積雪。該產品利用地面臺站積雪深度觀測數據進行驗證，結果表明其積雪分類精度為78.4%，去掉雪深≤3 cm 的樣本后，精度達到89.0%［19］。

（3）亞洲高山區MODIS 逐日無云積雪數據集（M*D10A1GL06）：數據來源于PANGAEA（https：//www.pangaea.de/）。該產品以MODIS 積雪產品MOD10A1 和MYD10A1 為基礎，使用8d MOYDGL06*產品（雙星8 天融合產品）分別代替Terra 和Aqua 雪產品中的云像素。然后與全球陸地冰空間測量計劃（Global Land Ice Measurements from Space，GLIMS）發布的RGI（Randolph Glacier Inventory）6.0 冰川編目數據產品進行融合，最終達到99.99%的去云效果［20］。本產品為二值產品，其編碼含義如表1所示。

表1 M*D10A1GL06積雪產品屬性定義Table 1 Attribute definition of M*D10A1GL06 snow product

1.3 Landsat-8

前期研究大多基于站點雪深觀測數據對MODIS積雪產品進行精度驗證，但由于二者空間尺度的不匹配，導致基于點-面驗證的結果存在較大的不確定性。因此，本研究采用空間分辨率為30 m的Landsat-8 陸地成像儀（Operational Land Imager，OLI）獲取的影像作為驗證數據。OLI 影像首先要選擇晴空狀態的數據，以避免云雪混淆造成驗證數據集的不確定性；影像應分布在青藏高原積雪分布的主要區域，且能包含高原主要的土地覆蓋類型；另外，為了避免破碎化積雪對OLI積雪制圖的影響，驗證時段選擇積雪分布較為穩定的時期。因此，選取時間范圍為2015 年12 月2 日至2016 年1 月23 日共34 景影像作為驗證MODIS 逐日積雪產品的參考數據。

采用混合像元分析法首先獲取OLI像元尺度的積雪面積比例（Fractional Snow Cover，FSC）。其中端元提取的規則如表2所示。將30 m的OLI積雪面積比例數據聚合成500 m，定義升尺度后的亞像元積雪面積比例閾值0.5，生成二值積雪數據。其中喜馬拉雅山中東段的18 景Landsat 積雪數據來源于ScienceDB：“2013—2020 年喜馬拉雅山中段和東段Landsat 8 積雪覆蓋范圍數據”。該數據集采用支持向量機（Support Vector Machine，SVM）分類方法，選取不同地形、陰影等條件的積雪特征訓練樣本進行積雪分類，結合地表水體等輔助數據及空間鄰域分析進行分類后處理。通過對比Sentinel-2 高分辨率積雪分類數據，在900 m×900 m的網格內，其積雪覆蓋率相關系數在0.95 以上，均方根誤差約0.1%。該數據集為二值影像，其編碼有1、2、100、200 和300，對應的編碼含義分別為積雪、陸地、河流、冰湖和湖泊［27］。

表2 Landsat-8 OLI端元提取規則［26］Table 2 Endmember extraction rules for Landsat-8 OLI［26］

1.4 土地覆蓋數據

MODIS 土地覆蓋產品（MCD12Q1）用來比較不同土地覆蓋類型下的MODIS 逐日無云積雪產品的精度。采用國際地圈生物圈計劃（International Geosphere-Biosphere Program，IGBP）分類方案確定的17個土地覆蓋類別，其中包括11個自然植被類別，3個開發和鑲嵌土地類別，以及3 個非植被土地覆蓋類別［28］。為了避免類型過多造成評價結果的復雜性，本文根據將IGBP 分類方案重新分類為七大類：森林、灌叢、草原、耕地、城市建設用地，冰川和裸地（圖1）［29］。

2 精度驗證方法

使用高分辨率遙感產品驗證低分辨率產品的精度是遙感產品驗證的常用手段［30-31］。本研究利用高分辨率的Landsat積雪圖對MODIS兩個版本的數據進行對比驗證，采用精度評估指標包括總體分類精度（overall accuracy，OA）、低估誤差（underestimated error，UE）、高估誤差（overestimated error，OE），積雪分類精度P和Kappa系數。其中Kappa一致性檢驗是評價兩幅圖像之間一致性的方法［32］。當Kappa 系數為0～0.20 時，表明兩幅圖像的一致性極低；當Kappa 系數為0.20～0.40 時，表明兩幅圖像一致性一般；當Kappa 系數為0.40～0.60 時，表明兩幅圖像一致性中等；當Kappa 系數為0.60～0.80 時，表明兩幅圖像一致性極好；當Kappa 系數為0.80～1.00時，表明兩幅圖像高度一致［33］。

精度評價指標計算公式如下所示：

式中：a為Landsat 和MODIS 均記錄為積雪的樣本數；b為Landsat記錄有雪而MODIS被分類為非雪的樣本數，即漏判數；c為Landsat 和MODIS 均為非雪的樣本數；d為Landsat記錄為非雪而MODIS被分類為積雪的樣本數，即為誤判數；Pe表示偶然一致性，即在偶然機會下預測結果與真實結果一致的概率。

3 結果與討論

3.1 精度評價

驗證結果表明，MODIS CGF SCE、M*D10A1 GL06 和MODIS_Dysno_Cloudfree 產品的總體精度（OA）分別為88.2%、88.2%和89.6%，表明三種積雪產品有著相對較高的精度。其中MODIS CGF SCE低估誤差（UE）最小，僅為9.4%，M*D10A1GL06和MODIS_Dysno_Cloudfree 積雪產品UE 分別為14.7%和18.6%。三套產品的高估誤差（OE）分別為12.2%、8.9%和11.2%。由于MODIS CGF SCE 產品積雪分類精度（P）最高，達到90.6%，其他兩種產品分別僅為85.3%和81.4%。Kappa 系數統計結果表明三種積雪數據與Landsat積雪數據的一致性均表現為極好，分別為0.642、0.629 和0.652，MODIS CGF SCE 和MODIS_Dysno_Cloudfree 與Landsat 積雪數據一致性最為接近，M*D10A1GL06產品略差。

M*D10A1GL06 和MODIS_Dysno_Cloudfree 產品均是基于MODIS V6 版本生產的去云積雪產品，兩種產品的用戶精度理論上是一致的，造成積雪分類誤差不一致的原因主要是云重分類誤差導致的，去云后M*D10A1GL06 的積雪分類精度優于MODIS_Dysno_Cloudfree，說明后者采用的云下信息恢復算法存在較大的誤分類現象，與站點驗證結論也有一定的差距［20］。MODIS CGF SCE 產品是基于大氣校正的地表反射率數據生成的，并且分土地覆蓋類型設置最優NDSI 閾值，在積雪分類時利用雪深等輔助數據進行插值，因此積雪分類精度最高；而MODIS V6產品是經過輻射校正的反射率數據計算NDSI。驗證結果表明，基于大氣校正的地表反射率數據在積雪信息提取方面更具優勢，積雪分類精度優于另外兩種去云產品。

圖2 MODIS逐日無云積雪產品驗證結果Fig.2 Validation results for three kinds of MODIS daily cloud-free snow cover products

3.2 土地覆蓋對MODIS積雪產品精度的影響

表3中統計了MODIS逐日無云積雪產品在不同土地覆蓋類型下的積雪分類精度。其中草原、裸地和森林的驗證樣本（N）分別占比67%（N=2 391 609），20%（N=729 769）和10%（N=370 443），其余驗證樣本分布在其他土地覆蓋類型，因樣本占比較小，未統計精度評價結果。結果顯示，三種去云積雪產品的總體分類精度在不同土地覆蓋類型下均大于80%，其中MODIS CGF SCE 在林區的總體分類精度最優，達到89.6%，MODIS_Dysno_Cloudfree在草原和裸地的總體分類精度最優，分別達到91.6%和84.0%。三種產品在裸地的總體分類精度均相對較差，介于81%～84%之間。MODIS CGF SCE 產品在草原的漏判和誤判誤差相當（～10%），在裸地存在嚴重的高估誤差（20.9%），但在林區漏判誤差較大（21.5%）。M*D10A1GL06 產品在裸地高估誤差大，在林區漏判誤差大，分別達到20.0%和27.2%。MODIS_Dysno_Cloudfree 產品僅在林區存在較為嚴重的漏判誤差（21.3%）。MODIS CGF SCE 產品的積雪分類精度在草原和裸地類型下均優于另外兩種產品，林區的積雪識別精度與MODIS_Dys-no_Cloudfree 產品相當，分別為78.2%和78.7%，M*D10A1GL06 產品在林區的積雪識別精度最低，為72.8%。與OLI 積雪數據相比，三種產品在草原和裸地類型下，一致性均較高，林區的一致性較差，但MODIS CGF SCE 產品在林區與OLI 一致性最高，為0.249。雖然該產品針對不同土地覆蓋類型分別設置NDSI 閾值，并結合了林區積雪指數（Normalized Difference Forest Snow Index，NDFSI）用于提高林區積雪分類精度，但驗證結果顯示，相比基于MODIS 標準逐日積雪產品生成的去云產品，對林區積雪識別精度的改進并不明顯，依然存在較高的漏判誤差。

表3 土地覆蓋類型對MODIS積雪產品精度的影響Table 3 Assessment of MODIS cloud-free snow products under different land cover types

3.3 FSC對MODIS積雪產品精度的影響

為了評估積雪覆蓋率（FSC）對MODIS 逐日無云產品積雪分類精度的影響，本研究基于OLI 獲取的FSC，對像元尺度的FSC 進行分級，驗證不同FSC級別MODIS 逐日無云積雪產品的積雪分類精度。結果如圖3 所示，三種積雪分類精度隨著積雪覆蓋率的增大在逐漸升高。其中MODIS CGF SCE 逐日無云積雪產品積雪分類精度在NDSI 大于10（原始值/100 代表實際值）均大于其他兩種產品，Dysno_Cloudfree 產品積雪分類精度相對較差。進一步說明基于大氣校正的地表反射率數據計算NDSI，進而對積雪信息進行提取，相比MODIS 標準逐日積雪產品具有一定的優勢。

圖3 不同積雪覆蓋率下的MODIS積雪分類精度Fig.3 Accuracy of snow classification in graded fractional snow cover

4 結論

近年來，MODIS積雪產品在水文、水資源管理、氣象和氣候變化等眾多領域得到了廣泛的應用。在青藏高原地區，受限于該地區較強的太陽輻射，積雪消融迅速且破碎化嚴重，使得MODIS 積雪產品在該地區的準確性明顯低于高緯度穩定積雪區。本文選取了青藏高原基于MODIS 數據發展的三套逐日無云積雪產品，其中MODIS CGF SCE 是基于MODIS 地表反射率數據生產的，其他兩套產品（M*D10A1GL06 和MODIS_Dysno_Cloudfree）是基于美國雪冰數據中心（National Snow and Ice Data Center，NSIDC）發布的MODIS 逐日積雪產品進行去云處理后獲取的。本文利用分辨率較高的Landsat-8 OLI 積雪影像數據，對比驗證了三套逐日無云積雪產品的精度，并分析了土地覆蓋類型及積雪覆蓋率對產品精度的影響，得出以下結論：

（1）三套無云積雪產品的總體精度均優于80%，MODIS_Dysno_Cloudfree 積雪產品總體分類精度最高，達到89.6%，Kappa 系數為0.652。MODIS CGF SCE 低估誤差最小，為9.4%，積雪分類精度也是最高的，達到90.6%。M*D10A1GL06 高估誤差最低（8.9%）。

（2）土地覆蓋類型對MODIS 積雪產品積雪分類的精度影響較大，尤其在林區。三種積雪產品在林區的積雪分類精度均較差，存在較高的漏判誤差。雖然MODIS CGF SCE 產品針對林區進行了NDSI閾值優化，但效果甚微。

（3）隨著積雪覆蓋率的增加，三套無云積雪產品的積雪分類精度也隨之增加，其中MODIS CGF SCE 產品在不同NDSI 等級（NDSI＞10）均優于其他兩套產品。

綜上，基于MODIS 地表反射率產品生產的逐日無云產品（MODIS CGF SCE），針對不同土地覆蓋類型進行NDSI 閾值優化，在積雪識別精度方面提升較為明顯，但是林區積雪識別誤差較大的問題依舊沒有得到有效解決。當積雪覆蓋率較低時，積雪識別精度相比未經過大氣校正的MODIS 標準積雪產品，精度相對較好，但提升不顯著。因此，在積雪破碎化嚴重的青藏高原地區，除了考慮NDSI 閾值的優化，應進一步結合該地區的積雪特性，包括環境特征，進一步發展適合青藏高原的積雪制圖算法。