北京山区泥石流的单沟预报
2023-02-28陈文鸿柳清文师春香
陈文鸿，余 斌，柳清文，路 璐，马 超，孙 帅，师春香
(1.成都理工大学 地质灾害防治与地质环境保护国家重点实验室，成都 610059； 2.福建省地质工程勘察院 自然资源部丘陵山地地质灾害防治重点实验室，福州 350002； 3.北京市地质矿产勘查院 北京市地质灾害防治研究所，北京 100120； 4.北京林业大学 水土保持学院，北京 100083； 5.中国气象局 国家气象信息中心，北京 100081)
1 研究背景
北京山区泥石流以暴雨类、低频率、高重度、低黏度、沟谷型泥石流为主，而以支沟泥石流集群暴发为其特点。其形成方式主要为沟床起动型泥石流[1-5]，其形成机理是在短历时强降雨的作用下，大量雨水在地表形成径流，汇集到主沟道并形成强大的山洪，继而侵蚀沟床并起动以前滑坡和崩塌堆积在沟床内的固体物质形成泥石流[6-10]。北京山区经常暴发泥石流，如2011年“7·24”和2012年“7·21”强降雨引发山洪泥石流导致密云区龙潭沟、房山区十渡镇受损严重，严重危害了当地人民的生命安全[11]。目前，北京市突发地质灾害专业监测预警工作已经全面展开，监测预警的最终目标是实现每条泥石流沟的单沟预警预报工作，因此，单沟泥石流预警模型是其关键技术与核心内容[12]。
目前针对北京山区泥石流灾害研究，谢洪等[13]认为应该根据泥石流单沟特征，通过制定减灾工程减轻和控制泥石流的危害。钟敦伦等[14]通过结合泥石流的影响因素，构建物源模型，在可拓学原理下，通过关联度计算得出区域泥石流预报和沟谷泥石流预报。周金星等[15]建立了北京暴雨泥石流发生当天的前3 d、5 d或15 d雨量与当日激发泥石流雨量的回归模型。倪化勇等[16]利用前期降雨量研究了降雨型泥石流预警预报模式，并给出相应的预警框架和建议。史明远[17]通过分析研究区泥石流类型、致灾能力、发展阶段和易发度特性，结合泥石流单沟前期有效降雨量和短期降雨量，提出未来6 h、24 h的气象预警雨量阈值。白利平等[18-19]基于GIS技术，提出北京地区泥石流危险度区划，并根据力学平衡条件，假设饱水状态下的沟床物质会受到重力、摩擦力、内聚力、浮力等的作用，推出松散固体物质起动判别式和雨强表达式。黄晓虎[20]对北京地区达摩小流域泥石流的启动机制进行了研究，他指出泥石流物源体启动的标志是出现切沟侵蚀，并计算了各种降雨情况下泥石流启动临界激发雨量。丁桂伶等[21]选取岩性、坡度、植被等评价指标，通过前期有效降雨量计算，提出北京地区的泥石流预警雨量阈值。涂剑等[22]用修正后的李氏法进行雨场分割，总结研究区泥石流雨量激发的条件。
马超等[23]通过野外观测的10 min降雨记录和沟道形态变化，分析了泥石流的触发条件和侵蚀过程，提出了10 min的降雨间隔强度灾害预测以及水泥沙量是控制侵蚀程度的首要因素。翟淑花等[24]通过绘制降雨安全等概率线，对北京山区泥石流进行预警。朱昳橙等[25]通过研究降雨量与降雨过程、地质灾害空间分布的对应关系，得出地质-降雨综合模型比地质-气象模型更适合怒江州滑坡地质灾害预警。李永超等[26]通过对79条冲沟的野外调查，采用信息值法、确定因子法和逻辑回归法对研究区泥石流易发性进行评价。蒋家沟泥石流预报模型是根据很多流域大量、长系列的观测数据提出的，虽然蒋家沟模型与北京山区目前预警模型具有一定的相似性，都使用降雨因子作为判断指标，但忽略了沟道内松散堆积物和利于泥石流发生的沟谷地势，这肯定会导致预报结果不精确[27]。谭炳炎等[28]提出了泥石流组合预报模型：Y=RM，其中Y表示泥石流发生的判别函数；R表示降雨条件函数；M表示地面环境动态函数。当Y<25时，不发生泥石流的机率为83%；Y>35时，发生泥石流的机率为85 %；Y=25～35时，介于发生与不发生之间，该模型主要运用于山区铁路沿线的暴雨泥石流预报。王毅等[29]用GIS技术与层次-信息量模型相结合的方法，确定在不同降雨下需要重点监测与防治的泥石流沟道。李大鸣等[30]通过分析泥石流形成条件，选取了泥石流7个危险因子，并应用模糊赋权法计算各个危险因子权重，从而建立了泥石流区域预报模型。
综上所述，泥石流预警模型成果很多，在区域预警预报中考虑了地形、降雨因子，在雨量指标上都用了前期累计雨量和泥石流激发雨量。但预报的准确性鲜有报道，原因在于大多数泥石流预警模型具有局地性，不能适用于北京山区泥石流预警。
目前，北京泥石流灾害研究多为区域性泥石流预报，以及泥石流危险区划和泥石流沟道侵蚀过程等方向的研究。关于单沟泥石流短期(10 min和1 h)预警的研究还很少。本文在前人研究基础上，收集近十几年北京山区单沟泥石流文献资料，通过研究影响北京山区单沟泥石流形成的地形条件、地质条件和水文条件等，得出适用于北京山区单沟泥石流的Ⅳ级预警模型，可有效预防和减轻灾害造成的人员伤亡和财产损失。
2 泥石流的形成条件
北京山区面积为10 417.5 km2，占北京市总面积的62%。北京山区地处燕山山系与西山太行山系交界处，在地质构造单元上属燕山沉降带。境内地形骨架形成于中生代燕山运动，西部山区与北部山区均为高山陡坡地形，岭谷相对高差均在300～600 m之间，河谷断面多呈“V”形。由于山区切割较深，相对高度较大，沟谷的纵坡大，沟道狭窄，崩塌时有发生，沟道和坡面碎屑物丰富，物源体主要为粗颗粒，为泥石流的发生提供了有利的地形条件[2]。泥石流的物源主要有两类，一是沟谷底部的松散堆积物，二是谷坡上的堆积物。北京山区发生泥石流的主要物源是沟谷底部的松散堆积物，沟道两侧斜坡物源少，虽然降雨会导致谷坡上堆积物汇入沟道，但对泥石流的形成影响不大[31-32]。
“消防水管效应”现象对沟床松散固体物质冲刷强烈，导致其失稳并参与泥石流活动[33]。根据Yu等[33]的研究，沟床起动类型的沟谷泥石流预报模型的基本概念和方法是基于泥石流形成的3个基本条件，即地形条件、地质条件和降水条件。通过收集泥石流的沟谷地形、地质和水文条件方面的特征，对它们之间的关系进行对比研究，从泥石流的形成机理上出发，对泥石流的发生作出预测。本文通过北京市气象局、区(县)防汛办，以及国家气象局在北京山区部分气象站点，共搜集了1989—2016年间26起泥石流事件(如图1)，研究区泥石流情况统计如表1。
表1 研究区泥石流情况统计Table 1 Statistics of debris flows in the mountainous area of Beijing
图1 研究区泥石流灾害点分布Fig.1 Locations of debris flow sites in the mountainous area of Beijing
根据余斌等[6，37]建立的贵州省望谟县2011年6月6日特大暴雨作用下的沟床起动类型泥石流预报模型，泥石流起动条件由地形条件、地质条件和降水条件组成，其表达式为
(1)
式中：P为预警值；R为泥石流形成区降水因子；T为泥石流形成区地形因子；G为泥石流形成区地质因子；Cr为模型判据临界值。1 h预报模型临界值(Cr)有2个，即Cr1=0.35，Cr2=0.47(分别对应模型中警报线和避难线)。10 min预报模型临界值(Cr)也有2个，即Cr3=0.078，Cr4=0.103(分别对应模型中警报线和避难线)。
2.1 地形条件
地形条件指有利于沟床物源起动的地形条件：形成区流域面积、沟道长度、沟床纵比降、形状系数等，主要注重泥石流形成区流域地形条件的作用，并用一个无量纲的地形条件来描述，沟床起动型泥石流无量纲的地形条件表达式为[37]
(2)
式中：J为泥石流形成区沟床纵比降；A是泥石流形成区流域面积(km2)；A0是单位面积，1 km2；L是形成区沟道长度(km)；A/L2为泥石流形成区形状系数；A/A0为无量纲化的泥石流形成区面积。
2.2 地质条件
在泥石流发育的过程中，泥石流流域岩性的平均坚固系数是重要的参数。它与泥石流的暴发频率和泥石流发生难易程度关系为：泥石流流域的岩石越坚固，其固体物源主要为崩塌产生，固体物质颗粒相对较大，积累时间较长，泥石流的暴发频率越低，泥石流发生越困难；反之，泥石流流域的岩石越软弱，其固体物源则主要是滑坡产生，固体物质颗粒相对较小，泥石流发生越容易[38]。
泥石流形成的地质条件可以表示为泥石流形成区流域岩性的坚固系数、断裂带(构造)、地震烈度、物理风化和化学风化等修正条件(通过地质构造、地震烈度、物理风化对岩石坚固系数进行修正，结合化学风化修正系数，得到综合地质因素G，G越大，越不利于泥石流的形成)。并用一个无量纲的地质条件来描述，沟床起动型泥石流无量纲的地质条件表达式为[37]
G=F0C1C2C3C4。
(3)
式中：F0为泥石流形成区平均坚固系数；C1、C2、C3、C4分别为形成区构造(断裂带)、地震烈度、物理风化和化学风化修正因子(具体取值见文献[37]—文献[39])。研究区泥石流地形、地质条件情况统计见表2(发生泥石流的沟道位置是查阅文献所得，未发生泥石流的沟道是选取发生泥石流的沟道附近的沟道)。
表2 泥石流地形、地质条件情况统计Table 2 Statistics of topography and geological conditions of debris flows
2.3 降雨条件
虽然极端降雨事件是一个低概率的事件，可是极端降雨事件在中国北部发生的次数最近这几年在逐渐增加，北京“7·21”就是一个低概率的暴雨事件。短历时(1 h或10 min时间段)强降雨是诱发沟床起动类型的沟谷泥石流的主要原因[9]。美国学者Kean等[40]认为激发这类泥石流，短时间的强降雨比总降雨量的作用更重要。这类泥石流的标准化无量纲临界降雨量可以表示为[38]
(4)
式中：R*为综合降雨临界值；R0为当地年平均降雨量(mm)；CV为当地10 min降雨变差系数(尽管本文采用的降水强度是1 h降水强度，但沟床起动类型泥石流与10 min降水强度相关性更好[41])；B为泥石流暴发前期累积降雨量(mm)；K为系数，在1 h模型中，k=12.5，在10 min模型中，k=8；I为泥石流激发1 h(1 h模型)或10 min(10 min模型)降雨量[38](mm)。
3 北京山区泥石流预报模型
3.1 模型分析
由于北京山区泥石流也是沟床起动型，因此适合该预报模型，为了与国内蓝、黄、橙、红Ⅳ级预警接轨，本文按照比例增加一个临界值范围，即：
Cr11=Cr1×0.8=0.28 ；
(5)
Cr101=Cr3×0.8=0.062 。
(6)
改进后的北京山区泥石流1 h预报模型判据临界值分别为Cr11=0.28、Cr12=0.35和Cr13=0.47， 表3 1 h预警降雨条件情况统计Table 3 Statistics of rainfall conditions for 1 h early-warning 图2 1 h预报模型临界线Fig.2 Critical lines of 1 h forecasting 将26 条泥石流沟和9条未发生泥石流的沟道代入式(1)得到1 h预报模型临界线(图 2)，其中有27条泥石流沟在黄色预警区以上(P≥0.28)；有27条泥石流沟在橙色预警区以上(P≥0.35)；有25条泥石流沟在红色预警区以上(P≥0.47)。9条未发生泥石流的沟道都在红色预警范围内，野外调查发现主要是因为未发生泥石流的沟道物源比较少，因此无法发生泥石流灾害。可见该预警模型只能适用于北京山区有物源的泥石流沟1 h预警。 2016年8月12日，密云县龙潭沟风景区以及毗邻查子沟村遭受暴雨袭击，当日位于景区内的16号曹庄子泥石流沟[34]、17号大木峪泥石流沟和18号艾洼峪泥石流沟[26]均发生泥石流。本文发生泥石流的降雨数据中，仅仅收集到龙潭沟流域周边的白龙潭气象站和查子沟气象站10 min降雨数据，艾洼峪沟沟内的气象站只有1 h降雨数据。由于白龙潭气象站距离龙潭沟流域较近(图1)，能更好地反映降雨情况，故选用白龙潭气象站数据作为龙潭沟流域3条泥石流10 min降雨数据，艾洼峪沟气象站数据作为龙潭沟流域3条泥石流1 h降雨数据。2016年8月12日6点50分之前，只有8月7日有累积雨量，为12 mm。曹庄子泥石流沟对面的栗村峪当地居民见到曹庄子泥石流的时刻大约在上午9时，而大木峪和艾洼峪泥石流的发生时刻没有报道[34]。2016年8月12日6：50—9：00时段的白龙潭气象站10 min降雨过程如图3(a)，艾洼峪气象站1 h降雨过程如图3(b)。 图3 降雨过程Fig.3 Rainfall process 图4(a)和图4(b)分别为曹庄子沟、大木峪沟和艾洼峪沟在2016年8月12日1 h和10 min模型预警指标P值与临界值关系。 图4 预报模型Fig.4 Forecast model 由龙潭沟流域1 h预警模型的预警指标P的变化过程(图4(a))可得：6：00的预警指标P值可判断为蓝色预警；7：00—11：00的预警指标P值可判断为红色预警。结合泥石流暴发的实际时间，1 h预警模式可以很好判断这次的3条沟泥石流预警，很可能可以很好地判断北京山区泥石流的发生。 由龙潭沟流域10 min预警模型的预警指标P的变化过程(图4(b))可得：7：20—9：00的预警指标P值可判断为黄色预警；8：30—8：40的预警指标P值可判断为橙色预警；7：50—8：00和8：20的预警指标P值可判断为红色预警。 结合泥石流暴发的实际时间，10 min预判断模式可以很好地判断这次的3条沟泥石流预警，很可能可以很好地判断北京山区泥石流的发生。其中当1 h 和10 min预警指标P值在不同的预警区域内时，按Ⅳ级预警模型中危险性区域最高的发布预警。 北京市2012年7月21日暴发特大山洪泥石流，其中暴雨中心所在的房山区暴发多处泥石流。北京山区大多数泥石流为沟床起动型泥石流，其形成机理是在短历时强降雨的作用下，大量雨水在地表形成径流，汇集到主沟道形成强大的山洪，继而侵蚀沟床并起动以前滑坡和崩塌堆积在沟床内的固体物质形成泥石流。短历时的强降雨是北京山区泥石流发生的最主要诱发因素。本文从北京近十几年的观测资料入手，在泥石流形成机理的基础上，结合泥石流形成的3大条件：地形、地质和降雨以及相互之间的关系，采用适用于北京山区泥石流预报模型P=RT0.2/G0.5得出预警临界值，进行Ⅳ级预警，得到以下结论： (1)当泥石流沟道的1 h预警指标P≥ 0.47时和10 min预警指标P≥0.103时，极可能发生泥石流，需进行红色预警。 (2)北京山区单沟泥石流短期预报模型可以较好地根据降雨的变化情况得出泥石流沟道不同时间段的泥石流的易发程度。 判断沟道泥石流的易发性是基于假设沟道内是有物源的前提下进行预警的，因为本文大多数泥石流沟数据来源于查阅大量的文献。沟道没有物源，即使预警值达到了危险红色预警范围，也不会发生泥石流。并且采矿区、道路建设等随意堆放的废矿石渣等特殊泥石流物源，以及强风化花岗岩地区的特别物源，都有可能影响预报模型的准确性。沟床起动型泥石流主要与沟床物质和沟床坡降有关，因此野外对沟床内物源颗粒粒径的大小、沟床物源堆积等特征的调查十分必要。虽然采用修正后的平均坚固系数表示地质条件有些不足，但这能较好且快速定量地给出沟道物源情况的指标。 因此，野外现场调查能提高预警的准确性。由于目前收集到的北京山区10 min降雨数据较少，北京山区单沟泥石流10 min预警模型还需进一步验证。
3.2 龙潭沟流域预报案例
4 结 论
