基于GRACE數據估計近年喜馬拉雅冰川質量變化

史紅嶺 陸 洋 高春春，4 朱傳東，4 杜宗亮，4

1 中國科學院測量與地球物理研究所大地測量與地球動力學國家重點實驗室，武漢市徐東大街340號，430077

2 地理空間信息工程國家測繪地理信息局重點實驗室，北京市蓮花池西路28號，100830

3 極地測繪科學國家測繪地理信息局重點實驗室，武漢市珞喻路129號，430079

4 中國科學院大學，北京市玉泉路甲19號，100049

喜馬拉雅山區冰川面積占亞洲冰川面積的1／4［1，2］。隨著該地區氣溫逐漸升高［1，3，4］，冰 川大量縮減消融，將形成嚴重的生態問題。精確評估喜馬拉雅山區冰川質量變化具有重要意義。隨著GIS在冰凍圈變化研究中的應用以及遙感和全球定位系統技術在冰川變化測量中逐步展開［6－7］，何春陽［6］利用表面數字高程模型計算各網點的高度變化，進而逐點推求整個冰川的變化。李德平［8］利用航片、Landsat TM 影像等資料對小冰期以來羌塘高原中西部冰川變化進行圖譜分析。但是，喜馬拉雅地區現代冰川數量龐大且氣候惡劣，地面監測其物質平衡存在諸多困難［9］。

2002年發射的重力恢復與氣候實驗衛星GRACE，開創了以高時空分辨率觀測全球重力場的新紀元，并廣泛應用于監測兩極冰蓋［10－15］、山地冰川［16－21］等物質平衡的時空變化。在青藏高原及附近地區，Rodell等［22］利用GRACE 衛星重力資料和土壤濕度模型分析印度西北部地下水儲量變化。Moiwo等［23］利用GRACE 數據和氣候模型分析青藏高原地區水儲量變化。Matsuo等［24］利用GRACE資料研究2003～2009年間亞洲地區冰川質量損失情況。Jacob 等［18］基于2003～2010年間GRACE Level2衛星重力場資料，采用冰川區域的Mascon數據計算除去極地冰蓋外面積大于100km2的山地冰川區域質量變化。Yi等［20］基于2003～2012年GRACE 數據得到亞洲區域冰川的質量變化。但是上述結果相差甚大，需要進一步研究分析。

本文基于CSR 最新發布的GRACE RL05時變重力場模型，利用點質量方法在頻域上估計喜馬拉雅地區2002－04～2014－03冰川質量變化，并通過假設檢驗和信息準則確定合理的最小二乘估計參數，討論喜馬拉雅地區冰川質量變化。

1 GRACE反演理論與后處理方法

1.1 GRACE反演理論

根據Wahr理 論［25］，GRACE時變重力場模型能夠反演得到地球表面質量薄層上任意一點的面密度變化：

其中，θ和φ分別為地心余緯和東經，a為地球平均半徑，ρave為地球平均密度，l和m為階和次，為正規化的締合勒讓德函數，kl為載荷勒夫數，ΔClm和ΔSlm分別為時變重力場模型的無量綱正規化球諧系數相對于均值的變化量。

質量變化的結果通常以等效水柱高Hw（θ，φ）表示：Hw（θ，φ）＝Δσ（θ，φ）／ρw，其中ρw為水密度。

1.2 數據及后處理方法

采用2002－04～2014－03GRACE RL05數據，球諧系數截斷到60階，并扣除非潮汐大氣、高頻海洋信號、各種潮汐、固體潮和固體極潮等影響［26］。由于GRACE數據中球諧系數C20項解算精度較低，本文采用SLR 測得的C20項對其作替換［27］。為了壓制GRACE 數據高階球諧系數中的較大噪聲以及相關誤差［25，28］，首先對球諧系數作去相關處理，對于次數小于等于10的球諧系數保持不變，次數大于10的奇偶項系數分別采用滑動窗口為7的3 次多項式擬合［29］。為了最大可能地利用數據原有的分辨率獲取小面積冰川質量變化，采用盡量小的平滑半徑以壓制較高階的誤差影響，采用平滑半徑為150km 高斯濾波對球諧系數作平滑處理［18，25］。考慮到GRACE 的計算結果無法分離出冰川均衡調整GIA 的影響，本文將Paulson 等［30］的GIA 模型展開成球諧系數并截取到60階，采用和GRACE數據相同的去相關和平滑處理。為去除水文信號對GRACE計算結果信號泄漏的影響，將GLDAS水文模型轉換為球諧系數并截取到60階，采用與GRACE數據相同的去相關和平滑處理［［31］。

為了得到喜馬拉雅冰雪質量變化，把喜馬拉雅區域分成多個小的任意定義的點質量（mascons）（圖1），然后通過擬合GRACE 時變重力場得到每個mascon的質量變化［18，32］。本文從倫道夫冰川目錄提供的冰川數據中提取喜馬拉雅地區冰川的分布范圍，并剔除面積小于100km2的冰川［33］。對于任意mascon（i），假設地球表面質量薄層任意一點的等效水柱高變化呈現均勻分布狀態，根據等效水柱高變化與球諧系數之間的關系，可求得單位等效水柱高分布的mascon對應的球諧系數［25］：

圖1 喜馬拉雅地區冰川范圍，不同顏色的方框代表不同的masconig.1 Ice－covered region in the Himalayas，each colored box－shaped region represents a mascon

假設任意mascon對應的未知等效水柱高變化為Mi（t），當E取得最小時可解得Mi（t）［18］：

2 時間序列分析

經過上述處理即可得到喜馬拉雅冰雪質量變化的時間序列。為了估算變化的趨勢，還需要對時間序列進行多項式擬合［34］：

式中，y（Δti）為ti時的剩余質量，Δti＝ti－tref，ti為時間點，tref為參考時間，xj（j＝0，…N）為待估多項式系數，Ck、Sk（k＝1，…M）分別為待估的余弦、正弦振幅，fk為頻率，εi為誤差項。

通過最小二乘（即εi平方和最小）便可對式（5）中未知參數進行估算。選取的待估參數不同，對冰雪質量變化趨勢的估算有較大影響［34］。回歸參數xj的選取是通過假設檢驗和信息準則AIC（akaike information criterion）和BIC（bayesian information criterion）進行評估確定的［34］。假設檢驗可以判斷回歸參數選取的統計顯著性。作原假設H0：xj＝0，設置信區間為95%，置信水平α＝5%。可以通過t檢驗和F檢驗判斷是否拒絕H0，從而檢驗參數xj（j＝1，2，3，4，…）回歸效果是否顯著。信息準則是根據后驗概率大小進行最優選擇的一種決策方法，假設檢驗和信息準則的計算模型為：

表1 回歸參數選擇的統計評價結果Tab.1 Results of regression parameter selection criteria evaluation

通過假設檢驗和信息準則可以確定N＝1（即線性擬合）為喜馬拉雅地區經過GIA 改正后的冰川質量變化時間序列的最佳擬合模型，t檢驗、F檢驗和信息準則確定線性擬合為最優。N＝2（即二項式擬合）為喜馬拉雅地區經過GIA 和GLDAS改正后的冰川質量變化時間序列的最佳擬合模型（表1），其中t檢驗同時支持線性和二項式擬合，F檢驗和信息準則確定二項式擬合為最優選擇。

圖2 喜馬拉雅地區冰川質量變化時間序列和變化趨勢Fig.2 The time series andits trend in Himalayas

3 結果分析

Jacob等［18］利用2003～2011期間GRACE數據估計亞洲區域的冰川質量變化，其中喜馬拉雅山和喀喇昆侖山區域的冰川質量變化約為－5±6 Gt／a。而Duan［35］利 用2002～2012年間GRACE時間序列得到的喜馬拉雅地區冰川質量變化約為－20±6Gt／a。Jacob［18］采用了CLM4和GLDAS／NOAH2.7.1的均值扣除水文影響，Duan［35］則是扣除GLDAS／NOAH 與WGHM 模型的均值。他們都考慮了LIA 模型的影響。根據Duan［35］的研究，喜馬拉雅地區LIA 的影響約為－1.5±1.1Gt／a，本文暫不考慮。基于133個月的GRACE月時變重力場數據，經過去相關和高斯濾波等相關后處理，基于劃分的macon反演得到喜馬拉雅地區2002－04～2014－03 冰川質量變化時間序列。扣除季節性的年變化，基于上述時間序列分析即可得到冰雪質量的變化趨勢和加速度。圖2 為2002－04～2014－03 喜馬拉雅區域的冰川質量變化時間序列和變化趨勢，可以看出喜馬拉雅地區冰川質量存在明顯的季節性變化，一般在春季和秋季分別達到最小值和最大值。扣除Paul的GIA 模型影響后［30］，喜馬拉雅地區在2002－04～2014－03冰川質量融化趨勢為－9.08±5.39Gt／a（圖2（a））。

在喜馬拉雅地區，由于冰川的地理分布特征和CSR月重力場模型分辨率的影響，很難精確地完全分離冰川質量變化信號和周圍的地下水衰竭信號，但是為了得到冰川質量變化的穩健估計，必須從GRACE觀測得到的質量變化趨勢中移除周圍地下水變化的影響［35］。進一步扣除GLDAS水文模型影響后的融化趨勢為－8.26±4.61Gt／a，加速度為－3.54±1.25Gt／a2（圖2b），消融具有加速趨勢。根據Yao的研究［5］，印度季風為青藏高原南部和喜馬拉雅山脈中東部帶來降雪，但是在過去的幾十年間，季風有所減弱，同時季風受到季節的影響，對冰川質量變化影響很大。此外降水是引起青藏區域水儲量變化的重要因素，同時也是水文模型模擬中的重要輸入因子，1979～2010年間喜馬拉雅的降水呈現下降趨勢［5］。但是目前對降水和季風的影響，尚未有清晰的認識。由于所采用的時間跨度、水文模型和反演方法不同，數值結果與Jacob和Duan等研究稍有不同，不過總體來說喜馬拉雅地區冰川主要呈現加速退縮趨勢，其中GIA 隆升和水文的影響有待于進一步詳細研究。

［1］張東啟，效存德，秦大河.近幾十年來喜馬拉雅山冰川變化及其對水資源的影響［J］.冰川凍土，2009（5）：885－895（Zhang Dongqi，Xiao Cunde，Qin Dahe.Himalayan Glaciers Fluctuation Over the Latest Decades and Its Impact on Water Resources［J］.Journal of Glaciology and Geocryology，2009，31（5）：885－895）

［2］王宗太，蘇宏超.世界和中國的冰川分布及其水資源意義［J］.冰川凍土，2003，25（5）：498－503（Wang Zongtai，Su Hongchao.Glaciers in the World and China：Distribution and Their Significance as Water Resources［J］.Journal of Glaciology and Geocryology，2003，25（5）：498－503）

［3］Hasnain S I.Status of the Glacier Research in the HKH Region［C］.ICIMOD，Kathmandu，Nepal，2000

［4］Ding Y，Griggs D J，Noguer M，et al.Climate Change 2001：the Scientific Basis［M］.Cambridge：Cambridge University Press，2001

［5］Yao T，Thompson L，Yang W，et al.Different Glacier Status with Atmospheric Circulations in Tibetan Plateau and Surroundings［J］.Nature Climate Change，2012，2（9）：663－667

［6］何春陽，丁永建，李新.冰川變化可視化計算方法研究［J］.冰川 凍 土，1999，21（2）：169－174（He Chunyang，Ding Yongjian，Li Xin.A Visualized Computation Method for Glacier Variation［J］.Journal of Glaciology and Geocryology，1999，21（2）：169－174）

［7］Li Z，Sun W，Zeng Q.Measurements of Glacier Variation in the Tibetan Plateau Using Landsat Data［J］.Remote Sensing of Environment，1998，63（3）：258－264

［8］李德平，王利平，劉時銀，等.小冰期以來羌塘高原中西部冰川變化圖譜分析［J］.冰川凍土，2009，31（1）：40－47（Li Deping，Wang Liping，Liu Shiyin，et al.Tupu Analysis of the Spatio－Temporal Glacier Variations in the Central and Western Qangtang Plateau Since the Little Ice Age［J］.Journal of Glaciology and Geocryology，2009，31（1）：40－47）

［9］Gardner A S，Moholdt G，Wouters B，et al.Sharply Increased Mass Loss from Glaciers and Ice Caps in the Canadian Arctic Archipelago［J］.Nature，2011，473（7 347）：357－360

［10］Velicogna I，Wahr J.Acceleration of Greenland Ice Mass Loss in Spring 2004［J］.Nature，2006，443（7 109）：329－331

［11］Velicogna I，Wahr J.Time Variable Gravity Observations of Ice Sheet Mass Balance：Precision and Limitations of the GRACE Satellite Data［J］.Geophysical Research Letters，2013，40（12）：3055－3063

［12］Chen J L，Wilson C R，Blankenship D，et al.Accelerated Antarctic Ice Loss from Satellite Gravity Measurements［J］.Nature Geoscience，2009，2（12）：859－862

［13］鄂棟臣，楊元德，晁定波.基于GRACE 資料研究南極冰蓋消減對海平面的影響［J］.地球物理學報，2009，52（9）：2 222－2 228（E Dongchen，Yang Yuande，Chao Dingbo.The Sea Level Change from the Antarctic Ice Sheet based on GRACE［J］.Progress in Geophysics，2009，52（9）：2 222－2 228）

［14］羅志才，李瓊，張坤，等.利用GRACE 時變重力場反演南極冰蓋的質量變化趨勢［J］.中國科學：地球科學，2012，42（10）：1 590－1 596（Luo Zhicai，Li Qiong，Zhang Kun，et al.Trend of Mass Change in the Antarctic Ice Sheet Recovered from the GRACE Temporal Gravity Field［J］.Science China Earth Sciences，2012，42（10）：1 590－1 596）

［15］鞠曉蕾，沈云中，張子占.基于GRACE衛星RL05數據的南極冰蓋質量變化分析［J］.地球物理學報，2013，56（9）：2 918－2 927（Ju Xiaolei，Shen Yunzhong，Zhang Zizhan.Antarctic Ice Mass Change Analysis based on GRACE RL05data［J］.Chinese Journal Geophysics，2013，56（9）：2 918－2 927）

［16］Luthcke S B，Arendt A A，Rowlands D D，et al.Recent Glacier Mass Changes in the Gulf of Alaska Region from GRACE Mascon Solutions［J］.Journal of Glaciology，2008，54（188）：767－777

［17］Riva R E M，Bamber J L，Lavallée D A，et al.Sea Level Fingerprint of Continental Water and Ice Mass Change from GRACE［J］.Geophysical Research Letters，2010，37（19）

［18］Jacob T，Wahr J，Pfeffer W T，et al.Recent Contributions of Glaciers and Ice Caps to Sea Level Rise［J］.Nature，2012，482（7 386）：514－518

［19］Baumann S，Menzel A，Seitz F.Estimating Glacier Mass Changes by GRACE Satellite Gravimetry in the Pamir and Tien－Shan mountains［C］.EGU General Assembly Conference Abstracts，2012

［20］Yi S，Sun W.Evaluation of Glacier Changes in High Mountain Asia Based on 10year GRACE RL05Models［J］.Journal of Geophysical Research：Solid Earth，2014，119（3）：2 504－2 517

［21］Chen J.Monitoring Terrestrial Water Cycle in Tibetan Plateau Using Satellite Gravimetry［J］.Acta Geologica Sinica－English Edition，2013，87（s1）：626－671

［22］Rodell M，Velicogna I，Famiglietti J S.Satellite－based Estimates of Groundwater Depletion in India［J］.Nature，2009，460（7258）：999－1 002

［23］Moiwo J P，Yang Y，Tao F，et al.Water Storage Change in the Himalayas from the Gravity Recovery and Climate Experiment（GRACE）and an Empirical Climate Model［J］.Water Resources Research，2011，47（7）

［24］Matsuo K，Heki K.Time－Variable Ice Loss in Asian High Mountains from Satellite Gravimetry［J］.Earth and Planetary Science Letters，2010，290（1）：30－36

［25］Wahr J，Molenaar M，Bryan F.Time Variability of the Earth’s Gravity Field：Hydrological and Oceanic Effects and their Possible Detection Using GRACE［J］.Journal of Geophysical Research：Solid Earth（1978－2012），1998，103（B12）：30 205－30 229

［26］Bettadpur S.CSR Level－2Processing Standards Document for Product Release 04［R］.Rep GRACE，2007

［27］Cheng M，Tapley B D.Variations in the Earth’s Oblateness During the Past 28Years［J］.Journal of Geophysical Research：Solid Earth（1978－2012），2004，109（B9）

［28］Swenson S，Wahr J.Monitoring Changes in Continental Water Storage with GRACE［J］.Space Science Reviews，2003，108（1－2）：345－354

［29］Swenson S，Wahr J.Post Pprocessing Removal of Correlated Errors in GRACE Data［J］.Geophysical Research Letters，2006，33（8）

［30］Paulson A，Zhong S，Wahr J.Inference of Mantle Viscosity from GRACE and Relative Sea Level Data［J］.Geophysical Journal International，2007，171（2）：497－508

［31］Rodell M，Houser P R，Jambor U，et al.The Global Land Data Assimilation System［J］.Bulletin of the American Meteorological Society，2004，85（3）：381－394

［32］Tiwari V M，Wahr J，Swenson S.Dwindling Groundwater Resources in Northern India，from Satellite Gravity Observations［J］.Geophysical Research Letters，2009（18）：252－260

［33］Pfeffer W T，Arendt A A，Bliss A，et al.The Randolph Glacier Inventory：a Globally Complete Inventory of Glaciers［J］.Journal of Glaciology，2014，60（221）：537－552

［34］Baur O.On the Computation of Mass－Change Trends from GRACE Gravity Field Time－Series［J］.Journal of Geodynamics，2012，61：120－128

［35］Duan J.Global Ice Mass Balance and its Contribution to Early Twenty－first Century Sea Level Rise［D］.Ohio：The Ohio State University，2013