2000–2020 年喀喇昆侖山洪扎河流域冰川高程變化數據集

吳坤鵬，劉時銀，朱鈺

1.中國科學院空天信息創新研究院可持續發展大數據國際研究中心，北京 100094

2.云南大學國際河流與生態安全研究院，昆明 650091

3.中國科學院西北生態環境資源研究院國家冰川凍土沙漠科學數據中心，蘭州 730000

引 言

作為冰凍圈系統的重要組成部分，冰川與區域氣候、自然環境演化、海平面變化、地球表面過程等有著密切聯系[1-6]。20 世紀90 年代以來，亞洲高山區冰川不斷退縮，冰川持續消融[7-8]，但空間差異顯著，主要表現為喜馬拉雅地區東部和念青唐古拉山脈的冰川強烈消融，喀喇昆侖山–帕米爾地區冰川呈穩定平衡狀態，甚至物質有所增加，稱之為“喀喇昆侖異?！盵9-10]?？錾奖ㄗ兓P系到下游地區的工農業、水力發電等的正常運作，其中西喀喇昆侖山是中國政府推動的“一帶一路”陸上絲綢之路關鍵區之一——中巴經濟走廊的必經之地，冰川表面高程變化關系到該地區水資源利用、基礎設施建設、災害防治預警，因此對該地區的冰川變化研究具有重要的社會意義。

當前對冰川變化研究運用最為廣泛的方法有傳統冰川學法、遙感監測法、模型模擬法等。傳統冰川學法通過采集足夠數量的花桿和雪坑數據估算冰川物質平衡，是監測單條冰川運用最為廣泛的方法；但冰川多發育在邊遠山區，自然環境惡劣，能夠開展野外實地監測的冰川有限[11]。模型模擬法可以模擬年際時間尺度的冰川物質平衡，并預測未來氣候變化情景下冰川變化趨勢[12-13]，但該方法需要詳細的冰川和氣象觀測資料驅動。遙感監測法通過比較不同時期遙感數據獲取的冰川規模、表面高程、運動速度等，分析冰川對氣候變化的響應[14-16]。相較于傳統冰川學法和模型模擬法，遙感能夠開展長時間序列、大空間尺度的冰川變化監測。隨著遙感技術的發展和海量數據的積累，遙感監測法已成為研究冰川變化的主要方法。

當前對喀喇昆侖山地區冰川變化研究大多數依賴國外遙感衛星數據，如KH9、ASTER、SPOT等，且與冰川高程變化相關的開源數據集主要為基于ASTER 提取的冰川高程變化[8-10]，包括2000–2016 年高亞洲冰川高程變化數據集[10]（空間分辨率30 m）和2000–2020 年全球冰川高程變化數據集[8]（空間分辨率100 m）。為凸顯我國遙感數據在科學研究中的實際價值，本文選用2020 年1–2 月獲取的資源3 號衛星立體像對數據，基于大地測量法[17]，通過DEM 誤差校正、冰雪穿透估算等[18]，獲取了洪扎河流域冰川表面高程變化，并以穩定地形區殘余高程差為統計特征，判別和驗證數據集產品的可靠性，為“喀喇昆侖異?！毖芯考霸摰貐^水資源利用、冰川災害防治等提供更精準的數據支撐。

1 數據采集和處理方法

1.1 數據源

2000–2020 年喀喇昆侖山洪扎河流域冰川高程變化數據集主要基于三種數據：2000 年2 月的數字高程模型（SRTM DEM）、2020 年1–2 月獲取的資源三號衛星立體像對數據和冰川編目數據（表1、圖1）。

表1 資源三號立體像對參數信息Table 1 Specifications of ZY3 stereo images

SRTM DEM 為美國地質調查局提供的、未填補“空洞”的 SRTM1 C 波段數據產品（https://earthexplorer.usgs.gov/），獲取時間為2000年2月，空間分辨率30 m，坐標系WGS84/EGM96。在90%的置信區間內，SRTM 垂直精度優于16 m，其高程精度受地形的影響較大，在平坦區域的高程精度可達10 m，且隨地表坡度變化[19]。SRTM 有兩套合成孔徑雷達系統，分別采用C 波段（波長5.6 cm，頻率2.2 Hz）和X 波段（波長3.1 cm，頻率8.8 Hz），為有效評估C 波段對積雪的穿透，本文還選取了覆蓋喀喇昆侖山地區的SRTM X 波段數據，通過C 波段和X 波段差值比較，估算C 波段對積雪的穿透深度[20]。

資源三號是我國首顆民用三線陣立體測繪衛星，于2012 年1 月9 日成功發射。資源三號衛星搭載1 臺多光譜和3 臺全色線陣光學相機。全色相機包括前視、正視和后視組成三線陣立體相機，其中前視和后視地面分辨率為3.5 m[21-22]。

冰川編目數據由國際冰雪數據中心（https://nsidc.org/data/glims）的RGI 6.0（Randolph Glacier Inventory 6.0）[23]，主要用于冰川區與非冰川區的劃分。由于喀喇昆侖山洪扎河流域存在大量躍動前進冰川，為了確保冰川表面高程變化的完整性，本文還采用了RGI 6.0[23]、中國第二次冰川編目[24]及GAMDAM（Glacier Area Mapping for Discharge from the Asian Mountains）[25]冰川編目的集合，以保證在研究時段冰川邊界的最大范圍。

1.2 數據處理方法和流程

冰川表面高程變化獲取主要由2 部分組成：基于立體像對提取2020 年數字高程模型，基于兩期數字高程模型獲取2000–2020 年的冰川表面高程變化。本文采用ENVI、ArcGIS 和Python 軟件平臺進行數據的處理和分析，其中ENVI 用于DEM 提取，ArcGIS 用于DEM 編輯與鑲嵌，Python 用于DEM 配準，主要處理流程如圖2。

圖2 ZY3 DEM 的提取及冰川表面高程變化生成流程示意圖Figure 2 The flow chart of ZY3 DEM extraction and glacier height changes

基于資源三號立體像對提取數字高程模型總體分為6 步，分別是輸入立體像對、定義地面控制點、定義連接點、設定DEM 提取參數、輸出并編輯DEM、鑲嵌DEM[26]。（1）輸入立體像對：左影像加載正視影像，右影像加載后視/前視影像，加載影像后，ENVI 軟件能夠自動識別RPC 文件；（2）定義地面控制點：ENVI 軟件提供不定義、交互式定義和讀取控制點文件3 種定義地面控制點的方式，其中不定義控制點提取的是相對高程，水平面為WGS84 基準面。由于獲取冰川表面高程變化需要進行兩期DEM 配準，因此提取DEM 可采用不定義控制點方式；（3）定義連接點：ENVI 軟件提供自動尋找、交互式定義和讀取連接點文件3 種定義連接點的方式，本文先選擇自動尋找連接點，然后進行交互式編輯，確保連接點誤差達到最??；（4）設定DEM 提取參數：設定投影參數為WGS84 UTM Zone 43N，輸出像元大小為10 m，背景值為0；（5）輸出并編輯DEM：生成DEM 并檢查由云量等造成的DEM 錯誤，通過手動編輯修改高程數據；（6）鑲嵌DEM：利用ArcGIS 軟件將提取的DEM 進行鑲嵌，得到研究區相對DEM（以下稱ZY3 DEM）。

兩期DEM 數據（SRTM DEM 和ZY3 DEM）之間的配準采用大地測量法。同一區域不同時間的兩幅DEM 數據，在沒有配準的情況下直接高程相減得到高程差值圖，而高程差值與坡度坡向之間存在如下余弦函數關系：

式中：dh為不同DEM 數據之間各個像元的高程差，a是平移矢量的大小，b是平移矢量的方向， 和φ分別為像元對應的坡度和坡向， ?????為不同DEM 數據之間的整體高程差異（垂直矢量的大?。?。從兩幅DEM 對應的像元屬性中可以獲得dh、坡度 和坡向φ。不同DEM 數據之間的高程差存在異常值，影響余弦曲線的擬合。本文對高程差樣本進行統計，認為低于5%或高于95%分位數的高程差是異常值，進行剔除。坡度的大小對DEM 數據匹配的精度有很大影響，因此要剔除地形平坦的區域和坡度太大的區域。本文剔除了地面坡度小于5°和大于80°的區域。然后通過最小二乘法擬合余弦曲線可以得到系數a、b、c，不同DEM 數據在X、Y 和Z 方向的平移量可表示為：

式中：α?為DEM 數據的平均坡度。求出shiftx，shifty和shiftz，并按照在X、Y 和Z 方向的平移量進行平移，但一次平移后的DEM 并不一定能夠完全消除空間匹配錯位造成的誤差，需要將平移后的DEM 再次迭代平移，重復上述步驟計算dh標準差和平移矢量。當dh的標準差減小幅度小于2%或平移距離小于0.5 m，即達到了兩幅DEM 水平方向的配準要求[17]。最后利用非冰川區高程差與最大曲率的關系，校正由不同分辨率造成的冰川區的高程差異殘差。最終得到的高程差值即為研究區高程變化信息。

2 數據樣本描述

本數據集的數據存儲格式為 GeoTIFF 格式，數據類型為 32 位浮點型，命名為DH_Karakoram_Hunza_20002020.tif，其中DH 表示為Difference Height，Karakoram_Hunza 表示為喀喇昆侖山洪扎河流域，20002020 表示為2000 年至2020 年時間段，數據單位是米，地理坐標系WGS84 UTM Zone 43N，空間分辨率30 m。樣本展示如圖3。

圖3 喀喇昆侖山洪扎河流域冰川表面高程變化Figure 3 The glacier height changes in Hunza Basin

3 數據質量控制和評估

由于SRTM 的C 波段能夠穿透積雪，在測定積雪覆蓋的冰川表面地形時得不到真實的地面信息。穿透深度隨著地面狀況的變化而不同，以新雪和粒雪的穿透深度最大，表磧物覆蓋的冰面穿透深度最小[27]，穿透深度能介于0–10 m 之間[20,28-29]。相比SRTM 的C 波段，X 波段對積雪的穿透效應要小[20]。利用配準和校正后的SRTM X 波段DEM 和SRTM C 波段DEM 相比較，得到SRTM C波段在喀喇昆侖山地區平均穿透1.6±0.05 m。SRTM DEM 和ZY3 DEM 在高海拔地區均存在一定“空洞”，但通過大地測量法獲取的冰川表面高程變化數據，在不同海拔均有變化信息，數據空洞占比較小，對冰川表面高程變化統計分析不會產生顯著影響。

基于SRTM DEM 和ZY3 DEM 獲取地面高程變化后，可利用穩定地形區高程差的均方根誤差（Root Mean Square Error,RMSE）或標準差（Standard Deviation,STDV）作為高程差不確定性的初步估計（圖4），但是會高估高程變化的不確定性，因為誤差在一個較大區域內平均后會減小[30-31]。根據空間鄰近相似性原理，相鄰的像元間存在自相關，因此需要除去自相關造成的影響。為此，引入標準平均誤差（Standard Error of the mean,SE）對高程差不確定性進行評估。

圖4 無冰區高程差及其直方圖Figure 4 Elevation changes and the histogram in off-glacier areas

式中：n為去自相關后的像元個數。除去自相關造成的影響，需要根據DEM 的分辨率計算空間自相關的長度?？臻g自相關距離可以通過自定義的方式確定，也可以通過莫蘭指數（Moran’s I）計算。研究表明，30 m 分辨率選取600 m 的自相關長度，10–20 m 分辨率選取400 m 的自相關長度較為合適[30]。因此，本文采用自定義方式，自相關長度選擇400 m。最終SRTM DEM 和ZY3 DEM 數據高程差的精度（σ）用無冰區高程差的均值MED和標準平均誤差SE計算：

結果顯示，誤差校正后，SRTM DEM 和ZY3 DEM 數據高程差的誤差減小到±1.1 m。

4 數據價值

相較于開源的冰川表面高程變化數據集，本數據集比2000–2016 年高亞洲冰川高程變化數據集[10]有更長的時間跨度、比2000–2020 年全球冰川高程變化數據集[8]具有更高的空間分辨率。更長的時間跨度與更高空間分辨率，能夠為喀喇昆侖山地區冰川變化研究提供更充實的基礎數據。

“喀喇昆侖異?！笔钱斍把芯繜狳c，主要表現為輕微的物質正平衡及大量躍動冰川的分布[15]。開源躍動冰川分布數據集多是根據冰川末端進退及冰川運動速度變化進行制備[32-33]，對冰川躍動過程中冰量的輸送研究較為欠缺，難以充分開展冰川躍動過程及其機理研究。本文基于SRTM DEM 和ZY3 DEM 獲取了喀喇昆侖山洪扎河流域21 世紀初冰川表面高程變化，為更完整的躍動冰川識別及典型冰川的躍動控制機理研究奠定了數據基礎。與此同時，洪扎河流域是中巴經濟走廊的必經之地，全面了解冰川表面高程變化、估算冰川物質平衡，為該地區水資源利用及冰川災害防治提供必要的基礎信息。

5 數據使用方法和建議

喀喇昆侖山洪扎河流域2000–2020 年冰川表面高程變化數據集通過國家冰川凍土沙漠數據中心NCDC（National Cryosphere Desert Data Center）網站提供全部數據下載服務。本數據集數據存儲格式為32 位浮點型GeoTIFF 柵格數據，地理坐標系WGS84 UTM Zone 43N?？梢栽贏rcGIS、SuperMap等GIS 平臺軟件直接調用，也可以使用MATLAB 等數據分析軟件編碼進行批量處理和分析。研究人員可利用該數據與其他時期冰川物質平衡相比較，或選擇該數據集作為冰川模型的輸入參數/驗證參數，開展冰川對氣候變化響應研究工作。

致 謝

感謝美國地質調查局（USGS）提供SRTM1 DEM 數據，感謝國際冰雪數據中心及國家冰川凍土沙漠科學數據中心提供冰川編目數據集。

數據作者分工職責

吳坤鵬（1990—），男，安徽省安慶市人，博士，副研究員，研究方向為冰川遙感。主要承擔工作：數據篩選，數據處理，論文撰寫。

劉時銀（1963—），男，河南省信陽市人，博士，研究員，研究方向為冰川水文。主要承擔工作：研究思路設計。

朱鈺（1992—），男，甘肅省平涼市人，博士，助理研究員，研究方向為冰川遙感。主要承擔工作：數據處理。