陶志勇，胡啟振，任曉奎

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

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

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

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

1 模型介紹

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

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

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

步驟 2 采用CEEMDAN 算法對(duì)數(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ù)測(cè)模型.

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

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

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

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

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

1.2 二層分解技術(shù)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

步驟 5 保存擁有最佳適應(yīng)度值的麻雀位置，若達(dá)到最大迭代次數(shù)時(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ù)來(lái)源 本文數(shù)據(jù)來(lái)源于2016 年10 月21日至2017 年11 月21 日江蘇省無(wú)錫市長(zhǎng)江水質(zhì)實(shí)時(shí)監(jiān)測(cè)站，監(jiān)測(cè)頻率為每日一次，共計(jì)397 個(gè)數(shù)據(jù)，經(jīng)預(yù)處理后，數(shù)據(jù)集如圖3 所示. 從圖3 中可以看出，DO 時(shí)間序列具有明顯的隨機(jī)性和不平穩(wěn)性.選 取2016 年10 月21 日 至2017 年10 月15 日 共計(jì)360 個(gè)數(shù)據(jù)作為訓(xùn)練集，剩余2017 年10 月16日至11 月21 日共37 個(gè)數(shù)據(jù)作為測(cè)試集.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表2 不同模型精度對(duì)比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 對(duì)DO 時(shí)間序列單次分解可以有效提高預(yù)測(cè)精度.

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

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

圖9 不同算法收斂曲線(xiàn)Fig. 9 Convergence curves of different algorithms

圖10 第2 組不同模型預(yù)測(cè)結(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ù)測(cè)精度有顯著提高.

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

3 結(jié)論

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