具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例提供一种基于多源互济供水模式的分布式水循环模型的设计方法，主要步骤包括嵌套式HRU划分、HRU属性设计、多源互济供水模块设计和SWAT模型改进四部分内容。多源互济供水模块设计就是设计并添加相关模块，增加SWAT模型的城乡多水源供水和多水源灌溉模拟功能，完善人工侧枝水循环，进而构建自然-人工分布式水循环模型。SWAT模型改进，即是修改SWAT模型相关模块，便于水源类型、水源个数、水源代码、取水优先序信息的读取以及多源互济模块的调用。所涉及的模块主要包括rchuse模块、res模块、watuse模块、irr_rch模块、irr_res模块、irrsub模块、subbasin模块、percmain模块、gwmod模块、surface模块、recday模块、recmon模块、main模块等。

本发明所述方法包括以下步骤：

步骤S1：嵌套式HRU划分；

步骤S2：HRU属性设计；

步骤S3：多源互济供水模块设计；

步骤S4：SWAT模型改进。

在所述的步骤S1中，采取基于“流域-水资源区-行政区-灌区-土地利用-土壤-坡度”多元属性的嵌套式坡面离散方法进行HRU划分。其主要步骤（图2）主要包括：（1）天然子流域划分；（2）天然HRU划分；（3）天然HRU叠加管理单元；（4）天然HRU叠加灌区。

（1）采用ArcGIS中的 Arc Hydro Tools 工具进行天然子流域划分。首先对 DEM底图进行填洼等修正操作，依据设定的生成子流域临界集水面积CSA（Critical SourceArea）阈值要求，确定格网水流流向，识别流域分水线，并依据流域坡度、坡向、坡长等河网特征参数划分天然子流域；

（2）按照子流域内的土地利用类型、土壤类型和坡度类型进行叠加，划分出若干个天然HRU。 HRU是指下垫面特征相对单一和均匀的区域，在这个区域中的下垫面具有相同的水文特性，每个HRU内只有一种植被、一种土壤和一种坡度类型的组合。

（3）采用选取HRU叠加地市/县行政区和水资源区的嵌套方式。采用GIS的“EditorToolbar”功能，将管理单元边界GIS图直接叠加到上一步所划分的天然HRU图上，在天然HRU与管理区划交界地区，依照边界线将天然HRU一分为二，并将划分的HRU附上管理单元属性；

（4）采用GIS的“Editor Toolbar”功能，在天然HRU与灌区交界地区，依照边界线将天然HRU一分为二，并将划分的HRU附上灌区属性，最终完成HRU的划分。最后，每一个HRU都具有子流域属性、水资源区属性、行政区属性、灌区属性、土地利用类型属性和土壤类型属性。

在所述步骤S2中，HRU属性的设计包括如下内容：

（1）HRU属性识别模块readattr构建：

构建HRU属性识别模块readattr，用于读取指定的HRU属性。设置的HRU属性分别是子流域属性、水资源区属性、行政区属性、灌区属性。输入文件格式见表1。

（2）HRU属性识别模块的调用

将构建的readattr模块（HRU属性识别模块）放在SWAT模型中的main主模块中，便于对readattr模块的调用。

在所述步骤S3中，多源互济供水模块Multi_sc的设计包括如下内容：（1）水源代码信息文件设计；（2）供水优先次序信息文件设计；（3）多源互济供水计算设计。

（1）水源代码信息文件设计：

用于读取指定的供水源信息。设置的供水源类型包括外调水、水库水、河道水、浅层地下水、深层地下水、坑塘水6种水源类型，并且设定单个HRU的同类型水源最多为5个，即每个HRU最多可以设置30个供水源。水源代码以输入文件读入的方式进行指定，并针对不同的HRU类型，居工地HRU以行政区为单位进行指定，同一行政区内的HRU供水源一致；灌溉用地HRU以灌区为单位进行指定，同一灌区内的HRU供水源一致。输入文件格式分别见表2和表3。

（2）供水优先次序信息文件设计：

用于读取指定的水源供水优先次序信息。每一个水源可以为多个HRU利用，每一个HRU也可以从多个供水源取水。根据每一个HRU的供水源类型和个数，设定供水次序。指定供水源的供水优先次序，供水优先次序以输入文件读入的方式进行指定，并针对不同的HRU类型，居工地HRU以行政区为单位进行指定，同一行政区内的HRU水源供水顺序一致；灌溉用地HRU以灌区为单位进行指定，同一灌区内的HRU水源供水顺序一致。优先次序以1到k（k为HRU的供水水源个数，k≤30）进行排序，1为最高优先序，对未指定的水源，序号设置为0。输入文件格式见表4和表5。

（3）取水控制量信息文件设计：

用于读取各行政区各年度的供水控制总量信息。每个HRU根据自身的行政区属性，识别所属行政区的年度地表水供水控制总量和地下水开采控制总量，用于后续的水源的取水量计算。输入文件格式见表6。

（4）多源互济供水计算设计：

多源互济供水模块Multi_sc的运行流程（图3）主要包括：1）HRU类型识别；2）供水源及其优先序识别；3）多源互济供水计算。

1）HRU类型识别

首先识别HRU的土地利用类型。如果是建设用地（城镇用地、农村用地），则程序进入城乡供水计算流程；如果是灌溉用地，则程序进入灌溉用水计算流程；如果是其他土地利用类型，则该HRU无供水要求，程序退出。

2）供水源及其优先序识别

多源互济供水模块通过识别HRU的取水源识别代码，调用相应的水源模块，河道、水库、浅层地下水、深层地下水、外调水和坑塘水的取水源识别代码分别是1、2、3、4、5、6。通过水源代码信息文件和供水优先次序文件的读取，识别每一个HRU的供水源个数、类型、取水先后顺序，根据取水源识别代码依次调用各水源模块。

3）多源互济供水计算

多源互济供水模块根据HRU的供水次序寻求各水源进行取水，直至满足该HRU的当日的生活、工业用水需求（建设用地HRU）或者农业灌溉用水需求（农业用地HRU）或者到最后一个供水源供水完毕为止。供水优先计算的步骤（图4）如下：

①明确目标HRU设定的当日需水量WD；

②明确目标HRU的供水水源个数k（k≤30）、水源代码及其供水优先次序；

③按次序调用水源模块，计算水源的取水量。水源模块包括rchuse模块（河道取水模块）、res模块（水库取水模块）、watuse模块（地下水、坑塘水的取水计算）。水源取水量取决于当日HRU的需水量和水源的可供水量。可供水量则取决于水源的可利用水量、取水工程（如引水渠道、供水管道、电井等）供水能力和取水控制量。计算公式如下：

式中，i表示HRU序号；j表示水源的供水优先次序号；WSP表示水源的当日实际取水量，m³；WD表示HRU的当日需水量，m³；Wsc表示水源的当日可供水量，m³；WF表示取水工程的供水能力，m³；WA表示水源的当日可利用水量，m³；WMX表示年度取水控制量，m³。WUM表示年用水控制总量，m³；WSM表示年度地表取水控制总量，m³；WGM表示年度地下水开采控制总量，m³。

④对供水优先次序为1的水源，则优先从该水源取水。如果水源的可供水量Wsc₁＞WD，那么该水源供水量为WSP₁=WD，供水程序结束，HRU的水源总供水量为WSP= WSP₁；否则，WSP₁=Wsc₁，并且HRU的需水量变为Wf =WD-Wsc₁，程序将继续寻找下一级供水源；

⑤对供水优先次序为j（j=2,…,k-1；k≤30）的水源，如该水源的当日可供水量Wsc_j＞Wf，那么该水源供水量为WSP_j=Wf，程序结束，HRU的水源总供水量为WSP=WSP+WSP_j；否则，WSP_j=Wsc_j，并且HRU的需水量变成Wf=Wf-Wsc_j，程序将继续寻找下一级供水源；

⑥对供水优先次序为k（k≤30）的水源，如果该水源的当日可供水量Wsc_k＞Wf，那么该水源供水量为WSP_k=Wf，程序结束，HRU的水源总供水量为WSP=∑WSP _i ；否则，WSP_k=Wsc_k，并且HRU的需水量变成Wf=Wf-Wsc_k，程序结束。

在所述步骤S4中，SWAT模型改进主要是将多源互济供水模块与SWAT模型中修改的模块进行衔接以实现数据实时交互，所述HRU根据输入的需水信息、供水水源类型信息、供用水优先序规则信息、水利工程信息和由SWAT模型提供的天然水文条件信息进行逐日水资源分配与调控，输出逐日“供水-用水-耗水-排水”的人工侧枝水循环信息，将人工侧枝水循环信息传递给SWAT模型。其中修改SWAT模型中的相关模块，便于水源类型、水源个数、水源代码、取水优先序信息的读取以及多源互济模块的调用。修改的模块主要有rchuse模块、res模块、watuse模块、irr_rch模块、irr_res模块、irrsub模块、subbasin模块、surface模块、gwmod模块、recday模块、recmon模块和main模块等。各模块主要功能如下：

rchuse 模块：河道供水模块，从指定的河道中取水，用于生活工业用水。

res模块：水库供水模块，从指定的水库中取水，用于生活工业用水。

watuse模块：地下水、坑塘供水模块，从指定的子流域内的浅层地下含水层、深层地下含水层和坑塘中取水，用于生活工业用水。

irr_rch 模块：河道灌溉模块，从指定的河道中取水，用于农田灌溉用水。

irr_res模块：水库灌溉模块，从指定的水库中取水，用于农田灌溉用水。

irrsub模块：地下水、外调水灌溉模块，从指定的子流域内的浅层地下含水层、深层地下含水层和外源水源中取水，用于农田灌溉。

percmain模块：渗漏损失模块，根据降雨入渗、灌溉入渗、坑塘入渗计算结果，模拟土壤水变化。

gwmod模块：地下水模块，根据地表入渗结果，模拟浅层、深层地下水变化以及地下径流。

subbasin模块：子流域模块，主要包括蒸发模拟、产流模拟、地下水模拟、植被生长模拟、农田管理模拟、营养物质迁移转化模拟等内容，通过对水文、水质过程相关模块的调用实现。

surface模块：地表径流模拟模块，模拟冠层截流、积雪融雪、地表产流等主要地表水文过程。

recday模块：日点源模块，通过输入逐日污染物排放量，模拟逐日点源变化过程。

recmon模块：月点源模拟，通过输入逐月污染物排放量，模拟逐月点源变化过程。

main模块：SWAT模型主模块，主要包括读入输入文件，分配数组大小、初始化参数、水文过程模拟等内容。

修改SWAT模型中的相关模块，具体包括：

（1）模块屏蔽

屏蔽水源模块，即rchuse 模块、res模块、watuse模块、irr_rch 模块、irr_res模块和irrsub模块，不再采取单一水源取水模式。将以上模块放在多源互济供水模块中进行调用。

（2）修改rchuse 模块、res模块

添加代码，将参数WSP _i （i=1,2）分别替换掉rchuse 模块、res模块中的参数waterrch、wuresn，实现多源互济供水模块Multi_sc与rchuse 模块、res模块的衔接和调用；

（3）修改watuse模块

①由于该模块中没有外调水供水功能，通过修改程序增加外调水取水以及外调水量限制约束，供水量不能超过调水限制总量，计算公式表达如下：

其中，waterout (i,j)表示第i天、第j个HRU的调水利用量，m³；WX5表示调水限制总量，m³。

②在watuse模块中添加代码，用参数WSP _i （i=3,4,5,6）分别替换掉watuse模块中的参数watershal、waterdeep、waterout和waterpnd，实现多源互济供水模块Multi_sc与watuse模块的衔接与调用。

（4）添加管网输水模拟功能

城市供水管网在向用水户输送水量的过程中，会因管网材质老化或破损等因素产生“跑、冒、滴、漏”等输水损失。因此有必要考虑城市供水损失。在rchuse 模块、res模块、watuse模块中添加如下公式的计算程序代码：

式中，pip表示供水管网漏损率。

（5）修改irrsub模块

由于模块中没有坑塘灌溉，添加坑塘灌溉功能；完善外调水灌溉功能，通过增加供水量限制约束，灌溉取水量不能超过调水限制总量。

其中，wirrout (i,j)表示第i天、第j个HRU的调水灌溉利用量，m³；WX5表示调水限制总量，m³。wirrpnt (i,j)表示第i天、第j个HRU的坑塘灌溉利用量，m³；WX6表示坑塘可供水量，m³。

（6）添加灌溉渠道模拟功能

①SWAT模型并不考虑灌溉渠道的输水损失对水循环的影响，灌溉渗漏视为系统损失。通过修改irr_rch模块、irr_res模块、irrsub模块的源代码，增加渠系输水过程模拟，包括渠道水量损失和渠道退水。其中，渠道水量损失包括渠道内的水分蒸发损失和渠道渗漏损失两部分。主要计算公式如下：

式中，ET _can 表示渠系蒸发损失量，mm；IRR _can 表示进入渠道的灌溉水量，mm；Ls _can 表示渠系渗漏损失量，mm；Surp _can 表示渠系退水量，mm；表示渠系水有效利用系数；表示渠系蒸发系数；表示渠系渗漏系数。

②渗漏损失水量进入上层土壤补给土壤水，参与土壤水分循环。通过修改percmain模块相关代码，添加渗漏损失计算程序。计算公式如下：

式中，表示第t天第一层土壤的土壤含水量，mm；表示第t-1天第一层土壤的土壤含水量，mm；inf _pcp 表示降雨入渗量，mm；inf _irr 表示灌溉入渗量，mm；inf _wet 表示湖库湿地渗漏量，mm；

（7）修改gwmod模块代码

供水管网漏失水量进入地下含水层，补给潜水。通过修改gwmod模块中的地下水补给代码，实现对管网漏损水量的模拟。计算公式如下：

式中，表示第t天的潜水补给量，mm；表示第t-1天的潜水补给量，mm；prc表示补给地下水的土壤渗漏水量，mm；GW_DELAY表示地下水补给延迟系数，mm；Area表示HRU的面积，m²；

（8）subbasin模块修改

在subbasin模块中添加多源互济供水模块，便于SWAT模型运行时的调用。通过读取指定的水源类型、水源个数、水源代码、取水量、取水时间，进行年内的水源动态互补供水运算，实现SWAT模型运行时的多水源联合供水模拟。

（9）surface模块修改

渠系退水是指渠道排泄超引灌溉水，渠道退水直接进入河道，参与河道汇流计算。通过修改surface模块相关代码，将叠加地表产流量。计算公式如下：

式中，表示渠道退水后的径流量，mm；表示渠道退水前的径流量，mm；

（10）点源模块修改

城市生活和工业污水产生后，通过排水管网系统直接排入河流，或者输送至污水处理厂进行处理。处理达标后的污水一部分直接排放进入河道，一部分通过深度处理达到再生水利用标准后回用于绿化、生活杂用以及生产等用途。SWAT模型采用其点源模块（recday模块和recmon模块）进行模拟。在recday模块和recmon模块中修改相关代码，将污染排放参数WDR分别替换参数floday、flomon。相关计算公式如下：

其中，WDR为城市污水产生量，m³；WP为污水排放量，m³；r表示用水消耗率；v表示污水处理厂的污水处理率；re表示再生水利用率。

（11）main主模块修改

将构建的readattr模块放在SWAT模型中的main主模块中，便于对readattr模块的调用。

用Fortran语言在Windows平台应用程序开发环境Visual Studio 2012中根据以上结构对SWAT源代码进行修改，实现了对多源互济供水模块的调用，具有了城市多水源供水和多水源农业灌溉模拟功能，从而构建了以SWAT模型为基础的、具有多源互济功能的分布式自然-人工水循环双向互馈模拟模型。

本发明选取白河流域作为实施案例，白河流域自然-人工水循环模拟过程介绍如下：

1、研究区概况

白河流域地处汉江流域中部，北起伏牛山地区，西与老鹳河为邻，东靠唐河流域，南抵汉江，大部分地区位于南阳盆地内，流域总面积为1.23万km²。地形呈北高南低态势，水系主要由白河及其支流湍河、刁河、严陵河等，多年平均水资源量20.8亿m³。流域内有河南省南阳市卧龙区、宛城区、新野县、湖北省襄阳市襄州区等3个地级市、共12个县级行政单元。白河流域境内有鸭河口灌区、引丹灌区等大型灌区4个，赵湾水库灌区、高丘水库灌区、陡坡水库灌区等中型灌区18个。为防洪抗旱、保障农业灌溉用水，流域内建有大型水库（鸭河口水库）1座，赵湾水库、虎山水库、陡坡水库、官沟水库等中型水库十几座。

2、基础数据收集

模型构建所需要的数据有空间数据，包括DEM（90m×90m），土地利用图（1:10万）、土壤分布图（1:100万）、行政区分布图、灌区分布图、以及水系图等；白河流域内及邻近的4个气象站的1990~2016年气象数据，包括日降水量、日最高及最低气温、相对湿度、日照时数风速等气象要素；水利工程数据，主要包括水库位置、死库容、调节库容、总库容等数据，灌溉渠道过水能力、电机井日提水能力等数据；新店铺水文站1991~2016年的逐月径流数据，用于模型的率定和验证；农业灌溉用水数据，通过查询2006~2016年南阳市水资源公报、襄阳市水资源公报等资料，获得历年农业灌溉用水数据；灌溉面积与种植结构数据，包括小麦、水稻、花生、油菜、芝麻、棉花、蔬菜、瓜果等各类作物种植面积。此外，还有播种、灌溉、施肥、收割等作物管理措施有关的时间、次数、单次（灌溉、施肥）数量等基本信息等。

3、建模过程

针对白河流域水资源分布、水利工程特点及各行业取用水情况，模型构建步骤如下：

（1）HRU划分：基于DEM数据、土地利用类型图、土壤类型图划分天然子流域和HRU，共获得34个天然子流域和507个天然HRU；在此基础上，根据行政区划及灌区分布情况，将天然HRU细分成1027个HRU，每个HRU均具有子流域属性、水资源区属性、行政区属性、灌区属性、土地利用类型属性和土壤类型属性。

（2）农业种植与管理信息输入：包括各行政区县的种植作物类型、作物种植面积与灌溉面积、轮作制度与灌溉制度等。

（3）供水规则信息输入：主要包括水源供水对象、分区供水原则、水源供水优先级、用水户分水原则、供水水源类型与个数、行业用水优先级等信息；此外，需要输入的信息还有水库、渠道、电机井等水利工程信息，渠系水有效利用系数、田间水有效利用系数等、耗水率等信息。

（4）气象数据信息输入：选择建立好的降水、气温、风速、辐射和相对湿度数据输入，再把所有属性数据及水库数据输入，模型构建完成后即可开始运行。

4、参数率定与模型校验

通过上述过程和分析，确定模型的主要参数值，参数调整的最终值见表7所示。通过对模型参数的调整，模型的径流模拟与实测径流过程对比结果见表8和图5所示。在率定期（1995~2005年），新店铺站月径流模拟值与实测值的相关系数为0.792，纳什效率系数为0.756；在验证期（2006~2016年），新店铺站月径流模拟值与实测值的相关系数为0.643，纳什效率系数为0.635。可以看出，月径流模拟值与实测值流量过程线拟合程度较好，模型的模拟精度达到了要求值。

对全部位于流域内的行政区（卧龙区、宛城区、镇平县）的模拟结果的偏差百分比见表9所示。可以看出，2010~2016年三个县（区）的供用水量偏差很小，其生活用水量、工业用水量、农业用水量、地表水供水量以及地下水供水量与实际水量的误差均在10%以内。可见，模型的模拟结果很好地反映了白河流域实际用水和供水情况。

5、模拟结果分析

通过对白河流域2016年水循环模拟汇总，获得流域水资源循环转化关系见图6所示。2016年全流域总降水量95.47亿m³，地表产流量22.73亿m³，年末土壤水蓄变量8.62亿m³、地表水蓄变量1.25亿m³、地下水蓄变量为0.65亿m³。全流域耗水总量为70.61亿m³，其中土壤蒸发29.90亿m³，植被散发33.09亿m³，截留蒸发6.25亿m³，积雪升华0.003亿m³，水面蒸发1.17亿m³，生活工业消耗0.45亿m³。考虑外调水利用，全流域的经济社会总用水量为9.53亿m³，其中，地表水利用量6.21亿m³（含外调水4.47亿m³），地下水利用量3.32亿m³；人工总耗水量为5.34亿m³，人工排水量为1.71亿m³，总出境水量20.30亿m³。

汇总2016年白河流域各水源供水情况如图7所示。2016年，白河流域水资源利用以地下水和外调水（丹江口水库供水）为主，全年地下水开发利用量占流域用水总量的34.83%（3.32亿m³），外调水利用量占流域用水总量的46.92%（4.47亿m³），河道水利用量最小，仅有0.18亿m³，占流域用水总量的1.95%。全年的用水峰值集中在3月和8月，3月份是冬小麦用水的关键的时期，外调水利用量为0.86亿m³，其次为蓄水供水（0.42亿m³），由于地表水资源量相对充足，地下水仅用于补充供水，补充供水量为0.26亿m³。8月份为玉米生长用水的关键的时期，外调水利用量为0.77亿m³，由于蓄水资源量不足，供水量仅为0.2亿m³，地下水起到了重要的补充供水作用，补充供水量为0.75亿m³，在数量上几乎与外调水供水量相等。此外，在10-11月、1月出现用水小高峰，主要是由于种植冬小麦后的灌溉用水。由于这段时间处于枯水季节，地表水资源量有限，供水主要由丹江口水库供水（外调水）和地下水补充灌溉。

表1

表2

表2中若无供水，代码以0表示。

表3

表3中若无供水，代码以0表示。

表4

表4中若无供水，代码以0表示。

表5

表5中若无供水，代码以0表示。

表6

表7

表7中，v表示调参值替代原参数值；r表示原参数值乘以调参值。

表8

表9

本发明以SWAT模型为基础，开发一种分布式自然-人工水循环模型，使其具有了自然-人工水循环动态互馈及水资源开发利用与调控一体化模拟功能，可完成基于城市多水源供水和多水源灌溉供水模式的流域自然-人工水循环过程模拟，在模型运行中，能够时刻保持自然水循环与人工侧枝水循环之间的动态互馈关系，不仅反映水循环过程对人工取用水的影响，同时还反映水资源开发利用与调控对水循环过程的实时干预作用，从而为深入了解强人类活动影响下的流域水循环机制，合理开发和利用水资源提供科学的参考依据。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何属于本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。