基于BMA多模型組合的疏勒河徑流預(yù)測(cè)研究

周 婷，溫小虎，馮 起，尹振良，楊林山

（1.中國科學(xué)院西北生態(tài)環(huán)境資源研究院，甘肅 蘭州 730000；2.中國科學(xué)院大學(xué)，北京 100049）

0 引言

準(zhǔn)確可靠的日徑流預(yù)測(cè)是水文研究的重要內(nèi)容，也是水資源管理的重要環(huán)節(jié)［1］，可為流域內(nèi)水資源的科學(xué)管理、合理利用以及水資源評(píng)估提供重要的科學(xué)依據(jù)［2］。其中，3日以上的中長(zhǎng)期日徑流預(yù)測(cè)對(duì)于流域內(nèi)的防洪、水資源調(diào)控、生態(tài)環(huán)境保護(hù)等具有十分重要的意義［3］。特別對(duì)于我國西北干旱半干旱地區(qū)，水資源已成為限制區(qū)域社會(huì)經(jīng)濟(jì)發(fā)展、影響生態(tài)安全的主要因素［4］。因此，準(zhǔn)確的徑流預(yù)測(cè)對(duì)于區(qū)域內(nèi)的水資源管理與可持續(xù)利用都具有非常重要的現(xiàn)實(shí)意義［5-6］。

分布式水文模型與機(jī)器學(xué)習(xí)方法是徑流模擬預(yù)測(cè)的常用工具［7-8］。分布式水文模型可以準(zhǔn)確地描述流域內(nèi)的水文物理過程，但通常需要大量的水文、氣象等數(shù)據(jù)作為輸入，而西北地區(qū)由于觀測(cè)站點(diǎn)稀疏，水文、氣象資料不足，使用分布式模型會(huì)使建模過程既耗時(shí)又昂貴［3，9］。機(jī)器學(xué)習(xí)方法僅從數(shù)學(xué)上考慮輸入與輸出之間的非線性關(guān)系，不需要描述水文過程的物理參數(shù)，通過挖掘數(shù)據(jù)本身的潛在規(guī)律對(duì)徑流進(jìn)行預(yù)測(cè)。近年來，機(jī)器學(xué)習(xí)由于算法簡(jiǎn)單、易于實(shí)現(xiàn)等優(yōu)點(diǎn)，在資料有限的條件下，展現(xiàn)出更優(yōu)的預(yù)測(cè)能力，已經(jīng)廣泛應(yīng)用于各種條件下的徑流模擬預(yù)測(cè)研究中［10］。其中，極限學(xué)習(xí)機(jī)（Extreme Learning Machine，ELM）模型、支持向量機(jī)（Support Vector Machine，SVM）模型、多元自適應(yīng)回歸樣條（Multivariate Adaptive Regression Spline，MARS）模型等在西北地區(qū)徑流預(yù)測(cè)研究中均取得了較好的預(yù)測(cè)結(jié)果［9-11］。雖然單一模型在特定的時(shí)空條件下都表現(xiàn)了較好的模擬預(yù)測(cè)性能。但是，由于徑流過程的高度非線性、形成機(jī)制的復(fù)雜性以及模型結(jié)構(gòu)自身的不確定性，單一模型的模擬預(yù)測(cè)結(jié)果存在較大不確定性［12］。

針對(duì)徑流預(yù)測(cè)研究通常采用單一方法進(jìn)行建模與預(yù)測(cè)，難以有效利用各預(yù)測(cè)模型優(yōu)勢(shì)的問題，有學(xué)者提出多模型融合的方法，即將多個(gè)單一模型進(jìn)行組合，通過降低單一模型的偏差，發(fā)揮各模型的優(yōu)勢(shì)，使預(yù)測(cè)值更加接近于實(shí)測(cè)值，進(jìn)而提高徑流預(yù)測(cè)結(jié)果的精確度與可靠性［13］。Zhang等［14］運(yùn)用奇異譜分析（SSA）與自回歸綜合移動(dòng)平均模型（ARIMA）結(jié)合對(duì)年徑流量進(jìn)行預(yù)測(cè)。結(jié)果表明組合模型具有更佳的預(yù)測(cè)能力。Wang等［15］將人工神經(jīng)網(wǎng)絡(luò)模型（ANN）與集成經(jīng)驗(yàn)?zāi)Ｊ椒纸猓‥EMD）結(jié)合用于中長(zhǎng)期徑流時(shí)間序列預(yù)測(cè)，其結(jié)果較單一ANN模型有了明顯的改進(jìn)。Shamseldin等［16］運(yùn)用簡(jiǎn)單平均法（SAM）、加權(quán)平均法（WAM）、神經(jīng)網(wǎng)絡(luò)方法（NNM）對(duì)5個(gè)水文模型進(jìn)行組合用于日徑流量預(yù)測(cè)，結(jié)果表明三組組合模型的預(yù)測(cè)精度均高于單一模型。以上方法均基于確定性的理論，集成不同模型的預(yù)測(cè)結(jié)果，能夠提供更精確的綜合預(yù)測(cè)結(jié)果，但無法定量評(píng)價(jià)模型結(jié)構(gòu)的不確定性［12］。基于貝葉斯理論發(fā)展起來的貝葉斯模型平均方法（Bayesian Model Averaging，BMA）通過計(jì)算單一模型的后驗(yàn)概率，并基于后驗(yàn)概率對(duì)組合模型成員進(jìn)行加權(quán)平均，判斷成員模型優(yōu)劣。因此，BMA可以有效處理組合模型成員的不確定性，在提供精度更高的預(yù)測(cè)結(jié)果的同時(shí)，還能提供可靠的預(yù)測(cè)概率，定量評(píng)價(jià)模型結(jié)構(gòu)不確定性對(duì)預(yù)測(cè)結(jié)果的影響［17］。近年來，BMA方法也被應(yīng)用于水文學(xué)領(lǐng)域，董磊華等［18］運(yùn) 用BMA方 法 綜 合SIMHYD、SMAR（Soil Moisture Accounting and Routing），以及新安江模型的預(yù)報(bào)結(jié)果，推求出可靠的預(yù)報(bào)值的概率分布；Jiang等［19］利用BMA方法集成多衛(wèi)星降水產(chǎn)品，對(duì)湘江日徑流進(jìn)行模擬；Xu等［20］以中國為研究區(qū)，研究了BMA在多模態(tài)水文預(yù)報(bào)后處理中的應(yīng)用。以上研究均表明BMA能夠提高水文預(yù)測(cè)的準(zhǔn)確性與可靠性，還能得到更為優(yōu)良的預(yù)測(cè)區(qū)間。然而，目前運(yùn)用BMA方法進(jìn)行中長(zhǎng)期日徑流預(yù)測(cè)的相關(guān)研究尚不多見。

疏勒河是河西走廊三大內(nèi)陸河之一，其上游出山徑流量制約著下游玉門、瓜州以及敦煌等地區(qū)的經(jīng)濟(jì)發(fā)展和生態(tài)環(huán)境變化［21］。因此，對(duì)疏勒河上游徑流進(jìn)行準(zhǔn)確地模擬預(yù)測(cè)對(duì)于加強(qiáng)區(qū)域內(nèi)水資源管理、實(shí)現(xiàn)水資源高效利用以及促進(jìn)人與生態(tài)可持續(xù)發(fā)展都具有重要的現(xiàn)實(shí)意義［22］。因此，本文將以疏勒河上游作為研究區(qū)，利用疏勒河上游出山口昌馬堡水文站2010—2017年實(shí)測(cè)日徑流數(shù)據(jù)，基于SVM、ELM、MARS模型對(duì)疏勒河上游流域未來1～7日徑流量進(jìn)行預(yù)測(cè)，運(yùn)用貝葉斯模型平均（BMA）方法對(duì)ELM、SVM、MARS模型的預(yù)測(cè)結(jié)果進(jìn)行組合，構(gòu)建徑流組合預(yù)測(cè)模型，并提供確定性預(yù)測(cè)與概率預(yù)測(cè)結(jié)果。

1 研究方法

1.1 極限學(xué)習(xí)機(jī)（Extreme Learning Machine，ELM）

ELM通過對(duì)輸入權(quán)值和偏置隨機(jī)賦值，并根據(jù)“最小二乘法”的原理，運(yùn)用Moor-Penrose偽逆矩陣計(jì)算輸出權(quán)值，在保證學(xué)習(xí)和訓(xùn)練精度的前提下，克服了傳統(tǒng)神經(jīng)網(wǎng)絡(luò)訓(xùn)練速度慢、容易陷入局部最優(yōu)等缺點(diǎn)。假設(shè)n、l、m分別是網(wǎng)絡(luò)輸入層、隱含層以 及 輸 出 層 的 節(jié) 點(diǎn) 數(shù)，對(duì) 于N個(gè) 樣 本(xi，ti)，i=1，2，…，N，其 中xi=[xi1，xi2，…，xin]T∈Rn，ti=[ti1，ti2，…，tim]T∈Rm，ELM模型可表示為：

其中：g（x）是第l個(gè)隱節(jié)點(diǎn)的閾值和非線性映射函數(shù)；βi是隱含層和輸出層的連接權(quán)值向量；ωi是輸入層和隱含層的連接權(quán)值向量；bi是隱含層神經(jīng)元閾值；oj是網(wǎng)絡(luò)輸出。此時(shí)，上式可以改寫為：

其中：H是隱含層的輸出矩陣；β是輸出權(quán)值矩陣；T是樣本集目標(biāo)矩陣。

實(shí)現(xiàn)基于ELM預(yù)測(cè)徑流量的過程如下。第一步：確定隱含神經(jīng)元個(gè)數(shù)，本研究根據(jù)均方誤差（Mean Squared Error，MSE）最小原則，設(shè)定隱藏神經(jīng)元個(gè)數(shù)。隨機(jī)設(shè)定輸入層與隱含層的連接權(quán)值ωi和隱含神經(jīng)元的閾值bi。第二步：選擇一個(gè)無限可微的函數(shù)作為隱含層神經(jīng)元的激活函數(shù)，進(jìn)而計(jì)算隱含層輸出矩陣H。通過交叉驗(yàn)證方法，選取誤差最小的sigmoid函數(shù)為激活函數(shù)。第三步：計(jì)算輸出層的連接權(quán)值矩陣β∶β=T+H。

1.2 支持向量機(jī)（Support Vector Machine，SVM）

支持向量機(jī)是根據(jù)統(tǒng)計(jì)學(xué)理論提出的一種通用學(xué)習(xí)方法，通過核函數(shù)把低維空間中的非線性回歸問題映射到高維特征空間中，在高維特征空間中求解凸優(yōu)化的問題，能很好地解決非線性、高維數(shù)、小樣本以及局部極小點(diǎn)等問題。主要轉(zhuǎn)換過程如下：

最終回歸函數(shù)為：

式中：?(x)為非線性映射函數(shù)；ω為超平面法向量；C為懲罰參數(shù)；b為超平面偏移量；ε為線性不靈敏損失函數(shù)；K(xi，yi)=?(xi)?(yi)是滿足Mercer條件的核函數(shù)；ai和a*

i為二次規(guī)劃中Lagrange乘子。SVM核函數(shù)的選擇直接影響模型徑流量預(yù)測(cè)結(jié)果好壞，高斯核函數(shù)是一種局部性極強(qiáng)的核函數(shù)，面對(duì)小樣本數(shù)據(jù)可以使用較少的參數(shù)取得較好的分類性能，因此本研究選取高斯徑向基核函數(shù)（radial）作為核函數(shù)。在此基礎(chǔ)上，通過交叉驗(yàn)證網(wǎng)格搜索算法確定最優(yōu)懲罰系數(shù)C與核函數(shù)參數(shù)γ。

1.3 多元自適應(yīng)回歸樣條（Multivariate Adaptive Regression Spline，MARS）

多元自適應(yīng)回歸樣條模型是一種非參數(shù)和非線性回歸方法，其原理是進(jìn)行局部回歸建模，通過樣條函數(shù)估算非線性模型的變化趨勢(shì)，將數(shù)據(jù)整體劃分為幾個(gè)子區(qū)域，在每個(gè)特定的子區(qū)域，響應(yīng)變量用線性回歸進(jìn)行擬合。MARS的優(yōu)點(diǎn)是能夠估算基函數(shù)的貢獻(xiàn)，并通過模擬變量之間的交互影響和加性來確定對(duì)應(yīng)變量。MARS把基函數(shù)作為基本單元，可表示為：

其中：hi(x)為第i個(gè)基函數(shù)，x為輸入變量；c為輸入?yún)?shù)的閾值向量。MARS的結(jié)構(gòu)為：

其中：f(x)為輸出結(jié)果；β0為初始常數(shù)；βi為第i個(gè)基函數(shù)的系數(shù)；m為基函數(shù)的個(gè)數(shù)；hi(x)為第i個(gè)基函數(shù)。

在建模過程中，模型使用常數(shù)值對(duì)目標(biāo)變量進(jìn)行估計(jì)，得到目標(biāo)變量的平均值，根據(jù)最大基函數(shù)的個(gè)數(shù)m，所有可能的基函數(shù)都會(huì)加入到模型中，從而導(dǎo)致結(jié)果出現(xiàn)過擬合現(xiàn)象。為避免這一局限，通過遵循廣義交互驗(yàn)證（Generalized Cross Validation，GCV）原則，識(shí)別重要性較低的基函數(shù)，并去除這些基函數(shù)，獲得最佳MARS模型。

式中：N為數(shù)據(jù)總數(shù)；H為基函數(shù)的個(gè)數(shù)；d為懲罰因子。

實(shí)現(xiàn)基于MARS預(yù)測(cè)徑流量的過程關(guān)鍵主要為懲罰因子d個(gè)數(shù)的選擇，懲罰因子數(shù)量決定了基函數(shù)的個(gè)數(shù)，d越小，生成的MARS模型中基函數(shù)越多。d越大，被放置的節(jié)點(diǎn)數(shù)越少，函數(shù)就會(huì)更加平滑，但同時(shí)被排除在外的基函數(shù)也就越多。本研究通過反復(fù)訓(xùn)練模型，根據(jù)MSE最小原則確定懲罰因子個(gè)數(shù)。

1.4 貝葉斯模型平均（Bayesian Model Averaging，BMA）

BMA方法對(duì)單個(gè)模型進(jìn)行加權(quán)平均，生成一個(gè)整體的概率分布函數(shù)（PDF），通過PDF集成不同模型的預(yù)測(cè)結(jié)果得到的更可靠的綜合預(yù)測(cè)值［23］。該方法的基本原理為：

設(shè)y為徑流模擬值，D=[d1，d2，…，dr]為徑流實(shí)測(cè)值，f=[f1，f2，…，fk]為所取K個(gè)水文模型組成的模型空間。根據(jù)總概率法則，BMA模擬變量y的概率密度函數(shù)可以表示為：

其中：p(fk|D)為模型所模擬序列fk的后驗(yàn)概率，即在給定的數(shù)據(jù)條件下模型為最優(yōu)模型的概率，反映了fk與實(shí)際徑流量的吻合程度。實(shí)際上，p(fk|D)是BMA的權(quán)重ωk，預(yù)測(cè)精度越高的模型分配到的權(quán)重越大，并且本質(zhì)上，BMA均值預(yù)測(cè)序列是對(duì)權(quán)重為ωi的不同模型的最優(yōu)預(yù)測(cè)序列進(jìn)行加權(quán)平均。pk(y|fk，D)是給定模型預(yù)測(cè)fk與數(shù)據(jù)D條件下的預(yù)測(cè)值Q的后驗(yàn)分布。BMA平均預(yù)測(cè)值與方差公式如下：

其中：σ2k為給定觀測(cè)數(shù)據(jù)D和模型fk的條件下模擬變量的方差。BMA預(yù)測(cè)變量方差包括模型間誤差和模型自身誤差，式（11）中

本研究運(yùn)用蒙特卡洛組合抽樣的方法生成BMA在任意時(shí)刻t的預(yù)測(cè)值的不確定性區(qū)間［24］。步驟如下：

（1）通 過 單 一 模 型 權(quán) 重［ω1，ω2，…，ωk］，在［1，2，…，K］中隨機(jī)生成整數(shù)k以抽選模型。具體步驟為：①設(shè)累計(jì)概率=0，計(jì)算；②在0～1之間隨機(jī)生成小數(shù)u；③若滿足，則選擇第k個(gè)模型。

（2）通過第k個(gè)模型在t時(shí)刻的概率分布隨機(jī)生成流量值是均值為f t k，方差為的正態(tài)分布。（3）對(duì)步驟（1）（2）重復(fù)M次，M是任意t時(shí)刻的樣本容量，本文M=10 000。

BMA在任意t時(shí)刻得到M個(gè)樣本后，將樣本從小到大排序，BMA的95%預(yù)測(cè)區(qū)間即為2.5%與97.5%分位數(shù)之間的部分。

BMA均值預(yù)測(cè)過程主要依賴于Zellner的g先驗(yàn)下的線性模型，為減少先驗(yàn)信息的主觀性，本文將參數(shù)先驗(yàn)概率分布設(shè)定為單位信息先驗(yàn)分布（UIP）。參數(shù)先驗(yàn)概率分布設(shè)定完成，還需要設(shè)置模型先驗(yàn)概率分布，本文選取隨機(jī)分布（Random Distribution）為先驗(yàn)概率分布。為減少對(duì)所有模型平均計(jì)算的繁雜過程，本文選取蒙特卡洛抽樣方法抽取后驗(yàn)概率較高的模型。

2 研究區(qū)概況與模型構(gòu)建

2.1 研究區(qū)概況

疏勒河發(fā)源于托勒南山與疏勒南山之間的沙果林那木吉木嶺。疏勒河上游（圖1）主要包含了出山口昌馬堡以上的山區(qū)（96.6°～99.0° E，38.3°～39.9°N），流域面積1.14×104km2，涵蓋了甘肅省酒泉市肅北蒙古族自治州縣、青海省西蒙古族自治州天駿縣等地區(qū)。地形上由托勒南山、疏勒南山以及疏勒河谷組成。河谷地區(qū)相對(duì)平坦，山區(qū)地勢(shì)高大險(xiǎn)峻、地形陡峭，海拔約在2 100～5 750 m，4 500 m以上廣泛分布著現(xiàn)代冰川，冰川資源豐富［25］。

圖1 研究區(qū)位置Fig.1 Location of the study area

疏勒河出山徑流的年內(nèi)分配差異較大，受氣候影響，徑流多集中在5月至9月，占全年徑流量的80%左右，最大日徑流可達(dá)415 m3·s-1，10月到次年4月徑流量則明顯偏少。其中，冬季降水常以固態(tài)形式存在，每年12月到次年3月日徑流量較少，甚至無徑流。4、5月氣溫逐漸回升，積雪和冰川融化使徑流量增加，隨之進(jìn)入一年中降水最為集中的季節(jié)，加上高山冰雪融水的增加，7、8月份出山徑流量達(dá)到最大值，10月之后隨著氣溫的降低和暖濕氣流影響的削弱，徑流量逐漸減少［26］。

2.2 數(shù)據(jù)獲取與處理

本文選用昌馬堡水文站2010—2017年日徑流量作為疏勒河出山徑流分析數(shù)據(jù)。其中，2010—2014年數(shù)據(jù)用于模型訓(xùn)練，2015—2017年數(shù)據(jù)用于模型測(cè)試。疏勒河上游出山口徑流量數(shù)據(jù)基本統(tǒng)計(jì)資料見表1。由表1可知，徑流量數(shù)據(jù)呈現(xiàn)明顯的高偏態(tài)分布（訓(xùn)練數(shù)據(jù)集為2.39，測(cè)試數(shù)據(jù)集為2.81），并且不同數(shù)據(jù)集日徑流數(shù)據(jù)的變異系數(shù)均大于1，這可能會(huì)對(duì)預(yù)測(cè)精度產(chǎn)生一定的影響。為了保證對(duì)模型的有效訓(xùn)練以及有效特征的提取，通常對(duì)輸入數(shù)據(jù)進(jìn)行歸一化處理。本文采用最小-最大歸一化方法對(duì)徑流數(shù)據(jù)時(shí)間序列進(jìn)行歸一化處理。公式為：

表1 昌馬堡日徑流量基本統(tǒng)計(jì)值Table 1 Basic daily runoff statistics of Changmapu

式中：Q*為經(jīng)標(biāo)準(zhǔn)化處理后的徑流序列；Q為原始徑流序列；Qmin為原始最小徑流量；Qmax為原始最大徑流量。

2.3 模型構(gòu)建

徑流時(shí)間序列能夠揭示其隨時(shí)間變化的規(guī)律，并將這種規(guī)律延伸到未來。因此，徑流時(shí)間序列可以作為徑流預(yù)測(cè)的基本輸入因素。本文通過計(jì)算2010—2017年實(shí)測(cè)日徑流的自相關(guān)函數(shù)（ACF）和偏自相關(guān)函數(shù)（PACF）來確定模型的輸入。如圖2所示，徑流時(shí)間序列的ACF與PACF均表現(xiàn)出拖尾特征，ACF逐漸衰減并在95%的置信區(qū)間內(nèi)表現(xiàn)出明顯的自相關(guān)性，特別是在時(shí)滯數(shù)小于5時(shí)，ACF大于0.8，表現(xiàn)出顯著的自相關(guān)性。同時(shí)，PACF顯著不為0的時(shí)滯數(shù)為1，2，3。因此，綜合ACF與PACF分析結(jié)果，確定Qt-3、Qt-2、Qt-1為模型輸入量。

圖2 實(shí)測(cè)徑流時(shí)間序列的ACF和PACFFig.2 ACF and PACF of measured runoff time series

構(gòu)建一個(gè)包含7列數(shù)據(jù)的表格，表格形式為（Qt-3，Qt-2，Qt-1，Qt+1，Qt+3，Qt+5，Qt+7），其中，Qt-3、Qt-2、Qt-1為輸入變量，分別表示t-3、t-2、t-1日的徑流量組合；Qt+1、Qt+3、Qt+5、Qt+7表示模型在t+1、t+3、t+5、t+7日的預(yù)測(cè)徑流量，即模型的輸出。

對(duì)于單一模型的構(gòu)建主要分為兩部分：第一部分為確定模型的輸入向量；第二部分為利用確定的輸入向量設(shè)置模型運(yùn)行的最優(yōu)基本參數(shù)，建立中長(zhǎng)期日徑流預(yù)測(cè)模型。對(duì)于BMA組合模型的構(gòu)建，主要是利用各模型后驗(yàn)概率將單一模型的預(yù)測(cè)結(jié)果組合起來得到集合預(yù)測(cè)結(jié)果，建立中長(zhǎng)期日徑流組合預(yù)測(cè)模型。此外，本文采用馬爾科夫蒙特卡洛抽樣對(duì)集合預(yù)測(cè)結(jié)果進(jìn)行10 000次抽樣，獲取集合預(yù)測(cè)結(jié)果的95%置信區(qū)間，對(duì)預(yù)測(cè)結(jié)果進(jìn)行不確定性分析。為驗(yàn)證BMA集合預(yù)測(cè)的有效性，在相同預(yù)見時(shí)間的條件下，與單一ELM、SVM、MARS模型進(jìn)行對(duì)比研究。本文主要研究技術(shù)路線如圖3所示。

圖3 研究技術(shù)路線Fig.3 Research technical route

2.4 模型評(píng)價(jià)指標(biāo)

模擬效率可體現(xiàn)模型在研究區(qū)的適用情況，本文選取相關(guān)系數(shù)R、納什效率系數(shù)NSE（Nash-Sutcliffe efficient）以及均方根誤差RMSE（Root-Mean-Square Error）對(duì)預(yù)測(cè)結(jié)果進(jìn)行評(píng)價(jià)，表達(dá)式如下：

式中：Qobs為實(shí)測(cè)值（m3·s-1）；Qf為預(yù)測(cè)值（m3·s-1）；為實(shí)測(cè)均值（m3·s-1）；為預(yù)測(cè)值均值（m3·s-1）；n為實(shí)測(cè)數(shù)據(jù)數(shù)量。

R是判斷實(shí)測(cè)值與預(yù)測(cè)值的線性相關(guān)程度。R取值范圍是-1至1，R的絕對(duì)值越接近于1，線性相關(guān)程度越高；NSE代表預(yù)測(cè)值與實(shí)測(cè)值在1∶1水平上的相似度，NSE在-∞至1.0之間取值，NSE≥0.5時(shí)預(yù)測(cè)結(jié)果可接受，NSE越接近1，擬合效果越好；RMSE反映預(yù)測(cè)值與實(shí)測(cè)值之間的偏差，該值越接近于0，模型的擬合度越高。當(dāng)R=1、NSE=1、RMSE=0時(shí)，認(rèn)為該模型為最佳模型［27］。

評(píng)價(jià)預(yù)測(cè)結(jié)果的不確定性區(qū)間的優(yōu)良性可以體現(xiàn)模型的可靠性。本文采用覆蓋率（CR）、平均帶寬（B）以及平均偏移幅度（D）三個(gè)主要指標(biāo)來分析比較BMA的不確定性區(qū)間［28］。

覆蓋率是評(píng)價(jià)模擬區(qū)間最常用的指標(biāo)，指預(yù)測(cè)區(qū)間所覆蓋的實(shí)測(cè)流量數(shù)據(jù)的比率。CR值越大，表明模擬區(qū)間的覆蓋率越高，預(yù)測(cè)結(jié)果包含的真實(shí)信息越多；平均帶寬B是在指定的置信水平，以及較高的覆蓋率保證的前提下，預(yù)測(cè)區(qū)間的平均帶寬越窄預(yù)測(cè)區(qū)間的不確定程度越小，反之，預(yù)測(cè)的不確定性程度越大；平均偏移幅度D是評(píng)價(jià)預(yù)測(cè)區(qū)間的中心線相較實(shí)際徑流過程線偏離程度的指標(biāo)，理論上，平均偏移幅度越小，表明預(yù)測(cè)區(qū)間的對(duì)稱性越好。

3 結(jié)果分析

3.1 單一模型預(yù)測(cè)精度分析

為了驗(yàn)證單一模型預(yù)測(cè)的有效性，表2～4給出了最優(yōu)參數(shù)下三個(gè)單一模型在不同時(shí)間的預(yù)測(cè)結(jié)果評(píng)價(jià)。由表2～4可知，隨著時(shí)間的增加，三個(gè)單一模型的預(yù)測(cè)精度逐漸降低。但總體上三個(gè)單一模型的預(yù)測(cè)值與實(shí)測(cè)值的R值均在0.75以上，表明預(yù)測(cè)值與實(shí)測(cè)值之間存在著很高的線性相關(guān)關(guān)系。NSE值均在0.55以上，說明預(yù)測(cè)值與實(shí)測(cè)值之間的擬合程度是可以接受的。因此，說明所構(gòu)建的單一模型的預(yù)測(cè)結(jié)果合理、有效，可以作為BMA集合預(yù)測(cè)的有效成員。

表2 ELM模型參數(shù)及預(yù)測(cè)結(jié)果評(píng)價(jià)Table 2 Parameter and evaluation of prediction results of ELM model

3.2 BMA均值預(yù)測(cè)與單個(gè)模型的預(yù)測(cè)值對(duì)比

本研究得到的BMA方法在測(cè)試期的（t+1）d、（t+3）d、（t+5）d與（t+7）d的預(yù)測(cè)結(jié)果如圖4與表5所示。通過水文過程曲線可以看出（圖4），BMA的整體預(yù)測(cè)效果比較好，預(yù)測(cè)值與實(shí)測(cè)值變化趨勢(shì)一致，洪峰出現(xiàn)時(shí)間同實(shí)際觀測(cè)情況一致。BMA對(duì)徑流低值的預(yù)測(cè)結(jié)果較好，但對(duì)徑流高值的預(yù)測(cè)有一定的誤差，并且隨著預(yù)測(cè)時(shí)間的增加，峰值預(yù)測(cè)誤差逐漸增大。這被認(rèn)為是目前徑流研究工作中的一個(gè)重要局限性，在其他徑流模擬預(yù)測(cè)研究中也出現(xiàn)了同樣的局限性［9-10］。盡管如此，從散點(diǎn)圖可以看出，BMA在（t+1）～（t+7）d上的中長(zhǎng)期日徑流的R均大于0.8，表明預(yù)測(cè)精度是可以接受的［29］。另外，由表5可知，BMA的預(yù)測(cè)精度隨著預(yù)測(cè)時(shí)間的增加而逐漸減小，但BMA在（t+1）～（t+7）d的R值均大于0.78，表明BMA的預(yù)測(cè)值與實(shí)測(cè)值之間存在較高的線性相關(guān)關(guān)系；NSE值均大于0.6，表明預(yù)測(cè)值與實(shí)測(cè)值之間的擬合效果較好。綜上，說明BMA可以用于（t+1）～（t+7）d的中長(zhǎng)期日徑流預(yù)測(cè)，并具有較高的預(yù)測(cè)精度。然而，BMA在（t+1）d后訓(xùn)練期的R值不如測(cè)試期的R值，而訓(xùn)練期的RMSE值卻優(yōu)于訓(xùn)練期。由表1可知，訓(xùn)練期的變差系數(shù)高于測(cè)試期，訓(xùn)練期數(shù)據(jù)離散程度大，導(dǎo)致訓(xùn)練期R值低于測(cè)試期；而測(cè)試期數(shù)據(jù)的偏度大于訓(xùn)練期，分布不均勻，產(chǎn)生的較大誤差。

表3 SVM模型參數(shù)及預(yù)測(cè)結(jié)果評(píng)價(jià)Table 3 Parameter and evaluation of prediction results of SVM model

表4 MARS模型參數(shù)及預(yù)測(cè)結(jié)果評(píng)價(jià)Table 4 Parameter and evaluation of prediction results of MARS model

表5 BMA在（t+1）～（t+7）d日徑流量預(yù)測(cè)結(jié)果評(píng)價(jià)Table 5 Evaluation of daily runoff predicted results of BMA in（t+1）～（t+7）days

由表2～5可知，BMA與組成BMA的3個(gè)單一模型均可適用于（t+1）～（t+7）d的中長(zhǎng)期日徑流預(yù)測(cè)，但進(jìn)一步對(duì)比分析可知，BMA的預(yù)測(cè)結(jié)果較單一模型具有更高的預(yù)測(cè)精度。以訓(xùn)練期和測(cè)試期的（t+1）d與（t+7）d預(yù)測(cè)結(jié)果為例，訓(xùn)練期（t+1）d時(shí)BMA的R值均為0.979，高于SVM（0.971）ELM（0.959）和MARS（0.972）的R值。BMA的NSE值為0.982，高于單一ELM（0.919）、SVM（0.919）、MARS（0.944）的NSE值。BMA的RMSE值較單一ELM、SVM、MARS的RMSE值 分 別 提 高 了28.35%、26.08%、13.77%。在訓(xùn)練期（t+7）d時(shí)，BMA與MARS的R值（0.786）相同，高于單一ELM（0.781）和SVM模型（0.778）。BMA的NSE值為0.618，分別高于三個(gè)單一模型的NSE值（ELM為0.610、SVM為0.595、MARS為0.612）BMA的RMSE值較單一ELM、SVM、MARS模型分別提高了4.79%、7.17%、4.58%。

在測(cè)試期（t+1）d時(shí)，BMA與SVM的R值相同（0.967），較ELM、MARS的R值 分 別 提 高 了1.83%、0.6%；NSE值較單一ELM、SVM、MARS模型分別提高了3.7%、3.1%、0.3%；RMSE值較單一ELM、SVM、MARS模 型 分 別 降 低 了33.7%、25.8%、3.8%。在測(cè)試期（t+7）d時(shí)，ELM、SVM、MARS的R值分別為0.781、0.783、0.781，BMA的R值 為0.800，高 于 單一模 型。BMA的NSE值為0.64，較單一模型ELM、SVM、MARS的NSE值分別提高了5.8%、9.1%、5.6%。另外，ELM、SVM、MARS的RMSE值分別為35.548、36.494、35.503，均大 于BMA的RMSE值33.878。因此，綜合R、NSE、RMSE的評(píng)價(jià)結(jié)果可以看出，BMA較單一模型無論是在（t+1）～（t+3）d的短期日徑流，還是在（t+5）～（t+7）d的中長(zhǎng)期日徑流預(yù)測(cè)均具有較高的預(yù)測(cè)精度，能夠提供更準(zhǔn)確的預(yù)測(cè)結(jié)果。

對(duì)不同等級(jí)的徑流量進(jìn)行預(yù)測(cè)對(duì)及時(shí)發(fā)布極端水文事件預(yù)警信息以及制定防治預(yù)案具有重要意義。本文參考董磊華等［18］研究成果，根據(jù)疏勒河的徑流量特征，將徑流量分為3個(gè)等級(jí)：流量從大到小排序，將前10%的大流量定為高水，中間50%的流量定為中水，后40%的小流量定為低水。對(duì)比分析BMA與單一模型ELM、SVM、MARS在不同等級(jí)徑流量的預(yù)測(cè)結(jié)果（表6）。從表6可以看出，BMA與單一模型對(duì)中水預(yù)測(cè)效果最好，BMA與單一模型均可預(yù)測(cè)出未來5～7 d中水徑流量；低水次之，除單一SVM模型可以預(yù)測(cè)未來3日的中水徑流量，其余模型均只能預(yù)測(cè)未來1日低水徑流量；高水最差，僅BMA與SVM模型能預(yù)測(cè)未來1日的高水徑流量。但相比之下，BMA較單一模型的預(yù)測(cè)效果更好，這是因?yàn)锽MA方法綜合發(fā)揮了單一模型各自的優(yōu)勢(shì)，在一定程度上克服了單一模型無法避免的局部極值問題［28］。

表6 BMA和組成它的3個(gè)模型的預(yù)測(cè)值在3個(gè)不同流量等級(jí)的統(tǒng)計(jì)結(jié)果Table 6 Statistical results of BMA and the forecast values of the three models that comprised it at three different flow levels

續(xù)表6

3.3 BMA徑流預(yù)測(cè)不確定性分析

對(duì)預(yù)測(cè)結(jié)果的不確定性分析可以判斷模型結(jié)構(gòu)的可靠性，對(duì)提高水資源管理的科學(xué)性、極端水文事件風(fēng)險(xiǎn)防范能力具有重要意義［29］。本文基于覆蓋度（CR）、平均區(qū)間寬度（B）以及平均偏移幅度（D）3個(gè)指標(biāo)對(duì)BMA的95%置信區(qū)間的徑流預(yù)測(cè)結(jié)果的不確定性進(jìn)行分析。從圖4可看出，大多數(shù)實(shí)測(cè)值都在不確定性區(qū)間之內(nèi)，測(cè)試期內(nèi)BMA在（t+1）～（t+7）d的95%置信區(qū)間的覆蓋度分別為達(dá)到了94.7%、94.1%、93.8%、92.9%，幾乎覆蓋了整個(gè)實(shí)測(cè)徑流序列，表明95%置信區(qū)間預(yù)測(cè)效果好，不確定性較小。

由于BMA對(duì)高水和低水的預(yù)測(cè)出現(xiàn)了較大偏差，存在較大的不確定性。因此，為進(jìn)一步研究模型結(jié)構(gòu)的不確定性，對(duì)每個(gè)等級(jí)徑流預(yù)測(cè)結(jié)果進(jìn)行評(píng)價(jià)。表7綜合了BMA在3個(gè)不同流量等級(jí)的預(yù)測(cè)區(qū)間優(yōu)良性，可以看出BMA預(yù)測(cè)區(qū)間覆蓋了整個(gè)低水部分，中水部分除在（t+7）d時(shí)刻上覆蓋率較低，其余時(shí)間均在95%以上。但是測(cè)試期與訓(xùn)練期在高水部分的覆蓋率均不足65%，其中測(cè)試期（t+7）d的高水部分最低覆蓋率僅43.52%，拉低了整體覆蓋率。同時(shí)，BMA 95%置信區(qū)間的平均帶寬在低水部分最小，在高水部分最大，表明BMA在高水部分的預(yù)測(cè)不確定性更大。就平均偏移幅度來說，BMA在高水部分的平均偏移幅度最大，中水部分最小，說明BMA在高水部分的預(yù)測(cè)值偏離實(shí)測(cè)值的程度最大，中水部分的預(yù)測(cè)值偏離實(shí)測(cè)值的程度最小。綜上，BMA的預(yù)測(cè)不確定性主要來源于高水部分，主要原因在于高水一般出現(xiàn)在汛期，徑流量受降水、融水影響較大，預(yù)測(cè)難度較大，因此高水預(yù)測(cè)的不確定性較大。另外，各流量等級(jí)的區(qū)間不確定性均隨著時(shí)間增加而增大。但總體而言，BMA的95%置信區(qū)間效果良好，對(duì)實(shí)測(cè)徑流量覆蓋率高。

表7 BMA 95%在3個(gè)不同流量等級(jí)的統(tǒng)計(jì)結(jié)果Table 7 Statistical results of BMA 95%confidence interval at three different traffic levels

4 討論

李洪源等［30］在運(yùn)用分布式水文模型SPHY對(duì)疏勒河上游徑流進(jìn)行了逐日與逐月徑流預(yù)測(cè)，并取得了較好的預(yù)測(cè)精度。但SPHY在率定期的日徑流NSE值為0.62，低于本文BMA在訓(xùn)練期（t+1）d的NSE值0.943。SPHY在測(cè)試期中的日徑流的NSE值為0.79，低于本文BMA在測(cè)試期（t+1）d的NSE值（0.935），表明BMA方法具有較高的精度，這可能是因?yàn)锽MA方法通過對(duì)各單一模型的后驗(yàn)概率進(jìn)行加權(quán)平均，有效解決了單一模型的不確定性，把多個(gè)單一模型的優(yōu)點(diǎn)組合起來，從而獲得了更準(zhǔn)確的預(yù)測(cè)值。

本文通過計(jì)算徑流的ACF與PACF確定模型的輸入，但如何利用徑流本身的自相關(guān)關(guān)系確定模型輸入暫無統(tǒng)一標(biāo)準(zhǔn)。對(duì)于自相關(guān)關(guān)系差的小流域來說，可以針對(duì)流域自身特點(diǎn)來確定模型的輸入。如于海姣等［31］在對(duì)小流域排露溝日降水-徑流模擬研究中，針對(duì)降水-徑流的非線性關(guān)系，選擇相對(duì)百分比誤差方法來確定模型的輸入。Gavin等［32］根據(jù)輸入與輸出之間的顯著線性關(guān)系，運(yùn)用互相關(guān)法確定模型輸入，對(duì)澳大利亞Murray河流的含鹽量進(jìn)行模擬。

BMA方法通過集成各個(gè)模型的自身優(yōu)勢(shì)，彌補(bǔ)單個(gè)模型的不足獲得較單一模型更高的預(yù)測(cè)精度。在疏勒河上游流域（t+1）～（t+7）d的中長(zhǎng)期日徑流預(yù)測(cè)中，BMA方法較單一模型的預(yù)測(cè)效果更好，精度更高，為缺乏資料的干旱半干旱地區(qū)的徑流預(yù)測(cè)提供了有效的預(yù)測(cè)工具。但是，BMA在中長(zhǎng)期日徑流預(yù)測(cè)中仍然存在著一定的不足。BMA在高徑流量和低徑流量預(yù)測(cè)上表現(xiàn)出較單一模型更高的精度，但是預(yù)測(cè)值仍然小于實(shí)測(cè)值，存在著較大的誤差。誤差原因可能是：疏勒河上游流域冰川積雪、凍土分布廣，由此形成了獨(dú)特的寒區(qū)水文過程，流域內(nèi)的溫度、蒸發(fā)、凍融作用都是影響高徑流量時(shí)期徑流變化的重要因素，所以，只將徑流數(shù)據(jù)作為模型輸入過于單一，導(dǎo)致預(yù)測(cè)序列出現(xiàn)偏差［31，33］。低徑流時(shí)期的徑流量數(shù)值較小且存在0值，經(jīng)歸一化處理后的徑流數(shù)據(jù)產(chǎn)生較多的0值，導(dǎo)致模型在運(yùn)行過程中出現(xiàn)一定的誤差［31］。因此，在日后的工作中可以嘗試從增加模型輸入、優(yōu)化模型參數(shù)計(jì)算等方面提高模型預(yù)測(cè)精度，更真實(shí)地反映徑流變化規(guī)律。

5 結(jié)論

本文利用BMA方法對(duì)疏勒河上游2010—2017年的（t+1）～（t+7）d中長(zhǎng)期日徑流進(jìn)行了預(yù)測(cè)，對(duì)BMA與組成BMA的ELM、SVM、MARS模型的預(yù)測(cè)結(jié)果進(jìn)行分析對(duì)比，并對(duì)BMA預(yù)測(cè)結(jié)果的不確定性進(jìn)行定量分析。得到以下結(jié)論：

（1）BMA在（t+1）～（t+7）d的徑流預(yù)測(cè)結(jié)果比單一ELM、SVM、MARS結(jié)果準(zhǔn)確度更高，不同時(shí)期徑流量的預(yù)測(cè)值與實(shí)測(cè)值的相關(guān)系數(shù)均在0.7以上，納什系數(shù)均在0.55以上，能夠再現(xiàn)徑流量的變化趨勢(shì)，說明BMA方法具有較好的實(shí)用性，在中長(zhǎng)期日徑流預(yù)測(cè)中可以提供更準(zhǔn)確的徑流量預(yù)測(cè)值。BMA對(duì)高水、低水的預(yù)測(cè)值與實(shí)測(cè)值的相關(guān)系數(shù)、納什系數(shù)在（t+1）d后均低于0.5，但仍較單一模型保持較高精度。

（2）貝葉斯模型加權(quán)平均（BMA）方法是一種通過綜合多個(gè)模型預(yù)測(cè)結(jié)果的后驗(yàn)分布來推斷預(yù)測(cè)結(jié)果的可靠概率分布分析的工具。BMA對(duì)（t+1）～（t+7）d的95%置信區(qū)間實(shí)現(xiàn)了對(duì)低水徑流量的100%全覆蓋，對(duì)中水徑流量的覆蓋率在70%以上，對(duì)高水徑流量的覆蓋率較低，但整體覆蓋率均在90%以上。表明BMA不僅可以提供更準(zhǔn)確的綜合預(yù)測(cè)結(jié)果，還能提供一個(gè)綜合預(yù)測(cè)區(qū)間來進(jìn)行模型結(jié)果不確定性分析，可成為資料有限條件下的徑流模擬預(yù)測(cè)的有效工具。