曹圣洁，夏 瑞，张 远，李正炎，任逸轩，塔 拉

1.中国海洋大学环境科学与工程学院，山东 青岛 266100 2.中国环境科学研究院水环境研究所，北京 100012 3.辽宁大学环境学院，辽宁 沈阳 110036

作为世界上最大的跨流域调水工程，南水北调中线工程(简称“中线工程”)自2014年12月12日正式通水以来，在很大程度上缓解了我国北方的缺水问题，是实现我国水资源优化配置和经济发展的重要战略性民生工程. 然而，如此规模宏大的水利工程在为供水地区带来经济生产力的同时，也在一定程度上伴随着对调水区下游生态系统的影响[1-5]. 近年来，国家领导人高度重视长江经济带水生态保护修复工作，汉江作为长江的最大支流，中线工程开通后对汉江下游水生态环境的影响成为我国环保和水利部门的重点关注. 汉江在近20年内多次暴发不同程度的硅藻水华事件，严重影响着沿江居民的饮用水安全状况[6-7]，尤其是2014年中线工程正式开通后，汉江下游武汉段已连续3年发生水华事件. 如何维系上游丹江口水库调水良好运行与下游河流生态健康的关系，不仅是国家河流生态健康管理的重要需求，也是当前国际上水文与生态交叉学科亟待解决的难题.

国外水利工程建设和投入相对较早，对水利工程的相关研究较多. 大量研究[8-11]表明，大型水利工程在建设、运行及移除过程中会对下游河流产生较大的物理和生态变化影响，并且比一般中小型水利工程的影响效应持续的时间更长. Matos等[12]基于历史数据分析了水利工程建设对非洲南部赞比西河水文情势的影响，发现水库下泄流量的改变对下游水环境要素存在较大的季节波动影响，导致河流生境发生显著改变；Bonacci等[13]针对意大利德瓦拉河近30年的水文情势开展研究，阐明了水库的修建和运行对河流下游水位、流量及输沙量的影响，从而引起流域多个支流磷负荷改变. 全球有超过一半的大型河流正遭受着高强度人类活动的干扰，河流水生态系统退化不再是单纯的河道内生态演变过程，而是水资源开发和水利工程建设影响下的复合水生态环境问题[14-15]. 由于国外人口数量和资源环境压力没有我国所面临的形势严峻，其水利工程远不及我国的南水北调、三峡和葛洲坝工程等规模巨大，因而在工程对河流水生态环境影响的案例、后果与对策上的研究相对较少[16].

近年来，许多国内专家学者针对中线工程对汉江中下游的影响从不同角度开展了研究. 李雨等[17]采用Mann-Kendall检验法和Spearman秩次相关检验法分析了中线工程调水前汉江中下游的水量变化趋势；彭聃等[18]对中线工程调水前汉江中下游的水质状况及时空变化趋势进行了分析评价；高永年等[19]基于综合评价指标体系方法研究了汉江下游水质、土壤、社会、水生物和水资源等要素与上游调水的响应关系；YANG等[20]利用GAM模型(广义相加模型)建立汉江中下游水华藻类生物量与环境因子的关系，探讨了通过调节丹江口水库蓄水量来控制下游水华暴发的可能性；殷大聪等[21-22]基于汉江下游实测数据和室内试验研究了水华硅藻的生物学特性，并提出预防汉江下游水华暴发的流量阈值范围. 然而，关于中线工程对汉江下游生态环境影响的已有研究多是基于情景假设和规划条件下的预判，针对当前实际调水条件下的研究成果较少. 随着中线工程于2014年底正式运行，科学识别调水前后汉江下游水生态环境特征和响应规律亟待深入开展.

该研究以汉江下游水华暴发最严重的武汉段为研究对象，选取2010—2017年调水前后汉江下游水文、气象、水质和水生态实测数据资料，重点围绕中线工程调水前后汉江下游多环境要素特征和响应规律开展研究，旨在阐明2个关键问题：①中线工程调水前后，汉江下游水文、气象、水质和藻密度发生了哪些变化？②如何科学分离并识别不同调水时期可能导致河流藻密度变化的关键影响要素及其贡献？从而为揭示汉江下游河流水华暴发成因提供依据，为科学识别大型水利工程影响下的河流水生态环境响应机制提供方法和手段.

1 材料与方法

1.1 数据来源

中线工程从汉江上游的丹江口水库引水，自流沿途供水并最终抵达京津地区. 汉江(106°E～114°E、30°N～34°N)流域面积15.1×104km2，干流全长 1 577 km，上、中、下游分别为丹江口以上河段、丹江口—钟祥河段和钟祥—武汉河段. 其中，汉江下游干流河长约380 km，流经天门市、潜江市、仙桃市、武汉市等地，多年平均气温15～17 ℃，气候温和湿润，雨量较为充沛. 汉江下游流域概况如图1所示.

图1 汉江下游流域概况Fig.1 Overview of the downstream of Hanjiang River

研究资料包括水文、气象、水质和藻密度数据. 水文数据选取汉江下游仙桃水文站2010—2017年流量逐日监测数据资料，分别以2010—2014年为调水前水文时间序列、以2015—2017年为调水后水文时间序列开展分析. 气象数据从中国气象数据网(http://data.cma.cn)下载，选取汉江下游武汉市气象站点(站位57494)2010—2017年日累计降雨量资料，其时间序列与水文资料时间序列保持一致. 水质和藻密度数据选取汉江下游武汉市白鹤嘴、琴断口、宗关3个主要供水水源地2010—2017年2—4月逐旬数据资料，在单因子变化特征分析中选取ρ(TP)、ρ(TN)、藻密度和ρ(Chla) 4个指标，其时间序列与水文和气象资料时间序列保持一致.

1.2 分析方法

1.2.1趋势检验方法

Spearman秩次相关检验法[23]用以对各项数据资料的年际变化趋势进行检验. 在置信水平α为0.95下，用|T|值(Spearman秩次相关检验值)来判断各时间序列的变化趋势是否显著. 若|T|≥1.64，表示通过检验，序列变化趋势显著；若|T|<1.64，表示未通过检验，序列变化趋势不显著.

1.2.2相关分析方法

GAM模型用以分析调水前后影响藻密度变化的关键因素及关键因子推移变化情况,能拟合响应变量与多个解释变量之间的非线性关系，确定各解释变量在方程中的重要程度[24-25]，其一般表达形式为

(1)

式中：s(μ)为响应变量的连续函数；μ为响应变量的期望值，该研究中μ为藻密度的对数值；a为恒定截距；fi为各解释变量的平滑函数，描述藻密度和第i个关键因子xi之间的关系.

利用R语言mgcv工具包中的gam函数实现对GAM模型的构建. 模型统计结果用P、edf、R2、Deviance explained等指标表征.P表示相关显著性水平，P≤0.05时，解释变量与响应变量相关性显著. edf表征被估计值的自由度，edf=1时，解释变量与响应变量呈线性关系；edf>1时，其值越大表示非线性影响力越强.R2为回归平方和与总离差平方和的比值，R2越大，回归效果越显著. Deviance explained为模型对响应变量总体变化的贡献率，贡献率越大，模型的拟合效果越好.

1.2.3流量分配计算方法

汉江下游流量分配的不均匀程度通过变异系数和流量变化幅度体现. 其中，变异系数〔见式(2)～(4)〕是反映河川流量年内分配不均匀性的指标，流量变化幅度〔见式(5)〕为河川月均流量最大值和最小值的比值，体现流量年内变化幅度的大小.

(2)

(3)

(4)

Cm=Rmax/Rmin

(5)

2 结果与分析

2.1 汉江下游水文变化特征

从仙桃断面流量变化趋势(见图2)可以看出，2010—2017年仙桃断面年均流量整体呈下降趋势. 中线工程调水前(2010—2014年)年均流量持续减少，并在2014年达到最低点；调水后(2015—2017年)流量有所回升，但仍低于调水前. 根据调水前后仙桃断面多年流量特征值分析结果(见表1)，调水前多年平均流量为 1 125.3 m3/s，调水后多年平均流量为996.08 m3/s，相比调水前下降了11.5%；对调水前后逐月流量进行Spearman趋势检验，发现调水前|T|值为4.72，大于1.64，流量呈显著下降趋势，调水后|T|值为0.09，小于1.64，流量上升趋势不明显.

图2 仙桃断面流量变化趋势Fig.2 Changes of flow rate at Xiantao cross-section

表1 调水前后仙桃断面多年流量特征值与Spearman检验结果

注：*表示在0.05置信水平上显著相关.

此外，枯水期(12月—翌年3月)调水后多年平均流量相比调水前下降了28.5%，流量在调水前后呈下降趋势，调水前|T|为1.22，接近显著水平1.64，调水后变化趋势不显著；丰水期(6—9月)调水后多年平均流量相比调水前下降了32.9%，调水前|T|为3.71，大于1.64，流量下降趋势显著，调水后变化趋势不显著，说明汉江下游流量在水量不同时期均呈现出较为显著的下降趋势.

调水前后仙桃断面流量分配计算结果表明，仙桃断面Cv(变异系数)由调水前的1.66升高为调水后的2.34，Cm(流量变化幅度)由调水前的4.41升高为调水后的5.95，说明调水后流量分配不均匀性增大，流量年内分配趋于不均匀. 综上，调水后汉江下游水文变异程度较大，引起水文情势变化的原因可能是人类活动的干扰.

2.2 汉江下游气象变化特征

从武汉市降雨量变化趋势(见图3)可以看出，调水前后年均降雨量在全年和枯水期呈现上升趋势，在丰水期有所下降. 根据调水前后多年平均降雨量分析结果(见表2)，调水后多年平均降雨量在全年和枯水期相比调水前分别上升了9.7%和24.7%，丰水期变化程度较小；降雨量Spearman检验结果显示，调水前后各个时期的|T|值均小于1.64，变化趋势不显著. 说明汉江下游气象变异程度较小，降雨量在调水前后呈现出不显著的上升趋势，同时期汉江下游径流过程可能受降雨等气候的影响较小.

图3 武汉市降雨量变化趋势Fig.3 Fluctuation of rainfall in Wuhan

2.3 汉江下游水质变化特征

选择ρ(TP) 和ρ(TN) 2个水质因子，分析其调水前后的变化特征和规律. 从水质因子变化趋势(见图4)可以看出，2010—2017年春季水华易发期间武汉市3个水源地的ρ(TP)、ρ(TN)在调水后均有不同程度的减小. Spearman检验结果(见表3)显示，调水前后2个时间序列的|T|值均小于1.64，说明调水前后汉江下游ρ(TP)、ρ(TN)的下降趋势不显著. 对汉江下游水质变化特征的分析结果表明，中线工程开通后汉江下游水质状况有所好转.

表2 调水前后武汉市多年平均降雨量变化与Spearman检验结果

图4 武汉市3个水源地水质因子变化趋势Fig.4 Water quality changes at three water sources in Wuhan

表3 水质及藻类因子变化与Spearman检验结果

注：*表示在0.05置信水平上显著相关.

2.4 汉江下游藻密度变化特征

以藻密度和ρ(Chla)作为藻类变化表征值，分析调水前后的变化特征和规律. 从藻密度变化趋势(见图5)可以看出，2010—2017年春季水华易发期间武汉市3个水源地的藻密度在调水前后略有降低；Spearman检验结果(见表3)显示，调水前3个水源地的藻密度|T| 值均大于1.64，藻密度下降趋势显著，调水后|T|值均小于1.64，藻密度呈现不显著的上升趋势.

从ρ(Chla)变化趋势(见图5)可以看出，ρ(Chla)在调水后大幅升高. 根据藻类因子变化结果(见表3)，白鹤嘴、琴断口、宗关的ρ(Chla)分别升高了80.5%、99.3%和96.2%；趋势检验结果显示，ρ(Chla)在调水前变化趋势不显著，在调水后呈显著上升趋势. 综上，藻类变化特征与水质变化特征差异较大，说明相比较水质变化对水生态带来的影响，汉江下游水华的发生对水文过程改变更加敏感，调水后汉江下游水生态恶化风险仍然较高.

2.5 影响藻密度变化的关键因子分析

为进一步识别调水前后可能导致汉江下游藻密度升高的关键要素，以藻密度为响应变量，利用GAM模型分析中线工程调水前后上述各水文、气象、水质因子与藻密度的相关性(见表4). 由表4可见：调水前，流量、ρ(TP)的P<0.05通过相关性检验，且R2较大，可以认为流量和ρ(TP)是影响藻密度变化的关键因子；调水后，流量、ρ(TP)的P<0.05通过相关性检验，且R2较大，流量和ρ(TP)仍然是影响藻密度变化的关键因子. 其中，流量的edf大于1，与藻密度非线性相关，ρ(TP)的edf等于1，与藻密度线性相关.

从表4也可以看出：调水前流量的R2最大，贡献率为27.7%，对藻密度变化的影响贡献最大；其次是ρ(TP)，R2为0.187，贡献率为20.5%. 调水后流量的R2为0.496，贡献率为65.4%，大于调水前的R2和贡献率，说明调水后流量对藻密度变化的影响比调水前更大；调水后ρ(TP)的影响贡献较小，且与调水前的影响程度持平. 综上，中线工程开通后，受上游调水的影响，汉江下游流量减小对藻密度变化的影响十分明显，而TP等营养盐浓度的影响相对减弱.

图5 武汉市3个水源地藻类变化趋势Fig.5 Algae trends of 3 water sources in Wuhan

表4 GAM模型拟合结果

注：** 表示在0.01置信水平上显著相关； *表示在0.05置信水平上显著相关.

3 讨论

通过研究分析，汉江水生态环境演变是一个受流域高强度人类活动和区域气候综合影响下的复合水生态环境问题. 在全球变暖的影响下，汉江流域暖干气候也在加强，近60年来汉江下游降雨量存在不明显的下降趋势[26]. 针对汉江下游仙桃流量与武汉市降雨量和丹江口水库下泄流量的相关性开展分析(见图6)，结果显示，汉江下游流量与降雨量相关性不显著，而与丹江口水库下泄流量呈现显著相关，说明汉江下游水文情势在中线工程调水前后受到降雨等自然气候条件影响较小，受水利工程等人类活动的影响较大. 受上游丹江口水库调水的影响，调水后汉江下游流量在全年和枯丰期均有不同程度的减小(见表1). 2012年丹江口水库大坝加高蓄水，为中线工程的开通做好了准备，汉江下游流量在2012—2014年呈显著下降趋势(见图2)，这与柏慕琛等[27]的研究结果一致；调水后在引江济汉工程(2014年建成并运行)的影响下，流量得到一定回升，汉江下游流量年内分配趋于不均匀，说明调水工程带来较大的流量波动，甚至影响着下游的流量分配.

图6 仙桃断面流量影响因素的GAM模型分析Fig.6 GAM model analysis of factors affecting the flow rate at Xiantao cross-section

中线工程调水前后汉江下游水质状况良好，调水后ρ(TP)和ρ(TN)均有下降，但ρ(Chla)大幅升高(见表3)，藻密度在调水后存在上升趋势(见图5). 水体中藻类生长超过一定阈值会导致水华的暴发[28]. 通常认为，过高的氮磷营养盐浓度是导致藻类形成的主要因素[29-30]，而汉江下游水质在调水后有所好转，ρ(TP)、ρ(TN)与藻密度呈负相关，造成这一现象的原因可能是支流对下游营养物变化的影响更大[31]，2015年以后丹江库区ρ(TP)和ρ(TN)下降[32]，调水后库区下泄流量将对下游水体产生稀释作用[20]，导致下游营养物浓度降低，而营养盐负荷仍保持在较高水平，引起藻类的增殖和聚集. 除此之外，低流速下营养盐对河流中硅藻生长的影响限制减弱，即使污染较低水华现象依然发生[33]. 2010—2017年的分析结果显示，汉江下游流量在调水后下降，这也可能是导致下游水质转好而水华暴发频次增加的原因.

GAM模型分析结果显示，汉江下游降雨量与藻密度变化相关程度较低(见表4)，气候条件的改变不是调水前后汉江下游水文情势和藻密度变化的关键影响因子. 调水前后影响汉江下游藻密度变化的关键因子是流量和ρ(TP)，其影响贡献依次减小，调水后流量对藻密度变化的贡献率较调水前大幅升高(见表4)，说明受上游调水工程的影响，河道下游水文情势减弱，致使河流呈现出类似湖库的水流缓慢流动的特征，水文要素对藻密度变化的影响显著升高[7,16]，流量等水文条件才是引起汉江下游水华发生的关键要素[34]. 由图7可见，当流量较小时，藻密度随流量的增大而下降，当流量较大时，河流的冲刷作用使上游浮游植物的在水体中的滞留时间减小，并被迅速带至河道下游[20]，可能导致河流下游藻密度随流量的增大而上升，呈现与调水前流量较大时类似的变化趋势. 水质因子中，ρ(TP)与藻密度在调水前后均呈线性正相关，ρ(TN)在调水前后对藻密度的影响均不大.

图7 调水前后流量对藻密度变化的GAM模型分析Fig.7 GAM model analysis of the effects of flow rate on algal density before and after water transfer

4 结论

a) 中线工程开通前，汉江下游水量充足，调水后受丹江口水库下泄流量减少的影响，汉江下游多年平均流量下降了11.5%，流量年内分配趋于不均匀，流量变化幅度增大；降雨量在调水后呈现不显著的上升趋势，汉江下游气象变异程度总体较小，相比自然气象条件，人类活动对汉江下游径流过程的影响更为显著.

b) 调水后汉江下游ρ(TP)、ρ(TN)减小，ρ(Chla)和藻密度显著上升，与中线工程开通前相比，调水后汉江下游水华暴发风险升高，说明调水后受上游来水量减少的影响，河道内水质状况虽未发生明显恶化，但汉江下游水生态恶化风险依然较高，汉江下游水华的发生对水文过程改变更加敏感.

c) 基于GAM模型的相关分析结果表明，调水前后影响汉江下游藻密度变化的关键因子是流量和ρ(TP)，调水前汉江下游流量和ρ(TP)对藻密度的贡献率分别为27.7%和20.5%，调水后贡献率分别为65.4%和20.5%，调水后汉江下游流量对藻密度变化的贡献率显著升高，说明上游调水引起的汉江下游流量减小对水华暴发的影响十分明显，而TP等营养盐浓度的影响相对减弱.