閆國輝，喬長錄，陳伏龍

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

0 引 言

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

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

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

1 模型構(gòu)建

1.1 經(jīng)驗?zāi)B(tài)分解法（EMD）

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

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

1.2 變分模態(tài)分解法（VMD）

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

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

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

1.3 馬爾可夫鏈預(yù)測的方法

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

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

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

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

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

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

1.4 EMD-WDD-MK 模型構(gòu)建

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

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

式中：?為年徑流量的預(yù)測值；Ai、Bi分別為狀態(tài)i區(qū)間的上限值和下限值。

④將IMF 分量及剩余分量的預(yù)測結(jié)果進(jìn)行重構(gòu)，重構(gòu)所得的量即為該預(yù)測模型的年徑流量預(yù)測值，然后進(jìn)行誤差分析與計算。

1.5 模型評價指標(biāo)

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

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

2 瑪河年徑流預(yù)測

2.1年徑流序列的分解

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

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

由圖1 可知，年徑流量時間序列的IMF 分量表現(xiàn)了該時間序列在頻率、振幅和波長方面的變化。其中IMF1 分量具有頻率最高、振幅最大、波長最短的特點，且由IMF1 分量到IMF5 分量，分量的振幅逐漸減小、頻率逐漸降低、波長逐漸變長。同時，瑪河年徑流量的IMF1-IMF5 分量分別具有準(zhǔn)3～6 a、準(zhǔn)6～9 a、準(zhǔn)8～10 a、準(zhǔn)11～15 a 和準(zhǔn)35 a 的波動周期。根據(jù)分解所得的趨勢分量可以看出瑪河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 分解所得分量表現(xiàn)由IMF1 分量到IMF5 分量的振幅逐漸變大、頻率逐漸增高和波長逐漸縮短的特點。通過與EMD分解得到的趨勢向量相比，VMD 分解得到的趨勢分量在反映徑流變化過程中更為精確。圖2中趨勢分量在2008年以后的趨勢情況與圖1中的趨勢分量大致相同，都表現(xiàn)為上升趨勢逐漸平緩并有下降趨勢。

2.2 小波閾值降噪

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

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

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

2.3 馬爾可夫鏈預(yù)測

2.3.1 計算均值與標(biāo)準(zhǔn)差

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

2.3.2 IMF1分量分級

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

表1 狀態(tài)分級表Tab.1 Status classification table

通過狀態(tài)分級表，將IMF1 分量序列中的數(shù)據(jù)進(jìn)行狀態(tài)劃分，結(jié)果如表2所示。

表2 IMF1分量序列及狀態(tài)Tab.2 IMF1 component sequence and status

2.3.3 建立狀態(tài)轉(zhuǎn)移概率矩陣

根據(jù)IMF1 分量序列中各年的狀態(tài)，可以得到不同步長的馬爾可夫鏈轉(zhuǎn)移概率矩陣。經(jīng)整理計算，各步長狀態(tài)轉(zhuǎn)移概率矩陣如下：

2.3.4 計算自相關(guān)系數(shù)

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

表3 各階自相關(guān)系數(shù)及歸一化Tab.3 Autocorrelation coefficients and normalization

2.3.5 狀態(tài)預(yù)測

根據(jù)IMF1分量1955-2012年狀態(tài)轉(zhuǎn)移矩陣對2013年IMF1分量狀態(tài)進(jìn)行預(yù)測，結(jié)果見表4。

表4 IMF1分量狀態(tài)預(yù)測Tab.4 IMF1 component state prediction

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

由表5可知max{Pi，i∈E}= 0.887，該權(quán)重和所對應(yīng)的i= 3，根據(jù)最大隸屬度原則，可以預(yù)測2014年IMF1 分量為狀態(tài)3。

表5 IMF1分量2014年狀態(tài)預(yù)測Tab.5 Status prediction of IMF1 component in 2014

2.3.6 預(yù)測值計算

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

2.4 預(yù)測結(jié)果對比分析

對3 種模型中不同分量進(jìn)行加權(quán)馬爾可夫鏈預(yù)測并重構(gòu)，得出2013 和2014年3 種模型不同的預(yù)測值，計算結(jié)果如圖4所示。

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

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

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

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

3 結(jié) 論

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

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

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

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