基于聚类分析法的地质灾害风险评价<br/>——以韩城市为例

基于聚类分析法的地质灾害风险评价
——以韩城市为例

2021-09-26李冠宇苏生瑞董永超

科学技术与工程 2021年25期

李冠宇，李鹏，郭敏，2，苏生瑞*，董永超

(1.长安大学地质工程与测绘学院，西安 710061；2.西安市地质环境监测站，西安 710007；3.陕西省水工环地质调查中心，西安 710068)

中国领土幅员辽阔，地形复杂多样，约65%为山地丘陵区，近年来，由于人类工程活动规模与强度不断增加，造成崩塌、滑坡、泥石流、地面塌陷等地质灾害频发[1]。地质灾害防治通报显示，2019年产生地质灾害数目、地质灾害造成失踪和死亡人数及带来的直接经济损失较前五年相比均有所减少[2]。尽管中国在防灾减灾上已取得了显著成效，但仍然面临着较为严峻的防灾减灾工作[3]。因此开展地质灾害风险评价和建立健全地质灾害综合防治体系迫在眉睫，识别和预测地质灾害发生尤为重要。

随着地理信息系统(geographic information system，GIS)技术、数学学科及信息全球化的发展，中外利用GIS对地质灾害进行风险性评价手段和方法逐渐成熟。由传统的专家打分法、模糊数学法、层次分析法、信息量法等评价模型，逐渐转变为结合数理统计和机器学习，支持向量机(support vector machine，SVM)、神经网络、聚类分析等多种评价模型相结合的评价方法[4-9]。像模糊数学法、指数法等传统的评价方法，需要采用专家打分法或层次分析法等方法确定各评价因子的权重，主观因素对评价结果的干扰较大[10]。而聚类分析法不采用专家打分法、层次分析法对评价因子进行权重计算，避免了主观因素的影响，其评价结果更为客观可靠。

最早的聚类分析概念由Tryon于1939年建立，早期应用于生物学和生态学研究[11]。因聚类分析具有计算方便、结果直观、不受其他约束条件的限制等优越性，21世纪以来已被广泛应用于地学与环境中大量存在的各种分类工作[12]。张海燕等[13]以吉林省灾害监测数据为基础，采用阈回归联合聚类分析方法对研究区进行地质灾害易损性评价，该方法对数据具有较高的拟合精度；王滨等[14]以重庆市北碚区为例，采用因子一聚类分析方法，对研究区泥石流进行危险性区划与实际情况吻合度较高；王欣[15]运用不确定性层次聚类分析方法对延安市宝塔区滑坡进行危险性评价，其精度较传统层次分析方法高；刘宏[16]采用模糊聚类分析对深圳地区场地建设进行适宜性评价，将场地分为适宜性差区、基本适宜区和适宜区。

在地质灾害风险评价中要考虑众多因素，其中主观干扰因素较多，由于K平均聚类分析方法在处理大数据集时，相对容易实现高效收敛，而且操作简单、分类结果较好，可以在一定程度上规避人为主观因素的影响。现以韩城市地质灾害风险评价项目为依托，利用GIS平台和SPSS统计分析相结合的手段，基于K平均聚类分析和系统聚类分析方法，以韩城市为例进行地质灾害风险评价，为韩城市防灾减灾提供一定的理论依据。

1 风险评价模型构建和聚类分析方法

1.1 风险评价概述

地质灾害风险性评价主要由易损性、危险性和风险性评价3部分组成[12]。现今普遍认同的灾害风险概念是由美国滑坡灾害专家Varnes提出的，即因地质灾害造成不好结果的可能性，包括地质灾害发生的可能性和其导致后果(损失)两个方面[17]。地质灾害风险性评价的基础是查清灾害发展历史及趋势、变化规律和形成条件，其主要是对灾害发展程度和其危害能力进行分析和判断。危险性评价是考虑地质灾害的活动程度和破坏能力，以及地质灾害发生概率的大小。易损性评价的主要对象是承灾体，是为了研究在已有经济技术条件下人类社会对地质灾害的抗御能力，得到不同社会经济要素的易损性参数。

1.1.1 评价单元选取

根据研究区的实际情况选取合适的评价单元[18]。规则栅格单元、自然斜坡和地貌单元和行政单元划分是划分单元格的主要方法。由于自然斜坡和地貌单元法受主观影响较大和行政单元划分对区域地质条件考虑欠缺的缺点，一般采取规则栅格单元法对研究区进行地质灾害危险性评价。而地质灾害易损性评价主要考虑的是经济因素，对区域地质条件考虑较少，一般采取行政单元划分法或规则栅格单元法对研究区进行评价。栅格单元划分的大小影响评价结果的精度，划分单元越小，评价结果精度越高，然得到栅格数据较多，后续处理任务较大，对电脑要求较高，反而不利于研究工作进行[18-19]。

1.1.2 评价因子的选取

根据研究区地质灾害分布规律和孕育条件，可将研究区致灾因子分为历史因子、基础地质因子和诱发因子三类。历史因子主要包括灾点密度；基础因子包括高程、坡度、地层岩性、地质构造和水系；诱发因子包括降雨量和人类活动，其中将人类活动分为修建公路和采矿两种。9个二级评价因子共同构成了研究区易发性评价指标体系。在易发性评价的基础上，考虑地质隐患点的规模大小以及稳定性两个影响因素，与易发性评价因子共同构成研究区危险性评价指标体系。易损性因子分为人口易损性、物质经济易损性和环境资源易损性，其中人口易损性为人口密度；物质经济易损性为公路交通密度和国内生产总值(gross domestic product，GDP)密度；环境资源易损性为耕地密度和矿产资源密度，5个因子构成地质灾害易损性评价体系。如表1所示。

1.2 K平均聚类分析

聚类分析就是根据多个指标进行数学分类的一种多元统计方法。聚类分析主要分为三类：①系统聚类主要是对小样本的样品间聚类以及对变量聚类；②动态聚类适用于样本量较大时样品间的聚类，常用K-平均聚类分析方法处理；③有序样品聚类是对有排列次序的样本的样品间聚类，要求必须是次序相邻的样品才能聚成一类。根据分类对象的不同，聚类分析可分为研究变量之间关系的R型和研究样品之间关系的Q型；根据数学原理和分类方法的不同，又将聚类分析分为动态聚类、模糊聚类、图论聚类、系统聚类和灰色聚类[11]。

表1 地质灾害风险性评价指标Table 1 Geological hazard risk evaluation index

采用动态分析中的K平均聚类分析方法进行危险性分析，其优点是：对于大数据集，K均值聚类分析方法相对是可伸缩和高效的，容易实现，收敛速度快；当结果簇是密集的，而簇与簇之间区别显著时，它的分类结果较好；主要调控的参数是簇数k。尽管无监督K平均聚类分析方法经常以局部最优结束，但一般情况下达到局部最优就已经可以满足聚类的需求[21-22]。具体步骤如下。

步骤一确定需要计算的数据集X。

步骤二确定k个聚类中心向量C1，C2，…，Cj，…，Ck。

步骤三对数据集中每一个数xi，计算其与聚类中心的距离(计算距离的方式有：欧式距离、曼哈顿距离、切比雪夫距离等，其中欧式距离应用较为广泛。)，若xi满足：

(1)

则xi∈Cj。

步骤四计算各聚类中心新的向量值C′j(1

(2)

式(2)中：m为数据属性个数；Nj为第j个聚类域sj中所包含的样本个数。

步骤五若聚类中心不再变动，迭代公式如式(3)所示，则终止计算；否则转到步骤三。

(3)

式(3)中：N=1,2,…,k；sj为中心为Cj的聚类域。

1.3 系统聚类分析

系统聚类分析是将所有单个样品(或变量)看作是多维变量(或样品)空间中矢量的一个端点，经由给定的公式运算，得出各“点”之间的相似性统计量，组成对称的相似性矩阵，在该矩阵的基础上进一步寻找各“点”之间，或“点”群之间的相似程度，依据相似程度大小，把样品(或变量)归类成群；分类结果由二维谱系图表示出来[12]。具体计算步骤如下。

步骤一合理选取变量和样品，构成原始数据矩阵，表达式为

X=[xij]n×m

(4)

式(4)中：xij为第i个样品的第j个指标的观测数据。

步骤二将原始数据矩阵进行数据变换处理，数据变换处理方法有标准化变换、对数化变换、规格化变换，根据实际需要选择变换处理方法，一般采取标准化变换处理。标准化变换公式为

郑大一附院形成了一种常态：出院时，由神经内科医护人员开出包括鼻饲进食及留置胃管的注意事项等出院处方，随后，由电话随访小组进行随访。

(5)

步骤三数据标准变换后，构建相关系数矩阵R=[rjk]m×m，其中相关系数rjk计算公式为

(6)

步骤四对相关系数矩阵R=[rjk]m×m采取一次聚类法或逐步聚类法，直到所有的变量和样品聚类完成为止。

一次聚类法：从相似性矩阵中选出最大的一对变量或样品，形成第1级分类，选出次大的一对变量或样品，形成第2级分类，继续操作，直到所有的变量或样品聚集成一类为止。

步骤五经与实际情况对比，选取合适的分类方案，整理联结表，绘制谱系图，并对分类结果进行地质解释。

2 研究区概况

韩城市位于陕西关中平原的东北隅，关中盆地与陕北黄土高原的过渡地带，北靠宜川，西临黄龙，南接合阳，东隔黄河与山西省乡宁、河津、万荣县相望。地理坐标东经110°07′19″～110°37′24″，北纬35°18′50″～35°52′08″，总面积1 621 km2。研究区属暖温带半干旱大陆性季风气候，年平均气温13.5 ℃以上，区内年平均降雨量为538.1 mm，但降水分布很不均，多集中于7、8、9月份。研究区水资源丰富，澽水、芝水、盘水、凿开河、白矾河等9条主要河流纵贯东西，汇入黄河。研究区地势西北高、东南低，自然地貌单元主要有中山区、低山丘陵区、黄土台塬区和河谷阶地区。研究区是陕西省的主要煤炭开发基地，采煤历史悠久，在开采煤矿时，常常形成一定规模的采空区，引发滑坡、崩塌、地面塌陷等地质灾害频发。根据韩城市地质灾害风险调查评价，研究区共有地质灾害隐患点175处，其中：崩塌85处(48.57%)，滑坡49处(28.00%)，地面塌陷34处(19.43%)，其他地质灾害6处(4.00%)。研究区地质灾害呈不均匀分布，在东北部的低山丘陵区及东南部的黄土台塬区发育较多，在同一镇办地质灾害发育情况也具有差异性，集中发育在桑树坪镇、芝川镇、西庄镇、板桥镇以及金城办，如图1所示。

图1 韩城市地质灾害隐患点分布图Fig.1 Distribution map of hidden danger points of geological disasters in Hancheng City

3 基于K平均聚类分析的危险性评价

3.1 危险评价体系

由于自然斜坡和地貌单元法受主观影响较大且行政单元划分对区域地质条件考虑欠缺的缺点，再结合研究区实际情况，采取规则栅格单元法，即50 m×50 m矩形进行等面积划分，共划分618 844个栅格，根据表1选取11个危险性二级评价因子共同构成了韩城市危险性评价指标体系，进行研究区的危险性评价。对于定性的因子，将研究区评价因子划分为：极低、低、中、高危险区4个等级[5，19-20]。

3.1.1 灾点密度

研究区灾害点沿南北方向呈带状分布，主要分布在低山丘陵区和黄土台源区，大都分布在河流附近，且受人类活动影响强烈，具有聚集性的特点，使用灾点密度分析实现灾点密度提取。

3.1.2 高程

通过研究区数字高程模型(digital elevation model，DEM)提取高程，发现随着高程的增大，灾点分布减少，主要分布在380～900 m，在1 100 m以上无灾害点分布。反而高程越低，人类工程活动越强烈，对环境破坏程度越大，采空区越多，造成地质灾害隐患点数量越多。

3.1.3 坡度

坡度是地质灾害发生的重要控制因素，坡度不同地质灾害发生的概率也不尽相同。通过研究区DEM进行表面分析获取研究区坡度，发现地质灾害分布在坡度小于40°以下。坡度越陡，受重力影响越大，斜坡越不稳定，但受人类工程活动相对较小。

3.1.4 地层岩性

地层岩性是地质灾害稳定性的重要因素，研究区共分布有四种不同软硬程度的岩土。根据岩性的力学性质，极低、低、中、高危险区分别为块状坚硬岩、层状中等坚硬岩、层状较坚硬岩、松散软弱岩。

3.1.5 构造

研究区主要包括8条断裂构造和两条褶皱构造，构造主要以南北向为主。地质构造的分布使周围岩土体发生错断和破碎，工程地质性质较差，以距构造距离对其进行缓冲区分析。极低、低、中、高危险区构造距离分别为>1 500、1 000～1 500、500～1 000、<500 m。

3.1.6 水系

研究区水资源丰富，有72处灾害点沿河流分布，因此对水系做距离分析。极低、低、中、高危险区水系距离分别为>900、600～900、300～600、<300 m。

3.1.7 降雨

降雨是诱发灾害发生的重要因素，降雨量从东南到西北逐渐减小，以降雨在空间上分布规律为依据分为6个区域，分别为>650、600～650、550～600、500～550、450～500、<450 mm。

3.1.8 公路

人类活动范围往往是沿着道路分布，以乡级及以上级别的道路做缓冲区分析。极低、低、中、高危险区水系距离分别为>3 000、2 000～3 000、1 000～2 000、<1 000 m。

3.1.9 采矿

采矿活动是研究区重要的人类活动，也是诱发地质灾害发生的重要因素，以矿产开采情况进行危险性分级。极低、低、中、高危险区分别为未开采区、闭坑煤矿区、计划开采区、开采区。

3.2 评价因子归一化

由于选取的评价因子中有定性和定量两种因子，且各因子间的单位不能统一，无法对地质灾害进行统一评价，因此需要对定性的危险性评价指标根据其对危险性的影响大小进行量化。评价因子定量后，在ArcGIS中利用栅格计算器根据式(7)实现数据的阈值法归一化处理，结果如图2所示。

图2 韩城市易发性评价指标归一化图Fig.2 Normalized map of Hancheng City susceptibility evaluation index

(7)

式(7)中：xi为样本数据；xmax为样本最大值；xmin为样本最小值；yi为经过归一化处理后的样本数据。

3.3 基于SPSS的K平均聚类分析

利用ArcGIS将评价区域的618 844个栅格转化为点，将9个因子归一化处理后的数据多值提取至点。提取的数据导入SPSS，为消除数据间的差异，先进行标准化处理，并对各评价因子进行K平均聚类分析，变量选取为标准化后的9个因子，聚类数为4，聚集中心如表2所示，经过18次迭代最终将618 844个栅格分为4类，如表3所示。

表2 聚集中心Table 2 Clustering center

表3 每个聚集中的观察值数目Table 3 Number of observations in each cluster

将SPSS统计分析的分类结果导入ArcGIS中，将点转化为栅格，根据各评价因子距离各聚集中心的距离大小，将韩城市地质灾害分为极低、低、中、高易发区4个等级，得到地质灾害易发性区划图，如图3所示。

3.4 危险性评价及区划

危险性反映的是地质灾害的活动程度和破坏能力，以及地质灾害发生概率的大小。在易发性区划的基础上，考虑地质隐患点的规模大小以及稳定性两个影响因素，对灾点做缓冲区分析，并根据灾点稳定性和破坏能力的大小分别对其赋值，并在ArcGIS中与易发性区划进行叠加运算，将计算结果通过自然断点法重分类后，结合实际情况与前人研究概况，对计算结果进行修饰、合并，将韩城市地质灾害危险性分为极低、低、中、高危险区4个等级，并根据危险性区划图得出各个等级危险区的面积，结果如表4、图4所示。

危险性高的区域主要分布在低山丘陵和黄土台源区；高危险区到极低危险区灾害点占比和灾点密度呈递减趋势，由此可见，基于K平均聚类分析方法得到的危险性区划效果良好，具有较高的准确度。

图3 韩城市地质灾害易发性区划图Fig.3 Geological susceptibility zoning map of Hancheng City

表4 地质灾害危险性区划统计表Table 4 Statistical table of geological disaster risk zoning

图4 韩城市地质灾害危险性区划图Fig.4 Geological hazard risk zoning map of Hancheng City

4 基于系统聚类的易损性评价

4.1 易损性评价指标体系

承灾体应是受地质灾害威胁的所有人与事物，但由于一些因子无法获取和量化，此外一些因子对整体性评价的影响很小。根据表1选取5个易损性二级评价因子共同构成了韩城市易损性评价指标体系。以韩城市8个镇办作为评价单元，通过对各评价指标统计、分析、计算得出各评价单元的评价因子数据如表5所示。

表5 韩城市各镇办易损性评价因子数据Table 5 Vulnerability assessment factor data of each town office in Hancheng City

4.2 基于SPSS的系统聚类分析

根据表5构建易损性评价因子原始数据矩阵X，在SPSS中采用标准变化对数据进行预处理后，进行系统聚类分析。对图5系统聚类结果，并结合实际情况分析，可以将韩城市易损性划分为极低、低、中、高4个区域，并得出各等级易损性分区面积，如表6和图6所示。其中高易损区的镇办人口密度大、乡镇化程度高、道路便利、经济较发达、人类工程活动强烈，易损性区划结果与实际情况相符。

图5 系统聚类分析结果图Fig.5 Hierarchical clustering analysis result graph

表6 地质灾害易损性区划统计表Table 6 Statistical table of geological disaster vulnerability zoning

5 地质灾害风险性评价

地质灾害风险性评价采用矩阵分析法，对地质灾害危险性和易损性区划结果利用ARCGIS栅格计算工具进行叠加运算，结果如图7所示，并对风险分区结果进行统计如表7所示。地质灾害风险等级划分建议如表8所示。从风险性区划图可以看出，风险性具有西低东高的特点，与危险性和易损性相符，且风险性分区和野外地质调查结果能够较好地相匹配。