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

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

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

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

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

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

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

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

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

1.1 GRACE反演理論

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

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

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

1.2 數(shù)據(jù)及后處理方法

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

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

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

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

2 時間序列分析

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

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

通過最小二乘（即εi平方和最?。┍憧蓪κ剑?）中未知參數(shù)進(jìn)行估算。選取的待估參數(shù)不同，對冰雪質(zhì)量變化趨勢的估算有較大影響［34］?；貧w參數(shù)xj的選取是通過假設(shè)檢驗和信息準(zhǔn)則AIC（akaike information criterion）和BIC（bayesian information criterion）進(jìn)行評估確定的［34］。假設(shè)檢驗可以判斷回歸參數(shù)選取的統(tǒng)計顯著性。作原假設(shè)H0：xj＝0，設(shè)置信區(qū)間為95%，置信水平α＝5%?？梢酝ㄟ^t檢驗和F檢驗判斷是否拒絕H0，從而檢驗參數(shù)xj（j＝1，2，3，4，…）回歸效果是否顯著。信息準(zhǔn)則是根據(jù)后驗概率大小進(jìn)行最優(yōu)選擇的一種決策方法，假設(shè)檢驗和信息準(zhǔn)則的計算模型為：

表1 回歸參數(shù)選擇的統(tǒng)計評價結(jié)果Tab.1 Results of regression parameter selection criteria evaluation

通過假設(shè)檢驗和信息準(zhǔn)則可以確定N＝1（即線性擬合）為喜馬拉雅地區(qū)經(jīng)過GIA 改正后的冰川質(zhì)量變化時間序列的最佳擬合模型，t檢驗、F檢驗和信息準(zhǔn)則確定線性擬合為最優(yōu)。N＝2（即二項式擬合）為喜馬拉雅地區(qū)經(jīng)過GIA 和GLDAS改正后的冰川質(zhì)量變化時間序列的最佳擬合模型（表1），其中t檢驗同時支持線性和二項式擬合，F(xiàn)檢驗和信息準(zhǔn)則確定二項式擬合為最優(yōu)選擇。

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

3 結(jié)果分析

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

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

［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］.地球物理學(xué)報，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 時變重力場反演南極冰蓋的質(zhì)量變化趨勢［J］.中國科學(xué)：地球科學(xué)，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衛(wèi)星RL05數(shù)據(jù)的南極冰蓋質(zhì)量變化分析［J］.地球物理學(xué)報，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，F(xiàn)amiglietti 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