基于NARX神經(jīng)網(wǎng)絡(luò)的千島湖藻類短期預(yù)測(cè)模型構(gòu)建

李 未，朱夢(mèng)圓，王裕成，朱廣偉

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

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

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

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

大數(shù)據(jù)及人工智能技術(shù)使得以數(shù)據(jù)驅(qū)動(dòng)的建模方式對(duì)時(shí)間序列進(jìn)行短期預(yù)測(cè)成為可能[14]，尤其是人工神經(jīng)網(wǎng)絡(luò)(artificial neural network，ANN)模型，因其具有良好的非線性映射和自學(xué)習(xí)能力，在非平穩(wěn)時(shí)間序列的高精度預(yù)測(cè)方面具有獨(dú)特的優(yōu)勢(shì). BP (back propagation，反向傳播)神經(jīng)網(wǎng)絡(luò)作為主流網(wǎng)絡(luò)之一，結(jié)構(gòu)簡(jiǎn)單、映射能力強(qiáng)，已在醫(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ò)，信息流動(dòng)是由輸入層單向流動(dòng)到輸出層，對(duì)于歷史輸入輸出缺乏記憶聯(lián)想，對(duì)時(shí)間序列時(shí)變特性的適應(yīng)能力不強(qiáng)，從而影響預(yù)測(cè)精度[24-25].

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

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

1 材料與方法

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

小金山國(guó)控監(jiān)測(cè)斷面(見圖1)距上游浙江省與安徽省交界的街口斷面約28.6 km，距下游三潭島斷面約15.0 km，水深約40 m. 布設(shè)在小金山國(guó)控監(jiān)測(cè)斷面的水質(zhì)自動(dòng)監(jiān)測(cè)系統(tǒng)為固定浮標(biāo)站(29°36′38″N、118°56′45″E)，浮標(biāo)懸掛有型號(hào)為YSI EXO2的多參數(shù)水質(zhì)儀(Yellow Springs Instrument Co.，美國(guó))，可自上而下地同步記錄水深(單位為m)和Chla濃度(單位為μg/L). 水深為0.1～10 m時(shí)，測(cè)量間隔為0.5 m；水深為10～40 m時(shí)，測(cè)量間隔為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ù)進(jìn)行分析.

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

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

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

圖 2 小金山國(guó)控?cái)嗝娓哳l監(jiān)測(cè)Chla濃度的剖面變化情況Fig.2 High frequency observed Chla profile at Xiaojinshan Staion

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

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

式中，ω為權(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)練時(shí)應(yīng)盡量防止網(wǎng)絡(luò)過度擬合；③仿真，其中如何初始化網(wǎng)絡(luò)和提高模型泛化能力對(duì)神經(jīng)網(wǎng)絡(luò)的構(gòu)建十分關(guān)鍵.

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

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

2 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

圖 6 兩類模型對(duì)未來6 d Chla濃度的預(yù)測(cè)值與觀測(cè)值對(duì)比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ù)測(cè)模型預(yù)測(cè)性能總體比較穩(wěn)定，預(yù)測(cè)值與實(shí)測(cè)值相關(guān)系數(shù)保持在0.8～0.9之間，均方誤差在15～30之間，可以較為準(zhǔn)確地預(yù)測(cè)未來0.5～7 d的Chla濃度值，說明NARX神經(jīng)網(wǎng)絡(luò)藻類預(yù)測(cè)模型對(duì)于千島湖Chla時(shí)序變化的預(yù)測(cè)是可行的.

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

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

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

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