陶志勇，胡啟振，任曉奎

(遼寧工程技術(shù)大學(xué) 電子與信息工程學(xué)院，遼寧 葫蘆島 125105)

河流中溶解氧（Dissolved Oxygen, DO）能反映水質(zhì)污染程度和水體自身凈化能力[1]，是衡量水質(zhì)優(yōu)劣的重要指標(biāo)之一. 當(dāng)河流受到有機(jī)物污染，水體的復(fù)氧速度趕不上分解有機(jī)物的耗氧速度時，會導(dǎo)致溶解氧含量快速下降，造成水中需氧生物和各種魚類的死亡[2]. 研究發(fā)現(xiàn)，水體復(fù)氧不利是河流水質(zhì)惡化的主要原因[3]. 因此，準(zhǔn)確預(yù)測未來河流中溶解氧的含量，掌握它的變化趨勢，可以及時做好河流污染防治工作，對生態(tài)環(huán)境的保護(hù)具有重要意義.

目前，溶解氧預(yù)測方法主要有統(tǒng)計(jì)預(yù)測、機(jī)器學(xué)習(xí)預(yù)測以及組合模型預(yù)測等. 統(tǒng)計(jì)預(yù)測的常用模型有差分自回歸移動平均模型（Autoregressive Integrated Moving Average Model，ARIMA）[4]等，如黃玥等[5]應(yīng)用ARIMA 模型對三峽庫區(qū)出口和入口端的DO 等水質(zhì)數(shù)據(jù)進(jìn)行預(yù)測. 隨著人工智能技術(shù)的發(fā)展，機(jī)器學(xué)習(xí)中的人工神經(jīng)網(wǎng)絡(luò)和支持向量機(jī)(Support Vector Machines，SVM)被廣泛應(yīng)用到水質(zhì)預(yù)測中，李鑫鑫等[6]用遺傳算法（Genetic Algorithm，GA）優(yōu)化神經(jīng)網(wǎng)絡(luò)的權(quán)值和閾值，建立GA-BP 神經(jīng)網(wǎng)絡(luò)模型對水質(zhì)參數(shù)進(jìn)行預(yù)測. 由于河流中溶解氧的質(zhì)量濃度參數(shù)具有非線性和不穩(wěn)定性的特征，較為復(fù)雜，難以進(jìn)行數(shù)學(xué)建模[7-8]，所以單一模型預(yù)測精度較低. 為解決上述問題，劉晨等[9]提出EEMD-FA-SVM 的組合預(yù)測模型，首先用集合經(jīng)驗(yàn)?zāi)B(tài)分解（Ensemble Empirical Mode Decomposition，EEMD）削弱DO 序列的復(fù)雜性，得到相對穩(wěn)定的子序列，然后將子序列分別用螢火蟲算法優(yōu)化的SVM 做預(yù)測，將各分量預(yù)測結(jié)果疊加，得到最終的預(yù)測結(jié)果. 盧毅敏等[10]提出CEEMDANSE-CS-Elman 組合預(yù)測模型，首先用自適應(yīng)噪聲的完整集成經(jīng)驗(yàn)?zāi)B(tài)分解（Complete Ensemble Empirical Mode Decomposition with Adaptive Noise，CEEMDAN）對DO 序列降噪，然后計(jì)算各分量的樣本熵值，將熵值相近的序列進(jìn)行合并，最后將各個序列分別進(jìn)行預(yù)測并將結(jié)果疊加，得到最終結(jié)果. 組合模型雖提升了DO 預(yù)測精度，但是EEMD和CEEMDAN 分解后的高頻率本征模函數(shù)（Intrinsic Mode Functions，IMF）易受到隨機(jī)因素的影響[11-12]，序列復(fù)雜性相對較高，導(dǎo)致模型預(yù)測精度降低.

在上述研究的基礎(chǔ)上，本文提出TDT-SSA-BP組合預(yù)測模型，首先應(yīng)用二層分解技術(shù)，將CEEMDAN 分解后的高頻率IMF 用變分模態(tài)分解（Variational Modal Decomposition，VMD）再次處理，從而充分削弱DO 序列的非線性和不平穩(wěn)性；然后提出麻雀搜索算法（Sparrow Search Algorithm，SSA）優(yōu)化神經(jīng)網(wǎng)絡(luò)的溶解氧預(yù)測模型，其中SSA 是2020 年提出的新型群智能優(yōu)化算法，具有更好的全局搜索和局部尋優(yōu)的能力以及更快的收斂速度[13]；接下來將二層分解得到的所有分量輸入到SSABP 模型中，得到各分量的預(yù)測結(jié)果，把各分量預(yù)測結(jié)果重構(gòu)，即得到最終的預(yù)測結(jié)果. 通過實(shí)驗(yàn)分析，本文所提模型相比于單一模型或其它組合模型擁有更高的預(yù)測精度.

1 模型介紹

本文所提TDT-SSA-BP 組合預(yù)測模型架構(gòu)如圖1 所示. 具體步驟如下：

圖1 TDT-SSA-BP 模型流程圖Fig. 1 Flow chart of TDT-SSA-BP model

步驟 1 采集DO 質(zhì)量濃度時間序列數(shù)據(jù)，對序列進(jìn)行預(yù)處理，形成數(shù)據(jù)集.

步驟 2 采用CEEMDAN 算法對數(shù)據(jù)集進(jìn)行第一層分解，得到模態(tài)分量IMF1，IMF2，···，IMFn-1和余量Residual.

步驟 3 用PE（Permutation Entropy，PE）計(jì)算步驟2 中第一層分解后各分量的排列熵值.

步驟 4 將排列熵值較高的模態(tài)分量采用VMD 算法進(jìn)行第二層分解，得到新的模態(tài)分量Mode1，Mode2，···，Modem.

步驟 5 使用SSA 算法優(yōu)化BP 神經(jīng)網(wǎng)絡(luò)的權(quán)值和閾值，建立SSA-BP 預(yù)測模型.

步驟6 把步驟2 和步驟4 中分解得到的所有模態(tài)分量輸入到SSA-BP 預(yù)測模型中，得到子序列的預(yù)測結(jié)果.

步驟7 重構(gòu)子序列的預(yù)測結(jié)果，得到最終的DO 質(zhì)量濃度預(yù)測結(jié)果.

1.1 數(shù)據(jù)預(yù)處理 由于設(shè)備故障或者人為等原因，采集到的原始DO 質(zhì)量濃度（mg/L）時間序列存在離群值和缺失值的情況. 若對離群值和缺失值不做處理，會造成時間序列數(shù)據(jù)起伏過大或者間斷，導(dǎo)致預(yù)測精度降低.

鑒于此，依據(jù)GB3838—2002《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》和箱線圖剔除離群值，因?yàn)槿笔?shù)據(jù)較少且短時間內(nèi)DO 質(zhì)量濃度不會有大幅度的波動，所以本文選擇中位數(shù)填充法對離群值和缺失值進(jìn)行填充.

此外，為提升預(yù)測模型的收斂速度和精度，在輸入神經(jīng)網(wǎng)絡(luò)前需對數(shù)據(jù)集進(jìn)行歸一化處理，公式如下：

1.2 二層分解技術(shù)

步驟 3 重復(fù)執(zhí)行步驟1 和步驟2，當(dāng)余量信號的極值點(diǎn)不超過2 個，此時余量無法繼續(xù)分解，最終DO 時序信號可表示為：

1.2.2 變分模態(tài)分解 相比于經(jīng)驗(yàn)?zāi)B(tài)分解，VMD[15]具有更完備的數(shù)學(xué)理論支撐，通過對排列熵值較高的IMF 分解，可有效降低時間序列的非線性和不穩(wěn)定性. 設(shè)DO 時間序列信號f(t)分解為有限帶寬的模態(tài)分量uk(t)且 中心頻率為 ωk，其算法步驟如下：

步驟 1 用Hilbert 對模態(tài)分量uk(t)進(jìn)行處理，得到解析信號和單邊頻譜.

步驟 2 將解析信號中混合中心頻率 e-jωkt，從而可將每個模態(tài)的頻譜轉(zhuǎn)換到基帶上.

步驟 3 通過高斯平滑計(jì)算uk(t)的帶寬.

因此，產(chǎn)生變分約束問題可表示為：

1.3 SSA-BP 預(yù)測模型及原理

1.3.1 麻雀搜索算法 SSA 算法通過模擬麻雀覓食的生物特性，將麻雀分為發(fā)現(xiàn)者、加入者和偵察者. 其中發(fā)現(xiàn)者負(fù)責(zé)尋找食物源，占種群的10%～20%；加入者則通過跟隨發(fā)現(xiàn)者來尋找食物；偵察者受到威脅時，位于種群邊緣的麻雀會向安全區(qū)移動，位于群體中央的麻雀會隨機(jī)移動，從而避免被捕食者攻擊.

發(fā)現(xiàn)者位置更新如下所示：

1.3.2 BP 神經(jīng)網(wǎng)絡(luò) BP 神經(jīng)網(wǎng)絡(luò)是多層前饋型神經(jīng)網(wǎng)絡(luò)，通常由輸入層、隱含層和輸出層3 層網(wǎng)絡(luò)結(jié)構(gòu)組成，如圖2 所示，各層神經(jīng)元之間是全互連連接，同一層中的神經(jīng)元之間互不相連，且傳遞函數(shù)一般采用sigmoid 函數(shù).

圖2 單隱含層BP 神經(jīng)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)Fig. 2 Topological structure of single hidden layer BP neural network

BP 網(wǎng)絡(luò)主要由信號的正向傳輸和誤差的反向傳遞兩個階段組成. 第一個階段中，輸入信號由輸入層進(jìn)入神經(jīng)網(wǎng)絡(luò)，然后經(jīng)過隱含層神經(jīng)元的處理 ，最后到達(dá)輸出層. 第二個階段中，若輸出結(jié)果和實(shí)際結(jié)果之間的誤差不在設(shè)定的范圍內(nèi)，則會進(jìn)行誤差的反向傳遞，由輸出層經(jīng)隱含層傳入到輸入層，并不斷地調(diào)整神經(jīng)元之間的權(quán)值和閾值使誤差降低，此過程反復(fù)循環(huán)，直到滿足精度要求或者達(dá)到最大迭代次數(shù)時，訓(xùn)練結(jié)束. 由于BP 神經(jīng)網(wǎng)絡(luò)具有較好的自學(xué)習(xí)和非線性映射能力，目前被廣泛應(yīng)用在各個領(lǐng)域.

1.3.3 麻雀搜索算法優(yōu)化神經(jīng)網(wǎng)絡(luò) BP 神經(jīng)網(wǎng)絡(luò)通過梯度下降法改變權(quán)值和閾值，但在訓(xùn)練過程中易陷入局部最優(yōu)解. SSA 算法具有更好的全局搜索和局部尋優(yōu)的能力，用來調(diào)整神經(jīng)網(wǎng)絡(luò)的權(quán)值和閾值，可有效提高預(yù)測精度.

SSA-BP 預(yù)測模型建立步驟如下：

步驟 1 選取樣本數(shù)據(jù)，對數(shù)據(jù)作歸一化處理，按一定比例劃分為訓(xùn)練集和測試集；確定神經(jīng)網(wǎng)絡(luò)的隱含層層數(shù)以及各層神經(jīng)元個數(shù).

步驟 2 設(shè)置麻雀種群規(guī)模為20，最大迭代次數(shù)為100，發(fā)現(xiàn)者比例為0.2，警戒值為0.8.

步驟 3 生成初始種群位置，用預(yù)測值和真實(shí)值的誤差矩陣的范數(shù)作為適應(yīng)度函數(shù)的輸出，計(jì)算所有個體的適應(yīng)度值，并根據(jù)適應(yīng)度值的大小排序.

步驟 4 將適應(yīng)度值較大的前4 個麻雀作為發(fā)現(xiàn)者，根據(jù)位置更新公式（8）進(jìn)行全局搜索；剩余的麻雀作為加入者，會跟隨發(fā)現(xiàn)者覓食，根據(jù)位置更新公式（9）進(jìn)行局部搜索；當(dāng)偵察者意識到危險(xiǎn)時，也會根據(jù)位置更新公式（10）進(jìn)行局部搜索.

步驟 5 保存擁有最佳適應(yīng)度值的麻雀位置，若達(dá)到最大迭代次數(shù)時，則迭代結(jié)束，否則重復(fù)步驟3 ～ 4，不斷更新最佳適應(yīng)度值.

步驟 6 將步驟5 中得到的最優(yōu)解作為神經(jīng)網(wǎng)絡(luò)的權(quán)值和閾值，建立SSA-BP 模型.

2 實(shí)驗(yàn)及分析

2.1 數(shù)據(jù)來源 本文數(shù)據(jù)來源于2016 年10 月21日至2017 年11 月21 日江蘇省無錫市長江水質(zhì)實(shí)時監(jiān)測站，監(jiān)測頻率為每日一次，共計(jì)397 個數(shù)據(jù)，經(jīng)預(yù)處理后，數(shù)據(jù)集如圖3 所示. 從圖3 中可以看出，DO 時間序列具有明顯的隨機(jī)性和不平穩(wěn)性.選 取2016 年10 月21 日 至2017 年10 月15 日 共計(jì)360 個數(shù)據(jù)作為訓(xùn)練集，剩余2017 年10 月16日至11 月21 日共37 個數(shù)據(jù)作為測試集.

圖3 原始DO 質(zhì)量濃度時序數(shù)據(jù)Fig. 3 Original DO concentration time series data

2.2 評價指標(biāo) 本文采用4 種評價指標(biāo)直觀地展示各個模型的預(yù)測效果，分別是平均絕對誤差（Mean Absolute Error，MAE）、均方根誤差（Root Mean Square Error，RMSE）、平均絕對百分比誤差（Mean Absolute Percentage Error，MAPE）和決定系數(shù)R2，計(jì)算公式如下：

2.3 實(shí)驗(yàn)仿真 首先使用CEEMDAN 分解DO 時間序列，得到IMF1，IMF2，···，IMF7和Residual，結(jié)果如圖4 所示. 由圖4 可知，IMF1具有較為復(fù)雜的波動特征，反映DO 序列受到隨機(jī)影響；IMF2～I(xiàn)MF7序列相對平穩(wěn)且具有一定的周期性；Residual時頻特征平緩下降，反映了DO 時間序列的長期趨勢.

圖4 CEEMDAN 分解DO 時間序列結(jié)果Fig. 4 DO time series results of CEEMDAN decomposition

接下來計(jì)算各分量的排列熵值，如圖5 所示.由圖5 可知，各個分量的排列熵值整體呈下降趨勢，說明序列的復(fù)雜程度和隨機(jī)性在逐漸減小. 其中IMF1的熵值遠(yuǎn)高于其他分量的熵值，可知IMF1時間序列的波動特征較為復(fù)雜，受到隨機(jī)因素的影響較大.

圖5 各分量排列熵值Fig. 5 Permutation entropy of each subcomponent

根據(jù)各序列的排列熵值以及考慮到預(yù)測效率的原因，本文僅對熵值較高的分量IMF1使用VMD進(jìn)行二次分解. 由于VMD 的模態(tài)分解個數(shù)K需要人為選擇，本文采用文獻(xiàn)[16]中的方法，其文獻(xiàn)是依據(jù)VMD 分解之后各模態(tài)分量的中心頻率來確定是否出現(xiàn)了過分解現(xiàn)象，若當(dāng)前模態(tài)分量的中心頻率差值遠(yuǎn)小于其它K值下的差值時，則表示VMD 出現(xiàn)過分解現(xiàn)象，此時的K-1 就是最佳分解個數(shù). 另外VMD 分解的其它初始參數(shù)中，懲罰因子設(shè)為1 000, 噪聲容忍度設(shè)為0. 不同分解個數(shù)K下的中心頻率如表1 所示.

表1 不同K 值下VMD 分解的中心頻率Tab. 1 VMD decomposition center frequency at different K values

從表1 中可以看出K=7 時，模態(tài)分量中心頻率的最小差值僅為81×10-6，遠(yuǎn)小于K=2～6 時中心頻率的差值，所以認(rèn)為K=7 時，VMD 出現(xiàn)了過分解現(xiàn)象. 另外，分解個數(shù)過多也會增大計(jì)算規(guī)模，所以最終選擇IMF1的VMD 模態(tài)分解個數(shù)為6，分解結(jié)果如圖6 所示.

圖6 VMD 分解IMF1 結(jié)果Fig. 6 VMD decomposition result of IMF1

至此，原始DO 時間序列經(jīng)過二層分解后，得到 分 量imf1，imf2，···，imf13，其 排 列 熵 值 如 圖7所示. 由圖7 可知，經(jīng)過VMD 的分解顯著降低了IMF1序列的復(fù)雜度.

圖7 二次分解后各分量排列熵值Fig. 7 Permutation entropy of each component after quadratic decomposition

將經(jīng)過二層分解得到的分量輸入到麻雀搜索算法優(yōu)化的神經(jīng)網(wǎng)絡(luò)模型中，在預(yù)測過程中，使用連續(xù)3 d 的DO 數(shù)據(jù)預(yù)測第4 d 的DO 數(shù)據(jù)，得到各分量的預(yù)測值，然后疊加各個分量的預(yù)測值，即可得到最后的DO 時間序列預(yù)測質(zhì)量濃度.

2.4 模型對比分析與評價 為驗(yàn)證本文所提TDTSSA-BP 模型的有效性，共設(shè)置兩組對比實(shí)驗(yàn)進(jìn)行分析與評價. 實(shí)驗(yàn)所用麻雀搜索算法和粒子群算法的種群數(shù)均設(shè)為20，最大迭代次數(shù)均設(shè)為100，神經(jīng)網(wǎng)絡(luò)隱藏層神經(jīng)元個數(shù)設(shè)為5，其它網(wǎng)絡(luò)參數(shù)均做相同設(shè)置.

第1 組對比試驗(yàn)：用本文所提模型與BP、SSABP、 CEEMDAN-SSA-BP、 VMD-SSA-BP、 TDTPSO-BP 模型做對比，結(jié)果如圖8 所示. 從圖8 中可以看出，BP 模型和SSA-BP 模型與實(shí)際值相比存在預(yù)測結(jié)果滯后的現(xiàn)象，采用CEEMDAN、VMD或者二層分解的模型能及時反映數(shù)據(jù)的變化，有效克服了預(yù)測滯后現(xiàn)象，本文所提模型的預(yù)測結(jié)果更貼近真實(shí)值，其各模型誤差對比如表2 所示.

圖8 第1 組不同模型預(yù)測結(jié)果Fig. 8 Prediction results of the first group of different models

表2 不同模型精度對比Tab. 2 Accuracy comparison of different models

由表2 可知，CEEMDAN-SSA-BP 模型相比于SSA-BP 模型MAE、RMSE 和MAPE 分別減少了62%、53%和62%，R2增加了105%，VMD-SSA-BP模型相比于SSA-BP 模型MAE、 RMSE 和MAPE分別減少了56%、61%和56%，R2增加了114%. 表明使用CEEMDAN 或VMD 對DO 時間序列單次分解可以有效提高預(yù)測精度.

本文所提預(yù)測模型與CEEMDAN-SSA-BP 模型對比，MAE、RMSE 和MAPE 分別減少了30%、36%和30%，R2增加了9%，與VMD-SSA-BP 模型相比，MAE、RMSE 和MAPE 分別減少了39%、24%和40%，R2增加了4%，表明二次分解技術(shù)比單次分解更顯著地提升了模型的預(yù)測精度.

與TDT-PSO-BP 模型對比，本文模型MAE、RMSE 和MAPE 分別減少了11%、8%和10%，R2增加了1%. 圖9 為兩種模型的Mode1序列預(yù)測的收斂曲線對比. 表明在做DO 時間序列的預(yù)測時，SSA 比PSO 收斂速度更快且具有更強(qiáng)的尋優(yōu)能力.第2 組對比實(shí)驗(yàn)：用本文所提模型與ARIMA[5]、GA-BP[6]和EEMD-PSO-SVM[9]溶解氧預(yù)測模型做對比. 結(jié)果如圖10 所示，由圖10 可知，TDT-SSABP 模型預(yù)測曲線更貼近真實(shí)值，誤差對比如表3所示.

圖9 不同算法收斂曲線Fig. 9 Convergence curves of different algorithms

圖10 第2 組不同模型預(yù)測結(jié)果Fig. 10 Prediction results of the second group of different models

由表3 可知，本文所提模型與ARIMA 模型相比MAE、RMSE 和MAPE 分別減少了71%、66%和71%，R2增加了74%，與GA-BP 模型相比MAE、RMSE 和MAPE 分別減少了74%、71%和74%，R2增加了134%，與EEMD-PSO-SVM 模型相比，MAE、RMSE 和MAPE 分別減少了63%、57%和63%，R2增加了32%，表明了本文模型與其它已報(bào)道模型相比，預(yù)測精度有顯著提高.

表3 本文模型與現(xiàn)有模型精度對比Tab. 3 The accuracy of this model is compared with the existing models

3 結(jié)論

為克服河流溶解氧質(zhì)量濃度序列的非線性和不平穩(wěn)性造成的預(yù)測精度難以提高的問題，本文提出了TDT-SSA-BP 預(yù)測模型. 首先應(yīng)用二層分解技術(shù)充分削弱序列的非線性和不穩(wěn)定性，與單次分解技術(shù)相比，預(yù)測精度得到了顯著提升；其次本文提出使用尋優(yōu)能力更強(qiáng)的SSA 算法優(yōu)化神經(jīng)網(wǎng)絡(luò)的權(quán)值和閾值并與PSO 算法做對比，所提SSA-BP 溶解氧預(yù)測模型取得了更高的預(yù)測精度且收斂速度更快；最后通過與現(xiàn)有的其它溶解氧預(yù)測模型做對比，TDT-SSA-BP 模型預(yù)測精度最高，進(jìn)一步驗(yàn)證了本文所提模型的有效性. 在后續(xù)的研究中，考慮融合二層分解技術(shù)和深度學(xué)習(xí)預(yù)測模型，進(jìn)一步提高預(yù)測精度.