聯(lián)合雙站InSAR與WoldView DEM 估算祁連山西段冰川多時段物質(zhì)平衡

李超，江利明，柳林，李濤，陳圓圓

1.中國科學院精密測量科學與技術(shù)創(chuàng)新研究院,武漢 430071;2.中國科學院大學 地球與行星科學學院,北京 100049;3.華中科技大學 物理學院引力實驗中心,武漢 430074;4.中國地質(zhì)大學 地理與信息工程學院,武漢 430074

1 引言

山地冰川是氣候變化的指示器，也是重要的淡水資源，近年來冰川消融對海平面上升的貢獻也越來越大（Houghton等，2001；Jacob等，2012）。隨著全球氣溫的不斷升高，高亞洲地區(qū)大多數(shù)冰川出現(xiàn)末端退縮和面積減小現(xiàn)象，冰川物質(zhì)虧損加?。˙olch 等，2012；Cogley，2016；Shean 等，2020；Bhattacharya 等，2021）。在此背景下，地處西北內(nèi)陸干旱區(qū)的祁連山冰川區(qū)域，同樣也經(jīng)歷了不同程度的冰川退縮（Yao 等，2012；張其兵等，2017；高永鵬 等，2018）。此外，冰川融水是中國西北干旱區(qū)重要的水資源，對社會經(jīng)濟可持續(xù)發(fā)展具有重要意義，但同時也會導致冰湖潰決、泥石流等災(zāi)害頻發(fā)。因此，加強祁連山冰川變化監(jiān)測以及冰川物質(zhì)平衡的估算，對區(qū)域水資源管理、生態(tài)環(huán)境保護和冰川災(zāi)害防治具有重要意義（謝自楚和劉潮海，2010；李新 等，2020）。

目前冰川物質(zhì)平衡的估計方法主要包括冰川學方法、積累區(qū)面積比率法、水文學法、重力測量法和大地測量方法等（Gardelle 等，2013）。受限于氣候和地形條件，冰川學方法實施困難且無法獲取大范圍冰川物質(zhì)平衡（楊大慶 等，1992）。重力測量法雖然可以通過重力衛(wèi)星觀測數(shù)據(jù)直接反演冰川物質(zhì)平衡，但無法獲取高空間分辨率的物質(zhì)平衡（Gardner等，2013；Matsuo和Heki，2010）。積累區(qū)面積比率法需要高精度的冰川物質(zhì)平衡和準確的冰川參數(shù)作為支撐。水文學法需要長時間序列的降水、蒸發(fā)、徑流等數(shù)據(jù)，其參數(shù)獲取較為困難。大地測量法主要是通過使用重復的、全覆蓋的數(shù)字高程模型（DEM）來比較冰面高程的變化，獲得高時空分辨率的冰川物質(zhì)平衡方法。本文特指采用DEM 差分法獲得冰厚及冰量變化的方法。近年來，隨著衛(wèi)星遙感技術(shù)的快速發(fā)展以及遙感衛(wèi)星的增多（周玉杉 等，2019；周建民 等，2021），高分辨率、高精度DEM 的日益增多，大地測量方法已成為獲取大范圍冰川物質(zhì)平衡的主流方法（K??b 等，2015；Neckel 等，2014；Toutin等，2011）。

位于祁連山西段的老虎溝12 號冰川是祁連山最大的冰川，為當?shù)亟?jīng)濟發(fā)展提供了重要的水資源，許多學者在此展開研究（Du 等，2016；Liu等，2019b；Qin等，2015；Sun等，2014；Wang等，2018；Zhang 等，2012）。然而，由于復雜的山地地形和惡劣的氣候條件，實測冰川物質(zhì)平衡記錄仍然很少，目前相關(guān)研究大多基于激光測高、SAR數(shù)據(jù)、重力數(shù)據(jù)和光學遙感數(shù)據(jù)。王玉哲等（2013）基于ICESat 數(shù)據(jù)利用差分方法估算了冰川高程隨時間的變化趨勢，進而獲得了2003 年—2009 年祁連山冰量變化。結(jié)果表明21 世紀初祁連山冰川處于物質(zhì)虧損狀態(tài)，但ICESat 數(shù)據(jù)在冰川上的光斑較少導致該結(jié)果有很大的不確定性。Gardner等（2013）基于ICESat 數(shù)據(jù)研究表明2003 年—2009 年整個祁連山的年平均高程變化-0.32±0.31 m/a。孫美平等（2015）基于兩次冰川編目數(shù)據(jù)對祁連山冰川變化進行了分析，研究表明近50 年間祁連山冰川面積和冰儲量分別減少約420.81 km2（-20.88%）和21.63 km3（-20.26%），且面積較小的冰川急劇萎縮是該區(qū)冰川面積減少的主要原因。隨著TanDEM-X 衛(wèi)星的成功發(fā)射，星載雙站SAR干涉測量開始應(yīng)用于冰厚變化監(jiān)測。目前TanDEM已成功應(yīng)用于青藏高原冰川區(qū)冰厚變化和物質(zhì)平衡的估算（Lin等，2017；Zhang等，2020）。Sun等（2018）基于TanDEM-X 雙站InSAR 數(shù)據(jù)，采用DEM 差分方法獲取了祁連山西段2000 年—2013 年的冰川表面高程變化結(jié)果，并與GPS 數(shù)據(jù)進行了對比分析，結(jié)果表明祁連山冰川處于消融狀態(tài)。高永鵬等（2018）基于校正后的SRTM數(shù)據(jù)與ASTER立體像對數(shù)據(jù)，利用大地測量法，對2000 年—2010 年祁連山地區(qū)冰川冰儲量變化進行應(yīng)用研究，結(jié)果表明2000 年—2010 年祁連山冰川厚度平均減?。?.68±2.76）m，平均冰川物質(zhì)平衡為（-0.48±0.23）m w.e./a，冰川儲量變化（-1.59±0.72）Gt。綜上所述，這些前期的工作一般側(cè)重于監(jiān)測冰川年代際的物質(zhì)平衡，缺乏年際物質(zhì)平衡研究，這在一定程度上限制了我們對冰川狀態(tài)的時間變化和冰川對區(qū)域氣候變化響應(yīng)的認識。

本研究旨在利用2013 年的TanDEM-X 雙站干涉合成孔徑雷達生成的DEM 數(shù)據(jù)（TDX DEM）以及2014 年—2015 年的高空間分辨率WorldView 光學DEM（WV DEM）進行差分處理獲得相應(yīng)時間段的冰厚變化，并結(jié)合2000 年的航天飛機雷達地形測繪任務(wù)（SRTM）DEM獲得了2000年—2015年祁連山西段的冰川高程變化。為了獲得高精度的冰川物質(zhì)平衡結(jié)果，本文充分考慮C波段和X波段雷達的穿透性和DEM 數(shù)據(jù)之間的季節(jié)性差異，獲得較為可靠的祁連山西段以及老虎溝12 號冰川2013 年—2015 年的年際冰川物質(zhì)平衡。最后，結(jié)合溫度和降水數(shù)據(jù)分析老虎溝12 號冰川物質(zhì)平衡變化的影響因素。

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

2.1 研究區(qū)概況

祁連山位于中國青海省東北部與甘肅省西部邊境，發(fā)育有冰川2859 條，冰川總面積1972.5 km2，估算儲水量為811.2×108m3。老虎溝12 號冰川位于祁連山西段北坡的大雪山地區(qū)，由東西兩支冰川組成，是祁連山區(qū)最大的極大陸型冷性復式山谷冰川，冰川表面較為平緩（坡度3°—6°）且積累區(qū)較寬，具有穩(wěn)定性冰川的特征，其冰川融水是昌馬河重要的徑流補給來源（劉宇碩 等，2013；于國斌 等，2014）。

祁連山冰川是中國西北干旱地區(qū)的“天然固體水庫”，也是西北山區(qū)內(nèi)各條河流的重要補給源頭。近年來，隨著全球不斷變暖，祁連山區(qū)的溫度也呈現(xiàn)出明顯上升趨勢，特別是20 世紀90 年代以后，氣溫上升幅度明顯增大。祁連山地區(qū)的年降水量也呈增加趨勢，其西段降水量增幅相對東段來說較大。祁連山區(qū)屬高原大陸性氣候，西北部主要受西風帶氣流控制，年平均氣溫5 ℃左右，降水主要集中在夏季，年降水量由東向西逐漸減少。該地區(qū)冰川融水是中國西北綠洲城市農(nóng)業(yè)灌溉和社會經(jīng)濟發(fā)展的重要水資源通道，因此研究該區(qū)域冰川物質(zhì)平衡對于甘肅省、青海省乃至全國都意義重大（張其兵 等，2017）。

2.2 數(shù)據(jù)

本研究主要采用3 種類型的DEM 數(shù)據(jù)進行處理：（1）一對TanDEM-X 雙站SAR 數(shù)據(jù)生成高分辨率高精度的 TDX DEM；（2）WorldView 光學影像生成的8 m 分辨率的DEM；（3）30 m 分辨率的SRTM DEM。在數(shù)據(jù)處理的過程中我們還相應(yīng)的采用RGI V6.0 冰川編目數(shù)據(jù)和植被積雪覆蓋數(shù)據(jù)。數(shù)據(jù)具體參數(shù)如下表1 所示，圖像覆蓋范圍如圖1所示。

表1 所用數(shù)據(jù)主要參數(shù)Table 1 Main parameters of DEM data

2.2.1 TDX DEM數(shù)據(jù)

TanDEM-X 是德國宇航局（DLR）的對地觀測雷達任務(wù)，由兩顆成像參數(shù)相同的X 波段雷達衛(wèi)星（TerraSAR-X、TanDEM-X）以近距離螺旋結(jié)構(gòu)緊密排列方式組成的雙衛(wèi)星系統(tǒng)進行觀測，消除了時間去相干和大氣去相干等去相干源（Breit等，2010）。本研究中采用的高分辨率高精度的TDX DEM是由2013年11月21日獲取的一對TanDEM-X雙站SAR 降軌數(shù)據(jù)生成，空間分辨率為10 m（Sun等，2018）。

2.2.2 SRTM DEM 數(shù)據(jù)

NASA航天飛機雷達地形任務(wù)（SRTM）數(shù)據(jù)集是由美國國家航空航天局（NASA），國家地理空間情報局（NGA）和德國航空航天中心（DLR）聯(lián)合開發(fā)的，其目的是使用雷達干涉測量法生成全球數(shù)字高程模型（DEM）。2000年2月11日至22日在STS-99 任務(wù)期間獲取了大約80%的陸地總面積（60°S—60°N）的雷達干涉數(shù)據(jù)（Neckel等，2013）。SRTM DEM 主要采用兩個波段傳感器分別為C 波段（波長5.6 cm）的SIR-C 傳感器和X 波段（波長3.1 cm）的X-SAR 傳感器，然而由于X 波段數(shù)據(jù)的條帶寬度較窄，僅呈相隔50 km 寬度交叉覆蓋，因此X 波段DEM 不連續(xù)，無法覆蓋整個研究區(qū)域（Farr等，2007）。

本文使用的是空洞填充后的1 弧秒C 波段和X波段SRTM DEM，可以分別從USGS 網(wǎng)站（https：//dds.cr.usgs.gov［2022-04-18］）和DLR 網(wǎng)站（https：//www.dlr.de［2022-04-18］）免費下載。C 波段和X波段SRTM DEM 具有相同的水平參考（WGS84 基準）；但二者的垂直參考基準不同（C 波段為EGM96大地水準面，X波段為WGS84橢球體），需在DEM 配準之前進行高程基準轉(zhuǎn)換。此外，兩個DEM 在數(shù)據(jù)質(zhì)量方面具有相同的水平精度（絕對誤差和相對誤差分別為20 m和15 m）。由于X波段數(shù)據(jù)的分辨率和信噪比稍高，因此X 波段數(shù)據(jù)的垂直精度（絕對精度16 m，相對精度6 m）優(yōu)于C波段數(shù)據(jù)（絕對精度16 m，相對精度10 m）（Berthier等，2006）。

2.2.3 光學DEM數(shù)據(jù)

本文采用美國冰雪數(shù)據(jù)中心提供的光學DEM數(shù)據(jù)，該數(shù)據(jù)主要是通過WorldView 衛(wèi)星獲取的Level-1B（L1B）全色（450—800 nm）影像，采用自動化的開源NASA艾姆斯立體聲管道方法（ASP）生成的數(shù)字高程模型（DEM）。主要參數(shù)如表1 所示，單個DEM 條帶內(nèi)的相對誤差通常小于1—2 m（Shean等，2016）。

2.2.4 冰川編目和氣象數(shù)據(jù)

本文使用RGI 6.0（Randolph Glacier Inventory，RGI）全球冰川編目數(shù)據(jù)確定冰川邊界。RGI 是在世界冰川編目（WGMS）的基礎(chǔ)上編制而成，并已在多個大規(guī)模的全球評估和建模研究中應(yīng)用（Vaughan 等，2013）。該數(shù)據(jù)庫覆蓋高亞洲地區(qū)的95537 個冰川，總冰川面積為97605 km2。研究區(qū)內(nèi)冰川覆蓋面積為93.8 km2。

中國區(qū)域地面氣象要素驅(qū)動數(shù)據(jù)集，是高時空分辨率近地表氣象數(shù)據(jù)集。該數(shù)據(jù)集是以國際上現(xiàn)有的再分析資料和衛(wèi)星遙感數(shù)據(jù)為背景場，融合了中國氣象局常規(guī)氣象觀測數(shù)據(jù)制作而成。其時空分辨率均較高，時間分辨率為3 h，水平空間分辨率為0.1°，精度介于氣象局觀測數(shù)據(jù)和衛(wèi)星遙感數(shù)據(jù)之間，可為中國區(qū)陸面過程模擬提供驅(qū)動數(shù)據(jù)。本文提取2010 年—2015 年祁連山區(qū)域日際溫度和降水數(shù)據(jù)，并分析其與冰川物質(zhì)平衡的聯(lián)系。

2.2.5 積雪和植被覆蓋數(shù)據(jù)

為獲得穩(wěn)定的基巖區(qū)域我們需要對非冰川區(qū)域進行冰川積雪和植被覆蓋區(qū)域的去除。積雪范圍數(shù)據(jù)使用的是全球8 d 500 m SIN 網(wǎng)格的MODIS/Terra 積雪覆蓋數(shù)據(jù)，此數(shù)據(jù)可以免費從國家冰雪數(shù)據(jù)中心（https：//nsidc.org/data/MOD10A2/versions/6［2022-04-18］）獲得。該積雪數(shù)據(jù)集主要是通過歸一化差異降雪指數(shù)（NDSI）獲得，此指數(shù)由8天的500 m 分辨率MOD10A1 瓦片生成，因積雪數(shù)據(jù)的有效時間是從2000 年2 月24 日開始，所以我們提取了2000年2月24日至3月24日的積雪范圍用于去除積雪區(qū)域?qū)RTM C/X 波段DEM 之間的垂直偏差的影響。

植被覆蓋數(shù)據(jù)主要采用的是中國向聯(lián)合國提供的首個全球地理信息公共產(chǎn)品。GlobeLand30 是一種30 m 分辨率的全球土地覆蓋數(shù)據(jù)產(chǎn)品，是通過使用各種多光譜圖像來開發(fā)和更新的WGS-84坐標系下的數(shù)據(jù)。已發(fā)布產(chǎn)品共包括10 個土地覆蓋類別，分別是：耕地、林地、草地、灌木地、濕地、水體、苔原、人造地表、裸地、冰川和永久積雪。在這項研究中，我們僅使用了3種類型的數(shù)據(jù)：耕地，人造地表和裸地（Jun等，2014）。

3 研究方法

隨著InSAR 技術(shù)和光學立體測繪技術(shù)的快速發(fā)展，多時相、高精度的DEM 日益增多。特別是DLR 的TanDEM-X 雙站InSAR 獲得的高分辨率、高精度DEM 和基于WorldView、SPOT 等光學影像的高精度DEM。這些高質(zhì)量、大范圍的DEM 數(shù)據(jù)產(chǎn)品，為利用多時相DEM 估算冰川物質(zhì)平衡提供了寶貴的基礎(chǔ)資料。本文主要是利用多時相DEM估算冰川物質(zhì)平衡，其數(shù)據(jù)處理流程主要包括以下幾個方面：首先是數(shù)據(jù)預(yù)處理，統(tǒng)一多源DEM的空間分辨率和坐標系統(tǒng)等（陳昊楠 等，2020）；其次是進行多時相DEM 配準與差值處理，與此同時進行奇異值去除（高程變化結(jié)果中的粗差）；然后是系統(tǒng)誤差的改正主要包括雷達穿透性改正和季節(jié)性改正；最后統(tǒng)計冰川研究區(qū)高程變化的平均值并將高程變化轉(zhuǎn)換到物質(zhì)平衡。具體實現(xiàn)流程如圖2所示。

圖2 數(shù)據(jù)處理技術(shù)路線Fig.2 Flowchart of data processing

3.1 DEM預(yù)處理

不同數(shù)據(jù)來源的DEM 往往會具有不同的時空分辨率、坐標系統(tǒng)以及高程基準等，例如SRTMX DEM 的空間分辨率為30 m，高程基準為WGS84橢球面，而SRTM-C DEM 的空間分辨率有30 m 和90 m 兩種，高程基準為EGM96 大地水準面。此外，近年來隨著星載光學攝影測量的不斷發(fā)展，光學影像的數(shù)據(jù)質(zhì)量也大有不同，利用光學影像生成的DEM 空間分辨率也不盡相同，所以在利用多時相DEM 估算冰川物質(zhì)平衡之前，需要對其進行預(yù)處理，基本處理步驟包括：（1）坐標系統(tǒng)和高程基準的轉(zhuǎn)換，將研究區(qū)的所有DEM 統(tǒng)一到一致的坐標系統(tǒng)和高程基準下，一般采用WGS84 坐標系和橢球面；（2）DEM 重采樣，將空間分辨率較低的DEM 進行插值，得到較高空間分辨率的DEM，使得所有DEM 具有一致的空間分辨率。本文首先將SRTM DEM 數(shù)據(jù)的坐標系統(tǒng)轉(zhuǎn)到統(tǒng)一的WGS84 系統(tǒng)，其次利用三次樣條插值法進行DEM重采樣，以WordView DEM 為基準，將所用到的DEM數(shù)據(jù)統(tǒng)一采樣到8 m分辨率。

Paul（2008）指出地形曲率較大的區(qū)域DEM重采樣會引入高程誤差。圖3為30 m SRTM 重采樣到8 m 的高程誤差與曲率的相關(guān)性散點圖，可以看出SRTM DEM的曲率大多分布在-0.01到0.02 m-1，相對應(yīng)的重采樣前后的SRTM DEM 差值在-0.1—0.2 m，兩者存在顯著的線性關(guān)系。本文利用獲得的線性關(guān)系方程式（y=9.617x-0.04）對高程誤差進行了改正，高程誤差改正平均值約為0.005 m。另外，TanDEM-X DEM 與WorldView DEM 分辨率相近，故未進行重采樣引起的高程誤差改正。

圖3 SRTM 30 m與8 m DEM 重采樣高程差與曲率的相關(guān)性Fig.3 Relation between curvature and the elevation differences computed from the original 30 m SRTM DEM and the one resampled from 8 m

3.2 DEM配準方法

采用Nuth 和K??b（2011）提出的DEM 高程差和坡度以及坡面朝向之間的關(guān)系，進行配準，獲取DEM之間的偏移總量和偏移總方向。公式如下：

式中，dh為不同DEM 數(shù)據(jù)間的高程差異；ψ和α分別為像元所對應(yīng)的坡向和坡度；為不同DEM數(shù)據(jù)間的整體高程差異。平面偏移會導致地面坡度產(chǎn)生偏差，所以我們需要對不同DEM 數(shù)據(jù)的高程差異基于坡度做歸一化處理，進而獲得高程差異與坡向之間的三角函數(shù)關(guān)系，公式如下所示：

式中，參數(shù)a，b和c可以通過回歸分析獲取，c代表垂直偏移量；a和b分別代表平面偏移的幅度和角度。不同DEM 數(shù)據(jù)間X，Y和Z方向上的偏移量校正公式如下：

式中，a 和b分別代表平面偏移幅度和角度，為DEM數(shù)據(jù)的平均坡度。

3.3 系統(tǒng)誤差改正

3.3.1 SRTM C波段穿透性改正

對于中低緯度地區(qū)的山地冰川而言，雷達信號對冰雪的穿透深度是SAR DEM 估算冰川物質(zhì)平衡的一個重要的系統(tǒng)誤差（Gardelle等，2012；Lin等，2017；Zhou 等，2018）。雖然X 波段雷達信號對雪或冰有一定的穿透能力，但是在先前青藏高原山地冰川的研究中基本被忽略，通常將冰川區(qū)SRTM X 波段和C 波段的高程差作為C 波段雷達信號在冰川區(qū)的穿透深度（Gardelle 等，2012）。然而，在冬季X 波段雷達信號對冰川也有一定的穿透性，其穿透深度是不可忽略的，所以本文將SRTM C 波段穿透性改正分為兩部分，一是SRTM C 波段與X 波段的穿透性深度差異，二是X 波段的穿透深度。

SRTM C 波段與X 波段的穿透性深度差異我們采用的是DEM 差值法。首先我們將高亞洲地區(qū)按1°×1°網(wǎng)格進行劃分，然后通過SRTM C 波段和X波段的高程做差，獲得該地區(qū)SRTM C 波段和X 波段的穿透性深度差異與高程的相關(guān)關(guān)系，最后利用此關(guān)系式求得整個冰川區(qū)的C波段和X波段的穿透性深度差異。主要步驟如下：（1）數(shù)據(jù)的預(yù)處理，將SRTM C 波段的高程系統(tǒng)由EGM96 轉(zhuǎn)為WGS84，實現(xiàn)兩個波段DEM 的高程基準的統(tǒng)一。（2）SRTM C 波段與X 波段DEM 進行配準，消除空間偏差后做差。（3）利用基巖區(qū)SRTM C 波段與SRTM X 波段高程差為零這一原則得到非冰川區(qū)的高程偏差，并從冰川區(qū)高程差中去除，從而消除高程偏差的影響。值得注意的是，此處基巖區(qū)指的是去除植被和積雪區(qū)的非冰川區(qū)。（4）在冰川區(qū)利用50 m 高程分段求取冰川區(qū)SRTM C 波段和X波段的穿透性深度差異與高程的關(guān)系。

因SRTM X 波段DEM 與TDX DEM 均為X 波段數(shù)據(jù)，所以在進行SRTM X 波段穿透深度改正時采用的是TanDEM-X 數(shù)據(jù)的穿透深度。因此，本文采用Liu 等（2019a）利用1 月和4 月的TanDEM-X數(shù)據(jù)獲得的普若崗日地區(qū)的X 波段的穿透性的平均值（0.61 m）進行X波段穿透深度改正。

3.3.2 季節(jié)性改正

采用Wang 等（2017）利用ICESat 數(shù)據(jù)得到的2003 年—2009 年的季節(jié)性改正數(shù)據(jù)的平均值進行季節(jié)性改正。因2013 年—2014 年的數(shù)據(jù)均在12 月不需要改正，只需將2015 年08 月的數(shù)據(jù)改正到12 月份，對獲得的季節(jié)性改正數(shù)據(jù)按月平均，獲得最后的改正結(jié)果為-0.23±0.10 m。

3.4 冰川高程變化平均值估算

在數(shù)據(jù)預(yù)處理、DEM 配準以及各項系統(tǒng)性誤差改正之后需要對獲取的冰川表面高程變化進行統(tǒng)計分析，為估算冰川物質(zhì)平衡提供基礎(chǔ)數(shù)據(jù)。但在進行統(tǒng)計工作之前，采用兩倍中誤差方法剔除冰川表面高程變化結(jié)果中的奇異值，以提高平均值估算的精度和準確性。為了抑制隨機誤差的影響我們采用50 m 高程分段統(tǒng)計冰川高程變化平均值，就是將冰川區(qū)的DEM 分割成多個連續(xù)的高程段（如：4100—4150 m、4150—4200 m和4200—4250 m 等）并且假定單個高程段內(nèi)（如：4100—4150 m）的像元具有相似的高程變化結(jié)果，然后分別計算各個高程段的冰川平均高程變化，最后根據(jù)各個高程段的冰川面積，采用面積加權(quán)方法（該高程段像素個數(shù)與總像素個數(shù)的比值），統(tǒng)計冰川區(qū)的高程變化平均值（式4）。

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

由于高亞洲地區(qū)氣候條件惡劣，人力很難到達，所以無法獲得DEM 獲取時間段內(nèi)實測的冰川區(qū)密度。在此背景下，通常采用密度模型估算的冰川密度進行冰川高程變化到冰川物質(zhì)平衡的轉(zhuǎn)換（Liu等，2020）。具體公式如下：

式中，dhpene和dhseason分別為穿透性改正和季節(jié)性改正，ρ為冰雪密度，在本研究中采用Huss 2013 年得到的適用于多種條件的冰川密度850±60 kg m-3（Huss，2013）。

3.6 不確定性分析

利用多時相DEM 估算冰川物質(zhì)平衡，其誤差來源主要包括以下幾個方面：（1）高程分段統(tǒng)計冰川高程變化均值的誤差；（2）C波段穿透性差異改正引起的誤差；（3）冰川季節(jié)性變化改正引起的誤差。本文最終采用下式（6）評估冰川物質(zhì)平衡不確定性：

式中，σ△h，σρ，分別為平均冰川高程變化值誤差、冰川密度誤差、雷達波段穿透差異誤差、季節(jié)性改正誤差。其中，σρ為冰川密度值的7%，即±60 kg/m3（Huss，2013），為季節(jié)性改正的數(shù)值，通過高程分段法計算的穿透深度利用誤差傳播定律獲得，σ△h計算公式如下：

式中，在某一高程區(qū)間內(nèi)，Ni為像元個數(shù)，為第i個高程段的高程變化均值誤差其計算公式如下：

式中，STDi為第i個高程段均值標準差，PS是所用DEM的像素分辨率；D是空間自相關(guān)的距離。Bolch認為D約為像素大小的20倍（Bolch等，2012）。

4 結(jié)果與討論

4.1 DEM配準前后結(jié)果分析

本文通過對2013年—2014年和2014年—2015年的DEM 配準前后高程變化結(jié)果的對比來評判配準方法的可行性以及可靠性。如圖4所示。

圖4 祁連山冰川配準前后高程變化，其中直方圖表示非冰川冰厚變化Fig.4 Elevation changes before and after co-registration of Qilian Mountains Glacier where the histogram shows the elevation changes in non-glacial area

由圖4 得：（1）因采用數(shù)據(jù)的時間間隔為1 年左右，故高程變化不會太大，配準之后的高程變化值的范圍減小，平均值減小，并向0 附近靠攏，表明配準后的結(jié)果精度更高。

（2）配準后的標準差整體趨勢是減小的，數(shù)據(jù)的離散程度減小，配準后的結(jié)果更可靠。

4.2 SRTM C波段穿透性深度

因為老虎溝12號冰川地區(qū)無SRTM X波段覆蓋，所以采用臨近1°×1°網(wǎng)格區(qū)域獲得的穿透性深度與高程的關(guān)系式進行外推獲得。采用臨近網(wǎng)格N39E095 的穿透深度差異數(shù)據(jù)，此網(wǎng)格SRTMC/X的穿透深度差異與高程呈現(xiàn)線性相關(guān)關(guān)系，方程式為Y=0.00409X-19.81（詳見Li 等，2021），具體如圖5所示。

4.3 祁連山西段冰川高程變化結(jié)果分析

（1）2013 年—2014 年冰厚變化。2013 年—2014年老虎溝12號冰川高程變化信號明顯（圖6），冰川末端高程變化為負，表明末端處于消融狀態(tài)且越靠近末端消融越快，而冰川上游高程變化為正處于積累狀態(tài)。祁連山西段冰川的高程變化均值為-0.35±0.034 m，老虎溝12 號冰川的高程變化均值為-0.32±0.024 m（表2），物質(zhì)平衡為-0.33±0.04 m w.e./a，表明2013 年—2014 年老虎溝12 號冰川處于虧損狀態(tài)，這與劉宇碩等實測的2014 年老虎溝12 號物質(zhì)平衡數(shù)據(jù)（-0.49 m w.e./a）相近（劉宇碩，2018）。

圖6 祁連山西段冰川2013年—2014年高程變化結(jié)果Fig.6 Glacier elevation changes in the West Qilian Mountains from 2013 to 2014

表2 2013年—2014年—2015年祁連山西段冰川和老虎溝12號冰川配準后統(tǒng)計結(jié)果（兩倍中誤差）Table 2 Statistical results after co-registration of Qilian Mountains glaciers and Laohugou No.12 glaciers from 2013 to 2014 to 2015（The error is twice）

（2）2014 年—2015 年冰厚變化。由圖7 可以看出祁連山西段地區(qū)的冰川高程變化在±1 m 波動，明顯比2013年—2014年的高程變化幅度要小，圖中祁連山西段冰川的高程變化均值為-0.03±0.004 m，老虎溝12 號冰川的高程變化均值在-0.04 m 左右（表2），從圖7 中可以明顯看出冰川末端區(qū)域處于消融狀態(tài)。我們以老虎溝12 號冰川為例，發(fā)現(xiàn)該冰川在2014 年—2015 年的高程變化的平均值與2013 年—2014 年相比要小，表明冰川厚度減薄變緩。

（3）2000 年—2015 年冰厚變化。本文通過對WorldView DEM 與SRTM 進行差值，經(jīng)過穿透性和季節(jié)性改正，最終獲得祁連山西段高程變化結(jié)果，如下圖8所示。從圖8中可以看出2000年—2015年高程變化集中在-15—15 m，老虎溝12號冰川的高程變化均值為-0.15 m/a，整體上冰川處于消融狀態(tài)。由于冰川末端消融較大，為了使結(jié)果更精確，本文未對2000 年—2015 年的冰川高程變化進行兩部中誤差去除。圖8明顯可以看出冰川末端處于消融狀態(tài)且越靠近末端消融越快。

圖8 祁連山西段冰川2000年—2015年高程變化結(jié)果Fig.8 Glacier elevation changes in the West Qilian Mountains from 2000 to 2015

4.4 祁連山西段冰川物質(zhì)平衡結(jié)果分析

根據(jù)大地測量法獲得祁連山西段冰川2013年—2014 年—2015 年的高程變化，結(jié)合冰雪密度進而獲得其冰川物質(zhì)平衡。結(jié)果表明祁連山西段在2013 年—2014 年和2014 年—2015 年的物質(zhì)平衡變化速率分別為-0.27±0.014 m w.e./a 和-0.024±0.084 m w.e./a，2000 年—2015 年的物 質(zhì)平衡為-0.098±0.01 m w.e./a。另外，老虎溝12號冰川的物質(zhì)平衡由2013 年—2014 年的-0.33±0.04 m w.e./a 減小到2014 年—2015年的-0.036±0.09 m w.e./a，其2000年—2015年的物質(zhì)平衡為-0.13±0.02 m w.e./a?？梢钥闯稣w上祁連山西段處于負物質(zhì)平衡狀態(tài)2014 年—2015 年物質(zhì)平衡相對2013 年—2014 年有所增加，物質(zhì)趨于平衡，尤其是老虎溝12號冰川。相比2013 年—2014 年結(jié)果，2014 年—2015 年冰川物質(zhì)平衡的不確定性較大，主要是由于季節(jié)性改正誤差的影響。由圖9、圖10和表3可以看出，在2015 年祁連山地區(qū)降水量顯著增加且溫度與前幾年相比變化不大（圖9、圖10），所以祁連山西段冰川在2014 年—2015 年物質(zhì)平衡有所增加的主要原因為降水量增加。綜合以上分析結(jié)果，祁連山冰川物質(zhì)平衡處于穩(wěn)定且略有虧損的階段，2014 年—2015 年物質(zhì)平衡相對2013 年—2014 年略有增加的主要原因為降水量增加。

圖9 2010年—2015年祁連山地區(qū)月均降水量（每年01-01—12-31）Fig.9 Average monthly precipitation in the Qilian Mountains from 2010 to 2015（January 1 to December 31 each year）

圖10 2010年—2015年祁連山地區(qū)溫度（每年01-01—12-31）Fig.10 The temperature in the Qilian Mountains from 2010 to 2015（from January 1 to December 31 each year）

表3 2010年—2015年祁連山地區(qū)溫度和降水統(tǒng)計Table 3 Statistics of temperature and precipitation in Qilian Mountains from 2010 to 2015

4.5 祁連山西段冰川高程變化與海拔之間的關(guān)系

同一區(qū)域，冰川整體上高程變化是類似的，但依然存在局部空間差異，往往與冰川海拔、冰川地形坡度、坡向和冰川動力等因素有關(guān)。本研究對老虎溝冰川高程變化與海拔之間的關(guān)系進行了剖面線分析，剖面線位置如圖11 綠線所示。由圖11 可以看出，隨著海拔的升高，冰川高程變化由負轉(zhuǎn)正，冰川高程變化結(jié)果與海拔之間呈負正相關(guān)關(guān)系，表明海拔是冰川表面高程變化的一種重要因素。

圖11 老虎溝12號冰川中心線剖面圖Fig.11 Centerline profile of Laohugou No.12 Glacier

5 結(jié)論

本文利用大地測量法，結(jié)合多時相高精度DEM 數(shù)據(jù)，充分考慮SRTM C波段的穿透深度和季節(jié)性改正兩類系統(tǒng)誤差的影響，獲得了祁連山西段川2000年—2013年和2013年—2015年的冰厚變化以及冰川物質(zhì)平衡，尤其是老虎溝12 號冰川的2000 年—2013 年和2013 年—2015 年多時段的冰川物質(zhì)平衡結(jié)果。

結(jié)果表明：2013 年—2015 年祁連山西段冰川冰厚減薄速率變緩，冰厚變化結(jié)果由2013 年—2014 年的-0.35±0.024 m 變?yōu)?014 年—2015 年的-0.03±0.004 m，物質(zhì)平衡變化速率也由-0.27±0.014 m w.e./a減緩至-0.024±0.084 m w.e./a，主要原因是2015年降水量的增加。其中就老虎溝12號冰川而言，冰川消融呈現(xiàn)相同的變化規(guī)律，即2013年—2015 年冰川消融速率降低。老虎溝12 號冰川的2013 年—2014 年和2014 年—2015 年的物質(zhì)平衡分別為-0.33±0.04 m w.e./a和-0.036±0.09 m w.e./a。結(jié)合2000 年—2015 年的冰川物質(zhì)平衡結(jié)果表明祁連山西段冰川尤其是老虎溝12 號冰川在2000 年—2015 年整體處于消融狀態(tài)。另外本文通過分析冰川中線位置處的高程變化，發(fā)現(xiàn)冰川高程變化與海拔之間呈正相關(guān)關(guān)系，隨著海拔的升高，冰川高程變化越小。

本研究初步驗證了光學立體測繪衛(wèi)星DEM 數(shù)據(jù)在求解冰川年際物質(zhì)平衡方面的可行性，但受限于數(shù)據(jù)獲取問題，難以獲得長時間序列的冰川年際變化信息。隨著衛(wèi)星遙感技術(shù)的不斷發(fā)展，高精度DEM 數(shù)據(jù)將在長時間序列的冰川物質(zhì)平衡估算中發(fā)揮更大的作用。

志 謝感謝美國地質(zhì)調(diào)查局（USGS）提供的SRTM-C DEM，德國航空航天中心（DLR）的提供SRTM-X DEM 和美國冰雪數(shù)據(jù)中心提供的World-View DEM 數(shù)據(jù)；感謝時空三級環(huán)境大數(shù)據(jù)平臺提供的地面氣象要素驅(qū)動數(shù)據(jù)；感謝匿名審稿人的寶貴意見和建議。