棉花灘水庫浮游植物群落結構演替特征及風險預警研究

崔玉潔，楊友德，鄭婉婷，成再強，陳天生，林曉芳

（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］。

由于物理模型的復雜性和局限性，單純基于數據挖掘算法和統計學理論，建立水華與驅動因素之間的關系的數據統計模型可避開其中的復雜環節，得以大量應用?？追毕璧龋?］利用水文和氣象監測，結合衛星遙感技術對太湖藍藻水華發生進行預測，準確率達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，推薦為水華風險預警優選組合。 □