基于地形指数的格陵兰冰面水系建模与评价

毛玮 陆欣 陆瑶 杨康,2,3

(1 南京大学地理与海洋科学学院,江苏 南京 210023; 2 江苏省地理信息技术重点实验室,江苏 南京 210023; 3 南方海洋科学与工程广东省实验室(珠海),广东 珠海 519080)

提要 每年消融期,冰面融水通过冰面水系输送至冰盖边缘,造成格陵兰冰盖物质损失。冰面水系的动态变化决定了冰面融水输送的方式与效率,显著影响格陵兰冰盖物质平衡。作为决定冰面水系动态变化的关键因素之一,冰面地形对冰面水系动态变化的控制程度有待研究。本研究选取格陵兰冰盖西南部典型区域(～1800 km²)作为研究区,以32 m 空间分辨率ArcticDEM 作为实验数据,利用地形湿度指数(Topographic Wetness Index,TWI)和最邻近河网高差指数(Height Above the Nearest Drainage,HAND)两种地形指数模拟冰面水系分布,结合2016—2019年消融期29 景10 m 空间分辨率Sentinel-2 和30 m 空间分辨率Landsat-8 影像提取冰面水系作为验证数据,通过模糊比较揭示了地形指数建模冰面水系的精度。研究结果表明: 地形指数能够较好地建模冰面水系分布(总体精度＞72%),冰面水系分布主要受地形因素控制,将地形指数引入冰面水文的研究,有助于加强对于冰面水文过程的理解。

0 引言

格陵兰冰盖是世界第二大冰盖,近年来,格陵兰冰盖物质损失加速,从而加快了全球海平面上升,对全球气候变化产生重要影响[1-2]。冰盖表面消融是格陵兰冰盖物质平衡的重要组成部分,对格陵兰冰盖物质损失的贡献约占50.3%[2]。每年消融期,大量冰面融水通过冰面水系(Supraglacial River Network)输送至冰盖内部或者汇入冰盖边缘海洋[3]。冰面水系的动态变化决定了冰面融水输送的方式和效率,对冰盖物质平衡具有重要影响[4]。冰面水系已成为格陵兰冰盖物质平衡研究的热点[3-8]。

冰面水系的动态变化受到多种因素的影响。冰面地形是决定冰面水系动态变化的关键因素。消融开始后,融水在冰面地形控制下汇流,首先在冰雪内部孔隙中汇集,当融水含量饱和后出露冰面形成冰面饱和湿雪区(Saturated Wet Snow Zone)[9]。随着消融进一步加强,冰面饱和湿雪区进一步发育,形成冰面水系[5]。由于地形对冰面水系动态变化的重要作用,基于数字高程模型(Digital Elevation Model,DEM)研究冰面水系动态变化的研究相继展开,例如Yang 等[10]利用经典的数字水文分析方法研究了DEM 数据模拟冰面水系分布的适用性,说明了DEM 水文建模结果可以在一定程度上反映冰面水系河道位置,但是无法反映冰面饱和湿雪等融水汇集特征。冰面地形并不是控制冰面水系动态变化的唯一因素。近年来,相关研究发现冰雪层密度、近表面不透水冰层(Ice Slab)的分布、吸光性杂质(Light-Absorbing Particles)和冰川藻类微生物(Glacier Algae)等因素共同影响融水汇集过程与冰面水系分布[11-14]。然而,目前针对冰面水系动态变化及其影响因素相关的研究较少,作为决定水系动态变化的关键因素,地形对冰面水系动态变化的控制程度有待研究。

地形对水文过程的影响,可通过地形指数来反映。地形指数综合了汇水面积、坡度等因素,能够较好地反映地表产汇流过程[15]和地表湿润区空间分布特征[16],量化地形对水文过程的控制作用[17]。地形指数已被广泛应用于流域水文特征、土壤含水量时空分布、植被分布格局以及灾害评估分析等相关研究中[15,17-21]。地形湿度指数(Topographic Wetness Index,TWI)[15]和最邻近河网高差指数(Height Above the Nearest Drainage,HAND)[22]是最为常用的两种地形指数。相较于陆地而言,冰面受到土壤、植被和人类活动等因素的干扰较少,地形对冰面水文过程的影响也更加直观。利用地形指数建模冰面水系,能够反映地形为单一变量控制下的冰面水系的分布特征,对比实际冰面水系的分布特征,能够反映地形对于冰面水系的控制程度,更好地揭示冰面水系的动态变化过程。

本研究选取格陵兰冰盖西南部典型区域作为研究区,以32 m 空间分辨率ArcticDEM 作为实验数据,利用TWI 和HAND 两种地形指数模拟冰面水系分布,结合2016—2019年消融期29 景10 m 空间分辨率Sentinel-2 和30 m 空间分辨率Landsat-8影像提取冰面水系作为验证数据,定量评价地形指数建模冰面水系的精度,揭示地形对冰面水系动态变化的控制程度。

1 研究区与数据

1.1 研究区

研究区位于格陵兰冰盖西南部,面积约1800 km2,地形平坦,平均坡度 0.57°,高程范围为 1700～ 1900 m,如图1 所示。过去20年,格陵兰冰盖西南部的物质平衡线(Equilibrium Line Altitude,ELA)约为1550 m,研究区原本属于物质积累区[23]。然而,近年来,随着格陵兰冰盖表面消融显著增强,该地区的物质平衡线快速上升至1700 m以上[24]。因此,研究区已由物质积累区变为消融区与积累区的过渡区域。每年消融期,研究区内分布着冰面湖(Supraglacial Lake)、冰面水系、冰面饱和湿雪区等大量冰面水文特征,可以作为研究冰面水文过程的典型区域。

图 1 研究区 Sentinel-2 遥感影像图.影像获取时间2016.7.25; 波段组合方式 Red (band 4),Green (band 3),Blue (band 2); 图中蓝色区域为冰面湖、冰面水系、冰面饱和湿雪区等冰面水文特征Fig.1.Sentinel-2 satellite image of the study area.Image was acquired on 25 July 2016,Red (band 4),Green (band 3),Blue (band 2).The blue area in the image indicates the supraglacial lake,supraglacial river networks and saturated wet snow zone

1.2 研究数据

本研究获取了2016—2019年消融期(6—8月)云量小于10%的19 幅Sentinel-2 和10 幅Landsat-8卫星影像,平均12.5 天覆盖一次,用于反映冰面实际水系的动态变化特征(表1)。Sentinel-2 是目前可以免费获取的空间分辨率最高的卫星遥感影像,可提供10 m 空间分辨率多光谱卫星遥感影像,在全球资源环境监测领域中有广阔的应用前景[25]。Landsat-8 可提供30 m 空间分辨率多光谱卫星遥感影像,重访周期为16 天[26]。Sentinel-2 和Landsat-8数据可通过美国地质调查局(United States Geological Survey,USGS,https://earthexplorer.usgs.gov/)免费获取。综合使用Sentinel-2 和Landsat-8两种卫星影像,可以实现对格陵兰冰面水文特征的高时空分辨率监测[27]。

表1 卫星影像和ArcticDEM 数据列表Table 1.Satellite images and ArcticDEM used in this study

此外,本研究选取了2 幅空间分辨率为32 m的镶嵌ArcticDEM 数据,用于建模冰面水系。ArcticDEM 是利用WorldView-1/2/3 等高分卫星遥感数据,通过立体像对处理生成的北极高分辨率数字地形模型,涵盖北纬60°以北的所有陆地区域[28]。目前,ArcticDEM 提供条带和镶嵌两类数据产品,条带ArcticDEM 是由输入的立体像对图像的重叠区域生成的条带; 镶嵌ArcticDEM 是从条带DEM 中选取出质量最好的数据编制而成,镶嵌ArcticDEM 被混合和羽化后,减少了数据空缺区域,应用较为广泛[28-29]。本研究所用的三种数据具有不同的空间分辨率,为便于对比分析,我们将三种数据的空间分辨率统一重采样至30 m。

2 研究方法

2.1 冰面水系建模方法

地形湿度指数TWI 由Beven 等[15]提出,被广泛应用于描述局部地形对土壤含水量空间分布和地下水位变化影响的相关研究[16,19,21,30-31]。目前的研究通常将TWI 表示为:

其中,α为单位等高线长度上的汇流面积(Contributing Area,单位m2·m-1),tanβ为该点处的局部水力坡降。TWI 是对流域中各点理论土壤含水量和径流产生潜力的量化,TWI 值越大,表明该区域的土壤含水量越容易达到饱和,该区域成为地表水体分布区的可能性越大。

最邻近河网高差指数HAND由Nobre等提出[22],被广泛应用于土壤含水量时空分布、灾害分析等相关研究中[32-33]。HAND 假设在流动的河流及其周围区域的土壤含水量应处于饱和状态,随着距河流的高差增大,区域重力势能增大,土壤含水量也随着相对重力势能差变小而降低,以此来表示土壤含水量的空间分布[22]。HAND 通过DEM获取河网(Drainage Network),计算每个栅格单元距最近河网的相对高差,HAND 值越低,表明该区域越接近河网,该区域土壤含水量越高,形成水体分布区的可能性越大。

本研究利用ArcGIS 和SAGA 两款常用的GIS分析软件处理ArcticDEM 数据,首先通过多流向算法(Multiple Flow Direction,MFD)计算流向、坡度与汇流面积栅格,再设定合理汇流面积阈值生成河网,进而计算TWI 和HAND。多流向算法认为水流应向邻域中所有高程较低的方向按照局部坡度分配,能较好地模拟水流在坡面等地形上的漫散流动,相较于D8 等单流向算法更符合流向的实际情况[34]。

2.2 冰面水系提取方法

本研究采用归一化水体指数(Normalized Difference Water Index,NDWI)提取冰面水系。NDWI是提取水体信息的常用方法,能够有效地提取冰面水系,在冰面水文的研究中应用较为广泛[27,35-36]。NDWI 基本思想是利用水体的光谱特性,通过波段运算构建水体指数,增强水体光谱信息,并选取合理阈值有效地区分冰面水体和背景,实现水体信息快速提取[37]。相较于陆地水系,冰面水系具有变化速度快、生命周期短的特点[7],单期遥感影像仅能反映特定消融时刻冰面水系的分布情况,难以反映整个消融期内冰面水系的动态变化。融合不同消融时期的多时相遥感影像,能够解决这一问题,较为完整的描绘冰面水系的动态变化。

本研究收集了研究区2016—2019年消融期(6～8月)的29 幅高质量遥感影像,包括19 幅Sentinel-2影像与10 幅Landsat-8 影像(表2)。研究分别计算不同时期遥感影像的NDWI,选取合理阈值提取水体信息。阈值选取的标准为能够识别细小冰面河,同时避免大量错提水体区域的出现[35],本研究选取的NDWI 阈值范围为0.12～0.20。最后,融合多期影像提取出的水体信息,得到水体分布的频次图,作为本研究的验证数据。遥感影像的获取以及NDWI 的计算均通过Google Earth Engine平台实现。

表2 2016—2019年高质量Sentinel-2 与Landsat-8 卫星遥感影像获取日期分布表Table 2.Acquisition date distribution of high-quality Sentinel-2 and Landsat-8 satellite images during 2016-2019

2.3 精度验证方法

本研究采用模糊比较与混淆矩阵相结合的方法进行精度验证。模糊比较的思想源于Hagen 等[38]基于考虑栅格的邻域效应,考虑分类数据之间类别与位置的不确定性进行相似性比较。本研究仅考虑位置的不确定性,基于距离函数,比较模拟与验证结果。通过对比DEM 和遥感影像提取出的河网,两者在部分区域存在1～3个栅格的位置偏差,为了使比较结果更加可靠,设定距离函数,在3个栅格的范围内,认为分类数据栅格单元之间存在相似性:

其中,Ri为真值数据栅格i到其邻近验证数据同类别栅格的距离,w为栅格单元大小,fi代表真值栅格i与模拟数据的相似性,即匹配程度。当距离小于一个栅格时,匹配度为1; 大于3 个栅格时,匹配度为0,其余由公式3 计算; 匹配度越高,代表两类数据越相似,最后,汇总该类别所有真值数据的匹配度,得到整体的匹配度[16]:

式中,F为整体的匹配度,fi代表真值栅格i与模拟数据的匹配程度,n代表真值栅格数量。

本研究中,将模拟和验证数据都划分为水体和非水体两类,结合混淆矩阵的相关指标进行精度验证,包括总体精度、漏分精度、错分精度和模糊Kappa 系数。四种精度均为将模糊比较应用于传统指标后拓展得到。其中,模糊Kappa 系数计算公式为:

其中,po为总体精度,pe为预测一致性,当只考虑位置的不确定时,与Kappa 系数计算方法一致[39]。模糊 Kappa 系数能够对分类不均衡问题进行一定的矫正[40],在研究中用于最优分类阈值的选取。

为了进一步分析地形指数建模结果,将精度验证结果可视化。按照匹配度将研究区划分为匹配、漏分、错分和过渡四个区域,匹配区域表示地形指数准确识别的水体,匹配度fi=1; 漏分区域和错分区域分别表示地形指数未提取和错误提取的水体,匹配度fi=0; 过渡区域则为介于匹配与不匹配之间的区域,包括漏分区域和错分区域的过渡,匹配度fi介于0 到1 之间。

3 结果

研究通过遥感影像提取实际冰面水系,通过地形指数建模冰面水系。遥感影像提取结果展示了完整消融期内冰面水系的演变过程,从6月中旬到8月底,冰面水系经历了从发育到旺盛,最后干涸的过程(图2)。融合29景多时相Sentinel-2与Landsat-8遥感影像水体提取结果得到冰面水系分布频次图(图3a),反映了冰面水系在消融期内的动态变化。在大部分区域,冰面水系出现频次小于两次,说明冰面水系并不是稳定地分布在同一位置,而是在消融的不同阶段发生形态和位置的变化,反映了冰面水系具有高度动态变化的特性[41],在主干冰面河、冰面湖以及高海拔(＞1800 m)冰面饱和湿雪区,水体出现频次高于5次,反映了在消融期大量冰面融水汇集到这些区域[42]。其中,部分冰面饱和湿雪区和冰面湖的出现频次大于10次,最高频次达到22次。通过ArcticDEM 处理得到的TWI 与HAND 两种地形指数结果与遥感影像提取的实际冰面水系整体分布较为一致(图3b 和3c)。

研究以水体分布频次图中频次大于两次的区域作为稳定实际冰面水系,设定不同 TWI 与HAND 阈值生成模拟冰面水系,通过模糊比较对比分析遥感影像提取的实际冰面水系与地形指数建模的冰面水系,发现在TWI 与HAND 的分类阈值分别取19.5 与1.0 时,模糊Kappa 系数最大(图4),此时模拟冰面水系与实际冰面水系匹配效果最好。因此,研究以这两个阈值作为最优阈值分类TWI 与HAND 两种地形指数,得到最优的模拟冰面水系(图5)。

研究对比分析了地形指数建模的冰面水系与遥感影像提取的实际冰面水系,发现两者具有较高的相似性(图5)。尤其是对于冰面饱和湿雪区,地形指数建模结果与实际冰面饱和湿雪区在分布和形态上具有很好的一致性,说明了融水在冰雪内部孔隙中汇集,融水含量饱和后出露冰面形成冰面饱和湿雪区,这一过程主要受到地形的控制。地形指数能够以较高的总体精度(＞72%)模拟实际冰面水系分布,错分区域较少,错分精度为78%(图6)。地形指数建模结果反映了以地形作为单一控制因素下的冰面水系空间分布,模拟冰面水系与实际冰面水系良好的对应关系,揭示了冰面水系的分布主要受到冰面地形的控制。

研究对比分析了TWI 和HAND 两种地形指数建模冰面水系的结果。在同样的错分精度(78%)下,TWI 的总体精度、漏分精度和模糊Kappa 系数均高于HAND,三个指标的差异分别为8%、17%和17%(图6),尤其是对于漏分精度,TWI 和HAND 差异较为显著,说明了相较于 HAND,TWI 能够更为准确的模拟冰面水系的空间分布。因此,TWI 建模冰面水系效果优于HAND。TWI和HAND 建模冰面水系结果的差异来自于两种地形指数表达水文过程的不同思路[15,22],相较于HAND 仅考虑栅格单元距离最邻近河网的相对高差,TWI 综合考虑了汇流面积与局部水力坡降,能够更好地反映冰面水系的形成过程。

研究进一步对比分析了地形指数建模不同冰面水系要素间的精度差异。研究区两个不同高程带内分布了不同类型的冰面水系要素。1700～ 1800 m 区域主要分布冰面河与冰面湖,1800～ 1900 m 区域主要分布冰面饱和湿雪区(图1),这也反映了研究区近年来正处于从积累区向消融区过渡的特征[24]。分高程带分析能够反映地形指数建模不同冰面水系要素精度的差异(图7)。TWI与HAND 两种地形指数在1800～1900 m 区域建模冰面水系的精度均明显高于1700～1800 m。TWI和HAND 在两个区域总体精度的差异均为10%,地形指数建模冰面饱和湿雪区的效果优于建模冰面河与冰面湖的效果,说明了相较于冰面河与冰面湖,冰面饱和湿雪区的空间分布受到地形的控制程度更为显著。

4 讨论

图2 实际冰面水系提取结果(仅展示可以代表不同消融阶段的12 幅提取结果)Fig.2.Results of supraglacial river networks derived from satellite images

图3 水体分布频次图对比地形指数结果.a) 水体分布频次图; b) TWI 指数结果; c) HAND 指数结果Fig.3.Comparison of surface meltwater frequency map and topographic index results.a) surface meltwater frequency map; b)TWI; c) HAND

图4 TWI(左)与HAND(右)不同阈值下的冰面水系地形指数结果精度对比Fig.4.Accuracy of topographic index modeling supraglacial river networks with different TWI (left) and HAND (right) thresholds

图5 实际冰面水系与模拟冰面水系对比.a) 实际冰面水系; b) TWI 模拟冰面水系; c) HAND 模拟冰面水系Fig.5.Comparison of supraglacial river networks derived from satellite images and modeled by topographic index.a) derived form satellite images; b) modeled by TWI; c) modeled by HAND

图6 地形指数建模精度Fig.6.Accuracy of topographic index modeling supraglacial river networks

地形指数建模结果存在一定的错分和漏分。错分区域和漏分区域主要集中在部分冰面河、冰面饱和湿雪区及干雪区与冰面饱和湿雪区的过渡区域(图8),这主要是由于冰面水系的分布受到多种因素控制,地形并不是唯一控制因素。地形指数建模结果仅反映以地形为单一控制因素下冰面水系的空间分布,未考虑其他环境变量。地形指数建模结果仅代表该区域具备生成水系的条件[15,22],是对实际水系分布的预测,而实际冰面水系的分布 往往还受到消融强度的影响,例如Cook 等[11]发现由于藻类微生物的大量繁殖降低了格陵兰冰盖西南部的冰面反照率,从而增强冰面消融并影响冰面融水分布; Tedesco 等[43]发现分布在冰面的吸光性杂质,同样能够降低冰面反照率,改变冰面消融强度与冰面融水分布; Hofer 等[44]发现夏季云层覆盖的减少能够增强区域冰面消融,同样会对冰面融水分布产生影响。此外,Macferrin 等[13]发现冰面融水由于再冻结作用会形成低渗透性的近表面不透水冰层,能够阻止融水下渗,从而产生更多的冰面水系。通过对比现有的近表面不透水冰层分布数据,证实本文研究区内存在近表面不透水冰层分布。近表面不透水冰层一旦形成,其厚度会持续增加并逐渐接近冰盖表面,且短时间内难以消失,持续影响着冰面水系的动态变化[45]。以上这些因素,都会对冰面水系的动态变化产生影响,从而造成了地形指数建模的冰面水系存在一定错分和漏分(图8)。

ArcticDEM 获取的时间与遥感影像获取的时间存在差异,可能导致地形指数建模的冰面水系存在一定错分和漏分。然而,2016—2019年研究区冰面消融带来的高程变化(约 1.5 m)[46]小于ArcticDEM 的垂直精度(约2 m)[47],这样的变化难以被ArcticDEM 捕捉。同时,冰面水系的分布取决于相邻DEM 像素的相对高程差异,研究区高程整体降低并不会带来相对高程差异的显著变化。因此,DEM 与遥感影像获取时间的差异对地形指数建模结果的影响较小。

本文通过设定阈值从TWI和HAND中提取模拟冰面水系,TWI和HAND的阈值大小代表着冰面水系分布的可能性。TWI较小的阈值(HAND较大的阈值)能够生成规模较大、分布广泛的模拟冰面水系,对应消融旺盛期的冰面水系分布特征; 而TWI较大的阈值(HAND较小的阈值)则可以生成规模较小、分布稀疏的模拟冰面水系,对应消融初期或末期的水系分布特征。因此,不同的地形指数阈值可以模拟冰面水系在消融期不同阶段的动态变化。本文通过对比遥感影像提取的实际冰面水系与模拟冰面水系,选取出地形指数最优分割阈值,此时模拟冰面水系能够最大程度地反映实际冰面水系的分布特征。本文将TWI和HAND地形指数引入冰面水系的研究中,通过精度验证与分析,证实了地形指数能够较好的模拟 与预测冰面水系的动态变化,存在很大的潜在运用价值,例如: (1)通过设定不同地形指数的阈值,可以模拟冰面水系在不同消融阶段的动态变化,不同的地形指数分割阈值,可以得到对应不同消融阶段的模拟冰面水系,能够模拟消融期内冰面水系的动态变化,将其与水文模型相结合,能够更为准确地建模冰面水系径流量,有助于提升对于冰面水文过程的理解; (2)以地形指数建模的冰面水系为基础,分析消融强度、藻类微生物、近表面不透水冰层等非地形因素对冰面水系动态变化的控制作用,进而揭示这些因素对冰面融水输送方式与效率的影响。在未来的研究中,进一步挖掘地形指数在冰面水系研究中的价值,并推广到整个格陵兰地区,有助于进一步提升对于冰面水文过程的理解。

图7 不同高程带地形指数精度对比Fig.7.Accuracy of topographic index modeling supraglacial river networks in different elevation zones

图8 地形指数建模冰面水系结果分析.a) TWI 建模结果; b) HAND 建模结果Fig.8.Analysis of topographic index modeling supraglacial river networks result.a) modeling result of TWI; b) modeling result of HAND

5 结论

本研究以32 m 空间分辨率ArcticDEM 为研究数据,利用TWI 和HAND 两种典型地形指数模拟冰面水系分布,结合2016—2019年消融期29景10 m 空间分辨率Sentinel-2 和30 m 空间分辨率Landsat-8 影像提取冰面水系作为验证数据,通过模糊比较的方法验证了地形指数建模冰面水系的精度,揭示了地形对冰面水系动态变化的控制程度。主要结论如下。

1.冰面水系动态变化主要受地形因素控制。TWI 与HAND 分别能以80%与72%的总体精度建模冰面水系分布,地形指数建模的模拟冰面水系反映了以地形为单一变量控制下的冰面水系分布特征,模拟冰面水系与实际冰面水系存在良好的对应关系。

2.TWI 建模冰面水系的效果优于HAND。TWI和HAND 均能较好建模冰面水系分布,但TWI 建模冰面水系的总体精度、漏分精度均高于HAND,两种精度的差异分别为8%与17%,TWI 对水文过程的表达,更能够反映冰面水系的形成过程。

3.地形指数建模冰面饱和湿雪区的效果优于建模冰面河与冰面湖的效果。地形指数的精度在主要分布冰面饱和湿雪区的区域均高于主要分布冰面河、冰面湖的区域,两个区域 TWI 和HAND 总体精度的差异均为10%。