孔 強(qiáng)

(石家莊市滹沱河生態(tài)工程運(yùn)維服務(wù)中心,河北 石家莊 050000)

在水文頻率計(jì)算方法上,一直面臨的一個(gè)問題是:基于不同樣本容量水文系列的頻率計(jì)算,會(huì)得出不同的結(jié)果,反映在實(shí)際中,就是很多已建水利工程,當(dāng)初的設(shè)計(jì)水文指標(biāo)值與現(xiàn)實(shí)情勢(shì)不符,須進(jìn)行重新核算。分析原因,較為普遍的觀點(diǎn)主要有:?自然因素和人類活動(dòng)使得水文事件產(chǎn)生條件發(fā)生變化,水文資料存在非平穩(wěn)、非一致性;?水文實(shí)測(cè)樣本容量偏小,不能滿足所需的充分樣本長(zhǎng)度;?選用線型參數(shù)本身存在誤差,影響了計(jì)算結(jié)果精度。 為此,眾多學(xué)者已進(jìn)行了多項(xiàng)研究,如由子系列組成混合分布分模型[1]以解決非一致條件的水文頻率問題,建立變參數(shù)PDS/GP 模型[2]以提高參數(shù)穩(wěn)定性;魯帆等[3]提出了貝葉斯MCMC 洪水頻率分析方法,用以提高分布參數(shù)的估計(jì)精度;戈立婷等[4]和黃華平等[5]分析了GEV 分布所需的充分樣本長(zhǎng)度和應(yīng)用歷史調(diào)查洪水資料保證充分樣本容量的方法。 這些研究對(duì)促進(jìn)水文頻率計(jì)算發(fā)展進(jìn)行了有益的探索。

本次研究通過計(jì)算年降水量樣本參數(shù)的估計(jì)區(qū)間,推斷水文總體正態(tài)分布的可能分布域,采用貝葉斯推斷,以數(shù)值法估計(jì)給定水文樣本的后驗(yàn)概率組合,計(jì)算一定概率下的年降水量值,旨在將正態(tài)分布作為水文事件的可能總體分布。

1 水文樣本序列的正態(tài)性檢驗(yàn)

1.1 頻率推斷與抽樣

水文頻率計(jì)算的基本假定是將水文事件作為隨機(jī)事件,認(rèn)為樣本序列是獨(dú)立隨機(jī)抽取自某一總體分布,其變化規(guī)律服從概率分布律[1]。 通常采用適線法,即選用對(duì)水文序列經(jīng)驗(yàn)點(diǎn)據(jù)擬合較好的線型,如極值Ⅰ型和Ⅱ型、廣義極值分布(GEV)、對(duì)數(shù)正態(tài)分布(LN)、皮爾遜Ⅲ型(P-Ⅲ)等,計(jì)算頻率設(shè)計(jì)值。 對(duì)于頻率線型存在的誤差,有一種觀點(diǎn)認(rèn)為:水文頻率推斷的非平穩(wěn)、非一致性假設(shè)已不滿足,一些“極端”水文事件可能產(chǎn)生系列突變[6]。 從大尺度歷史和未來趨勢(shì)來說,這些已觀測(cè)的“極端”事件,更是一種自然規(guī)律[7]。因此,以幾十個(gè)實(shí)測(cè)數(shù)據(jù)得出的水文系列長(zhǎng)周期突變結(jié)論,還有待跨學(xué)科深入探討。

從統(tǒng)計(jì)抽樣視角,由樣本估計(jì)的確定參數(shù)值,只是所有可能的參數(shù)取值的其中之一,反映在推斷結(jié)果上就是,即便對(duì)樣本擬合再好的線型,也是所有可能線型之一,只是概率不同,而取用基于樣本的某個(gè)確定的線型來代表整體分布,存在巨大理論缺陷,即便樣本擬合度再完美,也是不一定能正確反映總體分布。

1.2 水文事件總體分布正態(tài)性檢驗(yàn)

水文事件的產(chǎn)生過程受多重原因影響,包括自然環(huán)境、人類生產(chǎn)、星系運(yùn)動(dòng)等,對(duì)于水文事件總體分布是正態(tài)分布還是其他某種分布尚無(wú)定論,若將水文事件看成許多隨機(jī)因素疊加的結(jié)果,定性分析或可認(rèn)為水文事件應(yīng)是正態(tài)分布或近似正態(tài)分布,這也是對(duì)水文事件總體分布的合理猜測(cè)。 對(duì)于樣本的正態(tài)性檢驗(yàn)方法有很多,下面采用W 檢驗(yàn)(Shapiro-Wilk 檢)與D檢驗(yàn)(D’Agostino 檢驗(yàn))法[8],對(duì)地區(qū)年降水量等水文樣本進(jìn)行正態(tài)性檢驗(yàn)。 水文樣本容量(3≤n≤50)時(shí),采用W 檢驗(yàn):

設(shè)X1,X2,…,Xn是從總體X中抽取的容量為n的樣本,若X服從正態(tài)分布,應(yīng)滿足

統(tǒng)計(jì)量W為

對(duì)不同的n,系數(shù)αi(i=1,…,n)可查W 檢驗(yàn)統(tǒng)計(jì)量系數(shù)表得到;Wα是在給定的顯著性水平α下使P(W≤Wα)=α的臨界值。

假設(shè)H0:X服從正態(tài)分布,水文樣本容量(n＞50)時(shí), 采用D檢驗(yàn),檢驗(yàn)統(tǒng)計(jì)量D為

在H0為真時(shí),E(D) ≈0.28209479,令

D’Agostino 用隨機(jī)模擬法得到了Y的分位數(shù)表,在給定顯著性水平α后,用統(tǒng)計(jì)量Y進(jìn)行檢驗(yàn)的拒絕域?yàn)?,采用D 檢驗(yàn),由于樣本Y未落入拒絕域中,故在α= 0.05 時(shí),可認(rèn)為服從正態(tài)分布。

給定顯著性水平α=0.05 時(shí),由式(1)計(jì)算了資料中43 個(gè)城市樣本容量為22 的W值(n= 22,Wα=0.911),有41 個(gè)通過了正態(tài)檢驗(yàn),5 個(gè)水文站或入庫(kù)徑流量中,有3 個(gè)通過了正態(tài)檢驗(yàn);16 個(gè)城市樣本容量為36 的W值(n= 36,Wα= 0.935),有15 個(gè)通過了正態(tài)檢驗(yàn),5 個(gè)水文站或入庫(kù)徑流量中,有2 個(gè)通過了正態(tài)檢驗(yàn);由式(2)計(jì)算了7 個(gè)城市樣本容量為60 的D值(n=60,=1.13),全都通過了正態(tài)檢驗(yàn)。 可見降水量無(wú)論是W檢驗(yàn)還是D檢驗(yàn),都較好地滿足正態(tài)性。 而水文站或入庫(kù)徑流量正態(tài)性較差,可能是受人為影響因素較大,計(jì)算結(jié)果見表1。

表1 中,南昌、武漢、長(zhǎng)沙、?？?、重慶、成都、貴陽(yáng)、昆明、拉薩、西安、蘭州、西寧、烏魯木齊、天津、石家莊、太原、呼和浩特、沈陽(yáng)、哈爾濱、上海、杭州、合肥、福州、南京、銀川、北京、濟(jì)南、鄭州、廣州、南寧30個(gè)城市的1998—2020 年降水量數(shù)據(jù)均取自1998—2020 年的《中國(guó)統(tǒng)計(jì)年鑒》;榆林[9]、廬江[10]、永吉[11]、張家口[12]、石嘴山[13]、寶雞[14]、南京[15]、銀川[16]、咸陽(yáng)[17]、密云水庫(kù)入庫(kù)徑流[18]、灤河灤縣站年徑流[19]、北京[20]、濟(jì)南[21]、鄭州[21]、廣州[23]、南寧[24]、韶關(guān)[25]、隴南[26]、溫州[27]、洪家渡水文站構(gòu)皮灘水文站[28]20 個(gè)地區(qū)的年降水量或徑流量均來自公開發(fā)表論文或技術(shù)資料。

2 推斷方法

2.1 均值和方差的估計(jì)

對(duì)正態(tài)分布的均值和方差的估計(jì),采用極大似然估計(jì)量[8],即總體均值μ=X,總體方差樣本的極大似然值不代表所有可能取值,需要對(duì)、S2可能取值作區(qū)間估計(jì)。

2.1.1 總體均值μ的區(qū)間估計(jì)

由于總體方差σ2未知,可用σ2的無(wú)偏估計(jì)S2代替,則得隨機(jī)變量為含有未知參數(shù)μ的樞軸量[8],服從T～tn-1分布,由于t分布的概率密度函數(shù)是單峰對(duì)稱的,查t分布表可得自由度為n-1 的t分布的上側(cè)分位數(shù),使得

故μ的置信度為1-α的置信區(qū)間為

2.1.2 總體方差σ的區(qū)間估計(jì)2

樣本方差S2是σ2的無(wú)偏估計(jì),得隨機(jī)變量

由于χ2分布的概率密度函數(shù)圖形是單峰不對(duì)稱的,由χ2分布的上側(cè)分位數(shù)定義,有

故σ2的置信度為1-α的置信區(qū)間為

則標(biāo)準(zhǔn)差σ的置信度為1-α的置信區(qū)間為

取顯著性水平α=0.05,則置信區(qū)間為1-α=0.95。 對(duì)于正態(tài)分布的樣本均值置信區(qū)間取t分布,方差取χ2分布,分別計(jì)算樣本均值和方差置信區(qū)間。

2.2 可能的正態(tài)總體

正態(tài)分布樣本的均值和方差S2相互獨(dú)立[8],即值與S2值無(wú)關(guān)。 若水文樣本來自均值、方差估計(jì)區(qū)間的邊界值范圍內(nèi)的正態(tài)分布,則水文樣本分布必然會(huì)收斂于均值、方差估計(jì)區(qū)間內(nèi)的正態(tài)分布域。 因此樣本均值置信區(qū)間內(nèi)的所有可能值與方差置信區(qū)間內(nèi)所有可能值,組合成的是有界正態(tài)分布簇,其邊緣為置信區(qū)間上下限的4 個(gè)組合:(均值下限;標(biāo)準(zhǔn)差下限)、(均值下限;標(biāo)準(zhǔn)差上限)、(均值上限;標(biāo)準(zhǔn)差下限)、(均值上限;標(biāo)準(zhǔn)差上限),1 個(gè)最大似然正態(tài)分布是(樣本均值;樣本標(biāo)準(zhǔn)差)。

2.3 后驗(yàn)概率分布

2.3.1 貝葉斯公式

貝葉斯公式可以將先驗(yàn)知識(shí)和樣本信息結(jié)合起來進(jìn)行推斷,基本公式為

式中:θ代表參數(shù);x代表樣本;f(θ|x) 為后驗(yàn)概率密度函數(shù);π(θ) 為先驗(yàn)分布;f(x|θ) 為似然函數(shù)。

2.3.2 先驗(yàn)概率分布

前述已證實(shí),區(qū)域年降水量等非人為控制水文樣本總體趨向于正態(tài)分布,因此,可以水文事件正態(tài)分布作為貝葉斯推斷的先驗(yàn)分布:

式中:N(x;μ,σ) 為正態(tài)分布概率密度函數(shù);x為水文序列值;μ為均值;σ為標(biāo)準(zhǔn)差。

2.3.3 似然函數(shù)

由已知樣本擬合得到的分布曲線,如P-Ⅲ分布、GEV 分布等都是基于已知樣本的可能分布,其分布是總體分布的似然分布,其分布函數(shù)為似然函數(shù),若采用常用的P-Ⅲ型曲線作為擬合線型[29],則式(1)中似然函數(shù)f(x|θ) 可寫為

式中:Г(·) 為gamma 函數(shù);α0,β,α分別為分布的位置、尺度、形狀參數(shù),且有α≥α0、α ＞0、β ＞0,此3 個(gè)參數(shù)與總體參數(shù)Ex、Cv、Cs(期望值、變差系數(shù)、偏態(tài)系數(shù))關(guān)系如下:

2.3.4 后驗(yàn)概率分布

若先驗(yàn)分布π(θ) 和似然函數(shù)f(x|θ) 已知,則可以求出二者的乘積函數(shù)g(θ|x):

對(duì)乘積函數(shù)g(θ|x) 在Θ 上求解積分值k,則可由式(8)得出后驗(yàn)密度函數(shù)。 貝葉斯公式中,分母為兩個(gè)函數(shù)的積函數(shù)在積分域Θ 上,對(duì)參數(shù)θ的積分,理論上,先驗(yàn)分布π(θ) 和似然函數(shù)f(x|θ) 的積分域是[- ∞,+ ∞],一般采用馬爾可夫鏈蒙特卡羅(MCMC)方法求解貝葉斯公式。 而對(duì)于特定水文事件,定義積分域[- ∞,+ ∞] 是理論極限,在可遇見的未來,或存在實(shí)際的降水或徑流的上限和下限,故可取特定的計(jì)算區(qū)間[a,b], 下限為可能最小值,上限為可能最大值,選取計(jì)算步長(zhǎng),以兩個(gè)函數(shù)值乘積函數(shù)為分子目標(biāo)函數(shù)值,則式(5)可寫為

3 后驗(yàn)概率分布的數(shù)值求解

以南寧市1961—2020 年降水量樣本均值(μ=1302.14)、標(biāo)準(zhǔn)差(σ=227.241)生成正態(tài)分布隨機(jī)數(shù)為算例,對(duì)南寧市年降水量后驗(yàn)概率進(jìn)行數(shù)值求解。

3.1 確定正態(tài)分布組合

由式(3)和式(4)計(jì)算南寧市年降水量均值、方差的95%置信水平的取值范圍邊界值,可組合出5 個(gè)正態(tài)總體分布,見表2。

表2 南寧市降水量均值、方差置信區(qū)間邊界值正態(tài)分布組合

作為可能的正態(tài)先驗(yàn)密度函數(shù)簇邊界概率,由式(6)可計(jì)算5 個(gè)正態(tài)總體分布先驗(yàn)概率密度曲線,見圖1。

圖1 組合1 ～5 正態(tài)可能先驗(yàn)概率密度曲線

3.2 計(jì)算似然函數(shù)

對(duì)樣本數(shù)據(jù)以P-Ⅲ曲線進(jìn)行擬合,由式(7)求得P-Ⅲ擬合曲線密度函數(shù)(三參數(shù)gamma 分布)[29],作為貝葉斯公式的似然函數(shù),P-Ⅲ擬合頻率與樣本計(jì)算頻率擬合情況見圖2。

圖2 樣本頻率與P-Ⅲ擬合曲線對(duì)比

樣本頻率與P-Ⅲ擬合曲線擬合度及誤差計(jì)算見表3。

表3 南寧市年降水量P-Ⅲ擬合曲線擬誤差計(jì)算表

表4 傅里葉多項(xiàng)式參數(shù)

3.3 計(jì)算乘積函數(shù)

正態(tài)先驗(yàn)概率密度函數(shù)式(6)與似然函數(shù)式(7)相乘的解析解求解較為復(fù)雜,可采用數(shù)值法,積分域Θ取[400,2500],計(jì)算步長(zhǎng)設(shè)置為1mm,計(jì)算似然函數(shù)值與正態(tài)密度函數(shù)值的乘積值,以matlab 曲線擬合工具中8 參數(shù)傅里葉多項(xiàng)式對(duì)乘積值曲線進(jìn)行擬合,得出似然函數(shù)與組合1 ～5 可能正態(tài)先驗(yàn)密度函數(shù)乘積的函數(shù)擬合公式:

擬合公式的R-square 和Adjusted R-square 值均為1,SSE 和RSME 值見表5。 以組合1 正態(tài)先驗(yàn)值與P-Ⅲ似然函數(shù)數(shù)值計(jì)算,先驗(yàn)概率密度曲線、似然函數(shù)曲線、后驗(yàn)概率密度曲線見圖3,組合2 ～5 計(jì)算過程相同。

圖3 組合1 概率密度數(shù)值計(jì)算對(duì)比

表5 組合1 ～5 擬合優(yōu)度及k 值

3.4 計(jì)算k 值

在給定積分域Θ[400,2500]與乘積函數(shù)g(θ|x)曲線頂點(diǎn)值max(g(θ|x)) 圍成的矩形范圍內(nèi),采用MC 法計(jì)算5 個(gè)乘積函數(shù)g(θ|x) 在Θ 上的積分值,隨機(jī)取10 萬(wàn)個(gè)數(shù)值(1mm 單位),計(jì)算得g(θ|x) 所圍圖形面積,由式(8)可求得k值,見表5。

將k值代入式(9)可得5 個(gè)后驗(yàn)概率密度函數(shù)。組合1 ～5 后驗(yàn)概率密度曲線見圖4。

圖4 組合1 ～5 后驗(yàn)概率密度曲線對(duì)比

3.5 計(jì)算給定頻率的降水量值

給定頻率值,帶入組合1 ～5 的后驗(yàn)概率密度函數(shù),即可求出0.95 置信水平的降水量取值區(qū)間。 對(duì)圖4 中5 個(gè)后驗(yàn)概率密度取累計(jì)值,可得相應(yīng)的概率質(zhì)量函數(shù)曲線,見圖5。 對(duì)于給定的某個(gè)頻率值,其降水量取值為組合1 ～5 對(duì)應(yīng)的取值范圍:頻率0.10 =[1394,1492],頻率0.05=[1451,1562],頻率0.01=[1567,1693]。 由于正態(tài)分布總體均值和方差置信區(qū)間的對(duì)稱估計(jì),本方法數(shù)值計(jì)算結(jié)果相較于常用的水文頻率設(shè)計(jì)值的估計(jì)量抽樣分布的標(biāo)準(zhǔn)差計(jì)算方法[9]和含義都不同。

圖5 組合1 ～5 后驗(yàn)概率質(zhì)量曲線對(duì)比

4 結(jié) 論

4.1 推斷結(jié)果是1 個(gè)取值域

本方法將對(duì)水文事件的先驗(yàn)知識(shí)融入了頻率推斷過程,通過正態(tài)先驗(yàn)分布參數(shù)邊緣值組合,利用貝葉斯方法,靈活運(yùn)用數(shù)值模擬計(jì)算和MC 法,使水文頻率計(jì)算符合對(duì)事物的認(rèn)知,在水文頻率推斷方向上是正確的。 對(duì)給定樣本的水文序列頻率計(jì)算結(jié)果是1 個(gè)取值范圍,這個(gè)結(jié)果含義與教材中設(shè)計(jì)值的估計(jì)量抽樣分布的標(biāo)準(zhǔn)差計(jì)算方法和含義是完全不同的,在一定程度上解決了常用的擬合曲線可能很適合樣本數(shù)據(jù),但預(yù)測(cè)仍欠佳的問題,推斷結(jié)果所確定的值域有較好的概率理論依據(jù),豐富了水文頻率計(jì)算方法。

4.2 對(duì)樣本的擬合是誤差主要來源

對(duì)樣本數(shù)據(jù)的擬合有多種方法,每種方法會(huì)得出不同的擬合函數(shù),故而擬合函數(shù)對(duì)計(jì)算結(jié)果產(chǎn)生的誤差是不可避免的,而且是誤差的主要來源。 相較而言,計(jì)算過程中對(duì)組合1 ～5 的先驗(yàn)和后驗(yàn)的數(shù)值擬合函數(shù),以及MC 方法在設(shè)定積分域上對(duì)k值的模擬計(jì)算,所產(chǎn)生的誤差會(huì)隨著增加計(jì)算量而大幅減小。

4.3 采用數(shù)值法對(duì)貝葉斯公式直接求解是可行的

理論上,水文事件概率分布的兩端為[- ∞,+ ∞],根據(jù)可遇見的實(shí)際情況,將不能發(fā)生情況的概率認(rèn)定為0,對(duì)貝葉斯公式中分母積分域設(shè)置為固定值,采用MC 方法計(jì)算定積分,進(jìn)而求解k值。 采用數(shù)值法計(jì)算貝葉斯公式,直接模擬了貝葉斯公式分子乘積函數(shù),為貝葉斯公式在水文計(jì)算中的應(yīng)用提供了新方法。