1991—2021年天山融雪末期雪線高度遙感監(jiān)測研究

王靖文， 唐志光， 鄧 剛， 胡國杰， 桑國慶

（1.湖南科技大學(xué)測繪遙感信息工程湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 湘潭 411201；2.中國科學(xué)院西北生態(tài)環(huán)境資源研究院冰凍圈科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，甘肅 蘭州 730000）

雪線作為地球陸地表面積雪覆蓋區(qū)與非積雪覆蓋區(qū)之間的分界線［1］，其高度變化信息是冰雪水資源變化的直觀反映。在中緯度地區(qū)，季節(jié)變化能引起雪線位置的升降，這種臨時(shí)界線也稱為季節(jié)雪線（或瞬時(shí)雪線）［2-4］。而在冰川學(xué)中雪線特指消融末期積雪存在的下限［5］，等同于融雪末期雪線［6-7］。融雪末期雪線相比季節(jié)雪線更具有氣候意義，常用于指示冰川物質(zhì)平衡［8-11］，直接反映著冰川的進(jìn)退；研究融雪末期雪線高度變化有助于預(yù)測冰川和非季節(jié)性積雪未來的變化趨勢、理解區(qū)域和全球氣候變化［12］。因此，準(zhǔn)確獲取融雪末期雪線高度信息并實(shí)現(xiàn)其動(dòng)態(tài)監(jiān)測，不僅可為研究冰凍圈變化及其對氣候變化的響應(yīng)提供關(guān)鍵信息，對寒區(qū)水資源管理及其可持續(xù)利用也具有重要意義。

遙感技術(shù)在冰凍圈要素的時(shí)空變化研究中發(fā)揮著重要作用，已成為冰凍圈研究的重要技術(shù)手段［13-14］。MODIS 遙感數(shù)據(jù)因高時(shí)間分辨率的特點(diǎn)，在大范圍季節(jié)（或瞬時(shí)）雪線高度監(jiān)測中應(yīng)用廣泛［10,15-19］。Tang等［20］發(fā)展了基于MODIS積雪產(chǎn)品的融雪末期雪線高度遙感提取模型，提高了對亞洲高山區(qū)雪線高度的時(shí)空分布規(guī)律的宏觀認(rèn)識(shí)；但MODIS 數(shù)據(jù)因受空間分辨率的限制，僅適用于大尺度融雪末期雪線高度的宏觀遙感監(jiān)測。Landsat 系列衛(wèi)星提供了長時(shí)序和較高分辨率的遙感數(shù)據(jù)（近年發(fā)展的Sentinel-2衛(wèi)星時(shí)空分辨率更高），然而傳統(tǒng)計(jì)算機(jī)難以滿足處理海量遙感數(shù)據(jù)所需的存儲(chǔ)與計(jì)算資源，使得這些高分辨率的遙感大數(shù)據(jù)尚未被充分挖掘和應(yīng)用?，F(xiàn)有基于Landsat 系列遙感數(shù)據(jù)的雪線高度監(jiān)測研究僅限于小范圍局地區(qū)域［12,21-24］。近年來，遙感云計(jì)算平臺(tái)Google Earth Engine(GEE)的出現(xiàn)為遙感大數(shù)據(jù)挖掘提供了條件，并改變了傳統(tǒng)遙感數(shù)據(jù)處理和分析的模式，使得全球范圍長時(shí)間序列的分析和應(yīng)用成為可能［25］；GEE強(qiáng)大的存儲(chǔ)與云計(jì)算能力在積雪遙感領(lǐng)域也呈現(xiàn)出特有的優(yōu)勢［15,26-28］。此外，如何降低光學(xué)遙感中云遮蔽的影響，一直是積雪遙感監(jiān)測中面臨的重要問題，也是有效獲取雪線高度信息的關(guān)鍵。因此，綜合高分辨率遙感大數(shù)據(jù)的優(yōu)勢，降低云遮蔽的影響，發(fā)展遙感云計(jì)算平臺(tái)支持下的融雪末期雪線高度提取模型，是準(zhǔn)確監(jiān)測大范圍、長時(shí)序融雪末期雪線高度變化的重要途徑。

天山被譽(yù)為“中亞水塔”，橫亙于歐亞大陸腹地，由一系列高大山地、山間盆地和谷地組成，是世界上距離海洋最遠(yuǎn)、現(xiàn)代冰川發(fā)育最為廣泛的山系之一［29］。天山豐富的冰川、積雪是中亞干旱區(qū)的重要水資源，對氣候變化非常敏感。在過去的半個(gè)多世紀(jì)（1960—2016年），天山的年均氣溫呈顯著增加趨勢，增速為0.32 ℃·(10a)-1，比全球氣溫增加速率高出0.19 ℃·(10a)-1［30］。氣候變暖引起的天山冰雪變化直接影響河川徑流過程與水資源量的改變，導(dǎo)致中亞地區(qū)水資源的供需矛盾加劇，影響中亞區(qū)域國家之間關(guān)系以及“一帶一路”倡議［31］。為此，開展天山雪線高度變化的遙感監(jiān)測及其氣候響應(yīng)研究，可為研究氣候變化對天山山區(qū)的水資源影響提供重要信息，對科學(xué)管理水資源、服務(wù)于“一帶一路”倡議具有重要意義。

本文選取天山地區(qū)的4個(gè)山區(qū)流域作為研究對象，以Landsat 系列（TM、ETM+、OLI）衛(wèi)星數(shù)據(jù)為主要數(shù)據(jù)源，在遙感云計(jì)算平臺(tái)GEE 支持下，發(fā)展了融雪末期雪線高度提取模型，提取了4 個(gè)流域1991—2021年的融雪末期雪線高度數(shù)據(jù)，并分析了研究區(qū)融雪末期雪線高度的變化特征及其與氣候變化的關(guān)系。

1 研究區(qū)概況及數(shù)據(jù)

1.1 研究區(qū)概況

天山山脈橫跨中國、哈薩克斯坦、吉爾吉斯坦和烏茲別克斯坦等國家。天山地形復(fù)雜，高山區(qū)降水量大且氣溫低，冰川發(fā)育廣泛，其冰雪變化關(guān)系著諸多流域的水資源環(huán)境與水文過程。本文選擇瓊烏散庫什河流域、木扎提河流域、阿克牙孜河流域和瑪納斯河流域作為研究的典型流域（圖1），它們均發(fā)源于冰雪覆蓋區(qū)，冰雪融水在河川徑流中占有重要比例。瓊烏散庫什河和木扎提河發(fā)源于天山南坡，最終均匯入塔里木河流域，瓊烏散庫什河流域海拔4000 m 以上，主要是冰雪覆蓋，約占整個(gè)流域面積的35%，木扎提河流域內(nèi)發(fā)育有大量的現(xiàn)代冰川，約占整個(gè)流域面積的48%；阿克牙孜河流域約20%集中在海拔高程4000 m以上的高山帶；瑪納斯河發(fā)源于天山北坡，流入準(zhǔn)噶爾盆地腹地，流域內(nèi)海拔4000 m 以上常年被冰川和非季節(jié)性積雪覆蓋，約占整個(gè)流域面積的15%。

圖1 研究區(qū)概況圖Fig.1 Overview of the study area

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

研究采用了Landsat TM/ETM+/OLI、Sentinl-2 MSI、STRM DEM 和ERA5-LAND 再分析數(shù)據(jù)集，這些數(shù)據(jù)均通過GEE云計(jì)算平臺(tái)在線訪問和處理。

1.2.1 Landsat 數(shù)據(jù) 采用Landsat 系列數(shù)據(jù)中的陸表反射率（Land Surface Reflectance）數(shù)據(jù)，該數(shù)據(jù)已經(jīng)過輻射定標(biāo)和大氣校正，消除了大氣散射、吸收、反射引起的誤差，空間分辨率為30 m，時(shí)間分辨率為16 d。該數(shù)據(jù)提供了由CFMask 去云算法生成的質(zhì)量評估（QA）波段，包含了對云和云陰影的標(biāo)識(shí)。圖2 統(tǒng)計(jì)了各研究區(qū)逐年融雪期（6—9 月）共享的Landsat影像數(shù)量，以及研究區(qū)內(nèi)云量小于20%的影像頻次；利用Landsat 數(shù)據(jù)可以實(shí)現(xiàn)研究區(qū)1991 年至今的融雪末期雪線高度監(jiān)測。

圖2 研究區(qū)域融雪期覆蓋的影像數(shù)量Fig.2 Number of image covered by the study area at the melting season

1.2.2 Sentinel-2 MSI數(shù)據(jù) 歐洲航天局（ESA）發(fā)射的Senitnel-2 衛(wèi) 星 由Sentinel-2A 和Sentinel-2B 組成，分別于2015年、2017年發(fā)射，雙星重訪周期為5 d，均攜帶多光譜成像儀（MSI），共有13個(gè)通道，空間分辨率有10 m、20 m 和60 m。本文采用Sentinel-2 MSI 數(shù)據(jù)的Level-2A 產(chǎn)品對Landsat 數(shù)據(jù)的雪線高度提取結(jié)果進(jìn)行驗(yàn)證；該產(chǎn)品為經(jīng)過輻射定標(biāo)和大氣校正的地表反射率產(chǎn)品，選用數(shù)據(jù)波段的空間分辨率為10 m。

1.2.3 SRTM DEM數(shù)據(jù) 用于提取雪線高度信息的DEM 數(shù)據(jù)為SRTM DEM 數(shù)據(jù)（NASA 3.0 版）。該數(shù)據(jù)通過ASTER GDEM2、GMTED2010 和NED 數(shù)據(jù)對原始SRTM 數(shù)據(jù)的缺失值進(jìn)行填補(bǔ)，空間分辨率為30 m。該數(shù)據(jù)的全球垂直精度約為16 m，在中國地形平坦地區(qū)，絕對垂直精度約為14 m，地形起伏較大山區(qū)絕對垂直精度約為26 m［32］。

1.2.4 ERA5-LAND 再分析數(shù)據(jù)集 ERA5-LAND

再分析數(shù)據(jù)集來源于歐洲中期天氣預(yù)報(bào)中心（European Centre For Medium-Range Weather Forecasts, ECMWF），該數(shù)據(jù)集提供了1979 年1 月至今的全球降水、氣溫等變量數(shù)據(jù)。本文采用該數(shù)據(jù)集中的月平均氣溫和降水?dāng)?shù)據(jù)，空間分辨率為10 km，以探究研究區(qū)融雪末期雪線高度的變化與氣候因子之間的響應(yīng)關(guān)系。

2 研究方法

選取接近融雪末期的多時(shí)相Landsat TM/ETM+/OLI 影像為數(shù)據(jù)源，構(gòu)建雪線高度提取模型的技術(shù)路線如圖3所示。具體的方法主要包括積雪遙感識(shí)別、融雪末期“最小化”積雪范圍提取、區(qū)域雪線高度提取和精度驗(yàn)證。

圖3 總體技術(shù)路線Fig.3 Overall technical route

2.1 積雪遙感識(shí)別

采用SNOWMAP 算法［33］提取積雪范圍。該算法基于積雪在可見光波段，具有較高的反射率，在短紅外波段具有較強(qiáng)的吸收特征，建立歸一化積雪指數(shù)（Normalized Difference Snow Index,NDSI），設(shè)置合適的閾值，用以確定積雪覆蓋范圍。NDSI計(jì)算公式為：

式中：ρGreen和ρSWIR分別代表Landsat 數(shù)據(jù)中的綠光波段和短波紅外的反射率值。由于水體和積雪具有相似的波段反射特性，為了避免水體、陰影被誤識(shí)別為積雪，增加1 個(gè)判別條件即近紅外波段反射率值（ρNIR）＞0.11。相關(guān)研究證明了NDSI 閾值為0.29更適合于青藏高原地區(qū)的積雪制圖［34］。因此，本研究將NDSI＞0.29 &ρNIR＞0.11 的像元識(shí)別為積雪像元，不滿足上述條件的歸類為非雪像元，最終將識(shí)別得到的積雪、非雪與QA 波段識(shí)別得到的云合成得到積雪分類圖（雪、云、非雪）。

2.2 融雪末期“最小化”積雪范圍提取

通過瞬時(shí)雪線高度去云和時(shí)序融合，以減少云的影響并得到最能代表融雪末期的“最小化”積雪范圍。

2.2.1 瞬時(shí)雪線高度去云 對區(qū)域內(nèi)云和雪的比例進(jìn)行判斷，若云和雪面積比例≤1/2，則提取該區(qū)域的瞬時(shí)雪線高度以進(jìn)行去云處理。瞬時(shí)雪線去云策略：如果云像元對應(yīng)的地面海拔高度大于瞬時(shí)區(qū)域雪線高度，則歸類為積雪；反之，則歸類為陸地。因云覆蓋的比例過大，容易出現(xiàn)積雪像元極少或者積雪覆蓋區(qū)周圍均為云覆蓋，導(dǎo)致提取的瞬時(shí)雪線高度存在較大偏差或者無法提取。為了保證提取的瞬時(shí)區(qū)域雪線高度的有效性，本研究設(shè)定云和雪面積比例＞1/2時(shí)，則直接參與時(shí)序融合處理。

2.2.2 時(shí)序融合 對接近融雪末期的多時(shí)相積雪分類圖像（包含瞬時(shí)雪線高度去云后圖像）進(jìn)行時(shí)序融合，以完成積雪范圍的最小化合成，最終得到能代表融雪末期的“最小化”積雪范圍圖。時(shí)序融合策略：對每一像元，如果被任一時(shí)相識(shí)別為裸地，則歸類為裸地；如果50%以上的像元被識(shí)別為積雪，而其他時(shí)相都為云，則歸類為積雪；如果50%以上的像元被識(shí)別為云，而其他時(shí)相都為雪，則歸類為云。融合后得到的代表融雪末期的“最小化”積雪范圍圖，用于融雪末期雪線高度的提取。

2.3 區(qū)域雪線高度（RSLA）提取

區(qū)域雪線高度法的原理為：在目標(biāo)區(qū)域內(nèi)找出一個(gè)高程，使得該高程之上的非雪像元數(shù)量與該高程之下的積雪像元數(shù)量之和最小，那么這一高程就代表該區(qū)域內(nèi)的雪線高度（圖4），公式表示為：

圖4 區(qū)域雪線高度提取示意圖Fig.4 Sketch map of the extraction of regional snowline altitude

式中：sum為高程之下積雪像元數(shù)量（Cbelow）與高程之上非雪像元數(shù)量（Cabove）之和；RSLA 為區(qū)域雪線高度；elevmin為區(qū)域內(nèi)最低高程；elevmax為區(qū)域內(nèi)最高高程。

具體過程為：從區(qū)域最低高程開始，計(jì)算目標(biāo)高程之下的積雪像元數(shù)量與該高程之上的非雪像元數(shù)量之和，然后每隔1 m高程進(jìn)行1次該計(jì)算，直到區(qū)域內(nèi)最高高程。然后利用最小值函數(shù)（Min）來檢索數(shù)量和中的最小值sumRSLA，sumRSLA所對應(yīng)的高程即為目標(biāo)區(qū)域的雪線高度。該方法可以允許部分云覆蓋的存在，是實(shí)現(xiàn)目標(biāo)區(qū)域雪線高度的最優(yōu)估算，因而也適用于積雪圖像的去云處理。

2.4 雪線高度的精度評估

采用空間分辨率更高、重訪周期更短的Sentinel-2 影像，通過提取融雪末期“最小化”積雪范圍，來驗(yàn)證Landsat 數(shù)據(jù)提取融雪末期區(qū)域雪線高度的精度。具體驗(yàn)證方法：選取2019—2021年研究區(qū)接近融雪末期且云量低于10%的多時(shí)相Sentinel-2 遙感影像（表1），采用目視解譯法，對每一影像進(jìn)行積雪范圍制圖；然后將接近融雪末期的多時(shí)相積雪范圍圖進(jìn)行時(shí)序融合，得到Sentinel-2 提取的融雪末期“最小化”積雪范圍（作為融雪末期積雪范圍“真值”）。最后，基于Sentinel-2提取的融雪末期“最小化”積雪范圍圖，在雪線位置附近進(jìn)行驗(yàn)證樣本的隨機(jī)采樣，并結(jié)合DEM 數(shù)據(jù)和Landsat 提取的融雪末期區(qū)域雪線高度值，計(jì)算融雪末期雪線高度提取模型的總體精度（OA）、精確率（Pre）和召回率（Rec），計(jì)算公式如下：

表1 研究區(qū)域Sentinel-2影像覆蓋數(shù)量Tab.1 Number of Sentinel-2 image covered by the study area

式中：A為高于Landsat 融雪末期雪線高度且在Sentinel-2“最小化”積雪范圍中識(shí)別為積雪的點(diǎn)數(shù)；B為高于Landsat 融雪末期雪線高度且在Sentinel-2“最小化”積雪范圍中識(shí)別為非雪的點(diǎn)數(shù)；C為低于Landsat融雪末期雪線高度且在Sentinel-2“最小化”積雪范圍中識(shí)別為積雪的點(diǎn)數(shù)；D為低于Landsat融雪末期雪線高度且在Sentinel-2“最小化”積雪范圍中識(shí)別為非雪的點(diǎn)數(shù)。Pre＞Rec 時(shí)，代表雪線高度存在高估；Pre＜Rec時(shí)，代表雪線高度存在低估。

3 結(jié)果與分析

3.1 融雪末期雪線高度精度評估結(jié)果分析

通過建立混淆矩陣，計(jì)算研究區(qū)域融雪末期雪線高度提取結(jié)果的Pre、Rec、OA和Kappa系數(shù)，對提取結(jié)果的精度進(jìn)行定量評價(jià)（表2）。結(jié)果表明：4個(gè)流域融雪末期雪線高度提取結(jié)果具有較高的精度，OA 介于87.5%～94.4%之間，Kappa 系數(shù)介于0.86～0.93 之間；平均OA 和Kappa 系數(shù)分別為91.6%和0.9。然而，4 個(gè)流域Pre 均略低于Rec，說明相比時(shí)間分辨率更高的Sentinel-2融雪末期“最小化”積雪范圍，Landsat 提取的融雪末期雪線高度略有低估。以阿克牙孜河流域?yàn)槔?，OA 為92.7%，說明Landsat提取的融雪末期雪線高度與Sentinel-2提取的融雪末期“最小化”積雪范圍之間具有較高的一致性；精確率Pre為87.06%，說明在Landsat提取的融雪末期雪線高度之上且在Sentinel-2提取的“最小化”積雪范圍中被識(shí)別為非積雪的概率為12.94%；而該流域Rec為93.97%，說明在Landsat提取雪線高度之下且在Sentinel-2提取的“最小化”積雪范圍中被識(shí)別為積雪的概率僅為6.03%。

表2 Landsat融雪末期雪線高度提取結(jié)果精度評價(jià)Tab.2 Accuracy evaluation of snowline altitude extraction results at the end of melting season

3.2 融雪期雪線高度時(shí)空變化特征

獲取各流域6—9 月以16 d 為間隔的多年雪線高度數(shù)據(jù)（圖5）。雪線高度在6—9月呈現(xiàn)出強(qiáng)烈的季節(jié)性變化。從6月開始，隨著積雪的快速融化，研究區(qū)積雪覆蓋率逐漸降低、雪線高度逐漸升高，瞬時(shí)雪線高度的最大值主要出現(xiàn)在8 月，而8 月以后雪線高度又逐漸降低；且研究區(qū)瞬時(shí)雪線高度體現(xiàn)出較強(qiáng)的年際波動(dòng)（圖5）。以2020年阿克牙孜河流域的瞬時(shí)雪線高度變化為例（圖6），雪線高度在6—8 月逐漸升高，從6 月的3775 m 上升到8 月的3936 m，而9月雪線高度又下降到3458 m；經(jīng)模型提取的融雪末期雪線高度為3967 m，高于每一瞬時(shí)雪線高度。

圖5 6—9月雪線高度與平均積雪覆蓋率變化Fig.5 Variation of snowline altitude and snow cover percentage from June to September

圖6 瞬時(shí)雪線高度空間分布（以阿克牙孜河流域?yàn)槔〧ig.6 Spatial distribution of instantaneous snowline altitude(Take the Akyazi River basin for example)

3.3 融雪末期雪線高度年際變化特征

1991—2021 年研究區(qū)域融雪末期雪線高度總體上呈明顯的上升趨勢（圖7）。瓊烏散庫什河流域和瑪納斯河流域融雪末期雪線高度最高值均出現(xiàn)在2018年，最低值分別出現(xiàn)在2003年和1994年；木扎提河流域和阿克牙孜河流域融雪末期雪線高度最高值、最低值分別出現(xiàn)在2016年、2015年和1993年、1998年，4個(gè)流域近30 a融雪末期雪線高度振幅均不超過400 m。瓊烏散庫什河流域多年平均融雪末期雪線高度為4504 m，趨勢斜率為4.5 m·a-1，該流域近30 a 融雪末期雪線高度上升了約135 m；木扎提河流域多年平均融雪末期雪線高度為4207 m，趨勢斜率為3.1 m·a-1，近30 a 融雪末期雪線高度上升約93 m；瑪納斯河流域融雪末期雪線高度總體上與木扎提河流域接近，該流域多年平均融雪末期雪線高度為4212 m，趨勢斜率最高，為6.4 m·a-1，近30 a融雪末期雪線高度上升約192 m；阿克牙孜河流域多年平均融雪末期雪線高度為3955 m，趨勢斜率在4 個(gè)流域中最低，為2.8 m·a-1，近30 a 融雪末期雪線高度上升約84 m。此外，通過計(jì)算研究區(qū)融雪末期雪線高度的標(biāo)準(zhǔn)差（σ），來反映融雪末期雪線高度的年際波動(dòng)幅度（圖7）。結(jié)果表明：木扎提河流域融雪末期雪線高度標(biāo)準(zhǔn)差最?。é?59 m），說明該流域的融雪末期雪線高度年際波動(dòng)最小，瓊烏散庫什河流域融雪末期雪線高度標(biāo)準(zhǔn)差最大（σ=94 m），說明該流域?qū)?yīng)的融雪末期雪線高度年際波動(dòng)最大。

圖7 融雪末期雪線高度年際變化Fig.7 Interannual variation of snowline altitude at the end of melting season

3.4 融雪末期雪線高度變化與氣溫、降水之間的關(guān)系

雪線高度與氣候之間的關(guān)系非常密切。本文通過計(jì)算年平均和夏季平均氣溫、年總和夏季總降水量與融雪末期雪線高度之間的相關(guān)關(guān)系，分析融雪末期雪線高度與氣候變化之間的響應(yīng)關(guān)系（表3），其中，年降水量和年氣溫指的是由上一年10 月到當(dāng)年9月的總降水量和平均氣溫。瓊烏散庫什河流域年降水量和夏季降水量均呈微弱的增加趨勢，木扎提河流域、阿克牙孜河流域和瑪納斯河流域年降水量和夏季降水量均呈減少趨勢；各流域年均氣溫和夏季平均氣溫均呈上升趨勢（圖8）。其中，木扎提河流域多年年平均氣溫（-8.9 ℃）和夏季平均氣溫（1.9 ℃）最低，瑪納斯河流域多年年平均氣溫（-3.3 ℃）和夏季平均氣溫（8.7 ℃）最高；瓊烏散庫什河流域多年年平均降水量（485.8 mm）和夏季平均降水量（306.3 mm）最少，阿克牙孜河流域多年年平均降水量（930.6 mm）和夏季平均降水量（596.1 mm）最多。研究區(qū)融雪末期雪線高度與氣溫之間呈正相關(guān)關(guān)系，特別是與夏季氣溫之間呈顯著正相關(guān)關(guān)系（P＜0.05）；而融雪末期雪線高度與降水量之間呈較弱的負(fù)相關(guān)關(guān)系，這表明夏季氣溫是影響融雪末期雪線高度變化的主要因素。

圖8 融雪末期雪線高度與氣溫和降水量的年際變化Fig.8 Interannual variations of temperature and precipitation between snowline altitude at the end of melting season

表3 融雪末期雪線高度與氣溫、降水量之間的相關(guān)性Tab.3 Correlation between snowline altitude at the end of melting season,temperature and precipitation

為探究研究區(qū)融雪末期雪線高度對氣候變化的敏感性，采用多元線性回歸方法建立了融雪末期雪線高度與夏季氣溫和年降水量回歸模型：

式中：SLA為融雪末期雪線高度；λ為常數(shù)；T為夏季氣溫（℃）；P為年降水量（mm）；a、b為回歸系數(shù)。4個(gè)流域回歸系數(shù)如表4 所示，以阿克牙孜河流域?yàn)槔?dāng)夏季氣溫升高或降低1 ℃時(shí)，阿克牙孜河流域的融雪末期雪線高度將升高或降低76.2 m；當(dāng)年降水量增加或減少100 mm時(shí)，該流域的融雪末期雪線高度將降低或升高28.7 m。說明在阿克牙孜河流域，夏季氣溫升高1 ℃引起的融雪末期雪線高度上升需要增加265 mm 的降水量才能彌補(bǔ)。瓊烏散庫什河流域夏季氣溫回歸系數(shù)最大（82.9），意味著該流域融雪末期雪線高度變化對夏季氣溫響應(yīng)最為敏感?，敿{斯河流域年降水回歸系數(shù)最大（0.38），說明該流域融雪末期雪線高度變化對年降水量響應(yīng)最為敏感。

表4 融雪末期雪線高度回歸模型系數(shù)Tab.4 Regression model coefficients of snowline altitude at the end of melting season

4 討論

4.1 積雪制圖與雪線高度提取的不確定性

積雪制圖和雪線高度提取的不確定性主要來自以下4 個(gè)方面：（1）利用光學(xué)遙感影像提取山區(qū)積雪信息時(shí)，地形起伏產(chǎn)生的山體陰影會(huì)影響積雪范圍提取的準(zhǔn)確性，這種影響也被稱為地形效應(yīng)。通過地形校正［35-36］可以較好地消除地形效應(yīng)對山區(qū)積雪識(shí)別的影響。此外，如果是像元尺度雪線高度的提取，還應(yīng)定量分析坡向?qū)ρ┚€高度的影響。本研究專注的是區(qū)域性雪線高度的提取。該方法求解得到的唯一高程值代表了一個(gè)區(qū)域內(nèi)的綜合雪線高度。因而，本研究沒有考慮地形效應(yīng)以及坡向等微地形因子的影響。（2）云和云陰影的存在也是影響光學(xué)遙感中積雪信息提取的主要誤差來源。本文基于Landsat系列數(shù)據(jù)，采用瞬時(shí)雪線去云和時(shí)序融合的方法進(jìn)行去云處理。而Sentinel-1搭載了合成孔徑雷達(dá)（SAR），該數(shù)據(jù)具有全天時(shí)、全天候和穿透性的特點(diǎn)，不受天氣和太陽光照條件的影響；Sentinel-2數(shù)據(jù)為多光譜衛(wèi)星成像，該數(shù)據(jù)具有較高的時(shí)間分辨率（5 d）。在后續(xù)的進(jìn)一步研究中可以結(jié)合Sentinel-1 和Sentinel-2 數(shù)據(jù)，增強(qiáng)云下積雪識(shí)別的能力。（3）對于森林覆蓋區(qū)積雪的提取精度較差，樹木的冠層會(huì)擋住地面部分的積雪信息?；贜DSI、NDFSI（Normalized Difference Forest Snow Index）、NDVI(Normalized Difference Vegetation Index)的多指標(biāo)積雪識(shí)別方法［37-38］，可有效提取林下積雪信息。但融雪末期雪線高度要遠(yuǎn)高于森林植被區(qū)，所以本研究沒有考慮森林的影響。（4）因大部分裸露冰川表面不清潔（表磧覆蓋）反照率明顯低于積雪，因而不清潔的裸露冰川和非季節(jié)性積雪可以通過光學(xué)遙感（NDSI閾值）進(jìn)行區(qū)分。然而，融雪末期雪線高度的提取，也應(yīng)該注重區(qū)分非季節(jié)性積雪與表面清潔的冰川；例如可利用兩者在近紅外波段上所表現(xiàn)出的不同特征，根據(jù)近紅外波段影像直方圖分類出冰川和積雪［39-40］。特別是隨著Sentinel-1 等星載合成孔徑雷達(dá)的出現(xiàn)，使冰川和積雪的區(qū)分更具可能性。

本研究中的融雪末期雪線高度提取方法僅適合于區(qū)域（或者流域）尺度的雪線高度提取。但若注重分析像元尺度雪線高度的時(shí)空差異性（比如坡向、坡度等微地形因子的影響），應(yīng)采用像元尺度雪線高度提取方法［19,41-42］。

4.2 融雪末期雪線高度的空間分異

本研究發(fā)現(xiàn)，近30 a 來研究區(qū)域融雪末期雪線高度呈上升趨勢，這與天山物質(zhì)平衡處于負(fù)平衡狀態(tài)的研究結(jié)論是一致的［43］。從研究區(qū)域每年的融雪末期雪線高度來看，瓊烏散庫什河流域總體融雪末期雪線高度明顯高于其他3 個(gè)流域，阿克牙孜河流域融雪末期雪線高度處于最低水平。造成這種雪線高度分布特征的主要原因?yàn)椋涵倿跎焓埠恿饔蛭挥谔焐侥掀拢焐侥掀聻殛柶?，受到太陽輻射量較多，造成氣溫相對于北坡較高［44］；且該流域地勢總體比較高，受山體海拔影響，所以該地區(qū)融雪末期雪線高度最高；阿克牙孜河流域位于天山北坡，受大西洋和北冰洋的水汽影響［45］，降水量相對充沛，太陽輻射量少，氣溫相對較低，因而阿克牙孜河流域融雪末期雪線高度最低。

4.3 氣候變化對融雪末期雪線高度的影響

雪線作為一種氣候標(biāo)志線，其高度受氣候因素的影響同時(shí)也會(huì)反映氣候的變化?，F(xiàn)有的研究［42,46］表明：在全球氣候變暖的影響下冰川面積不斷萎縮，雪線高度持續(xù)升高，這與本文的研究相一致。本研究發(fā)現(xiàn)，夏季氣溫是影響融雪末期雪線高度變化的主要因素，以往有關(guān)高山區(qū)融雪末期雪線或者冰川平衡線變化的研究也得出類似結(jié)論：Tang等［20］基于MODIS 數(shù)據(jù)從30 km 格網(wǎng)尺度也發(fā)現(xiàn)了2001—2016 年亞洲高山區(qū)融雪末期雪線高度總體呈上升趨勢，特別是在天山地區(qū)和喜馬拉雅地區(qū)上升明顯，夏季氣溫是主要影響因素；陳安安等［47］、王寧練等［48］發(fā)現(xiàn)，夏季氣溫是影響冰川平衡線變化的主導(dǎo)因子。雪線高度與氣溫呈正相關(guān)關(guān)系，氣溫升高會(huì)促進(jìn)雪線高度的升高；雪線高度與降水量呈負(fù)相關(guān)關(guān)系，降水量增加會(huì)引起雪線高度的降低。在后續(xù)的研究中可以基于遙感提取的長時(shí)序雪線高度、氣溫和降水等氣象數(shù)據(jù)，利用機(jī)器學(xué)習(xí)方法構(gòu)建雪線高度變化的預(yù)測模型。

5 結(jié)論

本文基于GEE 云計(jì)算平臺(tái)和Landsat 影像，發(fā)展了融雪末期區(qū)域雪線高度遙感提取模型，提取了天山4個(gè)流域近30 a長時(shí)間序列的融雪末期雪線高度；并對4 個(gè)流域融雪末期雪線高度的時(shí)空變化特征及其與氣溫、降水之間的關(guān)系進(jìn)行分析，得出以下結(jié)論：

（1）模型提取的融雪末期雪線高度與Sentinel-2提取的融雪末期“最小化”積雪范圍具有較高的一致性，平均總體精度為91.6%；Kappa 系數(shù)平均值為0.9。利用該提取模型可準(zhǔn)確提取流域內(nèi)融雪末期區(qū)域雪線高度，可有效應(yīng)用于流域尺度、長時(shí)間序列雪線高度的遙感監(jiān)測。

（2）研究區(qū)1991—2021 年融雪末期雪線高度年際波動(dòng)較大，但總體上呈現(xiàn)出上升趨勢。阿克牙孜河流域融雪末期雪線高度上升速率最慢（2.8 m·a-1），瑪納斯河流域融雪末期雪線高度上升速率最快（6.4 m·a-1）。

（3）研究區(qū)融雪末期雪線高度主要受氣溫的影響，與夏季氣溫之間具有顯著的正相關(guān)關(guān)系，與降水量之間具有較弱的負(fù)相關(guān)關(guān)系。當(dāng)夏季氣溫升高（降低）1 ℃時(shí)，流域融雪末期雪線高度約升高（降低）50.7～82.9 m；當(dāng)年降水量增加（減少）100 mm 時(shí)，流域融雪末期雪線高度將約降低（升高）23.8～36.4 m。

（4）GEE 平臺(tái)強(qiáng)大的存儲(chǔ)與計(jì)算能力等，大幅提高了雪線高度遙感提取的效率，為長時(shí)序、大范圍遙感雪線高度監(jiān)測提供了強(qiáng)有力的計(jì)算資源保障。