Excel在求解水文頻率曲線中的應(yīng)用

李林波 魏 為 黃國芳

(云南省水利水電工程有限公司，云南 昆明 650000)

頻率曲線是推求設(shè)計降雨量、設(shè)計洪水等水文設(shè)計值的重要曲線[1]。該曲線一般繪制在頻率格紙(海森幾率格紙)上，但由于頻率格紙橫坐標(biāo)軸不均勻，傳統(tǒng)繪制方法多是在已繪好的頻率格紙上徒手繪制，不僅效率低，而且精度不高，存在人為誤差[2]。加之頻率曲線的求解過程計算量較大，給工作人員帶來諸多不便[3]。

前人在該領(lǐng)域的研究多為介紹頻率曲線求解原理，或講述求解步驟[4]，基于Excel軟件對此類問題進(jìn)行的研究多止于繪制頻率格紙，鮮有利用Excel軟件求解頻率曲線的研究，因此本文將著重解決這一問題，對其進(jìn)行研究。

1 基本概念與原理

1.1 水文頻率曲線

水文頻率曲線是水文要素值與其發(fā)生頻率之間的關(guān)系曲線，當(dāng)設(shè)計標(biāo)準(zhǔn)[設(shè)計頻率(P)或設(shè)計重現(xiàn)期(T)]已知時，可在該曲線上查得相應(yīng)的水文設(shè)計值(xp)，其意思是任何水文要素值x超過或等于xp的頻率為P，即

P(x≥xp)=P

(1)

水文頻率曲線有多種形式。《水利水電工程設(shè)計洪水計算規(guī)范》(SL 44—2006)規(guī)定應(yīng)使用皮爾遜Ⅲ型曲線(P-Ⅲ型曲線)，常見的有降雨量頻率曲線、洪峰流量頻率曲線[5-6]。

皮爾遜Ⅲ型曲線是一條一端有限一端無限的不對稱單峰、正偏曲線[7]，數(shù)學(xué)上常稱伽瑪分布，其概率密度函數(shù)為

(2)

式中：Γ(α)為α的伽瑪函數(shù)；α、β、a0為皮爾遜Ⅲ型分布的形狀尺度和位置未知參數(shù)，α>0，β>0。

(3)

(4)

(5)

水文計算中，一般需要求出隨機(jī)變量取值xp所相應(yīng)的指定頻率P，也就是對密度曲線進(jìn)行積分，即

(6)

實際計算時通常變換成如下形式進(jìn)行積分

(7)

(8)

式中：Φ為離均系數(shù)。

式(7)中被積函數(shù)只含有一個待定參數(shù)Cs，只需要假定一Cs值，便可通過積分求出P與Φ之間的關(guān)系。

對于若干給定的Cs值，當(dāng)Cs與Cv成倍數(shù)關(guān)系時，Φ和P有對應(yīng)數(shù)值表，可查閱“皮爾遜Ⅲ型頻率曲線的離均系數(shù)Φp值表”。

(9)

對于若干給定的Cs值，當(dāng)Cs等于Cv的一定倍數(shù)時，P-Ⅲ型頻率曲線的模比系數(shù)Kp也已制成表格，可查閱“皮爾遜Ⅲ型頻率曲線的模比系數(shù)Kp值表”。

(10)

1.2 頻率格紙

頻率格紙又叫海森幾率格紙，它是一種特殊的坐標(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%處為零，而海森幾率格紙在P=0.01%時的橫坐標(biāo)值為零，因此海森幾率格紙橫坐標(biāo)值可采用如下計算公式進(jìn)行計算：

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ù)。

2 水文頻率曲線求解方法與步驟

2.1 求解方法

(12)

(13)

2.2 求解步驟

a.點繪經(jīng)驗點據(jù)。如果用等概率公式P(x≥xi)=m/n的經(jīng)驗分布曲線估計總體分布曲線，存在不合理現(xiàn)象。當(dāng)m=n時，最末項的頻率為100%，即樣本最小值為總體中的最小值，不符合事實。因此，水文上用期望值公式估計頻率：

(14)

式中：m為按降水量遞減順序序列排列的項數(shù)；n為序列總項數(shù)。

利用式(14)求出經(jīng)驗頻率后，以經(jīng)驗頻率為橫坐標(biāo)、變量值為縱坐標(biāo)，點繪經(jīng)驗點據(jù)。

d.選擇一條與經(jīng)驗點據(jù)配合最佳的曲線作為采用曲線。該曲線的參數(shù)看作總體參數(shù)的估計值。

2.3 計算結(jié)果與表示

選取與經(jīng)驗點據(jù)配合最佳的理論頻率曲線為所求頻率曲線，通過讀圖可以得到不同頻率P對應(yīng)的xp值。

由于頻率比較抽象，為便于理解，常采用重現(xiàn)期表示。重現(xiàn)期是指在許多試驗中，某一事件重復(fù)出現(xiàn)的時間間隔的平均數(shù)，它是頻率的倒數(shù)。在工程水文中，重現(xiàn)期用T表示，一般以年為單位[10]：

(15)

如當(dāng)暴雨或洪水頻率為1%時，即重現(xiàn)期T=100年，稱此為百年一遇的暴雨或洪水。

3 計算實例

表1 興國縣2000—2017年年降水量 單位：mm

3.1 計算經(jīng)驗頻率曲線統(tǒng)計參數(shù)

a.在表格第一行依次輸入年份、年降水量R、序號等項目，見圖1。

b.做好表頭后，首先將A、B欄已知年份及其相應(yīng)的年降雨量輸入計算表格內(nèi)，然后在C欄填入序號1～18。之后在D欄對B欄的降水資料進(jìn)行由大到小的排序。

e.計算“Ki-1和(Ki-1)2。在F2處輸入“=E2-1”，回車；在G2處輸入“=F2^2”，回車。然后分別將光標(biāo)移至F2和G2的右下角，當(dāng)出現(xiàn)十字標(biāo)識后，向下拖動至需要計算數(shù)據(jù)行。

圖1 經(jīng)驗頻率曲線統(tǒng)計參數(shù)計算表

3.2 繪制頻率格紙

頻率格紙的橫向網(wǎng)格線為均勻分布，可直接由Excel軟件的圖表功能自動生成，而縱向網(wǎng)格線坐標(biāo)值是按照標(biāo)準(zhǔn)正態(tài)頻率曲線拉成一條直線的原理計算得出的，中間較密，兩端較稀，不能由Excel軟件的圖表功能自動生成。

因此，繪制頻率格紙的關(guān)鍵是繪制頻率格紙的縱向網(wǎng)格線，可以先用Excel軟件中的內(nèi)置函數(shù)NORMSINV(P)將不同頻率轉(zhuǎn)化為相應(yīng)的橫坐標(biāo)值，然后插入散點圖。函數(shù)NORMSINV(P)為返回標(biāo)準(zhǔn)正態(tài)累積分布函數(shù)的反函數(shù)值，其計算結(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.計算不同頻率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.計算頻率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.計算頻率格紙橫坐標(biāo)Lp值。由于海森幾率格紙在P=0.01%時的橫坐標(biāo)值為零，根據(jù)其橫坐標(biāo)值計算公式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軸最大值取大一點，本次繪制取Y軸范圍為0～4000。在D列依次輸入“0、4000、4000、0、0、4000、4000、0、0、4000、4000…”，直到P=99.99%對應(yīng)單元格為止。見圖2。

e.繪制網(wǎng)格線。?插入散點圖，選中C、D兩列，以C、D兩列數(shù)據(jù)分別作為散點圖的橫坐標(biāo)值和縱坐標(biāo)值，選擇插入圖表命令，在“所有圖表”中選擇“XY散點圖”，然后選擇“帶直線的散點圖”，回車，見圖2；?設(shè)置縱向網(wǎng)格線及橫坐標(biāo)刻度，刪除原縱向網(wǎng)絡(luò)線，選中所繪散點圖，將線條改為實線、黑色、0.5磅寬；將橫坐標(biāo)邊界設(shè)置為0～7.438，并刪除原橫坐標(biāo)刻度，選中剛繪制的散點圖，右擊鼠標(biāo)，在彈出的選項中選擇“添加數(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)格線格式”，修改線條為實線、黑色、0.5磅。至此便完成了海森幾率格紙的繪制。

圖2 頻率格紙橫、縱坐標(biāo)計算表與散點圖繪制

3.3 適線法配線

3.3.1 點繪經(jīng)驗點據(jù)

利用公式“NORMSINV(P′)”計算得到經(jīng)驗頻率P′對應(yīng)的標(biāo)準(zhǔn)正態(tài)分布分位數(shù)，再代入橫坐標(biāo)值計算公式“LP=-U0.01%+UP”可得到經(jīng)驗頻率P′對應(yīng)的橫坐標(biāo)值，見圖3。以其為橫坐標(biāo)，以降序排列的降雨量資料為縱坐標(biāo)，在圖表區(qū)添加系列，參照圖4設(shè)置系列名稱并框選X軸與Y軸系列值，然后回車。最后設(shè)置經(jīng)驗頻率散點圖為無線條散點圖，即完成了經(jīng)驗點據(jù)的繪制，見圖4。

圖3 經(jīng)驗點據(jù)橫、縱坐標(biāo)值計算表及其繪制

圖4 興國縣降水量頻率曲線

3.3.2 初定統(tǒng)計參數(shù)

3.3.3 選配理論頻率曲線

由圖4可知，理論頻率曲線(Cs=3Cv=0.75)的中段與經(jīng)驗頻率點據(jù)配合較好，但頭部偏于經(jīng)驗頻率點據(jù)下方，而尾部又偏于經(jīng)驗頻率點據(jù)上方，配合情況不理想。需調(diào)整參數(shù)進(jìn)行重新配線。

圖5 第一次配線計算結(jié)果與曲線繪制

b.改變參數(shù)進(jìn)行第二次配線，由第一次配線結(jié)果可知，需增大Cv值。取Cv=0.3，Cs=3Cv=0.9，參照選配理論頻率曲線第一步進(jìn)行配線，計算結(jié)果見圖6。由配線結(jié)果(圖4)可知，理論頻率曲線頭部配合較好，但尾部較經(jīng)驗點據(jù)偏低。

c.在第二次配線基礎(chǔ)上，為使尾部抬高一些與經(jīng)驗點據(jù)相配合，需增大Cs值，因此，取Cv=0.3,Cs=3.5Cv=1.05，再次計算理論頻率曲線，計算結(jié)果見圖6。由配線結(jié)果(圖4)可知，當(dāng)Cs=3.5Cv=1.05時，理論曲線與經(jīng)驗頻率點據(jù)配合較合適，讀圖得到保證率為20%、50%、75%及95%時年降雨量分別為1931.85mm、1504.31mm、1235.11mm、965.92mm(見表2)。

圖6 選配理論曲線計算結(jié)果

表2 興國縣不同保證率下年降水量對比情況

4 結(jié) 語

利用Excel軟件求解水文設(shè)計值、繪制頻率格紙和水文頻率曲線切實可行，該方法不僅能提高計算和制圖的效率，出圖效果更加整潔美觀，制圖精度相較于傳統(tǒng)的手繪曲線有所提高，為水文工作者提供了一種方便、簡潔、美觀的繪圖方法。