南极冰储量超过地表总冰储量的90%，其冰架物质平衡与全球气候变化、生态环境及人类社会发展等重大问题息息相关^［1］。作为衡量南极物质平衡状态与冰盖冰量损失速率的重要参数，南极冰架表面冰流速是全球气候变化建模中不可缺少的重要因子^［2，3］。同时，通过分析不同时期某一冰架表面结构特征（如冰裂隙、触地线等），可以进一步评判该冰架的稳定性及其未来的变化趋势^［4］。因此，进行长时间序列冰架表面冰流速和结构特征分析研究，可以揭示冰架的动态变化过程，进而预测海平面的变化趋势，对全球气候变化、极地冰川学的发展具有重要意义。

由于具有快速、大面积监测等优点，卫星遥感已广泛用于冰架表面冰流速和表面特征监测，且目前常用的遥感卫星包括合成孔径雷达（Synthetic Aperture Radar， SAR）和光学遥感。关于冰流速监测研究，基于SAR影像估算冰流速的方法包括差分雷达干涉（Differential Synthetic Aperture Radar， DInSAR）、多孔径干涉（Multiple Aperture Interferometric， MAI）和偏移量跟踪法等^［5］。DInSAR和MAI影像易受大气误差、去相关噪声等影响，致使测量精度降低^［6］。偏移量跟踪法虽弥补了该缺陷，但该方法需研究区域有一定的纹理特征，且精度受匹配精度影响较大^［7］。光学遥感监测基于图像特征追踪原理，通过计算不同时相遥感影像上同名点间的位移而计算相应的冰流速。相较于SAR监测，光学遥感监测成本相对较低，且其现存影像集时间跨度相对较大，更适用于长时间序列分析^［1］。国内外已有许多学者基于遥感数据开展南极冰架表面冰流速研究。Mouginot等^［8］基于Landsat-8、Sentinel-1及RADARSAT-2数据，生成了2005—2017年以年为单位的全南极冰流速图，该数据集空间分辨率为1 km，误差范围为1~17 m∙a^-1。Shen等^［9-10］基于Landsat-8数据，构建了2013—2019年以年为单位更高空间分辨率（105 m）的全南极冰流速度图。关于冰架表面特征监测，目前研究大多侧重于多源数据的融合，并在纹理提取与解译分析方法上进行改进。刘岩等^［11］提出一种基于ICESat-1数据提取冰裂隙位置及深度的方法，该方法可以探测冰裂隙峰值的应力点，进而追踪冰架最易崩解的地区。

位于南极西南方的Getz冰架（尤其是Getz F区域）是南极冰架冰量损失速率较快的地区之一。Getz冰架较阿蒙森海其他冰架对温度的升高更为敏感，而相关研究表明，冰川裂解与底部消融是南极冰架质量损失的主要部分，且这两个过程与全球温度升高密切相关^［12］。因此，开展Getz冰架表面冰流速研究显得尤为重要。Rosanova等^［13］基于Landsat系列数据测量得到南极西南极1973—1988年间平均冰流速约为200~800 m∙a^-1。Gardner等^［14］发现2008—2015年间Getz和Sulzberger冰架上游的冰川在触地线处的冰流速增加了90 m∙a^-1。

以上研究均基于遥感数据开展了南极冰架表面冰流速和结构特征分析，研究结果为揭示南极冰架稳定性及其变化趋势提供了一定参考。然而目前研究缺乏对南极冰架表面冰流速与表面高程关系探讨，而探讨冰流速与表面高程的关系，对深入探讨冰架表面冰流速的变化情况及其成因是极具参考价值的；同时，以往也未开展结合冰流速对Getz冰架稳定性进行评估的研究。Liu等^［15］提出了一种多尺度半自动化影像匹配算法重建冰架表面冰流速，该方法显著提高了计算速度，并对质量较差和虚假匹配点的修正效果较好，但其在南极Getz冰架表面冰流速重建的有效性还未曾探讨。

本文以Getz F为研究区，基于Landsat系列数据，采用多尺度图像半自动匹配算法重建了2000—2017年Getz F冰架表面冰流速，分析其时空变化特征及其随表面高程变化情况；进一步通过遥感影像增强处理及人工目视解译，提取了2000年和2017年的冰架表面结构特征（如冰裂隙），并结合冰流速变化趋势，评估Getz F冰架稳定性及其崩解趋势。

Getz冰架位于阿蒙森海中西部（115°~135° W），东邻别林斯豪森海，西邻罗斯海，面积约33 395 km²，宽度约650 km，是东南太平洋区域最大的冰架。受复杂的次冰地形及大量的冰流和出口冰川的补给，Getz冰架下存在复杂的次冰环流模式^［16］。Getz冰架前沿由Wright、Ducan等多座岛屿所固定，这些岛屿将该冰架分成了不同的部分，分别以A到F命名（图1）。本文研究区Getz F冰架位于Getz冰架最西部，是Getz冰架所有组成冰架中对温度升高最为敏感的一个冰架。

本研究基于Landsat系列数据重建Getz F冰架表面冰流速及结构特征，利用全南极陆地影像镶嵌图（Landsat Image Mosaic of Antarctica， LIMA）选取特征点进行影像配准，基于全南极冰流速度图进行精度对比验证，采用BedMachine2表面高程数据进行表面高程与冰流速的分区统计分析。

本研究采用的Landsat系列数据来源于美国地质勘探局（https：//earthexplorer.usgs.gov/），包括2000—2017年共10景Landsat-7 ETM+和Landsat-8 OLI影像（表1），数据产品等级为Level 1 GT（已经过辐射定标和几何粗校正），所有影像统一投影坐标系为“WGS 84/Antarctic Polar Stereographic”。Landsat-7 ETM+数据包含8个波段，Landsat-8 OLI数据包含11个波段，本文将选取其中空间信息最丰富的波段来提取Getz F冰架表面冰流速及结构特征（详细见3.1节）。

全南极陆地影像镶嵌图（Landsat Image Mosaic of Antarctica， LIMA）由美国地质勘探局、美国国家航空航天局和英国南极调查局基于1 073景Landsat-7 ETM+无云影像无缝拼接而成，覆盖了除82°~90° S的整片南极大陆，是现有空间分辨率最高的南极陆地影像镶嵌图，可较为清晰地展现南极冰架等的真实地表特征^［17］。本文使用其真彩色版本（空间分辨率为15 m）来选取特征点，进行影像的配准工作。

本文基于美国国家冰雪数据中心发布的全南极冰流速集（MEaSUREs Annual Antarctic Ice Velocity Maps 2005—2017， Version 1）进行精度验证^［18］。该数据集包括2005—2017以年为单位的共12景影像，空间分辨率为1 km，误差范围为1~17 m∙a^-1［8］，可于https：//nsidc.org/data/NSIDC-0720免费下载。

Bedmachine2表面高程数据来源于美国国家冰雪数据中心（https：//nsidc.org/data/nsidc-0756），由Morlighem等通过质量守恒、流线扩散等方法，并参考南极参考高程模型（Reference Elevation Model of Antarctica， REMA）制作而成^［19］，覆盖范围为70°~90° S，空间分辨率为500 m，投影坐标系为“WGS 84/Antarctic Polar Stereographic”。相较于2013年发布的Bedmap2数据，Bedmachine2表面高程数据具有更高的空间分辨率与更加丰富的细节信息。本文基于BedMachine2表面高程数据进行表面高程与冰流速的分区统计分析。

本文首先对Landsat系列影像进行预处理，再采用多尺度半自动影像匹配算法重建同名点的冰流速，最后使用普通克里金插值法提取冰流速二维图。

本文对Landsat系列影像预处理步骤包括：条带修复、几何精校正、波段选择、去云处理及影像增强。

条带修复：由于2003年Landsat-7 ETM+机载扫描行校正器（Scan Line Corrector，SLC）发生永久性故障，导致此后获取的ETM+影像（SLC-OFF影像）存在约22%数据缺失。本文基于经典的局部直方图匹配法（Local Linear Histogram Match，LLHM）进行Landsat-7 ETM+ SLC-OFF影像条带修复，主要思路如下：借助一幅覆盖同一地区的Landsat-7 ETM+影像（或Landsat其他系列影像）作为填充影像进行条带修复，若主影像数据缺失部分未被完整填充，则利用另一幅（或几幅）影像继续进行修复；若基于多幅填充影像进行修复后，主影像依旧存在数据缺失的部分，则该部分基于直接插值方法进行填充^［20］。

几何精校正：本文以LIMA全南极陆地影像镶嵌图作为基准影像，Landsat系列全色波段作为待配准波段，每幅影像选择10~20个地面控制点，基于二次多项式空间变换模型和三次卷积内插重采样法进行几何精校正^［4］。配准精度约0.5个像元（即误差约7.5 m），符合配准需求。

波段选择：本文以信息熵为主，标准差为辅对各波段所含的信息量进行评价。因具有丰富的空间信息和较高的空间分辨率（15 m），最终选取全色波段作为提取冰流速和冰架表面特征的最佳波段。

去云处理：（1）识别云层。基于归一化积雪指数和概率因子等获取潜在云层，采用形态学及光散射计算潜在云阴影层，利用归一化积雪指数识别潜在雪层；（2）云匹配。基于传感器视角、太阳方位角等预测云所在位置，同时将云层高度作为迭代条件，确保匹配步骤的合理性。

影像增强：基于空间域的Sobel滤波器来加强影像边缘信息，突出冰裂隙等要素的结构形态。

Bindschadler等^［21］提出了基于交叉互相关的匹配算法来提取两幅遥感影像上的同名点，其核心思想是在特征点的匹配窗口内遍历搜索灰度值相关性最大的点作为匹配点（图2）。该算法在冰流速重建研究中具有较高的应用价值，但仍存在一定局限性，如搜索时间慢，较难设置合适的搜索区域大小等。Liu等^［15］提出了一种多尺度半自动影像匹配算法对以上算法进行改进：

（1）影像金字塔的建立。将经过高斯低通滤波处理后的影像作为金字塔底层影像，然后将其分辨率压缩为原始影像分辨率的1/4，并基于双线性内插法使其恢复为原来的尺寸，从而建立影像金字塔。影像金字塔层数越高，分辨率越低。

（2）同名点搜索。设定配准窗口大小（一般为32×32或64×64）和搜索区域大小，进行同名点的搜索。

（3）异常点处理。基于冰架同一横截面上相近点间的运动速度和方向应相近等原理，对存在明显错误的异常点进行删除。

（4）多层次交互匹配。金字塔某层搜索结果可为其下一层搜索提供参考。由于完成删除后的结果是非均匀分布的，因此需基于反距离加权插值法估算缺失区域的数据。

从金字塔顶层开始，循环步骤2~4直到搜索至影像金字塔最底层为止。改进后的半自动影像匹配算法既能保证原始影像灰度信息的完整性，又可提高计算效率，较适用于大范围影像搜索研究。

半自动影像匹配算法计算结果为有相同特征的同名点对，而无明显特征区域则为数据空白区，需对该区域进行数据填充。本文使用克里金插值方法填补数据空白区。克里金方法基于自相关原理，具体插值公式如下：

式中：z0̑是点x0,y0处的估计值，即z0=zx0,y0，λi是权重系数，其设定需保证所求点x0,y0处的估计值z0̑与真实值z0间的差值最小，且需满足无偏估计条件Ez0̑-z0=0。

本文首先基于正态Q-Q（Quartile-Quartile）图来检测数据点分布是否满足正态分布，其次计算半变异函数来量化空间自相关性，并确定插值步长；随后根据半变异函数、采样位置及预测位置生成克里金矩阵，求解未知点的值；最后得到各未知区域预测值，并对预测值进行关联误差分析。

为分析Getz F冰架表面结构分布及其长期变化情况，选用冰流速研究时间范围内最早的2000年Landsat-7 ETM+和最新的2017年Landsat-8 OLI两幅影像，基于人工目视解译，绘制Getz F冰架表面结构（包括冰架前缘、冰裂缝、冰裂隙和纵向表面结构），各表面结构特征具体意义和识别标准如表2所示。

冰流速精度评估包括两方面：（1）基于多尺度半自动影像匹配算法重建同名点冰流速的精度评估；（2）基于普通克里金插值后的二维冰流速图的精度评估。

本文通过人工判读验证和与MEaSUREs全南极冰流速图的对比验证两种方法进行同名点冰流速精度评估。人工判读方法基于人工目视解译的127个均匀分布同名点进行验证，与MEaSUREs全南极冰流速图的对比验证法则随机选取2013—2014年和2014—2015年MEaSUREs全南极冰流速图上600个样本点进行精度评估。以本文重建同名点冰流速作为水平轴，以人工判读结果或MEaSUREs全南极冰流速同名点值作为纵轴，并拟合1∶1线进行分析（图3）。从图3（a）可以看出，与人工判读冰流速相比，本文重建冰流速均方根误差（RMSE）为12.25 m∙a^-1，标准偏差（Bias）为0.59 m∙a^-1；与MEaSUREs全南极冰流速图，本文重建冰流速均方根误差为24.47 m∙a^-1，标准偏差为14.91 m∙a^-1。两种方法均表明，本文基于多尺度半自动影像匹配算法重建的同名点冰流速结果较为可靠。

本文采用交叉验证法对普通克里金插值后的冰流速二维图精度进行评估，并拟合实际值和预测值的拟合线。若实际值与预测值具有较高的自相关性，且使用的克里金模型较为理想，则实际值和预测值的拟合线会接近于1∶1线^［4］。结果如图4所示，蓝色为实际值与预测值拟合线，可以发现其与1∶1线（黑色虚线）较为接近；同时，计算得该模型精度为±29 m，表明本文使用的普通克里金插值模型具有较高的自相关性，且得到的插值结果较为理想。

图5为基于普通克里金插值后2002—2017年间Getz F冰架冰流速二维图。可以发现，与其他年份冰流速相比，2000—2002年间冰架上游触地线处（与Berry冰川相接处）冰流速最低（500~550 m∙a^-1），较高冰流速区（600~900 m∙a^-1）面积最小，且极高冰流速区（900~950 m∙a^-1）面积也较小。同时，根据统计的2000—2017年各年Getz F冰架表面整体冰流速平均值可以发现（图6），Getz F冰架表面冰流速在2000—2002年间最低（538 m∙a^-1）。值得注意的是，2000—2002年间，Getz F冰架东部靠近Hanessian Foreland（HF）一侧尚有大片冰架存在，且其冰流速较高（600~650 m∙a^-1）。而该区域的冰架在之后的影像中消失，表明2002年后该处冰架可能发生断裂。这可能是由于该区冰架处于高冰流速区（650~950 m∙a^-1），与周围冰架冰流速差异较大，因此形成了冰裂隙，使海岸线变得崎岖，并与海洋间的摩擦增多，且同时受海洋应力和Hanessian Foreland山体的挤压而导致的。

由图5（b）可以看出，2002—2005年间Getz F冰架表面冰流速整体分布与2000—2002年［图5（a）］相似。冰架上游触地线处（与Berry冰川相接处）冰流速从500~550 m∙a^-1增加至600~650 m∙a^-1，极高冰流速区（900~950 m∙a^-1）面积较2000—2002年也有所增加，且极高冰流速区（900~950 m∙a^-1）向东部扩张明显；同时，2002—2005年间冰架整体冰流速平均值（图6）也有略微增加（从538 m∙a^-1提升至541 m∙a^-1）。

从图5（c）可以看出，2005—2007年Getz F冰架上游触地线处（与Berry冰川相接处）冰流速明显增大；极高冰流速区（900~950 m∙a^-1）与2000—2002年和2002—2005年相比，向西部扩张明显，致使2005—2007年该处海岸线向外明显突出；对比图5（b）与图5（c），中部偏东部的极高冰流速区（900~950 m∙a^-1）处的冰架在2005—2007年发生了断裂，产生了较多的冰裂隙，导致2005—2007年该处海岸线向里凹陷；区域西部整体冰流速有了一定提升，虽在靠近Berry冰川区域仍存在部分低冰流速区（0~250 m∙a^-1），但低冰流速区（0~250 m∙a^-1）面积远小于2002—2005年，且西部冰流速从上游触地线向海岸线逐渐增大，靠近Colburn Mount区域达450~500 m∙a^-1；而区域下游最东部冰流速也有了一定升高（从450~500 m∙a^-1上升至500~600 m∙a^-1）。结合冰流速平均值可以发现（图6），2005—2007年冰流速明显增大（612 m∙a^-1），是17年间冰流速最高的时间段。

由图5（d）可以看出，与2005—2007年相比，2007—2010年Getz F冰架上游触地线（与Berry冰川相接处）处冰流速明显减小，下游最东部冰流速再次回到300~350 m∙a^-1的较低水平，区域西部冰流速整体也呈下降趋势；同时，2007—2010年Getz F冰架海岸线也平滑了很多，仅剩冰架中部向海洋延伸，平均冰流速（565 m∙a^-1）较2007—2010年平均冰流速（612 m∙a^-1）大幅度下降，但依旧高于2005—2007年平均冰流速（514 m∙a^-1）。

从图5（e）中可以看出，与2007—2010年相比，2010—2013年Getz F区域极高冰流速区（900~950 m∙a^-1）面积明显缩小，仅在区域中部有小面积分布；区域中部的中高流速区向两边扩张明显，同时西部出现了一个500~550 m∙a^-1的中高冰流速中心，整体冰流速平均值从565 m∙a^-1上升至579 m∙a^-1（图6）。

图5（f）显示，2013—2014年极高冰流速区（900~950 m∙a^-1）向西部扩张明显，且海岸线有向外延伸的趋势；区域中部中高流速区继续向两边扩张，整体冰流速平均值也有了较大的提升（从579 m∙a^-1上升至611 m∙a^-1），是17年间的第二个平均冰流速峰值。

图5（g）与图5（f）总体冰流速分布情况相近。2014—2015年极高冰流速区（900~950 m∙a^-1）面积有所下降，西部冰流速整体呈下降趋势，而东部冰流速呈上升趋势，即2014—2015年冰架表面冰流速整体呈现向东部偏移的趋势。由图6（h）可知，2015—2017年海岸线中部向外延伸，补全了2014—2015年海岸线的缺口，且整体海岸线更加平滑。与2014—2015年相比，冰流速整体向东部偏移的趋势有所下降。

总的来说，从空间尺度上看，首先，Getz F冰架表面冰流速从上游触地线（与Berry冰川相接处）向下游海岸线逐渐增大，形成了一条明显的主流线，这表明Getz F冰架流量由Berry冰川补给较多，同时由于南极环极流（Antarctic circumpolar current， ACC）在135° W发展为环极深水（Circumpolar deep water， CDW），在向Getz F冰架流动的过程中发生了分流，其中一支分流发展为变性环极深水（Modified circumpolar deep water， MCDW），该分流从下方进入Getz F冰架，对冰架底部消融，使冰架变薄，加速了冰架的融解，且较薄冰架具有较低质量，使其受到的浮力小于后应力，加剧了冰流速从上游至下游的流动^［23］；其次，Getz F冰架表面冰流速从主流线向两边逐渐递减，且高冰流速区（850~950 m∙a^-1）逐渐向西部海岸线倾斜；从时间尺度看，Getz F冰架表面冰流速平均值在2000—2002年最低，随后冰流速值稳定增加，在2005—2007年达到峰值，然后又快速回落，而在2013—2014年再次达到峰值，此后冰流速虽缓慢下降，但依旧大于2000—2005年冰流速值；同时，2000—2017年间海岸线向外延伸较大。

为深入分析Getz F冰架表面冰流速的变化情况及其成因，本文统计了冰流速随表面高程变化情况。由图7可以看出，无论处于何种时间区间，Getz F冰架表面冰流速总体呈现随表面高程的增加而减少的趋势（除高程20~30 m区间）。一般情况下，温度随着海拔的升高而降低，而温度越低的区域，冰架融化速度相对较慢。除了温度的升高，海水温度升高同样会加剧冰架融化。当变性环极深水从下方进入Getz F冰架时，温暖的海水不断消融冰架底部，使冰架变薄，加速Getz F冰架底部融化。因此，Getz F冰架表面高程越高的区域，冰架表面冰流速反而越小。而其他研究在青藏高原地区得到了相反的结论，即青藏高原地区冰川表面冰流速随着表面高程的升高而增大^{［6，24-25］}。Zhang等^［24］研究指出，青藏高原冰川表面冰流速变化情况与表面高程有关，并推断可能的解释为高程较高地区接受到的太阳辐射较多，该区温度也相对较高，而温度的升高会对冰川内部和下方的热条件产生影响，如形成不稳定的冰下基，加速了冰架的融化。由于Getz F冰架所在地区表面高程（0~80 m）远低于青藏高原地区表面高程（一般为3 000~5 000 m），Getz F冰架受太阳辐射影响较小，且受变形环极深水这一特殊环流影响，更加加速了Getz F冰架低海拔区域的冰架融化。以往研究探讨冰流速与表面高程关系集中于青藏高原地区，对南极地区冰流速与表面高程关系的探讨还鲜有报告。本文为南极地区Getz F冰架表面冰流速与表面高程关系的探讨提供了一个实例，为进一步分析Getz F冰架表面冰流速变化情况及其成因具有一定参考价值。

图8为Getz F冰架在2000年和2017年的冰架表面结构。Getz F冰架纵向表面结构起源于内陆流向大陆架的支流冰川，冰裂隙在触地线及海岸线附近分布较多；由于冰流速相对较大，Getz F冰架东部主流线处冰裂隙分布也相对较多，但由于海岸线处缝合区的存在，防止了该区冰裂隙数量进一步增加；相比之下，Getz F冰架西部冰裂隙分布相对稀疏，但西部冰裂隙平均长度相对更长且深度相对更深。与2000年冰架表面结构相比，2017年Getz F冰架中西部下游冰裂缝数量明显增多，而东部中下游冰裂隙数量明显减少，且冰裂隙呈现由东部上游向东部下游移动的趋势。进一步结合Getz F冰架长时间序列冰流速分析（图5），推断以上现象可能与Getz F冰架流动方向有关。由于Getz F冰架向南方海岸线逐渐移动，且高冰流速区域（850~950 m∙a^-1）逐渐向西部海岸线倾斜，导致西部海岸冰裂隙数量逐渐增多，同时受变性环极深水影响，冰架融化速率加快，更加增加了该处的不稳定性。同时根据2000年冰架前缘存在的冰面融池及海岸凹陷推断，在2000年之前Getz F冰架前缘就已经发生过崩解。由于该处冰架较薄，随着冰裂隙及冰裂缝数量的不断增加，Getz F冰架前缘存在着较大不稳定性。

本文基于Landsat系列影像，采用多尺度半自动化影像匹配算法重建了2000—2017年Getz F冰架表面冰流速，利用普通克里金插值法提取冰流速二维图，通过遥感影像增强处理及人工目视解译，提取了2000年和2017年Getz F冰架表面结构特征，进一步探讨了Getz F冰架表面冰流速与表面高程关系，并结合冰流速对Getz F冰架稳定性进行了评估。

研究发现，Getz F冰架表面冰流速从上游触地线（与Berry冰川相接处）向下游海岸线逐渐增大，形成了一条明显的主流线，且冰流速从主流线向两边逐渐递减，高冰流速区（850~950 m∙a^-1）逐渐向西部海岸线倾斜；2000—2017年间Getz F冰架海岸线向外延伸较大，冰流速总体呈随表面高程的增加而减少的趋势。Getz F冰架表面冰流速平均值在2000—2002年最低，随后冰流速值稳定增加，在2005—2007年达到峰值，然后又快速回落，而在2013—2014年再次达到峰值，此后冰流速虽缓慢下降，但依旧大于2000—2005年冰流速值。

从2000—2017年Getz F冰架中西部下游冰裂缝数量明显增多，而东部中下游冰裂隙数量明显减少，且冰裂隙呈现由东部上游向东部下游移动的趋势。Getz F冰架流量由Berry冰川补给较多，且受变性环极深水消融影响较大；由于Getz F冰架向南方海岸线逐渐移动，且高冰流速区域（850~950 m∙a^-1）逐渐向西部海岸线倾斜，导致西部海岸冰裂隙数量逐渐增多。由于Getz F冰架前缘冰架较薄，随着冰裂隙及冰裂缝数量的不断增加，Getz F冰架前缘存在着较大不稳定性。

本文提出的冰架表面冰流速及结构特征分析方法和研究发现增加了全球变化背景下的冰架稳定性理解。后续的研究可以考虑增加冰流速与气象因子的相关性分析，进一步探讨冰架表面冰流速变化的原因。