CN108304668B - 一种结合水文过程数据和历史先验数据的洪水预测方法 - Google Patents

Abstract

本发明公开了一种结合水文过程数据和历史先验数据的洪水预测方法，利用加权聚类算法对影响降雨信息进行聚类分析，近似含水量信息；对数据进行清理，补充缺失的数据，修正错误数据；根据流域特征选取通用水文模型，获取水文过程中的变量之间的因果推理关系；对获取的一组组推理关系，采用有向图的形式进行表示，其中节点表示物理变量，边的方向表示因果推理关系，并以此拓扑结构构建贝叶斯网络；利用历史数据，结合通用水文模型，模拟产生中间结果数据，共同组成样本集；根据样本集获取历史先验知识，计算条件概率表，以实现对网络的训练，完成水文过程知识和历史先验知识的综合建模；输入测试数据，实现洪水预测。

Description

一种结合水文过程数据和历史先验数据的洪水预测方法

技术领域

本发明涉及一种结合水文过程数据和历史先验数据的洪水预测方法，属于洪水预报技术领域。

背景技术

山洪是山丘区中小河流由降水引发的突发性、暴涨暴落的洪水，而我国中小河流众多，各个流域具有地理地貌、气候的独特性，缺少统一便捷的建模手段。近年来，由于极端天气时间频多，地区持续性强降雨、突发性暴雨时常发生，而山丘地区下垫面呈现高陡坡、河流源短流急，在暴雨天气下极易发生山洪等山地灾害，对人民生命财产造成威胁。因此基于中小河流的洪水预报得到国家的高度重视。

在防洪领域中，常用的模型主要分基于物理(过程)驱动和基于数据驱动。物理驱动模型是基于水文概念模型的，该方法分析造成洪水的水文过程，建立适合某一流域的通用模型，模型参数需要依赖经验和人工交互不断迭代才能完成率定。模型预报的精度取决于建模者的知识和经验以及数据资料的完备情况。这类模型内包含的参数具备物理意义，具有良好的解释性，但中小河流复杂的环境因素加大了参数率定的难度。近年来，线性回归、神经网络、支持向量机等数据驱动模型被运用在洪水预报中，这类模型将水文过程视作黑箱子，不考虑系统内部的物理机制，通过建立输入输出样本间的映射关系以实现建模。由于无法获取精确的水文要素(降雨和流量)，也无法获取流域内水循环要素的可靠时空分布，造成智能模型在预报性能上有一定的局限性。水文预报中，各个因素存在不确定性，因素之间相互关联，这使得物理驱动模型和数据驱动模型均有各自的不足。

发明内容

发明目的：针对现有技术中存在的问题，本发明提供一种结合水文过程数据和历史先验数据的洪水预测方法。该方法选定适合的通用水文模型，分解水文过程构建贝叶斯网络，通过通用水文模型产生中间变量结果以丰富样本集，并用统计学方法获得变量之间的先验知识，完成水文过程的知识固化。本方法选用通用水文模型，避免了精确参数的率定，但同时利用贝叶斯网络的特点，能对不确定知识有良好的表达。

技术方案：一种结合水文过程数据和历史先验数据的洪水预测方法，包括如下步骤：

第一步：收集水文要素(历史降雨量数据、历史流量数据、历史蒸发量数据等)，根据降雨数据加权聚类近似推算土壤含水量，构成初始数据集。对数据进行清洗，对缺失数据进行补充；

第二步：根据流域特征选取适合该流域预报的水文模型，分析水文过程中的变量及变量之间的因果关系，收集水文过程知识，对水文过程知识进行建模，采用有向图的形式来表示水文过程，其中节点表示物理变量，边的方向表示因果推理关系，并以此拓扑结构构建贝叶斯网络；

第三步：利用历史数据，运行通用水文模型，模拟产生中间结果，历史水文要素数据和中间结果数据共同构成数据样本集，对样本集进行数理统计分析，分别计算每个变量的概率分布和变量之间的条件概率分布，以此获取历史先验知识，计算贝叶斯网络的每个节点的条件概率表，以实现对网络的训练，完成水文过程知识和历史先验知识的综合建模；

第四步：完成知识固化的贝叶斯网络可以替代水文模型，再输入新的水文要素，可以完成水文要素的预测。

第一步中加权聚类近似产生土壤含水量的步骤如下：

步骤1.1：收集整理历史降雨量；

步骤1.2：选取降雨特征，对土壤含水量计算点过去15天的降雨数据进行整理分析，抽取出多个时段的累计降雨量作为特征；

步骤1.3：根据经验对多个特征初步赋权，再进行k-means聚类，获得多个聚类质心(聚类质心为n位向量，a_i，i＝1,2,3…n)，计算相对湿润度A，计算公式为

其中w_i表示聚类质心的每个维度的权值；

步骤1.4：对相对湿润度A进行排序，对聚类质心进行定级，评定土壤的含水百分比；

步骤1.5：在进行新一轮土壤含水量预测时，根据已有的降水特征，计算与各个聚类质心的距离，选取较近的质心评级近似该土壤含水量。

第二步中获取水文过程知识的步骤如下：

步骤2.1：综合分析目标流域的水文特性、气候、土壤特性等，选定适合的通用模型；

步骤2.2：在选定通用模型之后，分析水文过程中的变量及变量之间的因果关系，收集水文过程知识(洪水的物理过程中各个物理量之间的因果关系)，具体的变量间因果关系如下所示：

(E_t,P_t,WM_t)→(R_t+1,WM_t+1)

(S_t,R_t+1)→(S_t+1,RS_t+1,RI_t+1,RG_t+1)

(TRS_t+1,TRI_t+1,TRG_t+1)→(TR_t+1)

(Q_t-1,Q_t,TR_t+1)→(Q_t+1,…,Q_t+h)

其中下标t表示时刻，E表示蒸发量，P表示降雨量，WM表示土壤含水量(张力水)，R表示产流量，S表示自由水含量，RS表示地表径流，RI表示壤中流，RG表示地下径流，TRS表示全流域的地表径流出流量，TRI表示全流域的壤中流出流量，TRG表示全流域的地下径流出流量，TR表示全流域出流，Q表示流量。

对水文过程知识进行建模，采用有向图的形式来表示水文过程，其中节点表示物理变量，边的方向表示因果推理关系；

步骤2.3：根据已有的网络拓扑图，构建贝叶斯网络，网络为有向图，有向边起点变量作为有向边终点变量的条件，节点保存变量的条件概率表。

第三步中利用贝叶斯网络以实现水文过程知识和历史先验知识的综合建模步骤如下：

步骤3.1：利用历史数据和通用的水文模型产生中间物理量的取值，需要产生的变量有R，RS，RI，RG，TRS，TRG，TRI，TRG，TR等，历史水文要素数据和中间结果数据共同构成数据集，对数据集进行离散化处理和规范化处理；R表示产流量，RS表示地表径流，RI表示壤中流，RG表示地下径流，TRS表示全流域的地表径流出流量，TRI表示全流域的壤中流出流量，TRG表示全流域的地下径流出流量，TR表示全流域出流。

步骤3.2：随机对数据集进行分组，其中2/3数据集构成训练集，1/3数据集构成测试集；

步骤3.3：针对网络中任意节点Y，其前置节点X₁,X₂,…,X_k，通过统计分析求解条件概率分布P(Y|X₁,X₂,…,X_k)，完成对节点参数的训练；为增加网络的自适应性和泛化能力，可以根据一定的规则对数据集进行扩充，如假定前置节点X_i符合某种分布，根据分布特性产生数据对数据集进行扩充。

第四步中利用综合模型进行预报的步骤如下：

步骤4.1：输入t时刻的水文要素X₁,X₂,…X_n的具体取值X₁＝x_t,1,X₂＝x_t,2,…X_n＝x_t,n；

步骤4.2：计算argmax_yP(Y＝y|X₁＝x_t,1,X₂＝x_t,2,…,X_n＝x_t,n)，此时y作为网络的输出，即径流取值；

步骤4.3：步骤4.2中同时也可以获得中间变量的节点，作为t+1时刻的输入，继续进行推理预测，完成结果的连续预测。

有益效果：与现有技术相比，本发明提供的结合水文过程数据和历史先验数据的洪水预测方法，对影响降雨信息聚类获得难以实际测量的土壤含水量信息，再通过分析通用水文模型，提取物理量之间的因果关系，建立贝叶斯网络，通过对历史样本的统计分析，实现对网络的训练，以实现水文过程知识的固化，通过贝叶斯网络对洪水进行预报。

附图说明

图1为本发明的方法流程图；

图2为是实施例中综合模型和通用水文模型对比图，(a)、(b)、(c)、(d)分别是预见期1～4小时的结果展示；

图3为模型拓扑结构。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，结合水文过程数据和历史先验数据的洪水预测方法，主要包括以下几步：

一是利用加权聚类算法对影响降雨信息量进行聚类分析，根据以往15天的降雨情况来表征土地的干旱程度，以此来近似含水量信息；二是对数据进行清理，补充缺失的数据，修正错误数据，采用近邻数据的均值来平滑缺失值或修正错误值；三是根据流域特征选取适合的通用水文模型，分析水文过程中的变量及变量之间的因果关系，收集水文过程知识；四是对水文过程知识进行建模，采用有向图的形式来表示水文过程，其中节点表示物理变量，边的方向表示因果推理关系，并以此拓扑结构构建贝叶斯网络；五是利用历史数据，结合通用水文模型，模拟产生中间结果数据，加入到样本集；六是对样本集进行统计分析，获取历史先验知识，计算条件概率表，以实现对网络的训练，完成水文过程知识和历史先验知识的综合建模；七是输入测试数据，实现洪水预测。

由于中小河流所处地理环境复杂，土壤类型多样，以往测量土壤含水量方法效果不佳，再加上部分地区并没有土壤含水量的资料信息，因此本文提出一种用影响降雨信息来近似推算土壤含水情况的方法。以下是加权聚类计算土壤含水量的算法步骤：

步骤1.1：收集整理历史降雨量；

步骤1.2：选取降雨特征，对土壤含水量计算点过去15天的降雨数据进行整理分析，抽取出多个时段的累计降雨量作为特征；

步骤1.3：根据经验对多个特征初步赋权，再进行k-means聚类，获得多个聚类质心(聚类质心为n位向量，a_i，i＝1,2,3…n)，计算相对湿润度A，计算公式为

其中w_i表示聚类质心的每个维度的权值；

步骤1.4：对相对湿润度A进行排序并定级，以此实现对聚类质心的评定(评定土壤的含水百分比)；

步骤1.5：在进行土壤含水量推算时，根据已有的降水特征，计算与各个聚类质心的距离，选取较近的质心评级近似该土壤含水量。

物理驱动模型在传统水文预报中有很多发展，因此也衍生了很多种适合不同流域特性的水文模型，水文模型中蕴含了很多水文过程的知识，对水文模型分析可以提炼水文过程知识，并采用贝叶斯网络进行表达，其具体算法步骤如下：

步骤2.1：综合分析目标流域的水文特性、气候、土壤特性等，选定适合的通用水文模型；

步骤2.2：在选定通用模型之后，分析水文过程中的变量及变量之间的因果关系，收集水文过程知识，具体如下所示：

(E_t,P_t,WM_t)→(R_t+1,WM_t+1)

(S_t,R_t+1)→(S_t+1,RS_t+1,RI_t+1,RG_t+1)

(TRS_t+1,TRI_t+1,TRG_t+1)→(TR_t+1)

(Q_t-1,Q_t,TR_t+1)→(Q_t+1,…,Q_t+h)

其中E表示蒸发量，P表示降雨量，WM表示土壤含水量(张力水)，R表示产流量，S表示自由水含量，RS表示地表径流，RI表示壤中流，RG表示地下径流，TRS表示全流域的地表径流，TRI表示全流域的壤中流出流量，TRG表示全流域的地下径流出流量，TR表示全流域出流，Q表示流量。

对水文过程知识进行建模，采用有向图的形式来表示水文过程，其中节点表示物理变量，边的方向表示因果推理关系；

步骤2.3：根据已有的网络拓扑图，构建贝叶斯网络，网络为有向图，有向边起点变量作为有向边终点变量的条件，节点保存变量的条件概率表。

完成对水文过程知识的获取之后，需要采用历史先验知识对该网络进行训练。由于网络中部分中间节点并没有数据，为解决这个问题，利用通用水文模型，基于历史水文要素数据，产生中间物理量结果数据，加入数据集，再基于数理统计分析，获取历史先验知识，对网络进行训练，其具体的算法步骤如下：

步骤3.1：利用历史数据和通用的水文模型产生中间物理量的取值，需要产生的变量有R，RS，RI，RG，TRS，TRG，TRI，TRG，TR等，历史水文要素数据和中间结果数据共同构成数据集，对数据集进行离散化处理和规范化处理；

步骤3.2：随机对数据集进行分组，其中2/3数据集构成训练集，1/3数据集构成测试集；

步骤3.3：针对网络中任意节点Y，其前置节点X₁,X₂,…,X_k，通过统计分析求解条件概率分布P(Y|X₁,X₂,…,X_k)，完成对节点参数的训练；为增加网络的自适应性和泛化能力，可以根据一定的规则对数据集进行扩充，如假定前置节点X_i符合某种分布，根据分布特性产生数据对数据集进行扩充。

在综合水文过程知识和历史先验知识建模后，可以用该模型进行水文预报，其具体的算法如下：

步骤4.1：输入t时刻的水文要素X₁,X₂,…X_n的具体取值X₁＝x_t,1,X₂＝x_t,2,…X_n＝x_t,n；

步骤4.2：计算argmax_yP(Y＝y|X₁＝x_t,1,X₂＝x_t,2,…,X_n＝x_t,n)，此时y作为网络的输出，即径流取值；

步骤4.3：步骤4.2中同时也可以获得中间变量的节点，作为t+1时刻的输入，继续进行推理预测，完成结果的连续预测。

为了验证本发明的效果，选取长湖汛期数据作为研究对象，昌化流域位于浙江省分水江流域上游，地势西北高东南高，属浙西山丘区，典型的中小河流水系。选取1998-2010年汛期场次洪水数据，数据时间间隔为1小时，其中1998年-2009年攻6790个数据作为训练样本，2010共671个数据作为测试样本。首先通过对雨量信息进行聚类，对土壤含水量进行近似估算，然后选定合适昌化流域的通用水文模型进行分析，建立贝叶斯网络模型，对样本数据进行数理统计分析，计算条件概率表，完成对网络的训练，最终完成水文过程的专家知识的固化，并投入预报。如图2所示，本文通过对比固化模型和通用水文模型，其确定性在不同预见期的对比结果如下：

预见期	固化模型	通用模型
			1h	0.939455	0.71
2h	0.981491	0.67
			3h	0.801938	0.67
4h	0.605704	0.66

实验结果表明：固化模型能对水文过程的专家知识进行良好的固化，并且具有良好的预测能力，这是因为固化模型不仅基于传统的通用水文模型，还借鉴了数据模型的特点，把严谨的数理推导转化成因果概率推导，弱化了计算性，降低了传统水文模型的应用难度。

Claims (8)

1.一种结合水文过程数据和历史先验数据的洪水预测方法，其特征在于，包括如下步骤：

步骤1，收集水文要素，根据降雨数据加权聚类近似推算土壤含水量，构成初始数据集，对数据进行清洗，对缺失数据进行补充；

步骤2，根据流域特征选取通用水文模型，通过分析通用水文模型，获取水文过程中的变量之间的因果推理关系，对每组因果推理关系进行建模，采用有向图的形式进行表示，其中节点表示物理变量，边的方向表示因果推理关系，并以此拓扑结构构建贝叶斯网络；

步骤3，利用所述初始数据集中的历史数据，运行通用水文模型，模拟产生中间结果，历史水文要素数据和中间结果数据共同构成数据样本集，对样本集进行统计分析，获取历史先验知识，计算贝叶斯网络的每个节点的条件概率表，以实现对网络的训练，完成水文过程知识和历史先验知识的综合建模；

步骤4，对于综合建模的模型，输入新的水文要素，完成水文要素的预测。

2.如权利要求1所述的结合水文过程数据和历史先验数据的洪水预测方法，其特征在于，步骤1中加权聚类近似产生土壤含水量的步骤如下：

步骤1.1：收集整理历史降雨量；

步骤1.2：选取降雨特征，抽取出多个时段的累计降雨量作为特征；

步骤1.3：对多个特征初步赋权，再进行k-means聚类，获得多个聚类质心，计算相对湿润度A；

步骤1.4：对相对湿润度A进行排序，对聚类质心进行定级，评定土壤的含水百分比；

步骤1.5：在进行土壤含水量预测时，根据已有的降水特征，计算与各个聚类质心的距离，选取较近的质心评级近似该土壤含水量。

3.如权利要求2所述的结合水文过程数据和历史先验数据的洪水预测方法，其特征在于，相对湿润度

其中w_i表示第i个聚类质心a_i的维度的权值，其中n表示聚类质心的维数。

4.如权利要求1所述的结合水文过程数据和历史先验数据的洪水预测方法，其特征在于，步骤2中获取水文过程知识的步骤如下：

步骤2.1：选定通用模型；

步骤2.2：在选定通用模型之后，分析水文过程中的变量及变量之间的因果关系，收集水文过程知识，对水文过程知识进行建模，采用有向图的形式来表示水文过程，其中节点表示物理变量，边的方向表示因果推理关系；

步骤2.3：根据已有的网络拓扑图，构建贝叶斯网络，网络为有向图，有向边起点变量作为有向边终点变量的条件，节点保存变量的条件概率表。

5.如权利要求4所述的结合水文过程数据和历史先验数据的洪水预测方法，其特征在于，收集水文过程知识，具体如下所示：

(E_t,P_t,WM_t)→(R_t+1,WM_t+1)

(S_t,R_t+1)→(S_t+1,RS_t+1,RI_t+1,RG_t+1)

(TRS_t+1,TRI_t+1,TRG_t+1)→(TR_t+1)

(Q_t-1,Q_t,TR_t+1)→(Q_t+1,…,Q_t+h)

其中，下标t表示时刻，h表示最大观测步长，E表示蒸发量，P表示降雨量，WM表示土壤含水量，R表示产流量，S表示自由水含量，RS表示地表径流，RI表示壤中流，RG表示地下径流，TRS表示全流域的地表径流出流量，TRI表示全流域的壤中流出流量，TRG表示全流域的地下径流出流量，TR表示全流域出流，Q表示流量。

6.如权利要求1所述的结合水文过程数据和历史先验数据的洪水预测方法，其特征在于，步骤3中利用贝叶斯网络以实现水文过程知识和历史先验知识的综合建模步骤如下：

步骤3.1：利用历史数据和通用的水文模型产生中间物理量的取值，需要产生的变量有R，RS，RI，RG，TRS，TRG，TRI，TR，历史水文要素数据和中间结果数据共同构成数据集，对数据集进行离散化处理和规范化处理；R表示产流量，RS表示地表径流，RI表示壤中流，RG表示地下径流，TRS表示全流域的地表径流出流量，TRI表示全流域的壤中流出流量，TRG表示全流域的地下径流出流量，TR表示全流域出流；

步骤3.2：随机对数据集进行分组，一组为训练集，一组为测试集；

步骤3.3：针对网络中任意节点Y，其前置节点X₁,X₂,…,X_k，通过统计分析求解条件概率分布P(Y|X₁,X₂,…,X_k)，完成对节点参数的训练。

7.如权利要求6所述的结合水文过程数据和历史先验数据的洪水预测方法，其特征在于，步骤3.3中为增加网络的自适应性和泛化能力，根据前置节点X_i符合的规则对数据集进行扩充。

8.如权利要求1所述的结合水文过程数据和历史先验数据的洪水预测方法，其特征在于，步骤4中利用综合模型进行预报的步骤如下：

步骤4.1：输入t时刻的水文要素X_1,t,X_2,t,…X_n,t的具体取值X_1,t＝x_t,1,X_2,t＝x_t,2,…X_n,t＝x_t,n；其中，x_t,i表示第i个观测点在时刻t具体的取值；

步骤4.2：计算argmax_yP(Y＝y|X_1,t＝x_t,1,X_2,t＝x_t,2,…,X_n,t＝x_t,n)，此时y作为网络的输出，即径流取值；

步骤4.3：步骤4.2中同时也可以获得中间变量的节点，作为t+1时刻的输入，继续进行推理预测，完成结果的连续预测。

Images