基于三种多变量不安定指数分析模型的汶川县地质灾害易发性评价

2022-03-10毕鸿基周天伦

中国地质灾害与防治学报 2022年1期

毕鸿基，聂磊，曾超，范晨，周天伦，刘晓

（1. 中国地质大学(武汉) 教育部长江三峡库区地质灾害研究中心, 湖北武汉 430074；2. 中交第二公路勘察设计研究院有限公司, 湖北武汉 430056）

0 引言

我国是地质灾害发育最为严重的国家之一，特大、重大地质灾害( 如滑坡、泥石流、崩塌等)时有发生[1]。2008年5月12日，四川省汶川县发生8.0级特大地震，汶川县地质环境遭受到了毁灭性的破坏。根据多位学者对汶川震后地质灾害发生规律的研究分析[2−6]，震后汶川灾区的地质灾害将持续活跃20～25 a。2017 年“6·24”茂县特大滑坡灾害显示，特大地震对后期地质灾害的影响还将持续更长的时间[7]。根据四川省地质灾害隐患点掌上查询系统，汶川县目前存在地质灾害隐患点2 152处，主要为崩塌、滑坡、泥石流和不稳定斜坡，对居民的人身财产安全存在极大的威胁。

近年来，多变量不安定指数法在地质灾害易发性评价中受到重视，该方法最早用于边坡稳定性评价研究中。1992年简李滨[8]提出将研究区域划分为网格，以不安定指数法进行区域性边坡稳定评估。多变量不安定指数分析模型以定量的对各评价因子进行赋值，排除了主观因素的影响，适应于大范围、地质条件复杂且对区域内地质灾害状况有一定了解的地区。

有许多学者以不同模型对汶川县内多种地质灾害危险性进行了研究。如许冲等[9]在GIS支持下，基于层次分析法对汶川地震区的滑坡易发性进行了评价；孟祥瑞等[10]运用因子分析法，在GIS平台下对研究区的地质灾害易发性进行评价；林虹宇等[11]在Matlab软件及GIS平台下，采用可拓模糊物元模型并结合层次分析法，对岷江上游典型泥石流进行了易发性评价。

虽然已有学者应用不同的方法对汶川县各种地质灾害易发性做出了评价，但大部分文献所使用的评价方法在评价因子的赋值方面仍然带有很大的主观性。为了弥补这一不足，本文基于多变量不安定指数分析法模型，并加以改进，开展汶川县地质灾害易发性评价，将汶川县已有的地质灾害普查数据作为数据源，实现定量地对各影响因子进行赋值，提高评估结果的合理性。

1 多变量不安定指数分析法模型建立

1.1 多变量不安定指数分析法：标准模型

多变量不安定指数分析法是在计量统计的基础上，于一种相对关系下求出不安定指数（Dt）用以描述该研究区域发生地质灾害的危险性的高低。

在用多变量不安定指数法评价地质灾害易发性时，以诱发地质灾害的影响因子为基础建立评价指标，各个因子共同影响评价不安定指数Dt。Dt可作为判定地质灾害易发性的指标，Dt越大，则意味着该范围内发生地质灾害的可能性就越大。建立模型的具体步骤如下：

首先，确定各个致灾因子中各个等级的评分值：通过ArcGIS，利用各专题图层与已有的地质灾害分布图进行叠加分析，记录各专题地图中不同等级的地质灾害栅格数量，再求得其各自所占百分比，用式（1）计算各致灾因子中各等级的评分值：

式中：i——影响因子序号；

Xi——各专题地图中不同等级的地质灾害点栅格数量所占总的地质灾害栅格数量的百分比；

Xi_max——最大百分比；

Xi_min——最小百分比。

然后，计算各个致灾因子的变异系数，式（2）为变异系数(Vi)的计算公式：

式中：µi——式（1）中Xi的平均值；

σi——Xi的标准差。

再根据各个致灾因子的变异系数通过公式（3）计算其影响权重比例Wi：

式中：i——影响因子序号；

n——影响因子个数；

Vi——各个致灾因子的变异系数。

最后将权重值（Wi）与评分值（di）代入多变量不安定指数（Dt）式（4）。

式中：n——影响因子个数；

dn——各个致灾因子各个分级的评分值；

Wn——各个致灾因子的权重。

根据（4）式，通过GIS中的栅格计算器将各专题图层代入进行计算，即可得到研究区各区域的不安定指数Dt，用以表示区域内地质灾害易发性的大小。

1.2 多变量不安定指数分析法：衍生模型

在诸多前人研究中，传统不安定指数法被广泛应用的同时，也有学者提出不同的数学模型求得不安定指数。

1.2.1 线性累加模型

各级评价因子权重值（Wi）以线性累加取代幂次相乘的方式，求得不安定指数Dt。如式（5）所示：

式中：i——影响因子序号；

n——影响因子个数；

Wi——各个致灾因子的权重；

di——各个致灾因子各个分级的评分值。

1.2.2 幂次累加模型

在吸取前人学者的研究经验后，本文提出以幂次累加方式求得Dt的模式，如公式（6）所示：

式中：i——影响因子序号；

n——影响因子个数；

Wi——各个致灾因子的权重；

di——各个致灾因子各个分级的评分值。

1.2.3 不同模型评价的实施过程

在地质灾害易发性评价中，分别应用式（5）或式（6）替换标准算法的式（4），通过ArcGIS中的栅格计算器将各专题图层代入进行计算，即可得到研究区各区域的不安定指数Dt，用以表示区域内地质灾害易发性的大小。

在本研究中，将分别利用幂次相乘、线性累加及幂次相加三种模式求得Dt，然后对各个模式的准确性进行对比，选取准确性最高的模式对研究区地质灾害易发性进行评价。

2 研究区概况

汶川县位于四川省阿坝藏族羌族自治州东南部的岷江两岸，东西宽84 km，南北长105 km，总面积4 084 km2。东北部贯穿龙门山脉，西南部为邛崃山系，地势上西北高东南低，西部多分布海拔3 000 m以上的高山，其中最高的四姑娘山海拔为6 250 m：东南部的岷江出口处海拔仅有780 m。汶川县降雨较多且集中、干湿季明显，雨量528.7～1 332.2 mm，全县各地区降雨量差距不大[12]，岷江由县北部入境，贯穿东部，长达88 km。杂谷脑河、寿江、草坡河为汶川县境内岷江三大支流，由西向东注入岷江。沿岷江及各大支流河谷、半山地带为全县主要经济活动区。汶川县地层发育比较完整，但奥陶系、志留系有大部分缺失，仅小范围出露，岩性变化较大，以第四系松散地层及强风化岩浆岩为代表的软弱岩土体分布广泛。由于汶川县处于九顶山新华夏构造带，地质构造较为复杂，断裂、褶皱较为发育。汶川县境内共有三条主要断裂带（青川-茂汶断裂带、北川-映秀断裂带、江油-灌县断裂带），三条断裂带皆呈北东-南西方向斜穿全县[13−14]（图1）。汶川县境内潜在的地质灾害隐患包括崩塌、滑坡、泥石流和不稳定斜坡，共有地质灾害隐患点2152处（图2）。

图1 汶川县地质概况Fig.1 Geological survey of Wenchuan County

图2 汶川县地质灾害隐患点分布图Fig.2 Distribution of potential geological hazards in Wenchuan

3 研究区地质灾害易发性评价

3.1 致灾因子的选取及分级

控制和影响地质灾害危险性程度的因子有多种。研究区地质灾害类型主要为滑坡、崩塌、泥石流及不稳定斜坡，因此这里选择对几种地质灾害都适用的6种评价因子：坡度、坡向、地层岩性、距水系距离、距断层带距离、植被覆盖率。

（1）坡度。坡度是影响滑坡、崩塌、泥石流和不稳定斜坡的重要影响因素。一方面坡度较缓的斜坡，会造成松散固体物质的堆积，为滑坡、泥石流的形成提供物质条件；另一方面斜坡的角度还会影响斜坡内应力的分布，当坡度较大时，会为斜坡两侧崩塌体的能量转换提供条件。因此，坡度对地质灾害的发生有较大影响。利用ArcGIS的空间分析功能对汶川县30 m分辨率的数字高程模型（Digital Elevation Model，DEM）进行空间分析，并采用自然断点法，得到坡度分级图见图3（a）。

图3 各评价因子分级图Fig.3 Classification chart of evaluation factors

（2）坡向。由于研究区处于高山地区，地形起伏较大，坡向不同的区域，通过太阳辐射、温度等途径造成坡面上自然地理要素的差异。并且研究区大部分范围内海拔都较高，导致上述差异在一定程度上被放大，从而加重了对岩石风化程度以及破碎程度的影响，进而影响了地质灾害的发生及规模。利用ArcGIS的空间分析功能对汶川县30 m分辨率的DEM进行空间分析，并采用自然断点法，得到坡向分级图，见图3（b）。

（3）地层岩性。岩土是地质灾害形成的基础，所以地层岩性对地质灾害的发生有着极大的影响。在致密坚硬、结构完整的岩土发育区，很少有崩塌、滑坡、泥石流活动；相反，在结构破碎或松散软弱岩土发育区，崩塌、滑坡、泥石流等地质灾害活动强烈。根据汶川县的岩土类型，可分为7类：震旦系砂岩、白云岩；志留系板岩、千枚岩；泥盆系石灰岩；石炭系灰岩；二叠系闪长岩；三叠系砂岩、千枚岩；侏罗系砂岩。地层岩性分级见图3（c）。

（4）植被覆盖率。植被起到护坡和防止水土流失的作用，并且植被根系具有加固土体，提高土体抗剪强度的能力，嵌入基岩或下部老黄土的根系还起到锚筋作用。所以植被的覆盖一定程度上可以减轻各种地质灾害发生的风险。采取归一化植被指数（Normalized Difference Vegetation Index，NDVI）来表现植被的生长状况及覆盖率，植被归一化指数是指遥感影像中，近红外波段的反射值与红光波段的反射值之差比上两者之和，可以通过ArcGIS中的栅格计算器在研究区的遥感影像中获得。如图3（d）所示，其中负值表示地面覆盖为云、水、雪等，对可见光高反射；0表示有岩石或裸土等；正值，表示有植被覆盖，且随覆盖度增大而增大，这里用自然断点法分为5级。

（5）距断层距离。断层极易诱发地质灾害，断层是岩层或岩体顺破裂面发生明显位移的构造。在地貌上，大的断层常常形成裂谷和陡崖，沿断裂带附近岩层因长期地质构造运动破坏严重，易发生崩滑现象，松散物源堆积较多后还可能造成较大规模泥石流。汶川县境内共有三条主要断裂带 (青川-茂汶断裂带、北川-映秀断裂带、江油-灌县断裂带 )，三条断裂带皆呈北东—南西方向斜穿全县。利用ArcGIS中的缓冲区分析功能获得研究区距断层距离分级栅格图见图3（e）。

（6）距水系距离。河流水系的冲蚀切割作用会造成更多的临空面，也会对土体斜坡的密实度造成影响，并且冲散的岩土也为滑坡、泥石流的发生提供松散物源。所以，距水系的距离也对地质灾害的发生有一定的影响作用。利用ArcGIS中的缓冲区分析功能获得研究区距水系距离分级栅格图，见图3（f）。

3.2 确定各致灾因子分级评分值及变异系数

利用ArcGIS中的叠加分析功能，将汶川县地质灾害隐患点分布图与各致灾因子分级图进行叠加分析。叠加分析后的地图中各分级的栅格数量所占总栅格量的百分比就相当于不同等级的地质灾害点数量所占总的地质灾害数量的百分比。

（1）利用式（1）计算出各因子评分值见表1—6。

表1 坡度评分值Table 1 Grading of slope

表2 坡向评分值Table 2 Grading of slope aspect

表3 地层岩性评分值Table 3 Score value of stratum lithology

表4 植被覆盖率评分值Table 4 Score value of vegetation coverage

表5 距断裂带距离评分值Table 5 Score value of distance from fault zone

表6 距水系距离评分值Table 6 Score value of distance from water system

（2）利用式（2）计算出各致灾因子的变异系数见表7。

表7 各致灾因子变异系数Table 7 Variation coefficient of disaster-inducing factor

3.3 地质灾害易性区划

利用式（3）确定的各致灾因子的权重值见表8。

表8 各致灾因子权重值Table 8 Weight value of disaster-inducing factors

应用传统的多变量不安定指数分析法模型进行评价，在对各个评价因子的各分级赋值时，多为位数较多的小数，而在使用ArcGIS软件对各评价因子的各个分级进行赋值时，value值字段只接受整数，若采用四舍五入取整，则会显著影响评价模型的精度。10kdn，将10kdn赋予各专题图层。然后利用改进后的传统多变量不安定指数（Dt）计算式（7），及改进后的衍生多变量不安定指数（Dt）计算式（8）和式（9）。

对模型进行改进，从而提高评价结果的精确度，选择将各评价因子各级的评分值小数点保留k位有效数字(本研究中k=1)，在利用ArcGIS进行重分类时，将第3.2节获得的各影响因子的各级评分值乘以10k，记为

其中，k=1，di∈[1,10]。由各公式的数学关系可以得出幂次相乘及线性累加模式下,Dt∈[1,10]；幂次累加模式下，Dt>n，这里即Dt>6。通过ArcGIS的栅格计算器，分别获得三种模式下研究区各区域的不安定指数Dt值，并分别用自然断点法将Dt值分为5级，即将研究区地质灾害易发性分为五个等级：极低危险区、低危险区、中危险区、高危险区、极高危险区（图4）。

图4 各模式下研究区地质灾害易发性分级图Fig.4 Classification of geological hazard susceptibility in the study area

如果缺少科学的验证，一切预测模型及结果在一定程度上毫无科学依据[15]。文中利用接受者操作特性曲线（Receiver Operating Characteristic Curve，ROC）[16−18]分析法对本研究中三种模式下的研究区地质灾害易发性区划结果进行验证，选取准确性最高的模型作为评价该研究区地质灾害易发性的最终模型，以预测地质灾害面积百分比累加与实际地质灾害隐患点面积百分比累加分别对三种模型构建检验折线（图5）。

图5 评价结果检验曲线Fig.5 Verification curve of evaluation results

图中，幂次相乘模式下，ROC曲线下与坐标轴围成的面积（Area Under Curve,简称AUC值）为75.29；线性累加模式下，AUC值为67.03；幂次累加模式下，AUC值为74.87。对比显示，幂次相乘的模式下准确性最高，最适合作为该研究区地质灾害易发性评价模型。

以幂次相乘模式下对研究区地质灾害易发性的区划结果为准，对研究区地质灾害易发性进行评价分析：

（1）1≤Dt<3.84，为地质灾害极低易发性区域，面积为718.5 km2，占研究区总面积17.6%，主要集中在汶川县西部的耿达、卧龙两镇内，未来几乎不会发生地质灾害。该区域距断裂带较远，水系较少，无居民，多为高山，不需要特地开展地质灾害防治工作。

（2）3.84≤Dt<5.42，为地质灾害低易发性区域，面积为787.6 km2,占研究区总面积19.3%，主要集中在汶川县西部及北部，主要分布在耿达、卧龙、龙溪三镇。该区域距断裂带较远，居民较少有较低可能发生地质灾害，但需要提高当地居民防范意识，及抗灾减灾意识，一旦发生地质灾害，可把损失降到最低。

（3）5.42≤Dt<6.76，为地质灾害中易发性区域，面积为784.8 km2,占研究区总面积19.2%，主要集中在汶川县中部及北部，该区域发生地质灾害的可能性较小。该区域距离断裂带有一定距离，有较多条水系支流，未来可能发生小规模滑坡、泥石流，可能会对该区域居民造成一定财产损失，但可以通过及时监测和一些工程措施加以避免。

（4）6.76≤Dt<8.24,为地质灾害高易发性区，面积为1 004.3 km2占研究区总面积24.6%，主要集中在汶川县中部与东部，分布在岷江左岸及紫坪铺水库上游地区，该区域发生地质灾害可能性较高，未来可能有中等规模的地质灾害发生，该区内居民居住较密集，人类活动频繁，地质灾害产生的破坏性也会较大。因此，需要加大地质灾害防治力度。

（5）8.24≤Dt≤9.96，为地质灾害极高易发性区域，面积为788.8 km2，占研究区总面积19.3%，主要集中在汶川县中部，岷江流域两岸与映秀-北川断裂带两侧，未来有极大可能发生地质灾害，且由于该区域人类活动频繁，水系复杂，岩层破碎，地质灾害的规模会较大，且该区域内以往都发生过大规模的滑坡、泥石流等，存在许多不稳定斜坡。因此，该区域内应以以防为主、防治结合、全面规划、综合治理的原则，严格开展地质灾害防治工作，以保护居民人身财产安全。

4 讨论

（1）在本研究中的各致灾因子中，距断层带距离是影响地质灾害易发性的主导因素。研究区汶川县为地震频发区，被三条断裂带贯穿，而断裂带附近地质灾害极易发生，尤其是“5·12”地震后，震发断裂带附近地质构造遭受严重破坏，近十几年来地质灾害频发，可以看出断裂带为影响该区域地质灾害易发性重要因素，而表8亦证明如此，这也就体现了对评价因子定量赋予权重的优势。

（2）在赋予6个评价因子权重时，以现有的地质灾害隐患点为基础，用定量的方法得到每个评价因子的权重值，相较于以往传统的专家打分法，更客观、准确。在权重确定过程中，由于该区域地质灾害多且大小不一，很难具体统计出各评价因子某一分级中灾害点的数量，所以这里将灾害点转化为栅格数据，以栅格数量表示灾害点数量，通过各致灾因子分级图层与现有灾害点栅格图层叠加分析，清楚的得出致灾因子各分级中灾害点的比例，相比于同样以定量方法获取各致灾因子权重的信息量法不需要繁琐的数据处理步骤。综上，多变量不安定指数分析法是一种可以对评价因子定量赋值的模型，可避免主观因素的影响，且易于操作，适用于已知研究区内部分地质灾害信息的大范围内的地质灾害易发性评价。

（3）从各影响因子权重来看，距断层带距离及地层岩性权重相似，同时为影响研究区地质灾害易发性的主要因素。分析其原因，主要是由于“5·12”震发断裂带附近大量分布着二叠系闪长岩及志留系板岩，其岩石结构在地震时受到严重破坏，所以地层也会出现与距断层带距离因子相似的现象。

（4）选用三种不同的数学模型分别获得研究区地质灾害易发性分区图，并用ROC曲线分别求得三种模式下评价结果的AUC值，比较后，选出最适合该研究区的评价模型，即以传统不安定指数分析法中的幂次相乘来求得不安定指数Dt。显然，传统的多变量不安定指数分析法模型在该研究区地质灾害易发性评价结果中较衍生出来的两种模式更为精确。