胡现振, 付少杰, 迟宏庆, 张辉军, 张志飞
(河北省地质环境监测院,河北省地质资源环境监测与保护重点实验室,河北 石家庄 050021)
地质灾害风险评价是我国目前正在开展的一项重要的区域性地质灾害防治工作[1-2]。近年来,国内外学者和工程技术人员针对地质灾害风险评价开展了一系列研究和实践,取得了一定成果。Ohlmacher等[3]采用多元逻辑回归方法对美国堪萨斯州东北部的山体滑坡进行了危险性评价,并确定了坡度是影响滑坡危险性的重要变量; Pradhan等[4]采用逻辑回归模型和人工神经网络对马来西亚槟岛地区的滑坡危险性进行了区划研究; Abbaszadeh等[5]采用人工神经网络模型选取地质、水文等14个因子对瑞典西南部滑坡进行了易发性评价; 张春山等[6]采用经验模型法和统计分析模型法对渭滨区地质灾害进行了风险评价; 齐信、徐继维等[7-8]对地质灾害风险评价的定义、发展历程、评价方法等进行了总结探讨; 刘洋等[9]采用综合信息量模型对通江县地质灾害危险性进行了评价; 甄艳等[10]利用GIS技术与信息量模型相结合的方法对都江堰虹口乡进行了风险评价; 刘小青等[11]通过构建权重线性组合(weighted linear combination,WLC)算子方法及模糊综合评判法对北川县进行了地质灾害风险评价; 闫成林等[12]采用模糊贝叶斯网络模型对辽东半岛东部海岸进行了风险评价。相关研究工作有力地推动了风险评价工作的发展,但我国地质灾害风险评价工作仍处于实践和探索阶段,针对不同的地质环境条件,尚未形成合适的方法和标准,风险评价结果与地质灾害风险精准防控尚存在一定差距[2,7,13]。基于地区实际情况的评价方法的合理选择以及评价成果如何更好地服务于风险管理,仍是需要进一步深入探讨的问题。
武安市地处太行山东麓,区内地质环境复杂,矿业开发及工程建设活动强烈,地质灾害多发频发,给当地居民带来了严重的安全隐患和经济损失,在该区开展地质灾害风险评价,加强地质灾害防治显得尤为迫切。本文以河北省武安市为研究区,以1∶5 万风险调查数据为基础,基于ArcGIS技术平台,利用层次分析(analytic hierarchy process,AHP)法与信息量相耦合的方法,针对不同灾种选取不同的影响因素,对地质灾害的易发性、危险性、易损性及风险进行全面评价,旨在为武安市地质灾害风险管控提供依据,同时为其他地区结合自身实际选用合理的评价方法提供思路和参考。
武安市位于河北省邯郸市西北部,东南距邯郸市区约29 km,是邯郸市辖的唯一县级市。研究区北与邢台沙河市接壤,东与邯郸市区、永年相邻,南接磁县、峰峰矿区,西倚涉县、山西左权县,总面积为1 806 km2,共辖22个乡镇,铁路、公路网络密集,交通十分便利。武安市属大陆性季风气候,受季风影响,四季较为分明,年平均气温为12.4~13.8 ℃。1956—2020年年平均降水量为 536.6 mm,受地理位置及地形等因素的影响,多年平均降雨量西北多东南少。区内发育子牙河和漳卫南运河两大水系,主要有北洺河、南洺河等5条河流,地表水较为发达。
研究区地处太行山东麓,大地构造属太行隆起与华北沉降带的过渡带,三面环山,山地丘陵起伏,地势NW高SE低,西北部平均海拔在1 350 m左右,东北部平均海拔在190 m左右。出露地层岩性主要有长城系砂岩、寒武系灰岩页岩、奥陶系白云岩灰岩和第四系砂砾石及亚砂土层。NNE向展布的新华夏构造类型最为发育,断裂结构面多以高角度冲断层形式出现,既有张性特征,又有强烈的挤压特点,断层面显示有逆扭和顺扭性质的擦痕。
武安市地形起伏大,构造发育,地质环境条件复杂,切坡修路建房、矿业开发等工程活动强烈,自然因素叠加人为因素使得地质灾害较为发育。据武安市1∶5万地质灾害风险调查结果,现有地质灾害隐患点共计69处(图1)。崩塌、滑坡、泥石流地质灾害主要发育在西北构造侵蚀中低山区,地面塌陷及伴生地裂缝主要分布于煤、铁资源较为丰富的中部丘陵及盆地区域。地质灾害隐患共威胁314户1 182人,房屋3 269间,公路5 445 m,铁路1 000 m,矿山15座,农田163.5亩,威胁财产达6 455.8万元。
图1 研究区地质灾害分布Fig.1 Geological hazard distribution in the study area
地质灾害风险评价主要包括易发性评价、危险性评价、易损性评价和风险评价4个部分[13-14]。地质灾害易发性评价是以影响地质灾害形成的内在控制因素为出发点,考察静态下区域内某空间中发生灾害的可能性大小[14-15]; 地质灾害危险性评价是在易发性评价基础上考虑降雨、地震等诱发因素作用下进而发生特定规模灾害的可能性; 地质灾害易损性评价是评价承灾体遭受地质灾害破坏的严重程度; 地质灾害风险评价是危险性和易损性评价叠加的结果。可见,地质灾害易发性评价是风险评价的基础和关键。地质灾害易发性评价方法主要有两类: 定性分析法和定量分析法[15]。定性分析主要根据经验进行评级分析,如层次分析法; 定量分析则更偏重于利用数学方法进行量化研究,如信息量法。为使评价结果更加客观科学,本文选取层次分析法与信息量法相耦合的评价模型进行易发性评价,有效解决了定量与定性、可靠度与效率之间的矛盾,使评价结果更接近实际情况[9]。其评价模型如图2所示。
图2 地质灾害风险评价模型Fig.2 Geological hazard risk assessment model
AHP法是一种定性与定量相结合的多目标决策方法[16],通过对评价目标进行拆解分析,将各因素分解成目标层、准则层和指标层,同层次的要素受到上层要素的约束,也对其下层要素具有支配作用,从而形成由上至下的递阶层次(图3)。
图3 AHP 法基本流程Fig.3 Basic flow of AHP
首先,根据建立的层次模型,构建两两判断矩阵,对选取的评价因子的相对重要性进行打分(表1)。
表1 判断矩阵的标度和含义Tab.1 Scale and meaning of judgment matrix
然后,采用和积法计算最大特征根与其对应的特征向量。在此基础上,进行一致性检验。公式为
(1)
式中:CR为随机一致性比率;CI为一致性指标;RI为同阶平均随机一致性指标(表2)。
表2 平均随机一致性指标值Tab.2 Average random consistency index values
当CR<0.1时,说明判断矩阵具有较好的一致性,可以通过检验; 反之,则要重新进行打分,直到通过检验。
地质灾害是在多种因素共同影响下形成的,信息量可反映一定地质环境下最容易导致灾害发生的因素与其细分后各区间的组合。通过对某个评价单元内在某种影响指标作用下灾害发生频数与区域内灾害发生频数进行比较来实现。每个评价单元均受较多因素的综合影响,各因素状态组合条件下地质灾害发生的信息量由公式(2)确定,其式为
(2)
式中:I为某个单元内地质灾害发生的总信息量,可作为地质灾害易发性指数,用来表示灾害发生可能性大小;Ni为某个因素i状态下地质灾害点数或面积,km2;Si为某个因素i状态下的分布面积,km2;N为评价区地质灾害总点数或总面积,km2;S为指评价区的总面积,km2。
评价时采用信息量法计算出各因素所占信息量,利用层次分析法计算出各因素所占权重,将两者进行结合,加权叠加得到总信息量I,计算公式为
I=∑Wi×Iij。
(3)
式中:Wi为层次分析法中计算得到的权重;Iij为第i个评价因子第j分级的信息量。
将计算得到的信息量值赋给对应的栅格图层,再利用ArcGIS栅格计算功能对各影响因子进行叠加,完成评价因子综合信息量的计算并生成叠加图层,根据自然断点法进行重分类,从而得到地质灾害易发性分级图。在易发性分区结果基础上,叠加降雨量图层与其进行栅格叠加运算,进一步形成危险性评价结果。将危险性评价和易损性评价结果采用矩阵分析法进行叠加,最终得到地质灾害风险评价结果。
地质灾害发育程度往往受到较多因素的控制和影响,区域不同灾害的严重程度则表现出明显的差异性和复杂性,所以在对研究区进行风险评价时首要任务就是确定评价单元。目前,国内外学者对评价单元的划分方法大致分为行政单元、自然斜坡或地貌单元、规则栅格单元等类型。根据各评价单元的特点并结合武安市实际情况,本文以高精度DEM数据为基础,选取规则栅格单元作为评价单元,分辨率为25 m×25 m,研究区划分栅格单元共计2 889 806个。
地质灾害作为一个庞大而复杂的复合非线性系统,不同的灾害类型受不同影响因子的制约。为确保易发性评价的准确性,同时结合武安市地质灾害孕灾特点,本研究对滑坡、崩塌、泥石流、地面塌陷、地裂缝分别选取了不同的评价因子,并利用ArcGIS进行了相关性分析,确保了评价因子选取的合理性。
崩塌、滑坡选取高程、坡向、归一化植被指数(normalized difference vegetation index, NDVI)、距水系距离、距断层距离、工程地质岩组、地貌类型、斜坡形态8项为评价因子; 泥石流选取流域面积、主沟纵坡降、沟壑密度、melton比率、水土流失侵蚀模数5项为评价因子; 地面塌陷、地裂缝则是根据矿山资料及形成状况采用定性方法进行了评价。
通过信息量计算模型,将各评价因子图层进行分级,计算得到各分级因子的面积及区域内分布的隐患点数量,再利用ArcGIS重分类功能,得出各指标的信息量,然后根据AHP原理计算得到各评价因子权重(表3,表4)。
表3 崩塌、滑坡地质灾害易发性评价因子信息量及权重Tab.3 Information content and weight for the susceptibility assessment factors of collapse and landslide
表4 泥石流地质灾害易发性评价因子信息量及权重Tab.4 Information content and weight for susceptibility assessment factors of debris flow
根据上述评价指标的信息量值进行评价因子的加权叠加计算,采用自然间断法将叠加计算的值分为非易发、低易发、中易发、高易发4个等级。将上述得到的各灾种评价结果归一化处理,再利用ArcGIS的叠加分析功能将其叠加,得到武安市地质灾害综合易发性评价结果(图4)。
图4 地质灾害易发性评价Fig.4 Geological hazard susceptibicity assessment
由图4可知: 地质灾害高易发区面积为168.48 km2,占研究区总面积的9.33%,主要分布于活水乡、管陶乡等中山区域; 中易发区面积为457.33 km2,占研究区总面积的25.32%,主要分布于贺进镇、冶陶镇、马家庄乡等低山区及地形坡度较大区域; 低易发区面积为393.74 km2,占研究区总面积的21.80%,主要分布于马家庄乡、矿山镇等丘陵平缓地带; 非易发区面积为768.45 km2,占研究区总面积的43.55%,主要分布于盆地及河谷平原区域。
地质灾害危险性是在易发性评价的基础上,考虑地质灾害的外在诱发因素,进一步刻画和预测地质灾害的影响范围和发生概率。评价结果见图5。
图5 地质灾害危险性评价Fig.5 Geological hazard danger assessment
以往地质灾害相关资料显示,降雨是研究区地质灾害发生的主要诱发因素,因此本次以降雨为诱发因子进行危险性评价。从降水量的空间分布来看,受地理位置及地形等因素的影响,多年平均降雨量东南少,降雨量在500 mm左右,中部、东部中等,降雨量在520~560 mm 之间,西北部降雨量相对较大,在560 mm以上。选取武安市各雨量站近10 a月累计平均降雨量,绘制了降雨量等值线图,并将其分为8个等级: <500 mm、(500,520] mm、(520,540] mm、(540,560] mm、(560,580] mm、(580,600] mm、(600,620] mm、>620 mm。不同降雨量对地质灾害的影响是不同的,利用降雨图层与易发性进行栅格叠加运算,采用自然间断法叠加计算形成危险性评价初步计算结果,结合实际调查结论,再进行平噪处理,得到危险性评价结果。
由图5可知: 高危险区面积为172.60 km2,占研究区总面积的9.56%; 中危险区面积为863.61 km2,占研究区总面积的47.82%; 低危险区面积为769.79 km2,占研究区总面积的42.62%。
易损性是地质灾害社会属性的反映,区域社会经济发展的程度也是其重要的影响因素。目前易损性评价尚没有统一的评价指标,笔者根据研究区地质灾害的分布规律以及承灾体所处的自然社会环境选择评价指标。本次在充分研究武安市地质灾害及其承灾体基本特征以及全面考虑研究区易损性评价数据的基础上,选择了人口、建筑物与道路3个影响因素,构建了适合研究区的易损性评价指标体系。并在此基础上,根据不同评价指标,结合研究区实际情况,构建各自的分级标准和赋值标准。所需数据来源于《武安市统计年鉴(2020)》[17]和第三次全国国土调查数据。
人口易损性是通过计算各区域人口密度得到人口密度分布图,并将其归一化处理,按照自然断点分级法将人口易损性分为3级从而得到人员易损性评价分区图层; 建筑物密度通过计算建筑物平面面积总和与区域面积的比值求得,采用对建筑物所占面积进行归一化的处理方法,取归一化值作为评价区内的基础易损性,利用ArcGIS软件分区处理得到建筑物易损性评价分区图层; 道路密度反映道路的密集程度,为区域内道路总长度与区域面积的比值,利用ArcGIS将道路密度归一化处理。根据交通设施易损性赋值原则,对不同设施类型和等级的交通设施进行了易损性赋值(表5),并叠加易损度栅格,形成交通设施易损性图层。在获得人口、建筑物和道路易损性后,运用层次分析法确定各评价因子权重(表6),采用ArcGIS栅格计算器分析功能进行叠加(综合易损指数=人口分布指数×0.54+建筑设施指数×0.30+交通设施指数×0.16),并按照自然断点法将易损性分为低易损、中易损、高易损3个等级,最终得到综合易损性评价图(图6)。
表5 交通设施易损性赋值Tab.5 Vulnerability assignment for traffic facilities
表6 易损性评价因子权重Tab.6 Weight value for vulnerability assessment factor
图6 地质灾害易损性评价Fig.6 Geological hazard vulnerability assessment
地质灾害风险评价为危险性评价和易损性评价结果的叠加。叠加运算采用矩阵分析法,首先对危险性和易损性等级由高到低分别赋值3、2、1,然后在ArcGIS中运用栅格计算器模糊叠加功能进行运算,最后根据《地质灾害风险调查评价技术要求(1∶50 000)》[18]风险等级划分矩阵(表7)得到风险栅格评价结果。
表7 地质灾害风险等级划分Tab.7 Classification of geological hazard risk degree
风险评价的目的在于服务地质灾害风险管控及国土空间规划等。然而,以栅格单元呈现的风险评价结果多数范围较小且较分散,在实际应用中存在诸多不便。因此,本文从管理应用角度出发并结合研究区实际情况,将以栅格单元呈现的评价结果按照临近性及相似性原则归并于斜坡单元(泥石流为整个流域),对不同风险等级的区域进了调整优化,最终得到研究区风险评价结果(图7)。
图7 地质灾害风险评价Fig.7 Geological hazard risk assessment
由图7可知: 高风险区主要分布于活水乡、管陶乡北部,面积为101.25 km2,占研究区总面积的5.61%; 中风险区主要分布于管陶乡、马家庄乡和矿山镇,面积为796.48 km2,占研究区总面积的44.10%; 低风险区主要分布于邑城镇、大同镇等地势平缓地带,面积为908.27 km2,占研究区总面积的50.29%。
经复核: 高风险区共有地质灾害隐患点30处,其中崩塌23处、滑坡4处、泥石流3处,威胁人口582人,威胁财产1 858.8万元; 中风险区共有地质隐患点39处,其中滑坡1处、崩塌7处、泥石流1处、地面塌陷27处、地裂缝3处,威胁人口600人,威胁财产4 597.0万元; 低风险区无地质灾害隐患点。地质灾害风险评价结果与地质灾害隐患发育情况较为一致。
(1)本文采用AHP-信息量耦合模型对研究区进行了风险评价,不仅考虑了每个因子所提供的信息量值和所占权重,而且将主观化的量值和客观化的量值进行结合。同时,根据研究区地质环境条件和地质灾害发育特征,针对不同的灾害类型选取了不同的评价因子,并从管理角度出发,将评价结果由栅格单元归并于斜坡单元,取得了较好的评价效果。
(2)武安市地质灾害高风险区主要集中在西北部侵蚀构造中低山区,面积为101.25 km2,占比5.61%; 中风险区主要分布于西北部低山区和中部丘陵区,面积为796.48 km2,占比44.10%; 低风险区主要分布于邑城镇、大同镇等地势平缓地带,面积为908.27 km2,占比50.29%。风险区分布与灾害发育情况基本吻合。
(3)将栅格单元归并于斜坡单元有利于风险区管理,但斜坡单元的划分精度会直接影响风险区的划分精度,因此,采用此方法时,应尽可能合理地划分斜坡(沟谷)单元。
(4)以栅格为单元进行风险评价,便于数据采集和空间分析,但在某种程度上会破坏地质单元的整体性,与地质环境条件的真实情况存在一定差距。今后可进一步探索斜坡单元和流域尺度的精细化风险评价,更好地指导地质灾害精准防控。