MODIS NDSI產(chǎn)品去云算法及最優(yōu)閾值選擇研究

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

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

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

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

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

1 引 言

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

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

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

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

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

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

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

2.1 研究區(qū)概況

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

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

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

2.2.1 MODIS積雪產(chǎn)品

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

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

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

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

表1 氣象站點(diǎn)位置及對(duì)應(yīng)的地表類型Table1 Location and land cover type of the meteorological stations

2.2.3 MODIS土地覆蓋產(chǎn)品

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

3 MODIS NDSI產(chǎn)品去云方法

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

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

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

如圖2 所示，假設(shè)M為t0時(shí)刻云遮擋的目標(biāo)像元，該像元在無云的t0+d時(shí)刻記為M′，假定SM′為M′的鄰近相似像元，對(duì)應(yīng)位置的SM為M的鄰近相似像元。那么，

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

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

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

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

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

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

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

3.3 去云流程

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

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

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

步驟3，采用基于鄰近相似像元的加權(quán)運(yùn)算，對(duì)剩余的云遮擋像元進(jìn)行去云處理，過程如下：

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

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

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

（4）移動(dòng)窗口，重復(fù)（1）—（3），對(duì)下一個(gè)有云像元進(jìn)行去云處理，直至整個(gè)影像的云覆蓋率為0。

3.4 去云結(jié)果與精度驗(yàn)證

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

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

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

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

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

表2 基于“云假設(shè)”的去云結(jié)果評(píng)價(jià)Table 2 Cloud removal result assessment based on “cloud assumption”

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

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

4.1 NDSI時(shí)序變化

圖5 為一些典型氣象站點(diǎn)所在像元的去云NDSI序列與站點(diǎn)雪深序列對(duì)應(yīng)圖（站點(diǎn)編號(hào)見圖1和表1）。站點(diǎn)所在的像元包括耕地、居民地、林地和混合像元，序列的日期從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）為氣象站點(diǎn)1 所在像元的NDSI 及雪深序列，該像元為純凈的居民地像元。穩(wěn)定積雪時(shí)期該像元的NDSI 值在0.6 附近，雪深在5 cm 以上。2017 年11 月10 日氣象站點(diǎn)觀測(cè)到雪深由0 變?yōu)? cm，對(duì)應(yīng)的NDSI 由負(fù)值上升為正值，之后一直在0.6 上下浮動(dòng)，直到2018 年3 月16 日融雪開始。融雪過程自2018 年3 月16 日持續(xù)到2018 年3 月28 日，期間氣象站點(diǎn)所測(cè)雪深由15 cm 逐漸減小，NDSI 值也逐漸降低。當(dāng)雪深降為0 時(shí)，即積雪完全融化，此時(shí)NDSI 變?yōu)樨?fù)值，融雪過程持續(xù)了12 d。

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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

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