2001—2017年尼泊尔中部地区滑坡变化及其影响因素*

孙国庆，陈方,3†，于博，王宁

(1 中国科学院遥感与数字地球研究所 中国科学院数字地球重点实验室， 北京 100094； 2 中国科学院大学， 北京 100049； 3 海南省地球观测重点实验室， 海南 三亚 572029)

滑坡是目前造成人类伤亡和财产损失的主要自然灾害之一[1]。山体滑坡主要由降雨、温度、地震和人类活动等因素造成[2]。深入开展降雨、温度和滑坡的关系研究，可以为滑坡机理和灾害预警分析提供重要的理论支撑。尼泊尔的地理环境主要由高山、丘陵和平原组成。其北部海拔最高达8 848 m，南部最低海拔仅70 m[3]。尼泊尔不仅地形复杂多样，气候也是世界上变化最为多样的国家之一[4]。每年雨季，短时间内的暴雨极易形成滑坡、泥石流等地质灾害。随着夏季气温升高，高山冰雪融化，水体渗入斜坡，边坡不稳定性增强，滑坡发生的概率增大[5]。自2000年以来的近17年内，尼泊尔因降雨引起多处滑坡，导致数千人丧生，对人民的生命财产造成巨大威胁。其中， 2004年7月，尼泊尔遭受15年来最为严重的雨灾，对公路等基础设施造成严重损坏[6]；2011年，尼泊尔因数月连续降雨引发的山体滑坡等自然灾害造成185人死亡[7]；2015年5月、7月和8月均因强降雨引发多起滑坡，先后造成80人死亡，数百人失踪[8]。降雨引起的滑坡已经成为导致尼泊尔重大经济损失和威胁人类性命的主要地质灾害。因此，深入研究滑坡与降雨和温度的关系，不仅有助于理解滑坡发生的机理，还可以为建立有效的滑坡预警系统提供可靠的理论支撑，也能有效地用于防灾减灾研究中。

目前国内外已有许多学者对滑坡与降雨和温度之间的关系进行探索，一些研究通过对降雨诱发滑坡的资料进行分析，总结滑坡与降雨和温度之间的关系[9]。李军和周成虎[10]研究滑坡大小与降雨的关系时，得出滑坡类型与前期降雨强度有关。此类研究多是基于历史记录数据，受研究区地理位置和区域的限制，而且多数发展中国家难以统计历史数据[11]。由于遥感影像具有成像周期短、覆盖范围广的特点，能够长期获取不同空间分辨率的地物信息，可以为识别滑坡信息提供长时间序列的数据源。目前，基于遥感影像获取滑坡信息主要通过人工目视解译、基于像元或面向对象的影像分析方法和机器学习构建滑坡提取模型。人工目视解译的过程虽然费时费力，但仍然是目前获取滑坡真实数据的可靠方法[12]。面向对象的影像分析方法，虽然可以充分利用遥感影像具有相似纹理和光谱特征的像素之间的几何特征，但是在图像分割和斑块合并的过程中需要大量的人工阈值设定，很难实现滑坡自动化提取。机器学习方法在滑坡提取方面目前应用比较广泛，但是需要大量的训练样本，构建的模型也多只适用于与训练样本具有一致分布概率的测试样本中，大大限制了模型的实用性。而且在构建机器学习模型时，需要人工选择合理的特征表达滑坡体，尽可能地区分滑坡与其他背景地物，对研究人员的知识背景要求比较高。

因此，本文使用多时相遥感数据对尼泊尔中部地区的滑坡信息进行初提取，再通过目视解译进行滑坡判读识别。此方法既充分利用了遥感影像的光谱纹理等信息，也实现了一定的自动化，降低了目视解译的成本。本研究通过遥感数据处理、GIS地理统计与数理统计方法，探究2001—2017年间尼泊尔中部地区的滑坡信息与降雨和温度等气象数据及DEM高程数据的关系。

1 实验区与数据

1.1 实验区概况

本研究以尼泊尔中部滑坡频发地区为研究区(经纬度范围为84°～85.6°E，26.5°～28.5°N，图1)。研究区北部为高原山地气候，中部为温带大陆性气候，南部为亚热带季风性气候。每年夏季，受南亚季风影响，在尼泊尔中部山地地区，湿暖空气遇冷，产生降雨，易诱发滑坡灾害。

图1 尼泊尔研究区的地理位置Fig.1 Geographical location of the study area in Nepal

1.2 数据获取与预处理

1.2.1 数据获取

本文采用空间分辨率为30 m的美国陆地卫星(Landsat)遥感数据对研究区的滑坡进行分析。其中，研究区所选择的影像序列为2000年1月至2017年12月，影像条带号为140/40、140/41、140/42、141/40、141/41、141/42、142/40、142/41、142/42。土地覆盖数据来源于ESA的全球土地覆盖数据集(https:∥www.esa-landcover-cci.org/?q=node/175)，空间分辨率为300 m，本研究将尼泊尔中部地区分为林地、草地、耕地和其他用地类型。地层岩性数据来源于ISRIC的SOTER数据集(https:∥www.isric.org/ explore /soter)。地震信息数据集来源于美国地震信息中心(https:∥earthquake.usgs.gov /earthquakes/)。

尼泊尔形成滑坡的因素有很多，本文主要研究滑坡与地形、降雨和温度之间的关系。其中，地形信息可通过SRTM(shuttle radar topography mission)DEM(digital elevation model)数据(https:∥earthexplorer. usgs.gov/)获取，空间分辨率为30 m。降雨数据来源于TRMM(tropical rainfall measuring mission)3B43产品(https:∥pmm.nasa.gov/ dataaccess/downloads/trmm)，空间分辨率为0.25°，时间分辨率为1 h，时间序列为2001年1月—2017年12月。气象温度数据来源于NOAA ESRL(Earth System Research Laboratory)的全球地表温度再分析数据集 GHCN CAMS(https:∥www. esrl.noaa.gov/psd/data/gridded/data.ghcncams.html)，空间分辨率为0.5°，时间分辨率为一个月，时间序列为2001年1月—2017年12月。

1.2.2 数据预处理

Google Earth Engine(GEE)是处理卫星影像和分析地理数据的可视化平台[13]，可以为本研究区提供Landsat卫星影像数据集和快速的影像预处理算法(https:∥developers. google.com/earthengine)，同时可以减少图像拼接、几何校正、辐射校正、大气校正、云掩膜等繁重的预处理工作。

本文基于Google Earth Engine平台，采用Yu等[13]提出的影像合成策略，将2000年1月—2017年12月的数据进行年尺度合成。其中，选择目标年份整年的Landsat5/7/8数据的6个波段(Landsat 8的波段2～波段7分别对应Landsat 5和Landsat 7的波段1、波段2、波段3、波段4、波段5和波段7)，在GEE平台中对相应的Landsat原始数据集按照以下步骤进行计算：

1)提取研究区内的Landsat年度原始影像数据集L。L={Lt1,Lt2,Lt3, …,Lti, …,Ltn}，其中，Lti为ti时间覆盖研究区的Landsat影像。

2)使用GEE平台上的Landsat云量计算算法和大气表观反射率计算算法，得到输入数据集的大气表观反射率数据和每个像元的云量得分。

3)以云量得分最低的像元重构目标年份的最小云量合成影像。该算法的例子可见https:∥code.earthengine.google.com/f80ead7606e2feeaa-3744cf3d03d761 d。

本文共使用5 493景Landsat TM/OLI影像数据(表1)，图1的遥感影像为2017年合成的Landsat影像。

表1 研究区Landsat卫星数据及数量表Table 1 Landsat satellite data of the study area

本研究处理的遥感数据为年尺度影像，因此将单位为mm/h的降雨数据转换为年均降雨数据(单位：mm/month)，以K/month为单位的温度数据转为以℃/month为单位，得到年均温度数据，便于后续分析。

2 研究方法

2.1 滑坡信息提取

2.1.1 选择特征初提取

本文基于随机森林的平均不纯度减少算法[14]进行特征选择，利用选取的特征进行滑坡信息初提取，不仅可以减少目视解译的工作量，还可以选择出区分滑坡和非滑坡的主要特征。随机森林由多个决策树组成，具有较好的稳定性和鲁棒性[15]。基于不纯度的方法，可以择优选出区分地物的重要特征。

我们选用遥感影像的蓝波段(blue)、绿波段(green)、红波段(red)、近红外波段(NIR)、短波红外1(SWIR1)波段、短波红外2(SWIR2)波段和NDVI(归一化差值植被指数)、NDBI(归一化建筑物指数)、MNDWI(改进的归一化差值水体指数)、SAVI(土壤调整植被指数)10个特征进行滑坡初提取实验。其中蓝波段、绿波段用于分辨土壤和植被，红波段和近红外波段是区分植被的最佳波段，SWIR1、SWIR2可以较好地减少大气层的影响，NDVI可以检测植被覆盖度，NDBI可以较好地识别城镇建设用地信息，MNDWI容易区分水体和阴影，SAVI可以减少土壤背景对植被指数的影响。本研究选用的4个遥感指数的计算公式如下，ρred为红波段的反射率，ρgreen为绿波段反射率，ρNIR为近红外波段反射率，ρSWIR1为短波红外1波段反射率，L是土壤调节系数，经过多次实验验证，L值取0.5时，效果最为明显[16]。

(1)

(2)

(3)

(4)

本研究基于GEE平台合成的年度遥感影像，参照Google Earth和高分一号卫星数据，通过目视解译勾选共约46万滑坡与非滑坡样本像素。在勾选样本时，采用均匀、代表性强的原则，同时通过对比相邻年份的影像信息，借助DEM数据，目视判断合成影像中的像素是否为滑坡。勾选的非滑坡样本包括林地、草地、耕地等地物。采用随机森林的平均不纯度减少算法，随机采样5 000个样本数据进行特征选择实验，计算每个特征的得分，重复1 000次，充分利用扰动数据集选择出具有鲁棒性的特征，得到累加得分(表2)。得分越高表示该特征越重要，经过统计，可以发现：NDBI、MNDWI、NDVI和红波段的得分较高，因此选用这4个特征对2000—2017年的遥感影像进行滑坡初提取。

表2 不同特征的累计得分Table 2 Cumulative scores for different characteristics

本文选用混淆矩阵评价特征选择初提取的结果。其中，通过Google Earth和高分一号卫星数据的目视解译，获得滑坡真值数据，并将其重采样为30 m的空间分辨率。然后，对初提取的滑坡结果进行精度评价，得到的精度均为70%～80%，主要是因为滑坡和裸土的光谱信息相似，难以保证提取精度。

2.1.2 目视解译精提取

基于滑坡初提取结果，本研究参照Google Earth和高分一号卫星数据，通过目视解译一一判读2000—2017年遥感影像中的滑坡信息。同时，删改补充勾画出尼泊尔中部地区的滑坡真实矢量图，并提取出每年新生的滑坡信息，以进行滑坡与其影响因素的研究，图2为2001—2017年滑坡解译结果分布图，图3(a)、3(b)为提取出的每年新生滑坡数目和面积的统计图。

2.2 滑坡数据格网

为了匹配降雨、温度、滑坡数据的分辨率，首先以研究区为对象，每年度创建一个数据格网，设定每个格网单元的尺寸，以研究每个数据格网内的滑坡数目、面积的变化及其影响因素。降雨数据(TRMM)的空间分辨率为0.25°，温度数据(GHCN CAMS)的空间分辨率为0.5°，因此本研究区设定每个数据网格单元的尺寸为0.05°×0.05°[17]。这样不仅可以实现最大空间分辨率的连续性，而且处理后的数据可直观显示出滑坡在研究区的分布变化情况。降雨、温度、DEM数据及土地覆盖数据均插值为0.05°的空间分辨率，以进行同尺度数据分析。

图2 2001—2017年尼泊尔中部地区滑坡解译结果分布图Fig.2 Distribution map of the landslide interpretation results in the study area from 2001 to 2017

2.3 偏相关分析

本研究利用数据格网分析年均降雨、温度与滑坡数目和面积的偏相关性，采用t检验法对偏相关结果进行显著性检验，计算每个数据格网的t值，且t0.1(1，15)=1.753、t0.05(1，15)=2.131。其中，当t值大于1.753时，表示在90%的置信水平下，呈显著正相关性。当t值大于2.131时，在95%的置信水平下，呈极显著正相关性。

图3 2001—2017年尼泊尔中部地区滑坡信息及其与海拔高度、地形坡度的关系图Fig.3 Landslide information and its relationships with altitude and slope in the study area from 2001 to 2017

(5)

(6)

(7)

3 滑坡变化特征及其影响因素分析

3.1 滑坡的时空变化特征

研究发现，尼泊尔中部地区的滑坡分布具有明显的地域特征和区域变化规律，中部山地地区是滑坡灾害的高发区，南部平原地区基本无滑坡分布。其中，巴格马蒂区的滑坡分布较其他地区密集，其北部最为密集，甘达基区的中东部和贾纳克布尔区也分布着较多的滑坡。然而，研究区内土地覆盖程度的不同，也较大地影响着滑坡的分布。其中，林地区域分布的滑坡数目最多，占总数目的58.33%，草地次之，占32.25%。

3.2 滑坡与地形高程的关系

由图3(c)可以看出：海拔在500～1 000 m的低山区，滑坡数占比为16.34%；1 000～2 500 m的中山区滑坡分布最多，占59.20%；2 500～3 500 m的高山区滑坡占18.24%。图3(d)表明：500～1 000 m海拔处的滑坡面积占比为23.98%，1 000～2 500 m的中山区占56.22%，2 500～3 500 m的高山区占14.98%。这主要是由于在500～3 500 m的山地地区，沟谷较多，斜坡坡面较大。海拔大于3 500 m的高山地区，冰川广泛分布，夏季易发生崩塌。小于500 m的平原、丘陵区，地势较为平坦，滑坡发生概率较小[19]。由此表明，滑坡主要分布在1 000～2 500 m海拔处。

3.3 滑坡与地形坡度的关系

滑坡灾害的空间分布与地形坡度也呈现着一定的关系。由图3(e)可以看出，坡度在10°～20°范围内，滑坡数占比为3.1%，20°～30°范围内占43.38%，30°～40°占41.72%，40°～50°占11.32%。图3(f)表明，10°～20°坡度范围内的滑坡面积占比为2.39%，20°～30°的滑坡面积占45.96%，30°～40°占40.03%，40°～50°占11.24%。

坡度在小于20°的平原或低丘陵处，不易形成滑坡灾害。坡度大于50°时，易发生崩塌。坡度一般在20°～40°的地段，多发生滑坡灾害，其中坡度在30°左右的滑坡数量最多，地形坡度较陡是引起滑坡灾害的重要因素之一。

3.4 滑坡与气候因素的关系

3.4.1 滑坡与降雨的偏相关分析

首先不考虑温度的影响，计算滑坡与年均降雨的偏相关系数，并采用t检验法进行显著性检验，得出显著性图4(a)、4(b)。结果表明，研究区的降雨量与滑坡数目呈正相关的面积占63.8%，呈显著正相关的面积占15.3%，呈极显著正相关占2.6%。如表3所示，滑坡数目与降雨的偏相关系数在林地区域内为0.671 5，草地为0.341 8，耕地为0.011 9，且其正相关性面积占比分别为74.2%、59.3%、48.4%。可见，海拔在2 500～3 500 m且坡度为40°～50°的林地(如博克拉城市附近)有较强的正相关，在海拔为1 000～2 000 m、坡度为25°～40°的林地(如拉梅奇哈普城市)也呈正相关。滑坡面积与降雨量呈正相关的面积占57.2%，呈显著正相关的面积占9.8%，呈极显著正相关占0.7%。其中，海拔1 500～2 000 m且坡度10°～30°的林地呈显著正相关，拉梅奇哈普城市附近正相关性较强。相比其他土地覆盖类型，林地中的滑坡变化与降雨的影响关系最强，其中坡度较大的林地，分布着较多的滑坡。因此，降雨量对尼泊尔中部地区的滑坡变化影响较强，而且，降雨对滑坡数目的影响相比滑坡面积更强。

表3 不同土地覆盖类型中滑坡信息与降雨的相关关系Table 3 Correlation between landslide information and rainfall for different types of land cover

3.4.2 滑坡与温度的偏相关分析

尼泊尔气候形态差异很大，北部山区冬季最低温度达-40 ℃，而南部平原夏季最高为45 ℃。不考虑降雨因素的影响，计算滑坡信息与温度的偏相关系数，并进行显著性检验，得出显著性图4(c)、4(d)。研究表明，滑坡数目和年均温度呈正相关的面积占54.4%，呈显著正相关的面积占8.5%，呈极显著正相关占1.2%。如表4所示，滑坡数目与温度的偏相关系数在林地区域内为0.261 7，草地为0.436 1，耕地为0.100 4，其正相关性面积占比分别为53.7%、65.1%、47.6%。其中，海拔3 000～5 500 m的草地呈显著正相关，在加德满都及西部地区的草地也呈现一定的正相关，此处海拔为800～2 000 m，坡度为25°～50°。滑坡面积与温度呈正相关的面积占51.3%，呈显著正相关的面积占6.7%，呈极显著正相关占0.5%。结果显示，只有在东北部高山地区的草地呈正相关性，大部门地区无相关性。因此，年均温度对尼泊尔中部地区的滑坡变化影响较降雨因素弱，且草地中的滑坡变化与温度的影响关系最强，均分布在高海拔地区。

图4 滑坡信息与降雨量、温度的显著性图Fig.4 The significance of landslide information with rainfall and temperature

表4 不同土地覆盖类型中滑坡信息与温度的相关关系Table 4 Correlation between landslide information and temperature for different types of land cover

3.5 滑坡与地层岩性的关系

除去降雨、温度等外部因素外，地层岩性也是造成滑坡的主要因素之一[19]。滑坡体通常发生在具有较多碎石，总体透水性较好的区域。本研究区的片麻岩、页岩、石英岩分布广泛[20]，夏季降雨频繁，极易引发滑坡灾害。

图5为每年新生滑坡数目、面积与研究区不同岩性之间的关系图。结果表明，滑坡多发生在片麻岩、页岩岩性地带。

图5 滑坡数目面积与地层岩性的关系图Fig.5 Relationships of lithology with the number(a)and the area(b)of landslide

3.6 滑坡与其他因素的关系

尼泊尔位于喜马拉雅强地震活动带，地质活动频繁，其中2011年9月18日在尼泊尔及印度锡金邦交界处发生6.9级地震，2015年4月25日尼泊尔发生8.1级地震[21]。由于地震是造成滑坡灾害的重要因素，因此2011、2015年的滑坡数明显高于其他年份。本研究在不同土地覆盖类型下，除去地震年份进行滑坡与降雨、温度的偏相关分析，得到表5所示的结果。去掉地震年份的滑坡信息与降雨、温度的偏相关系数比2001—2017年的值均有提升，表明地震对滑坡的影响较大。同时，本文选用2001—2017年尼泊尔中部地区发生的4.5级以上地震数据与滑坡数据进行叠加，图6(a)可直观显示，滑坡多分布在地震发生区域。

表5 不同土地覆盖类型中滑坡信息与降雨、 温度的偏相关系数(去除地震年份)Table 5 Partial correlation coefficients of landslide information with rainfall and temperature for different types of land cover(Remove the year of 2011 & 2015)

图6 2001—2017年尼泊尔中部地区4.5级以上地震、断裂带与滑坡的分布图Fig.6 The distribution maps of earthquakes above 4.5(a), fault zone(b) with landslides in central Nepal from 2001 to 2017

此外，断裂带也是诱发滑坡的重要因素，尼泊尔位于喜马拉雅山南麓，主要有喜马拉雅主中央断裂(MCT)、主边界断裂(MBT)和主前缘断裂(MFT)[22-23]，图6(b)为2001—2017年尼泊尔中部地区断裂带与其滑坡分布的显示图，可见，滑坡多分布在断裂带附近。滑坡的发生也受人类活动所影响[24]。其中，山区修路、搭桥和建设水利工程等均会破坏斜坡的稳定性，而且产生大量的弃石堆土。这些弃石堆土遇强降雨后极易引发滑坡灾害。

4 结论与讨论

尼泊尔中部地区的滑坡分布广泛，特别是每年雨季，易造成严重的人员伤亡和巨大的财产损失。本文充分分析滑坡信息与其高程、坡度、降雨、温度、土地覆盖类型等因素的关系。结果表明，滑坡主要发生在海拔1 000～2 500 m，坡度20°～40°处。而且降雨是诱发滑坡的主要因素，在降雨因素的影响下，林地发生滑坡灾害的概率最大。夏季高山地区冰雪融水，亦是诱发滑坡灾害的气象因素之一，其中，在温度因素的影响下，草地区域易发生滑坡灾害。岩土体是造成滑坡的物质基础，研究区片麻岩、页岩风化后，结构松散，加之断裂带、地震和人类活动的影响，遇雨易引发滑坡灾害。