王京達(dá)， 郝曉華， 和棟材， 王 建，3， 李弘毅， 趙 琴

（1.太原理工大學(xué)礦業(yè)工程學(xué)院，山西太原 030024； 2.中國科學(xué)院西北生態(tài)環(huán)境資源研究院，甘肅蘭州 730000；3.江蘇省地理信息資源開發(fā)與利用協(xié)同創(chuàng)新中心，江蘇南京 210023）

0 引言

積雪作為冰凍圈的重要組成部分，在地球表面分布較廣，大約98%的季節(jié)性積雪位于北半球，面積最大可達(dá)4.7×107km2[1］，且變化較為活躍，是影響氣候變化的重要因子[2-4］。在全球水循環(huán)中，積雪的積累和消融過程調(diào)節(jié)著水的年內(nèi)再分配，是干旱、半干旱地區(qū)春季最重要的淡水資源。積雪是氣候變化的指示器[5-6］，積雪范圍、積雪深度及雪水當(dāng)量等信息是水文模型和氣候模式的重要輸入?yún)?shù)[7-8］，全球尺度下的積雪參數(shù)監(jiān)測對于氣象、水循環(huán)、氣候變化以及能量平衡領(lǐng)域的研究具有重要的意義[9］。因此，準(zhǔn)確獲取北半球的積雪信息可為寒區(qū)水文循環(huán)、水資源管理、雪災(zāi)預(yù)警等重大研究提供必要可靠的參考資料[10］。

北半球積雪范圍和氣候變化研究需要長時間序列、高時空分辨率、高精度的積雪范圍數(shù)據(jù)。地面觀測積雪范圍資料是目前研究積雪長期變化特征最廣泛、最可靠的資料，具有時間尺度長、連續(xù)性好的優(yōu)點(diǎn)，但地面觀測資料無法準(zhǔn)確獲取積雪范圍[11-12］?？捎糜诒卑肭蜓芯康姆e雪范圍產(chǎn)品有微波產(chǎn)品和光學(xué)產(chǎn)品，微波積雪范圍產(chǎn)品主要來源于SMMR、SSM/I[13-14］和AMSR-E[15-16］傳感器。微波系列產(chǎn)品時間序列較長，但空間分辨率較低，目前可用的長時間序列微波積雪范圍產(chǎn)品的空間分辨率為25 km。 光學(xué)積雪范圍產(chǎn)品主要來源于AVHRR[17］傳感器、MODIS[18-19］傳感器，以及Landsat系列衛(wèi)星[20］和哨兵（Sentinel）系列衛(wèi)星[21］。MODIS積雪范圍產(chǎn)品（2000 年至今）時空分辨率較高，精度也較高，廣泛用于全球的水文和氣候變化研究，但對于氣候變化研究來說，時間序列僅有約20 年，較短的時間序列無法滿足研究；Landsat 時間序列長，空間分辨率高，但時間分辨率低（16 d），無法準(zhǔn)確監(jiān)測季節(jié)性積雪的變化；Sentinel系列衛(wèi)星雖然空間分辨率高，但時間序列短，時間分辨率較低，也無法滿足氣候變化研究需求。AVHRR 產(chǎn)品具有時間序列長（1981 年至今）、時空分辨率較高的特點(diǎn)（時間分辨率為1 d，空間分辨率為5 km），是研究北半球氣候變化的理想數(shù)據(jù)。

一些學(xué)者利用AVHRR 發(fā)展了積雪范圍產(chǎn)品。Hüsler 等[22］利用1 km 分辨率的AVHRR 原始數(shù)據(jù)集，開發(fā)了一種改進(jìn)的積雪范圍檢索算法來識別歐洲阿爾卑斯山的積雪范圍；Zhou 等[23］利用決策樹算法基于1 km 分辨率的AVHRR LAC 數(shù)據(jù)和HRPT數(shù)據(jù)，以中亞地區(qū)為例開發(fā)了一種長時間序列積雪范圍的識別算法。但此類研究主要集中在小區(qū)域的積雪識別，算法不具有全球普適性。日本宇航局（Japan Aerospace Exploration Agency，JAXA）利用AVHRR 和MODIS 結(jié)合發(fā)展的北半球1978—2015年5 km 分辨率的JASMES（JAXA Satellite Monitor?ing for Environmental Studies）積雪范圍產(chǎn)品，彌補(bǔ)了大尺度長時間序列逐日積雪范圍產(chǎn)品數(shù)據(jù)的缺失，但是該產(chǎn)品經(jīng)驗證發(fā)現(xiàn)在青藏高原精度較低[24］。Hao等[25］在中國區(qū)域發(fā)展了多級決策樹積雪識別算法，并生成了一套長時間序列積雪范圍產(chǎn)品，但算法中的閾值對于全球產(chǎn)品是否適用仍需要探索。

本研究的目的是利用AVHRR 數(shù)據(jù)，獲取多級決策樹積雪識別算法中適合北半球的閾值，生成1981—2000 年北半球歷史積雪范圍數(shù)據(jù)集，從而為氣候變化研究提供更長時間序列、更高精度的逐日積雪范圍數(shù)據(jù)集。研究主要解決兩方面問題：①如何通過有限的6個波段，結(jié)合研究區(qū)給定合理閾值；②如何填補(bǔ)AVHRR 在高緯度地區(qū)數(shù)據(jù)完全缺失和低緯度地區(qū)數(shù)據(jù)部分缺失問題，即如何填補(bǔ)數(shù)據(jù)空缺值。主要研究內(nèi)容包括：①以高空間分辨率的無云Landsat-5 TM 參考積雪范圍影像作為參考真值，獲取北半球積雪識別算法合理閾值；②填充積雪范圍數(shù)據(jù)中的空缺值，制作北半球AVHRR 逐日積雪范圍數(shù)據(jù)集；③利用北半球氣象臺站的雪深數(shù)據(jù)和Landsat-5 TM 作為參考真值，驗證數(shù)據(jù)集精度。該數(shù)據(jù)集將為北半球積雪變化研究提供重要支撐。

1 數(shù)據(jù)來源

1.1 AVHRR地表反射率數(shù)據(jù)

NOAA-AVHRR 地表反射率數(shù)據(jù)第四版（Sur?face Reflectance data Version 4，SR V4）被用作本研究的主要數(shù)據(jù)，來源于美國國家環(huán)境信息中心（https：//www. ncdc. noaa. gov/），它 是 由NOAA GAC（Global Area Coverage）L1 數(shù)據(jù)經(jīng)過一系列的定標(biāo)、輻射校正、大氣校正及幾何校正后得到的地表反射率數(shù)據(jù)和亮溫數(shù)據(jù)。該數(shù)據(jù)時間范圍為1981年6月24日至2019年5月16日，時間分辨率為1 d，空間分辨率為5 km。本研究使用1981 年6 月24日至2000年2月27日時間范圍內(nèi)的數(shù)據(jù)，數(shù)據(jù)被轉(zhuǎn)換投影坐標(biāo)為橢球體WGS84 Albers 投影坐標(biāo)，包含各個波段的地表反射率、亮溫和質(zhì)量控制標(biāo)記（quality control bit flags assessment，QA）等，具體如表1所示。

表1 NOAA-AVHRR地表反射率數(shù)據(jù)Table 1 NOAA-AVHRR surface reflectivity data

1.2 Landsat-5 TM數(shù)據(jù)

Landsat-5 TM 專 題 成 像 儀（thematic mapper，TM）是搭載在Landsat-5陸地觀測衛(wèi)星的傳感器，它于1984 年3 月發(fā)射成功直到2013 年6 月退役。本研究通過SNOWMAP 算法[18，26］獲取Landsat-5 TM二值積雪范圍，即雪和非雪。SNOWMAP 算法的核心是通過歸一化雪被指數(shù)（normalized difference snow index，NDSI）區(qū)分積雪和非雪，計算公式為

式中：NDSIL為Landsat-5 TM 的歸一化雪被指數(shù)；Band2 為Landsat-5 TM 影像的第二波段（0.52~0.60 μm）；Band5 為第五波段（1.55~1.75 μm）。算法設(shè)定NDSIL≥0.4 來區(qū)分雪和非雪，同時為了排除水體的干擾，設(shè)定第四波段（0.76~0.90 μm）的閾值大于0.11。利用該算法獲取Landsat-5 TM 空間分辨率為30 m 的二值積雪范圍后，通過最近鄰法將空間分辨率30 m 的Landsat-5 TM 二值積雪范圍影像重采樣為5 km 分辨率，轉(zhuǎn)換過程中設(shè)定5 km 像元中30 m 積雪像元占比≥50%作為積雪像元，否則為非雪像元，投影與AVHRR 保持一致，最終將獲取的5 km 分辨率Landsat-5 TM 二值積雪范圍影像作為參考真值。

Landsat-5 TM 參考積雪范圍影像分為兩組：第一組參考真值用作訓(xùn)練AVHRR 積雪識別算法閾值，選取1986年、1987年、1996年、1997年和1998年積雪季（11 月1 日開始至次年3 月31 日）的無云Landsat-5 TM 作為參考積雪范圍影像，共1 094 景，訓(xùn)練數(shù)據(jù)共包含積雪像元數(shù)為156 661 個，非雪像元數(shù)為696 296 個；第二組參考真值用于評估AVHRR 積雪范圍數(shù)據(jù)集，篩選1988 年、1989 年和1990 年積雪季Landsat-5 TM 無云且積雪占比大于50%的影像作為驗證數(shù)據(jù)集，共939景，驗證數(shù)據(jù)分布如圖1所示。

圖1 北半球氣象臺站和Landsat-5 TM驗證數(shù)據(jù)集的空間分布Fig. 1 Spatial distributions of meteorological stations and Landsat-5 TM validation data in the Northern Hemisphere

1.3 全球長時間序列逐日雪深數(shù)據(jù)集

本研究利用微波雪深數(shù)據(jù)集來填補(bǔ)AVHRR 中的空缺值，使用Che 等[27］采用SMMR 和SSM/I2 生成的全球長時間序列逐日雪深數(shù)據(jù)集。該微波雪深數(shù)據(jù)集，采用動態(tài)亮溫梯度算法，通過建立不同頻率亮溫差與實測雪深在空間和季節(jié)上的動態(tài)關(guān)系[28-29］獲取雪深。該數(shù)據(jù)集來源于時空三極環(huán)境大數(shù)據(jù)平臺（http：//poles. tpdc. ac. cn/zh-hans/）。本研究使用的時間范圍為1981 年6 月24 日至2000 年2月27日，數(shù)據(jù)集空間分辨率為25 km，利用最近鄰法將該數(shù)據(jù)集重采樣為5 km，并轉(zhuǎn)換為與AVHRR 相一致的投影。

1.4 地面雪深觀測數(shù)據(jù)

雪深數(shù)據(jù)用于對AVHRR 積雪識別算法和數(shù)據(jù)進(jìn)行評估，本研究使用的北半球雪深數(shù)據(jù)來自于NOAA 的國家環(huán)境信息中心網(wǎng)（http：//www. ncdc.noaa.gov/），是由世界氣象組織（World Meteorologi?cal Organization，WMO）提供的氣象臺站數(shù)據(jù)。本研究使用1981—1999 年2 546 個氣象臺站積雪季的實測雪深數(shù)據(jù)，使用的氣象臺站分布如圖1所示。

1.5 輔助數(shù)據(jù)

本研究共使用兩種輔助數(shù)據(jù)，一種是數(shù)字高程模型，另一種是氣候再分析數(shù)據(jù)-地表溫度產(chǎn)品。數(shù)字高程模型數(shù)據(jù)（Digital Elevation Model，DEM）為SRTM（Shuttle Radar Topography Mission）90 m數(shù)字高程數(shù)據(jù)庫第4 版本數(shù)據(jù)，是由美國國家航天局聯(lián)合美國國防部國家地理情報局通過航天飛機(jī)雷達(dá)測繪的數(shù)據(jù)產(chǎn)品（http：//srtm. csi. cgiar. org/），本研究中將STRM 數(shù)據(jù)用作云雪區(qū)分和積雪識別算法輔助數(shù)據(jù)，也用來作為消除積雪偽像元的輔助數(shù)據(jù)。

ERA5 氣候再分析數(shù)據(jù)-地表溫度逐日產(chǎn)品[30］是由歐盟資助歐洲中距離天氣預(yù)報中心（European Center for Medium-Range Weather Forecast，EC?MWF）執(zhí)行哥白尼氣候變化服務(wù)推出的適用于氣候研究的數(shù)據(jù)集。本研究中將ERA5數(shù)據(jù)集也作為輔助數(shù)據(jù)用來消除可能的積雪偽像元，空間分辨率為10 km，同樣采用最近鄰法將上述數(shù)據(jù)重采樣為5 km，并轉(zhuǎn)換為與AVHRR相一致的投影。

2 方法

本研究AVHRR 北半球積雪識別算法流程如圖2所示，主要分為以下四個步驟：

圖2 北半球AVHRR積雪識別算法流程Fig. 2 Flowchart of AVHRR snow discrimination algorithm in the Northern Hemisphere

（1）數(shù)據(jù)預(yù)處理

根據(jù)AVHRR SR V4 中自帶的質(zhì)量控制標(biāo)記剔除反射率無效值。

（2）云識別算法

針對冰云和積雪容易混淆的問題，結(jié)合云雪混淆區(qū)分技術(shù)，提高了積雪識別精度。

（3）積雪識別算法

利用大量的樣本進(jìn)行訓(xùn)練，針對北半球積雪研究區(qū)的特性，獲取多級決策樹積雪識別算法的最優(yōu)閾值，建立適合北半球的積雪識別算法，通過云、雪識別算法生成AVHRR L1級積雪范圍數(shù)據(jù)集。

（4）空缺值插補(bǔ)

針對AVHRR 高緯度地區(qū)數(shù)據(jù)完全缺失和低緯度地區(qū)部分?jǐn)?shù)據(jù)缺失問題，利用微波雪深數(shù)據(jù)集填充空缺數(shù)據(jù)。從而生成了一套北半球AVHRR 1981—1999 年逐日積雪范圍數(shù)據(jù)集，即AVHRR L2 級積雪范圍數(shù)據(jù)集。

2.1 數(shù)據(jù)預(yù)處理

根據(jù)AVHRR SR V4 中自帶的質(zhì)量控制標(biāo)記，如表2 所示，剔除AVHRR 各個通道的反射率無效值。有效值被作為下一步云、雪、非雪的區(qū)分，無效值被認(rèn)為是空缺值。

表2 AVHRR SR V4的質(zhì)量控制標(biāo)記描述Table 2 Descriptions of quality control bit flags assessment of AVHRR SR V4

2.2 云識別算法

關(guān)于云的識別算法，很多科學(xué)家也做了研究[10，31］，Hao 等[25］針對AVHRR SR V4 產(chǎn)品發(fā)展了一套云檢測算法，該算法可以有效解決云雪混淆問題。因此本研究結(jié)合Hao 等[25］生成的AVHRR 云雪混淆區(qū)分技術(shù)進(jìn)行云檢測，算法如表3所示。

表3 云檢測及其閾值Table 3 Cloud detection tests and their threshold values

2.3 積雪識別算法

由于Hao 等[25］開發(fā)的多級決策樹積雪識別算法是面對中國地區(qū)來提取積雪，在北半球算法的閾值還有待改進(jìn)。因此本研究采用該算法的方案（圖3），結(jié)合Landsat-5 TM 參考積雪范圍影像作為參考真值，來獲取北半球積雪識別算法最優(yōu)閾值。根據(jù)前人對AVHRR 積雪識別的分類研究[10，32］，涉及雪閾值算法的變量主要有Ref1、（Ref3/Ref2）和BT11三個一級閾值變量，NDVI 和（Ref3?Ref2）兩個二級閾值變量，NDSI 一個三級閾值變量，其中NDSI 的區(qū)分度最好。NDSI 通常使用綠色（大約0.50 μm波長）和短波紅外波段（大約1.60 μm 波長）計算。由于在AVHRR SR V4 中沒有1.60 μm 附近的短波紅外觀測，參考Hori 等[10］使用3.75 μm 波長的Ref3用 于NDSI 的 計 算，NDSI、NDVI 可 由 下 式 計 算得到。

（1）一級指標(biāo)閾值

如圖3 所示，一級指標(biāo)為Ref1、（Ref3/Ref2）、DEM 與BT11。一級指標(biāo)確立原則是排除大量的非雪數(shù)據(jù)，使盡可能多的雪被識別進(jìn)而再用下一級指標(biāo)進(jìn)行區(qū)分。一級指標(biāo)由于僅僅是初級篩選，因此一級指標(biāo)閾值采用Hao 等[25］算法的閾值。表4 展示了一級指標(biāo)和對應(yīng)的閾值，即當(dāng)Ref1＞0.14，（Ref3/Ref2）＜0.32、DEM＜1 300 m 并且BT11＜280 K 或者DEM≥1 300 m 并且BT11＜281.5 K 時，則認(rèn)為可能為積雪，否則為非雪，進(jìn)行下一級別區(qū)分。

表4 北半球積雪識別算法最優(yōu)閾值選取Table 4 Selection of optimal threshold for snow discriminant algorithm in the Northern Hemisphere

圖3 多級決策樹積雪識別算法流程Fig. 3 Flowchart of snow discriminant algorithm based on multi-level decision tree

（2）二級指標(biāo)閾值

如圖3所示，二級指標(biāo)是NDVI和（Ref3?Ref2）。將一級指標(biāo)篩選后可能為積雪的數(shù)據(jù)進(jìn)行分段統(tǒng)計，將非雪樣本1%各指標(biāo)對應(yīng)的反射率值作為二級指標(biāo)閾值，通過計算：獲取NDVI最優(yōu)閾值-0.11，即滿足NDVI＜-0.11 設(shè)定為第一類積雪（雪1）；獲取（Ref3?Ref2）最優(yōu)閾值-0.96，即滿足（Ref3?Ref2）＜-0.96 設(shè)定為第二類積雪（雪2），不符合條件的被認(rèn)為可能為積雪，進(jìn)行下一級別區(qū)分。

（3）三級指標(biāo)閾值

如圖3 所示，三級指標(biāo)是NDSI。利用二級指標(biāo)篩選后可能為積雪的數(shù)據(jù)，通過計算最佳總體精度的方法得出三級指標(biāo)，即逐步長方式計算最佳總體精度（步長設(shè)為0.01），最佳總體精度對應(yīng)的NDSI值即為三級指標(biāo)。最優(yōu)閾值的確定方法如圖4 所示。圖4（a）展示了積雪和非雪的NDSI 頻率分布，選擇具有最高總體精度OA（88%）時對應(yīng)的NDSI值作為最優(yōu)閾值（0.76）。計算結(jié)果顯示NDSI＞0.76，符合條件的設(shè)定為第三類積雪（雪3），不符合條件的被判別為非雪。北半球積雪識別算法的最優(yōu)閾值如表4所示。

圖4 三級指標(biāo)NDSI積雪和非雪頻率分布（a）和總體精度（b）Fig. 4 Frequency distribution of third-level indicator NDSI with and without snow（a）and overall accuracy（b）

2.4 空缺值插補(bǔ)

微波數(shù)據(jù)具有不受云干擾的特點(diǎn)，可對AVHRR 的空缺值進(jìn)行插補(bǔ)。利用Che 等[27］、Dai等[28］生成的全球被動微波數(shù)據(jù)逐日雪深數(shù)據(jù)集進(jìn)行插補(bǔ)。插補(bǔ)分為兩個步驟：

（1）將25 km 空間分辨率的雪深數(shù)據(jù)利用重采樣方法（最近鄰法）降尺度到5 km，投影和AVHRR一致。

（2）根據(jù)Che 等[27］研究結(jié)果，插補(bǔ)規(guī)則如下：①緯度處于0°~60° N 時，若AVHRR 的云像元對應(yīng)當(dāng)日雪深≥2 cm，則該云像元被分類為積雪像元；若雪深＜2 cm，則該云像元被分類為非雪像元。②緯度處于60°~90°N時，若所有像元對應(yīng)當(dāng)日雪深≥2 cm，則該像元被分類為積雪像元；若雪深＜2 cm，則該像元被分類為非雪像元。

3 結(jié)果與精度評估

3.1 結(jié)果

AVHRR L1 級和L2 級逐日積雪范圍數(shù)據(jù)集示意圖如圖5 所示。圖5（a）展示了通過云、雪識別后的北半球AVHRR L1 級逐日積雪范圍數(shù)據(jù)集示意圖，其中雪像元賦值為1，云像元賦值為255，陸地賦值為0。圖5（b）是北半球AVHRR L2級逐日積雪范圍數(shù)據(jù)集示意圖，其中雪像元賦值為1，云像元賦值為250，陸地賦值為0。

圖5 L1級（a）和L2級（b）數(shù)據(jù)集示意圖（選取日期為1998-01-13）Fig. 5 Schematic diagram of datasets at Level 1（a）and Level 2（b）（The selected date is 1998-01-13）

3.2 精度評估方法

3.2.1 基于氣象臺站雪深

利用1981—1999 年間北半球2 546 個氣象臺站的雪深數(shù)據(jù)分別對AVHRR L1 級和L2 級積雪范圍數(shù)據(jù)集進(jìn)行精度評估。設(shè)定氣象臺站雪深大于等于1 cm 時認(rèn)為該氣象臺站所在的像元為積雪像元，雪深小于1 cm 時認(rèn)為該氣象臺站所在的像元為非雪像元[33］。將氣象臺站積雪信息與AVHRR 各級積雪范圍數(shù)據(jù)集相應(yīng)位置的積雪信息進(jìn)行比對驗證。利用氣象臺站雪深數(shù)據(jù)轉(zhuǎn)化為雪和非雪的二值數(shù)據(jù)，并作為參考真值，然后分別構(gòu)建混淆矩陣分析總體精度及誤差。具體的SS、SN、NS 和NN 的定義參見表5，其中SS 指AVHRR 積雪范圍數(shù)據(jù)集和參考真值兩者都判別為雪的像元個數(shù)，NS 指AVHRR積雪范圍數(shù)據(jù)集判別為有雪且參考真值判別為非雪的像元個數(shù)，SN 指AVHRR 積雪范圍數(shù)據(jù)集判別為非雪且參考真值判別為雪的像元個數(shù)，NN 指兩者都判別為非雪的像元個數(shù)。

表5 精度驗證混淆矩陣的指標(biāo)含義Table 5 Index implication of confusion matrix for accuracy validation

分析數(shù)據(jù)集的總體精度（OA）、用戶精度（UA）、制圖精度（PA）和F-Score（FS）[34］。其中，OA 為雪和非雪像元被正確分類的總體精度；UA 為氣象臺站雪深數(shù)據(jù)被分類為積雪的概率，作為用戶精度；PA為AVHRR 積雪范圍數(shù)據(jù)集被分類為積雪的概率，作為制圖精度；FS 為UA 和PA 的調(diào)和平均值，是綜合考慮了UA 和PA 的結(jié)果，當(dāng)FS較高時說明UA 和PA較為均衡。

計算公式為

3.2.2 基于Landsat-5 TM

基于Landsat-5 TM 的驗證同樣使用3.2.1 節(jié)中構(gòu)建混淆矩陣的方法，進(jìn)一步評估AVHRR 積雪范圍數(shù)據(jù)集（L2 級）的精度。將AVHRR 積雪范圍數(shù)據(jù)集分別在北半球和中國、北美洲、歐亞大陸、西伯利亞、北極地區(qū)這五個北半球重要積雪區(qū)[35］進(jìn)行精度評估，分析數(shù)據(jù)集的OA、UA、PA 和FS。（Landsat-5 TM參考積雪范圍被當(dāng)作參考真值）

3.3 驗證結(jié)果與分析

3.3.1 基于氣象臺站雪深的精度評估

利用北半球1981—1999 年的共2 546 個實測雪深氣象臺站，分別對AVHRR L1 級和L2 級積雪范圍數(shù)據(jù)集進(jìn)行精度評估，計算結(jié)果如表6 所示，AVHRR L1 級數(shù)據(jù)集的總體精度為81.8%，AVHRR L2 級數(shù)據(jù)集的總體精度為82.2%，總體驗證精度較高。如圖6 所示，總體精度年際變化比較穩(wěn)定，都在80%左右，其中1996 年精度最高達(dá)到85.9%。AVHRR L1 級數(shù)據(jù)集的UA 為83.7%，PA為81.7%，AVHRR L2級數(shù)據(jù)集的UA為83.8%，PA為84.2%，表明AVHRR L1 級數(shù)據(jù)集漏分誤差略大于錯分誤差，尤其是在1993 年，由于AVHRR 大面積的數(shù)據(jù)空缺，導(dǎo)致漏分誤差較大，而AVHRR L2級數(shù)據(jù)集，制圖精度有所提高，表明雪深插值在保證精度的情況下，不僅可以填補(bǔ)空缺數(shù)據(jù)，而且可以有效降低漏分誤差，預(yù)測正確的積雪像元個數(shù)在AVHRR 識別的積雪像元中和氣象臺站記錄的積雪信息中的比例較高。FS 是平衡UA 和PA 的指數(shù)，AVHRR L1級和L2級數(shù)據(jù)集的FS 分別為82.9%和85.1%，F(xiàn)S 較高，表明漏分誤差和錯分誤差都比較均衡，算法可有效解決積雪識別問題。

圖6 L1級（a）和L2級（b）數(shù)據(jù)集的OA、UA和PA年際變化Fig. 6 Interannual variations of OA，UA and PA of datasets at Level 1（a）and Level 2（b）

表6 氣象臺站驗證評估結(jié)果Table 6 Verification and evaluation results of meteorological stations

3.3.2 基于Landsat-5 TM的精度評估

在空間上，利用1.2 節(jié)中的篩選方法，獲取了939 景Landsat-5 TM 參考積雪范圍影像作為參考真值，其中中國155 景、北美洲382 景、歐亞大陸455景、西伯利亞174 景及北極地區(qū)63 景。使用3.2.2節(jié)中的驗證方法，對待檢驗的AVHRR 積雪范圍數(shù)據(jù)集（L2級）進(jìn)行評估。

表7 展示了北半球和五個北半球重要積雪區(qū)（中國、北美洲、歐亞大陸、西伯利亞及北極地區(qū)）的驗證評估結(jié)果。結(jié)果表明，AVHRR 積雪范圍數(shù)據(jù)集的總體驗證結(jié)果OA、UA 和PA 分別為90.3%、90.2%和99.1%，該北半球積雪范圍數(shù)據(jù)集的總體精度、制圖精度和用戶精度都較高，數(shù)據(jù)具有一定可靠性。其中PA 最高達(dá)到99.1%，說明漏分誤差較低，預(yù)測正確的積雪像元個數(shù)在AVHRR 識別的積雪像元中和Landsat-5 TM 識別的積雪像元中的比例較高。FS 為93.4%，是UA 和PA 的平衡，較高則說明漏分誤差和錯分誤差較為均衡，積雪識別算法比較有效。對中國、歐亞大陸和西伯利亞地區(qū)的總體驗證結(jié)果OA、UA、PA 和FS 都在90.0%以上，北美洲地區(qū)總體精度OA 為85.2%，總體精度在中緯度地區(qū)較高，而在北極地區(qū)，OA 為73.7%，UA 為73.7%，PA 為99.9%，說明錯分誤差大于漏分誤差，且漏分誤差為0.1%，這是由于AVHRR 在高緯度地區(qū)數(shù)據(jù)的完全缺失，使用2.4 節(jié)中的數(shù)據(jù)和方法對其進(jìn)行插補(bǔ)，該數(shù)據(jù)集是通過微波數(shù)據(jù)和氣象臺站數(shù)據(jù)融合的數(shù)據(jù)集，說明該數(shù)據(jù)集錯分誤差較高，漏分誤差較低。

表7 Landsat-5 TM驗證評估結(jié)果Table 7 Verification and evaluation results of Landsat-5 TM

5 景獲取的Landsat-5 TM 驗證樣本的詳細(xì)信息如表8 所示，OA、UA、PA 和FS 的平均值分別為93.4%、93.7%、99.2%和96.2%，數(shù)據(jù)總體驗證精度較高，具有一定可靠性。每景驗證樣本對應(yīng)的Landsat-5 TM 參考積雪范圍、AVHRR 積雪范圍，以及AVHRR 和TM 的兩者對比如圖7 所示。說明在大范圍覆蓋的積雪區(qū)域，AVHRR 積雪范圍與Land?sat-5 TM 參考積雪范圍是較吻合的，而在積雪邊界區(qū)，AVHRR 積雪識別算法有未能識別出的Landsat-5 TM 參考積雪范圍影像的積雪像元，也有更多識別出的積雪像元，這可能是因為AVHRR 積雪范圍數(shù)據(jù)集對較低積雪范圍比例的像元識別能力較不準(zhǔn)確?？傮w而言，在積雪范圍識別上算法存在一定的漏分和錯分情況，但算法在保持漏分誤差和錯分誤差較為均衡的前提下，總體精度較高，數(shù)據(jù)具有可靠性。

圖7 Landsat-5 TM參考積雪范圍影像與AVHRR積雪范圍數(shù)據(jù)集對比Fig. 7 Comparison between Landsat-5 TM snow map and AVHRR snow cover extent datasets：China（a），North America（b），Eurasia（c），Siberia（d）and Arctic region（e）

表8 Landsat-5 TM驗證樣本信息Table 8 Verification sample information of Landsat-5 TM

4 結(jié)論

本研究基于AVHRR 地表反射率數(shù)據(jù)，獲取了基于多級決策樹積雪識別算法的各指標(biāo)最優(yōu)閾值，并結(jié)合云雪混淆區(qū)分技術(shù)，生成了北半球AVHRR 1981—1999 年L1 級逐日積雪范圍數(shù)據(jù)集。利用微波雪深數(shù)據(jù)插值填補(bǔ)了AVHRR 的空缺數(shù)據(jù)，生成了北半球AVHRR 1981—1999 年L2 級逐日積雪范圍數(shù)據(jù)集。最后利用氣象臺站實測雪深數(shù)據(jù)對L1級和L2級數(shù)據(jù)集分別進(jìn)行精度驗證，利用高空間分辨率遙感影像對L2級數(shù)據(jù)集進(jìn)一步進(jìn)行精度評估。

通過氣象臺站雪深數(shù)據(jù)對L1 級積雪范圍數(shù)據(jù)集驗證，評估結(jié)果表明：數(shù)據(jù)集總體精度為81.8%，用戶精度為83.7%，生產(chǎn)者精度為81.7%，數(shù)據(jù)集總體精度較高，漏分誤差比錯分誤差略高。FS 為82.9%，表明漏分誤差和錯分誤差比較均衡，整體上，北半球積雪識別算法精度較高。

通過氣象臺站雪深數(shù)據(jù)對L2 級積雪范圍數(shù)據(jù)集驗證，評估結(jié)果表明：數(shù)據(jù)集總體精度為82.2%，用戶精度為83.8%，生產(chǎn)者精度為84.2%，F(xiàn)S 為85.1%，L2 級數(shù)據(jù)集總體精度較高，漏分誤差和錯分誤差相近。同時結(jié)合Landsat-5 TM 參考積雪范圍影像，在面尺度上對L2 級數(shù)據(jù)集進(jìn)行進(jìn)一步驗證，結(jié)果表明：L2 級數(shù)據(jù)集總體精度為90.3%，用戶精度為90.2%，生產(chǎn)者精度為99.1%，在面尺度上，AVHRR L2 級積雪范圍數(shù)據(jù)集展現(xiàn)了更高的精度，表明插值方法可以有效解決AVHRR 在高緯度地區(qū)數(shù)據(jù)的完全缺失和低緯度地區(qū)數(shù)據(jù)的部分缺失問題。該算法將為北半球積雪范圍提取提供新的方案，L1 級和L2 級數(shù)據(jù)集也為北半球氣候變化研究提供了重要數(shù)據(jù)支撐。