棉花滩水库浮游植物群落结构演替特征及风险预警研究

崔玉洁，杨友德，郑婉婷，成再强，陈天生，林晓芳

（1.三峡库区生态环境教育部工程研究中心，湖北宜昌 443002；2.福建省龙岩水文水资源勘测分中心，福建龙岩 364000）

0 引 言

水库蓄水在发挥防洪、发电、航运等巨大综合效益的同时，也深刻的影响着水体的水动力、水环境特征［1］。随着蓄水，原河流水体水深加大，水体自净能力降低，易引发水体富营养化，流速趋缓也为漂浮型蓝绿藻提供了适宜的生境条件，引发水华现象。水华的暴发是适宜的气象条件、充足的营养盐和缓慢的水流共同作用的结果［2］。研究表明：随着水库运行年限延长，浮游植物优势藻种易由初期的河流型甲藻、硅藻向湖泊型蓝、绿藻演替，而蓝藻由于其释放的有毒藻毒素，一直是关注的焦点［3］。

水华预警是水华防控的重要手段，其预警模型可分为物理过程模型和数据统计模型两种，其中物理过程模型主要围绕水华形成的生源要素、气候水动力条件，以藻类增殖的过程来预测水华情势的强弱［4］。Kim 等［5］依托EFDC 构建韩国Han River的环境水力学模型，考虑三个优势藻种：硅藻、绿藻和蓝藻，较好的预测了Han River 叶绿素浓度；Zhang 等［6］针对日益严重的原甲藻水华，开发了集物理、营养循环和浮游生物生理学的三维预测模型，再现了美国Chesapeake 湾浮游植物种群越冬的模式。Kerimoglu 等［7］利用自养分类不同的浮游植物群，为法国Bourget 湖开发了一个蓝藻水华预测模型。由于浮游植物增殖影响因素众多，水动力学与藻类生理过程之间作用复杂，使用物理过程模型进行水华预警过程中，参数众多，且需要大量连续的基础数据支撑，如在CE-QUAL-W2 模型中，仅考虑单种浮游植物生长与温度的本构关系，就涵盖8个参数，复合藻种水华模拟参数个数将成指数增加［8］。

由于物理模型的复杂性和局限性，单纯基于数据挖掘算法和统计学理论，建立水华与驱动因素之间的关系的数据统计模型可避开其中的复杂环节，得以大量应用。孔繁翔等［9］利用水文和气象监测，结合卫星遥感技术对太湖蓝藻水华发生进行预测，准确率达80%。人工神经网络算法由于方法简单、无需建立复杂的数学计算模型、具有很强的适应性和容错性，在水华预报中得到广泛的应用。Shen 等［10］利用支持向量机较好地模拟和预测年际间藻华变化。覃苗等［11］用水温、溶解氧、电导率和气温作为输入变量，采用BP神经网络模型（Back Propagation）对东张水库叶绿素a 进行预测，预测结果较好，拟合度达0.83。神经网络虽以成功应用于水华预警中，但该模型不易收敛，需要反复多次调参才能趋于稳定。因此不仅考虑前一时刻的输入，而且赋予了网络对前期内容的“记忆”功能的改进后的RNN循环神经网络（Rerrent Neural Network）得以逐步推广，而长短期记忆网络LSTM（Long Short-Term Memory）则是在RNN 循环神经网络的基础上解决其存在的数据长期依赖问题，适合于处理和预测时间序列中间隔和延迟非常长的时间，在水质预测中进行了一定的应用［12］。

综上，本研究通过近两年对棉花滩水库的水文水生态跟踪监测，在分析浮游植物群落结构及生境变化的基础上，采用LSTM进行水华预警建模，以期为水华风险预警提供技术支撑。

1 研究区域与方法

1.1 研究区域

棉花滩水库大坝位于福建省永定区境内的汀江干流棉花滩峡谷河段中部，由汀江干流和黄潭河中下游河段组成，是以发电、防洪为主，航运、水产养殖为辅的大型调节型水库［13］。大坝工程于1998年开工，2002年12月全部完工，控制流域面积为7 907 km2，总库容20.35 m3［14］。库区属亚热带海洋性季风气候，多年平均气温为20.1 ℃。受温润暖湿气流影响，库区降水充沛，干湿季节分明，雨季相对集中，多年平均降水量1 400～1 800 mm，最大年平均降雨量近2 000 mm［15］。

1.2 监测方案

本次调查集中于棉花滩大坝库首区，样点布设见图1，其中汀江干流：上游源头自然河流段8 号、变动回水区7 号、回水淹没区6 号、5 号、交汇以后库区9 号、10 号、11 号、12 号。支流黄潭河：上游源头自然河流段4号、变动回水区3号、回水淹没区2号，1 号。2019年7-12月、2020年5-10月开展生态环境监测，共计34次。

图1 棉花滩水库采样点分布图Fig.1 Distribution of sampling sites in the Mianhuatan Reservoir

1.3 监测指标及方法

营养盐：取表层以下0.5 m 水体500 mL，用于测定水体中总氮（TN）、氨氮（NH4-N）、硝氮（NO3-N）、总磷（TP）、正磷酸盐（PO4-P）、叶绿素a（chlorophyll-a，简记为chl-a），现场浓硫酸低温保存，48 h 内进行水化学分析实验，分析测试方法依据《水和废水监测分析方法》［16］。

水环境参数：温度、溶解氧、电导率、pH、氧化还原电位，由水质多参数仪（Hydrolab DS5）现场测定表层水体获得。

藻类鉴定：取监测点表层水体1 L，每瓶加鲁哥试剂2.5 mL，低温保存，每静置8 h 以上采用医用输移管在微扰动模式下移除上层清液，反复几次，逐步沉淀浓缩至50 mL 后，在迅数藻类鉴定计数仪下对浮游藻类细胞进行分类镜检到属，其藻种鉴定参照《中国淡水藻类》［17］。

2 风险预报模型原理

2.1 模型介绍

长短期记忆（LSTM）网络最早于1997年被提出，它是递归神经网络（RNN）的改进版本［18］。LSTM 是一种时间递归神经网络，它继承了大多数RNN 模型的特点，解决了梯度反向传播过程中梯度消失的问题［19］。在RNN 的基础上，LSTM 增加了一个记忆“细胞”结构来判断信息是否有用。每个单元由一个输入门、一个忘记门和一个输出门组成。每条信息都进入了LSTM网络，并根据规则被认为是有用的，只留下符合规则的信息，不符合的信息通过遗忘门被遗忘，这对于具有长期序列依赖性问题的数据非常有效。

2.2 模型构建

图2 为LSTM 神经网络内部结构示意图，包括输入层、隐藏层和输出层［12］。

图2 LSTM神经网络结构图Fig.2 Neural network structure of LSTM

输入层：整个模型采用全连接结构，变量首先通过输入阀门进入神经网络结构，该阀门层的功能是决定更新的信息内容，包括sigmoid 函数层和tanh 函数层两部分，其层数与输入变量数量一致。其计算公式为：

隐藏层：在该模型隐藏层中设置50 个神经元，每个神经元细胞包含遗忘门和更新门。遗忘阀门主要用以决策，可自主决定舍去细胞状态中的信息。ht-1为上一层的输出结果，xt为现阶段的输入数据，输出结果的数值范围为［0，1］，Ct-1为各细胞的状态，其内部数字中1 意味着“完全保留”，0 意味着“完全舍弃”。更新阀门的功能为更新换代旧细胞和新细胞的状态，由状态Ct-1更新至状态Ct，该层的计算公式为：

输出层：该层功能主要为决策该层输出的类型和内容。这个输出过程为首先运行输入层中的sigmoid 层确定输出的哪些细胞状态，其次通过输入层中tanh函数层进行处理细胞状态，紧接着将其与sigmoid层的输出数据相乘，从而最终确定出输出内容。该层计算公式为：

本研究中的LSTM 神经网络模型使用Python 编程语言和其深度学习框架Keras 搭建。Keras 是一个用于快速构建深度学习网络的高级库，它的Sequential 模型可以快速堆叠不同类型的网络结构。由于输入输出参数量纲不一致，在进行运算之前需要进行数据标准化处理，计算公式如式（7）：

经过标准化处理的数据输入到LSTM 模型中。该神经网络模型中输入层设置一个单元，隐藏层设置50 个LSTM 神经元，输出层为一个单元，采用Dense（全连接）结构，该模型其余参数：迭代次数epochs 取值为1 000、分段长度batchsize 取值72、损失函数loss设置为mae，优化函数optimizer设置为adam。

2.3 误差评价指标

为分析各种参数预测性能的优劣，采用相对误差RE、平均绝对误差MAE和拟合度R2作为评价指标，计算方法如下：

式中：i表示样本脚标；n表示样本数量表示模拟值；yi表示实测值表示实测值的均值。

3 结果与分析

3.1 环境因子变化特征

表1 列举了棉花滩库区各环境因子平均值、变化范围以及变异系数。由表可知，其中平均水温27.63 ℃，平均气温29.35 ℃，波动幅度相对较小。由氧化还原电位可知，水体整体呈现出还原性，均值为-85.06，变异程度为所有参数最大。chla浓度变化范围为0.28～76.51 mg/m3，均值为16.16 mg/m3，离散程度仅次于氧化还原电位。溶解氧DO 含量变化范围位于5.04～10.74 mg/L，平均为7.70 mg/L，整体维持在较高水平。氨氮NH4-N含量位于0.34～0.90 mg/L之间，平均含量为0.06 mg/L，总磷TP 平均含量为0.02 mg/L，正磷酸盐PO4-P 含量偏低，总体来说，水质情况良好，处于地表水二类水质。

表1 棉花滩环境因子特征Tab.1 Characteristics of environmental factors in Mianhuatan Reservoir

3.2 浮游植物群落结构演替规律

棉花滩水库2019年、2020年共鉴定藻门6 门63 属（如表2），分别为绿藻门（21 属）、甲藻门（6 属）、硅藻门（19 属）、隐藻门（1 属）、蓝藻门（15 属）、黄藻门（1 属）。其中绿藻门属数最多，占浮游植物总属数的33.33%。小环藻为所有藻种中出现频率最高藻种，达到93.69%，为棉花滩水库常见优势藻种，其次为针杆藻、小球藻、衣藻。对比陈丽萍等在2010年9月紫金铜矿污染污染事件后棉花滩库区浮游植物群落结构构成为6 门35属［20］。

表2 棉花滩水库不同浮游植物藻属构成Tab.2 Composition of phytoplankton algae in Mianhuatan Reservoir

通过对棉花滩水库的浮游植物总藻密度可知，2019年7月总藻浓度达到最高值，达到0.35 亿个/L，2019年8月浓度大幅降低，9月以后藻类密度非常低。2020年6月总藻浓度达到最高值，达到0.286 亿个/L，7月浓度大幅降低，此后总藻密度一直呈上下波动状态。

图3 棉花滩水库浮游植物种类百分比Fig.3 Percentage of phytoplankton species in Mianhuatan Reservoir

3.3 风险预报模型结果分析

3.3.1 相关性分析

由于棉花滩水库共鉴定浮游植物藻种6 门63 属，且不同藻种生长关联要素众多，成因复杂，为开展水华风险预警，本文选择浮游植物总藻密度为预警目标，通过相关性分析，遴选棉花滩库区影响水华生长的关键因子，筛除关联性较弱因子，为后续LSTM 模型输入条件提供依据。此处主要采用棉花滩水库藻类生长的水质、气象因子与总藻密度开展Pearson 相关性分析，结果如表3 示。由表可知，总藻密度与水质因子中的水温、pH、氧化还原电位、溶解氧DO、五日生化需氧量、总磷、硝酸盐氮、高锰酸盐指数等水质参数存在显著相关关系（P＜0.01）。其中，藻密度与溶解氧的相关系数达0.57，符合浮游植物与溶解氧浓度成正相关关系的一般性规律［21］。水质因子中，营养盐（总磷、硝酸盐）对预测结果不敏感，说明主要影响棉花滩水库藻类繁殖的是其环境因子，营养盐不是藻类生长的限制因子。

表3 藻密度与各指标Pearson相关关系Tab.3 Correlation between algal density（units/L）and Pearson of each index

藻密度与水文气象条件中的气温（日均气温、最高气温、最低气温）、平均风速、风向、极大风速、入库及出库流量等存在显著相关性（P＜0.01）。

3.3.2 预测结果分析

依据相关性分析结果，选取显著相关（P＜0.01）的影响因素作为输入变量，分别以水质单因子、水文气象单因子和组合因子3种情况作为LSTM模型输入变量去预测藻密度浓度。

表4列举了气象因子中和水质因子逐一代入模型中后输出的模拟结果，结果表明，气象水文因子中入库流量模拟结果表现最好，实测值与预测值拟合度R2为0.68，出库流量次之，拟合度R2为0.60。水质因子中溶解氧为单一输入变量时的模拟结果时，实测值与预测值拟合度R2最好，可达0.65，平均绝对误差MAE为0.041 亿个/L，这与水华暴发时浮游植物增殖所产生的光合作用制氧关系密切［22］。其次为pH，拟合度R2为0.50，这与水华暴发时浮游植物增殖光合作用利用水中的CO2合成有机质，使得藻类的生物量增加，pH 升高有关［23］。水质因子中的溶解氧和pH 异常波动时水华暴发的结果，从一定程度上表明了预测结果的合理性。目前的研究中，有部分研究成果即以溶解氧DO浓度、pH来开展水华的预测预报工作［24，25］。

表4 单一因子作为输入变量模拟结果Tab.4 Simulation results of singal factor as input variables

图4 棉花滩水库总藻密度变化图Fig.4 Variation of total algae density in Mianhuatan Reservoir

图5（a）为仅以入库流量作为单一输入变量预测结果，由图可知，棉花滩水库藻密度预测与实测结果整体拟合度较好，特别是藻密度浓度峰值拟合效果较好，平均绝对误差MAE为0.037 亿个/L，这可能与棉花滩水库规律性季节调节方式有关。图5（b）为以溶解氧为输入变量时的模拟结果，虽平均误差整体相对其他水质因子较小，但是对藻类峰值时的模拟结果较差。

图5 分别以入库流量和溶解氧作为输入变量预测结果Fig.5 Forecast results of inflow flow and dissolved oxygen as input variables

当该模型输入变量为单因子时，拟合度均不超过0.70，模拟精度整体不高。为进一步提高模拟进度，此处对与水华存在相关性的16个因素分别进行了二因子组合，三因子组合和四因子组合，由于篇幅所限，此处仅列举预测模拟结果较好的情况（表6），由表可知，二因子“水温+入库流量”组合模式拟合度优于“水温+溶解氧”的二因子组合，R2可达到0.73。当增加为3因子时，即“水温+溶解氧+入库流量”的模式下，拟合度R2可进一步提高至0.75，当增加至四因子时，多因素组合下R2均可达到0.76，其中以“日均气温、水温、风向、入库流量”4 个指标作为模型输入变量，平均绝对误差MAE最低为0.031 亿个/L，同时由于这4 个因子可以结合天气预报提前获知，较容易实现水华的预测预报工作。

表5 组合因子作为输入变量模拟结果Tab.5 Simulation results of combining factors as input variables

在不同因子组合遴选过程中可知，棉花滩气温和出入库流量对藻类生长有较大的影响，气温升高会导致水库表层水温升高，营造利于藻类的生境条件，这与本人曾获取的淡水水库10大优势藻种的温度与藻类的本构关系曲线也可得到验证［8］。本研究中，预测结果对入库流量最为敏感，入库流量的大小可能会影响水库的水动力过程，改变水库垂向水温分层结构，影响藻类繁殖，其机制还有待进一步探究［26，27］。

对以上四因子组合（日均气温、水温、风向、入库流量）作为输入变量时，对实测值与预测值输出结果开展误差分析。从表6 中可以看出，当藻密度量级在2 500 万个/L 以下时，相对误差会出现较大值，此时预测值可信度较低，但此时为水华低风险时期，预测值不具有参考性；随着水华浓度的升高，实测值和预测值误差缩小，当实测藻密度高于107个/L 时，相对误差RE变化范围为0.02～0.73，整体处于可接受范围，能大体预测出水华风险趋势。

表6 实测值与预测值对比Tab.6 Comparison between Measured Values and Predicted Values

图6 为四因素组合下实测与预测结果对比图，由图6 可知，该模型整体拟合结果较好，虽然在藻密度浓度较低时，存在一定误差，但随着藻浓度升高，精度在可接受范围，加上这四个因子可以结合天气预报提前获知，较容易实现水华的预测预报工作，推荐作为棉花滩库区首选水华预测组合。

图6 以气温、水温、风向、入库流量作为输入变量预测结果Fig.6 Forecast results of air temperature，water temperature，wind direction and inflow flow as input variables

综合以上预测结果可知，气象水文单因子中，入库流量凭借其季节规律性和对降水的响应，可作为LSTM 模型预测棉花滩库区首选，水质单因子中，溶解氧可作为水华预测首选，但对藻密度峰值模拟精度稍差。多因素组合下，首选日均气温、水温、风向、入库流量四个指标作为模型输入变量，可较好的预测该区域水华。

4 结 论

通过对棉花滩水库近两年水生态环境连续监测，本文系统分析了该区域浮游植物群落结构特征，并通过相关性分析对影响水华关键因子进行筛选。依据相关性分析结果，选取显著相关（P＜0.01）的影响因素作为输入变量，分别以水质单因子、水文气象单因子和组合因子3 种情况作为LSTM 模型输入变量去预测总藻密度，主要结论如下。

（1）棉花滩水库属亚热带海洋性季风气候，气候温润多雨，整体水质处于地表水II 类水体，水质状况良好，浮游植物样品中共发现6 门63 属，主要优势藻种为小环藻、针杆藻、小球藻、衣藻，从属于绿藻和硅藻门，蓝藻水华整体风险较低；

（2）单因子中分别以气象水文因子中入库流量和水质因子中的溶解氧作为单一输入变量时的模拟结果最佳；多因子下以日均气温、水温、风向、入库流量的组合条件下，利用LSTM 模型预报浮游植物总藻种密度效果最佳，拟合度可达到0.76，特别当实测藻密度高于0.1 亿个/L 时，相对误差RE变化范围为0.02～0.73，推荐为水华风险预警优选组合。 □