联合SBAS-InSAR 和PSO-BP 算法的高山峡谷区地质灾害危险性评价

周定义，左小清※，喜文飞,2，肖 波,3，游 洪

（1. 昆明理工大学国土资源工程学院，昆明 650093； 2. 云南师范大学地理学部，昆明 650500；3. 云南交通职业技术学院公路学院，昆明 650500）

0 引 言

中国是地质灾害高发的国家，各类地质灾害给人民的生命和财产造成了巨大的损失[1]。滑坡、泥石流、崩塌作为主要的地质灾害，具有隐蔽性强、危害性大、突发性高等特点，广泛分布于中国山区和峡谷地带[2-4]。近年来，高山峡谷区地质灾害频频发生，给人民生命和财产安全带来严重威胁。因此，对高山峡谷地区地质灾害开展危险性评价研究，能够为防灾减灾事业及政府部门决策提供参考。

目前，国内外对于地质灾害危险性评价主要分为三类：经验模型[5-6]、统计分析模型[7-9]、机器学习模型[10-12]。以上三类方法证明了其在地质灾害评价中的有效性。但存在以下两个弊端：一是时效性差、识别不准确：基于收集的历史灾害数据，以前是滑坡或泥石流等地质灾害区域，可能已经被治理好的情况，如果再将其定义为地质灾害点，容易导致评价结果不准确。运用RS（Remote Sensing）技术解译得到的地质灾害点范围，无法识别潜在地质灾害，特别是在高山峡谷区。二是过多权值计算：三类评价模型都需要运用权重和信息量来进行地质灾害统计，这无疑加大专家参与的计算量。为此，亟需提出一种可行有效方法解决此弊端。

在以往研究中，笔者曾提出一种基于小基线集合成孔径雷达干涉（Small Baseline Subsets Interferometric Synthetic Aperture Radar，SBAS-InSAR）和粒子群优化-反向传播（Particle Swarm Optimization-Back Propagation，PSO-BP）神经网络算法的矿区地表沉降监测及预测模型[13]，最终证明了SBAS-InSAR 在矿区地表长时间沉降监测中的优势以及 PSO-BP 模型在矿区地表沉降预测中的有效性和合理性。近年来，众多学者运用小基线集合成孔径雷达干涉在地质灾害识别方面取得大量成功范例[14-15]，运用SBAS-InSAR 技术在地质灾害识别探测中的应用潜力已经得到国内外众多学者的认可。误差反向传播（Back Propagation，BP）学习算法作为人工神经网络（Artificial Neural Network，ANN）最为广泛运用的一种神经网络，已经成功应用于各类地质灾害危险性评价及预测等方面并证明了其有效性[16-18]。与传统的基于梯度下降的优化算法相比，粒子群（Particle Swarm Optimization，PSO）算法具有鲁棒性强、扩充性好、不易陷入局部最优等优点；而与进化规划、遗传算法等基于自然界进化过程的优化方法相比，粒子群算法的信息共享机制使得种群可以更快地收敛到最优值[19]。综上所述，利用SBAS-InSAR和PSO-BP 算法为地质灾害危险性评价提供了理论和实践基础。

为此，针对现有地质灾害危险性评价方法存在的弊端，提出联合SBAS-InSAR 和PSO-BP 相结合的方式对高山峡谷区地质灾害进行危险性评价。利用SBAS-InSAR 获取现有地质灾害点和潜在灾害的形变量，根据形变量值和引入高分辨率光学影像等作为辅助识别；利用PSO-BP 算法构建模型，通过学习、训练、测试、验证得到危险性指标，以期在减少专家参与进行大量统计计算的前提下，解决地质灾害评价时效性和准确性差的弊端。

1 研究方法与数据

1.1 研究区概况

东川小江流域为世界典型暴雨泥石流区，被称为“泥石流的天然博物馆”，该地区地质灾害类型主要以滑坡、崩塌、泥石流为主。泥石流的发生常伴随滑坡，属于暴雨型滑坡泥石流，即先滑坡，再经暴雨冲刷形成泥石流。东川区地质灾害活动由来已久，至今呈加剧发展趋势，截至2019 年4 月，东川区各类地质灾害隐患点数为269，潜在经济损失290 961 万元。本文选取小江沿线两侧高山峡谷作为研究区，如图1 所示。以小江为界，河谷凹陷，形成“V”字型，属于典型的深切割高山峡谷，东侧为牯牛寨山，最高峰海拔4 017.3 m；西部为拱王山，最高峰海拔4 344.1 m。该区域地势陡峭，独特的地形和地质构造致使局部区域暴雨多，土壤疏松，水土流失严重，进而导致滑坡、崩塌、泥石流等地质灾害频发。

1.2 研究数据

监测数据为2018 年至2019 年的40 期Sentinel-1A 雷达数据，每期1 景，用于监测研究区形变信息。Sentinel-1A的轨道类型为近圆形的太阳同步轨道，其轨道高度为693 km，倾角为98.18°。卫星重访周期为12 d，选用干涉宽模式（Interferometric Wide，IW）下的S1 TOPS-modes单视复数SLC（Single Look Complex）数据，幅宽250 km，升降轨数据的入射角分别为34.17°和39.35°。该数据为C 波段的升降轨数据，SLC 数据处理后的升降轨数据距离向分辨率为5 m，方位向分辨率为20 m。

辅助数据为谷歌卫星影像数据、坡度、坡向和曲率数据用于主被动遥感数据协同监测识别地质灾害；精密定轨星厉数据POD（Precise Orbit Ephemerides）用于提高雷达数据轨道精度；DEM（Digital Elevation Model）数据用于去除地形相位影响；降雨、断层、河流水系、地层岩性、土地利用和归一化植被指数（Normalized Difference Vegetation Index，NDVI）作为PSO-BP 模型的输入层评价因子。以雷达数据序列为基准数据集，数据源见表1。

表1 数据源一览表Table 1 List of data sources

1.3 研究方法

本研究的总体技术方案如图2 所示，主要技术流程为：1）地表形变信息获取及地质灾害识别；2）格网单元划分及评价因子选取；3）模型构建及危险指数评定。

1.3.1 地表形变信息获取及地质灾害识别

首先，依据欧空局下载的40 景升降轨道SCL 影像，利用精密定轨星厉数据POD 修正轨道信息，减小定时及轨道误差对干涉相位的影响。再次，选择同时位于时间基线中部及时序影像多普勒质心频率中心时段的数据作为主影像，对主影像和辅影像两两进行粗配准，再用增强谱分集算法进行精配准，通过设置临界基线和时间基线，生成像对，设置多视数Range Looks 和Azimuth Looks为4∶1，这样设置可以较好地抑制斑点噪声，采用Minimum Cost Flow 解缠方法和Goldstein 滤波方法做干涉工作流，最终生成干涉图，调整删除不理想的数据，对多视处理后的主从影像进行共轭相乘获得干涉图，利用DEM 对和精密定轨数据对干涉相位进行差分处理移除平地效应及地形效应，生成时间序列的干涉相位。具体包括感兴趣区选择、干涉处理、地形相位去除、差分模型构建4 个步骤。主辅影像干涉相位可表示为公式（1）[13]

式中φtopo为地形相位；φdef为形变相位；φatm为大气延迟相位；φflat为平地相位；φnoise为噪声相位。通过相位解缠提取有效的形变信息，运用最小二乘（Least Squares，LS）或奇异值分解（Singular Value Decomposition，SVD）算法得到变形速度，最后，对序列信息进行地理编码获得研究区域内升降轨下的雷达视线方向（Line of Sight，LOS）的形变信息。

利用不同轨道数据可以互补，使得监测结果更为准确全面，能够避免单一轨道带来的几何畸变等问题。研究区内植被覆盖度高，失相干严重，导致监测地面形变信息的能力较差，对形变监测精度的影响较大[22-23]，因此，引入归一化植被指数（Normalized Difference Vegetation Index，NDVI）对研究区植被覆盖数据进行分析，以剔除植被覆盖对形变监测精度的影响。同时，仅凭获取到的研究区形变监测结果，无法有效区分研究区形变监测结果是否为地质灾害，为此，引入高分辨率光学影像，即借助谷歌卫星影像对获取的形变区域进行人工识别辅助识别，结合形变范围、高程、坡度、植被覆盖和坡体是否具有地质灾害特征等进行识别，避免过度依赖形变结果导致的误判等问题。为验证地质灾害识别结果的准确性，选择野外实地勘察验证。

1.3.2 格网单元划分及评价因子选取

常用的地质灾害评价单元包括网格单元、斜坡单元、地形单元和地域单元4 种类型。对于高山峡谷地区，网格单元和斜坡单元是最好的选择。但斜坡单元具有不连续性，精细程度不高，划分的正确度没有保证。网格单元在评价过程中虽然对于地表起伏度不能很好显示，但是相较于斜坡单元具有更大优势，网格单元使整个研究区的评价结果过渡平稳，评价结果更加真实可信[24-25]。为此采用网格单元对研究区进行划分。网格单元大小划分根据薛凯凯等[26]采用DEM 试验进行回归分析而得到的经验公式（2）：

式中sG为网格单元大小，s为地形比例尺倒数。

查阅资料可知，用于地质灾害危险性评价的评价因子主要有高程、坡度、坡向、曲率、降雨量、地层岩性、土地利用类型、距断层距离、距河流水系距离[7,8,24-25]。本文在上述10 种评价因子的基础上引入NDVI 和升降轨数据获取的形变速率作为评价因子。根据识别结果选取部分地质灾害区域和选取部分为非地质灾害的区域，通过ArcGIS 多值到点工具，提取出地质灾害范围和非地质灾害的区域对应格网单元内的12 个评价因子，通过运用SPSSAU 软件中的灰色关联和因子分析，验证以上12 个因子作为评价因子的有效性和可靠性。

1.3.3 模型构建及危险指数评定

传统的BP 神经网络存在学习效率低，收敛速度慢等问题，利用粒子群算法来优化BP 神经网络的权值和阈值，可以减小训练误差，提高训练速度和精度。将识别出来的灾害和非灾害点集区域定量化处理，用“1”代表高危险区，“0”代表低危险区，运用ArcGIS 提取对应的评价因子，将其作为输入层，高低危险区作为输出层，构建PSO-BP 模型，其中隐含层和输入层激活函数分别选择tansig、purelin，训练函数选择trainlm。选取高危险区和低危险区共4 756 个点进行模型训练和验证，从4 756个点中随机选取其中的4 200 个点作为训练集，对随机抽取的556 个点进行验证，通过训练得到最优网络并保存网络模型，最终得到高、低危险区对应的输出范围。运用保存好的PSO-BP 模型对剩余的38 749 个点进行危险性指数预测，结合上述4 756 个点对应的危险性指数，将其得到的危险性指数导入到对应点属性表内，利用自然间断点分级法结合专家参与进行危险性分级，得到研究区危险性评价图。不论是经验模型、统计分析模型还是机器学习模型，都是在统计评价因子权重信息的基础上进行，同时均是基于收集的地质灾害数据进行危险性评价，为此，综合考虑，定量角度采用信息量和综合赋权法，信息量法对数据样本要求少，可解决多源数据差异性。地质灾害是在多种因素相互作用下而形成的，通过计算各个因素的信息量，并加权叠加各个因素的信息量，建立地质灾害易发性评价模型，信息量越大，地质灾害易发性越高，其计算方法如下：

式中I为评价区某单元信息量预测值；Ni为分布在因素Xi内特定类别内的灾害点单元数；N为研究区含有灾害点分布的单元总数；Si为研究区内含有评价因素Xi的单元个数；S为研究区总单元个数；Yi为致灾因子指标值。组合赋权法通过将主观赋权法与客观赋权法得到的权重值进行耦合，能够避免单一评价方法的缺点和不足。利用乘数归一法（公式（5））将层次分析法和熵权法得到的权重值进行耦合，得到各评价指标因子的综合权重值。

式中F总表示评价单元的综合权重值，Wj表示层次分析法求得的各评价因子的权重值；Vj表示熵权法求得的各评价因子的权重值。

最后，利用定性和定量的方法验证本文方法的可靠性和有效性，根据文献[21,27]通过ArcGIS 配准获取到研究区的地质灾害点，定性角度：选取姑海村、大白泥沟和大营盘3 个区域进行对比，定量角度：运用 AUC 值和准确率对3 种不同危险性评价方法得出的结果进行验证。ROC（Receiver Operating Characteristic Curve）曲线是目前应用最广泛的验证评价结果的方法之一，它表示拟合数据与实际数据的关系AUC（Area under Curve）可以反映和比较模型的评价预测精度，AUC 即ROC 曲线下的面积值介于0 到1 之间，AUC 越接近1 说明效果越好。通过SPSSAU 软件和Excel 计算得到3 种方法对应的AUC值和准确率。

2 结果与分析

2.1 地表形变信息及地质灾害识别

研究区域2018 年4 月25 日至2019 年4 月20 日雷达视线方向（Line of Sight，LOS）的形变速率图如图3、图4 所示。从图3 中可以看出，采用升轨获取的形变主要分布在小江东侧，最大形变速率为−162.05 mm/a，降轨获取的形变主要分布在小江西侧西北方向，最大形变速率为−120.43 mm/a。使不同轨道得到的形变结果不同，是由于升轨数据飞行方向大致从南到北，雷达视线方向位于右侧，能够很好地将峡谷两侧由西向东的地表形变监测出来，相反，降轨数据的飞行方向与之相反，能够将峡谷两侧由西向东的地表形变监测出来，利用不同轨道数据可以互补，使得监测结果更为准确全面，避免单一轨道带来的几何畸变等问题。

借助谷歌卫星影像对获取的形变区域进行人工识别辅助识别，结合形变范围、高程、坡度、植被覆盖和坡体是否具有地质灾害特征等进行识别，避免过度依赖形变结果导致的误判等问题。识别结果如图3、图4 所示。为更为直观地反映形变速率的实效性和准确性，将部分识别出来的地质灾害区形变速率与三维影像叠置（蓝色实线指向的图），从图中可以看到，利用升降轨对高山峡谷区地质灾害识别具有良好的互补优势，采用升轨识别的地质灾害主要分布于小江沿线东侧，利用降轨得到的结果相反。识别并不要求一次性将所有地质灾害识别出来，仅需要尽可能均匀地识别出地质灾害，为后期利用PSO-BP 模型做准备。本次识别的地质灾害包含了滑坡、泥石流、崩塌和潜在灾害点等。

通过实地勘察发现，识别出来的结果与实地勘察结果相吻合。如图3、图4 所示，灰色实线指向的图片为实地调研勘测到的，可以看出识别出的地质灾害多以滑坡、泥石流、危岩和崩塌为主。研究区地层岩性从元古界昆阳群至新生界第四系均有分布，包括沉积岩、岩浆岩、变质岩三大类。识别出来的地质灾害主要岩组类型为松散岩组。主要位于小江断裂带两侧，断裂两侧岩石挤压极为强烈，常形成数百米宽的碎粉岩、碎裂岩或糜棱岩带，两侧碎裂岩带宽约1 000～7 000 m，钙质胶结，断面倾向东，倾角为55°上下。小江深大断裂为境内地质主构造带，沿断裂挤压剧烈，岩层破碎，东支为重力负异常区，西支便于地应力集中，形成地壳脆弱地带，历史上具有长期活动的趋势，致使该区域地质灾害多发。

2.2 格网单元划分及评价因子选取

依据公式（2），经计算网格单元大小应选择30 m×30 m，考虑到制图需求以及自身电脑性能等问题，最终选择网格单元大小为100 m×100 m，共43 505 个网格单元。选取的部分地质灾害区和非地质灾害区，其中地质灾害区格网单元为2 094 个，非地质灾害区2 662 个，利用灰色关联分析，得到灰色关联排序，按关联度从大到小依次为降雨量（0.986）、地层岩性（0.982）、土地利用（0.981）、坡向（0.973）、距断裂带距离（0.972）、降轨形变速率（0.966）、高程（0.959）、距河流水系距离（0.953）、归一化植被指数（0.934）、升轨形变速率（0.926）、曲率（0.917）、坡度（0.902），关联度最高和最低的评价因子分别是降雨量和坡度，其关联度分别为0.986 和0.902。可以看出，选取的12 个评价因子与是否为地质灾害区具有高度关联。但还不能选取出最终所需要的评价因子，为此，运用SPSSAU 软件中的因子分析方法，通过因子分析来确定上述12 个因子之间是否存在共性，分别对12 个因子进行不同组合，当选取12 个因子进行组合时，其检验统计量（Kaiser-Meyer-Olkin，KMO）最佳。验证了选用以上12 个因子作为评价因子具有有效性和可靠性。

2.3 模型构建及危险指数评定

2.3.1 模型构建及训练验证

将识别出来的灾害和非灾害点集区域定量化处理，用“1”代表高危险区，“0”代表低危险区。根据多次试验得到隐含层为10 时结果最佳，构建的PSO-BP 模型训练参数如下：动量因子0.9、学习速率0.05、训练目标0.001、训练次数500 000、进化代数80、学习因子C1=2.8，C2=1.3，最大速度0.2、种群规模25。对随机抽取的556个点进行验证，其输出值与期望值如图5 所示。

从图5 中可以看到，验证结果总体可靠，虽然存在少量“异常”数据（紫色区域），并不影响其整体的预测[28]。利用构建的PSO-BP 模型得到高、低危险区的危险性指数[11]，在（0.85，1.02]为高危险区，在区间（−0.15，0.25]为低危险区。选用均方误差（Mean Squared Error，MSE）和平均绝对百分比误差（Mean Absolute Percentage Error，MAPE）作为模型的评价指标，最终计算出MSE和MAPE 分别为0.002 和1.01%，验证了该模型精度可靠，效果良好，保存训练好的网络。

2.3.2 地质灾害危险性评价

将保存好的网络运用于整个研究区，得到每个格网单元对应的危险性指数，按照文献[9-10]的危险性分级方式分为5 个等级，采用自然间断点分级法结合专家参与进行危险性分级处理，将其划分为[−0.92，−0.15]、（−0.15，0.25]、（0.25，0.85]、（0.85，1.02]、（1.02，1.5]共5个等级，分别对应极低危险区、低危险区、中危险区、高危险区、极高危险区。最终得到SBAS-InSAR 和PSO-BP模型下的研究区地质灾害危害性评价区划图如图6a 所示。利用信息量法和组合赋权法得到研究区地质灾害危险性评价，如图6b 和6c 所示。

选取姑海村、大白泥沟和大营盘3 个区域，如图6 所标示位置，采用本文方法、信息量法[7]和组合赋权法[27]进行定性比对可知姑海村分别为高危险区、高危险区、极高危险区；小白泥沟分别为极高危险区、极高危险区、低危险区；大营盘分别为低危险区、高危险区、极高危险区。大营盘属于村落，其危险性一般为低危险地区，但信息量法和组合赋权法得到的评价结果均为极高危险区，本文方法判别为低危险区。姑海村属于村落，其危险性一般为低危险地区，但信息量法和组合赋权法得到的评价结果均为极高危险区，本文方法判别为高危险区，致使本文方法判别为高危险区的原因，可能是由于姑海村位于白泥河泥石流沟对岸，所以其危险性较高。小白泥沟是典型的暴雨型泥石流沟，既有滑坡也有泥石流，其周围应为极高危险区，本文方法和信息量法均判为极高危险区，但组合赋权法为低危险区。

运用 AUC 值和准确率对3 种不同危险性评价方法得出的结果进行验证。如表2 所示。

表2 3 种方法AUC 值和准确率Table 2 AUC (Area Under Curve) values and accuracy of the three methods

从表2 可以看出，3 种方法均能有效的对地质灾害进行评价，相较于信息量法，组合赋权法和本文方法效果更优，本文方法优于其他两种方法。信息量法和组合赋权法计算量较大，给危险性评价增加了大量的专家参与与统计。本文方法从预测角度构建危险性评价，有效减少了计算量，其精度也较为可靠。

4 讨 论

目前，虽然国内外学者对于地质灾害危险性评价的研究已经取得大量成功范例，但是还存在地质灾害数据源时效性差以及需进行大量的专家参与统计等弊端。针对地质灾害数据源时效性差这一弊端，本文利用InSAR、谷歌卫星影像和其他辅助数据对高山峡谷地区地质灾害进行识别，利用升降轨获得的InSAR 形变速率，能够识别出高山峡谷区不同方向的已有活跃滑坡和潜在滑坡，利用谷歌卫星影像和其他辅助数据可以克服InSAR 识别结果不全和不准确。针对需大量专家参与统计这一弊端，本文提出利用PSO-BP 模型从预测角度，通过不同格网单元评价因子的输入来学习，得到对应格网单元的数值（危险性指数），有效的避免了大量的专家参与统计。

采用信息量法和已有学者利用组合赋权法获取的危险性评价图存在在部分评价错误的地方，导致评价错误的原因是，两种方法均是基于已有灾害点数据作为数据源进行危险性评价，其时效性差，存在以前是滑坡或泥石流等地质灾害区域，可能已经被治理好，如果再将其定义为地质灾害点，容易导致评价结果不准确，数据源的精确性直接影响着最终的评价结果。本文所用方法从定性和定量角度验证了其评价的有效性和准确性，所用方法从灾害点数据源获取角度就有本质区别，利用该方法可以识别出已有危险和潜在危险区域，提高了灾害数据源的时效性，同时运用PSO-BP 方法从预测角度构建危险性评价模型，避免了大量的评价因子权重和信息量统计，有效减少了专家参与的计算量。

5 结 论

本文通过分析现有地质灾害评价模型存在的弊端，提出利用SBAS-InSAR 和PSO-BP 算法对高山峡谷区地质灾害危险性进行评价，经过对研究区试验，得出以下结论：

1）与现有地质灾害数据收集方法相比，利用SBAS-InSAR 技术可以获取到研究区地表形变量，结合高分辨率影像等辅助信息，可有效识别出活跃的泥石流、滑坡、崩塌和潜在地质灾害。解决了现有地质灾害数据时效性差、不准确等弊端。

2）在高山峡谷地区，利用升降轨结合的方式对研究区地质灾害进行识别，其识别结果更为准确性和全面，避免了单一轨道存在SAR 成像几何畸变造成部分地质灾害不能识别、识别结果不全面等问题。同时，通过试验验证了升降轨形变速率可作为地质灾害危险性评价的重要评价因子。

3）PSO-BP 模型适用于地质灾害危险性评价，通过试验验证了该模型精度可靠，效果良好，其MSE 和MAPE分别为0.002 和1.01%。运用该模型避免了传统地质灾害评价需进行权值统计分级等弊端，从预测角度解决了过多的专家决策参与。

4）通过对比可知，本文方法相比信息量法和组合赋权法在地质灾害危险性评价中更具有准确性和有效性，利用本文方法能够为防灾减灾事业及政府部门决策提供一定参考价值。5）综合考虑选择信息量法和组合赋权法进行定量和定性比较，比较结果表明，本文方法有效的提高了地质灾害危险性评价的准确度，信息量法、组合赋权法和本文方法的AUC（Area Under Curve）值分别为0.694、0.721、0.785，准确率为73.3%、76.2%、79.8%。