本发明涉及气象学领域，特别是一种泥石流预警方法。
背景技术：
泥石流是山区常见的突发地质灾害之一，具有暴发突然、冲击力大、破坏力强等特点，一旦发生将会对人们生命财产带来巨大的损失。因此，泥石流预警预报是地质灾害主管部门及众多学者研究的焦点问题。
一直以来，针对降水在泥石流爆发中的触发作用公开了大量研究成果及工程实践经验。如，在统计分析研究区域已有泥石流灾害的基础上，选取多种形式的雨量组合作为预警指标，并采用线性拟合、Logistic模型或神经网络模型建立泥石流预警模型，并基于这些模型进行泥石流预测。
这些方法对泥石流的区域预警起了一定作用，但在适应性、可靠性方面还存在不足。
技术实现要素：
(一)要解决的技术问题
针对现有技术中泥石流预警方法存在的适应性较差、可靠性不高的问题，本发明提供了一种泥石流预警方法，旨在提高泥石流预警的普适性和可靠性。
(二)技术方案
为解决以上问题，本发明提供的技术方案为：
一种泥石流预警方法，包括：
获取泥石流待预警区域在观测时段t时各关键影响因子的具体数值；
利用多区域泥石流触发预警二叉树，确定所述泥石流待预警区域在观测时段t下对应的关键触发因子的预警数值。
具体地，所述的方法，还包括：
根据收集的有记载的泥石流发生区域在泥石流发生时的地质数据和降雨数据，采用降维方法确定泥石流触发时的关键地质因子和关键触发因子；
所述降维方法为以下任一种：采用主成分法、最小熵分析法、免疫进化算法和层次分析法。
具体地，所述的方法中，
所述地质数据包括：流域面积、主沟长度、沟床平均坡度、流域坡向、岩性、松散层厚度、松散物质比率和植被覆盖度；
所述降雨数据包括：10分钟雨强、30分钟雨强、小时雨强、降雨延时、最大小时雨强、当日雨量、前期有效雨量、总降雨量。
具体地，所述的方法中，
所述关键地质因子包括：流域面积、沟床平均坡度、松散物质比率；
所述关键触发因子包括：最大小时雨强与总降雨量的乘积。
具体地，所述的方法中，所述确定泥石流触发时的关键地质因子和关键触发因子的步骤之后，还包括：
根据收集的有记载的泥石流发生区域在泥石流发生时的所述关键地质因子和关键触发因子对应的值，利用基于二叉树的遗传规划方法确定描述所述关键地质因子与所述关键降雨因子之间函数关系的多区域泥石流触发预警二叉树。
具体地，所述的方法中，所述利用基于二叉树的遗传规划方法确定描述所述关键地质因子与所述关键降雨因子之间函数关系的多区域泥石流触发二叉树的步骤，包括：
步骤S10：获取收集的m个有记载的泥石流发生区域在泥石流发生时的所述关键地质因子和关键触发因子对应的值，确定m个所述关键地质因子对应的值为输入样本数据，m个所述关键触发因子对应的值为输出样本数据；
步骤S20：初始化：
获取终止符集；
获取函数集；
获取遗传规划各参数：群体规模n、交换概率α、突变概率β、最大突变深度γ、最大迭代次数T、模型误差δ、最大二叉树深度ζ、及终止准则，所述终止准则为：达到最大迭代次数T或不大于所述模型误差δ；
根据所述终止符集和所述函数集，采用完全法随机产生一个具有所述群体规模的二叉树群体，所述二叉树群体中的每一个个体均为一个不大于所述最大二叉树深度的二叉树；
步骤S30：确定二叉树群体内每一个二叉树的适应度，判断是否满足所述终止准则？
若是，则终止迭代，且转到步骤S40；
若否，则执行以下操作：
采用竞技选择法，从二叉树群体中，确定待执行复制操作的父代二叉树并执行复制操作；
根据所述交换概率α，从二叉树群体中，选择至少一个待执行交换操作的二叉树组，每个二叉树组内包括两个待执行交换操作的二叉树；
在每个待执行交换操作的二叉树组中，采用均匀分布的随机方法选取待执行交换操作的叶节点或中间节点并执行交换操作；
根据所述突变概率β及最大突变深度γ，从二叉树群体中，选择待执行突变的至少一个二叉树并对所述待执行突变的二叉树的叶节点或中间节点执行突变操作；
以增量为1的方式递增当前的迭代次数；
S40：若满足最大迭代次数而终止迭代，则确定当前群体中具有最大的适应度的二叉树为所述多区域泥石流触发二叉树；
若满足迭代误差不大于所述模型误差δ而终止迭代，则确定所述迭代误差不大于所述模型误差δ的二叉树为所述多区域泥石流触发二叉树。
具体地，所述的方法，在所述步骤S20中，
所述终止符集中的元素，包括：流域面积、沟床平均坡度、松散物质比率及至少一个随机数，其中，每个随机数在一个取值范围内均匀取值；
所述函数集中的元素，包括：四则运算符、三角函数、ex、及求平方根。
具体地，所述的方法中，所述根据所述终止符集和所述函数集，采用完全法随机产生初始群体的步骤，包括：
从所述函数集中选取元素做二叉树的根节点；
根据所述最大二叉树深度ζ，依次从所述函数集或所述终止符集中选取元素作为待定节点：
在待定节点的深度小于所述最大二叉树深度时，从所述函数集中随机选取元素作为所述待定节点；
在待定节点的深度等于所述最大二叉树深度时，从所述终止符集中选取元素作为所述待定节点。
具体地，所述的方法，在步骤S30中，
所述确定每一个二叉树的适应度的步骤，具体为：采用第一公式计算每一个二叉树的迭代误差，确定所述迭代误差的相反数为所述二叉树的适应度，所述第一公式为：
其中，pi为第i个输入样本数据对应于所述二叉树的预测输出， 1≥i≥m,oi为所述第i个输入样本数据对应的第i个输出样本数据。
(三)有益效果
本发明提供的泥石流预警方法根据表征不同沟谷泥石流降雨指数与其流域面积、松散物质率以及沟床平均坡度之间的非线性函数关系的多区域泥石流触发预警二叉树，实现了不同地质环境背景下的泥石流预警，精度高，可靠性高。
与现有技术相比，本发明提供的泥石流预警方法综合利用泥石流孕灾地质环境因子和降雨因子进行泥石流预警，普适性强、可靠性高，能够有效实现泥石流的预警。
附图说明
图1为本发明实施例的泥石流预警方法流程示意图；
图2为本发明实施例的泥石流预警方法中利用完全法生成的一个二叉树样本；
图3为本发明实施例中训练得到的临界降雨指数预测遗传规划模型图。
具体实施方式
为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。
泥石流是一种多发且频发的地质灾害，常对人民生命财产安全带来极大的威胁。有学者采用线性拟合来处理泥石流预警问题，结果显示，得到的预警模型精度不高；也有学者借助于人工神经网络较强的非线性映射函数的逼近能力建立泥石流预警模型并提高了预测精度，但基于人工神经网络建立的泥石流预警模型中，其结构形式的选择是不断训练和学习的过程，预测精度在很大程度上依赖于使用者的经验与水平，因此限制了模型推广和应用，造成目前泥石流预警的误报率和错报率比较高。
另一方面，目前泥石流的预警大多仅考虑降雨这一直观的泥石流触发因素。而并未涉及泥石流地形地貌、地质岩性等地质环境背景。
本发明提供的泥石流预警方法在获取的泥石流地质因子和降雨因子综合样本库的基础上，基于二叉树的遗传规划建立了泥石流的临界降雨指数预测的遗传规划模型，利用该遗传规划模型进行多个区域的泥石流预警，结果显示预警的精度高、可靠性强，能够达到及时预警的目的。
基于二叉树的遗传规划仿照生物界进化和遗传的过程，采用复杂多变的树状结构表达多个变量之间的线性或非线性数学关系，适用于处理大型、复杂的非线性问题。
为了更好地理解本发明，对基于二叉树的遗传规划的原理进行简单介绍。遗传规划一般包括的步骤为：1)选择一个初始群体；2)评价每个个体的适应度；3)进行选择(即复制)操作；4)进行交换和突变操作5)重复步骤2)-4)直至满足终止条件。在遗传规划中，通过对种群中个体的复制、交换和突变操作，能够搜索整个解空间，以取得最优解。
基于二叉树的遗传规划，随机产生一系列二叉树作为初始种群，种群中每个二叉树的节点值均取自终止集和函数集这两个集合。其中，终止集中的元素是要被操作的数据，函数集中的元素则是用来对数据进行操作的运算符。因此，遗传规划操作的对象就是整个二叉树群体。
如图1所示，本发明实施例的泥石流预警方法，包括：
步骤S100：获取的泥石流待预警区域在观测时段t时各关键影响因子的具体数值；
步骤S200：利用多区域泥石流触发预警二叉树，确定所述泥石流待预警区域在观测时段t下对应的关键触发因子的预警数值。
本发明提供的泥石流预警方法根据表征不同沟谷泥石流降雨指数与其流域面积、松散物质率以及沟床平均坡度之间的非线性函数关系的多区域泥石流触发预警二叉树，实现了不同地质环境背景下的泥石流预警，精度高，可靠性高。
进一步地，根据确定的关键触发因子的预警数值，结合该泥石流待预警区域当前的关键触发因子的值，就可以进行泥石流的预警。
具体地，所述的方法，还可以包括：
根据收集的有记载的泥石流发生区域在泥石流发生时的地质数据和降雨数据，采用降维方法确定泥石流触发时的关键地质因子和关键触发因子；所述降维方法为以下任一种：采用主成分法、最小熵分析法、免疫进化算法和层次分析法。
具体实施时，逐一利用上述降维方法根据收集的有记载的泥石流发生区域在泥石流发生时的地质数据和降雨数据，确定泥石流触发时的关键地质因子和关键触发因子；并对多种方法的模型复杂度和消耗的计算时间进行对比，并选择模型复杂度和计算时间均可以满足需求的方法作为最终使用的方法，并在后续步骤中使用利用该方法确定的关键地质因子和关键触发因子。
具体地，所述的方法中，所述地质数据可以包括：流域面积、主沟长度、沟床平均坡度、流域坡向、岩性、松散层厚度、松散物质比率和植被覆盖度；所述降雨数据可以包括：10分钟雨强、30分钟雨强、小时雨强、降雨延时、最大小时雨强、当日雨量、前期有效雨量、总降雨量。
需要说明的是，以上各术语均具有本领域内的一般含义。如，降雨延时是一次降雨过程的降雨的持续时间；前期有效雨量可以是前3天的有效雨量、前7天的有效雨量、或前15天的有效雨量。总降雨量则为前期有效降雨与当日降雨之和。
具体地，所述的方法中，所述关键地质因子包括：流域面积、沟床平均坡度、松散物质比率；所述关键触发因子包括：最大小时雨强与总降雨量的乘积。
泥石流是泥、砂、砾及巨石等固体物质与水的混合物，受重力作用后所产生的流动体。泥石流的爆发具有随机性和不确定性，受地形、地貌、地质以及人类工程活动等众多因子的影响。因此，泥石流预警模型因子的选择对预警精度具有重要影响，若选择因子太多，则会增加模型的负担，降低的模型的实际应用效果，若考虑的因子太少，又会因考虑不周而降低了模型的预测精度及可靠性。
鉴于此，在充分借鉴国内外已有研究成果的基础上，选择流域面积、松散物质比率、沟床平均坡度为影响泥石流临界雨量的自变量，以最大小时雨强和总累计雨量的乘积为因变量，以建立不同地质环境背景下的泥石流沟临界降雨指数预测模型。
1)流域面积：流域面积以流域出口为基准点，基准点以上的以分水岭为界的流域面积。流域面积的大小决定了汇水量的多少，流域面积也大，流域总汇水量越大。
2)松散物质率：松散物质的面积与泥石流发生与否具有较强的相关性，松散物质越多，固体物质越丰富，其爆发泥石流的可能性就越大。松散物质率是指流域内松散物质面积与流域面积的比值，其大小与岩层理化性质、地质构造、地下水的多寡、边坡坡度、植被状况等因素有关。
3)沟床平均坡度：泥石流发生流动需要一定的沟床坡度，若沟床坡度较缓，则无法提供泥石流运动所需的条件，沟床坡度越陡，泥石流动力越强。沟床平均坡度＝tan-1(高程差/主沟长度)，为达到单位的统一，采用平均坡度的正切值作为计算因子。
4)最大小时雨强与总累计雨量：基于泥石流触发的水理机制，直接激发雨量和前期降雨量是促使泥石流发生的重要参数，直接激发雨量常采用最大小时雨强，总累计雨量则为当日雨量与前期有效雨量之和，前期降雨量则采用前5天有效降雨量。为简化计算，以最大小时雨强(T) 和总降雨量(R)的乘积作为综合降雨驱动指标，并除以1000记作RTI 为降雨指数。
具体地，所述的方法中，所述确定泥石流触发时的关键地质因子和关键触发因子的步骤之后，还包括：
根据收集的有记载的泥石流发生区域在泥石流发生时的所述关键地质因子和关键触发因子对应的值，利用基于二叉树的遗传规划方法确定描述所述关键地质因子与所述关键降雨因子之间函数关系的多区域泥石流触发预警二叉树。
具体地，所述的方法中，所述利用基于二叉树的遗传规划方法确定描述所述关键地质因子与所述关键降雨因子之间函数关系的多区域泥石流触发二叉树的步骤，包括：
步骤S10：获取收集的m个有记载的泥石流发生区域在泥石流发生时的所述关键地质因子和关键触发因子对应的值，确定m个所述关键地质因子对应的值为输入样本数据，m个所述关键触发因子对应的值为输出样本数据；
步骤S20：初始化：
获取终止符集；
获取函数集；
获取遗传规划各参数：群体规模n、交换概率α、突变概率β、最大突变深度γ、最大迭代次数T、模型误差δ、最大二叉树深度ζ、及终止准则，所述终止准则为：达到最大迭代次数T或不大于所述模型误差δ；
根据所述终止符集和所述函数集，采用完全法随机产生一个具有所述群体规模的二叉树群体，所述二叉树群体中的每一个个体均为一个不大于所述最大二叉树深度的二叉树；
步骤S30：确定二叉树群体内每一个二叉树的适应度，判断是否满足所述终止准则？
若是，则终止迭代，且转到步骤S40；
若否，则执行以下操作：
采用竞技选择法，从二叉树群体中，确定待执行复制操作的父代二叉树并执行复制操作；
根据所述交换概率α，从二叉树群体中，选择至少一个待执行交换操作的二叉树组，每个二叉树组内包括两个待执行交换操作的二叉树；
在每个待执行交换操作的二叉树组中，采用均匀分布的随机方法选取待执行交换操作的叶节点或中间节点并执行交换操作；
根据所述突变概率β及最大突变深度γ，从二叉树群体中，选择待执行突变的至少一个二叉树并对所述待执行突变的二叉树的叶节点或中间节点执行突变操作；
以增量为1的方式递增当前的迭代次数；
S40：若满足最大迭代次数而终止迭代，则确定当前群体中具有最大的适应度的二叉树为所述多区域泥石流触发二叉树；
若满足迭代误差不大于所述模型误差δ而终止迭代，则确定所述迭代误差不大于所述模型误差δ的二叉树为所述多区域泥石流触发二叉树。
需要说明的是，在每次迭代中，通过复制操作、交换操作和突变操作，产生了对应于当前迭代次数的新的二叉树群体。
需要说明的是，每一个二叉树均代表一个由终止符集中的元素为自变量，由函数集中的元素为运算符的数学表达式。
需要说明的是，由完全法产生的初始个体，每一个叶节点的深度都等于给定的最大深度。二叉树的最大深度(也可以认为是层数)通常为10 以内，优选地取4～6。
所谓叶节点的深度，是指叶节点距树根的层数。用完全法产生初始群体时，首先从函数集中选取元素做根节点，然后根据给定的最大深度分别从函数集中及终止符集中选取元素。假若待定结点的深度小于给定的最大深度，则该节点从函数集中随机选取，以便算法树延伸。若待定节点的深度等于给定的最大深度，则该节点从终止集中选择，图3示出了一个二叉树。
具体地，在执行遗传操作(包括复制(即选择)、交换、和突变)，时，复制是将上一代优良个体未经任何改变地直接进入下一代群体，以体现进化过程中“优胜劣汰”原则。
具体地，采用“竞技选择法”进行个体复制，这种方法主要是每次从群体中随机选取一组包含k个个体的小群体，然后从中挑出适应度最高的个体作为复制的对象。被选中的个体经复制之后，仍保留在父代个体中，并且有可能再次进入随后迭代中随机选取的k个个体的小群体中。通过调整k值的大小，可以改变选择的力度。k值越大，被选出的个体越优，反之，被选出的个体只是相对较优而已。
交换操作，则是对选中的两棵二叉树中的中间节点或叶节点的位置及数值进行互换。
通常突变概率为小于0.1的数值，优选地，为0.05。突变概率过大则随机性太强，不利于解空间的搜索和解的收敛。突变在遗传规划中的作用比较小，通常，每代只有5％的个体发生突变，突变概率非常小。
交换是将两个父代个体的部分组份相互交换代替，形成两个新的个体，其方法步骤如下：
①从父代群体中选择两个优良个体作为交换对象，其选择方法同复制方法中任何一种择优选择法。
②在两个父代群体中用均匀分布的随机方法分别选取交换点(包括交换点以下的子树)成为交换段。
③交换两个父代个体的交换段，形成两个新的子代个体。
突变操作，通常将根据突变概率选取的二叉树个体中的一个节点作为突变点；删除突变点以下的子树或用随机产生的一棵新的二叉树更换突变点处及其下的子树。
根据需要，可以灵活选择在交换操作和突变操作中使用的二叉树群体，也即，该二叉树群体可以是在本次迭代中新产生的，也可以是上一次迭代中产生的。
复制操作不改变二叉树群体的个数；交换操作和突变操作也不改变二叉树群体的个数。
具体地，所述的方法，在所述步骤S20中，
所述终止符集中的元素，包括：流域面积、沟床平均坡度、松散物质比率及至少一个随机数，其中，每个随机数在一个取值范围内均匀取值；
所述函数集中的元素，包括：四则运算符、三角函数、ex、及求平方根。
经过多次尝试，并最终选择上述函数集中的元素，从而保证能在解空间内进行有效地搜索，又不会跳过最优解。
具体地，所述的方法中，所述根据所述终止符集和所述函数集，采用完全法随机产生初始群体的步骤，包括：
从所述函数集中选取元素做二叉树的根节点；
根据所述最大二叉树深度ζ，依次从所述函数集或所述终止符集中选取元素作为待定节点：
在待定节点的深度小于所述最大二叉树深度时，从所述函数集中随机选取元素作为所述待定节点；
在待定节点的深度等于所述最大二叉树深度时，从所述终止符集中选取元素作为所述待定节点。
具体地，所述的方法，在步骤S30中，
所述确定每一个二叉树的适应度的步骤，具体为：采用第一公式计算每一个二叉树的迭代误差，确定所述迭代误差的相反数为所述二叉树的适应度，所述第一公式为：
其中，pi为第i个输入样本数据对应于所述二叉树的预测输出， 1≥i≥m,oi为所述第i个输入样本数据对应的第i个输出样本数据。
将输入样本数据的值作为待评价的二叉树对应的数学表达式中自变量的值，则该数学表达式的输出即为预测得到的输出数据。
该输出数据与该输入样本数据对应的输出样本数据之间的误差，即为预测误差。该误差越大，则预测精度越低；也即该二叉树作为遗传中的一个个体，其适应度越差。
具体地，所述迭代误差的相反数或倒数均可以作为所述二叉树的适应度。这里优选采用相反数，可以避免在迭代中不同数值的量级相差太多，引起更大的计算误差。
与传统方法相比，本发明提供的泥石流预测方法利用基于二叉树的遗传规划来识别预测目标地区在触发泥石流时的临界降雨指数；基于训练样本，通过选择(复制)、交换和突变等操作，逐步迭代而逼近最优解；且不需要事先假定数据模型的种类。
特别地，与人工神经网络逼近最优解的方法相比，遗传规划能够输出各因素之间确定的非线性表达式，从而可以实现因变量和自变量关系的可视化。
本发明提供的泥石流预警方法确定了不同沟谷泥石流降雨指数与其流域面积、松散物质率以及沟床平均坡度之间的非线性函数关系，实现了不同地质环境背景下的泥石流预警。
在具体实施时，针对指定的区域，通过收集该区域的泥石流沟的流域面积、松散物质比率以及沟床平均坡度这些地质数据，就可以建立多区域泥石流触发二叉树；将最新观测的地质数据和降雨数据输入该多区域泥石流触发二叉树，就可以预测与之对应的泥石流爆发临界降雨指数，实现泥石流预警。
在确定了预测模型后，将测试样本带入模型，获取预测的临界降雨指数，并与实际降雨指数作比较，计算模型的预测精度。
当模型预测精度满足要求时，可将其他地区泥石流相关地质因子带入该模型，预测泥石流爆发临界降雨指数，实现不同地质环境背景下的泥石流预警。
以下以具体的实施例来详细说明本发明提供的泥石流预警方法。
根据泥石流发生地地质环境背景，在充分搜集历史泥石流灾害事件的基础上，综合分析确定泥石流爆发的主要影响因子，如流域面积、沟床平均坡度、松散物质比率以及小时雨强、当日雨量、前期有效雨量等。采用遥感解译，地面调查以及历史雨量站数据搜集等方式，并对主要因子进行量化。
如图3所示，通过搜集和建立历史上已经发生泥石流的灾害史，建立如表1所示泥石流灾害综合因子样本库；其中前31项样本作为训练样本，后5 个样本(第32-36个样本)作为测试样本对预测模型的可靠性进行验证。
表1泥石流爆发综合因子一览表
利用Matlab中的遗传规划工具箱进行模型训练及预测。设定遗传规划模型的参数如表2所列；选择表1中前31项的数据作为遗传规划模型的训练样本。经过二叉树学习和训练并满足终止准则后，最终得到的泥石流爆发降雨指数预测模型如图2中的二叉树所示，其中，变量X1、X2和X3分别代表区域流域面积、松散物质比率以及区域沟床平均坡度。
表2遗传规划的基本参数一览表
需要说明的是，图2中plus为"+"、times为"×"、minus为"-"、 mysqrt为mylog为exp为“ex”。
将第32-36个样本作为测试样本数据，将各测试样本中的各参数作为该泥石流爆发降雨指数预测模型的输入数据，运行该预测模型，即可得到各测试样本对应的降雨指数RTI预测值，具体如表3所列。
计算各预测值与泥石流灾害综合因子样本库中对应的泥石流爆发降雨指数的实测值的相对误差，如表3所列，全部的5个相对误差均在8％之内，满足预测精度要求。
表3基于遗传规划的泥石流爆发临界降雨指数预测误差一览表
最后应说明的是：以上所述的各实施例仅用于说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或全部技术特征进行等同替换；而这些修改或替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。