賀玉琪，王 棟，王遠坤

（南京大學 地球科學與工程學院，江蘇 南京 210023）

1 研究背景

受地形和氣候等諸多因素影響，降水時間序列多由趨勢項、周期項和隨機項混合疊加而成，兼有長短期變化規(guī)律[1]。傳統(tǒng)的時間序列預測方法如簡單平均法、指數(shù)平滑法等對數(shù)據(jù)的平穩(wěn)性有較高要求，不適用于預測非平穩(wěn)水文時間序列。神經(jīng)網(wǎng)絡（ANN）方法常用于水文時間序列分析[2-5]，但在設(shè)置復雜網(wǎng)絡參數(shù)和結(jié)構(gòu)的固定學習速率下，網(wǎng)絡收斂速度較慢[3]，訓練容易陷入局部最小值，導致學習失敗[4]，這些不足限制了ANN的應用。支持向量機回歸（Support Vector Regression，SVR）也常用于時間序列的研究[1，6-10]，SVR遵循結(jié)構(gòu)風險最小化原則[7]，其拓撲結(jié)構(gòu)由支持向量決定，避免了拓撲結(jié)構(gòu)難以確定的問題[8]，最終化為求解線性約束的凸二次規(guī)劃優(yōu)化問題[11]，因此可得全局最優(yōu)解。貝葉斯嶺回歸（Bayesian Ridge Regression，BRR）方法假設(shè)參數(shù)服從高斯分布，在估計過程中自動引入正則項，解決了過擬合問題[12]，且BRR適合解決小樣本高維度問題，明顯優(yōu)于ANN方法。BRR方法曾被用于降水資料時間序列的構(gòu)建，并體現(xiàn)一定的優(yōu)越性[13]。

小波變換可將時間序列分解到不同的時間尺度，反映數(shù)據(jù)在不同分辨率水平上的變化規(guī)律[15]，已被廣泛用于水文過程的多尺度分析研究[15-18]。以往研究中，小波分解后的各子序列采用相同模型，由于各序列的數(shù)據(jù)變化規(guī)律不同，在一些序列上SVR的精度較高，而BRR更適合另一些序列。本文針對這一點，提出不同子序列采用不同模型進行回歸，利用校驗數(shù)據(jù)從SVR和BRR中選擇對各子序列預測精度更高的模型，構(gòu)建BRR-SVR優(yōu)化模型，根據(jù)國內(nèi)外文獻調(diào)研結(jié)果，該優(yōu)化模型尚未用于水文及其他領(lǐng)域的研究。相比于單一的BRR和SVR模型，BRR-SVR模型能夠充分利用小波變換的優(yōu)勢，由于最終的預測結(jié)果是由不同序列的預測結(jié)果疊加而得，優(yōu)化模型保證了每個子序列上的預測誤差最小，使得整體預測精度最大。

2 基本方法

2.1 小波變換小波是一種長度有限、平均值為零的波形。小波函數(shù)通過伸縮與平移變換形成靈活可變的時頻窗進行局部化分析[5]。小波變換分為連續(xù)小波變換和離散小波變換，水文序列多為觀測的離散值[17]，因此采用離散小波變換。設(shè)離散信號f(n·Δt)為一個平方可積函數(shù)，其離散小波變換為：

式中：W f(a，b)為小波變換系數(shù)；a為尺度因子；b為平移因子；Ψ(t)為小波母函數(shù)；Ψ-- -- (-t-)為Ψ(t)的復共軛函數(shù)；Δt為采樣時間間隔；n為樣本總數(shù)。

對原始序列X，使用Mallat算法進行分解和重構(gòu)，Mallat算法是一種針對小波系數(shù)進行逐層分解和重構(gòu)的小波變換算法[19]。原序列經(jīng)第一次分解得到一個低頻序列a1和高頻序列d1，之后的每層分解均是將上一層的低頻信號aj-1分解成一個低頻序列aj和高頻序列dj，長度減半。對最大分解水平上的低頻序列ak和各高頻序列d（jj=1，2，…，k）通過小波逆變換進行重構(gòu)[5]，得到重構(gòu)序列Ak和D（jj=1，2，…，k）。各重構(gòu)序列的長度與原序列相等，原序列可由各重構(gòu)序列疊加而得，對于序列X有：

2.2 支持向量機回歸（SVR）支持向量機（SVM）是根據(jù)結(jié)構(gòu)風險最小化原理提出的一種機器學習算法，其基本思想是用核函數(shù)將原樣本空間中的非線性回歸轉(zhuǎn)化為高維空間中的線性回歸[9]。SVM 包括支持向量機分類和支持向量機回歸，SVR的原理如下。

給定樣本{(xi，yi)，i=1，2，…，m} ∈R n，m為樣本個數(shù)，n為樣本維度，回歸函數(shù)為：

式中：φ(x)為R n到高維特征空間的非線性映射；ω為超平面的權(quán)值向量；b為偏置。

根據(jù)結(jié)構(gòu)風險最小化原則，原問題轉(zhuǎn)化為在式（5）所示的約束條件下由式（4）求出最小值：

式中：常數(shù)C為誤差懲罰因子；m為樣本容量；ξi、ξ*i為松弛變量；ε為不敏感損失系數(shù)。

引入拉格朗日乘子，可以得到上述優(yōu)化問題的對偶問題：

式中：αi、αi*、αj、αj*為拉格朗日乘子；K(xi，xj)為核函數(shù)。由于最后化為一個線性約束的凸優(yōu)化問題，所以解具有全局最優(yōu)性和唯一性[20]。

2.3 貝葉斯嶺回歸（BRR）BRR是以貝葉斯理論為基礎(chǔ)的一種機器學習回歸算法[21]。貝葉斯線性回歸的函數(shù)見式（8），其目的是求使得損失函數(shù)式（9）最小的參數(shù)向量分布。

式中：n為樣本空間維度；m為樣本容量；ω為參數(shù)向量；Ψ(x)為輸入向量x的非線性函數(shù)，其中Ψ0(x)=1；ti為觀測值，ti =y(xi，ω)+ε，ε為噪聲，假定ε和ω分別服從高斯分布和則t是以y(x，ω)為均值的高斯分布，t的類條件概率密度函數(shù)為式（10），ω的先驗概率為式（11）：

依據(jù)貝葉斯規(guī)則有：

式中：p(ω|t)為后驗概率；p(t)為與ω無關(guān)的常數(shù)；c為常數(shù)，先驗概率對應嶺回歸中的L2 正則項，因此稱為貝葉斯嶺回歸。

BRR在估計過程中自動引入正則項，最終得到的是參數(shù)的后驗分布，避免了極大似然估計中的過擬合，得到更精確的參數(shù)估計。

3 模型構(gòu)建

模型的構(gòu)建分為數(shù)據(jù)預處理、模型的訓練和校驗，具體流程見圖1。在模型的校驗階段，需要評價模型的預測性能，本文采用確定系數(shù)R2、平均絕對百分比誤差MAPE和平均絕對誤差MAE這3項指標作為評價模型性能的依據(jù)，計算公式如下：

式中：yi為預測值；yi為觀測值；y為觀測值的平均值；m為樣本總數(shù)。R2越接近1，預測值與實測值越吻合，MAPE和MAE越小預測精度越高。

3.1 數(shù)據(jù)預處理不同月份的降水數(shù)據(jù)存在極端值和較多噪音，波動頻繁，經(jīng)過中心化處理可避免異常值和極端值的影響，減小誤差，加快模型訓練速度。中心化公式為：

式中：xt為t時刻降水量，mm；μ為整個時序降水均值，mm；σ為對應的標準差。

對中心化處理后的數(shù)據(jù)進行小波分解和重構(gòu)。本文選用Daubechies 小波族的“db3”正交小波和Mallat算法對數(shù)據(jù)做5層分解與重構(gòu)，得到一組低頻序列A5和5組高頻序列D（jj=1，2，…，5）。

為了適應支持向量機回歸和貝葉斯嶺回歸的分析方式[22]，對各子序列進行相空間重構(gòu)，將一維序列轉(zhuǎn)換為二維矩陣。將前m月的降水量作為輸入，下個月的降水量作為輸出。時間序列{xt，t=1，2，…，n}經(jīng)過相空間重構(gòu)后可得樣本：Xi ={xi，xi+1，…，xi+m-1}，Yi =xi+m，i=1，2，…，n-m，樣本量為（n-m）。

對6組子序列做相同的相空間重構(gòu)，將重構(gòu)后的樣本進行劃分，按時間先后順序分為訓練數(shù)據(jù)、校驗數(shù)據(jù)和測試數(shù)據(jù)。訓練數(shù)據(jù)用于模型學習歷史數(shù)據(jù)的規(guī)律；校驗數(shù)據(jù)用于參數(shù)的調(diào)整和選擇最優(yōu)模型；測試數(shù)據(jù)用于評估模型的預測性能。訓練數(shù)據(jù)和校驗數(shù)據(jù)都參與了模型的構(gòu)建，而測試數(shù)據(jù)對于模型是全新的數(shù)據(jù)。

3.2 模型的訓練和校驗由2.1節(jié)可知，最終的預測結(jié)果由6組子序列的預測結(jié)果疊加而得，若要整體的預測精度最大化，需要對6 組子序列分別建模，確保每個子模型的預測誤差最小。用校驗數(shù)據(jù)選取最佳模型，以R2、MAPE和MAE這3項指標來評估模型的精度，遵循少數(shù)服從多數(shù)的原則，當有兩項及以上指標的結(jié)果都說明一個模型的精度最高，該模型即為最佳模型。分析發(fā)現(xiàn)，用3項指標與僅用確定系數(shù)R2得到的結(jié)果基本一致，為了提升計算速率，僅用R2作為評價指標。3.2.1 BRR模型 貝葉斯嶺回歸中主要的超參數(shù)有正則化參數(shù)α1、α2和λ1、λ2，設(shè)置4 個超參數(shù)取值范圍均為{10-9，10-8，…，10-1}，不同取值的超參數(shù)構(gòu)成多組參數(shù)組合，用6組訓練數(shù)據(jù)集在各種參數(shù)組合下訓練模型，再用相應的校驗數(shù)據(jù)檢驗所有模型，以R2為評價指標，R2最大時對應的參數(shù)組合即為最佳參數(shù)，由此確定6個子序列對應的最佳BRR模型。3.2.2 SVR模型 支持向量機回歸的主要超參數(shù)有C、核函數(shù)和多項式的階數(shù)degree，設(shè)置C的取值范圍為{0.1，1，10，50，100}，degree的取值范圍為{1，2，3，…，10}，核函數(shù)的取值范圍為{‘rbf’，‘sigmoid’，‘poly’}。不同取值的超參數(shù)構(gòu)成多組參數(shù)組合。在各種參數(shù)組合下訓練模型，用校驗數(shù)據(jù)檢驗所有模型，R2最大時對應的參數(shù)組合即為最佳參數(shù)，由此確定6個子序列對應的最佳SVR模型。

3.2.3 BRR-SVR模型 在6組子序列已經(jīng)選取了對應的最佳SVR和BRR的基礎(chǔ)上，比較每組SVR和BRR在對應的校驗數(shù)據(jù)上的表現(xiàn)，發(fā)現(xiàn)對于高頻序列D1、D4、D5和低頻序列A5，BRR模型的精度往往高于SVR，而在序列D2和D3上，SVR的精度更高。因此，BRR-SVR的原理是，比較各序列對應的BRR 子模型和SVR子模型在校驗數(shù)據(jù)上的R2值，以R2值更大的模型作為最適合該序列的模型，得到6組子序列的最佳模型，從而構(gòu)建出BRR-SVR組合模型。由于R2表示預測值與實測值的逼近程度，R2越大代表預測值越接近真實值，因此選擇R2值更大的模型能夠確保每個子序列的預測值更逼近真實值，以此達到優(yōu)化的效果。

圖1 BRR-SVR模型流程圖

4 實例分析

以北京站、南京站和7 個太湖流域的雨量站1951—2015年間660 個月的降水資料為例，研究BRR模型、SVR模型和BRR-SVR模型對月降水量的預測性能。9個雨量站的降水數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)1951年1月1日至2015年12月31日的無空缺日降水數(shù)據(jù)，日降水數(shù)據(jù)經(jīng)過累加得到月降水數(shù)據(jù)。對月降水數(shù)據(jù)做中心化處理和五級小波分解重構(gòu)，得到6 組子序列。對所有站點的數(shù)據(jù)，令m=11，對各序列做相同的相空間重構(gòu)，重構(gòu)后每組樣本容量為649。為了保證能較好的訓練，防止模型對訓練數(shù)據(jù)過擬合，將樣本按時間順序分為訓練數(shù)據(jù)∶校驗數(shù)據(jù)∶測試數(shù)據(jù)=529∶60∶60。以上過程均通過PC機上的Python3.7編程實現(xiàn)。

4.1 模型參數(shù)校驗及優(yōu)化模型組合根據(jù)3.2節(jié)中的方法對SVR模型和BRR模型進行校驗，并選擇每組序列的最優(yōu)模型。不同站點不同序列所對應的SVR模型超參數(shù)不同，受篇幅限制，本文只列出北京站和南京站的SVR 超參數(shù)取值（見表1）。不同站點的6 組子序列對應的BRR模型超參數(shù)幾乎一致，故統(tǒng)一設(shè)定為α1=10-3，α2=10-8，λ1=10-8，λ2=10-3。最后利用校驗數(shù)據(jù)從最優(yōu)SVR和BRR中篩選出該序列的最優(yōu)模型，各站點的模型組合見表2。

表1 北京站和南京站不同序列的SVR超參數(shù)取值

表2 各站點不同序列的模型組合（1-BRR模型，2-SVR模型）

4.2 模型預測效果評價利用校驗后的模型對訓練期的數(shù)據(jù)進行擬合，并對校驗期和測試期的數(shù)據(jù)進行預測，3類模型對不同時期降水量的擬合和預測結(jié)果見表3、表4和表5。

由表3 可知，對于訓練數(shù)據(jù)，BRR-SVR的擬合精度明顯高于BRR和SVR，綜合3 項指標來看，SVR對訓練數(shù)據(jù)的擬合能力優(yōu)于BRR。由表5可知，對于測試數(shù)據(jù)，同樣是BRR-SVR的預測精度更高，但BRR的泛化能力總體上高于SVR。對比訓練數(shù)據(jù)和測試數(shù)據(jù)上的R2值可知，SVR的擬合精度基本在0.92左右，而預測精度大多在0.82～0.9之間；BRR的擬合精度和預測精度均在0.9左右。對比SVR和BRR 在訓練數(shù)據(jù)和預測數(shù)據(jù)上的表現(xiàn)可知，BRR的擬合能力稍遜于SVR，但面對新數(shù)據(jù)，BRR的泛化能力反而更強。

對比表4和表5可知，模型對不同時期的預測精度差異不大。從R2值上來看，BRR-SVR 大多高于0.89，BRR略低于BRR-SVR，對絕大多數(shù)站點，SVR模型的精度低于0.89，波動范圍較大。從平均絕對百分比誤差上看，大體上BRR-SVR＜BRR＜SVR，除了北京站大于2，其他站點基本處于0.2～0.6之間。從平均絕對誤差上看，總體滿足BRR-SVR＜BRR＜SVR，北京站的MAE值在8.32～12.82 之間，而其他站點的值處于18～30之間。由以上分析可知，3類模型對北京站的預測效果與長江中下游區(qū)域的水文站差異較大，北京市為國家重要的政治和經(jīng)濟中心，南京市為長江三角洲區(qū)域一體化發(fā)展等多個國家戰(zhàn)略交匯的中心城市，因此選擇具有代表性的南京站和北京站做進一步闡述。

4.2.1 北京站預測效果分析 各模型對北京站2011—2015年間的月降水量預測效果見圖2，圖中實線為觀測值，虛線為預測值。由圖2 可知，3 類模型都能較精確的預測下個月的降水量，從直觀上看，BRR模型的逼近能力不如另外兩類模型。

由表5可知，對于北京站，BRR-SVR與SVR的R2值均為0.966，高于BRR的0.932，從平均絕對百分比誤差來看，BRR-SVR的精度最高，其次是BRR，SVR的預測精度最低。從平均絕對誤差來看，SVR的精度略高于BRR-SVR，兩者均明顯高于BRR。由此可知，對于北京站，SVR的預測精度高于BRR，而BRR-SVR模型的預測精度與SVR基本相同，并未體現(xiàn)明顯的優(yōu)越性。

4.2.2 南京站預測效果分析 3類模型對南京站2011—2015年間月降水量的預測效果見圖3，由圖3可見，3類模型的預測值都能有效逼近觀測值，基本捕捉到了南京站月降水量的變化規(guī)律。除了2011年3月、2011年8月、2013年1月、2013年12月和2015年3月這幾個點的預測值有直觀上的差異外，3 者在其他點的預測值沒有明顯差別。此外，3 類模型對汛期降水量的預測效果要比非汛期更好，這對于城市的防洪建設(shè)極為重要。

表3 訓練期各站點BRR-SVR、BRR和SVR模型擬合精度計算結(jié)果

表4 校驗期各站點BRR-SVR、BRR和SVR模型預測精度計算結(jié)果

表5 測試期各站點BRR-SVR、BRR和SVR模型預測精度計算結(jié)果

由表5 可知，對于南京站，BRR-SVR模型的R2更接近1，平均絕對百分比誤差和平均絕對誤差都比SVR和BRR 小，可見BRR-SVR模型在南京站上更具優(yōu)勢。相比于SVR，BRR模型的R2更大，平均絕對百分比誤差更低，僅在平均絕對誤差上表現(xiàn)次于SVR，因此對于南京站，BRR的預測性能要優(yōu)于SVR。

4.2.3 北京站與南京站預測效果的對比 由4.1節(jié)中的分析可知，僅考慮R2值和平均絕對誤差的情況下，3類模型對北京站的預測精度高于長江中下游的水文站。以BRR-SVR模型為例，將北京站和南京站的降水量置于同一縱坐標水平下（圖4），由圖4可見，北京站每年11月至次年2月，月降水量曲線十分逼近橫坐標軸，這段期間降水量處于極低的水平，對應的MAPE表達式（15）中的分母幾乎為0，導致北京站的平均絕對百分比誤差偏大。因此，不考慮MAPE指標的情況下，可以認為模型對北京站的預測精度高于長江中下游。

圖2 BRR-SVR、BRR和SVR對北京站月降水量預測結(jié)果

圖3 BRR-SVR、BRR和SVR對南京站月降水量預測結(jié)果

4.3 討論本節(jié)以各模型對南京站的預測相對誤差分布（圖5）為例，探討B(tài)RR-SVR模型的精度高于BRR和SVR的原因，與yi分別表示各樣本的預測值與觀測值。由圖5可見，BRR-SVR的誤差整體處于較低水平，其在某一點的誤差值低于SVR和BRR或更接近兩者中的較低值。如2011年9月，BRR的誤差接近3，遠高于SVR，而BRR-SVR模型的誤差更接近SVR；2014年12月，SVR的誤差高于4，BRR的誤差也接近1，而BRR-SVR模型的相對誤差幾乎為0，2011年4月、2013年1月和12月都有類似情況，這說明BRR-SVR模型能夠?qū)深惸Ｐ汀叭￠L補短”，使模型的預測誤差整體上維持較低的水平。

圖4 BRR-SVR模型對北京站和南京站預測效果

圖5 BRR-SVR、BRR和SVR對南京站預測值相對誤差

5 結(jié)論與展望

本文構(gòu)建了BRR、SVR和BRR-SVR優(yōu)化模型，以北京站、南京站和太湖流域7個雨量站的數(shù)據(jù)來驗證模型的性能，結(jié)論如下：（1）對文中所應用的各站點，BRR-SVR、BRR和SVR3類模型都存在對汛期降水量的預測效果比非汛期更好的特點；3類模型對北京站的預測效果整體上優(yōu)于對長江中下游雨量站的預測效果。（2）對于北京站，訓練階段BRR-SVR為最優(yōu)模型，測試階段BRR-SVR的精度與SVR基本相同，未體現(xiàn)明顯優(yōu)越性；SVR的精度在訓練階段和測試階段均高于BRR。（3）對于長江中下游雨量站，BRR-SVR在訓練階段和測試階段均為最優(yōu)模型；BRR對訓練數(shù)據(jù)的擬合能力稍遜于SVR，但對新數(shù)據(jù)的泛化能力強于SVR。

整體上看，對于長江中下游雨量站，BRR-SVR能夠從BRR和SVR中選擇最優(yōu)模型，提升預測精度；對于北京站，由于SVR和BRR的精度已經(jīng)很高，故預測北京站降水使用單一SVR模型即可。本文依據(jù)經(jīng)驗選取Daubechies 小波為小波基，以后研究中將嘗試用其他類型小波。另外，對北京站的研究，若要得到更具代表性的結(jié)論，需要在以后的工作中獲取更多站點的降水數(shù)據(jù)來驗證。