水文频率曲线绘制的Matlab程序设计
李清富;闫鹏飞;孙静涛;程婷
【摘 要】基于目估适线法，选用皮尔逊Ⅲ型频率曲线作为水文频率分布线型，在Matlab平台上编写三个M文件，分别实现海森频率格纸的绘制、频率曲线参数x軃、Cv的求解以及水文频率曲线的绘制。通过多次的适配，到与实测水文资料配合较好的理想曲线的参数作为设计值并估算指定的水文特征值。并以某一水文测站24年的年径流实测资料为基础，检验了方法的正确性。为水文学中相关计算以及水文频率曲线的绘制提供了一种新的实现方法。%The P-III frequency curve is selected as the hydrological frequency distribution of linear,based on visual adaptability line. Three M-files was written in Matlab platform to draw the Hazen frequency grid paper,to solve the frequency curve parameters such as x、C v ,to draw the hydrological frequency curve,respectively. With the multiple adaptation,the ideal curve that fits well with the measured hydrological data is found. So we can estimate the value of the specified hydrological characteristics by the using of the ideal curve. Annual runoff measured data of a hydrological station,24 years is used to test the correctness of the method. It provid
es a new implementation for hydrological calculations and the drawing of the hydrological hydrology frequency curve.
【期刊名称】《河南科学》
【年(卷),期】2013(000)008
【总页数】5页(P1250-1254)
【关键词】Matlab;目估适线法;海森频率格纸;皮尔逊Ⅲ型频率曲线
【作 者】李清富;闫鹏飞;孙静涛;程婷
【作者单位】郑州大学水利与环境学院,郑州,450001;郑州大学水利与环境学院,郑州,450001;郑州大学水利与环境学院,郑州,450001;郑州大学水利与环境学院,郑州,450001
【正文语种】中 文
【中图分类】P333.9;TP312
我国水文频率计算一直采用皮尔逊型Ⅲ型频率曲线作为水文分布频率曲线.随着计算机技术的发展，水文频率分析软件应运而生.如由浙江大学唐启义教授研制的一款综合性的大型统计软件DPS（Data Processing System）统计软件，通用多功能数理统计和数学模型处理软件系统，将数值计算、统计分析、模型模拟以及画线制表等功能融为一体[1]；耿鸿江[2]，王双银[3]等人在研究了基于Excel的P-Ⅲ型分布曲线的绘制方法；赵培颖[4]等人研究了Visual Basic在绘制P-Ⅲ频率曲线中的应用.本文在已有学者的研究基础上，对如何应用Matlab绘制水文频率曲线进行了研究，并给出了利用Matlab绘制水文频率曲线的M文件，具有重复可操作性.
1 计算理论
1.1 目估适线法
目估适线法[5]是我国估计水文频率曲线统计参数最主要的方法.它是以经验频率点据为基础，选配一条拟合较好的理论频率曲线，利用频率曲线与经验频率点据相配合来估计水文要素总体的统计规律.具体步骤如下：
①绘制频率格纸；②将实测资料由大到小排列，计算各项的经验频率，在频率格纸上点绘经
验点据.③选用皮尔逊Ⅲ型频率曲线作为水文频率分布线型，采用矩法或其他方法估计出频率曲线参数的初估值、C v、Cs（凭经验初选为C v的某一倍数）；④根据拟定的¯、C v、Cs，查“皮尔逊Ⅲ型频率曲线的模比系数Kp值表”，得到模比系数Kp值，并计算xp值；⑤以xp为纵坐标，P为横坐标，即可得到频率曲线，将此线画在绘有经验点据的图上；⑥分析曲线与经验点据的配合情况，通过不断调整、C v、Cs值以及各、C v、Cs所对应的曲线，直到到理想的频率曲线为止；⑦求指定的水文变量变化值.
1.2 P－Ⅲ型曲线密度函数
英国生物学家皮尔逊提出了13种非正态的分布曲线，其中第Ⅲ型曲线被引入水文计算中，用途比较广泛.P-Ⅲ型曲线[5，8-9]是一条一端有限、一端无限的不对称单峰、正偏曲线，数学上成为伽马分布，其概率密度函数为：
式中：Γ（α）—α 的伽马函数；α、β、a0—PⅢ型曲线分布的形状、尺度和位置参数 a＞0，β＞0.
显然α、β、a0这3个参数求出来以后该密度函数便可以确定.3个参数与总体的3个参数x¯、C v、Cs具有如下关系：
在水文计算中，一般需要求出指定频率P所相应的随机变量X的取值xp，也就是通过对密度曲线进行积分[5]，即
2 程序设计
程序设计流程图如图1所示所示，主要包括5个步骤：①求解总体参数；②查“皮尔逊Ⅲ型频率曲线的模比系数Kp值表”确定Kp值；③绘制海森频率坐标格纸；④绘制P-Ⅲ型频率曲线；⑤选线.
根据流程图，利用M atlab[6-7]平台编制“zuobiao.m”、“zuobiao.m”、“huitu.m”共3个M文件，分别用来求解总体参数、绘制海森频率格纸以及绘制P-Ⅲ型频率曲线，并都保存在Matlab work子目录中，以便调用.
2.1 海森频率格纸的绘制
海森频率格纸是水文计算中常用的一种坐标格纸，其横坐标的分化就是按把标准正态频率曲线拉成一条直线的原理计算出来的.这种频率格纸的纵坐标仍是普通分格，但横坐标的分格是不相等的中间分格较密，越往两端越稀，中间分格在P=50%的两端是对称的[5].利用M atlab
平台编制绘制海森频率格纸的M文件，zuobiao函数的编程和变量解释如下：
图1 水文频率曲线计算流程图Fig.1  The calculation flowchart of hydrological frequency curve
2.2 求解总体参数的M文件编写
利用Matlab平台编制绘制P-Ⅲ型频率曲线的M文件，canshu函数的编程和变量解释如下：frequency函数计算频数
2.3 绘制P－Ⅲ型频率曲线的M文件编写
利用Matlab平台编制绘制P-Ⅲ型频率曲线的M文件，huitu函数的编程和变量解释如下：
图2 频率计算成果图Fig.2  The figure of the frequency calculation results
3 实例
某测站的年径流资料[5]，如表1所示.
表1 某站实测年径流资料Tab.1  Measured annual runoff data of a station年份  年径流深/mm
  年份  年径流深/mm  年份  年径流深/mm 1952  538.3  1960  964.2  1968  732.9 1953  624.9  1961  687.3  1969  1 064.5 1954  663.2  1962  546.7  1970  606.7 1955  591.7  1963  509.9  1971  586.7 1956  557.2  1964  769.2  1972  567.4 1957  998.0  1965  615.5  1973  587.7 1958  641.5  1966  417.1  1974  709.0 1959  341.1  1967  789.3  1975  883.5
首先，将表1中的年径流资料记为变量A，调入Matlab程序，在Command Window运行“canshu”得变差系数 C v＝0.263 31”，取 C v＝0.3，并分别取 Cs为 2C v、2.5C v、3C v，查“皮尔逊Ⅲ型频率曲线的模比系数 Kp值表”得到频率 P（%）0.01、0.1、0.2、0.33、0.5、1、2、5、10、20、50、75、90、95、99 时相对应的 Kp值，并在 Matlab 中分别记为变量 Kp 1、Kp 2、Kp 3.
然后，在 Matlab 的 Command Window 运行“zuobiao.m”、“huitu.m”子程序，得到三次配线结果（可根据原散点图与配线情况的配合情况适当修改程序确定配线的次数，直到到理想曲线）如图2所示.

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。