程兵芬 ,夏瑞,張遠*,張楠 ,張新飛 ,劉成建,
1.北京師范大學水科學研究院,北京 100875;2.中國環境科學研究院,北京 100012;3.西北大學城市與環境學院/陜西省地表系統與環境承載力重點實驗室,陜西 西安 710127
隨著全球水體富營養化問題的加劇,國內外多個大型河流出現了嚴峻水華問題,嚴重影響了當地河流的水體結構與水生態環境,對人類飲水安全造成了極大的威脅(Jung et al.,2011;竇明等,2002;張遠等,2017)。通常富營養化容易發生在水體流速相對靜止的湖泊和水庫,然而受多種因素的綜合影響,國內外諸多大型河流頻繁發生水華現象(Whitehead et al.,2015;梁開學等,2012;夏瑞等,2020)。營養鹽含量、水文條件和氣候條件的綜合作用是水華發生的主要因素(Carey et al.,2012);如營養鹽總氮和總磷作為藻類生長的必要條件(吳興華等,2017);流速與流量通過改變水體滯留時間影響藻類的生長,水位可以通過影響水體分層現象影響藻類密度的變化(Ji et al.,2017);降水的時空變化則會影響藻類的頻率與強度(Phlips et al.,2007);水溫則會影響不同藻類物種的耐受性,不同的溫度條件下,優勢藻類結構有很大差異等(Jung et al.,2009;Kim et al.,2014);如藍藻適宜溫度在30—35 ℃,綠藻適宜溫度均在20—35 ℃,藍綠藻適宜溫度均在25—35 ℃,硅藻在0—15 ℃都能夠的較好的生長,硅藻相比較藍綠藻更適合低溫條件(Chen et al.,2014)。綜合分析營養物質、溫度和水文條件的是探索河流藻華機制的必要條件(Kim et al.,2018)。在當前氣候變化和人類活動共同作用背景下,河流水華未來的演變趨勢及其主要影響因素也是當前亟需探索的科學問題。
漢江自 1992年第一次硅藻水華暴發以來,已前后暴發 10次不同程度的水華事件,嚴重時覆蓋河段達到河口以上500多千米,單次持續時間最長可達20天之久(梁開學等,2012;夏瑞等,2020)。針對近年漢江水華的連年爆發過程,較多國內外學者對水華的現狀及成因做了相關基礎研究,如梁開學等(2012)發現2000年以后,漢江中下游硅藻水華與1992—2000年相比,面積縮小,但頻率增加;景朝霞等(2019)指出2004—2014年漢江中下游非汛期富營養化嚴重,非汛期總磷和總氮濃度呈降低和增加趨勢;潘曉潔等(2014)發現漢江春季水華暴發期間,不同采樣點,硅藻種類占浮游植物比例在44.13%—51.29%之間;夏瑞等(2020)研究了識別了漢江河流型水華多影響要素,指出近50年來,漢江流域氣候逐漸呈現暖干化趨勢,且藻類密度的變化比營養鹽和水文因子更為敏感;吳興華等(2017)指出漢江硅藻水華發生的原因為適宜的氣候條件、高硅氮比和低流量;殷大聰等(2012)發現高營養鹽負荷、低流速、小流量、晴好的氣候條件是導致漢江水華主要因素。然而,受多種因素制約,這些研究缺乏在連續長時間序列尺度下,對人類和自然雙重影響下的河流水華變化特征及多影響要素綜合識別。
本研究以漢江武漢段為研究區域,基于STL趨勢分解法(Seasonal Trend Decomposition using LOSS)、MK突變檢驗法(Mann-Kendall)、冗余分析(Redundancy Analysis,RDA)等方法,主要分析了 2004—2017年 2—4月漢江武漢段藻類密度及其影響要素的突變規律,識別了藻類密度突變的主要影響因子,為漢江水華防治和預警預報提供技術支撐。
漢江發源于陜西省寧強縣潘冢山,流經陜西、湖北兩省,于武漢市進入長江,全長1577 km,落差1964 km,流域面積1.59×105km2。漢江流域屬北亞熱季風氣候區,氣候溫暖濕潤,夏季平均氣溫26 ℃,冬季平均氣溫在 0 ℃間,年平均降水量約為897 mm,降水集中在5—10月,多年平均徑流量為5.91×1010m3。漢江武漢段位于平原地區,全長50 km,是武漢市居民用水和工業用水的主要水源(殷大聰等,2017)。漢江武漢段徑流緩慢,既受到上游南水北調取水工程影響下的來水來沙條件額度制約,又受到長江頂托作用的影響,具有十分復雜的水文情勢,水質及水華問題關注度高(謝平等,2004)。
環境與生物監測數據來源于國家城市供水水質監測網武漢監測站提供的宗關、白鶴嘴與琴斷口3個站點的水質監測數據,主要監測指標為總氮(TN,mg·L?1)、總磷(TP,mg·L?1)、藻類密度(107cell·L?1),監測時段為 2004—2017 年 2—4 月每個自然月的上中下旬各10 d左右;這個時段漢江處于枯水期同時水華現象較為頻繁(梁開學等,2012;夏瑞等,2020)。在每個采樣點,采集1.0? L的表層水加魯哥氏液固定,靜置48 h后采用虹吸方法吸取上清液,留100 mL備用。藻類細胞的數量在24 ?h內分析測定,吸取0.1 mL滴入計數框,用視野法計數,每個樣計數3次,取平均值;采用分光光度法對葉綠素含量進行測定;采用鉬酸銨分光光度法和堿性鉀消解紫外分光光度法,對1? L水樣進行原位過濾并帶回實驗室進行營養成分分析,包括 TP和TN 等(Xia et al.,2019;國家環境保護總局,2002)。水文監測數據來源于長江水利委員會水文局(http://www.cjh.com.cn/),均為同時期仙桃水文站的流量(Q)、流速(v)、水位(WL1)以及水溫(t)。同時采用了漢口水文站的水位數據(WL2),仙桃水文站與漢口站水位差 ΔH。研究區域及觀測站點位置、分類見圖1。

圖1 研究區域及采樣點位分布Fig.1 Geographical location and distribution of monitoring station
研究聯合使用STL趨勢分解法、MK檢驗法、以及 RDA冗余分析法,旨在綜合識別與河流水華暴發密切關聯的藻類密度變化趨勢和突變點,辨析水華年和非水華年期間藻類密度主要驅動因子和演變規律,其中STL趨勢分解法主要識別藻類密度的時間周期變化趨勢及突變點,MK檢驗法則檢驗變化趨勢置信度水平,RDA冗余分析則識別影響藻類密度變化的主要驅動因子及貢獻。具體研究方法如下圖2。

圖2 河流水華拐點突變分析方法Fig.2 Mutation Test Method for the evolution and mutation attribution of water bloom
1.2.1 STL趨勢分解法
STL(A Seasonal and Trend Decomposition Using Loess)最早由Cleveland(1979)提出,它是以魯棒局部加權回歸作為平滑方法的時間序列分解方法。STL方法定義T(t)時刻原始數據Yt經STL分解為趨勢項Tt、季節項St和殘差項Rt:

式中,趨勢項Tt被認為是相對穩定的數據變化項,代表了低頻數據的變化趨勢,是對污染源和氣象條件的綜合反應;季節性成分St為高頻率變化的周期性穩定的數據擾動項,主要是由氣象條件的季節變化引起的;殘差Rt是不規則的隨機擾動引起的。STL方法最早應用于經濟學領域特別是對美國失業人口數據分析上,近年在環境科學、生態學、水質科學等學科中有著廣泛的應用(Silawan et al.,2008)。
1.2.2 MK突變檢驗法
Mann-Kendall(MK)統計檢驗方法由Mann and Kendall兩人在20世紀中葉發明并完善,它不需要樣本遵從一定的分布,且受異常值影響較小(Mann et al.,1945;Kendall et al.,1975)。MK 統計檢驗方法在降水、水文、氣候變化等時間趨勢研究上廣泛應用(Evans et al.,2000)。
設時間序列數據(x1,x2,…xn)是n個獨立的、隨機變量同分布的樣本,定義標準的正態統計變量S:

式中,sign()為符號函數。當Xi?Xj<0、=或>0 時,sign(Xi?Xj)為?1;當Xi?Xj=0 時,sign(Xi?Xj)為 0;當Xi?Xj>0 時,sign(Xi?Xj)為 1。

對于統計值Z來說,Z>0時,上升趨勢;Z<0時,下降趨勢。|Z|>1.28、|Z|>1.64、|Z|>2.32 分別表示時間變化趨勢通過了置信度α=90%、α=95%和α=99%的顯著性檢驗。
當采用非參數 Mann-Kendall法進行突變檢測時,檢驗統計量與上述Z有所不同。設有一時間序列如下:x1,x2, …,xn,構造秩序列ri,ri含義為xi>xj(1≤j≤i)的累計數。定義:

Sk均值E(Sk)與方差Var(Sk)的定義如下:

假定一組時間序列是隨機獨立的變量,則對這組時間序列的統計量UFk作如下定義:

式中,UFk為標準正態分布;給定顯著水平α,查詢正態分布表,得到臨界值Uα。通過統計序列UFk和反序列UBk交叉點對序列x變化趨勢進一步進行分析,同時明確突變時間,并指出突變區域。
1.2.3 RDA冗余分析法
20世紀50年代開始生態學家開始基于排序的方法研究植被群落的連續分布。冗余分析(Redundancy Analysis,RDA)是一種使用物種和環境因子組成數據的排序技術,是當前生態學領域較為常用的直接梯度變量約束排序分析方法(Ordination Analysis。Vonwehrden et al.,2009)。冗余分析最早由Vandenwollenberg A L發明,它的主要原理簡要如下:

式中,Z1代表獨立的內生變量n×r矩陣;Z2代表自變量或外生變量為n×t矩陣;W為t×d權重矩陣;A為d×r載荷矩陣;E為n×r殘差矩陣;式(9)為降秩回歸模型。
RDA計算結果可以將環境因子和物種因子等樣點投射到排序軸構成的二維平面上,并評價各因子之間的相關性和貢獻度(賴江山,2013)。
按照水華發生的統計標準(竇明等,2002),漢江武漢段在2008—2011年連續4年暴發水華事件,藻類密度均達到 1.0×107cell·L?1以上。由圖 3 可知,2004—2017年2—4月漢江武漢段宗關、白鶴嘴與琴斷口 3個觀測站藻類密度均值為 0.37×107cell·L?1;2008年水華事件中藻類密度達到最高值,即 5.3×107cell·L?1;2016 年藻類密度反彈上升,出現一次水華事件。從不同點位水華期間藻類密度比較上看,2008、2011、2016年水華暴發期間3個站點藻類密度峰值基本相近,2009年與2010年水華暴發期間 3個站點藻類密度峰值則呈現出宗關>琴斷口>白鶴嘴的趨勢。

圖3 2004—2017年漢江藻類密度(D)時間變化趨勢Fig.3 Variations of algae density (D) in Hanjiang River in 2004?2017
從STL分解的趨勢項變化上看(圖4、表1),2004—2017年,漢江3點位(宗關、琴斷口、白鶴嘴)藻類密度時間變化趨勢基本一致,2004—2008年呈現明顯的上升趨勢,2008—2017年整體為波動下降的趨勢。MK檢驗結果顯示,2004—2017年,藻類密度上升趨勢不顯著,2004—2008年藻類密度呈顯著上升趨勢(α=0.05),2008—2014年呈顯著下降趨勢(α≤0.1),2008年前后則是藻類密度變化的主要突變時段。

圖4 基于STL趨勢分解的2004—2017年漢江藻類密度趨勢項時間變化Fig.4 Trends of algae density in Hanjiang River in 2004?2017 based on STL method

表1 2004—2017年漢江藻類密度MK趨勢檢驗及拐點分析Table 1 Trend test and breakpoint analysis of algae density in Hanjiang River in spring 2004?2017 based on MK method
從季節變化上看(圖3,表2),3個站點的藻類密度變化具有明顯的季節性變化,2月中旬至 3月上旬為漢江武漢段水華易發期,且藻類密度呈現雙峰現象,其峰值分別在出現在2月中旬與3月上旬,說明這兩個時間較2月下旬比,更易發生水華;2008、2010、2016年峰值主要出現在 3月上旬,2009、2011年峰值出現在2月中旬;2004—2008年藻類密度最高值出現在3月上旬,2009—2017年其最高值出現在2月中旬,說明水華暴發時間有前移的趨勢。查看漢江中下游歷年水華發生情況統計數據,漢江水華主要發生在每年2—3月,主要優勢種均為適宜溫度較低的硅藻。2月上旬至3月下旬,漢江武漢段水體溫度均在0—15 ℃,均值溫度分別為6.8、9.0、7.7、9.8、13.0 ℃,且各旬溫度的最大值均超過 10 ℃,在營養鹽較充足的情況下,硅藻對低溫條件下也能有很好的適應性,這與辛小康等(2019)與夏瑞等(2020)研究結果一致,這一定程度上導致水華峰值呈前移趨勢。

表2 2004—2017年水華期間不同月旬漢江藻類密度峰值密度統計分布Table 2 Peaks of algae density in Hanjiang River during algae blooms in 2004?2017 in different month 107 cell·L?1
營養鹽的含量、溫度的適宜性以及水動力條件是影響藻類生長繁殖的主要影響因素,為解釋藻類時間變化趨勢和突變點,研究進一步分析了營養要素(TN、TP、N/P)、氣候要素(水溫t)、水文要素(流量Q、流速v、水位WL1與水位WL2,水位差ΔH)的時間變化趨勢與突變點分布。
從營養要素年變化趨勢上看(圖 5、表 3),2004—2017年2—4月TP處于 0.095—0.128 mg·L?1之間,與藻類密度的相關系數達到0.27(P<0.05),相關系數較小且相關性不明顯,這可能與TP變幅較小而藻類密度波動較大有關;但在所有的因子中TP仍與藻類密度的相關系數最大,且 2004—2017年TP時間變化趨勢不顯著。硅藻生長所需的TP為0.01 mg·L?1,TP>0.01 mg·L?1時,硅藻大量生長(Potapova et al.,2007);TP直接影響硅藻種群的穩定性,水華初期磷酸鹽濃度為0.06 mg·L?1,水華中期為 0.05 mg·L?1,而水華末期下降至 0.04 mg·L?1(Soininen et al.,2004)。2004—2017年2—4月TN處于 1.578—2.143 mg·L?1之間,且 2004—2017 年TN 時間變化趨勢不顯著。TN<0.2 mg·L?1時,low N硅藻大量生長,TN>3.0 mg·L?1時,high N 硅藻大量生長,TN 在 0.2—3.0 mg·L?1則為 low N 和 high N硅藻共同生長(Chen et al.,2004)。2004—2017年,漢江武漢段2—4月TN在low N和high N硅藻均可生長的濃度范圍,含量充足,這可能是造成藻類密度與TN含量沒有相關性的原因。2004—2017年2—4月N/P處于16.344—23.168之間,與藻類密度的相關性系數為?0.28(P<0.1),且 2004—2017年N/P無明顯上升趨勢(P>0.1)。在水華易發生的2月中旬到3月上旬,N/P均大于16,且隨著TP的變化出現相反的變化。研究表明漢江武漢段硅酸鹽 2月與 4月的含量均在 3 mg·L?1左右(殷大聰等,2017),其含量遠遠大于韓國Lower HanRiver硅藻水華(優勢種同樣為冠盤藻)發生時河段中硅的含量(0—1.5 mg·L?1)(Jung et al.,2009);硅負荷已不再是漢江水華的限制性因素;由此可見,與TN、N/P、硅酸鹽濃度相比,漢江武漢段水華的發生更易受到TP含量變化的影響,這與Xia et al.(2020)研究結果一致。

圖5 2004—2017年漢江藻類密度及相關因子時間變化趨勢Fig.5 Variations of algal density and influencing factors from 2004 to 2017

表3 2004—2017年漢江藻類密度及相關因子變化趨勢顯著性檢驗Table 3 Trends of algal density and influencing factors in Hanjiang River from 2004 to 2017 and their significance test with algae density
從氣候要素(環境因子)年變化趨勢上看,2004—2017年2—4月水溫處于9.593—14.437 ℃之間,與藻類密度的相關性系數達到?0.25(P<0.1),且2004—2017年無明顯的時間變化趨勢(P>0.1)。水華(Density>1.0×107cell·L?1)發生時,水體實測溫度范圍為10.011—13.031 ℃。漢江硅藻水華的主要優勢種一般為冠盤藻,其在低溫條件下也能有很好的適應性(夏瑞等,2020)。從水文要素年變化趨勢上看,2004—2017年2—4月水位、流量呈現明顯的下降趨勢,且流量與藻類密度相關系數為?0.21(P<0.1)。藻類密度與流速與水位均無明顯相關性。2004—2017年2—4月水位差范圍為7.714—10.056 m,與藻類密度的相關系數達到0.30(P<0.05);水位差對藻類密度的影響與漢江回水區的頂托對漢江下游河段的水動力條件變化密切相關(Qu et al.,2014)。由此可見,單獨水文要素的變化,并不足以引起漢江水華的發生,流量減少、流速下降,加之長江水位對漢口的頂托作用,易導致漢江武漢段出現水華。
從營養鹽含量、水文條件和氣候條件的變化趨勢與突變點時間分布上看(圖 6),TP、TN的 UF和 UB曲線交點主要出現在 2007—2008年和2015—2016年,水溫UF和UB曲線交點主要出現在2004—2007年和2014年,流量UF和UB曲線交點主要出現在2004年和2014年,水位差UF和UB曲線交點主要出現在2008—2009年和2014年。殷大聰等(2012)研究表明漢江武漢段TN和硅酸鹽等營養鹽負荷已經滿足春季暴發硅藻水華所需,TP的UF和UB曲線交點出現時間與藻類密度變化拐點(2008年藻類密度最高,2016年水華反彈)一致,表明TP與TN和硅酸鹽等營養鹽負荷相比,TP限制影響更大。夏瑞等(2020)發現,2008—2011年,漢江中下游主要河段年平均流量和水位均較1961—1999年明顯偏少,處于歷史上相對較低的一段時期。由此可見,漢江下游武漢段藻類密度的變化趨勢是多種因素共同作用的結果,營養負荷、氣候要素與水文要素綜合變化是漢江武漢段藻類密度2008年前后出現峰值及拐點,而2016年藻類密度反彈上升是出現一次水華事件的主要原因。

圖6 2004—2017年漢江環境因子與水文因子變化趨勢及突變分析Fig.6 Trends and mutation analysis of environmental and hydrological factors in Hanjiang River from 2004 to 2017
為對比分析水華年和非水華年影響因子的區別及影響,研究進一步選取 WL1為仙桃水文站的水位數據,代表漢江武漢段的水文情勢;WL2數據來源為漢口水文站水位數據,其位于漢口匯入長江口的下游,代表長江的水文情勢;水位差ΔH為漢口水文站與仙桃水文站水位差值,表征長江對漢江的頂托作用。
從2004—2017年2—4月漢江下游水華年與非水華年氣象水文要素統計結果上看(圖7),與非水華年相比,葉綠素濃度偏高14.4%,TN、TP均偏高1.9%左右,水溫偏低1 ℃,流速偏低0.05 m·s?1左右,仙桃、漢口平均水位分別偏低0.283、0.266 m,水位差偏小0.171 m,流量偏低15.0%。可見,水華發生流量小、氮磷負荷高且溫度偏低的年份。營養鹽含量已經不是河流藻類水華的主要限制因素(殷大聰等,2017);當河流中的總磷含量超過 0.030 mg·L?1時,營養鹽含量將不會成為主要限制因素(Sivapragasam et al.,2010)。Liu et al.(2012)研究表明,與湖泊富營養化相比,人類活動導致的河流水動力條件和水循環條件的變化對河流水華影響更大。Mitrovic et al.(2011)發現低流量提高了水的透明度,并創造了一個適當的光和溫度條件,以促進藻類生長和繁殖。盧大遠等(2000)指出快速水流會改變水平和垂直水交換率加速污染物的擴散和降解,最終降低河流中藻類繁殖率,從而抑制藻類生長。王紅萍等(2004)發現長江回水區的頂托對漢江下游河段的水動力條件也有較大影響,使得漢江武漢段水體滯留時間增加,加大了漢江水華發生的幾率;當長江水位較高時,假如漢江武漢段流量較小,漢江水位受到的頂托作用明顯,小流量對應高水位,當漢江流量較大時,長江頂托作用相對較小。盧大遠等(2000)指出,當漢江水流速度≤0.8 m·s?1、水流量≤500 m3·s?1、水溫≥10 ℃、pH>8.0、DO≥12 mg·L?1、色度≥15 度等是漢江硅藻水華指示因子的“預警值”;王紅萍等(2004)以藻類密度500×104cell·L?1為漢江春季水華預警值,提出宗關段水華預警流速為0.225 m·s?1。本研究中流量已經明顯大于500 m3·s?1;非水華年和水華年的流速均大于0.225 m·s?1;由此可見,低流速、小流量、高水位差等水動力條件易導致漢江硅藻生長和水華形成,近年漢江水華指示因子的“預警值”已經改變;漢江水華現象是多種因素綜合作用的結果,單一的限制因素分析已經滿足不了水華現象研究的需求。

圖7 2004—2017年2—4月漢江下游水華年與非水華年藻類密度及相關因子變化特征Fig.7 Characteristics of algal density and influencing factors in algal bloom year and non-algal bloom year in spring of 2004?2017 in Hanjiang River
為進一步厘清主要影響因子對藻類密度的影響和貢獻,先對水華年與非水華年漢江藻類密度與影響要素數據進行lg(x+1)轉換,然后基于RDA方法對其關系進行分析。由圖8可知,水華年前兩個排序軸對藻類密度變化的解釋率為40.9%,其中對藻類密度影響較大的影響要素為總磷與水位差,對葉綠素影響較大的影響要素為總磷,蒙特卡羅置換檢驗結果發現對漢江藻類變化趨勢的影響較為顯著的影響要素為水位差(F=7.4,P=0.010)與總磷(F=4.2,P=0.044),其他影響要素不顯著;非水華年前兩個排序軸對藻類密度變化的解釋率為35.1%,其中對藻類密度影響較大的影響要素為總磷、流量與水位,對葉綠素影響較大的影響要素為水位差與總氮,蒙特卡羅置換檢驗結果發現對漢江藻類變化趨勢的影響較為顯著的影響要素為流量(F=3.6,P=0.06)與總磷(F=2.8,P=0.088)。由此可見,磷是水華年與非水華年漢江藻類生長的主要限制營養元素,水華年水位差即長江與漢江水文交互作用對藻類密度生長影響顯著(Xia et al.,2020),非水華年漢江的流量對藻類密度影響顯著。

圖8 2004—2017年2—4月漢江下游水華年與非水華年RDA分析排序圖Fig.8 Results of RDA during algal and non-algal bloom years in Hanjiang River in 2004?2017
漢江武漢段是武漢市人民的飲用水水源地,水華的發生對飲用水安全造成極大的威脅,研究建立了基于拐點分析的河流水華診斷與歸因分析方法。綜合采用STL趨勢分解法、MK突變檢驗法、以及冗余分析等多源統計診斷與歸因分析方法,重點分析了2004—2017年2—4月漢江武漢段藻類密度的時間變化趨勢和突變點,識別了水華年和非水華年期間藻類密度主要影響因子。
(1)漢江武漢段水華年主要發生在2008—2011年和2016年,2月中旬至3月上旬藻類密度最大,均值為 0.79×107cell·L?1,2008 年水華事件為近 20年最嚴重。2004—2017年2—4月漢江武漢段宗關、白鶴嘴與琴斷口 3個觀測站藻類密度均值為0.37×107cell·L?1,峰值密度為 5.3×107cell·L?1;三站點藻類密度差異較小且時間變化趨勢一致,2004—2008年呈顯著上升趨勢,2008—2014年呈顯著下降趨勢,2008年前后為藻類密度突變點。
(2)營養負荷、氣候要素與水文要素綜合變化是導致 2004—2017年春季漢江水華現象的主要原因。近年漢江水華指示因子的“預警值”已經改變;總磷是水華年與非水華年漢江藻類生長的主要限制營養元素;低流速、小流量、低水位差等水動力條件易導致漢江硅藻生長和水華形成;與非水華年相比,總氮、總磷均偏高1.9%左右,水溫偏低1 ℃,流速偏低0.05 m·s?1左右,流量偏低15.0%。
(3)針對漢江水華的治理,建議除深入分析水華特征和影響因素外,應開展漢江總磷源解析,明確總磷來源,制定總磷污染防控對策;監控漢江武漢段的流量與漢口的水位,加強中上游丹江口水庫生態水量調度,預防“頂托”的形成。