本发明属于水文技术领域，具体涉及一种水文模型在地理空间上多流域并行率定的方法。
背景技术：
资料匮乏区的水文模拟问题是水文学研究的重点与难点，受到了广泛关注。相似流域间的模型参数移植是目前资料匮乏区进行模型应用和洪水预报普遍采取的手段，但仍存在参数不确定性大、移植性不稳定的问题。流域之间参数移植方法常因为参数过度拟合问题导致移植的不对称性，加上资料匮乏地区的参数移植效果通常难以进行验证，增加了无资料地区洪水预报的误差与不确定性。如何降低资料匮乏区模型参数推衍的不确定性是提高洪水预报精度的研究重点。
模型参数主要通过观测数据模拟率定和采用物理方法进行测量的方法获得，参数率定过程可能将数据和模型的不确定性转移到参数中，而能够直接采用物理方法测量的参数难以避免在测量过程中存在不确定性。资料匮乏区参数估算是水文学研究面临的重大挑战，也是我国中小河流洪水预报面临的一大难题。很大一部分中小河流观测设施建立时间短，有些站点仅有2～3年观测资料，模型参数估算的难度与不确定性尤为突出。当流域累积资料不足以率定水文模型参数时，需要通过移植相似流域的水文信息进行洪水预报作业。国际上通用做法是根据已知流域水文、地理特征进行水文相似性分析，确定与研究区域具有相似水文响应的流域，再根据相似度移植相似流域参数，实现资料匮乏区的模型模拟和洪水预报。然而，单一流域模型参数的优选在很大程度上依赖于率定期的水文气象信息，在优选过程中容易出现模型对率定期气象条件“过度拟合”问题，参数在时间和空间交叉验证时容易出现移植效果不一致、不对称的情况。参数移植的非对称现象严重影响了相似流域间信息移用的稳定性，加上难以客观验证资料匮乏区移植效果，加大了洪水预报的误差与不确定性。
针对以上不足，如何分析模型参数的不确定性和时空移植性，获取多个流域共享的模型参数，提高模型参数的可移植性，从而降低资料匮乏区参数推衍的不确定性，正是发明人需要解决的问题。
技术实现要素：
发明目的：为了获取多个流域共同适用的水文模型参数，为无资料地区水文模拟提供稳健可靠的参数，本发明提供了一种水文模型在地理空间上多流域并行率定的方法，通过增加水量平衡控制因子的方式，调整模型结构，采用多目标优化函数进行多个流域水文模型的同时率定，获取多个流域共同适用的模型参数组。有利于提高模型参数在无资料地区模拟的精度。
技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种水文模型在地理空间上多流域并行率定的方法，该方法包括以下步骤：
步骤1，对于现有的任意一种概念性降雨-径流模型，通过引入水量平衡控制参数η进行模型结构调整，将模型参数分为新增加的水量平衡控制参数和模型原有的动态平衡参数两类分别进行率定；
步骤2，选取多个研究流域，基于结构调整后的水文模型构建单个流域的水文模型，进行单个流域水文模型的率定，获取所有研究流域的水文模型参数率定的最优值；
步骤3，构建多目标优化函数，采用折衷优化算法同时率定多个流域水文模型参数，获取多个流域同时适用的模型参数。
进一步的，所述步骤1，具体包括：
步骤1.1，选取任意一种现有的概念性水文模型，对于选定的水文模型，引入水量平衡控制参数η进行模型结构调整，当土壤含水量与田间持水量的比值大于时η时，蒸散发量为潜在蒸散发能力，当土壤含水量与田间持水量的比值小于时η时，流域蒸散发受到土壤含水量的限制；
式中，η为水量平衡控制因子，eta为实际蒸散发，etp为蒸散发能力，sm为土壤含水量，fc为田间持水量；
步骤1.2，将所构建的水文模型参数分为控制水量平衡的参数η和控制动态状态的参数θ，θ即为选定的模型原有参数的集合，通过蒙特卡洛取样方法，根据选定模型参数的取值范围，生成n组模型动态参数组θ(i)，1≤i≤n，对每组参数都采用本文所选取的模型进行降雨-径流模拟，从而可以得出模拟的总径流为：
f(η)＝vim(η,θ)(2)
式中，vim为利用第i组动态状态控制参数及其对应的水量平衡参数进行模拟得到的总径流量，水量平衡参数η的取值在0-1之间，f(η)随着η值的增大而增大，若是模型能够达到水量平衡状态，则意味着在模型原有参数θ一定的情况下，必然有η在[0,1]之间搜索符合的值来满足以下公式：
vim(1,θ)＜vio＜vim(0,θ)(3)
式中，vio为观测的总径流，vim(1,θ)为在参数θ一定的情况下η取值为1时模拟得到的总径流，vim(0,θ)为在参数θ一定的情况下η取值为0时模拟得到的总径流；
针对特定的参数θ可以通过二分法搜索在[0,1]之间搜索能够满足模拟总径流等于观测总径流的η值：
vim(η(θ),θ)＝vio(4)
η(θ)为针对特定参数θ能够满足整个率定期内总水量平衡的η取值。
进一步的，所述步骤2中，选取多个研究流域，基于结构调整后的水文模型构建单个流域的水文模型，进行单个流域水文模型的率定，获取所有研究流域的水文模型参数率定的最优值，具体包括以下步骤：
步骤2.1，选择多个具有历史观测数据的研究流域，总流域个数为m，针对选取的研究流域，分别应用水文模型，单独进行每一个流域的模型参数率定；
步骤2.2，根据所选取的水文模型的原有参数(即本研究中的控制动态状态的参数θ)的可能取值范围，通过蒙特卡罗随机采样方法在取值范围内生成预定组数(设为n组)的随机动态参数组；
步骤2.3，利用公式(4),对第i组动态参数组θi，利用二分法在对η在[0,1]内进行搜索取值，通过模型模拟获取能够达到水量平衡的ηi值；
步骤2.4，基于参数ηi和θi，计算第i组参数的模拟目标函数值，选择水文模拟常用的纳什确定性作为本研究目标函数，计算公式如下；
式中，oi(ηi,θi)为由第i组参数进行模拟得到的纳什确定性，t为模拟的总时段数，qo(t)为t时刻的实测流量，qmi(t)为采用第i组参数在t时刻得到的模拟流量，
步骤2.5，对所有参数组进行模拟后，选取所有参数组中模拟效果最好的，即纳什确定性系数最大值作为该研究流域(设定流域编号为j)的最优模拟效果值
式中，
进一步的，所述步骤3中，构建多目标优化函数，采用折衷优化算法同时率定多个流域水文模型参数的，获取多个流域同时适用的模型参数，方法如下：
步骤3.1，基于折衷优化算法，构建多流域同时获得高精度模拟结果的多目标优化函数：
式中，oa(θi)为采用第i组模型动态参数组θi进行所有流域模拟得到的纳什确定性系数值；m为进行多流域同时率定的总流域个数；
步骤3.2，利用公式(7)，对所有通过蒙特卡罗方法随机生成的模型动态参数组θi，计算相应的水量平衡参数ηi和多个流域同时率定的最优目标函数值oa(θi)；
步骤3.3，对所有模拟的多目标函数值{oa(θ1),oa(θ2),...,oa(θn)}通过排序比较，选取最大值作为最佳目标函数值o′a，该目标函数其所对应的模型参数(θ′)即为率定得到的多个流域共同适用的模型参数；
步骤3.4，利用多流域共同适用参数(θ′)进行无实测资料流域的模拟。
进一步的，所述流域水文模型模拟的时段长为24小时。
有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：
本发明提供的一种水文模型在地理空间上多流域并行率定的方法，通过增加水量平衡控制因子，调整水文模型的结构，将模型参数分为水量平衡控制参数和动态过程控制参数，基于长期水量平衡关系实现对任一动态参数组所对应的水量平衡控制因子的率定。依据折衷优化算法，构建多流域同时模拟的多目标优化函数，实现对地理空间上多个流域水文模型的同时率定，获取稳健可靠的高移植性参数，解决资料匮乏区水文模型应用的难题。本方法充分考虑了流域之间参数移植方法常因为参数的过度拟合问题导致移植的不对称性，实现多流域共享参数的快速自动率定方法，既保证了不同流域之间共享参数移植的可靠性和稳定性，又解决了无资料地区水文模型参数推衍的难题，且方法中充分利用了多个流域的水文气象信息，增加了历史数据的丰富性，考虑了不同流域水文模型模拟精度的差异性，保证了结果的客观合理性，能够进一步促进无资料地区水文模型应用研究的深入发展。
附图说明
图1是本发明的计算流程示意图。
图2为本发明选取的有资料流域示意图。
图3为本发明中进行单个流域水文模型率定结果示意图。
图4为本发明中计算出的单个流域参数移植结果示意图。
图5为本发明中计算出的多个流域并行率定结果示意图。
图6为本发明中选取的假定无资料流域示意图。
图7为本发明计算的无资料流域参数移植结果示意图。
具体实施方式
下面结合附图和具体实施例对本发明作进一步描述。
如图1所示，本发明提供的一种水文模型在地理空间上多流域并行率定的方法，该方法包括以下步骤：
步骤1，对于现有的任意一种概念性降雨-径流模型，通过引入水量平衡控制参数η进行模型结构调整，将模型参数分为新增加的水量平衡控制参数和模型原有的动态平衡参数两类分别进行率定；
步骤2，选取多个研究流域，基于结构调整后的水文模型构建单个流域的水文模型，进行单个流域水文模型的率定，获取所有研究流域的水文模型参数率定的最优值；
步骤3，构建多目标优化函数，采用折衷优化算法同时率定多个流域水文模型参数，获取多个流域同时适用的模型参数。
进一步的，所述步骤1，具体包括：
步骤1.1，选取任意一种现有的概念性水文模型，对于选定的水文模型，引入水量平衡控制参数η进行模型结构调整，当土壤含水量与田间持水量的比值大于时η时，蒸散发量为潜在蒸散发能力，当土壤含水量与田间持水量的比值小于时η时，流域蒸散发受到土壤含水量的限制；
式中，η为水量平衡控制因子，eta为实际蒸散发，etp为蒸散发能力，sm为土壤含水量，fc为田间持水量；
步骤1.2，将所构建的水文模型参数分为控制水量平衡的参数η和控制动态状态的参数θ，θ即为选定的模型原有参数的集合，通过蒙特卡洛取样方法，根据选定模型参数的取值范围，生成n组模型动态参数组θ(i)，1≤i≤n，对每组参数都采用本文所选取的模型进行降雨-径流模拟，从而可以得出模拟的总径流为：
f(η)＝vim(η,θ)(2)
式中，vim为利用第i组动态状态控制参数及其对应的水量平衡参数进行模拟得到的总径流量，水量平衡参数η的取值在0-1之间，f(η)随着η值的增大而增大，若是模型能够达到水量平衡状态，则意味着在模型原有参数θ一定的情况下，必然有η在[0,1]之间搜索符合的值来满足以下公式：
vim(1,θ)<vio<vim(0,θ)(3)
式中，vio为观测的总径流，vim(1,θ)为在参数θ一定的情况下η取值为1时模拟得到的总径流，vim(0,θ)为在参数θ一定的情况下η取值为0时模拟得到的总径流；
针对特定的参数θ可以通过二分法搜索在[0,1]之间搜索能够满足模拟总径流等于观测总径流的η值：
vim(η(θ),θ)＝vio(4)
η(θ)为针对特定参数θ能够满足整个率定期内总水量平衡的η取值。
进一步的，所述步骤2中，选取多个研究流域，基于结构调整后的水文模型构建单个流域的水文模型，进行单个流域水文模型的率定，获取所有研究流域的水文模型参数率定的最优值，具体包括以下步骤：
步骤2.1，选择多个具有历史观测数据的研究流域，总流域个数为m，针对选取的研究流域，分别应用水文模型，单独进行每一个流域的模型参数率定；
步骤2.2，根据所选取的水文模型的原有参数(即本研究中的控制动态状态的参数θ)的可能取值范围，通过蒙特卡罗随机采样方法在取值范围内生成预定组数(设为n组)的随机动态参数组；
步骤2.3，利用公式(4),对第i组动态参数组θi，利用二分法在对η在[0,1]内进行搜索取值，通过模型模拟获取能够达到水量平衡的ηi值；
步骤2.4，基于参数ηi和θi，计算第i组参数的模拟目标函数值，选择水文模拟常用的纳什确定性作为本研究目标函数，计算公式如下；
式中，oi(ηi,θi)为由第i组参数进行模拟得到的纳什确定性，t为模拟的总时段数，qo(t)为t时刻的实测流量，qmi(t)为采用第i组参数在t时刻得到的模拟流量，
步骤2.5，对所有参数组进行模拟后，选取所有参数组中模拟效果最好的，即纳什确定性系数最大值作为该研究流域(设定流域编号为j)的最优模拟效果值
式中，
进一步的，所述步骤3中，构建多目标优化函数，采用折衷优化算法同时率定多个流域水文模型参数的，获取多个流域同时适用的模型参数，方法如下：
步骤3.1，基于折衷优化算法，构建多流域同时获得高精度模拟结果的多目标优化函数：
式中，oa(θi)为采用第i组模型动态参数组θi进行所有流域模拟得到的纳什确定性系数值；m为进行多流域同时率定的总流域个数；
步骤3.2，利用公式(7)，对所有通过蒙特卡罗方法随机生成的模型动态参数组θi，计算相应的水量平衡参数ηi和多个流域同时率定的最优目标函数值oa(θi)；
步骤3.3，对所有模拟的多目标函数值{oa(θ1),oa(θ2),...,oa(θn)}通过排序比较，选取最大值作为最佳目标函数值o′a，该目标函数其所对应的模型参数(θ′)即为率定得到的多个流域共同适用的模型参数；
步骤3.4，利用多流域共同适用参数(θ′)进行无实测资料流域的模拟。
进一步的，所述流域水文模型模拟的时段长为24小时。
以简单的降雨-径流模型hymod和美国20个研究流域为例，hymod是简单的概念性降雨-径流模型，包括融雪计算、蒸发、产流与汇流等模块，共有7个参数需要通过实测资料进行率定。研究流域的1970-1979年共10年的逐日降水、蒸发、流量和气温数据采用美国地质调查局(usgs)提供的历史观测数据。具体为：
步骤一、在hymod模型原有结构的基础上，进行水文模型结构的调整，引入水量平衡控制参数η，将模型参数分为水量平衡控制参数和动态平衡控制参数(即模型原有的所有参数集合)两类分别进行率定。具体为：
1)在hymod模型中引入水量平衡控制参数η，当土壤含水量与相当于田间持水量(在hymod中体现为最大蓄水能力cmax的比值大于时η时，蒸散发量为潜在蒸散发能力，当土壤含水量与田间持水量的比值小于时η时，流域蒸散发受到土壤含水量的限制；
式中：η为水量平衡控制因子，eta为实际蒸散发，etp为蒸散发能力，sm为土壤含水量，cmax为最大蓄水能力。
2)将hymod模型参数分为控制水量平衡的参数η和控制动态状态的参数θ，θ即为hymod模型原有参数的集合。通过蒙特卡洛取样方法，根据选定模型参数的取值范围，生成预设组数(n)的模型动态参数组θ(n)，从而可以得出利用本文构建模型进行流域降雨-径流模拟的总径流为：
f(η)＝vim(η,θ)(2)
式中，vim为利用第i套参数进行模拟得到的总径流量，水量平衡参数η的取值在0-1之间，f(η)随着η值的增大而增大，若是模型能够达到水量平衡状态，则意味着在模型原有参数θ一定的情况下，必然有η在[0,1]之间搜索符合的值来满足以下公式：
vim(1,θ)＜vio＜vim(0,θ)(3)
式中，vio为观测的总径流；vim(1,θ)为在参数θ一定的情况下η取值为1时模拟得到的总径流；vim(0,θ)为在参数θ一定的情况下η取值为0时模拟得到的总径流。
因此针对特定的参数θ可以通过二分法搜索在[0,1]之间搜索能够满足模拟总径流等于观测总径流的η值：
vim(η(θ),θ)＝vio。(4)
η(θ)为针对特定参数θ能够满足长期水量平衡的η取值。
步骤二、对于20个流域，选取10个作为有资料流域(如图2所示)，构建流域水文模型，进行单个流域水文模型的率定，获取所有研究流域的水文模型参数率定的最优值，具体为：
1)针对选取的研究流域，分别构建hymod模型，单独进行每一个流域的模型参数率定；
2)根据hymod模型参数的可能取值范围，通过蒙特卡罗随机采样生成10000组数的随机动态参数组；
3)对第i组动态参数组θi，利用二分法在对η在[0,1]内进行搜索取值，通过模型模拟获取能够达到水量平衡的ηi值；
4)基于参数ηi和θi，计算第i组参数的模拟目标函数值，本文采用纳什确定性系数作为目标函数进行参数优化，计算公式如下；
式中，oi(ηi,θi)为由第i组参数进行模拟得到的纳什确定性；t为模拟的总时段数；qo(t)为t时刻的观测流量；qmi(t)为采用第i组参数在t时刻得到的模拟流量；
5)对所有参数组进行模拟后，选取所有参数组中模拟效果最好的，即纳什确定性系数最大值作为该研究流域(设定流域编号为j)的最优模拟效果值
式中，
模拟结果如图3所示。同时，将每个流域率定的参数组移用到另外9个流域，进行流域之间参数移植效果的比较(图4)。
步骤三、构建多目标优化函数，采用折衷优化算法进行多个流域水文模型参数的同时率定，获取多个流域同时适用的模型参数，具体为：
1)基于折衷优化算法，构建所选定的10个流域同时获得高精度模拟结果的多目标优化函数：
式中：oa(θi)为采用第i套模型动态参数组θi进行所有流域模拟得到的纳什确定性系数值；m为进行多流域同时率定的总流域个数；
2)利用公式(7)，对所有10000组随机参数，计算相应的水量平衡参数ηi和多个流域同时率定的最优目标函数值oa(θi)；
3)对所有模拟的多目标函数值{oa(θ1),oa(θ2),...,oa(θn)}通过排序比较，选取最大值作为最佳目标函数值o＇a，该目标函数其所对应的参数组(θ′)即为率定得到的多个流域共同适用的模型参数。图5显示了10个流域并行率定的结果及其与单个流域率定、单个流域参数移用的模拟结果比较。从中可以看出，与单个流域参数移用相比，多流域共享参数的模拟效果更为稳定。
4)利用多流域共同适用参数(θ′)进行无实测资料流域的模拟。选取另外10个流域作为无资料流域(如图6所示)，进行多流域并行率定参数在无资料地区移植效果的检验。图7显示了10个用于参数率定流域的参数移用到无资料地区的模拟结果，分别是单个流域率定参数的移植和多流域并行率定参数的移植检验。图7同样显示了假定无资料流域的单个流域率定结果作为比对。结果显示了多流域并行率定参数在无资料地区移用的优势。
另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。
以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。
1.一种水文模型在地理空间上多流域并行率定的方法，其特征在于，该方法包括以下步骤：
步骤1，对于现有的任意一种概念性降雨-径流模型，通过引入水量平衡控制参数η进行模型结构调整，将模型参数分为新增加的水量平衡控制参数和模型原有的动态平衡参数两类分别进行率定；
步骤2，选取多个研究流域，基于结构调整后的水文模型构建单个流域的水文模型，进行单个流域水文模型的率定，获取所有研究流域的水文模型参数率定的最优值；
步骤3，构建多目标优化函数，采用折衷优化算法同时率定多个流域水文模型参数，获取多个流域同时适用的模型参数。
2.根据权利要求1所述的一种水文模型在地理空间上多流域并行率定的方法，其特征在于，所述步骤1，具体包括：
步骤1.1，选取任意一种现有的概念性水文模型，对于选定的水文模型，引入水量平衡控制参数η进行模型结构调整，当土壤含水量与田间持水量的比值大于时η时，蒸散发量为潜在蒸散发能力，当土壤含水量与田间持水量的比值小于时η时，流域蒸散发受到土壤含水量的限制；
式中，η为水量平衡控制因子，eta为实际蒸散发，etp为蒸散发能力，sm为土壤含水量，fc为田间持水量；
步骤1.2，将所构建的水文模型参数分为控制水量平衡的参数η和控制动态状态的参数θ，θ即为选定的模型原有参数的集合，通过蒙特卡洛取样方法，根据选定模型参数的取值范围，生成n组模型动态参数组θ(i)，1≤i≤n，对每组参数都采用本文所选取的模型进行降雨-径流模拟，从而可以得出模拟的总径流为：
f(η)＝vim(η，θ)(2)
式中，vim为利用第i组动态状态控制参数及其对应的水量平衡参数进行模拟得到的总径流量，水量平衡参数η的取值在0-1之间，f(η)随着η值的增大而增大，若是模型能够达到水量平衡状态，则意味着在模型原有参数θ一定的情况下，必然有η在[0，1]之间搜索符合的值来满足以下公式：
vim(1，θ)＜vio＜vim(0，θ)(3)
式中，vio为观测的总径流，vim(1，θ)为在参数θ一定的情况下η取值为1时模拟得到的总径流，vim(0，θ)为在参数θ一定的情况下η取值为0时模拟得到的总径流；
针对特定的参数θ可以通过二分法搜索在[0，1]之间搜索能够满足模拟总径流等于观测总径流的η值：
vim(η(θ)，θ)＝vio(4)
η(θ)为针对特定参数θ能够满足整个率定期内总水量平衡的η取值。
3.根据权利要求2所述的一种水文模型在地理空间上多流域并行率定的方法，其特征在于，所述步骤2中，选取多个研究流域，基于结构调整后的水文模型构建单个流域的水文模型，进行单个流域水文模型的率定，获取所有研究流域的水文模型参数率定的最优值，具体包括以下步骤：
步骤2.1，选择多个具有历史观测数据的研究流域，总流域个数为m，针对选取的研究流域，分别应用水文模型，单独进行每一个流域的模型参数率定；
步骤2.2，根据所选取的水文模型的原有参数，即控制动态状态的参数θ的可能取值范围，通过蒙特卡罗随机采样方法在取值范围内生成n组的随机动态参数组；
步骤2.3，利用公式(4)，对第i组动态参数组θi，利用二分法在对η在[0，1]内进行搜索取值，通过模型模拟获取能够达到水量平衡的ηi值；
步骤2.4，基于参数ηi和θi，计算第i组参数的模拟目标函数值，选择水文模拟常用的纳什确定性作为本研究目标函数，计算公式如下；
式中，oi(ηi，θi)为由第i组参数进行模拟得到的纳什确定性，t为模拟的总时段数，qo(t)为t时刻的实测流量，qmi(t)为采用第i组参数在t时刻得到的模拟流量，
步骤2.5，对所有参数组进行模拟后，选取所有参数组中模拟效果最好的，即纳什确定性系数最大值作为该研究流域的最优模拟效果值
式中，
4.根据权利要求1所述的一种水文模型在地理空间上多流域并行率定的方法，其特征在于，所述步骤3中，构建多目标优化函数，采用折衷优化算法同时率定多个流域水文模型参数的，获取多个流域同时适用的模型参数，方法如下：
步骤3.1，基于折衷优化算法，构建多流域同时获得高精度模拟结果的多目标优化函数：
式中，oa(θi)为采用第i组模型动态参数组θi进行所有流域模拟得到的纳什确定性系数值；m为进行多流域同时率定的总流域个数；
步骤3.2，利用公式(7)，对所有通过蒙特卡罗方法随机生成的模型动态参数组θi，计算相应的水量平衡参数ηi和多个流域同时率定的最优目标函数值oa(θi)；
步骤3.3，对所有模拟的多目标函数值{oa(θ1)，oa(θ2)，...，oa(θn)}通过排序比较，选取最大值作为最佳目标函数值o′a，该目标函数其所对应的模型参数θ′即为率定得到的多个流域共同适用的模型参数；
步骤3.4，利用多流域共同适用参数θ′进行无实测资料流域的模拟。
5.根据权利要求1-4任一所述的一种水文模型在地理空间上多流域并行率定的方法，其特征在于，所述流域水文模型模拟的时段长为24小时。
技术总结