表碛覆盖型冰川的提取方法及变化

薛 娇，姚晓军，张 聪，周苏刚，褚馨德

（西北师范大学地理与环境科学学院，甘肃 兰州 730070）

0 引言

冰川是指由降雪和其他固态降水积累，在重力作用下形成的具有流动特征的冰体［1］。作为冰冻圈的重要组成部分，冰川被视作气候变化的理想指示器［2］，对水资源平衡、生态系统安全和社会经济发展具有重要影响［3-5］。在气候变暖背景下，全球冰川普遍处于退缩状态，近60年全球参照冰川平均累积物质亏损达26.84 m w.e.［6］。然而，不同类型冰川对气候变化的响应程度差异明显［7］。

依据我国冰川所处区域的气候类型，可将其分为大陆型冰川和海洋型冰川（次分为季风海洋型冰川和西风海洋型冰川）［8］，前者广泛分布在昆仑山、天山、喀喇昆仑山和羌塘高原等地区，主要依靠低温严寒气候提供的冷储条件而存在，具有补给少、消融弱和运动速度慢等特点，占中国冰川面积的77.8%；后者主要分布在横断山、喜马拉雅山东段以及念青唐古拉山中东段，依靠较为充沛的降水条件而存在，具有温度高、消融强烈和运动速度快等特点，占中国冰川面积的22.2%［7，9］。由于受到发育环境和物理机制差异的影响，不同类型的冰川对气候变化的响应过程和程度存在较大的差异［9］，且部分冰川消融区覆盖着由冰川侵蚀、搬运和堆积作用而形成的表碛。与裸冰区相比，表碛覆盖区受到下垫面的反射率、颗粒物的大小及厚度等物理性质的影响［10］，呈现出与裸冰不同的响应特征［11-12］。因此研究不同类型表碛覆盖型冰川变化对正确认识气候变化规律及其影响程度意义重大。

冰川边界自动/半自动提取方法主要包括波段比值法［13］、雪盖指数法［14］、主成分分析法［15］、面向对象分析法［16］、合成孔径雷达干涉测量技术［17］和机器学习方法［18］等，但受到冰川区表碛和非冰川区地物光谱特征较为相似的限制，上述方法对表碛覆盖型冰川的提取效果较差，难以准确区分表碛的边界［19］。为提高冰川表碛区遥感识别准确性，Holobâcă等［20］使用基于合成孔径雷达相干图像和雪盖指数法的DebCovG-carto工具箱提取大高加索山脉两条冰川的表碛，该方法有效缩短了处理时间并降低了绘制错误；Xie等［18］分别采用雪盖指数法、大津阈值算法和机器学习算法提取巴基斯坦罕萨山谷冰川上的表碛，与人工目视解译相比，随机森林模型识别表碛准确度最高，为（95.5±0.9）%；宋波等［21］利用表碛与岩石的温差并结合地物光谱特征和DEM数据，实现了海洋型冰川表碛区的提取，但受影像质量影响，地表温度数据识别部分表碛区的有效性较低，导致分类中出现错分和漏分现象；上官冬辉［22］基于NDVI、波段比值、DEM数据和地表温度数据，利用知识库分类法对天山科其卡尔冰川的表碛覆盖区进行提取，精度达到90%；Shukla等［23］基于短波红外和热红外波段的特征差异提出归一化差异表碛指数（normalized difference debris index，NDDI）提取表碛覆盖型冰川的边界，该指数在ASTER影像上的提取效果较好，但受限于阈值影响，大范围冰川提取效果并不理想。为此，本文基于Landsat TM/ETM+/OLI遥感影像及ASTER DEM数据，利用提出的TDSI（temperature NDDI slope ice）方法提取天山、喀喇昆仑山和念青唐古拉山地区9条巨型表碛覆盖型冰川边界，分析2011—2020年冰川及其表碛面积的变化趋势，探讨大陆型冰川和海洋型冰川表碛覆盖面积变化的差异性及原因，从而为进一步探究表碛覆盖型冰川对气候变化的响应机制提供科学依据。

1 研究区概况

据中国第二次冰川编目［24］，中国共发育冰川48 571条，面积约5.18×104km2，其中表碛覆盖型冰川的数量和面积分别为1 723条和12 974.70 km2，占冰川总面积的25.05%［10］。为探究大陆型和海洋型表碛覆盖冰川变化的原因及差异，本文选取9条面积＞100 km2的表碛覆盖型冰川［24］作为研究对象（图1）。其中6条为大陆型冰川，分别为：①音苏盖提冰川（36°40′N，76°60′E），位于喀喇昆仑山脉乔戈里峰北坡，面积359.05 km2，长度约43 km，是中国目前已知最长最大的复式山谷冰川，中值海拔为5 475 m；②木斯塔冰川（35°51′N，76°16′E），位于音苏盖提冰川南侧，面积182.54 km2，海拔介于4 217～7 291 m，中值海拔为5 475 m；③托木尔冰川（41°55′N，79°59′E），位于天山西部托木尔峰西南坡，面积358.25 km2，长约39 km，中值海拔为4 744 m；④乌库尔冰川（42°16′N，80°20′E），位于天山西部汗腾格里峰木扎尔特河源区，多条不同规模的支冰川汇入，面积167.84 km2，长度超过30 km，中值海拔为4 462 m；⑤木扎尔特冰川（42°18′N，80°55′E），位于乌库尔冰川东侧，面积165.83 km2，长约37 km，中值海拔为4 655 m，是特克斯河地表水的主要源头；⑥穹特连冰川（41°57′N，80°90′E），位于托木尔冰川东侧，与其紧密相连，面积约140.95 km2，中值海拔为4 854 m。3条海洋型冰川分别为：⑦雅弄冰川（29°20′N，96°39′E），位于西藏东南部岗日嘎布山的然乌错流域，面积171.09 km2，长约33 km，海拔介于3 930～6 456 m，中值海拔为5 300 m，冰川末端与岗日嘎布湖相连；⑧夏曲冰川（30°30′N，94°28′E），位于念青唐古拉山东段南坡西侧，面积133.58 km2，长约22 km，中值海拔为4 969 m；⑨那龙冰川（30°33′N，94°50′E），位于夏曲冰川东侧，面积和长度分别为103.53 km2和20 km，中值海拔为5 291 m。

图1 典型表碛覆盖型冰川的分布和遥感影像Fig.1 Distribution and remote sensing images of typical debris-covered glaciers

2 数据与方法

2.1 数据源和预处理

本文所使用的遥感影像为2011—2020年Landsat TM/ETM+/OLI数据。为降低云雪影响，影像选取以冰川消融期无云无雪或少云少雪为主，从美国地 质 调 查 局（United States Geological Survey，USGS）网站（https：//earthexplorer.usgs.gov）共获取29景少云无季节性积雪的Landsat TM/ETM+/OLI遥感影像（表1）。为提高影像空间分辨率并保留其多光谱信息，对Landsat OLI影像的可见光波段与全色波段进行融合，融合后的影像空间分辨率为15 m［25］。对于Landsat ETM+影像首先利用ENVI中的Landsat_Gapfill插件进行坏损条带修复，然后对所有影像进行辐射定标（Radiometric Calibration模块）和大气校正（FLAASH Atmosphere Correction模块）等预处理，将影像中绿光（Green）、近红外（NIR）和短波红外（SWIR）波段的数字量化值转换为反射率，热红外数据转换为辐射率并计算出亮度温度和黑体辐射亮度，最后使用普朗克公式计算出地表温度［19，22］。

表1 本研究所用的Landsat卫星遥感影像信息Table 1 Information of Landsat satellite remote sensing images used in this study

式中：K1、K2为常数；Bs为黑体辐射亮度值。

用于提取山脊线和坡度的数字高程模型（DEM）数据为ASTER GDEM V2，从地理空间数据云（http：//www.gscloud.cn）下载获得。用于计算表碛区表面高程变化的DEM数据为SRTM V3和2017年9月 的ASTER Digital Elevation Model（AST14DEM）V003产品，分别从USGS网站和NASA EARTHDATA网 站（https：//earthdata.nasa.gov/）下载获得，空间分辨率均为30 m。ITS_LIVE（Inter-Mission Time Series of Land Ice Velocity and Elevation）数据从美国国家冰雪数据中心（National Snow and Ice Data Center，NSDIC）网 站（https：//nsidc.org/apps/itslive）下载获得，用来分析2013年、2016年和2019年9条冰川中流线上冰川流速的年际变化，空间分辨率为240 m［26］。气温和降水数据选择距研究冰川最近的阿克苏、皮山、波密和林芝四个气象站2000—2019年中国地面气候资料年和月数据集作为冰川区气候背景参考，该数据由国家气象科学数据共享服务平台（http：//data.cma.cn）提供。

2.2 冰川表碛识别及边界提取

虽然在一些研究中冰川表碛遥感识别取得了较好的效果，但尚未有普适性较好的表碛覆盖型冰川提取方法［27-28］。借鉴前人在冰川表碛遥感识别中采用的参数或指标，本文将预处理后的地表温度数据、NDDI指数、波段比值（NIR/SWIR1）和坡度数据相集成提出TDSI方法，尝试实现对冰川表碛的识别和边界的提取，具体流程如下：

（1）基于ASTER GDEM V2数据提取坡度，采用阈值24°［21］将冰川发育区分为表碛覆盖区（≤24°）和非表碛覆盖区两类［图2（a）］。

（2）计算冰川地表温度，其最佳阈值介于9～15℃，通过人工交互式方法得出一个最优阈值，分类得到表碛覆盖型冰川的二值图，阈值不同边界范围亦不同，图2（b）所示阈值为10℃。若使用的地表温度数据无法正确提取出表碛覆盖区，则采用NDDI指数进行提取［图2（c）］，阈值的取值范围为0.95～0.99。

式中：ρSWIR2为短波红外波段的反射率；εTIR为热红外波段的辐射率。

（3）对步骤（1）和（2）结果求交集，得到二值图［图2（d）为图2（a）和图2（b）的交集］。

（4）使用NIR/SWIR1和人工交互式方法提取裸冰区，最优阈值选取介于［0.0001，0.0003］［图2（e）］。

（5）对步骤（3）和（4）获取的结果求并集，得到表碛覆盖型冰川范围［图2（f）］，将其重分类为二值图后使用中值滤波器（窗口大小为3×3，单位为像元大小）剔除与冰川不相连的噪声像元［21］，然后将其转化为矢量多边形，剔除掉面积＜1 km2的多边形，得到表碛覆盖型冰川边界［图2（g）］；经平滑处理后利用山脊线自动提取方法［29］获得各冰川山脊线矢量数据，以此对冰川边界进行分割裁剪，从而得到单条冰川的矢量边界数据，并对提取的冰川边界进行人工目视修订［图2（h）］，消除阴影、积雪，以及其他干扰因素对冰川边界提取的影响，例如图2（g）所示将出水口部分误识为冰川，经人工目视修订后以出水口的源头作为冰川末端［图2（h）］；最后由图2（f）提取出表碛覆盖区域边界［图2（i）］。

图2 冰川及其表碛的提取流程Fig.2 Extraction procedure of glacier and its debris

2.3 冰川提取误差评估

冰川边界遥感解译结果的检验可以通过地面调查或甚高空间分辨率的遥感数据分类结果作为参考数据［22］。鉴于研究区地处自然条件较为恶劣的高海拔地区，实地考察难度较大，因此本文仅考虑由遥感影像空间分辨率所造成的误差［30-31］，可由下式计算得出。

式中：ε为影像空间分辨率造成的冰川面积误差；N为冰川轮廓的周长；A为半个像元的边长。其中，9条冰川由Landsat TM/ETM+/OLI（空间分辨率分别为30/30/15 m）遥感影像空间分辨率造成的冰川面积误差分别为75.60 km2、46.32 km2和31.16 km2，合计为153.08 km2，占冰川总面积的8.66%。

3 结果与分析

3.1 表碛覆盖型冰川提取精度评价

为验证基于TDSI方法提取的冰川及其表碛范围精度，选用2011—2020年同时段Landsat卫星假彩色合成影像对相应冰川的中国第二次冰川编目数据［24］进行人工修订，并采用冰川目视解译结果与提取结果比值的均值作为评价指标。对9条冰川试验结果表明，TDSI方法提取的冰川边界平均精度为91.23%，大陆型和海洋型冰川的提取精度分别为91.20%、90.97%；表碛覆盖区边界提取精度为91.87%，大陆型和海洋型冰川的提取精度分别为93.24%、87.22%（表2），表明TDSI方法在提取冰川表碛范围方面具有较高的准确度。

表2 大陆型和海洋型冰川及其表碛提取结果和精度Table 2 Results and accuracy of continental and maritime glaciers and their debris

3.2 2011—2020年典型表碛覆盖型冰川变化

基于TDSI方法提取2011—2020年9条表碛覆盖型冰川边界，对比显示2011—2020年典型表碛覆盖型冰川总体呈现微弱退缩趋势，而表碛所占面积比例均有所上升（表3）。6条大陆型冰川总面积从1 360.77 km2减少至1 358.00 km2，面积变化相对速率为-0.02%·a-1，平均每条冰川面积减少0.06%；其中裸冰区总面积减少22.82 km2，而表碛区总面积由168.19 km2增 加 到188.24 km2，增 加 幅 度 为11.92%。3条海洋型冰川总面积从408.2 km2减少为407.85 km2，共减少0.35 km2，平均每条冰川萎缩0.11%；其中裸冰总面积从361.05 km2减少为352.05 km2，表碛区总面积从47.15 km2增加到55.8 km2，增加幅度为18.35%。相比较而言，海洋型表碛覆盖型冰川面积和裸冰面积退缩均较大于大陆型表碛覆盖型冰川，且表碛区面积增加更为显著。

表3 2011—2020年大陆型和海洋型冰川变化Table 3 Variations of continental and maritime glaciers from 2011 to 2020

在这9条表碛覆盖型冰川中（图3），音苏盖提冰川和雅弄冰川面积减少最多，分别为0.45 km2和0.27 km2；那龙冰川无面积变化，其余冰川面积减少介于0.08～0.21 km2之间。其中，木扎尔特冰川、木斯塔冰川和音苏盖提冰川裸冰面积减少最多，分别为5.88 km2、4.61 km2和4.58 km2，穹特连冰川最少（-0.52 km2），其余冰川裸冰面积减少介于2.29～3.95 km2之间；木扎尔特冰川和音苏盖提冰川表碛面积增加最多，分别增加5.79 km2和4.13 km2；穹特连冰川表碛面积增加最少，仅0.42 km2，其他6条冰川表碛面积增加介于2.21～3.84 km2之间。不同时期的冰川表碛覆盖面积呈现波动变化趋势，如穹特连冰川2020年表碛覆盖面积最大（18.93 km2），2014年最小（15.50 km2）；木扎尔特冰川2019年表碛覆盖面积最大（25.33 km2），2014年面积最小（18.17 km2）。

图3 2011—2020年各冰川裸冰区和表碛区面积变化Fig.3 Area variations of clean ice and debris for each glacier from 2011 to 2020

3.3 2011—2020年典型表碛覆盖型冰川表碛区变化

以200 m高差为间隔对9条冰川海拔进行统计，结果显示2011—2020年大陆型冰川和海洋型冰川表碛覆盖面积变化在不同海拔上差异显著（图4）。2011年6条大陆型冰川表碛覆盖区介于海拔2 800～5 200 m之间，海拔3 800 m以下冰川表碛面积大于裸冰面积，尤其是海拔3 200 m以下冰川表碛面积占比达到100%［图4（a）］；随着冰川中值海拔由4 950 m上升至5 020 m，至2020年这6条大陆型冰川表碛覆盖范围的海拔区间扩大至2 800～5 800 m［图4（b）］，除海拔3 000 m以下表碛面积略有减少外，3 400 m以下冰川表碛面积占比达到100%，同时3 000～5 800 m各海拔区间冰川表碛面积均有所增加，尤其是5 400～5 800 m区域出现冰川表碛。与6条大陆型冰川其表碛主要分布在中值海拔以下区域不同的是，3条海洋型冰川的表碛在空间上分布范围更广，其中值海拔以上表碛面积占比（2020年为15.13%）明显高于大陆型冰川（2020年为2.74%）。由图4（c）～（d）可知，3条海洋型冰川表碛分布在海拔3 200～6 600 m之间，海拔4 000 m以下冰川表碛面积大于裸冰面积，表碛面积占比达到100%的海拔由2011年的3 400 m上升至2020年的3 600 m；同时海拔4 000 m以上地区表碛面积均有所增加，尤其是4 800 m以上表碛面积占比由2011年的15.85%增加为2020年的22.23%。

图4 2011年和2020年典型表碛覆盖型冰川裸冰与表碛覆盖区的海拔分布Fig.4 Elevation distribution of clean ice and debris of typical debris-covered glaciers in 2011 and 2020

4 讨论

4.1 典型大陆型和海洋型表碛覆盖冰川变化原因

气温和降水是影响冰川形成、发育和变化的主要气候因子，气温主要影响冰川消融，降水则在一定程度上决定冰川的物质补给量［1，32-33］。研究表明，夏季平均气温每上升1℃造成的冰川消融需要增加40%～50%的降水量才能弥补［34］，显然冰川受温度的影响更大，且气温升高导致降雪无法凝结成长时间的冰层，这促进了冰川的快速消融［35］。基于9条冰川附近气象站点的夏季平均气温和年降水量数据可知，2000—2019年大陆型和海洋型冰川区夏季气温均呈上升趋势（图5），升温速率分别为0.03℃·a-1和0.05℃·a-1；年降水量也呈增加趋势，分别为0.42 mm·a-1和2.05 mm·a-1。除那龙冰川较为稳定外，2011—2020年其他8条冰川均呈缓慢退缩态势，表碛面积增加了28.70 km2，这表明气温上升是典型冰川消融退缩和表碛面积增加的主要原因。其次，海洋型冰川区降水呈现更明显的增加趋势，冰川物质积累源增加更多，冰川顶部更易发生破碎和崩塌［18］，继而形成新的表碛覆盖物，因此表碛面积变化率较大陆型冰川更大。

图5 2000—2019年研究区夏季平均气温和年降水量变化Fig.5 Variations of average air temperature in summer and annual precipitation in the study area from 2000 to 2019

冰川运动产生的岩石碎屑物搬运受其表面速度控制［36］，因此冰川流速减缓在一定程度上也可增加表碛覆盖面积。2011—2018年大陆型冰川和海洋型冰川流速显著减缓，考虑两者差异，分别选取木斯塔冰川和夏曲冰川进行分析（图6）。其中，2011年木斯塔冰川中流线冰流速度最高，平均速度为39.68 m·a-1，至2016年 和2018年 分 别 减 缓 至30.56 m·a-1、23.06 m·a-1；夏曲冰川2011年平均速度为64.74 m·a-1，2016年平均速度为60.25 m·a-1，2018年减少为59.66 m·a-1。表碛覆盖型冰川上部积累区流速缓慢，但在消融区随着距冰川末端距离越近，流速由大逐渐减小，顺着冰川被运输的岩石碎屑物在消融区和末端大量堆积［37-38］。受流速减缓影响，输送的物质开始在距末端更远的区域积累，导致表碛覆盖面积逐渐向冰川上部扩张［18］。其次，海洋型冰川流速高于大陆型冰川，其上部输送的物质弥补了由末端终碛湖扩张所引起的表碛覆盖面积减少，使其表碛面积仍呈增加趋势。

图6 2011年、2016和2018年木斯塔冰川和夏曲冰川表面流速随末端距离变化Fig.6 Variations of surface velocities of Musita Glacier（a）and Xiaqu Glacier（b）with distance from glacier terminus in 2011，2016 and 2018

冰川表碛在向高海拔地区延伸同时，受冰川退缩影响，尤其是冰川终碛湖的扩张使得表碛面积增幅有所减少。湖水与冰川末端频繁热量交换加剧冰川退缩［39］，而发生剧烈变化的冰川消融区被表碛所覆盖，表碛物随着冰川末端消融或崩解入湖。在9条表碛覆盖型冰川中，雅弄冰川末端发育规模较大的终碛湖。人工目视解译结果显示，雅弄冰川的终碛湖处于持续扩张状态，2011年、2016年和2018年 面 积 分 别 为2.74 km2、2.82 km2和2.91 km2（图7）。与雅弄冰川表碛面积变化（2011—2018年为2.60 km2）相比，终碛湖扩张（0.17 km2）使得冰川表碛面积减少有限。

4.2 典型大陆型和海洋型表碛覆盖冰川高程变化差异

相较于面积、长度等参数，表碛覆盖型冰川表面高程对气候变化的响应更为敏感［40-41］。2011—2020年，托木尔冰川和那龙冰川面积变化甚小，二者表碛区面积占比亦为这6条大陆型和3条海洋型表碛覆盖冰川的最大值。为探讨大陆型和海洋型冰川表碛覆盖区高程变化差异，本文选择托木尔冰川和那龙冰川作为代表，首先对二者表碛分布区的SRTM DEM数据（采集时间为2000年2月）和AST14DEM数据（采集时间为2017年9月）进行配准处理［42］，然后通过栅格运算获得其表碛区高程变化。由图8可知，托木尔冰川表碛区表面高程变化介于-30～0 m，平均减薄（13.18±5.88）m；那龙冰川表碛区表面高程变化介于-50～0 m，平均减薄（23.14±6.26）m。与同一海拔高度带裸冰区相比，两种类型冰川的表碛区表面高程变化呈现不同特征，如托木尔冰川表碛区周围的裸冰表面高程变化明显大于表碛区，二者平均差值约2.53 m，而那龙冰川则相反，尤其在海拔4 600 m以下，表碛区高程变化幅度较裸冰区平均增加约8.34 m。这意味着表碛的存在对大陆型冰川消融起着一定的抑制作用，而对海洋型冰川则有助于加速消融。已有研究表明，与大陆型冰川相比，海洋型冰川的表碛厚度和热容量较小，内部热传导效率较高，从而导致其厚度减薄趋势更为显著［10，43-45］。因此，尽管表碛覆盖型冰川面积变化很小，但是其表面高程变化足以证明表碛覆盖型冰川正处于消融状态，且相同海拔高度带海洋型冰川消融更强于大陆型冰川。

图8 托木尔冰川和那龙冰川裸冰与表碛区高程变化Fig.8 Elevation variations of clean ice and debris in Tomur Glacier（a）and Nalong Glacier（b）

5 结论

基 于2011—2020年Landsat TM/ETM+/OLI遥感影像和ASTER DEM数据，利用TDSI方法提取中国境内9条面积＞100 km2的表碛覆盖型冰川，进而分析大陆型和海洋型表碛覆盖型冰川和表碛变化原因及两者差异性，得到以下主要结论：

（1）TDSI方法综合考虑了表碛光谱、地形和地表温度特征，能有效提取表碛覆盖型冰川，总体精度可达91.23%，大陆型和海洋型表碛覆盖型冰川的提取精度分别为91.20%和90.97%。

（2）2011—2020年9条表碛覆盖型冰川裸冰区面积减少2.05%，表碛面积增加13.33%。其中，6条大陆型表碛覆盖冰川末端均发生了轻微退缩，面积平均减少0.03%；3条海洋型表碛覆盖冰川中，除那龙冰川外其他两条冰川均在缓慢退缩，面积平均减少0.04%，表明表碛的存在可在一定程度上抑制冰川退缩。大陆型冰川表碛主要分布在其中值海拔以下，约占其总面积的95.37%（2020年），尤其在海拔3 400 m以下冰川表碛面积占比达到100%；而海洋型冰川表碛分布范围更广，2020年其中值海拔以上表碛面积占比为15.13%。近10年间二者均呈现向冰川上部扩张趋势，海洋型冰川表碛增加明显高于大陆型冰川。

（3）气温上升是导致冰川消融退缩和表碛增加的主要原因，同时冰川流速变化和终碛湖的演变也对表碛变化有一定影响。尽管表碛覆盖型冰川在面积变化上不明显，但其表碛区海拔降低足以证明表碛覆盖型冰川正处于消融状态，且与大陆型冰川相比，海洋型冰川厚度减薄更为明显，冰川消融更为强烈。

致谢：感谢美国地质调查局、地理数据空间云平台及美国国家冰雪数据中心分别提供了Landsat遥感影像、DEM数据和ITS_LIVE数据。