許陽(yáng)春，覃 苗，蘇玉萍*，林曉萍，蘇金洙

(1.福建師范大學(xué)環(huán)境科學(xué)與工程學(xué)院，福建 福州 350007；2.福建省污染控制與資源循環(huán)重點(diǎn)實(shí)驗(yàn)室，福建 福州 350007)

近年來(lái)，隨著我國(guó)經(jīng)濟(jì)的快速發(fā)展，生活及工業(yè)污水排放量增大，加重了水體的富營(yíng)養(yǎng)化[1]。水華是一種或多種浮游植物暴發(fā)性增殖的生態(tài)失衡現(xiàn)象。在我國(guó)和世界其他國(guó)家的多數(shù)湖泊和水庫(kù)中，常有關(guān)于淡水水華發(fā)生的報(bào)道[2-4]，水體富營(yíng)養(yǎng)化和水華暴發(fā)已成為全球性環(huán)境問(wèn)題[5]。大量營(yíng)養(yǎng)元素如氮和磷被認(rèn)為是促進(jìn)水華暴發(fā)的重要因素，研究表明，氮磷比值、水文和氣象條件也是誘導(dǎo)水華發(fā)生的主要因素[6]。Ndong等[7]構(gòu)建了基于風(fēng)向、氣溫、風(fēng)速的水華氣象指數(shù)模型，結(jié)果表明，在湖泊氮磷比滿足藍(lán)藻增殖的條件下，氣象條件對(duì)水華的發(fā)生產(chǎn)生較大的影響。水華發(fā)生機(jī)理的相關(guān)假說(shuō)也表明，在一定的水文、氣象條件下藻類會(huì)大量聚集[8-10]，進(jìn)一步說(shuō)明營(yíng)養(yǎng)鹽、微量元素、水文、氣象等條件都是影響水華發(fā)生的主要環(huán)境因子[11-14]。嚴(yán)重的水華會(huì)導(dǎo)致水體透明度降低，藻類死亡分解過(guò)程中還會(huì)消耗大量溶解氧，造成水體缺氧，影響水生生物的生活和生存。此外，部分有害水華藻類還會(huì)產(chǎn)生毒素和異味物質(zhì)，對(duì)水生態(tài)系統(tǒng)和人類健康造成極大威脅[15]。因此，為了有效治理水華，探討水華預(yù)警技術(shù)是當(dāng)務(wù)之急。水華預(yù)警模型是利用數(shù)學(xué)模型來(lái)描述藻類生長(zhǎng)與環(huán)境因子之間的關(guān)系，是研究水華變化過(guò)程較為有效的工具[16]。

國(guó)內(nèi)外學(xué)者已開(kāi)展了大量水華預(yù)警模型的研究。Dillon等根據(jù)湖泊中總磷和葉綠素濃度建立了第一個(gè)實(shí)證數(shù)據(jù)模型，對(duì)葉綠素濃度進(jìn)行預(yù)測(cè)[17]。Muhammetoglu等建立三維動(dòng)態(tài)水質(zhì)模型，該模型已成為測(cè)試湖泊富營(yíng)養(yǎng)化水平的重要工具[18]。隨著研究手段的快速發(fā)展和大量歷史數(shù)據(jù)的共享，更多的有效方法被用于建立數(shù)理模型，如多變量回歸模型等[19-20]。Shen等利用長(zhǎng)期觀測(cè)數(shù)據(jù)和支持向量機(jī)(LS-SVM)模型來(lái)模擬藻類增殖過(guò)程，結(jié)果表明，數(shù)據(jù)驅(qū)動(dòng)模型能夠模擬藻類實(shí)際增殖情況，具有良好預(yù)測(cè)能力；LS-SVM在藻類繁殖模擬中的成功應(yīng)用為有害藻華預(yù)測(cè)提供了新方法[21]。Recknagel等運(yùn)用日本霞浦湖12年的數(shù)據(jù)，以6種環(huán)境因子為自變量，葉綠素a(Chl-a)濃度為因變量，建立時(shí)間序列模型模擬葉綠素a濃度的變化[22]。此外，劉云翔等改進(jìn)CART決策樹(shù)，利用Fayyad邊界點(diǎn)判定定理，選出與水華發(fā)生相關(guān)性較大的影響因子，改進(jìn)后的水華預(yù)警模型能提高水華預(yù)測(cè)的正確率[23]。孔繁翔等利用水文和氣象監(jiān)測(cè)，結(jié)合衛(wèi)星遙感技術(shù)對(duì)太湖藍(lán)藻水華發(fā)生進(jìn)行預(yù)測(cè)，準(zhǔn)確率達(dá)80%[24]。張嬌等利用洱海藍(lán)藻水華發(fā)生的時(shí)空分布特征，結(jié)合水質(zhì)和氣象數(shù)據(jù)推測(cè)洱海藍(lán)藻水華暴發(fā)的主要影響因子，并探索基于多源遙感監(jiān)測(cè)數(shù)據(jù)的水華預(yù)警方法，其準(zhǔn)確率達(dá)89%[25]。

目前，神經(jīng)網(wǎng)絡(luò)已被廣泛應(yīng)用于水華的預(yù)測(cè)，在湖泊[26-27]、河流[28-29]和近海水域[30-31]中成效較好。BP(Back propagation)神經(jīng)網(wǎng)絡(luò)，是一種多層前饋型神經(jīng)網(wǎng)絡(luò)，主要特點(diǎn)是信號(hào)向前傳遞，誤差反向傳播，具有較高的構(gòu)建模型和數(shù)據(jù)擬合能力，非常適合處理非線性關(guān)系間內(nèi)部機(jī)制比較復(fù)雜的問(wèn)題[32]。目前基于BP神經(jīng)網(wǎng)絡(luò)的水華預(yù)警模型也日趨成熟。張娣等收集1986—2007年蘇州市吳縣的氣溫、氣壓、日照、相對(duì)濕度等氣象數(shù)據(jù)，選取其中8個(gè)年份的太湖藍(lán)藻暴發(fā)案例，分析了太湖藍(lán)藻暴發(fā)前一段時(shí)間的氣象條件，得出太湖藍(lán)藻水華暴發(fā)的主要限制因子為氣溫、降水、氣壓和相對(duì)濕度，以此為基礎(chǔ)建立的BP神經(jīng)網(wǎng)絡(luò)模型在太湖藍(lán)藻水華的模擬預(yù)測(cè)中取得較好的結(jié)果[32]。本研究以平潭三十六腳湖為研究對(duì)象，利用其2016年1月—2017年5月的歷史監(jiān)測(cè)數(shù)據(jù)進(jìn)行模型的預(yù)測(cè)研究，以期為平潭三十六腳湖的水華預(yù)警預(yù)報(bào)工作提供新的思路和科學(xué)依據(jù)。

1 材料和方法

1.1 模型演算數(shù)據(jù)來(lái)源

本文收集了三十六腳湖(圖1)自動(dòng)監(jiān)測(cè)站2016年1月至2017年5月期間記錄的3 500組水質(zhì)監(jiān)測(cè)數(shù)據(jù)及相應(yīng)時(shí)間的氣象數(shù)據(jù)(https：//rp5.ru)。水質(zhì)參數(shù)包括水溫(℃)、溶解氧(mg/L)、電導(dǎo)率(μs/cm)、pH(無(wú)量綱)、濁度(NTU)、總氮(mg/L)、總磷(mg/L)、氨氮(mg/L)和Chl-a(μg/L)。氣象參數(shù)包括平均氣溫(℃)、日照時(shí)數(shù)(h)、降水量(mm)、平均風(fēng)速(m/s)及風(fēng)向(°)。

1.2 BP神經(jīng)網(wǎng)絡(luò)演算模型

以水質(zhì)及氣象指標(biāo)的不同組合作為模型的輸入因子，以Chl-a作為輸出因子，對(duì)三十六腳湖2016年1月至2017年5月氣象及水質(zhì)歷史數(shù)據(jù)進(jìn)行模型的擬合演算及驗(yàn)證。

1.2.1 模型評(píng)價(jià)標(biāo)準(zhǔn)

為分析各種預(yù)測(cè)性能的優(yōu)劣，需建立統(tǒng)一、公正合理的評(píng)價(jià)體系。目前尚沒(méi)有比較完整的評(píng)價(jià)指標(biāo)可以使用直觀、統(tǒng)計(jì)的方法對(duì)預(yù)警演算模型的精度進(jìn)行評(píng)價(jià)。根據(jù)文獻(xiàn)[33]的研究，本文選取均方根誤差(RMSE)、標(biāo)準(zhǔn)偏差比(RSR)、可決系數(shù)(R2)作為性能評(píng)價(jià)指標(biāo)。

1)均方根誤差RMSE

(1)

2)標(biāo)準(zhǔn)偏差比RSR

(2)

3)可決系數(shù)R2

(3)

式(1)～(3)中：n為樣本數(shù)；Ti是第i個(gè)實(shí)測(cè)值；T′i是第i個(gè)預(yù)測(cè)值；RMSE反映了預(yù)測(cè)值與實(shí)測(cè)值的絕對(duì)偏離程度；RSR對(duì)數(shù)據(jù)的離散程度進(jìn)行估算，進(jìn)而評(píng)價(jià)模型精度；R2是對(duì)已制作好的演算模型進(jìn)行檢驗(yàn)，比較它們的預(yù)測(cè)結(jié)果與實(shí)際發(fā)生情況的吻合程度。

1.2.2 參數(shù)選取及模型演算

在模型演算之前，需對(duì)各個(gè)參數(shù)進(jìn)行篩選。因氨氮數(shù)值為區(qū)間范圍值，無(wú)法獲取點(diǎn)對(duì)點(diǎn)數(shù)據(jù)進(jìn)行模型演算，故將其剔除。作為營(yíng)養(yǎng)鹽參數(shù)代表的總氮、總磷與水華發(fā)生有著密切聯(lián)系，但總氮、總磷受外源影響較大，并且綜合前期研究結(jié)果發(fā)現(xiàn)，三十六腳湖湖區(qū)的總磷已經(jīng)達(dá)到了滿足湖區(qū)浮游植物生長(zhǎng)的較高值，在湖區(qū)水華頻發(fā)月份，該時(shí)間段的營(yíng)養(yǎng)鹽濃度不再是水華發(fā)生的主導(dǎo)因素，因而也將這兩個(gè)參數(shù)剔除。經(jīng)過(guò)篩選，最終選定6個(gè)水質(zhì)參數(shù)：水溫、溶解氧、電導(dǎo)率、pH、濁度、Chl-a，以及5個(gè)氣象參數(shù)：平均氣溫、日照時(shí)數(shù)、降水量、平均風(fēng)速、風(fēng)向，進(jìn)行模型運(yùn)算。其中設(shè)Chl-a為模型的輸出因子，即模型的預(yù)測(cè)結(jié)果，其余10個(gè)參數(shù)為輸入因子。

1.3 三十六腳湖BP神經(jīng)網(wǎng)絡(luò)演算模型驗(yàn)證

1.3.1 模型驗(yàn)證數(shù)據(jù)來(lái)源

模型驗(yàn)證數(shù)據(jù)來(lái)自三十六腳湖野外采樣數(shù)據(jù)，分別于2019年3月19日、4月3、13、18日及25日對(duì)三十六腳湖S1、S2、S3三個(gè)點(diǎn)位進(jìn)行了5次采樣(表1和圖2)。Chl-a濃度的實(shí)測(cè)值采用分光光度法[34]測(cè)定。

表1 庫(kù)區(qū)采樣點(diǎn)位圖情況

1.3.2 模型驗(yàn)證數(shù)據(jù)處理

為驗(yàn)證以水溫、電導(dǎo)率、平均氣溫3個(gè)指標(biāo)為輸入因子的模型演算效果，于2019年3月13日至2019年4月26日對(duì)平潭三十六腳湖進(jìn)行了5次野外采樣，將實(shí)驗(yàn)室檢測(cè)的Chl-a濃度與模型預(yù)測(cè)濃度進(jìn)行比對(duì)。

首先，對(duì)各個(gè)參數(shù)進(jìn)行歸一化處理：

(4)

Ti是第i個(gè)歷史數(shù)據(jù)樣本值，T′i是第i個(gè)歷史數(shù)據(jù)Ti的歸一化值，Tmin是歷史數(shù)據(jù)中的最小值，Tmax是歷史數(shù)據(jù)中的最大值。

1.3.3 模型驗(yàn)證過(guò)程

經(jīng)演算分析已經(jīng)確定模型的輸入因子組合為水溫、電導(dǎo)率、平均氣溫。利用MATLAB 2016b軟件進(jìn)行預(yù)警模型的運(yùn)行。將2016年1月至2017年5月這3個(gè)指標(biāo)及對(duì)應(yīng)的Chl-a的歷史數(shù)據(jù)在系統(tǒng)外部進(jìn)行歸一化處理，進(jìn)行擬合演算，獲取學(xué)習(xí)經(jīng)驗(yàn)，模擬輸入因子和輸出因子之間的關(guān)系。同樣，將2019年3月13日至2019年4月26日的水溫、電導(dǎo)率、平均氣溫組合作為輸入因子，在系統(tǒng)外部進(jìn)行歸一化處理，根據(jù)歷史數(shù)據(jù)的學(xué)習(xí)經(jīng)驗(yàn)運(yùn)行得到輸出因子的歸一化數(shù)據(jù)，即Chl-a濃度歸一化預(yù)測(cè)值。再進(jìn)行反向歸一化計(jì)算得到2019年3月13日至2019年4月26日的Chl-a濃度，進(jìn)而與野外水樣的實(shí)際檢測(cè)值進(jìn)行比對(duì)，確定模型預(yù)測(cè)精度。

2 結(jié)果與討論

2.1 模型演算與驗(yàn)證輸出結(jié)果

以篩選出來(lái)的輸入因子(5個(gè)水質(zhì)參數(shù)和5個(gè)氣象參數(shù))為基礎(chǔ)，變換不同參數(shù)組合作為輸入因子，在軟件MATLAB 2016b中進(jìn)行演算。表2選取了可決系數(shù)R2>0.5，即預(yù)測(cè)結(jié)果與實(shí)測(cè)值具有一定吻合程度的參數(shù)組合。而與日照時(shí)數(shù)和降水量相關(guān)的參數(shù)組合由于吻合程度過(guò)低，此處僅討論其余8個(gè)參數(shù)的參數(shù)組合作為輸入因子的模型演算結(jié)果。

表2 基于不同輸入因子組合的部分模型演算獲得的RMSE、RSR、R2結(jié)果

續(xù)表2

根據(jù)表2，在列舉出的誤差較小、擬合度較高的演算結(jié)果中，首先剔除RMSE及RSR值較大的組合，再以擬合度對(duì)演算結(jié)果進(jìn)行精度綜合評(píng)價(jià)。綜合所有演算結(jié)果，組合8與組合9的RMSE、RSR值幾乎一致，R2同為0.97，擬合度較好。但在增加輸入因子的條件下，模型的演算精度并未有大幅度的提高，再綜合考慮模型實(shí)際應(yīng)用的成本，應(yīng)在保證預(yù)測(cè)精度的前提下，盡量減少輸入因子，故最終確定水溫、電導(dǎo)率、平均氣溫(組合8)為模型的輸入因子。當(dāng)水溫、電導(dǎo)率、平均氣溫組合作為模型的輸入因子時(shí)，模型演算輸出的Chl-a結(jié)果如圖3，結(jié)果顯示實(shí)測(cè)值和預(yù)測(cè)值在所有數(shù)據(jù)點(diǎn)上十分接近，呈現(xiàn)很好的擬合結(jié)果。

利用BP神經(jīng)網(wǎng)絡(luò)預(yù)警模型演算獲得的2019年3月13日至4月26日每日間隔3 h的Chl-a濃度(圖4)。水樣的采集時(shí)間分別為3月19日17時(shí)、4月3日11時(shí)、4月13日11時(shí)、4月18日17時(shí)及4月25日11時(shí)，對(duì)應(yīng)的模型演算Chl-a濃度分別為13.31、13.81、14.96、16.70和20.12 μg/L。這5個(gè)模型預(yù)測(cè)的Chl-a濃度將與實(shí)測(cè)值的Chl-a濃度進(jìn)行比對(duì)，以確定模型的預(yù)測(cè)精度。

2.2 葉綠素a實(shí)測(cè)值

圖5為三個(gè)取樣點(diǎn)位在不同時(shí)間點(diǎn)的Chl-a濃度實(shí)測(cè)值。由于模型演算使用的水溫(℃)、電導(dǎo)率(μs/cm)和Chl-a(μg/L)濃度數(shù)據(jù)來(lái)源于2016年1月至2017年5月期間的三十六腳湖S3點(diǎn)位(在線監(jiān)控點(diǎn))。因此，選取實(shí)測(cè)的S3點(diǎn)位水樣的Chl-a濃度與模型演算的Chl-a濃度進(jìn)行比對(duì)。S3點(diǎn)位在這5日的Chl-a實(shí)測(cè)值分別為12.61、13.22、14.54、17.28和19.52 μg/L(圖5)。

2.3 模型輸出結(jié)果(預(yù)測(cè)值)與實(shí)測(cè)值誤差分析

表3整理了Chl-a濃度的實(shí)測(cè)值和模型演算的預(yù)測(cè)值比較。除4月18日外，其余日期模型預(yù)測(cè)值雖高于實(shí)測(cè)值，但很接近，表明模型能較好地為水華暴發(fā)前期的應(yīng)急工作提供一定的數(shù)據(jù)參考。模型使用的氣象數(shù)據(jù)為每間隔3 h記錄一次，即每日23時(shí)、2時(shí)、5時(shí)、8時(shí)、11時(shí)、14時(shí)、17時(shí)、20時(shí)，水質(zhì)數(shù)據(jù)為每隔4 h記錄一次，即每日0時(shí)、4時(shí)、8時(shí)、12時(shí)、16時(shí)、20時(shí)。時(shí)間點(diǎn)不完全一致，為使數(shù)據(jù)在時(shí)間點(diǎn)上盡量保持一致性，秉承時(shí)間點(diǎn)就近原則剔除了每日2時(shí)、14時(shí)的氣象數(shù)據(jù)，這樣每日有6個(gè)時(shí)間點(diǎn)的數(shù)據(jù)。因數(shù)據(jù)處理過(guò)程已有一定的偏差，故對(duì)模型預(yù)測(cè)結(jié)果的精度會(huì)產(chǎn)生一定的影響。

通過(guò)各評(píng)價(jià)指標(biāo)公式對(duì)模型的預(yù)測(cè)精度進(jìn)行判斷，計(jì)算得到模型誤差和擬合度。2019年5次野外水樣的實(shí)測(cè)值與模型預(yù)測(cè)值的擬合度(R2)高達(dá)0.97，標(biāo)準(zhǔn)偏差比(RSR)為0.24。說(shuō)明以氣溫、電導(dǎo)率、水溫3個(gè)指標(biāo)為輸入因子的BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型的Chl-a濃度預(yù)測(cè)值與實(shí)測(cè)值的擬合度較高，偏離程度較小，模型預(yù)測(cè)效果較好，達(dá)到了模型運(yùn)算期望。

表3 Chl-a濃度實(shí)測(cè)值與模型預(yù)測(cè)值比較

3 結(jié)論

本文運(yùn)用BP神經(jīng)網(wǎng)絡(luò)模型，通過(guò)組合不同水質(zhì)和氣象參數(shù)，比較模型的演算效果，確定水溫、電導(dǎo)率、平均氣溫組合為最佳輸入組合。進(jìn)一步以水溫、電導(dǎo)率、平均氣溫為輸入因子演算Chl-a濃度預(yù)測(cè)值，同時(shí)，用Chl-a濃度的實(shí)測(cè)值驗(yàn)證了BP神經(jīng)網(wǎng)絡(luò)演算預(yù)測(cè)的效果。以擬合度(R2)、標(biāo)準(zhǔn)偏差比(RSR)、均方根誤差(RMSE)三個(gè)指標(biāo)對(duì)模型進(jìn)行精度評(píng)價(jià)，發(fā)現(xiàn)模型擬合度(R2)高達(dá)0.97，標(biāo)準(zhǔn)偏差比(RSR)為0.24。因此，應(yīng)用BP神經(jīng)網(wǎng)絡(luò)，以水溫、電導(dǎo)率、平均氣溫為輸入因子演算三十六腳湖葉綠素a濃度具有可行性，該方法有望進(jìn)一步用于葉綠素a濃度和藍(lán)藻水華的預(yù)警。