基于多源數據的天山烏魯木齊河源1號冰川變化研究

李宏亮， 王璞玉，， 李忠勤， 王盼盼， 徐春海， 劉爽爽，金 爽，張正勇，徐麗萍

（1.中國科學院西北生態環境資源研究院冰凍圈科學國家重點實驗室/天山冰川觀測試驗站，甘肅蘭州730000；2.中國科學院大學，北京100049；3.石河子大學理學院，新疆石河子832003；4.西北師范大學地理與環境科學學院，甘肅蘭州730070）

0 引言

冰川是氣候變化的敏感指示器，也是全球水循環的重要組成部分［1］。在氣候變暖的影響下，全球大多數冰川呈現普遍變薄和退縮趨勢，對水資源、水循環和生態環境等都產生了顯著影響［2-3］，導致冰崩、冰湖潰決等災變風險的增加［4-6］。過去50年間全球氣溫升高，我國天山冰川面積和冰儲量均呈持續退縮趨勢，其面積和冰儲量年均退縮率分別為每年-0.7%和每年0.83%，且物質損失嚴重，減少了27%，同時1942—2014年間末端年均退縮速率為1.56 m·a-1［7-9］，短期內可使河川徑流量增加，但在長時間尺度上這種影響可能使河川徑流總量減少［10］，在冰川面積減少、厚度變薄及平衡線海拔升高等因素和機制的促進下，天山中部典型流域徑流量自20世紀90年代中期后遞減［11］，體現了氣候變化驅動下以冰川為代表的固態水體對流域水資源的調控機制。因此，開展天山地區冰川變化與氣候變化的聯系至關重要，對提高天山冰川變化的認識和水資源利用與管理具有重要意義。

基于各冰川參數的變化可評估氣候變化及管理水資源［12-13］。冰川物質平衡是表征冰川積累和消融量值的重要冰川參數之一，對氣候變化響應敏感。物質平衡及其動態變化是引起冰川規模和徑流變化的物質基礎，是連結冰川與氣候、冰川與水資源的重要紐帶。面積和長度（末端）是冰川的重要幾何形態參數，在冰川消融模擬、冰川幾何形態變化模擬、冰川水文等各類研究中，均充當不可或缺的基礎參數和邊界限制條件。通過定期對參照冰川的面積、長度（末端）等參數進行重復觀測，不僅可以精細研究冰川變化局部特征，捕捉變化過程，為冰川對氣候變化響應及其模擬研究提供數據基礎，而且是遙感資料不可或缺的驗證數據。

目前國內外在青藏高原及周邊地區開展了大量關于冰川變化時空格局及原因和對水資源及海平面變化的影響［14-20］等方面的研究。但仍存在一些問題。如：第一，大尺度冰川變化研究結果精度存在顯著差異且缺少實測驗證資料；第二，利用單條冰川觀測結果可進行實測驗證，但連續且詳細的監測冰川在全球范圍內非常少。因此，對現有定位觀測的參照冰川持續觀測尤為重要。烏魯木齊河源1號冰川（簡稱1號冰川）作為全球定位觀測的參照冰川之一，前人已經開展了大量研究［21-26］，隨著近年來各種高新技術手段的涌現，使得冰川觀測精度大幅提高，如三維激光掃描（Terrestrial Laser Scanner，TLS）靈活且更經濟，可獲得單條冰川年或季節內高時空分辨率冰面高程信息［27］；無人機航測技術（Unmanned Aerial Vehicle，UAV）具有機動性強、操作簡單，低空飛行可獲取冰面高分辨率正射影像（DOM）和數字高程模型（DEM），不僅能夠清晰解讀冰面微地貌特征，還在冰川地表模型構建方面具有明顯的優勢［28］；載波相位差分定位技術（Real Time Kinematic-Global Position System，RTK-GPS）是目前GPS測量中精度最高的一種定位方法，可實現密集點位覆蓋，適合于冰川表面地形高程測量等優勢［29］。前期觀測1號冰川主要利用RTK-GPS和TLS技術，其中，物質平衡觀測結果能夠與花桿/雪坑法得到的物質平衡很好的對應，但觀測結果僅由單一技術獲取的冰面高程得到［24-25］，結合RTK-GPS和TLS技術，本文首次引入UAV技術對1號冰川繼續觀測，同時試驗RTKGPS、TLS和UAV技術之間應用在冰川變化研究中的適用性。因此，本研究將基于不同時期RTKGPS、TLS和UAV技術獲取的冰川高程數據，分析1號冰川近期變化，比較三種觀測技術之間應用到冰川變化的適用性，同時提供近期準確的冰川變化信息，為下一步分布式冰川能量-物質模擬提供基礎。

1 研究區概況

1號冰川（43°06′N，86°49′E）是我國監測時間最長、觀測資料連續性最好的一條冰川，位于我國天山中部喀拉烏成山脈主脈北坡烏魯木齊河河源上游［圖1（a）］，屬雙支冰斗-山谷冰川，由東西兩支組成，呈東北走向［圖1（b）］。20世紀50年代至今，對該冰川進行連續的物質平衡、末端位置、表面高程等觀測，2012年冰川厚度為44.5 m，1980—2012年間冰川厚度減薄速率為0.34 m·a-1，同期末端退縮速率及物質平衡分別為4.4 m·a-1和-0.46 m w.e.·a-1，東西支末端變化差異顯著，在過去32年間呈加速退縮趨勢［23］。

2 數據與方法

2.1 冰川表面高程、面積和末端變化

為了開展近期1號冰川表面高程、面積和末端變化，本研究結合不同時期數據源開展了對比研究（表1）。

表1 本研究使用數據列表Table 1 Data information used in this study

2012年9月1日利用RTK-GPS獲得1號冰川表面高程和末端位置數據。GPS基站固定在臨近末端的位置［圖1（b）］，接收機在冰面進行同步觀測，測量間距為20 m，除冰川上部陡峭部分外，實現了整個冰川的觀測。

圖1 烏魯木齊河源1號冰川觀測示意圖Fig.1 The observation diagram on the Urumqi Glacier No.1：geographical location of the Urumqi Glacier No.1（a），observation sites（b）

2015年4月25日利用三維激光掃描技術獲取1號冰川表面點云數據，其中用于確定目標對象三維坐標信息的縱向（?）和橫向（φ）掃描角度可在地面激光掃描儀中完成，最終，獲得目標對象的三維坐標信息。同時，為了使三維激光掃描儀能夠最大可能的掃描到1號冰川全貌，將脈沖頻率分別設置為粗掃描（30 kHz）與精掃描（50 kHz），結合RTK-GPS確定的4個掃描站坐標信息［圖1（b）］，實現了點云配準以及重復掃描區域比例大于30%的要求［21］，最后進行坐標系統配準和點云數據分類及濾波。

結合RTK-GPS及三維激光掃描儀獲取的冰面高程數據，筆者于2018年4月24日使用大疆經緯M 200專業型四旋翼無人機對1號冰川進行了航空攝影測量以獲取冰川正射影像和數字表面模型（DSM）。但由于電池供電不足原因，此次航測只獲取到冰舌區冰面高程，具體航測范圍如圖1（b）所示，同時，在冰川基巖處和冰面花桿位置均勻的布設了8個地面控制點（GCPs）。解算軟件使用Pix4D mapper，得到高精度的DEM及正射影像。GCPs精度利用平均誤差（Emean）與均方根誤差（RMSE）來評估［30］，結果表明，無人機航測平均水平誤差為0.04 m，均方根誤差（RMSE）為0.06 m。

將上述不同時期的冰川觀測數據進行重采樣（5 m×5 m），并進行坐標歸一化處理，均采用統一的UTM投影和WGS84橢球體坐標系統。之后，通過不同時期的數據對比，開展冰川表面高程、面積和末端變化分析研究。

2.2 花桿/雪坑法冰川物質平衡

1號冰川物質平衡觀測始于1959年，每年4月底和8月末各觀測1次，研究期間均勻布設26～43個花桿，開展花桿/雪坑法冰川物質平衡觀測，以期和不同時期冰面高程數據變化，即大地測量法冰川物質平衡結果進行對比。測量包括研究時期內花桿高度、雪坑深、雪類型及密度。單點物質平衡計算方法見式（1），之后利用插值手段，結合等值線法與等高線法，得到觀測時段內整條冰川的物質平衡值。

式中：bs、bice、bsi分別為雪、冰川冰以及附加冰的物質平衡，具體可參見已有研究［31］。

花桿/雪坑法得到的冰川物質平衡誤差主要來源于野外觀測和計算整條冰川物質平衡時的插值，后者主要是將插值方法應用到冰川上未觀測的區域中［32］，由于缺少實測數據，本文依據Andreassen等［33］假設的插值造成誤差為±0.1 m w.e.·a-1作為本文插值造成的誤差。考慮以上因素，整條冰川的物質平衡年誤差（σglac）可以利用式（2）計算。

2.3 誤差分析

2.3.1 DEM配準及表面高程誤差

在計算多源DEM差值時，需要對原始高程數據進行配準和校正。本文采用Nuth等［34］提出的數據配準方法，以2015年DEM為參考，其余時期DEM為錯位DEM，選擇非冰川區作為感興趣區，提取兩期DEM非冰川區的高程差值dh、2015年DEM的坡度和坡向，通過余弦擬合求得平移矢量來校正DEM間水平和垂直誤差。在完成DEM數據配準后，由于空間分辨率的差異，不同DEM數據間還存在高程偏差，且這種偏差在冰川區和非冰川區保持一致，利用地面最大曲率與非冰川區高程差之間的關系進行糾正［35］，本研究基于2015年DEM數據，在ENVI 5.1軟件支持下計算像元最大曲率，建立了最大曲率與高程差的關系。

根據不同DEM數據間的高程殘差誤差滿足高斯分布的假設，表面高程誤差采用非冰川區高程殘差的均方根誤差或方差進行評估。DEM數據格網點的高程值存在較強的空間自相關，因此，對樣本數據選擇時，必須去除空間自相關，本文采用156 m去相關距離，基于Bolch等［36］研究DEM數據的誤差可用高程殘差平均值和標準誤差計算：

式中：N為像元數；SD與SE分別為標準差和標準誤差；MED為高程殘差平均值；E為誤差（表2）。

表2 感興趣區DEM校正前后誤差Table 2 The error of DEMs before and after correction in interesting area

2.3.2 冰川面積和末端變化誤差

基于RTK-GPS、TLS及UAV數據的冰川邊界，獲得1號冰川末端和面積變化，利用式（5）和式（6）進行冰川面積誤差評估［37-39］：

式中：UT與UA分別為末端及面積不確定誤差；λ為空間分辨率；ε為矯正誤差，保證在一個像元內。計算得到，2012—2015年和2015—2018年面積誤差分別為1.4×10-4km2和7×10-6km2。

2.3.3 大地測量法物質平衡估算

冰川物質平衡估算采用大地測量法，大地測量法物質平衡（Bgeo）是通過將冰川體積變化與平均密度相乘，然后除以面積得出的，該面積是兩期面積的平均值。Huss等［40］認為冰川密度是變化的且小于冰的密度，建議使用（850±60）kg·m-3較為合適。因此，本文采用（850±60）kg·m-3作為冰川體積-物質平衡轉換參數。

3 結果與分析

3.1 冰川表面高程及物質平衡變化

1號冰川表面高程在低海拔區域呈現明顯冰量損失，而在高海拔區域呈輕微減薄或變厚趨勢（圖2）。考慮到2018年無人機航測僅獲得1號冰川冰舌區表面高程，所以，2012—2018年間表面高程下降（1.34±0.88）m·a-1，對應物質平衡為（-1.13±0.18）m w.e.·a-1，累積物質虧損量（-6.78±1.08）m w.e.，東支冰量損失較西支明顯，而2015—2018年間冰川表面高程下降（2.03±0.96）m·a-1，物質平衡（-1.72±0.19）m w.e.·a-1，累 積 物 質虧 損 量（-5.16±0.57）m w.e.。2012—2015年間冰川表面高程下降（0.83±0.57）m·a-1，物質平衡（-0.71±0.17）m w.e.·a-1，累積物質虧損量（-2.13±0.51）m w.e.（表3）。研究期間表面高程減小區域主要位于冰舌區及東支右側，且東西兩支表面高程減薄速率隨海拔高度的增加而減小。2012—2015年表面高程在-6～2 m間變化，冰川上部呈積累狀態［圖2（a）］；2015—2018年冰舌區表面高程變化區間為-13～1 m，東支減薄速率大于西支［圖2（b）］，所以，2012—2018年冰舌區表面高程整體在-20～1 m間變化，東支末端減薄速率大于西支，海拔高度越高表面高程減薄越小，在西支表面高程變化上更為顯著［圖2（c）］。

圖2 烏魯木齊河源1號冰川表面高程變化Fig.2 Changes of surface elevation of the Urumqi Glacier No.1

表3 烏魯木齊河源1號冰川不同時期表面高程及物質平衡Table 3 Surface elevation and mass balance of the Urumqi Glacier No.1 in different periods

3.2 大地測量法與花桿/雪坑法物質平衡對比

Wang等［24］計算得到1981—2009年間大地測量法與花桿/雪坑法所得物質平衡的差異小于10%（表4）；隨著觀測儀器的發展和誤差評估方法的改進，Xu等［25］基于Zemp等［12］提出的關于大地測量法和花桿/雪坑法間的物質平衡減小差異（δ）的計算方法，得到了1981—2015年1號冰川物質平衡的減小差異 值（δ=0.53）通 過95%的置 信區 間檢 驗（|δ|＜1.96），大地測量法和花桿/雪坑法分別得到物質平衡值具有一致性。基于此，本文利用上述方法分析2012—2018年間兩種方法的差異性，其中，2012—2015年間1號冰川物質平衡的差異［（Bgeo-Bgla）/Bgla］為2.9%（表4），兩種方法的減小差異值（δ=0.78）通過95%的置信區間，這與前人研究結果一致。

表4 烏魯木齊河源1號冰川不同時期花桿/雪坑法與大地測量法累積物質平衡比較Table 4 Comparison of glaciological and geodetic cumulative mass balances of the Urumqi Glacier No.1 for given time periods

在整條冰川上利用大地測量法與花桿/雪坑法所得物質平衡的差異（δ）進行評估，而2018年4月無人機航測區域僅為1號冰川冰舌區域，為了對2012—2018年及2015—2018年間大地測量法與花桿/雪坑法獲取的物質平衡進行對比，提取對應花桿點的物質平衡。結果顯示，大地測量法和花桿/雪坑法的物質平衡兩者相關系數（R2）分別為0.93和0.91，均方根誤差（RMSE）為1.19 m w.e.和0.81 m w.e.（圖3），表明兩種方法得到的物質平衡能夠較好的對應。

圖3 烏魯木齊河源1號冰川冰舌區單點大地測量法與花桿/雪坑法物質平衡對比Fig.3 Geodetic versus glaciological mass balances for the ice tongue of Urumqi Glacier No.1

結合1號冰川物質平衡長期觀測結果發現，自1959年以來1號冰川年物質平衡和累積物質平衡呈減少趨勢，物質損失明顯［26］，基于線性統計方法按照斜率不同得到1號冰川物質平衡2個時期的變化特征：1985—1996年物質平衡為-0.27 m w.e.·a-1，物質損失顯著大于1959—1984年（-0.08 m w.e.·a-1）；從1997年開始，物質虧損更為強烈，使得1997—2018年間年物質平衡為-0.68 m w.e.·a-1，其中，2010年冰川物質平衡值低至-1.33 m w.e［42］。總的來說，2012—2018年間1號冰川物質損失（-0.64 m w.e.·a-1）大于1980—2012年（-0.47 m w.e.·a-1），2012—2014年物質平衡虧損減緩，而2015—2018年間物質平衡虧損增加，表明近期1號冰川物質平衡虧損仍在持續。將1號冰川年物質平衡與全球41條參照冰川物質平衡結果相比，其平均值變化趨勢相一致［42］，說明1號冰川物質平衡在一定程度上能代表全球山地冰川平均物質平衡變化。

3.3 冰川面積及末端變化

自1959年以來，基于野外觀測發現1號冰川退縮強烈［23］。考慮到2018年僅得到1號冰川3 785～4 185 m間冰川區，同時認為冰川上部面積保持不變的條件下，利用2015年4 185～4 485 m間冰川區得到2018年完整冰川區用于估算總面積。2012年，1號冰川總面積為1.59 km2，到2015年減少為1.56 km2，到2018年減少至1.52 km2，2012—2018年1號冰川面積整體減少了0.07 km2，年均變化率為-0.01 km2·a-1。其 中，2012—2015年 和2015—2018年均變化率分別為-0.01 km2·a-1和-0.02 km2·a-1，后者面積退縮強度大于前者，總的來說，2012—2018年1號冰川整體上呈現出持續退縮趨勢。

2012—2018年1號冰川末端變化退縮速率為6.28 m·a-1，東、西支退縮速率分別為7.64 m·a-1和4.93 m·a-1，其中，2012—2015年末端變化基本處于穩定狀態，在2016年東西支變化趨勢相反，2017—2018年退縮幅度顯著，但整體上東支退縮速率顯著大于西支（圖4）。在相同氣候背景下，末端變化差異不僅取決于物質平衡的變化，而且與冰川地形和熱力學參數有關的冰川動力學過程有關。首先，2012—2018年東支表面高程減薄速率大于西支［圖2（c）］，這與東西支末端退縮趨勢一致，體現了物質平衡對末端退縮的影響；其次，由于西支冰舌區坡向朝向為東南、東支受山體遮掩及西支面積和橫截面比東支小等因素的影響，在假設消融速率相同的情景下，導致2000年代之前西支退縮速率大于東支。此外，已知冰川運動速度的年際變化作用于末端位置，已有研究表明1號冰川2012年后東支運動速度顯著小于西支，低于1號冰川年均運動速度［41］，導致西支末端海拔較東支越來越高，且與西支相比東支末端受山體遮掩影響變小。最終，除2017年外西支退縮速率小于東支退縮速率。因此，通過分析物質平衡、面積及末端變化表明烏魯木齊河源1號冰川近期呈加速退縮趨勢。

圖4 1980—2018年烏魯木齊河源1號冰川末端變化Fig.4 Terminus changes of the Urumqi Glacier No.1 during 1980—2018

4 結論

本文基于RTK-GPS、TLS和UAV等資料分析了2012—2018年烏魯木齊河源1號冰川面積、末端和物質平衡變化，對比了近期烏魯木齊河源1號冰川變化，結果表明：

（1）2012—2018年期間，烏魯木齊河源1號冰川面積減少0.07 km2，年平均面積退縮率為-0.01 km2·a-1；烏魯木齊河源1號冰川末端在2012—2018年間呈退縮趨勢，變化率為6.28 m·a-1，東支退縮速率顯著大于西支，與2012—2015年相比，2015—2018年末端變化退縮幅度更為顯著。

（2）烏魯木齊河源1號冰川物質平衡在2012—2018年間表面高程下降（1.34±0.88）m·a-1，物質平衡為（-1.13±0.18）m w.e.·a-1，物質損失明顯，其中，2012—2015年間冰川表面高程下降（0.83±0.57）m·a-1，物質平衡為（-0.71±0.17）m w.e.·a-1；2015—2018年間冰川表面高程下 降（2.03±0.96）m·a-1，物質平衡為（-1.72±0.19）m w.e.·a-1，物質損失主要發生在消融區及東支右側，2012—2018年間物質平衡變化率大于1980—2012年間物質平衡變化率表明烏魯木齊河源1號冰川近期呈加速消融。

致謝：感謝西北研究院劉宇碩工程師、西北師范大學沈思民碩士研究生和蘭州大學鄭續碩士研究生在無人機等數據處理過程中的大力幫助！