基于Logistic回归和MCMC方法评价地震滑坡敏感性

2022-05-04贺倩,汪明,刘凯

水土保持研究 2022年3期

贺倩, 汪明, 刘凯

(1.北京师范大学地表过程与资源生态国家重点实验室, 北京 100875;2.北京师范大学减灾与应急管理研究院地理科学学部, 北京 100875;3.北京师范大学环境演变与自然灾害教育部国家重点实验室, 北京 100875)

地震滑坡是最常见的次生灾害之一，许多学者对地震滑坡的敏感性，危害和风险进行了研究[1]。在山区，地震引发的滑坡会掩埋村庄，阻塞道路和河流，阻碍救灾，其造成的损失可能比地震本身造成的损失更加严重[2]。因此，准确的评估地震滑坡的敏感性具有重要意义。我国西南地区地震多发，地震诱发的滑坡造成了交通阻塞，生命线中断，严重威胁着人民的生命财产安全[3]。

在滑坡敏感性评估中，Logistic回归是广泛应用的模型之一[1,4-5]。逻辑回归模型的参数估计以及参数不确定性的定量分析对于模型的构建来说十分重要。经典Logistic回归模型一般利用最大似然法来估计方程中的参数，最大似然是点估计，其对参数的不确定性估计较为困难[6-7]。马尔可夫链蒙特卡洛方法是在贝叶斯理论框架下，将蒙特卡洛模拟和马尔科夫链模型结合，实现抽样样本随模拟的进行而改变的动态模拟[8]，其可以获得参数估计的后验分布，从而利用参数后验分布来评估参数的不确定性。将MCMC方法应用于Logistic回归的参数估计能够结合马尔科夫模型长期预测的优势和Logistic模型模拟复杂系统变量和响应变量变化的能力，具有科学性和实用性[9]。MCMC方法在水文模型参数不确定性分析中应用广泛[10-12]，而在滑坡敏感性建模上的研究较少。此外，随着计算机技术的快速发展，基于MCMC方法估计Logistic模型中的回归系数易于实现[12-13]，特别是R语言的应用使得MCMC方法实际应用更加简便可行。

本文以2017-8-8九寨沟地震，2014-8-3鲁甸地震和2013-4-20芦山地震为案例进行研究。本文从地震动参数，地形以及地质三方面构建滑坡敏感性评价指标体系。基于MCMC和Logistic回归方法建立地震滑坡敏感性的评价模型，对模型的精度进行评价，然后分别绘制三次地震事件的滑坡敏感性图。通过MCMC获得的参数后验概率分布来量化参数估计值的不确定性。最后对区域影响因子的重要性进行评价。

1 研究区与数据源

1.1 研究区概况

中国西南地区地处地震带，该地区地壳活动频繁，地震多发，再加之地形起伏大，由地震引发的滑坡灾害严重。本文选择了位于该区域的三次地震事件作为研究案例。九寨沟地震发生于2017年8月8日，位于九寨沟国家公园，震级为7级，震源深度20 km，地震引发了各种类型的地质灾害，包括滑坡、溃坝和大坝毁坏等，共造成25人遇难[2]。鲁甸县位于云南省西北部，为中国地震最活跃的地区之一[14]。2014年8月3日，云南省昭通市鲁甸县发生6.5级地震，造成了617人死亡[15]，引发了大量滑坡[16]。芦山位于四川平原与青藏高原的交界处，横跨龙门山构造带和四川盆地两个构造单元。2013年4月20日，四川雅安芦山县发生7.0级地震，震中位于青藏高原最东缘，龙门山断裂的西南部，地震造成了217人死亡，诱发了数以千计的滑坡，造成了严重的破坏并严重阻碍了救援行动[17-18]。

1.2 数据源

1.2.1 滑坡编录数据滑坡编录数据对于评估区域规模的滑坡灾害或风险至关重要。本文的地震滑坡编录数据包括：九寨沟地震诱发的226处滑坡[19]，鲁甸地震诱发的滑坡1 024处[15]，芦山地震滑坡15 546处[20]。滑坡数据是研究学者根据地震前后的高分辨率的卫星影像进行系统准确的解译得到的。本文所使用的滑坡数据类型为点。本文采用Zhu等[21]提出的双缓冲区采样方法，在内部缓冲区(1 km半径)和外部缓冲区(15 km半径)之间的区域进行随机采样，分别获得每个事件等量的非滑坡样本点。然后将每个事件中滑坡和非滑坡点随机分成70%的训练样本用于模型的构建，30%的测试样本用于模型精度的验证。

1.2.2 指标选择参照前人的研究和数据的可用性[4, 22-25]，本文从地形，地质和地震动参数3个方面共选择了9个指标，包括高程，坡度，地形湿度指数(Topographic Wetness Index，TWI)，矢量稳健性测度(Vector Ruggedness Measure，VRM)，平面曲率，剖面曲率，距离断层的距离(简称断层距)，岩性，修正麦卡利烈度(Modified Mercalli Intensity，MMI)。地形数据为DEM，来源于SRTM (Shuttle Radar Topography Mission) DEM (http:∥srtm.csi.cgiar.org/)，空间分辨率为90 m×90 m。利用DEM数据和SAGA GIS平台提取得到研究区的高程、坡度、TWI，VRM、平面曲率和剖面曲率。地质数据包括断层和岩性数据。地震滑坡的空间分布与断层之间存在着较强的相关性[26]。本文所使用的断层数据来自GEM(Global Earthquake Model)全球活动断层数据库(the GEM Global Active Faults Database, GAF-DB)，GAF-DB汇总了许多地区的数据集，包含13 500个断层，空间完备性高[27]。岩性数据来源于高分辨率的全球岩性图(Global Lithology Map, GLiM)数据集[28]，该数据集的岩性类别分为3个级别，本文使用第一级的分类(16类)，包括松散沉积岩，碎屑沉积岩，火山碎屑岩等。地震动参数数据来自美国地质调查局(USGS)的ShakeMap(https:∥earthquake.usgs.gov/data/shakemap/)，ShakeMap提供了地震后地面运动和震动强度的近实时地图。本文使用ShakeMap中的修正麦卡利烈度产品。

2 方法与模型

2.1 Logistic回归

Logistic回归是分析滑坡敏感性的有效的方法[5,29]，其因变量是二分类变量(滑坡或者非滑坡)，解释变量(自变量)可以是离散的类别型，也可以是连续的数值型。Logistic回归的优势是变量不需要满足正态分布[30]。Logistic回归通过构建二元变量和多个自变量之间的回归关系来确定事件发生的概率，公式如下：

(1)

Z=β0+β1x1+β2x2+…+βkxk

(2)

式中:P表示单个栅格像素内滑坡发生的概率；xk表示自变量；β0表示算法的截距(常数项)；βk表示不同滑坡因子的回归系数；k为滑坡影响因子的个数。在本文中,Logistic回归的因变量指定为1,0,分别表示滑坡和非滑坡。模型输出在0～1之间表示滑坡发生的概率,即本文中的滑坡敏感性。

2.2 马尔可夫链蒙特卡罗

MCMC使用随机游动的Metropolis算法从Logistic模型的后验分布中进行模拟。其通过因变量和自变量构建的贝叶斯框架来解释先验信息。MCMC参数估计中包含3个关键部分：先验分布、似然函数和后验分布。Das等[7]建立了一个简单的频率Logistic模型的Bayes等价物，如下所示：

yi～Bernoulli(logit-1(η))

(3)

(4)

βj～N(0,0.00001)j=0,…,k

(5)

式中:yi表示因变量；βj是具有独立正态先验分布的系数;xj表示第j个变量的值；η是线性预测变量。利用贝叶斯理论,可以获得模型下参数的后验分布,利用贝叶斯公式进行扩展,其后验分布与先验分布和似然函数的乘积成比例关系[7,9]。

MCMC结合参数的先验分布通过迭代仿真得到参数估计值的后验分布和概率密度函数。本文基于R语言中的MCMCpack包进行模型参数的估计[31]。在每个事件中，使用70%的样本通过MCMClogit函数来拟合Logistic方程，其余30%来进行模型预测精度的验证。

2.3 模型精度评价

混淆矩阵是分类问题中模型精度评价的重要方法[32]。在二分类问题中，根据真实类别以及模型预测结果的组合划分为真正例(True positive，TP)，假正例(False positive，FP)，真反例(True negative，TN)以及假反例(False negative, FN)。在本研究中，正例为滑坡，反例为非滑坡，分别用1，0表示。TP和TN分别表示正确分为滑坡和非滑坡的数量，FP和FN分别表示错误地分为滑坡和非滑坡的数量。基于混淆矩阵计算出总精度(Accuracy)，精确率(Precision)和召回率(Recall)。精确率(Precision)在滑坡二分类问题中指的是在所有滑坡样本中预测正确的滑坡数比例，召回率(Recall)指的是预测为滑坡的样本点中真正的滑坡数比例，F1-score是基于召回率和精确率的调和平均定义的，计算公式见(6)—(9)。

(6)

(7)

(8)

(9)

此外，本文还计算了受试者工作特征曲线(Receiver Operating Characteristic，ROC)和曲线下面积(Area Under the Curve, AUC)。ROC反映了特异性和敏感性之间的关系，AUC反映模型的总精度，AUC范围为0.5～1，AUC值越大表明模型的预测能力越高。

3 结果与分析

3.1 指标共线性分析

指标之间的多重共线性会降低模型预测能力[33]。在建立模型前，本文利用容差(Tolerance，TOL)和方差膨胀因子(the variance inflation, VIF)[34-35]进行因子共线性的检验以确保因子之间相互独立。通常认为，当VIF值大于5时，因子之间存在强共线性，应该将该变量从模型中剔除[36-37]。VIF值越大，变量与其他变量之间就具有越强的共线性。TOL小于0.2表示因子之间存在多重共线性，小于0.1，表明存在严重的共线性。多重共线性分析表明所选择的指标因子的容差均在0.4以上，高于阈值0.1，并且其VIF值均小于3，低于阈值5(表1)。综合分析表明本文所选的因子之间不存在多重共线性问题，验证了评价指标的合理性。

表1 影响因子的多重共线性分析

3.2 指标与滑坡的关系分析

本文分析了滑坡数量占比与不同指标之间的关系。将各个连续的指标变量等分成30个区间，分析不同区间内滑坡数量的分布情况。由于篇幅有限，本文在此处呈现芦山地震滑坡的结果图。芦山的地震滑坡主要集中在高程值较小的区间上(图1A)，表明此次地震事件中区域内滑坡主要分布在低海拔地区。随着坡度的增大，芦山的滑坡比率与坡度呈现出正相关关系(图1B)。对于平面曲率和剖面曲率，随着曲率的增加，滑坡比率增加(图1C—D)。在TWI和VRM较大时，滑坡占比相对更高(图1E—F)，而在TWI和VRM较小时，滑坡数量较多。断层距越大，滑坡占比越小(图1G)，大部分滑坡点集中在断层距较小的区间内。地震参数MMI与滑坡比率呈现出明显的正相关关系(图1H)，随着MMI的增大，滑坡占比整体升高。

图1 芦山地震，滑坡发生比率与各指标之间的关系

3.3 MCMC的模型参数估计

本文对马尔可夫链的中burnin参数设置为2 000，表示弃用前2 000次的抽样消除初始值的效益影响。mcmc参数设置为40 000，表示使用后面的40 000次抽样进行后验参数的估计。此外，采样过程的收敛性检验至关重要，只有迭代收敛之后才能够用来获得参数的后验分布。本文利用样本的轨迹图和样本间的自相关性进行收敛性诊断。图2展示了芦山地震滑坡案例中的MMI和坡度参数估计的轨迹图和样本自相关图。通过样本路径图可以看出，采样过程基本稳定。从样本的自相关图可以看出最终样本达到平稳状态，收敛性较好。综合考虑样本采样轨迹图和自相关图，可以认为参数样本的分布已收敛到参数的后验分布。

图2 芦山地震中马尔可夫链迭代轨迹图和自相关图

MCMC参数估计的结果见表2。在三次滑坡事件中，9个因子对于滑坡的发生都有显著的影响(p<0.05)。 MCMC方法得到的是后验分布的均值(Mean)作为参数的估计值。为了比较不同参数的不确定性，利用变异系数(Coefficient of variation， CV)进行分析，可以看出，在九寨沟地震事件中，岩性因子的参数不确定性最高，其次为VRM，而其他因子的CV值都较低。对于鲁甸地震事件，剖面曲率的不确定性较高，其次为岩性和平面曲率因子。在芦山地震案例中，所有因子的不确定性均较低，只有平面曲率因子的CV值绝对值大于1。各个指标因子参数在芦山事件中的不确定性最小，分析认为这可能与样本数量有关(芦山地震事件中的滑坡样本数量是其他两个事件总和的12倍之多)。岩性因子在九寨沟和鲁甸两个事件中的不确定性都相对较高，分析认为可能是岩性数据的分辨率相对较低造成的，而九寨沟地震和鲁甸地震的区域面积小，一定程度上导致了其参数的不确定性较高。

表2 三次地震滑坡事件中基于MCMC的Logistic回归模型中的参数估计

3.4 模型预测精度评价

在每个事件中，利用30%的未参与模型训练的测试数据集来验证模型的预测精度。结果表明Logistic回归模型在3个区域都具有较高的预测精度。本文所使用的滑坡数据中，芦山地震中的滑坡点最多，其次是鲁甸，九寨沟滑坡最少。九寨沟地震，鲁甸地震和芦山地震的测试总精度分别为0.875，0.928，0.927，三次滑坡事件中精确率均在0.9以上，召回率均在0.8以上，F1-score均大于0.85(表3)。从各个预测结果的混淆矩阵图也可以看出(图3)，预测正确的滑坡(TP)和非滑坡(TN)占比最大。此外，三次地震滑坡事件中的AUC值均大于0.9。综合分析表明基于MCMC的Logistic回归模型具有较好的预测精度。

表3 模型测试精度评价参数

注：A，D代表九寨沟地震；B，E代表鲁甸地震；C，F代表芦山地震。

3.5 滑坡敏感性评价结果

利用得到的Logistic回归模型对研究区的每个像元进行预测，得到滑坡敏感性评价图(图4)。将滑坡点叠加在敏感性图上进行对比，通过对比发现3个事件中的地震滑坡点都基本上落在0.9～1的高敏感性区间内，仅有少部分的滑坡点落在敏感性较低的区域。在芦山地区，由于滑坡点数量多，分布密集，造成敏感性高的区域在空间上成片分布。

3.6 因子相对重要性评价

本文利用随机森林算法分别对3个案例中的因子重要性进行分析，对比因子相对重要性图(图5)发现，地震参数MMI，高程和断层距在三次地震事件中重要性排序均较高。在九寨沟和芦山地区，MMI为最重要的因子，而在鲁甸地震滑坡中，高程为最重要的影响因子。岩性在芦山地区的重要性相对较大，而在九寨沟和鲁甸地区重要性相对低，分析认为芦山地震滑坡分布的空间范围较大，区域内岩性种类较多，岩性指标对模型预测的滑坡敏感性的空间分布有较大影响。而九寨沟地震和鲁甸地震滑坡分布的范围较小，区域内岩性种类较为单一，对滑坡敏感性的空间分布影响不明显。

4 结论

本文采用基于MCMC方法的Logistic回归模型研究地震滑坡敏感性，利用MCMC方法确定模型参数，并获得参数的后验分布，对参数的不确定性进行了评价，并基于所得到的模型进行了区域的滑坡敏感性空间制图。在芦山地震事件中，参数的不确定性均比较低，在九寨沟事件中岩性的不确定性较高，而在鲁甸地震案例中，岩性，剖面曲率和平面曲率的不确定性较高。整体而言，三次事件中的大部分因子的回归系数具有较低的不确定性。本文所构建模型具有较高的预测精度，在三次地震滑坡事件中AUC值均大于0.9，表明了MCMC方法在Logistic回归模型中进行参数估计的有效性。滑坡影响因子的相对重要性在不同的区域有差异，其中MMI，高程以及断层距在3个案例的相对重要性都是最大的。