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

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

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

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

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

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

1 研究方法

1.1 地表質量變化反演

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

(1)

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

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

(2)

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

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

1.2 冰川質量變化計算

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

ΔICE=ΔTWS-ΔSM-ΔSWE

(3)

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

GLDAS(global land data assimilation system)水文模型由美國宇航局戈達德太空飛行中心和美國國家環(huán)境預報中心聯(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)、雪水當量和冠層水計算區(qū)域地表水儲量變化。

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

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

溫度數(shù)據(jù)采用美國國家海洋和大氣管理局提供的戈達德太空研究所對全球表面溫度分析數(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光學遙感影像，利用冰雪在可見光波段的高反射特性和在短波紅外的強吸收特性，通過波段比值法自動提取冰川邊界。波段比值法的關鍵在于閾值選取，經(jīng)過大量試驗得到閾值為2.3時可以更加有效區(qū)分冰川和陸地，但對水體和陰影錯分較多。為了獲得更加準確的結果，采用目視解譯法對錯分的水體和陰影進行修正，得到較為精確的冰川邊界。

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

1.4.2 冰川儲量估計

本文使用國家青藏高原科學數(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ù)制成，估計了第三極地區(qū)40多個典型冰川表面高程變化，空間分辨率為30 m。通過獲取念青唐古拉山脈西段冰川范圍內每個像元的高程變化值，結合式(4)可以得到2000—2015年間該地區(qū)的冰川體積變化ΔV

(4)

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

2 結果與分析

2.1 冰川質量變化

圖1為GRACE和GLDAS反演的念青唐古拉山脈陸地水儲量和地表水儲量變化序列。從圖1可看出，陸地水儲量變化和地表水儲量變化均表現(xiàn)出明顯的周期性，分別在冬季和夏季達到最小值和最大值。在時間尺度上，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

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

圖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℃，而冰川質量變化明顯減少。從整體看，兩條曲線具有明顯的相關性，其相關系數(shù)為-0.63。同時，計算了春季、秋季和冬季溫度異常與冰川質量變化曲線的相關系數(shù)，相關性均不明顯，其中秋季的相關系數(shù)較大也僅為-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 冰川面積變化

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

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

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

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

由第三極地區(qū)冰川表面高程變化數(shù)據(jù)產(chǎn)品獲取的念青唐古拉山脈西段冰川高程變化如圖6所示。由圖6可知，冰川的大部分冰舌區(qū)域均存在不同程度的表面高程減薄現(xiàn)象，其中北坡減薄現(xiàn)象更加明顯。經(jīng)統(tǒng)計，2000—2015年念青唐古拉山脈西段冰川平均減薄速率為38 cm/a，利用該數(shù)據(jù)產(chǎn)品計算的冰川儲量變化速率為-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為冰川平均密度參數(shù)，取850 kg/m3。

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

圖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 結 論

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

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

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

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

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