仕 玉 治， 彭 勇， 周 惠 成

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

（1）給定參數(shù)m、τ、ε的合理取值區(qū)間，對(duì)時(shí)間序列進(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ù)，并驗(yàn)證數(shù)據(jù)序列的混沌特性，在RVM內(nèi)循環(huán)中，利用式（8）、（9）迭代估計(jì)超參數(shù)αi、σ，給定一個(gè)αmax的值，將對(duì)應(yīng)的權(quán)重值設(shè)定為0，并剔除對(duì)應(yīng)的Xi，所剩余的X即為相關(guān)向量，對(duì)應(yīng)的權(quán)重向量

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

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

選取南方兩水庫入庫月徑流時(shí)間序列作為研究實(shí)例，分別為水庫1的51a（1953－01～2003－12）和水庫2的48a（1958－01～2005－12）入庫月徑流時(shí)間序列.水庫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.每一個(gè)樣本序列分成3個(gè)子樣本，對(duì)于水庫1，將前41a資料作為模型的訓(xùn)練樣本，中間包含豐、平、枯年份的5a資料作為檢驗(yàn)樣本，與前41 a配合確定合理的模型參數(shù)，剩余5a不參加確定模型參數(shù)，純粹用于檢驗(yàn)確定模型的外推預(yù)報(bào)能力.同樣，對(duì)于水庫2，用前38a的序列點(diǎn)作為訓(xùn)練樣本，中間5a作為檢驗(yàn)樣本，剩余5a序列點(diǎn)作為外推預(yù)測(cè)樣本.對(duì)數(shù)據(jù)進(jìn)行規(guī)格化處理，采用PSO算法優(yōu)化各方法中的參數(shù)，取相對(duì)誤差的絕對(duì)值（Emar）、相關(guān)系數(shù)（R）、確定性系數(shù)R2和合格率（定量）作為預(yù)報(bào)結(jié)果的評(píng)價(jià)指標(biāo).先以水庫1為例進(jìn)行計(jì)算分析，以水庫2作進(jìn)一步的驗(yàn)證.

4.1 實(shí)例分析

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

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

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

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

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

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

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

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

圖8 外推預(yù)測(cè)期實(shí)測(cè)流量與RVM預(yù)報(bào)區(qū)間對(duì)比圖Fig.8 Comparison between observed and predicted interval hydrograph during validated period

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

4.2 實(shí)例驗(yàn)證分析

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

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

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

5 結(jié) 論

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

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

（3）在進(jìn)行中長期徑流預(yù)報(bào)應(yīng)用時(shí)，相關(guān)向量機(jī)模型的不足之處是模型參數(shù)和樣本序列均以正態(tài)概率分布函數(shù)進(jìn)行推理，但從模型計(jì)算的結(jié)果來看可用于中長期徑流預(yù)報(bào)，下一步將以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] 王 文，許武成.對(duì)水文時(shí)間序列混沌特征參數(shù)估計(jì)問題的討論[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ù)報(bào)中的應(yīng)用[J].水利學(xué)報(bào)，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