基于相關(guān)向量機(jī)的中長期徑流預(yù)報模型研究

仕 玉 治， 彭 勇， 周 惠 成

（1.大連理工大學(xué) 水利工程學(xué)院，遼寧 大連 116024；2.山東省水利科學(xué)研究院，山東 濟(jì)南 250013）

0 引 言

中長期徑流預(yù)報對統(tǒng)籌安排防洪與抗旱、水庫調(diào)度與管理等事務(wù)，實現(xiàn)水資源的最大利用效益具有十分重要的實踐意義.由于水文系統(tǒng)本身的非線性及水文要素變化的不確定性，目前基于嚴(yán)密的物理方法還很難對徑流等水文現(xiàn)象進(jìn)行描述和預(yù)測，人們主要借助于成因分析法、水文統(tǒng)計法、模糊分析等方法來描述和預(yù)測水文過程.水文統(tǒng)計法依據(jù)水文資料的統(tǒng)計規(guī)律進(jìn)行預(yù)測，方法較為常用，它包括兩大類：一類是分析水文要素自身隨時間變化的統(tǒng)計規(guī)律，并建立模型進(jìn)行預(yù)測，如時間序列分析法；另一類是分析水文要素與多因子之間的相關(guān)關(guān)系，建立模型進(jìn)行預(yù)測，如多元回歸法.兩類均可直接利用原序列建立線性或非線性關(guān)系進(jìn)行預(yù)測，但其精度有時不能滿足工程需要.相空間重構(gòu)成為分析時間序列的一種嶄新的方法，通過挖掘或者恢復(fù)水文系統(tǒng)的多變量影響因子，重構(gòu)水文非線性動力系統(tǒng)，國外許多學(xué)者對短期徑流預(yù)報進(jìn)行研究，并取得了較好的應(yīng)用效果[1、2].時間滯時τ和嵌入維數(shù)m（文中也稱重構(gòu)參數(shù)）對時間序列的噪聲和數(shù)據(jù)量大小等影響因素比較敏感[3]，通常采用互信息法、關(guān)聯(lián)維數(shù)法、虛假近鄰法和Cao法等多種方法所得到的估計值差別較大，不利于獲得較好的預(yù)報精度.本文采用 Yu等[4]和Sivakumar[5]提出的方法，優(yōu)化得到重構(gòu)參數(shù).

一般線性方法難以描述水文系統(tǒng)非線性特征，許多新的方法逐步被引用到水文預(yù)報模型當(dāng)中，如貝葉斯理論、人工神經(jīng)網(wǎng)絡(luò)（ANN）、支持向量機(jī)（SVM）等[2、6、7]，進(jìn)一步發(fā)展了非線性徑流預(yù)報模型.2000年Tipping提出了一種新的稀疏概率模型相關(guān)向量機(jī)[8]（relevance vector machine，RVM），該方法用非線性核函數(shù)映射到高維空間，在高維空間進(jìn)行線性回歸，成功實現(xiàn)非線性向線性轉(zhuǎn)化，同時基于貝葉斯理論定義模型參數(shù)，不僅可以定量預(yù)報，而且能夠以概率分布的形式描述水文預(yù)報不確定度，可為水庫調(diào)度決策分析提供更多的可利用信息.目前，其已應(yīng)用到圖像分析[9、10]、信道均衡[11]等分類與回歸問題，獲得了較好的應(yīng)用效果.

綜上所述，確定自身前期影響因子和建立預(yù)報模型，是時間序列分析預(yù)測方法的關(guān)鍵，本文首先對徑流時間序列進(jìn)行相空間重構(gòu)，挖掘水文系統(tǒng)多變量因子；然后利用重構(gòu)后的時間序列建立RVM非線性徑流預(yù)報模型，并采用粒子群優(yōu)化（PSO）算法[12]辨識模型參數(shù)；最后應(yīng)用實例驗證本文模型的有效性.

1 混沌時間序列相空間重構(gòu)

相空間重構(gòu)是混沌理論的基礎(chǔ)，依據(jù)Takens理論[13]，對某一混沌時間序列｛xi：i＝1，2，…，n｝，只要適當(dāng)選取時間滯時τ和嵌入維數(shù)m，且嵌入維數(shù)滿足m≥2D＋1，其中D為飽和關(guān)聯(lián)維數(shù)，即可重構(gòu)與原未知動力系統(tǒng)具有相同幾何特征的m維相空間，則相空間中的點可以表示為Xi＝（xixi＋τxi＋2τ…xi＋（m－1）τ）T，由N個相點組成的延遲狀態(tài)向量表示為｛Xi：i＝1，2，…，N｝，其中N＝n－（m－1）τ，則相應(yīng)關(guān)聯(lián)積分表達(dá)式為

對于混沌時間序列，關(guān)聯(lián)積分C（r，m）與標(biāo)度尺度r近似成指數(shù)關(guān)系：C（r，m）∝rD，

2 相關(guān)向量機(jī)（RVM）徑流預(yù)報模型

已知相關(guān)向量｛Xi：i＝1，2，…，N｝，給定任意一個輸入向量X＊，則通過非線性映射到高維特征空間，然后在高維特征空間中進(jìn)行線性回歸得到預(yù)報輸出值，即RVM的模型輸出表示為

式中：K（·，·）為核函數(shù)，w為模型的權(quán)值.

式中：y＝ （y1y2…yN）；w＝ （w0w1…wN）；Φ為N×（N＋1）核函數(shù)矩陣，Φnn＝K（Xn，Xn），Φn1＝1.利用極大似然函數(shù)法估計w、σ時會導(dǎo)致嚴(yán)重的過擬合現(xiàn)象[8]，因此，為每個w定義高斯先驗概率分布函數(shù)

然后基于貝葉斯準(zhǔn)則計算權(quán)值的后驗概率分布，即

式中：后驗分布的協(xié)方差和均值分別為Σ＝

通過最大化邊緣似然分布函數(shù)

即可得超參數(shù)估計值，本文采用EM算法[8]內(nèi)循環(huán)迭代估計超參數(shù)α、σ，其αi和σ的迭代方程分別如下：

式中：γi＝1－αiΣii.

3 模型參數(shù)辨識及計算流程

RVM徑流預(yù)報模型需要解決兩個問題：（1）模型核函數(shù)選擇；（2）模型參數(shù)辨識過程中的目標(biāo)函數(shù)確立.對于核函數(shù)的選擇，線性核函數(shù)是徑向基核函數(shù)的特例，特定的sigmoid核函數(shù)功能上與徑向基核函數(shù)相同，核函數(shù)自身參數(shù)的個數(shù)太多不利于參數(shù)的選擇[7]，因此，本文選取徑向基核函數(shù)為核函數(shù)；其次，在模型參數(shù)識別過程中，通常選取訓(xùn)練樣本的擬合誤差最小為模型目標(biāo)函數(shù)，但是該方式下訓(xùn)練誤差收斂過程中會出現(xiàn)嚴(yán)重的過擬合現(xiàn)象，如圖1中擬合曲線1所示，訓(xùn)練階段擬合誤差非常小，幾乎接近于零，導(dǎo)致優(yōu)化參數(shù)不合理，外推期預(yù)報精度非常低.因此本文在訓(xùn)練過程中考慮具有豐、平、枯年份的檢驗樣本誤差來抑制過擬合，即綜合考慮訓(xùn)練樣本和檢驗樣本的誤差建立目標(biāo)函數(shù)，如下式：

式中：R1、R2分別為訓(xùn)練期、檢驗期的相對誤差絕對值的平均值，N1、N2分別為訓(xùn)練期、檢驗期的樣本個數(shù).其收斂過程如圖1中擬合曲線2所示，曲線2的收斂值要比曲線1的收斂值大，說明本文目標(biāo)函數(shù)具有抑制過擬合的能力.

最后，利用RVM模型進(jìn)行徑流預(yù)報的主要步驟如下：

（1）給定參數(shù)m、τ、ε的合理取值區(qū)間，對時間序列進(jìn)行相空間重構(gòu)；

圖1 不同目標(biāo)函數(shù)的訓(xùn)練誤差收斂過程Fig.1 Training error convergence process of different objective functions

（2）利用重構(gòu)后的訓(xùn)練樣本作為RVM的輸入進(jìn)行訓(xùn)練，采用PSO算法辨識方法參數(shù)，并驗證數(shù)據(jù)序列的混沌特性，在RVM內(nèi)循環(huán)中，利用式（8）、（9）迭代估計超參數(shù)αi、σ，給定一個αmax的值，將對應(yīng)的權(quán)重值設(shè)定為0，并剔除對應(yīng)的Xi，所剩余的X即為相關(guān)向量，對應(yīng)的權(quán)重向量

（3）依據(jù)優(yōu)化所得參數(shù)為模型的參數(shù)取值，給定任意一個輸入向量X＊，采用訓(xùn)練好的RVM進(jìn)行計算便可得到預(yù)報值的均值μ＊和方差，預(yù)報值服從均值μ＊和方差的后驗正態(tài)概率分布.

4 應(yīng)用實例分析

選取南方兩水庫入庫月徑流時間序列作為研究實例，分別為水庫1的51a（1953－01～2003－12）和水庫2的48a（1958－01～2005－12）入庫月徑流時間序列.水庫1的控制流域面積為11.45×104km2，多年平均月徑流量為1 251m3/s，變差系數(shù)為0.717，最大、最小月徑流量分別為5 000、248 m3/s；水庫2的控制流域面積為10.26×104km2，多年平均月徑流量為1 215m3/s，變差系數(shù)為0.856，最 大、最 小 月 徑 流 量 分 別 為 5 480、236 m3/s.每一個樣本序列分成3個子樣本，對于水庫1，將前41a資料作為模型的訓(xùn)練樣本，中間包含豐、平、枯年份的5a資料作為檢驗樣本，與前41 a配合確定合理的模型參數(shù)，剩余5a不參加確定模型參數(shù)，純粹用于檢驗確定模型的外推預(yù)報能力.同樣，對于水庫2，用前38a的序列點作為訓(xùn)練樣本，中間5a作為檢驗樣本，剩余5a序列點作為外推預(yù)測樣本.對數(shù)據(jù)進(jìn)行規(guī)格化處理，采用PSO算法優(yōu)化各方法中的參數(shù)，取相對誤差的絕對值（Emar）、相關(guān)系數(shù)（R）、確定性系數(shù)R2和合格率（定量）作為預(yù)報結(jié)果的評價指標(biāo).先以水庫1為例進(jìn)行計算分析，以水庫2作進(jìn)一步的驗證.

4.1 實例分析

RVM模型參數(shù)主要有相空間重構(gòu)參數(shù)（m、τ）、所選取的核函數(shù)自身參數(shù)及模型自動確定的超參數(shù)（α，σ）.給定一個較小的參數(shù)區(qū)間，進(jìn)行優(yōu)化計算，若優(yōu)化所得參數(shù)值為區(qū)間端點值，則進(jìn)一步擴(kuò)大區(qū)間重新計算，直至參數(shù)取值在區(qū)間范圍內(nèi)為止，即該區(qū)間為參數(shù)區(qū)間，由此確定徑向基核函數(shù)帶寬ε，混沌時間序列嵌入維數(shù)m、時間滯時τ的取值區(qū)間分別為[0.1，100]、[1，20]、[1，10]，另外對超參數(shù)初始化，取α（0）＝（0.25，0.25，…，

其次，對時間序列進(jìn)行相空間重構(gòu)，以重構(gòu)后的訓(xùn)練樣本作為徑流模型的輸入條件，采用PSO優(yōu)化模型參數(shù)（ε、m、τ），水庫1結(jié)果為（2.172 3，14，4），水庫2結(jié)果為（4.554，7，6）.根據(jù)優(yōu)化所得時間滯時，分別以嵌入維數(shù)m＝1～20繪制D－log2r折線斜率圖，如圖2所示，log2r在1～2，且m＞8時飽和關(guān)聯(lián)維數(shù)趨于穩(wěn)定值，即存在顯著標(biāo)度區(qū)，從而可以定性判斷兩水庫月徑流序列存在混沌特性，并由圖2（a）、（b）可以估計出1＜D＜3，重構(gòu)參數(shù)m滿足m≥2D＋1的條件，說明重構(gòu)參數(shù)是合理的.此外，在參數(shù)內(nèi)循環(huán)過程中，隨著超參數(shù)的迭代估計，邊緣似然分布函數(shù)值逐步趨于穩(wěn)定，如圖3所示；由圖4知超參數(shù)σ2收斂很快，迭代3次后基本達(dá)到最優(yōu)值，因此有的文獻(xiàn)也將超參數(shù)σ2作為一個固定值進(jìn)行內(nèi)循環(huán)計算.

圖2 水庫1和2月徑流關(guān)聯(lián)維數(shù)圖Fig.2 Correlation dimension of monthly runoff flow of Reservoirs One and Two

圖3 log2（p（y｜α，σ2））的收斂過程Fig.3 Convergence process of log2（p（y｜α，σ2））

圖4 超參數(shù)σ2的收斂過程Fig.4 Convergence process of hyper－parameterσ2

分析模型的模擬、檢驗、外推預(yù)報精度，并將本文模型（RVM）與應(yīng)用較為廣泛的最小二乘支持向量機(jī)模型（LSSVM），以及未考慮相空間重構(gòu)的（m＝12，τ＝1）相關(guān)向量機(jī)模型（RVM＊）和自動回歸滑動平均模型（ARMA（5，6））進(jìn)行對比分析.計算結(jié)果列于表1中，RVM方法對于訓(xùn)練期、檢驗期和外推預(yù)測期的預(yù)報結(jié)果見圖5～8.

總體而言，由表1知，考慮相空間重構(gòu)進(jìn)行預(yù)報時比未考慮相空間重構(gòu)時，RVM獲得比單一方法更優(yōu)的預(yù)報精度，與LSSVM和ARMA（5，6）的計算結(jié)果相比較，RVM的評價指標(biāo)值均優(yōu)于其相應(yīng)值，說明RVM具有較好的預(yù)報性能.按照《水文情報預(yù)報規(guī)范》（SL 250—2000）標(biāo)準(zhǔn)將徑流量劃分為枯、偏枯、平、偏豐、豐5個級別，對高流量（包括偏豐和豐流量）精度進(jìn)行了定量、定性比較分析，結(jié)果如表2所示.RVM在訓(xùn)練期、檢驗期及外推預(yù)測期的平均絕對相對誤差分別比LSSVM和ARMA的相應(yīng)值要小，但同時其預(yù)報精度均比預(yù)報總體時相應(yīng)值低，以多年變幅的20%為許可誤差，比較分析知，其定量合格率較本文所列其他方法有所提高.為提供更為充分的預(yù)報信息，本文對比分析了三階段高流量的定性預(yù)報合格率，除了在檢驗期RVM和LSSVM的合格率相同以外，其余兩階段RVM均獲得比其他方法更高的定性預(yù)報合格率，同樣說明RVM具有較強(qiáng)的高流量預(yù)報能力.

表1 水庫1月流量不同方法預(yù)測精度Tab.1 Prediction accuracy of monthly flow of Reservoir One resulting from various methods

圖5 訓(xùn)練期流量實測值與RVM預(yù)報值對比圖及相關(guān)圖Fig.5 Comparison and scatter plot between observed flow and predicted flow by RVM during training period

圖6 檢驗期流量實測值與RVM預(yù)報值對比圖及相關(guān)圖Fig.6 Comparison and scatter plot between observed flow and predicted flow by RVM during test period

圖7 外推預(yù)測期流量實測值與RVM預(yù)報值對比圖及相關(guān)圖Fig.7 Comparison and scatter plot between observed flow and predicted flow by RVM during validated period

圖8 外推預(yù)測期實測流量與RVM預(yù)報區(qū)間對比圖Fig.8 Comparison between observed and predicted interval hydrograph during validated period

進(jìn)一步考慮徑流預(yù)報的不確定性，以預(yù)報值的均值和方差為預(yù)報的后驗概率分布函數(shù)來描述預(yù)報值的不確定性，并討論了發(fā)生概率為80%的區(qū)間預(yù)報，其區(qū)間預(yù)報結(jié)果及實測流量過程如圖8所示.由圖8知，中低流量預(yù)報區(qū)間基本上可以包住實測流量，高流量區(qū)間上下限值對應(yīng)的級別能夠預(yù)報出實測值對應(yīng)的級別，概率區(qū)間預(yù)報是可靠的.

4.2 實例驗證分析

水庫2的統(tǒng)計參數(shù)與水庫1基本相同，但是變差系數(shù)較大，數(shù)據(jù)序列平穩(wěn)性相對較差，在同樣可行條件下，對水庫2進(jìn)行了計算，其預(yù)報結(jié)果的評價指標(biāo)列于表3.由表3知，RVM較LSSVM和ARMA（6，6）模型具有較高的預(yù)報精度，驗證說明了本文模型的有效性.

表2 水庫1高月流量不同方法預(yù)測精度Tab.2 Prediction accuracy of high monthly flow of Reservoir One resulting from various methods

表3 水庫2月流量不同方法預(yù)測精度Tab.3 Prediction accuracy of monthly flow of Reservoir Two resulting from various methods

5 結(jié) 論

（1）將混沌技術(shù)與相關(guān)向量機(jī)結(jié)合建立徑流預(yù)報模型，采用PSO算法辨識模型參數(shù)，優(yōu)化所得重構(gòu)參數(shù)滿足混沌理論條件，耦合方法比單一方法的預(yù)報精度有所提高，并對總體和高流量值進(jìn)行分析，取得比LSSVM和ARMA模型更優(yōu)的預(yù)報精度，說明本文模型的有效性.

（2）相關(guān)向量機(jī)為概率模型，能夠定量地、以概率分布的形式描述徑流預(yù)報不確定性，并給出指定發(fā)生概率下的區(qū)間預(yù)報.

（3）在進(jìn)行中長期徑流預(yù)報應(yīng)用時，相關(guān)向量機(jī)模型的不足之處是模型參數(shù)和樣本序列均以正態(tài)概率分布函數(shù)進(jìn)行推理，但從模型計算的結(jié)果來看可用于中長期徑流預(yù)報，下一步將以P－Ⅲ型概率分布函數(shù)進(jìn)行模型研究.

[1] SIVAKUMAR B. Chaos theory in hydrology important issues and interpretations[J].Journal of Hydrology，2000，227：1－20

[2] SIVAKUMAR B， JAYAWARDENA A W，F(xiàn)ERNANDO T M K G.River flow forecasting：use of phase－space reconstruction and artificial neural networks approaches [J].Journal of Hydrology，2002，265：225－245

[3] 王 文，許武成.對水文時間序列混沌特征參數(shù)估計問題的討論[J].水科學(xué)進(jìn)展，2005，16（4）：606－610

[4] YU X Y，LIONG S Y，BABOVIC V.EC－SVM approach for real time hydrologic forecasting [J].Journal of Hydroinformatics，2004，6（3）：209－223

[5] SIVAKUMAR B.Nonlinear determinism in river flow prediction as a possible indicator [J].Earth Surface Processes and Landforms，2007，32：969－979

[6] LIONG S Y，SIVAPRAGASAM C.Flood stage forecasting with SVM [J].Journal of the American Water Resources Association，2002，38（1）：173－186

[7] 林劍藝，程春田.支持向量機(jī)在中長期徑流預(yù)報中的應(yīng)用[J].水利學(xué)報，2006，37（6）：681－686

[8] TIPPING M E.The relevance vector machine[J].Advances in Neural Information Processing System，2000，12：652－658

[9] AGARWAL A，TRIGGS B.3Dhuman pose from Silhouettes by relevance vector regression[J].Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition，2004，2：882－888

[10] BOWD C，MEDEIROS F A，ZHANG Zuo－h(huán)ua，etal.Relevance vector machine and support vector machine classifier analysis of scanning laser polarimetry retinal nerve fiber layer measurements[J].Investigative Ophthalmology ＆ Visual Science，2005，46：1322－1329

[11] CHEN S，GUNN S R，HARRIS C J.The relevance vector machine technique for channel equalization application [J].IEEE Transactions on Neural Networks，2002，12（6）：1529－1532

[12] KENNEDY J，EBERHART R C.Particle swarm optimization[C]//Proceedings of IEEE Conference on Neural Networks.Piscataway：IEEE Press，1995：1942－1948

[13] KANTZ H，SCHREIBER T.Nonlinear Time Series Analysis [M].Cambridge：Cambridge University Press，1997