邢貞相 芮孝芳 馮 杰

(東北農(nóng)業(yè)大學(xué)水利與建筑學(xué)院1) 哈爾濱 150030)

(東北農(nóng)業(yè)大學(xué)農(nóng)業(yè)工程博士后科研流動(dòng)站2) 哈爾濱 150030)

(河海大學(xué)水文水資源學(xué)院3) 南京 210098)

BFS(貝葉斯概率預(yù)報(bào)系統(tǒng))是一個(gè)可與任一確定性水文模型協(xié)作進(jìn)行概率水文預(yù)報(bào)的通用理論框架[1],其理論基礎(chǔ)是貝葉斯公式.Krzysztofowicz相繼提出了線性正態(tài)假設(shè)、亞高斯轉(zhuǎn)換、概率定量降雨預(yù)報(bào)和概率河流水位預(yù)報(bào)的貝葉斯系統(tǒng)等[2],推動(dòng)了BFS的研究進(jìn)展.國(guó)內(nèi),張洪剛采用平穩(wěn)序列線性AR模型與線性擾動(dòng)模型(LPM)分別描述先驗(yàn)分布與似然函數(shù),在一定程序上降低了貝葉斯求解的復(fù)雜度[3];李向陽(yáng)等采用神經(jīng)網(wǎng)絡(luò)模型來(lái)描述先驗(yàn)分布與似然函數(shù),進(jìn)一步降低了貝葉斯求解過(guò)程的復(fù)雜程度[4].王建平等[5]將貝葉斯理論用于水質(zhì)模型的參數(shù)識(shí)別問(wèn)題,對(duì)復(fù)雜環(huán)境模型參數(shù)的不確定性進(jìn)行了研究.貝葉斯理論還是貝葉斯網(wǎng)絡(luò)模型的基礎(chǔ),它是一種不確定性知識(shí)的表達(dá)與推理模型,在建筑、經(jīng)濟(jì),環(huán)境等領(lǐng)域均有廣泛應(yīng)用[6].本文嘗試將貝葉斯理論與自適應(yīng)馬爾可夫鏈蒙特卡羅算法相結(jié)合來(lái)研究Nash模型參數(shù)的不確定性,并將其用于洪水概率預(yù)報(bào).

1 BFS的基本原理

BFS的理論依據(jù)就是下列貝葉斯公式

式中:π(θ|x)為參數(shù)的后驗(yàn)密度,它是在樣本 x給定條件下,參數(shù) θ的條件分布;π(θ)為θ的先驗(yàn)分布;p(x|θ)為似然函數(shù);Θ為θ的積分區(qū)間.

π(θ|x)集中了總體、樣本和先驗(yàn)等3種信息中有關(guān)θ的信息,是排除一切與θ無(wú)關(guān)的信息后所得的結(jié)果.基于后驗(yàn)分布 π(θ|x)對(duì) θ進(jìn)行統(tǒng)計(jì)推斷將更為有效、合理,稱之為貝葉斯統(tǒng)計(jì)推斷.

當(dāng)參數(shù)的先驗(yàn)密度與似然函數(shù)形式確定后,為獲得式(1)的后驗(yàn)密度解析式還需求得其右端分母的積分,而參數(shù)θ的積分區(qū)間只能靠實(shí)測(cè)資料估計(jì),無(wú)法獲得其真實(shí)的區(qū)間,所以,很難求得式(1)的解析式,為此,本文采用數(shù)值解法來(lái)獲得后驗(yàn)密度,即用馬爾可夫鏈蒙特卡羅隨機(jī)模擬的方法求其數(shù)值解.

作為隨機(jī)模擬方法的馬爾可夫鏈蒙特卡羅(MCMC)方法的關(guān)鍵是如何選擇推薦分布(轉(zhuǎn)移密度)使采樣更加有效.常用的采樣方法有Metropolis-Hastings算法、吉布斯(Gibbs)采樣和Adapative-Metropolis(AM)算法[7].這 3種方法中只有AM算法的不依賴于事先確定的推薦分布且可并行運(yùn)算,收斂速度快,故本文采用此算法.關(guān)于AM算法的具體過(guò)程和收斂判斷準(zhǔn)則及性能測(cè)試參見文獻(xiàn)[8-9].

2 Nash模型參數(shù)的不確定性

本文利用AM-MCMC算法將Nash模型參數(shù)k,n分兩種情況研究其不確定性:(1)將參數(shù)k視為隨機(jī)的,而參數(shù)n視為確定性的;(2)將2參數(shù)均視為隨機(jī)的.由于第一種情況是第二種情況的特例,故本文只介紹第二種情況的具體過(guò)程.Nash模型的輸入,即地面凈雨的計(jì)算采用斜線分割法.

選取長(zhǎng)江三峽沿渡河流域作為研究區(qū)域,共有30場(chǎng)洪水實(shí)測(cè)資料.該流域位于長(zhǎng)江三峽地區(qū),其水系流經(jīng)神農(nóng)架林區(qū)巴東縣.流域內(nèi)降水豐沛,流域多年平降雨量為1 337 mm,全年雨量以5～9月最多,約占全年68%.流域內(nèi)最大年降水量為2 448.2 mm,最小年降雨量為808.4 mm.

沿渡河流域面積601 km2,流域坡度較大,平均坡度為0.287%,高程垂直落差達(dá)2 800 m,山高坡陡,人類活動(dòng)影響較小,流域內(nèi)耕地面積占流域面積的10%左右,森林覆蓋率在70%以上.由于流域內(nèi)植被覆蓋良好,地表徑流中含沙量不大,除洪水期含沙量有所增大外,其余時(shí)間河水清澈.

2.1 先驗(yàn)分布的確定

1)參數(shù)n的先驗(yàn)分布的確定 假定其服從正態(tài)分布,根據(jù)地貌學(xué)的方法求得沿渡河流域的Nash模型參數(shù)n=3,可認(rèn)為是其分布的均值,設(shè)其先驗(yàn)方差為均值的10%,則得n的先驗(yàn)分布為n～N(3,0.3).

2)參數(shù)k先驗(yàn)分布的確定 首先率定參數(shù) ,選用該流域1981～1987年間28場(chǎng)洪水資料來(lái)率定參數(shù)k值.為保證計(jì)算精度,取計(jì)算時(shí)段長(zhǎng)為1 h.為避免異參同效現(xiàn)象的影響,令n=3保持不變,單獨(dú)率定k,方法采用矩法-優(yōu)選法.根據(jù)各場(chǎng)洪水k的率定結(jié)果求得k取值范圍為[1,1.96],其均值為1.19,方差為0.09,并假定其服從正態(tài)分布,即得k的先驗(yàn)分布為k～N(1.19,0.09).

2.2 似然函數(shù)的確定

由于有28場(chǎng)洪水資料參與率定,故采用多觀測(cè)擬合優(yōu)度的似然函數(shù)

式中:Q為流量;σei為第i個(gè)觀測(cè)與模型預(yù)報(bào)的誤差系列的標(biāo)準(zhǔn)差;N為實(shí)測(cè)序列的個(gè)數(shù),本文為N=28;其余符號(hào)意義同前.

2.3 AM-MCMC算法初始條件

AM-MCMC算法初始條件:初始協(xié)方差取為對(duì)角陣,初始化迭代次數(shù)為2 000,初始化階段次數(shù)為2 000,每次采樣為10 000次,算法并行運(yùn)行5次,這樣共將取樣(10 000-2 000)×5=40 000組(n,k)以用于沿渡河流域Nash模型參數(shù)的不確定性研究.

2.4 Nash模型參數(shù)不確定性分析

AM-MCMC運(yùn)算結(jié)束后,根據(jù)所抽40 000個(gè)樣本,統(tǒng)計(jì)得出Nash模型參數(shù)的邊緣后驗(yàn)密度分別為 k～N(2.03,0.09)(見圖 1),n～N(2.61,0.10)(見圖2),從圖1,圖2可見k和n的后驗(yàn)邊緣密度均近似服從正態(tài)分布,通過(guò)Kolmogorov-Smirnov假設(shè)檢驗(yàn)在顯著性水平為0.05時(shí)接受各自后驗(yàn)分布為正態(tài)分布的原假設(shè).圖3給出了2參數(shù)的后驗(yàn)均值迭代跡線、圖4給出了2參數(shù)的后驗(yàn)方差迭代跡線,從2圖中看出自第2 000次迭代后2參數(shù)的后驗(yàn)均值、后驗(yàn)方差均趨于穩(wěn)定,說(shuō)明所抽樣本已具有總體樣本的統(tǒng)計(jì)特性.圖5給出了兩參數(shù)的聯(lián)合后驗(yàn)概率密度,從圖中看出兩參數(shù)的聯(lián)合分布只有一個(gè)極值,其坐標(biāo)為兩參數(shù)的后驗(yàn)均值.圖6給出了兩參數(shù)樣本的散點(diǎn)圖,由圖可見,n與k之間存在著明顯的相關(guān)關(guān)系.

圖1 參數(shù)k的后驗(yàn)邊緣密度

圖2 參數(shù)n的后驗(yàn)邊緣密度

圖3 參數(shù)n,k的后驗(yàn)均值迭代跡線

圖6 參數(shù)n與k的散點(diǎn)圖

3 Nash模型的降雨徑流概率預(yù)報(bào)

芮孝芳[10]指出,產(chǎn)生水文模型的“異參同效”這一現(xiàn)象的原因至少有:目標(biāo)函數(shù)是多極值的;模型中包含的參數(shù)之間存在相互補(bǔ)償作用;模型參數(shù)具有隨機(jī)性.圖1和圖2雖給出了Nash模型兩參數(shù)的各自后驗(yàn)邊緣密度,但卻無(wú)法避免存在的“異參同效”現(xiàn)象,在實(shí)際水文預(yù)報(bào)時(shí),真正有意義的是兩個(gè)參數(shù)的組合,而不是單個(gè)參數(shù).為此,本文隨機(jī)選取AM-MCMC算法收斂后的10 000個(gè)參數(shù)組樣本分別對(duì)沿渡流域洪水進(jìn)行模擬,使某一場(chǎng)洪水的每個(gè)時(shí)段對(duì)應(yīng)所選取的不同參數(shù)組生成10 000個(gè)流量數(shù)值.用這些數(shù)據(jù)作為樣本來(lái)研究各時(shí)刻流量的統(tǒng)計(jì)特性,即可求得各時(shí)刻(包括洪峰時(shí)刻)流量的概率分布,其均值和方差及指定概率的置信區(qū)間.在作業(yè)預(yù)報(bào)時(shí)可采用每一時(shí)刻的預(yù)報(bào)流量樣本的均值作為其預(yù)報(bào)值.

表1中只給出了本文算法對(duì)沿渡河流域6場(chǎng)洪水(其他洪水限于篇幅未列出)的峰值概率預(yù)報(bào)及其80%的置信區(qū)間.在表1中,同時(shí)給出了當(dāng)參數(shù)k為隨機(jī)而n為確定時(shí)的相應(yīng)場(chǎng)次洪水的峰值預(yù)報(bào)結(jié)果(研究方案1為僅參數(shù)k為隨機(jī)的情況,方案2為參數(shù)k和n均為隨機(jī)的情況).通過(guò)對(duì)該流域30(其中的28場(chǎng)為參數(shù)率定過(guò)程所用過(guò)的洪水作為校核樣本,另810824和870827兩場(chǎng)為預(yù)報(bào)樣本)場(chǎng)洪水的預(yù)報(bào)結(jié)果可知,其中洪峰預(yù)報(bào)誤差在20%以內(nèi)的場(chǎng)次占總體的77%,洪峰誤差小于10%的場(chǎng)次占總體的60%.平均洪峰誤差為12.6%,所有洪峰滯時(shí)均在3 h以內(nèi),平均洪峰滯時(shí)為1.3,所有確定性系數(shù)均大于0.70,平均確定性系數(shù)為0.86;與單一參數(shù)k為隨機(jī)的模型預(yù)報(bào)結(jié)果相比,大部分洪水的洪峰誤差有所降低,確定性系數(shù)稍有提高;平均確定性系數(shù)相當(dāng),而平均洪峰滯時(shí)降低了58%.這說(shuō)明了Nash模型的確存在著較強(qiáng)的“異參同效”現(xiàn)象.兩場(chǎng)預(yù)報(bào)洪水的計(jì)算精度也較高.與僅k為隨機(jī)的情況下的預(yù)報(bào)結(jié)果相比,2參數(shù)均為隨機(jī)的計(jì)算洪峰均方差和80%的置信區(qū)間均有所增大,這說(shuō)明預(yù)報(bào)結(jié)果的不確定性增大了,這也正是由于增加了參數(shù)n的不確定性所致.綜述之,模型參數(shù)的不確定性對(duì)確定性系數(shù)影響較小,對(duì)洪峰誤差、洪峰滯時(shí)和置信區(qū)間影響較大.

圖7繪出了洪號(hào)為810714 a)和810824 b)2場(chǎng)洪水的洪峰后驗(yàn)密度直方圖及其極大似然估計(jì)的理論正態(tài)密度曲線,據(jù)圖看出各洪峰的密度直方圖與估計(jì)的理論正態(tài)密度曲線吻合較好.圖8繪出了這2場(chǎng)洪水的80%的置信區(qū)間與實(shí)測(cè)洪水的比較,據(jù)圖看出每場(chǎng)洪水的實(shí)測(cè)流量幾乎都包括在80%的置信區(qū)間內(nèi).圖9給出了這2場(chǎng)洪水基于AM-MCMC算法的Nash模型2參數(shù)均為隨機(jī)的BFS預(yù)報(bào)均值過(guò)程與實(shí)測(cè)過(guò)程的比較.由圖9可見,2場(chǎng)洪水的擬合精度都很高.

表1 沿渡河流域參數(shù)隨機(jī)的Nash模型的概率洪水預(yù)報(bào)成果表

圖7 洪峰后驗(yàn)密度直方圖及其理論密度曲線

圖8 概率預(yù)報(bào)的80%置信區(qū)間與實(shí)測(cè)過(guò)程比較

圖9 概率預(yù)報(bào)過(guò)程與實(shí)測(cè)過(guò)程比較

4 結(jié) 論

1)貝葉斯概率預(yù)報(bào)系統(tǒng)可與任一復(fù)雜的確定性水文模型協(xié)同工作,而無(wú)需附加任何假設(shè),是制定概率水文預(yù)報(bào)的通用理論框架.

2)AM算法采用并行抽樣,速度快,無(wú)需事先指定MCMC算法的推薦分布,且考慮所抽歷史樣本的信息,能準(zhǔn)確地獲得指定參數(shù)的總體分布特征,具有算法上的優(yōu)越性.

3)AM-MCMC算法能較好獲取Nash模型參數(shù)k,n的后驗(yàn)分布特征,Nash模型的兩個(gè)參數(shù)均存在較強(qiáng)的不確定性,沿渡河流域兩參數(shù)均近似服從正態(tài)分布.使模型的應(yīng)用不再受有限實(shí)測(cè)資料的制約.

4)貝葉斯概率洪水預(yù)報(bào)不僅可給出洪水各時(shí)刻的流量,而且能借助給出的各時(shí)刻的流量方差考慮洪水預(yù)報(bào)的不確定性,便于在實(shí)際應(yīng)用中估計(jì)各種防洪決策的風(fēng)險(xiǎn).

[1]Krzysztofowicz R.Bayesian theory of probabilistic via deterministic hydrologic model[J],Water Resour.Res.,1999,35(9):2 739-2 750.

[2]Krzysztofowicz R,Maranzano C J.Hydrologci uncertainty processor for probabilistic stage transition forecasting[J].Journal of Hydrology,2004,293(1-4):57-73.

[3]張洪剛.貝葉斯概率水文預(yù)報(bào)系統(tǒng)及其應(yīng)用研究[D].武漢:武漢大學(xué)水利水電學(xué)院,2006.

[4]李向陽(yáng),程春田,林劍藝.基于BP網(wǎng)絡(luò)的貝葉斯概率水文預(yù)報(bào)模型[J].水利學(xué)報(bào),2006,37(3):354-359.

[5]王建平,程聲通,賈海峰.基于MCMC法的水質(zhì)模型參數(shù)不確定性研究[J].環(huán)境科學(xué),2006,27(1):24-30.

[6]陳小佳,沈成武.既有橋梁的貝葉斯網(wǎng)絡(luò)評(píng)估方法.武漢理工大學(xué)學(xué)報(bào):交通科學(xué)與工程版,2006,30(1):132-135.

[7]Haario H,Saksman E,Tamminen J.An adaptive metropolis algorithm[J].Bernoulli,2001,7(2):223-242.

[8]Gelman A,Carlin J B,Stren H S,et al.Bayesian data analysis[M].London:Chapmann and Hall,1995.

[9]Gelman A,Rubin D B.Inference from iterative simulation using multiple sequences[J].Statistics Science,1992,7(4):457-511.