韓超，沈言龍，歐陽志棋，謝佩瑤，郭慧，陳思勇，王曉艷*

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

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

引 言

積雪是地表重要的覆蓋物，也是重要的淡水資源[1]，其較高的反照率和較低的導(dǎo)熱率特性，對(duì)地表輻射平衡起到了重要的影響，同時(shí)積雪時(shí)空變化對(duì)全球生態(tài)系統(tǒng)及區(qū)域氣候的變化起到至關(guān)重要的作用[2-7]。中國積雪主要分布于新疆，青藏高原和東北地區(qū)[2]。其中東北地區(qū)是中國重要的糧食生產(chǎn)基地。積雪對(duì)春季土壤墑情和作物生長的影響十分重要，會(huì)直接影響作物實(shí)際產(chǎn)量[8]。因而，研究東北地區(qū)積雪分布及其時(shí)空變化特性具有重要意義。

隨著遙感和地理信息技術(shù)的發(fā)展，人們逐步開展了一系列有關(guān)積雪產(chǎn)品的研發(fā)。MODIS積雪產(chǎn)品具有較高的時(shí)間分辨率，其空間分辨率為500 m，近年來被廣泛應(yīng)用于積雪時(shí)空變化分析研究。其中，歸一化差值積雪指數(shù)（Normalized Difference Snow Index，NDSI）在MODIS積雪面積制圖中被廣泛應(yīng)用[9-12]。我國東北地區(qū)地處高緯度，地形地貌相對(duì)復(fù)雜，森林分布廣泛[13-15]。受冬季太陽光照條件以及林地冠層的影響，當(dāng)前 MODIS積雪產(chǎn)品 V6版本在該地區(qū)存在過度云掩膜問題，MOD/MYD10A1的NDSI_Snow_Cover數(shù)據(jù)層存在大量的數(shù)據(jù)空缺，尤其在大興安嶺及小興安嶺森林地區(qū)該問題較為嚴(yán)重。本文對(duì)位于大興安嶺地區(qū)的一個(gè)站點(diǎn)（編號(hào)為50247，經(jīng)緯度為123.5579°E，52.0285°N）進(jìn)行了云覆蓋累計(jì)天數(shù)統(tǒng)計(jì)，發(fā)現(xiàn)該站點(diǎn)所處位置在MODIS積雪產(chǎn)品中月內(nèi)云覆蓋累計(jì)天數(shù)在1、2、3、11和12月份都高達(dá)26天以上。這導(dǎo)致該地區(qū)整個(gè)積雪季可用的NDSI數(shù)據(jù)極少，很大原因是MODIS積雪產(chǎn)品中將積雪林地錯(cuò)分為了云，從而限制了MODIS積雪產(chǎn)品在中國東北地區(qū)的應(yīng)用[16]。當(dāng)前的 MODIS積雪產(chǎn)品去云算法主要采用時(shí)空鄰域的像元信息進(jìn)行插值填補(bǔ)云下像元，這些方法對(duì)于大范圍連續(xù)的云覆蓋去云效果并不理想[17-23]。CHEN等[24]提出的時(shí)空自適應(yīng)的方法可以生成時(shí)空連續(xù)的 NDSI序列，但是當(dāng)存在過度云掩膜的情況下，該算法的效率和精度都會(huì)大大降低。另外，由于MODIS積雪產(chǎn)品中采用的low_NDSI screen，使得許多具有較低NDSI的森林積雪像元被判為無雪像元，造成森林積雪的漏分。

本文以MODIS V6版本MOD/MYD10A1和MOD09GA數(shù)據(jù)集為數(shù)據(jù)源，采用上下午星合成、決策樹判斷和時(shí)空自適應(yīng)去云算法，顯著提高了產(chǎn)品精度和計(jì)算效率，生成了2000-2020年中國東北地區(qū)積雪季逐日無云 NDSI序列產(chǎn)品。本數(shù)據(jù)可用于支持中國東北地區(qū)積雪分布時(shí)空變化、積雪物候等后續(xù)科研工作。

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

1.1 數(shù)據(jù)源

在MODIS V6積雪產(chǎn)品中，MOD10A1/MYD10A1產(chǎn)品包括Raw_NDSI、NDSI_snow_cover等數(shù)據(jù)集。兩個(gè)數(shù)據(jù)層分辨率一致，像元位置一一對(duì)應(yīng)，其中，Raw_NDSI數(shù)據(jù)集的數(shù)值范圍為-10000～10000，該值是每個(gè)像元真實(shí)NDSI值的10000倍。NDSI_snow_cover數(shù)據(jù)集在Raw_NDSI數(shù)據(jù)集的基礎(chǔ)上進(jìn)行了屬性分類，屬性編碼內(nèi)容如下表所示[25]。NDSI_snow_cover數(shù)據(jù)集中10-100的數(shù)值代表可能的積雪像元，該值是像元真實(shí)NDSI值的100倍。MOD09GA產(chǎn)品包括地面7個(gè)波段的地表反射率數(shù)據(jù)，本研究選取了MOD/MYD10A1中的Raw_NDSI和NDSI_snow_cover以及MOD09GA中的 sur_refl_bo4(band4)作為數(shù)據(jù)源（表1）。上述產(chǎn)品數(shù)據(jù)來源于 Google Earth Engine（https://code.eartengine.google.com）?？臻g分辨率為500 m，時(shí)間分辨率為1 d，數(shù)據(jù)格式為*.tiff。

表1 MODIS積雪產(chǎn)品V6版屬性定義Table1 Attribute definitions of MODIS Snow Product V6

1.2 數(shù)據(jù)處理方法

1.2.1 影像裁剪

在Google Earth Engine上首先對(duì)研究區(qū)積雪產(chǎn)品數(shù)據(jù)集進(jìn)行裁剪。由于所需數(shù)據(jù)集的空間分辨率都是500 m，故不用重采樣。數(shù)據(jù)統(tǒng)一處理后的空間范圍為北緯37°-54°30′，東經(jīng)114°-136°之間，范圍包含中國東北積雪區(qū)（如圖1），地理坐標(biāo)系為WGS84。

圖1 b2、b4和b6波段合成結(jié)果（2019-02-12）Figure 1 Synthetic results of Band2, Band4 and Band6（Feb.12, 2019）

1.2.2 影像合成

結(jié)合MOD/MYD10A1中的Raw_NDSI和NDSI_snow_cover兩個(gè)數(shù)據(jù)集，進(jìn)行上下午星合成，填充一些云像元，恢復(fù)其NDSI值[26]。算法規(guī)則為若原產(chǎn)品上午定義有云，下午定義無云，則將下午星NDSI值賦值給當(dāng)天；若上下午星都有云，則未被恢復(fù)，依然保留該像元為云像元；其余情況則采用上午星的NDSI值。經(jīng)過此步驟處理，可以一定程度地減小云覆蓋比例。以2018年為例，平均每景影像中云量占比由45%降至36%。

1.2.3 決策樹判別

采用決策樹判別方法，恢復(fù)MOD10A1中錯(cuò)分為云的像元NDSI值。此步驟只針對(duì)MOD10A1中的云像元進(jìn)行操作，決策樹流程見圖2。研究發(fā)現(xiàn)，在有雪季節(jié)，東北地區(qū)的云雪混淆現(xiàn)象非常明顯，主要表現(xiàn)為將積雪或者森林積雪錯(cuò)分為云。故此步驟針對(duì)每年11月至次年2月的積雪季數(shù)據(jù)實(shí)施。NSIDC（National Snow and Ice Data Center，美國國家冰雪數(shù)據(jù)中心）發(fā)布的MODIS全球積雪面積產(chǎn)品采用的NDSI閾值為0.4，可以有效進(jìn)行云雪的區(qū)分[27-28]。本文對(duì)東北地區(qū)大量的云樣本進(jìn)行了統(tǒng)計(jì)，結(jié)果表明云像元NDSI最大值不超過0.4。因此，首先對(duì)Raw_NDSI進(jìn)行均值濾波，算法規(guī)則是像元周圍9個(gè)像元的值直接平均賦值給中心像元，得到新的數(shù)據(jù)層Mean_NDSI。若MOD10A1中的云像元對(duì)應(yīng)的Mean_NDSI大于0.4，則認(rèn)為該像元不是云像元，恢復(fù)該像元的Raw_NDSI值。

圖2 決策樹流程圖Figure 2 The flow chart of decision tree

當(dāng) NDSI在0-0.4之間的像元，主要包含了森林積雪、積雪邊界和云像元。在假彩色合成影像（RGB=b2/b4/b6）中，積雪、林地積雪和云表現(xiàn)出不同的色彩，可以通過目視解譯直接分辨，如圖1所示。本文分別選取森林積雪像元和云像元進(jìn)行統(tǒng)計(jì)分析，發(fā)現(xiàn)二者第四波段地表反射率（sur_refl_bo4），即綠光波段具有明顯的差異。森林積雪像元在該波段值地表反射率大多分布在0-0.4之間；而云像元在綠光波段的反射率大于0.4，如圖3所示。故采用0.4作為閾值進(jìn)行判別，若綠光波段地表反射率小于或等于0.4，則該像元為非云像元，恢復(fù)其Raw NDSI值；若該波段地表反射率大于0.4，則該像元仍標(biāo)記為云。統(tǒng)計(jì)2018年數(shù)據(jù)，經(jīng)過此步驟圖中云量占比大約降至20%。此步驟消除了東北地區(qū)積雪產(chǎn)品中的過度云掩膜現(xiàn)象，可以大大提高后續(xù)計(jì)算效率，而且恢復(fù)的Raw NDSI為像元的NDSI真值，保證了生成的時(shí)空連續(xù)的NDSI序列的精度。

圖3 不同地表下b4波段地表反射率頻率曲線分布圖Figure 3 The curve distribution diagram of Band 4 reflectance frequency under different land surfaces

1.2.4 臨近日合成

一般來說，積雪的累積和消融會(huì)持續(xù)一段時(shí)間，故可以采用臨近日合成的方法，填補(bǔ)部分云下NDSI值。為了保證精度，本文只采用前后兩天的數(shù)據(jù)作為填補(bǔ)數(shù)據(jù)[30,31]。算法規(guī)則為若目標(biāo)像元前后天都是非云像元，則將前后天 NDSI值平均后賦值給當(dāng)天；若前后天中有一天為非云像元，則直接使用該非云像元的NDSI值賦值給當(dāng)天；若前后天都為云像元，則需尋找前后兩天像元繼續(xù)判別，算法規(guī)則同上；若前后兩天同樣也都是云像元，則目標(biāo)像元仍未得到恢復(fù)，計(jì)算結(jié)束，等待下一步去云操作。

1.2.5 時(shí)空自適應(yīng)去云

采用CHEN等[24]提出的時(shí)空自適應(yīng)去云算法，進(jìn)行全圖層徹底去云。其主要思路是判斷相似像元，基于相似像元的NDSI在自適應(yīng)時(shí)空中賦予相應(yīng)的權(quán)重去云。

其中，M(x0,y0,t0)是t0時(shí)刻待填充像元NDSI值，M(xk,yk,t0)是t0時(shí)刻相似像元NDSI值，M(x0,y0,t0+s)和M(xk,yk,t0+s)是t0+s時(shí)刻M(x0,y0)目標(biāo)像元和M(xk,yk)相似像元的NDSI值。Δt0和Δ(t0+s)表示兩個(gè)時(shí)刻相似像元和目標(biāo)像元各自的差異變化。在一定時(shí)期內(nèi)，我們可以認(rèn)為兩者差異變化相同，則得到公式（3）。

基于上述理論，通過時(shí)空相鄰相似像元的NDSI值插補(bǔ)云像元，可將云像元NDSI值得到恢復(fù)。在w×w局部窗口中選擇前n個(gè)NDSI差異最小的像元作為相似像元，其中相似像元差異定義為公式（4），Sk越小，表明兩者差異越小。公式（5）是填充目標(biāo)云像元的加權(quán)平均計(jì)算公式，公式（6）則為相似像元相應(yīng)的求權(quán)重公式，此處采用的是反距離權(quán)重算法。公式（7）中，Dk表示目標(biāo)像元和相似像元之間的相對(duì)距離，其中將Dk歸一化到1-w/2（w是窗口大?。１菊撐牡膮?shù)設(shè)置為n=20，w=15。

1.2.6 精度評(píng)價(jià)

采用云假設(shè)方法對(duì)最終去云結(jié)果進(jìn)行精度評(píng)價(jià)。未被云層覆蓋的真實(shí)地表 NDSI值是已知的，故可以作為驗(yàn)證數(shù)據(jù)。使用云假設(shè)將驗(yàn)證樣本進(jìn)行云掩膜并重新計(jì)算，對(duì)其計(jì)算得到的數(shù)據(jù)與驗(yàn)證數(shù)據(jù)進(jìn)行比較，評(píng)估去云算法的性能和數(shù)據(jù)產(chǎn)品質(zhì)量。本文采用2018年作為驗(yàn)證數(shù)據(jù)集，將積雪季各月云覆蓋率最小的圖像作為真實(shí)影像，使用各月云覆蓋率最大的圖像作為掩膜影像進(jìn)行驗(yàn)證，通過相關(guān)系數(shù)（r）、均方根誤差（RMSE）和平均絕對(duì)誤差（MAE）來評(píng)估產(chǎn)品數(shù)據(jù)的精度[24]。

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

本數(shù)據(jù)集包含2000-2020年中國東北地區(qū)積雪季（每年9月1日至次年4月30日）逐日無云的NDSI時(shí)間序列產(chǎn)品，均為柵格數(shù)據(jù)，空間分辨率為500 m。文件并命名為ChinaNE_ndsi_xxxxxxxx.tif（其中xxxxxxxx表示日期，例如20200101是指2020年1月1日）。其中，數(shù)據(jù)儲(chǔ)存格式為雙精度類型。數(shù)據(jù)集中數(shù)值范圍為-1～1，該值為去云后每個(gè)像元的NDSI值。圖4給出了數(shù)據(jù)集中2018年各月15日的東北積雪區(qū)無云NDSI的分布。需要注意的是，數(shù)據(jù)集中未對(duì)水體進(jìn)行掩膜，因此在使用該產(chǎn)品進(jìn)行積雪制圖時(shí)應(yīng)考慮水體要素。

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

MODIS V6版本從數(shù)據(jù)生產(chǎn)的每個(gè)環(huán)節(jié)進(jìn)行質(zhì)量控制，與之前V5相比，有了很大的改進(jìn)，對(duì)水、云、氣溶膠等進(jìn)行了大氣校正[25,31]。因此，本研究數(shù)據(jù)來源可靠。

從原始數(shù)據(jù)到最終產(chǎn)品的計(jì)算中，每個(gè)處理步驟都確保數(shù)據(jù)準(zhǔn)確性。上下午星合成去云是積雪產(chǎn)品去云常用的方法，主要是考慮同一天的地表具有穩(wěn)定性；決策樹中的閾值也是基于文獻(xiàn)或者統(tǒng)計(jì)特征選取，并且選取了大量執(zhí)行決策樹之后的結(jié)果與假彩色合成影像進(jìn)行了對(duì)比，該步驟可以有效恢復(fù)NDSI_snow_cover數(shù)據(jù)集中一些由森林積雪誤分成云的問題。最終采用CHEN等的方法進(jìn)行完全去云，該方法在原文中經(jīng)過云假設(shè)檢驗(yàn)證明了其精度，去云結(jié)果與原始數(shù)據(jù)的平均相關(guān)系數(shù)r為0.94，平均RMSE為0.12，平均MAE為0.09[24]。本文對(duì)最終產(chǎn)品也經(jīng)過云假設(shè)檢驗(yàn)，平均相關(guān)系數(shù)r、平均RMSE和平均MAE分別為0.95、0.10和0.06。最后對(duì)2018年逐日的數(shù)據(jù)產(chǎn)品均進(jìn)行了目視檢驗(yàn)，結(jié)果表明該數(shù)據(jù)集NDSI分布合理，可以作為東北地區(qū)積雪研究的基礎(chǔ)數(shù)據(jù)源。

圖4 積雪年2018年各月15日NDSI產(chǎn)品展示圖Figure 4 NDSI product display on the 15th day of each month in 2018 Snow Cover year

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

積雪是全球氣候變化的指示器，是全球氣候變化研究的關(guān)鍵變量，因此探究積雪時(shí)空變化受到眾多學(xué)者的廣泛關(guān)注。在現(xiàn)階段遙感積雪信息提取中，仍主要使用的是歸一化差值積雪指數(shù)（NDSI）方法。但當(dāng)前MODIS 積雪產(chǎn)品V6版本在東北地區(qū)，尤其是在東北森林地區(qū)存在過度云掩膜，使該地區(qū)有效的NDSI指數(shù)存在大范圍和長時(shí)間的數(shù)據(jù)值空缺，嚴(yán)重影響了MODIS積雪產(chǎn)品在該地區(qū)的實(shí)際應(yīng)用。本文針對(duì)MOD/MYD10A1的Raw_NDSI數(shù)據(jù)層中大量云像元帶來的數(shù)據(jù)空缺，采用多種有效的方法進(jìn)行數(shù)據(jù)恢復(fù)，生產(chǎn)了2000-2020年中國東北地區(qū)500 m分辨率積雪季逐日無云時(shí)間序列 NDSI數(shù)據(jù)。較已有的積雪產(chǎn)品，本數(shù)據(jù)集在去云和林區(qū)積雪識(shí)別方面有更好的精度，且有較長的時(shí)間跨度和較高的時(shí)間連續(xù)性，可為后續(xù)積雪二值提取，積雪混合像元分解等后續(xù)積雪產(chǎn)品的研究和生產(chǎn)工作提供數(shù)據(jù)源支撐。同時(shí)本數(shù)據(jù)集對(duì)東北地區(qū)積雪時(shí)空變化、森林水文過程、農(nóng)業(yè)生態(tài)、碳循環(huán)過程和氣候環(huán)境變化等研究具有重要的科學(xué)意義。

5 數(shù)據(jù)使用方法和建議

中國東北長時(shí)間序列逐日無云歸一化差值積雪指數(shù)（NDSI）數(shù)據(jù)集，保存為 TIFF格式，能夠在ArcGIS、ENVI及MATLAB等相關(guān)軟件中對(duì)數(shù)據(jù)進(jìn)行讀取、查看、編輯及后續(xù)的一系列統(tǒng)計(jì)分析工作。需要說明的是本數(shù)據(jù)集并未對(duì)水體進(jìn)行掩膜。此外，考慮到本數(shù)據(jù)集服務(wù)于后續(xù)的積雪產(chǎn)品生產(chǎn)工作，本文保留了東北地區(qū)邊界外的數(shù)值，即并未使用東北積雪區(qū)矢量邊界對(duì)數(shù)據(jù)集進(jìn)行裁剪工作。數(shù)據(jù)集中圖像所覆蓋的范圍完全包含了東北積雪區(qū)。由于本文采用的自適應(yīng)時(shí)空去云算法，需要用到前后日期的影像數(shù)據(jù)，因此導(dǎo)致本數(shù)據(jù)集中積雪季起始時(shí)期9月初和結(jié)束時(shí)期次年4月末的數(shù)據(jù)質(zhì)量存在一定問題，但是考慮9月初和4月末中國東北地區(qū)基本無積雪分布，因此不影響本數(shù)據(jù)集的使用。