一种考虑土地利用时空演变的径流变化归因识别方法与流程

本发明涉及水利工程领域中的流域径流变化归因识别技术，特别涉及一种考虑土地利用时空演变的径流变化归因识别方法。

背景技术：

近百年来，受全球气候变暖影响，水文循环加剧，极端径流事件频度、强度均有一定程度的增加。社会经济活动的持续增强和高密度水利工程群的建设，形成了复杂的人类活动影响链，研究表明众多人类活动中，土地利用变化改变了下垫面条件，从而改变了产汇流条件，对径流变化起到了尤为明显的作用。而长期以来，流域水资源开发利用所依据的历史径流序列和供需关系主要基于水文序列的平稳性假定，普遍缺少对气候变化和人类活动影响及其带来的不确定性研究。变化环境造成径流时空演变规律发生显著变异，并影响区域社会经济发展的水资源供需关系，从而给复杂水利水电系统的安全、经济运行带来风险。

径流变化的归因识别技术旨在定性和定量揭示变化环境下径流时空变化的成因，为流域水资源系统适应性调度和管理提供依据，通过定性或定量分析的归因识别结果做出具有针对性、便于操作的径流适应性调控决策。径流变化的归因识别不是一件简单的工作，目前还没有一个完善的模型或方法，主要采用数理统计分析和水文模型模拟等方法，水文模型模拟方法经历了从集总式模型向分布式模型的发展。由于分布式水文模型能够描述流域空间信息，被广泛应用于径流变化归因识别研究，通常的做法是用一期的土地利用资料代替一整段时期建模，无法考虑土地利用的时空演变，尤其是在人类活动剧烈的变化期。

技术实现要素：

发明目的：提供一种考虑土地利用时空演变的径流变化归因识别方法，以解决现有技术存在的上述问题。

技术方案：一种考虑土地利用时空演变的径流变化归因识别方法，包括以下步骤：

步骤1，根据径流突变点将研究时期划分为天然期和变化期；

步骤2，根据土地利用演变特征将变化期分成m段；

步骤3，建立流域swat模型；

步骤4，根据天然期资料率定流域swat模型参数；

步骤5，对流域径流变化逐段逐项滚动式归因定量计算。

根据本发明的一个方面：所述步骤1中，采用pettitt检验方法对流域径流时间序列进行突变检验，确定径流发生的变化高于预期值的时间，将整个研究时期划分为天然时期和变化时期，认定天然期受气候变化和人类活动的影响小于预期值，可反映该流域天然径流状况。

根据本发明的一个方面：所述步骤2中，统计不同时刻的各土地利用类型面积，反映各类型在数量上的变化，制作土地利用转移矩阵，反映不同土地类型在空间分布上的互相转换，根据土地利用的时空演变特征将变化期分成n段，分段归因以充分考虑不同阶段土地利用对径流的影响，n为自然数。

根据本发明的一个方面：通过arcgis平台下的arcswat模块建立流域swat模型，需要基于流域dem数据资料生成水系并进行子流域的划分，根据土地利用类型、土壤类型和坡度类型的空间分布情况生成水文响应单元，输入气象数据，运行模型。

根据本发明的一个方面：在进行参数敏感性分析后，排除众多参数中对模型结果影响低于预期值的参数，选取敏感性大于预定值的参数进行率定，具体的参数率定工作，利用swat-cup软件工具完成，得到了参数率定值之后，再对率定的效果进行评价。

根据本发明的一个方面：所述步骤5进一步分为：

步骤51、保持阶段i-1的下垫面条件不变，将阶段i的气象资料输入swat模型，则阶段i气候变化对径流的影响可以表示为：

δrc,i＝r(li-1,ci)-r(li-1,ci-1)(1)

式中，δrc,i为阶段i气候变化引起的径流变化；li-1为阶段i-1土地利用资料，ci-1和ci分别为阶段i-1和阶段i的气象资料；r(li-1,ci)和r(li-1,ci-1)分别表示固定阶段i-1的土地利用资料输入不变，分别采用阶段i和阶段i-1的气象资料，经swat模拟得到的径流；

步骤52、保持阶段i气候条件不变，模拟土地利用变化前后径流的变化，以表示阶段i土地利用对径流的影响：

δrl,i＝r(li,ci)-r(li-1,ci)(2)

式中，δrl,i为阶段i土地利用变化引起的径流变化，即间接人类活动对径流的影响；r(li,ci)和r(li-1,ci)分别表示固定阶段i的气象资料输入不变，分别采用阶段i和阶段i-1的土地利用资料，经swat模拟得到的径流；

步骤53、定量解析出气候变化和下垫面变化这种间接的人类活动对径流的影响后，剩下的径流变化近似归纳到流域内生活生产和流域外引水调水这种直接的人类活动：

δrd,i＝δrt,i-δrc,i-δrl,i＝ro,i-r(li,ci)(3)

式中，δrd,i为阶段i直接人类活动引起的径流变化，δrt,i为阶段i径流总变化，ro,i为阶段i实测径流。

一种考虑土地利用时空演变的径流变化归因识别方法，包括如下步骤：

步骤1，根据径流突变点将研究时期划分为天然期和变化期：采用如下检验方法找到径流时间序列突变点，判断突变点前后的累积分布函数是否存在显著的差异，构造统计量ut,n：

式中：xj，xi为时间序列；t为任一时间节点；sgn(·)为符号函数；n为时间序列长度；

记统计量ut,n的最大值为kτ，对应的样本τ就是突变点，显著性水平检验公示如下：

式中，kn＝max[|kτ|]；

确定径流发生显著变化的时间后，以突变点为界将整个研究时期划分为天然时期和变化时期，并认为天然期受气候变化和人类活动的影响小于预期值，可反映该流域天然径流状况；

步骤2，根据土地利用演变特征将变化期分成n段：搜集土地利用资料，统计不同时刻的各土地利用类型面积，反映各类型在数量上的变化，制作土地利用转移矩阵，反映不同土地类型在空间分布上的互相转换，根据土地利用的时空演变特征将变化期分成具有显著变化的n段，分段归因以充分考虑不同阶段土地利用对径流的影响，此时整个研究时期被划分成了n+1段；

步骤3，通过arcgis平台下的arcswat模块建立流域swat模型，具体分为以下几个子步骤：

步骤31、子流域划分：根据流域dem图层和人为设定的流域阈值，生成河网水系，以及与河网结构相关的若干出水口和入水口，在根据实际需要人为添加另外的出水口和入水口之后，确定流域总出水口，基于水力联系进行子流域划分，并分析计算各子流域的地形参数、几何参数和水流路径；

步骤32、水文响应单元生成：划分子流域后，需要对各子流域划分水文响应单元hru，反映不同下垫面条件下水文循环的特征，hru是一个子流域内具有相同的土地利用类型、土壤类型和坡度类型的区域，是模型运行中的最小地块单元，将土地利用、土壤、坡度类型的空间分布数据集载入模型，根据其组合和分布，可生成若干个hru；

步骤33、数据输入与模型运行：将各类土壤属性数据、气象数据等输入模型，在模型数据库中创建流域输入文件、子流域输入文件、hru输入文件、管理输入文件、土壤输入文件和地下水输入文件，然后进行模型的运行；

模型运行产生的结果主要包括流域内各hru、各子流域和各演算河段的汇总信息，分别存放在hru输出文件、子流域输出文件和主河道输出文件中；对模型运行结果进行分析计算，得出流域各部分的天然径流量；

步骤4，根据天然期资料率定流域swat模型参数，具体分为以下几个子步骤：

步骤41、参数敏感性分析：为了提高参数率定的效率，在进行参数敏感性分析后，排除众多参数中对模型结果影响较小的参数，选取敏感性较强的参数进行率定；

敏感性是一个无量纲的指数，反映了模型结果随参数改变而变化的程度，敏感性结果的计算形式如下：

式中：si为参数i的敏感性结果；o为模型输出结果，其数学表达式为o＝f(f1,f2,…,fn)；fi为影响模型结果的参数；n为参数个数；

敏感性结果通过绝对值、相对值或离差的形式进行表示；为对比不同参数的敏感性，选取相对值的形式，计算各参数的敏感性指数i，公式如下：

步骤42、参数率定：采用swat-cup软件工具完成，选择sufi-2算法，该算法考虑了模型结构、输入数据、参数和实测数据的不确定性，通过拉丁超立方随机抽样法随机生成一组参数，代入swat模型中进行模拟，并比较模拟值与实测值的差异，优选参数范围，通过多次迭代运算，不断减小参数的不确定性范围，使模拟值不断与实测值相接近；

步骤43、参数率定评价：选用相关系数r²和nashsuttcliffe模型效率系数nse评价swat模型参数的适用性；

相关系数r²用于评价模拟值和实测值之间的吻合程度，计算公式如下：

式中：qoi为第i时段的实测流量；为各时段实测流量的平均值；qsi为第i时段的模拟流量；为各时段模拟流量的平均值；n为实测数据个数，等于参数率定的时段数；

r²的值常介于0与1之间，愈接近1，说明实测值与模拟值线性相关程度愈高；r²大于0.8时，参数率定结果符合预期水平；

nse是一个整体综合指标，用于定量表示整个径流过程模拟值与实测值的拟合程度，计算公式如下：

式中各变量含义同上，nse的取值范围一般在0～1之间；当nse＝1时，模拟效果最好；当nse接近0或为负值时，模拟值可信度低；nse大于0.7时，参数率定结果符合预期；

步骤5，对流域径流变化逐段逐项滚动式归因定量计算，研究时期已经被分为n+1段，设i＝1,2,...,n，从变化时期第1段开始逐段逐项滚动式计算，具体分为以下几个子步骤：

步骤51、保持阶段i-1的下垫面条件不变，将阶段i的气象资料输入swat模型，则阶段i气候变化对径流的影响可以表示为：

δrc,i＝r(li-1,ci)-r(li-1,ci-1)(11)

式中，δrc,i为阶段i气候变化引起的径流变化；li-1为阶段i-1土地利用资料，ci-1和ci分别为阶段i-1和阶段i的气象资料；r(li-1,ci)和r(li-1,ci-1)分别表示固定阶段i-1的土地利用资料输入不变，分别采用阶段i和阶段i-1的气象资料，经swat模拟得到的径流；

步骤52、保持阶段i气候条件不变，模拟土地利用变化前后径流的变化，以表示阶段i土地利用对径流的影响：

δrl,i＝r(li,ci)-r(li-1,ci)(12)

式中，δrl,i为阶段i土地利用变化引起的径流变化，即间接人类活动对径流的影响；r(li,ci)和r(li-1,ci)分别表示固定阶段i的气象资料输入不变，分别采用阶段i和阶段i-1的土地利用资料，经swat模拟得到的径流；

步骤53、定量解析出气候变化和下垫面变化这种间接的人类活动对径流的影响后，剩下的径流变化可以近似归纳到流域内生活生产、流域外引水调水等直接的人类活动：

δrd,i＝δrt,i-δrc,i-δrl,i＝ro,i-r(li,ci)(13)

式中，δrd,i为阶段i直接人类活动引起的径流变化，δrt,i为阶段i径流总变化，ro,i为阶段i实测径流。

本发明达到的有益效果：首先，划分天然时期和变化时期，根据土地利用的演变特征将变化期细致划分，每段土地利用类型在面积和分布上具有显著变化，以充分考虑土地利用演变情况；其次，基于swat模型能够准确描述流域下垫面特征、有效模拟水文过程的特征，建立起流域分布式水文模型以描述水文过程；最后通过对流域研究时期径流变化逐段、逐项滚动式计算，获得的归因定量识别结果充分考虑了土地利用演变对径流的影响。

附图说明

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

图2为研究时期划分示意图。

具体实施方式

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

本发明提供一种考虑土地利用时空演变的径流变化归因识别方法，包括如下步骤：

步骤1，根据径流突变点将研究时期划分为天然期、变化期；

步骤2，根据土地利用演变特征将变化期分成m段；

步骤3，建立流域swat模型；

步骤4，根据天然期资料率定流域swat模型参数；

步骤5，对流域径流变化逐段逐项滚动式归因定量计算。

在进一步的实施例中，本发明采用pettitt检验方法对流域径流时间序列进行突变检验，确定径流发生显著变化的时间，从而将整个研究时期划分为天然时期和变化时期，并认为天然期受气候变化和人类活动影响较小，能够反映该流域天然径流状况。为了考虑土地利用情况的时空演变，尤其是变化时期，土地利用情况在时间和空间上的变化更加深刻，本发明根据土地利用的演变特征将变化期分成n段。通过arcgis平台下的arcswat模块建立流域swat模型，需要基于流域dem数据资料生成水系并进行子流域的划分，然后根据土地利用类型、土壤类型、坡度类型的空间分布情况生成水文响应单元，最后输入气象数据，运行模型。为了提高参数率定的效率，在进行参数敏感性分析后，排除众多参数中对模型结果影响较小的参数，选取敏感性较强的参数进行率定，具体的参数率定工作，利用swat-cup软件工具完成，得到了参数率定值之后，再对率定的效果进行评价。对流域径流变化逐段逐项滚动式归因定量计算，“逐段滚动”指分别将变化期内的每一段较其前一段进行归因定量计算。“逐项滚动”指针对每一段，首先保持下垫面条件不变，模拟气候条件变化前后的径流变化，以表示气候变化对径流的影响；再保持气候条件不变，模拟土地利用变化前后的径流变化，以表示土地利用对径流的影响；最后剩下的径流变化近似归纳到直接人类活动。

接下来，描述一种实施方案。如图1所示，一种考虑土地利用时空演变的径流变化归因识别方法，包括如下步骤：

步骤1，根据径流突变点将研究时期划分为天然期、变化期：

采用pettitt检验方法找到径流时间序列突变点，判断突变点前后的累积分布函数是否存在显著的差异，构造统计量ut,n：

式中：

记统计量ut,n的最大值为kτ，对应的样本τ就是突变点，显著性水平检验公示如下：

确定径流发生显著变化的时间后，以突变点为界将整个研究时期划分为天然时期和变化时期，并认为天然期受气候变化和人类活动影响较小，能够反映该流域天然径流状况。

步骤2，根据土地利用演变特征将变化期分成n段：

搜集土地利用资料，统计不同时刻的各土地利用类型面积，反映各类型在数量上的变化，制作土地利用转移矩阵，反映不同土地类型在空间分布上的互相转换。根据土地利用的时空演变特征将变化期分成具有显著变化的m段(图2)，分段归因以充分考虑不同阶段土地利用对径流的影响。如图2，此时整个研究时期被划分成了m+1段。

步骤3，通过arcgis平台下的arcswat模块建立流域swat模型，具体分为以下几个子步骤：

(1)子流域划分：根据流域dem图层和人为设定的流域阈值，生成河网水系，以及与河网结构相关的若干出水口和入水口。在根据实际需要人为添加另外的出水口和入水口之后，确定流域总出水口，基于水力联系进行子流域划分，并分析计算各子流域的地形参数、几何参数、水流路径等。

(2)水文响应单元生成：划分子流域后，需要对各子流域划分水文响应单元(hru)，反映不同下垫面条件下水文循环的特征。hru是一个子流域内具有相同的土地利用类型、土壤类型和坡度类型的区域，是模型运行中的最小地块单元。将土地利用、土壤、坡度类型的空间分布数据集载入模型，根据其组合和分布，可生成若干个hru。

(3)数据输入与模型运行：将各类土壤属性数据、气象数据等输入模型，在模型数据库中创建流域输入文件(.bsn)、子流域输入文件(.sub)、hru输入文件(.hru)、管理输入文件(.mgt)、土壤输入文件(.sol)、地下水输入文件(.gw)等输入文件，然后可以进行模型的运行。模型运行产生的结果主要包括流域内各hru、各子流域、各演算河段的汇总信息，分别存放在hru输出文件(output.hru)、子流域输出文件(output.sub)、主河道输出文件(output.rch)中。对模型运行结果进行分析计算，可以得出流域各部分的天然径流量。

步骤4，根据天然期资料率定流域swat模型参数，具体分为以下几个子步骤：

(1)参数敏感性分析：为了提高参数率定的效率，在进行参数敏感性分析后，排除众多参数中对模型结果影响较小的参数，选取敏感性较强的参数进行率定。

敏感性是一个无量纲的指数，反映了模型结果随参数改变而变化的程度。敏感性结果的计算形式如下：

式中：si为参数i的敏感性结果；o为模型输出结果，其数学表达式为o＝f(f1,f2,…,fn)；fi为影响模型结果的参数；n为参数个数。

敏感性结果可通过绝对值、相对值或离差的形式进行表示。为对比不同参数的敏感性，选取相对值的形式，计算各参数的敏感性指数i，公式如下：

(2)参数率定：采用swat-cup软件工具完成，选择sufi-2算法。该算法考虑了模型结构、输入数据、参数、实测数据等因素的不确定性。通过拉丁超立方随机抽样法(latin-hypercubesimulations)随机生成一组参数，代入swat模型中进行模拟，并比较模拟值与实测值的差异，优选参数范围，通过多次迭代运算，不断减小参数的不确定性范围，使模拟值不断与实测值相接近。

(3)参数率定评价：选用相关系数(r²)和nashsuttcliffe模型效率系数(nash-sutcliffeefficiencycoefficient，nse)2个指标来评价swat模型参数的适用性。

相关系数r²用于评价模拟值和实测值之间的吻合程度，计算公式如下：

式中：qoi为第i时段的实测流量；为各时段实测流量的平均值；qsi为第i时段的模拟流量；为各时段模拟流量的平均值；n为实测数据个数，等于参数率定的时段数。

r²的值常介于0与1之间，愈接近1，说明实测值与模拟值线性相关程度愈高。通常认为，r²大于0.8时，参数率定结果较为合理。

nse是一个整体综合指标，可以定量表示整个径流过程模拟值与实测值的拟合程度，计算公式如下：

式中各变量含义同上，nse的取值范围一般在0～1之间。当nse＝1时，模拟效果最好；当nse接近0或为负值时，模拟值可信度低。通常认为，nse大于0.7时，参数率定结果比较合理。

步骤5，对流域径流变化逐段逐项滚动式归因定量计算，如图2，研究时期已经被分为n+1段，设i＝1,2,...,n，从变化时期第1段开始逐段逐项滚动式计算，具体分为以下几个子步骤：

(1)保持阶段i-1的下垫面条件不变，将阶段i的气象资料输入swat模型，则阶段i气候变化对径流的影响可以表示为：

δrc,i＝r(li-1,ci)-r(li-1,ci-1)(21)

式中，δrc,i为阶段i气候变化引起的径流变化；li-1为阶段i-1土地利用资料，ci-1和ci分别为阶段i-1和阶段i的气象资料；r(li-1,ci)和r(li-1,ci-1)分别表示固定阶段i-1的土地利用资料输入不变，分别采用阶段i和阶段i-1的气象资料，经swat模拟得到的径流。

(2)保持阶段i气候条件不变，模拟土地利用变化前后径流的变化，以表示阶段i土地利用对径流的影响：

δrl,i＝r(li,ci)-r(li-1,ci)(22)

式中，δrl,i为阶段i土地利用变化引起的径流变化，即间接人类活动对径流的影响；r(li,ci)和r(li-1,ci)分别表示固定阶段i的气象资料输入不变，分别采用阶段i和阶段i-1的土地利用资料，经swat模拟得到的径流，其他项含义同上。

(3)定量解析出气候变化和下垫面变化这种间接的人类活动对径流的影响后，剩下的径流变化可以近似归纳到流域内生活生产、流域外引水调水等直接的人类活动：

δrd,i＝δrt,i-δrc,i-δrl,i＝ro,i-r(li,ci)(23)

式中，δrd,i为阶段i直接人类活动引起的径流变化，δrt,i为阶段i径流总变化，ro,i为阶段i实测径流，其他项含义同上。

以上详细描述了本发明的实施方式，在本发明的技术构思范围内，可以对本发明的技术方案进行多种跟她变换，这些等同变化均属于本发明的保护范围。