地球物理学报  2020, Vol. 63 Issue (11): 3981-3995   PDF    
CryoSat-2测高数据与Landsat 8光学数据融合提取南极Amery冰架接地线
李斐1,2, 郭云熙1, 张宇1, 朱婷婷2, 张胜凯1     
1. 武汉大学中国南极测绘研究中心, 武汉 430079;
2. 武汉大学测绘遥感信息工程国家重点实验室, 武汉 430079
摘要:本文利用CryoSat-2测高数据与Landsat 8光学数据,针对南极Amery冰架区域开展接地线提取研究.首先通过Landsat 8光学影像三次Hermite多项式插值处理,在坡度分析的基础上利用表面曲率改进方法获取接地线特征点;同时对CryoSat-2测高数据进行坡度分析和表面曲率计算,通过沿轨梯度分析方法提取接地点;最后将Landsat 8与CryoSat-2数据获取的接地点进行最小二乘融合得到融合接地线.实验结果表明,融合结果利用高空间分辨率光学数据不仅能保证接地线提取精度,同时测高数据还能够弥补光学数据受云遮挡导致的数据空缺,保持接地线的完整性.与MOA产品比较可以看出,融合数据点与MOA接地线平均距离为367 m,标准差为601 m,所有数据点中距离小于1 km的点占总数的93.19%,与MOA产品具有较好一致.本文提出的融合算法可以实现空间连续的接地线提取结果,对后续研究南极物质平衡、冰流速计算等具有重要的意义.
关键词: Landsat 8      CryoSat-2测高      坡度分析      表面曲率      最小二乘融合     
Grounding line extraction in the Amery ice shelf, Antarctic using fusion of CryoSat-2 altimetry data and Landsat 8 optical data
LI Fei1,2, GUO YunXi1, ZHANG Yu1, ZHU TingTing2, ZHANG ShengKai1     
1. Chinese Antarctic Center of Surveying and Mapping, Wuhan University, Wuhan 430079, China;
2. State Key Laboratory of Information Engineering in Surveying, Mapping and Remote sensing, Wuhan University, Wuhan 430079, China
Abstract: In this work, fusion of CryoSat-2 altimetry data and Landsat 8 data is applied to grounding line extraction in the Amery ice shelf area, Antarctic. Firstly, based on the 3-order Hermit polynomial interpolation processing of Landsat 8 optical images, the slope analysis is carried out to obtain the revised feature points using surface curvature. For the Cryosat-2 data, grounding point is extracted by using the along track gradient analysis after slope and surface curvature calculation. Finally, the least squares based fusion method is utilized to obtain the fusion grounding line. The results show that the fusion grounding line can guarantee the accuracy compared with optical data results, and the altimetry data are not affected by the cloud that can ensure the integrity of the grounding line. Compared with the Mosaic of Antarctica (MOA) product, the average distance between the fusion data points and the MOA grounding line is 367 m, and the standard deviation is 601 m. Points less than 1 km away from all data points account for 93.19% in total, which indicates the good agreement with the MOA product. The fusion algorithm maintains the grounding line characteristics extracted by Landsat 8, and the visual performance of overall results are also improved. It is of great importance for subsequent research on the Antarctic mass balance, ice velocity calculation and so on.
Keywords: Landsat 8    CryoSat-2 altimetry    Slope analysis    Surface curvature    Feature least squares fusion    
0 引言

接地线是指冰离开冰盖并进入海洋,直接作用于海平面的位置,是海洋通过浮动冰架底部融化对内陆流动产生最大影响的地方,其位置和结构的变化能够表明冰流的不稳定性或冰运动的变化.作为估算冰流速的数值模型输入,接地线不仅可提供极地内陆冰盖状态变化的早期预警,而且也为物质平衡及其估算提供有效参考.Schoof阐述了冰盖接地线稳定性,系统分析了接地线附近的冰通量和接地线的运动状态(Schoof, 2007).Hannes Konrad最新研究表明,通过卫星跟踪冰盖的接地线,南极洲冰下水域已经在温暖的海洋水流影响下融化,仅在2010年至2016年期间,南极洲冰盖面积损失达到1, 463±791 km2(Konrad et al., 2018).

目前,Bindschadler等发表的关于南极冰盖地面和自由浮动边界的高分辨率新地图,介绍了四种接地线提取方法(Bindschadler et al., 2011).Rignot等利用差分干涉雷达测量的方法(数据源包括多种SAR数据,时间跨度为1994—2007年),能获取最准确的接地线位置(Rignot et al., 2011, 2016; Milillo et al., 2017),但是由于南极极端环境条件,特别对于Amery冰架区域,即使短时间基线的SAR影像也存在严重的失相干问题,因而能重复运用区域较少,生成的干涉影像时间分辨率较低,且空间不连续,无法进行大范围和长时间序列分析.另外流体静力学确定的静水力学线是一条不连续的线,且南极整体静力学线长度远短于接地冰边界(Bindschadler et al., 2011),使得静水力学线存在不连续的特性.

许多学者已经利用CryoSat-2测高数据进行接地线的研究,包括CryoSat-2高程数据提取方法研究(Gourmelen et al., 2014),Hogg等通过CryoSat-2高程数据计算的冰盖表面坡度突变特点绘制冰架接地线,提取结果与合成孔径雷达差分干涉测量技术(DInSAR)和ICESat重复轨道法提取结果非常一致(Dawson and Bamber, 2017),Slater等基于CryoSat-2测高生成了分辨率1 km的南极数字高程模型,对进一步利用坡度分析确定南极接地线线具有重要的参考意义(Slater et al., 2018).Fricker等结合InSAR,可见光影像和重复轨道测高绘制Amery冰架(AIS)接地区(GZ)基准地图,有助于对冰架动态应力的理解(Fricker et al., 2009).虽然Landsat 8数据在南极地表温度,冰流速等领域已开展了一定的研究,但是目前通过Landsat 8光学数据进行南极接地线提取研究极少,因此如何利用Landsat 8光学数据和Cryosat-2测高数据获取接空间覆盖连续的接地线是本文研究重点.

综合目前研究表明,接地线是可以提供内陆冰盖不稳定性的关键指标,是冰流数值模型的边界条件,确定冰盖前进以及消退的指示器,确定接地线的位置对于量化分析冰流通量至关重要(Gladstone et al., 2012; Goldberg et al., 2009; Gladstone et al., 2010; Haseloff and Sergienko, 2018; De Rydt and Gudmundsson, 2016; Kingslake et al., 2018).多源遥感数据在南极冰架变化分析中发挥重要作用(张辛, 2014),基于遥感观测数据,国际上已经发布五种全南极接地线产品,详细参数信息如表 1所示.这五种产品无法每年提供南极完整的接地线,甚至有些产品已经停止更新,这在一定程度上限制了南极物质平衡的时间序列分析,而利用CryoSat-2测高数据和Landsat 8光学数据确定接地线可以在很大程度上解决这一问题.本文重点进行CryoSat-2测高数据与Landsat 8融合提取接地线方法研究.利用Landsat 8光学数据的纹理特征与地表坡度相关性,对Landsat 8影像进行三次Hermite多项式插值处理,并进行坡度和曲率计算,针对云层覆盖的影响以及纹理特征不明显均会导致表面数据缺失问题,利用覆盖全南极冰盖边缘的CryoSat-2测高数据不受云遮挡,且基底剪切的作用,冰面高程受地形变化的影响可以通过测高数据直接观测的优势,通过CryoSat-2高程数据坡度分析方法,利用坡度标准差和均值阈值,将地面坡度突变的位置作为接地线位置,从而将高程信息转化为接地线二维信息.最终对两种数据源获取的特征数据进行最小二乘融合,得到最终的接地线(Xie et al., 2016).

表 1 5种接地线产品参数 Table 1 Parameters of five grounding line products
1 实验数据

本文选取覆盖南极Amery冰架区域的Landsat 8 OLI影像进行研究,Landsat 8卫星包含陆地成像仪(Operational Land Imager, OLI)和TIRS热红外传感器(Thermal Infrared Sensor),OLI陆地成像仪包含9个波段,除全色光波段外(空间分辨率15 m),空间分辨率为30 m,成像宽幅为185×185 km.数据信息详细统计如表 2所示,其中LC81291102017329LGN00影像云覆盖率较高,但是云层没有遮挡关键区域.

表 2 实验中采用的Landsat 8光学影像信息 Table 2 The information of Landsat 8 optical image used in the experiment

CryoSat-2卫星是欧空局发射的用于监测全球冰盖和冰川的高度变化及海冰厚度变化的雷达测高卫星.SIRAL是CryoSat-2卫星的主要有效载荷,主要任务是观测冰盖的内部结构和范围,对海冰及其地貌进行研究,SIRAL能对冰盖的不规则、陡峭的边缘进行高精度测量.SIRAL的合成孔径雷达干涉测量模式(SARIn)对于地势复杂、险峻的冰盖区域,测量精度为1~3 cm(Wang et al., 2015; 张胜凯等,2015; 李斐等,2017).本文选取CryoSat-2卫星2017年8月到12月的SARIn模式L2i级别数据进行研究.此外,考虑到MOA接地线产品是唯一完整的覆盖全南极的接地线产品(Scambos et al., 2007),MOA接地线产品利用中分辨率成像光谱仪(MODIS)遥感数据进行坡度分析获取.其中MODIS是被动式成像分光辐射计,最大空间分辨率可达250 m,生成MOA产品的数据源时间跨度分别为2003—2004年、2008—2009年(见表 1),精度为±150~±200 m,包含空间覆盖连续的线段产品,较好的保持了接地线的完整性且具有较高的空间分辨率.因此,本文将选取Amery区域MOA产品作为本文的验证数据.

2 基于Cryosat-2测高数据与Landsat 8光学数据融合的接地线提取算法 2.1 Landsat 8光学影像坡度分析提取接地线

图 1所示,本文对Landsat 8光学影像进行三次Hermite多项式插值处理,通过比较坡度、表面曲率和高斯曲率等参数,筛选出遥感影像的纹理特征(张栋,2013),利用接地线区域坡度突变的特性,建立影像灰度值和坡度之间的映射关系,通过特征筛选的方法,获取遥感影像的接地线区域.

图 1 Landsat 8影像 黄色线为MOA接地线产品,黑色线为MEaSUREs接地线产品,即InSAR_GL_Antarctica_v2. Fig. 1 Landsat 8 OLI Yellow line is MOA grounding line, and the black line is MEaSUREs grounding lines product, InSAR_GL_Antarctica_v2.

由于接地线通常位于一个几公里甚至几十公里的区域内,Landsat 8光学影像空间分辨率为30 m,为了减小计算存储量以及提高计算效率,在研究区域设置100 m分辨率的空间格网,对Landsat 8光学影像进行三次Hermite多项式插值,获取光学影像中反映空间特征的像素矩阵.三次Hermite多项式插值满足插值函数在插值节点处与被插函数有相同的函数值,而且在某些节点处或全部节点上的导数值相同,可以有效地保持影像的纹理特征.

首先通过读取原始影像获取影像素点的空间位置{(x1, y1), (x2, y2), …, (xn, yn)}及对应的灰度值{DN1, DN2, …, DNn}, 构造二维的三次Hermite插值函数DN=f(x, y).对于插值影像图中(I, J)处的像素值,首先计算该点对应原图像中的坐标(i+v, j+u),最终(I, J)处插值结果是原图(i, j)处临近16个像素点卷积之和.如果用f(i, j)表示原始影像(i, j)处像素值,插值后对应坐标值为F(i+v, j+u),则其可用式(1)表示:

(1)

式中row,col代表行列,v代表行数偏差,u代表列数偏差,S(x)插值公式定义为(2):

(2)

由此Landsat 8遥感影像实现三次Hermite多项式插值流程如下:

(1) 影像预处理,采用像素标准差对区域进行拉伸增加影像对比度;

(2) 定义空间分辨率为100 m的插值经纬度坐标网格;

(3) 获取插值影像的像素坐标和定义的影像分辨率,在遥感影像每侧加载10个像素缓冲区便于插值,避免出错;

(4) 读取像素的行和列,并进行三次Hermite多项式插值,获得最终遥感影像插值矩阵.

对于光学数据而言,基于坡度突变导致的影像灰度值差异,本文在三次Hermite插值基础上,利用坡度信息和遥感影像表面曲率信息,用于表征遥感影像的纹理特征,包括影像边界、冰水边界、岩石以及冰裂缝等.其中,坡度的表示方法包含度数法、百分比法、密位法和分数法四种.而表面曲率包含高斯曲率和平均曲率,高斯曲率能实际反应拟合曲面的弯曲程度,高斯曲率变化比较大或者比较快时,表明曲面内部变化比较大.通过计算插值矩阵的表面曲率,获取曲率变化较大位置,采用曲率标准差和均值衡量.曲面上任意一点的高斯曲率K是该点主曲率k1k2的乘积,平均曲率H是空间曲面上某一点任意两个相互垂直的正交曲率的平均值,高斯曲率和平均曲率可以利用第二基本形式和第一基本形式表示,如式(3)、(4):

(3)

(4)

式中(E, F, G)表面的第一基本形式系数,(L, M, N)表面的第二基本形式系数,其中第一基本形式系数可以具体表示为式(5):

(5)

式中rxry分别表示xy方向上的向量,其组成由式(6)表示:

(6)

式中X, Y表示经纬度坐标,Z为影像像素值,x, y为变化方向.同时,L, M, N表示的表面第二基本形式系数可以定义为式(7):

(7)

式中rxxryyrxy为对应方向上的二阶导数构成的向量,可以用式(8)表示:

(8)

通过对遥感影像插值计算获取插值矩阵,结合插值矩阵以及经纬度空间格网,计算表面坡度和表面曲率,利用计算结果标准差以及均值确定坡度具有较大变化的位置,并结合原始影像的纹理特征,获取接地线结果.

2.2 CryoSat-2数据坡度分析提取接地线

CryoSat-2测高数据是最新的能覆盖全南极冰盖边缘的高时空分辨率数据源,其边界测量精度能达到1~3 cm,能准确反应冰盖边缘变化.在测高数据接地线提取时,首先对CryoSat-2 L2i级SARIn模式的原始数据进行预处理,利用测量置信度、测量质量、重定向器和校正错误标志来剔除质量差的数据.通过平面拟合获取研究区域的高程信息,利用坡度分析方法获取坡度变化大的位置,将三维高程数据转换为与位置相关的二维信息,便于后续与Landsat 8接地线结果进行融合.我们选取一部分CryoSat-2原始高程数据进行分级色彩显示,发现其分级边界与接地线位置具有相似性,如图 2所示.

图 2 CryoSat-2 SARIn模式提取的高程数据 图中黄色线为MOA接地线产品,红色线为MEaSUREs接地线产品,即InSAR_GL_Antarctica_v2. Fig. 2 Elevation data extracted by CryoSat-2 SARIn mode data and its corresponding backscattering intensity Yellow line is the MOA grounding line product, and the red line is the MEaSUREs grounding line product, InSAR_GL_Antarctica_v2.

在测高数据接地线提取时,首先采用二次曲面拟合的方法将重轨高程数据转换为连续曲面,然后采用3×3的网格进行坡度分析,如图 3,计算公式如式(9):

(9)

图 3 DEM坡度计算3×3窗口 Fig. 3 Digital elevation model slope calculation at 3×3 window size

式中,Slopewex坐标方向上的坡度,Slopesny坐标方向上的坡度.

可以通过几种不同的算法计算Slopewe和Slopesn,本文选取的算法如式(10):

(10)

式中的Cellsize是网格DEM的间隔长度,即DEM空间分辨率,e为计算窗口对应位置的高程值.通过计算拟合表面坡度,利用式(3)、(4)计算表面高斯曲率和平均曲率,确定坡度变化较大的区域,最后对测高数据进行沿轨分析,获取CryoSat-2高程数据接地点位置.

2.3 特征最小二乘加权融合

利用前述的光学数据接地点提取方法和CryoSat-2接地线提取方法,本文采用的最小二乘融合策略如图 4所示:

图 4 利用CryoSat-2高程数据与Landsat 8光学数据融合提取南极Amery冰架接地线流程 Fig. 4 The process flowchart of extracting grounding line from Amery Ice Shelf in Antarctica by fusion of CryoSat-2 elevation data and Landsat 8 optical data

在接地线提取过程中,对于Landsat 8遥感影像受云层等天气影响造成表面数据缺失,以及纹理特征不明显的区域,可以通过测高数据进行弥补,同时利用Cryosat-2单点定位精度高的特点,并进一步结合光学数据空间连续的优势,可优化接地线的纹理信息.本文在坡度分析基础上将二者的接地线提取结果转化为经纬度位置信息,采用最小二乘方法进行融合,如式(11):

(11)

式中X表示光学数据, 测高数据经纬度坐标x, y矩阵,Y为融合后的经纬度坐标矩阵,权矩阵W=DDT,其中

(12)

D-1左乘式(12),得到式(13):

(13)

用最小二乘法估计式(13),得到参数估计量,如式(14):

(14)

但是在实际计算中,只能通过原始模型采用最小二乘法得到随机误差项的近似估计量,以此得到权矩阵估计量,即:

(15)

上述过程表明,在Landsat 8地表的纹理特征提取基础上,利用Cryosat-2数据的坡度分析结果对非接地线区域的坡度突变以及光学数据无法反映的表面特征进行限制,减小非接地线区域坡度突变干扰,可对高空间分辨率的光学数据进行补充.由于Landsat 8获取的接地线数据点数量比CryoSat-2数据获取的点要多,故需要先在Landsat 8数据点集中搜索与每个CryoSat-2数据点距离最近的点后再进行融合,可以保持光学数据的纹理特征并保证接地线的空间连续性.

3 实验结果及分析 3.1 Landsat 8影像提取结果

本文对选取的遥感影像进行拼接处理,如图 5a所示,然后对拼接影像插值计算,结果如图 5b所示,图 5中可以看出光学影像纹理信息能直观反映地表变化特征,对于相邻像素之间无明显变化特征的位置,插值结果为空值,而有纹理变化特征的区域,通过对Landsat 8插值获取地表特征变化,并将其视为数字高程信息.从图 5b中可以看出,作为计算表面坡度和曲率的原始数据,大部分位置插值为0,与光学影像反应的地表平坦特征相符,剔除插值结果为空的位置,其统计直方图如图 6a所示.

图 5 Landsat 8原始影像与原始影像像素多项式插值计算结果 (a)多张Landsat8 OLI影像拼接生成的Amery区域影像;(b)原始影像像素插值生成的计算结果. Fig. 5 Landsat 8 original image and original image pixel polynomial interpolation calculation results (a) The Amery area image map generated by Landsat8 OLI images; (b) The calculation results generated by pixel interpolation of the original image.
图 6 Landsat 8表面插值及坡度计算结果 (a) Landsat8表面插值统计直方图;(b)插值坡度计算直方图;(c)坡度计算结果. Fig. 6 The results of Landsat8 surface interpolation and slope calculation (a) Landsat8 surface interpolation statistical histogram; (b) Interpolation slope calculation histogram; (c) Slope calculation result display.

其中插值结果最大值为284.395,最小值为-30.794,标准差为106.2554,均值157.4689,且插值变化较大区域和光学数据的纹理特征具有较好的一致性.此时对插值结果进行坡度计算,如图 6b,最大值为0.9372,最小值是坡度为0的无变化特征区域.从图 6c坡度分析表明,对于光学数据有纹理特征变化的区域,计算出的坡度值,与图 5插值结果一致.

同时,为了能更加直观反应Amery冰架表面特征,本文分析了插值后的平均曲率和高斯曲率,如图 7所示.其中平均曲率最大值为0.0033,最小值为-0.0034,高斯曲率最大值为1.051×10-5,最小值为-1.435×10-5,与光学影像表面纹理特征吻合.

图 7 Landsat 8表面插值曲率 (a)平均曲率;(b)高斯曲率. Fig. 7 The interpolation curvature of Landsat 8 surface (a) Mean curvature; (b) Gaussian curvature.

图 8给出了沿轨坡度值及对应的坡度分析结果.图中曲线分别对应插值结果,表面坡度结果,沿纬度方向坡度梯度以及表面平均曲率.由图 8插值结果可以看出,插值结果为0的区域为影像边界或者水体,此时计算坡度标准差,将坡度大于标准差的数据点作为初始结果,由此计算出Landsat 8光学影像接地线,如图 9所示.图中,蓝色线为Landsat 8提取结果,黄色线为MOA产品结果,黑色线为MEaSUREs接地线产品.由于Amery冰架区域接地线比较稳定,我们对提取的5975个数据点与MOA产品对比,统计结果如图 10所示,其中平均值为369 m,标准差为594 m,包含5598个距离小于1000 m的数据点,占总点数的93.69%,即光学提取结果与已有产品结果一致.

图 8 纬度分析示例 (a) Landsat8选取的数据点位置,图中黑色点沿纬度提取结果示例(Landsat8 1897),黄色线为MOA产品结果,黑色线为InSAR产品结果,即InSAR_GL_Antarctica_v2;(b)沿纬度分析曲线:影像插值结果,插值结果坡度计算,插值结果沿纬度梯度和平均曲率. Fig. 8 Example of latitude analysis (a) The location of the data points selected from Landsat 8, the black points in the figure are an example of extraction results along the latitude, the yellow line is the MOA grounding line product, and the black line is the MEaSUREs grounding lines product, ie InSAR_GL_Antarctica_v2; (b) Tthe analysis curve along the latitude, representing the image interpolation results, the slope calculation of the interpolation results, and the gradation gradient and the average curvature of the interpolation results.
图 9 Landsat 8光学数据提取结果 蓝色线为Landsat 8提取结果,黄色线为MOA接地线产品,黑色线为MEaSUREs接地线产品,即InSAR_GL_Antarctica_v2. Fig. 9 The results of Landsat8 optical data extraction The blue line is the result of Landsat 8 extraction, the yellow line is the MOA grounding line product, and the black line is the MEaSUREs grounding lines product, ie InSAR_GL_Antarctica_v2.
图 10 Landsat 8提取结果与MOA产品距离计算结果 Fig. 10 Distance between Landsat 8 extraction results and MOA products
3.2 CryoSat-2测高数据提取结果

CryoSat-2 SARIn测高数据轨道间距较大,即使对其进行重采样缩小轨道间距,DEM精度最高只达到1 km,而且在Amery冰架区域,由于地表较强的不规则变化,使得新生成的DEM也具有不确定性.为了解决上述问题,本文在对测高数据进行坡度分析分两步:(1)利用测高数据生成DEM进行坡度分析,确定坡度变化较大区域,如图 11所示;(2)为避免轨道间距过大,导致表面拟合时空间插值带来不确定性,需要针对每条轨道的高程变化特征,设定自适应的插值参数,图中11b展示了点数大于100个的轨道数据,11(c、d)为选取其中的3430轨道和186轨道,并利用一倍标准差和均值提取出梯度变化较大的数据点.Cryosat-2坡度和表面曲率计算结果如图 12所示.图 12a展示了坡度计算的结果,图中大部分红色区域计算结果为0或者趋近0,即地表无显著变化的区域,而深色区域对应坡度变化较大的区域,具有明显的边界纹理特征.图 12b统计结果可以看出,坡度最大值为0.3863,最小值为0,标准差0.0139.从图 12a看出,此时坡度在0.1弧度附近具有显著的边界特征.图 12c的平均曲率也能反映接地线区域的边界突变特性.通过比较坡度和平均曲率,本文选取坡度大小0.1±σ(σ即标准差,此处σ=0.0139)范围内的数据点和表面曲率大于标准差的数据点,作为接地点,且坡度0.1弧度与Hogg等给出的参数一致(Hogg et al., 2018).

图 11 CryoSat-2高程数据平面拟合及其沿轨数据 (a) CryoSat-2高程数据平面拟合,图中黄色线为MOA产品结果,黑色线为InSAR产品结果,即InSAR_GL_Antarctica_v2; (b)展示区域内坐标点数大于100的各轨道高程随纬度变化情况;(c)为选取出的3430轨道;(d) 186轨道. Fig. 11 The plane fitting of CryoSat-2 elevation data and CryoSat-2 along the track data (a) The plane fitting of CryoSat-2 elevation data, the yellow line is the MOA grounding line product, and the black line is the MEaSUREs grounding lines product, ie InSAR_GL_Antarctica_v2; (b) The variation of each orbital elevation along latitude in the region with more than 100 coordinate points; (c) Shows the orbit with the track number 3430; (d) The orbit 186.
图 12 CryoSat-2坡度和曲率计算 (a) CryoSat-2测高数据坡度计算结果;(b)坡度计算结果统计直方图;(c)高程平面平均曲率计算结果;(d)平均曲率计算结果统计直方图. Fig. 12 The slope and surface curvature calculation result of CryoSat-2 (a) Slope calculation result of CryoSat-2 altimetry data; (b) Statistical histogram of slope calculation result; (c) Elevation plane mean curvature calculation result; (d) Statistical histogram of mean curvature calculation result.

图 12看出,CryoSat-2测高数据空间分辨率较低,使得平面拟合获取高程信息时会抑制表面的坡度变化特征.为了改善坡度分析和平面曲率获得的接地点结果,本文对CryoSat-2测高数据进行沿轨分析,选取每一条轨道0.1±σ(标准差)范围内的数据点,对提取的接地点进行优化,最终获取4582个数据点.与MOA接地线距离平均值1252 m,标准差1172 m,统计结果如图 13a,其中误差在1000 m以内的点占总数的54.3%,图 13b展示了CryoSat-2坡度分析确定的接地点,与MOA和InSAR产品位置一致.

图 13 CryoSat-2提取结果与MOA接地点产品距离统计结果及其展示 (a)距离统计结果;(b) CryoSat-2坡度分析确定的接地点,绿色点为CryoSat-2坡度分析确定的点,黄色线为MOA产品结果,黑色线为InSAR产品结果,即InSAR_GL_Antarctica_v2. Fig. 13 The distance between CryoSat-2 and MOA grounding points (a) The statistical result of distance between CryoSat-2 and MOA grounding points; (b) The grounding points determined by CryoSat-2 slope analysis, the green point is determined by the CryoSat-2 slope analysis, the yellow line is the MOA grounding line product, and the black line is the MEaSUREs grounding lines product, ie InSAR_GL_Antarctica_v2.
3.3 融合结果

虽然坡度分析和平面曲率确定了Landsat 8光学影像和CryoSat-2测高数据位于接地线附近的数据点,但从图 6图 7可以看出,部分区域存在数据缺,且利用CryoSat-2测高确定的接地点数量相对有限.为了综合利用光学影像的纹理特征以及测高数据点提供的高程信息,采用最小二乘方法对两种提取结果进行融合,并将数据点转化为线,生成连续的接地线,如图 14所示.

图 14 利用CryoSat-2高程数据与Landsat8光学数据最小二乘融合提取南极Amery冰架接地线结果 红色线为融合结果,黄色线为MOA产品结果,黑色线为InSAR产品结果,即InSAR_GL_Antarctica_v2. Fig. 14 The fusion result of CryoSat-2 elevation data and Landsat8 optical date for the Antarctic Amery ice shelf grounding line by the least-squares constraints method The yellow line is the MOA grounding line product, and the black line is the MEaSUREs grounding lines product, ie InSAR_GL_Antarctica_v2.

图 15展示了两种数据在两个局部区域的提取结果,其中绿色点表示CryoSat-2测高数据提取结果,蓝色数据点表示Landsat 8插值坡度分析和曲率分析选取的结果,从图中可以看出融合结果能较好反映两种数据源与MOA接地线的位置关系.但是从图 16展示出的几个局部区域可以看出,上述区域Landsat 8影像存在纹理特征缺失的现象,导致得到的数据点较少,且插值结果无法反映影像表面特征,而通过高程坡度和平面曲率分析获取对应区域的特征点边界,能够对光学影像进行弥补,保证接地线的空间连续性.

图 15 融合结果细节示例一、二 图中绿色点为CryoSat-2数据获取的坡度分析点,蓝色点表示利用标准差选取的Landsat8光学影像坡度数据,红色线为融合结果,黄色线为MOA产品结果,黑色线为InSAR产品结果,即InSAR_GL_Antarctica_v2. Fig. 15 The details of fusion results The green points are calculter by the slope analysis using the CryoSat-2 data, the blue points are the slope data points of Landsat 8 optical image slope data selected by the standard deviation, and the red line is the fusion results, the yellow line is the MOA grounding line product, and the black line is the MEaSUREs grounding lines product, ie InSAR_GL_Antarctica_v2.
图 16 融合结果细节示例 图中绿色点为CryoSat-2测高数据坡度和平面曲率分析获取的数据点,蓝色点表示利用标准差选取的Landsat8光学影像坡度数据,黄色线为MOA接地线产品,黑色线为MEaSUREs接地线产品,红色线为融合结果. Fig. 16 The details illustration of fusion results The green points are the slope changing points from the CryoSat-2 data calculation, and the blue points are the slope data points of Landsat8 optical image selected by the standard deviation, the yellow line is the MOA grounding line product, the black line is the MEaSUREs grounding line product, and the red line is the fusion results.

此外,本文对融合结果与MOA产品结果进行定量分析,计算各融合数据点与MOA接地线距离,从表 3可以看出,融合结果与MOA距离平均值367 m,标准差601 m,所有数据点中距离小于1 km的点为5569个,占总数93.19%,绝大部分区域能够与MOA产品一致.定量分析结果如表 3所示,由图 10图 13可以看出融合结果不仅保持了Landsat 8提取的接地线的纹理特性,同时对CryoSat-2数据获取的接地线具有很大的提升,整体结果有一定的改善.

表 3 各提取结果与MOA产品距离 Table 3 Distance between extracted results and MOA products
4 结论

本文利用Landsat 8光学影像进行三次Hermite多项式插值计算,获取能表征地形突变的影像纹理特征,通过计算插值影像的坡度和表面曲率,获取冰盖地表起伏变化特征,利用坡度标准差和表面曲率标准差将计算结果反应的共同边界特征作为光学影像提取的接地线;同时,对CryoSat-2高程数据拟合生成DEM,进行坡度分析以及表面曲率计算,利用沿轨梯度分析方法,减小平面拟合空间插值带来的坡度分析不确定性.最后通过最小二乘算法将两个不同数据源获取的数据点边界进行融合得到融合接地线.

与MOA接地线产品比较发现,采用融合方法得到的Amery冰架接地线结果的平均距离误差为367 m,与MOA产品具有良好的一致性.同时融合产品能够有效的保持光学影像的纹理特征和反映测高数据的高程变化信息,测高数据同时能弥补光学影像数据缺失以及表面特征不明显的区域.与已有的产品(时间为2003—2004, 2008—2009)比较发现,在东南极地区,除了Frost, Denman和Recovery冰川平均退缩在19~45 m·a-1之间,Mertz, Budd and Shirase冰川和Slessor冰流平均增长速度在14~48 m·a-1之间之外,其他区域接地线移动变化主要集中在零附近.

Amery冰架接地线近10年内保持稳定,年平均变化约为±25 m,与Hannes Konrad的研究(Net retreat of Antarctic glacier grounding lines)具有一致性,也证明了本文算法的有效性.此外,通过融合获得Amery区域完整的接地线,能弥补DInSAR失相干导致的数据缺失,扩展接地线产品的时间覆盖范围,并对后续的物质平衡研究具有一定意义.

References
Bindschadler R, Choi H, Wichlacz A, et al. 2011. Getting around Antarctica:New high-resolution mappings of the grounded and freely-floating boundaries of the Antarctic ice sheet created for the International Polar Year. The Cryosphere, 5(3): 569-588. DOI:10.5194/tc-5-569-2011
Dawson G J, Bamber J L. 2017. Antarctic grounding line mapping from CryoSat-2 radar altimetry. Geophysical Research Letters, 44(23): 11886-11893. DOI:10.1002/2017GL075589
De Rydt J, Gudmundsson G H. 2016. Coupled ice shelf-ocean modeling and complex grounding line retreat from a seabed ridge. Journal of Geophysical Research:Earth Surface, 121(5): 865-880. DOI:10.1002/2015JF003791
Fricker H A, Coleman R, Padman L, et al. 2009. Mapping the grounding zone of the Amery Ice Shelf, East Antarctica using InSAR, MODIS and ICESat. Antarctic Science, 21(5): 515-532. DOI:10.1017/S095410200999023X
Gladstone R M, Payne A J, Cornford S L. 2010. Parameterising the grounding line in flow-line ice sheet models. The Cryosphere, 4: 605-619. DOI:10.5194/tc-4-605-2010
Gladstone R M, Payne A J, Cornford S L. 2012. Resolution requirements for grounding-line modelling:Sensitivity to basal drag and ice-shelf buttressing. Annals of Glaciology, 53(60): 97-105. DOI:10.3189/2012AoG60A148
Goldberg D, Holland D M, Schoof C. 2009. Grounding line movement and ice shelf buttressing in marine ice sheets. Journal of Geophysical Research:Earth Surface, 114(F4): F04026. DOI:10.1029/2008JF001227
Gourmelen N, Escorihuela M J, Shepherd A, et al. 2014. Fine Ice Sheet margins topography from swath processing of CryoSat SARIn mode data.//AGU Fall Meeting. San Francisco, California: AGU.
Haseloff M, Sergienko O V. 2018. The effect of buttressing on grounding line dynamics. Journal of Glaciology, 64(245): 417-431. DOI:10.1017/jog.2018.30
Hogg A E, Shepherd A, Gilbert L, et al. 2018. Mapping ice sheet grounding lines with CryoSat-2. Advances in Space Research, 62(6): 1191-1202. DOI:10.1016/j.asr.2017.03.008
Kingslake J, Scherer R P, Albrecht T, et al. 2018. Extensive retreat and re-advance of the West Antarctic Ice Sheet during the Holocene. Nature, 558(7710): 430-434. DOI:10.1038/s41586-018-0208-x
Konrad H, Shepherd A, Gilbert L, et al. 2018. Net retreat of Antarctic glacier grounding lines. Nature Geoscience, 11(4): 258-262.
Li F, Song G Y, Yang Y D, et al. 2017. Establishment of high precision DEM in Antarctic Dome A area with taking the waveform retracking, slope correction and the data fusion into account. Acta Geodaetica et Cartographica Sinica (in Chinese), 46(4): 403-410.
Milillo P, Rignot E, Mouginot J, et al. 2017. On the short-term grounding zone dynamics of Pine Island Glacier, West Antarctica, observed with COSMO-SkyMed interferometric data. Geophysical Research Letters, 44(20): 10436-10444. DOI:10.1002/2017GL074320
Rignot E, Mouginot J, Scheuchl B. 2011. Antarctic grounding line mapping from differential satellite radar interferometry. Geophysical Research Letters, 38(10): L10504. DOI:10.1029/2011GL047109
Rignot E, Mouginot J, Scheuchl B. 2016. MEaSUREs Antarctic grounding line from differential satellite radar interferometry, version 2. Boulder: NASA National Snow and Ice Data Center Distributed Active Archive Center.
Scambos T A, Haran T M, Fahnestock M A, et al. 2007. MODIS-based Mosaic of Antarctica (MOA) data sets:Continent-wide surface morphology and snow grain size. Remote Sensing of Environment, 111(2-3): 242-257. DOI:10.1016/j.rse.2006.12.020
Schoof C. 2007. Ice sheet grounding line dynamics:Steady states, stability, and hysteresis. Journal of Geophysical Research:Earth Surface, 112(F3): F03S28. DOI:10.1029/2006JF000664
Slater T, Shepherd A, McMillan M, et al. 2018. A new digital elevation model of Antarctica derived from CryoSat-2 altimetry. The Cryosphere, 12(4): 1551-1562. DOI:10.5194/tc-12-1551-2018
Wang F, Bamber J L, Cheng X. 2015. Accuracy and performance of CryoSat-2 SARIn mode data over Antarctica. IEEE Geoscience and Remote Sensing Letters, 12(7): 1516-1520. DOI:10.1109/LGRS.2015.2411434
Xie H, Chen L, Liu S, et al. 2016. A least-squares adjusted grounding line for the Amery ice shelf using ICESat and Landsat 8 OLI data. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 9(11): 5113-5122. DOI:10.1109/JSTARS.2016.2614758
Zhang D. 2013. Research on ice sheet feature extraction of Lambert Glacier drainage basin, Antarctica based on ICESat and ice radar data[Ph. D. thesis] (in Chinese). Nanjing: Nanjing University.
Zhang S K, Xiao F, Li F, et al. 2015. DEM development and precision analysis in two local areas of Antarctica, using CryoSat-2 altimetry data. Geomatics and Information Science of Wuhan University (in Chinese), 40(11): 1434-1439.
Zhang X. 2014. The monitoring of Antarctic snow and ice changes from the multiple-sources remote sensing data. Acta Geodaetica et Cartographica Sinica (in Chinese), 43(4): 437.
李斐, 宋国云, 杨元德, 等. 2017. 南极Dome A地区高精度DEM的建立——顾及波形重定、坡度改正及数据融合. 测绘学报, 46(4): 403-410.
张栋. 2013.基于ICESat和冰雷达数据的南极Lambert冰川流域冰盖特征提取研究[博士论文].南京: 南京大学.
张胜凯, 肖峰, 李斐, 等. 2015. 基于CryoSat-2测高数据的南极局部地区DEM的建立与精度评定. 武汉大学学报·信息科学版, 40(11): 1434-1439.
张辛. 2014. 多源遥感数据在南极冰雪变化监测中的应用研究. 测绘学报, 43(4): 437.