水華過程水質參數與浮游植物定量關系的研究
——以太湖梅梁灣為例

郭文景,符志友,汪 浩,吳豐昌 (1.北京師范大學水科學研究院,北京 100875；.中國環境科學研究院,環境基準與風險評估國家重點實驗室,北京 10001)

湖泊富營養化及其可能造成的水華暴發問題可造成湖泊生態系統結構失衡和功能退化,因而一直受到關注[1-4].水華暴發中藻類產生的藻類毒素可以威脅流域內人體健康和供水安全[5-6],并可影響到湖泊的漁業價值和旅游價值,因此水華暴發一直是湖沼學研究的重點.對于水華暴發機制的探討有助于深入理解水華暴發的過程并展開相應的控制和人為干預措施,但是對于水華暴發的發生過程以及湖泊水質參數對藻類增殖的影響一直有眾多不同的研究.目前研究水華發生過程中湖泊水質參數與浮游植物的定量關系的途徑主要有:基于氣象/水文/營養鹽指標響應的確定性的生態數學模型,對于該問題的研究從單一的營養鹽限制關系逐漸發展到結合熱力學或生物相互作用等的復雜系統[7-9],另一途徑是基于觀測數據的線性/非線性統計歸納方法,主要包括主成分分析、多元回歸分析和人工神經網絡等[10-13].已有研究的基本邏輯為:浮游生物指標存在對湖泊水質參數的及時響應,而不存在時間上的滯后性.隨著對水華暴發過程認識的逐漸深入,湖泊水質參數和氣候條件對于浮游植物增殖產生的滯后效應已有研究提及,研究[14]認為夏季常見的藍藻在水面大量聚集的現象可能源自于春季藻類在湖泊水面下的大量增殖,在合適的氣象和水溫條件下,大量藻類上浮聚集形成通常意義上的藍藻暴發.因此考慮到外界脅迫參數可能對藻類水華暴發產生的滯后效應,基于同時段檢測數據的相關性分析得出的結論可能會有所缺陷,因此需要一種考慮滯后效應的數據分析方法用以研究水華暴發的過程,但截至目前,只有較少的研究將這一滯后效應考慮到目前水華暴發過程的研究中[15].

時間序列分析方法普遍用于計量經濟學領域,此外近年來也常用于環境經濟學領域,如探討經濟增長和環境污染的相關關系等[16-18],其包括了格蘭杰因果關系檢驗和向量自回歸模型,其通過構建考慮到自變量滯后效應的回歸模型對時間變量進行擬合,可以避免無關變量之間分析導致的偽相關關系的出現,其結果具有更高的準確性和更豐富的實際意義,也可以為經濟學和環境領域的研究人員和政策制定者提供有效的政策建議和決策支持.

因此本研究利用時間序列方法對太湖梅梁灣區域 2000～2012年的監測數據進行分析,建立監測數據序列的格蘭杰因果關系模型和向量自回歸模型,探討水華暴發中湖泊水質參數對于浮游植物生物量的定量影響.由于上述時間序列方法對于長時間尺度的時間變量有著較為嚴苛的要求,首先通過平穩性檢驗和協整檢驗對獲得水質參數的時間序列進行篩選,滿足條件的時間序列經過格蘭杰因果關系檢驗驗證因果關系,并構建向量自回歸模型,建立水質參數和浮游植物生物量之間的定量關系.

1 方法和數據處理

1.1 研究區域

太湖(30°90′N～31°54′N,119°55.3′E～120°59.6′E)是我國第三大淡水湖,位于我國經濟發達的長江三角洲地區,湖泊水面面積為2427km2,平均深度為 2m,是我國富營養化程度較高的湖泊之一.太湖流域快速的經濟發展和人口增長帶來的污染物輸入導致了太湖從上世紀七八十年代的貧營養化退化為目前的富營養化狀態[19].富營養化引起的藻類尤其是藍藻水華暴發現象時有發生,其中在2007年發生在貢湖灣的藍藻暴發事件更是引起了公眾對于城市飲用水安全的關注[20].梅梁灣位于太湖的北部,湖區水域面積為100km2,水深 1.8～2.3m.作為無錫市重要的水源地和旅游區,梅梁灣在進入 21世紀以來已經呈現了嚴重的富營養化狀態,而水華暴發也在該區域頻頻出現.梅梁灣區域相對于中心湖區受太湖風浪作用的影響較小,因此該區域適合作為水華暴發的研究區域.

1.2 數據收集和預處理

監測數據來自于中國科學院太湖生態系統研究站的逐月監測數據(TLLER).梅梁灣湖區共設置 5個采樣點(圖 1),數據采樣時間為 2000～2012年,主要檢測參數為:透明度(SD),酸堿度(pH值),亞硝態氮(NO2--N),硝態氮(NO3--N),銨態氮(NH4+-N),總氮(TN),總磷(TP),水溫(WT).浮游植物藻類生物量用水體葉綠素 a(Chl-a)的含量來指示.為了防止水面風速、水氣交換和人為活動等因素的干擾,水樣的采集在水面下0.5m處進行.上述指標的檢測均按照湖泊生態系統觀測標準方法進行[21].

NOx-N表示亞硝態氮(NO2--N)和硝態氮(NO3--N)的總含量.為了消除梅梁灣區域內風浪流和外源污染物造成的水質參數的差異,所有梅梁灣的逐月監測數據進行均值處理.為了消除數據可能產生的異方差波動,也為了使數據具有統計學意義,所有平均后的監測數據進行自然對數處理(pH值除外),如時間序列TN自然對數處理生成的數據記為lnTN.共得到9個時間序列,每個時間序列有n=156個逐月監測數據.

圖1 太湖梅梁灣湖區監測點位置(共5個)Fig.1 Five sampling sites in Meiliang Bay of Taihu Lake

1.3 平穩性檢驗

Granger[22]和 Sims[23]將數據必須滿足的平穩性定義為:時間序列在時間尺度上的平均值、方差和自協方差無波動.由于對不滿足平穩性時間變量的相關性分析可能會導致偽相關關系的出現,并得出錯誤的結論,因此建立變量之間格蘭杰因果關系模型和向量自回歸模型的必要條件是參與的所有時間序列必須通過平穩性檢驗.

平穩性檢驗的常見方法是增廣的 Dickey-Fuller test(ADF test)[24-25],其基本原理是檢驗時間序列是否存在單位根,如果時間序列不存在單位根,則該時間序列具有平穩性,反之,則不具有平穩性.主要基于以下回歸方程進行檢驗:

式中: xt代表時間序列的變量; t為時間;Δ表示一階差分; ρ為單位根存在與否(賦值0或1); p代表模型的滯后階數; αi為系數;εt表示均值為0、方差為σ2的白噪聲.ADF檢驗的零假設設定為“該時間序列存在一個單位根”,當檢驗中的P值大于設定的顯著性水平值,零假設被接受,該時間序列存在一個單位根,時間序列不平穩,反之,則時間序列平穩.

1.4 協整檢驗

對于未通過平穩性檢驗的非平穩序列,Granger[22]和 Sims[23]認為,因為若干非平穩時間變量的線性組合可能會保持平穩,這些非平穩變量之間仍然可以進行因果關系檢驗并建立向量自回歸模型,但是這些非平穩變量須滿足:1)時間序列為同階不平穩,2)時間序列須通過協整性檢驗.其中,非平穩序列進行n次差分后的殘差通過平穩性檢驗,則稱該序列為n階單整,記為I(n),相同階數的非平穩序列稱同階不平穩.而協整性檢驗主要基于最小二乘法(OLS)進行,其方程為:

式中: t是指時間; μt表示隨機擾動,如果隨機擾動項μt表現出平穩性,則時間序列 yt和 xt之間具有協整性.關于協整檢驗的更多原理和模型方法參見 Engle[26]的文獻.目前常見的協整檢驗的方法是Johansen-Juselius co-integration檢驗(JJ檢驗),其主要方法是通過最大特征根檢驗和跡檢驗進行判斷[27].協整檢驗的零假設設定為“時間序列中最多存在n(n=0,1,2…)個協整方程”.

1.5 格蘭杰因果關系檢驗和向量自回歸模型

滿足平穩性要求或通過協整檢驗后,時間序列可以建立格蘭杰因果關系模型和向量自回歸模型.Granger[22]從概率論的角度解釋了格蘭杰因果關系檢驗的實質:時間序列X是Y的格蘭杰原因意味著在一定的滯后期內 X可以在一定程度上解釋和預測Y的變化.其建立的方程如下:

式中: αi, βi, λi和 δi為系數; μ1t和 μ2t為白噪聲; m 為最大滯后階數.格蘭杰因果關系檢驗的零假設設定為“X不是Y的格蘭杰原因”.關于格蘭杰因果關系的介紹參考文獻 [28].

通過了格蘭杰因果關系檢驗的時間變量構建穩定的解釋力較強的向量自回歸模型(VAR),其可以定量解釋時間序列之間在時間尺度上的定量關系,這一模型考慮了時間尺度上自變量對于因變量的滯后效應,因此更具有說服力.

1.6 數據分析

平穩性檢驗、協整檢驗、格蘭杰因果關系檢驗和向量自回歸模型的建立均通過常用的時間序列分析軟件—Eviews(Econometrics Views 6.0)進行.在構建平穩性檢驗、協整檢驗和格蘭杰因果關系模型時,零假設的顯著性臨界值設定為0.05.進行時間序列分析時,滯后階數對于模型的解釋能力十分重要,在本研究中,N=12確定為最大滯后階數,而最優滯后階數基于赤池信息準則(AIC)和施瓦茨準則(SC)確定[29-30].

圖2 梅梁灣湖區水體中Chl-a含量逐月平均數據的變化趨勢(2000～2012年)Fig.2 Long-term monthly variation of concentration of Chl-a in Meiliang Bay from 2000 to 2012

圖3 2000～2012年間梅梁灣區域若干水質參數逐月均值數據的變化趨勢Fig.3 Long-term monthly variations of several aquatic parameters and concentrations of nutrients in Meiliang Bay from 2000 to 2012

2 結果和討論

2.1 Chl-a濃度和水環境指標的時間變化趨勢

圖 2顯示的是梅梁灣湖區 Chl-a含量在 2000～2012年之間的月度變化,其最小值為0.56mg/m3(2001年 2月),最大值為 491.15mg/m3(2009年7月).Chl-a含量的逐月數據在表現出波動性的同時也表現出年際的增長性,其中在夏季達到峰值,而在冬季達到最低值.湖泊水質參數的變化趨勢見圖3,主要趨勢為在冬季和夏季達到極值.

2.2 平穩性檢驗結果

太湖梅梁灣區域時間序列的平穩性檢驗結果如表1所示,其中水質參數變量lnTN、lnSD、pH值、lnNH4+-N、lnNOx-N和lnNH4+-N/NOx-N平穩性模型的P值都小于0.05,因此零假設被拒絕,上述時間序列具有平穩性,Chl-a含量的時間變量(lnChl-a)由于零假設被接受,因此呈現出非平穩性.由于平穩性變量難以解釋非平穩性變量的變化趨勢,因此平穩性變量不認為是非平穩性變量Chl-a含量變化的格蘭杰原因,所以lnTN、lnSD、pH值、lnNH4+-N、lnNOx-N和lnNH4+-N/NOx-N在接下來的分析中被剔除.

對非平穩序列lnTP、lnN/P、lnWT和lnChl-a進行一階差分,并分別記為dlnTP、dlnN/P、dlnWT和 dlnChl-a.對上述序列的殘差進行平穩性檢驗,結果表明,其一階差分表現出平穩性,因此lnTP、lnN/P、lnWT和 lnChl-a為一階單整變量(記為I(1)),符合進行后續協整檢驗的要求.

表1 梅梁灣水質參數時間序列的ADF檢驗結果Table 1 Results of ADF unit root test for various aquatic parameters

2.3 協整檢驗結果

對梅梁灣水質參數和葉綠素變量 lnTP、lnN/P、lnWT和lnChl-a進行協整分析,使用方法為JJ協整分析,通過跡檢驗和最大特征根檢驗進行判斷是否存在協整關系.

表2 梅梁灣監測數據的JJ協整檢驗結果Table 2 Results of the Johansen-Juselius co-integration test for environmental parameters

水質參數和葉綠素變量lnTP、lnN/P、lnWT和lnChl-a的協整檢驗結果如表2所示,跡檢驗的統計值均大于 5%的臨界值,跡檢驗和最大特征根的結果表明上述變量在時間尺度上分別存在4個和 2個協整方程,即上述時間變量之間至少存在1個協整方程.

2.4 格蘭杰因果關系檢驗結果

對梅梁灣湖區的時間序列 lnTP、lnN/P、lnWT和lnChl-a呈現了同階不平穩,且通過了協整檢驗,滿足進行格蘭杰因果關系檢驗的條件,因此使用 Eviews軟件對水質參數 lnTP、lnN/P、lnWT和Chl-a含量lnChl-a的因果關系進行檢驗,結果如表3所示:

零假設“lnTP 不是lnChl-a的格蘭杰原因”的P值為0.0001,因此該假設被拒絕,湖泊水體中TP的變化是Chl-a變化的格蘭杰原因.另一個零假設“lnChl-a不是lnTP的格蘭杰原因”的P值為 0.8473,零假設被接受,因此,從統計學意義上湖泊水體中Chl-a的變化不是TP變化的格蘭杰原因,存在從湖泊總磷變化到 Chl-a含量變化的單向格蘭杰因果關系.格蘭杰因果關系的結果表明在滯后期n=2的時間尺度內,總磷的增加對浮游植物的增殖有促進作用,即表示研究區域冬季TP的變化可能會影響到春季浮游植物的生物量.

表3 Chl-a含量與湖泊水質參數變量的格蘭杰因果關系檢驗結果Table 3 Results of the Granger causality test for lnChl-a and various independent parameters

檢驗湖泊水體N/P和Chl-a含量之間關系時,其零假設“lnN/P不是 lnChl-a變化的格蘭杰原因”的P值為0.0098,因此lnN/P是lnChl-a變化的格蘭杰原因.而“lnChl-a不是lnN/P變化的格蘭杰原因”零假設的P值為0.6253,零假設被接受,lnChl-a不是lnN/P變化的格蘭杰原因,因此存在長期時間尺度內的從N/P到Chl-a的單向格蘭杰因果關系.

在檢驗湖泊水溫WT和Chl-a之間的因果關系時,“lnWT 不是lnChl-a的格蘭杰原因”的零假設被拒絕(P=0.0424),因此從統計學意義上水體溫度是Chl-a含量變化的格蘭杰原因,“lnChl-a不是 lnWT的格蘭杰原因”的零假設被接受(P=0.1053),因此可以認為lnChl-a不是lnWT的格蘭杰原因.因此同樣存在長期時間尺度內的從水溫WT到Chl-a的單向格蘭杰因果關系.格蘭杰因果關系結果表明在研究區域內,湖泊水溫的變化可以在長期尺度上影響 Chl-a濃度,因此研究區域冬季的溫度變化可能會影響次年秋季或冬季湖泊中水華暴發的強度.

2.5 向量自回歸模型建立

向量自回歸模型可以從時間尺度上定量檢驗自變量對于因變量的影響[31-32],同樣使用Eviews軟件建立格蘭杰因果關系模型中自變量TP、TP和WT與因變量Chl-a之間的向量自回歸模型(5)(括號內的數字為系數的標準差).

長時間尺度內的向量自回歸模型顯示,如果其他條件保持不變,水體中 TP的含量變化±1%,Chl-a的含量將會有±0.97%的變化,而1%的N/P變化只會造成0.078%Chl-a濃度的變化,1%水體溫度的變動會導致0.55%Chl-a含量的變化.即水體中TP濃度對于研究區域內浮游植物生物量的變化影響較大,水溫的影響次之,而N/P的影響較小.

2.6 解釋和討論

在方程(5)中,研究區域湖泊水體中水質參數TP、N/P、水溫和Chl-a含量呈現了長期均衡關系.關于營養鹽對于浮游植物增殖的影響一直存在爭議[4,33],本研究中格蘭杰因果關系檢驗的結果表明在研究時間段內存在單向的TP到Chl-a的格蘭杰因果關系,即梅梁灣湖區 TP的濃度對水體中浮游植物生物量存在較為顯著的促進作用,且對于浮游植物的影響遠大于其他水質參數.在太湖的相關研究中,Guo等[34]使用冗余分析和TITAN分析研究了太湖梅梁灣和五里湖地區營養鹽對浮游植物的影響,發現 TP可能是影響該區域藻類生長的限制因子.Xu等[35]通過原位實驗和野外采樣,估計了藻類暴發的營養鹽閾值為0.20mg TP/L.在本研究區域中,營養鹽的濃度為0.14mg TP/L (0.006～4.89mg/L).因此在該地區長時間尺度內富營養化過程的浮游植物增殖更多受到 TP的影響,且長時間的監測數據進一步表明TP的影響最大.對于N/P來說,格蘭杰因果關系檢驗表明N/P存在對Chl-a濃度變化的單向影響,但是向量自回歸模型的結果表示,1%的 N/P的變化只會造成 0.078%Chl-a含量的變化,相關性分析的結果也表明N/P與Chl-a濃度之間的相關關系不明顯(相關系數為0.32,P＞0.05,n=156,文章中未顯示).另外,研究區域內氮磷比為 35.07(0.03～127.86),高于 Redfield[36]對水體中藻類生長適宜氮磷比(N/P=16)的閾值.因此 N/P的變化并不是影響湖區中 Chl-a含量變化的主要因素.而對于水溫而言,格蘭杰因果關系檢驗的結果表明存在從湖泊水溫到浮游植物生物量的單向格蘭杰因果關系.溫度對藻類增殖的影響早有研究提及,Abrantes等[37]研究表明溫度是影響浮游植物生長的主要因素,LüRLING 等[38]在對藍藻增殖的研究中也得出了相似的結論.向量自回歸模型的結果進一步表明在該研究區域湖泊水溫對浮游植物生物量的影響要小于 TP的變化,但是遠大于N/P的變化.

TP、N/P、WT和Chl-a指標在平穩性檢驗中呈現了不平穩性(I(1)),而梅梁灣區域其他水質參數(lnTN、lnSD、pH、lnNH4+-N、lnNOx-N和 lnNH4+-N/NOx-N)的原始數據呈現平穩性(I(0)),對于不同階時間序列的相關性分析可能會導致偽回歸的產生,并且導致錯誤的結論,因此從統計學意義上,上述參數不作為水體中Chl-a含量的格蘭杰原因.對于氮源來說,梅梁灣區域研究時段內 TN的含量為 3.45(0.02～13.56mg/L),而在該區域的研究結果顯示,水華現象發生的TN閾值為0.8mg/L[35,39].因此梅梁灣區域內氮源對于浮游植物生長的限制作用可能十分有限.而對于pH值和SD來說,本研究的結果可以解讀為水體中pH值和SD的變化無法預測和影響水體中浮游植物生物量的變化,事實上水體中pH值和SD的變化可能是浮游植物增殖的結果,主要原因分別為浮游植物增殖消耗水體中CO2和HCO3-從而增加水體pH值[40]以及藻類增殖減少水體透明度.

在已有的對太湖區域藻類生物量進行模型預測工作時,陳宇煒等[10]也利用長時間的太湖水質參數和藻類生物量進行回歸統計,結果發現水體中 TP和水溫對于藻類生物量具有顯著影響,與本文結果較為相似.鄧建明等[41]利用廣義可加模型和太湖監測數據,分析自變量對于因變量的非線性關系,發現除總磷和溫度外,化學需氧量和電導率也是影響微囊藻生物量的關鍵環境因子,可能原因是模型選擇的差異.在少量時間序列分析方法在湖沼學領域的應用中,Lui等[15]利用向量自回歸模型擬合了香港吉澳島附近海域葉綠素熒光與相關水質參數的時間趨勢模型,該研究同樣也發現了時間變量之間的滯后效應.龔紹琦等[42]使用時間序列方法中的自回歸滑動平均模型(ARMA)建立了太湖不同湖區內葉綠素濃度的時間序列模型,并對一定時間內葉綠素濃度進行預測,取得了較好的預測精度.本文研究利用時間序列方法對太湖梅梁灣監測數據進行分析,發現了在一定滯后期內湖泊水質參數對于浮游植物生物量的影響,這一結論佐證了孔繁翔等提出的藍藻水華暴發的四階段假說-休眠、復蘇、生物量增加(生長)、上浮及聚集等四個階段[14].向量自回歸模型的結果定量的闡述了在該研究區域內不同水質參數的影響大小,梅梁灣湖區水華增殖過程中,TP的變化強烈導致了水體中Chl-a的變化,說明湖泊水體中 TP的輸入有利于水華暴發的形成,且 TP發揮的作用大于水體溫度和其他營養鹽,這可能和該流域內污染物的排放特征有關.而相對于其他營養鹽(如氮源),控制流域內的P輸入和湖泊內源P釋放對于抑制水華暴發可能會起到更大的作用.另外水溫的變化也導致了浮游植物生物量的變化,表明在長時間尺度內水溫的升高可能會提高水體中 Chl-a的含量.因此在全球氣候變暖的前提下水溫對藻類水華暴發的促進作用暗示了在氣候變暖的全球趨勢下,升高的水體溫度將提高該區域藻類水華暴發的概率和強度[43].

3 結論

3.1 太湖梅梁灣區域湖泊水質參數中TN、SD、pH、NH4+-N、NOx-N和NH4+-N/NOx-N不能解釋Chl-a含量變化的原因,而TP、N/P和WT與Chl-a之間均存在單向的格蘭杰因果關系.

3.2 研究區域內TP、N/P、WT和Chl-a含量之間存在長期均衡關系,在保持其他條件不變的基礎上,1%TP、N/P和WT的變化將會導致0.97%、0.078%和0.55%Chl-a含量的變化.水體中TP的濃度在水華暴發過程中的作用可能高于其他營養鹽(如不同形式的氮源),水溫的影響次之.

值得注意的是,時間序列方法較少在湖泊富營養化領域中應用,其表述的因果關系可以在數理統計學上對應用的場景進行解釋,同時模型的意義也需要結合實際的狀況進行闡述.另外本研究由于缺少對應的水文和氣象數據,因此未將上述因素考慮在內,但是水位、降水量和風力條件等因素對于水華暴發也有重要影響,因此在下一步的研究中應將這些因素考慮在內.

參考文獻：

[1]Heisler J, Glibert P M, Burkholder J M, et al. Eutrophication and harmful algal blooms: a scientific consensus [J]. Harmful algae,2008,8(1):3-13.

[2]Paerl H W, Otten T G. Harmful cyanobacterial blooms: causes,consequences, and controls [J]. Microbial Ecology, 2013,65(4):995-1010.

[3]Merel S, Walker D, Chicana R, et al. State of knowledge and concerns on cyanobacterial blooms and cyanotoxins [J].Environment International, 2013,59(3):303-327.

[4]馬健榮,鄧建明,秦伯強,等.湖泊藍藻水華發生機理研究進展[J]. 生態學報, 2013,33(10):3020-3030.

[5]Jochimsen E M, Carmichael W W, An J, et al. Liver failure and death after exposure to microcystins at a hemodialysis center in Brazil [J]. New England Journal of Medicine, 1998,338(13):873-878.

[6]De Figueiredo D R, Azeiteiro U M, Esteves S M, et al.Microcystin-producing blooms—a serious global public health issue [J]. Ecotoxicology and Environmental Safety, 2004,59(2):151-163.

[7]Vollenweider R A. Input-output models [J]. Schweizerische Zeitschrift für Hydrologie, 1975,37(1):53-84.

[8]Ye C, Shen Z, Zhang T, et al. Long-term joint effect of nutrients and temperature increase on algal growth in Lake Taihu, China [J].Journal of Environmental Sciences, 2011,23(2):222-227.

[9]李 博,孫冬梅,馮 平,等.于橋水庫浮游藻類時空分布規律的數值模擬研究 [J]. 中國環境科學, 2013,33(3):508-515.

[10]陳宇煒,秦伯強,高錫云.太湖梅梁灣藻類及相關環境因子逐步回歸統計和藍藻水華的初步預測 [J]. 湖泊科學, 2001,13(1):63-71.

[11]Rigosi A, Hanson P, Hamilton D P, et al. Determining the probability of cyanobacterial blooms: the application of Bayesian networks in multiple lake systems [J]. Ecological Applications A Publication of the Ecological Society of America, 2015,25(1):186-199.

[12]陳云峰,殷福才,陸根法.富營養化水體水華暴發的突變模型 [J].中國環境科學, 2006,26(1):125-128.

[13]鄭劍鋒,焦繼東.基于神經網絡的城市內湖水華預警綜合建模方法研究 [J]. 中國環境科學, 2017,37(5):1872-1878.

[14]孔繁翔,高 光.大型淺水富營養化湖泊中藍藻水華形成機理的思考 [J]. 生態學報, 2005,25(3):589-595.

[15]Lui G C, Li W K, Leung K M, et al. Modelling algal blooms using vector autoregressive model with exogenous variables and long memory filter [J]. Ecological Modelling, 2007,200(1):130-138.

[16]Bozoklu S, Yilanci V. Energy consumption and economic growth for selected OECD countries: Further evidence from the Granger causality test in the frequency domain [J]. Energy Policy, 2013,63(3):877-881.

[17]Wang S, Li Q, Fang C, et al. The relationship between economic growth, energy consumption, and CO2emissions: empirical evidence from China [J]. Science of the Total Environment, 2016,542(Pt A):360-371.

[18]Xiao F, Brajer V, Mead R W. Blowing in the wind: the impact of China's Pearl River Delta on Hong Kong's air quality [J]. Science of the Total Environment, 2006,367(1):96-111.

[19]Chen Y, Fan C, Teubner K, et al. Changes of nutrients and phytoplankton chlorophyll-a in a large shallow lake, Taihu, China:an 8-year investigation [J]. Hydrobiologia, 2003,506(1-3):273-279.

[20]Qin B, Zhu G, Gao G, et al. A drinking water crisis in Lake Taihu,China: linkage to climatic variability and lake management [J].Environmental Management, 2010,45(1):105-112.

[21]金相燦,屠清瑛.湖泊富營養化調查規范 [M]. 北京:中國環境科學出版社, 1990:138-207.

[22]Granger C W. Investigating causal relations by econometric models and cross-spectral methods [J]. Econometrica: Journal of the Econometric Society, 1969,37(3):424-438.

[23]Sims C A. Money, income, and causality [J]. The American economic review, 1972,62(4):540-552.

[24]Dickey D A, Fuller W A. Distribution of the estimators for autoregressive time series with a unit root [J]. Journal of the American statistical association, 1979,74(366a):427-431.

[25]Dickey D A, Fuller W A. Likelihood ratio statistics for autoregressive time series with a unit root [J]. Econometrica:Journal of the Econometric Society, 1981,49(4):1057-1072.

[26]Engle R F, Granger C W. Co-integration and error correction:representation, estimation, and testing [J]. Econometrica: Journal of the Econometric Society, 1987,55(2):251-276.

[27]Chang C-C. A multivariate causality test of carbon dioxide emissions, energy consumption and economic growth in China [J].Applied Energy, 2010,87(11):3533-3537.

[28]曹永福.格蘭杰因果性檢驗評述 [J]. 數量經濟技術經濟研究,2006,23(1):155-160.

[29]Akaike H. A new look at the statistical model identification [J].IEEE transactions on automatic control, 1974,19(6):716-723.

[30]Schwarz G. Estimating the dimension of a model [J]. The Annals of Statistics, 1978,6(2):461-464.

[31]高 山,黃 楊,王 超.貨幣政策傳導機制有效性的實證研究——基于我國利率傳導渠道的VAR模型分析 [J]. 財經問題研究, 2011,(7):50-58.

[32]高鐵梅.計量經濟分析方法與建模/EViews 應用及實例/數量經濟學系列叢書:EViews應用及實例 [M]. 清華大學出版社有限公司, 2006:267-314.

[33]Paerl H W, Xu H, McCarthy M J, et al. Controlling harmful cyanobacterial blooms in a hyper-eutrophic lake (Lake Taihu,China): The need for a dual nutrient (N & P) management strategy [J]. Water Research, 2011,45(5):1973-1983.

[34]Guo W, Zhu Y, Fu Z, et al. Influences of environmental factors on biomass of phytoplankton in the northern part of Tai Lake, China,from 2000 to 2012 [J]. Environmental Monitoring and Assessment, 2017,189(12):608-617.

[35]Xu H, Paerl H W, Qin B, et al. Nitrogen and phosphorus inputs control phytoplankton growth in eutrophic Lake Taihu, China [J].Limnology and Oceanography, 2010,55(1):420-432.

[36]Redfield A C. The biological control of chemical factors in the environment [J]. Science Progress, 1960,11(11):150-170.

[37]Abrantes N, Antunes S, Pereira M, et al. Seasonal succession of cladocerans and phytoplankton and their interactions in a shallow eutrophic lake (Lake Vela, Portugal) [J]. Acta Oecologica, 2006,29(1):54-64.

[38]L rling M, Eshetu F, Faassen E J, et al. Comparison of cyanobacterial and green algal growth rates at different temperatures [J]. Freshwater Biology, 2013,58(3):552-559.

[39]Xu H, Paerl H W, Qin B, et al. Determining critical nutrient thresholds needed to control harmful cyanobacterial blooms in eutrophic Lake Taihu, China [J]. Environmental Science &Technology, 2014,49(2):1051-1059.

[40]Oliver R L, Ganf G G. Freshwater blooms. In B. A. Whitton & M.Potts (Eds.), The ecology of cyanobacteria: their diversity in time and space [M]. Dordrecht: Springer Netherlands, 2002:149-194.

[41]鄧建明,秦伯強,王博雯.廣義可加模型在 R 中的快捷實現及藍藻水華預測分析中的應用 [J]. 生態學雜志, 2015,34(3):835-842.

[42]龔紹琦,黃家柱,李云梅,等.應用時間序列分析法對太湖葉綠素-α含量的動態研究 [J]. 海洋與湖沼, 2008,39(6):591-598.

[43]商兆堂,任 健,秦銘榮,等.氣候變化與太湖藍藻暴發的關系[J]. 生態學雜志, 2010,29(1):55-61.