基于突变理论的湟水河流域崩滑易发性评价

李 彬，周金喜，吴钶桥，马建全，祁汉文

（1.青海省青藏高原北部地质过程与矿产资源重点实验室，青海 西宁 810000；2.青海省有色第一地质勘查院（青海省隐伏矿勘查重点实验室），青海 西宁 810000；3.青海省有色第四地质勘查院，青海 西宁 810000；4.西安科技大学地质与环境学院，陕西 西安 710000）

0 引言

随着社会的不断发展，人类工程活动不断加强，地质灾害带来的不良影响越来越多，如何科学准确的评判地质灾害的危险性十分重要。崩滑灾害的发生极其复杂，受地形地貌、岩性组合、地质构造、降雨、地震和人类活动等多种因素的影响[1]。2014年宋盛渊等[2]采用突变级数理论的方法对岩桑树水电站库区发育的潜在单体滑坡进行危险性评价，取得了较好的效果；尚志海等[3]通过对梅州市地质灾害伤亡人数、直接经济损失和发生频率的统计，引入突变理论的突变级数法，建立了地质灾害灾度的评估指标体系和等级划分标准，并对梅州市地质灾害灾度进行了评估；2020年，刘晓宇等[4]利用突变理论选取岩性、风化程度、坡度等7 个指标，对西部山区威胁公路的10 个典型滑坡进行危险性评价，评价结果与现场调查结果具有较高的吻合度。目前，突变理论在灾害评价中的应用，仅限于单体灾害危险程度的研究或是某一地区灾情程度的研究上，对于区域易发性评价研究中尚未见相关报导[5]。本文旨在创建一种新的基于突变理论的崩滑灾害易发性评价方法，开发一种基于突变级数理论的崩滑灾害易发性评价平台。区别于前人研究成果的是对整个研究区域栅格单元进行突变级数赋值，进而对整个研究区崩滑灾害易发性进行评价，而不是单单对单体灾害的评价。

1 突变级数理论基本原理及易发性评价平台开发

1.1 突变级数法基本原理

将研究区内的评价因子按照倒立的树枝状结构进行多层次分解是突变级数法的基本思想[6]。按照“互补”或“非互补”原则，将最底层的控灾因子进行无量纲处理（归一化），将归一化结果代入相应的突变模型中，并计算出该层的突变级数。然后逐层向上计算各层的突变级数，并根据最高层的突变级数将所研究的系统分级[7 -8]。

（1）建立逐层结构模型

建立逐层结构模型是突变级数法的关键，要充分考虑各层因子之间的关系，逐层分解，选择相对应的突变函数，按先主后次排列成倒立的树枝状结构。常用的分解形式如图1 所示。图中，x为状态变量；u、v、w、t为控制变量。

图1 突变模型系统Fig.1 Mutation model system

（2）评价指标的无量纲化处理

各参评因子的取值范围和单位量纲均不相同，无法对各参评因子进行比较。因此，就有必要将各参评因子转变为0～1 之间的无量纲数据。值得注意的是，所有参评因子数据均应该转化为“越大越好”型数据，对“越小越好”型数据可采用倒数法进行标准化处理。对于正指标，采用式（1）进行标准化；对于负指标，采用式（2）进行标准化[2]。

式中：xi——评价指标值；

xmax、xmin——各评价指标的最大值和最小值。

（3）突变级数的计算

各参评因子无量纲处理后就可以选择对应的初等突变类型（表1）。

表1 初等突变类型Table 1 Primary mutation types

根据不同模型的归一化公式和最底层因子的无量纲数据，可计算出对应因子的突变级数值。在计算时，无相关性因子取最小值，有相关性因子取平均值。逐层向上计算突变级数，最终求出总突变级数。突变级数法流程见图2。

图2 突变级数法流程图Fig.2 Flow chart of catastrophe theory procedure

1.2 基于突变级数理论的易发性评价平台开发

开发基于突变级数理论模型的易发性评价平台采用的是MATLAB 语言，它是当今科研领域最常用的应用软件之一，具有强大的符号运算、矩阵计算及数据处理等功能[9]。MATLAB 有包罗万象的工具箱，“草稿纸式”的编程语言，易学易用，用户不仅可以调用其内部函数进行计算，还可以根据自己的算法进行扩展编程[9]。对突变级数理论模型进行编程使其程序化，最终开发出基于突变级数理论模型的崩滑灾害易发性评价平台，其平台开发流程图（图3）。

图3 基于突变级数理论的易发性平台开发流程图Fig.3 Development flow chart of susceptibility platform based on catastrophe series theory

2 研究区概况及地质灾害影响因素

2.1 研究区概况

湟水流域位于青藏高原与黄土高原的过渡地带，总体为高原地貌。地势起伏较大，西北高、东南低，亦有南北高、中间低的特点（图4）。最高海拔位于西北部的湟水干流源区，约4 350 m，最低海拔位于东南部的湟水干流河谷，约1 650 m。根据流域内实际，又可细分为中高山区、丘陵区和河谷区。由于湟水流域地形地貌以沟壑纵横的丘陵为主，气候气象具有干旱少雨但多集中于汛期的特点，出露岩性以黄土、泥质岩等软弱岩土体多见，植被覆盖较为稀疏且森林覆盖率低，断裂、褶皱等构造发育，因而湟水流域地质环境脆弱，水土流失严重，崩塌、滑坡等地质灾害多发易发；加之，流域内人口密集、人类工程活动强烈，对地质环境的改造明显和突出，人为因素导致的崩塌、滑坡等也越来越多。因此，湟水流域崩滑灾害多发，特别是中下游地区的丘陵区是青海省乃至全国地质灾害高易发区。

图4 湟水河流域崩滑灾害分布图Fig.4 Distribution map of landslide disasters in the Huangshui River Basin

2.2 因子选取及建立逐层机构模型

本文参考国内外研究学者经常选取的崩滑灾害影响因子，同时考虑所选评价因子容易量化程度[10]，最终选择具有代表意义的12 个因子作为崩滑灾害易发性评价指标。选取因子分别为：工程地质岩组、距离城镇、平面曲率、高程、剖面曲率、距离水系、地形地貌、距离公路、坡度、距离断层、地形湿度指数（TWI）、多年平均降水量。

将12 个控灾因子作为最底层因子指标进行多层次分解，先主后次排列成倒立的树枝状结构。逐层向下分解，直至分解到将灾害点危险度作为目标层，以地质因素和环境因素两个变量为准则层，地貌、地形、工程地质岩组、地质构造、水文气象、人类工程活动6 个变量作为中间层，高程、坡度、平面曲率、剖面曲率、工程地质岩组、地形湿度指数（TWI）、距离断层、年均降水量、地貌类型、距离水系、距离公路、距离城镇12 个控灾因子作为最底层变量，建立逐层结构模型。其中，目标层灾害点危险度即灾害点突变级数，各最底层变量即为控灾因子，各层之间根据变量的不同选用不同的突变模型。表2 为湟水河流域崩滑灾害危险度评价逐层结构模型。

表2 湟水河流域崩滑灾害危险度评价逐层结构模型Table 2 Layer-by-layer structure model for landslide hazard assessment in the Huangshui River Basin

3 研究区崩滑灾害评价

3.1 评价指标无量纲化

由于最底层评价指标之间的取值范围和单位量纲均不相同，各指标之间无法进行比较，需将其转变为0～1 之间的无量纲数据。利用归一化原则，将湟水河流域内1 171 处崩滑灾害进行无量纲处理。表3 为评价因子无量纲处理结果表。

表3 评价因子无量纲处理结果表Table 3 Dimensionless processing results of evaluation factors

3.2 突变级数的计算

在对流域内1 171 处崩滑灾害点的各类因子进行无量纲化后，利用突变模型的突变势函数计算各层评价指标的突变级数。按照表2 的突变模型，其中间层、准则层及目标层的突变级数计算公式为：

中间层：

根据式（3）—（11）计算得到湟水河流域崩滑灾害点的突变级数，即各灾害点危险度（表4）。

表4 湟水河流域崩滑灾害点危险度计算结果表Table 4 Calculation results of landslide hazard at various points in the Huangshui River Basin

3.3 崩滑灾害危险性评价

对于崩滑灾害危险性等级的划分，尚未形成统一的标准。因此，选用ArcGIS 重分类自然间断点法将研究区内1 171 处崩滑灾害点危险度进行分级，结合研究区崩滑灾害易发性实际特征，将崩滑灾害危险性划分为轻度危险、中度危险、重度危险和极度危险4 个等级，崩滑灾害危险性评价指标划分标准见表5。

表5 崩滑灾害危险性评价指标划分标准Table 5 Division standard of landslide hazard evaluation index

对表4 中研究区1 171 处崩滑灾害点突变级数按照划分标准表5 进行分类，其评价结果如表6 所示。

表6 崩滑灾害危险性评价结果Table 6 Risk assessment results of rockfall disasters

根据表6 危险性评价结果可知，研究区内1 171 处崩滑灾害中轻度危险点75 处、中度危险点201 处、重度危险点476 处、极度危险点419 处。根据评价结果，研究区内灾害点以重度和极度危险为主，与实际灾情状况一致。

按照突变级数法可以快速评价单体灾害危险性，其优点是操作简单，评价结果可靠，但其缺点是只能进行单体灾害危险性分析，不能进行全流域易发性评价。为了对湟水河流域全域易发性进行评价，必须开发一套基于突变级数理论的易发性平台，既满足单体灾害快速且准确率高的优点，也要能满足全流域易发性评价的要求。因此，利用MATLAB 编程，对突变级数理论模型进行升级，在单体灾害危险性评价基础上开发一套基于突变级数理论模型的易发性评价平台。

4 突变级数理论易发性平台开发

4.1 平台设计流程

（1）栅格单元赋值

通过突变级数计算获得了湟水河流域内1 171 处崩滑灾害的危险度，将崩滑灾害各指标中间层、准则层及目标层突变级数值导入ArcGIS，利用ArcGIS 栅格转出工具将数据转出为ASCII 格式。湟水河流域总面积为16 100.65 km2，按照30 m×30 m 栅格单元划分，可以将研究区划分为17 889 612 个栅格单元，利用MATLAB编程，对每个栅格单元进行突变级数赋值。该方法的优点是将整个流域离散化，具有比常规方法精度高、操作简单、评价结果与实际情况符合程度高等特点。具体步骤如下：

①构建评价因子栅格图。（加载研究区评价底图）

②对所有评价因子进行归一化处理，得到最底层D（运用式（1）（2）对选取的高程、坡度等12 个控灾因子进行归一化处理，使其可以在同一维度内进行对比，得到12 个控灾因子的归一化结果D）。

③根据论文中表2 中逐层结构模型及突变模型，最底层D的各个层进行栅格计算并叠加得到中间层C（选取不同的突变模型，运用式（3）-（8）对D进行栅格计算，得到地貌、地形、工程地质岩组、地质构造、水文气象、人类工程活动6 个中间层的突变级数C）。

④根据论文中表2 中逐层结构模型及突变模型，中间层C的各个层进行栅格计算并叠加得到中间层B（选取蝴蝶型和尖点型突变模型，运用式（9）（10），对C进行栅格叠加计算，得到地质和环境2 个准则层的突变级数B）。

⑤根据论文中表2 中逐层结构模型及突变模型，准则层B的各个层进行栅格计算并叠加得到目标层A，以此得到整个流域的灾害危险度突变级数（选取尖点型突变模型，运用式（11）对B进行栅格叠加计算，得到目标层灾害点危险度的突变级数A，即得到整个流域的灾害栅格单元突变级数）。

⑥根据论文中表5，将计算得到的突变级数栅格图进行重分类，得到易发性分区图5（基于表5 给出的灾害危险性评价指标划分标准，对流域内灾害危险度进行重分类并绘图，得到流域灾害易发性分区图，见图5）。

图5 湟水河流域崩滑灾害易发分区图Fig.5 Landslide susceptibility zoning map in the Huangshui River Basin

（2）预定义平台界面及属性

新建一个平台、定义平台的结构、对加载数据定义到函数句柄、定义输出及平台的调用句柄、定义及更新平台的调用句柄，输出界面见图6。

图6 预定义平台界面及属性Fig.6 Predefined platform interface and properties

（3）读取数据并绘图

（以“D1地貌”按钮为例）读取数据、剔除异常值、保存数据的文本格式、构建句柄定义的坐标区窗口、展示图像、设置颜色、保存数据、将更新后的数据和文本数据加载到句柄中、保存句柄到指定对象中、弹出确认提示，输出界面见图7。

图7 读取数据并绘图Fig.7 Data retrieval and plotting

（4）数据归一化并绘图

（以D1地貌为例“归一化”按钮）恢复对象中的句柄数据、定义“地貌归一化”值、构建句柄定义的坐标区、展示图像、设置颜色、设置透明度、设置图名、保存数据指定对象句柄中。设置背景色，选择对应突变函数，对所有因子进行归一化处理，对D1—D12重复进行步骤,输出界面见图8。

图8 湟水河流域栅格单元突变级数计算图Fig.8 Calculation of catastrophe series for grid units in the Huangshui River Basin

（5）突变级数计算并绘图

读取各因子归一化数据、分层计算输出目标值A、构建坐标区、目标值A归一化、展示A归一化值的图像、设置颜色、设置透明度、保存数据到指定对象句柄。

（6）计算结果输出

调用对象中保存得句柄信息、创建文件“YF.txt”、将dm 中的文本数据转换为字符串数组、将HC 中所有行数据输入YF.txt、将Ag 所有行数据输入YF.txt、保存结束，输出结果。

4.2 湟水河流域危险性评价结果

在地质灾害评价中地质灾害易发分区是一项重要内容，地质灾害易发分区是根据灾害的易发程度对研究区进行划分等级[11-13]。根据平台计算结果，将“YF.txt”文件导入ArcGIS，读取栅格单元突变级数，按照自然间断分级法将研究区栅格单元突变级数分为五个等级，即湟水河流域危险指数。分别是极低、低、中等、高和极高，得到基于突变级数理论的湟水河流域崩滑灾害易发分区图和易发分区统计表。图5 为湟水河流域崩滑灾害易发分区图。

根据表7 统计结果可知，湟水河流域崩滑灾害的极高易发区、高易发区主要分布于湟水河南北两岸以及一级支流两侧的低山丘陵区，面积分别为4 661.79 km2和5 511.44 km2，占整个研究区总面积的28.95%和34.23%；中易发区分布于湟水河南北两岸以及各大支沟两侧的构造侵蚀中山区，面积占研究区总面积的23.99 %；低易发区、极低易发区主要分布在湟水河河谷平原地带及全流域中高山区，面积占比分别占总面积的11.80 %和0.90%。根据评价结果，湟水河流域以中高易发区为主，占研究区总面积的63.18%（图9）。

表7 突变理论模型危险性分区统计表Table 7 Statistical table of hazard zoning using catastrophe theory model

图9 易发性评价因子归一化图Fig.9 Normalized evaluation factors for susceptibility assessment

4.3 评价结果精度验证

ROC 曲线为二分类模型的预测精度提供了定量化的研究方法[14-15]。评价过程简单，评价结果非常直观，以至于能仅凭肉眼进行预测效果初步判断（图10）。ROC 以曲线下面积判断预测精度好坏（AUC）。对于任何预测试验中，AUC值在0.5～1。特征曲线的基本原理是在崩滑危险性区划结果中随机选取了等量的崩滑点和随机点，来绘制ROC 曲线，并确定曲线下面积AUC值[16-17]。其评判标准如表8 所示。

表8 AUC 值评判依据Table 8 Criteria for AUC value evaluation

图10 基于突变理论的湟水河流域易发性模型ROC 精度验证结果Fig.10 ROC accuracy validation results for the susceptibility model in the Huangshui River Basin based on catastrophe theory

湟水河流域崩滑灾害分布特征与基于突变级数理论的易发性评价结果基本一致。为了验证这一结果，在ArcGIS 随机抽取30 %的灾害点并创建同等数量的随机点，通过SPSS 软件中的分析工具ROC 曲线，获得基于突变级数理论的湟水河流域易发性评价ROC 曲线。根据结果，基于突变级数理论的湟水河流域易发性模型AUC值为0.843，根据评判标准，表明该预测模型精度很好，这与流域内崩滑灾害灾情状况基本相一致，说明该模型满足评价要求。

5 结论

（1）本文以地质因素和环境因素两个变量为准则层，地貌、地形、工程地质岩组、地质构造、水文气象、人类工程活动 6 个变量作为中间层，高程、坡度、平面曲率、剖面曲率、工程地质岩组、地形湿度指数（TWI）、距离断层、年均降水量、地貌类型、距离水系、距离公路、距离城镇 12 个控灾因子作为最底层变量，建立了基于突变理论湟水河流域崩滑灾害危险度评价逐层结构模型。

（2）采用ArcGIS 重分类自然间断点法将研究区内1 171 处崩滑灾害点危险度计算结果划分为轻 度危险、中度危险、重度危险和极度危险4 个等级，结果显示，研究区内1 171 处崩滑灾害中轻度危险点75 处、中度危险点201 处、重度危险点476 处、极度危险点419处。研究区内灾害点以重度和极度危险为主，与实际灾情状况一致。

（3）开发了一种新的崩滑灾害易发性评价平台——突变级数理论平台对湟水河流域崩滑灾害易发性进行评价。结果表明，湟水河流域崩滑灾害的极高易发区、高易发区主要分布于湟水河南北两岸以及一级支流两侧的低山丘陵区，面积分别为4 661.79 km2和5 511.44 km2，占整个研究区总面积的 28.95% 和 34.23%；中易发区分布于湟水河南北两岸以及各大支沟两侧的构造侵蚀中山区，面积占研究区总面积的 23.99 %；低易发区、极低易发区主要分布在湟水河河谷平原地带及全流域中高山区，面积占比分别占总面积的 11.80 % 和 0.90%。

（4）通过ROC 曲线验证，基于突变级数理论的湟水河流域易发性模型ACU值为 0.843，表明该预测模型精度很好，其优点是操作简单，评价结果可靠，且弥补了突变级数只能进行单体灾害危险性分析的短板，为崩滑灾害防治提供理论依据。