基于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软件拟合精度最大提高一倍，若应用于工程实际，可以为工程计算提高精度和节约大量时间。