肖鵬峰，胡瑞，張正，秦棽，李震

1.南京大學(xué)地理與海洋科學(xué)學(xué)院，自然資源部國土衛(wèi)星遙感應(yīng)用重點(diǎn)實(shí)驗(yàn)室，江蘇省地理信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，南京 210023

2.中國科學(xué)院空天信息創(chuàng)新研究院，北京 100094

引 言

地表反照率指地表向各個方向反射的全部光通量與總?cè)肷涔馔康谋?。它量化了大氣與地表之間輻射的相互作用，是重要的氣候參數(shù)[1-2]，被廣泛用在地表能量平衡、中長期天氣預(yù)測和全球變化研究中[3]。積雪具有顯著的高反照率特性[4]，在局部或全球能量收支平衡中起著重要作用。近幾十年來北半球年平均輻射強(qiáng)迫每年增加約0.45 W/m2，其中積雪帶來的影響將近一半[5]。積雪反照率調(diào)節(jié)著積雪能量收支，進(jìn)一步控制積雪的消融過程，從而影響區(qū)域乃至全球的水循環(huán)過程。同時，積雪為世界上17%的人口提供飲用水及生產(chǎn)生活用水，在部分山區(qū)流域積雪融水占到年徑流量的50%以上[6]，是綠洲灌溉的寶貴水源[7]。此外，積雪污化物與環(huán)境污染相關(guān)，其濃度上升會導(dǎo)致積雪反照率顯著下降[8]。因此，積雪污化物的輻射強(qiáng)迫研究對人類污染排放與大氣活動有一定的指示作用[9]。

目前，國內(nèi)外主要有4種遙感反照率產(chǎn)品：美國宇航局16天合成地表反照率產(chǎn)品MCD43A3、逐日積雪反照率產(chǎn)品SAD、歐洲空間局八天合成地表反照率產(chǎn)品GlobAlbedo、北京師范大學(xué)八天合成地表反照率產(chǎn)品 GLASS。MCD43A3產(chǎn)品是中分辨率光譜成像儀（Moderate Resolution Imaging Spectroradiometer，MODIS）反照率產(chǎn)品中的一部分，采用的反演模型為半經(jīng)驗(yàn)核驅(qū)動模型（RossThick-LiSparse-Reciprocal，RTLSR），根據(jù)體散射核、幾何光學(xué)核和各向同性核的加權(quán)和來估算地表BRDF[10-11]。SAD產(chǎn)品是MODIS積雪產(chǎn)品MOD10A1/MYD10A1中的一部分，在球形雪粒假設(shè)下，利用離散縱坐標(biāo)輻射傳輸（Discrete Ordinates Radiative Transfer，DISORT）模型構(gòu)建等效光學(xué)粒徑為250 μm的各向異性響應(yīng)函數(shù)（Anisotropic Response Function，ARF）查找表，經(jīng)過窄寬轉(zhuǎn)換計算得到積雪黑空反照率[12]。GlobAlbedo產(chǎn)品的反演模型為半經(jīng)驗(yàn)核驅(qū)動模型RTLSR，經(jīng)過計算得到地表反照率[13-14]。GLASS產(chǎn)品采用直接估算方法，通過訓(xùn)練建立適用于全球不同區(qū)域、不同覆蓋類型、不同季節(jié)的格網(wǎng)查找表，直接估算得到地表反照率[15-17]。

然而，這些反照率產(chǎn)品應(yīng)用于積雪研究仍然存在問題。MCD43A3、GlobAlbedo和GLASS這三種地表反照率產(chǎn)品主要是基于植被的輻射傳輸特性生產(chǎn)的，不能充分反映積雪的輻射傳輸特性。專門針對積雪的全球反照率產(chǎn)品僅有SAD，但它在反演過程中未考慮雪粒徑和污化物濃度參數(shù)對積雪各向異性特性的影響，且缺乏白空反照率數(shù)據(jù)。同時，由于積雪反照率存在快速變化的特點(diǎn)，充分捕獲積雪反照率的時間變化特征需要高時間分辨率的積雪反照率產(chǎn)品。當(dāng)前的大部分產(chǎn)品采用了八天或16天合成的方法，如此低的時間分辨率無法表征積雪積累、消融過程中反照率快速變化的特征，導(dǎo)致這些反照率產(chǎn)品在積雪區(qū)域的精度欠佳[12,18]。而且，由于全球反照率產(chǎn)品驗(yàn)證數(shù)據(jù)主要分布在美洲、歐洲，在中國區(qū)域數(shù)量較少，因此這些產(chǎn)品在中國區(qū)域缺少系統(tǒng)性評價，產(chǎn)品選擇和使用的不確定性較大。

本文基于積雪的輻射傳輸模型發(fā)展了適應(yīng)中國區(qū)域特點(diǎn)的反照率反演算法，以 MODIS地表反射率產(chǎn)品作為輸入，生產(chǎn)了2000-2020年中國積雪反照率的逐日和八天合成去云的遙感產(chǎn)品。利用地面站點(diǎn)和樣方的積雪反照率觀測數(shù)據(jù)，對產(chǎn)品開展了精度驗(yàn)證分析，結(jié)果表明產(chǎn)品的均方根誤差、相關(guān)系數(shù)、平均偏差均優(yōu)于國際同類產(chǎn)品。產(chǎn)品可為我國積雪資源分布調(diào)查提供數(shù)據(jù)支持，為積雪類型劃分提供時空參數(shù)，為地表輻射研究和陸面過程模擬提供科學(xué)支撐，滿足我國輻射、氣候、環(huán)境、生態(tài)等相關(guān)研究的應(yīng)用需求。

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

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

本文主要基于MODIS地表反射率反演積雪反照率。因此，數(shù)據(jù)源包括MOD09GA地表反射率產(chǎn)品、積雪范圍、MCD12Q1土地覆蓋產(chǎn)品、SRTM數(shù)字高程數(shù)據(jù)、林區(qū)積雪反照率觀測數(shù)據(jù)共5類。

（1）MOD09GA地表反射率產(chǎn)品是對MODIS波段1到波段7的輻亮度數(shù)據(jù)進(jìn)行大氣校正處理，得到的500 m分辨率的反射率數(shù)據(jù)集和1 km分辨率的觀測幾何數(shù)據(jù)集[19]。本文使用2000-2020年逐日MOD09GA數(shù)據(jù)，最新版本為Version 6，它改進(jìn)了氣溶膠反演算法，構(gòu)建新的氣溶膠反演查找表，同時對數(shù)據(jù)內(nèi)部雪、云以及云影等檢測算法進(jìn)行了改進(jìn)。

（2）積雪范圍來自中國MODIS逐日無云500 m積雪面積產(chǎn)品[20]。它基于MOD09GA地表反射率產(chǎn)品，結(jié)合積雪野外調(diào)查數(shù)據(jù)、數(shù)字高程數(shù)據(jù)、土地覆蓋數(shù)據(jù)，通過決策樹分類方法生成了中國2000-2020年積雪面積產(chǎn)品，與MODIS積雪產(chǎn)品MOD10A1/MYD10A1相比具有更高的精度。

（3）MCD12Q1土地覆蓋產(chǎn)品提供每年的全球土地覆蓋類型[21]。它通過對MODIS地表反射率數(shù)據(jù)進(jìn)行監(jiān)督分類得到初始分類結(jié)果，利用先驗(yàn)知識和輔助信息進(jìn)行后處理，得到包含六種不同分類標(biāo)準(zhǔn)的土地覆蓋數(shù)據(jù)。本文使用最新Version 6版本的MCD12Q1數(shù)據(jù)，因目前僅有2001-2019年的逐年產(chǎn)品，故2000年的土地覆蓋用2001年代替，2020年的土地覆蓋用2019年代替。選擇其中的LC_Type1分類標(biāo)準(zhǔn)，提取森林范圍。

（4）數(shù)字高程數(shù)據(jù)采用2000年的SRTM數(shù)據(jù)[22]。SRTM雷達(dá)系統(tǒng)搭載于奮進(jìn)號航天飛機(jī)，采用干涉合成孔徑雷達(dá)技術(shù)生成了全球范圍內(nèi)的數(shù)字高程模型。本文使用經(jīng)過數(shù)據(jù)空白填補(bǔ)的最新Version 4版本，分辨率為90 m。

（5）林區(qū)積雪反照率數(shù)據(jù)來自長白山和小興安嶺森林生態(tài)監(jiān)測站的塔臺觀測。長白山森林生態(tài)監(jiān)測站位于吉林省安圖縣二道白河鎮(zhèn)，位置為128.47 E、42.40 N，海拔736 m。小興安嶺森林生態(tài)監(jiān)測站位于伊春市五營區(qū)，位置為129.20 E、48.23 N，海拔549 m。兩個站都架設(shè)有森林觀測塔，觀測塔下方為長勢均一的次生林，與周圍情況一致，具有較好的空間代表性。在觀測塔30 m高處架設(shè)四分量輻射表，自動獲取長時間序列的林上反照率連續(xù)觀測數(shù)據(jù)，時間范圍分別為 2003-2010年2017-2019年。

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

首先基于積雪范圍產(chǎn)品對反射率進(jìn)行掩膜，獲得積雪反射率；對積雪反射率進(jìn)行地形校正、各向異性校正，獲得窄波段反照率；經(jīng)過窄寬轉(zhuǎn)換，得到寬波段積雪反照率；基于土地覆蓋數(shù)據(jù)獲得林區(qū)范圍，對林區(qū)積雪反照率進(jìn)行校正；經(jīng)過異常值處理，獲得逐日積雪反照率產(chǎn)品；進(jìn)行八天均值合成，獲得去云積雪反照率產(chǎn)品。因此，數(shù)據(jù)處理方法和流程包括預(yù)處理、反照率計算、林區(qū)校正、異常值處理、任務(wù)運(yùn)行、后處理共六個部分（圖1）。

（1）預(yù)處理。基于MODIS地表反射率數(shù)據(jù)獲得積雪覆蓋范圍和云覆蓋區(qū)域，然后對反射率進(jìn)行掩膜，獲得積雪反射率。將數(shù)字高程數(shù)據(jù)和其他數(shù)據(jù)的投影轉(zhuǎn)換至MOD09GA的坐標(biāo)系統(tǒng)，并通過重采樣統(tǒng)一空間分辨率。利用數(shù)字高程數(shù)據(jù)，采用C校正模型對積雪反射率進(jìn)行地形校正，降低山區(qū)起伏地形對反照率反演的影響。

圖1 中國積雪反照率遙感產(chǎn)品生產(chǎn)流程Figure 1 Production process of remote sensing product of snow albedo over China

（2）反照率計算。采用基于積雪弱吸收特性建立的漸近輻射傳輸（Asymptotic Radiative Transfer，ART）模型[4]，它可以較好地模擬積雪的雙向反射特性[23]。因此，使用ART模型對積雪反射率進(jìn)行各向異性校正，獲得窄波段反照率。計算積雪的窄波段黑空和白空反照率的公式分別為：

式中，rb(θ0)和rw分別為積雪黑空和白空反照率，θ、θ0、φ分別為觀測天頂角、太陽天頂角、相對方位角，R(θ,θ0,φ)為積雪反射率，f(θ,θ0,φ)為積雪散射特征函數(shù)，K0(θ0)為半無限-無吸收介質(zhì)中光子溢出歸一化角分布函數(shù)，R0(θ,θ0,φ)為半無限-無吸收介質(zhì)的反射函數(shù)，計算公式分別如下：

式中，P(Ω)和Ω的計算公式分別為：

最后，采用針對積雪高反照率特性得出的窄寬轉(zhuǎn)換公式[24]，獲得寬波段積雪反照率：

式中，α為寬波段反照率，α1到α7分別為波段1到波段7的窄波段反照率。

（3）林區(qū)校正。林區(qū)積雪存在冠層截留、林下堆積、多次散射等復(fù)雜情況，反照率比純雪低得多。由于缺少實(shí)用的森林積雪輻射傳輸模型，利用線性擬合方式進(jìn)行林區(qū)積雪反射率校正。利用土地覆蓋數(shù)據(jù)獲得林區(qū)范圍，得到下墊面類型為森林的積雪區(qū)域?；谏稚鷳B(tài)監(jiān)測站的反照率觀測數(shù)據(jù)，通過格網(wǎng)搜索得到林區(qū)積雪反照率校正系數(shù)，利用該校正系數(shù)對所有的林區(qū)積雪反照率像元進(jìn)行校正。

（4）異常值處理。當(dāng)反射率產(chǎn)品中的云雪標(biāo)識錯誤時，誤判為雪的云會造成該像元反照率值顯著高于正確識別為雪的像元反照率值。因此，通過格網(wǎng)搜索確定閾值，檢測并剔除積雪反照率中異常高值的像元，得到正確的積雪反照率像元。

（5）任務(wù)運(yùn)行。產(chǎn)品生產(chǎn)采取云端（Google Earth Engine）和本地端（Python）交互生產(chǎn)的方式，以最大限度地提高生產(chǎn)效率。在本地端編寫腳本，提交生產(chǎn)任務(wù)至云端；在云端讀取源數(shù)據(jù)，生產(chǎn)積雪反照率產(chǎn)品，下載至本地進(jìn)行后處理。

（6）后處理。通過八天均值合成的算法，獲得去云的積雪反照率產(chǎn)品。最后將產(chǎn)品由tiff格式轉(zhuǎn)換為hdf格式，以提高存儲和使用的靈活性。并生成png格式的產(chǎn)品縮略圖，以供預(yù)覽。

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

基于上述數(shù)據(jù)處理方法和流程，生產(chǎn)得到的 2000-2020年中國積雪反照率遙感產(chǎn)品，命名為ChinaSA（China’s Snow Albedo）。數(shù)據(jù)示例如圖2，數(shù)據(jù)波段信息見表1。

產(chǎn)品的數(shù)據(jù)命名規(guī)則如下：PLT_MODIS_SAB_YYYYMMDD_RESOL_XXXXM_VXX.HDF，其中 PLT表示衛(wèi)星平臺，MODIS表示傳感器，SAB表示積雪反照率，YYYYMMDD表示年月日，RESOL表示時間分辨率，XXXXM表示以米為單位的空間分辨率，VXX表示版本號。

圖1 中國2020年1月1日積雪反照率遙感產(chǎn)品示意圖（審圖號：GS(2022)747號）Figure 2 Remote sensing product of snow albedo over China on January 1, 2020

表1 中國積雪反照率遙感產(chǎn)品波段信息Table 1 Band information of remote sensing product of snow albedo over China

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

3.1 驗(yàn)證方法

利用多種地面站點(diǎn)的積雪反照率觀測數(shù)據(jù)，結(jié)合積雪反照率樣方觀測數(shù)據(jù)，與遙感反演的積雪反照率產(chǎn)品進(jìn)行對比，計算均方根誤差（RMSE）、相關(guān)系數(shù)（R）、平均偏差（AD），驗(yàn)證積雪反照率產(chǎn)品的精度。

收集的地面站點(diǎn)觀測數(shù)據(jù)包括：中國氣象輻射基本要素逐小時曝幅量數(shù)據(jù)集、全球協(xié)調(diào)加強(qiáng)觀測計劃（CEOP）亞澳季風(fēng)青藏高原試驗(yàn)研究（CAMP-Tibet）數(shù)據(jù)集、中國陸地生態(tài)系統(tǒng)通量觀測研究聯(lián)盟（ChinaFLUX）數(shù)據(jù)集、中國典型積雪區(qū)積雪站積雪特性觀測數(shù)據(jù)集、森林生態(tài)監(jiān)測站數(shù)據(jù)集。利用高分辨率遙感影像對站點(diǎn)逐一進(jìn)行了空間代表性分析，發(fā)現(xiàn)擁有反射輻射數(shù)據(jù)的19個氣象站全部位于城區(qū)，觀測場的實(shí)測值高估了混合像元的反照率，因此主要選用另外三種數(shù)據(jù)進(jìn)行驗(yàn)證。CAMP-Tibet數(shù)據(jù)集包含2002-2004年7個站點(diǎn)的觀測數(shù)據(jù)，ChinaFLUX數(shù)據(jù)集包含2003-2010年10個站點(diǎn)的觀測數(shù)據(jù)，積雪站數(shù)據(jù)集包含2017-2020年6個站點(diǎn)的觀測數(shù)據(jù)。森林生態(tài)監(jiān)測站數(shù)據(jù)集包含2003-2019年2個站點(diǎn)的觀測數(shù)據(jù)，隨機(jī)抽取一半數(shù)據(jù)做校正，另一半數(shù)據(jù)做驗(yàn)證。站點(diǎn)數(shù)據(jù)的篩選規(guī)則為：（1）選取北京時間上午10:30左右的地面實(shí)測數(shù)據(jù)，以保證地面實(shí)測時間與衛(wèi)星過境時間一致；（2）地面實(shí)測與遙感反演的反照率值皆大于0且小于1，以去除無效數(shù)據(jù)；（3）站點(diǎn)所在像元當(dāng)天的積雪覆蓋率大于95%，保證觀測場的實(shí)測值基本能代表整個像元。經(jīng)過以上篩選，共得到242個站點(diǎn)驗(yàn)證數(shù)據(jù)。

站點(diǎn)對像元的驗(yàn)證仍然存在尺度效應(yīng)問題，進(jìn)一步通過樣方觀測提高驗(yàn)證的效度。針對中國主要積雪區(qū)，根據(jù)MODIS像元對應(yīng)的地面范圍，結(jié)合高分辨率影像提前確定樣方范圍，一般選擇下墊面較為平坦的農(nóng)田、戈壁、草原等開闊區(qū)域。設(shè)計便攜式反照率觀測裝置，分純雪像元樣方和混合像元樣方兩種情況進(jìn)行觀測。純雪像元樣方的反照率為中心點(diǎn)的實(shí)測值，混合像元樣方的反照率為不同地物的實(shí)測值與無人機(jī)獲取的面積比例進(jìn)行加權(quán)求和的結(jié)果。測得的樣方數(shù)據(jù)包括：2019年1月北疆地區(qū)67個樣方、2020年1月東北-內(nèi)蒙古地區(qū)73個樣方。樣方數(shù)據(jù)的篩選規(guī)則為：樣方所在像元當(dāng)天被產(chǎn)品標(biāo)識為積雪。利用該規(guī)則進(jìn)行篩選，共得到82個樣方驗(yàn)證數(shù)據(jù)。

篩選后的驗(yàn)證數(shù)據(jù)空間分布如圖3所示，覆蓋了東北-內(nèi)蒙古（A）、青藏高原（B）、北疆（C）三大典型積雪區(qū)。圖中的圓點(diǎn)表示站點(diǎn)驗(yàn)證數(shù)據(jù)，方形表示樣方驗(yàn)證數(shù)據(jù)，數(shù)字表示各站點(diǎn)、樣方用于驗(yàn)證的數(shù)據(jù)記錄條數(shù)。

圖3 中國積雪反照率遙感產(chǎn)品驗(yàn)證數(shù)據(jù)空間分布（審圖號：GS(2022)747號）Figure 3 Spatial distribution of validation data for remote sensing product of snow albedo over China

3.2 驗(yàn)證結(jié)果與精度比較

利用324個站點(diǎn)與樣方實(shí)測數(shù)據(jù)驗(yàn)證積雪反照率產(chǎn)品ChinaSA的精度，并與美國宇航局16天合成地表反照率產(chǎn)品 MCD43A3、逐日積雪反照率產(chǎn)品 SAD、歐空局八天合成地表反照率產(chǎn)品GlobAlbedo、北京師范大學(xué)八天合成地表反照率產(chǎn)品GLASS的精度進(jìn)行對比，結(jié)果如圖4所示。圖中左下方為林區(qū)實(shí)測數(shù)據(jù)驗(yàn)證結(jié)果，右側(cè)為非林區(qū)實(shí)測數(shù)據(jù)的驗(yàn)證結(jié)果，可見ChinaSA產(chǎn)品在非林區(qū)和林區(qū)都與實(shí)測數(shù)據(jù)具有更好的一致性。進(jìn)一步計算五種產(chǎn)品的精度指標(biāo)，如表2所示。ChinaSA產(chǎn)品與實(shí)測值的均方根誤差為0.114，相關(guān)系數(shù)為0.822，平均偏差為0.087，優(yōu)于國際同類產(chǎn)品。

圖4 五種反照率產(chǎn)品的積雪反照率與實(shí)測值對比圖Figure 4 Comparison of snow albedo of the five products and the in-situ data

表2 五種反照率產(chǎn)品的積雪反照率精度Table 2 Snow albedo accuracies of the five products

積雪反照率精度的提升，主要源于以下三方面的貢獻(xiàn)：（1）ChinaSA產(chǎn)品采用專門針對積雪弱吸收特性建立的ART模型，能夠更好地表達(dá)積雪的各向異性反射特性；其他地表反照率產(chǎn)品主要采用針對植被建立的輻射傳輸模型，在積雪區(qū)域的精度較低；積雪反照率產(chǎn)品SAD需要較多的輸入?yún)?shù)，限制了模型的適用性；（2）利用我國森林生態(tài)監(jiān)測站獲取校正系數(shù)，顯著提升了森林區(qū)域的積雪反照率反演精度；（3）通過地形校正處理，進(jìn)一步提升了起伏地表的積雪反照率反演精度。

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

本數(shù)據(jù)集為中國2000-2020年積雪反照率遙感產(chǎn)品，可為我國積雪資源分布調(diào)查提供數(shù)據(jù)支持，滿足我國輻射、氣候、環(huán)境、生態(tài)等相關(guān)研究的應(yīng)用需求。

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

數(shù)據(jù)可通過國家冰川凍土沙漠科學(xué)數(shù)據(jù)中心（www.ncdc.ac.cn）下載。數(shù)據(jù)下載后可用ArcGIS/ENVI等軟件或者相應(yīng)程序代碼讀取，數(shù)據(jù)中的波段與本文表1中的波段相對應(yīng)。

致 謝

感謝南京大學(xué)馮學(xué)智教授、張學(xué)良副教授和研究生葉李灶、馬威、劉豪、盛光偉、馬騰耀、宋依娜、李振世等參與北疆、青藏高原、東北-內(nèi)蒙古等地的積雪反照率野外觀測驗(yàn)證以及首席科學(xué)家王建研究員等專家對產(chǎn)品生產(chǎn)的指導(dǎo)、支持與幫助。