吳 娟，朱躍龍，金 松，楊 濤，馮 鈞，吳志勇，薛 濤，姜悅美

(1.太湖流域管理局水文局(信息中心)，上海 200434； 2.河海大學水文水資源學院，江蘇 南京 210098)

太湖(119° 52′ 32″E～120° 36′ 10″ E、30° 55′ 40″N～31°32′58″N)是我國第三大淡水湖，水面面積約2 338 km2，是典型的大型淺水湖泊。近幾十年來，太湖藍藻水華(以下簡稱“藻華”)呈現(xiàn)非線性、復雜性等特點[1]。藻華不完全是本地藍藻短時快速生長暴發(fā)所致[2]，而是已經(jīng)成為優(yōu)勢種群的藍藻在合適條件下，受物理(溫度[3]、光照[4]、水動力[5]等)、化學(氮、磷[6]等)[7]、生物(與其他藻類的競爭)[8]因素等共同影響，在較短的時間里改變水平與垂直位置形成“暴發(fā)”[9]，因此，藻華本質(zhì)是藻類生物量在水體中逐漸增加、上浮、聚集、遷移的過程[10]。藻華擴張分為時間擴張、空間擴張和生物量擴張。時間擴張用藻華持續(xù)時間、發(fā)生時間等表征，空間擴張用基于遙感影像的藻華面積表征，生物量擴張用葉綠素a、浮游植物總生物量等表征[11]。由于生物量統(tǒng)計主要以特定湖區(qū)采樣監(jiān)測、計數(shù)鑒定方法為主，受制于采樣站點的位置與頻次，獲取資料困難，因此，本文以研究藻華時間、空間擴張為主。

研究表明，氣候變化、水體富營養(yǎng)化是太湖藻華增多的重要原因[12]。大面積偶發(fā)性太湖藻華受氣象因子影響較大，而大面積頻發(fā)性藻華主要受營養(yǎng)鹽空間分布影響[13]。太湖處于較高的營養(yǎng)鹽水平，超過藻華發(fā)生的閾值[14]，但季節(jié)性營養(yǎng)鹽抑制過程[15]對藻華強度有一定的限制作用[16]，足以維持水華持續(xù)發(fā)生[17]。杭鑫等[18]認為氣溫、風、降水等氣象因子對太湖藻華影響較大，適宜的氣象因子為藻華暴發(fā)提供了有利條件。秦伯強等[19]研究表明：藻華面積形成過程較復雜，藍藻細胞生長階段，營養(yǎng)鹽濃度、溫度、光照等因素為藻華提供了物質(zhì)基礎(chǔ)；在藍藻水華暴發(fā)階段，藻華面積空間分布則受藍藻細胞團浮力作用與水動力湍流作用的共同影響，風速決定藻華垂向分布，風向和風速共同決定其水平分布。

藍藻生態(tài)動力學模型因過程復雜、參數(shù)眾多[20]，導致應(yīng)用在預(yù)測中存在一定的困難，而機器學習算法(如神經(jīng)網(wǎng)絡(luò)等)為太湖藻華預(yù)測提供了新的理論和方法[21]。機器學習是一門研究計算機模擬或?qū)崿F(xiàn)人類活動的學科[22]，也是知識發(fā)現(xiàn)、數(shù)據(jù)挖掘的重要基礎(chǔ)[23]。通過挖掘數(shù)據(jù)背后的深度價值和內(nèi)在聯(lián)系[24]，機器學習具有高度的非線性與靈活性[25]，應(yīng)用于短期降水預(yù)報[26]、徑流預(yù)測[27]、藻華預(yù)警[28]，有效提升了預(yù)報精度和效率。張艷會等[29]采用BP人工神經(jīng)網(wǎng)絡(luò)和模糊理論，建立了藻華發(fā)生的模糊風險評價方法；于家斌等[30]提出基于長短期記憶(LSTM)循環(huán)神經(jīng)網(wǎng)絡(luò)(RNN)藻華預(yù)測模型(LSTM模型預(yù)測精度較高,對樣本具有較好的適應(yīng)性)；羅曉春等[31]采用隨機森林機器學習算法分析同期氣象因子與藻華綜合指數(shù)的關(guān)系，定量評估影響藻華的主要氣象因子的貢獻率。上述研究為太湖藻華的預(yù)測預(yù)警與風險管理提供了參考。

太湖湖面寬廣，一般采用野外觀測和數(shù)值模擬研究湖流導致藻華遷移和堆積。風向分布具有不均勻性，并隨著風速減小差異有所擴大。在低風速條件下，無法形成穩(wěn)定環(huán)流，藍藻容易在太湖表面積聚，而此時各湖區(qū)之間風向差異較大，需要借助藍藻生態(tài)動力學模型深入研究。本研究側(cè)重采用全太湖、分湖區(qū)建模的思路，基于支持向量機(SVM)、長短記憶神經(jīng)網(wǎng)絡(luò)(LSTM)、極端梯度提升樹(XGBoost)模型，快速預(yù)判各湖區(qū)藻華面積，旨在為太湖藻華預(yù)測預(yù)警、風險管理提供技術(shù)支撐。

1 數(shù)據(jù)與方法

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

按照岸線、水質(zhì)和水下地形等因素，太湖可分為竺山湖、梅梁湖、貢湖等[32]，見圖1。近年來，中分辨率成像光譜儀(MODIS)數(shù)據(jù)因其光譜分辨率高、觀測周期短等特點，廣泛應(yīng)用于太湖藻華動態(tài)監(jiān)測[33]。考慮到太湖藻華每年4月開始零星暴發(fā)，12月以后基本無大范圍藻華暴發(fā)，因此，選取4—11月的 EOS /MODIS 影像數(shù)據(jù)(排除了云層、天氣質(zhì)量等對影像質(zhì)量影響較大的數(shù)據(jù))共761景，時間間隔從1 d到15 d不等，藻華面積來源于中國科學院南京地理與湖泊研究所解譯成果。本研究采用的藻華面積為絕對面積，指當像元內(nèi)存在水華時，只計算像元內(nèi)水華覆蓋部分的面積。太湖內(nèi)共布設(shè)33個監(jiān)測點，覆蓋了全部水域，監(jiān)測項目包括pH、總磷、總氮、高錳酸鹽指數(shù)、氨氮、葉綠素a、藍藻密度，僅每月上旬監(jiān)測1次，根本無法滿足藻華面積建模要求，因此本研究只能采用環(huán)湖10個水質(zhì)站、4個氣象站時間分布較精細的數(shù)據(jù)進行建模。水質(zhì)指標包括水溫、pH、電導率、溶解氧、氨氮、總磷、總氮、濁度、高錳酸鹽指數(shù)，氣象指標包括氣溫、風速、風向、降水。本研究的資料序列均為2014—2018年。由于東部沿岸區(qū)與東太湖沒有MODIS藻華資料，而貢湖、南部沿岸區(qū)均有水源地，因此，本研究采用全太湖、貢湖、南部沿岸區(qū)、中西北湖區(qū)(西部沿岸區(qū)+湖心區(qū)+竺山湖+梅梁湖)分別建模。南部沿岸區(qū)代表站包括幻溇、濮溇、湯溇，中西北湖區(qū)代表站包括長豐澗、大港、洪巷、官瀆、社瀆。

圖1 太湖分區(qū)與觀測站示意圖Fig.1 Sub-regions in Taihu Lake and observation stations

1.2 研究方法

太湖藻華受物理、化學、生物因素等共同影響，機理十分復雜，呈非線性變化特征。本研究需要先采用Pearson 相關(guān)系數(shù)去除不相關(guān)或重復變量，提取對各湖區(qū)藻華面積影響較大的主要氣象水文水質(zhì)因子。從時間、空間尺度上考慮到不同湖區(qū)水質(zhì)差異較大，MODIS成像時間為上午10：30左右，因此，從時間尺度上分成3個模型，模型1采用當天0：00—8：00之前的氣象、水文、水質(zhì)數(shù)據(jù)預(yù)測當天的藻華面積，模型2采用提前一天0：00—23：00的數(shù)據(jù)預(yù)測當天的藻華面積，模型3采用提前兩天0：00—23：00的數(shù)據(jù)預(yù)測當天的藻華面積。由于SVM、LSTM、XGBoost模型既可解決分類問題，也可解決回歸問題，并且處理非線性問題有一定的優(yōu)越性，在水質(zhì)評價與預(yù)測研究中有一定的應(yīng)用情景，因此，本研究分別開展兩類算法的應(yīng)用研究，共建立24個模型。訓練集、驗證集為2014—2018年序列數(shù)據(jù)，其中隨機抽取70%的樣本(533組)用于訓練，剩下30%的樣本(228組)作為驗證集。

1.2.1 數(shù)據(jù)預(yù)處理

由于輸入變量可能不在同一個數(shù)量級上，數(shù)據(jù)之間的差異性可能對模型學習能力存在一定的不利影響，為了減少數(shù)值差異、保證模型參數(shù)穩(wěn)定收斂，需要進行歸一化處理，確保各變量處于相同的量級。

(1)

式中：X——輸入變量；X′——變量X歸一化后的值；Xmin、Xmax——X的最小值、最大值。

1.2.2SVM

SVM是通過核函數(shù)，將低維輸入空間中線性不可分的點映射成高維特征空間中線性可分的點，通過劃分超平面，使所有的點到分類超平面的距離最大化[34]。其中，高維空間中距離分類超平面最近的那些點所對應(yīng)的低維空間點被稱為支持向量[35]。核函數(shù)的作用是接受兩個低維空間中的輸入向量，計算出經(jīng)過某種變換后二者在高維空間中的向量內(nèi)積值，具體見文獻[36]。

常用的核函數(shù)有線性、多項式、徑向基、sigmoid核函數(shù)。經(jīng)對比，本研究采用誤差較小、分類準確率較高的徑向基函數(shù)作為核函數(shù)[37]。SVM采用徑向基核函數(shù)，相應(yīng)參數(shù)為懲罰因子c和核參數(shù)g。為確定最佳的參數(shù)值，將c、g分別取以2為底的指數(shù)離散值，代入V折交叉驗證法中，選取平均誤差最小的c、g值；采用網(wǎng)格搜索法，減小取值范圍和間隔，確定最佳參數(shù)，根據(jù)計算，最佳參數(shù)為c=0.76，g=0.38，V=5。

1.2.3LSTM

RNN屬于深度學習方法[38]，用于處理非線性時間序列，但在實際建模中，存在梯度消失及梯度爆炸問題[39]。為此，Hochreiter等提出了LSTM，將RNN中隱含層中的神經(jīng)元替換為記憶體，實現(xiàn)信息保留與長期記憶。通過記憶細胞進行狀態(tài)信息存儲，門結(jié)構(gòu)負責細胞狀態(tài)的更新與保持[40]。每個記憶體包含一到多個記憶細胞和3種“門”，包括“遺忘門”、“輸入門”和“輸出門”。 LSTM訓練過程是尋找最優(yōu)參數(shù)，使模型收斂、損失函數(shù)達到最小[41]。算法見文獻[42]。LSTM 網(wǎng)絡(luò)層中隱藏神經(jīng)元數(shù)目是影響準確率的重要參數(shù)，為尋求模型的全局最優(yōu)解，采用基于梯度下降、具有相當魯棒性的ADAM自適應(yīng)學習率優(yōu)化算法；經(jīng)計算，初始學習率取0.001，隱藏神經(jīng)元數(shù)為256，可使全局最優(yōu)。

1.2.4XGBoost

XGBoost模型為樹結(jié)構(gòu)，是對梯度提升決策樹算法的改進，采用加法學習模型進行優(yōu)化[43]。通過不斷迭代生成新的樹，加入子樹使模型不斷逼近樣本分布，將多個分類準確率較低的樹模型組合成一個準確率較高的模型[44]。通過在Xgboost損失函數(shù)中加入正則項，既控制模型的復雜度，也降低模型的方差，避免模型過擬合[45]。算法見文獻[46]。XGBoost模型主要參數(shù)包括eta、subsample、max_depth。eta指學習速率，可以提高模型的魯棒性；subsample指訓練的實例樣本占總體樣本的比例；max_depth為樹的最大深度，用來避免過擬合。經(jīng)計算，最佳參數(shù)為eta=0.15，sub-sample=0.2，max_depth=8。

1.3 模型評價

常用評價指標包括：納什效率系數(shù)(NSE)、相對偏差(PB)，公式如下：

(2)

(3)

NSE取值范圍為負無窮至1，越接近1表示模擬效果越好、可信度越高；NSE接近0，表示模擬結(jié)果接近觀測值的平均值水平，即總體結(jié)果可信，但過程模擬誤差大。PB絕對值越小，模擬效果越好。

Pearson 相關(guān)系數(shù)(r)是描述兩個隨機變量線性相關(guān)的統(tǒng)計量，計算公式如下：

(4)

r取值范圍為-1.0～1.0，越接近1.0，正相關(guān)越顯著，越接近-1.0，負相關(guān)性越顯著。構(gòu)造統(tǒng)計量t檢驗相關(guān)性是否顯著，計算公式如下：

(5)

統(tǒng)計量t遵從自由度ν=n-2的t分布。給定顯著性水平α，查t分布表，若t>tα，認為相關(guān)系數(shù)是顯著的。

2 結(jié)果與分析

2.1 相關(guān)性分析

為了確定所選氣象、水文、水質(zhì)因子對太湖藻華面積模擬的有效性，提高藻華面積模型的準確性，采用Pearson相關(guān)系數(shù)去除相關(guān)性不顯著的變量。分別計算當天、提前1 d到提前30 d共計31種時間尺度條件下，全太湖、貢湖、南部沿岸區(qū)、中西北湖區(qū)藻華面積，以及氣象、水文、水質(zhì)數(shù)據(jù)之間的Pearson相關(guān)系數(shù)，發(fā)現(xiàn)當天、提前1 d、提前2 d時間尺度相關(guān)系數(shù)明顯高于其他時段。Pearson相關(guān)系數(shù)較顯著的因子如表1所示，全太湖藻華面積與水溫、風速、降水等氣象因子相關(guān)性較高，其中，與當天無錫2：00、8：00風速負相關(guān)性較顯著，分別達到-0.17、-0.23。不同時間尺度下的貢湖藻華面積與水溫、濁度、總氮、氨氮正相關(guān)性均較高，相關(guān)系數(shù)均可達0.15以上。南部沿岸區(qū)藻華面積與水溫呈顯著正相關(guān)性，相關(guān)系數(shù)均可達0.19以上，與當日東山站2：00、8：00風速呈顯著負相關(guān)性。中西北湖區(qū)藻華面積與當日宜興站2：00、8：00風速的負相關(guān)性較顯著。

表1 藻華面積與氣象、水文、水質(zhì)因子之間的Pearson相關(guān)系數(shù)Table 1 Pearson relative coefficient between cyanobacterial bloom area and factors including meteorological, hydrological and water quality factors

2.2 模型分析結(jié)果

將訓練集和驗證集分別代入訓練好的全太湖、貢湖、南部沿岸區(qū)、中西北湖區(qū)藻華面積模擬模型中，并對輸出值進行反歸一化處理。模擬值與觀測值之間的NSE與PB見表2。對于不同時間尺度(當天、提前1 d、提前2 d)下的全太湖藻華面積回歸模型訓練集與驗證集而言，XGBoost回歸模型NSE最高、PB最低，其次是LSTM回歸模型、SVM回歸模型；驗證期的XGBoost回歸模型，當天的模型1與提前1 d的模型2，相對優(yōu)于提前2 d的模型3。

表2 3種機器學習回歸模型評價結(jié)果Table 2 Evaluation results of three machine-learning regression methods

對于不同時間尺度下貢湖、南部沿岸區(qū)、中西北湖區(qū)的藻華面積回歸模型NSE，XGBoost模型訓練集與驗證集NSE普遍較高，訓練集LSTM模型的NSE高于SVM模型，驗證集SVM模型的NSE高于LSTM模型；驗證期的貢湖SVM模型：藻華面積模型1、南部沿岸區(qū)XGBoost藻華面積模型1、中西北湖區(qū)XGBoost藻華面積模型2的NSE高于相應(yīng)湖區(qū)的其他時間尺度模型，與全太湖藻華面積回歸模型表現(xiàn)也基本一致。

《太湖藍藻水華評價方法(試行)》將藍藻水華面積分為4類：<240 km2、[240，600)km2、[600，1 000]km2、>1 000 km2,相應(yīng)面積水域發(fā)生水華分別稱為零星湖區(qū)、局部湖區(qū)、區(qū)域、大范圍水華暴發(fā)[47]。以全太湖藻華面積過程線擬合效果為例，分析SVM、LSTM、XGBoost回歸模型模擬結(jié)果，見圖2～4。模擬結(jié)果表明：不同時間尺度下SVM、LSTM回歸模型訓練集的模擬值偏低，XGBoost回歸模型訓練集模擬值與實測值基本接近；盡管SVM、XGBoost回歸模型驗證集模擬值較實測值偏小，但有效模擬了零星湖區(qū)、局部湖區(qū)水華的發(fā)展趨勢；LSTM回歸模型驗證集的模擬值與實測值差異較大，不合理跳動較多，這也驗證了神經(jīng)網(wǎng)絡(luò)模型存在過學習、泛化能力差的固有缺陷。

圖2 基SVM于回歸模型的全太湖藻華面積訓練集與驗證集模擬Fig.2 Training set and verification set simulations of cyanobacterial bloom area with SVM regression model in Taihu Lake

圖3 基于LSTM回歸模型的全太湖藻華面積訓練集與驗證集模擬Fig.3 Training set and verification set simulations of cyanobacterial bloom area with LSTM regression model in Taihu Lake

圖4 基于XGBoost回歸模型的全太湖藻華面積訓練集與驗證集模擬Fig.4 Training set and verification set simulations of cyanobacterial bloom area with XGBoost regression model in Taihu Lake

從上述分析可知，回歸模型在驗證期模擬效果不太好，考慮到藻華面積預(yù)測與防控業(yè)務(wù)中主要關(guān)注藻華面積與相應(yīng)湖區(qū)面積比例的范圍，對精確數(shù)值關(guān)注較少。為進一步改進藻華面積模型，采用分類算法輸出離散的類別值。將藻華面積與相應(yīng)湖區(qū)面積比例按照<10%、[10%，25%)、[25%，50%)、>50%分成4類，分析不同時間尺度下SVM、LSTM、XGBoost分類模型的正確率，見表3。表3結(jié)果表明：全太湖各模型的訓練集正確率達到0.97以上，驗證集SVM、XGBoost分類模型正確率達到0.80以上，LSTM分類模型正確率略低，為0.66～0.86，誤差較大，見圖5。貢湖各模型的訓練集正確率達到0.89以上，驗證集正確率達到0.82以上；南部沿岸區(qū)各模型的訓練集正確率達到0.97以上，驗證集正確率達到0.80以上；中西北湖區(qū)訓練集正確率達到0.83以上，驗證集XGBoost模型正確率達到0.80以上，好于LSTM分類模型、SVM分類模型?？傮w而言，貢湖藻華面積SVM、LSTM、XGBoost分類模型驗證集正確率普遍優(yōu)于其他湖區(qū)，可能與貢湖風浪站相對于其他環(huán)湖站點位置更加偏向于貢湖內(nèi)部，氣象、水文、水質(zhì)數(shù)據(jù)代表性更強有關(guān)。

圖5 不同時間尺度下SVM、LSTM、XGBoost分類模型的藻華面積驗證集模擬Fig.5 Verification of cyanobacterial bloom areas with SVM, LSTM and XGBoost classification models in Taihu Lake

表3 3種機器學習分類模型評價結(jié)果Table 3 Evaluation results of three machine-learning classification methods

從上述研究可知，SVM、LSTM模型在藻華面積較低時，模擬與實測結(jié)果較接近；當藻華面積超過600 km2時，模擬與實測結(jié)果偏差較大?？赡茉蛴?個：一是太湖藻華面積高值相對較少，且水華從暴發(fā)到消退較快，在運用機器學習方法訓練模型時，如果目標值數(shù)量不均等、數(shù)值差異大，SVM、LSTM模型一般偏向于數(shù)量多、數(shù)值低的藻華面積的樣本。二是LSTM較適合大數(shù)據(jù)的機器學習，在數(shù)據(jù)量足夠多的情況下效果較優(yōu)，而本研究的數(shù)據(jù)量可能遠遠達不到大數(shù)據(jù)的要求，LSTM無法很好地保證模型的泛化能力，存在陷入局部最優(yōu)以及收斂時間長等問題[48]，導致模擬效果不佳。三是太湖風浪作用顯著，適度的間歇性風浪擾動可以促進氮磷等營養(yǎng)鹽釋放，刺激藍藻細胞生理特征改變，對太湖藻華遷移聚集影響較大，而模型只考慮了風場，流場數(shù)據(jù)的缺乏導致水動力輸移考慮可能不夠?？傮w而言，由于XGBoost模型在損失函數(shù)尋優(yōu)過程中加入了正則項來控制過擬合現(xiàn)象，模擬藻華面積與實測相比，精度較高、穩(wěn)健性較好。

3 結(jié) 論

a. XGBoost全太湖與分區(qū)的回歸模型訓練集與驗證集模擬效果均較好，LSTM回歸模型訓練集較好、但驗證集較差，SVM回歸模型介于XGBoost與LSTM回歸模型之間；SVM、XGBoost回歸模型的驗證集模擬值偏小，但有效模擬了藻華的發(fā)展趨勢；選取當天、提前1 d的氣象水文水質(zhì)因子作為太湖及分區(qū)藻華面積模型的輸入，模擬結(jié)果較好。

b. 在全太湖、中西北湖區(qū)藻華面積模擬中，XGBoost分類模型正確率較高；在貢湖、南部沿岸區(qū)藻華面積模擬中，SVM、LSTM、XGBoost分類模型正確率均較高。

c. SVM、LSTM模型在藻華面積較低時模擬與實測結(jié)果較接近，當藻華面積較高時模擬與實測偏差較大，可能與機器學習方法訓練模型一般偏向于數(shù)量多、數(shù)值低的藻華面積的樣本有關(guān)；XGBoost模型精度較高、穩(wěn)健性較好，適合在太湖藻華面積預(yù)測中推廣。