本发明属于水文技术领域，具体涉及一种基于临界雨量的分布式中小流域地质灾害及洪水预警方法。
背景技术：
随着自动雨量监测技术的发展，无人值守的自动雨量监测站以其测量快速、布设便利的特点站网密度正在不断加大，应用实时监测的降雨对山洪陡涨陡落和滑坡、泥石流等地质灾害频发的山区性小流域进行快速预警现已成为中小流域灾害控制管理的重要内容之一。在应用实时监测的降雨进行中小流域洪水和地质灾害预警时，其主要难点之一在于如何结合流域内土壤水含量情况，确定出易致灾害的临界降雨量，从而结合实时的降雨情况对流域内发生的洪水和地质灾害可能性进行快速评估预警。为了进一步促进流域洪水和地质灾害预警预报的发展，需要更深入研究临界雨量的计算方法。
在以临界雨量为标准，采用实时雨量与临界雨量相对比的方法进行流域洪水和地质灾害预警时，临界雨量的合理性对发布预警的合理性有着重要影响。然而，目前方法计算出来的临界雨量一般是依据历史降雨洪水资料，然而在实时降雨过程中，随着降雨情况不断发生变化，流域内土壤水含量也在实时发生变化，使得与土壤水含量有着紧密联系的临界雨量也在实时发生着变化。因此目前方法计算出来的临界雨量并不完全符合客观实际，从而限制了采用实时降雨数据在流域洪水和地质灾害预警中的应用，不利于国内中小流域实时洪水和地质灾害预警的发展。
针对以上不足，如何考虑流域土壤水含量的实时变化，结合实时变化的降雨量和土壤水含量计算出合理的临界雨量，并将其与实时的降雨量比较作为中小流域洪水和地质灾害快速预警的依据，正是发明人需要解决的问题。
技术实现要素：
为了解决现有技术中存在的不足，本发明提供了基于临界雨量的分布式中小流域地质灾害及洪水预警方法，具有数据来源稳定可靠、计算效率高、结果客观合理等优点，有利于中小流域洪水和地质灾害快速预警。
为解决上述问题，本发明具体采用以下技术方案：
基于临界雨量的分布式中小流域地质灾害及洪水预警方法，其特征在于，包括以下步骤：
步骤1，提取流域dem数据，将流域划分为若干个栅格；
步骤2，利用流域dem数据提取流域水系；
步骤3，依据各个站点的小时雨量数据，得到流域内的各个栅格上逐小时雨量数据；
步骤4，构建基于上层张力水含量和累计降雨量的临界雨量预警模型；
步骤5，实时洪水预报预警中，结合流域内的各个栅格上逐小时雨量数据和临界雨量预警模型，计算得到各个栅格单元内的临界雨量，并实现流域内分区分级预警；
步骤6，结合流域内的水系，依据流域内分区分级预警结果实现河流分段分级预警。
前述的基于临界雨量的分布式中小流域地质灾害及洪水预警方法，其特征在于，所述步骤1中提取流域dem数据，将流域划分为若干个栅格，具体包括以下步骤：
步骤1.1，填洼；
步骤1.2，计算流向；
步骤1.3，计算汇流流量，得到汇流栅格；
步骤1.4，确定流域出口站点；
步骤1.5，提取目标流域。
前述的基于临界雨量的分布式中小流域地质灾害及洪水预警方法，其特征在于，所述步骤2中利用流域dem数据提取流域水系，具体包括以下步骤：
步骤2.1，设置汇流阈值t；
步骤2.2，基于步骤1.3中的汇流栅格，依据汇流阈值t，汇流累计值低于阈值的栅格单元判断为坡地栅格，汇流累计值高于阈值的栅格单元判断为河道栅格，得到流域中的水系river；
步骤2.3，将2.2中得到的水系river与流域实际水系影像比较，若相差较大则修改阈值t再次提取水系，直到水系river与流域实际水系符合。
前述的基于临界雨量的分布式中小流域地质灾害及洪水预警方法，其特征在于，所述步骤3中依据各个站点的小时雨量数据，得到流域内的各个栅格上逐小时面雨量数据，具体包括以下步骤：
步骤3.1，对各个雨量站的日雨量资料和小时雨量资料进行整理，得到逐日的日雨量资料和逐小时的小时雨量资料；
步骤3.2，利用各个雨量站的日雨量资料对小时雨量资料进行校核；
式中：pd为日雨量数据；pht为1日中逐小时降雨量数据，t为逐小时的编码，从1到24；δp为日雨量和该日中逐小时降雨累计之和的差值，
步骤3.3，基于校核后各个雨量站小时雨量资料进行空间插值，得到流域上各个栅格单元内的逐小时雨量数据。
前述的基于临界雨量的分布式中小流域地质灾害及洪水预警方法，其特征在于，所述步骤4中构建基于上层张力水含量和累计降雨量的临界雨量预警模型，具体包括以下步骤：
步骤4.1，选择流域内洪峰较高、危害较大的典型洪水；
步骤4.2，设置不同的初始土壤水含量，并提前一段预热期，采用分布式的水文模型进行产流计算，得到每一个栅格单元内的逐小时上层张力水含量；
步骤4.3，在每一个栅格单元内，统计在不同初始土壤水含量的情况下，上层张力水含量达到饱和时所对应的累积降雨量accps，即该土壤水含量下的临界雨量；
步骤4.4，在每一个栅格单元内，统计能够使上层张力水含量达到饱和的最小累积降雨量出现时不同初始土壤水含量情况下所对应的各自上层张力水含量ws；
步骤4.5，在每一个栅格单元内，基于s个ws和accps点对，以ws为横坐标，accps点为纵坐标拟合曲线，得到曲线方程；
accps＝α×wsβ+μ(3)
式中：α为一个栅格单元内曲线方程的比例系数；β为一个栅格单元内曲线方程的次方系数；μ为一个栅格单元内的最大每小时可能降雨量；s为4.2中不同初始土壤水含量的种类数目，一般情况下s取6，分别对应的初始土壤水含量为饱和土壤水含量的0％，20％，40％，60％，80％，100％；
逐栅格确定各个栅格单元中的系数α、β和μ，从而确定每一个栅格单元内的临界雨量预警模型方程表达式，进而构建全流域的临界雨量预警模型。确定每一个栅格单元中的系数α、β和μ就能够确定每一个栅格单元中的临界雨量预警模型方程表达式，在每一个栅格单元的临界雨量预警模型方程表达式确定好之后就确定整个流域的临界雨量预警模型。
前述的基于临界雨量的分布式中小流域地质灾害及洪水预警方法，其特征在于，所述步骤5中在实时洪水预报预警中，结合流域内的各个栅格上逐小时面雨量数据和临界雨量预警模型，计算得到各个栅格单元内的临界雨量，并实现流域内分区分级预警，具体包括以下步骤：
步骤5.1，将各个站点的实时雨量插值为面雨量，并进行累加计算，得到实时累计降雨量空间分布栅格rastersump；
步骤5.2，以每一个栅格单元内的实时雨量为分布式水文模型的输入，计算得到每一个栅格单元内逐时刻上层张力水含量；
步骤5.3，将逐时刻上层张力水含量作为临界雨量预警模型的输入，结合步骤4.5中得到的流域中系数α的矩阵、系数β的矩阵和系数μ的矩阵计算出逐时刻的各个栅格单元的临界雨量，得到临界雨量空间分布栅格rastercrip；
步骤5.4，将5.1中得到的实时累计降雨量空间分布栅格减去5.3中得到的临界雨量空间分布栅格得到预警雨量判断栅格rasterexcp，若栅格单元中实时累计雨量高于临界雨量，则该栅格为需要进行预警的栅格单元，并结合流域实际情况对超过临界雨量的栅格单元进行分级；
在每一个栅格单元中：
rasterexcpe＝rastersumpe-rastercripe
式中：rasterexcpe为超过临界雨量的那部分雨量，rastersumpe为栅格单元内的累积雨量，rastercripe为栅格单元内的临界雨量，e为流域内栅格的编码，从1到sumraster，sumraster为流域内栅格单元的总数；
当0≤rasterexcpe<amt1时对编号为e的栅格发布蓝色预警；当amt1≤rasterexcpe<amt2时对编号为e的栅格发布黄色预警；当amt2≤rasterexcpe<amt3时对编号为e的栅格发布橙色预警；当amt3≤rasterexcpe时对编号为e的栅格发布红色预警，所述amt1、amt2以及amt3为流域地质灾害预警分级标准。一般来讲，将amt1设置为20，amt2设置为40，amt3设置为60。同时amt1、amt2以及amt3可依据各个流域多年实际情况进行相应调整。
前述的基于临界雨量的分布式中小流域地质灾害及洪水预警方法，其特征在于，所述步骤6中结合流域内的水系，依据流域内分区分级预警结果实现河流分段分级预警，具体包括以下步骤：
步骤6.1，访问预警雨量判断栅格中的每个栅格单元，并且计算每一个栅格单元n*n邻域范围内数值的均值ae均，滑块平均后的预警雨量判断栅格：
式中：ae均为e栅格单元n*n邻域范围内像元数据的均值，e为流域内栅格的编码，从1到sumraster，sumraster为流域内栅格单元的总数；i为e栅格单元n*n邻域范围内像元的编号，从1到n，n为邻域范围的边长；aei为每一栅格单元内的值；
步骤6.2，依据水系和滑块平均后的预警雨量判断栅格之间的相对位置关系，选择出水系中每一个河段经过的滑块平均后的预警雨量判断栅格的栅格单元；
步骤6.3，计算6.2中每一个河段经过的所有滑块平均后预警雨量判断栅格的栅格单元均值，并将该值赋予该河段；
式中：bj均为每一个河段经过的所有滑块平均后预警雨量判断栅格的栅格单元均值；j为河段编码，从1到r，r为水系中的河段总数；k为每一河段经过的栅格单元的编码，从1到m，m为每一河段经过的栅格单元总数；bjk为编码为j的河段经过的预警雨量判断栅格的栅格单元中的值，bjk就是预警雨量判断栅格的栅格单元中的值，类似于aei，亦即是每一条河段经过的栅格的aei，为以示河道预警与流域中面上的预警的区别，故单独命名为bjk以示区分；
步骤6.4，结合流域中的实际情况，依据6.3中河段被赋予的预警雨量的值对水系中所有河段进行分区分级预警，对各个河段发生山洪的可能性进行分级；
当
num1、num2以及num3为流域河段中山洪预警分级标准。一般来讲，将num1设置为10，num2设置为20，num3设置为30。同时num1、num2以及num3可依据各个流域多年实际情况进行相应调整。
本发明的有益效果：本发明提供的基于临界雨量的分布式中小流域地质灾害及洪水预警方法，根据土壤水蓄水容量的物理意义，结合流域实测降雨数据，建立了土壤含水量和实测雨量情况下与临界雨量之间的函数关系，推求了实时土壤含水量和实时雨量下的临界雨量，进而通过实测雨量与临界雨量的对比提取了超临界雨量警戒栅格的空间分布，结合流域内的水系，依据流域内分区分级预警结果实现河流分段分级预警。此方法从实时雨量资料和土壤水含量的出发，以实时计算得到的临界雨量作为发布预警的依据，既保证了预警发布的可靠性和实时性，又解决了无资料山区小流域的洪水和地质灾害预警问题，且方法中变量之间的函数关系明确，有利于临界雨量的快速自动生成，提高了山洪和地质灾害预警的自动化程度和实时性，同时通过分布式的计算方法充分利用了输入雨量数据的空间分布信息，考虑了临界雨量的空间变异性，保证了结果的客观合理性，可以进一步促进山区中小流域洪水和地质灾害预警的深入发展。
附图说明
图1是本发明的系统框图；
图2为本发明中提取出的流域河道示意图；
图3为本发明中通过空间插值得到的流域上各个栅格单元内的雨量值示意图；
图4为本发明中选择的典型洪水示意图；
图5为本发明中计算得到栅格单元内逐小时上层张力水含量变化过程的示意图；
图6为本发明中拟合得到的上层张力水含量和临界雨量预警模型示意图；
图7为本发明中计算得到的某时刻上层张力水含量空间分布示意图；
图8为本发明中计算得到的某时刻临界雨量空间分布示意图；
图9为本发明中计算得到的某时刻预警雨量判断栅格空间分布示意图；
图10为本发明中某时刻滑块平均后的预警雨量判断栅格空间分布示意图；
图11为本发明中判定的某时刻发布预警的河流空间分布示意图。
具体施方式
下面结合附图和具体实施例对本发明作进一步描述。
如图1所示，本发明提供的基于临界雨量的分布式中小流域地质灾害及洪水预警方法，包括以下步骤：
s1、提取流域数字高程模型(dem)数据，将流域划分为若干个栅格，具体包括以下步骤：
步骤1.1，填洼；
步骤1.2，计算流向；
步骤1.3，计算汇流流量，得到汇流栅格；
步骤1.4，确定流域出口站点；
步骤1.5，提取目标流域。
s2、利用流域dem数据提取流域水系，具体包括以下步骤：
步骤2.1，设置汇流阈值t；
步骤2.2，基于步骤1.3中的汇流栅格，依据汇流阈值t，汇流累计值低于阈值的栅格单元判断为坡地栅格，汇流累计值高于阈值的栅格单元判断为河道栅格，得到流域中的水系river，如图2所示；
步骤2.3，将步骤2.2中得到的水系river与流域实际水系影像比较，若相差较大则修改阈值t再次提取水系，直到水系river与流域实际水系基本符合。
s3、依据各个站点的小时雨量数据，得到流域内的各个栅格上逐小时雨量数据，具体包括以下步骤：
步骤3.1，对各个雨量站的日雨量资料和小时雨量资料进行整理，得到逐日的日雨量资料和逐小时的小时雨量资料；
步骤3.2，利用各个雨量站的日雨量资料对小时雨量资料进行校核；
式中：pd为日雨量数据；pht为1日中逐小时降雨量数据，t为逐小时的编码，从1到24；δp为日雨量和该日中逐小时降雨累计之和的差值，
步骤3.3，基于校核后各个雨量站小时雨量资料进行空间插值，得到流域上各个栅格单元内的雨量值，如图3所示。
s4、构建基于上层张力水含量和累计降雨量的临界雨量预警模型，具体包括以下步骤：
步骤4.1，选择流域内洪峰较高，危害较大的典型洪水，如图4所示；
步骤4.2，设置不同的初始土壤水含量，并提前一段预热期，采用分布式的水文模型进行产流计算，得到每一个栅格单元内的逐小时上层张力水含量，如图5所示；
步骤4.3，在每一个栅格单元内，统计在不同初始土壤含水量的情况下，上层张力水含量达到饱和时所对应的累积降雨量accps，即该土壤水含量下的临界雨量；
步骤4.4，在每一个栅格单元内，统计能够使上层张力水含量达到饱和的最小累积降雨量出现时不同初始土壤含水量情况下所对应的各自上层张力水含量ws；
步骤4.5，在每一个栅格单元内，基于s个ws和accps点对，以ws为横坐标，accps点为纵坐标拟合曲线，得到曲线方程，如图6所示；
accps＝α×wsβ+μ(3)
式中：α为一个栅格单元内曲线方程的比例系数；β为一个栅格单元内曲线方程的次方系数；μ为一个栅格单元内的最大每小时可能降雨量；s为4.2中不同初始土壤水含量的种类数目，一般情况下s取6，分别对应的初始土壤水含量为饱和土壤水含量的0％，20％，40％，60％，80％，100％。
逐栅格确定各个栅格单元中的系数α、β和μ，从而确定每一个栅格单元内的临界雨量预警模型方程表达式，进而构建全流域的临界雨量预警模型。确定每一个栅格单元中的系数α、β和μ就能够确定每一个栅格单元中的临界雨量预警模型方程表达式，在每一个栅格单元的临界雨量预警模型方程表达式确定好之后就确定整个流域的临界雨量预警模型。
s5、在实时洪水预报预警中，结合流域内的各个栅格上逐小时面雨量数据和临界雨量预警模型，计算得到各个栅格单元内的临界雨量，并实现流域内分区分级预警，具体包括以下步骤：
步骤5.1，将各个站点的实时雨量插值为面雨量，并进行累加计算，得到实时累计降雨量空间分布栅格；
步骤5.2，以每一个栅格单元内的实时雨量为分布式水文模型的输入，计算得到每一个栅格单元内逐时刻上层张力水含量，如图7所示；
步骤5.3，将逐时刻上层张力水含量作为临界雨量预警模型的输入，结合步骤4.5中得到的流域中系数α的矩阵、系数β的矩阵和系数μ的矩阵计算出逐时刻的各个栅格单元的临界雨量，得到临界雨量空间分布栅格，如图8所示；
步骤5.4，将5.1中得到的实时累计降雨量空间分布栅格减去5.3中得到的临界雨量空间分布栅格得到预警雨量判断栅格，若栅格单元中实时累计雨量高于临界雨量，则该栅格为需要进行预警的栅格单元，并结合流域实际情况对超过临界雨量的栅格单元进行分级，如图9所示。
在每一个栅格单元中：
rasterexcpe＝rastersumpe-rastercripe
式中：rasterexcpe为超过临界雨量的那部分雨量，rastersumpe为栅格单元内的累积雨量，rastercripe为栅格单元内的临界雨量，e为流域内栅格的编码，从1到sumraster，sumraster为流域内栅格单元的总数。
当0≤rasterexcpe<amt1时对编号为e的栅格发布蓝色预警；当amt1≤rasterexcpe<amt2时对编号为e的栅格发布黄色预警；当amt2≤rasterexcpe<amt3时对编号为e的栅格发布橙色预警；当amt3≤rasterexcpe时对编号为e的栅格发布红色预警。
amt1、amt2以及amt3为流域地质灾害预警分级标准。一般来讲，将amt1设置为20，amt2设置为40，amt3设置为60。同时amt1、amt2以及amt3可依据各个流域多年实际情况进行相应调整。
s6、结合流域内的水系，依据流域内分区分级预警结果实现河流分段分级预警的步骤包括：
步骤6.1，访问预警雨量判断栅格中的每个栅格单元，并且计算每一个栅格单元n*n邻域范围内数值的均值ae均，滑块平均后的预警雨量判断栅格，如图10所示：
式中：ae均为e栅格单元n*n邻域范围内像元数据的均值，e为流域内栅格的编码，从1到sumraster，sumraster为流域内栅格单元的总数；i为e栅格单元n*n邻域范围内像元的编号，从1到n，n为邻域范围的边长；aei为每一栅格单元内的值。
栅格中将要被访问的每个栅格单元称为待处理像元。待处理像元的值以及所识别出的邻域中的所有像元值都将包含在邻域统计数据计算中。同时各邻域可以重叠，因此一个邻域中的像元也可以包含在其他待处理像元的邻域中。
步骤6.2，依据水系和滑块平均后的预警雨量判断栅格之间的相对位置关系，选择出水系中每一个河段经过的滑块平均后的预警雨量判断栅格的栅格单元；
步骤6.3，计算6.2中每一个河段经过的所有滑块平均后预警雨量判断栅格的栅格单元均值，并将该值赋予该河段；
式中：bj均为每一个河段经过的所有滑块平均后预警雨量判断栅格的栅格单元均值；j为河段编码，从1到r，r为水系中的河段总数；k为每一河段经过的栅格单元的编码，从1到m，m为每一河段经过的栅格单元总数，bjk为编码为j的河段经过的预警雨量判断栅格的栅格单元中的值，bjk就是预警雨量判断栅格的栅格单元中的值，类似于aei，亦即是每一条河段经过的栅格的aei，为以示河道预警与流域中面上的预警的区别，故单独命名为bjk以示区分；
步骤6.4，结合流域中的实际情况，依据6.3中河段被赋予的预警雨量的值对水系中所有河段进行分区分级预警，对各个河段发生山洪的可能性进行分级，如图11所示。
当
num1、num2以及num3为流域河段中山洪预警分级标准。一般来讲，将num1设置为10，num2设置为20，num3设置为30。同时num1、num2以及num3可依据各个流域多年实际情况进行相应调整。
另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。
以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。