李林波 魏 為 黃國芳
(云南省水利水電工程有限公司,云南 昆明 650000)
頻率曲線是推求設(shè)計(jì)降雨量、設(shè)計(jì)洪水等水文設(shè)計(jì)值的重要曲線[1]。該曲線一般繪制在頻率格紙(海森幾率格紙)上,但由于頻率格紙橫坐標(biāo)軸不均勻,傳統(tǒng)繪制方法多是在已繪好的頻率格紙上徒手繪制,不僅效率低,而且精度不高,存在人為誤差[2]。加之頻率曲線的求解過程計(jì)算量較大,給工作人員帶來諸多不便[3]。
前人在該領(lǐng)域的研究多為介紹頻率曲線求解原理,或講述求解步驟[4],基于Excel軟件對此類問題進(jìn)行的研究多止于繪制頻率格紙,鮮有利用Excel軟件求解頻率曲線的研究,因此本文將著重解決這一問題,對其進(jìn)行研究。
水文頻率曲線是水文要素值與其發(fā)生頻率之間的關(guān)系曲線,當(dāng)設(shè)計(jì)標(biāo)準(zhǔn)[設(shè)計(jì)頻率(P)或設(shè)計(jì)重現(xiàn)期(T)]已知時(shí),可在該曲線上查得相應(yīng)的水文設(shè)計(jì)值(xp),其意思是任何水文要素值x超過或等于xp的頻率為P,即
P(x≥xp)=P
(1)
水文頻率曲線有多種形式?!端姽こ淘O(shè)計(jì)洪水計(jì)算規(guī)范》(SL 44—2006)規(guī)定應(yīng)使用皮爾遜Ⅲ型曲線(P-Ⅲ型曲線),常見的有降雨量頻率曲線、洪峰流量頻率曲線[5-6]。
皮爾遜Ⅲ型曲線是一條一端有限一端無限的不對稱單峰、正偏曲線[7],數(shù)學(xué)上常稱伽瑪分布,其概率密度函數(shù)為
(2)
式中:Γ(α)為α的伽瑪函數(shù);α、β、a0為皮爾遜Ⅲ型分布的形狀尺度和位置未知參數(shù),α>0,β>0。
(3)
(4)
(5)
水文計(jì)算中,一般需要求出隨機(jī)變量取值xp所相應(yīng)的指定頻率P,也就是對密度曲線進(jìn)行積分,即
(6)
實(shí)際計(jì)算時(shí)通常變換成如下形式進(jìn)行積分
(7)
(8)
式中:Φ為離均系數(shù)。
式(7)中被積函數(shù)只含有一個(gè)待定參數(shù)Cs,只需要假定一Cs值,便可通過積分求出P與Φ之間的關(guān)系。
對于若干給定的Cs值,當(dāng)Cs與Cv成倍數(shù)關(guān)系時(shí),Φ和P有對應(yīng)數(shù)值表,可查閱“皮爾遜Ⅲ型頻率曲線的離均系數(shù)Φp值表”。
(9)
對于若干給定的Cs值,當(dāng)Cs等于Cv的一定倍數(shù)時(shí),P-Ⅲ型頻率曲線的模比系數(shù)Kp也已制成表格,可查閱“皮爾遜Ⅲ型頻率曲線的模比系數(shù)Kp值表”。
(10)
頻率格紙又叫海森幾率格紙,它是一種特殊的坐標(biāo)系統(tǒng),其縱坐標(biāo)為均勻分格或?qū)?shù)分格的常規(guī)數(shù)學(xué)坐標(biāo),橫坐標(biāo)與頻率值(下側(cè)概率)的標(biāo)準(zhǔn)正態(tài)分布分位數(shù)有關(guān),兩端分格較稀而中間較密[8]。
由于標(biāo)準(zhǔn)正態(tài)分布分位數(shù)在P=50%處為零,而海森幾率格紙?jiān)赑=0.01%時(shí)的橫坐標(biāo)值為零,因此海森幾率格紙橫坐標(biāo)值可采用如下計(jì)算公式進(jìn)行計(jì)算:
LP=-U0.01%+UP
(11)
式中:LP為頻率P對應(yīng)的橫坐標(biāo)值;UP為頻率P對應(yīng)的標(biāo)準(zhǔn)正態(tài)分布分位數(shù);U0.01%為頻率P=0.01%對應(yīng)的標(biāo)準(zhǔn)正態(tài)分布分位數(shù)。
(12)
(13)
a.點(diǎn)繪經(jīng)驗(yàn)點(diǎn)據(jù)。如果用等概率公式P(x≥xi)=m/n的經(jīng)驗(yàn)分布曲線估計(jì)總體分布曲線,存在不合理現(xiàn)象。當(dāng)m=n時(shí),最末項(xiàng)的頻率為100%,即樣本最小值為總體中的最小值,不符合事實(shí)。因此,水文上用期望值公式估計(jì)頻率:
(14)
式中:m為按降水量遞減順序序列排列的項(xiàng)數(shù);n為序列總項(xiàng)數(shù)。
利用式(14)求出經(jīng)驗(yàn)頻率后,以經(jīng)驗(yàn)頻率為橫坐標(biāo)、變量值為縱坐標(biāo),點(diǎn)繪經(jīng)驗(yàn)點(diǎn)據(jù)。
d.選擇一條與經(jīng)驗(yàn)點(diǎn)據(jù)配合最佳的曲線作為采用曲線。該曲線的參數(shù)看作總體參數(shù)的估計(jì)值。
選取與經(jīng)驗(yàn)點(diǎn)據(jù)配合最佳的理論頻率曲線為所求頻率曲線,通過讀圖可以得到不同頻率P對應(yīng)的xp值。
由于頻率比較抽象,為便于理解,常采用重現(xiàn)期表示。重現(xiàn)期是指在許多試驗(yàn)中,某一事件重復(fù)出現(xiàn)的時(shí)間間隔的平均數(shù),它是頻率的倒數(shù)。在工程水文中,重現(xiàn)期用T表示,一般以年為單位[10]:
(15)
如當(dāng)暴雨或洪水頻率為1%時(shí),即重現(xiàn)期T=100年,稱此為百年一遇的暴雨或洪水。
表1 興國縣2000—2017年年降水量 單位:mm
a.在表格第一行依次輸入年份、年降水量R、序號等項(xiàng)目,見圖1。
b.做好表頭后,首先將A、B欄已知年份及其相應(yīng)的年降雨量輸入計(jì)算表格內(nèi),然后在C欄填入序號1~18。之后在D欄對B欄的降水資料進(jìn)行由大到小的排序。
e.計(jì)算“Ki-1和(Ki-1)2。在F2處輸入“=E2-1”,回車;在G2處輸入“=F2^2”,回車。然后分別將光標(biāo)移至F2和G2的右下角,當(dāng)出現(xiàn)十字標(biāo)識后,向下拖動至需要計(jì)算數(shù)據(jù)行。
圖1 經(jīng)驗(yàn)頻率曲線統(tǒng)計(jì)參數(shù)計(jì)算表
頻率格紙的橫向網(wǎng)格線為均勻分布,可直接由Excel軟件的圖表功能自動生成,而縱向網(wǎng)格線坐標(biāo)值是按照標(biāo)準(zhǔn)正態(tài)頻率曲線拉成一條直線的原理計(jì)算得出的,中間較密,兩端較稀,不能由Excel軟件的圖表功能自動生成。
因此,繪制頻率格紙的關(guān)鍵是繪制頻率格紙的縱向網(wǎng)格線,可以先用Excel軟件中的內(nèi)置函數(shù)NORMSINV(P)將不同頻率轉(zhuǎn)化為相應(yīng)的橫坐標(biāo)值,然后插入散點(diǎn)圖。函數(shù)NORMSINV(P)為返回標(biāo)準(zhǔn)正態(tài)累積分布函數(shù)的反函數(shù)值,其計(jì)算結(jié)果的精度可達(dá)到±3×10-7,該函數(shù)的詳細(xì)說明和用法可參考Excel軟件的幫助。
繪制頻率格紙的步驟具體如下:
a.新建Excel工作簿,在第一行依次輸入“頻率P”“P對應(yīng)的標(biāo)準(zhǔn)正態(tài)分布分位數(shù)”“X軸”“Y軸”,見圖2。
b.計(jì)算不同頻率P對應(yīng)的橫坐標(biāo)Lp值。在A欄中輸入頻率值0.01~99.99。在A2、A3處分別輸入“0.01”,在A4、A5處分別輸入“0.02”……依此類推,在A欄后續(xù)單元格中輸入縱向網(wǎng)絡(luò)線對應(yīng)頻率值,直至“99.99”止。
c.計(jì)算頻率P對應(yīng)的標(biāo)準(zhǔn)正態(tài)分布分位數(shù)Up。在B2單元格中輸入“=NORMSINV(A2%)”,回車,然后利用Excel填充柄向下填充至P=99.99%對應(yīng)單元格,就得到了頻率P對應(yīng)的標(biāo)準(zhǔn)正態(tài)分布分位數(shù)Up。
c.計(jì)算頻率格紙橫坐標(biāo)Lp值。由于海森幾率格紙?jiān)赑=0.01%時(shí)的橫坐標(biāo)值為零,根據(jù)其橫坐標(biāo)值計(jì)算公式LP=-U0.01%+UP,在C2處輸入“=-$B$2+B2”,回車,然后利用Excel填充柄向下填充,就得到了頻率P對應(yīng)的橫坐標(biāo)Lp值。見圖2中C列。
d.設(shè)置縱坐標(biāo)軸邊界。由降水資料可知最大年降水量Rmax為2265.8mm,為方便頻率曲線展布,可將Y軸最大值取大一點(diǎn),本次繪制取Y軸范圍為0~4000。在D列依次輸入“0、4000、4000、0、0、4000、4000、0、0、4000、4000…”,直到P=99.99%對應(yīng)單元格為止。見圖2。
e.繪制網(wǎng)格線。?插入散點(diǎn)圖,選中C、D兩列,以C、D兩列數(shù)據(jù)分別作為散點(diǎn)圖的橫坐標(biāo)值和縱坐標(biāo)值,選擇插入圖表命令,在“所有圖表”中選擇“XY散點(diǎn)圖”,然后選擇“帶直線的散點(diǎn)圖”,回車,見圖2;?設(shè)置縱向網(wǎng)格線及橫坐標(biāo)刻度,刪除原縱向網(wǎng)絡(luò)線,選中所繪散點(diǎn)圖,將線條改為實(shí)線、黑色、0.5磅寬;將橫坐標(biāo)邊界設(shè)置為0~7.438,并刪除原橫坐標(biāo)刻度,選中剛繪制的散點(diǎn)圖,右擊鼠標(biāo),在彈出的選項(xiàng)中選擇“添加數(shù)據(jù)標(biāo)簽”,選中出現(xiàn)的數(shù)據(jù)標(biāo)簽并設(shè)置為只顯示“X值”,拖動標(biāo)簽位置并修改標(biāo)簽內(nèi)容為相應(yīng)縱線對應(yīng)的頻率P值,最后刪除多余標(biāo)簽便完成了頻率格紙縱線的繪制;?設(shè)置橫向網(wǎng)格線及橫坐標(biāo)刻度,選中橫坐標(biāo)刻度,右擊鼠標(biāo),選擇“設(shè)置坐標(biāo)軸格式”,設(shè)置邊界范圍為0~4000,間距200,然后選中橫向網(wǎng)格線,從鼠標(biāo)右鍵菜單中選擇“設(shè)置網(wǎng)格線格式”,修改線條為實(shí)線、黑色、0.5磅。至此便完成了海森幾率格紙的繪制。
圖2 頻率格紙橫、縱坐標(biāo)計(jì)算表與散點(diǎn)圖繪制
3.3.1 點(diǎn)繪經(jīng)驗(yàn)點(diǎn)據(jù)
利用公式“NORMSINV(P′)”計(jì)算得到經(jīng)驗(yàn)頻率P′對應(yīng)的標(biāo)準(zhǔn)正態(tài)分布分位數(shù),再代入橫坐標(biāo)值計(jì)算公式“LP=-U0.01%+UP”可得到經(jīng)驗(yàn)頻率P′對應(yīng)的橫坐標(biāo)值,見圖3。以其為橫坐標(biāo),以降序排列的降雨量資料為縱坐標(biāo),在圖表區(qū)添加系列,參照圖4設(shè)置系列名稱并框選X軸與Y軸系列值,然后回車。最后設(shè)置經(jīng)驗(yàn)頻率散點(diǎn)圖為無線條散點(diǎn)圖,即完成了經(jīng)驗(yàn)點(diǎn)據(jù)的繪制,見圖4。
圖3 經(jīng)驗(yàn)點(diǎn)據(jù)橫、縱坐標(biāo)值計(jì)算表及其繪制
圖4 興國縣降水量頻率曲線
3.3.2 初定統(tǒng)計(jì)參數(shù)
3.3.3 選配理論頻率曲線
由圖4可知,理論頻率曲線(Cs=3Cv=0.75)的中段與經(jīng)驗(yàn)頻率點(diǎn)據(jù)配合較好,但頭部偏于經(jīng)驗(yàn)頻率點(diǎn)據(jù)下方,而尾部又偏于經(jīng)驗(yàn)頻率點(diǎn)據(jù)上方,配合情況不理想。需調(diào)整參數(shù)進(jìn)行重新配線。
圖5 第一次配線計(jì)算結(jié)果與曲線繪制
b.改變參數(shù)進(jìn)行第二次配線,由第一次配線結(jié)果可知,需增大Cv值。取Cv=0.3,Cs=3Cv=0.9,參照選配理論頻率曲線第一步進(jìn)行配線,計(jì)算結(jié)果見圖6。由配線結(jié)果(圖4)可知,理論頻率曲線頭部配合較好,但尾部較經(jīng)驗(yàn)點(diǎn)據(jù)偏低。
c.在第二次配線基礎(chǔ)上,為使尾部抬高一些與經(jīng)驗(yàn)點(diǎn)據(jù)相配合,需增大Cs值,因此,取Cv=0.3,Cs=3.5Cv=1.05,再次計(jì)算理論頻率曲線,計(jì)算結(jié)果見圖6。由配線結(jié)果(圖4)可知,當(dāng)Cs=3.5Cv=1.05時(shí),理論曲線與經(jīng)驗(yàn)頻率點(diǎn)據(jù)配合較合適,讀圖得到保證率為20%、50%、75%及95%時(shí)年降雨量分別為1931.85mm、1504.31mm、1235.11mm、965.92mm(見表2)。
圖6 選配理論曲線計(jì)算結(jié)果
表2 興國縣不同保證率下年降水量對比情況
利用Excel軟件求解水文設(shè)計(jì)值、繪制頻率格紙和水文頻率曲線切實(shí)可行,該方法不僅能提高計(jì)算和制圖的效率,出圖效果更加整潔美觀,制圖精度相較于傳統(tǒng)的手繪曲線有所提高,為水文工作者提供了一種方便、簡潔、美觀的繪圖方法。