0 引言

冰川是气候变化的指示器, 也是许多地区重要的水资源来源以及海平面上升的重要贡献者[1-3]。物质平衡是衡量冰川“健康”状况的重要指标, 但是由于山地冰川处于高海拔区域、 冰川类型多样、 通达性比较差, 传统的地面调查方法无法完成大范围冰川变化监测[4], 全球具有长时间尺度物质平衡观测的冰川非常少, 而长期的物质平衡观测也仅局限于一些较易到达的冰川区。随着对地观测技术的不断发展, 我们在高程数据的获取方面取得了很大的进步[5]。近年来, 高分辨率空间观测卫星(SPOT-5、 ALOS、 Cartosat-1等)已被广泛应用于冰川地形的测量, 为山地冰川物质平衡的估算[6-8]和冰湖体积的快速计算[9]提供了基础的DEM数据源。但由于受到技术的限制, 早期的遥感卫星大多不具备侧视成像功能, 相邻像对的基高比较小, 提取的DEM高程精度较低, 使当前大多数利用遥感卫星立体像对DEM数据获取冰川物质平衡的研究, 其时间跨度小于10年(主要集中在2000年以后), 因此较为缺乏早期的长时序的山地冰川物质平衡观测资料, 而这些资料对于认识和理解区域内冰川历史演变和气候变化具有非常重要的作用和意义[10]。

美国地质调查局(The US Geological Survey, USGS)在2002年解密了早期高分辨率照相侦察卫星KH-9“大鸟”, 其成像时间早(1971-1980年), 空间分辨率高, 地面覆盖范围广, 基高比大, 可为山地冰川长时序的物质平衡研究提供很好的基础数据源, 具有重要的研究价值[11-12]。目前, 利用KH-9 数据进行的冰川、 冰湖研究已经在整个青藏高原以及天山等区域取得较好的研究进展。上官冬辉等[9]结合KH-9、 SRTM、 SPOT-5和ALOS DEMs数据估算了天山默兹巴克湖泊2000-2008年的储水量, 建立了湖泊面积和湖泊体积的关系, 并利用该关系推算出了默兹巴克湖泊的补给率。Holzer等[13]利用1973年的KH-9数据和2013年的Pléiades数据获取了慕士塔格冰川近40年间的质量变化, 结果为(-0.01±0.30) m w.e.·a-1。表明该冰川在1973-2013年期间基本处于物质平衡状态。周玉杉等[14]利用1975年的KH-9数据和2000年SRTM1数据对帕米尔中部冰川的物质平衡进行了估计, 结果表明, 区域内冰川在1975-2000年期间的物质平衡为(0.04±0.04) m w.e.·a-1, 对应于(-0.08±0.08) Gt·a-1的物质流失, 整体上处于轻微的负物质平衡或平衡的状态。然而, 由于在偏远的山地冰川区缺乏高精度的地面高程数据, 很少有学者对生产的KH-9 DEM数据进行严格的精度评价、 精度影响分析和建立完整的KH-9 DEM提取技术流程。鉴于此, 本文利用KH-9立体像对, 在影像畸变校正、 影像拼接、 影像对比度增强、 DEM的提取与编辑的基础上, 建立了一套详细完整的DEM提取技术流程, 并利用ICESat GLAS高程数据对提取的DEM数据进行了精度评价与精度影响因素分析。

1研究区域

我们选取了两个研究区域, 以评估在不同地形条件下KH-9数据获取山地冰川DEM的能力。第一个研究区域位于西藏北部的普若岗日冰川(40 km×40 km； 中心经纬度： 33°45′ N, 89°22′ E； 图1左图黑框内), 该冰川所处区域的地形较为平坦, 高程变化范围为4 890～6 490 m。第二个研究区域位于西藏东南部的雅弄冰川(40 km×40 km； 中心经纬度： 29°10′ N, 96°29′ E； 图1右图黑框内), 该冰川周围多为高山地形, 高程变化范围为1 960～6 750 m。

2实验数据

本文所用数据包括： KH-9光学数据、 Landsat-7 ETM+ L1T遥感数据、 SRTM DEM和ICESat GLAS高程数据。所有影像数据的地理坐标系为WGS84, 采用UTM投影, 水准面模型为EGM96。数据详见表1。

2.1KH-9 Hexagon

KH-9 Hexagon是一种照相侦察卫星, 美国于1971年到1980年间共发射了20颗KH-9系列卫星(任务期编号： 1201～1220)。其携带的测图相机与美国宇航局在1984年使用的大像幅摄影机(Large Format Camera, LFC)相似, 为23 cm×46 cm的帧格式相机, 焦距为30.5 cm[11]。2002年USGS解密了KH-9卫星在12个任务期间(编号： 1205～1216)所拍摄的约29 000张影像。这些影像的像元大小为7 μm(3 600 dpi), 以8-bit TIFF格式存储, 存储大小为2.6 Gb, 拍摄高度为171 km, 影像中共覆盖了4个框标和1 058个十字标, 地面覆盖面积为250 km×125 km, 分辨率6～9 m[12]。本文共使用了6幅KH-9影像, 覆盖普若岗日冰川区的影像任务编号为1211, 覆盖雅弄冰川区域的影像任务编号为1209(表1, 图1)。

2.2Landsat-7 ETM+(L1T)

美国陆地7号卫星(Landsat-7)于1999年4月15日由美国航空航天局(NASA)发射升空, 其所携带的主要传感器为增强型主题成像仪ETM+ (Enhanced Thematic Mapper Plus)。2003年后, 由于Landsat-7的扫描行校正器发生了故障, 为了避免影像条带的影响且为了与SRTM获取的时间保持一致, 本文选择了获取于2000年前后的Landsat-7 L1T标准地形校正全色影像作为控制点的水平参考(表1)。

2.3SRTM DEM

SRTM DEM数据是在2000年2月由航天飞机雷达地形测绘任务(SRTM)、 美国国家航空航天局(NASA)、 美国国家测绘局(NIMA)和德国宇航中心(DLA)联合测量的。该数据覆盖了全球80%以上的陆地表面(60° N～56° S), 其水平精度优于45 m, 垂直精度优于10 m。可作为控制点的垂直参考(表1)。

2.4ICESat GLAS

ICESat(冰、 云、 陆地高程卫星)发射于2003年1月13日, 2009年10月11日停止运行。ICESat按近极地轨道飞行, 轨道高度为600 km, 其搭载的地学激光测高系统(GLAS)携带了3个激光测高计, 使用1 064 nm脉冲测量地面和浓云的高度, 发射频率为40 Hz, 光斑直径大约65 m, 相邻光斑间隔172 m[16-18]。本文使用了美国国家雪冰数据中心(NSIDC)最新发布的数据产品ICESat GLA14 (release 634), 利用NSIDC提供的NGAT程序提取了光斑的经纬度和高程信息, 其水平精度为(10.6±4.5) m, 垂直精度可达34 cm[19], 可作为控制点的垂直参考和DEM精度评价基准(表1)。

3数据处理

3.1KH-9数据的预处理

KH-9数据的预处理主要包括： 十字标识别、 畸变校正、 图像拼接、 图像增强等步骤(图2)。我们利用图像上覆盖的1 058个精密十字标来对影像的畸变进行校正, 重建影像在曝光时的几何结构。在该步骤中, 首先以亚像素精度获取了KH-9影像中所有十字标的中心坐标(图3a)。然后假设图像中心位置处的十字标没有畸变, 以1 cm的格网间距建立十字标的参考坐标集, 并计算出影像中的畸变十字标与对应位置处的参考十字光标的畸变矢量(图4, 表2)。然后, 基于二次多项式建立了影像十字标与参考十字标的转换关系, 并利用双线性内插算法对影像中的每个像素进行逐一校正, 最终使得所有影像的校正后误差均小于1个像素(表2)。在影像经过几何畸变校正后, 需要对影像中的所有十字标进行消除处理(图3b), 避免其对像对同名特征点的匹配造成干扰。在DEM提取之前, 还需对去除光标后的KH-9子影像进行拼接和对比度增强处理。我们通过特征点捕捉、 关键点自适应非极大值抑制、 基于尺度不变特征变换的关键点描述、关键点匹配、 图像合成等步骤实现图像的无缝拼接。通过51×51窗口的Wallis滤波和直方图均衡化实现了影像的对比度增强(图1)。

3.2DEM的提取

首先对影像进行内定向, 并使得所有影像的内定向中误差均小于0.3个像素(表2)。然后对影像进行绝对定向, 由于缺少RPC参数, 只能通过选取地面控制点来计算外方位元素。但由于在山地冰川区缺少高精度的地面实测数据, 使得控制点的获取变得异常困难。因此, 本文以Landsat ETM+、 SRTM DEM和ICESat遥感数据为参考在相对稳定和平坦的地形区域内选控制点。控制点的x和y坐标是以LandSat ETM+ L1T全色影像作为水平参考进行获取。而控制点的z坐标则主要通过ICESat GLA14高程数据进行获取, 在无ICEsat分布的区域内以SRTM高程数据补充。最终, 我们在普若岗日冰川区内共获取了32个控制点, 在雅弄冰川区域内获取了30个控制点, 并尽量使得所有的控制点均匀地分布在整景影像上(表3)。

在KH-9影像定向处理完成之后, 接着对KH-9立体像对进行空中三角测量解算。如表3所示, 两个研究区域内所有立体像对的总图像单位权RMSE均小于1个像素, 且控制点在水平和垂直方向的残差均小于11 m, 完全满足空中三角测量的精度要求。最后, 我们成功提取了两个研究区的KH-9 DEM, 并将其采样为空间分辨率30 m(图5)。

4DEM的精度评价与分析

4.1DEM的精度评价

DEM精度评定可通过两种不同的途径来进行, 一种是平面精度和高程精度分开评价, 另一种是两种精度同时评价。对于前者, 平面的精度结果可独立于垂直方向的精度结果而获得； 对于后者, 两种精度的获取必须同时进行。而在实际应用中, 由于平面精度比较容易控制, 一般只讨论DEM的高程精度, 并用高程精度表达DEM精度。由于ICESat GLA14为光斑数据, 因此本研究采用检查点法来对DEM进行精度评估。

选择位于非冰川区域且坡度小于15°的ICESat光斑作为检查点, 其分布如图6所示。获取每个检查点所在位置处DEM的高程值。为了更加准确地评价DEM的提取精度, 消除由于传感器不稳定引起的偏移误差, 首先对ICESat检查点和DEM进行了配准。采用Nuth提出的数据配准方法[20], 以ICESat为参考, DEM为错位高程, 提取了感兴趣区域内ICESat与DEM高程差值、 DEM的坡度、 坡向, 通过最小二乘余弦拟合求得平移矢量, 再根据平移矢量对DEM进行平移, 重复上述步骤直至x和y方向的平移量小于1 m。在完成数据的配准之后, 我们计算了所有检查点与DEM的高程差值, 并对其分布进行了统计和正态性分析。如图7所示, 左图为普若岗日区域的DEM误差分布直方图, 其总体分布的峰度为3.25(>0), 与正态分布(峰度=0)相比更加陡峭, 为尖顶峰, 偏度为0.22(>0)较正态分布(偏度=0)右偏, 由此可得, 普若岗日区域的DEM误差不近似服从正态分布； 右图为雅弄区域的DEM误差分布直方图, 其峰度为2.5, 偏度为-0.23, 同理雅弄冰川区域DEM误差分布也不近似服从正态分布。上述非正态性的误差分布主要是由KH-9 DEM中的异常值引起的。为了消除异常值对精度评价的影响, 我们采用了Joachim Höhle提出的DEM精度评价方法[21], 利用统计量： 中值(Med)、 Normalized Median Absolute Deviation (NMAD)和样本分位数作为评价指标对提取的DEM进行精度评价(表4)。其中, NMAD是一个健壮的统计量, 能够弹性的处理数据集中的异常值, 减少异常值对于数据集的影响, 提高DEM精度评价的准确性和可靠性。

NMAD=1.4826×Median(|Δh-mΔh|)(1)

式中： Δh表示每个检查点与DEM的高程差值； mΔh为所有Δh的中值。

由于受地形因素的影响, 两个研究区域的DEM提取精度相差较大(表4)。可以利用一种专门针对摄影测量DEM的概测法[22], 来对不同地形区域的DEM质量进行评估。该方法指出： 在地形平缓的区域, 提取的DEM垂直误差最大应不超过图像分辨率的1～3倍； 而在地形起伏较大的区域, 提取的DEM垂直误差最大应不超过为图像分辨率的5～7倍。因此, 综合考虑KH-9影像的空间分辨率和研究区地形条件, 普若岗日区域的DEM垂直误差应小于18 m, 雅弄区域的DEM垂直误差应小于42 m。很明显, 本文获取的两个区域的KH-9 DEM均完全符合精度要求。

4.2DEM精度影响分析

DEM精度评估结果主要受三方面因素的影响, 一是ICESat检查点的个数和采样方式, 二是控制点的精度和分布, 三是研究区地形特征。在检查点方面, 由于ICESat为极地卫星, 其在中低纬度的轨迹稀疏, 且在不同区域内光斑点的数量及其分布相差较大(图6, 表4)。因此, 在选取检查点时我们只能按照ICESat的轨迹选取, 而不能满足检查点的均匀分布, 这将不可避免的影响了DEM的精度评价结果； 在影像控制点的选取方面, 由于KH-9影像成像时间早(1971-1980年), 在影像拍摄时, 研究区域内没有道路交叉路口、 房屋等稳定的人工建筑物可作为控制点。我们只能根据不同的地形地貌特征, 在河流口或者是山脊处选取控制点。这导致了最终选取的控制点精度较差且分布不够均匀, 而DEM的质量与控制点的分布具有一定的相关性, 控制点所在范围及其附近地区生成的DEM较为完整光滑, 距离控制点越远, DEM的质量越差并且出现了很明显的异常值； 在研究区地形特征方面, 影响DEM精度的主要因素是由高差过大引起的阴影效应以及积雪覆盖引起的亮度饱和。在雅弄冰川区域, 较大的高差变化使得影像中存在着严重的阴影, 影响了立体像对的精确匹配。而严重的积雪覆盖在影像中表现出差的纹理特征和低的对比度, 导致了特征点的错误匹配, 从而在提取的DEM中引入了大量的误差(图5b)。

5结论与讨论

本文针对KH-9数据, 以普若岗日冰川区和雅弄冰川区为例, 在影像畸变校正、 拼接、 对比度增强、 DEM的提取与编辑等步骤的基础上建立了一套基于KH-9数据的山地冰川区DEM提取技术流程, 并利用ICESat高程数据, 通过检查点法对提取的DEM进行了精度评价。结果表明在普若岗日冰川区域, KH-9 DEM的NMAD为9.23 m； 而在雅弄冰川区域, 由于受到积雪覆盖和阴影的影响, KH-9 DEM的NMAD较普若岗日冰川区域大, 为13.86 m。虽然在不同地形条件下, 提取的DEM精度相差较大, 但与目前广泛应用于山地冰川长时间尺度物质量变化提取的SRTM DEM数据的垂向精度[23]进行对比发现, 我们提取的KH-9 DEM数据均可以很好的满足山地冰川物质平衡估算的精度要求。另外, KH-9数据具有质量稳定、 覆盖广泛、 成像时间早、 空间分辨率高的特点。因此, 基于KH-9 立体像对提取DEM的研究具有重要的实际应用价值和意义, 尤其是在估算山地冰川长时间尺度物质量平衡方面具有重要的作用。

地面控制点的数量、 分布和精度是影响DEM提取精度的最主要因素。增加控制点的数量会提高物方坐标的解算精度, 但是当控制点数量增加到一定程度时, 其对精度的影响会很小。鉴于KH-9影像年代久远且成像于偏远地区, 可作为控制点的特征地物较少, 控制点的选取较为困难, 建议控制点选取的数量为20～30个。控制点的分布与DEM的精度具有一定的相关性, 在选取控制点时应尽量将控制点立体、 均匀地分布在影像区内, 尤其应加强测区边沿的控制点布设。

另外, 应尽可能的利用高分辨率的影像数据作为控制点的坐标参考, 提高控制点的精度, 进而提高DEM的提取精度。

致谢： USGS为本研究提供了KH-9数据、 Landsat-7 ETM+ L1T数据、 SRTM DEM数据。NDSIC为本研究提供了ICESat GLA14高程数据。在此一并表示感谢。

参考文献：

[1] Yang Zhenniang. Glacier water resources of Qilian Mountains[J]. Journal of Glaciology and Geocryology, 1988, 10(1): 36-46. [杨针娘. 祁连山冰川水资源[J]. 冰川动土, 1988 , 10(1): 36-46.]

[2] Li Qiaoyuan, Xie Zhechu, Da Yanan, et al. The variation of glaciers in Nepal during 1980-2010[J]. Journal of Glaciology and Geocryology, 2017, 39(5): 935-948. [李巧媛, 谢自楚, 戴亚南, 等. 尼泊尔冰川1980-2010年的变化特征[J]. 冰川冻土, 2017, 39(5): 935-948.]

[3] Liu Kai, Wang Ninglian, Bai Xiaohua. Variation of glaciers in the Nubra basin, Karakoram Mountains, revealed by remote sensing images during 1993-2015[J]. Journal of Glaciology and Geocryology, 2017, 39(4): 710-719. [刘凯, 王宁练, 白晓华. 1993-2015年喀喇昆仑山努布拉流域冰川变化遥感监测[J]. 冰川冻土, 2017, 39(4): 710-719.]

[4] Zhang Wei, Hang Haidong. A review of the ice volume estimation of mountain glacier[J]. Journal of Glaciology and Geocryology, 2016, 38(6): 1630-1643. [张伟, 韩海东. 山地冰川储量估算方法研究进展[J]. 冰川冻土, 2016, 38(6): 1630-1643.]

[5] Jacob T, Wahr J, Pfeffer W T, et al. Recent contributions of glaciers and ice caps to sea level rise[J]. Nature, 2012, 482(7386): 514-518.

[6] Shangguan Donghui, Bolch T, Ding Yongjian, et al. Mass changes of southern and northern Inylchek Glacier, central Tian Shan, Kyrgyzstan, during～1975 and 2007 derived from remote sensing data[J]. Cryosphere, 2015, 9(2): 703-717.

[7] Bolch T, Pieczonka T, Benn D I. Multi-decadal mass loss of glaciers in the Everest area (Nepal Himalaya) derived from stereo imagery[J]. Cryosphere, 2011, 5(2): 349-358.

[8] Pieczonka T, Bolch T, Buchroithner M. Generation and evaluation of multitemporal digital terrain models of the Mt. Everest area from different optical sensors[J]. Isprs Journal of Photogrammetry & Remote Sensing, 2011, 66(6): 927-940.

[9] Shangguan Donghui, Ding Yongjian, Liu Shiyin, et al. Quick release of internal water storage in a glacier leads to underestimation of the hazard potential of glacial lake outburst floods from Lake Merzbacher in central Tian Shan Mountains[J]. Geophysical Research Letters, 2017, 44(19): 9786-9795.

[10] Berthier E, Arnaud Y, Vincent C, et al. Biases of SRTM in high-mountain areas: implications for the monitoring of glacier volume changes[J]. Geophysical Research Letters, 2006, 33(8): 153-172.

[11] Mollberg B H, Schardt B B. Mission report on the Orbiter Camera Payload System (OCPS) Large Format Camera (LFC) and Attitude Reference System (ARS)[R]. NASA; Washington, DC, United States, 1988.

[12] Day D. Eye in the sky: the story of the CORONA spy satellites[M]. Smithsonian Institution, 2015.

[13] Holzer N, Vijay S, Yao Tandong, et al. Four decades of glacier variations at Muztagh Ata (eastern Pamir): a multi-sensor study including Hexagon KH-9 and Pléiades data[J]. Cryosphere, 2015, 9(6): 2071-2088.

[14] Zhou Yushan, Li Zhiwei, Li Jia. Slight glacier mass loss in the Karakoram region during the 1970s to 2000 revealed by KH-9 images and SRTM DEM[J]. Journal of Glaciology, 2017, 63(238): 331-342.

[15] Farr T G, Kobrick M. The shuttle radar topography mission: a global DEM[J]. American Geophysical Union, 1999, 45(2): 37-55.

[16] Zwally H J, Schutz B, Abdalati W, et al. ICESat′s laser measurements of polar ice, atmosphere, ocean, and land[J]. Journal of Geodynamics, 2002, 34(3/4): 405-445.

[17] Schutz B E, Zwally H J, Shuman C A, et al. Overview of the ICESat mission[J]. Geophysical Research Letters, 2005, 32(21): 97-116.

[18] Magruder L A, Webb C E, Urban T J, et al. ICESat altimetry data product verification at white sands space harbor[J]. IEEE Transactions on Geoscience & Remote Sensing, 2006, 45(1): 147-155.

[19] Wang Yuzhe, Ren Jiawen, Qing Dahe, et al. Using satellite data to invert the change of glacier ice in the region: take Qilian Mountain as an example [J]. Journal of Glaciology and Geocryology, 2013, 35(3): 583-592. [王玉哲, 任贾文, 秦大河, 等. 利用卫星资料反演区域冰川冰量变化的尝试: 以祁连山为例[J]. 冰川冻土, 2013, 35(3): 583-592.]

[20] Nuth C, A. Kääb. Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change[J]. Cryosphere, 2011, 5(1): 271-290.

[21] Höhle J, H hle M. Accuracy assessment of digital elevation models by means of robust statistical methods[J]. Isprs Journal of Photogrammetry & Remote Sensing, 2009, 64(4): 398-406.

[22] Pellikka P, Rees W G. Remote sensing of glaciers: techniques for topographic, spatial and thematic mapping of glaciers[M]. CRC Press, 2009.

[23] Falorni G, Teles V, Vivoni E R, et al. Analysis and characterization of the vertical accuracy of digital elevation models from the Shuttle Radar Topography Mission[J]. Journal of Geophysical Research Earth Surface, 2005, 110(F2): 1-20.