2015–2017年中巴公路沿線山地冰川表面運動數據集

2015–2017年中巴公路沿線山地冰川表面運動數據集
作者：王磊　蔣宗立　劉時銀　張震
2018年11月9日

本作品收錄於《中國科學數據》

王磊, 蔣宗立, 劉時銀, 張震. 2015–2017年中巴公路沿線山地冰川表面運動數據集[J/OL]. 中國科學數據, 2018, 3(4). (2018-10-30). DOI: 10.11922/csdata.2018.0006.zh.

摘要&關鍵詞

摘要：中巴公路是一條連接中國西部與巴基斯坦及南亞、西亞的重要國際公路，對於沿線國家及高亞洲地區經濟文化有着非常重要的意義。而公路沿線的帕米爾高原和喀喇崑崙山脈是全球山嶽冰川分布最為密集的區域之一，冰川的運動時常對公路的安全運營造成威脅。本文獲取了2014年10月至2018年3月間覆蓋該區域（35°50′53″N–39°00′42″N，74°17′58″E–75°53′01″E）的Sentinel1A的SAR數據，使用特徵匹配方法計算得到了中巴公路沿線冰川這3年間的表面運動場。得益於Sentinel1A數據的IW模式，計算結果的空間分辨率約為100 m×100 m，靜止區域殘餘運動的統計分析表明總體精度約0.04 m/d。分析表明，中巴公路沿線北部和中部地區冰川表面流速相對較為穩定，而南部地區的大規模冰川流速較快，部分冰川積累區表面流速高達4 m/d以上。本數據集可作為高亞洲區域冰川物質平衡分析的基礎數據，支撐高亞洲區域水資源狀況評估以及氣候變化方面的研究。另外，山地冰川表面運動場的高精度監測結果對中巴公路及沿線居民區進行冰川災害的預警也有重要意義。

關鍵詞：中巴公路；特徵匹配方法；冰川表面運動；Sentinel1A

Abstract & Keywords

Abstract: None As an international highway, the Karakoram Highway (KKH) is of great importance to both China and Pakistan. Crossing one of the most glacier-concentrated regions – the Karakoram Mountains, KKH constantly faces threats from glacier motions along the way. This study obtained the displacement field of glacier surface along the KKH during October 2014 and March 2018 by feature-tracking method, using Sentinel1A SAR data of the research area of the same period. Results indicated that glacier surface velocities in the northern and central parts of the region are relatively stable, while parts of the large-scale glaciers in the southern part had higher velocities up to 4 m/d, particularly along the section of the highway near Batura Sar. This dataset can be used as basic data in research of glacial mass balance, as well as of hydrological resource dynamics and climate changes. In addition, the high-precision monitoring results of the displacement field are of great significance for glacier disaster warnings to KKH and the residents around.

Keywords: KKH; feature-tracking; glacier surface motion; Sentinel1A

數據庫（集）基本信息簡介

數據集名稱	2015–2017年中巴公路沿線山地冰川表面運動數據集
數據作者	王磊，蔣宗立，劉時銀，張震
數據通信作者	蔣宗立（jiangzongli@hnust.edu.cn）
數據時間範圍	2014年10月31日至2018年3月8日（時間間隔12/24天）
地理區域	35°50′53″N–39°00′42″N，74°17′58″E–75°53′01″E
空間分辨率	100 m×100 m
數據量	9.51 GB
數據格式	*.tif ( GeoTIFF, 32bit float )
數據服務系統網址	http://www.sciencedb.cn/dataSet/handle/628
基金項目	國家自然科學基金項目（41471067，41474014，41701087，41761144075)；科技部專項（2013FY111400）。
數據集組成	數據集由兩部分組成，一是冰川表面運動方向柵格數據集；二是冰川表面運動速度柵格數據集。文件夾名稱分別為Direction和Velocity。

Dataset Profile

Title	Dataset of glacier surface motion along KKH during 2015 – 2017
Data authors	Wang Lei, Jiang Zongli, Liu Shiyin, Zhang Zhen
Data corresponding author	Jiang Zongli (jiangzongli@hnust.edu.cn)
Time range	2014-10-31 to 2018-03-08（Interval 12/24day）
Geographical scope	35°50′53″N–39°00′42″N，74°17′58″E–75°53′01″E
Spatial resolution	100 m×100 m
Data volume	9.51 GB
Data format	*.tif (GeoTIFF, 32 bit float)
Data service system	http://www.sciencedb.cn/dataSet/handle/628
Sources of funding	National Natural Science Foundation of China (Nos. 41471067, 41474014, 41701087, 41761144075); Special project of the Ministry of Science and Technology (2013FY111400).
Dataset composition	The dataset consists of two data files: 「Direction」 stores raster data of glacier surface motion directions, and 「Velocity」 stores data of glacier surface motion speeds.

引 言

中巴公路（喀喇崑崙國際公路，Karakoram Highway，簡稱KKH）作為一條國際公路，是中巴友誼的重要見證，也是中國西部與中亞、西亞、南亞地區的交通要道，對中國和巴基斯坦以及周邊國家的經濟文化交流有着重要的意義。公路穿越的喀喇崑崙山系是全球低緯度冰川分布數量最多、規模最大的區域之一^[1]。公路經過的低海拔地區為典型的大陸性氣候，冬夏溫差大，降水稀少，而公路沿線的冰川所處的高海拔地區則為典型的高山氣候。受到減弱的印度洋季風和加強的西風的影響，部分高海拔冰川區年降水量可達2000毫米以上^[2]。

公路沿線50公里範圍內聚集了大小冰川2500餘條，其中包括發育於喬戈里峰的希斯巴冰川、發育於巴托拉峰的巴托拉冰川、帕蘇冰川等大型高物質平衡冰川^[3]。公路常年處於冰川災害的威脅之下，僅2008年一年固爾金冰川的加速前進就誘發了4次冰湖潰決型泥石流，導致公路嚴重受損^[4]。Amaury Dehecq在2015年使用光學影像獲取了喀喇崑崙山地區冰川表面流速，結果顯示區域內大部分冰川運動情況穩定，但存在部分較大規模冰川表面流速較快^[5]。許多研究表明該區域冰川運動存在異常^[6][7][8]，不同學者在不同時期的研究結果有一定出入^[9][5]，也說明了這一區域內冰川運動特徵的複雜性。

遙感手段憑藉大範圍長時間的對地觀測能力，已逐漸取代了GNSS（Global Navigation Satellite System，全球導航衛星系統定位）等傳統觀測手段對冰川等惡劣氣候環境區域的地面監測任務。而SAR（Synthetic Aperture Radar，合成孔徑雷達）技術經過幾十年的發展，已經具備了多波段、多極化、多模式的特性，能夠全天候、不受天氣影響地獲取較高時空分辨率的地面雷達反射數據^[10][11]，有利於開展大範圍的冰川表面運動監測。

由於間隔較長時間的影像對存在較嚴重的時間去相干，而偏移量跟蹤方法採用逐像素的匹配算法進行特徵點搜索，可以有效避免時間去相干帶來的影響。近年來，許多學者在全球許多區域使用偏移量跟蹤的方法進行了冰川表面流速的提取，Strozzi等人在2002年曾應用SAR偏移量跟蹤方法提取了冰川表面流速，並對方法的可行性與結果精度進行了論證，證明該方法是可行的^[12][13]。本文獲取了2014年10月至2018年3月間63期覆蓋中巴公路沿線區域的Sentinel1A數據，使用偏移量跟蹤方法對區域內的山嶽冰川進行了表面運動特徵分析，得到了58對連續的冰川表面運動場數據。

1 數據採集和處理方法

1.1 數據源

Sentinel1A是ESA（European Space Agency，歐洲航天局）研製的全天候、全球覆蓋的搭載了合成孔徑雷達的對地觀測成像系統。由圖1可以看出，研究區域南北狹長，需要縱向分布的南北中三景Sentinel1A影像才能完全覆蓋。本文選取了27號Path下的相鄰三個位置的IW模式、VV極化方式數據，每個位置獲取了2014年10月至2018年3月之間的63期數據，共計189景影像。原始數據具體情況如表1。

圖片

圖1 中巴公路及公路沿線冰川位置示意圖

表1 Sentinel1A數據使用情況

影像位置	影像對時間	影像數量	時間基線（天）	極化方式	處理模式
Path 27 Frame 121	2014-10-06至2016-12-31	30	24	VV	IW
	2017-01-24至2018-03-08	33	12
Path 27 Frame 117	2014-10-06至2016-12-31	30	24
	2017-01-24至2018-03-08	33	12
Path 27 Frame 112	2014-10-06至2016-12-31	30	24
	2017-01-24至2018-03-08	33	12

'註：VV：垂直同向極化方式；IW：Interferometric Wide swath, 干涉寬幅模式；'其中數據對20150710_20150827時間基線為48天，20150827_20151002的時間基線為36天。

另外，為了對SAR數據和計算結果進行基於地形的精校正，本文還獲取了NASA（美國宇航局）和USGS（美國地質調查與測繪局）發布的30 m分辨率SRTM（Shuttle Radar Topography Mission）數字高程模型，以獲取研究區域內的高程和坡度。

1.2 數據處理方法和流程

冰川具有流動性特徵，尤其是山嶽冰川，處於躍動期時流速較快，局部地區表面冰流速可達數百米每年。因此，冰川區域的兩幅SAR影像之間的去相干性較強。本文採用特徵匹配的方法獲取兩幅SAR影像的特徵點偏移量。

SAR特徵匹配方法的原理與光學圖像特徵匹配類似，根據設定的搜索窗口進行全局特徵搜索，得到特徵區域和特徵點的匹配對，根據設定的偏移距離閾值，將距離過遠的匹配對刪除。然後使用偏移量跟蹤的方法根據圖像紋理的最大相關係數計算得到特徵點在方位向和距離向的位移分量，再根據兩個方向上的位移量合成圖像特徵點的二位平面位移矢量^[14][15]，即圖像特徵點的偏移方向和偏移方向上的偏移量。配准偏移值包括：相同地表特徵在不同時間的運動位移量\({R}_{motion}\)，軌道偏移量\({R}_{orbit}\)（同一衛星獲取相同地物影像時重複軌道不夠精確導致的軌道偏移）和地形起伏引起的偏移\({R}_{dem}\)。配准偏移量\({R}_{offset}\)的構成可以表示為：

\({R}_{offset}={R}_{motion}+{R}_{orbit}+{R}_{dem}\) （1）

在基線較短和地形平坦的影像對中地形起伏引起的偏移量很小，可忽略不計。那麼\({R}_{offset}\)的公式就可以表示為：

\({R}_{offset}={R}_{motion}+{R}_{orbit}\) （2）

在本文研究區域內可以認為在短時間基線下除了冰川及冰雪覆蓋區域，其他地物是不存在位移的。為獲取冰川位移需要將軌道偏移量去除，對匹配窗口進行過採樣，並使用最小二乘函數擬合軌道偏移量：

\(R={a}_{0}+{a}_{1}x+{a}_{2}y+{a}_{3}xy\) （3）

其中\({a}_{0}，{a}_{1}，{a}_{2}，{a}_{3}\)是待定係數；\(x\)表示SAR坐標系下的距離向坐標值；\(y\)表示SAR坐標系下的方位向的坐標值。擬合根據內插確定的峰值來模擬相關函數。

本文採用瑞士GAMMA軟件平台和美國ESRI公司的ArcGIS軟件平台進行數據的處理和分析，如圖2所示。在計算之前，使用SRTM數據對SAR圖形進行基於地形的精配准，對SAR影像中由於地形起伏導致的誤差進行補償。而後，為了提高精度，首先目視解譯出SAR圖像的非冰川區域進行採樣，然後以128×256作為搜索窗口進行圖像特徵點和特徵區域的匹配，再使用偏移量跟蹤方法計算冰川區域的偏移向量。

圖片

圖2 計算冰川表面流速流程

數據處理流程主要包括預處理、影像配准、偏移量估計等步驟。預處理即將Sentinel1A的L1數據和SRTM DEM數據處理為適合GAMMA軟件平台的格式。影像配準是利用SAR影像配准原理獲取地表二維形變引起的像素偏移，為精確偏移量信息要藉助SRTM DEM進行地形校正編碼。然後基於軌道參數全局位移對主副影像進行粗配准，獲取初始偏移量，精度為0.5個像素。再根據初始偏移量依據像素級的配准原理設置一定尺寸的窗口和搜索步長在影像上滑動，通過計算相干性的方法獲取相同地表特徵的精確偏移量。本文使用的匹配窗口大小為128×256。

偏移量的估計是利用兩個時間段兩幅影像的特徵匹配結果，採用最小二乘法擬合相對軌道偏移的多項式，最後的標準差可以視為配准精度。將計算的配准偏移量去除擬合的軌道偏移量即為冰川表面的位移量。由於Sentinel1A的對地觀測方式為正攝測量，因此前述過程計算得到的結果可視為水平方向的偏移向量，在這裡引入坡度數據便可以計算得到冰川區域真實的三維方向運動矢量。

2 數據樣本

2015–2017年中巴公路沿線山地冰川表面運動數據集包含兩部分數據，一是研究區域內冰川表面運動方向數據，命名規則為YYYYMMDD_YYYYMMDD_direction.tif，表示兩個時間點間冰川表面的運動方向，以弧度值表示，存放在Direction文件夾中；二是研究區域內冰川表面運動速度數據，命名規則為YYYYMMDD_YYYYMMDD_velocity.tif，表示兩個時間點之間冰川表面的運動速度，單位為m/d，存放在Velocity文件夾中。樣本展示如圖3所示，圖3(a)是2017年5月帕米爾地區冰川流速，圖3(b)是2017年5月巴托拉峰地區冰川流速。

圖片（a）

圖片（b）

圖3 中巴公路沿線局部區域冰川表面運動分布圖

3 數據質量控制和評估

本數據集採用的Sentinel1A的IW模式數據分辨率為5m×20m，由分辨率導致的誤差是無法避免的。另外，在數據的計算過程中，例如影像配准也會產生誤差。可以認為誤差主要來自圖像對的偏移估計，其誤差大約等於距離向和方位向的誤差之和^[16]。本文使用的瑞士GAMMA雷達數據處理平台的影像配准算法誤差控制在0.01個像素範圍內，對應Sentinel1A距離向和方位向分別約為0.018 m和0.009 m，這個數量級的誤差對於冰川運動監測的影響有限。

本文採用的SAR圖像特徵匹配方法的誤差來源除了前文提到的GAMMA軟件系統誤差外，還包括影像配准誤差、地形引起地偏移、地理編碼導致的誤差以及軌道誤差等。為了量化誤差大小，假定非冰川區域在24天這樣較短的時間基線下是靜止的，因此結果中非冰川區域的偏移量就包含了上述所有的誤差源。本文對計算結果中的非冰川區域進行採樣統計，結果顯示整體非冰川區域的方差為0.0445 m/d，具體如表2所示。圖4則是採樣總體誤差的直方圖，可以看出本數據集的結果精度對於100 m×100 m的空間分辨率而言是足夠的。本文在進行採樣誤差分析時，假定非冰川區域均為靜止區域，然而事實上非冰川區域中還存在小規模水域的面積變化、冬春季節山腳積雪融化等導致的非冰川區域運動，因此通過對非冰川區域進行採樣統計的誤差可以認為是數據結果的最大誤差。

表2 非冰川區域誤差統計表（m/d）

數據對	方差	數據對	方差	數據對	方差
20141031_20141124	0.0540	20160809_20160902	0.0048	20170711_20170723	0.0835
20141218_20150111	0.0183	20160902_20160926	0.0042	20170723_20170804	0.1271
20150111_20150204	0.0414	20160926_20161020	0.0017	20170804_20170816	0.0744
20150312_20150405	0.0662	20161020_20161113	0.0012	20170816_20170828	0.0632
20150405_20150429	0.0917	20161113_20161207	0.0017	20170828_20170909	0.0606
20150429_20150523	0.0785	20161207_20161231	0.0013	20170909_20170921	0.0660
20150523_20150616	0.0680	20161231_20170124	0.0029	20170921_20171003	0.0531
20150616_20150710	0.0833	20170124_20170205	0.0202	20171015_20171027	0.0698
20150710_20150827	0.0300	20170205_20170217	0.0079	20171027_20171108	0.0630
20150827_20151002	0.0419	20170217_20170301	0.0101	20171108_20171120	0.0604
20151002_20151119	0.0406	20170301_20170313	0.0175	20171120_20171202	0.0864
20151213_20160106	0.0466	20170313_20170325	0.0167	20171202_20171214	0.1005
20160106_20160130	0.0022	20170325_20170406	0.0543	20171214_20171226	0.0622
20160130_20160223	0.0019	20170406_20170418	0.1159	20171226_20180107	0.0610
20160223_20160318	0.0032	20170418_20170430	0.0829	20180107_20180119	0.0558
20160318_20160411	0.0176	20170430_20170512	0.0406	20180119_20180212	0.0216
20160411_20160505	0.0346	20170512_20170524	0.0443	20180212_20180224	0.0676
20160505_20160529	0.0119	20170524_20170605	0.0206	20180224_20180308	0.0786
20160529_20160716	0.0041	20170605_20170629	0.0556	平均	0.0445
20160716_20160809	0.0015	20170629_20170711	0.0816

由於研究區域內極高山峰林立，因此存在部分雷達傳感器無法觀測到的疊掩和陰影區域，這部分區域的數據結果是不可信的。本文利用視域分析方法得到了疊掩和陰影區域的範圍，並將計算結果中這部分數據剔除。此外，春末夏初的部分結果中冰川的積累區出現了一些規模較小的數據空洞，是由於該區域圖像特徵點匹配對的偏移距離超過了設置的偏移量閾值而導致的失配准。春末夏初時地表氣溫變化劇烈，風力大，易發生雪崩，因此這部分數據在海拔落差較大的冰川積累區出現了小規模的數據空洞。

圖片

圖4 表面運動速度在非冰川區域的誤差

4 數據價值

Sentinel1A採用了C波段的合成孔徑雷達傳感器，對於大氣干擾有較好的屏蔽效果，加之Sentinel1A數據時間間隔較短，因此得到的冰川運動數據的精度和時間分辨率都比較高。本數據集還結合外部DEM數據進行了地形偏移補償，進而得到了中巴公路沿線冰川表面運動分布數據。冰川表面運動數據是研究冰川物質平衡，冰川運動特徵，以及冰川災害預警的重要基礎數據。對冰川表面運動數據集的充分利用，結合氣候、水文數據的綜合探討，可以對區域地理環境變化的特點和趨勢做出全面的分析；另外冰川物質平衡的分析結果，對區域乃至全球水資源分布和變化的研究有着重要意義。

5 數據使用方法和建議

2015–2017年中巴公路沿線山地冰川表面運動數據集通過CERN綜合中心數據資源服務網站提供全部數據下載服務。本數據集提供了ArcGIS軟件平台下生產的32位浮點型GeoTIFF柵格數據，除了*.tif數據本身之外，還有包括了影像坐標文件*.tfw和投影、金字塔信息*.tfw.aux.xml。其中流速數據以影像的像元值表示，單位為m/d；流向數據則以弧度值存儲在影像像元值中。可以在ArcGIS、SuperMap等GIS平台軟件直接調用，也可以使用Matlab等數據分析軟件編碼進行批量處理和分析。

致 謝

感謝ESA提供的Sentinel1A數據，感謝CGIAR-CSI（國際農業研究磋商組織）提供SRTM 30m DEM 數據，感謝GLIMS（全球陸地冰川空間測量）及中國科學院寒區旱區研究所科學數據中心提供冰川編目數據集。

參考文獻

1. 張祥松, 陳建明, 王文穎. 喀喇崑崙山巴托拉冰川的新近變化 [J]. 冰川凍土, 1996, 18(Special Issue): 13.
2. 施雅風, 張祥松. 喀喇崑崙山巴托拉冰川的近代進退歷史變化 [J]. 地理學報, 1978, 33(1): 16.
3. 施雅風 丁德文, 付連弟. 巴托拉冰川末端動態數值預報 [J]. 冰川凍土, 1982, 4(1): 10.
4. 廖麗萍, 朱穎彥, 楊志全, 等. 中國_巴基斯坦喀喇崑崙公路Ghulkin冰川百年進退變化 [J]. 冰川凍土, 2013, 35(6): 9.
5. ^ ^{5.0 5.1} DEHECQ A, GOURMELEN N, TROUVE E. Deriving large-scale glacier velocities from a complete satellite archive: Application to the Pamir–Karakoram–Himalaya [J]. Remote Sensing of Environment, 2015, 162: 55-66.
6. MUHAMMAD S, TIAN L. Changes in the ablation zones of glaciers in the western Himalaya and the Karakoram between 1972 and 2015 [J]. Remote Sensing of Environment, 2016, 187(505-12.
7. ASHRAF A, NAZ R, IQBAL M B. Altitudinal dynamics of glacial lakes under changing climate in the Hindu Kush, Karakoram, and Himalaya ranges [J]. Geomorphology, 2017, 283(72-9.
8. AKHTAR M, AHMAD N, BOOIJ M J. The impact of climate change on the water resources of Hindukush–Karakorum–Himalaya region under different glacier coverage scenarios [J]. Journal of Hydrology, 2008, 355(1-4): 148-63.
9. KHROMOVA T E, OSIPOVA G B, TSVETKOV D G, et al. Changes in glacier extent in the eastern Pamir, Central Asia, determined from historical data and ASTER imagery [J]. Remote Sensing of Environment, 2006, 102(1-2): 24-32.
10. 張震, 劉時銀, 魏俊鋒, 等. 新疆帕米爾躍動冰川遙感監測研究 [J]. 冰川凍土, 2016, 38(1): 11-20.
11. 楊魁, 楊建兵, 江冰茹. Sentinel_1衛星綜述 [J]. 城市勘測, 2015, 2): 24-7.
12. LI R, MA X, CHENG Y, et al., 2194-9034 [R]: ISPRS - International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2017.
13. SHANGGUAN D, LIU S, DING Y, et al. Characterizing the May 2015 Karayaylak Glacier surge in the eastern Pamir Plateau using remote sensing [J]. Journal of Glaciology, 2016, 62(235): 944-53.
14. 蔣宗立, 劉時銀, 許君利, 等. 應用SAR特徵匹配方法估計音蘇蓋提冰川表面流速 [J]. 冰川凍土, 2011, 33(3): 512-8.
15. 李佳, 李志偉, 汪長城, 朱建軍, 丁曉利. SAR偏移量跟蹤技術估計天山南依內里切克冰川運動 [J]. 地球物理學報, 2013, 56(4): 11.
16. QURESHI M A, YI C, XU X, et al. Glacier status during the period 1973–2014 in the Hunza Basin, Western Karakoram [J]. Quaternary International, 2017, 444: 125-36.

數據引用格式

王磊, 蔣宗立, 劉時銀, 張震. 2015–2017年中巴公路沿線山地冰川表面運動數據集[DB/OL]. Science Data Bank, 2018. (2018-07-02). DOI: 10.11922/sciencedb.628.