漢江中下游干流水華關(guān)鍵環(huán)境因子識(shí)別及閾值分析

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

(1.武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072；2.長(zhǎng)江水利委員會(huì)水文局長(zhǎng)江中游水文水資源勘測(cè)局，湖北 武漢 430010)

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

以往的研究認(rèn)為，漢江水華的發(fā)生是水文、氣象和營養(yǎng)鹽條件多方面綜合的結(jié)果[5]。國內(nèi)學(xué)者圍繞漢江水華已開展了大量的研究工作，謝平等[6]指出制約漢江水華發(fā)生的關(guān)鍵因子是流量和流速等水文因子，并非氮、磷等水質(zhì)因子和水溫等氣象因子；吳興華等[7]基于2015—2016年漢江硅藻水華的監(jiān)測(cè)數(shù)據(jù)，得出水華發(fā)生的成因是適宜的氣候條件、較高的硅氮比和低流量；王俊等[8]基于2018年漢江水華的應(yīng)急監(jiān)測(cè)數(shù)據(jù)，提出了嚴(yán)控污染輸入、優(yōu)化水量調(diào)度方案、完善管理機(jī)制等治理對(duì)策。總體來看，研究者基于不同時(shí)段、不同斷面的監(jiān)測(cè)數(shù)據(jù)得到的漢江水華的研究結(jié)論不盡相同。當(dāng)前，漢江中下游水華的預(yù)警和防治仍面臨著一些難點(diǎn)和挑戰(zhàn)[2,9]：①漢江中下游河段目前還未建立常規(guī)的水華監(jiān)測(cè)體系，應(yīng)急監(jiān)測(cè)期間的資料長(zhǎng)度較短且受到水量調(diào)度的影響，漢江水華的暴發(fā)成因，特別是預(yù)警研究尚未取得突破性進(jìn)展；②在氣候變化和人類活動(dòng)(如大型水利工程)的影響下，漢江中下游的環(huán)境因素復(fù)雜多變，以往研究多針對(duì)河流水文情勢(shì)變化對(duì)浮游植物生長(zhǎng)的單項(xiàng)研究，變化環(huán)境影響下的河流水華暴發(fā)與預(yù)測(cè)尚未取得好的效果；③某些重要環(huán)境因子一旦超過閾值，藻類在浮游植物群落的占比可能發(fā)生劇烈變化。因此，有必要探索漢江硅藻暴發(fā)的重要環(huán)境因子閾值，為治理漢江富營養(yǎng)化提供技術(shù)支撐。

本文基于漢江中下游干流硅藻水華發(fā)生前后的水文、水質(zhì)和水生態(tài)監(jiān)測(cè)數(shù)據(jù)，采用Pearson相關(guān)性檢驗(yàn)、分位數(shù)回歸和隨機(jī)森林模型方法分析水華發(fā)生的關(guān)鍵環(huán)境因子及閾值，以期為漢江中下游生態(tài)調(diào)度提供參考。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

研究區(qū)域?yàn)闈h江中下游皇莊至宗關(guān)水廠江段，以仙桃站為基本分析斷面。2021年1月18日，仙桃斷面疑似發(fā)生水華，長(zhǎng)江委水文局啟動(dòng)應(yīng)急響應(yīng)，于1月19日至2月2日開展應(yīng)急監(jiān)測(cè)和調(diào)度，分別在漢江中下游的皇莊、沙洋、興隆壩上、澤口、仙桃和宗關(guān)水廠6個(gè)漢江干流控制斷面開展跟蹤取樣監(jiān)測(cè)，斷面位置如圖1所示。

圖1 漢江中下游干流采樣斷面分布

1.2 數(shù)據(jù)來源

采用的水華監(jiān)測(cè)數(shù)據(jù)均來自長(zhǎng)江委水文局《漢江“水華”連續(xù)應(yīng)急監(jiān)測(cè)調(diào)查報(bào)告》及歷史水文資料。皇莊、沙洋、興隆壩上、澤口、仙桃和宗關(guān)水廠6個(gè)漢江干流控制斷面涉及的關(guān)鍵要素包括環(huán)境指標(biāo)(氣象、水質(zhì)和水文因子)和浮游植物(藻密度、葉綠素a質(zhì)量濃度)，其中，氣象數(shù)據(jù)來自國家氣象信息中心中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn)，其余數(shù)據(jù)均來自長(zhǎng)江委水文局。

2 研究方法

2.1 分位數(shù)回歸

浮游植物的生長(zhǎng)周期短，在環(huán)境因子變化較為劇烈時(shí)，其生物量會(huì)出現(xiàn)劇烈的時(shí)空變化，數(shù)據(jù)具有較強(qiáng)的時(shí)空異質(zhì)性，不能很好地滿足傳統(tǒng)最小二乘法中正態(tài)、等方差以及無自相關(guān)性等假設(shè)條件[10-11]。而分位數(shù)回歸方法適用的條件較為廣泛，可以很好地處理數(shù)據(jù)中的方差異質(zhì)性和異常值[12-13]。該方法由Koenker等[14]提出，可以根據(jù)響應(yīng)變量的條件分位數(shù)對(duì)預(yù)測(cè)變量進(jìn)行回歸，進(jìn)而得到所有響應(yīng)變量分位水平上的回歸模型。假設(shè)Y是一個(gè)連續(xù)型隨機(jī)變量，F(xiàn)Y(y)是其分布函數(shù)，一個(gè)特定的y值的發(fā)生概率為τ，那么第τ分位數(shù)就是指發(fā)生概率為τ時(shí)的y值，即：

(1)

如果p(Yyp)≤1-p，那么yp被稱為隨機(jī)變量Y的第p分位數(shù)。與普通最小二乘法把預(yù)測(cè)變量的期望平均值作為響應(yīng)變量的函數(shù)不同，分位數(shù)回歸可模擬一個(gè)或更多響應(yīng)變量的分位數(shù)[15]，分位數(shù)回歸的特征可以度量預(yù)測(cè)變量對(duì)響應(yīng)變量分布中心的影響, 還可以度量預(yù)測(cè)變量對(duì)響應(yīng)變量整體分布的影響。

2.2 隨機(jī)森林模型

由于河流水華的成因復(fù)雜，各影響因子之間的相互作用機(jī)制尚不明確[9]，因此大多數(shù)研究是基于黑箱模型[16-18]，采用數(shù)理統(tǒng)計(jì)方法或智能算法。在智能算法中，隨機(jī)森林模型是一個(gè)包含多個(gè)決策樹的分類器，并且其輸出的類別是由個(gè)別樹輸出的類別的眾數(shù)而定[19]。隨機(jī)森林模型可以產(chǎn)生高準(zhǔn)確度的分類器，在決定類別時(shí)可以評(píng)估變數(shù)的重要性，并且學(xué)習(xí)過程快速[20]。

選取偏差系數(shù)作為葉綠素a的擬合優(yōu)度指標(biāo)。偏差系數(shù)是對(duì)實(shí)際值與模型模擬值之間離差的一種度量，其計(jì)算公式為

(2)

式中：B為偏差系數(shù)；yi為因變量的實(shí)際值；y′i為因變量的模擬值。

3 結(jié)果與分析

3.1 監(jiān)測(cè)結(jié)果

3.1.1藻密度

2021年水華監(jiān)測(cè)期間，漢江干流不同監(jiān)測(cè)斷面的藻密度如圖2所示(1月30日無監(jiān)測(cè)數(shù)據(jù))。由圖2可知，從皇莊到宗關(guān)水廠斷面，不同監(jiān)測(cè)斷面之間具有顯著的空間差異性，其中宗關(guān)水廠的藻密度最大，仙桃斷面次之，皇莊斷面最小，表明漢江中下游的藻密度呈沿程增加的趨勢(shì)。隨著時(shí)間的推移，藻密度持續(xù)下降，1月24日各斷面間的藻密度差異最大。

圖2 2021年不同監(jiān)測(cè)斷面不同時(shí)間的藻密度

根據(jù)DB 44/T 2261—2020《水華程度分級(jí)與監(jiān)測(cè)技術(shù)規(guī)程》，各個(gè)斷面在不同時(shí)間的水華等級(jí)如表1所示。由此可見，1月23日水華程度最為嚴(yán)重，所有監(jiān)測(cè)斷面的水華等級(jí)均達(dá)到Ⅲ級(jí)輕度水華的標(biāo)準(zhǔn)；1月24日后，皇莊、沙洋和興隆壩上斷面的水華等級(jí)均降到Ⅱ級(jí)及以下；仙桃和宗關(guān)水廠斷面分別在1月28日和1月27日后保持Ⅱ級(jí)，無明顯水華。

表1 各監(jiān)測(cè)斷面的水華程度分級(jí)

3.1.2理化指標(biāo)

2021年漢江中下游干流水華發(fā)生期間，各監(jiān)測(cè)斷面的理化指標(biāo)均值和最大值如表2所示。可見，干流透明度均值基本呈遞減趨勢(shì)，在皇莊斷面最高，為90.25 cm，宗關(guān)水廠斷面最低，為43.60 cm。水溫的最大值在沙洋斷面，達(dá)到12.7 ℃，水溫均值在仙桃斷面最低，為9.08 ℃，在沙洋斷面最高，為10.77 ℃。溶解氧的質(zhì)量濃度均值和最大值均在興隆水庫壩上斷面最高, 仙桃斷面均值為12.65 mg/L，皇莊斷面均值最低，為11.35 mg/L，這是由于當(dāng)水體中充滿高密度的藻類細(xì)胞時(shí)，光合放氧會(huì)導(dǎo)致水中的溶解氧質(zhì)量濃度過高[7]。pH均值范圍為8.18～8.61，呈現(xiàn)沿程遞增趨勢(shì)，在宗關(guān)水廠斷面達(dá)到最大值。葉綠素a質(zhì)量濃度均值也呈沿程遞增趨勢(shì)，在宗關(guān)水廠斷面最高，為70.02 μg/L，皇莊斷面最低，為23.10 μg/L。高錳酸鹽指數(shù)的質(zhì)量濃度為2.67～4.20 mg/L。所有斷面的總磷質(zhì)量濃度均超過0.03 mg/L，總氮質(zhì)量濃度均值在1.5 mg/L左右，可見漢江下游河段較高的氮磷質(zhì)量濃度為水華的發(fā)生提供了充足的營養(yǎng)鹽條件。

表2 水華發(fā)生期間各監(jiān)測(cè)斷面的理化指標(biāo)均值和最大值

3.2 Pearson相關(guān)性

以Pearson相關(guān)性系數(shù)進(jìn)行雙側(cè)檢驗(yàn)，分別分析藻密度與監(jiān)測(cè)的氣象、水質(zhì)、水文因子的相關(guān)性。由于藻密度僅在水華發(fā)生后進(jìn)行監(jiān)測(cè)，無法體現(xiàn)水華的生消過程，因此選取葉綠素a來表征浮游植物的生物量。由于仙桃斷面水華暴發(fā)的頻率較高且該斷面設(shè)置了自動(dòng)監(jiān)測(cè)站，因此重點(diǎn)對(duì)該站進(jìn)行分析。圖3為仙桃站葉綠素a質(zhì)量濃度與藻密度的相關(guān)關(guān)系，可見，葉綠素a質(zhì)量濃度與藻密度呈顯著線性正相關(guān)(R2=0.797)，因此可以用葉綠素a質(zhì)量濃度來表征浮游植物的生物量。基于水華程度分級(jí)標(biāo)準(zhǔn)，當(dāng)藻密度小于1×107個(gè)/L時(shí)無明顯水華，該藻密度值為水華的臨界值。因此依據(jù)藻密度和葉綠素a質(zhì)量濃度的擬合公式，葉綠素a質(zhì)量濃度的臨界值為76.47 μg/L，可以作為水華發(fā)生的一項(xiàng)判別標(biāo)準(zhǔn)。

圖3 仙桃站葉綠素a質(zhì)量濃度與藻類密度的相關(guān)關(guān)系

采用仙桃站2021年1月1日至31日的數(shù)據(jù)，分別判斷葉綠素a與氣象、水文、水質(zhì)因子的相關(guān)性。氣象因子與葉綠素a的Pearson相關(guān)性分析結(jié)果表明，葉綠素a與平均氣溫在0.01的顯著性水平上呈正相關(guān)關(guān)系(R2=0.639)，與平均風(fēng)速的相關(guān)關(guān)系(R2=-0.104)和日照時(shí)數(shù)的相關(guān)關(guān)系(R2=-0.279)均不顯著。

葉綠素a與水位和流量均呈顯著負(fù)相關(guān)關(guān)系，R2分別為-0.723和-0.592，與以往的研究結(jié)果一致，即漢江水華的發(fā)生與枯水期低流量有關(guān)[7]。圖4為葉綠素a質(zhì)量濃度與流量的關(guān)系，1月20日前，流量呈下降趨勢(shì)且整體偏低，此時(shí)葉綠素a的質(zhì)量濃度逐漸上升并在20日達(dá)到最大值，之后葉綠素a質(zhì)量濃度持續(xù)下降。這是由于在本次漢江水華發(fā)生后，自1月19日開始興隆水庫按“沖蓄結(jié)合”方式進(jìn)行應(yīng)急調(diào)度。興隆樞紐以上河段的水華影響通過丹江口、王甫洲等梯級(jí)水庫的應(yīng)急調(diào)度解決，1月24日開始丹江口水庫日均下泄流量由620 m3/s加大至800 m3/s，并持續(xù)6 d。興隆樞紐以下漢江河段的水華影響在維持興隆水庫出入庫平衡的基礎(chǔ)上，通過引江濟(jì)漢工程應(yīng)急調(diào)度解決。

圖4 仙桃站葉綠素a質(zhì)量濃度與流量的關(guān)系

葉綠素a與水質(zhì)理化指標(biāo)的Pearson相關(guān)性檢驗(yàn)表明，葉綠素a與水溫、pH和溶解氧均有顯著的正相關(guān)性(P<0.01)，相關(guān)性系數(shù)分別為0.388、0.526和0.407。表3為1月21日漢江中下游斷面總氮、總磷質(zhì)量濃度變化情況，可以看出：①支流氮、磷營養(yǎng)鹽水平明顯高于干流，表明支流匯入帶來的營養(yǎng)鹽負(fù)荷可能對(duì)漢江干流水質(zhì)造成一定影響；②干流斷面的總氮、總磷質(zhì)量濃度雖與非水華期無明顯差異，但總氮、總磷的質(zhì)量濃度較高，總氮質(zhì)量濃度普遍超過1.3 mg/L，總磷質(zhì)量濃度普遍超過0.03 mg/L，說明此次漢江的硅藻水華是多種環(huán)境因子綜合作用的結(jié)果，即適宜的氣象和水文條件加上較高的營養(yǎng)鹽所致。

表3 1月21日漢江中下游斷面總氮、總磷質(zhì)量濃度變化

3.3 分位數(shù)回歸結(jié)果

水華發(fā)生會(huì)影響水體的pH值和溶解氧質(zhì)量濃度，藻密度增加后，水中CO2被利用而減少，導(dǎo)致pH值上升，同時(shí)導(dǎo)致溶解氧增加。因此，pH值和溶解氧因子作為水華發(fā)生后的結(jié)果，不適宜用來分析其對(duì)藻類的影響。圖5為仙桃站氣溫和流量的參數(shù)估計(jì)結(jié)果隨著葉綠素a豐度分位數(shù)的變化，包括分位數(shù)回歸斜率及其95%的置信區(qū)間，普通最小二乘法得出的斜率及其95%的置信區(qū)間。由圖5可見，氣溫在葉綠素a的所有豐度上均具有顯著的正斜率；流量在葉綠素a的所有豐度上均具有顯著的負(fù)斜率，各個(gè)因素對(duì)葉綠素a的影響作用都是非線性的。以流量因子為例對(duì)結(jié)果進(jìn)行說明，最小二乘法得到的斜率始終在-0.2附近，表明是一種均值回歸；而分位數(shù)回歸結(jié)果(圖6(b))表明，在0.05、0.25、0.50、0.75豐度分位數(shù)下，直線斜率分別為0.03、-0.3、-0.18和-0.16，可見在不同的葉綠素a分布條件下，流量對(duì)葉綠素a的影響不同(圖6)。

(a)氣溫

(a)氣溫

3.4 隨機(jī)森林模型模擬結(jié)果

基于實(shí)測(cè)氣溫和流量數(shù)據(jù)，通過隨機(jī)森林模型預(yù)測(cè)葉綠素a質(zhì)量濃度，圖7為葉綠素a質(zhì)量濃度的實(shí)測(cè)值與模擬值對(duì)比。由圖7可見，模擬值的相對(duì)誤差為-3.13%，說明隨機(jī)森林結(jié)果具有可信度。從影響因子的重要性來看，流量和氣溫因子對(duì)葉綠素a質(zhì)量濃度影響的重要程度分別為55.8%和44.2%。

圖7 葉綠素a質(zhì)量濃度的實(shí)測(cè)值與模擬值對(duì)比

3.5 基于水華防控的環(huán)境閾值分析

確定水華開始出現(xiàn)臨界藻密度作為臨界預(yù)警時(shí)刻。識(shí)別臨界藻密度對(duì)應(yīng)的水質(zhì)、氣象和水文因子，作為控制水華的環(huán)境因子閾值。

前文分析可知，氣溫是顯著影響葉綠素a含量的因子之一。為了探究氣溫的閾值，統(tǒng)計(jì)2021年1月3 d、5 d、7 d、9 d、11 d滑動(dòng)積溫，分別計(jì)算其與葉綠素a質(zhì)量濃度的相關(guān)關(guān)系，結(jié)果見圖8(a)，發(fā)現(xiàn)葉綠素a質(zhì)量濃度與7 d滑動(dòng)積溫的相關(guān)性最大，R2達(dá)到0.909。擬合葉綠素a質(zhì)量濃度與7 d滑動(dòng)積溫得到方程如圖8(b)所示，7 d滑動(dòng)積溫的閾值為56.02 ℃，即當(dāng)7 d滑動(dòng)積溫超過56.02 ℃時(shí)，將是漢江中下游硅藻水華暴發(fā)的危險(xiǎn)時(shí)間。

(a)葉綠素a質(zhì)量濃度與滑動(dòng)積溫的相關(guān)關(guān)系

圖9(a)為仙桃站1992—2021年1—3月的平均流量過程與所有水華發(fā)生期間的藻密度關(guān)系。可見，大多數(shù)的水華事件發(fā)生在平水年和枯水年，并且對(duì)應(yīng)1—3月的平均流量頻率均大于37%。圖9(b)為歷年1—3月的最小7 d平均流量，2008年以前最小7 d平均流量低于550 m3/s的年份都發(fā)生了水華，但2008年后發(fā)生水華的最小7 d平均流量有了明顯上升。有水華年份的最小7 d平均流量明顯低于無水華年份，當(dāng)最小7 d平均流量高于800 m3/s時(shí)，從未發(fā)生過藻類水華，當(dāng)其高于700 m3/s時(shí)，有時(shí)會(huì)發(fā)生水華(2009年)。仙桃站1992—2021年1—3月的最小7 d平均流量的范圍為364～1 164 m3/s。以50 m3/s為步長(zhǎng)，統(tǒng)計(jì)在不同流量區(qū)間內(nèi)發(fā)生水華的次數(shù)，據(jù)此得到不同流量下發(fā)生水華的概率，見圖9(c)。可見，流量與不發(fā)生水華概率的擬合相關(guān)性較好，R2達(dá)到0.922 3。通過擬合關(guān)系推求得到了不同最小7 d平均流量下不發(fā)生水華的概率，基于1992—2021年的樣本數(shù)據(jù)，仙桃站在60%、70%、80%、90%不發(fā)生水華的概率下，對(duì)應(yīng)的最小7 d平均流量應(yīng)分別高于310 m3/s、470 m3/s、610 m3/s、780 m3/s。

(a)歷年1—3月平均流量與藻密度關(guān)系

4 結(jié) 論

a.葉綠素a質(zhì)量濃度分別與氣溫、流量因子有顯著的正相關(guān)性和負(fù)相關(guān)性，但各個(gè)因子對(duì)葉綠素a質(zhì)量濃度的影響都是非線性的。

b.基于篩選出的環(huán)境因子，葉綠素a質(zhì)量濃度模擬值與實(shí)測(cè)值的相對(duì)誤差值為-3.13%，具有較好的擬合優(yōu)度。流量和氣溫因子對(duì)葉綠素a影響的重要程度分別為55.8%和44.2%。

c.基于仙桃站水華防控的環(huán)境閾值分析表明，硅藻水華暴發(fā)的水質(zhì)、氣溫和流量臨界預(yù)警時(shí)刻為：葉綠素a質(zhì)量濃度高于76.47 μg/L，7d滑動(dòng)積溫超過56.02 ℃，最小7 d平均流量低于780 m3/s。