基于地形起伏度的冰湖溃决隐患研究——以希夏邦马峰东部为例

贺鹏， 童立强， 郭兆成， 涂杰楠， 王根厚

(1.中国地质大学(北京)地球科学与资源学院，北京 100083；2.中国自然资源航空物探遥感中心，北京 100083)

0 引言

在全球气候变暖背景下，近年来青藏高原暖湿化特征显著，冰川活动性增强，冰湖溃决事件频发[1]。西藏地区是我国冰湖溃决灾害高发区，以喜马拉雅山区尤甚[2]。自20世纪30年代以来，该区有文献记载的溃决冰湖33处，冰湖溃决事件37次，对我国西藏南亚大通道和川藏铁路等重要交通干线建设构成严重威胁[3-6]。冰湖溃决的形成机制复杂，诱发成因多样，其中冰崩、冰滑坡及相关组合诱因占70%[7]。冰崩、冰滑坡等灾害多发于斜坡处，由冰川断裂崩解形成[8]。目前，冰湖溃决危险性判别及评价流程已形成较为完整的体系[9]，其中斜坡地质灾害敏感性是定量评价过程中的重要指标之一[10]，而地形起伏度又是实现准确评价的前提[11]。

作为描述区域地形特征的宏观指标，地形起伏度是指某一确定面积中最高点与最低点海拔的差值[12]，可以真实反映地表起伏特征变化，与斜坡地质灾害，特别是冰崩灾害的发生显著相关[13]。近年来，基于地貌发育的基本理论，前人在不同地区和尺度下展开研究： 涂汉明等[12]利用国家数字地形数据库，基于最大高差法得出中国地形起伏度最佳统计单元为21 km2； 刘爱利等[14]基于中国1∶100万数字高程模型(digital elevation model,DEM)，探讨DEM定量化、自动化的地貌分类应用，得出在1∶100万比例尺下最佳统计单元为6.5 km2； 韩海辉等[15]应用均值变点法提取青藏高原地区地形起伏特征并分级分析，认为最佳统计窗口为1.17 km2； 张军等[16]以新疆1∶25万DEM数据为基础，利用邻域分析及均值变点法得出2.56 km2的窗口尺度为最佳统计单元； 毕晓玲等[17]基于人工判读法得出四川省地形起伏度的最佳统计单元为2.25 km2，并进一步探讨了其在区域地质灾害评价中的应用； 千琳勃等[18]利用均值变点法获取了六盘山区地形起伏度，确认最佳统计单元为1.96 km2，并分析了研究区地形空间变化特征； 张竞等[19]针对京津冀地区不同地貌样本计算了最佳分析窗口，认为地貌起伏在400 m内以4.64 km2为佳，而起伏度大于400 m的以5.35 km2更为适合； 王岩等[20]运用均值变点法，在青海德贵地区开展1∶5万比例尺下地形起伏度研究，论证了最佳分析窗口为0.16 km2。以往研究中以中低尺度下全国及西部地区的居多，以均值变点法的应用频率最高、效果最好[21]。然而，在我国藏南冰湖集中区1∶5万比例尺下的研究较少，将地形起伏度应用于冰湖溃决隐患识别方面的分析案例还较少。在当前冰湖溃决灾害频发的严峻形势下，亟须开展研究，以丰富并完善我国冰湖溃决隐患分析评价方法。本文以希夏邦马峰东部冰湖集中区为研究区，以国产高分二号(GF-2)卫星影像及ASTER GDEM V3为数据源，利用均值变点法确定地形起伏度最佳提取单元，结合冰湖分布特征开展溃决隐患相关性分析。

1 研究区概况

研究区位于聂拉木县以北，地处喜马拉雅山中段希夏邦马峰东南部区域，在行政区划上隶属西藏自治区日喀则市聂拉木县，范围涵盖波绒乡、门布乡、亚来乡及聂拉木镇等4个乡镇，总面积约730 km2(图1)。地理坐标介于E28.17°～28.5°，N85.75°～86.00°之间。

图1 研究区遥感影像

受青藏高原新构造运动影响，形成整体以高落差基岩山地夹持纵深河谷的地貌特征，区内冰湖分布广泛、重力地貌发育、各类地质灾害频发[22]。受近年暖湿化气候影响，喜马拉雅山区冰川退缩明显，冰湖显著扩张，溃决风险加剧，溃后引发链式次生灾害事件也随之增多[23]。据前人研究，希夏邦马峰东部及周边地区存在多个高溃决风险冰湖，其中嘉龙错、次仁玛错和帕曲错等已发生溃决，对下游居民的生命财产安全造成巨大损失[24-25]。

2 研究数据与方法

2.1 冰湖信息提取

采用人机交互解译方法，基于国产GF-2卫星正射级产品数据，并配合图像增强提取处理等方法，实现对区内冰湖信息的准确获取。数据由中国自然资源航空物探遥感中心提供。

2.2 地形起伏度计算

2.2.1 地形数据

本文采用ASTER GDEM V3版地形数据。该数据于2019年8月发布，覆盖了从N83°～S83°的所有陆地表面，其全球空间分辨率为30 m，垂直分辨率约为10 m，较V2版主要减少了高程值空白区域及水域数值异常[26]。数据来源于美国国家航空航天局(https: //earthdata.nasa.gov/)。

2.2.2 均值变点法

地形起伏度一般基于DEM数据通过窗口分析法求取，其关键是确定最佳统计窗口大小。以某一起伏地形上任意点为中心计算地形起伏度，随着统计窗口尺度的扩大，相对高程差与地形起伏度的变化率呈对数曲线特征，其中必定存在一个变化率由大变小(即曲线由陡转缓)的“拐点”，该点有且唯一，对应的统计窗口大小即适宜计算尺度。本文通过均值变点法确定研究区最佳统计单元，具有客观性，相关计算原理参见相关文献[15-20]，在此不再赘述。

3 适宜计算尺度确定

3.1 主要步骤

以n×n的矩形作为分析窗口，对希夏邦马峰东部冰湖集中区进行窗口尺度递增的起伏度计算，规定起止矩形窗口尺度分别为3×3个像元和65×65个像元，单次增幅为2个像元(表1)。

表1 地形起伏度分析窗口尺度设置

基于ArcGIS软件的Spatial Analyst模块下neighborhood statistics工具对区内DEM栅格数据进行不同尺度窗口下的遍历计算，得到各窗口的栅格最大及最小值； 利用Raster Calculator工具计算起伏度值，并将分析窗口大小与该尺度下起伏度值进行关联(表1)； 运用Band Collection Statistics工具进行栅格数据统计分析，并将结果属性表输出为表格； 利用VBA编程采用均值变点法计算统计量，其计算公式为：

(1)

(2)

将表1中各项数值代入式(1)—(2)中，实现对窗口面积与起伏度的曲线拟合与均值变点法中S-Sj的计算。变点的存在会使S和Sj的差距增大，S-Sj的最大值对应的分析窗口大小即为最佳分析窗口。

3.2 地形起伏度拟合曲线

前人研究表明地势起伏度与分析窗口大小的变化曲线为逻辑斯蒂曲线[12]，因此，对地形起伏度和统计单元面积(分析窗口)进行曲线拟合，如图2所示。

图2 窗口大小与地形起伏度拟合曲线

在图2的拟合曲线中，决定系数R2为0.950 5，拟合度良好，通过统计学检验。地形起伏度随统计单元尺度的扩大而增加，变化曲线表现为对数函数或幂函数特征。统计单元面积小于32.49×104m2(窗口尺度19×19)时，其增长趋势迅速； 统计单元面积在32.49×104～47.61×104m2区间(窗口尺度在19×19～23×23)时，地形起伏度的增速由急转缓； 而统计单元面积超过47.61×104m2(窗口尺度23×23)之后，地形起伏度的增速放缓并趋于平稳。

3.3 最佳窗口尺度确定

如图2所示，拟合曲线内存在一个增速逐渐趋缓的点(非数学拐点)，即所谓“变点”，该点对应的统计面积即为最佳分析窗口尺度。将起伏度计算结果逐一代入统计量公式，得到非线性系统数列样本j的S与S-Sj值，其中S-Sj最大值所对应的分析窗口大小，即所求的适宜计算尺度，计算结果见表2。

表2 均值变点法分析统计结果

根据表2数据，构建S-Sj与非线性系统数列样本j的变化拟合曲线(图3)，可以看到该曲线为近似抛物线型，在j=10时出现S-Sj的最大值，可知此处为所求的变点，即最佳分析窗口的尺度，其窗口尺度为21×21，统计单元面积为39.69×104m2。

图3 S-Sj值的变化趋势

4 地形起伏度对冰湖溃决的影响分析

4.1 基于遥感的冰湖提取结果

基于GF-2卫星数据的冰湖信息提取结果显示(图4)，区内现存冰湖多达1 020个，总面积17.67 km2。从规模上来看，可分为大型(≥1 km2)、中型(0.1～1 km2)及小型(≤0.1 km2)3类； 就类型而言，包括冰蚀湖、冰碛阻塞湖、冰斗湖及冰川槽谷湖4类。区内冰蚀湖仅3处，分别是位于希夏邦马峰东部地区的岗西错、贡错及嘎龙错，其总面积达12.21 km2，是全区仅有的3个大型冰湖； 冰碛阻塞湖18处，主要是分布在希夏邦马峰东南部地区的嘉龙错和达热措等冰湖，其中中型冰湖8处、小型冰湖10处，总面积为3.12 km2； 冰斗湖30处，均为小型冰湖，以中南部地区分布居多，总面积0.41 km2； 冰川槽谷湖969处，规模上除2处中型冰湖外其余均为小型，多分布于冰川侵蚀作用形成的谷地内。

图4 冰湖提取结果及分类

4.2 地形起伏度分级

依据上述统计分析结果，采用21×21的网格大小(39.69×104m2)对研究区内地形起伏度进行计算。为便于后续分析，进一步对区内起伏度进行分级，具体分为6个等级： 起伏度在[0,30) m为平原，[30,70) m为台地，[70,200) m为丘陵，[200,500) m为小起伏山地,[500,1 000) m为中起伏山地； [1 000,2 500) m为大起伏山地。

结果显示(图5)，区内平原分布相对较少，整体以小起伏山地为主，在拉布日峰至俄热村之间、希夏邦马峰东南侧多有分布； 丘陵及中起伏山地所占比例相近，前者多见于北部纳朵热至大来村、岗西错东部等地，后者主要分布于中部及南部冰湖集中分布区域； 大起伏山地相对较少，主要位于中西部的希夏邦马峰及西南部的多尔雷山附近。

图5 研究区地形起伏度分级

4.3 地形起伏度对冰湖溃决的影响分析

如前文所述，地形起伏度是定量化分析斜坡地质灾害敏感性的重要因子之一，更是冰湖溃决危险性评价体系中不可或缺的指标，对研究全区冰湖溃决隐患具有重要的指示意义。本文以冰湖周边地形起伏组合特征为依据，叠加冰湖空间分布特征，以分析区内地形起伏度对冰湖溃决的影响。

4.3.1 地形起伏度与冰湖溃决的关系

通过地形起伏度与冰湖空间分布特征的叠加分析发现就冰湖后缘及两侧起伏度特征而言，主要可分为单一起伏等级型、多起伏等级组合型及跨起伏等级型3类冰湖(以下简称单一型、组合型和跨级型)。

单一型冰湖周缘起伏度类型一致，多发育于冰川槽谷区，湖泊类型上以小型冰碛湖为主。区内地形起伏度较小、地势变化缓、稳定性高，不具备溃决条件。从影像特征及实地验证结果看，冲堆普流域(希夏邦马峰与多尔雷山之间)的典型冰湖(图6(a))特征与分析结果一致，作为冰川退缩消融过程中的残留湖泊，多处于同一地形起伏等级(图6(b))，周边地势平缓(图6(c))，冰湖间多由常年地表径流相连，由于湖面较浅、流动性好，整体溃决隐患极低。

(a) 单一型冰湖影像 (b) 单一型冰湖起伏特征 (c) 单一型冰湖实地调查照片

(d) 组合型冰湖影像 (e) 组合型冰湖起伏特征 (f) 组合型冰湖实地调查照片

(g) 跨级型冰湖影像 (h) 跨级型冰湖起伏特征 (i) 跨级型冰湖实地调查照片

图6 典型冰湖与地形起伏度叠加分析图及实地照片

组合型冰湖周边起伏度类型具显著差异，湖泊后缘多与冰川相接，湖泊类型以中-小型冰斗湖和冰碛阻塞湖为主。此类湖泊所处地质环境条件相对复杂，周边各地势起伏等级单元均有分布，后缘冰川易崩落入湖，形成涌浪可能诱发溃决，溃决概率较高，以嘎龙错较为典型。遥感影像及实地验证结果显示，嘎龙错周边地形特征(图6(d))与分析结果相符，属组合型(图6(e))，冰湖前缘终碛堤两侧及后缘部分区域地形高陡(图6(f))，局部具备冰滑坡、冰崩等成灾条件，可能引发冰湖溃决灾害。

跨级型冰湖周边起伏度的等级变化不连续，多为冰蚀湖、冰碛湖，规模均在中型以下，以多尔雷山东北部的嘉龙错最为典型。从遥感影像(图6(g))及实地调查结果看，该湖东侧地势较缓，而西北侧起伏度变化明显(图6(h))，与分析结果一致。冰湖侧碛垄与后缘相接处可见显著地形起伏度跨级界线(图6(i)红线标识位置)，湖面与后缘冰川具有一定高差，上部冰川易失稳发生大规模滑塌，进而诱发冰湖溃决灾害，因此这类冰湖发生溃决的可能性极高。

4.3.2 对冰湖溃决隐患的评价

基于上述分析，可进一步依据地形起伏度对冰湖溃决隐患进行分级评价(表3),其中单一型冰湖规模不一、以小型为主，由于周边地形平缓，不具备成灾条件，故整体溃决隐患低； 组合型冰湖规模以小型居多，具一定成灾条件，溃决隐患中等； 跨级型冰湖规模适中，后缘冰川所处地形落差较大，具有显著的成灾条件，溃决隐患高。区内属低溃决隐患的冰湖993处，占比97.35%，在规模上涵盖小、中、大3种类型冰湖，较典型的有大型冰湖贡错及中型冰湖达热措等； 中溃决隐患的冰湖18处，占比1.77%，规模上包括小型及大型2类，以大型冰湖岗西错、嘎龙错等最为典型； 高溃决隐患的冰湖9处，占比0.88%，包括小型及中型2种规模的冰湖，其中属中型冰湖嘉龙错溃决隐患最高(图7)。

表3 冰湖溃决隐患评价结果

图7 地形起伏度对冰湖溃决隐患的影响分析

从评价结果中不同隐患等级冰湖的数量来看，全区近98%的冰湖溃决隐患极低，仅有不足2%的冰湖存在溃决的可能； 但就湖泊规模而言，中-高隐患的冰湖面积达11.37 km2，占到全区冰湖总面积的64.31%，其对下游河道两侧牧民以及聂拉木县城构成了极大的威胁。

值得注意的是，冰湖溃决造成的损失与其规模大小并非绝对正相关。例如2014年那隆藏布支沟发生的冰湖溃决事件[6]，该冰湖溃前面积仅0.24 km2，按本文湖泊规模划分方案属中型，由于冰湖后缘冰川与常年积雪区边坡高陡，在连续降雨条件下冰川发生崩塌后冲入湖中，进而诱发溃决； 在2016年7月聂拉木县樟木镇上游发生的次仁玛错右侧支沟冰湖溃决事件中[6]，冰湖溃决前面积仅0.01 km2，属小型规模，但溃后形成洪水导致樟木口岸被淹，造成较大经济损失。

因此，本次依据地形起伏度评价出的9处高溃决隐患冰湖虽规模偏小，但在当前全球气候暖湿化的背景下，由冰川退缩消融、冰崩/冰滑坡、埋藏冰融化、强降水等诱发溃决的可能性极高，潜在威胁较大。特别是位于多尔雷山东北处冲堆普流域内的嘉龙错曾发生过2次溃决事件[27]，再次发生灾害的概率较高，需重点关注。

5 结论

本文基于GIS技术对西藏聂拉木县境内希夏邦马峰东南部冰湖集中分布区域开展地形起伏度研究，从冰湖地形起伏度组合特征的角度出发，对区内冰湖溃决隐患进行分析评价，得出如下结论：

1)利用均值变点法实现区内地形起伏度的提取，经有效计算及准确性检验，该区在1∶5万比例尺下最佳采样单元为21×21(采样单元面积为0.39 km2)，提取结果可准确反映地形变化特征。

2)研究区整体地势以小起伏山地为主，中起伏山地及丘陵次之，大起伏山地、台地及平原分布较少，多见高落差基岩山地夹持纵深河谷的地貌特征。

3)区内冰湖周边地形起伏度特征可分为单一起伏等级型、多起伏等级组合型、跨起伏等级型3种组合类型，对应溃决隐患等级依次升高，以跨起伏等级型冰湖溃决隐患最大。

4)全区1 020个冰湖中，溃决隐患为低、中、高级别的数量分别占总体的97.35%，1.77%和0.88%，中-高隐患的冰湖面积占到全区总体的64.31%，高溃决隐患冰湖中以嘉龙错特征最为典型，对下游河道两侧牧民以及聂拉木县城构成威胁。

5)结果显示，以地形起伏度作为冰湖溃决隐患评价分析指标，对完善现有冰湖溃决危险性判别方法和体系具有实际意义，可有效提高评价结果精度，为同类地区开展相关研究提供理论参考。