漢江中下游干流水華關鍵環境因子識別及閾值分析

田 晶，郭生練，王 俊，張德兵，羅春艷

(1.武漢大學水資源與水電工程科學國家重點實驗室，湖北 武漢 430072；2.長江水利委員會水文局長江中游水文水資源勘測局，湖北 武漢 430010)

水華指淡水水體中藻類大量繁殖的一種自然生態現象，是水體富營養化的一種特征。與出現在湖泊、水庫等靜水生態系統中的藍藻、綠藻水華相比，河流水華并不常見，但卻具有影響范圍廣、暴發成因復雜、控制難度大等特點[1-2]。漢江中下游近年來水華暴發頻次增加，引起廣泛關注。自1992年漢江中下游首次暴發硅藻水華以來，水華的發生頻次持續上升，持續時間也有所延長，以往發生的河段多為興隆以下河段，并有向支流蔓延的趨勢[3]，興隆低水頭閘壩型水庫建成后，2018年在興隆庫區也出現了水華現象。水華發生時水體呈棕褐色，散發腥味，對漢江中下游生產生活產生一定的影響[4]。因此，探究影響漢江水華發生的關鍵因子，精準防控水華的暴發，是目前迫切需要解決的問題。

以往的研究認為，漢江水華的發生是水文、氣象和營養鹽條件多方面綜合的結果[5]。國內學者圍繞漢江水華已開展了大量的研究工作，謝平等[6]指出制約漢江水華發生的關鍵因子是流量和流速等水文因子，并非氮、磷等水質因子和水溫等氣象因子；吳興華等[7]基于2015—2016年漢江硅藻水華的監測數據，得出水華發生的成因是適宜的氣候條件、較高的硅氮比和低流量；王俊等[8]基于2018年漢江水華的應急監測數據，提出了嚴控污染輸入、優化水量調度方案、完善管理機制等治理對策?？傮w來看，研究者基于不同時段、不同斷面的監測數據得到的漢江水華的研究結論不盡相同。當前，漢江中下游水華的預警和防治仍面臨著一些難點和挑戰[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)?？梢姡髁颗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。