南极冰盖表面消融影响冰盖物质平衡和全球能量平衡[1]。一方面,表面消融形成的融水以三种方式影响南极冰盖的物质平衡[2,3],融水汇集后从冰盖边缘流出[4];渗入裂隙后加剧裂隙的扩张,破坏冰架稳定性甚至导致冰架的大范围的快速崩解;渗入冰盖内部,促进冰盖内部通道的发育,抵达冰盖底部后减弱冰盖底部运动阻力,导致冰盖底部滑动[5,6]。另一方面,南极地表97.8%为雪覆盖[7],干雪具有高反照率,而蓝冰、湿雪、融水反照率较低[8],其吸收的太阳辐射约为干雪吸收的三倍[9,10,11],形成“反照率低-吸收太阳辐射多-加剧融化”的正反馈机制[12,13],因此,消融范围及持续时间的变化将影响全球的辐射收支[5,14]。

尽管机制有所差异,但蓝冰、湿雪和融水都是冰雪消融的产物[3,15,16],本文统称它们为消融区。在大多数年份南极冰盖表面消融区面积占比都小于20%,远低于占比常年超过40%的格陵兰冰盖,但南极冰盖表面融化的面积是格陵兰的两倍多[9]。目前,格陵兰冰盖消融研究体系较完善[17,18],南极冰盖消融研究相对缺乏[1,19]。南极冰盖消融研究更多为大尺度冰盖整体的研究,区域性精细研究较少且多集中于南极半岛[20]。从消融的观测对象上看,冰面融池及河流水系的研究相对深入[21,22],而湿雪和蓝冰监测有待更多的开展。

遥感技术是冰盖表面冰雪消融大范围观测的最佳手段。湿雪的范围及持续时间的变化可以通过微波辐射计、散射计、SAR以及光学数据监测。前三者为微波数据不受云的影响而光学数据易受到云的影响。微波辐射计数据对于消融的监测主要依赖于干湿雪亮温的显著差异[23,24],散射计数据对于消融的监测主要依赖于干、湿雪后向散射系数的显著差异[25],两种数据主要监测大尺度区域冰雪表面消融起止和持续时间、范围随时间的变化,无法满足对消融特征空间细节方面的描述[26,27,28]。SAR数据对于湿雪的监测也主要依赖于干、湿雪后向散射系数的显著差异[29,30],其高空间分辨率可描述消融空间细节特征,可用于精确监测湿雪范围[31,32],但重访周期较长,适用于精确的年际变化分析,无法满足季节性分析[16,26]。光学数据对于湿雪、融水及蓝冰的监测主要依赖于水、冰和雪光谱特征的差异。例如利用融水和冰雪在蓝光和红光波段的光谱差异使用冰雪归一化水体指数开展冰面融水提取[4,22];基于光学数据对南北极水系水文特征的分析[4,22,33,34];基于Landsat-7影像数据蓝冰与湿雪裸岩的光谱反射率差异对整个南极蓝冰区的提取分析[35]。光学数据用于湿雪的分析却鲜少开展。综上所述,冰雪消融的监测难以同时满足高时、空间分辨率,高时间分辨率的被动微波数据难以满足消融细节特征的描述,高空间分辨率的SAR及光学数据难以实现季节性变化的相关分析,而两类数据很难融合。因此,目前亟需高时、空分辨率冰盖消融的季节性分析。

本文选用Landsat-8数据开展东南极达尔克冰川消融区监测,其多光谱信息可用于湿雪、蓝冰和融水的提取,其较高的空间分辨率可满足对冰盖表面消融细节特征的监测。此外,达尔克冰川刚好处于Landsat-8不同轨道影像的重叠区,重访周期为2—7天,可满足融化季节内多时间序列消融观测。本文提出并验证了基于Landsat-8影像的MNDWIice消融区自动监测方法,以此为基础,对2016—2017年东南极达尔克冰川消融区的时空分布及季节性特征进行了初步分析。

1、研究区与数据

1.1研究区概况

达尔克冰川位于东南极拉斯曼丘陵东侧(69°26′S,76°27′E),邻接埃默里冰架(图1a),距离中山站直线距离不超过3km,为距离中山站最近的冰川。达尔克冰川流域起源于南极冰盖内部(图1b),流入接地线东部的冰源为自上游600多千米的高原,而流入接地线西部的冰源为自上游仅不到10km的冰盖边缘(图1c),总面积超过1900km[2][36]。达尔克冰川南部接地线距离冰川前缘约15km,漂浮部分即冰架面积约为20km[2](图1),是典型的冰盖溢出型冰川。冰川前端入海口宽度为3km,随着冰川的不断运动而崩解[37]。

1.2研究数据

本研究采用2016年9月1日—2017年4月30日覆盖达尔克冰川无云或少云的Landsat-8影像,此时间段内覆盖达尔克冰川的数据有41景,目视解译可用于研究的少云或无云影像共有18景影像(表1)。数据来自USGS网站:https://glovis.usgs.gov。影像由于采集时间的不同,太阳天顶角不同,日地距离也可能存在差异,为了降低这些差异的影响,使得不同时相影像数据具有可比性,将Landsat-8Level-1数据的像元灰度值(DN值)转换为大气表观反射率(TOA反射率)。反射率计算采用Kang等[38]的方法,计算公式来自Landsat手册(http://Landsat.usgs.gov/),计算所需参数从每个Landsat-8数据集的头文件中读取。此外,考虑到高纬度地区的同一幅影像,太阳高度角从中间到边缘的差异较大,计算采用逐像元计算的太阳高度角替代原始使用的单一中心点太阳高度角[38]。下文中所做的计算及制图全部基于反射率数据。

2、研究方法

2.1消融区的光谱特征

利用多光谱光学数据识别目标地物的思路是利用目标地物与背景地物在不同波长上的波谱反射特性差异,选择最佳的波段进行比值运算方式,通过选择合适的阈值最大程度地区分目标地物与背景地物。本文利用Landsat-8影像识别的消融区,包含裸露的蓝冰(无论是否发生融化)、湿雪和融水。

图1达尔克冰川概况.

蓝冰生成过程中,冰内的气泡被挤出,使得冰晶变大,蓝冰呈现蓝色。而湿雪和融水由于水的存在也显示蓝色,三者在遥感真彩色波段上均显示蓝色,区别于干雪(图2)。Landsat-8影像不同地类的光谱信息如图3所示,蓝冰、湿雪和融水在近红外波段(波段5)表现出强烈的吸收而干雪表现出强烈的反射,在该波段消融区(蓝冰、湿雪和融水)与干雪、阴影和云差异最为明显,在蓝光波段(波段2)各地类差别最小。基于此,本文选择利用近红外波段(波段5)和蓝光波段(波段2)进行比值运算,对消融区进行提取。

表1本研究阈值试验所用Landsat-8影像

*其中有云指原始数据云量>20%,少云指0%<云量≤20%

图2蓝冰、湿雪、融水、干雪、阴影样本对比图.底图皆为2016年12月21日Landsat-8反射率数据,显示方式皆为波段(5,3,2)对应(红绿蓝)标准差拉伸增强图像

2.2改进的冰雪归一化消融指数

本文采用的改进的冰雪归一化消融指数MNDWIice是基于归一化水体指数NDWI及冰雪归一化水体指数NDWIice提出的。归一化水体指数NDWI是对绿光波段(G)与近红外波段(NIR)进行归一化差值处理,以凸显陆地表面的开阔水域[39],其表达式为:

冰雪归一化水体指数(NDWIice)是对蓝光波段(B)与红光波段(R)进行归一化差值处理,以凸显冰盖表面融水[22],其表达式为:

NDWIice最早提出是用于对格陵兰冰盖表面融水的提取[22],目前已用在南极表面融水深度的确定[4]。本文基于消融区在Landsat-8数据的近红外波段(波段5)和蓝光波段(波段2)的光谱特征的强响应,提出改进的冰雪归一化消融指数MNDWIice为:

式中,Bblue和Bnir分别为Landsat-8蓝光波段(波段2)和近红外波段(波段5)的反射率。

基于2016年12月21日获取的两期影像,选取融水、蓝冰、湿雪、干雪、云和阴影6类地物的纯像元样本,计算各地类对应NDWI、NDWIice和MNDWIice的指数值,如图3(b)所示。试验显示MNDWIice对于消融区的区分效果最好,其他两种指数除对融水的区分度较明显外,其余几类地物的指数值相差较小,很难进行区分。

图3各地物类别反射率及指数特征.

2.3基于F值法的阈值获取和消融区提取精度评估

为选取最优阈值,引入机器学习评估方法——精度F值法(F-Measure),用F值评估分类结果精度。F值由准确率(Precision,P)和召回率(Recall,R)两个度量值计算得到,F值越高表明提取效果越好。准确率即为影像分类中的用户精度,召回率即为影像分类中的生产者精度(制图精度)。目标像元分为四类:TP,TN,FP,FN。TP是指消融区被正确划分为消融区的像元数,TN是指干雪被正确划分为干雪的像元数,FP是指干雪被错误划分为消融区的像元数,FN是指消融区被误判为干雪的像元数[40,41]。准确率P用于评估被正确提取的消融像元占所有被提取为消融像元的比例,召回率R用于评估被正确提取的消融像元占应该被提取为消融像元的比例。F值法是准确率和召回率的加权调和平均。

根据召回率与准确率对于研究的重要性给β赋值,当两者同样重要时,β值取1;当准确率更重要,β值取小于1的值;召回率更重要时,β值取大于1的值[40]。本研究需同时考虑消融区被准确分类为消融区的概率和分类结果影像中消融区像元对应原始影像为消融像元的概率,因此定义β值取1。

根据MNDWIice阈值选取和分类结果精度评估需要,我们分别选取了试验样本和验证样本(如图1所示),为保证阈值选取结果的准确性,试验样本选择位于消融区和干雪区交界处的容易发生误分的区域,图1(d)中7个白框即为试验所用的试验样本区(包含干雪区及消融区),试验样本的面积占研究区面积的9.94%,各期影像的干雪及消融区试验样本即基于各期影像反射率数据对白框区域采用目视解译获取;验证样本选择的是30个空间随机分布的纯消融像元或者纯非消融区像元,每个验证样本像元大小为5×5,像元需满足纯样本条件,即像元区域为纯消融区像元或者纯非消融区像元,由于消融初期消融区面积小,20160907期影像难以选择一定数量的消融区纯样本,考虑到验证样本的一致性,故从20160914期影像选取30个纯消融样本和30个纯非消融样本,后续影像由于阴影等其他影响不满足纯像元条件的,选择其他纯像元做补充,保证验证样本为纯消融或纯非消融样本。试验样本与验证样本存在一定比例的重合,但各期影像中重合面积不超过试验样本面积的0.95%,不超过验证样本面积的22.34%。

我们旨在利用单一的MNDWIice阈值进行不同时相的消融区的分类,单一阈值TMNDWIice计算方法如下:

式中,TMNDWIice为18期影像利用F值法选取到的每一期最优阈值,Fmax(i)为最优阈值对应的F值。每一期的最优阈值TMNDWIice获取方法:根据经验设定MNDWIice阈值范围为[0.01,0.2],以0.001为步长,利用试验样本数据计算每个阈值分类结果对应的F值,获取MNDWIice阈值-F值曲线,选择F值最大值Fmax(i)点处的阈值为最优阈值(图4)。

图420161221期影像阈值选取实验

本文基于Landsat-8level1数据计算反射率,利用反射率数据计算得到各时相MNDWIice数据,利用18期影像的MNDWIice结合样本区目视解译结果确定各时相最优阈值及单一阈值,最后利用目视解译的验证样本结果,对消融区自动提取结果进行精度评估,具体实现流程详见图5。

3、结果

3.1基于MNDWIice消融区提取效果评估

基于试验样本的阈值试验,18期影像消融区提取MNDWIice最优阈值如图6a所示,各时相MNDWIice阈值分布范围为0.1—0.17,根据每一期的最优阈值和对应的F值,依据公式7计算了18期影像消融区提取的单一阈值TMNDWIice为0.136。采用各期影像最优阈值对其余18期影像试验样本的提取结果进行F值法精度评估(图6a),Fmax(i)值整体优于65.9%,准确率R优于56.4%,召回率P优于65.3%;采用单一阈值0.136对其余18期影像试验样本的提取结果进行F值法精度评估(图6b),Fmax(i)值整体优于63.4%,准确率R优于47.2%,召回率P优于59.4%。

图5基于Landsat-8数据的MNDWIice消融区自动监测方法技术流程图

图618期影像基于试验样本的消融区提取结果.

a)各时相最优阈值及对应F值;b)各时相单一阈值提取的召回率、准确率及F值

利用17期影像(18期影像除去20160907期影像)的验证样本,分别对采用单一阈值及最优阈值进行消融区提取的17期影像提取结果进行F值法精度评估(表2)。采用单一阈值0.136的提取精度为,Fmax(i)值整体优于67.7%,17期影像的Fmax(i)均值高于81.5%,准确率R优于56.8%,召回率P优于66.4%。采用最优阈值的提取精度为,Fmax(i)值整体优于71.3%,17期影像的Fmax(i)均值高于83.5%,准确率R优于56.8%,召回率P优于85.9%。

表217期影像验证样本分类精度评估

图7分别显示了利用单一阈值提取消融区分类结果与试验样本目视解译结果对比,F值最低的两期影像和分类结果F值最高的一期影像及对应的分类结果图。分类精度较低的分别为20160907及20170104两期影像,其中20160907期影像精度最低为19.0%,消融初期太阳高度角低,受阴影影响大,因此存在阴影误分为消融区现象;融化像元较难目视勾选,试验样本在目视解译时也存在错误分类的现象。20170104期影像的分类结果的F值为47.2%,据参与该时间段南极科考的人员回顾,一月份刚刚下完新雪,对提取结果可能存在影响,且同样受到云及阴影的影响。但两幅影像的分类结果,目视效果良好(图7)。20170311期影像提取结果最好,分类结果对应的F值达到了78.6%。

图7单一阈值提取结果对比:提取精度最低两期影像与提取精度最高一期影像对比图.

3.2达尔克冰川消融区季节性空间分布特征

利用MNDWIice采用单一阈值0.136自动提取的2016年9月—2017年3月达尔克冰川的消融区时空分布结果如图8、9所示。从时间分布来看,消融区面积呈现出较为明显的季节特征(图8),消融区从9月份或之前开始出现,面积整体呈上升趋势,在南极的正夏1月份达到面积最大值,之后消融区面积开始下降,与中山站自动气象站温度数据变化趋势基本一致。阈值提取的消融区面积,12月份的4幅影像单一阈值提取的面积明显大于最优阈值提取的面积,其余12幅影像单一阈值提取的面积则小于或约等于最优阈值提取的面积。

消融区MNDWIice均值随时间变化总体上呈现与消融面积相似的变化趋势,消融前期均值呈上升趋势,在2016年12月28日达到最大,最大为0.119,随后均值下降,呈现出明显的季节性变化特征。这表明消融前期消融面积较小,且MNDWIice均值较小,随着消融的发生,消融面积与MNDWIice均值增大,到1月份前后均达到最大值,随后消融面积与MNDWIice均值降低。消融发生的时间不晚于Landasat-8数据观测到的最早时间(9月7日),且此时气象观测数据显示温度低于零摄氏度。

从空间分布来看(图9),消融主要发生在冰盖前缘,由海岸线向内扩展的十几公里内。消融最早发生在坡度大的东部接地线处。消融持续时间最长的区域主要分布在东部接地线处以及冰架前端,由表面纹理可以看出,这两处区域都覆盖了大量表面裂隙,表面粗糙度较其他处更大。

4、结论与讨论

本文根据Landsat-8多光谱影像显示出的干雪与消融区(蓝冰、湿雪和融水)的光谱特征差异,构建了改进的冰雪归一化消融指数MNDWIice,利用18期影像获取了用于提取南极季节性消融区的MNDWIice阈值(0.136)。经试验证明,摒弃太阳高度角度很低、阴影影响很大的一期数据,采用单一阈值0.136自动提取消融区的平均精度为81.5%,最低精度为67.7%。另外,本研究验证样本在随机布点的基础上,考虑样本均匀分布的同时,特意保留了干雪验证样本阴影区等易发生误分区域的覆盖,这会一定程度上造成精度的低估。因为阈值选取所用数据覆盖了9月到次年4月份,太阳高度角从7.49°到35.62°,兼顾到消融季节光学数据的复杂性;试验数据18景,样本量具备统计学意义;试验样本及验证样本均匀分布,满足试验的随机性。因此,本文认为在精度要求不高于80%时,基于MNDWIice采用单一阈值0.136可用于整个南极冰盖季节性消融区提取。消融区的错分和漏分主要来自于阴影及云的影响,存在一定程度的干雪阴影区被误提为消融区像元的现象,后续工作将增加云及阴影指数,最大程度降低云及阴影的影响,从而进一步提高提取精度。

图8消融区的季节变化特征.

图9基于18期影像达尔克冰川消融季节变化.图例由左至右表示观测到的发生融化的时间越来越早

消融区提取结果显示,达尔克冰川的消融呈现出明显的季节分布,随着气温升高,消融区面积呈现上升趋势,在1月份达到最大值,说明太阳辐射很可能是达尔克冰川季节性表面消融的主导因素。早期的消融区主要分布在达尔克冰川的东部坡度较大的接地线区,很可能是坡度大加剧了下降风,而下降风的侵蚀会促使蓝冰的形成并会带走蓝冰表面的降雪。蓝冰的反照率低,会加剧表面融化。消融持续时间最长区域为裂隙富集区,裂隙增加了表面粗糙度,可能是消融发生的有利条件。由于研究仅局限于达尔克冰川区域且仅做了单年的分析,要想获取整个南极的消融特征,需要更大范围和更长时间序列的观测。

参考文献:

[3]杨康.格陵兰冰盖表面消融研究进展[J].冰川冻土,2013,35(1):101-109.

璩榆桐,程晓,刘岩.基于Landsat-8的东南极达尔克冰川季节性表面消融信息提取[J].极地研究,2020,32(02):164-176.

基金:国家重点研发计划项目(2016YFA0600103);国家自然科学基.