何 妹，趙華榮，姚 越，金 鑫

(桂林理工大學環(huán)境科學與工程學院，廣西 桂林 541006)

P-Ⅲ型頻率曲線適用于中國大部分地區(qū)的降雨徑流序列，成為水文計算中常用的分布頻率曲線。計算機技術(shù)在計算、繪圖等方面都具有較強大的功能，在水文與水資源方面的應用也越來越廣泛，由于P-Ⅲ型頻率曲線計算過程和圖形繪制較復雜，可以充分發(fā)揮計算機的優(yōu)勢。對水文頻率曲線的繪制和計算主要有：車國文[1]利用Excel繪制皮爾遜型曲線，其結(jié)果滿足工程要求；李清富等[2]利用MATLAB編寫M文件實現(xiàn)P-Ⅲ型曲線的繪制；趙培穎等[3]研究了Visual_Basic繪制P-Ⅲ頻率曲線的應用；許義和等[4]基于Matlab對P-Ⅲ型曲線繪制軟件的研發(fā)與應用；俞超鋒等[5]基于MATLAB對P-Ⅲ型頻率曲線參數(shù)估算的優(yōu)化；李傳博等[6]實現(xiàn)了基于C語言的P-Ⅲ型曲線繪制方法；謝子波等[7]利用基于R軟件實現(xiàn)水文頻率計算適線繪制。

上述P-Ⅲ型曲線計算和繪制過程中，大多采用商業(yè)軟件，部分軟件也不是針對統(tǒng)計計算設計的，操作過程比較復雜。如EXCEL、Visual Basic、MATLAB和C語言，在實現(xiàn)參數(shù)計算和圖形繪制時，需要較強的計算及編程能力。文獻[7]采用R軟件繪制P-Ⅲ型曲線過程中對最優(yōu)參數(shù)的確定需要人工輸入查找。為了彌補上述問題，本文采用R軟件在文獻[7]的基礎(chǔ)上，增加了自動查找最優(yōu)適線參數(shù)的功能。同時根據(jù)P-Ⅲ型頻率曲線在實際工程應用中可能會強調(diào)較大值或較小值擬合精度的情況，如對特大洪水頻率計算需要考慮較大值擬合精度，用水保證率則需考慮較小值擬合精度。因此，將水文序列劃分為較大值和較小值，增加了較大值和較小值擬合最優(yōu)參數(shù)的查找并繪制對應水文頻率曲線，可根據(jù)工程實際情況選擇使用。

1 R軟件對水文頻率計算的實現(xiàn)

1.1 R軟件運用原理及優(yōu)勢

R軟件是由奧克蘭大學的Ross Ihaka和Rontleman兩位統(tǒng)計學家基于S語言開發(fā)的一個面向?qū)ο蟮哪_本語言，該語言免費開源，并以兩者首字母“R”命名。它是一種比較適合統(tǒng)計計算的語言，可進行數(shù)據(jù)輸入、輸出、分支和循環(huán)，具備用戶可自定義函數(shù)等功能[8]。

R語言可以在命令窗口進行統(tǒng)計計算、預測分析和繪制精美圖形，實現(xiàn)數(shù)據(jù)挖掘和可視化，完成較復雜的數(shù)據(jù)分析等工作。此外，R語言有數(shù)千個的程序包，用戶能通過選擇相應的程序包快速將R應用到相關(guān)領(lǐng)域[9]。相比其他軟件繪制P-Ⅲ型水文頻率曲線，R語言更具有優(yōu)勢。R語言與其他編程語言的對比見表1。

表1 P-Ⅲ型曲線繪制常用軟件比較

P-Ⅲ型水文頻率分析計算過程包括參數(shù)估計、理論頻率、經(jīng)驗頻率的計算及圖形的繪制。以上內(nèi)容在R軟件中較容易實現(xiàn)。本文主要使用openxlsx和e1071程序包，其中openxlsx包是用于讀取excel表格中水文序列值；e1071程序包有P-Ⅲ型曲線計算過程中所需要的_gamma函數(shù)。在水文頻率計算分析中常用函數(shù)及對應R程序包見表2。

表2 水文頻率計算分析常用函數(shù)類型

由于R軟件自身的圖形編程界面不友好，一般都采用第三方圖形界面。比較常用的編輯器有Rstudio和Notepad++。本次采用帶NPPtoR插件的Notepad++編輯器編寫代碼，將編寫好的程序代碼通過F8功能鍵推送到R軟件進行調(diào)試或運行。

P-Ⅲ型曲線在中國水文頻率計算分析中應用較廣泛，該線性的不確定性影響是較低的[10]。P-Ⅲ型水文頻率計算中常用的參數(shù)估計方法有矩法、適線法、圖解法、權(quán)函數(shù)法、概率權(quán)重法和極大似然法等[11]。本文采用適線法推求水文參數(shù)，適線法可分為目估適線法和優(yōu)化適線法，目估適線受到人為因素的影響，較難找到最佳擬合曲線。優(yōu)化適線法是通過計算縱向離差平方和查找到最佳擬合曲線。水文頻率分析計算過程見圖1。

圖1 水文P-Ⅲ曲線分析計算流程

1.2 水文頻率計算原理

P-Ⅲ型曲線的概率密度函數(shù)為：

(1)

其中Γ(α)為α的伽馬函數(shù)；α、β、a0分別代表P-Ⅲ型的形狀、尺度和位置參數(shù)，并且α、β均大于零[12]。α、β、a0通過均值(Ex)、變差系數(shù)(Cv)和偏態(tài)系數(shù)(Cs)3個統(tǒng)計參數(shù)關(guān)系可求出，關(guān)系如下：

(2)

(3)

(4)

其中根據(jù)矩法計算參數(shù)估計：

(5)

(6)

(7)

使用R軟件mean函數(shù)計算均值Ex，sd函數(shù)計算標準差σ，skewness函數(shù)計算偏斜度Sc，計算出以上參數(shù)后，該概率密度函數(shù)也隨之確定。由式(8)計算出對應的一組經(jīng)驗頻率值：

(8)

式中m——水文序列從大到小排列的序數(shù)；n——實測序列年數(shù)；p——經(jīng)驗頻率。

在P-Ⅲ型頻率計算中，需要求出某一個數(shù)值x所對應的概率P，故將概率密度函數(shù)進行積分得出大于等于xp頻率P，見式(9)：

P=P(x≥xp)=

(9)

通過以上參數(shù)可以變換得出式(10)：

xp=Ex(1+Cvφ)

(10)

(11)

根據(jù)上述參數(shù)之間的關(guān)系可求出離均勢數(shù)φ和水文變量值xp，從而可以計算出理論頻率P。

2 R軟件在P-Ⅲ型水文頻率計算的應用

2.1 繪制海森概率格紙

繪制水文頻率P-Ⅲ型曲線之前，先繪制機率格紙，即海森概率格紙。海森概率格紙的橫坐標數(shù)值是非均勻分布的數(shù)字坐標，而縱坐標數(shù)值則是均勻排列的數(shù)字格式，其中橫坐標與頻率值(下側(cè)概率)的標準正態(tài)分布分位數(shù)P=50%有關(guān)，p=0.01%時橫坐標值為0。

海森概率格紙繪制方法如下：首先，將(0.01，0.5，1，5，10，20，30，40，50，60，70，80，90，99.99，99.99)賦值到海森概率格紙的橫坐標，利用qnorm函數(shù)求出對應的分位數(shù)值，并調(diào)用axis函數(shù)將所求的分位數(shù)替換橫坐標刻度單位值；調(diào)用qgamma函數(shù)計算對應水文頻率值作為縱坐標值，根據(jù)經(jīng)驗頻率值P用qnorm函數(shù)計算對應x軸數(shù)值；最后調(diào)用abline函數(shù)完成海森概率格紙的繪制。

R軟件繪制機率格紙的相關(guān)函數(shù)[7]包括伽馬函數(shù)分布函數(shù)(pgamma)、分位數(shù)函數(shù)(qgamma)、正態(tài)分布函數(shù)(qnorm)。函數(shù)調(diào)用方式為pgamma(q,hape,rate,scale=1/rate)；qgamma(p,shape,rate,scale=1/rate)；qnorm(p,mean,sd)。R軟件中調(diào)用pgamma、qgamma函數(shù)的程序語句為：P=1-pgamma(x-a0，a，β)；x=qgamma(1-P，a，β)+a0。

2.2 水文頻率計算適線過程

由于目估適線法受到人為因素的影響。本文采用縱向離差平方和作為判斷適線優(yōu)劣的標準，通過計算機自動“捕捉”最佳適線參數(shù)，使適線更加準確。依照設計洪水頻率曲線適線離差平方和最小準則，當縱向離差平方和最小，即∑(xi-xp)2=Δmin適線結(jié)果最佳，反之，較差。

采用矩法計算的參數(shù)估計值時，Cs偏差較大，通常將矩法所算出的參數(shù)值作為適線調(diào)整初值。適線過程不計算具體的Cs值，而是計算Cs為Cv某一倍數(shù)變化。在完成繪制海森概率格紙后，用qgamma函數(shù)和pgamma函數(shù)計算全序在Cs與Cv不同倍數(shù)下的理論頻率和對應的水文變量值(即縱坐標)。通過qnorm將序列轉(zhuǎn)變?yōu)檎龖B(tài)函數(shù)的分位數(shù)(即橫坐標)。最后，分別計算全序、較大值及較小值在Cv和Cs不同倍數(shù)下的縱向離差平方和，求出最小縱向離差平方和及對應參數(shù)作為最佳參數(shù)組合，再調(diào)用lines函數(shù)在概率格紙上繪制頻率曲線。

3 實例計算

3.1 驗證過程

本次實例采用文獻[13]中某站1949—2000年的實測水文序列驗證，具體步驟如下：①利用R軟件讀取Excel表格中的年徑流量序列，使用相關(guān)函數(shù)計算序列長度(n)，均值(Ex)、變差系數(shù)(Cv)和偏態(tài)系數(shù)(Cs)等參數(shù)；②繪制海森概率格紙，根據(jù)式(8)計算經(jīng)驗頻率值并點繪在概率格紙上；③將年徑流量分成較大值序列與較小值序列；④調(diào)用for函數(shù)計算P-Ⅲ型曲線參數(shù)，將矩法估算Cv和Cs值作為初值，將Cv設定范圍從Cv-0.1到Cv+0.1；Cs設定范圍在Cv～4Cv之間，兩者均按照0.001的步長增加；⑤確定最佳適線參數(shù)。根據(jù)經(jīng)驗點與頻率曲線曲線的縱向離差平方和，查找全序、較大值和較小值縱向離差平方和的最小值，保存最小縱向離差平方和對應參數(shù)，退出循環(huán)。⑥圖形繪制。根據(jù)最小縱向離差平方和對應參數(shù)調(diào)用lines函數(shù)繪制P-Ⅲ曲線(圖2)；⑦在R窗口顯示部分水文頻率所對應的水文變量值(表3)。

圖2 水文P-Ⅲ曲線適線結(jié)果

表3 部分水文頻率下的水文變量

3.2 驗證結(jié)果

本實例中經(jīng)驗頻率、矩法和適線法結(jié)果見圖2，其中使用矩法計算的Cv和Cs值分別為0.358和0.519，其縱向離差平方和為11 138，偏差較大，精確度不高，一般不建議用該法擬合水文頻率曲線。適線法對全序、大值及小值進行擬合，結(jié)果為：較大值最佳適線參數(shù)為Cv=0.394、Cs=1.383×Cv，其最小縱向離差平方和為2 989；較小值最佳適線參數(shù)為Cv=0.383、Cs=2.368×Cv，最小縱向離差平方和為1 559。全序列最佳適線參數(shù)為Cv=0.383、Cs=1.992×Cv，其最小縱向離差平方和為5 776，全序列適線比文獻[7]的縱向離差平方和少186表明擬合結(jié)果更佳。

為了檢驗R優(yōu)化成果，對部分文獻的實例數(shù)據(jù)用R軟件再次進行參數(shù)擬合，并與原文獻結(jié)果進行比較。由于原文獻中沒有給出縱向離差平方和，因此采用原參數(shù)通過本文的方法進行計算，結(jié)果見表4。采用本文適線方法所得的縱向離差平方和比原參數(shù)計算結(jié)果都要小一個數(shù)量級，說明本文所采用方法可以大大提高擬合精度。

表4 不同軟件適線結(jié)果與本文適線結(jié)果比較

4 結(jié)論

a) R軟件結(jié)合Notepad++具有操作簡單，適用于小規(guī)模的軟件開發(fā)等功能。其軟件包統(tǒng)計函數(shù)多，編程及繪圖方便，計算精度較高，實用性強，易于操作，便于實際應用推廣。

b) 通過實例驗證結(jié)果可知R軟件不僅可以精準地計算出水文序列的特征參數(shù)，繪制海森概率格紙，而且能快速獲取最小縱向離差平方和及對應參數(shù)，并根據(jù)最佳適線參數(shù)組合(Cv=0.383、Cs=1.992×Cv)繪制P-Ⅲ曲線的最優(yōu)曲線，該法相比同軟件繪制所得的縱向離差平方和少了186，達到優(yōu)化效果。

c) 通過不同軟件適線結(jié)果比較可知R軟件擬合精度最大提高一倍，若應用于工程實際，可以為工程計算提高精度和節(jié)約大量時間。