CN109408848B - 一种考虑径流演变时空异质性的分布式归因方法 - Google Patents

Abstract

本发明公开了一种考虑径流演变时空异质性的分布式归因方法，包括如下步骤：建立分布式水文模型，对研究区域按子流域进行天然径流还原计算；根据径流演变特征划分研究时期；分析天然径流时空分布与影响因子相关性，构建重要影响因子集；对天然径流演变逐时段、逐子流域、逐因子滚动式归因；对实测径流演变逐时段、逐片区、逐因子滚动式归因。本发明根据径流演变的时程特征，将研究时期细分为多段，根据空间地形地貌特征将研究区域细分到子流域，分别从天然和实测径流入手，将径流演变在时间和空间的维度上归因到影响因子，获得的归因结果充分考虑了径流演变的时空异质性。

Description

一种考虑径流演变时空异质性的分布式归因方法

技术领域

本发明设计水利工程领域中的径流演变归因技术，特别涉及一种考虑径流演变时空异质性的分布式归因方法。

背景技术

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

径流演变的归因技术旨在定量分析径流时空变化的成因，为预测未来径流、制定水资源管理的适应性对策提供依据，针对不同的成因和影响程度做出具有针对性、便于操作的适应性调控决策。国内外变化环境下径流变化的归因研究取得了长足的发展，但从当前见诸文献的成果看，目前的归因技术通常将研究时期分为天然期和变化期两段，给出整个研究区域的归因结果，对径流演变的时空异质性考虑尚显不足。

发明内容

发明目的：提供一种考虑径流演变时空异质性的分布式归因方法，以解决现有技术存在的上述问题。

技术方案：一种考虑径流演变时空异质性的分布式归因方法，包括以下步骤：

步骤(1)建立分布式水文模型，对研究区域按子流域进行天然径流还原计算；具体为:

依据研究区域地形地貌数据生成水系、划分子流域，基于地形地貌数据、土地利用覆被数据、土壤类型数据生成水文响应单元，从而建立起分布式物理水文模型SWAT模型结构。采用LH-OAT(Latin Hypercube One-factor-at-a-time)敏感性分析技术，在模型众多参数中选出重要敏感的参数。通过率定敏感参数，获得满足评价要求的模型，进行天然径流还原计算。

步骤(2)根据径流演变特征划分研究时期；具体为:

初始的长时间序列往往是由多段具有明显变化特征的中短时间序列组合而成，而总体的时程演变规律并不明显。因此，在传统分析的基础上对初始时间序列进行分段，将其划分为若干个径流子序列，然后分段分析各径流序列的时程演变特征。

步骤(3)分析天然径流时间、空间分布与影响因子时间、空间分布相关性，构建重要影响因子集；具体为:

根据气象站气象因子时间序列计算子流域气象因子面平均序列；从土地利用/覆被资料中提取子流域的土地利用类型面积。对各个时段，分析影响因子绝对值与变化值的空间分布状况，采用独立性检验的方法评价天然径流时空变化与影响因子的相关程度。参照独立性检验的结果，确定径流时空分布与影响因子之间的相关性程度，选取高相关性的因子组成重要影响因子集。

步骤(4)对天然径流演变逐时段、逐子流域、逐因子滚动式归因；具体为:

述步骤4中，对天然径流演变逐时段、逐子流域、逐因子滚动式归因。将全流域共划分为n个子流域，研究期共包含m个时段，重要影响因子集中共有n_F个因子，将第j个子流域(j＝1,2,...,n)在第i个时段(i＝1,2,...,m-1)内的天然径流量记为

在第i+1个时段(i＝1,2,...,m-1)内的天然径流量记为

令第i个时段重要影响因子集中其他n_F-1个影响因子保持不变，将影响因子F替换为第i+1个时段的相应数值，输入流域SWAT模型模拟天然径流，将第j个子流域的径流量记为

则可计算第j个子流域从第i个时段到第i+1个时段的天然径流变化量

影响因子F对天然径流变化的贡献

第j个子流域从第i个时段到第i+1个时段的天然径流变化程度

影响因子F变化引起天然径流变化程度

为：

各时段内的流域径流特性存在空间差异，各子流域的径流量分属不同的大小级别，各级别的径流在时段间演变时所受各影响因子的影响不尽相同，展开不同大小级别的径流演变归因分析；各时段间的流域径流演变特性也存在空间差异，各子流域的径流变化程度分属不同级别，各变化程度的径流在演变时所受各影响因子的影响不尽相同，展开不同变化程度的径流演变归因分析。

步骤(5)对实测径流演变逐时段、逐片区、逐因子滚动式归因,具体为:

流域实测径流在数量上等于天然径流扣除人类直接取用水的值，其重要影响因子为天然径流重要影响因子集和人类直接取用水。依据研究区域N个片区的水文监测资料计算出每个片区的实测径流。天然径流演变按照步骤4逐时段、逐片区、逐因子滚动式归因，人类取用水因素对片区实测径流演变的贡献

和贡献度

计算如下：

有益效果：通过建立基于物理机理的分布式水文模型，用以还原计算研究区域天然径流；依据地形地貌划分子流域，依据径流时程演变规律划分子时段，后续的径流演变归因都基于此；分析天然径流时空分布与影响因子相关性，构建了重要影响因子集；通过对天然径流演变逐时段、逐子流域、逐因子滚动式归因、对实测径流演变逐时段、逐片区、逐因子滚动式归因，在时间和空间两个维度给出了归因结果，充分考虑了径流演变的时空异质性。

附图说明

图1为径流演变时空异质性的示意图。

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

图3为研究时期划分时段方法的流程图。

图4为泰森多边形构建的示意图。

图5为子流域K-均值聚类的流程图。

具体实施方式

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

本发明提供一种考虑径流演变时空异质性的分布式归因方法，包括以下步骤：

步骤(1)建立分布式水文模型，对研究区域按子流域进行天然径流还原计算；具体步骤如下：

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

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

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

步骤(14)模型参数敏感性分析：采用LH-OAT(Latin Hypercube One-factor-at-a-time)敏感性分析技术，将模型众多参数按照敏感度进行排序；首先将P个参数的原始范围分成N_par段，对其进行N_par次拉丁超立方抽样，通过扰动抽样点P次，计算敏感度；对于第a个参数在第b个抽样点处的敏感度S_a,b(％)可以由下式计算：

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

步骤(15)率定敏感性参数集：选取满足E_NS＝max{E_NS(t)}，即E_NS最大的一次作为模型的最佳模拟，E_NS如下式计算：

式中，Q_ot为第t时刻的实测径流，

为各时段实测径流的平均值，Q_st为第t时刻的模拟模拟，T为序列长度；

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

式中，

为各时段模拟径流的平均值，Q_ot为第t时刻的实测径流，

为各时段实测径流的平均值，Q_st为第t时刻的模拟，T为序列长度。

步骤(2)根据径流演变特征划分研究时期；具体步骤如下：

在传统分析的基础上对初始时间序列进行分段，将其划分为若干个径流子序列，然后分段分析各径流序列的时程演变特征；对时间序列进行多层次的分段。将初始时间序列作为一级时间序列，进行Pettitt突变检验，找到一级突变点，于是初始时间序列被划分为突变点前后的两个二级时间序列，以此类推，不断运用Pettitt法寻找逐级时间序列的突变点，在突变点处将时间序列进行分段；通过Mann-Kendall趋势分析的显著性检验，当某时间序列具有显著的趋势变化特征，或序列长度L小于给定的长度阈值L_m，或序列级别数C大于给定的级别阈值C_m时，将不再对该时间序列进行分段；对于时间序列X，序列长度为T，构造时刻t时的Pettitt检验统计量U_t,T：

式中x_d和x_c为时间序列X的要素，记统计量U_t,N的最大值为k_τ＝max{U_t,T}，对应的时间τ就是突变点，显著性水平P检验公式如下：

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

T为序列长度；

对于时间序列X，构造Mann-Kendall检验统计量S和Z_c：

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

步骤(3)分析天然径流时间、空间分布与影响因子时间、空间分布相关性，构建重要影响因子集；具体为：

步骤(31)确定子流域影响因子特征值：根据各气象站气象因子值求解各子流域气象因子值，首先，根据各气象站的位置点据，构建泰森多边形；每个气象站点据分别落在一个泰森多边形中，认为泰森多边形覆盖的面积即为该气象站所控制的面积，由此可得各气象站的控制范围；用加权平均的方法，由各站气象因子值计算各子流域气象因子值：

式中，C_j为第j个子流域的气象因子值；A_j为第j个子流域的面积；A_js为第s个气象站对应的泰森多边形在第j个子流域中的面积；C_s为第s个气象站的气象因子值；S为流域内的气象站个数；n为流域内的子流域个数；

各子流域土地利用类型可从全流域土地利用数据空间分布图中直接提取各子流域的土地利用类型构成，从而得到提取子流域的土地利用类型面积，根据各子流域在各时段内的影响因子值，可计算各时段间各影响因子的变化值；

步骤(32)采用K-均值聚类法，将子流域划分为确定数目的类。基本思想为：随机选取N_K个初始的聚类中心，将子流域中的每个影响因子特征值分配到与之距离最近的聚类中心所属的类中，然后不断重新计算各个类的聚类中心并重合进行聚类，直到聚类中心和聚类结果收敛为止，即最终的聚类结果满足总聚类平方偏差和最小；在根据聚类结果选取新的聚类中心时，采用的方法是计算各类子流域的均值，公式如下：

式中，N_K为类的数目；D为子流域特性的维数；N_l为第l个类中的子流域数目；y_l为第l个类的聚类中心，y_l,r为其第r维的值；b_lj为第l个类中的第j个子流域，b_lj,r为其第r维的值；根据y_l,r的计算值，可以得到各聚类中心的表达形式y_l；

基于流域内各子流域的影响因子特征值，依照各因子在各时段内和时段间的特性进行子流域K-均值聚类；根据各类子流域的位置分布情况，可概括各时段内影响因子的空间分布特征以及各时段间影响因子变化程度的空间分布特征；

步骤(33)通过独立性检验，分析径流时空演变和影响性子的相关性：在分析径流空间分布与影响因子的相关性时，基于各时段内的静态分布格局和各时段间的动态变化格局分别展开独立性检验，将流域内的全部子流域作为样品，将径流和影响因子作为分类变量，设研究期共包含m个时段，独立性检验的对象包括：(1)第i个时段(i＝1,2,...,m)内径流特性与影响因子特性的独立性检验；(2)第i个时段与第i+1个时段(i＝1,2,...,m-1)间径流变化特性与影响因子变化特性的独立性检验；根据上述n_F×(2m-1)组独立性检验的结果，可对径流空间分布与各影响因子的相关性程度做出综合评价。

步骤(4)对天然径流演变逐时段、逐子流域、逐因子滚动式归因；具体如下：

步骤(41)对天然径流演变逐时段、逐子流域、逐因子滚动式归因：将全流域共划分为n个子流域，研究期共包含m个时段，重要影响因子集中共有n_F个因子，将第j个子流域(j＝1,2,...,n)在第i个时段(i＝1,2,...,m-1)内的天然径流量记为

在第i+1个时段(i＝1,2,...,m-1)内的天然径流量记为

令第i个时段重要影响因子集中其他n_F-1个影响因子保持不变，将影响因子F替换为第i+1个时段的相应数值，输入流域SWAT模型模拟天然径流，将得到的模拟径流量记为

则可计算第j个子流域从第i个时段到第i+1个时段的天然径流变化量

影响因子F对天然径流变化的贡献

第j个子流域从第i个时段到第i+1个时段的天然径流变化程度

影响因子F变化引起天然径流变化程度

为：

各时段内的流域径流特性存在空间差异，各子流域的径流量分属不同的大小级别，各级别的径流在时段间演变时所受各影响因子的影响不尽相同，展开不同大小级别的径流演变归因分析，设在第i个时段(i＝1,2,...,m-1)内各子流域的径流特性共分为Z类，各类之间的区别在于径流量大小不同，对Z类径流分别进行演变格局中所呈现特征的归因分析，以受影响因子F贡献最大的子流域的面积和为度量指标，确定影响因子F对其影响程度，计算公式如下：

其中，

式中，

为第i时段内第z类大小级别的径流在接下来的径流演变过程中受F因子主导的比例；A_j为第j个子流域的面积；n为流域内子流域的个数；m为研究期内包含的时段数；Z为第i个时段内基于径流大小级别的子流域聚类数目；

各时段间的流域径流演变特性也存在空间差异，各子流域的径流变化程度分属不同级别，各变化程度的径流在演变时所受各影响因子的影响不尽相同，展开不同变化程度的径流演变归因分析；设从第i个时段到第i+1个时段(i＝1,2,...,m-1)各子流域的径流变化特性共分为G类，各类之间的区别在于径流量变化程度的不同；对G类径流分别进行演变格局中所呈现特征的归因分析，以受影响因子F贡献最大的子流域的面积和为度量指标，确定影响因子F对其影响程度，计算公式如：

其中：

式中，

为从第i时段到第i+1时段第g类变化程度的径流在其演变过程中受影响因子F主导的比例；A_j为第j个子流域的面积；n为流域内子流域的个数；m为研究期内包含的时段数；G为从第i时段到第i+1时段基于径流变化程度的子流域聚类数目。

步骤(5)对实测径流演变逐时段、逐片区、逐因子滚动式归因；具体如下：

对实测径流演变逐时段、逐片区、逐因子滚动式归因：流域实测径流在数量上等于天然径流扣除人类直接取用水的值，其重要影响因子为天然径流重要影响因子集和人类直接取用水，根据流域内代表性水文站的控制范围，将整个流域划分为N个不重叠的片区，每个片区的出口处各有1个水文站，每个片区实测径流的计算公式如下所示：

式中，

为第p个片区在第i时段内的实测径流；

和

为第p和k个水文站在第i时段内的实测径流；A_p和A_k为第p和k个水文站控制区域的面积；m为研究期内包含的时段数；N为选定的水文站个数，等于片区数；

于是，从第i时段到第i+1时段，第p个片区实测径流的变化量

可按照下式进行计算：

为衡量人类直接取用水对实测径流演变的贡献，需计算在不受人类取用水影响时各片区天然径流的变化量，基于由流域SWAT模型输出的各子流域天然径流，按照下式计算各片区的天然径流：

式中，

为第p个片区在第i时段内的天然径流；

为第j个子流域在第i时段内的天然径流；

为第j个片区的面积；A_pj为第j个子流域在第p个片区内的面积；m为研究期内包含的时段数；N为片区个数；n为子流域个数。

从第i时段到第i+1时段，第p个片区天然径流的变化量

可按照下式进行计算：

在步骤4中已经将天然径流演变归因到重要影响因子集，计算了各重要影响因子对天然径流演变的贡献，将其折算到各片区，计算如下：

式中：

为从第i时段到第i+1时段第p个片区内由影响因子F引起的天然径流变化量；其余变量的含义与上文一致；

人类取用水因素对片区实测径流演变的贡献

和贡献度

计算如下：

Claims (7)

1.一种考虑径流演变时空异质性的分布式归因方法，其特征在于，包括如下步骤：

步骤(1)建立分布式水文模型，对研究区域按子流域进行天然径流还原计算；

步骤(2)根据径流演变特征划分研究时期；

步骤(3)分析天然径流时间、空间分布与影响因子时间、空间分布相关性，构建重要影响因子集；

步骤(4)对天然径流演变逐时段、逐子流域、逐因子滚动式归因；

步骤(5)对实测径流演变逐时段、逐片区、逐因子滚动式归因；

所述步骤(4)具体为：

步骤(41)对天然径流演变逐时段、逐子流域、逐因子滚动式归因：将全流域共划分为n个子流域，研究期共包含m个时段，重要影响因子集中共有n_F个因子，将第j个子流域，j＝1,2,…,n，在第i个时段，i＝1,2,…,m-1，内的天然径流量记为

在第i+1个时段，i＝1,2,…,m-1，内的天然径流量记为

令第i个时段重要影响因子集中其他n_F-1个影响因子保持不变，将影响因子F替换为第i+1个时段的相应数值，输入流域SWAT模型模拟天然径流，将得到的模拟径流量记为

则计算第j个子流域从第i个时段到第i+1个时段的天然径流变化量

影响因子F对天然径流变化的贡献

第j个子流域从第i个时段到第i+1个时段的天然径流变化程度

影响因子F变化引起天然径流变化程度

为：

各时段内的流域径流特性存在空间差异，各子流域的径流量分属不同的大小级别，各级别的径流在时段间演变时所受各影响因子的影响不尽相同，展开不同大小级别的径流演变归因分析，设在第i个时段，i＝1,2,…,m-1，内各子流域的径流特性共分为Z类，各类之间的区别在于径流量大小不同，对Z类径流分别进行演变格局中所呈现特征的归因分析，以受影响因子F贡献最大的子流域的面积和为度量指标，确定影响因子F对其影响程度，计算公式如下：

其中，

式中，

为第i时段内第z类大小级别的径流在接下来的径流演变过程中受F因子主导的比例；A_j为第j个子流域的面积；n为流域内子流域的个数；m为研究期内包含的时段数；Z为第i个时段内基于径流大小级别的子流域聚类数目；

各时段间的流域径流演变特性也存在空间差异，各子流域的径流变化程度分属不同级别，各变化程度的径流在演变时所受各影响因子的影响不尽相同，展开不同变化程度的径流演变归因分析；设从第i个时段到第i+1个时段，i＝1,2,…,m-1，各子流域的径流变化特性共分为G类，各类之间的区别在于径流量变化程度的不同；对G类径流分别进行演变格局中所呈现特征的归因分析，以受影响因子F贡献最大的子流域的面积和为度量指标，确定影响因子F对其影响程度，计算公式如：

其中：

式中，

为从第i时段到第i+1时段第g类变化程度的径流在其演变过程中受影响因子F主导的比例；A_j为第j个子流域的面积；n为流域内子流域的个数；m为研究期内包含的时段数；G为从第i时段到第i+1时段基于径流变化程度的子流域聚类数目。

2.根据权利要求1所述的一种考虑径流演变时空异质性的分布式归因方法，其特征在于：所述步骤(1)具体包括如下步骤：

步骤(11)依据研究区域地形地貌数据生成水系，划分子流域；

步骤(12)基于地形地貌数据、土地利用覆被数据、土壤类型数据生成水文响应单元，从而建立起分布式物理水文模型SWAT模型结构；

步骤(13)输入数据试运行模型；

步骤(14)采用LatinHypercube One-factor-at-a-time敏感性分析技术，在模型众多参数中选出重要敏感的参数；

步骤(15)通过率定敏感参数，获得满足评价要求的模型，进行天然径流还原计算。

3.根据权利要求1或2所述的一种考虑径流演变时空异质性的分布式归因方法，其特征在于：所述步骤(1)具体为：

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

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

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

步骤(14)模型参数敏感性分析：采用LatinHypercube One-factor-at-a-time敏感性分析技术，将模型众多参数按照敏感度进行排序；首先将P个参数的原始范围分成N_par段，对其进行N_par次拉丁超立方抽样，通过扰动抽样点P次，计算敏感度；对于第a个参数在第b个抽样点处的敏感度S_a,b由下式计算：

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

步骤(15)率定敏感性参数集：选取满足E_NS＝max{E_NS(t)}，即E_NS最大的一次作为模型的最佳模拟，E_NS如下式计算：

式中，Q_ot为第t时刻的实测径流，

为各时段实测径流的平均值，Q_st为第t时刻的模拟径流，T为序列长度；

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

式中，相对误差Re为百分比，

为各时段模拟径流的平均值，Q_ot为第t时刻的实测径流，

为各时段实测径流的平均值，Q_st为第t时刻的模拟，T为序列长度。

4.根据权利要求1所述的一种考虑径流演变时空异质性的分布式归因方法，其特征在于：所述步骤(2)具体为：

在传统分析的基础上对初始时间序列进行分段，将其划分为若干个径流子序列，然后分段分析各径流序列的时程演变特征；对时间序列进行多层次的分段；将初始时间序列作为一级时间序列，进行Pettitt突变检验，找到一级突变点，于是初始时间序列被划分为突变点前后的两个二级时间序列，以此类推，不断运用Pettitt法寻找逐级时间序列的突变点，在突变点处将时间序列进行分段；通过Mann-Kendall趋势分析的显著性检验，当某时间序列具有显著的趋势变化特征，或序列长度L小于给定的长度阈值L_m，或序列级别数C大于给定的级别阈值C_m时，将不再对该时间序列进行分段；对于时间序列X，序列长度为T，构造时刻t时的Pettitt检验统计量U_t,T：

式中x_d和x_c为时间序列X的要素，记统计量U_t,N的最大值为k_τ＝max{U_t,T}，对应的时间τ就是突变点，显著性水平P检验公式如下：

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

T为序列长度；

对于时间序列X，构造Mann-Kendall检验统计量S和Z_c：

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

5.根据权利要求1所述的一种考虑径流演变时空异质性的分布式归因方法，其特征在于：所述步骤(3)具体为：

步骤(31)根据气象站气象因子时间序列计算子流域气象因子面平均序列；从土地利用/覆被资料中提取子流域的土地利用类型面积；将子流域气象因子和子流域的土地利用类型面积作为径流变化的影响因子；

步骤(32)对各个时段，依照影响因子时空分布特征对各子流域进行K-均值聚类，分析影响因子绝对值与变化值的空间分布状况；

步骤(33)采用独立性检验的方法评价天然径流时空变化与影响因子的相关程度；参照独立性检验的结果，确定径流时空分布与影响因子之间的相关性程度，选取高相关性的影响因子组成重要影响因子集。

6.根据权利要求1或5所述的一种考虑径流演变时空异质性的分布式归因方法，其特征在于：所述步骤(3)具体为：

步骤(31)确定子流域影响因子特征值：根据各气象站气象因子值求解各子流域气象因子值，首先，根据各气象站的位置点据，构建泰森多边形；

每个气象站点据分别落在一个泰森多边形中，认为泰森多边形覆盖的面积即为该气象站所控制的面积，由此得到各气象站的控制范围；用加权平均的方法，由各站气象因子值计算各子流域气象因子值：

式中，C_j为第j个子流域的气象因子值；A_j为第j个子流域的面积；A_js为第s个气象站对应的泰森多边形在第j个子流域中的面积；C_s为第s个气象站的气象因子值；S为流域内的气象站个数；n为流域内的子流域个数；

各子流域土地利用类型是从全流域土地利用数据空间分布图中直接提取各子流域的土地利用类型构成，从而得到提取子流域的土地利用类型面积，根据各子流域在各时段内的影响因子值，计算各时段间各影响因子的变化值；

步骤(32)采用K-均值聚类法，将子流域划分为确定数目的类；思想为：随机选取N_K个初始的聚类中心，将子流域中的每个影响因子特征值分配到与之距离最近的聚类中心所属的类中，然后不断重新计算各个类的聚类中心并重合进行聚类，直到聚类中心和聚类结果收敛为止，即最终的聚类结果满足总聚类平方偏差和最小；在根据聚类结果选取新的聚类中心时，采用的方法是计算各类子流域的均值，公式如下：

式中，N_K为类的数目；D为子流域特性的维数；N_l为第l个类中的子流域数目；y_l为第l个类的聚类中心，y_l,r为其第r维的值；b_lj为第l个类中的第j个子流域，b_lj,r为其第r维的值；根据y_l,r的计算值，得到各聚类中心的表达形式y_l；

基于流域内各子流域的影响因子特征值，依照各因子在各时段内和时段间的特性进行子流域K-均值聚类；根据各类子流域的位置分布情况，概括各时段内影响因子的空间分布特征以及各时段间影响因子变化程度的空间分布特征；

步骤(33)通过独立性检验，分析径流时空演变和影响性子的相关性：在分析径流空间分布与影响因子的相关性时，基于各时段内的静态分布格局和各时段间的动态变化格局分别展开独立性检验，将流域内的全部子流域作为样品，将径流和影响因子作为分类变量，设研究期共包含m个时段，独立性检验的对象包括：(1)第i个时段，i＝1,2,…,m，内径流特性与影响因子特性的独立性检验；(2)第i个时段与第i+1个时段，i＝1,2,…,m-1，间径流变化特性与影响因子变化特性的独立性检验；根据n_F×(2m-1)组独立性检验的结果，对径流空间分布与各影响因子的相关性程度做出综合评价。

7.根据权利要求1所述的一种考虑径流演变时空异质性的分布式归因方法，其特征在于：所述步骤(5)具体为：

对实测径流演变逐时段、逐片区、逐因子滚动式归因：流域实测径流在数量上等于天然径流扣除人类直接取用水的值，其重要影响因子为天然径流重要影响因子集和人类直接取用水，根据流域内代表性水文站的控制范围，将整个流域划分为N个不重叠的片区，每个片区的出口处各有1个水文站，每个片区实测径流的计算公式如下所示：

式中，

为第p个片区在第i时段内的实测径流；

和

为第p和k个水文站在第i时段内的实测径流；A_p和A_k为第p和k个水文站控制区域的面积；m为研究期内包含的时段数；N为选定的水文站个数，等于片区数；

于是，从第i时段到第i+1时段，第p个片区实测径流的变化量

按照下式进行计算：

为衡量人类直接取用水对实测径流演变的贡献，需计算在不受人类取用水影响时各片区天然径流的变化量，基于由流域SWAT模型输出的各子流域天然径流，按照下式计算各片区的天然径流：

式中，

为第p个片区在第i时段内的天然径流；

为第j个子流域在第i时段内的天然径流；

为第j个片区的面积；A_pj为第j个子流域在第p个片区内的面积；m为研究期内包含的时段数；N为片区个数；n为子流域个数；

从第i时段到第i+1时段，第p个片区天然径流的变化量

按照下式进行计算：

在步骤4中已经将天然径流演变归因到重要影响因子集，计算了各重要影响因子对天然径流演变的贡献，将其折算到各片区，计算如下：

式中：

为从第i时段到第i+1时段第p个片区内由影响因子F引起的天然径流变化量；其余变量的含义与上文一致；

人类取用水因素对片区实测径流演变的贡献

和贡献度

计算如下：

Images