基于深度学习和高分遥感影像的青藏高原东部融冻泥流阶地提取研究

王习敏， 黄荣刚， 徐志达， 焦志平， 江利明

（1.中国科学院精密测量科学与技术创新研究院大地测量与地球动力学国家重点实验室，湖北 武汉 430077；2.中国科学院大学地球与行星科学学院，北京 100049）

0 引言

融冻泥流阶地是一种冰缘区域斜坡受到反复冻融和重力作用形成的台阶状堆积微地貌［1］。全新世以来青藏高原边缘地带（如东部地区）近地表经历了多次冻融循环的复杂历史演化过程［2-3］，古融冻泥流阶地分布广泛，其空间分布对于重建第四纪古气候环境和多年冻土分布范围均具有重要意义［4-6］。

目前，针对青藏高原融冻泥流阶地观测较少，大范围制图研究基本空白。在先前的研究中，受到卫星遥感影像分辨率过低的限制，实地调绘和航空遥感是融冻泥流阶地制图的两种主要方式［7-11］。例如，1984—1990年，郭东信等［11］对青藏公路风火山垭口盆地融冻泥流阶地先后进行了两次地面调绘，现场测量还在很多泥流阶地站点开展过实验，如格陵兰东北部［12］、斯瓦尔巴群岛［13］、瑞典阿比斯库［14］、科罗拉多前沿山脉［15］等。自20世纪90年代开始，航空遥感与摄影测量技术逐渐用于泥流阶地提取和制图。Walsh等［10］使用ADAR-5500航空影像绘制了美国冰川国家公园的泥流阶地；Ridefelt等［7-8］利用航空摄影测量手段在瑞典阿比斯库山进行了多期融冻泥流阶地识别与分析。结合实地勘测调查和航空影像解译，陈游东［16］获取了中天山典型地区融冻泥流阶地的分布信息，并分析了其空间分布规律及对公路工程的灾害性影响。但是，以上手段只针对地块和小范围的泥流阶地，难以满足大范围制图的需求。

近年来，随着遥感技术的发展，高分辨率卫星遥感技术在融冻泥流阶地大范围识别中具有良好的应用潜力。然而，受融冻泥流阶地几何形态不一、表面覆盖多样（如裸地、草甸和稀疏植被）等因素影响，融冻泥流阶地纹理复杂，基于遥感影像的自动提取依然面临很大的挑战。快速发展的深度学习方法，如U-Net［17］，SegNet［18］和DeepLab［19］，在遥感影像地物分类和目标识别等方面取得了突破性的进展，在地球科学和遥感领域得到了广泛的应用［20-23］。其中，DeepLab V3+［24］深度学习模型，能获取多尺度上下文信息和清晰的边界，已用于青藏高原热融滑塌制图［25-26］，有望为复杂环境下基于遥感影像的大范围泥流阶地自动提取提供一种新途径。

因此，本文提出了一种基于DeepLab V3+和高分辨率遥感影像的融冻泥流阶地提取方法，并选择甘孜州新都桥周边为研究区，在青藏高原开展融冻泥流阶地的遥感自动识别研究。首先，利用研究区Google Earth高分辨率遥感影像，通过人工勾勒方式绘制了泥流阶地几何边界，建立深度学习的样本数据集；然后，构建DeepLab V3+泥流阶地深度学习模型，并将训练好的模型用于高分辨率遥感影像进行泥流阶地自动提取；最后，对提取结果进行精度评定，并分析了泥流阶地空间分布特征及其与地形因子的关系。本研究为融冻泥流阶地的大范围遥感制图提供了新思路。

1 研究区与数据

研究区位于四川省甘孜藏族自治州，如图1，总面积为900 km2。该区域平均海拔3 945 m；年降水量约553 mm，降水量主要集中于6—7月，日降雨量最高达到50 mm；温度范围为-15～28℃，最高温度分布于7月，最低温度主要在1月和12月。该区域地表覆盖主要是草地、少量灌木和农田等。区域内融冻泥流阶地分布较多，主要集中在318国道、215省道沿线。

图1 研究区位置Fig.1 Location of the study area

本文选用Google Earth遥感影像作为实验数据，空间分辨率为0.59 m，数据采集时间为2016—2019年。此 外，选 用2008年 的12.5 m分 辨 率ALOS PALSAR DEM数据，用于泥流阶地的地形因子分析。2020年30 m分辨率的地表覆盖数据用于分析融冻泥流阶地与地表下垫面的关系，数据来源于国家基础地理信息中心全球地表覆盖数据产品服务网站。

2 研究方法

2.1 DeepLab V3+模型

本文采用最近提出的DeepLab V3+网络，开展融冻泥流阶地遥感自动提取，该深度学习语义分割模型将深度卷积神经网络（DCNN）与概率图模型结合，取得了比DCNN更好的检测精度，并将在空间金字塔池化结构（ASPP）和解码器部分（Encoder）引入深度可分离卷积降低网络的计算复杂度［27］，提高了模型学习的效率，已在遥感图像分割和目标识别中取得了广泛应用［27-28］。

DeepLab V3+网络结构如图2所示，采用下采样编码和上采样解码结构。在编码阶段，通过骨干网络，如Xception等［29］，结合多尺度空洞卷积，实现多尺度深度特征感知。在解码阶段，通过融合骨干网络低层次特征和空洞卷积高层次特征，并采用简单的上采样操作，逐步提高语义分割的分辨率，实现与输入图像一致大小的分割结果。其流程主要包括：采用骨干网络（如Xception65［29］）进行深度可分离卷积运算，然后通过ASPP模块获取多尺度特征（空洞率分别为1，6，12，18），并将ASPP空洞卷积层特征与全局池化层进行特征的融合，以此得到多尺度的高级语义特征输出结果；将从主干网络中提取的低级语义特征和编码阶段提取的高层次多尺度特征图进行融合，并采用3×3的卷积运算和4倍上采样，提高其空间分辨率。

图2 DeepLab V3+网络结构Fig.2 Structure of DeepLab V3+

2.2 基于DeepLab V3+的泥流阶地自动提取

基于DeepLab V3+的泥流阶地自动提取处理流程如图3所示。首先，训练样本准备，制作泥流阶地数据集。然后，设置参数，对深度学习网络Deep-Lab V3+进行训练。其次，利用训练好的模型对研究区域Google Earth高分辨率影像进行预测、后处理。最后，对提取的泥流阶地进行精度评估。

图3 基于DeepLab V3+的泥流阶地提取流程图Fig.3 Flowchart of solifluction terraces extraction based on DeepLab V3+

2.2.1 多类型泥流阶地正负样本制作

综合泥流阶地明暗相间的纹理特征、长条状几何形态、集群分布等特点，结合Google Earth三维地形分析，采用ArcGIS矢量化工具勾勒泥流阶地前沿陡坎边缘线作为样本。在区域A影像上勾勒泥流阶地陡坎边界，共计4 752条，主要包括高陡坎强纹理的泥流阶地、低陡坎弱纹理的泥流阶地、植被覆盖的泥流阶地和裸地覆盖的泥流阶地4种类型，如图4（a）～4（d）所示。此外，本文还勾勒出易混淆的其他目标作为负样本，如山脊线、山谷线、田地间陡坎、冲沟等线状目标区域，如图4（e）～4（h）所示，共计110个多边形。最终勾勒的泥流阶地和非泥流阶地样本分布如图1所示。基于人工勾勒的正负样本矢量多边形，结合对应的Google Earth高分辨率影像，生成大小为513×513的正负样本切片数据集，所得正负样本切片数量分别为3 046和1 863。

图4 泥流阶地人工勾勒正负样本示例Fig.4 Examples of manually delineated positive and negative samples：solifluction terraces with high risers and strong textures（a）；solifluction terraces with low risers and weak textures（b）；solifluction terraces covered by bare soil（c）；solifluction terraces covered by vegetation（d）；ridge lines（e）；valley lines（f）；raised paths through farmlands（g）；gullies（h）

2.2.2 DeepLab V3+模型训练与泥流阶地提取

基于生成的泥流阶地数据集，对DeepLab V3+模型进行迭代训练，获取其模型参数。在训练过程中，参数设置如下：批处理量为12；迭代次数为6×104次；初始学习率为1×10-3；权重衰减系数设置为4×10-5；采用的优化器是Adam；骨干网络为Xception65；ASPP结构中空洞率分别为［1，6，12，18］，ASPP模块可有效提取不同尺度下的语义特征，一定程度上有助于提高不同尺度泥流阶地的提取效果。

在模型训练好后，将其用于整个实验区的融冻泥流阶地预测提取。首先，将实验区所有Google Earth高分辨率影像切割成513×513大小的切片；然后，利用训练好的模型对各切片进行泥流阶地提取，将切片分割成泥流阶地陡坎和其他背景像素。最后，使用GDAL进行镶嵌处理，输出完整的预测结果。

为了剔除错误提取的泥流阶地，结合几何特征、集群存在特征、与等高线相对平行的走势关系等先验知识，对预测结果进行后处理。后处理的规则主要包括：（1）剔除面积小于40 m2、长度小于10 m、宽度小于4 m的提取目标；（2）删除条数少于10条的零散分布提取目标；（3）去除个别与等高线走势相悖的目标。

2.3 精度评估

为评估泥流阶地提取方法的性能，采用人工勾绘的结果作为参考，分别根据公式（1）、（2）和（3）计算以下三个统计指标对提取结果进行精度评估［30］。

式中：P（Precision）和R（Recall）分别为精确度、召回率，F1是综合两者的度量指标，TP、TP和FN分别正确提取、错误提取、漏提取的泥流阶地陡坎面积。

3 结果与分析

3.1 泥流阶地提取结果

图5为基于深度学习的融冻泥流阶地深度学习提取结果，共9 203条，总面积为2.12 km2。据提取结果发现，在本研究区内，泥流阶地沿山谷线两侧分布，主要位于国道等不同等级道路两侧，呈现了聚集分布。

图5 基于DeepLab V3+模型的泥流阶地自动提取结果Fig.5 Results of solifluction terraces based on DeepLab V3+

对于不同覆盖类型具有明显特征的泥流阶地，本文方法在目标提取率和边界吻合程度等方面均取得很好的效果，泥流阶地的边界提取效果较好，但类型间也存在差别。如图6所示，植被覆盖的泥流阶地，条状形态、纹理不大明显的部分，提取结果有待提高；陡坎较高的泥流阶地，边界位置相对精准；陡坎纹理不清的泥流阶地，深度学习模型仅能提取到部分结果。此外，相比较于形态尺度较小的泥流阶地陡坎，形态尺度较大的泥流阶地陡坎提取效果更好。

为了定量评估泥流阶地深度学习提取质量，本文分别统计了融冻泥流阶地正确提取、错误提取以及漏提取的面积，并根据公式（1）、（2）和（3）计算深度学习自动提取结果的精度、召回率和F1（表1）。在训练区域A范围内，精度F1达到0.738；在预测区域B范围内，精度F1达到0.68。综合图5图6和表1看，本文方法总体提取精度较高，但也存在较多的错检和漏检，主要原因有：（1）算法容易错把少量梯田、植被、冲沟当成泥流阶地陡坎；（2）由于样本有限以及自动后处理面积阈值剔除等情况，在小面积泥流阶地、弱纹理泥流阶地（特别草甸覆盖区域）等区域存在漏提取。

图6 泥流阶地提取结果与人工判读比对Fig.6 Comparison between automated extraction solifluction terraces and the interpretation results：solifluction terraces covered by vegetation（a）；solifluction terraces with high risers（b）；solifluction terraces covered by meadow（c）；solifluction terraces covered by bare soil（d）

表1 融冻泥流阶地提取结果精度评估Table 1 Evaluation of extracted result of solifluction terraces

综上，泥流阶地提取方法性能分析表明，Deep-Lab V3+深度学习模型能够有效提取研究区泥流阶地陡坎的边界，且拥有良好的泛化能力，在形态、纹理相似的泥流阶地区域（如G318国道沿线）具有推广应用潜力。

3.2 实地考察验证

为进一步验证提取的结果，2020年夏季进行了实地考察［图7（a）］，考察点位于G318国道沿线，如图1中黑色定点标记所示。该处泥流阶地陡坎特征不明显，对其提取具有较大难度。但是通过自动提取结果和人工解译结果对比［图7（b）］表明，深度学习提取方法能够较好地提取大部分泥流阶地，但由于部分泥流阶地陡坎特征纹理太弱，导致部分泥流阶地无法有效提取［图7（b）］。

图7 提取结果的实地验证Fig.7 Field validation of extraction results：on-site photos（a）；results comparison（b）

4 讨论

4.1 地形因子统计分析

为分析泥流阶地的空间分布特征，本文讨论了研究区内泥流阶地与地形因子的统计关系。图8为区域A、区域B的地形因子（坡度、坡向、高程）统计图。结果显示：

（1）在研究区A内，泥流阶地坡向主要集中在西北方向［图8（a1）］。在高程统计方面，平均值为3 614.420 m，标准差为81.483 m，且主要集中在3 550～3 650 m区间，占研究区A内泥流阶地的43%

［图8（b1）］。其中，最高的泥流阶地，其高程为3 875 m，位于318国道沿线的水桥村北侧920 m附近；最低的泥流阶地高程值为3 438 m，位于营官寨村西北部2.5×103m附近。泥流阶地坡度主要集中分布在10°～30°区间，占全部泥流阶地数量的79%。

（2）在研究区B内，泥流阶地坡向主要集中在正西、西北和东北等阴坡方向，与研究区A内的泥流阶地具有类似的坡向分布特征［图8（a2）］。在高程分布方面，泥流阶地主要集中在3 650～3 750 m分布范围内，占研究区B内泥流阶地的43%；区域内泥流阶地平均高程为3 674.694 m，标准差为86.240 m。在坡度上，集中分布在10°～30°，并且在20°～25°区间的泥流阶地数量最多。

图8 泥流阶地的地形因子统计分析Fig.8 Statistical analysis of topographic factors of solifluction terraces：aspect（a），elevation（b），slope（c）of Area A and B

总体而言，本研究区域内泥流阶地以西北、北等阴坡方向为主；在坡度上，集中分布在10°～30°之间；在高程上，3 650～3 750 m区间为主要分布范围。

4.2 几何因子统计分析

此外，本文还对泥流阶地陡坎的几何因子（面积、周长、长度、宽度）进行了统计分析和讨论，如图9所示。统计结果显示：

图9 泥流阶地几何因子统计Fig.9 Statistics of geometric factors of solifluction terraces：area（a）；perimeter（b）；length（c）；width（d）

（1）在研究区A中，泥流阶地面积平均值为369 m2，标准差为330.79 m2。绝大部分泥流阶地面积小于1 000 m2，占总数的95%。其中，小于500 m2的泥流阶地占79%。泥流阶地周长平均为128.46 m，标准差为79.59 m，最大值和最小值分别为890.65 m和20.50 m。长度、宽度大部分小于500 m，均值分别为54 m、26 m。

（2）在研究区B中，泥流阶地面积均值463 m2，标准差为412.46 m2，小于1 000 m2的占91%；周长平均为146.84 m，标准差为86.36 m，100～200 m区间最为集中；长度、宽度均值分别为62 m、29 m。

总体而言，研究区A和研究区B内的泥流阶地在几何形态上具有较高的相似性，但也存在一定差异。例如，研究区B的平均长度和宽度均大于研究区A。

4.3 地表覆盖

研究区域的现今地表覆盖空间分布如图10所示，主要有草地、森林、耕地、水体、人造地表、积雪和冰川共6种。对各类型地表覆盖的泥流阶地进行统计分析，如图11所示。结果表明：泥流阶地主要分布在草地、森林、耕地三种地表覆盖类型处。区域A中的草地类型占94%；区域B中的草地类型占95%。

图10 泥流阶地地表覆盖空间分布图Fig.10 Spatial distribution map of surface cover on solifluction terraces

图11 泥流阶地地表覆盖类型统计Fig.11 Statistics on surface coverage types of solifluction terraces

5 结论

基于DeepLab V3+深度学习网络和Google Earth高分辨率卫星影像，本文提出了一种大范围融冻泥流阶地遥感自动提取的方法，并在四川省甘孜州新都桥地区开展了实验研究，分析了该区泥流阶地分布特征及其与地形因子、几何因子、地表覆盖等因素的关系。主要结论如下：

（1）与人工解译结果对比，本方法提取结果的综合精度达到0.68以上，并经野外调查验证了其有效性。陡坎较高的泥流阶地提取精度最高，对于弱纹理的泥流阶地精度较低但能实现其定位。

（2）共识别了9 203条融冻泥流阶地，主要分布在新都桥镇附近的山谷两侧，海拔高程主要分布在3 650～3 750 m区间，且主要位于20°～25°的阴坡上。

（3）泥流阶地面积和周长分别在1 000 m2、150 m以内，地表覆盖类型主要是草地，占比高达95%。

致谢：本研究使用的数据来自于Google Earth平台提供的遥感影像数据、日本JAXA提供的ALOS PALSAR DEM数据、国家基础地理信息中心全球地表覆盖数据产品中心提供的地表覆盖数据、第二次青藏高原综合科学考察项目提供的野外照片资料，本文训练样本制作中武汉工程大学王慧妮老师给予了大力支持，中国科学院精密测量科学与技术创新研究院陈圆圆同学提供了绘图调色指导，在此表示衷心感谢。