王曉艷， 陳思勇， 郭慧， 謝佩瑤， 王建， 郝曉華

1. 蘭州大學 資源環(huán)境學院, 蘭州 730000;

2. 南京大學 地理與海洋科學學院, 南京 210008;

3. 中國科學院西北生態(tài)環(huán)境資源研究院, 蘭州 730000;

4. 江蘇省地理信息資源開發(fā)與利用協(xié)同中心, 南京 210023

1 引 言

積雪是全球氣候變化的指示器，對氣候系統(tǒng)具有顯著的正反饋作用，是全球氣候變化研究的關鍵變量（秦大河 等，2007；Choi等，2010）。積雪在可見光和近紅外波段的高反照率特征會影響地表能量平衡（Robinson 等，1986；Huang 等，2017），積雪融化消耗能量會有效延緩溫室效應（Robinson和Kukla，1985）。與此同時，積雪融化導致地表反照率下降，使得地面吸收更多的熱量，促進溫度升高（楊梅學 等，2000）。在干旱和半干旱地區(qū)，季節(jié)性融雪是一種主要的水源，有超過世界六分之一的人口依賴冰川和積雪融水維持生命（Barnett等，2005；車濤和李新，2005；王建和李碩，2005；趙曉艷 等，2022）。積雪范圍的定量化研究對氣候變化、地表能量平衡和水文研究都具有重要意義（王建 等，2018）。

傳統(tǒng)的氣象站和積雪測線可以提供精確的地表積雪信息，但是由于站點的有限性并不能用于大范圍積雪監(jiān)測。衛(wèi)星遙感技術可以提供大范圍的積雪信息面上監(jiān)測，彌補了常規(guī)觀測在時間和空間上的不連續(xù)問題。中等分辨率成像儀MODIS（Moderate-resolution Imaging Spectroradiometer） 包含積雪識別所需的可見光和短波紅外波段，并且具有很好的時間和空間分辨率，在積雪檢測中得到廣泛的應用（Hall等，2002；Frei等，2012）。

大量研究表明，基于MODIS的積雪產品在晴空下具有較高的精度（Klein和Barnett，2003；Huang等，2011；Yang等，2015；Simic等，2004；Parajka和Bl?schl，2006）。然而，云遮擋導致MODIS產品大量的數(shù)據(jù)缺失并限制了MODIS 產品的應用，如何高精度地恢復云下像元、生成時空連續(xù)的高精度MODIS 積雪產品成為近年來研究的熱點問題。目前MODIS 積雪產品去云方法可以歸納為空間濾波方法、時間濾波方法、時空聯(lián)合濾波方法和多源數(shù)據(jù)融合方法4 大類（Li等，2019a）?？臻g濾波是當前常用的去云方法，該方法通常在空間鄰域上選擇無云像元來估計云覆蓋像元，對于大片的云遮擋區(qū)域不能獲得滿意的結果（Tran 等，2019；Gafurov等，2015）；采用Terra和Aqua一天合成以及多天合成的時間濾波方法，在連續(xù)多天的云覆蓋情況下無法恢復云下數(shù)據(jù)，同時在積雪快速變化的階段，多天合成會導致很大的誤差（Wang等，2014；Hall等，2010；Dozier等，2008）；時空聯(lián)合的方法綜合利用積雪的時空特征，在不同區(qū)域開展的實驗表明，時空聯(lián)合可有效提高積雪產品去云的精度（Li等，2017；Da Ronco和De Michele，2014；Parajka和Bl?schl，2008；邱玉寶 等，2017）；多源融合方法充分利用光學觀測、微波觀測和地面觀測等不同數(shù)據(jù)源的互補信息，也是一種有效的云去除方法。光學與微波數(shù)據(jù)相結合，可以用于較粗分辨率的積雪產品，但是產品的精度依賴于輸入的遙感產品精度（Li等，2019b；Huang等，2016；Yu等，2016；Deng等，2015；Gao等，2010；Liang等，2008；Chen等，2018）。近年來，出現(xiàn)了一些基于新技術的去云算法。Dong 和Menzel（2016）在實地雪深觀測的基礎上，利用條件概率插值對MODIS 積雪圖上的剩余云量進行重新分類，對德國西南部云覆蓋和積雪變化較高地區(qū)的積雪產品進行了去云處理；Huang等（2018）采用隱形馬爾科夫隨機場框架內集成MODIS光譜信息、時空上下文信息和環(huán)境關聯(lián)的技術，對2006年—2008年積雪季上里約熱內盧盆地的積雪產品去云；Hou等（2019）采用時空濾波結合機器學習的技術，生成了2001 年—2016 年中國北疆地區(qū)逐日無云MODIS 積雪面積比例產品。以上的MODIS 去云算法綜合MODIS 數(shù)據(jù)的時空特性并結合多源數(shù)據(jù)的優(yōu)勢特征，取得了較好的效果。

當前的積雪研究中常用的MODIS 積雪產品包括版本V005和V006，以上對積雪產品去云的研究主要針對V005 版本提供的二值積雪面積以及積雪面積比例產品。V005 版本二值積雪面積產品是基于歸一化差值積雪指數(shù)NDSI（Normalized Difference Snow Index）的特定閾值生成的，積雪面積比例產品也是通過與NDSI 建立線性回歸關系產生（Hall等，2002；Hall和Riggs，2007；Salomonson和Appel，2004）。然而，采用統(tǒng)一算法生成的MODIS積雪產品在不同地形、不同地表條件下的精度不盡相同，許多地區(qū)的積雪產品并不能達到令人滿意的精度（黃曉東 等，2007；曹云剛 等，2007；唐志光 等，2013；于小淇 等，2017；王曉艷 等，2017；Riggs等，2017）。NASA（National Aeronautics and Space Administration）于2016 年發(fā)布的MODIS 積雪產品V006 版本不再提供二值積雪面積和積雪面積比例產品，而是直接提供歸一化差值積雪指數(shù)NDSI 產品。與傳統(tǒng)的二值積雪產品和積雪面積比例產品相比，NDSI 產品可以保留更多的信息，從而更加靈活地應用于水文和生態(tài)過程模擬（Riggs 和Hall，2015）。然而，由于云的遮擋，當前的NDSI 數(shù)據(jù)集存在大量的數(shù)據(jù)缺失，影響了數(shù)據(jù)的使用。因此，亟需開展相關研究，生成高精度逐日無云NDSI數(shù)據(jù)集以更好地滿足積雪遙感監(jiān)測的需求。

Jing等（2019）提出了一種通過高斯核函數(shù)和誤差校正的融合框架生成無云MODIS NDSI 產品的方法，在青藏高原地區(qū)取得了較高的精度。該方法重構的NDSI 值為0 和10—100 的整數(shù)，并未考慮0—10的NDSI低值。然而在森林地區(qū)，由于冠層的遮擋，森林積雪通常具有較低的NDSI 值（Rittger等，2013；Hall 等，2007）。因而，恢復云下像元的NDSI原始值，包括低值的NDSI是非常必要的。

本文發(fā)展了一種基于鄰近相似像元的云去除方法。在NDSI 影像中，對于t0時刻某一個云覆蓋的目標像元，依次在前后日期無云的影像上尋找到與目標像元NDSI 值相近的像元作為鄰近相似像元，然后用t0時刻的這些鄰近相似像元進行加權計算得到目標像元的NDSI 值。本文以東北積雪區(qū)為實驗區(qū)，選用一個積雪季（2017-10-01—2018-04-30）的NDSI 影像進行去云實驗，并且將去云的NDSI 序列與氣象站點雪深序列進行對比用于確定二值積雪制圖中最優(yōu)的NDSI閾值。

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

2.1 研究區(qū)概況

東北是中國3大主要積雪區(qū)之一，包括黑龍江、吉林、遼寧和內蒙古東部地區(qū)，這里冬季漫長寒冷，積雪豐富，研究區(qū)包括森林、草原、耕地、居民地等多種地表類型（圖1）。東北平原被大興安嶺山脈、長白山脈和小興安嶺山脈所環(huán)繞，平均海拔為443 m。該區(qū)域的植被分布類型主要為溫帶針闊葉混交林、寒溫帶針葉林和溫帶草原3 類，植被覆蓋度平均在70%以上，其中森林覆蓋率達到40%（Che 等，2016）。選取該區(qū)域為研究區(qū)，獲取逐日無云NDSI 時間序列，結合氣象站點的雪深數(shù)據(jù)，可以揭示不同地表類型下NDSI 序列與積雪分布的變化規(guī)律，從而增強NDSI 在積雪識別中的應用。

圖1 研究區(qū)及氣象站點分布圖Fig. 1 Map of study area and the location of the meteorological stations

2.2 數(shù)據(jù)及預處理

2.2.1 MODIS積雪產品

為了減小積雪監(jiān)測的低估誤差和高估誤差，以便在全球尺度上獲得精確的積雪分布，MODIS積雪產品V006（MOD10A1，MYD10A1）不再提供V005 中的二值積雪面積和積雪面積比例產品，而是采用保守的積雪監(jiān)測方法，只提供MODIS NDSI和MODIS NDSI SNOW COVER 產品。其中，MODIS NDSI 的中像元的值在-10000—10000 之間，是將像元的原始NDSI 值擴大10000 倍的結果；MODIS NDSI SNOW COVER 中保留了可能是積雪的像元的NDSI 值（0—100），并且已經識別出云（250）、內陸水體（237）、海洋（239）等地物類型，同時還標識出了缺數(shù)據(jù)（200）、無精度（202）、夜間（211） 和探測器飽和（254） 等（Riggs 和Hall，2015）。MODIS NDSI SNOW COVER 中，NDSI 值在0—0.1 時均被置為0，10—100 為原始NDSI 值擴大100倍的結果（Riggs等，2016）。

東北地區(qū)森林覆蓋度高達40%，由于冠層的遮擋作用，森林積雪通常具有較低的NDSI 值（王曉艷 等，2017）。由于MODIS NDSI SNOW COVER將0—0.1 的NDSI 值均置為0，可能會丟失一些森林積雪的信息。因此我們選擇MODIS NDSI 序列進行去云處理，以獲取森林地區(qū)NDSI 的時序變化規(guī)律。本文選用了2017-10-01—2018-04-30 的MODIS Terra 和MODIS Aqua 的每日積雪產品（MOD10A1和 MYD10A1 V006），數(shù)據(jù)下載于美國國家雪冰數(shù)據(jù) 中 心（http：//nsidc.org/［2020-06-03］）。使 用MRT 工具將數(shù)據(jù)轉換為WGS84 投影。采用MODIS NDSI SNOW COVER 數(shù)據(jù)層的云類別對無云的MODIS NDSI 進行云掩膜，再將去云結果與MODIS NDSI實際值進行對比來評價該去云算法的精度。

2.2.2 氣象站點雪深數(shù)據(jù)

本文選取了研究區(qū)內均勻分布的24 個國家氣象局氣象站點，站點分布見圖1。以站點的雪深數(shù)據(jù)作為積雪分布的真值，與同時期的NDSI 進行對比分析。站點數(shù)據(jù)采用以cm 為單位的整數(shù)記錄了各個站點的積雪厚度，氣象站點的雪深數(shù)據(jù)存在個別異常值，在使用前進行了剔除。表1為各站點的位置信息及對應的地表類型。

表1 氣象站點位置及對應的地表類型Table1 Location and land cover type of the meteorological stations

2.2.3 MODIS土地覆蓋產品

為了確定不同地表覆蓋類型下NDSI 的最優(yōu)閾值，本研究采用2017 年空間分辨率為500 m 的土地覆蓋數(shù)據(jù)MCD12Q1。該產品是年度土地覆蓋分類產品，有5 種分類標準，本文選用Land_cover_type5 分類標準，在使用前將土地覆蓋類型合并林地、草地、耕地、水體和其他。該數(shù)據(jù)來源于美國航空航天局（https：//lpdaac.usgs.gov［2020-06-03］）。

3 MODIS NDSI產品去云方法

為了充分利用無云的MODIS NDSI 數(shù)據(jù)信息，在使用基于鄰近相似像元去云方法之前，首先使用MOD10A1/MYD10A1 結合和鄰近時間數(shù)據(jù)濾波去除部分云像元。MOD10A1/MYD10A1 結合和鄰近時間數(shù)據(jù)濾波的去云操作已經被廣泛采用（Paudel和Andersen，2011；Gao等，2010；Hall 等，2010），因此接下來的部分將重點闡述基于鄰近相似像元的去云方法。

3.1 基于鄰近相似像元去云算法

地表的地物分布具有一定的相關性，在空間上并不是相互獨立的，通常鄰近像元具有很大的光譜相似性，被稱為鄰近相似像元，這些鄰近相似像元具有兩個基本特征：（1）近似相等的光譜特征；（2）在多時相影像上具有相同的改變（Liu等，2019；Zhu 等，2010；Cheng 等，2017）。

如圖2 所示，假設M為t0時刻云遮擋的目標像元，該像元在無云的t0+d時刻記為M′，假定SM′為M′的鄰近相似像元，對應位置的SM為M的鄰近相似像元。那么，

圖2 鄰近相似像元示意圖Fig.2 Sketch map of similar pixels.

根據(jù)對鄰近相似像元的假設，從t0到t0+d，目標像元和鄰近相似像元的改變應該是一致的，也就是說Δ和Δ′相等。于是可以得到，

因此，選擇已知的M′、SM和SM′，即可預測出云覆蓋像元M 的NDSI 值。然而，由于一段時間內地表類型的改變，可能會使得某些鄰近相似像元的變化并不相同。選取當前窗口下n個鄰近相似像元，并對每個鄰近相似像元賦予權重來計算目標像元的值，可以有效避免由于單個鄰近相似像元的突變導致的較大預測誤差。這里，與目標像元的NDSI 差值最小的前n個像元被認為是鄰近相似像元。那么，目標像元的計算公式如下：

式中，xi，yi為第i個鄰近相似像元的x，y坐標值，x0，y0為目標像元的x，y坐標值，L為當前窗口的大小。式（6）將Di約束在1到1+ 2 范圍。

3.2 窗口大小與鄰近相似像元個數(shù)的選擇

為了確定式（4）中鄰近相似像元的個數(shù)n和式（6）中窗口的大小L，本文首先將無云的NDSI影像進行云掩模，然后采用不同L和n值進行云去除，采用本文3.3 節(jié)的去云算法預測云下的NDSI值，然后計算預測的NDSI值與NDSI實際值的相關系數(shù)，結果如圖3 所示。當L=15，n=20 時，相關系數(shù)r達到95.73%，繼續(xù)增加L和n，r并無明顯變化，綜合考慮預測精度和計算耗時，本文選擇的窗口大小為15×15，鄰近相似像元個數(shù)為20。

圖3 預測的NDSI與真值的相關系數(shù)變化曲面Fig. 3 The surface diagram of the correlation coefficient between predicted NDSI and true NDSI

3.3 去云流程

步 驟1，對MOD10A1 和MYD10A1 的NDSI 數(shù)據(jù)層進行合成得到MOYD。優(yōu)先采用MOD10A1 的有效數(shù)據(jù)，即如果MOD10A1 為無云像元，則取該值作為MOYD 的值；若MOD10A1 為云遮擋，MYD10A1 無云，則采用MYD10A1 為MOYD 的值；若均為云遮擋，則該像元仍然為云像元，等待下一步去云操作。早期的MYD10A1 使用B7 代替B6進行積雪判別具有較低的精度（Riggs 等，2016）。Gladkova 等（2012）發(fā)展的QIR（Quantitative Image Restoration）算法被用來恢復Aqua MODIS 的B6 波段用于積雪監(jiān)測，結果表明積雪監(jiān)測的精度與Terra MODIS 的MOD10A1相當。目前，QIR 算法已集成到V006 的MYD10A1 中（Riggs 等，2016）。因此，聯(lián)合V006 中MOD10A1 和MYD10A1 可以充分利用每天的無云信息。

步驟2，利用臨近時間的數(shù)據(jù)進行前后時間段的去云處理。本文選用了前后兩天的數(shù)據(jù)進行融合，結果記為ATC（Adjacent Temporal Composite）。如果選擇的時間間隔太長，可能會由于積雪的變化導致去云結果存在較大的誤差。如果僅采用前后一天的數(shù)據(jù)，則剩余的云量太多導致最后一步去云運算量過大。東北地區(qū)積雪相對穩(wěn)定，采用前后兩天數(shù)據(jù)融合是合理的。具體融合規(guī)則如下：對于當前的云遮擋像元，若前一天和后一天該位置像元均無云，則用這兩天的像元平均值來替代云遮擋像元；若只有一天為無云，則采用該無云像元的值替代云遮擋像元；若前一天和后一天該位置像元均被云遮擋，則進一步考慮前兩天和后兩天的影像，云像元的替換規(guī)則與前后各一天的替換規(guī)則一致。

1.2 家庭農場在現(xiàn)代農業(yè)生產的優(yōu)勢 家庭農場的內在市場化運作方式和家庭經營的穩(wěn)定性決定了在市場經濟下有較高的競爭力，在發(fā)達國家，以家庭農場為經營方式的已存在近200年，今天依然為農業(yè)的最重要經營方式。我國家庭農場是在家庭承包經營基礎上發(fā)展起來的，它保留了家庭承包經營的穩(wěn)定性特點，同時又吸納了現(xiàn)代農業(yè)市場化運作方式，其優(yōu)勢有以下幾點。

步驟3，采用基于鄰近相似像元的加權運算，對剩余的云遮擋像元進行去云處理，過程如下：

（1）對于t0時刻云遮擋的目標像元M，如果該目標像元在t0+d時刻是無云的，即M′已知，那么以該目標像元為中心，在15×15 窗口范圍內選取前20 個與目標像元具有最小NDSI 差異的像元SM′作為它的鄰近相似像元，同時還需滿足SM也是無云的，即M′、SM′和SM均為無云像元。

（2）依據(jù)式（6）計算每個鄰近相似像元到中心像元的距離Di，然后由式（5）計算出該鄰近相似像元的權重Wi，最后根據(jù)式（4）計算出云覆蓋的目標像元M的NDSI值。

（3）d 從1 開始，即先分別在前后一天的影像上尋找滿足條件的鄰近相似像元，如果前后一天都能找到滿足條件的鄰近相似像元，那么計算出的兩個M平均值作為最終的預測值；若只有一天滿足，則用這一天的M作為最終值；如果d=1時沒找到符合條件的20 個鄰近相似像元，則d=d+1 繼續(xù)以上操作，直到找到滿足條件的鄰近相似像元，計算出M。

（4）移動窗口，重復（1）—（3），對下一個有云像元進行去云處理，直至整個影像的云覆蓋率為0。

3.4 去云結果與精度驗證

本研究對東北地區(qū)2017-10-01—2018-04-30的MODIS NDSI 數(shù)據(jù)進行了去云處理，得到了無云的NDSI 序列。圖4 給出了部分影像逐步去云的結果，圖中CF為云覆蓋比例。

圖4 最后一列從上到下分別為2018-01-01、2018-02-01 和2018-03-01 云去除后的NDSI 影像，影像上NDSI 的分布符合東北積雪分布規(guī)律。在這些日期，大興安嶺、小興安嶺以及長白山地區(qū)都被積雪覆蓋，而森林地區(qū)積雪NDSI 較低，在0—0.3 之間；高緯度地區(qū)的草地耕地居民地積雪具有較高的NDSI 值，在0.5 以上；西南部草原耕地少雪，NDSI為負值。

圖4 MODIS NDSI去云影像結果Fig. 4 Cloud removal results for the MODIS NDSI images

為了進一步定量驗證去云的精度，本文采用“云假設”的方法?！霸萍僭O”的思想是在無云的影像上添加云掩膜，然后對添加云掩膜的影像進行去云處理，再將去云結果與真實的影像進行對比來評價算法的好壞。為了使得添加的云接近真實的云分布，本文選取了每個月中云量最小的兩景影像作為NDSI 影像的真值，將本月中云量最大的云覆蓋到該NDSI 影像上，這里采用的是MODIS NDSI SNOW COVER中的云類別（值為250的像元）對MODIS NDSI進行云掩膜。

MOD10A1和MYD10A1的合成以及ATC 已被廣泛采用并證明具有較高的精度（Paudel和Andersen，2011；Dong和Menzel，2016）。這里只討論基于鄰近相似像元的去云結果，表2中的樣本比例即為采用基于鄰近相似像元進行去云的像元占總像元的比例，時間窗口是指完全去除云遮擋需要的天數(shù)。使用相關系數(shù)r、均方根誤差RMSE和平均偏差MAE這幾個指標（式（7）—（9））對去云結果進行精度評價。

表2 基于“云假設”的去云結果評價Table 2 Cloud removal result assessment based on “cloud assumption”

表2 給出了10 景實驗影像在基于鄰近相似像元去云過程中云像元的比例、時間窗口的大小以及r、RMSE、和MAE 的值。較高的r值和較低的RMSE 以及較低的MAE 代表反演的結果更加接近真值。10景實驗影像去云結果中，相關系數(shù)r的平均值為0.95，均方根誤差RMSE 的均值為0.08，平均絕對誤差MAE 的均值為0.05。時間窗口從22 天到61 天不等，時間窗口的大小與云覆蓋比例并無顯著相關，對去云結果的精度也無明顯影響。

4 NDSI 的時序特性及其在積雪識別中的最優(yōu)閾值

4.1 NDSI時序變化

圖5 為一些典型氣象站點所在像元的去云NDSI序列與站點雪深序列對應圖（站點編號見圖1和表1）。站點所在的像元包括耕地、居民地、林地和混合像元，序列的日期從2017 年10 月1 日—2018 年4 月30 日共212 天，圖中橫軸day1 代表2017年10月1日。

圖5 典型地表的NDSI序列與雪深變化曲線（橫軸的序列日期從2017-10-01—2018-04-30； 第一天代表2017-10-01）Fig.5 The curve of NDSI and snow depth at typical surface （The date on the x-axis is from October 1， 2017 to April 30， 2018， and the first day is October 1， 2017）

圖5（a）為氣象站點1 所在像元的NDSI 及雪深序列，該像元為純凈的居民地像元。穩(wěn)定積雪時期該像元的NDSI 值在0.6 附近，雪深在5 cm 以上。2017 年11 月10 日氣象站點觀測到雪深由0 變?yōu)? cm，對應的NDSI 由負值上升為正值，之后一直在0.6 上下浮動，直到2018 年3 月16 日融雪開始。融雪過程自2018 年3 月16 日持續(xù)到2018 年3 月28 日，期間氣象站點所測雪深由15 cm 逐漸減小，NDSI 值也逐漸降低。當雪深降為0 時，即積雪完全融化，此時NDSI 變?yōu)樨撝?，融雪過程持續(xù)了12 d。

圖5（b）對應的氣象站點2所在像元為純耕地像元，該像元的NDSI 時序變化與氣象站點所測的雪深序列也有很好的對應關系。氣象站點觀測到雪深由0 變?yōu)?5 cm 的同時，NDSI 由負值驟升到0.6。之后漫長的冬季，NDSI 一直保持在0.6 左右，站點觀測的雪深在10—20 cm。2018 年3 月18 日開始，氣象站點所測雪深逐漸降低，NDSI 值也隨之降低。至2018 年3 月26 日積雪完全融化，雪深降為0，NDSI 同時降為負值，積雪融化過程大約持續(xù)8天。

圖5（c）對應的氣象站點4所在像元為耕地居民地的混合像元，該像元的NDSI 變化與雪深變化也有很好的一致性。在穩(wěn)定的積雪季，像元NDSI值在0.6上下波動，雪深在5—10 cm。從2018年3月26日開始至2018年3月29日，雪深由5 cm 降為0，NDSI值也由0.6降為負值。因雪深較淺，融雪過程僅持續(xù)3天。

圖5（d）對應氣象站點19，該站點所在的像元為大興安嶺的落葉林。林地NDSI 的變化與站點所測雪深具有很好的對應關系。2017年10月初有一次5 cm 降雪，因為積雪快速融化，NDSI 值有響應但并未能升至大于0。接下來10月下旬和11月初的降雪NDSI值都有相應的變化響應。2017年11月10日至2018年3月18日為穩(wěn)定的積雪期，雪深在15 cm左右，此時NDSI 的值在0.2—0.4 之間，低于積雪期耕地和居民地的NDSI值。從2018年3月9日開始NDSI值持續(xù)下降，至2018年3月25日NDSI降為負值，雪深由10 cm降為0，該融雪過程持續(xù)16天。

圖5（e）對應的氣象站點22 位于小興安嶺的針闊葉混交林。自2017年11月10日站點觀測到的雪深由0變?yōu)? cm，接下來的積雪季，雪深逐步增加，2018年3月份雪深增至25 cm，而NDSI的值始終在0—0.2之間。融雪開始于2018年3月28日，至2018 年3 月31 日雪深變?yōu)?，同時NDSI 降至小于0。因該像元融雪開始時間較晚，較高的氣溫使得融雪過程僅持續(xù)了3 天時間。

圖5（f）為緯度較低的氣象站點15所在像元，從該站點的積雪觀測數(shù)據(jù)來看，在2018年3月12日有一場明顯的降雪，NDSI 的時序變化圖也準確地反映了此次降雪過程。

圖5（a）—5（e）站點均位于高緯度地區(qū)，冬季漫長寒冷，有明顯的積雪穩(wěn)定期。NDSI 與氣象站點所測雪深有很好的對應關系。在氣象站點觀測的雪深由0 變?yōu)榉橇銜r，NDSI 會有相應的躍升。隨著積雪融化，觀測的雪深逐漸減小為0，相應的NDSI 值也會變?yōu)樨撝?。圖5（f）表明，NDSI 也能準確反映短暫的降雪過程。

4.2 積雪提取中NDSI的最優(yōu)閾值選擇

圖6 表明，NDSI 時序變化與觀測的雪深具有很好的對應關系，因此NDSI 是進行積雪識別的有效指數(shù)。接下來討論NDSI 的閾值與積雪識別精度的關系。

圖6 積雪識別精度隨NDSI閾值的變化曲線Fig. 6 The curve of snow recognition accuracy with NDSI threshold

通常，站點雪深數(shù)據(jù)大于等于1 cm 時，認為該站點所在的像元為有雪像元；站點雪深為0，則認為該站點所在的像元為無雪像元（Parajka和Bl?schl，2006；Ault 等，2006）。那么根據(jù)站點雪深數(shù)據(jù)，可以將站點所對應的像元分為積雪像元和無雪像元，作為積雪分布的真值。選取一定的NDSI 閾值，同樣可以將像元分為積雪像元（NDSI>給定閾值）和無雪像元（NDSI<給定閾值），得到二值積雪分布。以雪深的積雪分類作為真值，對采用不同NDSI閾值制作的二值積雪分布圖進行精度評價，評價指標包括總體精度OA （Overall Accuracy）、低 估誤差UE （Underestimation Error）和OE 高估誤 差（Overestimation Error）（Parajka 和Bl?schl，2008），見式（10）—（12），其中，a、b、c、d含義如表3 所示。因為森林和非森林地區(qū)的NDSI值有較大差異，這里分別對森林和非森林地區(qū)進行評價。

表3 基于雪深和NDSI的二值積雪識別的混淆矩陣Table 3 Confusion matrix comparing the snow depth（SD） and NDSI binary snow cover

（1）非森林像元。采用了圖1中的非森林站點1—17，雪深大于等于1cm 則記為snow。隨著NDSI閾值的變化，總體精度的變化如圖6（a）：當NDSI=0.1時總體精度最高，達到95.6%，此時，高估誤差和低估誤差分別為1.3%和4.7%。

（2）森林的像元。通常站點安置在森林附近的道路旁或者居民區(qū)，站點所在的像元為混合像元。為了統(tǒng)計純森林像元的NDSI 變化與雪深的關系，參照Google Earth 選取站點附近（<2 km）的純森林像元，用站點的雪深數(shù)據(jù)代表周圍森林像元的雪深。本文選取了站點18—24 周圍的18 個純森林像元，對比森林像元的NDSI 序列與站點雪深的變化規(guī)律，統(tǒng)計不同NDSI 閾值下積雪識別的精度。在森林地區(qū)，隨著NDSI 閾值的變化，森林像元的積雪識別的精度變化如圖6（b）：當NDSI=0時，精度最高可達到93.5%，低估誤差和高估誤差均為3.8%。選取的NDSI 閾值大于0.1，低估誤差會急劇增加，總體精度也相應快速降低。

5 結 論

本文提出了一種基于鄰近相似像元的去云方法，對MODIS 積雪產品NDSI V006 進行去云處理，得到逐日無云NDSI 序列。選取中國東北地區(qū)2017-10-01—2018-04-30 一 個 積 雪 季 的MODIS NDSI 產品進行去云實驗并采用“云假設”方法進行精度評價。結果表明，本文提出的基于鄰近相似像元的去云方法可以很好的恢復云下像元的NDSI 值，預測值與真值的相關系數(shù)達到0.95，是一種有效的去云方法。

無論在林地還是非林地、高緯度還是低緯度地區(qū)，無云的NDSI 序列與氣象站點觀測到的雪深序列都具有很好的一致性，NDSI 可以準確地反映站點所在像元的積雪的變化。無雪時的林地和非林地像元的NDSI 基本為負值，在有雪的情況下，林地像元和非林地像元的NDSI 存在差異。非林地積雪像元具有較高的NDSI值，約為0.6。由于森林冠層的影響，林地積雪的NDSI 值較低，大興安嶺落葉林具有較高的透過率，林地積雪的NDSI 值在0.3 左右，而小興安嶺的針闊葉混交林透過率較低，林地積雪的NDSI值在0—0.2之間。

將雪深大于等于1 cm 記為有雪，作為地表積雪分布的真值，統(tǒng)計不同NDSI 閾值下積雪識別的精度。結果表明，在森林和非森林地區(qū)的NDSI 的閾值分別取0 和0.1 時，可以得到最高的積雪識別精度95.6%和93.5%。Zhang 等（2019）基于中國境內的200多個站點的日雪深數(shù)據(jù)，推薦中國境內對MODIS V6 產品使用NDSI 閾值0.1 來判定積雪。由于采用的森林站點數(shù)量有限，Zhang 等得出的NDSI=0.1 的閾值主要是指在非森林地區(qū)，本研究的結論與此一致。

Hall 等（2002）給出的MODIS 全球積雪面積產品NDSI閾值為0.4；郝曉華等（2008）選取祁連山中部山區(qū)的3 個積雪子區(qū)進行積雪識別并采用Landsat 數(shù)據(jù)進行精度評價，結果表明，NDSI 平均閾值為0.33 時積雪識別的總精度達到最高。本文給出的NDSI 閾值與以上研究的結果有較大差異，主要是因為去云之后的NDSI 產品，在積雪提取中不需要再考慮云雪混淆問題，只需區(qū)分有雪和無雪的情形，故NDSI 閾值適當降低，減少低估誤差的同時不會引入將云錯分為雪的高估誤差。

圖6 表明，非森林像元的NDSI 閾值從-0.1 到0.4，積雪識別的精度均在90%以上，NDSI=0 時的積雪識別精度為93.2%，略低于NDSI=0.1 時的95.6%。而森林像元的NDSI 閾值大于0.1 之后，低估誤差急劇增加，積雪識別的總精度明顯下降。因此，如果在積雪識別過程中存在地表分類圖，則可以在森林和非森林地區(qū)分布采用NDSI 閾值為0 和0.1，若沒有相應的地表分類圖，則可以統(tǒng)一使用NDSI=0作為積雪識別的閾值。

本文提出的MODIS NDSI 產品去云方法，可以高精度地恢復云下像元，獲取逐日無云的NDSI 產品，使其更好地應用于積雪識別中?；谡军c雪深選取的東北地區(qū)森林和非森林地區(qū)NDSI 最優(yōu)閾值，為NDSI 在積雪識別中的應用提供了有力的支撐。