龐書劍，柯長青，周興華，張其兵，范宇賓，喻薛凝

1.南京大學地理與海洋科學學院,南京 210023;

2.自然資源部第一海洋研究所,青島 266000;

3.湖南文理學院資源環(huán)境與旅游學院,常德 415000

1 引言

冰川作為冰凍圈的重要組成部分，對氣候變化極為敏感，了解和掌握其動態(tài)變化，對于指示氣候具有重要的作用（Bojinski等，2014；Qin等，2018；楊瑞敏等，2019）。中國青藏高原地區(qū)分布有大量的山地冰川，由于溫室效應(yīng)不斷積累，其物質(zhì)平衡整體呈減少的趨勢（王寧練等，2019），但是不同地區(qū)的冰川也存在差異。自1970年以來，青藏高原南部和東南部冰川出現(xiàn)嚴重退縮，西部喀喇昆侖山、西昆侖山地區(qū)冰川則保持了相對穩(wěn)定甚至前進的狀態(tài)，冰川物質(zhì)變化整體較?。╖hou等，2018）。2000年—2016年整個亞洲高山區(qū)除了西昆侖的冰川處于異常積累狀態(tài)，其余地區(qū)冰川均處于負物質(zhì)平衡狀態(tài)，東昆侖山與青藏高原內(nèi)部地區(qū)（-0.14±0.08 m·w·e/a）消融較為劇烈（Brun等，2017）。其中，位于東昆侖山西段的昆侖山第二大冰川分布區(qū)木孜塔格峰1972年—2011年冰川呈物質(zhì)虧損狀態(tài)（-0.06±0.01 m·w·e/a）（蔣宗立等，2019），而2000年—2015年東昆侖山中心區(qū)域的7條主要冰川以馬蘭山和布喀達坂峰的冰川物質(zhì)損失最為顯著，其余冰川則處于接近平衡或微弱負平衡狀態(tài)，冰川的消融對周圍的勒斜武擔湖與可可西里湖等產(chǎn)生了一定的補給效應(yīng)（Zhou等，2019）。

冰川物質(zhì)平衡作為反應(yīng)冰川消融積累狀態(tài)的重要指標，通常根據(jù)冰川表面的高程變化對其進行估算，而冰面高程監(jiān)測的遙感手段主要包括合成孔徑雷達干涉測量、激光雷達測高技術(shù)以及光學立體像對（Zhang等，2019）。過去主要采用光學立體像對獲取冰川表面高程，如KH-9（Keyhole-9）解密航天數(shù)據(jù)，由于成像時間久遠并位于偏遠地區(qū)，控制點的選取較為困難，影像畸變難以糾正，生成的DEM精度較差（Pieczonka和Bolch，2015）。近年，也有很多學者利用ASTER（Advanced Spaceborne Thermal Emission and Reflection Radiometer）和SPOT立體像對數(shù)據(jù)對青藏高原地區(qū)冰川物質(zhì)平衡進行研究，得到的DEM結(jié)果精度有了很大提高（Pieczonka等，2013；Zhou等，2017），但光學影像仍受到云層與復雜天氣狀況的影響，存在一定的局限性。相比于光學方法，合成孔徑雷達能夠穿透云層，不受天氣條件限制對冰川進行全天時全天候觀測，并且具有更高的精度（Li和Lin，2017），ERS1/2由于最短重返周期可達1天，過去常被用于冰面高程提取，但是1天間隔內(nèi)冰川流動的殘余相位難以剔除（黃丁發(fā)和劉國祥，2009），隨著德國宇航局在2007年與2010年分別發(fā)射了TerraSAR-X與TanDEM-X兩顆并列飛行的X波段雷達衛(wèi)星，由于其雙站飛行測量模式幾乎完全消除了由冰川表面復雜的運動帶來的時間失相干，大幅提高了雷達干涉測量的精度，并且基于該數(shù)據(jù)的雙站干涉測量方法還解決了常規(guī)干涉測量解纏相位跳躍、不連續(xù)，以及軌道基線誤差引起的長波趨勢面等問題（孫亞飛等，2016）。Liu等（2019）等利用TanDEM-X數(shù)據(jù)研究了青藏高原普若崗日冰川在2011年—2016年冰川損失加劇的情況，蔣宗立等（2018）也采用雙站干涉測量的方法研究了黃河源區(qū)阿尼瑪卿山典型冰川物質(zhì)平衡。激光雷達高度計相比于其他遙感測高方法，具有更高的測高精度和更加廣泛的覆蓋范圍，也常被用于冰川物質(zhì)平衡監(jiān)測，K??b等（2012）利用ICESat估算了喜馬拉雅山區(qū)域2003年—2009年冰川物質(zhì)平衡，Wang等（2017a）也根據(jù)同樣的數(shù)據(jù)研究了亞洲高山區(qū)冰川厚度的大規(guī)模季節(jié)性變化。

馬蘭山位于昆侖山中部，常年被冰川覆蓋，總儲水量2.2×108m3，平均年融水量6.8×106m3，周圍有許多大的湖盆，如太陽湖、科考湖、可可西里湖、飲馬湖等（姜珊等，2012），其變化對于周圍的氣溫、降水以及附近湖泊河流的融水補給具有重要影響。本文采用TerraSAR-X/TanDEM-X合成孔徑雷達數(shù)據(jù)進行雙站干涉測量，并結(jié)合ICESat-2激光雷達高度計數(shù)據(jù)與SRTM DEM，研究了馬蘭山地區(qū)2000年—2020年冰川表面高程變化與物質(zhì)平衡，分析了兩種高程測量方法的精度誤差和有效性，并探討該區(qū)域主要冰川近20年來的表面高程變化與物質(zhì)平衡分布特征，以及冰川對氣候變化的響應(yīng)關(guān)系。

2 研究區(qū)與數(shù)據(jù)

2.1 研究區(qū)

馬蘭山地處青藏高原可可西里北部，屬于東昆侖山系，其北部是太陽湖，東南是可可西里湖，分布范圍介于35°46＇—35°52＇N，90°32＇—90°45＇E（姜珊等，2012），平均海拔5790 m，最高可達6056 m，地理位置和高程如圖1所示。馬蘭山一帶身處內(nèi)陸，主要受到西風環(huán)流的影響（姜珊，2012），同時冰川發(fā)育。根據(jù)GLIMS（Global Land Ice Measurements from Space）項目發(fā)布的全球冰川輪廓的完整目錄RGI 6.0（Randolph Glacier Inventory 6.0）顯示，馬蘭山整個山體被冰川覆蓋，共由41條冰川組成，面積達188 km2，形成了一個很大的冰帽冰川群，屬極大陸型冰川，雪線海拔在5340—5540 m之間；其山脊偏向北側(cè)，呈東西向分布，東西長約38 km，南北寬約15 km（蒲健辰等，2001）。馬蘭山南坡地勢平緩，北坡地勢陡峭，南坡冰川面積遠大于北坡，以冰帽型冰川和寬大平緩的山谷冰川為主體，冰舌短小而拱圓；北坡雖然冰川作用面積較小，但冰川作用高差較大，可達800—1000 m，其冰川類型主要為山谷冰川，冰帽整體規(guī)模較?。ㄖx自楚等，2000）。

圖1 馬蘭山地區(qū)高程與冰川分布（圖中黑色邊界為RGI亞洲高山區(qū)山脈劃分界限）Fig.1 Map showing elevation of Malan Mountain and glacier distribution（The black boundary in the image marks the boundaries of the RGI Asian High Mountain Range）

2.2 數(shù)據(jù)

2.2.1 TerraSAR-X/TanDEM-X合成孔徑雷達

TerraSAR-X與TanDEM-X是德國宇航中心分別于2007年和2010年發(fā)射的兩顆X波段雷達衛(wèi)星組，兩顆衛(wèi)星以螺旋軌道緊密編隊同步飛行，工作頻率為9.65 GHz，波長3.1 cm，入射角為22°—55°，共有3種觀測模式，分別是雙站模式（Bistatic mode）、單站模式（Monostatic mode）和交替雙站模式（Alternating bistatic mode），重返周期為11天，極化方式包括單、雙、全極化。本研究使用兩對經(jīng)過配準的雙站模式干涉數(shù)據(jù)對（CoSSC，Coregistered Single look Slant range Complex）進行干涉，極化方式為HH極化，空間分辨率約為12 m×12 m，此模式使用TSX（TerraSAR-X）或TDX（TanDEM-X）作為發(fā)射器，然后由兩個衛(wèi)星同時記錄散射的信號，避免了由于時間去相關(guān)和大氣干擾而可能產(chǎn)生的誤差，數(shù)據(jù)列表如表1所示。

表1 TerraSAR-X/TanDEM-X SAR影像參數(shù)Table 1 Parameters of TerraSAR-X/TanDEM-X Synthetic Aperture（SAR）images

2.2.2 ICESat-2激光雷達高度計

ICESat-2是美國國家航空航天局（NASA）于2018年9月15日發(fā)射的一顆激光雷達測高衛(wèi)星，重返周期為91天，軌道傾角為92°，主要用于測量冰蓋和冰川物質(zhì)變化、海冰厚度，以及全球森林和其他生態(tài)系統(tǒng)的植被高度（Smith等，2021）。其搭載的高級地形激光測高儀系統(tǒng)ATLAS（Advanced Topographic Laser Altimeter System）以500 km固定軌道高度10 kHz頻率沿軌每隔0.7 m發(fā)射綠色（532 nm）激光脈沖，每個脈沖點在地表的直徑約為14 m。每個脈沖會被衍射光學元件分成3對發(fā)射波束，每對波束間的距離為3.3 km，每對波束的左右波束的跨軌距離為90 m。本文使用ATL06（Advanced Topographic Laser Altimeter System 06）陸冰高度數(shù)據(jù)集，空間分辨率為20 m，以WGS84橢球為空間參考，其高度值為沿軌40 m分段的地表高程平均值。為了讓ICESat-2足跡點在冰川區(qū)覆蓋更加均勻，采用的數(shù)據(jù)時間為2018年11月至2019年4月、2019年11月至2020年4月，跨越春、冬兩個季節(jié)，與SRTM采集時間（2月）稍有差異，足跡點空間分布如圖2所示，冰川區(qū)足跡點數(shù)量分布圖3所示。

圖3 冰川區(qū)ICESat-2足跡點數(shù)量隨海拔變化（柱狀圖為ICESat-2在不同海拔高度足跡點數(shù)量，曲線為SRTM隨海拔分布）Fig.3 Changes in the number of ICESAT-2 footprints on glacier with elevation.The histogram shows the number of footprints of ICESAT-2 at different altitudes，and the curve is the distribution of SRTM with altitude

2.2.3 SRTM數(shù)字高程模型

SRTM（Shuttle Radar Topography Mission）是由美國、德國、意大利等航天機構(gòu)合作于2000年2月實施航天飛機雷達地形測繪任務(wù)，其攜帶了兩套合成孔徑雷達系統(tǒng)，分別采用C波段（波長5.6 cm，2.2 GHz）和X波段（波長3.1 cm，8.8 GHz）（Zhang等，2019）。由于X波段在中緯度地區(qū)存在許多空洞（Farr等，2007），對馬蘭山覆蓋缺失，文中使用SRTM-C DEM作為馬蘭山2000年時的冰面高程，并采用SRTM-X DEM對C波段進行冰雪穿透深度校正以獲得更加精確的冰面高程。SRTM-C DEM以EGM96坐標系為垂直空間參考，SRTM-X DEM以WGS84坐標系為垂直空間參考，二者均以WGS 84坐標系為水平參考，空間分辨率約為30 m×30 m，高程精度約為10 m。

2.2.4 Landsat-7陸地資源衛(wèi)星

Landsat-7是美國的陸地衛(wèi)星計劃（Landsat）中的第7顆，于1999年4月15日在加利福尼亞范登堡空軍基地用DeltaⅡ火箭發(fā)射。衛(wèi)星攜帶增強型專題制圖儀ETM+（Enhanced Thematic Mapper）傳感器，共有8個波段，除了長波紅外與全色波段，其余波段分辨率為30 m（Barsi等，2016）。本文所用的數(shù)據(jù)為經(jīng)過輻射校正和幾何校正的2級產(chǎn)品，時間為2000年、2007年—2012年冬季，用于提躍動取冰川末端輪廓線。根據(jù)不同年份的假彩色合成影像，手工勾畫躍動冰川末端輪廓線，并以不同輪廓線間沿主流線方向的平均距離估算冰川末端前進距離。

2.2.5 RGI冰川編目與氣象再分析資料

全球陸地冰川太空測量計劃（GLIMS）起源于ASTER科學團隊的一個項目，主要使用光學衛(wèi)星儀器的數(shù)據(jù)來監(jiān)測全球冰川，獲取冰川的圖像（RGI Consortium，2017）。RGI是GLIMS數(shù)據(jù)庫的補充數(shù)據(jù)集，是由IPCC推動產(chǎn)生的全球冰川輪廓編目（RGI Consortium，2017）。采用GLIMS于2017年7月發(fā)布的6.0版本，來確定冰川的邊界以及冰川相關(guān)信息的提取。

全球降水氣候中心（GPCC）氣象再分析資料是基于全球67200個記錄時間超過10年的氣象站的質(zhì)量控制數(shù)據(jù)生成的，采用V2018版本每月降水數(shù)據(jù)集，空間分辨率為1°×1°，獲得研究區(qū)2000年—2020年降水量變化。馬蘭山地區(qū)同期氣溫變化通過全球歷史氣候網(wǎng)絡(luò)（GHCN）每月地表氣溫再分析資料獲得，采用結(jié)合氣候異常監(jiān)測系統(tǒng)數(shù)據(jù)的GHCN_CAMS版本，空間分辨率為0.5°×0.5°。

3 方法

3.1 雙站干涉測量提取冰川表面高程變化

雙站干涉測量方法用于提取2000年—2012年馬蘭山冰川表面高程變化，并生成2012年3月的冰川表面DEM。該方法源自于合成孔徑雷達干涉測量二軌差分法，由于TSX與TDX對冰川表面的同時觀測，時間基線為零，其干涉相位之中不存在由冰川運動產(chǎn)生的形變相位，將其與SRTM DEM反演的地形相位進行差分（Liu等，2019），即可獲得從SRTM采集時間到TanDEM-X采集時間內(nèi)的冰川表面形變相位，進而計算冰川表面高程變化。主要步驟如下（蔣宗立等，2019）：

（1）地形相位模擬與配準。將SRTM-C DEM作為參考DEM模擬地形相位并利用互相關(guān)系數(shù)配準到TSX/TDM SAR坐標系。

（2）差分干涉與相位解纏。將差分干涉得到的形變相位進行濾波并以0.3為相干性閾值掩膜掉低相干性區(qū)域后，解纏后得到真實形變相位轉(zhuǎn)換為高程差。

（3）誤差糾正與DEM生成。利用無冰區(qū)地形平坦區(qū)域的點雙線性插值面對基線引起的系統(tǒng)誤差予以整體改正，將優(yōu)化的高程差與SRTM疊加得到TDX DEM。

馬蘭山冰川被兩對跨軌干涉數(shù)據(jù)對覆蓋（圖1），左右兩對干涉數(shù)據(jù)對的間隔時間為11天（表1），期間冰面高程變化可忽略不計。

3.2 基于ICESat-2與TanDEM估算冰川表面高程變化

采用ICESat-2激光雷達高度計對2012年—2020年、2000年—2020年冰面高程變化進行估算。ICESat-2足跡點在冰川區(qū)空間分布均勻（圖2），海拔覆蓋范圍為4900—5800 m，僅缺失海拔在5800—6000 m的冰川區(qū)域，其面積為2.04 km2，僅占馬蘭山冰川總面積的1.07%（圖3），所以采用該激光雷達高度計對冰川高程變化監(jiān)測是可行的。

將ICESat-2與SRTM、TDX DEM進行配準，采用雙線性內(nèi)插將DEM值提取至ICESat-2足跡點，進行高程差值計算，剔除差值超過±150 m的點。然后將冰川按高程劃分為相等間隔的分組，在每個分組內(nèi)將高程差的中值作為冰川表面高程變化值，按照高程分組的面積權(quán)重計算整個冰川表面的高程變化（Wang等，2017b）。

為了減少在低海拔與高海拔高程分組內(nèi)足跡點數(shù)量分布較少帶來的誤差，分別采用100 m、200 m、300 m的高程間隔分組對冰川高程變化進行重復計算，最終采用3次計算結(jié)果的平均值作為最終冰川高程變化值?？紤]到足跡點空間分布的均勻性會對估算結(jié)果產(chǎn)生比較大的影響，只對馬蘭山區(qū)域整體和面積大于2 km2的12條主要冰川進行高程變化監(jiān)測，篩選排除ICESat-2足跡點缺失高程分段的面積超過10%的冰川，最終共對9條冰川高程變化進行了統(tǒng)計。

3.3 不同高程數(shù)據(jù)配準

將SRTM-C DEM轉(zhuǎn)換到WGS84坐標系后以其為基準分別與TanDEM、ICESat-2進行配準。根據(jù)兩個不同高程數(shù)據(jù)間的偏差與坡度、坡向間存在的三角函數(shù)關(guān)系進行配準（Nuth和K??b，2011）：

式中，dh為高程偏差，α為坡度，ψ為坡向，a、b、c為3個參數(shù)，其中a的振幅為偏移矢量的大小，b是偏移矢量方向，c是兩DEM間的平均偏差除以平均坡度的正切值，在無冰區(qū)利用最小二乘法可以對其進行求解。

利用高程偏差的5%和95%分位數(shù)剔除異常值后，通過不斷迭代直到標準偏差小于2%或者位移量小于0.5 m。

3.4 偏差校正

3.4.1 季節(jié)校正

TSX/TDX干涉數(shù)據(jù)對的采集時間為3月，ICESat-2數(shù)據(jù)采集時間為11月到4月，跨越了冬季和春季，為了更準確的估計冰川物質(zhì)平衡，與SRTM的采集時間（2月）相對應(yīng)，需要對季節(jié)性冰川物質(zhì)平衡變化進行校正。馬蘭山區(qū)域在冰川物質(zhì)平衡線海拔高度附近的年降雨量（降雪量）為400—500 mm，其中90%的降雪都處于6月份到9月份（謝自楚等，2000），因此馬蘭山冰川屬于夏季積累型冰川。降雪、升華、再凍結(jié)和雪、冰融化決定了冰川物質(zhì)平衡狀態(tài)（Maussion等，2014）。由于主要的消融與降雪都集中在夏季，在寒冷干燥的非融雪季節(jié)（包括冬季和春季），由升華導致物質(zhì)虧損基本被少量降雪所彌補，少量的融雪被再凍結(jié)過程所抵消（Zhang等，2020），馬蘭山冰川在冬季和春季沒有明顯的季節(jié)性物質(zhì)平衡變化，所以季節(jié)校正值設(shè)定為0。

3.4.2 穿透校正

SRTM DEM由C（5.6 cm）波段雷達生成，由于C波段對冰蓋穿透可達1—2 m，對干雪的穿透最高可達10 m，會對2000年冰川表面高程產(chǎn)生一定的低估。本研究利用SRTM-X與SRTM-C的差值來估算C波段的穿透深度，以100 m高程間隔分段計算穿透深度，得到隨海拔變化的穿透深度差曲線，然后利用該曲線對冰川區(qū)域的每個像素進行校正（Gardelle等，2013）。由于SRTM-X在中緯度地區(qū)存在大量空洞，并未覆蓋馬蘭山區(qū)域，為了估算冰川區(qū)全部高程范圍的穿透深度，減少冰川間不同積雪覆蓋間引起的穿透深度差異，將東昆侖山全部冰川穿透深度的平均值作為馬蘭山冰川的穿透深度。

以X波段與C波段的穿透深度差異來估算C波段的穿透深度，并沒有考慮X波段對雪的穿透，所以經(jīng)過校正的冰面高程實際為X波段DEM冰面高程。TanDEM-X與SRTM-X波長相同，穿透效應(yīng)一致，由于冬季馬蘭山區(qū)域幾乎沒有降雪，積雪覆蓋在2000年和2012年冬季可以認為幾乎相同，因此經(jīng)過穿透校正的SRTM-C與TanDEM-X在對冰川表面的高程的低估上可以抵消，所以對2000年—2012年冰面高程變化的結(jié)果無需考慮X波段的穿透。在利用ICESat-2激光雷達高度計估算2000年—2020年與2012年—2020年冰面高程變化時，ICESat-2為綠光（532 nm）不會對冰雪表面產(chǎn)生任何穿透效應(yīng)，而利用SRTM-X DEM校正后的冰面高程仍然包含X波段相對于激光測高的穿透深度，因此這會導致2000年與2012年的冰面高程相對于激光雷達測高產(chǎn)生一定低估，從而導致物質(zhì)平衡估算結(jié)果偏高。目前還沒有有效方法對X波段的穿透深度進行精確評估，根據(jù)Dehecq等（2016）的研究，X波段對干雪的穿透范圍在2—6 m，平均穿透深度約為4 m。基于X波段只對干雪有穿透這一假設(shè)，Liu等（2019）基于冬季與春季的積雪覆蓋差異，利用2012年4月份與1月的份TanDEM-X相減得到了普若崗日冰川X波段平均穿透深度為0.61±0.06 m；Zhou等（2019）以平均4 m雪厚和40%的冬季干雪覆蓋率估算了昆侖山中心區(qū)域的X波段穿透約為1.5 m。由于根據(jù)氣候條件對積雪覆蓋的假設(shè)并不準確，X波段的穿透是一個明顯高估的值，所以沒有將其考慮進物質(zhì)平衡估算結(jié)果。

除研究李杜外，聞一多研究的唐代詩人還有李商隱和韓愈，甚至構(gòu)想了“唐代六大詩人之研究”，可惜他的宏大研究計劃至死未竟。

3.5 冰川物質(zhì)平衡計算

由于缺乏對馬蘭山區(qū)域冰川近表層密度的實測數(shù)據(jù)，根據(jù)Huss（2013）的研究結(jié)果，采用850±60 kg/m3作為冰川體積物質(zhì)轉(zhuǎn)換的密度，其中60 kg/m-3作為冰川物質(zhì)平衡估算結(jié)果誤差進行計算。最終，冰川物質(zhì)平衡的不確定度采用以下公式進行計算（Pieczonka和Bolch，2015）：

式中，um為不確定度，Δh為冰川的表面高程變化值，ρI為冰的密度（850 kg/m3），Δρ為冰密度的不確定度（60 kg/m3），ρw為水的密度，t為監(jiān)測時間跨度，E為高程精度。

3.6 精度驗證

根據(jù)配準之后的結(jié)果，利用非冰川區(qū)高程差值對精度進行驗證。為了排除DEM像元間的自相關(guān)，根據(jù)Bolch等（2011）等定義的方式，結(jié)合研究區(qū)地形，按照500 m為間隔對無冰區(qū)裸地進行抽樣統(tǒng)計。對最終冰川表面高程變化的不確定度（e）采用下式進行評估：

式中，SE（Standard Error of the Mean）為標準平均誤差，STDVnoglac為無冰區(qū)標準差，N為無冰區(qū)采樣點個數(shù)，MED（Mean Elevation Difference）為無冰區(qū)高程殘差均值。

ICESat-2的精度驗證采用同樣的方式，考慮到軌道間存在自相關(guān)，對不同軌道之間按照5 km間隔進行篩選，同軌的足跡點仍按照500 m間隔進行抽樣統(tǒng)計，去除絕對誤差在3倍標準差范圍外的點，最終高程變化不確定度采用下式進行評估：

式中，ebias為陸冰算法得到的足跡點傳播誤差均值，傳播誤差包括一階光子偏差和采樣誤差。

4 結(jié)果

4.1 穿透深度

利用SRTM-X與SRTM-C相減估算的穿透深度如圖4所示，在東昆侖山區(qū)域，隨著海拔的升高，C波段對冰雪的穿透深度逐漸增加，呈波動上升趨勢，在4900—6000 m范圍內(nèi)，其穿透深度最小為1.32 m（4900—5000 m），最大為5.14 m（5600—5700 m），平均穿透深度為3.91±0.86 m，略高于Zhou等（2019）的穿透深度估算結(jié)果（3.26 m），這是由于所用方法的差異造成的，并利用馬蘭山東部臨近的五雪峰冰川的穿透深度通過線性擬合外推來估算馬蘭山區(qū)域的穿透深度，此方法無法準確估算高海拔地區(qū)的穿透深度。本文為了得到全部海拔范圍的穿透深度，利用東昆侖山的全部冰川估算，由于馬蘭山位于東昆侖山中心區(qū)域最南端，相同海拔更高緯度地區(qū)冰川積雪覆蓋更厚，覆蓋面積更廣，所以可能會對結(jié)果產(chǎn)生一定的高估。

圖4 X波段與C波段SRTM穿透深度差異（100 m高程分段，誤差棒為標準差）Fig.4 Penetration depth difference between SRTM-X and SRTM-C at each 100 m altitude band.The error bars indicate the standard deviation of the estimated penetration depth difference

4.2 冰川表面高程變化與物質(zhì)平衡

通過雙站差分得到的2000年—2012年冰川表面高程變化如圖5所示。從圖5中可以看出，除躍動冰川（Ⅲ）外，馬蘭山地區(qū)的所有冰川在末端都有不同程度的減薄，在頂端會有一些增厚區(qū)域的出現(xiàn)。12條主要冰川（表2）均呈負物質(zhì)平衡狀態(tài)，但其物質(zhì)虧損速率也存在一定的空間差異性，其中Ⅶ和Ⅸ兩條冰川在12年間減薄超過7 m，物質(zhì)平衡為-0.52±0.04 m·w·e/a，是整個研究區(qū)消融最劇烈，冰量損失最嚴重的冰川，其余冰川減薄都在2.5—6.0 m范圍內(nèi)。

圖5 2000年—2012年馬蘭山冰川表面高程變化Fig.5 Thickness changes of glacier in Malan Mountain from 2000 to 2012

根據(jù)ICESat-2估算2012年—2020年冰川高程變化和物質(zhì)平衡結(jié)果如表2所示，可以發(fā)現(xiàn)在2012年之后，馬蘭山冰川整體仍呈強烈的消融狀態(tài)，區(qū)域內(nèi)所有冰川平均減薄2.08±1.07 m，物質(zhì)平衡為-0.22±0.11 m·w·e/a，其中Ⅺ與Ⅷ在8年間減薄最多，分別為-2.98±1.07 m和-2.20±1.07 m，其余冰川減薄則小于2 m。2012年—2020年馬蘭山9條主要冰川消融速度要低于2000年—2012年，但Ⅶ、Ⅸ在2012年前后物質(zhì)平衡差異較大，在2012年后消融速度驟減，Ⅻ甚至變?yōu)檎镔|(zhì)平衡，這主要是ICESat-2足跡點在這3條面積較小的冰川內(nèi)空間分布不均勻?qū)е碌模捎谧阚E點只覆蓋了這3條冰川西側(cè)，且分布于發(fā)生積累區(qū)域的足跡點數(shù)量要多于消融區(qū)，因此這3條冰川在2012年—2020年物質(zhì)平衡估算結(jié)果偏高。

表2 馬蘭山地區(qū)冰川物質(zhì)平衡Table 2 Glacier mass balance in Malan Mountain

根據(jù)2000年—2020年冰川物質(zhì)平衡的結(jié)果可以發(fā)現(xiàn)（圖6），冰川表面高程變化隨著海拔的升高波動上升，物質(zhì)虧損逐漸減少，在海拔高于5800 m的范圍內(nèi)已經(jīng)基本呈平衡或弱正平衡狀態(tài)。整體而言，2012年前后兩段時間內(nèi)冰面高程變化趨勢大體相同，都處于較強的消融狀態(tài)，物質(zhì)持續(xù)虧損，但根據(jù)表2中主要冰川物質(zhì)平衡估算結(jié)果與圖6中冰川區(qū)高程變化曲線，2012年—2020年除了在5100—5300 m海拔范圍內(nèi)冰川消融速度加快，整體冰川消融速率相較于2000年—2012年有所降低。

圖6 冰川區(qū)高程變化曲線（100 m高程分段）Fig.6 Glacier thickness change curve with an altitude bin of 100 m

5 討論

5.1 冰川對氣候變化的響應(yīng)

馬蘭山地區(qū)氣溫與降水變化趨勢如圖7所示。自2000年以來，氣溫的逐年升高（0.07℃/a）是冰川消融的主要因素，導致冰川表面整體高程降低。降水主要受到西風環(huán)流控制，同期內(nèi)降水量波動變化，總體呈上升趨勢。由于年降水的微弱增加并不足以彌補由氣溫升高產(chǎn)生消融導致的冰量損失，所以馬蘭山冰川在此20年間呈明顯負平衡狀態(tài)（-0.24±0.06 m·w·e/a），物質(zhì)虧損嚴重。根據(jù)Yao等（2012）對青藏高原冰川對氣候環(huán)流響應(yīng)的研究以及王盛等（2011）基于度日模型對同樣受到西風環(huán)流控制的祁連山七一冰川對氣候變化敏感性的研究結(jié)果，氣溫是影響冰川物質(zhì)平衡最主要的因素。馬蘭山地區(qū)2000年—2012年間年均氣溫（0.1℃/a）與夏季氣溫（0.03℃/a）都呈明顯上升趨勢，2012年—2020年間則升溫速度減緩（0.06℃/a），夏季氣溫甚至出現(xiàn)了微弱的波動下降趨勢（-0.03℃/a），而馬蘭山冰川的消融主要發(fā)生在夏季，這導致了消融速度的減緩。雖然夏季降水與年均降水量在2012年之后的增速要小于2012年之前，但根據(jù)GPCC降雨再分析數(shù)據(jù)可知2000年—2020年這20年間馬蘭山地區(qū)年均降水量僅為306.16 mm，較為稀少，所以降水小幅度變化以及增速差異對整個冰川補給效應(yīng)的影響不大，而且2000年—2012年均降水量為303.53 mm，2012年—2020年均降水量為314.4 mm，雖然2012年后波動上升趨勢減緩，但實際對冰川物質(zhì)的補給量要稍多。氣溫變化的差異主導了冰川的變化狀態(tài)，降雨量雖然存在一定波動變化差異，但總體對冰川補給影響較小，所以2012年之后冰川的物質(zhì)損失速度稍有減慢。

圖7 2000年—2020年馬蘭山地區(qū)氣溫和降水變化Fig.7 Temperature and precipitation changes in Malan Mountain from 2000 to 2020

5.2 躍動冰川

2000年—2012年馬蘭山地區(qū)所有冰川末端都呈現(xiàn)出不同程度的減薄，只有CN5Z122A0003冰川末端明顯增厚，異常的增厚導致了在5100—5300 m海拔范圍內(nèi)2012年后冰川消融速度快于2012年之前，這與整個馬蘭山冰川在2012年前后消融快慢的整體趨勢相反。該冰川末端異常增厚的原因是2007年—2012年發(fā)生了躍動（Zhang等，2020）。

如圖8所示，通過對比不同年份Landsat-7影像發(fā)現(xiàn)該冰川末端在2000年—2007年退縮了約47 m，處于正常的消融狀態(tài)，而在2007年—2012年前進了251 m；其中，在2007年—2009年前進了117 m，2009年—2010年前進了93 m，2010年—2012年前進了41 m，前進速度最快的時間在2009年—2010年間。

圖8 2000年—2012年CN5Z122A0003（躍動冰川）末端輪廓線（背景為2017年12月Landsat-8假彩色合成圖像）Fig.8 Glacier outlines between 2000 and 2012 for Glacier CN5Z122A0003（surge-type glacier）.The background is the Landsat 8 false-color image acquired in December 2017

根據(jù)冰川在2000年—2012年的平均消融速度，估算出冰川在2007年時表面高程，TanDEM則可以認為是冰川發(fā)生躍動后的冰面高程，其在躍動前后的高程變化如圖9所示。結(jié)果表明，躍動從5200—5250 m海拔范圍開始，積累區(qū)面積約為6.53 km2，平均表面高程下降9 m，最大下降29 m，體積約減少5.9×107m3，接收區(qū)面積約為1.94 km2，表面高程平均上升27 m，最大上升為59 m，體積約增加5.2×107m3，大量冰體通過冰崩、物質(zhì)重分配向冰川末端轉(zhuǎn)移，導致冰川末端異常增厚。

圖9 2007年—2012年躍動冰川沿主流線高程變化Fig.9 Elevation change of surge-type glacier along the main flow line between 2007 and 2012

6 結(jié)論

采用TerraSAR-X/TanDEM-X合成孔徑雷達、SRTM DEM、ICESat-2激光雷達高度計等，基于雙站差分干涉測量，對馬蘭山冰川表面高程變化和物質(zhì)平衡變化進行了分析，結(jié)果表明：2000年—2020年馬蘭山冰川平均表面高程變化量為-5.64±0.96 m，物質(zhì)平衡為-0.24±0.06 m·w·e/a，呈明顯負平衡趨勢，冰量損失嚴重。結(jié)合GPCC與GHCN_CAMS氣象再分析資料發(fā)現(xiàn)，冰川減薄主要是由于氣溫升高帶來的消融，降雨量的微弱增加僅彌補了少部分物質(zhì)損失，2012之后的升溫變慢，夏季均溫波動下降導致了冰川消融速度減緩。其中，位于馬蘭山南坡的CN5Z122A0003為躍動冰川，其在2007年—2012年發(fā)生躍動，導致冰川末端在此期間發(fā)生多次前進，根據(jù)Landsat-7影像研究發(fā)現(xiàn)，2009年—2010年間前進速度最快，在一年內(nèi)前進了93 m，然后逐漸減緩，冰川末端共前進約為251 m，大量冰體向下游轉(zhuǎn)移，導致下游的接收區(qū)發(fā)生異常增厚。

用于提取冰面高程變化的兩種測量方法之間存在一定的偏差，ICESat-2激光雷達高度計為不連續(xù)的點數(shù)據(jù)，而干涉測量生成的SRTM、TanDEM為30 m×30 m空間分辨率的柵格數(shù)據(jù)，僅以一個點周圍臨近像元線性插值的結(jié)果作為該點處的高程值，這種由面到點的轉(zhuǎn)換存在一定的誤差，尤其是在地形起伏較大、坡度較陡的區(qū)域。此外，由于激光雷達不會對冰雪產(chǎn)生穿透效應(yīng)，而X波段雷達對干雪存在一定的穿透，干涉測量的得到高程在積雪覆蓋較多較厚的高海拔區(qū)域會有輕微低估，一定程度上影響了物質(zhì)平衡估算結(jié)果的精度。