基于BPCC的流域水沙效应与气候波动及覆被变化

气候波动和覆被变化是引起流域径流和泥沙输移在较长时期内趋势性变化的最主要因素。土地利用/覆被变化影响流域内降雨的截留、下渗、蒸发等水文要素及其产汇流过程，从而对流域的水文过程产生直接的影响。在一定的条件下，土地利用/覆被变化对流域水量与水质造成的影响非常巨大，可能引发或加剧洪涝、干旱、河流与地下水位异常变化以及水质恶化[9]，加大洪涝灾害、水源地污染等事件发生的频率和强度。

对土地覆被变化在流域范围内引发的水文效应，目前的研究方法主要有传统的水文分析方法（实验流域法、时间序列分析法[10]）、水文模型方法及模型试验与数值模拟相结合的方法。随着计算机技术的进步，研究流域水文过程及土地利用变化产生的水文效应的数值模型有了长足进步。集总式模型如HBV、CSC、CHARM、新安江、陕北模型等，模型参数的物理意义不强，且将整个流域视为一个单元而没有考虑参数的变化特性和流域的空间差异，适用于土地利用/覆被变化单一的小尺度流域。分布式水文模型更能准确地描述水文循环的时空变化过程，能有效地利用地理信息系统技术、遥感技术和测雨雷达技术提供的大量空间信息，及时地模拟人类活动或下垫面因素的变化对流域水文循环过程的影响[11]。陈军锋等[12]用CHARM模型和SWAT模型模拟了大渡河上游梭磨河流域气候波动和覆被变化对水文的影响，得出气候波动造成的径流变化占3/5～4/5、由覆被变化引起的径流变化占1/5的结论。陈利群和刘昌明[13]将SWAT和VIC模型应用于黄河源区，得出气候变化对黄河源区的径流影响在95%以上的结论。

由于在建立分布式水文模型时对各种参数的影响机理尚难以完整描述，模型参数有很大的不确定性，致使模型的计算结果存在不确定性；另外，相对于影响径流的另一主要因素降雨而言，覆被变化对径流到底起到多大作用迄今并没有较为明确的结论。本节将通过参数敏感性分析消除模型中各参数的不确定性对计算结果的影响，并在此基础上对比分析了覆被变化与气候波动的作用效果。

4.2.2.1 BPCC模型原理

1.模型综述

分布式水文及泥沙侵蚀模型BPCC（Basic Pollution Calculated Center）从流域DEM高程数据开始，以地理信息系统提供的气象数据和下垫面数据（植被、土壤）等作为输入条件，采用TOPAZ模型划分DEM得到子流域。每个子流域对应唯一河段，且被划分为左、右、源三个坡面。当子流域划分足够细密时，三个坡面能够逼近自然坡面的划分方式。由于同一坡面各单元的下垫面条件（土壤类型、植被和土地利用方式等）不一定相同，需根据土地利用方式、土壤类型和植被类型的各种组合将坡面归类为单一的植被、土壤、土地利用方式，这样的坡面单元构成了模型基本计算单元。坡面单元和每个河段作为“元流域”进行产流计算。其中，坡面单元是水文响应过程的核心，由植被截留、地表填洼、地表径流、壤中流和地下水[14]等模块构成，分别用以求解降雨蒸散发、下渗、地表径流、壤中流、潜水出流等水文过程，并在此基础上考虑雨滴溅蚀与细沟侵蚀作用，作为坡面泥沙侵蚀模块进行坡面产沙计算，然后以坡面的产流和产沙过程作为输入，在连接这些子流域的沟道内进行汇流和泥沙输移计算，最后得到流域出口的各质量源输出过程。这样，将流域产汇流过程概化为“坡面-沟道”系统，可以反映流域降雨及下垫面条件的空间变化。采用具有一定物理机理的数学方程来描述产流和汇流过程，使模型既得到了简化又提高了精度，同时保持了分布式模型的优点。计算水流条件的水文模块与坡面产沙-沟道输沙模块相耦合，可以同时计算流域出口流量过程和含沙量过程。模型中各子模型的结构关系如图4.2-1所示。

图4.2-1　分布式水沙模型的结构

2.坡面水文过程描述

坡面单元是流域水文响应过程的最小单元，需要完成产流、产沙计算，是模型的核心部分。坡面单元模型主要由植被截留模块、地表水模块、壤中流模块和土壤侵蚀模块等构成，水文过程涉及散发、入渗、地表径流、壤中流和潜水出流等，水力过程涉及坡面的侵蚀产沙和沟道的泥沙输移。以下主要对坡面径流模型、坡面产输沙模型、沟道汇流模型及沟道泥沙侵蚀模型进行详细论述。

（1）植被冠层截留模型。

植被冠层是影响降水传输的第一个作用层，降水通过林冠层后形成冠层截留和穿透雨，改变降雨特性，削减降雨动能，对土壤水分收支、地表径流、河川径流调节有重大影响。目前冠层截留模型，如Horton模型、Rutter模型和Gash解析模型，多为机理性模型，参数较多，应用受限。在分布式水沙模型BPCC中，穿透雨量由降雨总量和植被冠层的截留容量共同影响，仅考虑了水量平衡而未考虑降雨在叶面上的运动，节约了计算空间。截留容量代表覆被冠层对降水的最大截留能力，受季节、植被种类和叶面积指数等因素的影响，其确定方法主要有浸水-叶面积法、基于野外实验数据的回归法和基于微波衰减技术的遥感法等[15]。

（2）蒸散发模型。

冠层截留的水分、开敞的水面、裸露土壤中的孔隙水或土壤水经植物根系至叶面气孔处的水分等，可转化为水蒸气返回大气中，发生蒸散发。对蒸散发过程的正确理解及蒸散发量的准确估算，是认识气候变化条件下水循环特征的关键问题之一。

潜在蒸发和实际蒸发，是计算蒸散发模型的两个关键问题。经验公式法、水汽扩散法、能量平衡法和综合法等，是计算潜在蒸发的主要方法，本研究选用世界上应用最广的Penman公式法。对于实际蒸散发的计算，包括传统的水量平衡法、波温比和涡度相关法、依据潜在蒸散发量进行换算的经验公式法和基于蒸发互补原理的一系列方法[16]。本研究实际蒸散发由植被冠层截留水分蒸发率、裸土蒸发和植物蒸腾率三部分组成。

（3）饱和-非饱和土壤水运动模型。

经植被截留后的穿透雨到达地表，在毛管力和重力的共同作用下，渗入地表并在土壤孔隙中运移。土壤水分通过土壤空隙的吸收、保持和传递作用，经降雨、蒸发、渗漏等，进行重新分布，当土壤达到饱和时，一部分水通过侧向排水作用形成壤中流，另一部分则受重力作用控制形成地下径流。本研究将饱和带和非饱和带（潜水层）的水分运动统一考虑，同时考虑蒸发、蒸腾、入渗以及水的再分配，确定潜水面位置和压力水头。

（4）坡面径流模型。

当降雨持续时间较长或降雨强度较大，即降雨使土壤达到饱和或降雨强度超过土壤的实际入渗能力时，地表渗透不会发生，多余的降雨会首先填充地表洼地，而后形成沿坡面流动的坡面流。一般情况下，坡面流是形成流域洪峰的主要部分。同圣维南方程相比，运动波模型是一种更好的数学描述方式，是目前坡面流模拟中最常用的方法。Woolhiser和Ligget[17]认为在运动波波数K＞10时，运动波模型可以很好地描述坡面流运动，其后，陈力等证实自然界坡面流的运动波波数一般远大于10。

3.沟道汇流过程描述

经坡面单元产生的坡面流、壤中流和地下水在坡脚出流，后汇入沟道，经过沟道的逐级输移，在流域出口形成质量源的输出过程，即流域的径流过程。沟道汇流的演算可分为集总式和分布式两类[18]，前者仅考虑某个断面水流的时间函数，而后者则可取到沿沟道的若干断面，描述水文要素在空间的分布。汇流演算模型可分为水文学模型和水力学模型两类，前者一般只是时间的函数，而后者考虑了空间因素。

水流在自然沟道中输移和汇聚，往往会因地形地势等边界条件的改变而发生流态的变化，缓流、临界流及急流交替甚至同时出现，在地形跨度较大的流域，甚至伴随间断流的发生，因此需采用考虑空间因素的分布式水力学模型。考虑到沟道的复杂性，河道断面信息的不足，以及计算本身的稳定性要求，在汇流计算时采用了改进的扩散波方法[19]，扩散波方程为

式中：Q为流量；t为时间；x为纵坐标；C为波速系数；D为扩散系数。

Muskingum-Cunge演算公式形式为

其中

式中：C 1、C 2、C 3为马斯京根法流量系数。

马斯京根法的槽蓄系数K和流量比重因子ε分别为

其中

式中：C为波速。

当ε≤0.5时，达到稳定条件。柯朗数Cr=CΔt/Δx越接近1，收敛性越好。

式（4.2-4）适用于任何断面的河槽，断面平均流速U由曼宁公式推求。

4.坡面产输沙模型

土壤侵蚀是流域地貌形态演变的主要过程之一，是造成土壤退化、生态环境破坏、河道萎缩等诸多自然灾害的根源之一[20]。

土壤侵蚀是一个极其复杂的能量耗散过程。具有一定动能的雨滴撞击土壤颗粒，破坏土壤结构，降雨形成的地表径流，在势能作用下携带泥沙沿坡面流向坡下或沟道，形成汇流过程。沟道中，水流由势能转化为动能，引起土壤冲刷，水流条件改变则发生河床淤积。水流以推移或悬移的方式将泥沙一起向下游输送，形成流域产沙。

土壤侵蚀和泥沙输移过程与水力条件和水流过程密不可分，因此，模型的泥沙侵蚀模块应与水流输移模块相一致，即在坡面产流和沟道汇流的基础上，计算坡面的产沙和沟道的输沙过程，最终得到流域出口的含沙量过程。

坡面水力侵蚀产沙过程可分为雨滴击溅侵蚀过程及坡面流的冲刷和输移两大子过程。雨滴击溅主要起破坏土壤结构作用，为坡面流输移提供物质来源。坡面流过程包括冲刷、输移和沉积三个子过程，而侵蚀方式可分为片流侵蚀和细沟侵蚀。坡面流的水力特性是决定其侵蚀产沙过程的最主要因素，此外，还受降雨、地形（坡度、坡长等）、土壤特性、植被和人类活动等因素的制约。

（1）雨滴击溅侵蚀。

降雨是引起流域土壤侵蚀的主要能量来源之一，其对坡面侵蚀产沙的影响主要为：第一，决定坡面径流量，影响坡面侵蚀方式的演变及产沙过程；第二，打击土壤表面，分离土壤颗粒，为坡面流提供泥沙来源；第三，加强坡面流紊动动能，提高径流输移能力。一般地，雨强是影响径流量和击溅侵蚀量的主要因素，击溅侵蚀量还受土壤特性（强度、结皮、矿物成分）、前期含水量和坡度的影响。

（2）坡面流侵蚀。

坡面流侵蚀主要包括片流侵蚀和细沟侵蚀。片流侵蚀是指沿坡面运动的薄层水流对坡面土壤的分散和输移过程，主要发生在坡面上部无细沟区和下部细沟间区，是沟间地泥沙输运的主要动力。坡面水流本身只能输送颗粒较小的悬移质，经雨滴击溅后，推移质方能被坡面流输运。细沟侵蚀是指汇集成股流后的坡面流对土壤的冲刷和搬运过程。细沟中的水流集中，流速及水深增大，侵蚀特性发生本质变化，侵蚀量明显增加。但是，细沟形成的过程和临界状态具有很强的随机性，细沟流态也不稳定，一般很难形成较为成熟的理论，目前仅依据观测和实验数据等得到经验公式。

自然界中，降雨影响片流侵蚀量，而片流侵蚀可直接传送到细沟边壁而引起细沟侵蚀，加之试验资料缺乏，在实际应用中，很难将片流侵蚀和细沟侵蚀划分出明显的界限。因此，模型中采用“侵蚀-沉积”理论[21]，及水流侵蚀能力是水流所耗费能量的函数，而与之携带的泥沙数量无关。水流耗费能量主要来自于水流与坡面之间的剪切作用以及水流的紊动动能。

5.沟道产输沙模型

泥沙在沟道中输移过程中受到水流冲刷，发生浅沟、切沟侵蚀，随水流条件改变，也可在河道中淤积。天然情况下，水流和泥沙在自然输移过程中呈现出强烈的非恒定、非均匀的特性，河道的槽蓄和河床的冲淤变化，使得洪峰和沙峰在传播过程中表现出衰减和恢复等动态特性。模型采用非恒定悬移质不平衡输沙方程[22]，用以客观描述自然河道中水流、泥沙的动态演进过程。

6.模型参数说明

分布式水沙模型BPCC涉及参数众多，但多出现在用以描述水文、水力及泥沙动力过程的数学物理方程中，因而具有明确的物理意义。根据物理过程描述，可分为地形参数、植被参数、土壤水分参数以及土壤侵蚀参数等四大类。在参数的选取过程中，本研究尽量参照国内外已有的数据库参数和已发表的研究成果，同时，对于无法直接获取的参数，则根据经验和模型率定予以确定。

4.2.2.2 模型率定与验证

选取典型代表年的日降雨-径流过程，进一步应用分布式水沙模型BPCC揭示和了解镇江关流域径流在不同年份的变化规律，分析径流和其他水文要素与气候、覆被等影响因子之间的相互关系。

将1993年、1996年和1998年作为率定期，2000年和2001年作为验证期，时间步长1800s，在人工试错法的基础上对模型参数进行自动优化处理。采用Nash-Sutcliffe效率系数（E NS）、相关系数（R 2）判断模型率定和验证结果的精度。如果日地表径流的E NS≥0.5，R 2≥0.6，则认为率定后的参数符合要求。(www.xing528.com)

率定期和验证期的日径流过程、日平均泥沙浓度过程和日输沙率过程如图4.2-2～图4.2-6所示，图中雨强为流域的平均雨强，模型精度评价结果见表4.2-1。

从率定和验证的结果来看，模拟镇江关流域出口日流量及日平均输沙浓度、日输沙率变化过程的相关系数和效率系数均符合要求，模拟精度较高且变幅不大，说明模型参数稳定，模拟结果合理有效。

4.2.2.3 模型敏感参数分析

图4.2-2　分布式水沙模型BPCC率定期计算值与实测值对比（1993年）

图4.2-3　分布式水沙模型BPCC率定期计算值与实测值对比（1996年）

图4.2-4　分布式水沙模型BPCC率定期计算值与实测值对比（1998年）

图4.2-5　分布式水沙模型BPCC验证期计算值与实测值对比（2000年）

图4.2-6　分布式水沙模型BPCC验证期计算值与实测值对比（2001年）

表4.2-1　分布式水沙模型BPCC日径流及日平均输沙浓度、日输沙率过程评价参数

为消除各种参数变化对模型计算带来的误差，在实际论证覆被和气候变化对流域水文的影响时，需首先通过对模型参数的敏感性分析找出最敏感参数并加以优化，在固定优化参数的基础上进行下一步计算。BPCC模型中与覆被相关的主要参数包括地表最大填洼量、冠层截留指数及土壤饱和渗透率。选取1996年序列进行参数的敏感性分析，采用扰动分析法，即以最优参数为基准分别变化±15%和±30%，以考察模型输出年径流量Q（主要考虑水量平衡）、Q max（最大径流）、Nash-Sutcliffe效率系数（E NS）、相关系数（R 2）的变化情况。各参数变化后的径流过程见图4.2-7（为显示变化规律图中仅列出参数变化±30%的结果）和表4.2-2～表4.2-4。

图4.2-7　BPCC模型系数变化对应的径流过程

表4.2-2　BPCC模型填洼系数改变后的径流变化情况

表4.2-3　BPCC模型截留系数改变后的径流变化情况

表4.2-4　BPCC模型饱和渗透系数改变后的径流变化情况

考虑流域土地利用/覆被变化对填洼的影响，当雨强超过土壤下渗能力时，净雨开始填注地表洼地，拦蓄的水量即为填洼量。流域地形、土地利用及覆被变化影响洼地的容积、数量、面积，从而影响填洼量[23]。由图4.2-7（a）和表4.2-2中可见，填洼系数主要影响的是峰值，填洼量越小峰值越大，即最大径流量越大，但对径流总量的影响不显著。

降水经过林冠后，林冠拦截部分雨量，削减降雨动能，改变降雨分布格局，影响林下土壤水分配及营养物质的循环，同时，有学者研究认为植被的截留作用只有在对地表径流速度产生影响的小暴雨过程中才表现出来[24]。结合图4.2-7（b）和表4.2-3可见，截留系数不属于敏感参数，不过，在水量平衡研究中，截留起着不可忽视的作用，有研究发现树冠拦截10%～40%的雨量，一般为10%～20%，因地表覆被类型、密度、雨强、蒸发等多种因素而不同[23]。本研究中植被最大截留量等于叶面积指数（LAI）乘以一个特定的存储值，虽然截留系数不是敏感参数，但是变化值较大时仍能在一定程度上影响径流总量，是模型中重要的影响参数。

由图4.2-7（c）和表4.2-4很容易看出，饱和渗透率是影响计算精度的最敏感参数。饱和渗透率越大，相应产生的壤中流和基流也越大，由水量平衡可知，地表径流的峰值越小，年径流量越大。在变化30%情况下总径流改变量达21.3%，是模型中最主要的影响参数。

4.2.2.4 覆被变化与气候波动的流域水沙效应

气候与覆被状况（下垫面条件）是影响洪水径流及流域产沙的两个最主要条件，但多数研究只关注其中一个影响因子。本研究模拟了镇江关流域20世纪80年代以来的径流变化过程，同时研究了下垫面条件的变化与气候波动对径流和泥沙的影响以及各自的贡献，初步揭示了该流域径流变化的基本规律。

在BPCC模型中，与气候相关的因素包括降雨和温度，覆被条件包括土地利用方式、土壤属性和NDVI（normalized difference vegetation index）值（由于模型三种因素相互作用，模型将其综合考虑为覆被条件）。为了避免采用单一年份计算引起的偶然性，模型统计了1980—2003年期间镇江关流域的降雨、温度和NDVI值与1980年相比较的变化趋势（图4.2-8）。图4.2-9（彩图4）给出了1980年和2000年的土地利用方式，由于土壤属性变化不大，不考虑其对径流改变的作用。按照波动趋势得到24年后的结果，见表4.2-5，表中变幅记为Fi。

图4.2-8　1980—2003年镇江关流域降雨、气温及NDVI值的变化

图4.2-9　镇江关流域土地利用方式对比

图4.2-9（彩图3）与表4.2-5数据显示，1980—2003年期间镇江关流域降雨有所减少，温度上升，植被覆盖增加，由图可见，土地利用中农田略有减少（0.1%），草地向灌木和森林转变（约1%）。这说明经过20多年的森林保护和管理工作，覆被条件有所好转。

表4.2-5　1980—2003年镇江关流域波动前后各因子取值

以1980年资料作为比较标准，按照趋势线计算出24年后的降雨量、温度及下垫面数据。首先计算每个因子变化后的径流及泥沙过程，其次计算两种因子组合变化后的径流和泥沙过程，最后计算三种因子共同变化后的径流和泥沙过程，将每一种工况结果与1980年径流与泥沙过程比较，得到每种工况对径流深改变的“贡献率”和每种影响因子“单位贡献率”，以此作为流域水文和泥沙改变的定量化指标。采用上述方法的计算结果见表4.2-6和表4.2-7。

表4.2-6　BPCC模型不同模拟工况下计算的流域径流深变化

表4.2-7　BPCC模型不同模拟工况下计算的流域输沙量变化

表4.2-6中数据以1980年工况作为变化前的基数，以三种因子均发生波动的工况为变化后的结果，即D 0=-32.8mm，取各种组合工况的径流深变化量Di与D 0绝对值的比值为该工况的“贡献率”，用以表述影响因子在相同时期内对径流的贡献；降雨、气温和覆被在单独改变工况下的“贡献率”与变幅（实际变化率）的绝对值的比值称为该影响因子的“单位贡献率”，用以表述影响因子在单位变幅条件下对径流的贡献。由表4.2-6得到以下结论：

（1）岷江上游流域1980—2003年的24年来降雨、气温和覆被变化的综合影响使径流深减少了11.5%。

（2）降雨、气温和下垫面三种因素单独改变对径流深变化的“贡献率”分别为39.3%、45.4%、15.9%，气温的贡献率较降雨稍大，覆被的贡献率最小，说明24年来气温升高对流径流改变的作用最大。

（3）三种因素单独改变对径流深变化的“单位贡献率”分别为9.4%、2.4%、12.2%，说明如果发生相同的变幅（波动率），地表覆被对径流的影响是最大的，如果短时间内地貌发生突变（如地质灾害、森林火灾、人类砍伐森林开垦土地等大规模高强度的经济活动），都会剧烈影响流域的径流过程，导致生态环境的恶化。

（4）如果将降雨和温度看作是气候波动，覆被条件看作人为因素，则气候波动和人为因素对径流量的贡献率分别为84.5%、15.9%，人为因素占到总体变化的近1/6，说明通过植被维护、退耕还林等工作改善了该流域的水文条件。

同样，表4.2-7中数据以1980年工况作为变化前的基数，以三种因子均发生波动的工况为变化后的结果，即T 0=-48.18万t，取各种组合工况的产沙变化量T i与T 0绝对值的比值为该工况的“贡献率”，用以表述影响因子在相同时期内对产沙的贡献；降雨、气温和覆被在单独改变工况下的“贡献率”与变幅（实际变化率）的绝对值的比值称为该影响因子的“单位贡献率”，用以表述影响因子在单位变幅条件下对产沙的贡献。由表4.2-7得到以下结论：

（1）岷江上游流域1980—2003年的降雨、气温和覆被变化的综合影响使产沙量减少了62.8%。

（2）降雨、气温和下垫面三种因素单独改变对产沙量变化的“贡献率”分别为44.9%、44.2%、25.6%，降雨的贡献率较气温稍大，覆被的贡献率最小，说明24年来降雨减少对产沙量改变的作用最大，而径流变化的趋势略有不同。

（3）三种因素单独改变对径流深变化的“单位贡献率”分别为10.7%、2.3%、19.7%，与对径流的影响相比，覆被对产沙的影响更为显著。

（4）如果将降雨和气温看作是气候波动，覆被条件看作人为因素，则气候波动和人为因素对径流量的贡献率分别为80.2%、25.6%，人为因素占到总体变化的近1/3，说明人类活动对流域产沙的影响高于对径流的影响，通过植被维护、退耕还林等工作改善了该流域的产沙条件。