結(jié)合時(shí)間序列分解和神經(jīng)網(wǎng)絡(luò)的河流溶解氧預(yù)測(cè)

盧毅敏，張紅

(1. 福州大學(xué) 空間數(shù)據(jù)挖掘與信息共享教育部重點(diǎn)實(shí)驗(yàn)室，福建 福州 350108；2. 福州大學(xué) 地理空間信息技術(shù)國(guó)家地方聯(lián)合工程研究中心，福建 福州 350108；3. 數(shù)字中國(guó)研究院(福建)，福建 福州 350003)

河流中的溶解氧(DO)是反映水質(zhì)狀況及自凈能力的重要指標(biāo)[1].高質(zhì)量濃度溶解氧有利于降解河流中各類污染物，并有效控制底泥釋放氮、磷和有機(jī)物，當(dāng)河流復(fù)氧速度遠(yuǎn)低于耗氧速度時(shí)，低質(zhì)量濃度溶解氧將導(dǎo)致需氧生物死亡及水質(zhì)惡化[2].因此，針對(duì)河流溶解氧預(yù)測(cè)的研究具有重要意義，為水質(zhì)管理和污染預(yù)警提供決策支持.

基于機(jī)理模型的河流溶解氧預(yù)測(cè)方法需要大量的基礎(chǔ)資料作為支撐[3]，而小流域資料貧乏，隨著機(jī)器學(xué)習(xí)與智能傳感器技術(shù)的快速發(fā)展，基于數(shù)據(jù)驅(qū)動(dòng)的非機(jī)理組合預(yù)測(cè)模型逐漸興起[4-7].但河流的動(dòng)態(tài)性、不確定性，以及繁雜性使得河流溶解氧隨著時(shí)間隨機(jī)變化，呈現(xiàn)非線性、非平穩(wěn)性特征[8]，預(yù)測(cè)精度難以提高.文獻(xiàn) [9-10]采用小波分解法進(jìn)行平穩(wěn)化處理及降噪，并取得較好的預(yù)測(cè)效果，但這些方法未對(duì)DO序列的不同時(shí)頻特征進(jìn)行深入挖掘，并且小波基函數(shù)依賴于人為選擇，給預(yù)測(cè)結(jié)果帶來(lái)一定主觀影響.經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)法很好地解決了這個(gè)問(wèn)題[11]，僅根據(jù)數(shù)據(jù)自身極值特點(diǎn)進(jìn)行分解,而具有自適應(yīng)噪聲的完整集成經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMDAN)法解決了EMD模態(tài)混疊問(wèn)題[12]，以及集合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)存在大量集成平均計(jì)算次數(shù)的問(wèn)題[13].因此，本文提出一種結(jié)合CEEMDAN分解和Elman動(dòng)態(tài)神經(jīng)網(wǎng)絡(luò)的河流溶解氧預(yù)測(cè)模型.

1 研究方法

1.1 CEEMDAN分解法

DO時(shí)序數(shù)據(jù)具有明顯的非線性和非平穩(wěn)性特征，任一時(shí)間可具有多種波動(dòng)模式，為了提取DO時(shí)序數(shù)據(jù)潛在的變化特性、周期特征，以及長(zhǎng)期趨勢(shì)，實(shí)現(xiàn)對(duì)DO序列時(shí)頻特征的充分挖掘、平穩(wěn)化處理及降噪，引入CEEMDAN分解法，使噪聲殘留引起的重構(gòu)誤差在分解階段疊加抵消.

進(jìn)行I次實(shí)驗(yàn)，通過(guò)對(duì)每次分解后的余量序列添加白噪聲，并進(jìn)行EMD分解，可得第k個(gè)模態(tài)分量IMFk(t)及第k個(gè)余量Rk(t)為

(1)

Rk(t)=rk-1(t)-IMFk(t).

(2)

式(1)，(2)中：I為500；εk-1取值為ε0[std((rk-1(t))/std(Ek-1(wi(t)))]，ε0取0.2，std(·)為標(biāo)準(zhǔn)差算子，能使每個(gè)分解過(guò)程具有適當(dāng)?shù)男旁氡龋籈(·)為分解算子.

當(dāng)余量序列Rk(t)的極值點(diǎn)個(gè)數(shù)小于2時(shí)，結(jié)束分解，共得到k個(gè)代表DO序列不同時(shí)頻特征的模態(tài)分量，以及代表DO序列長(zhǎng)期趨勢(shì)的最終余量R(t).原始DO時(shí)序數(shù)據(jù)S(t)可表示為

(3)

1.2 樣本熵

若對(duì)分解后的每個(gè)DO時(shí)頻特征都進(jìn)行建模，會(huì)帶來(lái)預(yù)測(cè)誤差累積;若僅根據(jù)頻率特征將各特征重組為高頻、中頻、低頻3組，會(huì)丟失部分隱含的DO時(shí)序變化信息.樣本熵(SE)是通過(guò)計(jì)算時(shí)間序列的復(fù)雜度來(lái)衡量信號(hào)產(chǎn)生新模式的概率，抗噪能力強(qiáng)，采用較少的數(shù)據(jù)段即可得到穩(wěn)定的熵值，可充分挖掘水環(huán)境系統(tǒng)復(fù)雜性[14].因此，文中以SE衡量DO各時(shí)頻特征的自相似性，熵值越大,時(shí)頻特征越復(fù)雜，特征序列的自相似性越小，所包含的特征信息與變化細(xì)節(jié)越重要，對(duì)DO預(yù)測(cè)結(jié)果影響越大，建模時(shí)越應(yīng)當(dāng)保留.將各DO時(shí)頻特征依序分別組成維數(shù)為m和m+1的向量序列Sm={s1,s2,…,st-m}和Sm+1={s1,s2,…,st-m}，其中，si={ui,ui+1,…,ui+m-1}(i=1,2,…,t-m)為IMFk(t)中第i個(gè)數(shù)據(jù)開(kāi)始，連續(xù)m+1個(gè)數(shù)據(jù)組成的向量,m取4.

樣本熵定義為

(4)

各時(shí)頻特征的樣本熵值計(jì)算式為

(5)

式(4),(5)中：r取0.1～0.2std(IMFk(t))具有較合理的統(tǒng)計(jì)特性，文中取0.15std(IMFk(t));Bm(r)，Bm+1(r)分別為向量Sj與模板向量Si的匹配概率.

將熵值近似的分量進(jìn)行合并，得到n個(gè)新的重組序列Xt={x1,x2,…,xt}(t=1,2,…,T).

1.3 CS優(yōu)化的Elman神經(jīng)網(wǎng)絡(luò)

DO時(shí)序數(shù)據(jù)具有延續(xù)性及自相關(guān)性，歷史時(shí)序數(shù)據(jù)對(duì)未來(lái)指標(biāo)值的變化具有很大影響.Elman動(dòng)態(tài)神經(jīng)網(wǎng)絡(luò)作為時(shí)間序列預(yù)測(cè)模型的一種，承接層可以記憶一定程度的DO歷史數(shù)據(jù)，并與當(dāng)前時(shí)刻DO值共同成為隱含層的當(dāng)前輸入.這可以表達(dá)歷史DO數(shù)據(jù)與未來(lái)DO數(shù)據(jù)間的時(shí)間延遲，捕捉DO數(shù)據(jù)的時(shí)間變化特征，具有適應(yīng)時(shí)變特性的能力[15-17].新興的布谷鳥(niǎo)搜索(CS)算法只有兩個(gè)參數(shù)(待優(yōu)化的初始權(quán)值與閾值組數(shù)n=25，新解的概率p=0.25)，通用性好，搜索速率快，不易陷入局部最小值，能有效克服神經(jīng)網(wǎng)絡(luò)存在的問(wèn)題[18-19]，并消除預(yù)測(cè)過(guò)程中隨機(jī)因素的干擾，提高預(yù)測(cè)精度.將誤差值定義為適應(yīng)度值fit，依據(jù)萊維(Levy)飛行原理尋找新解，采用偏好隨機(jī)游動(dòng)法替換該解，即

(6)

式(6)中:α>0為步長(zhǎng)縮放因子；?為逐點(diǎn)乘積運(yùn)算；levy(λ)為隨機(jī)搜索路徑，可以產(chǎn)生隨機(jī)步長(zhǎng)；R為(0,1)的隨機(jī)數(shù)；xt，p,xt，q表示t代的兩個(gè)隨機(jī)解.最終得到最小適應(yīng)度值fitmin及最優(yōu)解xbest.

1.4 組合預(yù)測(cè)模型的構(gòu)建

基于CEEMDAN分解、SE和CS-Elman動(dòng)態(tài)神經(jīng)網(wǎng)絡(luò)構(gòu)建的河流溶解氧預(yù)測(cè)方法流程圖，如圖1所示.首先，使用CEEMDAN方法分解原始DO時(shí)序數(shù)據(jù)，實(shí)現(xiàn)對(duì)非線性、非平穩(wěn)性DO序列的平穩(wěn)化處理及降噪，充分挖掘DO時(shí)序信號(hào)中不同時(shí)間尺度的特征信息及噪聲，提取DO時(shí)序隨時(shí)間變化的波

圖1 河流溶解氧預(yù)測(cè)方法流程圖Fig.1 Flow chart of dissolved oxygen prediction in rivers

動(dòng)特征、周期特征及長(zhǎng)期趨勢(shì).然后，通過(guò)計(jì)算樣本熵值衡量各時(shí)頻特征的自相似性，將熵值近似的特征重新組合為新序列，以減小計(jì)算量和誤差累積.同時(shí)，保留對(duì)DO預(yù)測(cè)結(jié)果具有重要影響的特征信息與變化細(xì)節(jié).最后，對(duì)新序列分別構(gòu)建CS優(yōu)化的Elman預(yù)測(cè)模型，將預(yù)測(cè)值疊加，得到最終預(yù)測(cè)結(jié)果.

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

2.1 數(shù)據(jù)預(yù)處理

以福建省的晉江流域?yàn)檠芯繀^(qū)域，該流域位于北緯24°49′36″～25°35′13″，東經(jīng)117°41′13″～118°41′49″，流域總面積為5 629 km2，河長(zhǎng)為182 km，河道平均坡降為1.9%，是福建省第3大河流，泉州市境內(nèi)第一大河流，也是福建省經(jīng)濟(jì)最發(fā)達(dá)地區(qū)之一.晉江流域水資源是泉州市重要的飲用水源地，其水質(zhì)優(yōu)劣與泉州市經(jīng)濟(jì)發(fā)展和人民生活水平密切相關(guān).

實(shí)驗(yàn)數(shù)據(jù)來(lái)源于福建省環(huán)保廳地表水質(zhì)自動(dòng)監(jiān)測(cè)站點(diǎn).晉江流域共有泉州石礱、安溪南英和南安秋陽(yáng)3個(gè)地表水質(zhì)監(jiān)測(cè)站點(diǎn)，選取各站點(diǎn)2017年3月27日至2019年1月2日的DO監(jiān)測(cè)數(shù)據(jù)作為研究數(shù)據(jù).每日自00:00開(kāi)始監(jiān)測(cè)，每間隔4 h采樣一次.由于監(jiān)測(cè)設(shè)備故障、網(wǎng)絡(luò)傳輸錯(cuò)誤等問(wèn)題，建模前必須對(duì)原始DO監(jiān)測(cè)數(shù)據(jù)進(jìn)行預(yù)處理.首先，依據(jù)GB 3838-2002《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》和箱線圖剔除異常值；其次，考慮到實(shí)驗(yàn)數(shù)據(jù)集缺失值很少，且DO數(shù)據(jù)一般短時(shí)間內(nèi)波動(dòng)較小，所以使用中值插補(bǔ)法對(duì)缺失值進(jìn)行插補(bǔ)；最后，計(jì)算日均值，分別得到的647條日監(jiān)測(cè)數(shù)據(jù)構(gòu)成完整的DO時(shí)間序列S(t)={s1,s2,…,st}.

2.2 河流溶解氧質(zhì)量濃度預(yù)測(cè)

DO時(shí)間序列CEEMDAN分解結(jié)果,如圖2所示.由圖2可知：原始DO時(shí)序數(shù)據(jù)表現(xiàn)出明顯的隨機(jī)性和非線性，IMF 1～I(xiàn)MF 3時(shí)頻特征起伏變化明顯，表明DO質(zhì)量濃度受到隨機(jī)影響；IMF 4～I(xiàn)MF 8時(shí)頻特征具有明顯的周期特征，可知DO時(shí)間序列具有季節(jié)性變化；余量時(shí)頻特征較為平緩，表明DO時(shí)間序列的長(zhǎng)期趨勢(shì)，各IMF分量迭代次數(shù)逐漸降低至0，重構(gòu)后相對(duì)百分比誤差達(dá)到10-13數(shù)量級(jí)，表明DO原始數(shù)據(jù)得到完全分解.

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

計(jì)算各IMF分量樣本熵，如表1所示.表1中：Δ為差值.由表1可知：各分量的樣本熵值整體表現(xiàn)為遞減趨勢(shì)，說(shuō)明各分量隨著波動(dòng)頻率降低，序列的復(fù)雜程度越小，隨機(jī)性也越?。籌MF 1熵值最大，和相鄰分量IMF 2的樣本熵差值為0.230，在所有相鄰樣本熵差值中為第2大，說(shuō)明IMF 1與IMF 2之間的相似性較小，兩者所包含的特征信息與變化細(xì)節(jié)的差異較大，所以IMF 1對(duì)DO預(yù)測(cè)結(jié)果影響較大，建模時(shí)將IMF 1單獨(dú)作為一個(gè)新的重組序列；IMF 2～I(xiàn)MF 4熵值較大，且IMF 2與IMF 3，IMF 3與IMF 4之間的樣本熵差值在所有差值中較大，說(shuō)明IMF 2～I(xiàn)MF 4所包含的特征信息與變化細(xì)節(jié)對(duì)預(yù)測(cè)結(jié)果均有一定影響，也分別單獨(dú)作為一個(gè)新的重組序列；IMF 5～I(xiàn)MF 7熵值較小，序列的復(fù)雜程度較低，且IMF 4與IMF 5之間的樣本熵差值在所有差值中最大，IMF 5與IMF 6，IMF 6與IMF 7之間的樣本熵差值最小，說(shuō)明IMF 4與IMF 5所包含的特征信息的相似性很小，不應(yīng)疊加，而IMF 5～I(xiàn)MF 7這3個(gè)分量具有一定的相似性，所以將這3個(gè)分量疊加作為一個(gè)新的重組序列；同理將余量與IMF 8疊加作為一個(gè)新的重組序列.

表1 各IMF分量樣本熵Tab.1 Sample entropy of each IMF component

根據(jù)文獻(xiàn)[20]計(jì)算各重組序列的平均周期，各重組序列的分量組成及其平均周期結(jié)果，各IMF分量重組結(jié)果，如表2所示.表2中：T為平均周期.由表2可知：各重組序列的平均周期呈現(xiàn)遞增趨勢(shì)；重組序列1～3為分解出的高頻序列部分，平均周期分別為6.55，8.44，12.27 d，表明短期內(nèi)DO影響因子隨機(jī)作用引起的DO質(zhì)量濃度波動(dòng)情況；重組序列4,5為分解出的中頻序列部分，平均周期分別為23.28，51.92 d，表明河流DO質(zhì)量濃度的月變化、季節(jié)性變化；重組序列6為分解出的低頻部分，平均周期為337.50 d，表明從長(zhǎng)期看，晉江流域DO的周期性是按年度變化的.

表2 各IMF分量重組結(jié)果Tab.2 Results of recombined of IMF components

通過(guò)計(jì)算時(shí)間序列的自相關(guān)系數(shù)確定模型結(jié)構(gòu)，當(dāng)滯后階數(shù)為N時(shí)，時(shí)間序列自相關(guān)系數(shù)較小，可以認(rèn)為前N-1天的DO質(zhì)量濃度對(duì)第N天的DO質(zhì)量濃度影響最大，確定Elman模型輸入神經(jīng)元個(gè)數(shù)為N-1，輸出神經(jīng)元個(gè)數(shù)為1.對(duì)于重組序列Xt，循環(huán)將連續(xù)N天的數(shù)據(jù)分為一組，用前N-1天的DO數(shù)據(jù)預(yù)測(cè)第N天的DO質(zhì)量濃度，即前N-1天的DO時(shí)序數(shù)據(jù)作為輸入，第N天的數(shù)據(jù)作為輸出，得到樣本P={p1,p2,…,pT-N+1}，其中,pi=[xi,xi+1,…,xi+N-1]T(i=1,2,…,T-N+1)，共T-N組，將前3(T-N+1)/4組作為訓(xùn)練數(shù)據(jù)，后(T-N+1)/4組作為測(cè)試數(shù)據(jù).N為3，共得到645組實(shí)驗(yàn)數(shù)據(jù)，將前500組作為模型的訓(xùn)練數(shù)據(jù)，后145組作為測(cè)試數(shù)據(jù).

將訓(xùn)練樣本P={p1,p2,…,pT-N+1}歸一化,并輸入CS-Elman模型，將最優(yōu)初始權(quán)值和閾值xbest賦予Elman神經(jīng)網(wǎng)絡(luò)，使用誤差反向傳播和梯度下降法對(duì)各個(gè)隱含層神經(jīng)元的權(quán)值系數(shù)進(jìn)行修正，直到訓(xùn)練誤差小于閾值.為保證模型高效、穩(wěn)定運(yùn)行，將Elman神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)速率設(shè)為0.1，訓(xùn)練目標(biāo)最小誤差設(shè)為0.000 1，隱含層輸出采用激活函數(shù)tansig處理，模型訓(xùn)練結(jié)束后，將測(cè)試樣本的DO數(shù)據(jù)輸入到訓(xùn)練好的模型中進(jìn)行預(yù)測(cè)，經(jīng)過(guò)線性傳遞函數(shù)purelin處理并反歸一化，得到各重組特征序列的預(yù)測(cè)結(jié)果，如圖3所示.由圖3可知：隨著DO數(shù)據(jù)波動(dòng)趨于平穩(wěn)，預(yù)測(cè)值與實(shí)際值越相吻合.

(a) 重組序列1 (b) 重組序列2

(c) 重組序列3 (d) 重組序列4

(e) 重組序列5 (f) 重組序列6圖3 重組序列單步預(yù)測(cè)結(jié)果Fig.3 Single-step prediction results of recombined sequences

圖4 DO時(shí)間序列單步預(yù)測(cè)Fig.4 Single-step prediction of DO time series

2.3 模型對(duì)比分析與精度評(píng)價(jià)

為評(píng)估文中模型(CEEMDAN-SE-CS-Elman)的有效性，將其與其他傳統(tǒng)時(shí)間序列預(yù)測(cè)模型進(jìn)行對(duì)比，結(jié)果如圖4所示.由圖4可知：文中模型的預(yù)測(cè)曲線更貼合真實(shí)曲線，預(yù)測(cè)精度更高.

精度評(píng)價(jià)結(jié)果,如表3所示.表3中：EMA為平均絕對(duì)誤差；EMPA為平均絕對(duì)百分誤差、ERMS為均方根誤差；R2為可決系數(shù).由表3可知：相較于其他經(jīng)典的時(shí)間序列預(yù)測(cè)模型(LSTM，ARIMA，Elman動(dòng)態(tài)神經(jīng)網(wǎng)絡(luò)預(yù)測(cè))，文中模型的精度更高，更適合作為晉江流域河流溶解氧預(yù)測(cè)的基準(zhǔn)模型；相較于未分解的單一Elman模型，文中模型的EMA提高0.17，EMPA提高2.60%，ERMS提高0.26，R2提高0.197 5，說(shuō)明采用分解法對(duì)溶解氧時(shí)序數(shù)據(jù)進(jìn)行平穩(wěn)化處理及降噪，提取溶解氧不同時(shí)頻特征，能顯著提高預(yù)測(cè)精度，且CEEMDAN分解比EMD，EEMD分解更有效；相較于CEEMDAN-CS-ELman,文中CEEMDAN-SE-CS-Elman模型EMA提高0.04，EMPA提高0.49%，ERMS提高0.05，R2提高0.031 2，說(shuō)明將樣本熵值近似的溶解氧時(shí)頻特征重組，能減小誤差累積，保留重要信息，有效提高河流溶解氧預(yù)測(cè)精度；與采用遺傳算法優(yōu)化的Elman神經(jīng)網(wǎng)絡(luò)模型(CEEMDAN-SE-GA-Elman)進(jìn)行對(duì)比，文中模型EMA提高0.02，EMPA提高0.30%，ERMS提高0.02，R2提高0.016 4，證明CS優(yōu)化算法的優(yōu)越性，能進(jìn)一步提高河流溶解氧預(yù)測(cè)精度.綜合分析，文中構(gòu)建的模型較其他模型預(yù)測(cè)效果更好，誤差評(píng)價(jià)指標(biāo)均為最優(yōu).

表3 不同預(yù)測(cè)模型精度評(píng)價(jià)

2.4 模型應(yīng)用結(jié)果

為了驗(yàn)證文中模型的實(shí)用性，分別對(duì)晉江流域安溪南英和南安秋陽(yáng)站點(diǎn)同期DO時(shí)序數(shù)據(jù)進(jìn)行預(yù)測(cè)，預(yù)測(cè)的結(jié)果，如圖5所示.由圖5可知：預(yù)測(cè)值與實(shí)際值擬合較好，安溪南英站點(diǎn)EMA為0.20，EMPA為2.50%，ERMS為0.31，R2為0.928 1；南安秋陽(yáng)站點(diǎn)EMA為0.17，EMPA為2.39%，ERMS為0.25，R2為0.902 5.

實(shí)驗(yàn)結(jié)果表明：文中模型具有一定的實(shí)用性，但對(duì)于突變數(shù)據(jù)(如安溪南英站點(diǎn)第133 d)，雖能準(zhǔn)確預(yù)測(cè)變化趨勢(shì)，但存在一定預(yù)測(cè)誤差，在今后的建模中，將針對(duì)研究提高突變數(shù)據(jù)的預(yù)測(cè)精度，考慮加入其他DO影響因素，使模型更加準(zhǔn)確、穩(wěn)定.

(a) 安溪南英 (b) 南安秋陽(yáng)圖5 不同監(jiān)測(cè)站點(diǎn)的預(yù)測(cè)結(jié)果Fig.5 Prediction results of different monitoring sites

3 結(jié)論

針對(duì)小流域基礎(chǔ)資料少和溶解氧指標(biāo)值的隨機(jī)波動(dòng)性與非平穩(wěn)性造成預(yù)測(cè)精度難以提高的問(wèn)題，提出了基于CEEMDAN分解、樣本熵和CS-Elman動(dòng)態(tài)神經(jīng)網(wǎng)絡(luò)的組合預(yù)測(cè)模型.以福建省晉江流域溶解氧數(shù)據(jù)作為實(shí)例驗(yàn)證，得到以下2個(gè)結(jié)論.

1) 設(shè)計(jì)CEEMDAN和樣本熵相結(jié)合的溶解氧時(shí)間序列分解方法，實(shí)現(xiàn)對(duì)溶解氧數(shù)據(jù)的平穩(wěn)化處理及降噪，充分挖掘溶解氧數(shù)據(jù)隱含的不同時(shí)間尺度特征，相較傳統(tǒng)分解方法，既保留了細(xì)節(jié)信息，又減小誤差累積，有助于解讀溶解氧指標(biāo)隨時(shí)間變化的內(nèi)在機(jī)理.

2) 對(duì)提取的溶解氧時(shí)頻特征，分別構(gòu)建Elman動(dòng)態(tài)神經(jīng)網(wǎng)絡(luò)模型進(jìn)行訓(xùn)練和預(yù)測(cè)，采用全局搜索能力較強(qiáng)的CS算法優(yōu)化模型，進(jìn)一步提高了河流溶解氧預(yù)測(cè)精度.

結(jié)果表明，提出的CEEMDAN-SE-CS-Elman模型與其他時(shí)間序列預(yù)測(cè)模型相比，EMA，EMPA，ERMS及R2均有所提高，可為數(shù)據(jù)驅(qū)動(dòng)的小流域短時(shí)河流溶解氧預(yù)測(cè)提供新方法.