胡騰飛，施勇，毛勁喬，欒震宇，陳煉鋼，陳黎明，金秋，徐祎凡，戴會超

(1. 南京水利科學(xué)研究院水文水資源研究所，江蘇 南京 210029; 2. 水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029； 3. 河海大學(xué)水利水電學(xué)院，江蘇 南京 210098)

河流和湖泊水位建模常用技術(shù)包括具有物理機(jī)制的水動力模型[1-2]和數(shù)據(jù)驅(qū)動技術(shù)[3-5].與機(jī)理模型相比，數(shù)據(jù)驅(qū)動技術(shù)能從有限的訓(xùn)練樣本中學(xué)習(xí)輸入輸出間的內(nèi)在聯(lián)系，因此模型輸入通常只包含所關(guān)注變量及其關(guān)聯(lián)要素的時間序列[6].數(shù)據(jù)驅(qū)動技術(shù)不需要對物理過程的全面認(rèn)識，且能更快地提供預(yù)測結(jié)果.河湖水位預(yù)測常用的數(shù)據(jù)驅(qū)動技術(shù)有支持向量回歸(support vector regression，SVR)和多種人工神經(jīng)網(wǎng)絡(luò).基于統(tǒng)計(jì)學(xué)習(xí)理論的SVR在水位預(yù)測領(lǐng)域被認(rèn)為是比人工神經(jīng)網(wǎng)絡(luò)更好的選擇[7].

輸入變量選擇(input variable selection，IVS)目標(biāo)是確定SVR模型的最優(yōu)結(jié)構(gòu)形式.常用的IVS方法可分為獨(dú)立于模型的Filter技術(shù)和基于模型的Wrapper技術(shù)以及Embedded技術(shù)[8].Filter技術(shù)通過統(tǒng)計(jì)方法分析可用數(shù)據(jù)以衡量變量或變量組合的重要性，選擇結(jié)果對不同數(shù)據(jù)驅(qū)動方法具有通用性；然而，此類技術(shù)忽略了底層模型預(yù)測能力對選擇結(jié)果的響應(yīng).Wrapper技術(shù)通過搜索變量的所有(或大部分)組合來確定對應(yīng)數(shù)據(jù)驅(qū)動方法最大預(yù)測能力的輸入變量集合[9]；然而，由于底層模型沒有針對候選變量進(jìn)行充分率定，其真實(shí)預(yù)測能力可能被掩蓋[10].Embedded技術(shù)可使不相關(guān)或冗余變量的影響在訓(xùn)練中逐漸消失；然而，此類技術(shù)僅面向特定的數(shù)據(jù)驅(qū)動方法.

為更好地解決河湖水位SVR建模中的輸入變量選擇問題，文中提出基于進(jìn)化算法的SVR模型結(jié)構(gòu)-參數(shù)同步優(yōu)化方法.為驗(yàn)證所提方法的有效性，該方法被應(yīng)用于長江流域2個不同區(qū)域的實(shí)際案例中.

1 研究方法與數(shù)據(jù)

1.1 研究區(qū)域

文中選取荊江和洞庭湖區(qū)域作為SVR模型結(jié)構(gòu)-參數(shù)同步優(yōu)化方法的驗(yàn)證區(qū)域，具體涉及洞庭湖城陵磯(七里山)站和荊江陳二口站水位模擬預(yù)測.該區(qū)域河湖水流既相互補(bǔ)給又相互制衡，水流運(yùn)動復(fù)雜，適合檢驗(yàn)水位建模方法的適用性及合理性.洞庭湖位于中國長江中游地區(qū)荊江南岸，是中國面積第二大淡水湖.洞庭湖在調(diào)蓄長江洪水方面發(fā)揮著重要作用[11]，也因此呈現(xiàn)“漲水是湖，退水為洲”的獨(dú)特景觀格局.長江洪水通過荊江南岸的松滋、太平和藕池三口分流入湖(曾有第四口調(diào)弦口，于1959年封堵)；洞庭湖還承接其4條主要支流—湘、資、沅、澧的來水，入湖水流經(jīng)湖泊調(diào)蓄后在城陵磯匯入長江.

1.2 研究方法

基于SVR構(gòu)建復(fù)雜河湖水系的水位模型時，為應(yīng)對考慮關(guān)聯(lián)要素時滯后輸入變量搜索空間的高維特性，提出了基于進(jìn)化算法的SVR模型結(jié)構(gòu)-參數(shù)同步優(yōu)化方法.首先簡要介紹SVR及核函數(shù)的基本理論，再重點(diǎn)介紹文中提出的SVR模型結(jié)構(gòu)-參數(shù)同步優(yōu)化方法.

1.2.1 支持向量回歸SVR

支持向量機(jī)的基礎(chǔ)為VAPNIK[12]創(chuàng)建的統(tǒng)計(jì)學(xué)習(xí)理論.近年來，支持向量機(jī)憑借其優(yōu)良特性和強(qiáng)大的泛化能力受到研究者青睞[13-14].

1) 支持向量回歸函數(shù).SVR通過數(shù)學(xué)方法歸納樣本輸入輸出數(shù)據(jù)中隱含的依賴關(guān)系，這一依賴關(guān)系獲得后便可對未知輸出進(jìn)行預(yù)測.對于數(shù)據(jù)集DS={(x1,y1), (x2,y2),..., (xl,yl)}，xi∈Rn為輸入變量，yi∈R1是相應(yīng)的輸出變量，數(shù)據(jù)集DS的輸入輸出映射可通過以下回歸函數(shù)估計(jì)

f(x)=ωTφ(x)+b,

(1)

式中：ω為權(quán)重向量；φ為將x映射到高維特征空間的非線性映射；b為截距.

基于ε不敏感損失函數(shù)(ε>0為誤差閾值)，可通過優(yōu)化以下函數(shù)得到回歸函數(shù)[12]，即

(2)

由于權(quán)重向量ω的高維特性，通常求解式的對偶問題，引入拉格朗日乘子得到

(3)

借助核函數(shù)，高維特征空間內(nèi)積φ(xi)Tφ(xj)可直接在輸入空間通過訓(xùn)練樣本計(jì)算得出，不需要進(jìn)行輸入向量映射到特征空間的計(jì)算，核函數(shù)的引入大大減小了計(jì)算量.

通過求解式(3)，回歸函數(shù)可最終表示為

(4)

2) 核函數(shù).在水文水資源領(lǐng)域，絕大多數(shù)應(yīng)用都采用徑向積函數(shù)并獲得了理想的模型精度[15].因此文中河湖水位建模采用徑向積函數(shù)，其表達(dá)式為

K(xi,xj)=exp(-γ‖xi-xj‖2)，

(5)

式中：γ> 0為控制核寬度的參數(shù).

1.2.2 SVR模型結(jié)構(gòu)-參數(shù)同步優(yōu)化

SVR模型的輸入輸出結(jié)構(gòu)一般可表示為

(6)

式中：fSVR為基于SVR構(gòu)建的輸入輸出映射；變量xi，(i=1, 2, …,N)為目標(biāo)變量y的關(guān)聯(lián)要素;t為時間;mi和ni，(i=0, 1, …,N)為對應(yīng)變量的最小和最大時滯.

若模型輸入不考慮變量y的前期值，則將式(6)簡化為

(7)

優(yōu)化SVR模型的輸入輸出結(jié)構(gòu)，即是要確定對目標(biāo)變量預(yù)測有重要作用的關(guān)聯(lián)要素及其時滯.在目標(biāo)變量關(guān)聯(lián)要素眾多、各要素時滯不同且變化范圍較大時，潛在可行的模型輸入變量及時滯組合極多，從中選擇最優(yōu)組合比較困難.鑒于基于種群的進(jìn)化算法可很好地應(yīng)對高維搜索空間，文中河湖水位SVR模型結(jié)構(gòu)優(yōu)化將基于進(jìn)化算法開展.另一方面，為保證輸入變量及時滯組合對應(yīng)的預(yù)測能力得以準(zhǔn)確體現(xiàn)，在結(jié)構(gòu)優(yōu)化的同時對模型參數(shù)也進(jìn)行同步優(yōu)化.為避免模型出現(xiàn)過擬合問題，SVR模型訓(xùn)練采用n折交叉驗(yàn)證方式.基于進(jìn)化算法的SVR模型結(jié)構(gòu)-參數(shù)同步優(yōu)化方法中，種群染色體適應(yīng)度評估流程見圖1.

圖1 SVR模型結(jié)構(gòu)-參數(shù)同步優(yōu)化中染色體的評估流程

由圖1可知，SVR完全嵌入了結(jié)構(gòu)-參數(shù)同步優(yōu)化流程，因此該方法屬于基于模型的IVS方法.對于常規(guī)的僅支持浮點(diǎn)數(shù)編碼的進(jìn)化算法，SVR模型結(jié)構(gòu)-參數(shù)同步優(yōu)化中染色體評估按以下步驟執(zhí)行：① 將染色體分為2個部分，一部分提供變量時滯信息，另一部分提供模型參數(shù)；② 將指示每個變量時滯的2個浮點(diǎn)數(shù)向上取整以獲取該變量的最大最小時滯；根據(jù)所有變量的最大最小時滯信息組織模型輸入；③ 使用n折交叉驗(yàn)證訓(xùn)練SVR模型；④ 將交叉驗(yàn)證均方誤差用作所評估染色體的適應(yīng)度.

需要特別說明的是，進(jìn)化算法中變量時滯的搜索范圍需設(shè)為略大于其合理范圍.在步驟②的取整操作后，若變量時滯超出其合理范圍則在組織模型輸入時徹底舍棄該變量.這種設(shè)定可靈活放棄對目標(biāo)變量預(yù)測作用不大的關(guān)聯(lián)變量.

文中SVR模型結(jié)構(gòu)-參數(shù)同步優(yōu)化方法具有以下特點(diǎn)：① 在SVR模型結(jié)構(gòu)確定過程中同步優(yōu)化模型參數(shù)，可準(zhǔn)確反映候選變量及時滯的真實(shí)預(yù)測能力，有助于發(fā)掘SVR的應(yīng)用潛力；② 模型訓(xùn)練采用了n折交叉驗(yàn)證，有利于避免過擬合問題；③ 所提方法中進(jìn)化算法的引入，可幫助應(yīng)對復(fù)雜水系考慮關(guān)聯(lián)要素時滯后輸入變量搜索空間的高維特性.

1.3 研究數(shù)據(jù)

為構(gòu)建洞庭湖城陵磯(七里山)站水位SVR模型，收集了2014—2017年城陵磯站逐日平均水位、洞庭湖四水(湘江湘潭站、資水桃江站、沅水桃源站和澧水石門站)以及長江宜昌站逐日平均流量數(shù)據(jù)；為構(gòu)建荊江陳二口站水位SVR模型，收集了2014—2017年陳二口站、上游枝城站、下游馬家店站的逐日平均水位，以及同期松滋河新江口站和沙道觀站逐日平均流量數(shù)據(jù).

2 結(jié)果與討論

基于文中提出的結(jié)構(gòu)-參數(shù)同步優(yōu)化方法，構(gòu)建了洞庭湖城陵磯水位SVR模型和荊江陳二口水位SVR模型.文中介紹2個模型的構(gòu)建過程、模型結(jié)構(gòu)參數(shù)及模型精度，從而驗(yàn)證論文提出方法的有效性.

2.1 洞庭湖城陵磯水位SVR模型

2.1.1 模型輸入輸出結(jié)構(gòu)

洞庭湖城陵磯水位SVR模型的輸入為考慮時滯的長江和四水流量，以及城陵磯前期水位.由于城陵磯水位同時受長江和四水來水影響，且兩者時滯年內(nèi)和年際波動較大，因此在模型輸入中考慮自身站位的前期水位有望大幅提升模型預(yù)測能力.

通過SVR構(gòu)建的城陵磯逐日水位與其影響因素間的映射可表示為

(8)

式中：Lclj為城陵磯逐日平均水位；Di(i=1,…,5)分別為湘江湘潭、資水桃江、沅水桃源、澧水石門以及長江宜昌逐日平均流量.各變量時滯mi和ni(i=0, 1, …,5)單位為d，其最小容許值均為1 d.

2.1.2 模型結(jié)構(gòu)-參數(shù)同步優(yōu)化

文中使用2016—2017年相關(guān)數(shù)據(jù)優(yōu)化城陵磯水位SVR模型的結(jié)構(gòu)和參數(shù)并最終訓(xùn)練模型，而2014—2015年數(shù)據(jù)則用于模型驗(yàn)證.其中，訓(xùn)練數(shù)據(jù)中城陵磯水位最高值稍大于驗(yàn)證數(shù)據(jù)，從而保證模型的泛化能力.使用遺傳算法進(jìn)行城陵磯水位SVR模型的結(jié)構(gòu)和參數(shù)優(yōu)化.遺傳算法參數(shù)設(shè)置：浮點(diǎn)數(shù)編碼；最大進(jìn)化代數(shù)300，種群規(guī)模300；參數(shù)ε,C,γ搜索下邊界1.0×10-6、上邊界1.0×106；變量時滯搜索下邊界-0.5、上邊界11，若變量時滯取整得到0，則在組織模型輸入時舍棄對應(yīng)變量；所有模型輸入都線性歸一化至[0, 1]區(qū)間，確保各個輸入在模型訓(xùn)練時具有一致權(quán)重；模型結(jié)構(gòu)-參數(shù)同步優(yōu)化統(tǒng)一采用5折交叉驗(yàn)證的訓(xùn)練方式.

優(yōu)化程序執(zhí)行完畢后，得出城陵磯水位SVR模型參數(shù)值ε=0.006 8，C=33 026.465 2，γ=0.002 1，選中的模型輸入變量時滯見圖2(灰色色塊).對選中的變量時滯進(jìn)行了敏感性分析：① 基于模型訓(xùn)練數(shù)據(jù)，將每個時滯對應(yīng)的輸入數(shù)據(jù)改變±10%而其余輸入保持不變；② 計(jì)算2種輸入對應(yīng)水位預(yù)測結(jié)果的差異的絕對值；③ 取差異絕對值的中位數(shù)作為該時滯對水位變化的影響強(qiáng)度；④ 將所有時滯的影響強(qiáng)度一起排序，圖2中色塊越黑意味著該時滯的影響強(qiáng)度越大.

圖2 城陵磯水位SVR模型輸入變量時滯

由圖2可知，對城陵磯水位預(yù)測影響最大的模型輸入為同站點(diǎn)前期水位，且前1天至前5天水位信息的影響強(qiáng)度逐步弱化.由于河流流量存在自相關(guān)性，圖2中各條河流流量時滯近似還原了水流運(yùn)動時間，而敏感性分析結(jié)果恰當(dāng)?shù)胤从沉瞬煌恿鲗Τ橇甏壦蛔兓绊懙牟町愋?長江來水對位處洞庭湖與荊江交匯處的城陵磯水位有著最為顯著的影響，且由于長江流量年內(nèi)波動大，其時滯變化范圍在所有河流中也最大.洞庭湖四水中，湘江影響最強(qiáng)，澧水影響最弱，資水、沅水介于二者之間.

2.1.3 模型精度評估

文中采用多重指標(biāo)對構(gòu)建的城陵磯水位SVR模型進(jìn)行精度評估，指標(biāo)包括均方根誤差(RMSE)、決定系數(shù)(R2)和平均絕對誤差(MAE).

訓(xùn)練期和驗(yàn)證期的城陵磯水位SVR模型精度評估結(jié)果見表1，城陵磯水位預(yù)測值與實(shí)測值對比見圖3.

由表1可知，2個時段模型RMSE均小于0.060 m，R2均大于0.999 0，MAE均小于0.040 m，取得了理想的模擬預(yù)測精度.模型在訓(xùn)練期和驗(yàn)證期的精度差異很小，未出現(xiàn)過擬合問題，證明模型結(jié)構(gòu)-參數(shù)同步優(yōu)化采用的5折交叉驗(yàn)證發(fā)揮了預(yù)期作用.從圖3可以看出，2個時段的水位預(yù)測值與實(shí)測值吻合較好，特別是在驗(yàn)證期高水位時，模型精度沒有下降，這一結(jié)果可歸因于模型在訓(xùn)練期有機(jī)會學(xué)習(xí)更嚴(yán)重的洪水情景.訓(xùn)練期和驗(yàn)證期模型的最大高估誤差分別為0.17 m和0.29 m，而最大低估誤差分別為-0.32 m和-0.29 m.在訓(xùn)練期和驗(yàn)證期，分別有96.3%和94.2%的水位預(yù)測誤差絕對值不超過0.10 m.

表1 城陵磯水位SVR模型精度評估結(jié)果

圖3 城陵磯水位預(yù)測值與實(shí)測值對比圖

2.2 荊江陳二口水位SVR模型

陳二口水位連續(xù)監(jiān)測開始時間較晚，當(dāng)前有向歷史年份外延陳二口水位數(shù)據(jù)的實(shí)際需求.文中擬基于陳二口上、下游臨近站點(diǎn)的水位和流量信息構(gòu)建陳二口水位SVR模型，因此輸入變量不考慮自身站位的前期水位.

2.2.1 模型輸入輸出結(jié)構(gòu)

荊江陳二口水位SVR模型的輸入為考慮時滯的枝城站和馬家店水位，以及松滋河新江口和沙道觀流量.

通過SVR構(gòu)建的陳二口逐日水位與其影響因素間的映射可表示為

(9)

式中：Lcek為陳二口水位；L1和L2分別為枝城站水位和馬家店水位；D1和D2分別為新江口流量和沙道觀流量.各變量時滯mi和ni，i= 1, 2, …, 4,單位為d，其最小容許值均為0.

2.2.2 模型結(jié)構(gòu)-參數(shù)同步優(yōu)化

為實(shí)現(xiàn)陳二口水位序列向歷史年份外延，文中使用2014—2016年相關(guān)數(shù)據(jù)優(yōu)化陳二口水位SVR模型的結(jié)構(gòu)和參數(shù)并最終訓(xùn)練模型，2017年數(shù)據(jù)則用于模型驗(yàn)證.陳二口水位一般低于枝城站水位、高于馬家店水位，且更接近后者.然而，枯水期陳二口水位與馬家店水位差異顯著，進(jìn)入汛期后兩者差異變小，因此簡單地考慮站點(diǎn)間距進(jìn)行陳二口水位插值難以獲得準(zhǔn)確結(jié)果.

陳二口水位SVR模型結(jié)構(gòu)和參數(shù)優(yōu)化中，遺傳算法設(shè)置與2.1.2節(jié)類似.優(yōu)化得到的陳二口水位SVR模型參數(shù)值ε= 0.095 8，C=274 302.672 1，γ=0.000 1，選中的模型輸入變量時滯見圖4.對選中的變量時滯進(jìn)行了敏感性分析，方法與2.1.2節(jié)類似.

圖4 陳二口水位SVR模型輸入變量時滯

由圖4可知，對陳二口水位預(yù)測作用最大的變量為枝城站和馬家店水位.對于td陳二口水位預(yù)測，上述站點(diǎn)td水位的作用明顯大于t-1 d水位，這與3個站點(diǎn)距離較近有關(guān)(枝城站至馬家店距離不超過35 km).此外，新江口流量(t和t-1 d)和沙道觀流量(td)也被選入模型輸入，說明松滋河流量對陳二口水位的預(yù)測發(fā)揮了積極的輔助作用.

2.2.3 模型精度評估

陳二口水位SVR模型的精度評估仍然基于RMSE,R2和MAE，評估結(jié)果見表2，陳二口水位預(yù)測值與實(shí)測值對比見圖5.表2中2個時段模型RMSE均小于0.070 m，R2均大于0.998 0，MAE均小于0.050 m，說明優(yōu)化后模型取得了理想的模擬預(yù)測精度.模型訓(xùn)練期和驗(yàn)證期的精度未出現(xiàn)過大差異，不存在過擬合問題.由圖5可知，2個時段的水位預(yù)測值與實(shí)測值均可較好吻合，偏離理想線的情況極少.訓(xùn)練期和驗(yàn)證期模型的最大高估誤差分別為0.15 m和0.20 m，最大低估誤差分別為-0.15 m(但存在一個異常點(diǎn)：-1.03 m)和-0.18 m，最大低估幅度與最大高估幅度非常接近，說明模型沒有出現(xiàn)系統(tǒng)性偏差.

表2 陳二口水位SVR模型精度評估結(jié)果

圖5 陳二口水位預(yù)測值與實(shí)測值對比

3 結(jié) 論

1) 基于結(jié)構(gòu)-參數(shù)同步優(yōu)化的水位預(yù)測模型可準(zhǔn)確反映不同影響因素對水位預(yù)測的作用大小，城陵磯水位預(yù)測最主要的外部因素為長江來水和湘江來水，陳二口水位預(yù)測則為枝城站和馬家店水位.

2) 結(jié)構(gòu)-參數(shù)同步優(yōu)化可充分發(fā)掘SVR潛力，城陵磯水位模型(RMSE<0.060 m，R2>0.999 0，MAE<0.040 m)和陳二口水位模型(RMSE<0.070 m，R2>0.998 0，MAE<0.050 m)均取得了理想精度.

3) 提出方法不僅能實(shí)現(xiàn)較高的模型精度，其采用的n折交叉驗(yàn)證方式還可有效避免模型過擬合問題，構(gòu)建的水位預(yù)測模型在訓(xùn)練和驗(yàn)證期的精度差異不明顯.