基于EMD-WDD-MK模型的瑪納斯河年徑流預測

閆國輝，喬長錄，陳伏龍

（石河子大學水利建筑工程學院，新疆石河子832000）

0 引 言

徑流預測是進行流域水資源優化配置的前提，但是目前還難以用物理成因分析法準確預報未來某一時段的徑流值。徑流過程是一個十分復雜的過程，涉及水文與氣象等多個因素，受隨機因素的影響，過程表現出非線性且非平穩的特點［1］。徑流預報的主要方法有傳統統計方法、概率論方法、模糊集理論等方法［2-4］。隨著計算機技術的進步，新的預測方法也逐漸出現，例如支持向量機、人工神經網絡以及模糊數學等。同時對徑流內部變化規律的掌握，以及對徑流序列中噪聲的處理，也有助于提高模型預測精度。目前常用的徑流序列分析方法有主成分分析法、小波分解方法、經驗模態分解法（Empirical Mode Decomposition，EMD）和變分模態分解法（Variational Modal Decomposition，VMD）［5-7］。新方法的產生，并不代表著舊方法的落伍。尋求不同方法的優點，并將不同方法進行耦合建模，構建具有更高預測精度的預測模型是目前徑流預測的一個新的研究方向。陳旭等［8］提出了基于EMD 和PSO 的耦合模型，對汾河上游的年徑流進行了預測，得到了較為滿意的結果。呂晗芳等［9］構建了VMD-LSSVM 模型，對上靜游等水文站的月徑流進行了預測，并與單一LSSVM 模型預測結果進行對比，提高了模型的預測性能。梁浩等［10］構建了多種模型耦合的預測模型，對渭河流域的徑流進行了預測，體現了"分解-合成"預測方法的適用性。但目前對EMD方法應用時，多數未考慮到分解得到的高頻分量中含噪聲信號過多，未對高頻分量進行處理。且目前基于EMD 與VMD 構建的耦合模型，多用于內陸河流，在西北干旱區河流的適用性和可行性需要進一步驗證。

在水資源配置與規劃越來越嚴格的今天，對規劃年河流來水量的預測變得更加重要，在缺水的西北干旱區這種情況尤為明顯。由于徑流預測研究中以西北干旱區內陸河為對象的研究較少，且西北地區較為缺水，需要對規劃年的徑流量有更高精度的預測。為了提高干旱區河流年徑流預測精度，使農業及城市配水有更可靠的依據，本文將EMD 與VMD 分別與加權馬爾可夫鏈進行耦合建模，并引入小波降噪，構建EMD-MK、EMD-WDD-MK 和VMD-MK 模型。以位于西北干旱區的瑪納斯河（以下簡稱瑪河）為例，進行該方法的應用研究，以期建立一種適合干旱區河流的年徑流預測模型，為瑪河流域水資源配置與規劃提供可供參考的依據。

瑪河位于新疆準噶爾盆地南部，干流全長324 km，流域面積33 500 km2，地處干旱區內陸，屬于典型的溫帶大陸性氣候，平均氣溫為4.7～5.7 ℃。該河發源于天山北坡，是天山北麓徑流量最大的河流，主要以冰川融水及降雨為補給，冰雪融水補給占46%，降雨占26%，降雪18%，最后流入尾閭瑪納斯湖［11，12］。

1 模型構建

1.1 經驗模態分解法（EMD）

Huang 等［13］于1998年提出了經驗模態分解法，該方法對處理非線性或非平穩時間序列方面具有較強自適應性。本質上，EMD 可以對非平穩時間序列進行處理，并逐漸將序列分解為不同頻率的曲線和趨勢，從而獲得不同特征尺度序列的一種方法。每個序列稱為本征模式函數（Intrinsic Mode Function，IMF），原始序列的趨勢或剩余量通常是最低頻率的IMF 分量。本征模態函數必須滿足2 個條件：①相對于時間軸來說數據表現出局部對稱性，即在任何時間點局部平均值都為零［14］；②數據系列中，極點和零點在數量上必須相等或者最多相差1 個。經驗模式分解的結果可表示為［15］：

式中：Ci(t)為對應原始序列中不同特征尺度的序列分量；rn為剩余分量，可以表現原始序列X（t）的趨勢變化。

1.2 變分模態分解法（VMD）

Dragomiretskiy 等［16］于2014年提出了一種處理非平穩時間序列的新方法，即變分模態分解法（VMD）。VMD在分解序列過程中，采用非遞歸和變分分解對序列進行處理，將序列分解成幾個分量和一個趨勢向量。

VMD將變分問題引入到時間序列的處理中，其關鍵內容是變分問題構建完之后，對問題的求解。在將原始徑流序列輸入算法中后，經過算法運算，將輸出K個模態函數uk（t），其目的是為了將模態的估計寬帶之和最小。分解中變分約束模型如下式所示［17］：

式 中：{uk}={u1，u2，…，uK}是模態函數的集合；{ωk}={ω1，ω2，…，ωK}是對應模態函數中心頻率的集合；δ（t）為狄克拉函數；?為內卷運算符。

1.3 馬爾可夫鏈預測的方法

馬爾可夫過程是一類隨機過程，其最基本的特征是“馬氏性”，也稱“無后效性”，即在系統（隨機過程）的“現在狀態”已知的條件下，其“將來狀態”與“過去狀態”無關。馬爾可夫鏈定義如下：

考慮隨機過程{xn，n∈T}，若對任意整數n∈T和任意i0，i1，…，in∈T，條件概率滿足：

也就是說此時狀態的概率只與前一個時刻所處的狀態相關，與其他時刻所處的狀態無關，那么過程{xn，n∈T}為馬爾可夫鏈［18］。

一般考慮齊次馬爾可夫鏈，即對任意的m，k∈T有：

式中：Pij(m；k)為系統在m時刻處在狀態i，經過k步轉移至狀態j的概率；Pij(k)為系統從狀態i，經k步轉移至狀態j的概率，此時轉移概率與初始時刻無關，k取1時，Pij(1)記為Pij。

1962年，Thomas和Fiering［19］第一次利用一階馬爾可夫鏈來模擬河川徑流，之后更多的學者將馬爾可夫鏈用于徑流預測中。McMillan 等［20］通過馬爾可夫鏈-蒙特卡羅采樣方式，建立了徑流預測模型。隨著科學技術的進步，以及計算機的產生與發展，馬爾可夫鏈開始與智能算法相結合進行徑流預測。Aksoy 等［21］基于馬爾可夫鏈和合成數據提出條件人工神經網絡模型對某干旱區月降水量進行了預測。國內也有很多學者進行馬爾可夫鏈預測徑流的研究工作，例如，韓貽強［22］基于傳統馬爾可夫鏈對大沽夾河年徑流狀態進行預測。沈冰等［23］基于改進的加權馬爾可夫鏈徑流預測模型對寶雞市的渭河年徑流量進行預測。藍永超等［24］結合時間序列分析和馬爾可夫鏈建立了時間序列-馬爾可夫鏈耦合模型，并對黃河上游的年徑流進行多年的預測。劉新有等［25］以狀態預測概率為權重，將傳統馬爾可夫鏈徑流預測模型得出的預測狀態加權求和，實現了徑流數值預測并對怒江干流進行了預測。對于時間離散且狀態離散的時間序列，隨機過程預測方法的運用變得越加廣泛，馬爾可夫鏈預測就是其中最重要的一種［26］。

1.4 EMD-WDD-MK 模型構建

EMD-WDD-MK 模型的預測方法步驟：①首先利用EMD 方法將年徑流序列進行分解，可分解得到幾個IMF 分量和一個剩余分量；②由于EMD 分解得到的IMF1-IMF3 分量為高頻含噪聲分量，所以使用小波閾值法進行降噪處理；③使用馬爾可夫鏈對高頻降噪后的IMF 分量以及低頻IMF 分量進行預測，預測中引入級別特征值，定義級別特征值H為：

式中：η為概率作用系數，取值通常為2 或4，取值越大，大概率的影響就會越大。依據概率最大隸屬度原則判斷分量的分級狀態后，由式（6）計算分量的預測值［25］：

式中：?為年徑流量的預測值；Ai、Bi分別為狀態i區間的上限值和下限值。

④將IMF 分量及剩余分量的預測結果進行重構，重構所得的量即為該預測模型的年徑流量預測值，然后進行誤差分析與計算。

1.5 模型評價指標

為了對比EMD-MK、EMD-WDD-MK 和EMD-VMD 三種模型的預測結果精度，本文選取了預報合格率（QR）、平均絕對誤差（MAE）、平均相對誤差（MAPE）和均方根誤差（RMSE）4 種指標判定模型的預測性能。計算公式如下：

式中：n為預測結果合格的個數；m為總預測結果個數；yi為實測年徑流量模型預測得出的年徑流量。

2 瑪河年徑流預測

2.1年徑流序列的分解

利用EMD 與VMD 對瑪河肯斯瓦特水文站1956-2012年的年徑流序列進行分解，分解結果如圖1 和圖2 所示。模型以兩種分解方式分解得到的1955-2012年分量為模型的率定期，以2013和2014年為模型預測結果的驗證期，通過對比3種模型在驗證期的預測結果精度，來判斷模型的預測水平。

圖2 瑪河1955-2012年年徑流量序列VMD分解Fig.2 VMD decomposition of annual runoff series of Mahe River from 1955 to 2012

由圖1 可知，年徑流量時間序列的IMF 分量表現了該時間序列在頻率、振幅和波長方面的變化。其中IMF1 分量具有頻率最高、振幅最大、波長最短的特點，且由IMF1 分量到IMF5 分量，分量的振幅逐漸減小、頻率逐漸降低、波長逐漸變長。同時，瑪河年徑流量的IMF1-IMF5 分量分別具有準3～6 a、準6～9 a、準8～10 a、準11～15 a 和準35 a 的波動周期。根據分解所得的趨勢分量可以看出瑪河1955-2008年的年徑流量一直處于上升趨勢，2008年以后，上升趨勢逐漸平緩。

圖1 瑪河1955-2012年年徑流量序列EMD分解Fig.1 EMD decomposition of annual runoff series of Mahe River from 1955 to 2012

與EMD 分解所得向量不同，圖2 中VMD 分解所得分量表現由IMF1 分量到IMF5 分量的振幅逐漸變大、頻率逐漸增高和波長逐漸縮短的特點。通過與EMD分解得到的趨勢向量相比，VMD 分解得到的趨勢分量在反映徑流變化過程中更為精確。圖2中趨勢分量在2008年以后的趨勢情況與圖1中的趨勢分量大致相同，都表現為上升趨勢逐漸平緩并有下降趨勢。

2.2 小波閾值降噪

由于人為因素以及自然因素的影響，采集整編得到的徑流資料往往含有隨機噪聲，這會造成徑流序列的信噪比降低，所以需要采用一些方法對原始徑流資料消去噪聲處理，提高資料序列信噪比，進而提高預測精度。

在EMD分解得到的年徑流序列的各分量中，高頻分量中所含的噪聲較多，所以選擇IMF1-IMF3分量進行小波閾值降噪處理。處理過程中，通過對比不同處理方法所得去噪序列的信噪比（SNR）和均方根誤差（RMSE）。在硬閾值、軟閾值和固定閾值降噪處理3 種處理方法中，選取SNR較高，RMSE較低的方法進行降噪處理。對IMF1～IMF3 分量的降噪處理結果如圖3所示。

圖3 IMF1-IMF3分量降噪結果圖Fig.3 Denoising results of IMF1-IMF3 components

2.3 馬爾可夫鏈預測

2.3.1 計算均值與標準差

以EMD 分解得到，并通過降噪后的IMF1 分量的預測為例進行預測求解。根據IMF1 分量的數據，用式（11）和（12）分別計算出徑流序列的均值和均方差，計算結果序列均值為xˉ=-6×106m3、序列標準差為S=8.8×107m3。

2.3.2 IMF1分量分級

如表1 所示，根據計算得到的均值和均方差，把IMF1 分量序列劃分低、偏低、平、偏平、高等5個狀態。

表1 狀態分級表Tab.1 Status classification table

通過狀態分級表，將IMF1 分量序列中的數據進行狀態劃分，結果如表2所示。

表2 IMF1分量序列及狀態Tab.2 IMF1 component sequence and status

2.3.3 建立狀態轉移概率矩陣

根據IMF1 分量序列中各年的狀態，可以得到不同步長的馬爾可夫鏈轉移概率矩陣。經整理計算，各步長狀態轉移概率矩陣如下：

2.3.4 計算自相關系數

利用式（13）、式（14）分別計算得IMF1 分量序列的各階自相關系數和歸一化權重，并將各階自相關系數歸一化后作為各步長的馬爾可夫鏈權重，結果如表3所示。

表3 各階自相關系數及歸一化Tab.3 Autocorrelation coefficients and normalization

2.3.5 狀態預測

根據IMF1分量1955-2012年狀態轉移矩陣對2013年IMF1分量狀態進行預測，結果見表4。

表4 IMF1分量狀態預測Tab.4 IMF1 component state prediction

由表4可得max{Pi，i∈E}= 0.884，該權重和所對應的i= 3，根據最大隸屬度原則，可以預測IMF1分量2013年值的狀態為狀態3。之后將IMF1分量2013年的預測值加入IMF1序列中，以1955-2013年的序列資料預測2014年的IMF1 分量狀態，狀態預測的結果見表5。

由表5可知max{Pi，i∈E}= 0.887，該權重和所對應的i= 3，根據最大隸屬度原則，可以預測2014年IMF1 分量為狀態3。

表5 IMF1分量2014年狀態預測Tab.5 Status prediction of IMF1 component in 2014

2.3.6 預測值計算

根據模糊集理論中的級別特征值，根據式（5）和式（6）分別求出級別特征值和該年的預測值。為了突出較大概率隸屬度的作用，本文中η取2，把表3和表4中的數據分別代入（5）和式（6）中，分別計算出H 和x?。經過計算2013年與2014年EMD 分解得到的IMF1分量的預測值別分別0.33和0.33，對每個分量進行重復計算，即可得出各分量的預測值，通過重組可得預測年份的年徑流量。

2.4 預測結果對比分析

對3 種模型中不同分量進行加權馬爾可夫鏈預測并重構，得出2013 和2014年3 種模型不同的預測值，計算結果如圖4所示。

圖4 瑪河2013和2014年年徑流量預測結果Fig.4 Annual runoff forecast results of Mahe River in 2013 and 2014

由圖4 可知，EMD-MK 模型所得2013年和2014年瑪河年徑流量預測結果分別為14.53 億m3和13.08 億m3，相對誤差分別為18.81%和15.04%。VMD-MK 模型的預測結果分別為13.58 億m3和15.33 億m3，相對誤差分別為11.04%和34.83%。EMD-WDD-MK 模型在2013年和2014年的年徑流量預測值與實測值的相對誤差分別為18.32%和11.87%。

結合圖4 與表6 可以看出，雖然VMD-MK 模型在2013年的預測結果相對誤差相較于其他模型較小，但是2014年的預測結果誤差達到了34.83%，使其合格率低于70%，超出了可以參考的誤差范圍。在合格率方面，EMD-MK 模型與EMD-WDD-MK模型在2013年和2014年的預測值與實測值相對誤差均小于30%，預測結果都為合格。在平均絕對誤差的對比中可以發現，引入小波降噪的EMD-WDD-MK 模型比EMD-MK 模型的平均絕對誤差減小了0.21。在平均相對誤差和均方根誤差的對比中，也體現出EMD-MK 模型優于VMD-MK 的耦合模型，且EMD-WDD-MK 模型的表現最優，其精度也是3 種模型中最高的。

表6 瑪河各模型年徑流預測結果Tab.6 Annual runoff prediction results of each model of Mahe River

3 結 論

（1）利用本文構建的EMD-MK、EMD-WDD-MK 和VMDMK 模型對瑪河2013年和2014年的年徑流量進行了預測。通過對比預測結果的誤差，可以得出EMD-WDD-MK 模型相較于其他模型預測更為準確的結論。

（2）利用EMD 與VMD 分解方法對年徑流量序列進行分解有助于了解徑流周期變化的特征，并且能將非線性、非平穩年徑流序列轉化為平穩時間序列，有助于提高預測的準確性。同時使用小波閾值法，對高頻含噪聲的IMF1-IMF3分量進行降噪處理，降低噪聲對預測結果的影響，進而提高預測精度。

（3）與普通的馬爾可夫預測相比，加權馬爾可夫鏈以歸一化各步長的自相關系數為權重，更能充分發揮資料系列自身的相關性，且具有更好的精度。同時文中使用均方差將資料序列進行分級，更能體現資料序列的自身特點，分區更合理。

（4）隨著實測數據的增加，時間序列不斷延長。自相關系數、狀態轉移矩陣和權重有更多資料的支撐，會發生某些變化，這種變化使得預測模型不斷完善，預測精度也在不斷地提高。所以在序列中加入新的實測數據是十分重要的，這樣可以為徑流預測的結果提供更多的資料支持。 □