基于SWAT模型的流域天然径流计算方法与流程

本发明属于天然径流计算领域，尤其涉及一种基于swat模型的流域天然径流计算技术。

背景技术：

近年来，随着经济社会的发展以及城市化进程的加快，许多流域的下垫面发生了变化。工农业取水、生活取水、地下水开采等人类活动使得实际测得的径流已然无法代表无人类活动影响下的天然径流。计算流域天然径流是流域水资源评价、水资源配置的基础，只有了解流域真实的天然产流，才能谈及合理利用。

计算天然径流的方法有很多，利用分布式水文模型模拟天然径流是方法之一。swat模型是考虑了流域整个产汇流物理过程的分布式水文模型，通过输入流域的气候数据、土壤分布数据和土地利用数据，模拟整个流域的水流运动。swat模型可以很好的应用于流域天然径流的模拟，有效解决上述下垫面变化以及人类取用水活动影响监测径流的问题。

目前，利用swat模型模拟天然径流存在以下不足：(1)模型自带参数值确定功能效率低。swat模型自带的参数值确定及敏感性分析功能效率较低、运算时间长。(2)参数值的不确定性大。由于天然径流难以测得，相应的历史水文数据缺乏，导致模型模拟的参数值存在很大的不确定性，进而导致模拟精度低下，模拟的径流依然无法代表天然径流。

技术实现要素：

针对现有技术的不足，本发明通过搜集dem、土地利用和土壤等资料，建立基础数据库，在此基础上，对目标流域建立了流域swat模型。对流域中的天然子流域确定了相关模型参数值并验证其合理性，然后建立了天然子流域参数与其他子流域参数的相关关系，最后基于参数值确定后的模型计算流域的天然径流，提出了一种基于swat模型的流域天然径流计算方法。

为解决上述技术问题，本发明采用如下技术方案：一种基于swat模型的流域天然径流计算方法，包括以下步骤：

步骤1，建立流域swat模型，构建基础数据库；

步骤2，基于sufi-2方法对天然子流域确定参数值并验证其正确性；

步骤3，建立其余各子流域参数与天然子流域参数相关关系；

步骤4，输入参数，通过swat模型计算流域天然水资源量。

根据本发明的一个方面：所述步骤2的天然子流域为受人类活动影响小于预期值、具备历史实测径流资料的子流域，参数确定及验证的范围仅为天然子流域。所述步骤3的相关关系是将其余个子流域的下垫面情况以及是否与天然子流域接壤指标化，通过指标计算其余各子流域参数值。步骤4中的所述输入参数是代表天然状况下模型的参数。

进一步地，提供一种基于swat模型的流域天然径流计算方法，包括以下步骤：

步骤1，建立流域swat模型和基础数据库；所述基础数据库中的数据包括空间数据资料、土地利用资料、土壤空间分布资料、气象资料和水文资料；上述数据经数据清理、格式统一化处理和投影变换后添加到模型中；其中，空间数据资料、土地利用资料和土壤空间分布资料需要将其图层本身的地理坐标系统转换为投影坐标系统；

基础数据库建立完成后，利用dem图提取流域水系；根据设定的子流域面积范围与子流域进出口的位置等进行子流域的划分；根据不同的下垫面情况，进一步划分水文响应单元hrus；

完成水文响应单元hrus的划分后，将气象数据加载进模型，利用模型进行水文模拟；在模型的运行中，需要对径流模拟、气象模拟、潜在蒸发量模拟和河道演算的方法作出选择；利用以日降水观察为基础的scs径流曲线法模拟径流；利用hargreaves法模拟潜在蒸散发；利用变动存储系数法或马斯京根法进行河道演算；

步骤2，基于sufi-2方法对天然子流域确定参数值并验证其正确性；

模型参数确定是指调整模型参数值、初始边界条件以及限制条件的过程，目的是使模型的模拟计算结果与历史实测数据相吻合；利用swat模型模拟天然径流，需要历史实测天然径流数据；在确定模型参数前，先选取流域内下垫面改变小于预期值、受人类活动影响小于预定值、能反映流域天然径流情况以及有部分历史径流资料的子流域作为天然子流域，基于sufi-2方法确定天然子流域的参数并验证其合理性；

在sufi-2方法中，首先给参数一个相对较大的范围，使实测数据尽量多的被包含在95％不确定性范围内，然后逐步地缩小不确定性的区间范围，同时改变参数值的范围；每一次参数范围的改变，都将重新进行敏感性矩阵和协方差矩阵的计算；更新参数，再进行下一步模拟，使模拟值逐渐接近实测值；其中，敏感性矩阵和协方差矩阵的计算方程如下：

j＝1,2,3,...,m(3)

式中，j为敏感性矩阵；ji,j为敏感性矩阵元素；i为组数；j为列数；δbj为第j个要率定的参数；δgi为参数的敏感度；为敏感性矩阵的行数；c为协方差矩阵；为目标函数计算结果的偏差，m为自然数；

在swat模型中关于径流的参数有26个，先利用sufi-2算法率定天然子流域的26个参数；率定时选定预热期、率定期和验证期，利用sufi-2算法对典型子流域进行参数值的确定，调整参数，使得模拟值与实测值吻合，最终确定参数最佳适应值；

确定天然子流域参数值后利用历史实测径流数据进行验证，为了评价参数值是否可靠，选取纳什效率系数ens和相关系数r²作为评价指标；ens取值范围是0～1，取值越接近1代表整个模拟径流过程和实测径流数据的拟合程度越好；相关系数r²的取值范围是0～1，取值越接近1，模拟值与实测值的线性相关程度越高；ens和r²的计算公式如下：

其中：qp为实测值；q0为模拟值；qavg为实测平均值；n为实测数据个数；

其中：q0i为实测流量；为实测流量平均值；qsi为模拟流量；为模拟流量平均值；

当ens<0.6，r²<0.6时，认为模拟值接近实测值，参数率定合理；若ens和r²的值不在合理范围内，利用sufi-2方法重新确定参数值；

步骤3，建立其余各子流域参数与天然子流域参数相关关系

受人类活动影响大于预期值或缺少历史实测径流资料的子流域模拟的径流无法代表天然径流，因此本发明建立天然子流域参数与其他子流域参数的相关关系，将天然子流域参数移植或倍比移植到其余的子流域；

步骤3具体包括如下子步骤：

步骤31，将其他子流域分为两部分，与天然子流域接壤的其他子流域为一部分，与各天然子流域都不接壤的其他子流域为另一部分；

步骤32，与天然子流域接壤的其他子流域直接沿用与之接壤的天然子流域的参数，如果一个子流域与多个天然子流域接壤，计算此子流域与各天然子流域接壤长度占比，将各天然子流域参数值按占比分配给此子流域，具体公式如下：

其中：ea,j是a号子流域第j个参数的参数值；n是与a号子流域接壤的天然子流域数；xi是第i个天然子流域与a号子流域的接壤长度；ei,j是与a号子流域接壤的第i个天然子流域第j个参数的参数值；

步骤33，与各天然子流域都不接壤的子流域先将其土地利用和土壤分布情况指标化；具体做法是：计算子流域不同土地利用类型面积占比和不同土壤分布面积占比，然后将土地利用数据和土壤分布数据做归一化处理，归一化公式如下：

其中：ti,j是第i个子流域第j个土地利用类型/土壤类型的归一化数据；gi,j是第i个子流域第j个土地利用类型/土壤类型的面积占比；gi,min是第i个子流域中土地利用类型/土壤类型面积占比的最小值；gi,max是第i个子流域中土地利用类型/土壤类型面积占比的最大值；

计算与各天然子流域都不接壤的其他子流域的土地利用和土壤分布归一化数据与各天然子流域土地利用和土壤分布归一化数据的欧式距离；计算子其他子流域与各天然子流域欧氏距离占比，将各天然子流域参数值按占比分配给其他子流域，公式如下：

其中：ea,j是a号子流域第j个参数的参数值；n是天然子流域数；oi是第i个天然子流域与a号子流域的欧式距离；ei,j是第i个天然子流域第j个参数的参数值；

步骤4，输入参数，通过swat模型计算流域天然水资源量

天然子流域和其余子流域的参数都确定后，将参数输入swat模型，以月尺度模拟目标流域在选定时间段内的径流，并根据流域边界上事先添加的控制点，用流出目标流域的水量减去流入该流域的水量，即得目标流域天然径流系列。

与现有技术相比，本发明具有以下优点和有益效果：

1.本发明充分利用了能反映流域天然径流情况的子流域的历史径流资料，建立天然子流域参数与受人类活动影响较剧烈的子流域参数的相关关系，进而将天然子流域确定性较高的参数移植到其他子流域，提高整个目标流域天然径流模拟值的可靠性。

2.本发明提出的模型参数确定方法效率高，与物理实际相结合，为天然径流的计算提供了新方法、新思路。

附图说明

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

图2为流域swat模型建立流程图。

具体实施方式

具体地，本发明为了解决现有技术存在的问题，提供了一种基于swat模型的流域天然径流计算方法，该方法主要包括以下步骤：

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

步骤2，基于sufi-2方法对天然子流域确定参数值并验证其正确性；

步骤3，建立各子流域参数与天然子流域参数相关关系；

步骤4，输入参数，通过swat模型计算流域天然水资源量。

步骤1中基础数据库所需要的资料有：研究区域的空间数据资料、土地利用资料、土壤空间分布资料、气象资料和水文资料等。在建模的过程中，根据流域实际情况设定子流域面积范围以及流域径流进出口的位置，基于此进行子流域的划分。

步骤2中选取流域内受人类活动影响较小、能大致反映流域天然径流情况以及有部分历史径流资料的子流域作为天然子流域。

步骤3中相关关系的建立以是否与天然子流域相接壤为评判依据，将子流域分为两部分。每部分的子流域以不同的指标进行参数的移植。与天然子流域接壤的子流域沿用天然子流域的参数或者根据接壤长度占比分配。与天然子流域不接壤的子流域根据各自的下垫面条件与各天然子流域下垫面条件的相似性进行参数的移植。

即，在本案中，通过输入多种数据资料，建立流域swat模型；选定天然子流域，确定天然子流域参数值并验证其合理性；在受人类活动影响较大和历史水文资料缺失的其他子流域，根据与天然子流域的相关关系进行参数的移植；最后通过参数确定后的swat模型计算流域天然径流。

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

一种基于swat模型的流域天然径流计算方法，包括以下步骤：

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

如图2所示，swat模型需要多种数据资料的支持，包括：空间数据资料、土地利用资料、土壤空间分布资料、气象资料和水文资料等资料。这些资料来源不同，形式各异，需要经过统一的处理和投影变换后再添加到模型中。其中，空间数据资料、土地利用资料、土壤空间分布资料需要将其图层本身的地理坐标系统转换为投影坐标系统。

基础数据库建立完成后，先利用dem图提取流域水系；再根据设定的子流域面积范围与子流域进出口的位置等进行子流域的划分；根据不同的下垫面情况，进一步划分水文响应单元(hrus)。

完成hrus的划分后，将气象数据加载进模型，然后利用模型进行水文模拟。在模型的运行中，需要对径流模拟、气象模拟、潜在蒸发量模拟、河道演算等的方法作出选择。在swat模型中，有三种模拟径流的方法可供选择：日降水数据/scs径流曲线法，以日为时间单位模拟径流；小时降水数据/green&ampt法，以日为时间单位模拟径流；小时降水数据/green&ampt法，以小时为时间单位模拟径流。根据不同的资料情况和模拟时间尺度，可以选择不同的方法。在本例中选用以日降水观察为基础的scs径流曲线法。潜在蒸散发有三种计算方法：hargreaves、priestley-taylor和penman-monteith。本例选用hargreaves法。在swat模型中，可利用变动存储系数法(variablesrorage)和马斯京根法(muskingum)这两种方法进行河道演算，且二者模拟的精度相差不大。本例选用马斯京根法。

步骤2，基于sufi-2方法对天然子流域确定参数值并验证其正确性；

模型参数的确定是指调整模型参数值、初始边界条件以及限制条件的过程，目的是使模型的模拟计算结果与历史实测数据相吻合。所以，利用swat模型模拟天然径流，需要历史实测天然径流数据。因此，在确定模型参数前，先选取流域内下垫面改变较少、受人类活动影响较小、能大致反映流域天然径流情况以及有部分历史径流资料的子流域作为天然子流域，基于sufi-2方法确定天然子流域的参数并验证其合理性。

sufi-2算法是一种把尽量多的实测数据包含在尽量窄的不确定性带中的算法。算法首先给参数一个比较大的范围，使实测数据尽量多的被包含在95％不确定性范围内，然后逐步地缩小不确定性的区间范围，同时改变参数值的范围。每一次参数范围的改变，都将重新进行敏感性矩阵和协方差矩阵的计算，然后更新参数，再进行下一步模拟，使模拟值更加接近实测值。敏感性矩阵和协方差矩阵的计算方程如下：

j＝1,2,3,...,m(3)

其中：j为敏感性矩阵；ji,j为敏感性矩阵元素；i为组数；j为列数；δbj为第j个要率定的参数；δgi为参数的敏感度；为敏感性矩阵的行数；c为协方差矩阵；为目标函数计算结果的偏差。

swat模型中关于径流的参数有26个，本例先利用sufi-2算法率定天然子流域的26个参数。率定时选定预热期、率定期和验证期，利用sufi-2算法对典型子流域进行参数值的确定，调整参数，使得模拟值与实测值吻合。最终确定参数最佳适应值。

确定天然子流域参数值后利用历史实测径流数据进行验证，为了评价参数值是否可靠，选取纳什效率系数ens和相关系数r²作为评价指标。ens取值范围是0～1，取值越接近1代表整个模拟径流过程和实测径流数据的拟合程度越好。相关系数r²的取值范围是0～1，取值越接近1，模拟值与实测值的线性相关程度越高。ens和r²的计算公式如下：

其中：qp为实测值；q0为模拟值；qavg为实测平均值；n为实测数据个数。

其中：q0i为实测流量；为实测流量平均值；qsi为模拟流量；为模拟流量平均值。

当ens<0.6，r²<0.6时，认为模拟值接近实测值，参数率定合理。若ens和r²的值不在合理范围内，利用sufi-2方法重新确定参数值。

步骤3，建立各子流域参数与天然子流域参数相关关系；

受人类活动影响大或缺少历史实测径流资料的子流域模拟的径流无法代表天然径流。因此本发明建立天然子流域参数与其余各子流域参数的相关关系，将天然子流域参数移植或倍比移植到其余的子流域。步骤3包含的具体子步骤如下：

步骤(1)，将子流域分为两部分，与天然子流域接壤的子流域为一部分，与各天然子流域都不接壤的子流域为另一部分。

步骤(2)，与天然子流域接壤的子流域直接沿用与之接壤的天然子流域的参数。如果一个子流域与多个天然子流域接壤，计算此子流域与各天然子流域接壤长度占比，将各天然子流域参数值按占比分配给此子流域。具体公式如下：

其中：ea,j是a号子流域第j个参数的参数值；n是与a号子流域接壤的天然子流域数；xi是第i个天然子流域与a号子流域的接壤长度；ei,j是与a号子流域接壤的第i个天然子流域第j个参数的参数值。

步骤(3)，与各天然子流域都不接壤的子流域先将其土地利用和土壤分布情况指标化。具体做法是：计算子流域不同土地利用类型面积占比和不同土壤分布面积占比，然后将土地利用数据和土壤分布数据做归一化处理，归一化公式如下：

其中：ti,j是第i个子流域第j个土地利用类型/土壤类型的归一化数据；gi,j是第i个子流域第j个土地利用类型/土壤类型的面积占比；gi,min是第i个子流域中土地利用类型/土壤类型面积占比的最小值；gi,max是第i个子流域中土地利用类型/土壤类型面积占比的最大值。

然后计算与各天然子流域都不接壤的子流域的土地利用和土壤分布归一化数据与各天然子流域土地利用和土壤分布归一化数据的欧式距离。计算子流域与各天然子流域欧氏距离占比，将各天然子流域参数值按占比分配给子流域。公式如下：

其中：ea,j是a号子流域第j个参数的参数值；n是天然子流域数；oi是第i个天然子流域与a号子流域的欧式距离；ei,j是第i个天然子流域第j个参数的参数值。

步骤4，输入参数，通过swat模型计算流域天然水资源量。

各个子流域的参数都确定后，将参数输入swat模型，以月尺度模拟目标流域在选定时间段内的径流，并根据流域边界上事先添加的控制点，用流出目标流域的水量减去流入该流域的水量，即得目标流域天然径流系列。