基于WEP的横断山垂直带水文过程模拟研究

苏辉东，贾仰文，牛存稳，倪广恒，刘佳嘉，刘 欢，杜军凯

（1.中国水利水电科学研究院 流域水循环模拟与调控国家重点实验室，北京 100038；2.清华大学 水利水电工程系，北京 100084）

1 研究背景

我国是一个山地多而人均水资源少的国家。我国70%的陆域面积为山地，40%多的人口生活在山区。而我国又是一个人均水资源短缺的国家，人均水资源量仅为全球1/4的水平，且绝大部分水资源分布在山区[1]。横断山是长江、怒江、澜沧江等江河重要的水源地，仅长江上游流域（宜昌以上）多年平均径流量约为4515亿m3，约占长江平均总水量的46%。以横断山为代表的亚洲高山区为“亚洲的水塔”，对数亿人的持续供水做出了重大贡献[2-3]。开展山地水文过程模拟及水资源演变规律研究，可为我国水资源安全保障提供支撑，具有重要科学意义。

山地水文过程及演变规律研究需要用到水文模型。水文模型可以分为集总式模型（如：TANK模型[4]、新安江模型[5]）和分布式水文模型（如：SHE[6]、TOPMODEL[7]、SWAT[8]、VIC[9]、WEP（Water and Energy transfer Processes）等）。山地水文过程模拟面临很多复杂问题，主要包括时空尺度问题和不均介质问题。时空多尺度一直以来是山地水文过程研究的重点问题，也是难点问题。在不同水文尺度上，其水循环和水文要素表现不一。针对这个问题，国内外学者进行了许多研究。Samuel和Sivapalan[10]在澳洲的三个不同流域分析了场次降雨、季节降雨、年内降雨、年际降雨等不同尺度对洪水频率曲线的影响规律。Liu等[11]使用EEMD法从3年、5年、12年和46年等周期分析了雅鲁藏布江流域径流波动和它对气候变化的响应，结果表明年际振荡对雅鲁藏布江的径流变化起着重要作用。徐飞[12]从年季月多时间尺度模拟分析了五台山清水河流域、崇陵小流域、拒马河流域多空间尺度的水循环过程，结果揭示了不同时空尺度下流域水文过程对气候变化和土地利用的响应不同。不均介质问题研究包括大孔隙优先流、喀斯特溶洞、岩石裂隙、基岩凹凸面、土石介质等方面研究。另外动植物影响[13]、干湿/冻融交替作用[14-15]、人为扰动（如耕地等）等因素会增加介质的不均质性。土石二元介质对入渗产流影响具有不确定性，碎石一方面会减少入渗，而碎石与土壤接触部分又会增加孔隙度，从而增加入渗[16]。例如：片麻岩介质和石灰岩介质透水性和持水性各异，对入渗产流会造成不同的影响[17]。受岩溶裂隙[18]影响，喀斯特山地[19]土层薄、蓄水能力弱，补径排迅速，水位流量季节变化剧烈。此外，土壤斥水性[20]也影响土壤入渗产流，斥水性越大，入渗率就会减小。

国内外已在这方面开展了一些相关实验和分析研究，集中在对野外小流域典型观测[21]，通过分析土壤含水率、蒸散发量等因素，研究山地流域不同山坡结构、地质和地形对径流的影响。就山坡水文过程与产汇流机制研究总体而言，需要综合考虑气候、地形和地质条件约束。目前，尽管国内外学者提出了众多分布式水文模型，但针对山区尤其是山坡产汇流、无资料地区水文模拟等问题的研究仍较薄弱。分布式流域水文模型对山地水文过程垂直地带性与演变特征的研究不够，尤其是受山地复杂地形、气候、下垫面条件的限制，给参数方案和模拟结果带来了很大不确定性，相关研究尚不足。

横断山海拔落差大，导致水-热-植被分布具有显著的垂直地带性，有“一山四季”和“一山四带”的现象。气象、植被、土壤、水分具有显著垂直地带性特点，山地垂直地带性导致参数不确定性很大，气候条件恶劣，地形复杂，位置偏僻，导致数据缺乏（1万km2不足一个气象观测站点）。因此，研究横断山垂直地带性对水文过程的影响具有挑战性与科学价值，同时又有很重要的实践意义。本文以横断山垂直带水文过程模拟研究为目的，采用分布式流域水文模型WEP，提出垂直带参数化方案以更好地刻画垂直带水文过程，提高山地垂直带的水文模拟精度。

2 研究区概况与研究方法

2.1 研究区概况横断山（24°35′—33°34′N，96°56′—104°30′E）位于我国西南部，青藏高原东南部，面积约50万km2，其中90%以上为山地。横断山位于我国第一、第二阶梯的交界地带，新构造运动活跃。横断山海拔高程差值大，导致水分、温度、植被垂直地带性显著。横断山区自西向东有怒江、澜沧江、金沙江、雅砻江、大渡河、岷江等河流，其分布如图1所示。横断山地区山高谷深，水流丰沛，主要河流均发源于青藏高原地区，大江大河主要为自北向南流。

图1 横断山流域边界及水系分布

2.2 研究方法WEP分布式流域水文模型[22]，能够用于模拟流域内水、热（通量）的运移过程。采用子流域套等高带作为基本计算单元[23]、并开发了流域二元水循环模型[24]、嵌套了跨区域水资源配置耦合模型[25]，使得模型日趋完善。WEP模型先后在日本、韩国、中国等各种尺度、复杂地质气候条件下的流域应用，取得了良好模拟效果。WEP模型结构如图2所示，水循环、能量过程的模拟方法详见参考文献[20]。

2.3 数据来源横断山分布式流域水文模型的基础数据主要包括以下五大类：（1）水文气象数据主要来自观测站点；（2）地理高程及地形数据；（3）河网水系数据；（4）土壤信息及水文地质等数据；（5）土地利用及植被覆盖数据。其数据来源汇总如表1所示。

图2 WEP模型的水平结构与垂直结构

表1 数据资料描述及其来源汇总

3 WEP模型构建与垂直带参数化方案

3.1 横断山WEP模型的构建

3.1.1 流域边界确定 横断山区只是一个地理单元，并不是一个完整的流域。要构建分布式流域水文模型，需要构建一个或者若干个闭合的流域，使得这些流域能够全部或者大体上囊括横断山区。根据横断山的三级水资源分区，需要适当扩大模拟域的边界。横断山区东北角的小部分区域为白河水系，其属于黄河流域，南部有小部分区域属于元江—黑水河，这两部分没有进入本次的研究区范围。本次构建的横断山流域面积约85万km2，主要包含“横断六江”和“长江源区”。横断山流域分布如下图1所示。

3.1.2 河网提取 WEP河网包括实际河网和虚拟河网。实际河网由测量所得，包括五级河流（如图3所示）。虚拟河网提取主要包括以下步骤：对DEM进行裁剪、重分类、填洼等处理；根据实测的河网水系，对DEM进行降高程处理；计算“流向”、“汇”、确定积水阈值，确定流域水库、湖泊、流域出口等位置，制取分割点；提取计算得到虚拟河网。横断山流域的河网如图3所示。

图3 基于WEP提取的实际——虚拟水系

3.1.3 计算单元划分 为既满足垂直带的模拟分析要求，也不造成过重的计算负担，将横断山流域85万km2的面积划分成了具有拓扑关系的3089个子流域，再按照WEP的等高带划分规则，细分成14 995个基本计算单元。

3.1.4 输入数据制备 WEP输入的数据主要包括DEM、径流数据、气象数据、土地利用和土壤数据，其中DEM和气象站点分布如图4所示。

图4 横断山DEM和气象站点分布

3.2 垂直带参数化方案的改进垂直带参数化方案主要是针对山地气象、植被、土壤参数的垂直带的空间分异而进行处理的，旨在加强对山地垂直带的刻画和模拟，达到提高模拟精度和提升水文过程刻画的准确度。本文提出的垂直带参数化的步骤包括：划分获得山地垂直带谱；根据WEP及三级水资源区划对参数进行分区，通过ArcGIS将参数分区和垂直带嵌套，获得参数编号；通过卫星遥感、站点观测、野外试验、模型经验等方法确定垂直带参数值。

3.2.1 横断山垂直带划分 根据横断山的子流域形心点的高程，对横断山区进行垂直带划分。横断山（291～7556 m）一共划分成6个垂直带，分别为：0～1000 m带，1000～2000 m带，2000～3000 m带，3000～4000 m带，4000～5000 m带及＞5000 m带。

3.2.2 垂直带参数分区 根据横断山所属的二、三级水资源分区和横断山6个垂直带对横断山流域进行参数化分区，构建WEP模型，将横断山区分为9个基本参数区，其分区结果如图5所示。

图5 横断山流域海拔高程（DEM）及子流域划分

3.2.3 垂直带参数值的确定

（1）气温的垂直带修正。WEP采用中国地面气候资料数据集中的降水、气温、日照、风速、相对湿度等5类数据进行空间插值。受高程差大的影响，简单的不考虑海拔变化的插值，不再适合横断流域气象数据空间插值展布，因此需要进行高程的修正。降水的修正比较复杂，本文主要针对温度进行修正。温度的垂直带修正步骤如下：

将气象站点的温度通过气温～高程的关系（海拔每上升100 m下降0.6℃）将气温进行换算到参考的高程点温度。本文将海平面作为参照点时，温度修正为：

根据RDS反距离平均法计算插值点的参考温度。结合RDS法计算公式如式（2）和式（3）：

将得到的气温换算到子流域形心点的气温，将该气温作为输入温度，展布到整个子流域：

（2）植被参数的垂直带方案。为了得到植被垂直带的参数方案值，本文先用遥感数据进行分析，然后结合实测数据进行校正、规律分析总结，得到概化后的垂直带植被参数值。

本文基于遥感（Modis）解译的横断山流域多年平均月叶面积指数LAI和NDVI进行分析，以100 m为等高带，分别统计该等高带内的平均LAI和叶面积指数veg。其结果分别如图6和图7所示。其中Modis遥感数据没有植被覆盖度（veg）产品，需要根据NDVI数据估算区域植被覆盖度[26]。植被覆盖度veg的计算公式：

图6 基于Modis的叶面积指数（LAI）随高程变化结果

图7 基于Modis的植被覆盖度（veg）随高程变化结果

式中：veg为像元内的植被覆盖度；N为像元NDVI数值；Nveg为像元全由植被所覆盖下的 NDVI；Nsoil表示像元全由土壤覆盖下的NDVI，Nsoil会随时间变化，一般在-0.1～0.2之间[27]。根据NDVI制成的横断山流域多年平均月植被覆盖度（veg）分布如图7所示。这里需指出，概化后的垂直带参数化表仅供参考，在率参时还需适当调整参数。

最终，得到了横断山垂直带植被叶面积指数（LAI）参数化的结果（表2）和垂直带植被覆盖度（veg）的参数化结果（表3）。其按照高程分成了6个参数化带，且每个带都概化了参考的1—12月际变化值。

表2 横断山垂直带植被叶面积指数（LAI）参数

表3 横断山垂直带植被覆盖度（veg）参数

（3）土壤参数的垂直带方案。土壤具有垂直地带性，土壤参数对水文模型的模拟也影响很大。土壤层厚度和饱和含水率是高敏感性影响因子，饱和导水系数和湿润锋土壤吸力是中敏感性影响因子。横断山流域土壤类型分布见图8。

图8 横断山流域土壤类型分布

垂直带土壤参数确定的基本步骤是：确定垂直带上土壤的主要质地组成，其方法如式（6）所示；建立土壤质地和土壤水力特征参数（如土壤层厚度和饱和含水率等）之间的联系或者一一对应关系；再结合野外实测的土壤电导率、土壤温度、湿度数据进行校验分析；按照本文提出的垂直带参数分区方法，确定每个垂直带上的土壤参数的平均值：

微生物接种剂在农业领域应用技术已逐步成熟，广泛应用于多种牧草、粮食作物、蔬菜和中草药等领域。研究表明，大豆植株接种根瘤菌和假单胞菌，大豆磷含量、铁离子含量分别提高了88.9%、115.7%，且其全碳、全氮和豆血红蛋白的含量以及根瘤数目均有增加[1]。Abbasi等[2]从小麦根际分离出具有分泌IAA功能的细菌并接种于小麦根际，发现该菌对小麦促生效果明显，给辣椒接种芽孢杆菌(Bacillus)后辣椒产量和次生代谢产物(包括吲哚乙酸、铁载体和几丁质酶)含量均可增加[3]。

式中：M，N分别为垂直带上土壤栅格的行和列数；Sandij、Siltij和Clayij分别为砂粒、粉粒和黏粒的百分比含量；θsand、θsilt和θclay砂土、粉土和壤土的各类土壤参数。

再结合野外实测的土壤电导率、土壤温度、湿度数据进行分析。综合上述分析，得到了垂直带土壤的参数化结果，其结果如表4所示。

表4 垂直带土壤水力特性参数初值

3.3 方案验证设计本文选择的验证时间段为1960—2015年逐月天然径流过程（其中参数率定期为1960—1990年，验证期为1991—2015年），选择横断山流域的14个水文站点进行验证对比。分别运用三种参数方案运行WEP模型，进行模拟结果对比分析。14个站点分别为巴塘、直门达、奔子栏、石鼓、泸宁、雅江、甘孜、屏山、夹江、彭山、三头坝、洼里、泸定、华弹，分布如图9所示。

图9 垂直带参数化化方案选择验证的站点分布

为验证垂直带参数化方案的模拟效果，本文设置三种不同的参数化方案进行对照。三种参数化方案分别为：

“垂直带参数化方案”。如前文所述，将流域先划分成若干个参数区，根据高程——空间插值、遥感——野外实测结合确定气象、植被、土壤的概化垂直带参数值，即为垂直带参数化方案。在垂直带参数化方案中，植被变量被当成随高程变化的概化参数，在模型中可以进行调节。

“普通参数化方案”。气象只采用RDS空间插值、土地利用、土壤、植被都是变量，不作参数处理，在模型中不能调节。其获取方法为直接采用Modis数据等进行统计获得、不做参数处理。每个等高带 （计算单元） 具有一个单独的植被、土壤、水文参数。

“均一参数方案”。整个流域不进行水资源分区，全流域用一个相同的参数。对整个流域Modis数据根据土地利用统计获取森林、草地、作物三种类型的LAI和veg，所有等高带具有相同的植被、土壤参数。

4 结果与分析

4.1 横断山垂直带水文过程模拟验证采用水文站控制断面月平均径流量对整个模型进行参数率定和验证，时间段为1960—2015年共56年的连续模拟（部分站径流数据有缺失）。垂直带参数化方案验证站点、验证时段设定和验证效果，见表5和图10。

表5 基于垂直带参数化方案的横断山流域径流模拟结果

图10 基于垂直带参数方案的径流过程模拟结果

在此基础上，利用改进后的WEP模型对横断山区域各水文站点径流过程进行模拟。14个水文站点的模拟月径流与还原径流在率定期的NS效率系数介于0.73～0.84，相对偏差RE范围介于-12.1%～11.3%；验证期的NS效率系数介于0.71～0.85，相对偏差RE范围介于-3.1%～8.7%。因此，总得来说，模型精度尚可，垂直带参数化方案较好地模拟了横断山区垂直带的水文过程。

4.2 垂直带参数化方案对模拟改进效果分析本研究通过设置三种不同的植被参数化方案进行对比，分析植被参数随高程变化对模型模拟结果的影响。各方案分别按月份设置不同的参数集，将这些参数集带入模型中分别进行模拟，不同方案下不同水文站模拟的NS效率系数和RE结果如图11和图12所示。

图11 三种方案的NS效率变化对比

图12 三种方案的相对偏差RE变化对比

从图11中可以看出，“普通参数化”方案和“垂直带参数化”方案NS效率系数基本相同，但前者相对误差略高于后者；“均一参数化”方案的效果要低于其他两种方案。由此可知，通过考虑植被参数的高程影响，能够提高模型模拟效率。

从图11及图12中可以看出，垂直带参数化方案比普通方案的NS效率系数和RE偏差都有所改进，NS提高0.03，RE减少1%。总体来说改进的并不明显，但是也有部分改进比较大的站点。而相比均一化方案，垂直带参数化方案NS提高0.3、RE减少19.6%。垂直带参数化方案的优点：（1）相对于固定的计算，垂直带参数化方案具有更多的灵活性，提高精度；（2）在没有连续年份Modis数据的情况下，可通过高程分组获取森林、草地、作物三种类型的植被参数；（3）参数简化，易于获取；（4）更好地反应了气象、植被、土壤、水分等信息的垂直地带性，弥补了以往水文模型对参数的垂直变化描述不足。

高原山地区域在气候、地形以及垂直梯度变化和地质等方面存在着差异，植被垂直带状分布明显，而原有WEP模型对植被参数的垂直变化描述不足。通过考虑气象、植被、土壤参数随高程变化的空间变异性，可提高山区水文模拟精度。

5 结论

本文针对横断山海拔落差大的特点，构建了分布式流域水文模型WEP，根据横断山气象—植被—土壤的垂直带特征提出了参数化方案，评估了径流模拟精度。主要结论如下：

（1）14个水文站点的模拟月径流与还原径流在率定期NS效率系数介于0.73～0.85，相对偏差RE范围介于-12.1%～11.3%；验证期的NS效率系数介于0.73～0.88，相对偏差RE范围介于-9.3%～11.3%。总得来说，垂直带参数化方案较好地模拟了横断山区的径流过程。

（2）“垂直带参数化方案”比“普通的参数方案”的NS提高0.03、RE降低1%，比“均一化参数方案”的NS提高0.3、RE降低19.6%。通过考虑参数的垂直变化影响，能够提高模型模拟效果。

（3）垂直带参数方案具有以下特点：相对于固定普通参数化方案的计算，垂直带参数化方案具有更多的灵活性，可提高精度；在没有连续年份Modis数据的情况下，可通过高程分组获取森林、草地、作物三种类型的植被参数；参数简化，易于获取；更好地反映了气象、植被、土壤、水分等信息的垂直地带性，弥补以往水文模型对参数的垂直变化描述不足。