Excel在求解水文频率曲线中的应用

李林波 魏 为 黄国芳

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

频率曲线是推求设计降雨量、设计洪水等水文设计值的重要曲线[1]。该曲线一般绘制在频率格纸(海森几率格纸)上，但由于频率格纸横坐标轴不均匀，传统绘制方法多是在已绘好的频率格纸上徒手绘制，不仅效率低，而且精度不高，存在人为误差[2]。加之频率曲线的求解过程计算量较大，给工作人员带来诸多不便[3]。

前人在该领域的研究多为介绍频率曲线求解原理，或讲述求解步骤[4]，基于Excel软件对此类问题进行的研究多止于绘制频率格纸，鲜有利用Excel软件求解频率曲线的研究，因此本文将着重解决这一问题，对其进行研究。

1 基本概念与原理

1.1 水文频率曲线

水文频率曲线是水文要素值与其发生频率之间的关系曲线，当设计标准[设计频率(P)或设计重现期(T)]已知时，可在该曲线上查得相应的水文设计值(xp)，其意思是任何水文要素值x超过或等于xp的频率为P，即

P(x≥xp)=P

(1)

水文频率曲线有多种形式。《水利水电工程设计洪水计算规范》(SL 44—2006)规定应使用皮尔逊Ⅲ型曲线(P-Ⅲ型曲线)，常见的有降雨量频率曲线、洪峰流量频率曲线[5-6]。

皮尔逊Ⅲ型曲线是一条一端有限一端无限的不对称单峰、正偏曲线[7]，数学上常称伽玛分布，其概率密度函数为

(2)

式中：Γ(α)为α的伽玛函数；α、β、a0为皮尔逊Ⅲ型分布的形状尺度和位置未知参数，α>0，β>0。

(3)

(4)

(5)

水文计算中，一般需要求出随机变量取值xp所相应的指定频率P，也就是对密度曲线进行积分，即

(6)

实际计算时通常变换成如下形式进行积分

(7)

(8)

式中：Φ为离均系数。

式(7)中被积函数只含有一个待定参数Cs，只需要假定一Cs值，便可通过积分求出P与Φ之间的关系。

对于若干给定的Cs值，当Cs与Cv成倍数关系时，Φ和P有对应数值表，可查阅“皮尔逊Ⅲ型频率曲线的离均系数Φp值表”。

(9)

对于若干给定的Cs值，当Cs等于Cv的一定倍数时，P-Ⅲ型频率曲线的模比系数Kp也已制成表格，可查阅“皮尔逊Ⅲ型频率曲线的模比系数Kp值表”。

(10)

1.2 频率格纸

频率格纸又叫海森几率格纸，它是一种特殊的坐标系统，其纵坐标为均匀分格或对数分格的常规数学坐标，横坐标与频率值(下侧概率)的标准正态分布分位数有关，两端分格较稀而中间较密[8]。

由于标准正态分布分位数在P=50%处为零，而海森几率格纸在P=0.01%时的横坐标值为零，因此海森几率格纸横坐标值可采用如下计算公式进行计算：

LP=-U0.01%+UP

(11)

式中：LP为频率P对应的横坐标值；UP为频率P对应的标准正态分布分位数；U0.01%为频率P=0.01%对应的标准正态分布分位数。

2 水文频率曲线求解方法与步骤

2.1 求解方法

(12)

(13)

2.2 求解步骤

a.点绘经验点据。如果用等概率公式P(x≥xi)=m/n的经验分布曲线估计总体分布曲线，存在不合理现象。当m=n时，最末项的频率为100%，即样本最小值为总体中的最小值，不符合事实。因此，水文上用期望值公式估计频率：

(14)

式中：m为按降水量递减顺序序列排列的项数；n为序列总项数。

利用式(14)求出经验频率后，以经验频率为横坐标、变量值为纵坐标，点绘经验点据。

d.选择一条与经验点据配合最佳的曲线作为采用曲线。该曲线的参数看作总体参数的估计值。

2.3 计算结果与表示

选取与经验点据配合最佳的理论频率曲线为所求频率曲线，通过读图可以得到不同频率P对应的xp值。

由于频率比较抽象，为便于理解，常采用重现期表示。重现期是指在许多试验中，某一事件重复出现的时间间隔的平均数，它是频率的倒数。在工程水文中，重现期用T表示，一般以年为单位[10]：

(15)

如当暴雨或洪水频率为1%时，即重现期T=100年，称此为百年一遇的暴雨或洪水。

3 计算实例

表1 兴国县2000—2017年年降水量 单位：mm

3.1 计算经验频率曲线统计参数

a.在表格第一行依次输入年份、年降水量R、序号等项目，见图1。

b.做好表头后，首先将A、B栏已知年份及其相应的年降雨量输入计算表格内，然后在C栏填入序号1～18。之后在D栏对B栏的降水资料进行由大到小的排序。

e.计算“Ki-1和(Ki-1)2。在F2处输入“=E2-1”，回车；在G2处输入“=F2^2”，回车。然后分别将光标移至F2和G2的右下角，当出现十字标识后，向下拖动至需要计算数据行。

图1 经验频率曲线统计参数计算表

3.2 绘制频率格纸

频率格纸的横向网格线为均匀分布，可直接由Excel软件的图表功能自动生成，而纵向网格线坐标值是按照标准正态频率曲线拉成一条直线的原理计算得出的，中间较密，两端较稀，不能由Excel软件的图表功能自动生成。

因此，绘制频率格纸的关键是绘制频率格纸的纵向网格线，可以先用Excel软件中的内置函数NORMSINV(P)将不同频率转化为相应的横坐标值，然后插入散点图。函数NORMSINV(P)为返回标准正态累积分布函数的反函数值，其计算结果的精度可达到±3×10-7，该函数的详细说明和用法可参考Excel软件的帮助。

绘制频率格纸的步骤具体如下：

a.新建Excel工作簿，在第一行依次输入“频率P”“P对应的标准正态分布分位数”“X轴”“Y轴”，见图2。

b.计算不同频率P对应的横坐标Lp值。在A栏中输入频率值0.01～99.99。在A2、A3处分别输入“0.01”，在A4、A5处分别输入“0.02”……依此类推，在A栏后续单元格中输入纵向网络线对应频率值，直至“99.99”止。

c.计算频率P对应的标准正态分布分位数Up。在B2单元格中输入“=NORMSINV(A2%)”，回车，然后利用Excel填充柄向下填充至P=99.99%对应单元格，就得到了频率P对应的标准正态分布分位数Up。

c.计算频率格纸横坐标Lp值。由于海森几率格纸在P=0.01%时的横坐标值为零，根据其横坐标值计算公式LP=-U0.01%+UP，在C2处输入“=-$B$2+B2”，回车，然后利用Excel填充柄向下填充，就得到了频率P对应的横坐标Lp值。见图2中C列。

d.设置纵坐标轴边界。由降水资料可知最大年降水量Rmax为2265.8mm，为方便频率曲线展布，可将Y轴最大值取大一点，本次绘制取Y轴范围为0～4000。在D列依次输入“0、4000、4000、0、0、4000、4000、0、0、4000、4000…”，直到P=99.99%对应单元格为止。见图2。

e.绘制网格线。ⓐ插入散点图，选中C、D两列，以C、D两列数据分别作为散点图的横坐标值和纵坐标值，选择插入图表命令，在“所有图表”中选择“XY散点图”，然后选择“带直线的散点图”，回车，见图2；ⓑ设置纵向网格线及横坐标刻度，删除原纵向网络线，选中所绘散点图，将线条改为实线、黑色、0.5磅宽；将横坐标边界设置为0～7.438，并删除原横坐标刻度，选中刚绘制的散点图，右击鼠标，在弹出的选项中选择“添加数据标签”，选中出现的数据标签并设置为只显示“X值”，拖动标签位置并修改标签内容为相应纵线对应的频率P值，最后删除多余标签便完成了频率格纸纵线的绘制；ⓒ设置横向网格线及横坐标刻度，选中横坐标刻度，右击鼠标，选择“设置坐标轴格式”，设置边界范围为0～4000，间距200，然后选中横向网格线，从鼠标右键菜单中选择“设置网格线格式”，修改线条为实线、黑色、0.5磅。至此便完成了海森几率格纸的绘制。

图2 频率格纸横、纵坐标计算表与散点图绘制

3.3 适线法配线

3.3.1 点绘经验点据

利用公式“NORMSINV(P′)”计算得到经验频率P′对应的标准正态分布分位数，再代入横坐标值计算公式“LP=-U0.01%+UP”可得到经验频率P′对应的横坐标值，见图3。以其为横坐标，以降序排列的降雨量资料为纵坐标，在图表区添加系列，参照图4设置系列名称并框选X轴与Y轴系列值，然后回车。最后设置经验频率散点图为无线条散点图，即完成了经验点据的绘制，见图4。

图3 经验点据横、纵坐标值计算表及其绘制

图4 兴国县降水量频率曲线

3.3.2 初定统计参数

3.3.3 选配理论频率曲线

由图4可知，理论频率曲线(Cs=3Cv=0.75)的中段与经验频率点据配合较好，但头部偏于经验频率点据下方，而尾部又偏于经验频率点据上方，配合情况不理想。需调整参数进行重新配线。

图5 第一次配线计算结果与曲线绘制

b.改变参数进行第二次配线，由第一次配线结果可知，需增大Cv值。取Cv=0.3，Cs=3Cv=0.9，参照选配理论频率曲线第一步进行配线，计算结果见图6。由配线结果(图4)可知，理论频率曲线头部配合较好，但尾部较经验点据偏低。

c.在第二次配线基础上，为使尾部抬高一些与经验点据相配合，需增大Cs值，因此，取Cv=0.3,Cs=3.5Cv=1.05，再次计算理论频率曲线，计算结果见图6。由配线结果(图4)可知，当Cs=3.5Cv=1.05时，理论曲线与经验频率点据配合较合适，读图得到保证率为20%、50%、75%及95%时年降雨量分别为1931.85mm、1504.31mm、1235.11mm、965.92mm(见表2)。

图6 选配理论曲线计算结果

表2 兴国县不同保证率下年降水量对比情况

4 结 语

利用Excel软件求解水文设计值、绘制频率格纸和水文频率曲线切实可行，该方法不仅能提高计算和制图的效率，出图效果更加整洁美观，制图精度相较于传统的手绘曲线有所提高，为水文工作者提供了一种方便、简洁、美观的绘图方法。