汉江中下游干流水华关键环境因子识别及阈值分析

田 晶，郭生练，王 俊，张德兵，罗春艳

(1.武汉大学水资源与水电工程科学国家重点实验室，湖北 武汉 430072；2.长江水利委员会水文局长江中游水文水资源勘测局，湖北 武汉 430010)

水华指淡水水体中藻类大量繁殖的一种自然生态现象，是水体富营养化的一种特征。与出现在湖泊、水库等静水生态系统中的蓝藻、绿藻水华相比，河流水华并不常见，但却具有影响范围广、暴发成因复杂、控制难度大等特点[1-2]。汉江中下游近年来水华暴发频次增加，引起广泛关注。自1992年汉江中下游首次暴发硅藻水华以来，水华的发生频次持续上升，持续时间也有所延长，以往发生的河段多为兴隆以下河段，并有向支流蔓延的趋势[3]，兴隆低水头闸坝型水库建成后，2018年在兴隆库区也出现了水华现象。水华发生时水体呈棕褐色，散发腥味，对汉江中下游生产生活产生一定的影响[4]。因此，探究影响汉江水华发生的关键因子，精准防控水华的暴发，是目前迫切需要解决的问题。

以往的研究认为，汉江水华的发生是水文、气象和营养盐条件多方面综合的结果[5]。国内学者围绕汉江水华已开展了大量的研究工作，谢平等[6]指出制约汉江水华发生的关键因子是流量和流速等水文因子，并非氮、磷等水质因子和水温等气象因子；吴兴华等[7]基于2015—2016年汉江硅藻水华的监测数据，得出水华发生的成因是适宜的气候条件、较高的硅氮比和低流量；王俊等[8]基于2018年汉江水华的应急监测数据，提出了严控污染输入、优化水量调度方案、完善管理机制等治理对策。总体来看，研究者基于不同时段、不同断面的监测数据得到的汉江水华的研究结论不尽相同。当前，汉江中下游水华的预警和防治仍面临着一些难点和挑战[2,9]：①汉江中下游河段目前还未建立常规的水华监测体系，应急监测期间的资料长度较短且受到水量调度的影响，汉江水华的暴发成因，特别是预警研究尚未取得突破性进展；②在气候变化和人类活动(如大型水利工程)的影响下，汉江中下游的环境因素复杂多变，以往研究多针对河流水文情势变化对浮游植物生长的单项研究，变化环境影响下的河流水华暴发与预测尚未取得好的效果；③某些重要环境因子一旦超过阈值，藻类在浮游植物群落的占比可能发生剧烈变化。因此，有必要探索汉江硅藻暴发的重要环境因子阈值，为治理汉江富营养化提供技术支撑。

本文基于汉江中下游干流硅藻水华发生前后的水文、水质和水生态监测数据，采用Pearson相关性检验、分位数回归和随机森林模型方法分析水华发生的关键环境因子及阈值，以期为汉江中下游生态调度提供参考。

1 研究区概况与数据来源

1.1 研究区概况

研究区域为汉江中下游皇庄至宗关水厂江段，以仙桃站为基本分析断面。2021年1月18日，仙桃断面疑似发生水华，长江委水文局启动应急响应，于1月19日至2月2日开展应急监测和调度，分别在汉江中下游的皇庄、沙洋、兴隆坝上、泽口、仙桃和宗关水厂6个汉江干流控制断面开展跟踪取样监测，断面位置如图1所示。

图1 汉江中下游干流采样断面分布

1.2 数据来源

采用的水华监测数据均来自长江委水文局《汉江“水华”连续应急监测调查报告》及历史水文资料。皇庄、沙洋、兴隆坝上、泽口、仙桃和宗关水厂6个汉江干流控制断面涉及的关键要素包括环境指标(气象、水质和水文因子)和浮游植物(藻密度、叶绿素a质量浓度)，其中，气象数据来自国家气象信息中心中国气象数据网(http:∥data.cma.cn)，其余数据均来自长江委水文局。

2 研究方法

2.1 分位数回归

浮游植物的生长周期短，在环境因子变化较为剧烈时，其生物量会出现剧烈的时空变化，数据具有较强的时空异质性，不能很好地满足传统最小二乘法中正态、等方差以及无自相关性等假设条件[10-11]。而分位数回归方法适用的条件较为广泛，可以很好地处理数据中的方差异质性和异常值[12-13]。该方法由Koenker等[14]提出，可以根据响应变量的条件分位数对预测变量进行回归，进而得到所有响应变量分位水平上的回归模型。假设Y是一个连续型随机变量，FY(y)是其分布函数，一个特定的y值的发生概率为τ，那么第τ分位数就是指发生概率为τ时的y值，即：

(1)

如果p(Yyp)≤1-p，那么yp被称为随机变量Y的第p分位数。与普通最小二乘法把预测变量的期望平均值作为响应变量的函数不同，分位数回归可模拟一个或更多响应变量的分位数[15]，分位数回归的特征可以度量预测变量对响应变量分布中心的影响, 还可以度量预测变量对响应变量整体分布的影响。

2.2 随机森林模型

由于河流水华的成因复杂，各影响因子之间的相互作用机制尚不明确[9]，因此大多数研究是基于黑箱模型[16-18]，采用数理统计方法或智能算法。在智能算法中，随机森林模型是一个包含多个决策树的分类器，并且其输出的类别是由个别树输出的类别的众数而定[19]。随机森林模型可以产生高准确度的分类器，在决定类别时可以评估变数的重要性，并且学习过程快速[20]。

选取偏差系数作为叶绿素a的拟合优度指标。偏差系数是对实际值与模型模拟值之间离差的一种度量，其计算公式为

(2)

式中：B为偏差系数；yi为因变量的实际值；y′i为因变量的模拟值。

3 结果与分析

3.1 监测结果

3.1.1藻密度

2021年水华监测期间，汉江干流不同监测断面的藻密度如图2所示(1月30日无监测数据)。由图2可知，从皇庄到宗关水厂断面，不同监测断面之间具有显著的空间差异性，其中宗关水厂的藻密度最大，仙桃断面次之，皇庄断面最小，表明汉江中下游的藻密度呈沿程增加的趋势。随着时间的推移，藻密度持续下降，1月24日各断面间的藻密度差异最大。

图2 2021年不同监测断面不同时间的藻密度

根据DB 44/T 2261—2020《水华程度分级与监测技术规程》，各个断面在不同时间的水华等级如表1所示。由此可见，1月23日水华程度最为严重，所有监测断面的水华等级均达到Ⅲ级轻度水华的标准；1月24日后，皇庄、沙洋和兴隆坝上断面的水华等级均降到Ⅱ级及以下；仙桃和宗关水厂断面分别在1月28日和1月27日后保持Ⅱ级，无明显水华。

表1 各监测断面的水华程度分级

3.1.2理化指标

2021年汉江中下游干流水华发生期间，各监测断面的理化指标均值和最大值如表2所示。可见，干流透明度均值基本呈递减趋势，在皇庄断面最高，为90.25 cm，宗关水厂断面最低，为43.60 cm。水温的最大值在沙洋断面，达到12.7 ℃，水温均值在仙桃断面最低，为9.08 ℃，在沙洋断面最高，为10.77 ℃。溶解氧的质量浓度均值和最大值均在兴隆水库坝上断面最高, 仙桃断面均值为12.65 mg/L，皇庄断面均值最低，为11.35 mg/L，这是由于当水体中充满高密度的藻类细胞时，光合放氧会导致水中的溶解氧质量浓度过高[7]。pH均值范围为8.18～8.61，呈现沿程递增趋势，在宗关水厂断面达到最大值。叶绿素a质量浓度均值也呈沿程递增趋势，在宗关水厂断面最高，为70.02 μg/L，皇庄断面最低，为23.10 μg/L。高锰酸盐指数的质量浓度为2.67～4.20 mg/L。所有断面的总磷质量浓度均超过0.03 mg/L，总氮质量浓度均值在1.5 mg/L左右，可见汉江下游河段较高的氮磷质量浓度为水华的发生提供了充足的营养盐条件。

表2 水华发生期间各监测断面的理化指标均值和最大值

3.2 Pearson相关性

以Pearson相关性系数进行双侧检验，分别分析藻密度与监测的气象、水质、水文因子的相关性。由于藻密度仅在水华发生后进行监测，无法体现水华的生消过程，因此选取叶绿素a来表征浮游植物的生物量。由于仙桃断面水华暴发的频率较高且该断面设置了自动监测站，因此重点对该站进行分析。图3为仙桃站叶绿素a质量浓度与藻密度的相关关系，可见，叶绿素a质量浓度与藻密度呈显著线性正相关(R2=0.797)，因此可以用叶绿素a质量浓度来表征浮游植物的生物量。基于水华程度分级标准，当藻密度小于1×107个/L时无明显水华，该藻密度值为水华的临界值。因此依据藻密度和叶绿素a质量浓度的拟合公式，叶绿素a质量浓度的临界值为76.47 μg/L，可以作为水华发生的一项判别标准。

图3 仙桃站叶绿素a质量浓度与藻类密度的相关关系

采用仙桃站2021年1月1日至31日的数据，分别判断叶绿素a与气象、水文、水质因子的相关性。气象因子与叶绿素a的Pearson相关性分析结果表明，叶绿素a与平均气温在0.01的显著性水平上呈正相关关系(R2=0.639)，与平均风速的相关关系(R2=-0.104)和日照时数的相关关系(R2=-0.279)均不显著。

叶绿素a与水位和流量均呈显著负相关关系，R2分别为-0.723和-0.592，与以往的研究结果一致，即汉江水华的发生与枯水期低流量有关[7]。图4为叶绿素a质量浓度与流量的关系，1月20日前，流量呈下降趋势且整体偏低，此时叶绿素a的质量浓度逐渐上升并在20日达到最大值，之后叶绿素a质量浓度持续下降。这是由于在本次汉江水华发生后，自1月19日开始兴隆水库按“冲蓄结合”方式进行应急调度。兴隆枢纽以上河段的水华影响通过丹江口、王甫洲等梯级水库的应急调度解决，1月24日开始丹江口水库日均下泄流量由620 m3/s加大至800 m3/s，并持续6 d。兴隆枢纽以下汉江河段的水华影响在维持兴隆水库出入库平衡的基础上，通过引江济汉工程应急调度解决。

图4 仙桃站叶绿素a质量浓度与流量的关系

叶绿素a与水质理化指标的Pearson相关性检验表明，叶绿素a与水温、pH和溶解氧均有显著的正相关性(P<0.01)，相关性系数分别为0.388、0.526和0.407。表3为1月21日汉江中下游断面总氮、总磷质量浓度变化情况，可以看出：①支流氮、磷营养盐水平明显高于干流，表明支流汇入带来的营养盐负荷可能对汉江干流水质造成一定影响；②干流断面的总氮、总磷质量浓度虽与非水华期无明显差异，但总氮、总磷的质量浓度较高，总氮质量浓度普遍超过1.3 mg/L，总磷质量浓度普遍超过0.03 mg/L，说明此次汉江的硅藻水华是多种环境因子综合作用的结果，即适宜的气象和水文条件加上较高的营养盐所致。

表3 1月21日汉江中下游断面总氮、总磷质量浓度变化

3.3 分位数回归结果

水华发生会影响水体的pH值和溶解氧质量浓度，藻密度增加后，水中CO2被利用而减少，导致pH值上升，同时导致溶解氧增加。因此，pH值和溶解氧因子作为水华发生后的结果，不适宜用来分析其对藻类的影响。图5为仙桃站气温和流量的参数估计结果随着叶绿素a丰度分位数的变化，包括分位数回归斜率及其95%的置信区间，普通最小二乘法得出的斜率及其95%的置信区间。由图5可见，气温在叶绿素a的所有丰度上均具有显著的正斜率；流量在叶绿素a的所有丰度上均具有显著的负斜率，各个因素对叶绿素a的影响作用都是非线性的。以流量因子为例对结果进行说明，最小二乘法得到的斜率始终在-0.2附近，表明是一种均值回归；而分位数回归结果(图6(b))表明，在0.05、0.25、0.50、0.75丰度分位数下，直线斜率分别为0.03、-0.3、-0.18和-0.16，可见在不同的叶绿素a分布条件下，流量对叶绿素a的影响不同(图6)。

(a)气温

(a)气温

3.4 随机森林模型模拟结果

基于实测气温和流量数据，通过随机森林模型预测叶绿素a质量浓度，图7为叶绿素a质量浓度的实测值与模拟值对比。由图7可见，模拟值的相对误差为-3.13%，说明随机森林结果具有可信度。从影响因子的重要性来看，流量和气温因子对叶绿素a质量浓度影响的重要程度分别为55.8%和44.2%。

图7 叶绿素a质量浓度的实测值与模拟值对比

3.5 基于水华防控的环境阈值分析

确定水华开始出现临界藻密度作为临界预警时刻。识别临界藻密度对应的水质、气象和水文因子，作为控制水华的环境因子阈值。

前文分析可知，气温是显著影响叶绿素a含量的因子之一。为了探究气温的阈值，统计2021年1月3 d、5 d、7 d、9 d、11 d滑动积温，分别计算其与叶绿素a质量浓度的相关关系，结果见图8(a)，发现叶绿素a质量浓度与7 d滑动积温的相关性最大，R2达到0.909。拟合叶绿素a质量浓度与7 d滑动积温得到方程如图8(b)所示，7 d滑动积温的阈值为56.02 ℃，即当7 d滑动积温超过56.02 ℃时，将是汉江中下游硅藻水华暴发的危险时间。

(a)叶绿素a质量浓度与滑动积温的相关关系

图9(a)为仙桃站1992—2021年1—3月的平均流量过程与所有水华发生期间的藻密度关系。可见，大多数的水华事件发生在平水年和枯水年，并且对应1—3月的平均流量频率均大于37%。图9(b)为历年1—3月的最小7 d平均流量，2008年以前最小7 d平均流量低于550 m3/s的年份都发生了水华，但2008年后发生水华的最小7 d平均流量有了明显上升。有水华年份的最小7 d平均流量明显低于无水华年份，当最小7 d平均流量高于800 m3/s时，从未发生过藻类水华，当其高于700 m3/s时，有时会发生水华(2009年)。仙桃站1992—2021年1—3月的最小7 d平均流量的范围为364～1 164 m3/s。以50 m3/s为步长，统计在不同流量区间内发生水华的次数，据此得到不同流量下发生水华的概率，见图9(c)。可见，流量与不发生水华概率的拟合相关性较好，R2达到0.922 3。通过拟合关系推求得到了不同最小7 d平均流量下不发生水华的概率，基于1992—2021年的样本数据，仙桃站在60%、70%、80%、90%不发生水华的概率下，对应的最小7 d平均流量应分别高于310 m3/s、470 m3/s、610 m3/s、780 m3/s。

(a)历年1—3月平均流量与藻密度关系

4 结 论

a.叶绿素a质量浓度分别与气温、流量因子有显著的正相关性和负相关性，但各个因子对叶绿素a质量浓度的影响都是非线性的。

b.基于筛选出的环境因子，叶绿素a质量浓度模拟值与实测值的相对误差值为-3.13%，具有较好的拟合优度。流量和气温因子对叶绿素a影响的重要程度分别为55.8%和44.2%。

c.基于仙桃站水华防控的环境阈值分析表明，硅藻水华暴发的水质、气温和流量临界预警时刻为：叶绿素a质量浓度高于76.47 μg/L，7d滑动积温超过56.02 ℃，最小7 d平均流量低于780 m3/s。