一种基于Budyko理论的水文随机插值径流预测方法

摘要

本发明公开了一种基于Budyko理论的水文随机插值径流预测方法，包括以下步骤：确定待预测流域，所述流域分为多个子流域；利用Budyko理论对每个子流域多年平均年径流深进行空间预测，得到年径流预测值，作为径流的确定性成分；利用水文随机插值方法对多个子流域多年平均年径流预测偏差进行插值，作为径流的随机性成分；利用所述确定性成分与随机性成分耦合，得到最终预测径流。本发明用于预测无资料地区的径流深并解释径流的空间分布特征，提高了径流预测精度，为径流区域化提供了一种新的思路和方法。

说明书

技术领域

本发明涉及一种基于Budyko理论的水文随机插值径流预测方法，尤其是涉及无资料地区的径流预测、径流空间分布特征及规律分析领域，属于水文学及水资源技术领域。

背景技术

采用空间插值技术的径流图化研究可用于揭示流域的水资源空间分布特征、预测无资料地区的多年平均年径流，这是水文学的基本问题之一，也是水资源规划管理的研究热点。与基于点观测的降水和潜在蒸散发不同，流域出口的径流观测受到整个流域降水和下垫面特征的共同作用，而非仅代表某点的径流特征，因而径流的空间插值不仅受实测径流值控制，也与河流水系的区域水文特性有关。

从水文学的确定性角度来说，流域的径流过程是区域水文循环的重要组成部分，分布式降雨径流模型能模拟径流对降水的响应，模型的径流模拟结果也解释了其区域空间分布规律。然而这种方法很大程度上依赖于数据集的可用性。考虑到水文气象要素和下垫面特征的径流区域化的插值方法一定程度上可弥补分布式模型的不足，再如建立径流和气象、土地利用和地形等控制要素的经验关系来进行空间插值。但经验关系仅适用于特定流域并不能不经验证而直接应用于其他类似流域。Budyko理论作为一种相对简单的半经验性公式，常用来估计多年尺度下的陆面的水量平衡问题，定量分析区域降水转化为蒸散发和径流的比例，预测无资料地区的可利用水量和径流变化特征。然而Budyko理论能较好描述大流域尺度上的水文-气象特征，对特定流域或地区的径流或蒸发预测还存在较大偏差。为提高Budyko区域水文预测方法的精度，很多研究尝试建立曲线参数与土地利用方式/植被、地形、气候变异性等的关系，然而尚未将Budyko预测值与流域实测径流的偏差作为区域统计特性纳入考虑。

地统计插值方法(克里金插值)基于区域变量的统计特征相似性，其特定位置的变量估计值是所有样本区域变量值的加权组合。传统的地统计插值方法的变量空间相似性由基于空间欧氏距离的点对之间的变异函数来描述。克里金插值作为一种最优线性无偏估计方法，使预测值和实测值的方差最小。根据对确定性趋势值的不同处理，克里金插值可以分为不同的类型。如普通克里金插值(Ordinary Kriging)以未知常量作为确定性趋势值，泛克里金插值(Universal Kriging)以空间坐标的回归方程作为确定性趋势，块克里金(block-kriging)不同于普通克里金插值等常规克里金方法，估计的是块状区域变量值而非针对点估值，因而点与点的变异函数由点与方块变异函数代替。

随机插值方法是水文气象领域区域变量的常用插值方法，尤其是降水、蒸发这样的空间点过程。如将该类方法应用于径流插值那就意味着忽略了流域的横向径流形成过程，将空间连续的径流过程作为空间点特性，因而丢失了径流的空间变化特征，可能会导致径流的过高估计。

沿河流流动的水流受嵌套水系结构的影响，是空间连续的整体过程，有嵌套关系的流域上下游集水区域与相邻的无水力关系的流域插值过程有所差异，因而传统的地统计插值方法不适用于径流插值过程。也就是，传统插值过程中采用的欧式距离多数情况下不能很好的适用于径流插值过程中的嵌套流域的空间距离的定义。考虑到水系结构的影响，Gottschalk开发了一套径流水文随机插值方法，同时考虑了径流的空间联系整体特性以及流域分级水系结构的共同作用，重新定义了流域空间距离的计算，点对之间的变异函数由流域的变异函数代替，并增加了水量平衡约束。但是该统计方法并未考虑到区域多年平均年径流受气候和下垫面影响的确定性趋势组分，同时前人的工作也表明同时考虑区域变量的随机和确定性规律，并将其同时应用于水文气象地统计领域有较好的适用性。

发明内容

发明目的：为了解决现有技术存在的问题，预测无资料地区的径流深并分析径流的空间分布特征，提高径流预测精度，本发明提供一种基于Budyko理论的水文随机插值径流预测方法。

技术方案：一种基于Budyko理论的水文随机插值径流预测方法，包括以下步骤：

(1)确定待预测子流域及多个已知子流域；通过水文站观测每个已知子流域实际多年平均年径流深R(x)的；利用Budyko理论对每个已知子流域多年平均年径流进行空间预测，得到多年平均年径流深预测值R_d(x)，作为待预测子流域径流的确定性成分；

(2)将R(x)与R_d(x)的差值作为多年平均年径流预测偏差；利用水文随机插值方法对多个已知子流域多年平均年径流预测偏差进行插值，得到已知子流域径流深偏差插值结果R_s^*(x)，作为待预测子流域径流的随机性成分；

(3)利用所述确定性成分与随机性成分耦合，得到待预测子流域的最终预测径流：

R^*(x)＝R_d(x)+R_s^*(x)(1)

式中，R^*(x)表示待预测子流域x处的最终预测径流。

优选的，步骤(1)中，分别统计多年平均降水量、实际蒸散发量、潜在蒸散发量，利用Budyko理论对多年平均年径流进行空间预测：

P-E＝R_d(x)(3)

由公式(2)和(3)得公式(4)：

其中，P、E、E₀、R_d(x)分别是多年平均年降水量、实际蒸散发量、潜在蒸散发量和budyko预测的径流深(mm)；ω是公式参数，且ω取值范围是(1,∞)，ω值越大，降水转化为径流的比例越低。

优选的，步骤(2)中，利用水文随机插值方法对多年平均年径流预测偏差进行插值包括：在block kriging插值方法的基础上，考虑流域水系结构特征的流域距离，同时增加流域水量平衡约束，预测的径流深偏差作为空间上的连续过程，预测公式为：

式中，R_s^*(A₀)是待预测子流域径流深偏差，R_s(A_i)是第i个面积为A_i的已知子流域计算得到的径流深偏差，Λ是权重矩阵，R_s是径流深偏差矩阵，R_s^*(A₀)是R_s(A₀)的估计值，其中A₀是预测的单元面积，根据最优线性无偏估计要求E[R_s^*(A₀)-R_s(A₀)]＝0计算权重矩阵Λ：

C(u_i,u_j)是每对已知子流域的拟合协方差函数值，C₀(u_i,u₀)是待求的u₀位置和每个已知子流域u_i处的拟合协方差函数值，λ_i为待预测子流域x处对应的各已知子流域的权重值，μ为拉格朗日系数，计算得到权重矩阵后代入式(3)得到径流深偏差的预测值。

优选的，权重矩阵Λ的计算方法：

将式(4)表示为矩阵形式：CΛ＝C₀，因而权重矩阵Λ表示为：

Λ＝C^-1C₀.>

其中：

其中，矩阵C是每对已知子流域中的点的协方差函数矩阵，C₀是每对已知子流域中的点与待预测子流域中的点的协方差函数矩阵，u_i为已知子流域中的点的位置，u₀为待预测子流域中的点的位置，μ为拉格朗日系数。

优选的，式(5)、(6)是一个待预测区域的径流深偏差表达式，如果插值过程同时计算M个子流域的径流深偏差，式(5)不变，最优权重矩阵为：

且

其中L_i是第i个已知子流域的对于每个已知子流域的权重值，μ^*和μⁱ为拉格朗日系数，所述C和C₀为：

其中，矩阵K表示为：

矩阵V和G为：

其中ΔA_i是第i个非嵌套子流域的集水面积，n_i为每个非嵌套子流域ΔA_i所含的网格个数，上述矩阵方程组包括水量平衡约束，水量平衡约束即流域出口处的实测径流量应等于所有子流域待插值径流量之和，表示为：

R_T是流域出口的实测径流量。

优选的，将每个非嵌套已知子流域ΔA_i划分为n_i个面积为a的网格，每个ΔA_i的径流偏差预测值是权重系数与已知子流域计算得到径流深偏差的线性组合，将式(16)进一步表示为：

其中，n_T是基本网格的个数，r_T是流域出口的实测径流深；

由Gottschalk定义流域距离，即流域内所有网格点对的距离的期望值，表示为：

其中，A、B表示两个已知子流域；A₁和A₂分别是已知子流域A、子流域B的流域面积，u₁、u₂分别为已知子流域A、已知子流域B的位置；

根据流域距离，绘制出实验协方差函数与流域距离的关系图，基于流域支撑A₁和A₂，理论协方差函数Cov(A,B)为：

其中，Cov_p是点协方差函数。

优选的，由试错法率定最优数学模型，绘制实验协方差函数与流域距离的关系图时仅采用两两独立的已知子流域。

优选的，所述子流域至少有40个。

有益效果：本发明提供的一种基于Budyko理论的水文随机插值径流预测方法，通过增加水文随机插值外部漂移函数，将Budyko径流预测偏差应用水文随机插值方法，将两种方法耦合形成一种新的基于Budyko理论的径流水文随机插值方法，并提高了径流的估计精度。该方法可应用于水文循环中径流要素的区域化，给出特定流域内径流的空间分布特征并预测无资料地区的径流深，相比较传统的Budyko理论和一般的水文随机插值方法提高了径流的预测精度，为径流图化研究提供了一种新的思路和方法，为更好地进行水资源规划管理提供了科学依据。

附图说明

图1为本发明基于Budyko理论的水文随机插值径流预测方法的流程图；

图2为淮河流域蚌埠以上地形水系图；

图3为淮河流域蚌埠以上水文站及子流域图；

图4为淮河流域Budyko曲线E/P～E₀/P关系图；

图5为淮河子流域径流深的实验协方差和理论协方差函数；

图6为淮河子流域径流深偏差的实验协方差和理论协方差函数；

图7(a)为傅抱璞公式实测与预测径流深关系；

图7(b)为水文随机插值方法的实测～预测径流深的交叉验证结果；

图7(c)为耦合方法的交叉验证结果(图中虚线为1:1线)；

图8(a)为Budyko法预测多年平均年径流空间分布图；

图8(b)为水文随机插值方法预测多年平均年径流空间分布图；

图8(c)为基于Budyko理论的水文随机插值径流预测方法预测多年平均年径流空间分布图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。

本实施例以淮河流域蚌埠以上40个子流域为实验样本，淮河流域是中国第六大水系，位于中国南北交界的气候过渡区，从东部的暖温带变化到西部的湿润气候。流域人口密度大，是中国的主要农业区之一，每年将要消耗大量水资源以维持工农业及生活用水。人均水资源和单位面积的水资源不足全国平均水平的1/5，更重要的是超过一半的水资源被过度开采，远高于国际内陆河流的推荐利用水平(30％)。加之很高的降雨集中度，意味着全年降水有相当大的比例集中在短短几个月内，使得流域极易受到旱涝灾害的侵袭。频繁的旱涝灾害也睁大了水资源利用和防洪措施实施的难度。地形图见图2，水文站及子流域见图3，蚌埠以下淮河流域面积121,000km²，河网水系由国家地理信息系统提供的数据包得到，该流域上游共40个子流域，水文气象资料系列从1961年到2000年，其中独立的子流域为27个，及13个嵌套的子流域。该研究区北部的子流域如中牟、砖桥、告成、周口的干旱系数(E₀/P)相对较高超过了1.3，而南部各子流域较为湿润如梅山水库、横排头、潢川等的干旱系数小于0.8，多年平均年径流深约400mm，最小为北部接近黄河流域的90mm到南部山区的1000mm，淮河流域径流的时间空间变异性北部较南部更大。

本实施例针对1961年到2000年共40年的淮河流域，通过实验对比Budyko理论预测、水文随机插值方法以及本发明的基于Budyko理论的水文随机插值径流预测方法三种模型的预测结果，从而证明本方法的有益效果。如图1所示，包括四个部分，利用Budyko理论对多年平均年径流进行空间预测，水文随机插值方法预测多年平均年径流，基于Budyko理论的水文随机插值耦合方法，交叉验证评价不同模型预测结果。

1.傅抱璞公式径流预测，各流域的实际蒸散发根据多年平均年水量平衡公式(E＝P-R)求得，各已知子流域1961-2000年多年平均水量平衡要素见表1，根据公式(1)Budyko曲线E/P～E₀/P关系图由图4所示，图中两条边界线分别为干旱地区的水量平衡控制(E＝P)和湿润地区的能量平衡控制(E＝E₀)。淮河流域的Budyko曲线形状由参数ω决定，各子流域的ω值根据公式(2)求得并列于表1，ω值的范围从黄尾河流域的1.43到蒋家集流域的3.16，均值为2.32。

所有已知子流域的ω值根据非线性拟合方法(使蒸散发的Budyko预测值和多年平均水量平衡得到的实测值二者的平均标准误差最小即MAE最小)求得的全流域ω值如图4所示，该拟合的ω值为2.213，同40个流域的ω均值十分接近。将ω＝2.213代入公式2，得到：

该公式所表达的即为径流的确定性成分，干旱系数(E₀/P)越小，相应的流域的径流深越大，图4中的北部较大的径流深说明了这些子流域相较其他流域更为干旱。多年平均年降水、实际/潜在蒸散发及公式(20)表示的Budyko预测径流深以及预测偏差、平均标准误差等均在表1中给出，对预测结果的评价由表2给出。Budyko预测径流深的平均标准误差MAE为94mm，均方根误差RMSE为112mm，其中黄尾河子流域的绝对误差最大为328.03mm，息县子流域的误差最小为23.77mm，相对误差最大的为新郑子流域91mm，约站流域实测径流深的81.6％，最小为响洪甸水库的36.94mm，占实测径流深的4.99％。

表1淮河各子流域水文气象要素及预测径流深

表2三种方法在淮河流域的径流预测结果

2.径流深的水文随机插值方法，作为对比，将水文随机方法直接应用于淮河各子流域的径流预测，实验协方差由下式计算：

其中，表示所有已知子流域的实测径流深(mm)，d是两两已知子流域之间的地统计距离。为求出不同已知子流域间的地统计距离，淮河流域被划分为40行×50列的网格，由公式(18)所有已知子流域对(即已知子流域A和B，共计820对)的地统计距离由已知子流域A和B所含的两两网格之间距离的平均值得到，根据求出的两两已知子流域之间的距离和相应的各流域的实测径流深，由公式(21)求出实验协方差C(d)，进而绘出流域距离和实验协方差的关系图，该步骤中采用的是两两独立的流域的实验协方差函数和相应的流域距离并于每50km间距求得实验协方差函数的平均值绘于图中，如图5所示，实验协方差函数的最佳拟合函数为：

C(d)＝600000×exp(-d/28.62). (22)

式(22)所示拟合的指数函数用于计算公式(19)中的理论协方差函数，矩阵C,C₀,K,V及G采用matlab程序也可以相应计算得到，权重系数矩阵由此求出，进而得到淮河流域的径流插值结果，预测误差由表1给出。直接利用水文随机插值方法预测径流深的平均标准误差MAE为134mm，均方根误差RMSE176mm，最大绝对误差是黄尾河流域(448mm)，最小绝对误差是响洪甸水库流域(3mm)。由此可见，插值误差较Budyko理论的预测误差要大，说明了径流的确定性成分更大的决定了实测径流的空间分布，而这也影响了其空间插值精度。

3.基于Budyko理论的水文随机插值方法，通过去除径流空间的确定性成分利用空间插值方法对空间自相关的随机误差可以对径流预测进行校正。首先将傅抱璞公式用来表示径流的的确定性成分即外部漂移函数R_d^*(x)，径流实测值与预测值的偏差作为随机成分，利用水文随机插值方法进行预测，插值结果见表1。

两两已知子流域之间的实验协方差R_s^*(x)与相应流域距离的关系图见图6，如下指数函数作为其最佳拟合方程：

C(d)＝3000×exp(-d/48.34). (23)

通过式(23)以及公式(11)～(15)利用matlab程序计算得到矩阵C,C₀,K,V及G，然后计算得到权重系数矩阵进而求出径流偏差的预测值。该插值方法针对径流的空间偏差，将傅抱璞公式的径流预测值叠加该径流偏差的插值预测结果最终得到淮河流域的径流预测结果。该径流的耦合预测方法结果仍由表1和2给出，并针对不同方法的径流预测精度进行了对比。该方法在40个流域应用结果为平均标准误差MAE为47mm，均方根误差RMSE69mm，最大绝对误差是黄尾河流域(236mm)，最小绝对误差是蒋家集流域(1.5mm)。

4.三种方法的径流空间预测对比，如表2所示，本发明提出的基于Budyko理论的水文随机插值径流预测方法综合考虑了空间径流的确定性和随机成分，可用于系统性降低空间径流的预测误差，相较Budyko理论方法和直接利用水文随机插值两种方法，本发明提出的方法得到的径流预测平均绝对误差MAE和均方根误差RMSE都有显著降低。如黄尾河流域的最大误差由Budyko理论方法预测的328mm和直接利用水文随机插值预测的448mm降低到了236mm。表2给出的交叉验证结果决定系数由0.81和0.54提高到了0.93。

径流预测和实测值的相关关系图由图7给出，本发明方法计算绘出的二者相关关系达到0.95，与1:1线的偏离度最小。与之相比的，Budyko理论方法和直接应用水文随机插值方法得到的径流实测～预测相关关系较差(相关系数R2分别为0.58和0.82)，特别是后两种方法明显高估了径流低值，低估了径流高值(由图7可以看出)。

由上述三种预测方法(分别是Budyko理论、水文随机插值及本发明提出的耦合方法)绘出的淮河流域的径流空间分布图如图8所示，该三种方法的淮河流域径流分布图有明显差别，前两种方法与本发明方法相比明显低估了北部较干旱地区的径流深，也说明了本发明的合理性。本发明提出的方法综合考虑了径流的空间整体性特征和河网水系对径流的影响，同时综合了径流的确定性和不确定性成分，特别是将径流的随机特征加至预测方法的考虑范围，有利于提高空间径流的预测精度，也为其他水文变量的预测提供了一种新的思路和方法。