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

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

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

0 引言

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

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

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

1 研究區(qū)概況

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

2 數(shù)據(jù)與方法

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

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

表1 本研究使用數(shù)據(jù)列表Table 1 Data information used in this study

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

圖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日利用三維激光掃描技術(shù)獲取1號冰川表面點(diǎn)云數(shù)據(jù)，其中用于確定目標(biāo)對象三維坐標(biāo)信息的縱向（?）和橫向（φ）掃描角度可在地面激光掃描儀中完成，最終，獲得目標(biāo)對象的三維坐標(biāo)信息。同時(shí)，為了使三維激光掃描儀能夠最大可能的掃描到1號冰川全貌，將脈沖頻率分別設(shè)置為粗掃描（30 kHz）與精掃描（50 kHz），結(jié)合RTK-GPS確定的4個(gè)掃描站坐標(biāo)信息［圖1（b）］，實(shí)現(xiàn)了點(diǎn)云配準(zhǔn)以及重復(fù)掃描區(qū)域比例大于30%的要求［21］，最后進(jìn)行坐標(biāo)系統(tǒng)配準(zhǔn)和點(diǎn)云數(shù)據(jù)分類及濾波。

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

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

2.2 花桿/雪坑法冰川物質(zhì)平衡

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

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

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

2.3 誤差分析

2.3.1 DEM配準(zhǔn)及表面高程誤差

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

根據(jù)不同DEM數(shù)據(jù)間的高程殘差誤差滿足高斯分布的假設(shè)，表面高程誤差采用非冰川區(qū)高程殘差的均方根誤差或方差進(jìn)行評估。DEM數(shù)據(jù)格網(wǎng)點(diǎn)的高程值存在較強(qiáng)的空間自相關(guān)，因此，對樣本數(shù)據(jù)選擇時(shí)，必須去除空間自相關(guān)，本文采用156 m去相關(guān)距離，基于Bolch等［36］研究DEM數(shù)據(jù)的誤差可用高程殘差平均值和標(biāo)準(zhǔn)誤差計(jì)算：

式中：N為像元數(shù)；SD與SE分別為標(biāo)準(zhǔn)差和標(biāo)準(zhǔn)誤差；MED為高程殘差平均值；E為誤差（表2）。

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

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

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

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

2.3.3 大地測量法物質(zhì)平衡估算

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

3 結(jié)果與分析

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

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

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

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

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

Wang等［24］計(jì)算得到1981—2009年間大地測量法與花桿/雪坑法所得物質(zhì)平衡的差異小于10%（表4）；隨著觀測儀器的發(fā)展和誤差評估方法的改進(jìn)，Xu等［25］基于Zemp等［12］提出的關(guān)于大地測量法和花桿/雪坑法間的物質(zhì)平衡減小差異（δ）的計(jì)算方法，得到了1981—2015年1號冰川物質(zhì)平衡的減小差異 值（δ=0.53）通 過95%的置 信區(qū) 間檢 驗(yàn)（|δ|＜1.96），大地測量法和花桿/雪坑法分別得到物質(zhì)平衡值具有一致性?；诖耍疚睦蒙鲜龇椒ǚ治?012—2018年間兩種方法的差異性，其中，2012—2015年間1號冰川物質(zhì)平衡的差異［（Bgeo-Bgla）/Bgla］為2.9%（表4），兩種方法的減小差異值（δ=0.78）通過95%的置信區(qū)間，這與前人研究結(jié)果一致。

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

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

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

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

3.3 冰川面積及末端變化

自1959年以來，基于野外觀測發(fā)現(xiàn)1號冰川退縮強(qiáng)烈［23］?？紤]到2018年僅得到1號冰川3 785～4 185 m間冰川區(qū)，同時(shí)認(rèn)為冰川上部面積保持不變的條件下，利用2015年4 185～4 485 m間冰川區(qū)得到2018年完整冰川區(qū)用于估算總面積。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，后者面積退縮強(qiáng)度大于前者，總的來說，2012—2018年1號冰川整體上呈現(xiàn)出持續(xù)退縮趨勢。

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

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

4 結(jié)論

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

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

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