CN108897977A - 一种基于大区域水文模拟的径流演变不确定归因方法 - Google Patents

Abstract

本发明公开了一种基于大区域水文模拟的径流演变不确定归因方法，该方法考虑了模型参数不确定性演化和归因路径多样性而导致不确定性的时空径流演变归因，包括如下步骤：建立数据集；描述历史径流演变规律；建立大区域水文模型；模型参数敏感性分析；模型参数率定及不确定性分析；径流演变不确定归因分析。本发明能够利用有限的径流监测资料获得全区域径流变化归因结果，并且充分考虑了模型参数不确定性演进和归因路径多样而导致不确定性。

Description

一种基于大区域水文模拟的径流演变不确定归因方法

技术领域

本发明设计水利工程领域中的径流演变归因技术，特别涉及一种基于大区域水文模拟的、考虑了模型参数不确定性演化和归因路径多样性而导致不确定性的时空径流演变归因方法。

背景技术

径流是水文循环中最重要的组成部分之一，在当前的变化环境下，理解径流的产生、变化及变化的潜在原因，对进行高效的水资源管理有着重要意义。径流过程与大气环流、气候变化、流域内下垫面和人类社会经济等诸多要素密切相关，径流变化是这些要素共同作用、交织影响的综合结果，从而显示出复杂多变、难以预测的特征。伴随人口的快速增长，水资源的供需矛盾日益激烈，厘清气候变化和人类活动等诸多因子分别对径流的影响作用，对预测未来水资源情势、水资源管理适应性决策有着至关重要的作用。

径流演变的归因技术旨在定量分析径流时空变化的成因，为预测未来径流、制定水资源管理的适应性对策提供依据，针对不同的成因和影响程度做出具有针对性、便于操作的适应性调控决策。国内外变化环境下径流变化的归因研究取得了长足的发展，但从当前见诸文献的成果看，径流变化归因技术有待进一步深化和改进，表现为：对于大区域的径流演变归因，多采用数据分析或概念性方法，笼统给出整个区域的归因结果；即使能够分子流域进行径流演变归因，对水文监测数据有着极高的要求，往往也受限于此；归因分析的不确定性来源很多，而不确定性评估尚显不足。

发明内容

发明目的：提供一种基于大区域水文模拟的、考虑了模型参数不确定性演化和归因路径多样性而导致不确定性的时空径流演变归因方法，以解决现有技术存在的上述问题。

技术方案：一种基于大区域水文模拟的径流演变不确定归因方法，包括以下步骤：

步骤1，建立数据集；

步骤2，描述历史径流演变规律；

步骤3，建立大区域水文模型；

步骤4，模型参数敏感性分析；

步骤5，模型参数率定及不确定性分析；

步骤6，径流演变不确定归因分析。

所述步骤1中，建立包含水文气象数据、空间数据和地理特征数据的数据集。更详细地，包括：水文站径流监测数据、气象站降水、最高气温、最低气温、风速、相对湿度、太阳辐射数据，地形地貌数据，土地利用与覆被数据，土壤类型数据，水库数据，取用水数据。

所述步骤2中，不同的统计方法对同一个时间序列检验可能得到不一致的结果，因此采用Mann-Kendall非参数检验和Spearman’s Rho检验两种方法描述历史径流的演变特征并相互佐证；采用Pettitt检验和有序聚类方法确定径流发生统计意义上显著突变的时间，划分基准期和变化期，认定该时间点前的基准期径流受变化环境影响统计意义上不显著，可反映天然情况。

所述步骤3中，依据大区域地形地貌数据生成水系、划分子流域，基于大区域地形地貌数据、土地利用覆被数据、土壤类型数据生成水文响应单元，从而建立起大区域物理水文模型SWAT模型结构，输入数据试运行。

所述步骤4中，采用LH-OAT(Latin Hypercube One-factor-at-a-time)敏感性分析技术，在模型众多参数中选出重要敏感的参数。为了反映大区域模型参数的空间异质性，对每个观测控制点以上的子流域集合采用一套敏感性参数集，整个大区域采用多套敏感性参数集。

所述步骤5中，针对各流域(每个观测控制点以上的子流域集合)的敏感参数集进行参数率定和不确定分析。大区域的水文模拟包含了很多不确定性和复杂性，这些不确定性和复杂性可以反映在模型参数上。本发明采用贝叶斯理论进行参数不确定性分析，采用贝叶斯形式，模型参数的后验分布可以表述为：

式中，θ为模型参数，为观测值，为模型参数的后验分布，p(θ)为模型参数的先验分布，为似然函数，为常数。进而模型参数的后验分布可以表述为：

假定残差是不相关的，那么似然函数可以写成：

式中，f_a(b)指基于b对a概率密度函数的估计。进一步假定残差服从正态分布那么似然函数可以写成：

本发明认定模型参数的先验分布在其初始范围中服从均匀分布，一旦先验分布和似然函数定义完成，就可以推断模型参数的后验分布。然而大区域的水文模型是高度非线性系统，无法使用解析的方法获得后验分布。此时，蒙特卡罗方法(MC)方法可以被用来生成后验分布的样本，然而标准的MC方法在计算上效率不高，因此本发明采用马尔可夫链蒙特卡罗方法(MCMC)来获得后验分布的样本，通过应用DREAM(DiffeRential EvolutionAdaptive Metropolis)算法实现。

所述步骤6进一步分为：

步骤61、筛选径流演变影响因子。大区域水文模型的输入包括，气象输入：降水、降水、最高气温、最低气温、风速、相对湿度、太阳辐射数据，人类活动输入：土地利用与覆被数据，水库数据，取用水数据。分析突变时间前后，输入因子的统计特征值，筛选出发生统计意义显著变化的因子作为归因分析的影响因子。

步骤62、根据筛选出的影响因子生成情景集。

步骤63、将情景集输入模型，各子流域敏感参数集分别采用最佳模拟、不确定性上限，不确定性下限，分别模拟对应径流。

步骤64、计算各因子对径流变化的贡献及参数不确定性演化的归因不确定性区间。

为了更清楚地阐述步骤62-64，给出以下2因子、4因子、N因子算例。当仅含2项影响因子F₁、F₂时，每项因子具有2个状态，分为“-1”基准期状态，“1”变化期状态，情景组合方案可以设置为情景S₁至S₄：

为了书写简便记因子F₁对径流的影响量可以表示为：

或

因子F₂对径流的影响量可以表示为：

或

这就是归因路径的多样性问题，不同的路径计算得到的因子贡献量是不相等的，因为真实情况下自然界的影响因子间是相互不独立的，在不同状态的F₂下，F₁对径流的影响量是不一样的。为了规避归因路径的多样性带来的不确定性，本发明定义影响因子F₁、F₂对径流变化的贡献量分别为：

将影响因子对径流变化的贡献量定义为各路径的平均值，这样求得某因子的贡献量为其他因子平均状态下的贡献量。因此4因子情形下，情景组合方案可以设置为：

影响因子F₁、F₂、F₃、F₄对径流变化的贡献量按前述定义分别为：

进一步推广到N因子情形，第q个因子对径流变化的贡献量可以表述为：

式中，α_q,p为权重系数，表示因子F_q在情景p下的状态，等于-1或1，具体取值可以在情景组合方案表中查得。

有益效果：首先建立起的大区域物理水文模型，利用有限的水文观测资料率定合格后，能够从机理上准确地描述细致到水文响应单元的水文进程；筛选出重要的影响因子后，构造多因子情景组合方案集，通过大区域物理水文模型模拟，能够获得全空间的归因分析结果；考虑了参数不确定性演化和归因路径多样性引起的不确定性，给出了归因结果的不确定性区间。

附图说明

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

具体实施方式

下面结合附图，通过实施例对本发明的技术方案做进一步具体描述。

在进一步的实施例中，本发明建立包含水文气象数据、空间数据和地理特征数据的数据集。更详细地，包括：水文站径流监测数据、气象站降水、最高气温、最低气温、风速、相对湿度、太阳辐射数据，地形地貌数据，土地利用与覆被数据，土壤类型数据，水库数据，取用水数据。采用Mann-Kendall非参数检验和Spearman’s Rho检验两种方法描述历史径流的演变特征并相互佐证；采用Pettitt检验和有序聚类方法确定径流发生统计意义上显著突变的时间，划分基准期和变化期，认定该时间点前的基准期径流受变化环境影响统计意义上不显著，可反映天然情况。依据大区域地形地貌数据生成水系、划分子流域，基于大区域地形地貌数据、土地利用覆被数据、土壤类型数据生成水文响应单元，从而建立起大区域物理水文模型SWAT模型结构，输入数据试运行。采用LH-OAT敏感性分析技术，在模型众多参数中选出重要敏感的参数。为了反映大区域模型参数的空间异质性，对每个观测控制点以上的子流域集合采用一套敏感性参数集，整个大区域采用多套敏感性参数集。针对各流域(每个观测控制点以上的子流域集合)的敏感参数集进行参数率定和不确定分析。筛选径流演变影响因子，根据筛选出的影响因子生成情景集，将情景集输入模型，各子流域敏感参数集分别采用最佳模拟、不确定性上限，不确定性下限，分别模拟对应径流，计算各因子对径流变化的贡献及参数不确定性演化的归因不确定性区间。

具体的，如图1所示，一种基于大区域水文模拟的径流演变不确定归因方法，包括如下步骤：

步骤1，建立包含水文气象数据、空间数据和地理特征数据的数据集。更详细地，包括：水文站径流监测数据、气象站降水、最高气温、最低气温、风速、相对湿度、太阳辐射数据，地形地貌数据，土地利用与覆被数据，土壤类型数据，水库数据，取用水数据。

步骤2，采用Mann-Kendall非参数检验和Spearman’s Rho检验两种方法描述历史径流的演变特征并相互佐证。对于时间序列X，样本容量为n，构造Mann-Kendall检验统计量S和Z_c：

式中，x_d和x_c为时间序列X的要素，当Z_c＞0时，表明时间序列X随时间有增加趋势，Z_c＜0时表明时间序列X随时间有减少趋势，当|Z_c|＞Z_1-α/2时，Z_1-α/2为标准正态离差，α为显著性水平，表明时间序列变化趋势统计意义显著。

构造Spearman’s Rho检验统计量ρ_s：

式中，x_i为样本，t_i为时间，为样本均值，为时间均值。当样本容量n>30时，统计量ρ_s为正态分布，记当Z_c＞0时，表明时间序列X随时间有增加趋势，Z_c＜0时表明时间序列X随时间有减少趋势，当|Z_c|＞Z_1-α/2时，变化趋势统计意义显著。

采用Pettitt检验和有序聚类方法确定径流发生统计意义上显著突变的时间。构造时间t时的Pettitt检验统计量U_t,n：

记统计量U_t,N的最大值为k_τ＝max{U_t,n}，对应的时间τ就是突变点，显著性水平检验公式如下：

P＝2exp{-6(k_τ)²/(n³+n²)} (21)

有序聚类发则将原径流时间序列X分成两组样本{xh}(k＝1,2,…,τ)和{xh}(k＝τ+1,2,…,n)，通过样本的统计特性识别序列的突变点，实质是在保持原序列次序的基础上推求最优分割点。最优分割点即使得分割后两样本内部的离差平方和较小，两样本间的离差平方和较大的序列点。记分割后前后两组样本内部的离差平方和为：

式中，V_τ和V_n-τ为离差平方和，和为两组样本均值，总离差平方和为：

S_n(τ)＝V_τ+V_n-τ (24)

满足条件S_n＝min{S_n(τ)}的τ为最优分割点；

按照突变时间前后将研究时期划分为基准期和变化期，认定该时间点前的基准期径流受变化环境影响统计意义上不显著，可反映天然情况。

步骤3，建立大区域水文模型，具体分为以下几个子步骤：

步骤31、生成水系，划分子流域：依据地形地貌数据生成大区域河网水系结构，包括河网水系、连接点、出水口、入水口。根据实际情况修整河网水系、连接点，添加另外的出入水口或删除不需要的出入水口。设定子流域面积阈值，根据水力联系划分子流域，计算子流域几何参数、地形参数和水流路径，生成子流域报告。

步骤32、生成水文响应单元：载入坡度数据、土地利用覆被数据、土壤类型数据，设定归类阈值，将阈值内的坡度、土地利用覆被、土壤类型认定为同一类，同时具有相同类坡度、土地利用覆被、土壤类型的单元为水文响应单元，是模型结构的最小单元，也是水文平衡计算的参照单元。生成水文响应单元报告。

步骤33、输入数据试运行：创建模型数据库，包括配置文件、土壤数据、天气发生数据、子流域数据、水文响应单元数据、主要河道数据、地下水数据、取用水数据、管理数据、水库数据、流域数据，试运行模型。

步骤4，模型参数敏感性分析：采用LH-OAT(Latin Hypercube One-factor-at-a-time)敏感性分析技术，在模型众多参数中选出重要敏感的参数。首先将P个参数的原始范围分成N段，对其进行N次拉丁超立方抽样，通过扰动抽样点P次，计算敏感度。对于第a个参数在第b

个抽样点处的敏感度S_a,b(％)可以由下式计算：

式中，M(·)代表模型方程，f_a代表参数扰动，e_a,b代表第a个参数在第b个抽样点处的取值。参数的敏感性可以依据每个参数的敏感度均值排序，选取最敏感的。为了反映大区域模型参数的空间异质性，对每个观测控制点以上的子流域采用一套敏感性参数集，整个大区域采用多套敏感性参数集。

步骤5，针对各流域(每个观测控制点以上的子流域集合)的敏感参数集进行参数率定和不确定分析。大区域的水文模拟包含了很多不确定性和复杂性，这些不确定性和复杂性可以反映在模型参数上。本发明采用贝叶斯理论进行参数不确定性分析，采用马尔可夫链蒙特卡罗方法(MCMC)来获得敏感参数集后验分布的样本，通过应用DREAM(DiffeRential Evolution Adaptive Metropolis)算法实现。

本发明采用贝叶斯理论进行参数不确定性分析，采用贝叶斯形式，模型参数的后验分布可以表述为：

式中，θ为模型参数，为径流量观测值，为模型参数的后验分布，p(θ)为模型参数的先验分布，为似然函数，为常数，进而模型参数的后验分布可以表述为：

假定残差是不相关的，那么似然函数可以写成：

式中，y₁(θ)、y₂(θ)、y_n(θ)为模拟径流，为中元素，f_a(b)指基于b对a概率密度函数的估计，进一步假定残差服从正态分布那么似然函数可以写成：

模型参数的先验分布在其初始范围中服从均匀分布，一旦先验分布和似然函数定义完成，就可以推断模型参数的后验分布；采用马尔可夫链蒙特卡罗方法(MCMC)来获得后验分布的样本，通过应用DREAM(DiffeRential Evolution Adaptive Metropolis)算法实现；具体分为以下几个子步骤：

步骤51、认为参数在原始范围内服从均匀分布分布，生成候选样本：

式中θ为模型参数，和为迭代过程中生成的候选样本，k为迭代过程的某一次，γ(δ)、Z₁(a)、Z₂(b)为DREAM算法定义的参数，残差ε～N_d(0,b)，δ、a、b为DREAM算法定义的参数，b为一个小的常数；

步骤52、根据交叉概率CR替换候选样本：

式中，U∈[0,1]服从均匀分布，CR为交叉概率，d_eff为DREAM算法定义的参数；

步骤53、计算候选样本的Metropolis接受概率：

式中Y_obs为径流观测值。

步骤54、如果接受则向前移动，不接受仍保持原位；

重复步骤51-54，直到达到所需规模。

步骤55、采用Gelman收敛诊断指标判断抽样序列是否收敛，当时，收敛，并认定收敛后的样本分布稳定。

步骤56、选取收敛后的敏感参数样本作为敏感参数集的后验分布抽样，采用构建好的模型模拟所以敏感参数样本取值下的径流量。选取满足E_NS＝max{E_NS(i)}，即E_NS最大的一次作为模型的最佳模拟，E_NS如下式计算：

式中，Q_oi为第i时段的实测径流，为各时段实测径流的平均值，Q_si为第i时段的模拟模拟。

E_NS用以描述模型的拟合程度。除E_NS外，采用确定性系数R²评价模拟径流和实测径流的线性相关程度，采用相对误差Re(％)评价模拟结果的偏差，计算公式分别为：

式中，为各时段模拟径流的平均值，其他含义同上。

引入95％不确定性区间描述不确定性，定义为敏感参数集按后验分布抽样升序或降序排列，由2.5％和97.5％分位点所夹区间的参数样本模拟得到的径流区间。采用P-factor和R-factor两个指标来模拟的不确定性。P-factor定义为实测径流被包含在95％不确定性区间中的百分比；R-factor定义为95％不确定性区间平均宽度比标准差，计算公式为：

式中Q_U和Q_L分别为95％不确定性区间的上界和下界，σ为实测径流的标准差。P-factor越接近1，R-factor越接近0，结果越好，当R-factor小于1时，可以认为模拟结果满足要求。

步骤6，径流演变不确定归因分析，具体分为以下几个子步骤：

步骤61、筛选径流演变影响因子。大区域水文模型的输入包括，气象输入：降水、降水、最高气温、最低气温、风速、相对湿度、太阳辐射数据，人类活动输入：土地利用与覆被数据，水库数据，取用水数据。分析突变时间前后，输入因子的统计特征值，筛选出N个发生统计意义显著变化的因子作为归因分析的影响因子；

步骤62、根据筛选出的影响因子生成N情景集；

步骤63、将情景集输入模型，各子流域敏感参数集分别采用最佳模拟、不确定性上界，不确定性下界，分别模拟对应径流；

步骤64、分别计算选取最佳模拟、不确定性上界，不确定性下界的参数集，各因子对径流变化的贡献：

绘制各个子流域的各因子对径流演变贡献空间分布图，描述大区域径流演变归因空间规律。对于每个子流域、每个因子，不确定性上界，不确定性下界的参数集分别求得的贡献度所夹区间，即为参数不确定性演化的归因不确定性区间。

Claims (8)

1.一种基于大区域水文模拟的径流演变不确定归因方法，其特征在于，包括以下步骤：

步骤(1)建立数据集；

步骤(2)描述历史径流演变规律；

步骤(3)建立大区域水文模型；

步骤(4)模型参数敏感性分析；

步骤(5)模型参数率定及不确定性分析；

步骤(6)径流演变不确定归因分析。

2.根据权利要求1所述的一种基于大区域水文模拟的径流演变不确定归因方法，其特征在于：所述步骤(1)具体为：

建立包含水文气象数据、空间数据和地理特征数据的数据集，包括：水文站径流监测数据、气象站降水、最高气温、最低气温、风速、相对湿度、太阳辐射数据，地形地貌数据，土地利用与覆被数据，土壤类型数据，水库数据，取用水数据。

3.根据权利要求1所述的一种基于大区域水文模拟的径流演变不确定归因方法，其特征在于：所述步骤(2)具体为：

采用Mann-Kendall非参数检验和Spearman’s Rho检验两种方法描述历史径流的演变特征并相互佐证；采用Pettitt检验和有序聚类方法确定径流发生统计意义上显著突变的时间点，划分基准期和变化期，该时间点前的基准期径流受变化环境影响统计意义上不显著，可反映天然情况。

4.根据权利要求3所述的一种基于大区域水文模拟的径流演变不确定归因方法，其特征在于：所述步骤(2)具体过程为：

采用Mann-Kendall非参数检验和Spearman’s Rho检验两种方法描述历史径流的演变特征并相互佐证，对于时间序列X，样本容量为n，构造Mann-Kendall检验统计量S和Z_c：

式中，x_d和x_c为时间序列X的要素，当Z_c＞0时，表明时间序列X随时间有增加趋势，Z_c＜0时表明时间序列X随时间有减少趋势，当|Z_c|＞Z_1-α/2时，Z_1-α/2为标准正态离差，α为显著性水平，表明时间序列变化趋势统计意义显著；

构造Spearman’s Rho检验统计量ρ_s：

式中，x_i为样本，t_i为时间，为样本均值，为时间均值；当样本容量n>30时，统计量ρ_s为正态分布，记当Z_c＞0时，表明时间序列X随时间有增加趋势，Z_c＜0时表明时间序列X随时间有减少趋势，当|Z_c|＞Z_1-α/2时，Z_1-α/2为标准正态离差，变化趋势统计意义显著；

采用Pettitt检验和有序聚类方法确定径流发生统计意义上显著突变的时间，构造时间t时的Pettitt检验统计量U_t,n：

记统计量U_t,N的最大值为k_τ＝max{U_t,n}，对应的时间τ就是突变点，显著性水平检验公式如下：

P＝2exp{-6(k_τ)²/(n³+n²)}

有序聚类发则将原径流时间序列X分成两组样本{xh}(h＝1,2,…,τ)和{xh}(h＝τ+1,2,…,n)，通过样本的统计特性识别序列的突变点，实质是在保持原序列次序的基础上推求最优分割点；最优分割点即使得分割后两样本内部的离差平方和较小，两样本间的离差平方和较大的序列点；记分割后前后两组样本内部的离差平方和为：

式中，V_τ和V_n-τ为离差平方和，和为两组样本均值，总离差平方和为：

S_n(τ)＝V_τ+V_n-τ

满足条件S_n＝min{S_n(τ)}的τ为最优分割点；

按照突变时间前后将研究时期划分为基准期和变化期，认定该时间点前的基准期径流受变化环境影响统计意义上不显著，可反映天然情况。

5.根据权利要求1所述的一种基于大区域水文模拟的径流演变不确定归因方法，其特征在于：所述步骤(3)具体为：

步骤(31)生成水系，划分子流域：依据地形地貌数据生成大区域河网水系结构，包括河网水系、连接点、出水口、入水口；修整河网水系、连接点，添加另外的出入水口或删除不需要的出入水口；设定子流域面积阈值，根据水力联系划分子流域，计算子流域几何参数、地形参数和水流路径，生成子流域报告；

步骤(32)生成水文响应单元：载入坡度数据、土地利用覆被数据、土壤类型数据，设定归类阈值，将阈值内的坡度、土地利用覆被、土壤类型认定为同一类，同时具有相同类坡度、土地利用覆被、土壤类型的单元为水文响应单元，是模型结构的最小单元，也是水文平衡计算的参照单元，生成水文响应单元报告；

步骤(33)输入数据试运行：创建模型数据库，包括配置文件、土壤数据、天气发生数据、子流域数据、水文响应单元数据、主要河道数据、地下水数据、取用水数据、管理数据、水库数据、流域数据，试运行模型。

6.根据权利要求1所述的一种基于大区域水文模拟的径流演变不确定归因方法，其特征在于：所述步骤(4)具体为：

模型参数敏感性分析：采用LH-OAT(LatinHypercube One-factor-at-a-time)敏感性分析技术，在模型众多参数中选出重要敏感的参数；首先将P个参数的原始范围分成N段，对其进行N次拉丁超立方抽样，通过扰动抽样点P次，计算敏感度，对于第a个参数在第b

个抽样点处的敏感度S_a,b(％)可以由下式计算：

式中，M(·)代表模型方程，f_a代表参数扰动，e_a,b代表第a个参数在第b个抽样点处的取值；参数的敏感性可以依据每个参数的敏感度均值排序，选取最敏感的；为了反映大区域模型参数的空间异质性，对每个观测控制点以上的子流域采用一套敏感性参数集，整个大区域采用多套敏感性参数集。

7.根据权利要求1所述的一种基于大区域水文模拟的径流演变不确定归因方法，其特征在于：所述步骤(5)具体为：

所述步骤5中，针对各流域(每个观测控制点以上的子流域集合)的敏感参数集进行参数率定和不确定分析；大区域的水文模拟包含了很多不确定性和复杂性，这些不确定性和复杂性可以反映在模型参数上；本发明采用贝叶斯理论进行参数不确定性分析，采用贝叶斯形式，模型参数的后验分布可以表述为：

式中，θ为模型参数，为径流量观测值，为模型参数的后验分布，p(θ)为模型参数的先验分布，为似然函数，为常数，进而模型参数的后验分布可以表述为：

假定残差是不相关的，那么似然函数可以写成：

式中，y₁(θ)、y₂(θ)、y_n(θ)为模拟径流，为中元素，f_a(b)指基于b对a概率密度函数的估计，进一步假定残差服从正态分布那么似然函数可以写成：

模型参数的先验分布在其初始范围中服从均匀分布，一旦先验分布和似然函数定义完成，就可以推断模型参数的后验分布；采用马尔可夫链蒙特卡罗方法(MCMC)来获得后验分布的样本，通过应用DREAM(DiffeRential Evolution Adaptive Metropolis)算法实现；具体分为以下几个子步骤：

步骤51、认为参数在原始范围内服从均匀分布分布，生成候选样本：

式中θ为模型参数，和为迭代过程中生成的候选样本，k为迭代过程的某一次，γ(δ)、Z₁(a)、Z₂(b)为DREAM算法定义的参数，残差ε～N_d(0,b)，δ、a、b为DREAM算法定义的参数，b为一个小的常数；

步骤52、根据交叉概率CR替换候选样本：

式中，U∈[0,1]服从均匀分布，CR为交叉概率，d_eff为DREAM算法定义的参数；

步骤53、计算候选样本的Metropolis接受概率：

式中Y_obs为径流观测值；

步骤54、如果接受则向前移动，不接受仍保持原位；

重复步骤51-54，直到达到所需规模；

步骤55、采用Gelman收敛诊断指标判断抽样序列是否收敛，当时，收敛，并认定收敛后的样本分布稳定；

步骤56、选取收敛后的敏感参数样本作为敏感参数集的后验分布抽样，采用构建好的模型模拟所以敏感参数样本取值下的径流量，选取满足E_NS＝max{E_NS(i)}，即E_NS最大的一次作为模型的最佳模拟，E_NS如下式计算：

式中，Q_oi为第i时段的实测径流，为各时段实测径流的平均值，Q_si为第i时段的模拟模拟；

E_NS用以描述模型的拟合程度，除E_NS外，采用确定性系数R²评价模拟径流和实测径流的线性相关程度，采用相对误差Re(％)评价模拟结果的偏差，计算公式分别为：

式中，为各时段模拟径流的平均值，其他含义同上；

引入95％不确定性区间描述不确定性，定义为敏感参数集按后验分布抽样升序或降序排列，由2.5％和97.5％分位点所夹区间的参数样本模拟得到的径流区间，采用P-factor和R-factor两个指标来模拟的不确定性，P-factor定义为实测径流被包含在95％不确定性区间中的百分比；R-factor定义为95％不确定性区间平均宽度比标准差，计算公式为：

式中Q_U和Q_L分别为95％不确定性区间的上界和下界，σ为实测径流的标准差；P-factor越接近1，R-factor越接近0，结果越好，当R-factor小于1时，可以认为模拟结果满足要求。

8.根据权利要求1所述的一种基于大区域水文模拟的径流演变不确定归因方法，其特征在于所述步骤6进一步分为：

步骤61、筛选径流演变影响因子，大区域水文模型的输入包括，气象输入：降水、降水、最高气温、最低气温、风速、相对湿度、太阳辐射数据，人类活动输入：土地利用与覆被数据，水库数据，取用水数据；分析突变时间前后，输入因子的统计特征值，筛选出发生统计意义显著变化的因子作为归因分析的影响因子；

步骤62、根据筛选出的影响因子生成情景集；

步骤63、将情景集输入模型，各子流域敏感参数集分别采用最佳模拟、不确定性上限，不确定性下限，分别模拟对应径流；

步骤64、计算各因子对径流变化的贡献及参数不确定性演化的归因不确定性区间；

为了更清楚地阐述步骤62-64，给出以下2因子、4因子、N因子算例；当仅含2项影响因子F₁、F₂时，每项因子具有2个状态，分为“-1”基准期状态，“1”变化期状态，情景组合方案可以设置为情景S₁至S₄：

为了书写简便记因子F₁对径流的影响量可以表示为：

或

因子F₂对径流的影响量可以表示为：

或

这就是归因路径的多样性问题，不同的路径计算得到的因子贡献量是不相等的，因为真实情况下自然界的影响因子间是相互不独立的，在不同状态的F₂下，F₁对径流的影响量是不一样的；为了规避归因路径的多样性带来的不确定性，本发明定义影响因子F₁、F₂对径流变化的贡献量分别为：

将影响因子对径流变化的贡献量定义为各路径的平均值，这样求得某因子的贡献量为其他因子平均状态下的贡献量；因此4因子情形下，情景组合方案可以设置为：

影响因子F₁、F₂、F₃、F₄对径流变化的贡献量按前述定义分别为：

进一步推广到N因子情形，第q个因子对径流变化的贡献量可以表述为：

式中，α_q,p为权重系数，表示因子Fq在情景p下的状态，等于-1或1，具体取值可以在情景组合方案表中查得。