基于R軟件優化繪制皮爾遜Ⅲ型頻率曲線

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

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

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

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

1 R軟件對水文頻率計算的實現

1.1 R軟件運用原理及優勢

R軟件是由奧克蘭大學的Ross Ihaka和Rontleman兩位統計學家基于S語言開發的一個面向對象的腳本語言，該語言免費開源，并以兩者首字母“R”命名。它是一種比較適合統計計算的語言，可進行數據輸入、輸出、分支和循環，具備用戶可自定義函數等功能[8]。

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

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

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

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

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

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

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

1.2 水文頻率計算原理

P-Ⅲ型曲線的概率密度函數為：

(1)

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

(2)

(3)

(4)

其中根據矩法計算參數估計：

(5)

(6)

(7)

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

(8)

式中m——水文序列從大到小排列的序數；n——實測序列年數；p——經驗頻率。

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

P=P(x≥xp)=

(9)

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

xp=Ex(1+Cvφ)

(10)

(11)

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

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

2.1 繪制海森概率格紙

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

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

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

2.2 水文頻率計算適線過程

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

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

3 實例計算

3.1 驗證過程

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

圖2 水文P-Ⅲ曲線適線結果

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

3.2 驗證結果

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

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

表4 不同軟件適線結果與本文適線結果比較

4 結論

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

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

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