梁 浩，黃生志，孟二浩，黃 強(qiáng)

（西安理工大學(xué) 西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710048）

1 研究背景

河川徑流作為水循環(huán)的關(guān)鍵環(huán)節(jié)，對(duì)其進(jìn)行準(zhǔn)確預(yù)報(bào)十分重要［1］。當(dāng)前水資源系統(tǒng)受到諸多因素的影響，徑流序列表現(xiàn)出多時(shí)間尺度變化特性和非平穩(wěn)特性，這給徑流序列的預(yù)測(cè)帶來了難度。尤其是在極端天氣頻發(fā)的背景下，徑流序列呈現(xiàn)更大的波動(dòng)性［2］。當(dāng)前，徑流預(yù)測(cè)模型已有不少，從不同的模型中優(yōu)選可靠的徑流預(yù)報(bào)模型，獲得更高的徑流預(yù)測(cè)精度，是水文預(yù)報(bào)工作者迫切開展的工作。

為了獲得未來幾天、幾個(gè)月甚至幾年的徑流預(yù)測(cè)值，水文預(yù)報(bào)人員已經(jīng)提出了多種方法。這些水文預(yù)報(bào)模型方法可以大致分為過程驅(qū)動(dòng)和數(shù)據(jù)驅(qū)動(dòng)兩大塊［3］。過程驅(qū)動(dòng)的模型以基于物理過程的水文模型為主，如王中跟等［4］通過SWAT模型利用地理信息系統(tǒng)和遙感的數(shù)據(jù)模擬流域上復(fù)雜的水文過程。數(shù)據(jù)驅(qū)動(dòng)則主要通過針對(duì)水文站、氣象站等觀測(cè)的數(shù)據(jù)，結(jié)合統(tǒng)計(jì)方法進(jìn)行徑流模擬預(yù)測(cè)，如運(yùn)用廣泛的多元線性回歸（MLR）［5-6］、人工神經(jīng)網(wǎng)絡(luò)模型（ANN）［7-8］和支持向量機(jī)模型（SVM）［9-10］等諸多基于人工智能的模型［11］。水文模型要涉及較為復(fù)雜的水文過程，需要大量準(zhǔn)確的水文過程資料，需要考慮不同時(shí)空尺度下水文系統(tǒng)規(guī)律的差異性，存在參數(shù)不確定性以及水文模型普適性較差等問題［12］。數(shù)據(jù)驅(qū)動(dòng)模型不需要考慮水文過程的物理機(jī)制，對(duì)時(shí)間序列進(jìn)行數(shù)學(xué)分析，在解決數(shù)據(jù)預(yù)測(cè)問題以及時(shí)間序列分析等方面有獨(dú)特的優(yōu)勢(shì)［13-14］。模型對(duì)不同數(shù)據(jù)有其適用性，本文首先優(yōu)選多元線性回歸、人工神經(jīng)網(wǎng)絡(luò)和支持向量機(jī)3種數(shù)據(jù)驅(qū)動(dòng)模型中對(duì)渭河流域徑流序列適用性最強(qiáng)的單一模型。

徑流過程本身具有一定的周期性、趨勢(shì)性以及隨機(jī)性特征，給徑流預(yù)測(cè)帶來難度，所以徑流的多時(shí)間尺度分析和認(rèn)識(shí)對(duì)提高徑流預(yù)測(cè)精度十分必要。分解技術(shù)可以將徑流序列分解為一系列相對(duì)穩(wěn)定的子序列。近年來研究者基于分解技術(shù)進(jìn)行的徑流預(yù)測(cè)研究已經(jīng)有很多，對(duì)徑流進(jìn)行分解重構(gòu)，最終實(shí)現(xiàn)徑流預(yù)測(cè)的方法不斷被應(yīng)用，有經(jīng)驗(yàn)?zāi)B(tài)分解法（EMD）、集合經(jīng)驗(yàn)?zāi)B(tài)分解法（EEMD）和小波分解（WD）等方法。張洪波等［15］探究了對(duì)EMD 分解預(yù)測(cè)后的分量采用不同的重構(gòu)方式及對(duì)于徑流預(yù)報(bào)精度的影響。Tan 等［16］和杜懿等［17］采用EEMD 對(duì)徑流進(jìn)行分解，更好地得到徑流的有用信息，提高預(yù)測(cè)精度。Mohammad 等［18］將EEMD 應(yīng)用于韓國(guó)和伊朗河流的模擬預(yù)測(cè)，也驗(yàn)證了其良好性能。周婷等［19］探討了基于小波分解的SVM 徑流模擬預(yù)測(cè)，徑流預(yù)測(cè)精度得到了提高。以上研究表明了基于分解技術(shù)的徑流研究對(duì)徑流預(yù)測(cè)精度的提高是有效的。然而，大多數(shù)研究者聚焦于采用單一分解方法的混合模型進(jìn)行徑流預(yù)測(cè)研究，缺乏對(duì)基于多種分解方法的混合模型進(jìn)行綜合對(duì)比分析及其適應(yīng)性研究。

以徑流時(shí)間序列為研究對(duì)象的數(shù)據(jù)驅(qū)動(dòng)型徑流預(yù)測(cè)模型，一般情況下模型的輸入項(xiàng)為降雨等氣象因子，同時(shí)考慮歷史的徑流序列。以輸入因子為依據(jù)的預(yù)測(cè)模型可以大致分為3類：基于氣象因素的徑流預(yù)測(cè)模型；基于歷史徑流序列以及氣象因素的徑流預(yù)測(cè)模型；僅考慮歷史徑流序列的徑流預(yù)測(cè)模型。白云天等［20］在對(duì)尼爾基水庫(kù)的入庫(kù)徑流預(yù)報(bào)時(shí)選取的預(yù)報(bào)因子為前期徑流以及降水。變化環(huán)境下水文序列發(fā)生了改變，給徑流預(yù)測(cè)研究帶來了困難。為得到較高的徑流預(yù)測(cè)精度，豐富輸入因子也是關(guān)鍵途徑之一。孟二浩等［21］在徑流預(yù)測(cè)時(shí)豐富了預(yù)報(bào)因子，加入北大西洋年代際震蕩、北極震蕩等因子對(duì)徑流預(yù)測(cè)精度提高起到了作用。楊開斌等［22］探究了在汛期時(shí)ENSO 事件對(duì)于徑流預(yù)測(cè)的影響。相關(guān)的研究表明大氣環(huán)流異常因子對(duì)徑流具有一定影響。值得一提的是，當(dāng)前的研究主要集中于大氣環(huán)流異常因子對(duì)于徑流整體序列的影響，而沒有考慮分解技術(shù)下大氣環(huán)流異常因子對(duì)徑流子序列的影響。

基于此，本文對(duì)比分析了基于EMD、EEMD和WD方法構(gòu)建的多種混合模型進(jìn)行徑流預(yù)報(bào)的精度和效果，并融合大氣環(huán)流異常因子的信息，以期揭示不同混合模型對(duì)河川徑流的預(yù)測(cè)能力，為變化環(huán)境下徑流預(yù)報(bào)精度的提高提供新思路。

2 單一預(yù)測(cè)模型

2.1 多元線性回歸多元線性回歸（MLR）是用于擬合多個(gè)自變量和因變量之間的線性關(guān)系，通過訓(xùn)練確定得出模型需要的參數(shù)，從而確定具體的MLR 方程，之后便可通過此方程來預(yù)測(cè)因變量［23］。其一般形式為：

式中：Y為因變量；βi為回歸系數(shù)，i=0,1,…,k；Xi為自變量，i=1,2,…,k；μ為隨機(jī)誤差項(xiàng)，。回歸系數(shù)β的求解采用最小二乘法，使誤差項(xiàng)的平方和最小。

2.2 人工神經(jīng)網(wǎng)絡(luò)人工神經(jīng)網(wǎng)絡(luò)（ANN）是一種通過模仿動(dòng)物大腦神經(jīng)的行為特征，從而進(jìn)行大規(guī)模分布式并行信息處理的數(shù)學(xué)模型。ANN 由3部分組成：輸入項(xiàng)、包含激勵(lì)函數(shù)的隱含層和輸出項(xiàng)。其依據(jù)的主要思想是：通過對(duì)輸入輸出數(shù)據(jù)的訓(xùn)練，建立能夠以任意精度逼近最佳輸出的非線性函數(shù)的模型［24］。

ANN以神經(jīng)元作為基本處理單元。簡(jiǎn)單的ANN是多輸入一輸出的，其表達(dá)式為：

式中：Ii為隱含層的值；yi為輸出值；xj(j=1,2,…,n) 為不同的輸入信號(hào)；θi為閾值；wji為從神經(jīng)元j到神經(jīng)元i的權(quán)值；f（·）為傳遞函數(shù)。

ANN 通過不斷地調(diào)整權(quán)值和閾值的大小，從而形成一種由輸入到隱含層的固有模式，然后通過傳遞函數(shù)由隱含層到輸出層。合理選擇隱含層數(shù)和神經(jīng)元個(gè)數(shù)對(duì)ANN十分重要。

2.3 支持向量機(jī)支持向量機(jī)（SVM）通過引入核函數(shù)，將原低維空間中非線性不可分的數(shù)據(jù)映射到高維特征空間。將兩類數(shù)據(jù)分類問題轉(zhuǎn)換為在高維特征空間中求最優(yōu)的線性分類面，即找到的分類平面與最靠近的數(shù)據(jù)點(diǎn)之間距離最大［25］。SVM 進(jìn)行模擬預(yù)測(cè)主要思路：通過訓(xùn)練期對(duì)輸入輸出數(shù)據(jù)的訓(xùn)練，建立結(jié)構(gòu)風(fēng)險(xiǎn)最小化的模型，然后給定輸入數(shù)據(jù)就可以得到相應(yīng)的輸出值。通用的SVM 表達(dá)式如下：

式中：xi為輸入值；yi為與xi對(duì)應(yīng)的輸出值；sgn（·）為符號(hào)函數(shù)；為核函數(shù)；為調(diào)控SVM確定最優(yōu)分類平面的參數(shù)。

核函數(shù)的選取是SVM 應(yīng)用的關(guān)鍵之處。常用的核函數(shù)有：線性核函數(shù)、多項(xiàng)式核函數(shù)以及徑向基核函數(shù)等。核函數(shù)計(jì)算公式如下：

式中?函數(shù)將數(shù)據(jù)從低維空間映射到高維特征空間。

SVM 用于回歸分析，能夠?qū)τ?xùn)練數(shù)據(jù)進(jìn)行學(xué)習(xí)，并實(shí)現(xiàn)其余數(shù)據(jù)的模擬預(yù)測(cè)，很好的解決了非線性擬合問題。對(duì)于樣本數(shù)據(jù)，xi、yi分別為輸入變量和預(yù)測(cè)變量，n為數(shù)據(jù)量［26］。其回歸函數(shù)為：

式中：φ(x) 為輸入空間到高維空間的非線性映射；ω為權(quán)值向量；b為偏差。

SVM 基于結(jié)構(gòu)風(fēng)險(xiǎn)最小化原理，而不是經(jīng)驗(yàn)風(fēng)險(xiǎn)最小化，提高了模型的泛化能力，求得風(fēng)險(xiǎn)最小的f(x)：

式中：ε為損失系數(shù)；ξi、為松弛變量；C為懲罰因子。

3 分解技術(shù)

3.1 經(jīng)驗(yàn)?zāi)B(tài)分解經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）簡(jiǎn)言之就是將一段時(shí)間序列依據(jù)指定的條件分解為多個(gè)本征模態(tài)分量和一個(gè)趨勢(shì)項(xiàng)，使得每個(gè)本征模態(tài)分量的極值點(diǎn)和過零點(diǎn)數(shù)相等或者最多相差1，且相對(duì)于局部零值點(diǎn)是對(duì)稱的［27］。其具體步驟如下：

（1）輸入需要處理的原始時(shí)間序列y(t) ；

（2）找出序列中所有的極大值點(diǎn)和極小值點(diǎn)，根據(jù)三次樣條函數(shù)分別擬合出上包絡(luò)線yp(t)和下包絡(luò)線并求出上下包絡(luò)線的均值序列(t) 。計(jì)算公式如下：

（4）計(jì)算剩余項(xiàng)z1(t)，以z1(t)作為原始序列，重復(fù)第（2）和第（3）步。如果最后的剩余項(xiàng)zi(t)滿足單調(diào)函數(shù)或至多具有一個(gè)局部極值點(diǎn)，此時(shí)zi(t)=q(t)為趨勢(shì)項(xiàng)，完成整個(gè)分解。否則繼續(xù)進(jìn)行第（1）—（3）步。

式中r'i(t)為滿足本征模態(tài)分量條件的序列。

EMD分解原始序列y(t)得到n個(gè)本征模態(tài)分量和趨勢(shì)項(xiàng)?？杀硎緸椋?/p>

式中：ri(t)為原始時(shí)間序列的本征模態(tài)分量，i=1,2,…,n；q(t)為原始時(shí)間序列的趨勢(shì)項(xiàng)。

3.2 集合經(jīng)驗(yàn)?zāi)B(tài)分解集合經(jīng)驗(yàn)?zāi)B(tài)分解（EEMD）與EMD 的目的是相同的，將原有的時(shí)間序列分解得到若干不同尺度的本征模態(tài)分量IMF 以及一個(gè)趨勢(shì)項(xiàng)。Huang 等［28］發(fā)現(xiàn)EMD 分解得到的某些本征模態(tài)分量IMF 包含了不止一個(gè)尺度下的信號(hào)，有些信號(hào)出現(xiàn)在不同的IMF 中。通過加入白噪音改進(jìn)EMD 得到的EEMD 解決了模態(tài)混疊的問題，能夠更精確的分解得到不同尺度的本征模態(tài)分量和趨勢(shì)項(xiàng)。EEMD分解時(shí)間序列的步驟如下：

（1）輸入需要處理的原始時(shí)間序列y(t)；

（2）產(chǎn)生m組白噪聲序列nj(t) (j=1,2,…,m)，要求m組白噪聲的均值為零。分別加入原始序列中，序列變?yōu)閦j(t)：

（3）對(duì)zj(t)進(jìn)行EMD分解，得到n個(gè)本征模態(tài)分量rij(t)以及趨勢(shì)項(xiàng)qj(t)：

（4）由第（3）步可以得到m組n個(gè)本征模態(tài)分量以及n個(gè)趨勢(shì)項(xiàng)，求各項(xiàng)均值為序列分解的最終分量；

可以看出，EEMD同樣將時(shí)間序列y(t)分為了若干的本征模態(tài)分量和趨勢(shì)項(xiàng)。公式如下：

式中：r'i(t)為原始時(shí)間序列的若干本征模態(tài)分量；q'(t)為原始時(shí)間序列的趨勢(shì)項(xiàng)。

3.3 小波分解小波分析（WD）是一種對(duì)時(shí)間序列在時(shí)域和頻域上進(jìn)行逐級(jí)分析的方法［29］。小波是均值為0，長(zhǎng)度有限的波形。滿足如下條件的函數(shù)為小波母函數(shù)：

其中：ψ(ω)為小波母函數(shù)；ψ(t)為在頻率ω處的Fourier變換；ψ*(ω)為ψ(ω)的共軛函數(shù)。

小波變換是在時(shí)間和頻率兩域?qū)π盘?hào)的局部特征進(jìn)行表征。在R上平方可積函數(shù)構(gòu)成的函數(shù)空間L2(R)，f(t)∈L2(R)。則關(guān)于母小波函數(shù)ψ的連續(xù)小波變換公式為：

對(duì)于離散的信號(hào)，有離散小波變換，離散連續(xù)小波變換中的伸縮因子a和平移因子b得到，通常取。公式如下：

1989年，塔式算法融入小波變換中，實(shí)現(xiàn)信號(hào)的分解重構(gòu)，即著名的Mallat 算法。分解和重構(gòu)過程表達(dá)式為：

式中：cj為原始時(shí)間序列；h、g分別為低通、高通濾波器；cj-1、dj-1分別為分解尺度j-1下的低頻信號(hào)與高頻信號(hào)；分別為h、g的共軛反轉(zhuǎn)；的卷積；為對(duì)卷積的二次下抽樣；Ucj-1為對(duì)cj-1的二元上抽樣；U、D為二元上、下抽樣算子。

分解、重構(gòu)后得到與原時(shí)間系列同長(zhǎng)度的子序列，便于對(duì)時(shí)間序列進(jìn)行預(yù)測(cè)。WD 方法應(yīng)用于時(shí)間序列預(yù)測(cè)的步驟如下：（1）選擇合理的小波函數(shù)，并選擇合適的時(shí)間尺度（或小波分解水平）；（2）通過濾波器對(duì)時(shí)間序列進(jìn)行分解，時(shí)間序列分解得到一個(gè)逼近信號(hào)序列和若干個(gè)細(xì)節(jié)信號(hào)序列；（3）利用小波重構(gòu)方法可以得到與原序列長(zhǎng)度相同子序列，最后進(jìn)行時(shí)間序列模擬預(yù)測(cè)等工作。

結(jié)合以上模型方法，本文通過多站優(yōu)選MLR、ANN、SVM 模型，選取最優(yōu)的單一模型并基于此模型結(jié)合多種分解技術(shù)（EMD、EEMD、WD）構(gòu)建混合預(yù)測(cè)模型，對(duì)比分析不同站點(diǎn)采用混合模型的徑流預(yù)測(cè)精度高低；探究變化環(huán)境下豐富輸入因子（大氣環(huán)流異常因子）對(duì)于混合模型徑流預(yù)測(cè)精度的影響。主要技術(shù)路線如圖1所示。

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

4.1 研究區(qū)域概況以渭河流域?yàn)檠芯繀^(qū)域，選取北道、林家村、華縣3 個(gè)水文站作為研究對(duì)象。

圖1 技術(shù)路線

其中，北道位于渭河上游，林家村位于渭河中游，華縣站位于渭河下游。北道位于甘肅省內(nèi)，地處秦嶺北麓，集水面積1846 km2。渭河支流牛頭河匯集湯浴河、后川河等諸多小河于天水市北道區(qū)注入渭河。林家村水文站位于陜西省寶雞市，其集水面積30 661 km2，多年平均含沙量63.7 kg/m3，是渭河干流的控制站，屬于國(guó)家重點(diǎn)水文站。華縣站位于陜西省渭南市，集水面積106 498 km2，多年平均含沙量49.3 kg/m3，是國(guó)家重點(diǎn)水文站，同時(shí)也是渭河流域干流的控制站。渭河流域示意圖如圖2。

圖2 渭河流域

采用的數(shù)據(jù)為北道、林家村、華縣3 個(gè)水文站1960年1月至2010年12月徑流，3 個(gè)水文站所在區(qū)域內(nèi)相應(yīng)時(shí)間的降雨、蒸發(fā)以及大氣環(huán)流異常因子（北極濤動(dòng)（Arctic Oscillation，AO）、太平洋十年濤動(dòng)（Pacific Decadal Oscillation，PDO）、Nino3.4、太陽(yáng)黑子）的月尺度數(shù)據(jù)。其中，以徑流作為預(yù)測(cè)目標(biāo)，以降雨、蒸發(fā)和前期的徑流作為模型的輸入項(xiàng)，將大氣環(huán)流異常因子作為輸入項(xiàng)增加項(xiàng)。

本文收集的4 種具有代表性的大氣環(huán)流異常因子月尺度數(shù)據(jù)，其中，AO 數(shù)據(jù)從美國(guó)國(guó)家海洋和大氣管理（National Oceanic and Atmospheric Administrator，NOAA）國(guó)家氣象數(shù)據(jù)中心獲取，網(wǎng)址為（http：//www.ncdc.noaa.gov/teleconnections/ao.php）；PDO 數(shù)據(jù)從東京氣候中心獲取，網(wǎng)址為（http：//ds.data/jma.go.jp/tcc/products/elnino/decadal/pdo.html）；Nino3.4 從NOAA 的地球系統(tǒng)研究實(shí)驗(yàn)室獲得，網(wǎng)址為（http：//www.esrl.noaa.gov/psd/data/correlation/nina34.data）；太陽(yáng)黑子數(shù)據(jù)從NOAA的國(guó)家地球物理數(shù)據(jù)中心獲得，網(wǎng)址為（https：//www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/SUNSPOT/）。

由于降雨蒸發(fā)等因素對(duì)徑流的影響存在滯后效應(yīng)，不同因子對(duì)徑流的影響時(shí)間存在差異，故本文基于相關(guān)系數(shù)對(duì)模型輸入項(xiàng)進(jìn)行優(yōu)選。通過對(duì)于不同滯時(shí)下，各因子與目標(biāo)因子的相關(guān)性計(jì)算，分析選取合理的因子，這對(duì)徑流的預(yù)測(cè)是至關(guān)重要的。由于在徑流預(yù)測(cè)中，考慮了滯時(shí)為一個(gè)月的徑流作為預(yù)測(cè)下一月徑流的輸入項(xiàng)，所以徑流預(yù)測(cè)的預(yù)見期為1個(gè)月。

4.2 徑流序列突變分析以及時(shí)間尺度特性分析采用Mann-Kendall 變異診斷法對(duì)各研究站點(diǎn)徑流序列進(jìn)行趨勢(shì)分析以及變異點(diǎn)診斷，結(jié)果如圖3所示。

由圖3（a）北道站的UF曲線可知：自1970年9月起北道站徑流量有明顯的減少趨勢(shì)，且在1986年8月以后這種減少趨勢(shì)超過了0.05的顯著性水平，意味著北道站徑流量下降趨勢(shì)顯著。同時(shí)可以看到北道站UF與UB曲線交點(diǎn)位于1988年1—4月，不在信度線之間，即徑流量沒有發(fā)生顯著的突變。由圖3（b）林家村站的UF曲線可知：自1963年11月起林家村站徑流量有減少趨勢(shì)，在1989年11月后減少趨勢(shì)超過顯著性水平0.05的臨界線，并且自1992年10月開始趨勢(shì)出現(xiàn)更明顯的下降。林家村站UF與UB曲線交點(diǎn)位于1989年1月不在信度線之間，即徑流量于此月發(fā)生突變，但交點(diǎn)已很接近信度線。由圖3（c）華縣站的UF曲線可知：自1963年11月起華縣站徑流量有減少趨勢(shì)，在1980年7月至1984年11月之間有增加趨勢(shì)，之后徑流量為持續(xù)減少趨勢(shì)且減少趨勢(shì)超過顯著性水平0.05的信度線。華縣站UF與UB曲線交點(diǎn)位于1971年11月且交點(diǎn)位于信度線之間，即此月徑流量發(fā)生由多到少的突變。

較長(zhǎng)時(shí)期的河川徑流變化受多種因素的影響，是氣候變化和人類活動(dòng)因素（如修建大量水利工程等）綜合作用產(chǎn)生的結(jié)果［30］。渭河流域上的3個(gè)研究站點(diǎn)月徑流序列均有趨勢(shì)性變化，且華縣站月徑流序列發(fā)生了變異。渭河流域有大量灌區(qū)和水庫(kù)水電站，主要灌區(qū)集中在關(guān)中平原，關(guān)中九大灌區(qū)依渭河橫貫關(guān)中平原［31］。灌區(qū)灌溉用水量波動(dòng)不大，且大多灌區(qū)建成時(shí)間較早。渭河流域的大中型水庫(kù)修建時(shí)間為1960—1980年代［32］。理論上認(rèn)為模型可以在訓(xùn)練期捕捉并學(xué)習(xí)到灌溉工程用水以及水庫(kù)運(yùn)行等人類活動(dòng)信息。選取此3個(gè)研究站點(diǎn)，以期探究不同混合模型對(duì)變化環(huán)境下徑流預(yù)測(cè)的效果及適用性，并且可以更真實(shí)的認(rèn)識(shí)大氣環(huán)流異常對(duì)于分解后徑流子序列的影響。

圖3 各站徑流突變點(diǎn)檢驗(yàn)UF、UB過程線

考慮滯后效應(yīng)以及徑流變化可能對(duì)模型效果的影響，故劃分本文的訓(xùn)練期為1960年1月至2000年12月，驗(yàn)證期為2001年1月至2010年12月。在模型的訓(xùn)練期，模型將捕捉并學(xué)習(xí)到氣候變化和人類活動(dòng)因素的部分信息，并相應(yīng)的體現(xiàn)在模型的參數(shù)變化上。在訓(xùn)練期實(shí)現(xiàn)對(duì)于模型參數(shù)的確定，驗(yàn)證期為預(yù)測(cè)驗(yàn)證樣本進(jìn)行徑流預(yù)測(cè)［33］。在徑流模型評(píng)價(jià)體系建立中，選取了納什系數(shù)和均方根誤差為衡量徑流預(yù)測(cè)精度的評(píng)價(jià)指標(biāo)。

4.3 選用MLR、ANN、SVM 模型進(jìn)行徑流預(yù)測(cè)本文基于相關(guān)系數(shù)選取合理的輸入因子，分別計(jì)算北道、林家村、華縣3個(gè)水文站徑流、降雨、蒸發(fā)與徑流的相關(guān)系數(shù)分析，結(jié)果如圖4。

圖4 基于相關(guān)系數(shù)的輸入因子優(yōu)選

基于相關(guān)系數(shù)優(yōu)選輸入因子（降雨、蒸發(fā)、前期的徑流），相關(guān)系數(shù)大于0.3即滿足95%置信度檢驗(yàn)，滿足置信度檢驗(yàn)的序列作為預(yù)選輸入項(xiàng)［25］。可以看出，華縣站徑流自滯時(shí)為1 時(shí)，相關(guān)系數(shù)大于0.3；降雨與徑流的滯時(shí)為0、1、2、7、8個(gè)月的相關(guān)系數(shù)絕對(duì)值均大于0.3；蒸發(fā)與徑流的滯時(shí)為1、2、3、7、8、9 個(gè)月的相關(guān)系數(shù)絕對(duì)值大于0.3。徑流預(yù)報(bào)模型輸入項(xiàng)需要高信息量，強(qiáng)相關(guān)性的因子。一些因子雖然與徑流相關(guān)性好，但它們的作用是重復(fù)的。這樣的因子并不能同其他因子配合共同對(duì)徑流起作用，結(jié)果可能使預(yù)測(cè)精度降低，故需要剔除［34］。本文通過窮舉法［35］獲得最佳的輸入項(xiàng)。如華縣站的最佳輸入因子有：滯時(shí)為1 的徑流，滯時(shí)為0、1、2 個(gè)月的降雨以及滯時(shí)為1、2、3、7 個(gè)月的蒸發(fā)。同理優(yōu)選得到北道站、林家村站的最佳輸入因子，然后以輸入因子、目標(biāo)因子建立單一模型（MLR、ANN、SVM）。本文選用ANN 中的BP 三層神經(jīng)網(wǎng)絡(luò)，模型訓(xùn)練函數(shù)采用trainlm［36］。對(duì)于隱含層神經(jīng)元個(gè)數(shù)，本文在訓(xùn)練過程預(yù)先選定一個(gè)較小的值，然后依次增加神經(jīng)元數(shù)，取驗(yàn)證數(shù)據(jù)誤差達(dá)到最小時(shí)的隱層神經(jīng)元數(shù)為最終隱層單元數(shù)。由于徑向基核函數(shù)具有穩(wěn)定性和去記憶性，故SVM 模型構(gòu)建中，對(duì)比選用了徑向基核函數(shù)，不敏感損失系數(shù)ε取0.001，采用網(wǎng)格尋優(yōu)方法確定懲罰因子C和核函數(shù)參數(shù)g［37］。在預(yù)測(cè)驗(yàn)證期，將相應(yīng)滯時(shí)下的前期徑流、降雨、蒸發(fā)作為輸入項(xiàng)，輸入已訓(xùn)練好的模型，可得到預(yù)測(cè)徑流?；诓煌瑔我荒Ｐ偷膹搅黝A(yù)測(cè)結(jié)果，如表1所示。

通過表1 可以得出，對(duì)于不同站點(diǎn)，采用MLR、ANN、SVM 模型預(yù)測(cè)徑流精度有較大差距。林家村站體現(xiàn)尤為明顯，采用MLR 模型驗(yàn)證期其納什系數(shù)僅為0.150，均方根誤差為0.733；采用ANN模型，驗(yàn)證期其納什系數(shù)為0.332，均方根誤差為0.642。據(jù)調(diào)查可知，林家村站上游有大型灌區(qū)水庫(kù)，需要取水調(diào)水，人為因素對(duì)徑流的影響很強(qiáng)，使得林家村站徑流預(yù)測(cè)模型中，MLR 模型的泛化能力弱，相比采用描述線性關(guān)系的MLR 模型，采用模擬非線性關(guān)系的ANN 模型在處理徑流及其影響因素之間復(fù)雜模糊的關(guān)系中更具有優(yōu)勢(shì)，泛化能力更強(qiáng)，預(yù)測(cè)時(shí)能得到更高的精度。采用SVM 模型預(yù)測(cè)林家村站徑流，驗(yàn)證期其納什系數(shù)為0.603，相較ANN 模擬結(jié)果提高了81.6%；相應(yīng)的均方根誤差為0.495，降低了22.9%，其模擬效果優(yōu)于ANN的模擬效果，說明SVM相比ANN模型在預(yù)測(cè)徑流過程中更具有優(yōu)勢(shì)。從模型原理探究原因，ANN 模型通過訓(xùn)練后建立了基于經(jīng)驗(yàn)風(fēng)險(xiǎn)最小化的網(wǎng)絡(luò)模型，存在著局部極小值以及具有不穩(wěn)定性等缺點(diǎn)，而SVM 模型引用核函數(shù)，并且以結(jié)構(gòu)風(fēng)險(xiǎn)最小化為原則，最終得到確定的唯一解，可以彌補(bǔ)ANN 的一些不足，所以模擬預(yù)測(cè)效果優(yōu)于ANN。北道站以及華縣站的徑流預(yù)測(cè)結(jié)果分析亦可得到相同的結(jié)論。

表1 基于MLR、ANN和SVM模型的徑流預(yù)測(cè)

4.4 基于混合模型的徑流預(yù)測(cè)本文對(duì)北道、華縣和林家村站的徑流分別采用EMD、EEMD 和WD技術(shù)處理。EMD分解方法自適應(yīng)的將北道站、華縣站和林家村站分解得到的子序列個(gè)數(shù)分別為7、6、7。子序列分別包含各站徑流的高頻部分、低頻部分以及趨勢(shì)項(xiàng)。EEMD分解方法將北道站、華縣站、林家村站分解后，得到子序列數(shù)分別為7、8、8。由于徑流序列是離散的，所以采用離散小波分解徑流。對(duì)于WD方法，小波函數(shù)以及分解層數(shù)的選擇極為關(guān)鍵。通過對(duì)比優(yōu)選，本文選取的小波函數(shù)為Daubechies小波系中的db3［38］，小波分解層數(shù)為5層，子序列有1個(gè)逼近信號(hào)（a5）和5個(gè)細(xì)節(jié)信號(hào)（d1、d2、d3、d4、d5）。

耦合不同分解方法與SVM 建立混合預(yù)測(cè)模型，即根據(jù)不同分解方法將徑流序列分解為多個(gè)子序列；將各子序列分別建立預(yù)測(cè)模型進(jìn)行單獨(dú)模擬，最后將預(yù)測(cè)的各子序列值進(jìn)行累加，從而得到最終的預(yù)測(cè)徑流序列。以降雨、蒸發(fā)以及前期徑流為輸入因子，考慮各因子與各徑流子序列的滯后性。對(duì)于各徑流子序列進(jìn)行模擬預(yù)測(cè)時(shí)，采用窮舉法選取有效滯時(shí)下的影響因子作為模型的輸入因子。將預(yù)測(cè)精度最高的徑流子序列累加合成完整預(yù)測(cè)徑流序列。同樣，對(duì)于混合模型的預(yù)測(cè)徑流精度通過納什系數(shù)和均方根誤差來評(píng)價(jià)?；诨旌夏Ｐ偷膹搅黝A(yù)測(cè)結(jié)果見表2和圖5。

通過表2基于混合模型的徑流預(yù)測(cè)結(jié)果，可以看出不同研究站點(diǎn)的徑流序列經(jīng)過不同的方法分解后，子序列通過SVM 預(yù)測(cè)合成后的徑流精度有明顯差別。EEMD-SVM 分解預(yù)測(cè)的精度要普遍高于EMD-SVM 分解預(yù)測(cè)的精度。北道、華縣、林家村徑流預(yù)測(cè)精度提高幅度分別為21.6%、1.9%、4.8%。EEMD 方法改進(jìn)了EMD 方法所存在的模態(tài)混雜問題，在時(shí)間序列分解上更有優(yōu)勢(shì)。WD-SVM分解預(yù)測(cè)的精度要高于EEMD-SVM 分解預(yù)測(cè)精度。但是不同站點(diǎn)區(qū)別很大，EEMD-SVM 模型在華縣站徑流預(yù)測(cè)中泛化能力較弱，采用WD-SVM模型進(jìn)行分解預(yù)報(bào)相較EEMD-SVM 模型徑流精度提高幅度為31.0%。而北道站，WD-SVM 與EEMD-SVM 徑流預(yù)測(cè)精度相差很小。姜浩等［39］的研究發(fā)現(xiàn)：依據(jù)通量矩陣可得WD 方法在分解能力上優(yōu)于EMD。由以上結(jié)果可知，不同徑流序列對(duì)于分解方法有著不同的適應(yīng)性，但一般來說預(yù)測(cè)精度高低符合如下表達(dá)式：WD-SVM＞EEMD-SVM＞EMD-SVM。

表2 基于混合模型的徑流預(yù)測(cè)

施能［40］的研究表明大氣環(huán)流是造成我國(guó)冬季氣候變化的重要原因之一。本文研究區(qū)域?yàn)槲己恿饔?，海陸熱力的差異?huì)給當(dāng)?shù)貧夂蛟斐奢^大影響，同時(shí)可能會(huì)對(duì)當(dāng)?shù)貜搅鳟a(chǎn)生影響。所以本文的輸入因子在降雨、蒸發(fā)、前期徑流的基礎(chǔ)上，考慮大氣環(huán)流異常因子（AO、PDO、Nino3.4、太陽(yáng)黑子）作為徑流預(yù)測(cè)的輸入因子，并對(duì)增加大氣環(huán)流異常因子前后徑流預(yù)測(cè)精度進(jìn)行了對(duì)比分析。原始徑流采用WD 方法分解得到子序列，通過相關(guān)系數(shù)計(jì)算大氣環(huán)流異常因子相對(duì)于各分解序列的滯時(shí)，優(yōu)選大氣環(huán)流異常因子作為輸入項(xiàng)，采用SVM 模型進(jìn)行預(yù)測(cè)。表3 為基于WD-SVM 混合模型的大氣環(huán)流異常因子篩選結(jié)果。

如表3 所示，對(duì)于研究站點(diǎn)大氣環(huán)流異常因子的篩選結(jié)果可知：4 個(gè)大氣環(huán)流異常因子中，Nino3.4 對(duì)渭河徑流影響明顯，而AO、PDO 和太陽(yáng)黑子的影響不明顯；北道、林家村和華縣站依次由西向東分布，季風(fēng)對(duì)中國(guó)東部和中國(guó)西部的影響強(qiáng)度不同，且季風(fēng)也會(huì)對(duì)Nino3.4 產(chǎn)生影響，故Nino3.4 對(duì)于不同水文站徑流的影響會(huì)不同［41］；對(duì)北道站經(jīng)小波分解后的徑流子序列d3 進(jìn)行預(yù)測(cè)時(shí)，增加滯時(shí)為3、4、5、8、9、10 的Nino3.4 作為模型輸入項(xiàng)；累加分解預(yù)測(cè)后的徑流子序列，徑流精度得到了提高。針對(duì)小波分解后的各徑流子序列通過支持向量機(jī)進(jìn)行預(yù)測(cè)時(shí)，支持向量機(jī)中懲罰因子C和核函數(shù)參數(shù)g結(jié)果如表4所示。

圖5 基于分解的混合模型及融合大氣環(huán)流異常因子的WD-SVM模型的徑流預(yù)測(cè)

最終徑流預(yù)測(cè)結(jié)果如表5 所示。驗(yàn)證期混合模型以及融合大氣環(huán)流異常因子的WD-SVM 模型預(yù)測(cè)徑流與實(shí)測(cè)徑流對(duì)比情況見圖5。

由表5 基于WD-SVM 模型豐富輸入因子后的徑流預(yù)測(cè)，可知相比輸入因子只有降雨、徑流以及前期徑流，增加大氣環(huán)流異常因子對(duì)3站的徑流預(yù)測(cè)精度均有提高。北道站徑流預(yù)測(cè)中納什系數(shù)提高了0.25%，均方根誤差降低了0.78%；林家村站納什系數(shù)提高了0.23%，均方根誤差降低了1.0%；華縣站納什系數(shù)提高了0.7%，均方根誤差降低了16%。其次，對(duì)各研究站點(diǎn)徑流極大值點(diǎn)（10 個(gè)左右）的預(yù)測(cè)精度進(jìn)行了分析。WD-SVM 混合模型融入大氣環(huán)流異常因子后，北道站其均方根誤差減小0.9%；林家村站其均方根誤差減小1.3%；華縣站其均方根誤差減小4.7%。結(jié)果印證基于分解技術(shù)下大氣環(huán)流異常因子對(duì)于徑流子序列有影響，豐富輸入因子可以提高徑流預(yù)測(cè)精度，尤其是對(duì)于徑流極大值點(diǎn)的預(yù)測(cè)精度提高更為明顯。

表3 基于WD-SVM混合模型的大氣環(huán)流異常因子篩選結(jié)果

表4 基于小波分解的徑流子序列支持向量機(jī)模型參數(shù)值

表5 基于WD-SVM混合模型融合大氣環(huán)流異常因子前后的徑流預(yù)測(cè)

綜上所述，豐富輸入因子后的WD-SVM 模型對(duì)研究站點(diǎn)的徑流預(yù)測(cè)精度均高于80%。林家村站和華縣站訓(xùn)練期以及驗(yàn)證期的徑流預(yù)測(cè)精度不低于北道站（水利工程影響較?。Ｐ藿ǖ拇罅克こ滩⑽磳?duì)徑流預(yù)測(cè)造成較大影響。其他學(xué)者也做了不少有關(guān)渭河流域徑流模擬預(yù)測(cè)的研究，左德鵬等［42］基于SWAT 模型和SUFI-2 算法的研究結(jié)果中，林家村站和華縣站月徑流模擬結(jié)果為：林家村站校準(zhǔn)期和驗(yàn)證期納什系數(shù)為0.23 和0.75、華縣站校準(zhǔn)期和驗(yàn)證期納什系數(shù)為0.72 和0.78；Zhao 等［43］基于分布式水文模型DTVGM 得到咸陽(yáng)站（咸陽(yáng)站位于渭河干流，林家村站和華縣站之間）月模擬徑流和觀測(cè)徑流納什系數(shù)在0.6～0.7 之間；黎云云等［44］基于SWAT 模型的渭河流域徑流模擬中，林家村站校準(zhǔn)期和驗(yàn)證期納什系數(shù)為0.58 和0.80、華縣站校準(zhǔn)期和驗(yàn)證期納什系數(shù)為0.77 和0.82。而本文基于WD-SVM 混合模型徑流預(yù)測(cè)結(jié)果中，林家村站校準(zhǔn)期和驗(yàn)證期納什系數(shù)0.880 和0.868、華縣站校準(zhǔn)期和驗(yàn)證期納什系數(shù)為0.890 和0.855。通過單一預(yù)測(cè)模型和混合預(yù)測(cè)模型預(yù)測(cè)結(jié)果、以及與其他學(xué)者相關(guān)研究結(jié)果對(duì)比可知，基于分解技術(shù)來提高徑流預(yù)測(cè)精度的有效性以及變化環(huán)境下融合大氣環(huán)流異常因子會(huì)進(jìn)一步提高徑流預(yù)測(cè)精度。

5 結(jié)論

（1）本文采用MLR、ANN 和SVM 模型對(duì)渭河流域北道、林家村、華縣站徑流進(jìn)行預(yù)測(cè)，結(jié)果表明，依次采用MLR、ANN、SVM 模型進(jìn)行徑流預(yù)測(cè)，預(yù)測(cè)精度不斷提高。表明渭河流域上徑流與其影響因素之間復(fù)雜的非線性關(guān)系占較大比重，在對(duì)徑流預(yù)測(cè)過程中基于結(jié)構(gòu)風(fēng)險(xiǎn)最小化的SVM 要優(yōu)于基于經(jīng)驗(yàn)風(fēng)險(xiǎn)最小化的ANN。

（2）混合模型的徑流預(yù)測(cè)研究中，3 種混合模型的徑流預(yù)測(cè)精度由高到低排序?yàn)椋篧D-SVM＞EEMD-SVM＞EMD-SVM。

（3）WD-SVM 混合模型中，預(yù)報(bào)因子加入大氣環(huán)流異常因子，各站點(diǎn)的模型預(yù)測(cè)精度均得到了提高，尤其是對(duì)于徑流極大值點(diǎn)的提升較為明顯。最終，北道、林家村、華縣站驗(yàn)證期徑流的納什系數(shù)分別為0.818、0.868和0.855。

（4）在對(duì)渭河流域不同站點(diǎn)徑流預(yù)測(cè)研究中，相較于直接對(duì)徑流序列預(yù)測(cè)，基于分解的混合模型預(yù)測(cè)精度更高。采用徑流分解技術(shù)后可更好地呈現(xiàn)出徑流的周期、趨勢(shì)等本質(zhì)特征，對(duì)提高徑流預(yù)測(cè)精度起促進(jìn)作用。在接下來的研究工作中，將有3個(gè)突破點(diǎn)：①大氣環(huán)流異常因子對(duì)徑流影響機(jī)理的深入研究，實(shí)現(xiàn)作為預(yù)報(bào)因子的有效辨別和預(yù)處理；②徑流預(yù)測(cè)的模型進(jìn)一步優(yōu)化，提高其學(xué)習(xí)、泛化能力；③分解方法的優(yōu)化，對(duì)徑流序列進(jìn)行更為本質(zhì)的剖析，提高徑流的預(yù)測(cè)精度。