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

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

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

引 言

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

目前，國(guó)內(nèi)外主要有4種遙感反照率產(chǎn)品：美國(guó)宇航局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ū)動(dòng)模型（RossThick-LiSparse-Reciprocal，RTLSR），根據(jù)體散射核、幾何光學(xué)核和各向同性核的加權(quán)和來(lái)估算地表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)過(guò)窄寬轉(zhuǎn)換計(jì)算得到積雪黑空反照率[12]。GlobAlbedo產(chǎn)品的反演模型為半經(jīng)驗(yàn)核驅(qū)動(dòng)模型RTLSR，經(jīng)過(guò)計(jì)算得到地表反照率[13-14]。GLASS產(chǎn)品采用直接估算方法，通過(guò)訓(xùn)練建立適用于全球不同區(qū)域、不同覆蓋類型、不同季節(jié)的格網(wǎng)查找表，直接估算得到地表反照率[15-17]。

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

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

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

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

（3）MCD12Q1土地覆蓋產(chǎn)品提供每年的全球土地覆蓋類型[21]。它通過(guò)對(duì)MODIS地表反射率數(shù)據(jù)進(jìn)行監(jiān)督分類得到初始分類結(jié)果，利用先驗(yàn)知識(shí)和輔助信息進(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)號(hào)航天飛機(jī)，采用干涉合成孔徑雷達(dá)技術(shù)生成了全球范圍內(nèi)的數(shù)字高程模型。本文使用經(jīng)過(guò)數(shù)據(jù)空白填補(bǔ)的最新Version 4版本，分辨率為90 m。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 驗(yàn)證方法

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

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

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

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

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

圖4 五種反照率產(chǎn)品的積雪反照率與實(shí)測(cè)值對(duì)比圖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)品采用專門針對(duì)積雪弱吸收特性建立的ART模型，能夠更好地表達(dá)積雪的各向異性反射特性；其他地表反照率產(chǎn)品主要采用針對(duì)植被建立的輻射傳輸模型，在積雪區(qū)域的精度較低；積雪反照率產(chǎn)品SAD需要較多的輸入?yún)?shù)，限制了模型的適用性；（2）利用我國(guó)森林生態(tài)監(jiān)測(cè)站獲取校正系數(shù)，顯著提升了森林區(qū)域的積雪反照率反演精度；（3）通過(guò)地形校正處理，進(jìn)一步提升了起伏地表的積雪反照率反演精度。

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

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

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

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

致 謝

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