时间序列冰面湖水深反演及储水量变化监测

杨晓彤,谷祥辉,赵彬如,张峰,邢喆

(国家海洋信息中心,天津 300171)

0 引言

格陵兰冰盖物质损失是导致全球海平面上升加速的重要因素[1],其中冰盖表面融水的流失和冰川崩解是格陵兰冰盖物质损失的主要途径[2-3]。作为格陵兰冰盖表面融水的重要存储单元[4],冰面湖会引起水力压裂事件,使融水突然流入冰盖底部,导致格陵兰冰盖物质损失[5-7]。对消融期内冰面湖的形成、演变和消失过程开展变化监测,尤其是掌握冰面湖能够存储多少融水、冰面湖储水量的动态变化情况,为掌握格陵兰冰盖表面融水的存储、传输和释放机制发挥着至关重要的作用[8]。

格陵兰冰面湖实地调查难度大。鉴于冰面湖分布在格陵兰冰盖的大部分地区,利用无人机等新型测绘手段会极大耗费人力和物力成本[9],卫星监测则成为获取冰面湖信息最快速、最有效的手段[10-14]。基于卫星影像计算冰面湖储水量,首先需要获取冰面湖的水深信息。但目前星载格陵兰冰面湖水深反演技术主要受到两方面制约。一是缺少实测水深数据导致常用的各类回归分析模型无法得到应用。即使能够获取少量水深实测数据,由于实测点过于稀疏,也只能用于少数冰面湖水深反演精度的验证[15-16]。如何对星载冰面湖水深反演结果开展精度验证仍有待研究。二是目前大多极地冰面湖水深反演所用数据源的空间分辨率无法满足日益增加的精细化冰面湖变化监测需求。当前基于卫星遥感的冰面湖水深反演的常用数据源包括中低分辨率的MODIS[17-18]、ASTER[19-20]、Landsat8[21-23]和Sentinel2[24-25]等,高分辨率影像由于受可见光波段较少和影像幅宽较窄的限制,一些物理参数无法从影像上直接获取,因此应用较少。尽管已有研究显示应用数字高程模型(digital elevation model,DEM)进行冰面湖等水文要素的识别需要设定合理的阈值[26],但是在无法获取实测水深的情况下,根据格陵兰冰面湖的形成原理,DEM对于获取冰面湖水深信息仍将发挥积极作用,已成为格陵兰冰面湖水深反演的重要支撑数据[27]。对于多时相冰面湖水深反演,难点在于获取多时相DEM数据与各时相影像相对应,而Sundal等[28]指出冰面湖具有年际重现的特征,且在不同年份同一位置出现的冰面湖的形态基本保持不变。Leeson等[29]的研究表明虽然格陵兰冰盖表面高程处于动态变化中,但冰面地形洼地的形态在较长时期内基本保持稳定,揭示了将单一时相DEM数据产品和多时相卫星遥感影像用于冰面湖水深反演的可能性。

综上所述,本研究基于覆盖2017年整个消融期的9幅4波段WorldView多光谱影像和Arctic digital elevation model(ArcticDEM)数据,以格陵兰东北部一个研究区为例,研究了一套兼顾工作效率和提取精度的时间序列冰面湖水深反演方法,完成了覆盖一个完整消融期的冰面湖储水量计算,并开展了研究区内冰面湖面积、最大水深和储水量的变化分析,以达到为格陵兰冰盖表面融水的存储、输送和释放机制研究提供有效技术支撑和科学参考的目的。

1 研究区与数据

1.1 研究区概况

研究区位于格陵兰岛东北部,经纬度范围为23°4′56″W～23°35′8″W,78°59′20″N～79°4′8″N,总面积为163.2 km2。该区域位于一个快速冰流区,在2017年6月至8月内共有7个冰面湖在此形成、发育并消失,选择该研究区可以验证本文提出的方法对于不同发展阶段冰面湖储水量计算的适用性,并开展一个消融期内时间序列冰面湖变化监测。因此选择该区域开展冰面湖储水量计算方法研究,具有一定的典型性。

图1为2017年8月1日的WorldView2遥感影像图,显示了在盛夏时期该区域内的7个冰面湖的形态特征及分布情况,为了便于后续分析,本文将研究区内的冰面湖进行了编号。

图1 研究区内冰面湖分布图

1.2 数据及数据预处理

研究数据包含9幅4波段WorldView2/3卫星影像,覆盖整个研究区的ArcticDEM数据,以及研究区附近实测站点气温数据。

9幅WorldView2/3影像的获取日期标识分别为20170614、20170624、20170706、20170715、20170720、20170801、20170810、20170819和20170826,基本覆盖了格陵兰冰盖的一个消融期,同时相邻日期影像的时间间隔不超过10 d,适用于监测冰面湖储水量在整个消融期内的变化发展趋势。影像预处理步骤主要包括辐射校正、FLAASH大气校正和基于RPC模型的正射校正、影像之间的相对校正等,并将全部影像设置为WGS84 UTM 27N投影。

ArcticDEM是利用光学立体影像制作的高分辨率、高质量的北极数字高程模型。本研究选取的ArcticDEM产品能够覆盖研究区内全部冰面湖,产品生产所用星源为WorldView2,所用影像时相为20170504,早于本研究卫星影像获取时间,且在格陵兰消融期之前。产品精度为-0.041 m,能够为开展冰面湖储水量计算及变化监测提供数据支撑。

气温数据为美国国家海洋和大气管理局发布的格陵兰岛DANMARKSHAVN气象观测站点记录的2017年实测气温数据。该站点位于76.767°N,18.667°S,为距离研究区最近的气象站点。数据包含观测日期、观测小时点、气温、海平面气压、风向、风速等信息。由于冰盖表面消融是一个累积过程,仅收集与影像日期相对应的单日气温数据不能代表区域内整体的气温变化趋势,因此本研究从2017年的观测数据中提取6月1日至8月31日的每小时气温数据,计算得到日平均数据,用于近似代表研究区的日平均气温变化走势。

2 研究方法

本研究技术方案见图2,研究步骤主要包括冰面湖边界高精度提取、冰面湖水深反演及精度验证、冰面湖储水量计算和时间序列冰面湖变化分析等。

图2 时间序列冰面湖储水量计算方法技术路线图

2.1 冰面湖边界提取

冰面湖边界提取采用基于形态学滤波的时间序列冰面湖信息提取方法,实现方法及提取结果详见前期研究成果[30]。应用该方法提取的冰面湖精度平均优于2个像元,时间序列冰面湖提取总精度达88.37%,满足后续水深反演和储水量计算精度要求。

2.2 冰面湖水深反演

1)基于物理模型的水深计算法。根据比尔-朗伯定律,Philpot[31]提出了反演冰面湖水深的单波段法,计算如式(1)所示。

z=g-1[ln(Ad-R∞)-ln(Rw-R∞)]

(1)

式中:z为冰面湖水深;Ad为冰面湖底部的反射率;R∞为深水处的反射率;Rw为观测到的冰面湖水体在某一波段的表观反射率;g为双向光谱衰减值[32]。一旦参数(Ad,g,R∞)被确定,利用该公式和影像表观反射率即可计算水深。Ad的数值可通过读取冰面湖边缘像元的反射率值获得,如果同景影像上有海洋覆盖,则可通过读取海水的反射率作为R∞的数值,g则根据经验或参考他人研究选取。

针对本研究,由于研究区远离格陵兰冰盖边缘,且WorldView影像幅宽较窄,无法直接获得R∞。g的数值会随着同类型遥感影像、不同特征研究区、不同气象和水文条件等发生变化,尤其对于WorldView影像,能够参考获取g值的研究较少。因此,本研究无法直接应用物理模型计算水深。

2)基于数字高程模型的形态拟合法。格陵兰冰面湖是由冰面融水在地形相对地洼的区域汇集而形成的,而DEM数据可以反映包括地形洼地在内的冰面高程信息,因此学者们提出了用DEM拟合洼地三维形态,计算冰面湖储水量的方法。假设DEM数据获取时,冰面地形洼地均未被融水填充,则DEM数据能够较好地刻画洼地三维形态,且已有研究证明了单时相DEM数据对于反演不同时相的冰面水深的可靠性。因此,本研究利用ArcticDEM和冰面湖边界信息,采用形态拟合方法反演9期影像上所有冰面湖的水深。

3)水深反演精度交叉验证。由于无法获取实测水深数据,因此将精度验证常用的交叉比对方法引入格陵兰冰面湖水深反演精度评价。由于基于物理模型的水深计算方法的可靠性在南北极冰面湖水深反演研究中得到了大量验证[33-35],本研究对形态拟合法得到的水深反演结果与基于物理模型的水深计算结果进行比对。若二者一致性较高,则认为基于形态拟合法的水深反演结果具有较高的精度,交叉比对方法如下。

选择盛夏时期、冰面湖发育较为完整的时相为20170801的水深形态拟合结果,从当期影像和形态拟合水深反演结果上提取全部像元的波段表观反射率(Rw) 与水深(z)的对应关系,即点对(Rw,z),分别从蓝光、绿光、红光波段上提取3组点对。

分别从3组点对中随机抽取50%的样本点作为输入值,代入水深反演的物理模型中(式(1)),基于Levenberg Marquardt算法[36-37]进行非线性最小二乘法拟合,得到物理模型3个系数Ad、R∞和g的拟合结果,通过拟合结果决定系数R2,从蓝、绿、红波段中确定最优的水深反演物理模型计算波段。

用最优波段中剩余50%的点对作为验证点,将拟合系数代入物理模型计算每个验证点的水深值,与形态拟合得到的水深结果进行比对,以验证两种方法的结果一致性。将基于物理模型的水深计算结果看做实测水深,精度评价指标有决定系数R2、均方根误差(RMSE)和平均相对误差(MRE)。

基于拟合系数计算其他8期影像上的冰面湖水深,与形态拟合法的水深反演结果进行比对,计算R2、RMSE和MRE,验证二者的结果一致性。

2.3 冰面湖储水量计算

以冰面湖边界提取结果为范围约束,采用像元迭代积分方法计算每个冰面湖在不同时相上的体积,即可得到冰面湖体积,即储水量信息。冰面湖储水量V的计算如式(2)所示。

(2)

式中:s为冰面湖像元数量;hs为冰面湖第s个像元的水深深度;As为像元面积。

3 结果及分析

3.1 非线性拟合结果分析

从8月1日影像和基于形态拟合法的水深反演结果上共提取了1 080 524个点对,从中随机抽取50%用于物理模型参数的非线性拟合。蓝、绿、红波段的非线性拟合结果见表1,绿光波段的拟合决定系数R2最大为 0.837。图3为蓝、绿、红波段反射率与形态拟合法水深反演结果的拟合结果图,图中散点颜色由蓝到红表示点的分布密度由低到高,可以看出绿光波段的拟合效果最好(图3(b))。

表1 基于日期为20170801的影像的非线性拟合结果

图3 可见光波段与水深反演结果的非线性拟合散点图

综上,确定绿光波段为物理模型水深反演最优波段。基于绿光波段的水深反演物理模型的拟合系数分别为:Ad=0.791,R∞=0.066,g=0.104。将拟合系数代入式(1),Rw取绿光波段的反射率值,即可得到基于物理模型的水深计算结果。

3.2 水深反演精度评价及分析

利用剩余50%的点对作为验证点,开展反演结果一致性交叉比对。经计算,两个水深反演结果的R2达到0.803,RMSE为0.853 m,MRE为18.09%。

采用相同的方法,开展其他时相影像上基于物理模型的水深反演结果与基于形态拟合的水深反演结果的交叉比对。经计算,R2最小为0.581,最大为0.936,平均为0.728,具有较强的统计相关性;RMSE最大为1.743 m,平均优于0.690 m;MRE最大为23.80%,平均优于20.43%(表2)。水深反演的总体精度能够达到80%,冰面湖水深反演的精度可以满足后续冰面湖储水量计算的要求。

表2 水深反演精度交叉验证结果

3.3 储水量计算结果及变化分析

研究区内7个冰面湖的时间序列储水量计算结果见表3,研究区6月1日至8月31的日平均气温变化情况见图4,冰面湖储水量、冰面湖面积和最大水深随时间的动态变化情况见图5。

表3 7个冰面湖的时间序列储水量计算结果(×104) m3

研究区内冰面湖的总储水量呈先升后降的变化规律,与气温变化走势一致(图4)。研究区内冰面湖的总储水量于8月1日达到最大,为1 753.79×104m3,7个冰面湖的总储水量累计可达5 331.35×104m3。

图4 研究区日平均气温变化图

虽然相隔较近但7个冰面湖的消长变化特征却有所不同,储水能力也有较大差异。4个冰面湖(4号至7号)的变化趋势总体与气温变化同步,即在8月1日气温最高时,储水量、面积和最大水深均达最大(图5(d)～图5(g))。在相同气温条件下,1号冰面湖直到8月19日一直处于扩张状态,而后迅速结冰,符合格陵兰冰面湖快速演进的变化特点[38-41]。2号冰面湖和3号冰面湖则在气温仍处于上升趋势的过程中,即出现面积和储水量减小的现象。

储水能力方面,储水能力最强的1号冰面湖最大储水量达到993.39×104m3。1号冰面湖最大水深可达14.11 m,最大面积可达138.46×104m2(图5(a)),均是7个冰面湖中最大的,且该湖水深和面积处于持续增大的状态,表明该冰面湖所在地形洼地的深度和范围均是研究区内最大的,且该湖所在洼地的地势较低,持续有冰面融水向此处汇集。储水能力最弱的5号冰面湖,最大储水量仅为1.54×104m3,该湖的最大面积和最大水深均是7个冰面湖中最小的(图5(e)),即其所在地形洼地的深度和范围均从根本上限制了融水在此汇集。

注:图中对冰面湖的最大水深和最大面积进行了数值标注。图5 研究区内7个冰面湖时间序列消长变化图

虽然2号冰面湖的面积最大可达105.5×104m2,小于7号冰面湖的最大面积124.08×104m2,大于4号冰面湖的最大面积86.94×104m2,但2号冰面湖的储水能力却显著偏弱,最大仅为16.71×104m3。7号冰面湖储水量最大可达475.07×104m3,4号冰面湖为397.27×104m3。原因有二:一是该湖的最大水深仅0.84 m,即所在洼地的深度较浅,在面积条件相当的情况下,洼地深度限制了融水存储;二是该湖在7月15日即出现缩减趋势,远远早于其他冰面湖。冰面湖快速干涸往往作为辅助判断锅穴存在的依据。据此,推测2号冰面湖湖底可能有锅穴存在,导致冰面融水向冰盖内部输送。7月15日至8月1日2号冰面湖向冰盖内部输送的水量约1.79×104m3。

类似地,监测到3号冰面湖的储水量在8月1日即降低。经影像判读,该湖通过冰面径流与外界连通,推测发生了融水输送现象。经计算,7月20日至8月1日3号冰面湖在冰盖表面输送的融水约1.4×104m3。此外,2号和3号冰面湖发生储水量降低时,其最大水深却并未下降,亦表明二者储水量的下降不是由于洼地地势、深度,或气温变化等原因导致的,而是发生了冰盖表面融水输送现象。

对于4～7号冰面湖,其储水量与最大面积和最大水深的消长趋势基本保持一致(图5(d)至图5(g)),亦未发生储水量的异常降低,表明这些冰面湖不存在导致融水“流失”的要素。尽管如此,5号和6号冰面湖由于面积和深度均较小,其储水能力也较小,弱于发生融水输送现象的2号冰面湖和3号冰面湖。

综上,冰面湖的储水能力由其所在洼地的地形条件、冰面湖底部是否存在锅穴、冰面湖的连通性等条件共同决定,其中地形条件指洼地所处的地势条件、洼地范围大小以及洼地的深度等,是决定冰面湖储水量的决定性因素。

4 结束语

本研究开展了9期4波段WorldView2/3影像上的冰面湖水深反演方法研究,通过精度交叉比对验证了方法的适用性和有效性,并以此计算每个冰面湖的储水量。通过监测时间序列影像上的冰面湖面积、最大水深和储水量的变化,分析了研究区内7个冰面湖在一个完整的消融期内的时空变化特征和冰面湖的形成、发育和消亡特点,得到以下结论。

1)对于缺少实测水深数据的格陵兰区域,可以选择合适的单一时相Arctic DEM数据,开展基于形态拟合法的时间序列冰面湖水深反演,水深反演精度满足要求。通过交叉比对验证了基于形态拟合法和基于物理模型法的水深反演结果的一致性,二者决定系数R2平均达到0.728,具有较强的统计相关性。在9期影像上,基于形态拟合法的水深反演RMSE平均优于20.43%,MRE平均优于0.69 m,证明了方法的可行性和有效性。

2)针对高空间分辨率的4波段WorldView影像,应优先选择绿光波段开展基于物理模型的水深反演。针对本研究,基于单期水深信息形态拟合结果与波段反射率拟合得到的物理模型拟合系数(Ad=0.791,R∞=0.066,g=0.104)可以直接应用于其他时相影像的冰面湖水深反演。经交叉验证,反演结果与形态拟合法所得结果一致性较好。

3)本研究从冰面湖储水量、面积和最大水深3方面开展了冰面湖变化监测。在整个消融期内,研究区内冰面湖的总储水量最大可达1 753.79×104m3,总储水量累计可达5 331.35×104m3。虽然7个冰面湖在空间分布上相隔较近,但各自消长变化趋势和储水能力却存在较大差异。冰面湖的储水能力由其所在洼地的地形条件(地势、范围、深度)、冰面湖底部是否存在锅穴以及冰面湖的连通性共同决定,其中地形条件是决定性因素。

4)通过监测冰面湖储水量和最大水深的变化规律,能够辅助判断湖底锅穴的存在。针对本研究,2号和3号冰面湖的面积和储水量峰值的出现日期均早于研究区内其他冰面湖,推测原因为通过湖底锅穴和冰面径流,发生了融水输送现象。经估算,7月15日至8月1日,2号冰面湖通过湖底锅穴向冰盖内部输送的水量累计约1.79×104m3,3号冰面湖通过冰面河在冰盖表面输送的融水约1.4×104m3。

后续研究拟将本文研究方法应用至更大空间范围,分析冰盖表面融水输送现象的影响因素和时空规律、冰面融水的年际变化特征等,以更好揭示格陵兰冰盖表面融水的存储、输送和释放机制。