基于GAM模型的鄱陽湖葉綠素a與水質因子相關性分析

袁偉皓 ，王華 *，夏玉寶 ，曾一川 ，鄧燕青，李媛媛 ，張心悅

1. 河海大學環境學院，江蘇 南京 210098；2. 河海大學淺水湖泊綜合治理與資源開發教育部重點實驗室，江蘇 南京 210098；3. 江西省水文局，江西 南昌 330000

鄱陽湖是中國長江中游典型的通江湖泊，也是中國第一大淡水湖，位于江西省，地處九江、上饒、南昌三市（115°50′—116°44′E，28°25′—29°45′N），上游承接贛江、撫河、信河、饒河、修河5條主要河流來水，經湖區調蓄后由湖口注入長江，是一個季節性較強的吞吐型湖泊，在維系區域水量平衡與生態安全方面發揮著重要作用（陳煉鋼等，2020；俞珊妮等，2020；陳旻坤等，2021）。近年來，隨著鄱陽湖流域地區社會經濟的快速發展，大量氮、磷等營養物質排入鄱陽湖水體，導致富營養化水平呈逐年上升趨勢，營養鹽總體維持在中營養水平，但是上升速度不斷加快，一些區域甚至已經呈現富營養化水平（劉靜等，2020；溫春云等，2020）。

氮磷是藻類增殖的物質基礎，其與藻類生物量之間的關系一直是研究湖泊富營養化的重點之一（劉靜，2020；錢奎梅等，2021）。葉綠素a（Chl-a）是藻類重要的表征指標，水體中Chl-a含量的多少反映了浮游植物生物量的高低，其濃度與水環境質量密切相關，是水體理化性質動態變化的綜合反映指標，在水體富營養化評價中起關鍵作用（陳金紅等，2020）。目前，眾多研究關注湖泊水質因子與Chla之間的相關性，但得出的結論都是線性、對數和正負相關等（孫越峰等，2020；朱廣偉等，2020；曾一川等，2020），由于Chl-a與各水質因子之間的關系非常復雜，同時各類湖泊具有其獨特的性質，尤其是針對具有典型通江特性的鄱陽湖。

常用的線性回歸分析可能會忽略解釋變量（Explanatory variable）與響應變量（Response variable）之間的非線性關系；同時，傳統的線性分析方法不能適用多變量之間的統計學分析，而廣義加性模型（Generalized Additive Model，GAM）是一種非參數的分析模型，能直接處理響應變量和多個解釋變量之間的非線性關系，不需要假定數據的分布，具有很高的靈活性，近年來在環境等眾多領域具有非常廣泛的應用（Peng et al.，2020；賀祥等，2017；丁隆強等，2020）。但是，關于GAM模型在鄱陽湖水環境研究中的應用相對較少。程新等（2017）運用GAM模型對鄱陽湖表層水體中藻類豐度與環境因子之間的關系進行了研究，研究表明：2012年10月—2013年7月，鄱陽湖藻類豐度與水溫和葉綠素a呈極顯著的正相關關系，與氨氮、總磷呈極顯著負相關關系。因此，深入探究鄱陽湖不同區域水質惡化和富營養化的原因顯得尤為重要，構建GAM模型對鄱陽湖較長時間序列的Chl-a與水質因子進行相關性研究也十分必要。

本文運用非線性的分析方法，探討了鄱陽湖不同區域的Chl-a變化的限制性水質指標與主要的可能因素，給出了相應的可能有效解決途徑，具有十分重要的理論和實踐意義，也為鄱陽湖長遠的富營養化防控和有效科學的管理提供了依據。

1 材料與方法

1.1 樣品采集與測定

本研究選取鄱陽湖主湖區3個典型點位，由西向東分別為蚌湖（115°58′49″E，29°18′22″N）、都昌（116°11′35″E，29°14′36″N）和蛇山（116°22′57″E，29°05′14″N），于 2016 年 1 月—2020 年 12 月開展了近5年的連續野外監測，采樣頻次為每月1次，每次用 2 L有機玻璃采水器采集距離湖水表層0.5 m的3組平行水樣，采集的水樣用500 mL聚乙烯塑料瓶收集，現場測定水溫（WT）、透明度（SD）等參數，后及時帶回實驗室測定葉綠素a（Chl-a）、總氮（TN）、總磷（TP）、氨氮（NH3-N）和高錳酸鹽指數（CODMn）。其中，葉綠素a采用丙酮萃取分光光度計法測定（楊玉珍等，2011）；總氮和總磷分別采用堿性過硫酸鉀消解-紫外分光光度法（HJ 636—2012）及鉬酸銨分光光度法（GB 11893-89）測定；氨氮采用納氏試劑分光光度法測定（HJ 535—2009）；高錳酸鹽指數則采用高錳酸鹽指數測定法（GB 11892-89）；透明度采用塞式盤法測定；水溫等則為采樣現場儀器直接測定。

1.2 數據分析方法

GAM模型是廣義線性模型（Generalized Linear Model，GLM）的半參數拓展，可直接擬合因變量與多個自變量之間的非線性關系，對不同的函數進行加和，找出其中的規律，適應于各種不同類型分布的函數分析（郭亮等，2017；楊佳星等，2020）。其一般表達形式如下：

式中：

g(y)——連接函數；

s(x)——連接解釋變量的光滑函數；

φ——隨機殘差項。

GAM模型的運算通過R語言實現，數據分析通過R中的mgcv軟件包實現。運用GAM模型對鄱陽湖3個研究點位Chl-a與各環境因子分別進行相關性分析，分析步驟如下：（1）基于chl-a與環境因子的Pearson相關系數，確定2個主要解釋變量；（2）基于解釋變量間的 Pearson相關系數大小，進行共曲線性診斷（當相關系數大于0.5時，可認為2個解釋變量之間存在共曲線性，需剔除1個變量；反之，則不存在，無需剔除）；（3）基于變量預分析結果，確定GAM模型的擬合方程，代入得相關性（r2）、自由度（edf）、顯著性取值（P值）、解釋率（Deviance Explained，D-E）和擬合圖像等結果信息（當edf=1時，變量間呈現線性關系，其值越大則代表非線性影響能力越強；P值代表統計結果的顯著性水平，與評估各因子間的相關性水平類似，“*”越多顯著性越強。“***”對應P<0.01，表示極顯著，結果可信度高；D-E為模型對變量總體變化的解釋率）（張智淵等，2018）。但樣品檢測結果存在一定量的異常值（有異常過大及過小的現象出現），為保證較好的代表性，異常的數據結果在后面結果中不予呈現。

2 結果與討論

2.1 葉綠素a含量變化特征

總體而言，研究期間（2016年1月—2020年12月）3個研究點位Chl-a質量濃度的波動范圍為0.0002—0.0510 mg·L?1，鄱陽湖各點位間 Chl-a 質量濃度隨時間的變化明顯，夏季Chl-a質量濃度普遍高于其他季節，具體如圖1所示。蚌湖點位質量濃度高于其他測點，5年內整體呈現先下降后上升的趨勢，且于2018—2020年夏季出現多次峰值，最大值出現在 2018年 10月，質量濃度達 0.0510 mg·L?1，5 年均值為 0.0124 mg·L?1；都昌點位 Chla 質量濃度介于 0.0002—0.0219 mg·L?1之間波動，整體趨勢較為穩定，分別于2016—2017年和2020年夏季出現3次峰值，最大值出現在2020年10月，5年均值為 0.0062 mg·L?1；蛇山點位Chl-a質量濃度介于 0.0006—0.0234 mg·L?1之間波動，整體趨勢也相較穩定，分別于2017年和2020年夏季出現兩次峰值，最大值出現在2017年8月，5年均值為0.0057 mg·L?1。

圖1 研究點位葉綠素a含量變化Fig. 1 Variation of Chl-a concentration at the research sites

就蚌湖和都昌點位而言，其位于贛江下游流域，較為直接地受到贛江水質的影響。盛文濤等（2021）研究發現，近年來贛江的營養鹽濃度及TN/TP顯著高于其他4河，同時贛江年徑流量位于5河之首，豐水期贛江攜帶大量營養鹽通量進入鄱陽湖，在其下游附近逐漸出現滿足藻類增殖條件的區域，因此該區域內Chl-a多會在豐水期結束后的時期（10月左右）達到峰值。而對于蛇山點位而言，其位于鄱陽湖東北部，李云良等（2016）研究發現鄱陽湖于高洪水位季節，污染物（營養鹽等）因東北部湖灣區存在順時針方向環流而發生長時間滯留和富集；因此，蛇山點位于8月左右達藻類增殖最佳條件（流速、營養鹽和水溫等），這也解釋了蛇山點位Chl-a質量濃度峰值多出現在8月的原因。

2.2 環境因子變化特征

2016年1月—2020年12月鄱陽湖3個研究點位的環境因子隨時間的變化見圖2。就3個點位的水溫而言，點位之間的差異不大，5年內呈現夏季高，冬季低的特點，夏季水溫介于33.4—34.4 ℃之間波動；冬季水溫則介于2.1—2.7 ℃之間波動。就透明度而言，3個點位的SD介于0.1—1.5 m之間波動，隨時間變化明顯，3個點位SD 5年均值由西向東分別 0.44、0.48、0.45 m。就總氮質量濃度而言，3個點位TN質量濃度呈現整體下降趨勢，點位間差異不大；蚌湖點位TN質量濃度介于0.27—3.51 mg·L?1之間波動，5 年均值為 1.24 mg·L?1，于2016—2018年冬季出現多次峰值；都昌點位TN質量濃度介于0.19—2.45 mg·L?1之間波動，5年均值為1.25 mg·L?1；蛇山點位TN質量濃度介于0.34—2.78 mg·L?1之間波動，5 年均值為 1.43 mg·L?1。就總磷而言，3個點位整體質量濃度介于 0—0.125 mg·L?1之間波動，趨勢相對穩定，蛇山點位于2016年3月和2017年3月出現了兩次峰值，分別達到了 0.206 mg·L?1和 0.228 mg·L?1；3 個點位 5 年均值分別為 0.053、0.058 和 0.061 mg·L?1，自西向東呈現遞增趨勢。就氨氮而言，其變化趨勢整體上與TN相似，于2016—2018年內波動顯著，后兩年趨向于穩定，各點位多年均值分別為0.236、0.157、0.188 mg·L?1。就高錳酸鹽指數而言，蚌湖的CODMn顯著大于其他點位，5年來整體質量濃度趨向于穩定，介于 1.0—3.0 mg·L?1之間波動。

圖2 研究點位環境因子隨時間的變化Fig. 2 Variation of environmental factors at the research sites

2.3 葉綠素a與環境因子相關性

2.3.1 蚌湖葉綠素a與環境因子的相關性

對于蚌湖點位，根據Chl-a與環境因子間的相關系數判斷，TN和CODMn與Chl-a之間的相關系數較高，因此選取 TN和 CODMn為解釋變量代入GAM模型；同時，TN和CODMn間的相關系數為0.18，故排除共曲線性可能，無需剔除變量。根據變量預分析結果，將TN和CODMn代入如下模型方程擬合：

式中：

x1——蚌湖點位 TN 質量濃度（mg·L?1）；

x2——蚌湖點位 CODMn質量濃度（mg·L?1）。得到r2為0.435，TN和CODMn的自由度分別為4.79和2.40，解釋率為52.3%。

從圖3中可以看出，蚌湖點位Chl-a與TN和CODMn存在非線性關系，TN的非線性影響能力更強，Chl-a與TN的關系更為復雜，當TN質量濃度為0—1 mg·L?1時，TN對Chl-a的影響最為顯著，兩者呈現明顯正相關關系；當TN的質量濃度為1—2.6 mg·L?1時，Chl-a質量濃度會隨著TN的變化而略微減少；當TN的質量濃度大于2.6 mg·L?1時，二者又再次表現較明顯的正相關關系。CODMn與Chl-a呈現先增大后減小的變化趨勢，當CODMn的質量濃度為0—3.3 mg·L?1時Chl-a質量濃度會隨著CODMn的變化而顯著增加；當 CODMn的質量濃度大于3.3 mg·L?1時，Chl-a質量濃度會隨著 CODMn的增加而減少。

圖3 蚌湖TN和CODMn與Chl-a質量濃度的關系Fig. 3 The relationship between TN, CODMn and Chl-a in Banghu

2.3.2 都昌葉綠素a與環境因子的相關性

對于都昌點位，根據 chl-a與環境因子間的相關系數判斷，SD和TN與Chl-a之間的相關系數較高，因此選取SD和TN為解釋變量代入GAM模型；同時，SD和TN間的相關系數為?0.51，可能存在共曲線性，因此需剔除1個變量。

根據變量預分析結果，將SD和TN分別作為解釋變量構建GAM模型，方程如下：

式中：

x1——都昌點位SD（m）；

x2表示都昌點位 TN 質量濃度（mg·L?1）。基于表1結果可得，解釋變量SD的顯著性和解釋率均較好，擬合程度較高，因此選取SD的模型為最優模型。

表1 GAM模型分析結果Table 1 Analysis results of GAM model

從圖4中可以看出，都昌點位Chl-a與SD存在較強的非線性關系，當SD為0.2—0.6 m時，Chl-a質量濃度基本保持不變，維持在相對較低水平；當SD為0.6—0.8 m時，Chl-a質量濃度增加明顯，而當介于0.8—1.2 m時，Chl-a質量濃度逐漸減少；但當SD大于1.2 m時，Chl-a質量濃度又顯著增加。

圖4 都昌SD與Chl-a質量濃度的關系Fig. 4 The relationship between SD and Chl-a in Duchang

2.3.3 蛇山葉綠素a與環境因子的相關性

對于蛇山點位，根據Chl-a與環境因子間的相關系數判斷，TN和CODMn與Chl-a之間的相關系數較高，因此選取 TN和 CODMn為解釋變量代入GAM模型；同時，TN和CODMn間的相關系數為?0.30，故排除共曲線性可能，無需剔除變量。根據變量預分析結果，將TN和CODMn代入如下模型方程擬合：

式中：

x1——蛇山點位 TN 質量濃度（mg·L?1）；

x2——蛇山點位 CODMn質量濃度（mg·L?1）。得到r2為0.627，TN和CODMn的自由度分別為5.05和5.64，解釋率為71.0%，解釋率較高。

從圖5中可以看出，蛇山點位Chl-a與TN和CODMn存在較強非線性關系，CODMn的非線性影響能力更強，Chl-a與CODMn的關系更為復雜。當TN質量濃度為0—1 mg·L?1時，TN對Chl-a的影響最為顯著，兩者呈現顯明顯的負相關關系，當TN質量濃度大于1 mg·L?1時，Chl-a質量濃度隨著TN的增加而保持相對的穩定，后出現略微減小趨勢。對于Chl-a與CODMn而言，當CODMn質量濃度介于0—2 mg·L?1時，Chl-a保持相對穩定；當 CODMn質量濃度為 2.0—2.4 mg·L?1時，Chl-a質量濃度略有增長；當CODMn質量濃度大于2.4 mg·L?1時，Chla與其呈現明顯的負相關關系。

圖5 蛇山TN和CODMn與Chl-a質量濃度的關系Fig. 5 The relationship between TN, CODMn and Chl-a in Sheshan

2.4 討論

鄱陽湖3個研究點位藻類增殖與水質因子的關系存在較為顯著的空間差異，位于西部的蚌湖、中部的都昌和東部的蛇山點位Chl-a質量濃度分別與TN和CODMn、SD及TN和CODMn呈現非線性關系，所得的對應解釋率分別為 52.3%、45.1%和71.0%。更多地，蚌湖與蛇山點位由于地理位置的差異，Chl-a質量濃度的伴隨變化也不盡相同，需要給出不同點位（區域）的限制性水質指標與主要的可能因素，并探討相應的有效解決途徑。

蚌湖屬鄱陽湖營養鹽質量濃度較高的點位，通過GAM模型模擬得出Chl-a與TN和CODMn存在非線性關系，TN的非線性影響能力更強，Chl-a與TN的關系更為復雜，其并未呈現隨TN和CODMn質量濃度增加而增加的趨勢。蚌湖點位的Chl-a質量濃度主要受到氮含量與有機物含量的限制，對于該點藍藻水華防控而言，應當加強其TN的監控，并盡量保持相對較低水平，雖出現TN質量濃度增加而Chl-a質量濃度減少的階段，但其可能原因是非適應性的藻類細胞死亡，而伴隨著的是較強適應性的藻種成為優勢藻，藻類組成發生變化（朱偉等，2008），因此會伴隨該藻種的增殖與TN呈現明顯的正相關關系，Chl-a質量濃度顯著升高的現象；同時，也應當重點關注蚌湖點位的CODMn，盡量將有機物含量控制在相對較小的水平，盡管出現CODMn質量濃度增加而Chl-a質量濃度減小的現象，但湖泊中過多的有機物含量會間接影響水體溶解氧和微生物等因子，對水體健康不利（石玉等，2021）。

通過GAM模型擬合得出，都昌點位Chl-a與SD呈現較強且復雜的非線性關系，原因可能在于：鄱陽湖河湖相通，泥沙輸入、風浪攪動、航運采砂等眾多因素都會導致水體透明度下降（程時長等，1993；賴普文等，2015；李海軍，2017；劉同宦等，2020）。伴隨著藍藻等浮游生物的繁殖，水體中的有機物質被利用，同時浮游生物的繁殖還可能減小風浪，降低泥沙攪動，具體表現為，在一定范圍內，水體的SD與Chl-a含量之間的正相關關系。對于都昌點位的Chl-a調控，應當主要關注透明度指標，使其保持較為合理的區間，過高和過低的水體透明度對于湖泊生態環境均有一定程度的影響；同時，也應當加強都昌點位的營養鹽因子調控，雖然本研究中氮磷等營養物質對于都昌點位并非關鍵性限制因子，但是TN/TP也可能是潛在的限制因子。

通過GAM模型擬合得出，蛇山點位Chl-a與TN和CODMn存在較強非線性關系，CODMn的非線性影響能力更強，Chl-a與CODMn的關系更為復雜，與蚌湖點位類似，蛇山也未呈現隨TN和CODMn增加而增加的趨勢。蛇山點位的Chl-a質量濃度主要受到有機物含量與氮含量的限制，與蚌湖點位不同，整體而言Chl-a質量濃度伴隨TN和CODMn質量濃度變化的幅度不大。對于蛇山點位的Chl-a調控和藍藻水華防控應當充分注重水質指標的綜合控制，保持水體氮含量與有機物含量于合理區間內，維護該區域的水生態平衡。

本文基于GAM模型對鄱陽湖Chl-a與水質因子進行了相關性分析，但也存在一定不足：蚌湖和都昌點位的最優模型解釋率中等，分別為52.3%和45.1%，原因可能為缺乏考慮風場和流場等環境因子對于Chl-a質量濃度的影響。在對鄱陽湖風場和湖流的綜合研究結果顯示，鄱陽湖中部和東部湖區在風場和湖流形態的雙重作用下，其水流呈現“環流”現象，這也導致了水體更新時間相較鄱陽湖其他區域更長（Yuan et al.，2020），營養鹽與有機物在此范圍內會停留更長的時間，也可能是導致該區域近年來藍藻水華現象頻發的主要原因，過多的氮含量與有機物含量無法被利用或降解，伴隨 Chl-a質量濃度與環境因子的逐漸惡化，亦可能是導致鄱陽湖水環境惡化的重要原因（Zhang et al.，2015）。可在今后的研究中加強對于鄱陽湖典型點位的風場和流場的連續監測，分析考慮綜合環境因子（水文過程、氣象因素和水質因子）對于鄱陽湖營養化的影響。

3 結論

鄱陽湖Chl-a的關鍵性水質指標與主要的可能因素的判別，對于永葆鄱陽湖“一湖清水”和水環境永續可持續發展具有重要理論和實踐意義，本文的主要結論如下：

（1）鄱陽湖3個研究點位（西部蚌湖、中部都昌和東部蛇山）Chl-a與水質因子的關系存在較為顯著的空間差異，其 Chl-a質量濃度分別與 TN和CODMn、SD及TN和CODMn呈現非線性關系，模型解釋率分別為52.3%、45.1%和71.0%；蚌湖與蛇山點位由于地理位置的差異，Chl-a質量濃度的伴隨變化也不盡相同。

（2）運用GAM模型擬合得出，蚌湖Chl-a與TN和CODMn存在非線性關系，TN的非線性影響能力更強，Chl-a與TN的關系更為復雜，蚌湖點位的Chla質量濃度主要受到氮含量與有機物含量的限制，對于該點藍藻水華防控應注意控制較低的氮含量和合理范圍的有機物含量；都昌點位Chl-a與SD呈現較強且復雜的非線性關系，對于都昌點位的 Chl-a控制，應當主要關注透明度指標，但TN/TP也可能是潛在的限制因子；蛇山點位Chl-a與TN和CODMn存在較強非線性關系，CODMn的非線性影響能力更強，Chl-a與CODMn的關系更為復雜，對于蛇山點位的藻類或藍藻水華防控，需要充分注重水質指標的綜合控制，不僅僅局限于營養鹽與有機物含量。

（3）鄱陽湖不同區域Chl-a質量濃度與各水質因子的相關程度不盡相同，區別于太湖等湖泊，鄱陽湖因其特殊的性質，具有復雜的風場和湖流形態特征，在今后的研究中應當加強對于鄱陽湖風場和流場的連續監測，分析考慮綜合環境因子（水文過程、氣象因素和水質因子等）對于鄱陽湖營養化的影響。