Landsat 8全色影像阴影的扁平冰山出水高度提取

关真富，程晓，2，3，4，刘岩，3，4，璩榆桐，李腾

1.北京师范大学 全球变化与地球系统科学研究院遥感科学国家重点实验室,北京 100875;2.中山大学测绘科学与技术学院,珠海 519082;3.南方海洋科学与工程广东省实验室(珠海), 珠海 519082;4.中国高校极地联合研究中心,北京 100875

1 引 言

随着全球气候变暖，冰架崩解事件的发生愈益频繁（刘岩 等，2013；Liu 等，2015）。冰架崩解产生冰山，冰山作为南大洋中重要的移动淡水源，随着洋流漂移和融化，影响上层海洋层次划分的稳定（Gladstone 等，2001）。冰山数目的增加或者冰架崩解频率的增加，会使得更多的淡水进入海洋进而影响海洋的洋流循环（Silva 等，2006），从而对南大洋的水文、生态系统产生重要的影响（Schwarz 和Schodlok，2009）。因此，准确测量冰山体积对精确评估冰山对海洋的淡水输入量非常重要。

遥感是测量冰山面积和冰山厚度进而计算冰山体积的主要手段（冯准准等，2013）。冰山厚度可通过冰雷达直接测量，也可通过测量冰山出水高度间接获取。使用具有冰穿透能力的机载雷达是直接测量冰山结构和厚度的一种重要方式，但这种方法有强烈的时空上的局限性，难以大范围实施（Blankenship 等，2002）。ICESat、ICESat-2、CryoSat-2 和Jason 系列测高卫星数据是精确测量冰山出水高度的另一种重要手段（Zwally 等，2002；Dibarboure 等，2011；王显威 等，2013；Markus等，2017；Tournadre 等，2018），例如，Wang 等（2014）首次使用ICESat/GLAS 时间序列数据生成默茨冰舌崩解冰山的出水高度信息图；Tournadre等（2015）利用Jason-1/2 雷达高度计波形信息计算了南极冰山平均出水高度；Li 等（2018）利用CryoSat-2 数据计算C28A 和C28B 扁平冰山出水高度变化。但是，由于受到测高卫星仅沿轨间隔采样特征的限制，很难覆盖到南极所有的冰山，特别是较小的冰山。

南极冬季沿岸海面覆盖大量海冰，海冰表面非常平坦可以近似为平面，此时太阳高度角很低，被海冰包围的冰山在海冰表面可形成很长的阴影，这使得利用阴影测量冰山出水高度成为一种可能，Li 等（2019）在利用UAV 数据评价海冰表面特征时简单讨论了低太阳高度角下冰山阴影与其高度的关系，但国内外对此尚无系统全面的研究。阴影测高方法目前主要用于建筑物高度提取，其利用建筑物及阴影成像几何模型测量建筑物高度，所采用的数据以高分辨率影像为主。国内外许多学者对此进行了研究，例如Irvin 和McKeown（1989）使用航空影像建立阴影与建筑物的关系；Cheng 和Thiel（1995）、Shettigara 和Sumerling（1998）及何国金等（2001）利用SPOT 全色图像中的建筑物阴影直接估算建筑物的高度取得了较好的结果；谢军飞和李延明（2004）分析了IKONOS卫星图像上建筑物阴影与实际高度的几何关系；冉琼等（2008）利用“北京一号”小卫星影像建立了建筑物高度测量的模型，测高精度达到亚像元级；张晓美等（2011）建立了ALOS卫星阴影测高模型。

冰山与建筑物等固定不动的物体有所差别，由于消融的影响，冰山的出水高度是变化的。例如，Li 等（2018）利用多期CryoSat-2 雷达高度计数据测量了南极默茨冰山出水高度的变化率，夏季约为0.7 m/月，冬季约为0.8 m/月，夏季与冬季的出水高度变化率很低且相近。因此，由于冰山出水高度的季节性差异小，冬季冰山出水高度的测量结果仍具有较好的代表性。

本文根据Landsat 8 卫星成像过程中太阳、卫星、冰山、阴影之间的几何关系，建立阴影测高模型，利用多期影像，使用阈值法和灰度梯度差异提取冰山阴影，精确计算阴影长度，测量冰山出水高度，采用交叉验证的方法评估了Landsat 8全色15 m影像提取冰山出水高度的能力。

2 研究方法

2.1 阴影测高模型

Landsat 8 卫星星下点成像，其Level-1T 数据产品默认卫星高度角为90°（来源：https：//www.usgs.gov［2019-10-31］）即卫星可以观测到冰山的全部阴影，因此可采用最简单的阴影测高模型：

式中，L为阴影长度，H为扣除海冰厚度的冰山出水高度，θ为太阳高度角。

根据上述模型，冰山出水高度测量精度取决于阴影长度的测量精度和太阳高度角计算精度。阴影长度的测量精度与影像像元大小和太阳高度角大小直接相关。单个影像像元测量的理论精度为ΔH：

ΔL为影像像元大小。图1 显示不同太阳高度角不同分辨率影像单个纯阴影像元的理论测量精度图。

图1 不同分辨率图像不同太阳高度角下单个纯阴影像元的理论测高精度Fig.1 The theoretical measurement presions of a single shade pixel in different resolutions and different solar elevation angles

2.2 冰山出水高度测量

冰山出水高度测量过程主要包括影像正影、阴影提取、阴影长度计算、成影点太阳高度角计算和冰山出水高度计算，具体实施流程如图2所示。

图2 冰山出水高度测量流程Fig.2 Processing steps of the iceberg freeboard measurement

2.2.1 正影

正影是通过图像旋转将冰山影子调整为图像的正南方向（如图3 正影影像图），旋转角度由太阳方位和影像成图坐标体系决定，其目的是实现阴影长度的自动测量。Landsat 8 卫星Level-1T 数据产品采用南极极方位立体投影（标准纬线南纬71°），其正影旋转角度由太阳方位角和冰山的经度决定，计算公式如下：

式中，γ为影像的旋转角度，α为冰山所在位置太阳方位角，β为冰山所在位置的经度。

图像旋转角度依赖于太阳的方位角和冰山的位置，不同时相不同位置的冰山其对应的旋转角度也不同。为了确保获取旋转后图像的空间几何位置信息，同时避免复杂的坐标之间的转化，正影前同时获取原始影像的（x，y）坐标影像，同步对（x，y）坐标影像进行旋转。

2.2.2 阴影提取

由于冰雪反射率高，冰山的阴影（阴影区）与冰山表面和海冰表面（非阴影区）灰度差异明显，并且阴影区边缘的灰度梯度变化明显（图3）。根据这两个特征，本文联合阈值分割法和灰度梯度差异提取阴影。通过目视解译获取阴影阈值，可初步提取阴影。由阈值分割法获取的阴影边缘可能与灰度梯度峰谷值所在位置不一致，因为边缘点是一阶导数（梯度）的局部最值所对应的像素点，因此我们沿成影方向以阈值分割法获取阴影边缘点为中心搜索最近的灰度梯度的峰值点和谷值点，替代原始的边缘点，更新阴影区。

图3 冰山示例和冰山阴影剖面示意图Fig.3 Iceberg example and schematic diagram of iceberg shadow profile

2.2.3 阴影长度自动计算

沿成影方向，遇到的第一个阴影区边缘点为冰山成影点，另一个边缘点为阴影外边缘点。获取冰山成影点和阴影外边缘点的（x，y）坐标值，依据两个点坐标值计算南极极方位立体投影下的阴影长度L。

2.2.4 太阳高度角计算

Landsat 8 影像产品仅给出每景影像中心点的太阳高度角和太阳方位角，每景影像的不同位置太阳高度角相差可能达到几度，所以精确计算每个冰山的太阳高度角是非常必要的。本文采用美国国家可再生能源实验室太阳辐射研究实验室测量和仪表团队开发的太阳位置算法SPA（NREL’s Solar Position Algorithm），根据成影点的经纬度坐标计算其太阳高度角和方位角，精度可以达到±0.0003°（Iqbal，1983；Meeus，1991）。因此，由太阳高度角造成的阴影测高误差可忽略。

2.2.5 冰山出水高度测量

冰山出水高度计算工作如下：

式中，Hiceberg为冰山出水高度，θ为太阳高度角，Hseaice为海冰出水高度，设其为常数。下文中所提取的冰山出水高度仅为扣除海冰出水高度后的高度。

2.3 测量精度评估

由于无实测数据用于结果精度验证，本文利用不同时相同名成影点进行交叉精度验证，用均方根误差（RMSE）、平均绝对误差（MAE）、平均误差（AE）和决定系数（R2）评估测量精度。

3 试验数据

试验选择了南极普利兹湾沿岸2016年8月29日、9月7日、9月16日共3期Landsat 8全色15 m Level-1T地形校正影像数据（表1）。对应3期影像内涵盖了大量的冰山，在覆盖时间内冰山位置基本未动，试验选择5 个扁平冰山为例进行3 期影像的冰山出水高度的提取和交叉验证（图4）。

表1 本研究试验所用Landsat 8影像Table 1 Landsat 8 image used in the test of this study

图4 东普利兹湾2016年9月7日Landsat 8全色影像Fig.4 Landsat 8 panchromatic image of east Prydz Bay on September 7，2016

4 结果分析

如图5所示，试验结合阈值法和灰度梯度差异提取了5 个冰山3 期影像的阴影，获取了冰山上每个成影点的阴影线（图5 蓝色和红色线），通过目视判读，它们与影像阴影非常吻合。通过严格的位置匹配，3 期数据一共获得52 个落在一个像元内的同名成影点（图5黄色点）。

图5 不同时相试验冰山的阴影长度提取结果示意图Fig.5 Schematic diagram of shadow length extraction results of test icebergs in different time

表2 和图6 详细显示了52 个同名成影点的3 期出水高度测量结果、3期相互交叉验证结果以及3期分别与3 期均值的交叉验证结果。52 个同名成影点测量冰山的出水高度在21.9—59.2 m，不同期测量结果之间交叉验证的绝对误差范围在0—6.8 m，R²均高于0.968，MAE 均小于1.5 m，RSME 均低于2.0 m。

图6 冰山出水高度交叉验证结果图Fig.6 Cross validation result of icebergs freeboard

表2 同名成影点冰山出水高度验证结果Table 2 Validation of icebergs freeboard with the matching points

3 期影像试验冰山区域太阳高度角的平均值分别为5.25°、8.15°和11.23°，单个阴影像元测高的理论精度为1.38 m、2.15 m 和2.98 m。图7 所示，与三期均值比较时，78.85%以上的冰山出水高度的绝对误差小于1.38 m，仅不到3.8%的绝对误差大于2.98 m；8月29日、9月7日、9月16日共3期出水高度估算的绝对误差在单个阴影像元测高理论精度内的比例分别是92.31%、 98.08%、96.15%，说明阴影长度测量精度在一个像元以内。

图7 冰山出水高度估算误差统计Fig.7 Error statistics of estimation of icebergs freeboard

续表

5 结 论

本文提出了一种快速提取冰山在海冰上的阴影、高精度测量成影点太阳高度角、自动计算阴影长度，并利用阴影测高模型测量冰山出水高度的方法。阴影测高的误差取决于阴影长度测量误差和太阳高度角的计算误差。此方法采用逐个冰山高精度计算太阳高度角，精度优于0.0003°，由此太阳高度角造成的误差可忽略。阴影长度测量误差与影像像元大小和太阳高度角大小密切相关，研究使用南极冬季影像中心的太阳高度角不高于11.01°的3期Landsat 8全色15 m 影像进行试验，以扁平冰山为例，验证此方法阴影长度测量误差优于一个像元，即使用此方法，当Landsat 8全色15 m影像太阳高度角低于11.31°时的冰山出水高度的测量精度优于3 m。虽然不及Wang 等（2014）利用ICESat/GLAS 数据提取的默茨冰舌崩解冰山±0.5 m的精度，但接近无人机和WorldView 系列卫星影像使用立体像对摄影测量在极地地区生成数字高程模型的垂直精度。

此冰山出水高度测量方法采用的阴影测高模型简单、人工干预少自动化程度高、对影像空间分辨率要求低、测量精度高。同样也具有一定的局限性：冰山聚集时阴影易受遮挡；冰山形态各异表面地形起伏，可能会造成的阴影的重叠和遮挡。另外，本文所提出的方法具有一定的时间局限性，仅能获取冬季冰山的出水高度。尽管冰山出水高度季节性变化率较小（Li 等，2018），但在精细测量中仍要考虑其季节差异，需结合激光高度计和雷达高度计数据对不同季节的冰山出水高度进行测量。

本试验结果表明利用中高分辨率光学卫星影像，基于阴影测高方法在估算南极冰山出水高度方面有着巨大的应用潜力。在后续研究中，将针对此方法的局限，考虑不同类型的冰山、冰山之间阴影的遮挡，同时考虑海冰出水高度以及冰山阴影的全自动提取。