王 瑞，萬定生

(河海大學(xué)計算機(jī)與信息學(xué)院，南京 211100)

中小河流的水文數(shù)據(jù)有時序特點和復(fù)雜性特點，水文時間序列預(yù)測前提是需要提取符合時序數(shù)據(jù)的特征，要能夠?qū)W習(xí)中小河流復(fù)雜的數(shù)據(jù)特征。為了確保水文時間序列預(yù)測結(jié)果更為準(zhǔn)確，水文時間序列特征提取成為預(yù)測工作的關(guān)鍵步驟。支持向量機(jī)[1](support vector machine,SVM)在結(jié)構(gòu)風(fēng)險最小化理論保證了其具備良好的推廣能力。近年來，支持向量機(jī)模型與其他方法相結(jié)合逐漸成為研究熱點[2]。

時間序列的特征提取方法有4種，分別是基于統(tǒng)計特征的分類特征提取、基于構(gòu)建模型的分類特征提取、基于變換的分類特征提取以及基于分形理論的特征提取?；诮y(tǒng)計量的特征提取方法是最直接的特征提取方法，提取時間序列數(shù)據(jù)在統(tǒng)計學(xué)上的特征構(gòu)成特征向量，統(tǒng)計特征有兩種：時間域與頻率域。時間域特征包括均值、峰值、方差、均方根和峰值因子等[3]；基于構(gòu)建模型的分類特征提取方法是將提取時間序列特征等價于提取模型因子，首先分析數(shù)據(jù)特點，并根據(jù)已完成分析的不同的數(shù)據(jù)特點構(gòu)建不同的模型[4]，如針對相對穩(wěn)定的時間序列數(shù)據(jù)，可以通過自回歸滑動平均模型提取特征，針對相對不穩(wěn)定的時間序列可以對數(shù)據(jù)進(jìn)行神經(jīng)網(wǎng)絡(luò)建模，利用神經(jīng)網(wǎng)絡(luò)模型演化求解代表時間序列數(shù)據(jù)變化的特征；基于變換的分類特征提取方法，將時間序列特征數(shù)據(jù)在不同域中映射變換，時域和頻域的變換是比較常見的域變換，使得顯著特征在特定維度凸顯，典型的域變換有傅里葉變換和小波變換[5]；基于分形理論的分類特征提取方法是針對符合遞歸生成原則和自相似結(jié)構(gòu)特點的現(xiàn)象，如自然界中廣泛分布的分形現(xiàn)象如水位線，河流的流向等進(jìn)行理論分析，利用分形理論對時間序列數(shù)據(jù)進(jìn)行特征提取時，將時間序列數(shù)據(jù)轉(zhuǎn)換成信號，分析信號在特定尺度下的可分形特征[6]。

然而，在實際展開水文時間序列預(yù)測研究中，基于構(gòu)建模型的分類特征提取方法應(yīng)用居多，采用經(jīng)典神經(jīng)網(wǎng)絡(luò)建模提取特征[7]；或者直接選用距離預(yù)報前時刻4 h或距離預(yù)報前時刻6 h作為基本特征[8]，缺少對水文時間序列特征提取的更深層更細(xì)致研究。

水文時間序列特征選擇是典型的迭代優(yōu)化問題，神經(jīng)網(wǎng)絡(luò)的參數(shù)訓(xùn)練也是迭代問題，同時神經(jīng)網(wǎng)絡(luò)參數(shù)訓(xùn)練是一個比較耗時的過程，如果使用神經(jīng)網(wǎng)絡(luò)計算特征選擇過程的誤差，求解這組特征組合等價于雙層神經(jīng)網(wǎng)絡(luò)優(yōu)化問題，時間復(fù)雜度較高。

因此，基于構(gòu)建模型的分類特征提取方法，現(xiàn)提出一種基于支持向量回歸和高斯過程回歸的水文時間序列特征提取方法，首先列出水文時間序列候選特征，并隨機(jī)組合特征，不同特征組合下，分別對訓(xùn)練數(shù)據(jù)進(jìn)行支持向量回歸建模與高斯過程回歸建模，利用遺傳算法演化求解一組目標(biāo)特征組合，使得支持向量回歸和高斯過程回歸輸出誤差同時最小。以屯溪流域水文時間序列數(shù)據(jù)進(jìn)行實驗分析。

1 相關(guān)研究

1.1 支持向量回歸

給定i個水文時間序列樣本數(shù)據(jù)集{xi,yi}，xi為訓(xùn)練樣本輸入，yi為訓(xùn)練樣本輸出，φ(xi)為xi通過非線性映射函數(shù)將低維特征映射到高維特征空間H的特征向量，f(xi)對應(yīng)的最優(yōu)超平面為

f(xi)=wTφ(xi)+b

(1)

式(1)中：w為法向量；b為偏置量。

當(dāng)且僅當(dāng)f(xi)與yi完全相同時，記為損失為零，SVM認(rèn)為此時學(xué)習(xí)到的超平面最優(yōu)。與此不同，支持向量回歸(support vector regression，SVR)[9]引入松弛因子ε，并定義當(dāng)f(xi)與yi之間的差別絕對值大于ε時才計算損失，當(dāng)f(xi)與yi之間的差別絕對值小于ε時認(rèn)為學(xué)習(xí)到的超平面為最優(yōu)。SVR問題形式化為

(2)

約束條件為

(3)

最終得到函數(shù)表達(dá)式為

(4)

1.2 高斯過程回歸

高斯過程回歸(gaussian process regression,GPR)[10]基于高斯過程先驗對數(shù)據(jù)進(jìn)行回歸分析。高斯過程回歸模型包括回歸殘差和高斯過程先驗兩個內(nèi)容。假設(shè)一個存在殘差的回歸模型。其計算公式為

Y=f(X)+ξ

(5)

(6)

(7)

式中：In為n維單位矩陣；K(x,x)=(xij)為對稱正定協(xié)方差矩陣，xij通過核函數(shù)度量xi與xj之間的相關(guān)性；K(x*,x)=K(x,x*)T為訓(xùn)練集x和測試集x*之間的協(xié)方差矩陣；K(x*,x*)為測試集本身的協(xié)方差矩陣。平方指數(shù)核函數(shù)具有無窮可微的特點，使高斯過程回歸平滑，因此選用平方指數(shù)核[11]。其計算公式為

(8)

式(8)中：p1、p2為高斯函數(shù)中的可調(diào)參數(shù)。預(yù)測值y的后驗分布為

(9)

(10)

1.3 遺傳算法

遺傳算法(genetic algorithm,GA)[12]是一種自適應(yīng)全局搜索最優(yōu)解的迭代算法，基于生物進(jìn)化論自然選擇和遺傳學(xué)機(jī)理的生物進(jìn)化理論，模擬自然進(jìn)化過程，即保留遺傳過程中優(yōu)勢特征，淘汰劣勢特征。

遺傳算法中問題的每一個候選解被編碼成染色體個體，若干個個體相互排列組合構(gòu)成若干個群體，即所有可能解。遺傳算法在迭代開始時，首先隨機(jī)產(chǎn)生一個初始解，根據(jù)設(shè)定的目標(biāo)函數(shù)對單個個體進(jìn)行評估，計算此時的適應(yīng)值，根據(jù)適應(yīng)值對個體進(jìn)行輪盤賭選擇策略、交叉運算、選擇操作等操作產(chǎn)生新一代種群，對產(chǎn)生的新一代種群進(jìn)行適應(yīng)度計算評估，直至達(dá)到目標(biāo)函數(shù)優(yōu)化終止條件。

遺傳算法過程圖解如圖1所示。

圖1 遺傳算法過程圖Fig.1 Genetic algorithm process diagram

2 基于支持向量回歸和高斯過程回歸的水文時間序列特征提取方法

基于支持向量回歸和高斯過程回歸的水文時間序列特征提取方法架構(gòu)總體分為兩部分。

模塊一：從屯溪流域屯溪水文站獲取原始水文數(shù)據(jù)，將水位時間數(shù)據(jù)按照相等小時時間間隔排列，構(gòu)造待插補數(shù)據(jù)輸入序列，構(gòu)建Stacking集成學(xué)習(xí)法，Stacking集成學(xué)習(xí)小波神經(jīng)網(wǎng)絡(luò)和支持向量機(jī)對數(shù)據(jù)分析的結(jié)果，將待插補水位數(shù)據(jù)序列輸入已構(gòu)造好的Stacking學(xué)習(xí)模型中，輸出插值補充好的完整的水文時間序列。

模塊二：將插補好的水文時間序列作為特征提取的輸入，首先對數(shù)據(jù)集逆行歸一化處理，然后分別構(gòu)建支持向量回歸模型和高斯過程回歸模型，利用遺傳算法迭代求解最優(yōu)特征組合。

基于支持向量回歸和高斯過程回歸的水文時間序列特征提取方法架構(gòu)圖如圖2所示。

圖2 特征提取方法架構(gòu)圖Fig.2 Feature extraction method architecture diagram

2.1 水文時間序列特征提取方法原理分析

水文時間序列的連續(xù)性是把某一時間段內(nèi)水文數(shù)據(jù)按整點時間的連續(xù)性進(jìn)行排列，須保證整點時刻均為有效數(shù)據(jù)。中小河流原始水文數(shù)據(jù)缺失情況普遍存在，提出采用Stacking集成學(xué)習(xí)法，融合小波神經(jīng)網(wǎng)絡(luò)(wavelet neural network，WNN)[13-14]模型和支持向量機(jī)(support vector machine，SVM)模型，預(yù)測待插補水文數(shù)據(jù)空缺值。具體實現(xiàn)算法如算法1所示。

算法1

輸入：待插補水文時間序列觀測數(shù)據(jù)集D={(x1,y1),(x2,y2),…,(xm,ym)}；

輸出：預(yù)測值result

1：輸入X(m)

2：數(shù)據(jù)歸一化

/*本實驗采用0-1標(biāo)準(zhǔn)歸一化*/

3：fori=1,2,…,mdo

4：WNN預(yù)測；

5：end for

6：fori=1,2,…,mdo

7：SVM預(yù)測

8：i++；

9：end for

10：Stacking融合WNN模型&SVM模型

11：做出預(yù)測結(jié)果

12：輸出待插補預(yù)測數(shù)據(jù)

在特征提取過程中，將水文時間序列數(shù)據(jù)輸入，列出水文時間序列候選特征，置為0或1兩個狀態(tài)，0代表舍棄當(dāng)前特征，1代表選擇當(dāng)前特征，列出特征組合的所有解，將求解最優(yōu)特征組合問題等價于0-1規(guī)劃問題，并對水文時間序列數(shù)據(jù)分別構(gòu)建支持向量回歸模型和高斯過程回歸模型，利用遺傳算法迭代求解一組最優(yōu)特征組合。具體實現(xiàn)算法如算法2所示。

算法2

輸入：水位時間序列數(shù)據(jù)X(l)；

輸出：最優(yōu)特征組合Fj；

Fj表示距離當(dāng)前預(yù)報時刻的第j個時刻

1：輸入X(l)

2：水位時間序列數(shù)據(jù)進(jìn)行數(shù)據(jù)歸一化

/*本實驗采用0-1標(biāo)準(zhǔn)歸一化*/

3：GA種群初始化

4：輪盤賭選擇、交叉和變異

5：GA演化

6：計算適應(yīng)度fitness

7：SVR建模&GPR建模

8：if fitness=Min(RMSE)

9：else go to 6

10：輸出最優(yōu)特征組合Fj

3 實驗

3.1 實驗準(zhǔn)備

以屯溪流域為中小河流代表流域進(jìn)行研究，屯溪流域地勢西高東低，氣候濕潤。屯溪位于新安江上游，新安江上游多為峽谷，河形彎曲，右岸靠山，左岸河谷平原，經(jīng)新安江上游率水、橫江兩溪匯流以后成為屯溪。屯溪流域地形圖如圖3所示。

中小流域的水文數(shù)據(jù)受水文數(shù)據(jù)觀測儀器精度，地理環(huán)境和氣候環(huán)境因素等影響，中小河流數(shù)據(jù)記錄缺失和數(shù)值缺失現(xiàn)象普遍存在。

采用Stacking集成學(xué)習(xí)法，融合小波神經(jīng)網(wǎng)絡(luò)模型和支持向量機(jī)模型預(yù)測水文數(shù)據(jù)空缺值，預(yù)測水文缺失數(shù)據(jù)，并將完整的水文時間序列結(jié)果作為輸入數(shù)據(jù)。

圖3 屯溪流域地形圖Fig.3 Tunxi basin topographic map

選取屯溪流域屯溪水文站2020年5月1日00：00—6月2日1：00的水位時間序列數(shù)據(jù)，時間間隔1 h。前620條用于訓(xùn)練，提取水文時間序列特征；后154條數(shù)據(jù)用于測試，將提取的水位特征應(yīng)用于水位預(yù)測，測試本方法所提取特征在水文預(yù)測中的表現(xiàn)能力。水文時間序列數(shù)據(jù)的輸入格式如表1所示。

表1 輸入數(shù)據(jù)格式Table 1 Input data format

模型參數(shù)的設(shè)定是模型成功的重要因素。通過暴力破解算法窮舉法，設(shè)初始種群個數(shù)為20，交叉率為0.5，變異率為0.5，種群迭代次數(shù)為100次。

3.2 實驗準(zhǔn)備

(1)均方根誤差。采用均方根誤差(root mean squared error，RMSE)作為特征選擇評價指標(biāo)。即RMSE越小，則本次特征提取越具有代表性。其計算公式為

(11)

式(11)中：n為測試樣本總數(shù)；yi為預(yù)測值；Yi為真實值。

(2)決定系數(shù)。決定系數(shù)R2度量自變量解釋比例，反映回歸方程擬合優(yōu)程度。其中。R2越接近1，說明數(shù)據(jù)擬合效果就越好。其計算公式為

(12)

3.3 實驗結(jié)果與分析

對于選定訓(xùn)練數(shù)據(jù)集進(jìn)行水文時間序列特征提取，將距離預(yù)報時刻前10 h水位特征進(jìn)行提取，被提取特征結(jié)果為[1，2，5，6]，表示將預(yù)報時刻前第一個小時，第二個小時，第五個小時，第六個小時水位數(shù)據(jù)作為數(shù)據(jù)特征，應(yīng)用于水文時間序列預(yù)報。水文時間序列特征提取結(jié)果如表2所示。

表2 水文時間序列特征篩選結(jié)果表Table 2 Table of hydrological time series feature extraction

將預(yù)報時刻前第一個小時、第二個小時、第五個小時、第六個小時的水位數(shù)據(jù)作為特征輸入數(shù)據(jù)，未來1 h的水位Yt+1作為輸出進(jìn)行預(yù)報，即

Yt+1=f(Xt-6,Xt-5,Xt-2,Xt-1,Yt)

(13)

為了使實驗對比更明顯，利用CNN單一神經(jīng)網(wǎng)絡(luò)方法提取的歷史水文時間序列特征,所篩選特征表示為Xt-4,Xt-3,Xt-2,Xt-1，結(jié)合LSTM應(yīng)用于水文時間序列預(yù)測，其預(yù)測結(jié)果如圖4所示。

基于支持向量回歸和高斯過程回歸建模，并用遺傳算法演化出一組水文時間序列特征組合，將篩選出的這組特征與水文時間序列作為輸入數(shù)據(jù)，通過相空間重構(gòu)將低維時間序列轉(zhuǎn)換成高維時間序列，并結(jié)合長短期記憶(long short-term memory，LSTM)神經(jīng)網(wǎng)絡(luò)對其進(jìn)行水文時間序列預(yù)測。其預(yù)報結(jié)果如表3所示。

圖4 CNN特征提取方法預(yù)測結(jié)果Fig.4 Prediction results of CNN feature extraction method

表3 支持向量回歸和高斯過程回歸方法的水位預(yù)測結(jié)果Table 3 Table of hydrological time series method based on SVR_GPR

為了使實驗效果更清晰，選出101條預(yù)測結(jié)果進(jìn)行展示，經(jīng)支持向量回歸和高斯過程回歸特征提取方法應(yīng)用于水位預(yù)測結(jié)果如圖5所示。

將基于CNN提取的水文時間序列特征與基于支持向量回歸和高斯過程回歸提取的特征用于水文時間序列預(yù)測，并將其預(yù)測結(jié)果進(jìn)行對比，同時選擇RMSE和R2對上述兩種提取方法進(jìn)行模型參數(shù)檢驗，得到的參數(shù)統(tǒng)計結(jié)果如表4所示。

圖5 支持向量回歸和高斯過程回歸特征提取方法預(yù)測結(jié)果Fig.5 Predicted results of SVR-GPR method

表4 不同特征提取方法預(yù)測結(jié)果對比Table 4 Table of comparison of prediction results of different feature extraction methods

4 結(jié)論

提出的基于支持向量回歸和高斯過程回歸的水文時間序列特征提取方法，對屯溪流域水文時間序列數(shù)據(jù)進(jìn)行應(yīng)用分析。基于支持向量回歸和高斯過程回歸的水文時間序列特征提取方法，避免神經(jīng)網(wǎng)絡(luò)優(yōu)化遺傳算法不斷尋參的過程，采用支持向量回歸模型和高斯過程回歸模型對水文時間序列進(jìn)行建模，將雙層網(wǎng)絡(luò)優(yōu)化簡化為單層優(yōu)化問題，大大提高了水文時間序列特征篩選的時間效率，在水文時間序列預(yù)測上具有實際的應(yīng)用意義。

實驗結(jié)果說明，基于支持向量回歸和高斯過程回歸的水文時間序列特征提取方法依賴于水文時間序列變化趨勢，不同于直接選用預(yù)報前幾個小時的數(shù)據(jù)作為特征輸入數(shù)據(jù)，通過這種提取水位特征方法，能及時有效地捕捉影響當(dāng)前水位的高度相關(guān)特征，并將提取的特征結(jié)果結(jié)合神經(jīng)網(wǎng)絡(luò)應(yīng)用于水文時間序列預(yù)測，大大提高了水文時間序列預(yù)測的精準(zhǔn)度，同時為水利信息化建設(shè)提供一些實質(zhì)性的方法捕捉水文時間序列特征。

研究將從對水文時間序列進(jìn)行汛期與非汛期不同周期出發(fā)，研究汛期與非汛期的水文時間序列特征提取，找出水文時間序列預(yù)測中的更具普適性的水文時間序列特征。