基于BP神經(jīng)網(wǎng)絡(luò)的瀟河流域水質(zhì)預(yù)測

張秀菊， 王柳林， 李秀平， 王靈生

(1.河海大學(xué) 水文水資源學(xué)院， 江蘇 南京210098； 2.晉中市水利局， 山西 晉中 030600；)

1 研究背景

優(yōu)質(zhì)的水環(huán)境質(zhì)量是提升城鄉(xiāng)居民生活環(huán)境及人文內(nèi)涵的必要條件，改善水環(huán)境是各流域機(jī)構(gòu)與地方政府長期且艱巨的任務(wù)，而水質(zhì)預(yù)測則是水環(huán)境管理與保護(hù)中的重要環(huán)節(jié)。水質(zhì)預(yù)測是利用現(xiàn)有實測資料，通過建立水質(zhì)模型預(yù)測未來水環(huán)境狀況的一種方法?？茖W(xué)準(zhǔn)確的水質(zhì)預(yù)測有助于認(rèn)識水體環(huán)境的變化規(guī)律和發(fā)展態(tài)勢，可以為水環(huán)境保護(hù)與水污染防治工作提供技術(shù)支撐，提高管理部門的決策主動性。

國內(nèi)外在水質(zhì)預(yù)測方面的研究始于20世紀(jì)20年代。初期學(xué)者們主要關(guān)注機(jī)理性水質(zhì)模型，此類模型從物理、化學(xué)、生物等多方面因素研究水質(zhì)的變化機(jī)理，受到諸多學(xué)者的青睞[1]。機(jī)理性水質(zhì)模型先后經(jīng)歷了一維、二維及三維幾個發(fā)展階段，一維水質(zhì)模型有S-P水質(zhì)模型[2]、MIKE11模型[3]、 QUAL2K模型[4]等，隨后MIKE21水動力-水質(zhì)模型[5]、CE-QUAL-W2水質(zhì)模型[6]等二維水質(zhì)模型得到發(fā)展，二維模型考慮了污染物濃度在河流縱向和橫向上的演變，預(yù)測內(nèi)容較一維模型明顯增加。后來學(xué)者們研發(fā)了EFDC模型[7]、WASP模型[8]等三維水質(zhì)模型，可以模擬更加復(fù)雜的水質(zhì)狀況。

上述機(jī)理性水質(zhì)模型較為復(fù)雜，在實際運(yùn)用中需要大量的基礎(chǔ)資料作為支撐[9]，例如利用水動力-水質(zhì)耦合模型進(jìn)行水質(zhì)預(yù)測時需要有實測的河道斷面尺寸、河道流量、水位及水質(zhì)資料，且需要掌握研究河段范圍內(nèi)的時段降雨與蒸發(fā)量，以及精確的污染源、取排口等資料，基礎(chǔ)資料成為模擬與預(yù)測精度的主要約束。故而一些學(xué)者開始轉(zhuǎn)向研究非機(jī)理性水質(zhì)模型，如：陳昭明等[10]運(yùn)用多元回歸模型對漢豐湖水質(zhì)進(jìn)行了預(yù)測；黃玥等[11]采用時間序列模型對三峽水庫出庫水質(zhì)進(jìn)行了預(yù)測；Zhou等[12]利用灰色預(yù)測模型對水庫水質(zhì)進(jìn)行了預(yù)測；Liu等[13]利用支持向量機(jī)對水質(zhì)進(jìn)行了預(yù)測。上述非機(jī)理性水質(zhì)模型在一些流域得到了實際運(yùn)用，然而這些模型自身也存在諸多缺陷。例如回歸分析模型相對簡單，但其對統(tǒng)計數(shù)據(jù)的要求較高，需要大樣本且數(shù)據(jù)須有較好的分布規(guī)律；時間序列模型具有相對完善的理論基礎(chǔ)，但其預(yù)測的準(zhǔn)確性較差[14]；灰色預(yù)測模型適用于歷史數(shù)據(jù)較少且不連續(xù)的情況，但模型容易受到不穩(wěn)定數(shù)據(jù)的影響，導(dǎo)致預(yù)測誤差較大；支持向量機(jī)適用于小樣本，但其對參數(shù)和核函數(shù)的選擇較為敏感。

基于以上研究方法的不足，考慮到人工神經(jīng)網(wǎng)絡(luò)模型具有很強(qiáng)的自適性、自學(xué)習(xí)和高度容錯能力，本文以瀟河流域晉中段為研究區(qū)域，在Python環(huán)境下建立BP神經(jīng)網(wǎng)絡(luò)模型，通過模擬流域內(nèi)重鉻酸鹽指數(shù)(CODCr)和高錳酸鹽指數(shù)(CODMn)的時空變化規(guī)律，預(yù)測未來一定時段的水質(zhì)態(tài)勢，為水資源管理與水環(huán)境保護(hù)提供參考依據(jù)。

2 數(shù)據(jù)來源與研究方法

2.1 研究區(qū)概況

瀟河發(fā)源于山西省晉中市昔陽縣，流經(jīng)昔陽、壽陽、和順、榆次、小店、清徐6個縣區(qū)，最終匯入汾河，河長為147 km，總流域面積為3 894 km2，其中晉中段流域面積為3 781 km2。瀟河流域建有2座中型水庫，5座小型水庫，目前可以供水的水庫為松塔水庫、蔡莊水庫與鄭家莊水庫。流域內(nèi)有6條一級支流，現(xiàn)設(shè)有松曲、馬坊、蔡莊水庫、冷泉寺、東趙、郝村6個水質(zhì)監(jiān)測斷面。瀟河流域晉中段水系分布及各水質(zhì)監(jiān)測斷面位置如圖1所示。

圖1 瀟河流域晉中段水系分布及各水質(zhì)監(jiān)測斷面位置

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

水質(zhì)數(shù)據(jù)取自瀟河流域6個監(jiān)測斷面的水質(zhì)資料，數(shù)據(jù)來源為晉中市水利局。本次選取的數(shù)據(jù)長度為2017年1月-2020年8月，由于部分月份河流有斷流、結(jié)冰或缺測情況，故導(dǎo)致這些月份水質(zhì)數(shù)據(jù)缺失，水質(zhì)數(shù)據(jù)缺失情況統(tǒng)計見表1。為了確保水質(zhì)數(shù)據(jù)的連續(xù)性，本文運(yùn)用拉格朗日插值法補(bǔ)充缺失值[15]，拉格朗日插值多項式如下：

表1 研究區(qū)水質(zhì)指標(biāo)數(shù)據(jù)缺失情況統(tǒng)計表

(1)

式中：th、t、th+1為相鄰的3個月份；xh、L1(t)、xh+1為相鄰3個月份對應(yīng)的水質(zhì)指標(biāo)(CODCr或CODMn)濃度，mg/L，其中L1(t)為數(shù)據(jù)缺失月份水質(zhì)指標(biāo)濃度插值。

2.3 研究方法

2.3.1 BP神經(jīng)網(wǎng)絡(luò)原理 1957年，Rosenblatt基于大腦神經(jīng)元工作原理發(fā)明了感知器，它是最簡單的人工神經(jīng)網(wǎng)絡(luò)，僅能解決線性問題。1974年，Paul Werbos發(fā)明了反向傳播算法，使得利用多層神經(jīng)網(wǎng)絡(luò)求解非線性問題成為可能，由此產(chǎn)生了BP(back propagation)神經(jīng)網(wǎng)絡(luò)。圖2為BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖。

圖2 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖

圖3 BP神經(jīng)網(wǎng)絡(luò)模型流程圖

2.3.2 數(shù)據(jù)處理 由于部分監(jiān)測斷面的實測水質(zhì)數(shù)據(jù)波動范圍較大，可能影響模型的收斂速度，因此在訓(xùn)練前需要對數(shù)據(jù)進(jìn)行預(yù)處理。本研究選擇改進(jìn)的最大最小值歸一化方法，可以有效避免出現(xiàn)0值，從而降低原始數(shù)據(jù)對模型訓(xùn)練的影響[17]。公式(2)與(3)是對數(shù)據(jù)進(jìn)行預(yù)處理和反預(yù)處理的方法。

數(shù)據(jù)預(yù)處理方法的計算公式為：

(2)

反預(yù)處理方法的計算公式為：

(3)

式中：x為CODCr或CODMn的原始數(shù)值；x′為經(jīng)數(shù)據(jù)處理后的數(shù)據(jù)； min(x)和max(x)分別為CODCr或CODMn數(shù)值的最小值和最大值。

2.3.3 模型構(gòu)建 BP神經(jīng)網(wǎng)絡(luò)模型構(gòu)建主要包括神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)設(shè)置、激活函數(shù)的選取、求解器的選取3個部分。

(1)BP神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)設(shè)置。理論上，對于一個3層的BP神經(jīng)網(wǎng)絡(luò)，只要隱含層節(jié)點數(shù)足夠多，則可以逼近任意非線性函數(shù)[18]。本研究確定BP神經(jīng)網(wǎng)絡(luò)的總層數(shù)為3層。

BP神經(jīng)網(wǎng)絡(luò)在訓(xùn)練過程中，輸入層節(jié)點數(shù)設(shè)定為k個，輸出層節(jié)點數(shù)為1個。輸出層，即第n個節(jié)點的CODCr或CODMn的濃度為前k個節(jié)點的歷史數(shù)據(jù)的函數(shù)，其預(yù)測動態(tài)系統(tǒng)[19]為：

yn=f(xn-1,xn-2,…,xn-k)

(4)

式中：yn為對應(yīng)編號為n的CODCr或CODMn；k為滯后編號數(shù)。

這樣整個系統(tǒng)的輸入為某編號為n的前k個編號的CODCr或CODMn，而輸出則是第n個節(jié)點的CODCr或CODMn。若節(jié)點數(shù)太少，則網(wǎng)絡(luò)訓(xùn)練性能會變差，而節(jié)點數(shù)太多，容易出現(xiàn)過擬合的狀況。因此，本研究在程序設(shè)計時，設(shè)定了相應(yīng)區(qū)間和步長，利用計算機(jī)自動尋優(yōu)確定節(jié)點數(shù)量，最終確定的節(jié)點數(shù)見表2。以表2中的松曲斷面CODCr為例，其對應(yīng)的輸入層節(jié)點數(shù)為6個，在模型訓(xùn)練期間，首先輸入2017年1月-2017年6月的CODCr濃度實測值，通過模型訓(xùn)練輸出2017年7月的CODCr濃度模擬值，然后依次滑動，最后輸入2019年11月-2020年4月的CODCr濃度實測值，通過模型訓(xùn)練輸出2020年5月的CODCr濃度模擬值。

(2)激活函數(shù)的選取。常見的激活函數(shù)有identity函數(shù)、logistic函數(shù)、tanh函數(shù)和relu函數(shù)。其中identity函數(shù)為線性函數(shù)，當(dāng)計算問題為非線性時，該激活函數(shù)不適于單獨使用；logistic函數(shù)和tanh函數(shù)計算量大且容易出現(xiàn)梯度消失問題；relu函數(shù)可以加快收斂速度，且解決了梯度彌散問題，故本研究選用relu函數(shù)，即f(x)=max(0,x)。

(3)求解器的選取。求解器可以用來優(yōu)化權(quán)重，常見的求解器有l(wèi)bfgs、sgd、adam。其中sgd求解器每步接受的信息有限，造成目標(biāo)函數(shù)收斂不穩(wěn)定甚至不收斂；adam求解器主要適用于數(shù)據(jù)集較多的情況；lbfgs為準(zhǔn)牛頓方法家族中的優(yōu)化程序，對于小型數(shù)據(jù)集，lbfgs可以收斂得更快并且性能更好，故本研究選用lbfgs求解器。

3 結(jié)果與分析

3.1 模型訓(xùn)練

以瀟河流域2017年1月-2020年5月監(jiān)測的水質(zhì)指標(biāo)CODCr和CODMn逐月濃度數(shù)據(jù)作為訓(xùn)練集，利用構(gòu)建好的模型對6個監(jiān)測斷面的水質(zhì)數(shù)據(jù)進(jìn)行擬合，模型的最大學(xué)習(xí)次數(shù)設(shè)為2 000，學(xué)習(xí)速率為0.001。模型訓(xùn)練結(jié)果見圖4和5，其中，圖4為各水質(zhì)監(jiān)測斷面CODCr濃度模擬結(jié)果，圖5為各水質(zhì)監(jiān)測斷面CODMn濃度模擬結(jié)果。由于各斷面BP神經(jīng)網(wǎng)絡(luò)輸入節(jié)點數(shù)不一致(見表2)，故圖4和5存在各斷面起始模擬值的時間不一致的情況。

圖4 研究區(qū)各水質(zhì)監(jiān)測斷面CODCr濃度模擬結(jié)果

圖5 研究區(qū)各水質(zhì)監(jiān)測斷面CODMn濃度模擬結(jié)果

表2 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)表

由圖4、5可以看出，瀟河流域各監(jiān)測斷面的實測水質(zhì)數(shù)據(jù)在時間尺度上呈現(xiàn)非線性變化，而BP神經(jīng)網(wǎng)絡(luò)模擬值與水質(zhì)數(shù)據(jù)實測值的變化趨勢基本一致，表明BP神經(jīng)網(wǎng)絡(luò)經(jīng)訓(xùn)練后能較好地反映水質(zhì)變化趨勢。

現(xiàn)將各斷面水質(zhì)指標(biāo)的實測值與模擬值進(jìn)行對比分析，用相關(guān)性系數(shù)和平均相對誤差來衡量模型擬合的效果，表3為水質(zhì)指標(biāo)擬合效果表。

由表3可以看出，6個監(jiān)測斷面的CODCr濃度相關(guān)性系數(shù)均超過了0.98，為強(qiáng)相關(guān)。相關(guān)性系數(shù)由大到小排序為：東趙>冷泉寺>郝村>馬坊>蔡莊水庫>松曲，相關(guān)性系數(shù)最大的為東趙斷面，其相關(guān)性系數(shù)為0.999 8，相關(guān)性系數(shù)最小的為松曲斷面，相關(guān)性系數(shù)為0.989 2。6個監(jiān)測斷面的CODCr濃度平均相對誤差的平均值約為3%，平均相對誤差由小到大排序為：東趙<郝村<冷泉寺<馬坊<蔡莊水庫<松曲，平均相對誤差最大的為松曲斷面，對應(yīng)的平均相對誤差為6.59%，平均相對誤差最小的為東趙斷面，對應(yīng)的平均相對誤差為0.63%。整體而言，各斷面CODCr濃度的擬合效果較好，其中東趙斷面的擬合效果最佳。從空間位置上看，下游監(jiān)測斷面的CODCr濃度模擬效果好于上游斷面。6個監(jiān)測斷面的CODMn濃度相關(guān)性系數(shù)均超過了0.99，屬于強(qiáng)相關(guān)。相關(guān)性系數(shù)由大到小排序為：冷泉寺>松曲>蔡莊水庫>東趙>郝村>馬坊，相關(guān)性系數(shù)最大的為冷泉寺斷面，對應(yīng)的相關(guān)性系數(shù)為0.999 9，相關(guān)性系數(shù)最小的為馬坊斷面，對應(yīng)的相關(guān)性系數(shù)為0.993 1。6個監(jiān)測斷面的CODMn濃度平均相對誤差的平均值為1.5%，平均相對誤差由小到大排序為：冷泉寺<松曲<蔡莊水庫<東趙<郝村<馬坊，平均相對誤差最大的為馬坊斷面，對應(yīng)的平均相對誤差為4.63%，平均相對誤差最小的為冷泉寺斷面，對應(yīng)的平均相對誤差為0.15%。整體而言，各斷面CODMn濃度的擬合效果較好，其中冷泉寺斷面的擬合效果最佳。從空間位置上看，上游監(jiān)測斷面CODMn濃度的模擬效果好于下游斷面。

表3 研究區(qū)各斷面水質(zhì)指標(biāo)擬合效果表

3.2 模型驗證

利用模型模擬期優(yōu)化得到的權(quán)值與閾值，對瀟河流域2020年6-8月的水質(zhì)進(jìn)行驗證。具體檢驗結(jié)果見表4和5。其中表4為各斷面CODCr濃度模型驗證結(jié)果，表5為各斷面CODMn濃度模型驗證結(jié)果。

表4 2020年6-8月研究區(qū)各水質(zhì)斷面CODCr濃度模型驗證結(jié)果 mg/L

表5 2020年6-8月研究區(qū)各水質(zhì)斷面CODMn濃度模型驗證結(jié)果 mg/L

根據(jù)表4中的驗證結(jié)果，蔡莊水庫斷面在2020年6月份、東趙斷面在2020年8月份、郝村斷面在2020年7月份的CODCr濃度相對誤差較大，各斷面其余月份的CODCr濃度的相對誤差均未超過20%。因此，構(gòu)建的BP神經(jīng)網(wǎng)絡(luò)模型通過檢驗，可以用于預(yù)測瀟河流域的CODCr濃度；根據(jù)表5中的驗證結(jié)果，除冷泉寺斷面在2020年7月份的誤差較大外，其余斷面各月份的CODMn濃度驗證結(jié)果較好。因此，構(gòu)建的BP神經(jīng)網(wǎng)絡(luò)模型通過檢驗，可以用于預(yù)測瀟河流域的CODMn濃度。

4 討 論

BP神經(jīng)網(wǎng)絡(luò)對于歷史數(shù)據(jù)具有很好的適應(yīng)性，能較好地適用于水質(zhì)濃度變化緩慢的情況，且模擬與驗證精度較高。但若由于人為因素影響，個別月份水質(zhì)指標(biāo)濃度發(fā)生較大變化時，BP模型的精度會降低?？傮w而言，本文在模擬期各斷面水質(zhì)指標(biāo)的平均相對誤差均小于7%，模擬精度較高。

根據(jù)相關(guān)文獻(xiàn)，驗證期的水質(zhì)預(yù)測誤差是不易控制的，例如董瑞瑞等[20]利用MIKE11水動力-水質(zhì)模型對漢江中下游CODMn濃度進(jìn)行預(yù)測時，其驗證期的平均相對誤差為14.5%；王瞾瑞等[21]利用一維水質(zhì)預(yù)測模型對CODCr濃度進(jìn)行預(yù)測，驗證期CODCr濃度的平均相對誤差為24%。本研究中，在驗證期，月最大相對誤差未超過30%(東趙斷面)，各斷面大多數(shù)月份的誤差小于15%。其中，各斷面CODCr濃度的平均相對誤差均小于18%，CODMn濃度的平均相對誤差均小于12%，驗證精度優(yōu)于其他模型。因此，本文驗證精度較優(yōu)，其預(yù)測結(jié)果可以接受。

為剖析蔡莊水庫、東趙、郝村、冷泉寺斷面部分月份水質(zhì)指標(biāo)驗證誤差偏大的原因，對該4個斷面水質(zhì)指標(biāo)的實際檢測濃度進(jìn)行比較，2017年1月-2020年8月該4個斷面的水質(zhì)指標(biāo)濃度變化情況見圖6。

由圖6(a)可以看出，蔡莊水庫斷面CODCr濃度在2017年1月-2020年4月平穩(wěn)振蕩，2020年4月之后，CODCr實測濃度明顯降低，較其他年份同時期CODCr濃度值小，由此導(dǎo)致2020年6月CODCr濃度模擬值大于實測值，相對誤差較大。根據(jù)瀟河實測流量資料，蔡莊水庫在2020年4月-2020年8月的河道流量偏小，故推測CODCr濃度降低的原因應(yīng)該是當(dāng)?shù)厮块T實施了嚴(yán)格的污染物入河管控措施，污染物濃度出現(xiàn)減小趨勢。

由圖6(b)可以看出，東趙斷面CODCr濃度在2017年1月-2020年6月基本保持平穩(wěn)振蕩，但2018年4月和2020年8月存在兩個峰值，其濃度遠(yuǎn)高于歷史平均數(shù)值，例如2020年8月CODCr的實測濃度是往年相同月份的2倍。經(jīng)核對，東趙斷面在2020年8月份的河道流量較往年大，故CODCr濃度大幅度增大的原因可能是存在臨時性或突發(fā)性的入河污染物增多的問題。

由圖6(c)可以看出，郝村斷面CODCr濃度在2017年1月-2020年8月平穩(wěn)振蕩，而郝村斷面在2020年7月的CODCr濃度相對誤差較大，其原因應(yīng)是由于模型自身存在誤差，BP神經(jīng)網(wǎng)絡(luò)出現(xiàn)了陷入局部極值的現(xiàn)象。

由圖6(d)可以看出，冷泉寺斷面在冬季1或2月份均出現(xiàn)檢測峰值，但總體上CODMn濃度呈現(xiàn)減小的態(tài)勢，且2020年4月之后，CODMn濃度與其他年份同時期相比顯著偏小，由此導(dǎo)致2020年7月驗證的相對誤差變大。由于冷泉寺斷面在流域各斷面中水質(zhì)相對較差，近些年流域機(jī)構(gòu)加強(qiáng)了對周邊入河污染物的管控，這也是該斷面CODMn濃度下降的原因。

圖6 2017年1月-2020年8月4個斷面的實測水質(zhì)指標(biāo)濃度變化情況

BP神經(jīng)網(wǎng)絡(luò)模型只考慮了瀟河流域各水質(zhì)指標(biāo)的歷史數(shù)據(jù)，而外部環(huán)境的改變對河流水質(zhì)的影響較大，會對神經(jīng)網(wǎng)絡(luò)訓(xùn)練過程產(chǎn)生干擾，從而影響模型精度。如何降低外部因素的干擾，亦或考慮將影響水質(zhì)的影響因素加入模型中，還有待進(jìn)一步研究。

5 結(jié)論與展望

本文分析了傳統(tǒng)水質(zhì)預(yù)測方法的優(yōu)缺點，通過構(gòu)建BP神經(jīng)網(wǎng)絡(luò)模型，利用瀟河流域晉中段實測水質(zhì)資料進(jìn)行水質(zhì)模擬與權(quán)值優(yōu)化，進(jìn)而對瀟河流域6個水質(zhì)監(jiān)測斷面的CODCr和CODMn指標(biāo)未來3個時段的濃度進(jìn)行預(yù)測，主要結(jié)論與展望如下：

(1)BP神經(jīng)網(wǎng)絡(luò)能較好地反映水質(zhì)數(shù)據(jù)的變化趨勢，訓(xùn)練集的各斷面水質(zhì)指標(biāo)平均相對誤差均小于7%，相關(guān)性系數(shù)均超過了0.98；驗證集的各斷面水質(zhì)指標(biāo)平均相對誤差均低于18%，預(yù)測精度高，可用于預(yù)測瀟河流域的CODCr和CODMn濃度，為水資源管理與水環(huán)境保護(hù)提供參考依據(jù)。

(2)本研究在選取神經(jīng)網(wǎng)絡(luò)節(jié)點數(shù)量時，在程序中設(shè)置了相應(yīng)區(qū)間和步長，對訓(xùn)練集的平均相對誤差和相關(guān)性系數(shù)進(jìn)行了約束，利用計算機(jī)自動尋優(yōu)，降低了模型的不確定性，節(jié)省了參數(shù)率定的時間。

(3)為了提高收斂速度與預(yù)測精度，以及改善過擬合的問題，將人工神經(jīng)網(wǎng)絡(luò)與其他算法的結(jié)合應(yīng)是今后研究的重點。