李 未，朱夢圓，王裕成，朱廣偉

1. 中國科學(xué)院南京地理與湖泊研究所，湖泊與環(huán)境國家重點實驗室，江蘇 南京 210008

2. 杭州市生態(tài)環(huán)境局淳安分局，浙江 杭州 311700

千島湖又名新安江水庫，是長三角最大的戰(zhàn)略水源地及杭州市的飲用水源，以水質(zhì)優(yōu)異而著稱[1-2]. 千島湖設(shè)有4個國控監(jiān)測斷面，其中，街口斷面位于安徽省和浙江省交界處，水質(zhì)考核標(biāo)準(zhǔn)為GB 3838-2002《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》河道Ⅱ類；庫體內(nèi)有3個監(jiān)測斷面，分別為小金山、三潭島和大壩前，水質(zhì)考核標(biāo)準(zhǔn)均為GB 3838-2002湖庫Ⅰ類(除小金山TP濃度≤0.025 mg/L外). 近年來，隨著流域經(jīng)濟的發(fā)展，千島湖承接的外源污染負(fù)荷增大，加上水文氣象過程的影響，庫區(qū)水質(zhì)出現(xiàn)不穩(wěn)定性，局部地區(qū)開始出現(xiàn)富營養(yǎng)化和藻類異常增殖現(xiàn)象[3-7]. 尤其是小金山斷面，地處自上游狹長河道往開敞型水面過渡的區(qū)域，極易受到上游河流斷面的影響，導(dǎo)致水質(zhì)狀況無法穩(wěn)定達(dá)標(biāo)，成為千島湖面臨的水環(huán)境保護難題[8-9].

為了快速捕捉水質(zhì)與藻類變化，自2011年起千島湖開始在國控監(jiān)測斷面建設(shè)集水溫、pH、溶解氧濃度、濁度、葉綠素a濃度等參數(shù)于一體的剖面高頻監(jiān)測系統(tǒng). 在獲取高頻時間序列對水質(zhì)狀況進行實時自動監(jiān)測的同時，嘗試用科學(xué)手段挖掘數(shù)據(jù)特征和規(guī)律，在千島湖沉降通量、熱成層和溶解氧分層的穩(wěn)定性以及藻類垂向分布規(guī)律方面取得了新進展[9-13]. 但是，目前通過高頻時間序列實現(xiàn)對千島湖水質(zhì)水華的短期預(yù)測，從而對藻類異常增殖帶來的水質(zhì)風(fēng)險做出預(yù)判的研究較少.

大數(shù)據(jù)及人工智能技術(shù)使得以數(shù)據(jù)驅(qū)動的建模方式對時間序列進行短期預(yù)測成為可能[14]，尤其是人工神經(jīng)網(wǎng)絡(luò)(artificial neural network，ANN)模型，因其具有良好的非線性映射和自學(xué)習(xí)能力，在非平穩(wěn)時間序列的高精度預(yù)測方面具有獨特的優(yōu)勢. BP (back propagation，反向傳播)神經(jīng)網(wǎng)絡(luò)作為主流網(wǎng)絡(luò)之一，結(jié)構(gòu)簡單、映射能力強，已在醫(yī)學(xué)[15]、大氣[16-17]、水文[18]、土壤[19-20]、生物[21]、環(huán)境[22-23]等多學(xué)科領(lǐng)域廣泛應(yīng)用. 但是，BP神經(jīng)網(wǎng)絡(luò)屬于靜態(tài)神經(jīng)網(wǎng)絡(luò)，信息流動是由輸入層單向流動到輸出層，對于歷史輸入輸出缺乏記憶聯(lián)想，對時間序列時變特性的適應(yīng)能力不強，從而影響預(yù)測精度[24-25].

NARX (nonlinear auto-regressive with exogenous inputs，帶有外部輸入的非線性自回歸)神經(jīng)網(wǎng)絡(luò)是動態(tài)神經(jīng)網(wǎng)絡(luò)的一種，其輸入結(jié)果是當(dāng)前外部輸入和歷史輸出結(jié)果的非線性函數(shù)，可看作是有時延輸入的BP神經(jīng)網(wǎng)絡(luò)加上輸出到輸入的延時反饋連接. 由于網(wǎng)絡(luò)結(jié)構(gòu)中存在時延和反饋，因此其可反映系統(tǒng)的歷史狀態(tài)信息，是一種有記憶功能的神經(jīng)網(wǎng)絡(luò)，可以模擬時間序列長期動態(tài)特征[26-30]. 近年來，NARX神經(jīng)網(wǎng)絡(luò)在環(huán)境空氣質(zhì)量[31-32]、水文預(yù)報[33]、大壩變形[34]、地表沉降[35]、農(nóng)作物水分蒸騰[36]、股票擇時[37]等預(yù)測中均取得了較好的效果. 研究[27,34]表明，NARX神經(jīng)網(wǎng)絡(luò)可用于預(yù)測時間序列，且通常可以保留信息的時間是常規(guī)遞歸神經(jīng)網(wǎng)絡(luò)的2～3倍.

鑒于此，該研究以千島湖小金山國控監(jiān)測斷面的水質(zhì)自動監(jiān)測站葉綠素a (Chla)高頻監(jiān)測數(shù)據(jù)為研究對象，構(gòu)建基于NARX神經(jīng)網(wǎng)絡(luò)的藻類預(yù)測模型，分析該模型在千島湖Chla濃度時序變化預(yù)測的效果，探討最優(yōu)預(yù)見期，以期為構(gòu)建以數(shù)據(jù)驅(qū)動的千島湖水華監(jiān)測預(yù)警系統(tǒng)提供科學(xué)依據(jù).

1 材料與方法

1.1 研究點概況與浮標(biāo)數(shù)據(jù)

小金山國控監(jiān)測斷面(見圖1)距上游浙江省與安徽省交界的街口斷面約28.6 km，距下游三潭島斷面約15.0 km，水深約40 m. 布設(shè)在小金山國控監(jiān)測斷面的水質(zhì)自動監(jiān)測系統(tǒng)為固定浮標(biāo)站(29°36′38″N、118°56′45″E)，浮標(biāo)懸掛有型號為YSI EXO2的多參數(shù)水質(zhì)儀(Yellow Springs Instrument Co.，美國)，可自上而下地同步記錄水深(單位為m)和Chla濃度(單位為μg/L). 水深為0.1～10 m時，測量間隔為0.5 m；水深為10～40 m時，測量間隔為2 m. 浮標(biāo)在2016年每3 h (分 別 為00:00、03:00、06:00、09:00、12:00、15:00、18:00、21:00)記錄一次數(shù)據(jù)，2017年以后每4 h(分別為00:00、04:00、08:00、12:00、16:00、20:00)記錄一次數(shù)據(jù). 為保持?jǐn)?shù)據(jù)的一致性，該研究取用每日00:00與12:00的數(shù)據(jù)進行分析.

圖 1 千島湖庫體及國控監(jiān)測斷面分布Fig.1 Qiandaohu Reservoir and location of state-controlled sections

1.2 數(shù)據(jù)處理

小金山國控監(jiān)測斷面逐日高頻剖面觀測數(shù)據(jù)的監(jiān)測時間為2016年9月8日-2019年12月2日.從Chla濃度隨時間的剖面變化〔見圖2(a)〕可見：Chla濃度剖面存在明顯的季節(jié)性變化特征，冬季、春季水柱垂向差別不明顯，呈幾乎完全混合的狀態(tài)；夏季、秋季上層Chla濃度明顯高于下層. 隨機取夏季某一次監(jiān)測剖面數(shù)據(jù)，分析Chla濃度隨水深的變化〔見圖2(b)〕發(fā)現(xiàn)：水深10 m以上，Chla濃度較高且變化較大；水深10 m以下，Chla濃度較低且趨于穩(wěn)定.因此，對每日00:00和12:00 0～10 m的Chla濃度監(jiān)測數(shù)據(jù)運用梯形求和公式，得到隨時間變化的Chla濃度在水深10 m以上的沿深平均值，時間間隔為0.5 d.

圖 2 小金山國控斷面高頻監(jiān)測Chla濃度的剖面變化情況Fig.2 High frequency observed Chla profile at Xiaojinshan Staion

高頻監(jiān)測期間，由于儀器維護或惡劣天氣等人為不可抗拒原因，出現(xiàn)了一些時長為半天至十幾天的測量間斷. 缺測值采用窗口長度為7的移動中位數(shù)進行線性填充，即以缺測值位置為中心，向前、向后各移動3個位置，形成長度為7的移動窗口；若沒有足夠的數(shù)據(jù)填滿窗口，則窗口自動在前或向后繼續(xù)移動或直至端點處截斷. 當(dāng)窗口被截斷時，根據(jù)窗口內(nèi)的數(shù)據(jù)計算出中位數(shù)，替代缺測值.

1.3 NARX神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)及構(gòu)建

NARX神經(jīng)網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu)分為輸入層、隱含層、輸出層3個層次以及輸出到輸入的延時(見圖3).神經(jīng)網(wǎng)絡(luò)的輸出延時保存后，通過外部反饋引入輸入層，與輸入樣本共同學(xué)習(xí)，按一定的訓(xùn)練標(biāo)準(zhǔn)計算網(wǎng)絡(luò)的實際輸出值與期望輸出值的誤差，不斷進行誤差反向傳播，從而調(diào)整網(wǎng)絡(luò)各層權(quán)重，使誤差達(dá)到最小，完成學(xué)習(xí)目的. 隱含層神經(jīng)元的數(shù)目根據(jù)經(jīng)驗公式〔見式(1)〕進行確定[32,35].

圖 3 NARX神經(jīng)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)Fig.3 Topology architecture of NARX neural network

式中，h為隱含層神經(jīng)元數(shù)目，n、m分別為輸入層、輸出層單元數(shù)，a為1～10之間的常數(shù).

NARX神經(jīng)網(wǎng)絡(luò)訓(xùn)練采用列文伯格-馬夸爾特(Levenberg-Marquardt，LM)算法[38-39]，計算公式：

式中，ω為權(quán)重閾值參數(shù)，J為雅克比矩陣，JT為J的轉(zhuǎn)置，μ為學(xué)習(xí)常數(shù)，I為單位矩陣，e為誤差向量.

NARX神經(jīng)網(wǎng)絡(luò)建模的基本步驟主要包括：①初始化網(wǎng)絡(luò)，包括網(wǎng)絡(luò)參數(shù)的選擇和設(shè)定；②訓(xùn)練，訓(xùn)練時應(yīng)盡量防止網(wǎng)絡(luò)過度擬合；③仿真，其中如何初始化網(wǎng)絡(luò)和提高模型泛化能力對神經(jīng)網(wǎng)絡(luò)的構(gòu)建十分關(guān)鍵.

1.4 NARX神經(jīng)網(wǎng)絡(luò)藻類預(yù)測模型性能評價

為了評價NARX神經(jīng)網(wǎng)絡(luò)藻類預(yù)測模型的預(yù)測性能，采用均方誤差(MSE)和相關(guān)系數(shù)(R)對預(yù)測值和實際觀測值進行對比分析. 均方誤差反映訓(xùn)練輸出值與目標(biāo)值之間的誤差，其值越小表示模擬效果越好；相關(guān)系數(shù)范圍在0～1之間，其值越大說明模擬數(shù)據(jù)和實測數(shù)據(jù)相關(guān)度越高，模型的模擬精度也越高[32,34].

2 結(jié)果與討論

為了探討NARX神經(jīng)網(wǎng)絡(luò)藻類預(yù)測模型在千島湖Chla濃度時序變化的預(yù)測情況，分別采用連續(xù)3 d的Chla濃度觀測值以及連續(xù)7 d的Chla濃度觀測值作為輸入樣本，未來0.5～7 d的Chla濃度作為輸出，建立輸入-輸出之間的響應(yīng)模型共28個.

模型采用3層神經(jīng)網(wǎng)絡(luò)，即輸入層、隱含層和輸出層各1層，其中將利用連續(xù)3 d Chla濃度預(yù)測未來0.5～7 d Chla濃度的預(yù)測模型定義為第1類模型，輸入層為6個單元，輸入向量Xi=(Ci,Ci+1, …,Ci+5)，輸出層為1個單元，輸出值Yi=Ci+j(j=6,7,…,19)；將利用連續(xù)7 d Chla濃度預(yù)測未來0.5～7 d Chla濃度的預(yù)測模型定義為第2類模型，輸入層為14個單元，輸入向量Xi=(Ci,Ci+1, …,Ci+13)，輸出層為1個單元，輸出值Yi=Ci+j(j=14,15,…,27). 該研究指定隱含層神經(jīng)元個數(shù)為10，延時階數(shù)為3 d. 初始化權(quán)值及閾值取[-1, 1]之間的隨機數(shù).

經(jīng)插補后，小金山監(jiān)測斷面2016年9月8日-2019年12月2日的Chla濃度數(shù)據(jù)共2 361個，按占總樣本數(shù)量75%、15%和10%的比例，分為訓(xùn)練集、驗證集和測試集. 其中，訓(xùn)練集樣本在模型訓(xùn)練過程中輸入到網(wǎng)絡(luò)中，神經(jīng)網(wǎng)絡(luò)在完成初始化之后，根據(jù)輸出值與標(biāo)注值之間的誤差不斷進行權(quán)值和偏置值的調(diào)整；驗證集樣本不直接參與上述訓(xùn)練調(diào)整，主要用于測度在訓(xùn)練過程中網(wǎng)絡(luò)泛化能力的表現(xiàn)，在泛化能力停止改進時停止訓(xùn)練，從而防止神經(jīng)網(wǎng)絡(luò)訓(xùn)練中發(fā)生過擬合現(xiàn)象，導(dǎo)致泛化能力的下降；測試集樣本對訓(xùn)練過程不施加影響，而在訓(xùn)練期間及訓(xùn)練后作為獨立于訓(xùn)練的樣本數(shù)據(jù)，對神經(jīng)網(wǎng)絡(luò)的性能進行測試、分析及評價.

經(jīng)過多次試驗對比發(fā)現(xiàn)，兩類模型在訓(xùn)練過程中訓(xùn)練集、測試集和驗證集的試驗效果均十分理想，相關(guān)系數(shù)在0.953～0.986之間，均方誤差在1.387～4.483之間，模型訓(xùn)練性能較好，可用于預(yù)測.

以2019年1-7月數(shù)據(jù)作為預(yù)測數(shù)據(jù)輸入已訓(xùn)練好的28個NARX神經(jīng)網(wǎng)絡(luò)進行計算，輸出未來0.5～7 d的Chla濃度預(yù)測值，并與實際觀測值進行比較，從而評價模型預(yù)測效果(見表1). 兩類模型的預(yù)測性能總體比較穩(wěn)定，相關(guān)系數(shù)保持在0.8～0.9之間，均方誤差在15～30之間. 隨著預(yù)見期的增加，兩類模型的均方誤差均呈上升趨勢，相關(guān)系數(shù)呈下降趨勢，但趨勢平緩，符合神經(jīng)網(wǎng)絡(luò)的預(yù)測性能，即隨著預(yù)測時間的增加，預(yù)測精度逐漸下降[23,32]. 相比而言，在未來0.5～4 d的預(yù)測中，第1類模型的預(yù)測精度優(yōu)于第2類模型，即使用連續(xù)3 d的Chla濃度作為輸入的預(yù)測模型精度較高；而在未來4.5～7 d的預(yù)測中，第2類模型的預(yù)測精度優(yōu)于第1類模型.

表 1 兩類模型預(yù)測性能對比Table 1 Prediction performance of two NARX models

圖4～6分別為第1類模型和第2類模型對于未來0.5、3和6 d Chla濃度預(yù)測值與觀測值的對比. 由圖4可見，對于未來0.5 d的預(yù)測，兩類模型預(yù)測值與觀測值的趨勢變化情況基本相似，吻合程度較高，但在第2類模型中存在一些時刻預(yù)測值略小于實測值，整體預(yù)測程度不如第1類模型. 由圖5可見，對于未來3 d的預(yù)測，兩類模型均存在某些預(yù)測值小于實測值的情況，但第1類模型與實測值吻合程度良好，而第2類模型明顯存在預(yù)測趨勢平緩的現(xiàn)象，尤其是沒有預(yù)測出6月Chla濃度的高值區(qū)，7月的Chla濃度存在明顯的相位差. 由圖6可見，對于未來6 d的預(yù)測，第2類模型的預(yù)測值與實測值吻合程度高于第1類模型，尤其是對5月后高低值的振蕩趨勢模擬較好，而第1類模型沒有預(yù)測出6月以后Chla濃度的高值區(qū).

在已構(gòu)建的預(yù)測千島湖水體富營養(yǎng)化的BP神經(jīng)網(wǎng)絡(luò)模型[22]中，為了預(yù)測Chla濃度的周尺度變化，采用了水溫、pH、Chla濃度、透明度和總氮5個輸入變量，并且隨著輸入變量減少為1個，模型訓(xùn)練過程中的相關(guān)系數(shù)降至0.744. 而NARX神經(jīng)網(wǎng)絡(luò)藻類預(yù)測模型僅用Chla濃度作為輸入變量，得到了比BP神經(jīng)網(wǎng)絡(luò)高的預(yù)測精度，進一步說明NARX神經(jīng)網(wǎng)絡(luò)藻類預(yù)測模型在時間序列預(yù)測中的優(yōu)勢，以及減少輸入?yún)?shù)可以降低輸入?yún)?shù)帶來的模型不確定性[23,40].

圖 4 兩類模型對未來0.5 d Chla濃度的預(yù)測值與觀測值對比Fig.4 Performance of two NARX models forecasting the dynamics of Chla in the future 0.5 day

圖 5 兩類模型對未來3 d Chla濃度的預(yù)測值與觀測值對比Fig.5 Performance of two NARX models forecasting the dynamics of Chla in the future 3 day

圖 6 兩類模型對未來6 d Chla濃度的預(yù)測值與觀測值對比Fig.6 Performance of two NARX models forecasting the dynamics of Chla in the future 6 day

3 結(jié)論

a)基于輸入樣本不同而定義的兩類NARX神經(jīng)網(wǎng)絡(luò)藻類預(yù)測模型預(yù)測性能總體比較穩(wěn)定，預(yù)測值與實測值相關(guān)系數(shù)保持在0.8～0.9之間，均方誤差在15～30之間，可以較為準(zhǔn)確地預(yù)測未來0.5～7 d的Chla濃度值，說明NARX神經(jīng)網(wǎng)絡(luò)藻類預(yù)測模型對于千島湖Chla時序變化的預(yù)測是可行的.

b)隨著預(yù)見期的變化，模型性能不盡相同. 建議在實際應(yīng)用過程中，采取二者相結(jié)合的方法，在未來0.5～4 d的預(yù)測中，使用連續(xù)3 d的Chla濃度作為輸入的預(yù)測模型；在未來4.5～7 d的預(yù)測中，使用連續(xù)7 d的Chla濃度作為輸入的預(yù)測模型.

c)相比于基于機理過程的藻類預(yù)測數(shù)值模型，基于數(shù)據(jù)驅(qū)動的NARX神經(jīng)網(wǎng)絡(luò)藻類預(yù)測模型結(jié)構(gòu)簡單，運行成本較低，具有較強的實用性和時效性. 如果采用其他點位的監(jiān)測數(shù)據(jù)對所構(gòu)建模型進行訓(xùn)練、測試和驗證，同樣可以實現(xiàn)對其他點位Chla時序變化的預(yù)測，因此又具有較好的可移植性.

d)目前預(yù)測獲得的時間序列平滑性有待進一步提升，這主要是因為原始時間序列為非平穩(wěn)序列，不僅包含較平穩(wěn)的有用信息，還包含影響數(shù)據(jù)可靠性和預(yù)測精度的隨機噪聲. 在接下來的研究中，將針對時間序列的特性，繼續(xù)構(gòu)建基于小波分析的NARX神經(jīng)網(wǎng)絡(luò)藻類預(yù)測模型，進一步優(yōu)化模型參數(shù)，提高Chla濃度的預(yù)測精度.

e)總體來講，基于NARX神經(jīng)網(wǎng)絡(luò)的藻類預(yù)測模型是可靠的，可以用于千島湖重點水域藻類生長的短期預(yù)測，為水華的監(jiān)測預(yù)警系統(tǒng)構(gòu)建提供借鑒與依據(jù).