李曉東，李林，史飛飛,4*，蘇文將，肖建設(shè)，李紅梅,5

1.青海省氣象科學(xué)研究所，西寧 810008

2.青海省防災(zāi)減災(zāi)重點(diǎn)實(shí)驗(yàn)室，西寧 810008

3.青海省氣象局，西寧 810008

4.青海師范大學(xué)地理科學(xué)學(xué)院，西寧 810008

5.青海省氣候中心，西寧 810008

引 言

第六次IPCC評估報(bào)告指出，隨著全球變暖的加劇，強(qiáng)降水、高溫?zé)崂说雀哂绊憵夂蚴录⒏鼮轭l繁[1]。被譽(yù)為世界“第三極”的青藏高原，對全球氣候變化的反饋具有一定的“前瞻性”。近年來青藏高原部分地區(qū)異常降雪過程增多，積雪狀況也發(fā)生了顯著而又獨(dú)特的變化[2-3]。青藏高原作為全球典型的高海拔積雪分布區(qū)以及中國和周邊諸多河流的發(fā)源地，積雪消融及其變化狀況對該地區(qū)乃至全球尺度的水循環(huán)、氣候變化進(jìn)程將產(chǎn)生較大影響，因而準(zhǔn)確獲積雪變化信息具有十分重要意義[4-6]。遙感具有宏觀、高效和高頻次等優(yōu)勢，能夠快速實(shí)現(xiàn)積雪的動態(tài)監(jiān)測。在近五十余年的積雪遙感監(jiān)測過程中，國內(nèi)外學(xué)者利用 SMMR、SSM、AMSR-E等微波傳感器，以及 AVHRR、MODIS、Suomi NPP和FY-3 MERSI等光學(xué)傳感器研制了多種全球和區(qū)域尺度的積雪覆蓋、積雪深度和雪水當(dāng)量產(chǎn)品[7-10]，但受青藏高原云系旺盛、積雪斑塊化嚴(yán)重以及積雪消融較快等因素影響，現(xiàn)有的全球積雪產(chǎn)品精度在高原地區(qū)普遍不夠理想[11]。為提高積雪產(chǎn)品在青藏高原的監(jiān)測精度，黃曉東等[12]、于金媛等[13]和邱玉寶等[14]均基于MODIS研制了逐日無云積雪產(chǎn)品，唐志光等[15]發(fā)展了MODIS積雪面積比例產(chǎn)品，上述積雪產(chǎn)品的監(jiān)測精度均得到較大程度提高并在積雪評估研究中得到廣泛應(yīng)用。

青藏高原雪層普遍偏淺且呈非連續(xù)性分布，在反映積雪類型分區(qū)和分布特征方面積雪日數(shù)相比其他積雪監(jiān)測數(shù)據(jù)更具有代表性[16-17]，但目前積雪日數(shù)多以地面氣象站的積雪觀測資料為主，受青藏高原地形復(fù)雜和氣象站點(diǎn)稀少等因素影響其對積雪信息描述存在較大不確定性。本數(shù)據(jù)集利用2002-2021年MODIS V6版積雪產(chǎn)品中的歸一化積雪指數(shù)（Normalized Difference Snow Index，NDSI）等數(shù)據(jù)，通過逐日積雪覆蓋提取、積雪覆蓋范圍校正、雪季內(nèi)積雪日數(shù)合成、冰川區(qū)積雪日數(shù)修正和氣象站點(diǎn)數(shù)據(jù)校正等步驟，制備了青海省積雪日數(shù)遙感監(jiān)測數(shù)據(jù)集，將為青藏高原地區(qū)的積雪分布及雪水儲量評估、積雪變化分析和雪災(zāi)風(fēng)險(xiǎn)區(qū)劃等研究工作提供數(shù)據(jù)支撐。

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

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

青海省積雪日數(shù)數(shù)據(jù)集的制作流程主要分為8個(gè)步驟：獲取MODIS積雪產(chǎn)品、積雪覆蓋提取、積雪覆蓋范圍校正（包括：上下午數(shù)據(jù)雪蓋合成、臨近三天雪蓋合成）、雪季內(nèi)積雪日數(shù)合成、冰川區(qū)積雪日數(shù)修正、基于氣象站點(diǎn)的積雪日數(shù)校正、精度驗(yàn)證和構(gòu)建積雪日數(shù)數(shù)據(jù)集，具體流程如圖1所示。其中輸入的原始數(shù)據(jù)為2002年10月至2021年5月的MODIS V6版積雪數(shù)據(jù)，精度驗(yàn)證采用地面氣象臺站的人工觀測數(shù)據(jù)，最終輸出各積雪季內(nèi)青海省的積雪日數(shù)數(shù)據(jù)。

1.2 數(shù)據(jù)獲取

積雪日數(shù)數(shù)據(jù)集的構(gòu)建選用 MODIS積雪產(chǎn)品，該產(chǎn)品在青藏高原地區(qū)的積雪研究中廣泛使用[4]，其500 m的空間分辨率以及逐日的觀測頻次能夠滿足高原地區(qū)積雪特性監(jiān)測的實(shí)際需求。MODIS積雪產(chǎn)品已更新至第六版，本文結(jié)合Terra和Auqa雙星積雪數(shù)據(jù)，其產(chǎn)品編號分別為MOD10A1和MYD10A1，在谷歌地球引擎（Google Earth Engine，GEE）中獲取2002-2021年青海省的逐日積雪數(shù)據(jù)。

圖1 數(shù)據(jù)處理流程Figure 1 Data processing flow

1.3 積雪覆蓋提取方法

此前在MODIS V5版積雪產(chǎn)品中是利用NDSI進(jìn)行積雪覆蓋的二值判識，其閾值設(shè)定為0.4，但在MODIS V6版積雪產(chǎn)品中不再提供積雪覆蓋的直接判識結(jié)果，轉(zhuǎn)而提供NDSI指數(shù)。ZHANG等[18]和高楊等[19]研究指出在對青藏高原進(jìn)行積雪覆蓋判識時(shí)將NDSI的閾值設(shè)定為0.4會存在嚴(yán)重漏判，而將NDSI的閾值設(shè)定為0.1較0.4更為合理。本文在此研究基礎(chǔ)上，通過前期對MODIS影像中積雪覆蓋目視判識和進(jìn)行野外觀測對比，獲取了青海省積雪判識時(shí) NDSI的經(jīng)驗(yàn)性閾值，并最終采用分段閾值的方法開展青海省積雪覆蓋范圍提取，具體方法見公式（1）。在GEE平臺下通過編寫腳本實(shí)現(xiàn)了逐日積雪覆蓋提取，同時(shí)利用MODIS積雪產(chǎn)品中的云（250）、內(nèi)陸水體（237）和無數(shù)據(jù)（200）等質(zhì)量控制編碼值對每日MOD/MYD10A1數(shù)據(jù)中的非積雪信息進(jìn)行掩膜。

式中：B2為0.55-0.57 um可見光波段的反射率，B1為1.62-1.65 um短波紅外波段的反射率，NDSIth為判定閾值，通過實(shí)驗(yàn)將青海省內(nèi)10-12月的NDSIth設(shè)定為0.1，1-3月的NDSIth設(shè)定為0.13，4-5月的NDSIth設(shè)定為0.08。

1.4 積雪覆蓋范圍校正

上下午數(shù)據(jù)雪蓋合成是指對MODIS/Terra上午星積雪數(shù)據(jù)（MOD10A1）判識的積雪覆蓋范圍和MODIS/Aqua下午星積雪數(shù)據(jù)（MYD10A1）判識的積雪覆蓋范圍進(jìn)行最大合成，具體是將一日內(nèi)兩幅積雪判識后的二值柵格圖像進(jìn)行疊加，將每個(gè)柵格單元值取兩幅積雪判識結(jié)果中的最大值，最后合成一日內(nèi)具有最大積雪覆蓋范圍的二值數(shù)據(jù)。

臨近三天雪蓋合成是指在較短時(shí)間內(nèi)積雪未發(fā)生快速消融，在當(dāng)天積雪監(jiān)測結(jié)果中有柵格像元受云影響無法判識是否存在積雪時(shí)，如前天和后天該像元上有積雪存在，則認(rèn)定當(dāng)天依舊有雪覆蓋。

1.5 積雪日數(shù)合成

積雪覆蓋日數(shù)合成是將每日的積雪覆蓋提取結(jié)果依據(jù)使用要求進(jìn)行旬/月/季/年或任意時(shí)段的合成，具體是將每日內(nèi)積雪覆蓋的二值柵格數(shù)據(jù)進(jìn)行空間疊加，將每個(gè)柵格單元值進(jìn)行累加后，便可生成統(tǒng)計(jì)時(shí)段內(nèi)的累計(jì)積雪日數(shù)柵格數(shù)據(jù)。本數(shù)據(jù)集是按積雪季進(jìn)行合成，其中當(dāng)年積雪季一般從上年10月開始至當(dāng)年5月結(jié)束。

1.6 冰川區(qū)積雪日數(shù)修正

青藏高原冰川分布區(qū)海拔高且常年有積雪覆蓋，但上述區(qū)域經(jīng)常受云遮擋導(dǎo)致部分地區(qū)積雪日數(shù)的遙感監(jiān)測結(jié)果偏低，針對該問題選用中國第二次冰川數(shù)據(jù)集[20]獲取青海省內(nèi)冰川分布區(qū)，并對上述區(qū)域內(nèi)的積雪日數(shù)進(jìn)行校正，即整個(gè)雪季內(nèi)冰川區(qū)均有積雪覆蓋。

1.7 站點(diǎn)積雪日數(shù)校正

利用 MODIS積雪數(shù)據(jù)初步獲取的積雪日數(shù)會受到云的干擾造成一定誤差，可利用地面觀測的積雪日數(shù)對其進(jìn)行校正[21]。其中積雪日數(shù)觀測數(shù)據(jù)來源于地面氣象臺站，從中國氣象局綜合氣象信息共享平臺（China Integrated Meteorological Information Sharing System，CIMISS）獲取并已進(jìn)行質(zhì)量控制，觀測員于每日8時(shí)定時(shí)開展觀測，如在后期有降雪過程則在14時(shí)和20時(shí)進(jìn)行補(bǔ)測，規(guī)定在積雪觀測場內(nèi)積雪覆蓋范圍超過一半，且多次測量的積雪深度均≥2 cm時(shí)則記錄當(dāng)天積雪日數(shù)為1。在校正方法的具體實(shí)施時(shí)，首先按3:2比例隨機(jī)將全省氣象臺站劃分為建模站點(diǎn)（31個(gè)）和驗(yàn)證站點(diǎn)（20個(gè)），其次需獲取整個(gè)積雪季內(nèi)（上年10月至當(dāng)年5月）建模氣象臺站觀測的積雪日數(shù)和與氣象臺站空間位置對應(yīng)處MODIS提取的遙感積雪日數(shù)，對上述兩組時(shí)序數(shù)據(jù)按公式(2)進(jìn)行擬合，以獲得校正系數(shù)a和b。最后當(dāng)其擬合效果較好時(shí)則按照公式(3)對MODIS積雪數(shù)據(jù)提取的積雪日數(shù)進(jìn)行逐柵格像元校正。

式中：X為利用MODIS積雪數(shù)據(jù)監(jiān)測的積雪日數(shù)，Y為站點(diǎn)觀測的積雪日數(shù)，單位均為天（d），a、b為校正系數(shù)，Y1為MODIS積雪數(shù)據(jù)監(jiān)測的積雪日數(shù)，Y2為校正后的MODIS積雪數(shù)據(jù)監(jiān)測的積雪日數(shù)。

1.8 精度驗(yàn)證

將積雪日數(shù)數(shù)據(jù)集與驗(yàn)證站點(diǎn)（20個(gè)）的積雪日數(shù)觀測數(shù)據(jù)進(jìn)行對比，以評價(jià)數(shù)據(jù)的準(zhǔn)確性，按下列公式計(jì)算決定系數(shù)(R2)和相對誤差(Yr)，一般R2越接近1，Yr越小則認(rèn)為制備的積雪日數(shù)數(shù)據(jù)精度較高。

式中：Xi為站點(diǎn)積雪覆蓋日數(shù)，Yi為遙感監(jiān)測的積雪覆蓋日數(shù)，Y?為遙感監(jiān)測的積雪覆蓋日數(shù)平均值，X?為站點(diǎn)積雪覆蓋日數(shù)的平均值，n為驗(yàn)證站點(diǎn)個(gè)數(shù)。

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

青海省積雪日數(shù)遙感監(jiān)測數(shù)據(jù)集共由19個(gè)雪季年的柵格數(shù)據(jù)文件組成，為便于數(shù)據(jù)在通用遙感軟件下處理與展示，將數(shù)據(jù)存儲為tif格式，其中文件名以雪季年（上年10月至當(dāng)年5月）命名，空間分辨率為500 m，屬性值為雪季內(nèi)的積雪日數(shù)，單位為天（d），覆蓋范圍為青海省，圖2展示了2017和2019雪季內(nèi)青海省積雪日數(shù)的空間分布狀況。

圖2 2017雪季（a）和2019雪季（b）積雪日數(shù)空間分布圖Figure 2 The spatial distribution of snow cover days in 2017 (a) and 2019 (b) snow season

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

以參與驗(yàn)證的地面臺站人工觀測的積雪日數(shù)作為“真值”，逐積雪季對本文制作的積雪日數(shù)數(shù)據(jù)集進(jìn)行精度檢驗(yàn)，結(jié)果表明在各個(gè)雪季內(nèi)遙感提取的積雪日數(shù)與臺站觀測相關(guān)程度較高，其R2在0.75-0.93，平均相對誤差25%。其中2017和2019雪季的積雪日數(shù)及其精度狀況如圖2、圖3所示，積雪日數(shù)高值區(qū)主要分布在三江源腹地、可可西里地區(qū)和祁連山地區(qū)，而積雪日數(shù)低值區(qū)主要分布在柴達(dá)木盆地，其較好地展示了積雪在空間尺度上的分布差異，并與實(shí)際積雪分布現(xiàn)狀較為相符；同時(shí)，針對異常降雪年份例如2019年，強(qiáng)降雪日數(shù)異常偏多、氣溫普遍較低、積雪持續(xù)時(shí)間長，在青海省內(nèi)瑪多縣、雜多縣、甘德縣和玉樹市等地雪災(zāi)等級為歷史最重，本文積雪日數(shù)數(shù)據(jù)也能較好刻畫積雪覆蓋在空間和時(shí)間尺度上的積累狀況，并對雪災(zāi)發(fā)生的風(fēng)險(xiǎn)也具有很好的指示意義。

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

本數(shù)據(jù)以MODIS V6版積雪產(chǎn)品為基礎(chǔ)，通過逐日積雪覆蓋提取、積雪覆蓋范圍校正、雪季內(nèi)積雪日數(shù)合成和氣象站點(diǎn)數(shù)據(jù)校正等步驟制作了2003-2021雪季內(nèi)青海省的積雪日數(shù)數(shù)據(jù)集，將對開展青藏高原積雪狀況評估、積雪演變特征分析、雪災(zāi)風(fēng)險(xiǎn)評估及預(yù)測等研究工作提供數(shù)據(jù)支撐。

圖3 2017雪季（a）和2019雪季（b）MODIS積雪日數(shù)與站點(diǎn)積雪日數(shù)相關(guān)性圖Figure 3 The diagram of correlation between snow cover days in MODIS and site snow days in 2017 (a) and 2019 snow seasons (b)

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

本數(shù)據(jù)集為tif文件格式，可使用ArcGIS和ENVI等圖像處理軟件進(jìn)行查看、統(tǒng)計(jì)和分析等，各數(shù)據(jù)以積雪季命名，數(shù)據(jù)起始時(shí)間為2003雪季（即2002年10月至2003年5月），數(shù)據(jù)結(jié)束時(shí)間為2021雪季（即2020年10月至2021年5月）。