吳坤鵬，劉時銀，朱鈺

1.中國科學(xué)院空天信息創(chuàng)新研究院可持續(xù)發(fā)展大數(shù)據(jù)國際研究中心，北京 100094

2.云南大學(xué)國際河流與生態(tài)安全研究院，昆明 650091

3.中國科學(xué)院西北生態(tài)環(huán)境資源研究院國家冰川凍土沙漠科學(xué)數(shù)據(jù)中心，蘭州 730000

引 言

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

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

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

1 數(shù)據(jù)采集和處理方法

1.1 數(shù)據(jù)源

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

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

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

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

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

1.2 數(shù)據(jù)處理方法和流程

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

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

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

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

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

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

2 數(shù)據(jù)樣本描述

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

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

3 數(shù)據(jù)質(zhì)量控制和評估

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

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

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

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

結(jié)果顯示，誤差校正后，SRTM DEM 和ZY3 DEM 數(shù)據(jù)高程差的誤差減小到±1.1 m。

4 數(shù)據(jù)價值

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

“喀喇昆侖異常”是當(dāng)前研究熱點，主要表現(xiàn)為輕微的物質(zhì)正平衡及大量躍動冰川的分布[15]。開源躍動冰川分布數(shù)據(jù)集多是根據(jù)冰川末端進退及冰川運動速度變化進行制備[32-33]，對冰川躍動過程中冰量的輸送研究較為欠缺，難以充分開展冰川躍動過程及其機理研究。本文基于SRTM DEM 和ZY3 DEM 獲取了喀喇昆侖山洪扎河流域21 世紀(jì)初冰川表面高程變化，為更完整的躍動冰川識別及典型冰川的躍動控制機理研究奠定了數(shù)據(jù)基礎(chǔ)。與此同時，洪扎河流域是中巴經(jīng)濟走廊的必經(jīng)之地，全面了解冰川表面高程變化、估算冰川物質(zhì)平衡，為該地區(qū)水資源利用及冰川災(zāi)害防治提供必要的基礎(chǔ)信息。

5 數(shù)據(jù)使用方法和建議

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

致 謝

感謝美國地質(zhì)調(diào)查局（USGS）提供SRTM1 DEM 數(shù)據(jù)，感謝國際冰雪數(shù)據(jù)中心及國家冰川凍土沙漠科學(xué)數(shù)據(jù)中心提供冰川編目數(shù)據(jù)集。

數(shù)據(jù)作者分工職責(zé)

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

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

朱鈺（1992—），男，甘肅省平?jīng)鍪腥耍┦?，助理研究員，研究方向為冰川遙感。主要承擔(dān)工作：數(shù)據(jù)處理。