变化环境下流域水资源评价方法及二元水循环过程模拟

（一）基础数据准备和处理

模型的输入数据包括水文气象、地表高程信息、河网信息、土地利用／覆盖、土壤信息、水文地质、水利水保工程、社会经济及供用水信息，见表8-1。

表8-1　模型输入的数据

续表

1.水文气象

水文气象资料主要包括降水、日照、气温、相对湿度和风速等气象信息以及实测和还原径流资料。

（1）降水。采集源于水文和气象两个部门的水文气象站点长系列过程数据，具体信息特征如下。

1）水文部门雨量信息参数：选用伊洛河流域1956～2000 年45 年系列雨量站点逐日降水信息，日内过程选用雨量站点降水要素摘录信息。

2）气象部门雨量信息参数：选用伊洛河流域1956～2000 年45 年系列气象站点逐日降水信息。

（2）其他气象信息。收集整理了1956～2000 年逐日气象要素信息，统计项目包括日照、气温、相对湿度和风速。

（3）径流资料。来源于黄河水利委员会水文局，包括实测和还原径流。采集并整理了流域黑石关、白马寺和龙门镇3 个水文站1956～2000年45 年系列逐日的实测流量和逐月的还原流量信息。

2.地表高程信息

本次研究采用的伊洛河流域DEM 来自于美国地质调查局（USGS）EROS数据中心建立的全球陆地DEM（也称GTOPO30）。GTOPO30可直接从互联网上下载，网址是http：／／edcdaac．usgs．gov／gtopo30／gtopo30．asp。GTOPO30为栅格型DEM，它涵括了全球陆地的高程数据，采用WGS84基准面，水平坐标为经纬度坐标，水平分辨率为30 弧秒，整个GTOPO30数据的栅格矩阵为21600行、43200列。详细高程图见图8-8。

3.河网信息

河网信息主要包括从DEM提取的实测河网、模拟河网以及河道断面的统计信息。

（1）实测河网。取自于全国1∶25 万地形数据库。

（2）模拟河网。利用GIS软件，参照实测的水系图，从全流域栅格型DEM上提取出来，使模拟河网与实测水系比较一致。

（3）河道断面。模型计算中需要河道断面形状参数，涉及整个流域子流域的河道。实际上不可能获得如此多的实测资料，所以在分析了大量河道断面实测数据的基础上，用统计等方法来推算河道断面形状参数。

图8-8　伊洛河流域地表高程图

首先是推算实测断面形状参数。实测资料取自于水文年鉴的实测大断面成果表，为获得最大的过水断面，尽可能取不同年份中实测序列最长的实测大断面成果。将最大过水断面概化为梯形，如图8-9 所示。依据实测资料推算梯形的上底宽、下底宽和高度。具体计算如下所述。

图8-9　河道断面示意图

（1）河道过水断面的上底宽。依据资料数据，以测量过程中的左右岸的起点距之差为过水断面的上底面宽。

（2）河道过水断面的高。下底高程的计算分两种方法，对于规则断面，根据实测数据确定底面高程；对于不规则断面，首先利用线性回归的方式确定过水断面的两个侧面位置，取剔除特异点后的相关系数较大的实测点为侧面，并认为其他实测点是在下底面上，以这些实测点高程的平均值作为过水断面的下底高程。过水断面上下底面的高程差即为河道过水断面的高。

（3）河道过水断面的下底宽。利用侧面的回归方程和下底面的高程，插值确定下底面的左右岸位置，从而确定下底面的宽度。

（4）推算225 段河道的断面形状参数。

4.土地利用／覆盖

本次采集到的地表环境信息主要包括以下信息。

（1）土地利用信息。包括经国家相关部门审查批准生产的1986 年、1996 年和2000年三个时段的1∶10万土地利用图。土地利用源信息为各时段的TM 数字影像，波段为4、3、2；地表空间分辨率为30m。土地利用类型的分类系统采用国家土地遥感详查的两级分类系统，累计划分为6 个一级类型和31 个二级类型。通过地表抽样调查，遥感解译精度为93．7%。

（2）植被指数。包括1980～2000年21 年逐旬NOAA／AVHRR 影像，地表分辨率为8km。在该源信息的基础上，依次提取出植被指数（NDVI）、植被盖度（VEG）和叶面积指数（LAI）等有关植被时数信息。

5.土壤信息

土壤及其特征信息采用全国第二次土壤普查资料，其中土壤分布图为比例尺分别为1∶100万和1∶10万两套，土层厚度和土壤质地均采用《中国土种志》上的“统计剖面”资料。为进行分布式水文模拟，根据土层厚度对机械组成进行加权平均，采用国际土壤分类标准进行重新分类。

6.水文地质

（1）主要水文地质参数。伊洛河流域水文地质参数分布（μ值、K 值）均采用全国水资源规划地下水评价中的相关资料。

（2）岩性分区。采用《中国水文地质分布图》的分区资料。

（3）含水层厚度。采用《中国水文地质分布图》的分区资料。

7.水利水保工程

主要的水利工程包括水库和灌区。

（1）水库。本研究重点考虑了截止到2000 年流域内已起用的大型水库，主要包括水库的空间定位与属性数据两个方面。

水库的空间定位是指确定水库坝址处的空间位置，定位后才能进一步确定水库控制的汇流范围。空间定位依据的资料主要有全国1∶25 万地形数据库、《黄河流域地图集》、全国1∶10万土地利用数据及搜集的各种文字资料。以地形数据库为基础得到大多数水库初步的空间位置，再利用其他资料对初步结果进行补充和修正。

水库的属性数据包含的内容较多，主要有水库起用日期、水位—库容—面积曲线、特征库容、特征水位、淤积状况、时间系列蓄变量、供水目标等。表8-2为流域内主要水库。

表8-2　伊洛河流域大型水库列表

（2）灌区分布。为了研究农业灌溉用水情况，进行了灌区数字化工作。主要是确定了灌区的空间分布范围（见图8-10），收集并整理了灌区的各类属性数据。在灌区数字化过程中，主要参考了国家基础地理信息中心开发的“全国1∶25 万地形数据库”（包括其中的水系、渠道、水库、各级行政边界、居民点分布等）、中科院地理所开发的1∶10 万土地利用图以及黄河水利委员会勘测规划设计院编写的《黄河灌区资料简编》和黄河水利委员会编制的《黄河流域地图集》等资料。

8.社会经济及供用水信息

（1）社会经济信息。主要来源于全国水资源规划、水资源开发利用调查评价部分的成果，以水资源三级区和地级行政区为统计单元，收集整理了1980 年、1985 年、1990 年、1995 年、2000年等5 年与用水关联的主要经济社会指标。2000 年按水资源四级区和地级行政区填报，其余4个年份按水资源二级区和省级行政区填报。

图8-10　伊洛河流域灌区分布图

（2）供、用、耗水信息。主要来源于全国水资源规划、水资源开发利用调查评价部分的成果，以水资源三级区和地级行政区为统计单元，收集整理了1980 年、1985 年、1990年、1995 年、2000年等5 个典型年份不同用水门类的地表水、地下水供、用、耗水信息。2000年按水资源四级分区和县级行政区填报，其余4 个年份按水资源三级区和地级行政区口径填报。

（3）灌溉制度。伊洛河P＝75%的灌溉制度。

（4）种植结构。现状年流域各种作物播种面积。

（二）参数初值确定

WEP L模型的参数可分为三类：第一类是地表面及河道系统参数，包括坡面及河道的曼宁糙率、河床材料的覆盖厚度及导水系数、城市土地利用的不透水率及地表洼地最大储留深；第二类是植被参数，包括植被覆盖度、叶面积指数（LAI）、植被高度、根系深度、最小叶孔阻抗以及地表覆盖的空气动力学参数等；第三类是土壤与含水层参数，包括土壤层厚度、土壤空隙度、土壤入渗湿润峰吸力、饱和土壤导水系数、土壤水分—吸力特征曲线参数、土壤水分—导水系数关系参数、含水层厚度、含水层导水系数等。以上参数均具有物理意义，可根据观测实验数据和遥感数据估算，但由于参数在每个计算单元内仍具有空间变异性，模拟计算时往往使用其计算单元内平均参数或称为有效参数，因此，一些关键参数需要结合模型检验，根据流量过程线和地下水水位的观测结果进行调整。

1.地表面及河道系统参数

地表面以及河道系统参数影响流域汇流，包括坡面汇流与河道汇流。河道汇流中考虑了河水与地下水的水量交换。曼宁糙率反映边界表面的粗糙程度对水流阻力作用影响的参数，表面越粗糙，糙率值越大。坡面汇流的曼宁糙率取计算单元内各类土地利用的曼宁糙率面积调和平均值，河道曼宁糙率根据各断面实测洪水反演值（水文年鉴）并参照王国安、李文家《水文设计成果合理性评价》设定。

河水与地下水的水量交换与河床材料特性密切相关，模型将河床材料导水系数除以河床材料厚度的商作为调试参数，首先根据各地河床材料和厚度选取初值，再作模型率定。

城市土地利用的不透水率和地表洼地最大储留深影响产流过程，前者根据遥感影像进行确定，后者则根据各类土地利用进行初值设定，然后进行模型率定。

2.植被等地表覆盖物参数

植被覆盖度是区域内植被垂直投影面积占地表面积的百分比，是影响地表植被蒸腾、土壤水分蒸发损失以及能量交换的一个重要的控制因子，因而是在流域水循环模型模拟中需要考虑的一个重要参数。

叶面积指数（LAI）是陆面水文过程中十分重要的结构参数，是表征植被冠层结构最基本的参量之一，它控制着植被的许多生物、物理过程，如光合、呼吸、蒸腾、碳循环和降水截获等，为模型定量化描述植被冠层截留量提供了有力的支持。

植被高度和根系深度作为植被的主要参数，影响着水分循环。植被根系深度控制土壤蓄水量的深度，该深度内的水分是参与水分循环的活跃分子，是水循环模拟中需定量描述的参数。

最小叶孔阻抗以及地表覆盖的空气动力学参数影响水文过程的蒸散发，也是进行蒸散发计算过程中的关键性参数。

植被等地表覆盖物参数根据NOAA影像或者GMS影像进行确定，本书主要是采用国家重点基础研究（973）发展规划项目“黄河流域水资源演化规律与可再生性维持机理”中选用的参数，然后在模型中进行率定。模型参数见表8-3 和表8-4。

表8-3　植被参数

续表

表8-4　地表覆盖的空气动力学参数(www.chuimin.cn)

3.土壤参数确定

土壤层厚度是和植被根系深度互相制约的参数，其值的确定直接影响土壤持水率，进而影响蒸散发。对于集总式概念模型，土壤层厚度所带来的影响隐含在模型平衡中，而对于具有物理机制的分布式模型中，土壤层厚度则需要作为一个可变参数进行考虑。土壤层厚度的初值可以采用两种方法确定：一种是假定土壤层厚度的分布受植被分布的影响，即在同一单元内，土壤层厚度等于最大根系深度；另一种是根据浅层土壤的分类来确定。本书土壤层厚度设定主要采用第二类方法，利用《中国土种志》上的“统计剖面”资料，按照国际土壤分类标准进行重新分类，分为四类十二种，四类为砂土类、壤土类、黏壤土类以及黏土类，然后按照土壤分类确定土壤层厚度，并将最大根系深度作为限制性条件。

土壤水分特性参数影响水分循环，在WEP L模型中，根据黄土高原土壤水分运动实验研究成果，对土壤水分特性参数进行了概化计算，推定出四类土壤的水分特性参数（见表8-5）。为了考虑冰冻对山区土壤导水系数的影响，假定土壤导水系数在0℃以下随温度绝对值的增加呈指数递减。

另外，在蒸发计算及热循环计算中，用到的土壤及其他介质的热力学特性参数，根据土壤成分构成和土壤含水率，通过建立经验公式考虑土壤的热力学特性参数随土壤湿度的变化（见表8-6）。

表8-5　土壤水分特性参数

表8-6　土壤及其他介质的热力学特性参数

4.地下水含水层参数

平原区地下水含水层的渗透系数、传导系数和给水度等参数参照《全国水资源综合规划》确定，然后进行模型率定。而山区饱和土壤层的参数结合径流模拟验证推定：渗透系数取表层土壤饱和透水系数×横向不均匀系数（一般取5），传导系数定为渗透系数×土壤层厚度（m），给水度定为0．05。

（三）参数率定

1.灵敏度分析

参数不确定性依赖于模型结构，而气候、土壤、土地利用、水文以及地理等空间信息获取的难度增加了这种不确定性。目前，预测模型不确定性最普遍的方法是灵敏度分析。参数灵敏度分析的目的是分析系统参数对模型输出的影响因子，从而衡量参数对物理过程的重要性，对系统参数进行筛选，便于模型校正和参数评估。

模型灵敏度分析是建立、改进、检验以及校正水文模型最有效的方法。Lane与Richards以及Rabitz指出灵敏度分析的作用：①能够鉴别模型中最灵敏的参数，进而简化并促进模型的校正，或者使未来研究或者现场测试更有针对性；②说明模型对参数值以及边界条件的典型改变的反应是否是现实的；③证明模型的概念十分敏感，能够代表现实自然系统的行为；④将模型缩减到最精炼的结构。Saltelli指出的灵敏度分析能够增强模型的可信度以及预报功能。可见，灵敏度分析对于模型的校正以及分析是非常必要的。

在模型计算过程中，没有必要也不可能对所有参数都加以考虑，只需通过对敏感性较大的主要参数进行准确评估，而对于敏感性较弱的参数，可以根据资料以及经验粗略得出的取值范围选出，即可建立合理可靠的系统模型。

在复杂系统中，灵敏度指标（一般取为一阶灵敏度系数，即系统输出对系统参数的一阶导数）常常无法直接计算得到。常用的简化计算方法是扰动分析法，即对某一系统参数进行微小摄动，同时固定其他参数取值，进行系统计算，得到相应系统输出，然后采用差分计算得到灵敏度大小。灵敏度计算公式如下：

本章主要对WEP L模型的三类参数进行灵敏度分析。模型中的输出项比较多，选取与水资源量有直接关系的分项进行研究，包括河川径流量、地下水补给量以及蒸发量等，研究模型不同参数变化对各水资源量的影响。根据参数灵敏度因子的相对大小，本章将参数按高、中、低灵敏度进行区分。高灵敏度参数为灵敏度因子大于0．1 的参数；中灵敏度参数为灵敏度因子介于0．01 和0．1 之间的参数；而低灵敏度参数为灵敏度因子小于0．01 的参数。

根据模型计算可以得出：低灵敏度参数包括植被覆盖度、叶面积指数、植被高度、最小叶孔阻抗、空气动力学参数等地表覆盖物参数，河道和坡面曼宁糙率、河道形态、土壤初始含水率、田间持水率、地下水初始水位以及土壤饱和导水系数等；中灵敏度参数包括植被根系深度、土壤水分—吸力特征曲线参数、土壤入渗湿润峰吸力、土壤水分—导水系数关系参数、地下含水层厚度、含水层的给水度以及渗透系数；高灵敏度参数包括降水、温度、各种土地利用类型的洼地最大储留深、土壤最大含水率、土壤层厚度以及河床材料导水系数（见表8-7）。

表8-7　WEP-L分布式模型参数相对灵敏度分析表

续表

针对以上分析，本章选取灵敏度较高的参数分析它们对模型输出的影响。含水层渗透系数属于地下水模拟中比较重要的参数，虽其灵敏度处于中等，但也选取了此参数。结果表明，各参数对模型输出量的影响不尽相同，而且同一参数对不同输出分项的影响也存在差异，如图8-11 所示。

图8-11 （a）反映的是地表洼地最大储留深与各输出分项的关系，横坐标为地表洼地最大储留深变化率，纵坐标为各输出项的变化率与地表洼地最大储留深变化率的比值。从图中可以看出，地表洼地最大储留深与河川基流量、地下水补给量以及蒸发量呈正相关关系，而与地表径流量及河川径流量呈负相关关系。其中，对地表径流量的影响最大，而对地下水补给量影响最小。可见，地表径流量变化对该参数较为敏感，而地下水补给量相对较弱。地表洼地最大储留深的增加使得地表截留量增加，相应减少了地表径流量；同时延长截留水量的滞留时间，促进了截留水量的入渗和蒸发，增加河川基流量。

图8-11（b）反映的是土壤最大含水率与各输出分项的关系。土壤最大含水率的增加有利于土壤持水量的增加，则地表水向土壤水垂向补给增加，而土壤水向地下水的垂向补给减少，使得地表径流和河川基流减少，而土壤水的蒸发增加。可见，土壤最大含水率与蒸发量呈正相关关系，而与其他分项呈负相关关系，其中，对地表径流量的影响最大，而对蒸发量的影响最小。

图8-11　伊洛河模型主要参数对输出项的影响

（a）洼地储留深；（b）土壤最大含水率；（c）河床材料渗透系数；（d）潜水含水层的渗透系数；（e）土壤层厚度；（f）年平均降水深；（g）年平均温度

图8-11（c）和（d）反映的是河床材料及含水层渗透系数与各输出分项的关系。渗透系数的增加有利于水量垂直向下补给，相应河川基流量增加，蒸发量和地表径流量减少。河床材料和含水层的渗透系数虽然对河川径流量、地表径流量、河川基流量以及蒸发量的影响相同，但在影响强度上存在较大的差别；地下水补给量与河床材料渗透系数呈正相关关系，而对含水层的渗透系数响应不大。

图8-11（e）反映的是土壤层厚度与各输出分项的关系。土壤层厚度与蒸发量呈正相关关系，而与其他分项呈负相关关系，其中，对地表径流量的影响最大，而对蒸发量的影响最小。

图8-11（f）表示降水与各输出分项的关系。从图中可以看出降水与各输出分项呈现正相关关系，其中对地表径流量的影响最大，对河川基流量的影响最小。可见，降水作为水文循环的主要环节，对各水资源量的变化起着关键性作用。

图8-11（g）反映的是温度与各输出分项的关系。从图中可以看出，温度与蒸发量呈正相关关系，而与其他分项呈负相关关系，其中，对地表径流量的影响最大，而对蒸发量的影响最小。可见地表径流量对温度的变化敏感性较强，而蒸发对温度的变化敏感性相对较弱。温度的增加有利于蒸发，包括地表截留蒸发以及土壤蒸发，在整个水循环中，蒸发增加使得其他水量相应减少，与温度呈负相关关系。

在灵敏度较高的因子中，降水和温度受气候、区域位置以及区域地理特征等因素的影响，可以根据水文站资料和气象站资料进行确定。土壤层厚度不确定性较高，在本模型中土壤层厚度根据土壤类型的分布确定，同时根据植被根系深度进行约束。因此本模型参数的率定不考虑降水、气温和土壤层厚度这三个参数，而仅对地表洼地最大储留深、含水层的渗透系数、土壤最大含水率、河床材料渗透系数等四个参数进行参数优选。

模型参数的率定指提供给模型研制具有代表性的输入、输出资料、调整参数，确定一组最优化的参数，使模型拟合实测资料最好，达到最优化。参数率定的方法有很多，常用的有人工优选和自动优选。自动优选往往采用遗传算法、罗森布瑞克法、单纯型法以及SEC法。尽管自动优选法具有不依赖于参数初始值、能够在较短的时间内达到全局最优点、精度较高以及弥补了工作人员缺乏经验的不足等优点，但是自动优选法的基本思想认为在特定模型结构下只有唯一一组最佳参数与之对应，然而由于模型结构的复杂性与数据的不确定性，自动优选法通常不能寻优到模型的唯一真值，并且也无法判断算法是否达到全局最优，不能为深入研究复杂模型提供有效途径，因此它还不能完全替代人工优选方法，需要与人工优选方法结合使用。对于本模型来说，系统庞大，计算时间长，采用自动优选参数速度慢，难以满足计算的要求。因此采用人工优选法进行模型参数的率定，并充分借鉴前人成果以及综合规划成果中确定的参数成果，弥补个人在参数选择中的随意性，以增加模拟结果的客观性和可信度，同时大大提高计算速度。

2.率定效果评价指标

为评价模型验证的好坏，选择的主要评价标准包括：径流量误差、Nash Sutcliffe效率以及模拟流量与统计流量的相关系数。

（1）径流量误差。径流量误差是整个模拟期模拟径流量与实测径流量差值百分比的绝对值，径流量误差绝对值越小越好。

（2）Nash Sutcliffe效率。Nash与Sutcliffe在1970年提出了模型效率系数（也称确定性系数）来评价模型模拟结果的精度，它更直观地体现了实测过程与模型模拟过程拟合程度的好坏，公式如下：

（3）相关系数。相关系数是对两个变量间关系的量度，考查两个事物之间的关联程度。相关系数的绝对值越大，相关性越强，相关系数越接近于1；相关度越弱，相关系数越接近于0。

通常情况下，相关系数0．8～1．0为极强相关，0．6～0．8为强相关，0．4～0．6 为中等相关，0．2～0．4为弱相关，小于0．2为极弱相关或无相关。其计算公式如下：

式中：rXY为相关系数；N 为系列的样本数；X、Y 分别代表实测系列和模拟系列的数值。

校正准则包括：模拟期平均年均径流量误差尽可能小；Nash Sutcliffe效率尽可能大；模拟流量与观测流量的相关系数尽可能大。

（四）模型验证

为进行模型验证，在1956～2000年共45 年历史水文气象系列及相应下垫面条件下进行连续模拟计算。取1956～2000年为模型校正期，主要校正的参数包括河床材料透水系数、土壤最大含水率、各种土地利用的洼地最大储留深以及地下水含水层的渗透系数四个参数。

选取伊洛河黑石关水文控制站为验证站，分别对有、无人工取用水两种模拟结果进行校验。其中，无人工取用水的校验是将黑石关模拟计算的径流过程与“还原”后的径流过程进行对比；有人工取用水的校验则是将其模拟计算的径流过程与实测的径流过程进行对比。

从无人工取用水模拟的校验结果来看（见表8-8），黑石关1956～2000 年多年平均径流量误差为6．3%，Nash系数为0．84，相关系数达0．93。从各支流校验结果来看，龙门镇径流量误差为6．1%，Nash 系数为0．78，相关系数达0．94；白马寺径流量误差为6．3%，Nash系数为0．84，相关系数达0．92。

表8-8　模型验证结果

从有人工取用水模拟的校验结果来看（见表8-8），黑石关1956～2000 年多年平均径流量误差为2．9%，Nash系数为0．80，相关系数达0．92。从各支流校验结果来看，龙门镇径流量误差为2．6%，Nash 系数为0．75，相关系数达0．90；白马寺径流量误差为7．1%，Nash系数为0．80，相关系数达0．90。

从水资源评价和水资源演变规律研究的角度出发，这个精度验证结果是比较理想的。两种情况下的月径流拟合结果见图8-11～图8-13 和图8-15～图8-18，实测日径流过程见图8-14。

图8-12　黑石关水文站实测月径流量校验结果

图8-13　黑石关水文站还原月径流量校验结果

图8-14　黑石关水文站实测日径流量校验结果

图8-15　白马寺水文站实测月径流量校验结果

图8-16　白马寺水文站还原月径流量校验结果

图8-17　龙门镇水文站实测月径流量校验结果

图8-18　龙门镇水文站还原月径流量校验结果