都偉冰， 張世瓊， 李均力， 包安明， 王雙亭， 史寧可，許琳娟， 高 鑫， 馬丹丹， 鄭巖超

（1.河南理工大學(xué)測繪與國土信息工程學(xué)院，河南焦作 454003；2.中國科學(xué)院新疆生態(tài)與地理研究所荒漠與綠洲生態(tài)國家重點(diǎn)實(shí)驗(yàn)室，新疆烏魯木齊 830011；3.新疆遙感與地理信息系統(tǒng)應(yīng)用重點(diǎn)實(shí)驗(yàn)室，新疆烏魯木齊 830011；4.黃河水利委員會黃河水利科學(xué)研究院，河南鄭州 450003）

在全球變暖背景下，氣候變化加劇，中亞冰川變化表現(xiàn)出明顯區(qū)域差異。青藏高原中部和西北部冰川相對穩(wěn)定，而高原周邊山區(qū)的冰川物質(zhì)虧損嚴(yán)重［1］?？龅貐^(qū)冰川異常，冰川躍動現(xiàn)象延伸到昆侖西部和帕米爾，與中亞其他地區(qū)的負(fù)質(zhì)量平衡形成鮮明對比［2］。目前，喀喇昆侖冰川質(zhì)量增加暫緩，且開始消融［3］。冰川變化形式復(fù)雜多樣，冰川積累或消融必須考慮表面高程變化因素。了解當(dāng)前和未來冰川質(zhì)量的變化，以此預(yù)測干旱區(qū)冰川融水徑流的未來變化情況和自然災(zāi)害的發(fā)生，提前提出防治措施，對于淡水資源安全與經(jīng)濟(jì)可持續(xù)發(fā)展等有重要意義［4］。

目前，測定冰川變化主要有花桿測量法、衛(wèi)星定位監(jiān)測方法、激光測高等［5］。冰川區(qū)坡度大、冰磧覆蓋和冰前湖發(fā)育特征使冰川變化的局部復(fù)雜化，增加了冰川高程變化監(jiān)測的難度［6］。隨著遙感技術(shù)的興起，由于其多時相大范圍覆蓋的特點(diǎn)，被廣泛應(yīng)用于測定冰川高程變化。

本文采用2003—2009 年的ICESat 衛(wèi)星激光測高數(shù)據(jù)，構(gòu)建適合中亞高山區(qū)冰川表面高程數(shù)據(jù)的處理模型，獲取中亞地區(qū)2003年各腳點(diǎn)每年的擬合高程，最后根據(jù)擬合高程結(jié)合降水、溫度以及地形等因素分析研究區(qū)冰川表面高程變化特征及驅(qū)動因素。

1 原理與方法

1.1 點(diǎn)云去噪及其精度優(yōu)化方法

ICESat衛(wèi)星在中低緯度地區(qū)沒有精密航天器指向控制，造成單一重復(fù)軌道在中低緯度地區(qū)并不完全匹配，在軌道上相隔幾百米甚至幾公里［7］。本文通過建立去噪算法和多項(xiàng)式模型等提高ICESat 激光測高數(shù)據(jù)在空間和時間上的一致性。

本文以SRTM 高程數(shù)據(jù)作為基準(zhǔn)面，采用標(biāo)準(zhǔn)差法剔除由于受云層等因素影響的ICESat 高程異常值，從而保證實(shí)驗(yàn)數(shù)據(jù)的可靠性。所有激光腳點(diǎn)與基準(zhǔn)DEM的差值服從正態(tài)分布時，設(shè)激光光束腳點(diǎn)與基準(zhǔn)DEM 的高程差值為D=[d1,d2,d3,···，dN-1,dN]，其平均值為-d，標(biāo)準(zhǔn)差為e，di表示第i個高程腳點(diǎn)，N為激光腳點(diǎn)的總個數(shù)。則激光腳點(diǎn)與基準(zhǔn)DEM的高程差值落在(-3e,+3e)的概率P如下式所示。

如果di服從正態(tài)分布，在3倍標(biāo)準(zhǔn)差的方法中，異常值就可定義為與平均值偏差超過3倍標(biāo)準(zhǔn)差的值；如果di不服從正態(tài)分布，可以采用遠(yuǎn)離平均值的多倍標(biāo)準(zhǔn)差來描述。后續(xù)在2.3.2節(jié)中，具體說明如何根據(jù)正態(tài)分布結(jié)果，采用遠(yuǎn)離1倍標(biāo)準(zhǔn)差，來剔除異常值。

1.2 大范圍冰川高程變化時序重建

以SRTM高程為自變量，采用公式（1）剔除置信度低的高程腳點(diǎn)，并建立區(qū)域高程擬合多項(xiàng)式模型，確定每年激光雷達(dá)衛(wèi)星ICESat 腳點(diǎn)高程與SRTM高程之間的關(guān)系。對每組給定的數(shù)據(jù)擬合點(diǎn)(xi,yi)(其中，i≤m，m為腳點(diǎn)數(shù)量；xi為中腳點(diǎn)i對應(yīng)的重采樣SRTM高程；yi為ICESat激光光束腳點(diǎn)i的高程)，為了獲取xi與yi的n（n≤m）次多項(xiàng)式函數(shù)關(guān)系Pn(xi)。

式中：k為第k（k≤n）多項(xiàng)式；ak為擬合多項(xiàng)式的系數(shù)。當(dāng)擬合值Pn(xi)與yi越接近，表示激光腳點(diǎn)i獲取的高程值越精確。當(dāng)所有激光腳點(diǎn)的擬合高程Pn(xi)與ICESat 腳點(diǎn)高程yi差值的平方和I值最小時，由系數(shù)ak構(gòu)成的多項(xiàng)式模型也就越精確。

當(dāng)擬合函數(shù)為n次多項(xiàng)式時，滿足式（3）的Pn(xi)稱為最小二乘擬合多項(xiàng)式；當(dāng)n=1時稱為線性擬合，I為系數(shù)a0,a1,···，ak(k=0,1,···，n)的多元函數(shù)。上述問題即為求I=I(a0,a1,···，ak)的極值問題。對公式（3）以aj為未知數(shù)求導(dǎo)數(shù)為0 時(j=0,1,···，n)，得多元函數(shù)求極值的必要條件公式如下：

將公式（4）按照關(guān)于aj(j=0,1,···，n)展開，用矩陣表示為式（5）的法方程組。

公式（5）的系數(shù)矩陣是對稱正定矩陣，存在唯一解。當(dāng)n=3時，求解得到的擬合高程Pn(xi)與基準(zhǔn)高程xi相關(guān)系數(shù)R2≥0.95，選n=3 得到區(qū)域的擬合高程，即為激光測高某個時間激光腳點(diǎn)i的擬合高程值與基準(zhǔn)DEM的擬合多項(xiàng)式如下：

式中：xi為基準(zhǔn)高程值；Pn(xi)為待校正高程腳點(diǎn)數(shù)據(jù)；a3，a2，a1，a0分別為擬合多項(xiàng)式的系數(shù)。

2 實(shí)驗(yàn)

2.1 研究區(qū)概況

中亞山區(qū)位于亞洲中部海拔在3500 m以上，冰川面積約1.3×105km2。為研究中亞高山冰川區(qū)內(nèi)部的冰川變化差異，本文根據(jù)地貌、地形等生態(tài)地理區(qū)域系統(tǒng)將研究區(qū)劃分為3個區(qū)域（圖1）。I區(qū)主要包括西藏及青海南部地區(qū)，該區(qū)東南部冰川大部分屬于海洋性冰川，內(nèi)部冰川大多屬于大陸性冰川；II區(qū)主要包括新疆地區(qū)及青海北部和甘肅北部的部分區(qū)域，該區(qū)冰川以大陸性冰川為主，但天山、昆侖山東段以及喀喇昆侖山北坡等部分地區(qū)冰川屬于亞大陸性冰川；III區(qū)在中國境外，與中國相鄰，地勢東高西低，大部分為溫帶大陸性氣候。

圖1 中亞冰川分區(qū)概況示意圖Fig.1 Overview of glacier zoning in Central Asia

2.2 數(shù)據(jù)與預(yù)處理

本文使用數(shù)據(jù)主要包含以下3部分：（1）ICESat數(shù)據(jù)來源于冰、云和陸地測高衛(wèi)星，衛(wèi)星于2003年1月13 日成功發(fā)射，是首顆激光測高衛(wèi)星，搭載地學(xué)激光測高系統(tǒng)GLAS［8］；ICESat數(shù)據(jù)產(chǎn)品在美國國家冰雪數(shù)據(jù)中心(http://www.nsidc.org)免費(fèi)獲取。根據(jù)研究目的，本文采用2003—2009 年的ICESat GLAH14數(shù)據(jù)，為便于時序重建，避免降雪和融化等季節(jié)間冰川高程變化的影響，選擇每年的同一月份數(shù)據(jù)［9］；根據(jù)統(tǒng)計(jì)的數(shù)據(jù)量，選擇數(shù)據(jù)量最多且最完整的每年3 月的數(shù)據(jù)進(jìn)行處理分析。（2）本文選擇能夠免費(fèi)獲取的90 m 分辨率包含研究區(qū)范圍的SRTM3數(shù)據(jù)作為參考DEM。（3）冰川編目數(shù)據(jù)采用2017年7月發(fā)布的RGI 6.0(Randolph Glacier Inventory 6.0)數(shù)據(jù)庫，RGI 6.0 以2010 年左右的夏季無云Landsat TM/ETM+為主要影像來源，用于獲取研究區(qū)冰川邊界。

2.3 數(shù)據(jù)處理

首先從ICESat 數(shù)據(jù)中提取腳點(diǎn)的坐標(biāo)、高程，進(jìn)行坐標(biāo)轉(zhuǎn)換和地理配準(zhǔn)；再剔除受地形、云層等因素影響產(chǎn)生的異常值，本文根據(jù)正態(tài)分布檢驗(yàn)結(jié)果選定異常值并將其剔除。利用去除異常值后的ICESat 坐標(biāo)數(shù)據(jù)，在SRTM DEM 數(shù)據(jù)中重采樣，獲取ICESat 腳點(diǎn)坐標(biāo)對應(yīng)的SRTM DEM 高程，建立每年的多項(xiàng)式擬合函數(shù)，以解決ICESat 數(shù)據(jù)重復(fù)軌道不完全匹配的問題，并計(jì)算相關(guān)系數(shù)R2，作為評估擬合效果參數(shù)。以2003 年的SRTM DEM 高程數(shù)據(jù)作為基準(zhǔn)，獲取2003年每個腳點(diǎn)數(shù)據(jù)每年的擬合高程，實(shí)現(xiàn)時序重建。

2.3.1 數(shù)據(jù)預(yù)處理 由于本文始終以SRTM DEM 為基準(zhǔn)進(jìn)行時序重建，因此需將Topex/Poseidon 高程轉(zhuǎn)換為WGS84高程，如下式。

式中：H84為轉(zhuǎn)換后點(diǎn)云WGS84高程；hglas為ICESat點(diǎn)云相對于Topex/Poseidon 參考橢球的高程（m）；N是大地水準(zhǔn)面和參考橢球面的距離（m）；數(shù)值0.7的單位為m，為2 個參考橢球間的高程差異［10］。本文在30 km×30 km 尺度范圍內(nèi)選取最高點(diǎn)作為SRTM DEM 山頂點(diǎn)，對比參考橢球?yàn)閃GS84 高分辨率四維地球影像篩選均勻分布可靠的山頂點(diǎn)，記錄其作為控制點(diǎn)的150 個山頂點(diǎn)的經(jīng)緯度坐標(biāo)，對SRTM DEM 進(jìn)行配準(zhǔn)，提高不同坡度、坡向的ICESat 腳點(diǎn)精度［11］。

2.3.2 剔除異常值 根據(jù)獲取的冰川邊界篩選出研究區(qū)內(nèi)ICESat 數(shù)據(jù)高程點(diǎn)。以Ⅰ區(qū)的2003 年數(shù)據(jù)為例，檢驗(yàn)ICESat 腳點(diǎn)高程是否符合正態(tài)分布。如表1所示，顯著性為0.00＜0.05，所以該區(qū)域高程值不服從正態(tài)分布，因此，選定遠(yuǎn)離平均值1倍標(biāo)準(zhǔn)差的值為異常值并將其剔除，確保數(shù)據(jù)的準(zhǔn)確性。

表1 I 區(qū)的2003年ICESat腳點(diǎn)高程正態(tài)分布檢驗(yàn)Tab.1 Normal distribution test of ICESat foot elevation in area I in 2003

2.3.3 多項(xiàng)式擬合函數(shù) 以SRTM DEM高程為自變量，與ICESat 高程進(jìn)行多項(xiàng)式擬合，建立多項(xiàng)式擬合函數(shù)如式（6），并計(jì)算相關(guān)系數(shù)R2。如圖2 所示，通過對比Ⅰ區(qū)內(nèi)2003 年點(diǎn)云數(shù)據(jù)去除異常值前后擬合效果，發(fā)現(xiàn)去除異常值的R2更接近1，擬合函數(shù)的擬合效果更好，更符合冰川表面高程的實(shí)際情況。

圖2 Ⅰ區(qū)2003年數(shù)據(jù)去除異常值前后對比Fig.2 Comparison of data in zone I before and after removing outliers in 2003

3 結(jié)果與分析

3.1 多項(xiàng)式擬合模型評估

由于衛(wèi)星數(shù)據(jù)的處理過程帶有保密性，難以給予物理模型對數(shù)據(jù)源進(jìn)行校正。本文通過統(tǒng)計(jì)學(xué)的方法，采用確定系數(shù)R2評估多項(xiàng)式擬合效果，表示模型對現(xiàn)實(shí)數(shù)據(jù)擬合的程度，R2的計(jì)算公式如下。

式中：SSR 為預(yù)測數(shù)據(jù)與原始數(shù)據(jù)均值之差的平方和；SST 即原始數(shù)據(jù)和均值之差的平方和；SSE 為擬合數(shù)據(jù)和原始數(shù)據(jù)對應(yīng)點(diǎn)的誤差的平方和；R2的取值范圍為［0,1］。當(dāng)擬合數(shù)據(jù)和原始數(shù)據(jù)對應(yīng)點(diǎn)的誤差越小時，SSE越小，擬合數(shù)據(jù)越可靠，此時R2的取值越接近1。因此，當(dāng)R2的取值越接近1，擬合效果越好。

由表2 所示，整體冰川區(qū)、I 區(qū)、II 區(qū)、III 區(qū)的R2都大于0.98，說明整體擬合效果很好。冰川區(qū)由物質(zhì)平衡線分為消融區(qū)和積累區(qū)，在物質(zhì)平衡線以下為消融區(qū)，物質(zhì)平衡線以上為積累區(qū)，冰川表面高程變化可總體歸為這2 種模式，因此利用三次多項(xiàng)式更能符合2 種模式的冰川表面高程變化現(xiàn)象，且符合冰川表面高程自身變化規(guī)律。

表2 研究區(qū)2003—2009年3月多項(xiàng)式擬合模型的R2Tab.2 R2of the polynomial fitting model from 2003 to March 2009 in the study area

3.2 冰川表面高程變化分析

在數(shù)據(jù)處理的過程中由于系統(tǒng)誤差會產(chǎn)生部分異常值，本研究采用式（9），根據(jù)正態(tài)性檢驗(yàn)結(jié)果，將數(shù)據(jù)中與變化平均值的偏差超過2 倍標(biāo)準(zhǔn)差的高程變化平均值作為異常值進(jìn)行剔除。各區(qū)不同時間段內(nèi)冰川表面高程平均變化量見表3。

表3 各區(qū)不同時間段內(nèi)冰川表面高程平均變化量Tab.3 Average variation of glacier surface elevation in different periods of time in each region

式中：E為誤差=高程變化值-高程變化均值。

由圖3a可知，整體區(qū)域即整個中亞地區(qū)所有的冰川區(qū)域，總的冰川表面高程下降9.59±1.89 m，在2004年高程下降了5.89±3.76 m，次年回升到與2003年相當(dāng)?shù)母叨?，之后逐年緩慢下降；圖3b中，Ⅰ區(qū)冰川表面高程在2003—2004 年間下降5.78±2.13 m，2005 年略有回升，之后高程逐年緩慢下降，變化斜率為-0.581 m·a-1；圖3c 中，II 區(qū)其冰川表面平均高程在2003—2009 年間下降了7.87±5.03 m，2005 年冰川表面高程出現(xiàn)大幅上升，甚至超過2003年的高程值，但整體仍處于下降趨勢，變化斜率為-1.226 m·a-1；圖3d所示，與其他分區(qū)相比，Ⅲ區(qū)每年的冰川表面平均高程變化波動較大。

圖3 各個研究區(qū)冰川表面高程變化時序重建Fig.3 Time series reconstruction of glacier surface elevation change in each study area

如圖4所示，整個研究區(qū)內(nèi)天山、帕米爾高原和青藏高原東南部地區(qū)的冰川高程下降劇烈，青藏高原南部地區(qū)冰川高程呈輕微下降趨勢。中部的昆侖山西段地區(qū)部分冰川高程呈現(xiàn)升高狀態(tài)，西段以及昆侖山東部的祁連山地區(qū)都呈下降趨勢。北部的天山地區(qū)冰川高程變化呈現(xiàn)出明顯的區(qū)域差異，東天山、北天山和西天山南部地區(qū)冰川高程下降明顯，西天山北部冰川高程下降速率相對較慢。整個研究區(qū)內(nèi)除昆侖山部分地區(qū)冰川高程增厚外，其他地區(qū)整體呈下降趨勢，只是下降速率存在地區(qū)差異。與2019 年Treichler 等［2］的研究相比，各個區(qū)域的冰川表面高程變化趨勢基本符合。但由于本文只采用了每年3 月的數(shù)據(jù)，數(shù)據(jù)量小，分布稀疏，部分區(qū)域缺少數(shù)據(jù)支持且表面高程變化特征表現(xiàn)性差，例如在Treichler等［2］的研究中，中天山東南部地區(qū)、昆侖山大部分地區(qū)以及青藏高原中部的部分地區(qū)冰川呈輕微增厚狀態(tài)，但本文中在該部分地區(qū)數(shù)據(jù)稀少，無法確認(rèn)冰川變化情況。且本文中顏色柵格表示小區(qū)域范圍內(nèi)激光點(diǎn)的變化率，因此，與Treichler 等［2］研究中大分區(qū)冰川表面高程變化率存在部分差異。

圖4 各區(qū)冰川表面變化速率分布Fig.4 Distribution map of glacier surface change rate in each region

根據(jù)2019年Treichler等［2］對MERRA-2和ERAInterim氣候再分析數(shù)據(jù)的處理，顯示在2003—2009年間，氣溫呈上升趨勢，除昆侖山以及周圍少數(shù)地區(qū)外，亞洲冰川表面高程整體呈下降趨勢。青藏高原夏季降水量，整體呈現(xiàn)自東南向西北遞減的分布規(guī)律［12］，東北部干旱區(qū)雨季降水量也呈增加趨勢［13］。圖3 中，在2004—2005 年間各區(qū)的冰川表面高程變化劇烈。據(jù)相關(guān)研究［14］顯示，2003 年12 月至2004年2月，研究區(qū)冬、春季大部分地區(qū)氣溫較常年同期明顯偏高，降水偏少，造成冰川消融。2004 年新疆南部等地夏季高溫日數(shù)普遍達(dá)到10～20 d，局部地區(qū)達(dá)到40 d［15］。2005 年西部地區(qū)伴有降水降雪的增加，造成各區(qū)的冰川表面高程在2005年都呈上升趨勢［16］。其中，新疆東南部和西部以及青藏高原北部等地較常年同期偏高3倍［17］。

地形上，冰川在低海拔坡度小于30°時，冰川厚度減薄隨著坡度的減小而加劇，冰川表面高程降低較快；中高海拔地區(qū)，溫度逐漸降低，受降水的影響較大，易導(dǎo)致冰川表面高程增高；而在高海拔地區(qū)平均坡度＞35°，易發(fā)生冰崩、雪崩等現(xiàn)象，使得冰川物質(zhì)減少［18］。結(jié)合地形因素，降水隨海拔增高而增多，但有一個最大降水量高度，超過此高度，山地降水不再隨高度遞增，隨氣候干濕而異［19］。I區(qū)中，青藏高原平均海拔4000 m以上，海拔高度對氣溫的影響已超過緯度位置作用。因此，青藏高原除南部和東南部的冰川退縮較嚴(yán)重外，大部分區(qū)域冰川高程變化相對緩慢［20］。II 區(qū)南部有昆侖山和阿爾金山，這2個地區(qū)海拔高，干旱少雨，冰川高程變化相對穩(wěn)定，但天山西北部受降水量影響較大的外山脈冰川融化速度較快。祁連山地表平均溫度在海拔高的地區(qū)變暖趨勢明顯，因此，祁連山地區(qū)冰川表面高程緩慢下降［21］。Ⅲ區(qū)中西帕米爾高原山脈交錯復(fù)雜，降水豐富，該地區(qū)冰川高程變化強(qiáng)烈。東帕米爾高原東南邊緣高山阻隔了來自印度洋、太平洋的暖濕氣流，氣候干燥，冰川高程緩慢下降［22］。就下降趨勢而言，研究區(qū)內(nèi)境外區(qū)域冰川表面高程下降趨勢最快，為-1.663 m·a-1，新疆地區(qū)次之，青藏高原地區(qū)最為緩慢。

4 結(jié)論

本文針對ICESat 衛(wèi)星測高數(shù)據(jù)，考慮到重復(fù)軌道在中低緯度地區(qū)空間匹配性弱的問題，選取SRTM DEM 為高程基準(zhǔn)面，建立冰川區(qū)點(diǎn)云去噪及其精度優(yōu)化算法和多尺度冰川區(qū)表面高程時間序列分析模型，為點(diǎn)云類遙感監(jiān)測冰川表面高程變化提供參考。主要結(jié)論如下:

（1）采用正態(tài)分布顯著性檢驗(yàn)，進(jìn)行粗差剔除，在不服從正態(tài)分布情況下，選定遠(yuǎn)離線性回歸1 倍標(biāo)準(zhǔn)差的高程值為異常值，能夠有效剔除冰川區(qū)ICESat腳點(diǎn)高程中的異常值。

（2）利用三次多項(xiàng)式模型，對中亞山區(qū)整體和分區(qū)域的不同尺度冰川表面高程進(jìn)行擬合，模型R2都在0.98 以上，表明ICESat 數(shù)據(jù)和SRTM 高程數(shù)據(jù)的三次多項(xiàng)式關(guān)系在該區(qū)具有較強(qiáng)的普適性。

（3）亞洲高山區(qū)的冰川表面高程2003—2009年整體下降幅度在10 m以內(nèi)，表現(xiàn)出明顯的區(qū)域差異，對于2004 年前后表現(xiàn)出強(qiáng)烈變化，由于環(huán)流和季風(fēng)等因素的影響，造成的溫度以及降水異常，從而導(dǎo)致冰川異常波動。

本研究方法對冰川區(qū)點(diǎn)云類高程腳點(diǎn)監(jiān)測具有一定的通用性，但會使點(diǎn)云數(shù)據(jù)更稀疏，另對基準(zhǔn)DEM的依賴度較高。后續(xù)，研究將進(jìn)一步加入其他冰川區(qū)點(diǎn)云數(shù)據(jù)，研究時間序列建立所適用的空間尺度，建立更靈活多變的區(qū)域冰川表面高程時序監(jiān)測模型，對影響冰川表面高程變化的地形、環(huán)境等因子進(jìn)行更有效全面的分析。