不同土地利用方式下漓江流域降水-径流特征变化分析研究

阳 扬，翟禄新，贾艳红，薛开元

（1.广西师范大学环境与资源学院，广西桂林 541004；2.珍稀濒危动植物生态与环境保护教育部重点实验室，广西桂林 541004）

0 引 言

气候和下垫面变化是影响流域水文循环过程改变的主要因素，其中径流变化对不同时空尺度下气候和下垫面变化的响应研究是水资源规划及其管理的基础［1，2］。对于描述自然流域的水文模型，模型参数大多反映了流域土壤和植被等的水文物理特性［3-5］，与流域的土壤、植被、地形地质条件等相关。由于以往假定流域特征是处于稳定状态，即土壤、地形、地貌等下垫面在长时间尺度上是不变的，因而通常率定后的模型参数随时间固定不变。但实际情况是，人类活动通过水利工程建设、改变土地利用和植被覆盖等方式直接或间接影响流域对大气降水的分配以及蒸散发等，从而影响地表径流和地下径流的形成过程，表现出非一致性特征［6-8］。因此，水文学中传统的静态参数假设已不能反映变化环境下流域的特征，水文模型参数应随流域特征产生动态变化。Wagener T 在探讨如何模拟环境变化对水文过程的影响时，提出模型参数与流域特征具有复杂的相关关系，应当是流域多种特征的函数，随着流域特征的变化而变化［9］。Merz R 等对奥地利273 个流域的水文状况进行模拟，发现降雪融雪过程的相关参数、土壤特性参数和产流过程相关参数等均随时间发生显著变化［10］。由此看来，土地利用和植被覆盖与模型参数的动态变化存在相关性，即提供一种思路：基于土地利用指标与模型参数的相关性，可以分析人类活动影响下流域特征变化对降水量-径流机制的影响，从而实现流域变化的动态定量感知。孙瑜对汉江流域和潮河流域开展研究，认为不同模型的同类参数与环境要素呈现出同向的相关关系，如蒸发过程参数均与植被覆盖条件正相关，产流过程参数均与植被覆盖条件负相关［4］。熊立华［11］等考虑植被和人类活动影响的水文模型参数时变特征，认为最大蓄水能力Sc与流域植被覆盖条件、土地利用变化、水利工程修筑等具有较强的相关关系。

本文选用集总式水文模型MODHYDROLOG 对漓江流域径流开展模拟。该模型结构简单清晰、便于率定和计算，而且由于集总式水文模型对基础资料要求较低，可在研究区基础地理数据并非特别丰富的情况下应用，降低了开展水文模拟的难度［12-18］。同时该模型考虑了超渗产流和蓄满产流两种产流机制，因此在干旱、湿润流域都能得以广泛应用。据此，本文在不同时段人类活动和气候变化情景下，通过时变参数分析流域降水-径流特征的变化情况，探索剧烈人类活动影响下流域下垫面变化动态定量感知，为流域非稳态背景下的水文模拟提供新的思路，以提高定量分析环境变化对于流域水文变量影响的精度，揭示变化环境下整个流域气候变化、人类活动、土壤和植被等变化对河川径流的影响，从而为流域防洪减灾、流域生态系统的稳定性评价、水源涵养功能提升等提供科学依据。

1 研究区概况

漓江起源于猫儿山，流经兴安县、灵川县、临桂县、桂林市区、阳朔县、平乐县城之后与恭城河交汇，流域全长约214 km，面积约为1.23 万km2，地理坐标为24°38′10″ N～25°53′59″ N、110°07′39″E～110°42′57″E（图1）［19］。整体地势北高南低，以漓江为中轴线，呈西北东南向分布，北部和东部为碎屑岩中低山地貌，中南部主要为低山、丘陵、岩溶地貌，具有典型的喀斯特地貌特征［20，21］。本文研究区为桂林水文站以上流域范围，集水面积约为2 762 km2。漓江桂林水文站以上流域水利工程较多，已建成的水库主要有：青狮潭水库、川江水库、小溶江水库、斧子口水库等［19］。

图1 漓江流域范围Fig.1 Scope of Lijiang River

2 数据与方法

2.1 数据收集

本文主要收集桂林水文局的桂林水文站2007-2016年日尺度径流数据及基本气象数据。研究还收集了流域内的图形数据、DEM 数据以及相关的规划图等，主要有2020年桂林市行政区划图、漓江流域30 m 分辨率数字高程，使用2008、2010、2013以及2015年共4年桂林市1∶1 万土地利用现状图，其他年份采用年内线性插值方法推求获取。土地利用分类参照表1。

表1 中国土地利用/土地覆盖遥感监测数据分类系统Tab.1 Classification system of land use/land cover remote sensing monitoring data in China

2.2 MODHYDROLOG 模型简介

MODHYDROLOG 模型是一个复杂的地下水补给模型，最初应用于澳大利亚［21，22］。该模型有5 个储存库和15 个参数。其结构及主要参数分别如图2、表2 所示。MODHYDROLOG 模型模拟的径流产出由3个部分组成：超渗地表径流、蓄满壤中流和基流。这3 种径流在MODHYDROLOG 模型中的主要模拟过程为：首先降雨受到截留损失，超过最大截留储量后剩余的降水一部分入渗进入土壤，另一部分形成超渗地表径流和洼地蓄水。入渗的土壤水又会经过3个分配过程：补充地下水量，补充土壤含水量和蓄满壤中流。以上过程中，截留量、洼地容量、土壤含水量都会受到蒸发作用，前两部分蒸发按潜在蒸发能力计算，而在实际计算中，一般用实测蒸发量代替而作为模型的一个输入量。模型原理及说明参见文献［21，22］。

表2 MODHYDROLOG 模型参数及取值范围Tab.2 Parameters and values range of MODHYDROLOG

图2 MODHYDROLOG 模型结构Fig.2 Structure of MODHYDROLOG

2.3 参数优选

上述模型用Nelder-Mead 算法进行参数优化。Nelder-Mead算法是由Nelder和Mead于1962年在Spendley等［23］工作的基础上设计出新的单纯形搜索方法，该方法为求解无约束优化问题的局部搜索算法，并无须目标函数的任何导数信息。对于D个变量的函数最小化问题，Nelder-Mead 单纯形法使用反射、扩张、收缩、压缩等操作，通过比较单纯形的D+1 个顶点的目标函数值，用新的点取代目标函数值最大的顶点，逐步迭代并不断更新，最终逼近问题的最优解。本文在利用系列数据优化评价的基础上，对模型进行逐年优化，推求不同年份的模型参数，以体现参数的时变特征。

2.4 模型评价指标

验证模型准确性和评价模拟效果，可以通过以下统计学参数来反映。

主要有决定系数（R2），纳什效率系数（NSE）［24］，标准偏差（Dv）［25］、克林效率系数（KGE）［26］：

式中：Qs和Qo分别代表第i 个样本的模拟值和实测值；μs、μo分别为模拟值和实测值的平均值；σs、σo分别为模拟值和实测值的标准偏差；r为模拟值和实测值的相关系数。

两种评价因子R2、NSE与KGE数值越接近1，则径流的模拟值与实测值的差值越小，模拟效果越好。通常情况下，R2、NSE与KGE评价因子数值都在0.5 以上的径流模拟就可表明模型精度较好［27］。据Donigian［28］的标准，认为Dv在观测值的10%范围内的模拟已经达到了“非常好”的标准，越接近于0 则表明模拟效果越好。

3 结果与分析

3.1 参数优化及精度分析

率定期选择漓江流域2008.1.1-2016.12.31 的降水、潜在蒸散发、气温等气象数据，将2007全年数据年作为模型的预热期。应用Nelder-Mead 算法对MODHYDROLOG 模型进行参数优化选择，设置待调参数上下边界条件，以克林效率系数（KGE）为率定指标。模型优化结果显示，模型的多年KGE值接近于0.8，通过模拟效率系数检验；将最优参数结果代入MODHYDROLOG模型，对2008-2016年漓江流域径流过程进行模拟，比较分析决定系数（R2），纳什效率系数（NSE），标准偏差（Dv），表明模型的决定系数R2、纳什效率系数NSE均超过0.7，标准偏差Dv均在0.2%以内，达到了“非常好”的标准。因此，该模型对漓江流域日径流模拟精度较高。图3 给出了MODHYDROLOG 模型日径流模拟结果，由图3可知，两种模型的日径流过程与实测径流过程具有很好的一致性，模型在枯水年、丰水年都能够较好的再现日径流变化过程，同时，洪峰的量级和出现时间基本一致，模拟效果令人满意。

图3 MODHYDROLOG 模型径流模拟结果Fig.3 Runoff simulation results of MODHYDROLOG

3.2 蒸散发和径流成分模拟结果

MODHYDROLOG 模型将蒸散发计算分为3层：植被截留蒸散发、土壤蒸散发和洼地蒸散发，总蒸发量为三层之和，计算得到的年蒸散发量见表3。

表3 MODHYDROLOG 模型模拟蒸散发情况 mmTab.3 Simulation of evapotranspiration by MODHYDROLOG

径流成分是由水源划分决定的，MODHYDROLOG 是三水源模型，由地表径流、壤中流和基流3 种径流成分构成。MODHYDROLOG 模型模拟所得各径流量如图4所示。2008-2016年地下水基准面、深层渗漏量、基流变化如图5 所示，该图体现了模拟期内地下水基准面年均值、深层渗漏量和基流年总值之间的相互关系。

图4 模拟径流成分比较Fig.4 Comparison of simulated runoff components

图5 地下水基准面、深层渗漏、基流变化情况比较Fig.5 Comparison of changes of groundwater datum，deep leakage and base flow

3.3 土地利用变化与模型参数相关分析结果

将土地利用相关指标进行统计，通过综合对比模型的相关参数来揭示径流的变化特征，从而反映漓江流域的特征变化情况。如图6 所示，将插值后的土地利用类型及变化与模型参数进行相关性分析。结果表明，不同土地利用类型的改变与参数的相关性差异较大。其中参数sub与草地为极显著相关（p≤0.01）。参数dlev与城乡，k1 与耕地和水域呈显著相关（p≤0.05）。此外，参数crak与林地和草地、em与水域、dlev与耕地、林地、水域，k1与城乡，k3与耕地、林地、城乡相关（p≤0.1）。相关参数表明，土地利用变化通过影响不同径流成分的比例影响产汇流过程，最终作用于模拟径流量。主要是地下水储量及其与河道之间的水量交换等相关参数与土地利用的变化密切相关，与植被作用相关的模型参数也有微弱相关性。而渗透（coeff、sq、vcond）和洼地的相关参数（dsc、ads、md）并无明显显著性。

图6 土地利用变化与模型参数的相关性Fig.6 Correlation between land use change and model parameters

3.4 土地利用方式变化对流域降水-径流特征的影响

3.4.1 耕地面积变化对流域降水-径流特征的影响

图7 表明，耕地面积与参数dlev 表现为显著负相关，相关系数为-0.63，说明耕地会抑制地下水储量基准面抬升，使得地下水深层渗漏量（与外流域的水量交换量）增多，不利于径流的产生；同时，耕地面积与k1、k3 呈显著正相关，相关系数分别为0.72 和0.59，其中与k1 为极显著，说明耕地会促进地下水含水层与河流的水量交换作用，地下水补给河道水量，从而产生更多的径流。此外，耕地面积与insc和em 分别呈微弱的正相关和负相关，反映出耕地面积增加有促进冠层截流和减缓植被蒸腾的作用，一定程度上增加土壤蓄水量和径流量。

图7 耕地与参数的相关性Fig.7 Correlation between cultivated land and parameters

3.4.2 林地面积变化对流域降水-径流特征的影响

图8 表明，林地面积与sub 呈微弱的正相关，因而林地对壤中流的形成起着微弱的促进作用，这是因为林地促进土壤非均质性增加更易形成壤中流［3］。林地面积与crak 呈显著负相关，即表明林地会抑制入渗水量对地下水的优先补给，即在土壤水分达到饱和前进入地下水的水量将会减少。林地面积对dlev、k1、k3 的影响与耕地相同，林地会降低地下水基准面，使得地下水深层渗漏增加，且更易于补给河道，有增加径流的作用。

图8 林地与参数的相关性Fig.8 Correlation between forest and parameters

3.4.3 草地面积变化对流域降水-径流特征的影响

图9 表明与耕地和林地不同的是，草地并没有体现出与地下水相关参数的显著性。唯一的是草地与crak 为正相关，相关系数为0.62，也就是说草地一定程度上促进了在水分入渗过程中补给地下水的比例。同时，草地与参数sub表现出极显著负相关，相关系数为-0.84，即减小了壤中流的形成比例，这与增加地下水补给相一致，也体现了草地与林地的差别。另外，草地与土壤最大储水量呈微弱正相关，显著性系数为0.13，即一定程度上会促进土壤最大储水量的增加。因此，草地较全面的体现了入渗水的3个分配方式的相关性，促进了土壤的储水能力，增加了地下水的补给比例，降低了壤中流形成比例。

图9 草地与参数的相关性Fig.9 Correlation between grassland and parameters

3.4.4 水域面积变化对流域降水-径流特征的影响

图10 表明水域面积与参数dlev呈显著正相关，相关系数为0.59，说明水域面积会促进地下水基准面的抬升，减少地下水的深层渗漏量，有增加径流的作用。与参数k1 和k3 呈显著负相关，相关系数分别为-0.71 和-0.58，其中与k1 为极显著相关性，显著系数为0.03，表现出含水层基流交换减弱，并且补给河道的基流量减小，不利于径流形成。此外，水域面积与参数insc表现出微弱的负相关，而与em呈显著正相关，表明水域面积增加的同时会减小植被面积的比例，从而影响整个流域的截留能力，提供更多的水分供流域蒸发。

图10 水域与参数的相关性Fig.10 Correlation between water areas and parameters

3.4.5 城乡面积变化对流域降水-径流特征的影响

图11 表明城乡面积与dlev呈极显著正相关，相关系数为0.67，即城乡面积有利于抬升地下水的基准面，减少地下水的深层渗漏量，增加河道径流量。城乡面积与参数k1、k3 呈显著负相关，表现出弱化含水层基流对河道的补给，不利于径流形成。城乡面积与insc呈反比，说明对冠层截留有抑制作用，可适当增加净雨，进而形成较多的径流量。

图11 城乡面积与参数的相关性Fig.11 Correlation between urban and rural areas and parameters

总体来说，5 种土地利用类型中除草地外，其余4 种土地利用类型方式都对地下水补给施加不同程度的影响，进而影响基流的水量交换，通过参数dlev、k1、k3 从深层渗漏和地下水补给河流得以表现。土壤渗透的相关参数（coeff、sq、vcond）以及洼地的相关参数（dsc、ads、md）与流域降水-径流特征相关性不显著，即土地利用方式的改变并不会显著地通过改变土壤渗透性和流域内洼地的面积来影响径流形成。因此，土地利用类型对不同水文过程的影响程度不同，因而单一组分对产汇流的影响也存在差异，表4归纳了土地利用类型对相关参数及径流的影响。

表4 土地利用类型对相关参数及径流的影响Fig.4 Effects of land use types on relevant parameters and runoff

4 讨 论

研究结果表明，模型参数经率定后，MODHYDROLOG 模型能较好地模拟漓江流域的日径流过程，NSE 和R2均超过0.7，Dv在0.2%以内。主要由于MODHYDROLOG 模型考虑了超渗产流和蓄满产流两种产流机制，具有一定的物理基础，在枯水期和丰水期都能很好的再现日径流的变化过程，因而能在处于湿润地区的漓江流域适用性较好，一定程度上反映了研究区漓江流域的产汇流特征。但由于漓江流域全年降雨量较大且集中降雨较多，这使得模型主要以超渗产流的机制为主，因此在枯水期存在一定误差。

当耕地、林地面积的减小，草地、水域面积的增加时，与之对应的是漓江流域的地下水基准面的抬升，深层渗漏减少，同时基流交换愈加强烈，地下水大量补充河道产生径流。如图5所示，由于地下水基准面抬升作用，深层渗漏量有减小趋势；同时，地下水与河道间的基流交换有增强趋势，产生更多径流。林地的土壤类型多为不饱和雏形土和疏松岩性土，土质较为疏松，经土壤入渗补充地下水的能力较强，所以地下水补给量高，具有较强的水土保持和涵养水源的功能，也因根系作用易于形成壤中流。草地的土壤类型多为饱和黏性土，此类土表层疏松，但深层的土层黏粒含量相对较高，这使降水通过入渗来补给地下水变得困难，而且草地的根系分布较浅，仅能从浅层土壤中吸收水分，因而与地下水有关参数的相关性并不显著。研究时段内草地面积的变化仅微弱上升并不显著，因此并不能作为反映流域特征变化的依据［29］。另外，水域面积的增加造成更多的水量直接汇入河道，同时加大了下渗后介于地下水和河道之间的基流交换；与之同样效果的是，城乡建设增加了不透水面积，因此更多的降水直接进入河道形成径流。

模拟结果表明，模型中蒸散发和壤中流等相关参数与大部分土地利用指标相关性较弱，一方面源于植被的作用，另一方面也体现了土地利用方式的改变对漓江流域的蒸散发和壤中流等影响不大。此外，土地利用指标与模型参数的相关性分析结果表明，土地利用方式与渗透和洼地的相关参数（coeff、sq、vcond、dsc、ads、md）的显著性水平较低（p＞0.2），并无表现出明显的相关性，与Francis Chiew 认为这类参数与流域特征相关性极弱的研究结果相一致［3］。

通过土地利用变化与模型参数的相关分析，针对不同土地利用类型作用于模型的时变参数的情况，最终解释土地利用变化对流域产汇流过程的影响，即流域降水-径流特征。土地利用变化对流域径流量的影响取决于流域内土地利用变化的程度，由于本文所采用的时间序列较短，土地利用指标的精度不高，可能存在参数与流域特征指标显著性不高的情形。另外，对于参数如k1、k3 等虽与土地利用指标存在显著性，但由于其作为方程的系数且物理意义不太明确，故在解释流域的具体特征变化上略显无力。此外，模型对于地表水与地下水的相互作用能较好的概括，但不一定对流域特征较好地表达；流域中有新建水库参与降水径流过程，但参数中体现不足，可能与运行时间较短有关；在分析土地利用方式变化时，没有考虑不同方式之间的转化，仅单一地考虑某一利用土地面积的变化造成的降水-径流特征的影响。

模型中与地下水相关的参数与土地利用类型指标变化有较强的相关性，而这些参数又将直接作用于降雨转化为径流的整个过程，因而一定程度上解释了漓江流域降雨-径流特征。具体而言，土地利用类型的变化主要与模型中补给地下水、地下水基准面及基流（crak、dlev、k1、k3）有较强的相关性。在概念性水文模型中，如dlev 一类的参数，不仅反映地下水的深层渗漏，也对模拟流域水量平衡方面发挥补偿调节作用［30］，即有些土地利用改变，会造成流域降水径流平衡误差增大时，模型以dlev 参数来进行调节补偿。另一方面，研究区中下游平原区多石灰岩分布，因此可以通过地下暗河发生地下水深层渗漏。

5 结 论

（1）MODHYDROLOG 模型能很好的模拟漓江流域的径流量变化情况，但由于漓江流域降雨充沛且分配不均，主要以超渗产流为主、蓄满产流为辅，同时在降雨强度较小的枯水期容易高估径流量，因而产生偏差。

（2）相关分析表明，模型参数主要是地下水相关参数dlev、k1、k3 与绝大部分土地利用类型变化关系密切，个别土地类型的变化与植被相关参数sub、crak、em 有着不同程度的相关性。主要表现为耕地和林地面积与参数dlev 负相关，与k1、k3 正相关，水域和城乡反之。此外，林地和草地与参数crak 相关，水域与参数em 也达到了相关性。而渗透（coeff、sq、vcond）和洼地的相关参数（dsc、ads、md）并无明显显著性。

（3）对模型参数和下垫面指标进行相关性分析，分析土地利用对产汇流径流成分的影响，主要体现在耕地和林地抑制地下水储量基准面抬升，使得地下水深层渗漏量增多，同时促进了含水层与河道之间的基流交换作用，有利于地下水对河道的补给。草地可促进水分入渗补给地下水的比例，减小了壤中流的形成比例，促进土壤的储水能力。水域和城乡与耕地和林地的作用机制相反，其对基流交换的抑制作用减缓了地下水对河道的补给量。