基于高分三号影像的2019–2020 年高亚洲地区典型冰川表面流速数据集

张齐民，吕明阳，闫世勇*

1.中国矿业大学，环境与测绘学院，江苏徐州 221116

2.中国科学院空天信息创新研究院，数字地球重点实验室，北京 100094

引 言

随着全球气候变暖的进程日益加剧，作为气候变化的指示器之一，冰川已成为气候变化研究的关键指标[1-4]。相比于两极冰川，青藏高原地区的冰川对于气候变化的响应更为敏感，其升温速度是同期全球其他区域均值的两倍，冰川加速融化、冰川融水径流量显著增加，整体呈现出失衡的状况，因此研究青藏高原地区冰川流速时空变化对于认识该区域冰川动力学特征及其对气候变化的响应具有重要的意义[5-7]。冰川流速是冰川变化的一个重要参数，能够反映一段时间内冰川的活动性及其物质平衡情况。因而，冰川表面运动时空分布特征对于研究该区域的物质平衡具有重要指示作用[8-11]。

山地冰川所处区域往往气候条件恶劣，实地冰川运动监测成本高且难以在广泛开展。随着科技的进步和遥感技术的发展，基于遥感影像的冰川运动监测已成为研究冰川活动性的最为重要的手段，具有获取数据快、监测范围广、重访周期短等特点[12-15]，主要有光学遥感和雷达遥感两种观测方式。光学影像覆盖范围广，存档数据多，但大气和天气条件等因素对其适用性有较大的限制，高海拔冰川区常年积雪覆盖，易导致光学影像上出现过饱和现象，破坏冰川监测的完整性[9]。相比较而言，合成孔径雷达（SAR）具备全天时全天候的成像能力，不受天气影响且拥有较高的时空分辨率，在冰川活动性监测中有着较强的适用性[16]。

当前青藏高原冰川运动信息的获取绝大多数依赖国外遥感卫星数据，国产卫星遥感影像在冰川运动广域监测中的应用寥寥无几，导致其在冰川活动性研究中的价值难以体现。为了能够及时准确地了解青藏高原地区典型山地冰川的运动状况，凸显我国遥感数据在气候研究中的实际价值，本文选用2019–2020 年获取的高分三号（GF-3）遥感数据影像，基于改进的并行化Pixel-Tracking（PT）技术[17-20]，通过复杂全局形变拟合估计、几何畸变校正和降噪滤波等处理[21-22]，获取了高亚洲地区典型冰川的表面流速，并以非冰川区残余位移量作为统计特征，判别和验证数据集产品的可靠性，最终实现基于国产高质量遥感影像的冰川流速数据集产品生成与应用推广，为青藏高原气候变化研究和冰川运动相关灾害预测提供更精准的数据支持。典型冰川分布如图1 所示，其中冰川边界参照中国第二次冰川编目数据[23]。

图1 典型冰川分布示意图

1 数据采集和处理方法

1.1 数据源

GF-3 卫星是中国自主研制的首颗民用C 波段全极化星载卫星，重访周期为29 天。因其全天时、全天候以及高空间分辨率的特点，在冰川运动监测中有着重要优势。GF-3 卫星具有12 种常规成像模式，最高分辨率可达1 m，最大幅宽约为650 km，可满足不同领域的需要[24]。本文选用2019–2020年覆盖研究区精细条带1 模式(FSI)的SAR 影像进行冰川运动提取，原始数据信息参数见表1，同时采用30 m 分辨率的ASTER DEM 数字表面模型数据辅助SAR 影像配准和地理编码。

表1 GF3 雷达影像对参数信息

1.2 数据处理方法和流程

本研究利用覆盖研究区典型山地冰川的多期GF-3 遥感影像，采用改进的PT 算法高精度监测冰川表面运动的空间分布。基于非冰川区在监测周期内未发生移动的假设，精确地估算出全局轨道相关偏移量、地形起伏相关偏移量和噪声偏移量等无关信息，并对其进行补偿，最终得到冰川区较为可靠的流速结果。PT 算法利用影像强度信息的归一化互相关算法（Normalized Cross Correlation，NCC）进行影像配准，通过计算主副影像模板的相关程度来实现亚像素级配准。该方法精度高，但计算效率低，巨大的运算量对于青藏高原广泛分布的典型冰川区流速提取工作来说十分不利[17]。为此本文利用图形处理器（Graphics Processing Unit，GPU）良好的性能实现了PT 算法的并行化改进，极大提高了运行效率，为冰川表面运动的时序监测提供了技术保障。本研究技术路线如图2 所示。

图2 基于GF-3 影像的冰川流速提取流程示意图

1.2.1 基于RANSAC 的全局形变估计

两次或多次成像过程中，卫星姿态及位置差异会引入轨道相关的全局性偏移量，一般采用二次多项式函数模型进行拟合和补偿[25]。为了降低匹配噪声点和冰川区匹配点给全局模型参数估计造成的不利影响，本文假设非冰川区域在观测周期内未发生位移且面积占比较大，利用随机抽样一致性算法（Random Sample Consensus，RANSAC），通过不断迭代确定符合模型的匹配点，并重新计算变换参数，最终在满足阈值条件时终止迭代计算[26]，由此可以获得变换模型参数并去除轨道相关的偏移量。该方法能够在提高计算效率的同时保证全局形变参数估计的精确性，从而准确补偿因卫星轨道及传感器姿态差异而导致的整体性位移。

1.2.2 基于外部DEM 的地形校正

研究区地形复杂多变，冰川分布区高差大，卫星对冰川区成像时的位置差异会引入与地形相关的偏移量，进而降低冰川运动监测的准确性[27]。一般而言，方位向的地形偏移量主要与所在区域的地形起伏、成像时轨道交角及卫星入射角有关，而距离向地形起伏导致的偏移量则与垂直基线长度呈正比[21]。基于外部DEM 数据和两景SAR 影像的空间基线情况，能够准确获取地形起伏在距离向和方位向所产生的偏移，对全局形变校正后的偏移量进行精确的地形补偿，可进一步提高冰川表面运动估计的准确性。

1.2.3 冰川运动分布滤波与去噪

由于数据质量和处理方法等各因素的限制，经全局形变补偿和地形效应校正后得到的偏移量分布仍然不同程度地受到误配点噪声影响，因此有必要开展冰川运动信息滤波。在冰川运动分布场中，斑点噪声属于高频信息，可以通过低通滤波手段来削弱噪声。然而冰川的边缘信息也属于高频信号，为了保证冰川边缘位移细节信息的完整性和最终结果的准确性，本文采取自适应中值滤波（Rankedorder based Adaptive Median Filter，RAMF）来削弱噪声的影响，借助滤波窗口可变的特点，RAMF既能保证冰川偏移量信息的精度，又能适当地过滤噪声，在冰川运动提取中具有较好的应用价值[17]。

2 数据样本描述

2019–2020 年高亚洲地区典型冰川表面运动分布数据集所用数据为GF-3 FSI 精细条带模式数据，冰川运动结果以GeoTiff 格式分别存储在4 个不同文件夹中，存储格式为32 位浮点型。文件命名由三个部分组成。其中前两个部分形式如下：AREA_YYYYMMDD-YYYYMMDD，AREA 表示冰川名称或区域，YYYYMMDD-YYYYMMDD 表示研究时间段。第三部分表示冰川运动结果的不同类型：VEL 表示冰川表面的平均流速，单位为cm/d，存放在Glacier Velocity 文件夹中；DIS 表示相应冰川表面位移量数据，单位为m，存放于Glacier Motion 文件夹中；DRan 表示相应冰川距离向位移量数据，单位为m，存放于Range Motion 文件夹中；DAzi 表示相应冰川方位向位移量数据，单位为m，存放于Azimuth Motion 文件夹中。图3 为冰川表面运动数据样本展示，(a)(b)(c)(d)分别为South Inylchek 冰川的距离向位移量、方位向位移量、合成方向累积位移量及平均流速示意图。累计位移量单位为m，正负号表示方向；平均流速单位为cm/d，黑色轮廓为冰川边界。

图3 典型冰川表面运动分布示意图

3 数据质量控制和评估

为保证冰川表面运动采样点的可靠性，偏移量跟踪算法中匹配模板和搜索步长等参数的选择和设定需要充分考虑研究区冰川规模、运动特征及其地表后向散射的稳定性。根据研究区的冰川规模大小，分别设定为：模板窗口128×128 像素，采样窗口5×5 像素。搜索窗口的选择需要根据不同区域冰川的运动量来进行相应的调整，从而保证匹配结果具有较高可靠性，进而获取冰川区完整的位移量分布信息。

为了得到更为可靠的偏移量结果，本文以模板匹配过程中的统计特性对GF-3 影像对中目标点进行质量评估，采用归一化互相关系数和协方差作为判别条件，反映相应区域的匹配质量。此外，本文基于RANSAC 算法的全局形变估计，有效改正了因轨道和卫星成像姿态差异带来的偏差。同时，给出基于外部DEM 的地形误差校正方法和相应数据处理流程，以此来弥补复杂地形带来的误差。高亚洲地区地形复杂，气候多变，常年受积雪覆盖和降水等因素的影响，偏移量跟踪结果不可避免地受到噪声干扰，本文采用RAMF 算法进行滤波处理，可以有效削弱噪声对于结果的影响，提升冰川运动监测结果的质量。

冰川区气候条件恶劣，高寒高海拔交通不便，难以广泛开展冰川流速野外验证。因此，为进一步验证冰川运动结果的可靠性，本文基于非冰川区在监测周期内稳定的合理假设，利用偏移量跟踪结果中非冰川区的残余位移量进行精度分析[24,28]。本文不同研究区内非冰川覆盖区域的残差统计直方图如图4 所示，均值及方差见表2。非冰川区残差统计分析表明，在研究周期内冰川运动总体精度约为0.5 m。该精度已基本满足实际需要，具有一定的可靠性和准确性。

图4 非冰川区残余误差直方图

表2 非冰川区误差统计表（m）

4 数据价值

随着对地观测技术的不断发展，遥感已经成为研究冰川性质及其动态变化特征的重要技术手段。当前，青藏高原冰川运动信息数据集多来自国外遥感卫星影像，我国高质量遥感影像在该区域的应用十分稀少，未能充分挖掘和体现国产遥感数据在冰川活动性研究中的应用价值。本文基于C 波段GF-3 雷达影像的FSI 模式数据提取了青藏高原地区典型山地冰川表面2019–2020 年运动分布，一定程度上填补了青藏高原地区没有国产遥感卫星广域冰川流速数据集的空白，凸显出我国遥感数据在气候研究中的应用价值。此外，通过全面了解冰川运动情况，能更好地开展我国冰川资源的调查，为合理利用青藏高原地区冰川资源提供必要的基础信息。

5 数据使用方法和建议

2019–2020 年高亚洲地区典型冰川表面运动分布结果数据存储格式为32 位浮点型 GeoTiff 栅格数据，地理坐标系为WGS-1984，ArcMap 等GIS 软件或ENVI、PIE 等软件均支持该数据的读取和操作，其中冰川流速数据单位为cm/d，累计位移量数据单位为m。研究人员可根据需要灵活选择数据，也可根据当前或后期的冰川运动监测结果，开展青藏高原典型冰川区表面应力分布分析与研究工作。

致 谢

感谢中国科学院空天信息创新研究院（AIR）提供的GF-3 数据，美国航空航天局（NASA）网站提供的ASTER 30 m DEM 数据，中国科学院寒区旱区环境与工程研究所提供的中国第二次冰川编目数据集。