結合時間序列分解和神經網絡的河流溶解氧預測

盧毅敏，張紅

(1. 福州大學 空間數據挖掘與信息共享教育部重點實驗室，福建 福州 350108；2. 福州大學 地理空間信息技術國家地方聯合工程研究中心，福建 福州 350108；3. 數字中國研究院(福建)，福建 福州 350003)

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

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

1 研究方法

1.1 CEEMDAN分解法

DO時序數據具有明顯的非線性和非平穩性特征，任一時間可具有多種波動模式，為了提取DO時序數據潛在的變化特性、周期特征，以及長期趨勢，實現對DO序列時頻特征的充分挖掘、平穩化處理及降噪，引入CEEMDAN分解法，使噪聲殘留引起的重構誤差在分解階段疊加抵消.

進行I次實驗，通過對每次分解后的余量序列添加白噪聲，并進行EMD分解，可得第k個模態分量IMFk(t)及第k個余量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(·)為標準差算子，能使每個分解過程具有適當的信噪比；E(·)為分解算子.

當余量序列Rk(t)的極值點個數小于2時，結束分解，共得到k個代表DO序列不同時頻特征的模態分量，以及代表DO序列長期趨勢的最終余量R(t).原始DO時序數據S(t)可表示為

(3)

1.2 樣本熵

若對分解后的每個DO時頻特征都進行建模，會帶來預測誤差累積;若僅根據頻率特征將各特征重組為高頻、中頻、低頻3組，會丟失部分隱含的DO時序變化信息.樣本熵(SE)是通過計算時間序列的復雜度來衡量信號產生新模式的概率，抗噪能力強，采用較少的數據段即可得到穩定的熵值，可充分挖掘水環境系統復雜性[14].因此，文中以SE衡量DO各時頻特征的自相似性，熵值越大,時頻特征越復雜，特征序列的自相似性越小，所包含的特征信息與變化細節越重要，對DO預測結果影響越大，建模時越應當保留.將各DO時頻特征依序分別組成維數為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個數據開始，連續m+1個數據組成的向量,m取4.

樣本熵定義為

(4)

各時頻特征的樣本熵值計算式為

(5)

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

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

1.3 CS優化的Elman神經網絡

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

(6)

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

1.4 組合預測模型的構建

基于CEEMDAN分解、SE和CS-Elman動態神經網絡構建的河流溶解氧預測方法流程圖，如圖1所示.首先，使用CEEMDAN方法分解原始DO時序數據，實現對非線性、非平穩性DO序列的平穩化處理及降噪，充分挖掘DO時序信號中不同時間尺度的特征信息及噪聲，提取DO時序隨時間變化的波

圖1 河流溶解氧預測方法流程圖Fig.1 Flow chart of dissolved oxygen prediction in rivers

動特征、周期特征及長期趨勢.然后，通過計算樣本熵值衡量各時頻特征的自相似性，將熵值近似的特征重新組合為新序列，以減小計算量和誤差累積.同時，保留對DO預測結果具有重要影響的特征信息與變化細節.最后，對新序列分別構建CS優化的Elman預測模型，將預測值疊加，得到最終預測結果.

2 實驗分析

2.1 數據預處理

以福建省的晉江流域為研究區域，該流域位于北緯24°49′36″～25°35′13″，東經117°41′13″～118°41′49″，流域總面積為5 629 km2，河長為182 km，河道平均坡降為1.9%，是福建省第3大河流，泉州市境內第一大河流，也是福建省經濟最發達地區之一.晉江流域水資源是泉州市重要的飲用水源地，其水質優劣與泉州市經濟發展和人民生活水平密切相關.

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

2.2 河流溶解氧質量濃度預測

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

圖2 DO時間序列CEEMDAN分解結果Fig.2 CEEMDAN decomposition results of DO time series

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

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

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

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

通過計算時間序列的自相關系數確定模型結構，當滯后階數為N時，時間序列自相關系數較小，可以認為前N-1天的DO質量濃度對第N天的DO質量濃度影響最大，確定Elman模型輸入神經元個數為N-1，輸出神經元個數為1.對于重組序列Xt，循環將連續N天的數據分為一組，用前N-1天的DO數據預測第N天的DO質量濃度，即前N-1天的DO時序數據作為輸入，第N天的數據作為輸出，得到樣本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組作為訓練數據，后(T-N+1)/4組作為測試數據.N為3，共得到645組實驗數據，將前500組作為模型的訓練數據，后145組作為測試數據.

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

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

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

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

圖4 DO時間序列單步預測Fig.4 Single-step prediction of DO time series

2.3 模型對比分析與精度評價

為評估文中模型(CEEMDAN-SE-CS-Elman)的有效性，將其與其他傳統時間序列預測模型進行對比，結果如圖4所示.由圖4可知：文中模型的預測曲線更貼合真實曲線，預測精度更高.

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

表3 不同預測模型精度評價

2.4 模型應用結果

為了驗證文中模型的實用性，分別對晉江流域安溪南英和南安秋陽站點同期DO時序數據進行預測，預測的結果，如圖5所示.由圖5可知：預測值與實際值擬合較好，安溪南英站點EMA為0.20，EMPA為2.50%，ERMS為0.31，R2為0.928 1；南安秋陽站點EMA為0.17，EMPA為2.39%，ERMS為0.25，R2為0.902 5.

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

(a) 安溪南英 (b) 南安秋陽圖5 不同監測站點的預測結果Fig.5 Prediction results of different monitoring sites

3 結論

針對小流域基礎資料少和溶解氧指標值的隨機波動性與非平穩性造成預測精度難以提高的問題，提出了基于CEEMDAN分解、樣本熵和CS-Elman動態神經網絡的組合預測模型.以福建省晉江流域溶解氧數據作為實例驗證，得到以下2個結論.

1) 設計CEEMDAN和樣本熵相結合的溶解氧時間序列分解方法，實現對溶解氧數據的平穩化處理及降噪，充分挖掘溶解氧數據隱含的不同時間尺度特征，相較傳統分解方法，既保留了細節信息，又減小誤差累積，有助于解讀溶解氧指標隨時間變化的內在機理.

2) 對提取的溶解氧時頻特征，分別構建Elman動態神經網絡模型進行訓練和預測，采用全局搜索能力較強的CS算法優化模型，進一步提高了河流溶解氧預測精度.

結果表明，提出的CEEMDAN-SE-CS-Elman模型與其他時間序列預測模型相比，EMA，EMPA，ERMS及R2均有所提高，可為數據驅動的小流域短時河流溶解氧預測提供新方法.