基于冰川流速的喀喇昆仑地区典型冰川厚度反演与冰储量估算

王纪硕, 张齐民, 孙一丹, 张豪磊, 张子彦, 闫世勇

(1.中国矿业大学 环境与测绘学院,江苏 徐州 221116; 2.自然资源部国土环境与灾害监测重点实验室,江苏 徐州 221116; 3.山东省国土测绘院 基础测绘中心,山东 济南 250102)

在当前全球变暖的背景下[1],以青藏高原为核心的“亚洲水塔”正在发生以失衡为特征的剧烈变化,青藏高原山地冰川加速退缩,对人类社会文明和生态安全的影响不可低估[2]。青藏高原喀喇昆仑山是全球山岳冰川最发达的高大山脉,也是中纬度地区现代冰川最发育的地带,而冰川作为当地重要的淡水资源库,其厚度与储量的变化对周边区域的生态平衡起着至关重要的作用[3]。

获取冰川厚度的方法主要有实地测量和间接估计2类。采用人工钻孔法进行实地测量,会受到地理环境限制,无法进行大范围作业。根据厚度-面积比例关系式可间接估计冰川的平均厚度,但由于不同冰川各自的属性不同,使得该方法针对区域性冰川有效,对于单一冰川并不适用。20世纪80年代初期,我国首次采用探地雷达(ground penetrating radar,GPR)对天山乌鲁木齐河源1号冰川进行冰厚测量并取得一些成果,随着GPR技术的不断发展,该方法被广泛应用在我国青藏高原山地冰川的冰厚提取与冰储量估算研究中,但GPR探测冰厚的准确性受到电磁波穿透冰体能力和雷达系统图像处理能力的制约,并且冰层的细微变化会分散GPR的信号,从而降低图像清晰度,当前GPR技术多应用于小型冰川厚度提取,在青藏高原大规模冰川厚度测量中的应用很少[4]。卫星遥感技术的发展使得大规模监测冰川成为可能,在山地冰川应用较多的是美国航天飞机雷达地形测绘使命(Shuttle Radar Topography Mission,SRTM) 数字高程模型(digital elevation model,DEM)数据和Landsat卫星光学影像。文献[5]根据浅冰近似(shallow ice approximation,SIA)理论,采用SRTM DEM数据进行冰川厚度反演。SIA理论是基于冰川中流线平均基底剪切应力经验公式和冰川中流线平均坡度进行冰川厚度反演,然而该方法只能求得冰川的平均冰厚和总体积,并且由于DEM数据获取时间具有单一性,难以获取冰川坡度的年际变化,同时冰川表面运动、温度、降水等因素也都严重影响冰川厚度的变化,因而不能为冰川厚度的持续性变化提供参考。文献[6]采用Landsat卫星影像获取青藏高原东南部冰川的表面状态,根据质量平衡分布数据、冰川表面速度和冰流原理估算出冰川冰厚的空间分布;文献[7]采用Landsat 7、Sentinel-2卫星影像对喜马拉雅山脉冰川的流速进行反演,探索冰川动态变化与质量守恒的相互作用机制。但光学影像在青藏高原冰川监测中的应用极大地受限于恶劣的天气状况[8],而雷达影像具有穿透力强、覆盖范围广、不受云雨天气影响等优点,利用不同时期的雷达影像,基于冰川流速对冰川厚度进行估算,能够简便高效地通过数值模拟的方法实现冰川厚度的反演。本文基于Sentinel-1A卫星干涉宽幅(Interferometric Wide swath,IW)模式下的单视复数影像和SRTM DEM数据,采用像素跟踪算法和层流理论进行冰川厚度反演与冰储量估算,以认识喀喇昆仑地区冰川厚度空间分布及冰厚与气候之间的响应关系;将估算结果与法国南部比利牛斯天文台数据服务中心的全球冰川厚度数据产品进行对比分析,从而验证本文方法的可行性。研究结果可为研究喀喇昆仑地区冰川厚度变化提供参考。

1 研究区概况和数据来源

喀喇昆仑山是全球山地冰川最发达的高大山脉,从阿富汗最东部向东南延伸约480 km,平均海拔超过6 000 m,山区终年低温,年平均气温低于0 ℃,由于冬春两季受到西风大环流的影响,冰川成冰区域年降水量超过1 000 mm,独特的地理环境使得喀喇昆仑山有超过37%的面积被冰川覆盖。在当前气候变暖的形势下,全球绝大多数冰川出现萎缩现象,而喀喇昆仑地区冰川对气候变化响应不够显著,在某些情况下出现末端前进、面积增加的现象,被称为“喀喇昆仑异常”[9],本文选取喀喇昆仑地区却哥隆玛冰川、锡亚琴冰川和巴尔托洛冰川3条规模较大且具代表性的冰川作为研究对象。

为研究喀喇昆仑地区冰川厚度分布,采用6景Sentinel-1A卫星降轨单视复数影像,时间跨度为2019年1月—2020年1月,影像距离向和方位向空间分辨率分别约为5、20 m,观测范围覆盖研究区目标冰川,影像数据全部来自欧洲航天局(European Space Agency,ESA) (https://search.asf.alaska.edu),具体影像信息见表1所列。另外,美国地质调查局90 m空间分辨率的SRTM DEM数据主要用于研究区地形位移量补偿、地理编码及坡度计算;第6版兰多夫冰川编目 (Randolph Glacier Inventory 6.0,RGI 6.0)用于提取研究区冰川矢量边界,数据来源于全球陆地冰空间测量计划(http://www.glims.org/RGI/);冰川厚度数据产品(https://doi.org/10.6096/1007)[10]作为验证数据用于结果的精度评定。

表1 研究区影像信息

2 研究方法

为获取目标冰川流速,并在其基础上获取冰川厚度信息,本研究采用像素跟踪算法和层流理论进行冰川流速提取与冰川厚度的反演,通过与全球冰川厚度产品进行对比分析,验证该方法的反演精度,相应技术流程如图1所示。图1中,NCC表示归一化互相关(normalized cross-correlation)。

图1 冰川厚度反演技术流程

2.1 像素跟踪算法

对于流速较快的山地冰川,像素跟踪算法具有良好的适用性[11],本研究利用覆盖研究区的典型冰川多期Sentinel-1A卫星遥感影像,采用改进的像素跟踪算法进行位移量计算,该方法以NCC算法为基础,通过主、辅影像的非冰川区稳定像素的亚像素级配准结果精确解算主、辅影像之间的几何变换模型,进而估计冰川表面运动速度分布。NCC的计算准则公式为:

(1)

通过上述复杂的匹配过程,能够获取研究区总体位移量,总体位移量可分为冰川运动位移量和非冰川位移量,其计算公式为:

D=Dglacier+Dnonglacier

(2)

其中:D为研究区总体位移量;Dglacier为冰川运动位移量;Dnonglacier为非冰川位移量,其中包含合成孔径雷达(synthetic aperture radar,SAR)影像轨道不同引起的全局位移量、电离层位移量、噪声位移量和地形高低起伏产生的位移量。为去除SAR影像轨道不同引起的全局位移量,采用随机抽样一致性算法来估算变换模型,获得变换模型参数去除轨道相关位移量;由于喀喇昆仑地区处于中低纬度,而电离层引起的位移量对于高纬度地区较为明显,因此电离层影响可以忽略不计;为去除图像配准过程中产生的噪声点,在图像后处理过程中采用自适应中值滤波进行滤波处理;冰川区地形陡峭引起的位移量受到空间基线长度的影响,采用DEM数据和基线信息对产生的地形位移量进行补偿。在考虑非冰川位移量的基础上进行位移合成和地理编码,即可得到冰川区表面流速分布。

2.2 层流理论

采用层流理论进行冰川厚度的反演[12],该方法的基本原理是综合考虑冰川坡度、基底滑动、自身重力等因素对冰川厚度的影响,根据输入参数反演冰川厚度,如图2所示。

图2 层流理论示意图

图2中:H为冰厚;τ为基底剪应力;F为冰体沿斜面方向向下滑动的力,即重力沿斜面的分力;Flim为冰体静止时的最大摩擦限度;m为冰体质量。随着冰体质量增加,当斜面方向向下滑动的力超过摩擦限度,就会产生运动。

假设冰体底面积为S,冰体在斜面固定的力为τS,基底剪应力等于冰体重力沿斜面方向的分力,即

τS=mgsinα,m=ρHS

(3)

其中:g为重力加速度;α为冰川表面坡度[6];ρ为冰密度。

冰-岩界面条件变化产生的底部摩擦力不可忽略,需要增加形状因子f′,则有:

τ=f′ρgHsinα

(4)

其中,f′与山谷冰川横断面结构有关,山谷冰川f′取值范围[13]为0.5～0.9。

根据层流理论,冰川表面运动速率us计算公式为:

(5)

其中:ub为冰川基底运动速率;A为蠕变参数,与环境温度、冰体杂质等因素相关[14];n为格伦定律指数,在1.5～4.2范围内变动,与冰川所受的应力有关[15]。

结合式(4)、式(5)可以得出:

(6)

其中,k为基底滑动系数。

3 数据处理与结果分析

3.1 冰川流速提取与结果分析

利用2019年1月—2020年1月的6景降轨SAR数据,采用像素跟踪算法进行冰川流速提取,设置主影像模板窗口的大小为96,辅影像搜索窗口的大小为32,步长为2,可以在冰川区达到较高的配准精度。为削弱噪声的影响,采用自适应中值滤波进行降噪处理,设置最小滤波半径为7,最大滤波半径为15,得到目标冰川平均流速如图3所示。图3中,冰川运动方向如红色箭头所示。

图3 3条冰川流速分布图

由于冰川底部基岩本身的差异及冰川受山谷拖滞作用影响,冰川的形态呈现多样性、复杂性。从图3可以看出:却哥隆玛冰川和巴尔托洛冰川最大速率位于冰川上游(约为38 cm/d),随着冰川坡度不断减缓,冰川流速逐渐放缓,但均呈上游积累区相对于下游消融区流速较快、冰川中心相对于冰川边缘流速较快的特点,这是冰体下滑过程中与周边山体发生摩擦导致速度减缓,另外,表面冰碛物的阻热特性一定程度上减缓了冰川的消融、阻碍冰川运动;而锡亚琴冰川末端出现流速激增的情况(约为40 cm/d),相对于上游流速明显增大,期间可能发生冰川跃动,这也是喀喇昆仑地区跃动冰川常见现象[16],此外冰川流速变化符合“上游向下游逐渐减缓、中央向两侧逐渐减缓”的运动规律。

3.2 冰川厚度反演与结果分析

为分析冰川厚度分布,在获取流速的基础上,采用层流理论进行冰川厚度反演,其中蠕变系数A=2.4×10-24/(Pa3·s),滑动系数k=0.2,冰密度ρ=917 kg/m3,重力加速度[10]g=9.81 m/s2;综合考虑喀喇昆仑山脉所处气候条件及冰川结构等影响因素,将格伦定律指数n取值为2.6,形状因子f′取值为0.6。目标冰川厚度分布如图4所示。图4中,a1—a2、b1—b2、c1—c2表示冰川剖面线始末位置。

图4 3条冰川厚度反演图

由图4可知:却哥隆玛冰川的冰厚最大值出现在支流与干流交汇口,约为380 m,并且大部分区域厚度为220～250 m,占冰川面积的60.04%;锡亚琴冰川最大冰厚出现在冰川干流下游,约为380 m,大部分区域冰厚为180～220 m,占冰川面积的42.35%;而巴尔托洛冰川冰厚最大值集中在冰川干流上游,随着海拔不断升高,冰厚在210～300 m范围内不断增加,且冰川大部分区域冰厚为170～200 m,占冰川面积的57.09%。总体而言,3条冰川均呈从上游到下游逐渐减薄、支流向干流不断增厚的趋势。

为进一步分析地形因素对冰川厚度的影响,本文沿着图4中的剖面线提取3条冰川的冰厚并进行拟合,得到剖面图,如图5所示。

图5 3条冰川干流高程和冰厚变化曲线

从图5可以看出,却哥隆玛冰川冰厚存在2个峰值,分析其原因,有以下2点:① 可能是冰川支流推动作用的影响,使得上游和支流的物质在支流交汇处释放,造成支流交汇处冰体较厚;② 可能是由于冰川的搬运作用,使得冰碛物在冰川下游拐角处堆积,一定程度上阻碍了冰体的运动和消融,促进物质积累。

锡亚琴冰川由于冰下水压力的增加导致冰川末端基底滑动增强与冰形变的加速[16],致使上游积累区物质向下游推进,出现冰川下游相对上游较厚的现象。

巴尔托洛冰川为典型的树状型冰川,主干呈东西走向,冰川整体坡度较缓,厚度自东向西逐渐减薄,与高程变化一致,冰川下游海拔较低、气温较高等因素促进冰雪消融并抑制冰雪积累[17]。

3.3 冰储量估算与统计分析

为进一步估算冰川储量,将冰川厚度栅格影像重采样至50 m分辨率,与冰川边界矢量叠置分析,得到完整的冰川区,利用RGI 6.0全球冰川编目数据,通过统计冰川区的栅格像元数目,将栅格像元大小与每个像元厚度值的乘积之和作为冰川区冰储量。却哥隆玛冰川、锡亚琴冰川和巴尔托洛冰川统计面积分别为283.08、1 053.88、781.06 km2,平均厚度值分别为237.27、229.66、172.39 m,却哥隆玛冰川、锡亚琴冰川、巴尔托洛冰川的冰储量估算结果分别为67.42、242.18、149.72 km3。

为分析上述结果的可靠性,采用文献[10]中的法国南部比利牛斯天文台数据服务中心的冰川厚度产品作为验证数据,该产品经过精度检验,已被应用于全球冰川流速与厚度特征分析等研究中。本文反演结果与厚度产品对比见表2所列。

表2 3条冰川厚度、冰储量本文结果与厚度产品对比

从表2可以看出,却哥隆玛冰川、锡亚琴冰川、巴尔托洛冰川的平均厚度与厚度产品的误差分别为14.05%、13.86%、11.31%,而冰储量的误差分别为14.05%、14.30%、13.85%,表明本文冰厚反演方法具有较好的可行性,可用于辅助分析青藏高原冰川厚度时空特征差异。

4 讨 论

本文采用层流理论进行冰川厚度反演,该方法对于光学影像和雷达影像的流速结果都具有良好的适用性,但模型中不同参数的组合易导致“异参同效”现象的发生。本研究根据冰川表面流速、坡度、形状因子等模型参数的不确定性,考虑冰川不同区域的气象条件差异,结合不同山谷冰川所处气候类型、外界应力及横断面结构等属性的差异,设置合理的参数,尽可能减少“异参同效”现象,使得该方法获得较为理想的结果。

模型反演存在不确定性,其中输入数据、模型参数是造成反演结果不确定性的主要来源。

1) 根据式(6)可得,若冰川表面流速增加10倍,则会对冰厚高估超过75%[10],会导致冰厚模拟值在冰川的活跃阶段与静止阶段之间不一致,这进一步增加了厚度值的不确定性。

2) 对于山地冰川,可变参数τ和形状因子f′的取值分别与冰川谷壁横向阻力和冰川横截面的形状有关,但实测值与理论值往往有很大差异[18]。在区域尺度上,由于冰川类型、气候、地形等因素的显著差异,且统一特征参数难以确定以及缺少有关物质平衡的区域性数据,该方法在区域研究中具有较大的局限性。

3) 虽然本文中冰川流速有具体的时间信息,但由于DEM数据的匮乏,无法获取目标时间段的DEM数据,采用了2003年公开发布的SRTM DEM数据进行坡度提取,时间的不一致最终也会直接导致冰厚反演结果的误差[19]。

目前对于冰体流动的相关研究中,关于冰体在给定应力作用下的变形规律尚认识不清,且冰体的变形速率与所处温度和外界施加的压力密切相关。虽然众多研究者设计了尽可能接近冰川的冰厚实验来建立冰厚模型,但是实验条件均为理想状态,相对于自然条件存在一定差别。

影响冰川厚度变化的因素有很多,降水和温度在气象因子中占有重要地位[20]。喀喇昆仑地区降水较为充沛,但3条冰川所在的区域年均降水量有较大差异,却哥隆玛冰川年均降水量约为1 000 mm,而锡亚琴冰川、巴尔托洛冰川年均降水量分别约为600、700 mm,并且却哥隆玛冰川、锡亚琴冰川、巴尔托洛冰川年平均气温分别为-13.79、-22.93、-20.61 ℃,降水量和气温的差异会影响冰川消融与积累的转化效率,进而影响冰川的补给与冰储量的变化,增加了冰厚反演结果的不确定性。

5 结 论

本文利用6景Sentinel-1A卫星影像和90 m 空间分辨率SRTM DEM数据,应用层流理论模型对却哥隆玛冰川、锡亚琴冰川、巴尔托洛冰川3条冰川进行厚度反演及冰储量估算,分析冰厚空间分布主要特点及其原因,结论如下:

1) 基于冰川表面流速和层流理论能够反演复杂地形区山地冰川的厚度信息,通过与法国南部比利牛斯天文台数据服务中心的冰川厚度数据产品对比验证了结果的有效性。本文方法可为青藏高原地区山地冰川厚度反演与冰储量估算研究提供参考。

2) 却哥隆玛冰川、锡亚琴冰川、巴尔托洛冰川的平均厚度分别为237.27、229.66、172.39 m,同一地区不同冰川厚度略有差异,冰川的搬运作用和堆积作用与冰厚的空间分布存在着必然的联系,但冰川整体上呈从上游到下游逐渐减薄、支流向干流不断增厚的趋势。

3) 层流理论物理意义明确,相对于传统经验公式更具有稳健性,基于层流理论进行冰厚估算时,由于模型过于简化,估算结果误差偏大[15],实测资料的稀缺一定程度上限制了该方法结果的验证,如何有效利用丰富的冰川编目资料及其他辅助数据改进厚度估算模型是以后进一步研究的方向。