河流生态系统是多样化生态系统中最重要的组成之一, 为人类和其他生物的生存与发展提供了重要的生态系统服务^[1]。随着人类活动对河流的干扰日益加剧, 河流生态系统退化已经成为全球性的生态环境问题^[2]。河流生态系统健康的恢复和维持, 已成为流域管理的重要策略, 是实现流域可持续发展的基础^[3]。

水文情势的自然变化对水质、水温、溶解氧、营养盐等物理化学过程以及河流地貌形成过程起到了主导作用^[4—5], 是维持和保护本土物种和生态完整性的关键因素^[6]。因此, 定量描述水文变量和生态变量之间的关系成为国内外研究的热点问题之一^[7—8]。水文改变指标(Indicators of Hydrological Alteration, IHA)方法不仅可以描述流量变化特征, 还建立起了不同水文变量与水生生境、栖息地环境等生态特征的联系^[9]。为了定量评估生态水文改变程度, Richter等^[4]进一步发展了基于IHA的变化范围法(Range of variability approach, RVA)。Yang等^[10]基于遗传算法, 构建了IHA指标与描述生物多样性指标香农指数之间的拟合关系, 定量描述了水文指标与鱼类种群丰度的联系。然而, 33个IHA指标之间存在着自相关性和冗余性^[11], 使得其在水资源管理实践中难以操作和管理^[12—13]。此外, 气候变化耦合人类活动影响引起的环境变化, 使得水文指标的变化更为显著, 导致了严重的河流生态系统退化^[14—16]。因此, 优选表征生态水文信息的关键水文变量, 了解其变化特征和潜在驱动因素, 对流域水资源的有效管理及河流生态系统的保护具有重要意义^[17—18]。

本文以老哈河流域为例, 利用可变下渗容量模型(Variable Infiltration Capacity, VIC)重建天然流量序列；采用主成分分析(Principal component analysis, PCA)方法从IHA指标中筛选生态最相关水文指标(The most ecologically relevant hydrological indicators, ERHIs)；基于“观测-模拟”对比分析法, 定量评估气候变化和人类活动对河道径流、ERHIs等水文变量变化的影响, 以期为河流适应性管理提供基础依据。

1 研究区概况

老哈河发源于燕山山脉北部的光头山, 全长约430公里, 是西辽河的主要源头之一。本研究选取兴隆坡水文站(42°19′N, 119°26′E)以上集水面积作为研究区。研究区面积为18112km², 呈不对称扇形, 介于41°—42°45′N, 117°15′—120°E(图 1)。流域地形以山地为主, 地势起伏较大, 自西南向东北高程由2017m下降为427m。流域地处中温带半干旱大陆性季风气候区, 夏季高温多雨, 多年平均降水量约为411.7mm, 约88%的年降水发生在5—9月份。近50年流域平均气温为7.0℃, 且表现出缓慢上升的趋势。流域内水利工程众多, 包括2座大型水库(三座店、打虎石)、1座中型水库(二道河子)以及30多座小型水库, 总蓄水量超过5亿m³, 主要满足流域内农业灌溉的需要^[19]。其中, 三座大中型水库的基本信息见表 1。

2 数据及方法 2.1 研究数据

本文收集的资料主要分为地理信息数据、水文气象数据和社会经济数据。地理信息数据包括：美国地质调查局(USGS)网站提供的30m分辨率的全球数字高程模型数据(https://www.usgs.gov/)；联合国粮食与农业组织(FAO)提供的土壤类型数据；美国马里兰大学提供的1km分辨率的全球土地覆被类型数据。水文气象数据包括：从内蒙古水文局获取的1964—2016年17个雨量站日降雨数据和兴隆坡水文站的日流量数据, 从国家气象局获取的老哈河流域内及周边6个气象台站1964—2016年数据(http://data.cma.cn/), 包括日最高和最低气温、10m平均风速, 并采用反距离平方加权法对站点数据进行空间插值。社会经济数据主要为中国科学院资源环境科学与数据中心提供的人口和国内生产总值(GDP)数据(http://www.resdc.cn)以及地方统计局提供的粮食产量数据。

2.2 研究方法 2.2.1 天然径流重建

本文采用可变下渗容量模型(VIC)重建天然径流序列。作为大尺度分布式水文模型, VIC模型模拟径流被广泛应用于径流演变归因^[20]、干旱的演变及影响^[21—22]和生态相关径流指标的评估^[23]等方面。模型主要通过读取区域土壤类型和土地利用类型数据, 利用降水、气温和风速等水文气象数据驱动计算产流并汇流至流域出口, 依据基准期的实测流量数据率定和验证模型参数, 进而还原河流天然流量序列。因此, 本文应用Mann-Kendall(MK)趋势检验^[24]、Pettitt突变点检测^[25]和降水径流双曲线法分析流域降水、潜在蒸散发和径流序列的变异特征, 并将研究期划分为基本不受人类活动影响的基准期和受人类活动影响的变化期。其中, 潜在蒸散发量由Penman-Monteith公式^[19]计算得到。

2.2.2 基于PCA的生态最相关指标优选及其改变度计算

水文变化指标(IHA)包含了33个水文变量, 按照其特征可分为5组(表 2), 涵盖了径流量的量级、出现时间、历时、频率和变化率等水文特征。利用IHA 7.1软件, 可计算各水文变量的统计特征。各水文指标更加详细的说明和具体计算方法请参见IHA软件使用手册^[26]。

主成分分析(PCA)是通过降维的方法将一系列相关变量转化为少数几个不相关变量, 并较大程度保留原始数据的关键信息^[27]。本文中, IHA共有33个指标, 且各指标的量纲与数值大小存在明显差异, 因此首先要对数据进行标准化处理。基于特征值≥1和累计贡献率达到70%—90%两个原则, 可确定主成分(PC)的个数, 并将主成分因子载荷矩阵中绝对值最大的变量作为该主成分的解释变量, 从而优选出生态最相关水文指标(ERHIs)。

变化范围法通常被用来评估水文指标的改变程度。该方法以基准期各指标发生概率为75%和25%对应的值作为上下阈值, 可估算单一指标的水文改变度, 并通过计算多个指标的综合改变度反映河流受影响程度。其计算公式如下：

(1)

(2)

式中, D_i为第i个水文指标改变度；D₀为多个指标综合改变度；N_{0, i}为第i个指标在变化期水文指标落在RVA阈值范围内的实际年数；N_e为变化期水文指标落于RVA阈值范围内的预期年数；0≤|D|＜0.33为低度改变；0.33≤|D|＜0.67为中度改变；0.67≤|D|≤1为高度改变。

2.2.3 环境变化对水文指标影响的定量计算

“观测-模拟”对比分析法基于水文模型还原天然径流, 被广泛应用于定量刻画环境变化对径流的影响^[28]。假定气候变化和人类活动是相互独立的, 则气候变化和人类活动对水文变量影响的定量分离如下：

(3)

(4)

(5)

(6)

式中, ΔX_t和|ΔX_t|分别表示水文变量的总变化值和其绝对值；ΔX_c和ΔX_h分别表示气候变化和人类活动引起的水文变量变化量；X_{1, obs}和X_{2, obs}分别为基准期和变化期实测序列水文变量多年平均值；X_{1, sim}和X_{2, sim}分别为基准期和变化期模拟序列水文变量多年平均值；I_c和I_h分别表示气候变化和人类活动对水文变量变化的贡献率。

3 结果与分析 3.1 水文气象要素变异分析与天然径流重建

老哈河流域1964—2016年水文气象要素演变特征如表 3和图 2所示。该流域年降水量和潜在蒸散发量呈不显著下降趋势, 其MK检测量Z值分别为-0.81和-0.70。年径流量下降趋势明显, 其MK检测量Z值为-5.04, 特别是2000年之后, 年径流量下降显著。Pettitt突变点检测和降水径流双累积曲线结果一致表明：流域年降水量和潜在蒸散发量无明显突变特征, 而年径流量存在两次明显的突变, 分别为1979年和1999年。因此, 整个研究期可分为3个阶段：1964—1979年, 基准期；1980—1999年, 变化期Ⅰ；2000—2016年, 变化期Ⅱ。

利用VIC模型还原天然径流, 首先将基准期(1964—1979年)划分为模型率定期(1964—1974年)和验证期(1975—1979年)。利用基准期确定的模型参数, 模拟了流域月平均径流, 其与实测月平均径流的对比如图 3。本文采用相关系数E_CC、Nash-Sutcliffe效率系数E_NSC和相对偏差E_BIAS三个指标评估模型的模拟精度。在率定期, E_CC、E_NSC和E_BIAS分别为0.91、0.83、1.2%；在验证期, E_CC、E_NSC和E_BIAS分别为0.90、0.80、1.6%。以上结果表明, 水文模型模拟的径流序列满足精度要求^[29], 可被用来进行水文指标的定量归因分析, 其中变化期内实测径流与模拟径流的差别可认为是由人类活动引起的。

3.2 生态最相关指标优选及演变特征分析

利用PCA方法计算的33个水文变量特征值及累计贡献率如图 4所示。其中, 第1主成分(PC1)的特征值为15.13, 累计贡献率为46%；前7个主成分的特征值均大于1, 累计贡献率为83%。根据2.2.3中提出的主成分提取原则, 确定PC1—PC7为所需主成分。表 4给出了上述7个主成分的因子载荷矩阵。依据解释变量的确定原则, 最终优选出7个生态最相关水文指标, 依次为最大7日流量、2月流量、上升率、低流量年内发生次数、高流量年内平均历时、6月流量、年最小流量出现时间。

IHA变量之间的高度自相关性, 使得信息冗余问题显著, 影响了IHA在水资源管理实践中的应用^[13]。如表 4所示, 主成分的高载荷值在一定程度上出现了聚类现象, 例如, PC1中, 最大1、3、7、30、90日流量均有较高的载荷值。图 5表明, 同组IHA变量的变异度之间差异较小, 例如, 变量31—33在变化期Ⅰ和变化期Ⅱ的变异度均较大, 从而影响IHA变量的综合变异评估。因此, 优选ERHIs, 减少信息冗余性, 成为准确评估河流水文特征、构建生态水文联系的关键^[30]。本文选出了7个ERHIs, 其中, 年最小流量出现时间在多项研究中被选为ERHIs^[13], 而年极端流量的出现时间正是水生生物进入生命循环周期的重要信号^[6]。此外, IHA的5个组别中, 均有至少一个变量被选为ERHIs。因此, 可以认为本文基于PCA选取的ERHIs是合理的, 并且有效反映了原始数据集的关键信息。

图 6反映了7个ERHIs的时间变化特征, 除年最小流量出现时间呈增长趋势外, 其余6个ERHIs均表现为下降趋势。其中, 6月流量、最大7日流量、高流量年内平均历时和上升率下降明显, 表明流域高流量脉冲在量级、历时和频率方面均减少；2月流量下降趋势不显著, 年最小流量出现时间略有提前, 低流量年内发生次数显著减少, 表明流域低流量现象有所改善。这种变化特征可能是因为流域内众多水库“蓄洪补枯”的运行调度模式。此外, 在变化期Ⅰ, 各指标波动特征明显, 综合改变度为0.45；在变化期Ⅱ, 各指标主要表现为显著降低, 综合改变度为0.74。如图 5所示, 在两个变化期, 相比于模拟数据集, 实测数据集IHA各变量改变度D_i绝对值增加明显。这表明在变化期Ⅱ, 人类活动对IHA变量的影响更显著。

3.3 水文变量变异驱动因素定量归因分析

变化环境下, 河流水文特征发生了显著的改变。如表 5和图 2所示, 流域年径流量下降明显, 特别是在变化期Ⅱ(2000—2016)。相比于变化期Ⅰ, 尽管变化期Ⅱ内气候因素导致的径流下降有所增加, 从-1.0mm减少到-8.6mm, 但人类活动导致的径流下降变化量更大, 从-6.2mm减少到-20.9mm。此外, 变化期Ⅱ人类活动对径流下降的贡献率小于变化期Ⅰ, 主要有两个原因：一是变化期Ⅰ整体径流减少量较小；二是变化期Ⅱ多年平均降水量为385.4mm, 与基准期425.2mm和变化期Ⅰ423.3mm有较大差距, 导致了气候变化贡献率的相对增加。整体而言, 人类活动对径流下降的影响更强烈, 其全变化期平均贡献率为74.3%, 远高于气候变化的贡献率25.7%。

为进一步分析环境变化对流域水文情势的影响, 计算了7个ERHIs的变化量及气候变化和人类活动的贡献率。如表 6所示, 在变化期Ⅰ, 气候变化对2月流量、年最小流量出现时间、高流量年内平均历时和上升率影响显著, 贡献率分别为-105.2%、-536.6%、-125.3%和-128.6%, 对应的人类活动贡献率分别为205.2%、636.6%、225.3%和28.6%；在变化期Ⅱ, 受气候变化影响显著的有2月流量和年最小流量出现时间, 其气候变化贡献率分别为-170.6%和-316.6%, 相应的人类活动贡献率为70.6%和416.6%。其中, 由于研究区内气候变化和人类活动对生态最相关水文指标的影响并不总是同向的, 从而使得共同影响下的变化量小于单个因子影响下的变化量, 造成了贡献率超过100%, 贡献率的绝对值反映了因子对水文指标的影响程度。两变化期对比, 气候变化的贡献率整体降低；全变化期内, 人类活动的贡献率整体高于气候变化的贡献率, 说明人类活动仍是影响流域水文情势的主要驱动因素。

4 讨论

本文优选出7个ERHIs, 其中, 月平均流量(2月和6月流量)主要影响水生生物栖息地的可利用性；最大7日流量, 对重塑河道形态和自然栖息地有重要作用；年最小流量出现时间则反映了生物体生命循环的兼容性以及对生存压力的可预见性；高流量平均历时主要通过洪泛区影响水生生物的栖息地和水流有机质富集；上升率对缓解缺水造成的水生生物胁迫有重要意义。该7个ERHIs与以往研究的结果不完全一致, 顾玉娇等^[31]在老哈河流域三座店水库所在子流域的研究中选取了6月流量、最小30日流量和零流量天数三个主要指标；程俊翔等^[13]在洞庭湖的研究选取了3月流量、6月流量、最小3日流量、最大7日流量、年最小流量出现时间、低流量年内平均历时和流量逆转次数7个指标；Gao等^[7]统计分析了美国189个水文观测点的数据, 选取了2月流量、3月流量、6月流量、11月流量、最小30日流量、最大7日流量、高流量年内平均历时和上升率, 与本文的结果有较大的相似。此外, Gao等^[7]的研究还表明不同流域、不同数据集分析结果存在一定的差异, 但往往包含三个共同的要素：至少一个月均流量统计数据, 至少两个极端事件统计数据代表高和低极端, 以及至少一个与低流量和高流量频率相关的统计数据。这与本文选取的生态最相关指标有很好的一致性。

以往研究表明, 人类活动是影响流域河道径流量^[20]、水文干旱^[32]和水文健康^[33]的主要因素, 综合贡献率为80%左右, 这与本文中人类活动对径流量和水文情势影响特征一致。此外, 本文结果表明, 人类活动对部分生态最相关水文指标变化的影响更为强烈, 远高于气候变化的影响, 特别是2000年以后。这主要与流域内社会经济的发展有关。流域内生产活动主要为农业、畜牧业和采矿业, 粮食产量相较基准期增加约400万t^[33], 灌溉耕地面积增加约31万hm², 且在2000年之后呈现出持续增长的趋势(图 7)；国内生产总值GDP快速增长, 特别是2000年之后, 增幅明显。研究指出, 流域范围内, 人类取用水比例近10年整体呈上升趋势, 平均为70%左右^[33]。此外, 为满足流域用水需求, 三座店水库于2003年动工建设, 总蓄水量为3.69亿m³, 大规模的农业灌溉显著减少了河道内径流量。

变化环境下, 水文指标变化的驱动因素尤为复杂。Yong等^[19]研究表明, 老哈河流域1970—1979年和1990—1999年降水量较大, 其中1990—1999年年平均降水量为475.4mm, 大于1970—1979年的449.4mm；2000年以后, 年降水量较小, 平均为390.2mm。因此, 在变化期Ⅰ, 2月流量、6月流量和最大7日流量变化不大, 甚至有小幅增加。2000年之后, 区域社会经济高速发展, 伴随着降水量的减少, 6月流量和最大7日流量显著降低, 同时人类活动对以上两个指标的贡献率明显增加。此外, 水库“蓄洪补枯”的运行调度模式, 增加了2月流量, 减少了6月流量、最大7日流量和低流量年内发生次数, 改变了年最小流量出现时间。由于社会经济的发展, 人类生产、生活取用水增加, 也影响了水库下泄的总流量, 进而减少了高流量年内平均历时。人类取用水减少了河道内径流量, 未严格遵循生态控制指标的水库调度破坏了天然水文情势, 造成了流域生态赤字的累积^[34], 显著影响了河流生态健康, 使得区域水资源开发与利用以及河流生态保护与修护面临挑战。

5 结论

(1) MK检验表明, 年降水量和潜在蒸散发量呈不显著下降趋势, 年径流量下降趋势明显；Pettitt突变点检测和降雨径流双累积曲线表明, 年径流量序列在1979和1999年出现突变。1980年和2000年均为经济增长速率明显增大的起始点, 表明流域内经济活动显著影响了河道径流。基于VIC模型的径流重建, 在率定期和验证期, 相关系数E_CC均大于0.9, Nash-Sutcliffe效率系数E_NSC均大于0.75, 相对偏差E_BIAS均小于10%, 模拟精度符合研究要求。

(2) 基于PCA方法, 确定了7个主成分, 累积贡献率为83%, 其解释变量分别为2月流量、6月流量、最大7日流量、年最小流量出现时间、低流量年内发生次数、高流量年内平均历时和上升率, 即为7个生态最相关水文指标。在时间尺度上, 除年最小流量出现时间外, ERHIs主要表现为下降趋势。变化期Ⅰ和变化期Ⅱ7个ERHIs综合改变度分别为0.45和0.74, 表明水文情势变化呈增长的趋势。

(3) 归因分析表明, 人类活动是导致径流量下降和流域水文情势变化的主要因素。其中, 全变化期内, 人类活动对径流下降的平均贡献率为74.3%。对于7个ERHIs, 变化期Ⅱ比变化期Ⅰ变化显著, 人类活动贡献率相对增加, 对上述7个变量的贡献率依次为70.6%、-118.8%、-64.7%、416.6%、-100.3%、-82.2%和-90.9%, 表明人类活动在变化期Ⅱ对水文情势的影响比变化期Ⅰ更为剧烈。

本文假设气候变化和人类活动是独立的, 但气候变化和人类活动往往存在一定的关联, 如雨量条件会影响人类取用水特征：地区少雨干旱时, 会使得人类取水量增大, 进一步影响了河道的水文要素。因此, 本文计算的气候变化和人类活动对水文情势影响的贡献率可能与实际值存在一定差异。