2015–2017年中巴公路沿线山地冰川表面运动数据集 中华文库

摘要：中巴公路是一条连接中国西部与巴基斯坦及南亚、西亚的重要国际公路，对于沿线国家及高亚洲地区经济文化有着非常重要的意义。而公路沿线的帕米尔高原和喀喇昆仑山脉是全球山岳冰川分布最为密集的区域之一，冰川的运动时常对公路的安全运营造成威胁。本文获取了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: 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

中巴公路（喀喇昆仑国际公路，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 数据源

Sentinel1A是ESA（European Space Agency，欧洲航天局）研制的全天候、全球覆盖的搭载了合成孔径雷达的对地观测成像系统。由图1可以看出，研究区域南北狭长，需要纵向分布的南北中三景Sentinel1A影像才能完全覆盖。本文选取了27号Path下的相邻三个位置的IW模式、VV极化方式数据，每个位置获取了2014年10月至2018年3月之间的63期数据，共计189景影像。原始数据具体情况如表1。

图1 中巴公路及公路沿线冰川位置示意图

表1 Sentinel1A数据使用情况

'注：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的对地观测方式为正摄测量，因此前述过程计算得到的结果可视为水平方向的偏移向量，在这里引入坡度数据便可以计算得到冰川区域真实的三维方向运动矢量。

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 中巴公路沿线局部区域冰川表面运动分布图

本数据集采用的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）

由于研究区域内极高山峰林立，因此存在部分雷达传感器无法观测到的叠掩和阴影区域，这部分区域的数据结果是不可信的。本文利用视域分析方法得到了叠掩和阴影区域的范围，并将计算结果中这部分数据剔除。此外，春末夏初的部分结果中冰川的积累区出现了一些规模较小的数据空洞，是由于该区域图像特征点匹配对的偏移距离超过了设置的偏移量阈值而导致的失配准。春末夏初时地表气温变化剧烈，风力大，易发生雪崩，因此这部分数据在海拔落差较大的冰川积累区出现了小规模的数据空洞。

图4 表面运动速度在非冰川区域的误差

Sentinel1A采用了C波段的合成孔径雷达传感器，对于大气干扰有较好的屏蔽效果，加之Sentinel1A数据时间间隔较短，因此得到的冰川运动数据的精度和时间分辨率都比较高。本数据集还结合外部DEM数据进行了地形偏移补偿，进而得到了中巴公路沿线冰川表面运动分布数据。冰川表面运动数据是研究冰川物质平衡，冰川运动特征，以及冰川灾害预警的重要基础数据。对冰川表面运动数据集的充分利用，结合气候、水文数据的综合探讨，可以对区域地理环境变化的特点和趋势做出全面的分析；另外冰川物质平衡的分析结果，对区域乃至全球水资源分布和变化的研究有着重要意义。

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（全球陆地冰川空间测量）及中国科学院寒区旱区研究所科学数据中心提供冰川编目数据集。

王磊, 蒋宗立, 刘时银, 张震. 2015–2017年中巴公路沿线山地冰川表面运动数据集[DB/OL]. Science Data Bank, 2018. (2018-07-02). DOI: 10.11922/sciencedb.628.