磨英飞,潘宏坚,邓荫万
(广西壮族自治区地质环境监测站,广西 南宁 530201)
广西北流市是典型的桂东南花岗岩丘陵区,滑坡、崩塌地质灾害多发,开展地质灾害易发分区,在地质灾害预警和防治规划、国土空间规划中具有重要应用。目前,地质灾害易发性评价方法有综合指数法、证据权法、信息量法、多元统计分析法、逻辑回归法等,综合指数法、信息量法、多元统计分析法等在选取评价因子和权重中存在较大主观性和人为因素,不同人员选取评价因子和赋予的权重可能会有差别,同时评价模型可靠度直接取决于地质灾害原始数据的精度和发育数量,存在评价因子组合状态比较多、样本需求大、实际统计数量受限等缺点,具有一定的局限性。证据权法是采用贝叶斯条件概率原理、综合各种证据层来预测某种事件发生概率的一种定量方法,在地质灾害易发性评价中,充分考虑了地质灾害与评价因子之间的相关关系以及评价因子之间的相互关系,评价结果较其他方法具有客观性和可靠性。国内学者和工作人员应用证据权法[1-6]开展地质灾害易发性评价,均取得了较好的效果。如胡燕等[7]采用证据权法开展巴东县城滑坡灾害易发性评价,将巴东县城划分为滑坡极高易发区、高易发区、中易发区、低易发区与极低易发区5类;罗可[8]应用证据权法进行大埔县崩塌滑坡地质灾害易发性评价分区,98.37%的灾害点位于中、高易发区内。本文采用证据权法对广西北流市南部石窝镇、六靖镇、清湾镇[9]等地区进行地质灾害易发性评价和分区。
研究区位于广西东南的玉林市北流市南部,面积397.9 km2,多年平均降雨量约2 200 mm,降雨多集中在4—9月,占年均降雨量80%。河流发育有5条小河流。
研究区属构造剥蚀丘陵地貌,分为高丘、低丘2种类型。
(1)高丘区山顶高程+300.0~+616.8 m,自然坡度一般30°~40°,局部45°~50°,沟谷切割深200~300 m,多呈“V”形,分布面积约192.6 km2。
(2)低丘区山顶高程+100.0~+300.0 m,自然坡度一般20°~30°,局部35°~40°,沟谷切割深50~150 m,多呈“U”形,分布面积205.1 km2。
植被以速生桉树为主,覆盖率较高,村庄多分布在低丘区的丘陵坡脚。
研究区地层岩性有侵入岩和变质岩(图1)。
侵入岩侵入时期为寒武纪—志留纪、白垩纪,岩性以黑云(钾长)二长花岗岩为主,少量黑云花岗片麻岩、片麻状黑云二长花岗岩。
变质岩包括寒武系浅变质岩及晚元古代平政组花岗变粒岩。其中,寒武系浅变质岩岩性为千枚岩、片岩、板岩,由含砾长石石英砂岩、不等粒砂岩、粉砂岩、页岩浅变质形成;平政组花岗变粒岩岩性为黑云变粒岩、黑云二长变粒岩、黑云斜长变粒岩及少量黑云石英片岩、含硅线堇青黑云片岩等。
花岗岩、花岗变粒岩区丘陵山坡分布第四系坡残积粉质黏土、砂质黏土,厚2~3 m;下部为全风化层砂土,以砾砂、粗砂为主,局部含黏性土、中砂、粉砂,呈散体状,遇水易软化崩解,厚度10~30 m。浅变质岩区丘陵山坡第四系坡残积黏土、粉质黏土,厚度一般在0.5~3.0 m,局部厚5.0~8.0 m。
研究区构造以南南西—北北东向为主,主要发育有陆川米场—石窝断裂、大坡外—六靖断裂、上珍断裂、中龙断裂、侯山断裂、六仑断裂、黄田断裂7条次级断裂。
研究区破坏地质环境条件的人类工程活动主要有城镇工程建设、居民临坡建房、公路建设形成高切坡、高填方等。集镇及村屯建设挖、填方形成高5~15 m、坡度50°~80°的高陡人工边坡,道路建设高挖低填形成高5~20 m、坡度45°~70°的高陡边坡,多数人工边坡无支护、排水等防护措施,强降雨作用下易引发崩塌、滑坡地质灾害。经调查发现,研究区有553处高陡人工边坡。
研究区发育地质灾害点155处。其中,崩塌30处,滑坡17处,不稳定斜坡108处。崩塌滑坡物质以残坡积土体和全风化花岗岩为主,花岗岩、花岗变粒岩区共发育地质灾害106处,占灾害点总数的68.4%,浅变质岩区共发育地质灾害49处,占灾害点总数的31.6%。地质灾害以小型为主。滑坡体厚1~5 m,崩塌体厚2~3 m。地质灾害主要发育在居民屋后人工边坡和道路两侧的人工边坡。
证据权法是一种条件独立假设的前提下,基于贝叶斯条件概率的定量预测方法。证据权法最初是作为医疗诊断支持的方法,20世纪80年代末,加拿大数学地质学家Bonham-Carter G F等[10-15]将该方法引入到矿产资源定量预测与评价。该方法作为一种人工智能模型,已被国内外学者广泛用于多元信息综合评价,近年来又被引入到地质灾害的评价中。其基本原理为:假设在一特定的研究区内,已知有n种二值证据图层,1代表因子对灾害发生的证据存在,0代表不存在,通过证据权模型给出该二值化的证据因子图层的权重,最终叠加多元图层,实现地质灾害易发性评价。证据权法的分析流程如下。
图1 研究区地质简图Fig.1 Geological sketch map of the study area
2.1.1 权重计算
(1)
式中,W+为证据因子存在区的权重值;W-为证据因子不存在区的权重值,其大小表示证据因子与地质灾害发生和不发生的关系密切程度。证据因子和灾点正相关表示为W+>0,W-<0,负相关为W+<0,W->0,不相关时权重为0。P为不同条件发生的条件概率。
已知地质灾害点的先验概率为Po=D/T。证据因子权重由落入特定证据因子图层的灾点数和全部灾点数之比与证据因子图层面积和调查区总面积之比的比值决定。为表示证据因子对于地质灾害发生和不发生的区分能力,可计算相对系数C=W+-W-,用来度量证据因子和地质灾害之间的相关性大小,C值越大越有利于地质灾害的出现,C值越小越不利于地质灾害的发生。
2.1.2 证据综合
在权重值计算及分析的基础上,通过证据层的优选,选择权重较大、与地质灾害关系密切的证据层,剔除权重较小、与地质灾害关系不密切的证据层;进一步进行证据因子相对灾点的条件独立性检验,剔除地质灾害权重相对较小而与其他证据因子相关性大的证据层。对最终筛选出的n个关于地质灾害点条件独立的证据因子,根据贝叶斯法则,研究区任一单元K发生地质灾害的可能性,即对数后验概率见式(2)。
(2)
式中,O为D的概率,O(D)=D/(T-D);D为存在地质灾害的单元网格数;Bi为第i个证据层;K(i)为在第i个证据因子层存在时为+,不存在时为-;Wi为第i个证据因子存在或不存在的权重。
最后计算后验概率:
(3)
后验概率值的大小表示易发性的高低,其值在0~1。后验概率值越大,表示易发性越高;后验概率值越小,表示易发性越低。
不同地域地质环境条件和地质灾害发育特征不同,选取的地质灾害易发评价因子[16-19]存在差异。研究区地貌类型单一,岩组较单一,植被茂盛且类型较单一,地质灾害主要形成于第四系残坡积土体和全风化花岗岩中。根据地质灾害发育特征和形成机理,选择地形地貌、地层岩性、斜坡结构类型、土体厚度、构造、人类工程活动等证据因子进行地质灾害易发性评价,各证据因子细分不同的证据层。
2.2.1 地形地貌
地形地貌选取坡度、坡面曲率等证据因子,将坡度划分为 ≤10°、(10°,15°]、(15°,20°]、(20°,25°]、(25°,35°]、(35°,45°]、(45°,55°]、>55° 共8个证据层;将坡面曲率划分为≤-7、(-7,-3]、(-3,3]、(3,7]、>7共5个证据层。
2.2.2 地层岩性
研究区地质灾害主要发育在丘陵斜坡的第四系残坡积粉质黏土和全风化花岗岩中,沟谷及河岸分布的第四系全新统冲洪积层地质灾害不发育,各地层分布面积和地质灾害发育情况如图2所示。易崩易滑地层为晚志留世扶新单元(S3F)、六杨单元(S3Ly)侵入黑云(钾长)二长花岗岩、寒武纪古桑单元(∈G)、平政组(Pt3p)黑云变粒岩、片岩、寒武系黄洞口组第一段(∈h1)浅变质长石石英砂岩与页岩,共划分为11个证据层。
图2 研究区各地层岩性分布面积及发育灾害数量柱状图Fig.2 Histogram of lithologic distribution area and development disaster of each stratum in the study area
2.2.3 斜坡结构类型
研究区的浅变质岩由沉积岩变质形成,继承了沉积岩的岩层层理,按斜坡坡向与岩层层理倾向间的夹角,划分为顺向坡(夹角≤30°)、同向斜交坡(30°,60°]、横向坡(60°,120°]、反向斜交坡(120°,150°]、逆向坡(150°,180°]共5个证据层,侵入岩及平政组花岗变粒岩按顺向坡考虑。
2.2.4 土层厚度
土层厚度考虑斜坡残坡积层和全风化层的厚度,划分为≤ 3 m、(3,6 m]、(6,9 m]、(9,12 m]、>12 m共5个证据层。
2.2.5 构造
包括与断裂的距离和岩体节理裂隙等结构面对地质灾害的相关性。
与断裂的距离划分为[0,200 m]、(200,500 m]、(500,1 000 m]、>1 000 m共4个证据层。结构面考虑不同倾向的结构面对地质灾害发育的控制或影响,划分为[0°,30°]、(30°,60°]、(60°,90°]、(90°,120°]、(120°,150°]、(150°,180°]、(180°,210°]、(210°,240°]、(240°,270°]、(270°,300°]、(300°,330°]、(330°,360°] 共12个证据层。研究区结构面倾向在120°~130°、170°~190°、225°~245°、265°~280°等方位区段上有比较明显的优势,与120°~140°、170°~210°、225°~270°等方向的地质灾害发育数量较多基本吻合(图3),在0°~100°、300°~360°方位区段上地质灾害分布零散。
图3 研究区结构面倾向与地质灾害分布统计Fig.3 Structural inclination and distribution statistics of geological hazards in the study area
2.2.6 人类工程活动
人类工程活动考虑丘陵坡脚房屋建设强度、切坡强度,由于研究区内道路沿线未发育地质灾害点,同时切坡强度也考虑了道路的人工边坡分布,因此不考虑交通建设强度指标。
丘陵坡脚房屋建设强度可以采用房屋面积密度表示,选用最新的1∶1万测绘信息和遥感影像资料,统计房屋面积密度,利用自然间断点法分为4个等级,划分为[0,23 780]、(23 780,74 312]、(74 312,157 540]、>157 540共4个证据层。
切坡强度采用高陡斜坡和地质灾害的分布点密度表示,研究区调查了553处高陡斜坡和155处地质灾害点,大部分地质灾害点都是建房切坡引发的,高陡斜坡在强降雨等因素作用下可能失稳发生崩塌滑坡,因此高陡斜坡和地质灾害的分布点密度划分为[0,11.41]、(11.41,28.29]、(28.29,55.67]、>55.67共4个证据层。
坡度、坡面曲率、斜坡结构类型、土层厚度、与断裂的距离、结构面、丘陵坡脚房屋建设强度、切坡强度等均为连续型数据,地层岩性为分类型数据。对于连续型证据因子数据,采用计算累计证据权重判断拐点的方法进行分级。先按照较小的级距将连续性证据因子分级,然后在ArcGIS平台下初步计算W+、W-累计权重和相对系数C,初步判断累计权重值曲线的拐点作为分级的断点值,将证据因子划分为若干级,重新计算各个因子等级的权重,再绘制累计权重曲线,判断曲线拐点和分级的合理性,并进行多次迭代,直至分类趋于合理;分类型数据可直接分级。最终各评价因子的证据层权重及后验概率计算结果见表1—表9。
表1 坡度证据层权重及后验概率计算结果Tab.1 Calculation results of slope evidence layer weight and posterior probability
表2 坡面曲率证据层权重及后验概率计算结果Tab.2 Calculation results of slope curvature evidence layer weight and a posteriori probability
证据权法的前提是各证据因子相互之间应满足条件独立性假设。为确保评价指标的相互独立性,需对评价指标进行相关性分析,此次采用卡方检验方法进行检验。首先,将各因子图层二值化,即将相对系数C<0和C=0对应的栅格赋值为0,将C>0的栅格单元赋值为1,在GIS软件内将地质灾害点与二值模式栅格图层叠加分析;然后,用地质灾害点图层汇总二值的数量,计算期望值和卡方值(二联表计算过程简略),各证据因子相互的卡方值见表10。
表3 地层岩性证据层权重及后验概率计算结果Tab.3 Stratum lithology evidence layer weight and posterior probability calculation results
表4 斜坡结构证据层权重及后验概率计算结果Tab.4 Calculation results of evidence layer weight and posterior probability of slope structure
表5 土层厚度证据层权重及后验概率计算结果Tab.5 Calculation results of evidence layer weight and posterior probability of soil layer thickness
表6 与断裂的距离证据层权重及后验概率计算结果Tab.6 Calculation results of weight and posterior probability of evidence layer of distance from fault
根据统计学概率,当P(K2≥K)=0.001时,查表得卡方临界值为10.828,即当卡方值小于10.828时(可信度99.9%),两因子之间不存在相关性,此时2组证据因子层可以同时参与易发性指数计算。
表7 结构面倾向证据层权重及后验概率计算结果Tab.7 Calculation results of evidence layer weight and posterior probability of structural plane tendency
表8 房屋建设强度证据层权重及后验概率计算结果Tab.8 Calculation results of evidence layer weight and posterior probability of housing construction intensity
表9 切坡强度证据层权重及后验概率计算结果Tab.9 Calculation results of evidence layer weight and posterior probability of slope cutting strength
根据表10计算结果,地形坡度、坡面曲率、地层岩性、土层厚度、切坡强度5个证据因子相互之间的卡方值小于临界值10.828,满足各证据因子相互之间独立性假设,满足基于贝叶斯条件概率的证据权计算方法。研究区地质灾害主要发育于第四系残坡积土体和全风化花岗岩、花岗变粒岩土体中,斜坡结构类型、结构面对地质灾害的控制和影响仅属个别现象,断裂对地质灾害的影响微弱。因此,选择地形坡度、坡面曲率、地层岩性、土层厚度、切坡强度5个证据因子组合进行易发性评价,也基本符合研究区的地质灾害发育现状。
研究区易发性评价最终选取地形坡度、坡面曲率、地层岩性、土层厚度、切坡强度等5个证据因子,采用证据权法计算各证据因子的权重和后验概率,利用ArcGIS提取各证据层所在斜坡单元的后验概率值,按斜坡单元对各证据因子的后验概率值进行叠加计算,采用自然间断法,将研究区划分为地质灾害极高易发区、高易发区、中易发区、低易发区4个等级(图4),其中极高易发区面积62.16 km2、高易发区135.22 km2、中易发区171.23 km2、低易发区29.29 km2。
表10 各证据因子相互卡方值计算结果Tab.10 Calculation results of mutual chi-square values of various evidence factors
图4 北流市南部易发性评价分区Fig.4 Zoning map of vulnerability assessmentin the south of Beiliu City
根据易发性分区结果,极高易发区分布有69处地质灾害点,高易发区分布有45处,中易发区分布有23处,低易发区分布有18处,极高、高易发区分布的灾害点占73.5%。通过绘制ROC曲线(图5),计算AUC值为 0.718,说明证据权法在北流市南部易发性评价的可信度较高。
采用证据权法开展北流市南部地质灾害易发性评价,根据北流市南部地质灾害发育和分布特征,选择合理的证据因子,通过检测各证据因子独立性,选择地形坡度、坡面曲率、地层岩性、土层厚度、切坡强度5个证据因子组合作为评价指标,划分出地质灾害极高、高、中、低易发区。经验证,73.5 %的地质灾害点落在极高、高易发区中,划分结果具有较高的可信度。
图5 易发性分区与地质灾害数量ROC曲线Fig.5 Prone zoning and ROC curve of geological disaster quantity