衛星重力與光學遙感組合的念青唐古拉山脈冰川變化分析

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

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

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

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

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

1 研究方法

1.1 地表質量變化反演

利用重力場模型的球諧系數變化計算地球表面質量變化的模型為

(1)

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

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

(2)

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

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

1.2 冰川質量變化計算

念青唐古拉山脈為冰川聚集區，該地區質量變化主要由冰川消融主導。因此，根據水量平衡原理，冰川質量變化可表示為

ΔICE=ΔTWS-ΔSM-ΔSWE

(3)

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

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

1.3 氣象數據處理

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

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

1.4 冰川信息提取

1.4.1 冰川面積提取

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

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

1.4.2 冰川儲量估計

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

(4)

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

2 結果與分析

2.1 冰川質量變化

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

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

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

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

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

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

圖3 念青唐古拉山脈夏季冰川質量變化和溫度異常曲線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 冰川面積變化

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

利用兩幅Landsat影像獲取了念青唐古拉山脈西段冰川邊界，經統計2003年冰川總面積為619.305 km2，2017年時為547.003 km2，這14年間冰川總面積退縮了72.302 km2，減少比例約為11.67%，年均變化速度為-5.16 km2/a。另外，從文獻[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數字高程模型(digital elevation model，DEM)數據生成200 m等高線，統計相鄰兩等高線間冰川面積變化，見表1。由表1可以看出，海拔在5600～6200 m區間內分布的冰川數量最多，占比91%以上。從年均變化率上看，低海拔區域冰川年均變化率最大，退縮趨勢明顯，這與低海拔地區溫度較高有關。隨著海拔的升高，溫度降低，冰川年均變化率逐漸減小。在海拔5600～5800 m區間內，冰川面積變化最大，達到-34.56 km2(圖5)，該區間內包含有大量冰川邊界，屬于易消融區域。在6200 m以上高海拔地區，溫度較低有利于冰川的發育，冰川面積不減反而緩慢增加。圖5中白色網格表示海拔在5600～5800 m區間范圍，底圖為2017年冰川區Landsat影像。

表1 念青唐古拉山脈西段不同海拔冰川面積變化統計

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

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

(5)

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

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

圖5 5600～5800 m區間冰川分布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 結 論

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

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

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

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

(4) 衛星重力和遙感技術的組合應用可以優勢互補，能夠充分發揮衛星重力在大、中尺度的快速探測能力和遙感技術精確的冰川信息提取能力，實現對變化劇烈冰川的快速定位和精準監測，但在冰川儲量估計方面仍存在一定缺陷，有待進一步研究。