盧毅敏，張紅

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

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

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

1 研究方法

1.1 CEEMDAN分解法

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

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

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

(3)

1.2 樣本熵

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

樣本熵定義為

(4)

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

(5)

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

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

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

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

(6)

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

1.4 組合預測模型的構建

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

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

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

2 實驗分析

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

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

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

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

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

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

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

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

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

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

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

(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)的有效性，將其與其他傳統(tǒng)時間序列預測模型進行對比，結果如圖4所示.由圖4可知：文中模型的預測曲線更貼合真實曲線，預測精度更高.

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

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

2.4 模型應用結果

為了驗證文中模型的實用性，分別對晉江流域安溪南英和南安秋陽站點同期DO時序數(shù)據(jù)進行預測，預測的結果，如圖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.

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

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

3 結論

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

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

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

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