论文部分内容阅读
摘 要：山洪是山区中影响范围广，人员伤亡较大的一类地质灾害。针对山洪预警模型，国内外学者从多个角度对其进行过研究，研究区域多为大流域，方法多为雨量临界值预警。为了建立一套适用于北京山区小流域的山洪预警模型，本文研究了小流域单沟洪水模型预警在系统，此系统是基于充分的野外地质调查，在获得详细的沟域地质及土地利用数据后，利用合理的数学模型对流域内的降雨所产生的洪水产汇流过程进行全过程模拟，并对汇流过程中可能产生的灾害进行预警预报的一个系统。此单沟洪水模型具有精度高，计算结果准确的特点，其相应的预警系统能够对小范围居民区等重要建构筑物进行灾害预报，并预估人员撤离的安全时间，对实际的防灾减灾工作具有重要意义。
关键词：山洪;小流域;预警;模型;水流流速;北京
中图分类号：P694 文献标识码：A 文章编号：1007-1903（2019）02-0018-07
Abstract： Mountain flood is a kind of common geological disaster in mountainous area， which has a wide influence and heavy casualties. As for the mountain flood warning model， domestic and foreign scholars have studied it from multiple perspectives. Most of the research objects are focused on the large watersheds， and the methods are adopted as rainfall critical value warning. In order to establish a set of flash flood warning model suitable for Beijing mountain watershed， this paper studies the single-trench flood early-warning model of small watershed. This system is established as a system of early warning based on plenty of field geological investigation， the detailed field geological and land use data， using a reasonable mathematical model of the basin rainfall flood aquatic confluence process generated by the whole process of simulation. The single-trench flood model has the characteristics of high accuracy and accurate calculation results， and its corresponding early-warning system can predict the disaster of important structures such as small residential areas， and estimate the safe time of evacuation， which is of great significance to the actual disaster prevention and reduction work. The system has a good demonstration effect on the calculation and three-dimensional simulation of mountain flood runoff in Beijing area， and provides decision-making basis for related departments in disaster prediction， personnel evacuation and standard formulation.
Keywords： Mountain flood; Small watershed; Early warning， Model; Flow velocity; Beijing
0引言
北京地區河流均属海河流域，主要有永定河、潮白河、北运河、大清河和蓟运河五大水系。除北运河外，其余各河均发源于北京市境外。北京市属温带半干旱、半湿润季风区，多年平均降水量585mm。降水量年际变化较大，年内分布不均，汛期降水约占全年降水量的85%（王亚娟等，2016）。洪涝灾害是北京历史上发生频率既高而又危害最严重的自然灾害，北京洪灾大致可分为大河洪灾、城市洪灾、涝渍灾害以及山洪泥石流灾害和泥沙灾害。据不完全统计，1949—1990年，北京大小洪涝灾害累计受灾面积2770万亩次。全市因洪涝灾害造成农业损失1.82亿元（按当年价计算），死亡人口246人，倒塌房屋13.63间（王金如等，1998）。2012年“7·21”特大暴雨灾害中，山洪泥石流灾害尤为严重，全市多条已干涸多年的山洪、泥石流沟道洪水瞬时达到历史最大值;截至2012年7月31日，全市受灾人口127.48万人，紧急转移9.27万人，直接经济损失161.57亿元（霍风霖等，2014）。山区小流域因山洪沟坡陡、流域面积小，具有突发性强、破坏性大、预见期短、暴涨暴落、难预报和难预防等特点。所以，只有通过及时有效的山洪灾害预报预警才可获得较长预见期，以减少人员伤亡和财产损失（马传波，2016）。因而，在北京地区因地制宜地选用恰当的模型，建立完善的山洪预警系统，提高预警准确性势在必行（宋瑶等，2016）。 为了对山区暴雨和洪水做出有效的预报和预警，国内外学者近年来从多个不同的领域进行了许多有意义的研究（王礼先等，2001）。由于山洪和强对流降水的密切关系，气象工作者主要着眼于一些改进触发山洪的强降水监测预警方面的研究（李致家等，2004;张火青等，1996;张廷治等，1996）。水文防汛部门则采用传统的水文学临界雨量阈值法，以及一些较为稳定的水文产汇流模型对山洪进行预报，根本上也极其依赖于气象预报（芮孝芳，2004）。我国从20世纪60年代起，以各省、市、自治区为单位，进行大规模的水文调查和研究，各省分别制定了水文手册和暴雨洪水查算手册（邵利萍，2009），为各地的小流域水文计算提供了依据。总的来说，小流域洪水计算的核心思想是利用各地区的暴雨洪水观测资料，研究其规律并建立其与流域特征之间的经验关系（刘川佳仔，2013）。在小流域产汇流计算方法的研究上，一般利用实测资料拟合建立暴雨强度、汇流时间、降雨损失和流域特性之间的数量关系，可根据流域土质、植被情况分别确定相应的参数并做成表格，以便日后查用，也可利用迭代法求解（罗东，2009;徐德龙等，2000;潭惠平，1999;姜典，1984;陈家琦等，1983）。
区别于以往单纯采用雨量作为预警判别条件来对山洪进行预报，本文研究的北京山区单沟洪水预警模型采用水动力学原理，利用降雨作为输入条件，考虑地形、地质条件、土地利用类型等因素，全面准确的对山区小流域洪水的流量及流速进行预估，为相关部门提供有效的防灾减灾的科学依据。
1 系统设计思路
模型采用的计算方法是基于水动力学的基本方程，包括流体力学的连续性方程和运动方程，进行修正后利用Riemann通量格式计算。整个模型包括洪水产汇流部分以及预警预报部分，其中产汇流部分包含三个模块，地形数据处理、下渗率及糙率赋值、汇流及模型计算。针对某特定山洪单沟，结合图1所示的预警流程，本研究利用所采集的高程数据、降雨数据、土地利用类型数据以及地质岩性数据，代入模型计算，获得系统预警所需的结果。
（1）地形数据处理
地形数据处理是模型计算中最基础的模块，处理地形数据，首先需要根据单沟沟域的山脊、沟道等明显的地形特征在图片上勾勒出泥石流沟流域的边界。其次在Arcgis中对DEM数据进行配准。接下来依据已勾勒的边界剖分网格，选取网格中心点，计算出坐标并插值得出中心点高程数据。
（2）下渗率及糙率赋值
该模块主要用来计算不同地质条件的下渗和植被条件的汇流情况，利用得到的沟域地质图以及植被分布图，依据行业内已有的标准，对相应的地质岩性赋予渗透性系数，对相应的植被类型赋予糙率系数，在计算中分别对应于土壤的下渗以及汇流过程，前者主要影响水流的下渗量，后者影响水流的流速及流量。
（3）汇流及模型计算
汇流模块是模型计算的核心模块，输入标准化的降雨数据，利用前期处理过的地形数据与植被和地质数据，代入模型进行计算。此模型重点考虑陡坡汇流情況，在原有的二维浅水波方程的基础上做了对连续方程的源汇项以及动量方程的降雨动量做了修正，采用有限体积Riemann通量计算格式，建立基于网格的水动力陡坡汇流计算方法。
（4）系统预警部分
预警预报部分是系统应用的最重要的一个模块，基于前述产汇流部分的计算结果，洪水单沟流域的任意一点在某一降雨过程中的流量与流速均可以得出。有别于以往根据降雨量等级判定洪水级别的预警方法，本系统采用更为直接的利用流量和流速进行联合预警。两者的关系可以是同时达到预警值预警或其中一者达到预警值预警。为了能够有效的对村庄、农田、重要建构筑物进行保护，我们选择沿着汇流方向，距离保护区域300m处设置预警点，预警等级分为蓝、黄、橙、红4个等级，与4个预警等级相对应的流量与流速值可以进行人工设定，通常依据经验或者实验模拟得出。
2 数学模型设计
（1）产汇流基本方程
北京市山洪预警系统主要应用于北京市西山以及北山地区的单沟小流域，为了有效的模拟实际山洪的产汇流过程，本系统利用的数学模型是在溃坝大坡度陡变水流计算模式的基础上，采用修正的二维浅水波方程，基本方程如式（2-1）。
式中：h为水深;u为局部坐标x方向的流速;qe为源汇项;α为地形坡度;φ为降雨角度;g为重力加速度;S0为底坡比降;Sf为阻力项;σ为风应力;ρ为空气密度;p为降雨强度;vm为雨滴降落速度;d为下渗强度;i为蒸发强度。
由于沟域边界一般为山脊线，故边界处不考虑外界水流流入，只以降水量作为输入水量。初始时间除降水量之外，各个物理量均为0。
模型算法的基本思路是，基于修正后的流体力学连续性方程和运动方程，对研究区域采用适当宽度的距离进行网格剖分，计算方式采用水位与流量时间交错（李大鸣等，2011），考虑到沟域内不同通道情况，使用不同的通道流量公式。
连续方程：
运动方程：
式中：h为水深，z为水位，z= z0+h，z0为底高程，u、 v分别为x和y方向的平均流速，q为源汇项，g为重力加速度，n为糙率。
经过对两个方程的计算可得出流域汇流的计算结果。
（2）地面渗流计算公式
地面渗流是计算地表水经过土地下渗后的地表产流，在土壤含水率达到饱和前地表产流强度时，由于土壤入渗的速率低于降水强度，产生的地表径流由下式计算：
式中：rs为地表径流，i为降水强度，fm为最大入渗速率，f为透水层平均入渗速率。若土壤入渗率进入平稳期，则根据霍顿入渗公式，采用下式进行计算：
其中，n =2.5;Rp为时段降水量;h为时段地表径流。
模型在参考溃坝陡坡汇流模型时，由于陡坡或水层汇流时，划分单元较小，坡降较大，易出现水流量小但计算流量偏大的情况，这种情况我们称之为虚假流量现象。为了解决这一问题，需适当的减少计算时的时间步长，或对计算模式有所修正。
（3）地表糙率计算公式
土地利用类型数据以及地质数据对模型的汇流过程具有重要影响。土地利用类型是通过地表不同的区域的糙率系数来对汇流过程的流速产生影响，进而对最终的流量和流速产生作用。其取值一般依据国家防汛抗旱总指挥部编写的《洪水风险图编制导则》提供的经验系数，以及实验室以往建立水力学模型的经验进行修正。地质数据主要是考虑到水流下渗过程中，不同岩层渗透性系数的不同，从而在不同区域造成下渗量的差异，导致汇流流量的变化。其赋值依据《水利水电工程地质勘察规范》（GB 50287-99）规定，按照不同的岩体特征和土类划分渗透性等级和系数。
糙率主要影响谢才系数C，进而影响到通道流速和流量，其中两者之间的关系是通过谢才公式建立的。计算通道上的流速和流量时，采用谢才公式，如下：
式中：v为断面平均流速（m/s），R为水力半径（m），J为水力坡度，C为谢才系数。
本次模型主要采用巴甫洛夫斯基公式来计算C值，即：
式中：n为糙率。
3 预警功能实现
预警模块是模型的应用模块，是实现模型计算的最重要的功能。预警一般以小流域为单元，研究确定符合当地实际的动态预警指标，提高预警精准度，延长预见期，增强预警时效（孙东亚，2019）。基于上一小节模型的产汇流部分的计算结果，特定沟域内任意一点的流量与流速均可以得出。不同于以往单纯利用雨量大小及强度进行沟域山洪预警值设定，这里我们利用流量与流速值进行联合预警。首先确定预警点位置，在山洪的沟域内，村庄等人居或者重要建构筑物，通常是在沟域下游至沟口的位置，这里是山洪汇流的最终聚集点。为了有效的对村庄等区域进行保护，预留人员撤离时间，我们一般将预警点设置在沿着汇流方向距离村庄300m位置处。预警值的设定基于流量和流速共同确定，两者可以同时达到设定值预警，也可其中之一达到设定值预警。预警级别分四级，分别为蓝、黄、橙、红，每一级的设定值均由人为设定，可根据历史经验或实验分析来获得。
系统中将预警功能与沟域山洪汇流展示结合在一起，当沟域出现暴雨山洪时，水流流经预警点，当流量及流速触及某一预警等级的设定值时，预警警报响起，预警点处闪烁相应级别的旗帜，同时在展示页面展示自第一次预警触发到最后一次预警解除时，预警点的流速及流量的数值表，并显示相应预警级别的颜色。
4 结果验证
4.1 孙胡沟地质背景
孙胡沟，位于怀柔区中部的琉璃庙镇境内。全村总面积约18735亩，村内共有人口686人，辖9个自然村。整个村南北长、东西短，东窄西宽呈月牙斧状。东与琉璃庙镇崎峰茶村、西台子交界，南与怀柔区雁栖镇交界，西与琉璃庙镇河北村和鱼水洞村毗邻，北与崎峰茶村接壤。属华北平原军都山系，年平均气温13℃～15℃，主要农作物为玉米，黄豆和高粱等小杂粮，盛产板栗、核桃、大扁仁。
孙胡沟历史上曾发生过多次山洪泥石流灾害，最近一次发生于1969年8月10日。在该次泥石流发生前，当地连续降雨达20天，至8月10日则以枣树林为中心的军都山南北坡普降大雨到暴雨，降雨量高达297mm/d，最大降雨强度达132mm/h。多日连续降雨导致斜坡土体饱和，加之当地坡体植被几乎全被破坏，当晚孙胡沟暴发大型山洪，主沟及几条支沟均有人员伤亡，全村总计死亡31人，造成了重大人员财产损失。
（1）岩性
孙胡沟出露的地层由老至新有太古界、中新元古界、中生界和第四系，其岩性主要为侵入岩和白云岩。其中，侵入岩主要包括两期，一期为太古代侵入的英云閃长岩，主要分布在孙胡沟的中部;另一类为中生代侵入的二长花岗岩，主要分布在孙胡沟南端。白云岩主要为中元古界长城系白云岩，主要分布在孙胡沟北侧。
孙胡沟内的第四系主要分布于各主要支沟沟床和少量残坡积物，多数为全新统冲积、冲洪积碎石、砂卵石、亚粘土沉积物，厚度0～20m，砾径数厘米至数米，分选性差。
（2）地形地貌
孙胡沟地区海拔高度多在430～ 1160m，南部高，北部低。总的地貌特点是坡陡峰峻，地形起伏大，相对高差在300～500m。
（3）土地利用类型
根据孙胡沟的遥感影像，解译后得到土地利用类型可分为林地、庄稼地、农村宅基地、裸地、果园、草地等。
4.2 模型准确性验证
我们利用沟内的流速仪数据对模型准确性进行验证。首先在系统中输入计算所需的各类数据（图2）。实测点位于计算区域内，主沟偏上游位置，图3中“☆”位置（图中线条为地面等高线）。
选取2018年8月17日22：00：00至2018年8月18日10：30：00的降雨数据作为模型的验证样本。通过模型计算，计算并输出实测点所在位置网格（编号：6246）的流速，并与实测值进行比较。
流速仪于验证时间段内测得的数据列于表1第二列中，观察整体数据，发现在未降雨的情况下，仪器仍然有一初始流速——经计算大约为1.33m/s，故本场降雨所产生的流速应当比仪器显示流速小1.33m/s，整理于表1第三列。
以21：55：00为零点，每5min为一个单位，绘制校正后实测流速随时间的变化，如图4所示。橙色的点为流速的仪器实测值，蓝色曲线为流速的模型计算值。
从图4中可以看出，模型计算值与实测值趋势基本相同，实测值皆落于计算曲线附近。将计算值与实测值对应数据列于表2当中。 从表2中可以看出，验证的7个时间点中，绝对误差均小于0.08m/s;最小误差为0.004073m/s，出现在2018/8/18 4：50：00;实测值与计算值的平均值相差0.017313m/s。
产生误差的原因分析：①由于仪器测量位置以及模拟计算的精度，实测点与模拟的计算通道不完全重合，存在系统误差;②实际情况复杂，存在诸多计算模型中无法完全考虑到的情况，如是否存在暗渠、较小的坑塘洼地等;③实测数据精度受实测环境和器材精度影响，存在机械误差，如当地风速对流速仪旋桨的影响。
综上所述，计算值与实测值走势相同，绝对误差相差较小，模拟效果良好，计算值可靠。
图5为沟域在某次降雨过程中所计算出的积水深度，其中细红色线条为沟域的边界，蓝色表示积水淹没区，绿色表示未被积水淹没区，沟域内的粗红线为人为设置的预警横断线。此图表示在实际发生的某次的降雨过程中，人为设置的三条预警横断线分别出现了黄色、橙色、红色预警（以相应颜色的旗帜表示）。浅蓝色箭头表示洪水流出方向。
5 结论
（1）本文建立的小流域单沟山洪产汇流计算模型，选取示范沟并对其历史降雨过程进行了模拟，模拟效果良好，误差较小，说明该模型适用于北京山区小流域的洪水模拟。模型的特点，一是基于山区地形，以陡坡汇流的方式对原洪水汇流过程的计算方法进行了修正，二是考虑了小流域的地质岩性以及土地利用类型，使得模型计算更加符合特定沟域的自身条件，更具针对性。
（2）系统在基于模拟状况良好的模型基礎上，在沟域内设置特定预警点，以流速和流量联合预警的方式对山洪进行预报预警。这种方法具有预警失误率低，预报时间准确以及山洪影响面积及程度可度量的特点，并且能够为下游人员撤离预留充足时间。该系统适用于复杂环境下居民区及重要区域的预警预报，对应急避险具有很强的实际效用。
（3）系统采用C#语言编写，模型的各个模块以及集成的运行效率直接影响到系统运行的效率，故在保证原有程序不出现bug的基础上，增加对程序语言的优化能够提升计算效率并减少对计算资源的占用，是模型改进的一个重要的方向。
参考文献
陈家琦，张恭肃，1983.推理公式的应用与改进[J].水文， （1）：2-9.
霍风霖，张力，2014.北京山洪灾害防治体系现状与建议[J].中国防汛抗旱， 24（3）：13-15.
姜典，1984.推理公式法在暴雨时深关系不连续的应用[J].水文， （1）：11-12+62.
李致家，刘金涛，葛文忠，等，2004.雷达估测降雨与水文模型的耦合在洪水预报中的应用[J].河海大学学报（自然科学版）， （6）：601-606.
李大鸣，吕会娇，2011.山区暴雨泥石流预报数学模型的研究[J].中国农村水利水电， （6）：24-28+35.
刘川佳仔，2013. 北京小流域产流特点及暴雨洪水计算研究[D].清华大学.
罗东，2009.推求中小流域暴雨洪水方法探讨[J].气象水文海洋仪器， 26（2）：158-161.
马传波， 2016.山区小流域洪水预警临界雨量值确定方法的研究与应用[J].吉林水利，（6）：50-52+62.
芮孝芳， 2004.水文学原理[M].北京：中国水利水电出版社.
宋瑶，程金花，胡晓静，等， 2016.基于垂向混合产流模型的北京山区山洪预警模型构建[J].水土保持通报， 36（3）：353-357.
邵利萍， 2009. 无资料小流域洪水叠加计算方法初探[D].浙江大学.
孙东亚，2019.中国山洪灾害监测预警体系建设[J].中国防汛抗旱， 29（2）：4-5.
潭惠平，1999.特小流域设计洪水计算方法简析[J].河北水利水电技术， （2）：23-24.
王亚娟，赵小伟，臧敏，等， 2016.北京市“7·20”特大暴雨洪水分析[J].北京水务， （5）：1-6.
王金如，顾小明，1998. 北京的水害及防治[C]//中国灾害防御协会，北京减灾协会.中国减灾与新世纪发展战略——首届“中国21世纪安全减灾与可持续发展战略高级研讨会”论文集：5.
王礼先，于志民，2001.山洪及泥石流灾害预报[M].北京：中国林业出版社：1-23.
徐德龙，肖华， 2000.小流域设计洪水推理公式计算方法探讨[J].人民长江， （7）：13-14.
张火青，魏文秋，1996.天气雷达在洪水预警报中的应用[J].水利水电科技进展， （3）：21-25.
张廷治，李守智，李祥云，等，1996.诱发山洪泥石流特强暴雨的特征[J].气象， （5）：43-47.
关键词：山洪;小流域;预警;模型;水流流速;北京
中图分类号：P694 文献标识码：A 文章编号：1007-1903（2019）02-0018-07
Abstract： Mountain flood is a kind of common geological disaster in mountainous area， which has a wide influence and heavy casualties. As for the mountain flood warning model， domestic and foreign scholars have studied it from multiple perspectives. Most of the research objects are focused on the large watersheds， and the methods are adopted as rainfall critical value warning. In order to establish a set of flash flood warning model suitable for Beijing mountain watershed， this paper studies the single-trench flood early-warning model of small watershed. This system is established as a system of early warning based on plenty of field geological investigation， the detailed field geological and land use data， using a reasonable mathematical model of the basin rainfall flood aquatic confluence process generated by the whole process of simulation. The single-trench flood model has the characteristics of high accuracy and accurate calculation results， and its corresponding early-warning system can predict the disaster of important structures such as small residential areas， and estimate the safe time of evacuation， which is of great significance to the actual disaster prevention and reduction work. The system has a good demonstration effect on the calculation and three-dimensional simulation of mountain flood runoff in Beijing area， and provides decision-making basis for related departments in disaster prediction， personnel evacuation and standard formulation.
Keywords： Mountain flood; Small watershed; Early warning， Model; Flow velocity; Beijing
0引言
北京地區河流均属海河流域，主要有永定河、潮白河、北运河、大清河和蓟运河五大水系。除北运河外，其余各河均发源于北京市境外。北京市属温带半干旱、半湿润季风区，多年平均降水量585mm。降水量年际变化较大，年内分布不均，汛期降水约占全年降水量的85%（王亚娟等，2016）。洪涝灾害是北京历史上发生频率既高而又危害最严重的自然灾害，北京洪灾大致可分为大河洪灾、城市洪灾、涝渍灾害以及山洪泥石流灾害和泥沙灾害。据不完全统计，1949—1990年，北京大小洪涝灾害累计受灾面积2770万亩次。全市因洪涝灾害造成农业损失1.82亿元（按当年价计算），死亡人口246人，倒塌房屋13.63间（王金如等，1998）。2012年“7·21”特大暴雨灾害中，山洪泥石流灾害尤为严重，全市多条已干涸多年的山洪、泥石流沟道洪水瞬时达到历史最大值;截至2012年7月31日，全市受灾人口127.48万人，紧急转移9.27万人，直接经济损失161.57亿元（霍风霖等，2014）。山区小流域因山洪沟坡陡、流域面积小，具有突发性强、破坏性大、预见期短、暴涨暴落、难预报和难预防等特点。所以，只有通过及时有效的山洪灾害预报预警才可获得较长预见期，以减少人员伤亡和财产损失（马传波，2016）。因而，在北京地区因地制宜地选用恰当的模型，建立完善的山洪预警系统，提高预警准确性势在必行（宋瑶等，2016）。 为了对山区暴雨和洪水做出有效的预报和预警，国内外学者近年来从多个不同的领域进行了许多有意义的研究（王礼先等，2001）。由于山洪和强对流降水的密切关系，气象工作者主要着眼于一些改进触发山洪的强降水监测预警方面的研究（李致家等，2004;张火青等，1996;张廷治等，1996）。水文防汛部门则采用传统的水文学临界雨量阈值法，以及一些较为稳定的水文产汇流模型对山洪进行预报，根本上也极其依赖于气象预报（芮孝芳，2004）。我国从20世纪60年代起，以各省、市、自治区为单位，进行大规模的水文调查和研究，各省分别制定了水文手册和暴雨洪水查算手册（邵利萍，2009），为各地的小流域水文计算提供了依据。总的来说，小流域洪水计算的核心思想是利用各地区的暴雨洪水观测资料，研究其规律并建立其与流域特征之间的经验关系（刘川佳仔，2013）。在小流域产汇流计算方法的研究上，一般利用实测资料拟合建立暴雨强度、汇流时间、降雨损失和流域特性之间的数量关系，可根据流域土质、植被情况分别确定相应的参数并做成表格，以便日后查用，也可利用迭代法求解（罗东，2009;徐德龙等，2000;潭惠平，1999;姜典，1984;陈家琦等，1983）。
区别于以往单纯采用雨量作为预警判别条件来对山洪进行预报，本文研究的北京山区单沟洪水预警模型采用水动力学原理，利用降雨作为输入条件，考虑地形、地质条件、土地利用类型等因素，全面准确的对山区小流域洪水的流量及流速进行预估，为相关部门提供有效的防灾减灾的科学依据。
1 系统设计思路
模型采用的计算方法是基于水动力学的基本方程，包括流体力学的连续性方程和运动方程，进行修正后利用Riemann通量格式计算。整个模型包括洪水产汇流部分以及预警预报部分，其中产汇流部分包含三个模块，地形数据处理、下渗率及糙率赋值、汇流及模型计算。针对某特定山洪单沟，结合图1所示的预警流程，本研究利用所采集的高程数据、降雨数据、土地利用类型数据以及地质岩性数据，代入模型计算，获得系统预警所需的结果。
（1）地形数据处理
地形数据处理是模型计算中最基础的模块，处理地形数据，首先需要根据单沟沟域的山脊、沟道等明显的地形特征在图片上勾勒出泥石流沟流域的边界。其次在Arcgis中对DEM数据进行配准。接下来依据已勾勒的边界剖分网格，选取网格中心点，计算出坐标并插值得出中心点高程数据。
（2）下渗率及糙率赋值
该模块主要用来计算不同地质条件的下渗和植被条件的汇流情况，利用得到的沟域地质图以及植被分布图，依据行业内已有的标准，对相应的地质岩性赋予渗透性系数，对相应的植被类型赋予糙率系数，在计算中分别对应于土壤的下渗以及汇流过程，前者主要影响水流的下渗量，后者影响水流的流速及流量。
（3）汇流及模型计算
汇流模块是模型计算的核心模块，输入标准化的降雨数据，利用前期处理过的地形数据与植被和地质数据，代入模型进行计算。此模型重点考虑陡坡汇流情況，在原有的二维浅水波方程的基础上做了对连续方程的源汇项以及动量方程的降雨动量做了修正，采用有限体积Riemann通量计算格式，建立基于网格的水动力陡坡汇流计算方法。
（4）系统预警部分
预警预报部分是系统应用的最重要的一个模块，基于前述产汇流部分的计算结果，洪水单沟流域的任意一点在某一降雨过程中的流量与流速均可以得出。有别于以往根据降雨量等级判定洪水级别的预警方法，本系统采用更为直接的利用流量和流速进行联合预警。两者的关系可以是同时达到预警值预警或其中一者达到预警值预警。为了能够有效的对村庄、农田、重要建构筑物进行保护，我们选择沿着汇流方向，距离保护区域300m处设置预警点，预警等级分为蓝、黄、橙、红4个等级，与4个预警等级相对应的流量与流速值可以进行人工设定，通常依据经验或者实验模拟得出。
2 数学模型设计
（1）产汇流基本方程
北京市山洪预警系统主要应用于北京市西山以及北山地区的单沟小流域，为了有效的模拟实际山洪的产汇流过程，本系统利用的数学模型是在溃坝大坡度陡变水流计算模式的基础上，采用修正的二维浅水波方程，基本方程如式（2-1）。
式中：h为水深;u为局部坐标x方向的流速;qe为源汇项;α为地形坡度;φ为降雨角度;g为重力加速度;S0为底坡比降;Sf为阻力项;σ为风应力;ρ为空气密度;p为降雨强度;vm为雨滴降落速度;d为下渗强度;i为蒸发强度。
由于沟域边界一般为山脊线，故边界处不考虑外界水流流入，只以降水量作为输入水量。初始时间除降水量之外，各个物理量均为0。
模型算法的基本思路是，基于修正后的流体力学连续性方程和运动方程，对研究区域采用适当宽度的距离进行网格剖分，计算方式采用水位与流量时间交错（李大鸣等，2011），考虑到沟域内不同通道情况，使用不同的通道流量公式。
连续方程：
运动方程：
式中：h为水深，z为水位，z= z0+h，z0为底高程，u、 v分别为x和y方向的平均流速，q为源汇项，g为重力加速度，n为糙率。
经过对两个方程的计算可得出流域汇流的计算结果。
（2）地面渗流计算公式
地面渗流是计算地表水经过土地下渗后的地表产流，在土壤含水率达到饱和前地表产流强度时，由于土壤入渗的速率低于降水强度，产生的地表径流由下式计算：
式中：rs为地表径流，i为降水强度，fm为最大入渗速率，f为透水层平均入渗速率。若土壤入渗率进入平稳期，则根据霍顿入渗公式，采用下式进行计算：
其中，n =2.5;Rp为时段降水量;h为时段地表径流。
模型在参考溃坝陡坡汇流模型时，由于陡坡或水层汇流时，划分单元较小，坡降较大，易出现水流量小但计算流量偏大的情况，这种情况我们称之为虚假流量现象。为了解决这一问题，需适当的减少计算时的时间步长，或对计算模式有所修正。
（3）地表糙率计算公式
土地利用类型数据以及地质数据对模型的汇流过程具有重要影响。土地利用类型是通过地表不同的区域的糙率系数来对汇流过程的流速产生影响，进而对最终的流量和流速产生作用。其取值一般依据国家防汛抗旱总指挥部编写的《洪水风险图编制导则》提供的经验系数，以及实验室以往建立水力学模型的经验进行修正。地质数据主要是考虑到水流下渗过程中，不同岩层渗透性系数的不同，从而在不同区域造成下渗量的差异，导致汇流流量的变化。其赋值依据《水利水电工程地质勘察规范》（GB 50287-99）规定，按照不同的岩体特征和土类划分渗透性等级和系数。
糙率主要影响谢才系数C，进而影响到通道流速和流量，其中两者之间的关系是通过谢才公式建立的。计算通道上的流速和流量时，采用谢才公式，如下：
式中：v为断面平均流速（m/s），R为水力半径（m），J为水力坡度，C为谢才系数。
本次模型主要采用巴甫洛夫斯基公式来计算C值，即：
式中：n为糙率。
3 预警功能实现
预警模块是模型的应用模块，是实现模型计算的最重要的功能。预警一般以小流域为单元，研究确定符合当地实际的动态预警指标，提高预警精准度，延长预见期，增强预警时效（孙东亚，2019）。基于上一小节模型的产汇流部分的计算结果，特定沟域内任意一点的流量与流速均可以得出。不同于以往单纯利用雨量大小及强度进行沟域山洪预警值设定，这里我们利用流量与流速值进行联合预警。首先确定预警点位置，在山洪的沟域内，村庄等人居或者重要建构筑物，通常是在沟域下游至沟口的位置，这里是山洪汇流的最终聚集点。为了有效的对村庄等区域进行保护，预留人员撤离时间，我们一般将预警点设置在沿着汇流方向距离村庄300m位置处。预警值的设定基于流量和流速共同确定，两者可以同时达到设定值预警，也可其中之一达到设定值预警。预警级别分四级，分别为蓝、黄、橙、红，每一级的设定值均由人为设定，可根据历史经验或实验分析来获得。
系统中将预警功能与沟域山洪汇流展示结合在一起，当沟域出现暴雨山洪时，水流流经预警点，当流量及流速触及某一预警等级的设定值时，预警警报响起，预警点处闪烁相应级别的旗帜，同时在展示页面展示自第一次预警触发到最后一次预警解除时，预警点的流速及流量的数值表，并显示相应预警级别的颜色。
4 结果验证
4.1 孙胡沟地质背景
孙胡沟，位于怀柔区中部的琉璃庙镇境内。全村总面积约18735亩，村内共有人口686人，辖9个自然村。整个村南北长、东西短，东窄西宽呈月牙斧状。东与琉璃庙镇崎峰茶村、西台子交界，南与怀柔区雁栖镇交界，西与琉璃庙镇河北村和鱼水洞村毗邻，北与崎峰茶村接壤。属华北平原军都山系，年平均气温13℃～15℃，主要农作物为玉米，黄豆和高粱等小杂粮，盛产板栗、核桃、大扁仁。
孙胡沟历史上曾发生过多次山洪泥石流灾害，最近一次发生于1969年8月10日。在该次泥石流发生前，当地连续降雨达20天，至8月10日则以枣树林为中心的军都山南北坡普降大雨到暴雨，降雨量高达297mm/d，最大降雨强度达132mm/h。多日连续降雨导致斜坡土体饱和，加之当地坡体植被几乎全被破坏，当晚孙胡沟暴发大型山洪，主沟及几条支沟均有人员伤亡，全村总计死亡31人，造成了重大人员财产损失。
（1）岩性
孙胡沟出露的地层由老至新有太古界、中新元古界、中生界和第四系，其岩性主要为侵入岩和白云岩。其中，侵入岩主要包括两期，一期为太古代侵入的英云閃长岩，主要分布在孙胡沟的中部;另一类为中生代侵入的二长花岗岩，主要分布在孙胡沟南端。白云岩主要为中元古界长城系白云岩，主要分布在孙胡沟北侧。
孙胡沟内的第四系主要分布于各主要支沟沟床和少量残坡积物，多数为全新统冲积、冲洪积碎石、砂卵石、亚粘土沉积物，厚度0～20m，砾径数厘米至数米，分选性差。
（2）地形地貌
孙胡沟地区海拔高度多在430～ 1160m，南部高，北部低。总的地貌特点是坡陡峰峻，地形起伏大，相对高差在300～500m。
（3）土地利用类型
根据孙胡沟的遥感影像，解译后得到土地利用类型可分为林地、庄稼地、农村宅基地、裸地、果园、草地等。
4.2 模型准确性验证
我们利用沟内的流速仪数据对模型准确性进行验证。首先在系统中输入计算所需的各类数据（图2）。实测点位于计算区域内，主沟偏上游位置，图3中“☆”位置（图中线条为地面等高线）。
选取2018年8月17日22：00：00至2018年8月18日10：30：00的降雨数据作为模型的验证样本。通过模型计算，计算并输出实测点所在位置网格（编号：6246）的流速，并与实测值进行比较。
流速仪于验证时间段内测得的数据列于表1第二列中，观察整体数据，发现在未降雨的情况下，仪器仍然有一初始流速——经计算大约为1.33m/s，故本场降雨所产生的流速应当比仪器显示流速小1.33m/s，整理于表1第三列。
以21：55：00为零点，每5min为一个单位，绘制校正后实测流速随时间的变化，如图4所示。橙色的点为流速的仪器实测值，蓝色曲线为流速的模型计算值。
从图4中可以看出，模型计算值与实测值趋势基本相同，实测值皆落于计算曲线附近。将计算值与实测值对应数据列于表2当中。 从表2中可以看出，验证的7个时间点中，绝对误差均小于0.08m/s;最小误差为0.004073m/s，出现在2018/8/18 4：50：00;实测值与计算值的平均值相差0.017313m/s。
产生误差的原因分析：①由于仪器测量位置以及模拟计算的精度，实测点与模拟的计算通道不完全重合，存在系统误差;②实际情况复杂，存在诸多计算模型中无法完全考虑到的情况，如是否存在暗渠、较小的坑塘洼地等;③实测数据精度受实测环境和器材精度影响，存在机械误差，如当地风速对流速仪旋桨的影响。
综上所述，计算值与实测值走势相同，绝对误差相差较小，模拟效果良好，计算值可靠。
图5为沟域在某次降雨过程中所计算出的积水深度，其中细红色线条为沟域的边界，蓝色表示积水淹没区，绿色表示未被积水淹没区，沟域内的粗红线为人为设置的预警横断线。此图表示在实际发生的某次的降雨过程中，人为设置的三条预警横断线分别出现了黄色、橙色、红色预警（以相应颜色的旗帜表示）。浅蓝色箭头表示洪水流出方向。
5 结论
（1）本文建立的小流域单沟山洪产汇流计算模型，选取示范沟并对其历史降雨过程进行了模拟，模拟效果良好，误差较小，说明该模型适用于北京山区小流域的洪水模拟。模型的特点，一是基于山区地形，以陡坡汇流的方式对原洪水汇流过程的计算方法进行了修正，二是考虑了小流域的地质岩性以及土地利用类型，使得模型计算更加符合特定沟域的自身条件，更具针对性。
（2）系统在基于模拟状况良好的模型基礎上，在沟域内设置特定预警点，以流速和流量联合预警的方式对山洪进行预报预警。这种方法具有预警失误率低，预报时间准确以及山洪影响面积及程度可度量的特点，并且能够为下游人员撤离预留充足时间。该系统适用于复杂环境下居民区及重要区域的预警预报，对应急避险具有很强的实际效用。
（3）系统采用C#语言编写，模型的各个模块以及集成的运行效率直接影响到系统运行的效率，故在保证原有程序不出现bug的基础上，增加对程序语言的优化能够提升计算效率并减少对计算资源的占用，是模型改进的一个重要的方向。
参考文献
陈家琦，张恭肃，1983.推理公式的应用与改进[J].水文， （1）：2-9.
霍风霖，张力，2014.北京山洪灾害防治体系现状与建议[J].中国防汛抗旱， 24（3）：13-15.
姜典，1984.推理公式法在暴雨时深关系不连续的应用[J].水文， （1）：11-12+62.
李致家，刘金涛，葛文忠，等，2004.雷达估测降雨与水文模型的耦合在洪水预报中的应用[J].河海大学学报（自然科学版）， （6）：601-606.
李大鸣，吕会娇，2011.山区暴雨泥石流预报数学模型的研究[J].中国农村水利水电， （6）：24-28+35.
刘川佳仔，2013. 北京小流域产流特点及暴雨洪水计算研究[D].清华大学.
罗东，2009.推求中小流域暴雨洪水方法探讨[J].气象水文海洋仪器， 26（2）：158-161.
马传波， 2016.山区小流域洪水预警临界雨量值确定方法的研究与应用[J].吉林水利，（6）：50-52+62.
芮孝芳， 2004.水文学原理[M].北京：中国水利水电出版社.
宋瑶，程金花，胡晓静，等， 2016.基于垂向混合产流模型的北京山区山洪预警模型构建[J].水土保持通报， 36（3）：353-357.
邵利萍， 2009. 无资料小流域洪水叠加计算方法初探[D].浙江大学.
孙东亚，2019.中国山洪灾害监测预警体系建设[J].中国防汛抗旱， 29（2）：4-5.
潭惠平，1999.特小流域设计洪水计算方法简析[J].河北水利水电技术， （2）：23-24.
王亚娟，赵小伟，臧敏，等， 2016.北京市“7·20”特大暴雨洪水分析[J].北京水务， （5）：1-6.
王金如，顾小明，1998. 北京的水害及防治[C]//中国灾害防御协会，北京减灾协会.中国减灾与新世纪发展战略——首届“中国21世纪安全减灾与可持续发展战略高级研讨会”论文集：5.
王礼先，于志民，2001.山洪及泥石流灾害预报[M].北京：中国林业出版社：1-23.
徐德龙，肖华， 2000.小流域设计洪水推理公式计算方法探讨[J].人民长江， （7）：13-14.
张火青，魏文秋，1996.天气雷达在洪水预警报中的应用[J].水利水电科技进展， （3）：21-25.
张廷治，李守智，李祥云，等，1996.诱发山洪泥石流特强暴雨的特征[J].气象， （5）：43-47.