基于NARX神經網絡的千島湖藻類短期預測模型構建

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

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

2. 杭州市生態環境局淳安分局，浙江 杭州 311700

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

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

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

NARX (nonlinear auto-regressive with exogenous inputs，帶有外部輸入的非線性自回歸)神經網絡是動態神經網絡的一種，其輸入結果是當前外部輸入和歷史輸出結果的非線性函數，可看作是有時延輸入的BP神經網絡加上輸出到輸入的延時反饋連接. 由于網絡結構中存在時延和反饋，因此其可反映系統的歷史狀態信息，是一種有記憶功能的神經網絡，可以模擬時間序列長期動態特征[26-30]. 近年來，NARX神經網絡在環境空氣質量[31-32]、水文預報[33]、大壩變形[34]、地表沉降[35]、農作物水分蒸騰[36]、股票擇時[37]等預測中均取得了較好的效果. 研究[27,34]表明，NARX神經網絡可用于預測時間序列，且通?？梢员Ａ粜畔⒌臅r間是常規遞歸神經網絡的2～3倍.

鑒于此，該研究以千島湖小金山國控監測斷面的水質自動監測站葉綠素a (Chla)高頻監測數據為研究對象，構建基于NARX神經網絡的藻類預測模型，分析該模型在千島湖Chla濃度時序變化預測的效果，探討最優預見期，以期為構建以數據驅動的千島湖水華監測預警系統提供科學依據.

1 材料與方法

1.1 研究點概況與浮標數據

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

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

1.2 數據處理

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

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

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

1.3 NARX神經網絡結構及構建

NARX神經網絡的拓撲結構分為輸入層、隱含層、輸出層3個層次以及輸出到輸入的延時(見圖3).神經網絡的輸出延時保存后，通過外部反饋引入輸入層，與輸入樣本共同學習，按一定的訓練標準計算網絡的實際輸出值與期望輸出值的誤差，不斷進行誤差反向傳播，從而調整網絡各層權重，使誤差達到最小，完成學習目的. 隱含層神經元的數目根據經驗公式〔見式(1)〕進行確定[32,35].

圖 3 NARX神經網絡拓撲結構Fig.3 Topology architecture of NARX neural network

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

NARX神經網絡訓練采用列文伯格-馬夸爾特(Levenberg-Marquardt，LM)算法[38-39]，計算公式：

式中，ω為權重閾值參數，J為雅克比矩陣，JT為J的轉置，μ為學習常數，I為單位矩陣，e為誤差向量.

NARX神經網絡建模的基本步驟主要包括：①初始化網絡，包括網絡參數的選擇和設定；②訓練，訓練時應盡量防止網絡過度擬合；③仿真，其中如何初始化網絡和提高模型泛化能力對神經網絡的構建十分關鍵.

1.4 NARX神經網絡藻類預測模型性能評價

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

2 結果與討論

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

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

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

經過多次試驗對比發現，兩類模型在訓練過程中訓練集、測試集和驗證集的試驗效果均十分理想，相關系數在0.953～0.986之間，均方誤差在1.387～4.483之間，模型訓練性能較好，可用于預測.

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

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

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

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

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

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

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

3 結論

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

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

c)相比于基于機理過程的藻類預測數值模型，基于數據驅動的NARX神經網絡藻類預測模型結構簡單，運行成本較低，具有較強的實用性和時效性. 如果采用其他點位的監測數據對所構建模型進行訓練、測試和驗證，同樣可以實現對其他點位Chla時序變化的預測，因此又具有較好的可移植性.

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

e)總體來講，基于NARX神經網絡的藻類預測模型是可靠的，可以用于千島湖重點水域藻類生長的短期預測，為水華的監測預警系統構建提供借鑒與依據.