一种改善水文模型径流模拟的方法

1.本发明属于水文模型参数优化技术领域，具体涉及一种改善水文模型径流模拟的方法。


背景技术：

2.近几年，洪水和山洪事件的数量和强度不断增加，不仅造成了环境问题，还导致了巨大的人员和经济损失。其中，土壤水分是水文学中的关键状态变量，因为它控制着从土地中渗透、流失或蒸发的降雨比例，对水文预报起到重要作用。但是，由于土壤质地的空间异质性，土壤水分在空间分布上差异显著，想要得到精确区域乃至全球范围内的土壤水分，仅仅依靠地面站点的观测是远远不可能实现的。因此，土壤水分的遥感观测成为必要手段，通过遥感监测可以较为方便地获得大尺度上的土壤水分观测数据。然而遥感监测在时间和空间分辨率上同样存在缺陷，且遥感监测数据与地面测量数据相比误差较大。为此，科学家们将数据同化方法引入到水文和陆面过程研究中，以期解决土壤水分、土壤温度等关键变量的估计精度问题。
3.数据同化作为数据处理的关键技术，已日趋成为生态水文过程研究与遥感反演研究的前沿与热点问题，这项技术能够考虑模型结构及模型输入数据的不确定性，将其与遥感信息结合，引入遥感观测数据，在卫星过境时刻实时修正模型状态变量，可以充分发挥陆面过程模型和遥感数据的优势，为土壤含水量的连续模拟提供新途径。且其在协调、融合遥感技术与生态建模方面起到了重要的桥梁作用，能够有效地融合多源遥感数据以及模型的模拟结果，进而提高土壤水分的预测精度。
4.水文模型同样也是获得土壤水分条件准确估计值的一种有效方法，并已在全球范围内使用。然而，水文模拟是“不完美的”，因为它们包含不确定性，这些不确定性主要与
①
用于驱动模型的水文数据的质量和数量以及由于规模不兼容造成的代表性误差有关；
②
真实物理过程的简单表示导致不可避免的假设和简化，从而不可避免地不完美地逼近复杂的现实和
③
参数估计的误差，这可能导致模型输出的巨大误差。因此，近几年，一些科学家提出可以将数据同化技术应用到水文模型当中，实时更新水文模型中的土壤水分的状态变量，从而提高水文模型的预报精度，改善水文模型在研究区域的应用。


技术实现要素：

5.解决的技术问题：针对上述技术问题，本发明提供了一种改善水文模型径流模拟的方法，能够提高径流模拟结果的精度，从而为水文洪水预报提供帮助。
6.技术方案：一种改善水文模型径流模拟的方法，包括步骤如下：
7.s1.采用misdc模型作为水文模型，进行参数率定；
8.s2.手段1，采用直接插入法替换水文模型的土壤水分数据，将实际的土壤水分数据用于水文模型模拟径流；手段2，采用集合卡尔曼滤波器法将实际的土壤水分数据与水文模型模拟的土壤水分数据相结合，得到的状态变量值用于模型模拟径流；
9.s3.由s2中两种手段得到水文模型径流模拟的优化结果，分别比较两种优化结果的ns系数和优化前后的标准化均方根误差nrmse，采用ns系数较大者和nrmse值较小者对应的手段得到径流参数模拟的优化结果。
10.优选的，所述步骤s1中参数率定的标准通过如下计算确定：
[0011][0012]
式中，q
obs
为实测流量，q
sim
为模型模拟得到的流量，上标“—”表示平均值，当计算的err值到达最小值时的参数即为最优参数。
[0013]
优选的，所述步骤s1包括：
[0014]
s1.1通过水量平衡方程将入渗、蒸发以及排水过程联系起来求得土壤蓄水量，此过程为misdc模型中的swb模型部分；
[0015]
s1.2将土壤蓄水量数据输入misdc模型，输出径流模拟结果；
[0016]
其中，所述水量平衡方程如下：
[0017][0018]
式中，t为时间；w(t)为土壤蓄水量；w
max
为土壤最大蓄水量；f(t)为降雨入渗到土壤当中的过程；e(t)为蒸散率；g(t)为由于重力导致土壤排水的过程。
[0019]
优选的，所述入渗过程采用green-ampt方程，其表示如下：
[0020][0021]
式中，r(t)表示降雨过程；ks表示土壤饱和导水率；ψ表示湿润锋处的吸力；f表示降雨事件开始时的累积入渗深度；l表示土层的厚度。
[0022]
优选的，所述蒸发过程表示如下：
[0023][0024]
式中，etp表示潜在蒸散发。
[0025]
优选的，所述潜在蒸散发的计算方法如下：
[0026]
et
p
＝-2+b[ξ(0.46ta(t)+8.13)]
[0027]
式中，ta(t)表示气温，单位为℃；ξ表示当天日照小时数与全年的日照小时数(365
×
12h)的比值；b为待率定的参数。
[0028]
优选的，所述排水过程的计算方法如下：
[0029]
g(t)＝ks[w(t)/w
max
]
3+2/m
[0030]
式中，ks表示土壤饱和导水率；m表示土层的孔隙大小分布指数。
[0031]
优选的，所述步骤s1.2包括：
[0032]
s1.2.1雨量过剩ε的计算，计算公式如下：
[0033][0034]
式中，r为降雨强度；r为降雨事件开始时的降雨深度；λ为待率定的参数；s是连接土壤蓄水量w与雨量过剩的一个参数，其计算公式如下：
[0035]
s＝a[w
max-w(t)]
[0036]
式中，a为待率定的参数；
[0037]
s1.2.2径流量计算，运用直接径流水文过程线的方法计算流域的径流量，此过程为misdc模型中的降雨-径流模型部分；流域的平均滞后时间г的计算公式如下所示：
[0038]
γ＝η1.19a
0.33
[0039]
式中，η是待率定的参数；a是流域的面积，单位为km2。
[0040]
优选的，所述手段1的替换方法为：将实际的土壤含水量转换为实际的土壤蓄水量，并替换模型中计算得到的土壤蓄水量；其中，土壤含水量与土壤蓄水量的转换关系如下：
[0041][0042]
式中，w(t)为土壤蓄水量；θ(t)为土壤含水量；θr为土壤剩余含水量。
[0043]
优选的，所述手段2实现的具体目标函数如下：
[0044][0045]
式中，上标a和f分别表示分析和预报的过程；x
i,k+1
表示在k+1时刻，第i个成员的值；z表示观测向量；h表示观测算子；v
i,k
是均值为0，方差为rk的高斯白噪声；kk表示卡尔曼增益矩阵，其计算公式如下：
[0046][0047]
式中，上标t表示矩阵的转置，表示k+1时刻的先验信息误差方差预测矩阵，计算公式如下，
[0048][0049]
式中，n表示集合成员的数量。
[0050]
有益效果：本发明使用两种手段将土壤含水量同化到水文模型当中，提高了水文模型的预报精度，可以为洪水预报报警提供帮助。
[0051]
手段1为直接插入法(di)，其主要思想为将实测的土壤水分数据直接插入替换模型模拟得到的土壤水分数据，而这土壤水分数据是模型模拟径流的中间状态变量，决定着模型最后的径流模拟，故这一手段将实际的土壤水分数据与模型相结合，提高了模型的预报精度。
[0052]
手段2为集合卡尔曼滤波器法(enkf)，其重要思想为将实测的土壤水分数据与模型模拟得到的土壤水分数据相结合，计算得到一个两者之间最优的状态估计值，这种方法既考虑了模型模拟的不确定性，同时也考虑了实际情况的不确定性，提高了模型的预报精度。
附图说明
[0053]
图1是misdc模型在王勿桥流域模型模拟的结果；
[0054]
图2是集合卡尔曼滤波器在misdc当中实现过程的流程图；
[0055]
图3是采用di方法使用cldas 0-5cm处的土壤水分数据同化到模型当中的结果；
[0056]
图4是采用enkf方法使用cldas 0-5cm处的土壤水分数据同化到模型当中的结果。
具体实施方式
[0057]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用于解释本发明，并不用于限定本发明。
[0058]
实施例1
[0059]
一种改善水文模型径流模拟的方法，包括步骤如下：
[0060]
s1.采用misdc模型作为水文模型，进行参数率定；参数率定的标准通过如下计算确认：
[0061][0062]
式中，q
obs
为实测流量，q
sim
为模型模拟得到的流量，上标“—”表示平均值。当计算的err值到达最小值时的参数即为最优参数；
[0063]
s2.手段1，采用直接插入法替换水文模型的土壤水分数据，将实际的土壤水分数据用于水文模型模拟径流；手段2，采用集合卡尔曼滤波器法将实际的土壤水分数据与水文模型模拟的土壤水分数据相结合，得到的状态变量值用于模型模拟径流；
[0064]
s3.由s2中两种手段得到水文模型径流模拟的优化结果，分别比较两种优化结果的ns系数和优化前后的标准化均方根误差nrmse，采用ns系数较大者和nrmse值较小者对应的手段得到径流模拟的优化结果。
[0065]
其中，所述步骤s1包括：
[0066]
s1.1通过水量平衡方程将入渗、蒸发以及排水过程联系起来求得土壤蓄水量，此过程为misdc模型中的swb模型部分；
[0067]
所述水量平衡方程如下：
[0068][0069]
式中，t为时间；w(t)为土壤蓄水量；w
max
为土壤最大蓄水量；f(t)为降雨入渗到土壤当中的过程；e(t)为蒸散率；g(t)为由于重力导致土壤排水的过程。
[0070]
所述入渗过程采用green-ampt方程，其表示如下：
[0071][0072]
式中，r(t)表示降雨过程；ks表示土壤饱和导水率；ψ表示湿润锋处的吸力；f表示
降雨事件开始时的累积入渗深度；l表示土层的厚度。
[0073]
所述蒸发过程表示如下：
[0074][0075]
式中，et
p
表示潜在蒸散发。
[0076]
所述潜在蒸散发的计算方法如下：
[0077]
et
p
＝-2+b[ξ(0.46ta(t)+8.13)]
[0078]
式中，ta(t)表示气温，单位为℃；ξ表示当天日照小时数与全年的日照小时数(365
×
12h)的比值；b为待率定的参数。
[0079]
所述排水过程的计算方法如下：
[0080]
g(t)＝ks[w(t)/w
max
]
3+2/m
[0081]
式中，ks表示土壤饱和导水率；m表示土层的孔隙大小分布指数。
[0082]
s1.2将土壤蓄水量数据输入misdc模型，输出径流模拟结果；
[0083]
其中，所述步骤s1.2包括：
[0084]
s1.2.1雨量过剩ε的计算，计算公式如下：
[0085][0086]
式中，r为降雨强度；r为降雨事件开始时的降雨深度；λ为待率定的参数；s是连接土壤蓄水量w与雨量过剩的一个参数，其计算公式如下：
[0087]
s＝a[w
max-w(t)]
[0088]
式中，a为待率定的参数；
[0089]
s1.2.2径流量计算，运用直接径流水文过程线的方法计算流域的径流量，此过程为misdc模型中的降雨-径流模型部分；其中，由于流域上的降雨从落地点到流域出口断面是需要一定时间的，而且流域上各处的降雨汇聚到出口断面的时间各不相同，但水文学家已证明流域滞时是与平均流域汇流时间等价的，因此流域的平均滞后时间г的计算公式如下所示：
[0090]
γ＝η1.19a
0.33
[0091]
式中，η是待率定的参数；a是流域的面积，单位为km2。
[0092]
其中，所述步骤s2涉及的两种手段如下：
[0093]
所述手段1的替换方法为：将实际的土壤含水量转换为实际的土壤蓄水量，并替换模型中计算得到的土壤蓄水量；其中，土壤含水量与土壤蓄水量的转换关系如下：
[0094][0095]
式中，w(t)为土壤蓄水量；θ(t)为土壤含水量；θr为土壤剩余含水量。
[0096]
所述手段2实现的具体目标函数如下：
[0097][0098]
式中，上标a和f分别表示分析和预报的过程；x
i,k+1
表示在k+1时刻，第i个成员的值；z表示观测向量；h表示观测算子；v
i,k
是均值为0，方差为rk的高斯白噪声；kk表示卡尔曼增益矩阵，其计算公式如下：
[0099][0100]
式中，上标t表示矩阵的转置，表示k+1时刻的先验信息误差方差预测矩阵，计算公式如下，
[0101][0102]
式中，n表示集合成员的数量。
[0103]
具体地，在本实施例中：
[0104]
(1)研究区域选择的是淮河流域上的王勿桥流域，流域面积约为200km2。
[0105]
(2)步骤s1中，所述的水文模型为misdc模型，使用cldas产品的2015年到2016年的降雨与气温数据作为模型的输入数据，率定参数，结果如图1所示，图中q
obs
表示实测流量，q
sim
表示模型模拟的流量。
[0106]
(3)步骤s2中，采用di法使用cldas 0-5cm处的土壤水分数据同化到水文模型当中，结果如图3所示，图中，qa表示采用数据同化方法后得到的径流量。
[0107]
(4)步骤s2中，采用enkf法同样使用cldas0-5cm处的土壤水分数据同化到水文模型当中，过程如图2所示，结果如图4所示。
[0108]
通过上述过程，得到了两个手段的结果(图3和图4)，可以看到，采用di法，模型模拟的结果的ns系数从0.326提高到了0.374，提高了约14.72％，nrmse的值为0.964，小于1，表示同化的结果是改善的。采用enkf法，模型模拟的结果的ns系数从0.326提高到了0.385，提高了约18.10％，nrmse的值为0.962，小于1，表示同化的结果是改善的。并且对比这两个结果，采用enkf法的结果改善较高。
[0109]
本实例仅仅以一个流域使用0-5cm处的土壤水分数据进行同化研究，但此方法不仅仅只限于此。其同样可以比较使用不同层的土壤水分数据对模型径流模拟结果的影响，找到一个最优的手段以及使用土壤水分数据的最优深度。
[0110]
通过采用本发明公开的上述技术方案，得到了如下有益的效果：
[0111]
本发明分别采用两种手段将实测的土壤水分数据同化到水文模型当中，改善水文模型径流模拟的精度。该方法的优势在于只需要降雨和气温数据作为模型的输入数据，同化的数据是土壤水分数据。这三种数据的获取途径有许多种，能够较为方便的运用于各个流域地区，并且避免了大量的数据处理工作。而且，对于本发明而言，采用的两种手段同化土壤水分数据对于洪水预报的精度是有所改善的。
[0112]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。