周 苗，常曉濤，朱廣彬，瞿慶亮，劉 偉

1. 自然資源部國土衛(wèi)星遙感應(yīng)用中心，北京 100048； 2. 山東科技大學(xué)測繪與空間信息學(xué)院, 青島 266590

冰川是冰凍圈的主要組成部分，對氣候變化極為敏感，被認(rèn)為是氣候變化的最佳天然指示器之一[1]。冰川主要分布在極區(qū)和高海拔地區(qū)，隨著氣候的周期性變化處于一種動態(tài)平衡當(dāng)中。隨著全球變暖，大多數(shù)冰川呈加速消融趨勢[2-4]，增加了發(fā)生冰湖潰決、泥石流等冰川災(zāi)害的可能[5]。念青唐古拉山脈是我國第二大冰川聚集區(qū)，冰川消融量巨大，2000—2015年間念青唐古拉山脈冰川總面積減少了463.36 km2,消融量高居亞洲其他山脈之首。開展念青唐古拉山脈冰川變化監(jiān)測對區(qū)域生態(tài)保護(hù)、災(zāi)害預(yù)警等具有重要的現(xiàn)實(shí)意義。

隨著空間對地觀測技術(shù)的發(fā)展，冰川變化監(jiān)測技術(shù)更加多樣化，監(jiān)測效率更加高效。近年來，衛(wèi)星重力和光學(xué)遙感技術(shù)在冰川研究中廣泛應(yīng)用。衛(wèi)星重力可直接獲取區(qū)域質(zhì)量變化，在研究冰川變化機(jī)制和大、中尺度質(zhì)量變化方面具有一定優(yōu)勢。文獻(xiàn)[6]利用GRACE、GRACE-FO時變重力場模型解算的等效水高序列評估了2002年4月至2019年9月期間世界(不包括南極和格林蘭島)冰川和冰蓋的質(zhì)量平衡。結(jié)果表明，世界冰川和冰蓋的平均損失速率為281.5±30 Gt/a。文獻(xiàn)[7]基于GRACE衛(wèi)星重力資料指出，2003—2015年間青藏高原冰川質(zhì)量一直減少。文獻(xiàn)[8]利用GRACE衛(wèi)星重力數(shù)據(jù)在頻域上計(jì)算了2002—2014年喜馬拉雅山冰川質(zhì)量變化，發(fā)現(xiàn)該區(qū)域冰川質(zhì)量變化整體上呈現(xiàn)加速消融趨勢。光學(xué)衛(wèi)星遙感技術(shù)能很好地反映冰川變化的細(xì)節(jié)，是提取冰川面積和冰面儲量變化的最快速、最有效手段。文獻(xiàn)[9]利用Landsat TM/ETM+/OLI影像數(shù)據(jù)，結(jié)合比值閾值法與目視解譯提取冰川邊界，對1990—2015年喜馬拉雅山冰川的消融情況進(jìn)行了研究。分析認(rèn)為，近25年來喜馬拉雅山冰川整體呈退縮趨勢，冰川面積由23 229.27 km2減少到20 676.17 km2，退縮率為10.99%，并有加速退縮的趨勢。文獻(xiàn)[10]利用資源三號立體影像和SRTM DEM數(shù)據(jù)對念青唐古拉山脈2000—2017年的冰川物質(zhì)平衡變化進(jìn)行估計(jì)。研究發(fā)現(xiàn)，念青唐古拉山脈整體冰川物質(zhì)平衡為負(fù)值，其中念青唐古拉山脈西段2013—2017年間冰川儲量損失最快。文獻(xiàn)[11]基于地形圖和Landsat TM/OLI遙感影像，利用目視解譯和波段比值法研究了南迦巴瓦峰地區(qū)1980—2015年的冰川變化。研究發(fā)現(xiàn)南迦巴瓦峰地區(qū)冰川面積持續(xù)減小并呈加速退縮的趨勢。

國內(nèi)外學(xué)者在利用衛(wèi)星重力和遙感技術(shù)實(shí)現(xiàn)冰川質(zhì)量、面積變化方面進(jìn)行了較為廣泛的研究，為冰川質(zhì)量變化監(jiān)測做出了重要貢獻(xiàn)，但兩種技術(shù)協(xié)同應(yīng)用方面仍有一些不足，具有一定的研究空間。本文結(jié)合GRACE、GRACE-FO衛(wèi)星重力數(shù)據(jù)和遙感資料對念青唐古拉山脈冰川變化展開研究，并利用氣象數(shù)據(jù)對冰川的變化機(jī)制及影響因素進(jìn)行分析，為該區(qū)域冰川水資源變化監(jiān)測及生態(tài)保護(hù)提供參考。

1 研究方法

1.1 地表質(zhì)量變化反演

利用重力場模型的球諧系數(shù)變化計(jì)算地球表面質(zhì)量變化的模型為

(1)

本文采用美國得克薩斯大學(xué)空間研究中心提供的GRACE、GRACE-FO RL06版本時變重力場數(shù)據(jù)，模型的最大階數(shù)為96，時間跨度為2002年8月至2020年4月，共有33個月的數(shù)據(jù)缺失，其中2017年6月至2018年6月這11月的數(shù)據(jù)空白是由于衛(wèi)星任務(wù)更替造成的，對除了這11個月的其他缺失數(shù)據(jù)采用3次樣條插值方法補(bǔ)齊。由于GRACE、GRACE-FO衛(wèi)星對低階項(xiàng)不敏感，本文使用文獻(xiàn)[12]提供的一階項(xiàng)和衛(wèi)星激光測距(satellite laser ranging，SLR)獲得的C20項(xiàng)[13]代替球諧系數(shù)中的相應(yīng)項(xiàng)。另外，從2016年10月開始，GRACE和GRACE-FO衛(wèi)星上僅有一顆加速度計(jì)數(shù)據(jù)用于地球時變重力場的計(jì)算，這會造成該段時間GRACE、GRACE-FO得到的C30項(xiàng)精度較低，本文使用SLR獲得的C30對其進(jìn)行替換[14]。同時，為了削弱高階次噪聲和奇偶階相關(guān)誤差的影響[15]，對位系數(shù)進(jìn)行300 km高斯濾波和滑動窗去相關(guān)濾波[16]。最后，結(jié)合式(1)得到更精確的地表質(zhì)量變化反演公式

(ΔClmcos(mλ)+ΔSlmsin(mλ))

(2)

式中，w(l)為高斯平滑函數(shù)。

經(jīng)過以上處理得到的結(jié)果中仍包含冰后回彈(glacial isostatic adjustment,GIA)和泄露誤差的影響，本文利用ICE-6G-C模型[17]進(jìn)行GIA改正。同時，采用文獻(xiàn)[18]提出的空域法對泄露信號進(jìn)行恢復(fù)?？沼?yàn)V波法的原理是假定地表質(zhì)量變化是由一系列位置已知、質(zhì)量未知的質(zhì)量塊(Mascon)引起的，根據(jù)Mascon的球諧展開和綜合結(jié)果與GRACE結(jié)果的線性關(guān)系，可推求每個Mascon對應(yīng)的系數(shù)，此系數(shù)序列即為各個Mascon對應(yīng)的區(qū)域平均結(jié)果。

1.2 冰川質(zhì)量變化計(jì)算

念青唐古拉山脈為冰川聚集區(qū)，該地區(qū)質(zhì)量變化主要由冰川消融主導(dǎo)。因此，根據(jù)水量平衡原理，冰川質(zhì)量變化可表示為

ΔICE=ΔTWS-ΔSM-ΔSWE

(3)

式中，ΔICE為冰川儲量變化；ΔTWS為陸地水儲量變化；ΔSM為土壤水儲量變化；ΔSWE為雪水當(dāng)量變化。其中，ΔTWS為GRACE、GRACE-FO計(jì)算結(jié)果，主要包括地下水、土壤水、積雪水、地表徑流、植被水等。ΔSM和ΔSWE均可在GLDAS水文模型中得到。

GLDAS(global land data assimilation system)水文模型由美國宇航局戈達(dá)德太空飛行中心和美國國家環(huán)境預(yù)報(bào)中心聯(lián)合建立，包括NOAH、VIC、CLM、MOSAIC 4種陸地表面模式。本文采用2003年1月至2019年12月Noah陸面模式的Version1數(shù)據(jù)[19]，時間分辨率為1個月，空間分辨率為1°×1°。選取該數(shù)據(jù)的4層土壤水(0～200 mm)、雪水當(dāng)量和冠層水計(jì)算區(qū)域地表水儲量變化。

1.3 氣象數(shù)據(jù)處理

降水?dāng)?shù)據(jù)采用美國國家宇航局和日本空間發(fā)展局提供的TRMM(tropical rainfall measuring mission)全球降雨數(shù)據(jù)產(chǎn)品[20]，空間分辨率為0.25°×0.25°，時間分辨率為1個月。該產(chǎn)品提供的是月降雨率(mm/h)數(shù)據(jù)，為了反映降水和水儲量變化的關(guān)系，將其轉(zhuǎn)換成對應(yīng)月份的降雨量。為了保持與GRACE數(shù)據(jù)空間分辨率的一致，采用三次樣條插值方法將降雨數(shù)據(jù)插值成1°×1°，然后扣除各個格網(wǎng)點(diǎn)上的平均值，得到中國區(qū)域1°×1°的月降雨異常。

溫度數(shù)據(jù)采用美國國家海洋和大氣管理局提供的戈達(dá)德太空研究所對全球表面溫度分析數(shù)據(jù)集[21](GISS Surface Temperature Analysis，GISTEMP)，數(shù)據(jù)的空間分辨率為2°×2°，時間分辨率為1月，覆蓋范圍是89.0N—80.0S,1.0E—359.5E，該數(shù)據(jù)集提供了1980年至今的所有數(shù)據(jù)，具有較好的時間連續(xù)性。同降雨數(shù)據(jù)的處理步驟類似，得到中國區(qū)域的月溫度異常。

1.4 冰川信息提取

1.4.1 冰川面積提取

冰川面積提取主要借助Landsat光學(xué)遙感影像，利用冰雪在可見光波段的高反射特性和在短波紅外的強(qiáng)吸收特性，通過波段比值法自動提取冰川邊界。波段比值法的關(guān)鍵在于閾值選取，經(jīng)過大量試驗(yàn)得到閾值為2.3時可以更加有效區(qū)分冰川和陸地，但對水體和陰影錯分較多。為了獲得更加準(zhǔn)確的結(jié)果，采用目視解譯法對錯分的水體和陰影進(jìn)行修正，得到較為精確的冰川邊界。

冰川研究中遙感影像篩選應(yīng)滿足少云、無積雪覆蓋、消融末期3項(xiàng)原則，基于此原則經(jīng)過大量篩選，最終選取了2003-12-20(TM)、2017-10-07(OLI)兩景Landsat影像進(jìn)行預(yù)處理，主要包括輻射定標(biāo)、大氣校正和影像裁剪等。

1.4.2 冰川儲量估計(jì)

本文使用國家青藏高原科學(xué)數(shù)據(jù)中心(http:∥data.tpdc.ac.cn/zh-hans/)提供的第三極地區(qū)冰川表面高程變化數(shù)據(jù)產(chǎn)品[22]。該數(shù)據(jù)基于2000年的SRTM(shuttle radar topography mission)和2015年前后ASTER立體相對數(shù)據(jù)制成，估計(jì)了第三極地區(qū)40多個典型冰川表面高程變化，空間分辨率為30 m。通過獲取念青唐古拉山脈西段冰川范圍內(nèi)每個像元的高程變化值，結(jié)合式(4)可以得到2000—2015年間該地區(qū)的冰川體積變化ΔV

(4)

式中，Δhi表示第i個像元的高差；sp為單個像元面積；n為像元的個數(shù)。

2 結(jié)果與分析

2.1 冰川質(zhì)量變化

圖1為GRACE和GLDAS反演的念青唐古拉山脈陸地水儲量和地表水儲量變化序列。從圖1可看出，陸地水儲量變化和地表水儲量變化均表現(xiàn)出明顯的周期性，分別在冬季和夏季達(dá)到最小值和最大值。在時間尺度上，GRACE的結(jié)果和GLDAS的結(jié)果具有較好的一致性，但GRACE結(jié)果的振幅較大，這是因?yàn)镚RACE的結(jié)果是綜合各種因素后的水儲量變化，而GLDAS水文模型只包含表層土壤水變化，不包括深層土壤水、冰川、地下水、地表河流及湖泊等水儲量變化[23]。對GRACE和GLDAS的結(jié)果進(jìn)行最小二乘線性擬合，得到他們的變化速率分別為-2.76 cm/a和-0.48 cm/a，這說明地表水儲量變化對念青唐古拉山脈整體質(zhì)量虧損影響較小。

圖1 念青唐古拉山脈陸地水儲量變化和地表水儲量變化曲線Fig.1 Change curve of land water reserve and surface water reserve in Nyainqentanglha Mountain

念青唐古拉山脈不存在大型河流和湖泊，地下水和深層土壤水在凍土層的保護(hù)和補(bǔ)給也相對穩(wěn)定，因此，計(jì)算的冰川質(zhì)量變化結(jié)果受其他因素的影響較小。圖2給出了念青唐古拉山脈冰川質(zhì)量變化序列，其中藍(lán)線表示原始冰川質(zhì)量變化、紅線表示經(jīng)過泄露誤差改正的結(jié)果。圖2中藍(lán)線的變化速率為-2.48 cm/a，對總質(zhì)量變化的貢獻(xiàn)達(dá)到了89%，表明念青唐古拉山脈2003—2019年冰川質(zhì)量虧損嚴(yán)重，尤其是2013年之后念青唐古拉山脈冰川有加速消融趨勢，這一結(jié)果與文獻(xiàn)[10]的結(jié)論一致。此外，在2015年1—3月間念青唐古拉山脈地區(qū)存在一個明顯的質(zhì)量變化異常信號，這一異常在地表水儲量變化序列中未出現(xiàn)，結(jié)合2015年3月前后研究區(qū)附近地震的分布情況，認(rèn)為該質(zhì)量異常可能與地球內(nèi)部的物質(zhì)遷移有關(guān)。

圖2 冰川質(zhì)量變化和空域法改正曲線Fig.2 The curve of glacier mass change and spatial domain method correction

對比圖2兩條曲線能夠看出(圖中坐標(biāo)軸不一樣)，兩條曲線變化趨勢基本一致但振幅相差較大，這與空域法建立的Mascon起到的質(zhì)量約束作用有關(guān)。圖2紅線為泄露改正后的念青唐古拉山脈冰川儲量變化曲線，其變化速率為-15.53 cm/a，同時計(jì)算了念青唐古拉山脈西段和東段兩個Mascon的質(zhì)量變化速率分別為-20.23 cm/a和-10.82 cm/a，其中念青唐古拉山脈東側(cè)變化速率結(jié)果與文獻(xiàn)[24]使用空域法的結(jié)果基本一致。

為探究冰川質(zhì)量變化及影響機(jī)制，本文提取了夏季溫度異常和降水異常的變化曲線(圖3、圖4)。由圖3可以看出，2003—2019年間夏季溫度呈緩慢上升趨勢，研究時段內(nèi)平均溫度大約升高了1℃，而冰川質(zhì)量變化明顯減少。從整體看，兩條曲線具有明顯的相關(guān)性，其相關(guān)系數(shù)為-0.63。同時，計(jì)算了春季、秋季和冬季溫度異常與冰川質(zhì)量變化曲線的相關(guān)系數(shù)，相關(guān)性均不明顯，其中秋季的相關(guān)系數(shù)較大也僅為-0.36。由圖4可以看出，念青唐古拉山脈降水有緩慢減少的趨勢，表明溫度升高和降水減少是致念青唐古拉山脈冰川質(zhì)量虧損的主要原因。

圖3 念青唐古拉山脈夏季冰川質(zhì)量變化和溫度異常曲線Fig.3 The glacier quality change and temperature abnormal curves in the Nyainqentanglha Mountains in summer

圖4 念青唐古拉山脈夏季降水異常Fig.4 Abnormal precipitation in the Nyainqentanglha Mountain in summer

2.2 冰川面積變化

念青唐古拉山脈東段為海洋型冰川，冰川活動劇烈，常年云量較多，很難獲取同一尺度、覆蓋整個范圍的光學(xué)影像。因此，在冰川面積及儲量變化研究中以念青唐古拉山主峰所在區(qū)域的冰川為研究對象，該冰川位于念青唐古拉山脈西段(89.89°E—90.84°E，29.83°N—30.57°N)，分布獨(dú)立，靠近納木錯湖泊，平均海拔約5500 m，是較典型的亞大陸型冰川[25]。

利用兩幅Landsat影像獲取了念青唐古拉山脈西段冰川邊界，經(jīng)統(tǒng)計(jì)2003年冰川總面積為619.305 km2，2017年時為547.003 km2，這14年間冰川總面積退縮了72.302 km2，減少比例約為11.67%，年均變化速度為-5.16 km2/a。另外，從文獻(xiàn)[26]得知，2010年念青唐古拉山脈西段冰川總面積為571.81 km2，就2003—2010年和2010—2017年2個時段而言，存在明顯的變化差異，2003—2010年冰川面積的變化速率為-6.785 km2/a，而2010—2017年為-3.544 km2/a。

為探究不同海拔冰川面積變化，利用90 m分辨率的SRTM數(shù)字高程模型(digital elevation model，DEM)數(shù)據(jù)生成200 m等高線，統(tǒng)計(jì)相鄰兩等高線間冰川面積變化，見表1。由表1可以看出，海拔在5600～6200 m區(qū)間內(nèi)分布的冰川數(shù)量最多，占比91%以上。從年均變化率上看，低海拔區(qū)域冰川年均變化率最大，退縮趨勢明顯，這與低海拔地區(qū)溫度較高有關(guān)。隨著海拔的升高，溫度降低，冰川年均變化率逐漸減小。在海拔5600～5800 m區(qū)間內(nèi)，冰川面積變化最大，達(dá)到-34.56 km2(圖5)，該區(qū)間內(nèi)包含有大量冰川邊界，屬于易消融區(qū)域。在6200 m以上高海拔地區(qū)，溫度較低有利于冰川的發(fā)育，冰川面積不減反而緩慢增加。圖5中白色網(wǎng)格表示海拔在5600～5800 m區(qū)間范圍，底圖為2017年冰川區(qū)Landsat影像。

表1 念青唐古拉山脈西段不同海拔冰川面積變化統(tǒng)計(jì)

2.3 冰川儲量變化與不確定性分析

由第三極地區(qū)冰川表面高程變化數(shù)據(jù)產(chǎn)品獲取的念青唐古拉山脈西段冰川高程變化如圖6所示。由圖6可知，冰川的大部分冰舌區(qū)域均存在不同程度的表面高程減薄現(xiàn)象，其中北坡減薄現(xiàn)象更加明顯。經(jīng)統(tǒng)計(jì)，2000—2015年念青唐古拉山脈西段冰川平均減薄速率為38 cm/a，利用該數(shù)據(jù)產(chǎn)品計(jì)算的冰川儲量變化速率為-0.145 km3/a，相當(dāng)于-0.123 Gt/a。為了與遙感估算的冰川儲量變化進(jìn)行對比，盡可能減小觀測時段不一致帶來的誤差，通過計(jì)算得知，2002年8月至2015年12月念青唐古拉山脈西段冰川質(zhì)量變化速率為-12.78 cm/a，利用式(5)轉(zhuǎn)化為冰川體積變化為-0.088 km3/a，相當(dāng)于-0.075 Gt/a

(5)

式中，ΔV為冰川體積變化；ΔH為GRACE計(jì)算的等效水高的變化；S表示兩期冰川面積的平均值(583.154 km2)；ρg為冰川平均密度參數(shù)，取850 kg/m3。

衛(wèi)星重力和遙感估計(jì)的冰川儲量變化存在較大差異，其原因在于，GRACE重力信號中包含凍土融水、地下水、湖泊水等其他信號，這些信號的不確定性影響了GRACE觀測結(jié)果的精度，使冰川儲量的估計(jì)結(jié)果存在較大的誤差[27]。此外，冰川表面高程變化的不確定性對遙感估計(jì)結(jié)果影響較大，本文采用的冰川儲量估計(jì)方法依賴每個像元的高程變化，無論是SAR還是光學(xué)立體像對在地形起伏較大的區(qū)域獲取的DEM精度普遍不高，即便使用了DEM配準(zhǔn)方法減少了相位誤差和一些系統(tǒng)差的影響，但殘余誤差對結(jié)果的影響依然很大。

圖5 5600～5800 m區(qū)間冰川分布Fig.5 Distribution of glaciers in the 5600～5800 m interval

圖6 念青唐古拉山脈西段2000—2015年間冰川表面高程變化Fig.6 Changes of glacier surface elevation in the western Nyainqentanglha Mountain from 2000 to 2015

3 結(jié) 論

利用GRACE、GRACE-FO重力衛(wèi)星數(shù)據(jù)，結(jié)合GLDAS水文模型，反演了念青唐古拉山脈冰川質(zhì)量變化，借助氣象資料說明念青唐古拉山脈的常年質(zhì)量虧損和區(qū)域溫度升高有密切聯(lián)系。此外，以念青唐古拉山脈西段冰川為重點(diǎn)研究對象，結(jié)合Landsat光學(xué)遙感數(shù)據(jù)和表面高程產(chǎn)品對冰川的面積和儲量變化進(jìn)行詳細(xì)分析。主要結(jié)論如下。

(1) 2003—2019年間念青唐古拉山脈冰川質(zhì)量整體虧損嚴(yán)重，經(jīng)過空域法泄露誤差改正后，念青唐古拉山脈冰川質(zhì)量變化速率為-15.53 cm/a，其中念青唐古拉山脈西段和東段冰川質(zhì)量變化速率分別為-20.23 cm/a和-10.82 cm/a。對夏季的溫度變化和降水異常分析結(jié)果表明，念青唐古拉山脈的冰川退縮和降水減少、溫度升高有關(guān)。

(2) 念青唐古拉山脈西段冰川在不同海拔處冰川面積分布有所不同。海拔5400 m以下分布大量小型冰川，這部分冰川消融速度最快；海拔大于6200 m的區(qū)域溫度較低，有利于冰川的發(fā)育，冰川面積有所增加。海拔在5600～5800 m區(qū)間的冰川面積變化最大，其值為-34.56 km2，占總面積變化的48%，該區(qū)間內(nèi)分布大量冰川邊界和冰舌區(qū)域，在冰川中最易消融。

(3) 利用衛(wèi)星重力和遙感技術(shù)獲得的念青唐古拉山脈西段冰川儲量變化速率分別為-0.075 Gt/a和-0.12 Gt/a，二者變化趨勢相同，但數(shù)值差異較大，其主要原因是GRACE后處理誤差和冰川表面高程數(shù)據(jù)精度不高導(dǎo)致的。

(4) 衛(wèi)星重力和遙感技術(shù)的組合應(yīng)用可以優(yōu)勢互補(bǔ)，能夠充分發(fā)揮衛(wèi)星重力在大、中尺度的快速探測能力和遙感技術(shù)精確的冰川信息提取能力，實(shí)現(xiàn)對變化劇烈冰川的快速定位和精準(zhǔn)監(jiān)測，但在冰川儲量估計(jì)方面仍存在一定缺陷，有待進(jìn)一步研究。