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

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

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

0 引言

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

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

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

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

1 研究方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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ū)相對平坦，山區(qū)地勢高大險峻、地形陡峭，海拔約在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ù)用于模型測試。疏勒河上游出山口徑流量數(shù)據(jù)基本統(tǒng)計資料見表1。由表1可知，徑流量數(shù)據(jù)呈現(xiàn)明顯的高偏態(tài)分布（訓(xùn)練數(shù)據(jù)集為2.39，測試數(shù)據(jù)集為2.81），并且不同數(shù)據(jù)集日徑流數(shù)據(jù)的變異系數(shù)均大于1，這可能會對預(yù)測精度產(chǎn)生一定的影響。為了保證對模型的有效訓(xùn)練以及有效特征的提取，通常對輸入數(shù)據(jù)進(jìn)行歸一化處理。本文采用最小-最大歸一化方法對徑流數(shù)據(jù)時間序列進(jìn)行歸一化處理。公式為：

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

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

2.3 模型構(gòu)建

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

圖2 實測徑流時間序列的ACF和PACFFig.2 ACF and PACF of measured runoff time series

構(gòu)建一個包含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ù)測徑流量，即模型的輸出。

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

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

2.4 模型評價指標(biāo)

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

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

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

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

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

3 結(jié)果分析

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

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

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

3.2 BMA均值預(yù)測與單個模型的預(yù)測值對比

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

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

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

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

由表2～5可知，BMA與組成BMA的3個單一模型均可適用于（t+1）～（t+7）d的中長期日徑流預(yù)測，但進(jìn)一步對比分析可知，BMA的預(yù)測結(jié)果較單一模型具有更高的預(yù)測精度。以訓(xùn)練期和測試期的（t+1）d與（t+7）d預(yù)測結(jié)果為例，訓(xùn)練期（t+1）d時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時，BMA與MARS的R值（0.786）相同，高于單一ELM（0.781）和SVM模型（0.778）。BMA的NSE值為0.618，分別高于三個單一模型的NSE值（ELM為0.610、SVM為0.595、MARS為0.612）BMA的RMSE值較單一ELM、SVM、MARS模型分別提高了4.79%、7.17%、4.58%。

在測試期（t+1）d時，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%。在測試期（t+7）d時，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的評價結(jié)果可以看出，BMA較單一模型無論是在（t+1）～（t+3）d的短期日徑流，還是在（t+5）～（t+7）d的中長期日徑流預(yù)測均具有較高的預(yù)測精度，能夠提供更準(zhǔn)確的預(yù)測結(jié)果。

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

表6 BMA和組成它的3個模型的預(yù)測值在3個不同流量等級的統(tǒng)計結(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ù)測不確定性分析

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

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

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

4 討論

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

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

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

5 結(jié)論

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

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

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