第31卷第8期2013年8月
文章编号：1004-3918（2013）08-1250-05
河南科学
HENAN SCIENCE
Vol.31No.8Aug. 2013
水文频率曲线绘制的M ａｔｌａｂ程序设计
李清富，
摘
闫鹏飞，孙静涛，程婷
（郑州大学水利与环境学院，郑州４５０００１）
要：基于目估适线法，选用皮尔逊Ⅲ型频率曲线作为水文频率分布线型，在M ａｔｌａｂ平台上编写三个Ｍ文件，分
軃、C v 的求解以及水文频率曲线的绘制. 通过多次的适配，找到与实测水别实现海森频率格纸的绘制、频率曲线参数x
文资料配合较好的理想曲线的参数作为设计值并估算指定的水文特征值. 并以某一水文测站２４年的年径流实测资料为基础，检验了方法的正确性. 为水文学中相关计算以及水文频率曲线的绘制提供了一种新的实现方法. 关键词：M ａｔｌａｂ；目估适线法；海森频率格纸；皮尔逊Ⅲ型频率曲线中图分类号：Ｐ３３３．９；ＴＰ３１２
文献标识码：A
ＤｅｖｅｌｏｐｍｅｎｔｏｆM ａｔｌａｂＰｒｏｇｒａｍｆｏｒＨｙｄｒｏｌｏｇｉｃａｌＦｒｅｑｕｅｎｃｙＣｕｒｖｅＤｒａｗｉｎｇ
Li Qingfu ，Ｙan Ｐｅｎｇｆｅｉ，
Ｓun Ｊｉｎｇｔａｏ，
Ｃheng Ｔｉｎｇ
（ＳｃｈｏｏｌｏｆＥｎｖｉｒｏｎｍｅｎｔａｎｄＷａｔｅｒＣｏｎｓｅｒｖａｎｃｙＥｎｇｉｎｅｅｒｉｎｇ，ＺｈｅｎｇｚｈｏｕＵｎｉｖｅｒｓｉｔｙ，Ｚｈｅｎｇｚｈｏｕ４５０００１，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：The Ｐ－ＩＩＩｆｒｅｑｕｅｎｃｙｃｕｒｖｅｉｓｓｅｌｅｃｔｅｄａｓｔｈｅｈｙｄｒｏｌｏｇｉｃａｌｆｒｅｑｕｅｎｃｙｄｉｓｔｒｉｂｕｔｉｏｎｏｆｌｉｎｅａｒ，ｂａｓｅｄｏｎｖｉｓｕａｌ
ａｄａｐｔａｂｉｌｉｔｙｌｉｎｅ．ＴｈｒｅｅＭ－ｆｉｌｅｓｗａｓｗｒｉｔｔｅｎｉｎＭａｔｌａｂｐｌａｔｆｏｒｍｔｏｄｒａｗｔｈｅＨａｚｅｎｆｒｅｑｕｅｎｃｙｇｒｉｄｐａｐｅｒ，ｔｏｓｏｌｖｅ軃、C v ，ｔｏｄｒａｗｔｈｅｈｙｄｒｏｌｏｇｉｃａｌｆｒｅｑｕｅｎｃｙｃｕｒｖｅ，ｒｅｓｐｅｃｔｉｖｅｌｙ．Ｗｉｔｈｔｈｅｔｈｅｆｒｅｑｕｅｎｃｙｃｕｒｖｅｐａｒａｍｅｔｅｒｓｓｕｃｈａｓx
ｍｕｌｔｉｐｌｅａｄａｐｔａｔｉｏｎ，ｔｈｅｉｄｅａｌｃｕｒｖｅｔｈａｔｆｉｔｓｗｅｌｌｗｉｔｈｔｈｅｍｅａｓｕｒｅｄｈｙｄｒｏｌｏｇｉｃａｌｄａｔａｉｓｆｏｕｎｄ．Ｓｏｗｅｃａｎｅｓｔｉｍａｔｅ
ｔｈｅｖａｌｕｅｏｆｔｈｅｓｐｅｃｉｆｉｅｄｈｙｄｒｏｌｏｇｉｃａｌｃｈａｒａｃｔｅｒｉｓｔｉｃｓｂｙｔｈｅｕｓｉｎｇｏｆｔｈｅｉｄｅａｌｃｕｒｖｅ．Ａｎｎｕａｌｒｕｎｏｆｆｍｅａｓｕｒｅｄｄａｔａｏｆａｈｙｄｒｏｌｏｇｉｃａｌｓｔａｔｉｏｎ，２４ｙｅａｒｓｉｓｕｓｅｄｔｏｔｅｓｔｔｈｅｃｏｒｒｅｃｔｎｅｓｓｏｆｔｈｅｍｅｔｈｏｄ．Ｉｔｐｒｏｖｉｄｅｓａｎｅｗｉｍｐｌｅｍｅｎｔａｔｉｏｎｆｏｒｈｙｄｒｏｌｏｇｉｃａｌｃａｌｃｕｌａｔｉｏｎｓａｎｄｔｈｅｄｒａｗｉｎｇｏｆｔｈｅｈｙｄｒｏｌｏｇｉｃａｌｈｙｄｒｏｌｏｇｙｆｒｅｑｕｅｎｃｙｃｕｒｖｅ．
Ｋｅｙｗｏｒｄｓ：M ａｔｌａｂ；ｔａｒｇｅｔａｄａｐｔａｂｉｌｉｔｙｌｉｎｅ；ｔｈｅＨａｚｅｎｆｒｅｑｕｅｎｃｙｇｒｉｄｐａｐｅｒ；Ｐ－ＩＩＩｆｒｅｑｕｅｎｃｙｃｕｒｖｅ我国水文频率计算一直采用皮尔逊型Ⅲ型频率曲线作为水文分布频率曲线. 随着计算机技术的发展，水文频率分析软件应运而生. 如由浙江大学唐启义教授研制的一款综合性的大型统计软件ＤＰＳ（ＤａｔａＰｒｏｃｅｓｓｉｎｇＳｙｓｔｅｍ）统计软件，通用多功能数理统计和数学模型处理软件系统，将数值计算、统计分析、模型模拟以及画线制表等功能融为一体[1]；耿鸿江[2]，王双银[3]等人在研究了基于Ｅｘｃｅｌ的Ｐ－Ⅲ型分布曲线的绘制方法；赵培颖[4]等人研究了ＶｉｓｕａｌＢａｓｉｃ在绘制Ｐ－Ⅲ频率曲线中的应用. 本文在已有学者的研究基础上，对如何应用Ｍａｔｌａｂ绘制水文频率曲线进行了研究，并给出了利用Ｍａｔｌａｂ绘制水文频率曲线的Ｍ文件，具有重复可操作性.
１计算理论
目估适线法
目估适线法[5]是我国估计水文频率曲线统计参数最主要的方法. 它是以经验频率点据为基础，选配一条拟合较好的理论频率曲线，利用频率曲线与经验频率点据相配合来估计水文要素总体的统计规律. 具体步骤如下：
收稿日期：2013-01-10
基金项目：国家自然科学基金项目（５０９７８２３４）
作者简介：李清富（１９６６－），男，河南林州人，教授，博士，从事交通运输工程领域的研究
闫鹏飞（１９９１－），男，河南商丘人，主要从事道路桥梁方面的研究．
１．１
2013年8月李清富，等：水文频率曲线绘制的M ａｔｌａｂ程序设计
-１２５１-
①绘制频率格纸；②将实测资料由大到小排列，计算各项的经验频率，在频率格纸上点绘经验点据. ③选軃、C v 、C s 用皮尔逊Ⅲ型频率曲线作为水文频率分布线型，采用矩法或其他方法估计出频率曲线参数的初估值x 軃、C v 、C s ，查（凭经验初选为C v 的某一倍数）；④根据拟定的x “皮尔逊Ⅲ型频率曲线的模比系数K p 值表”，得到模比系数K p 值，并计算x p 值；⑤以x p 为纵坐标，P 为横坐标，即可得到频率曲线，将此线画在绘有经验点据的軃、C v 、C s 值以及各x 軃、C v 、C s 所对应的曲线，直到找到图上；⑥分析曲线与经验点据的配合情况，通过不断调整x 理想的频率曲线为止；⑦求指定的水文变量变化值.
１．２Ｐ－Ⅲ型曲线密度函数
英国生物学家皮尔逊提出了１３种非正态的分布曲线，其中第Ⅲ型曲线被引入水文计算中，用途比较广泛. Ｐ－Ⅲ型曲线[5，8鄄9]是一条一端有限、一端无限的不对称单峰、正偏曲线，数学上成为伽马分布，其概率密度函数为：
βa -1-β（x-a ）
x-a 0） e ，（ｆx ）=Γ（α）
α
式中：Γ（α）—α的伽马函数；α、β、a 0—ＰⅢ型曲线分布的形状、尺度和位置参数ａ＞０，β＞０.
軃、C v 、C s 具有显然α、β、a 0这３个参数求出来以后该密度函数便可以确定. 3个参数与总体的3个参数x 如下关系：
α=42；
C s
β=2；軃C v C s x
軃a 0=x （1-2C v
. C s
在水文计算中，一般需要求出指定频率Ｐ所相应的随机变量X 的取值x p ，也就是通过对密度曲线进行积分[5]，即
β
P =P （x ≥x p ）=
Γ（α）
α
乙
x p
∞
（x-a 0） e
a -1-β（x-a 0）
d x .
２程序设计
程序设计流程图如图１所示所示，主要包括５个步骤：①求
解总体参数；②查“皮尔逊Ⅲ型频率曲线的模比系数K p 值表”确定K p 值；③绘制海森频率坐标格纸；④绘制Ｐ－Ⅲ型频率曲线；⑤选线.
根据流程图，利用M ａｔｌａｂ[6鄄7]平台编制“ｚｕｏｂｉａｏ．ｍ”、“ｚｕｏｂｉａｏ．ｍ”、
“ｈｕｉｔｕ．ｍ”共3个Ｍ文件，分别用来求解总体参数、绘制海森频率格纸以及绘制Ｐ－Ⅲ型频率曲线，并都保存在Ｍａｔｌａｂｗｏｒｋ子目录中，以便调用. ２．１海森频率格纸的绘制
海森频率格纸是水文计算中常用的一种坐标格纸，其横坐标的分化就是按把标准正态频率曲线拉成一条直线的原理计算出
图１水文频率曲线计算流程图来的. 这种频率格纸的纵坐标仍是普通分格，但横坐标的分格是
Ｆｉｇ．１Ｔｈｅｃａｌｃｕｌａｔｉｏｎｆｌｏｗｃｈａｒｔｏｆｈｙｄｒｏｌｏｇｉｃａｌ
不相等的中间分格较密，越往两端越稀，中间分格在Ｐ＝５０％的两
ｆｒｅｑｕｅｎｃｙｃｕｒｖｅ
[5]
端是对称的.
利用M ａｔｌａｂ平台编制绘制海森频率格纸的Ｍ文件，ｚｕｏｂｉａｏ函数的编程和变量解释如下：ｆｕｎｃｔｉｏｎｚｕｏｂｉａｏ％定义坐标函数
ｘ＝［０．０１：０．０１：．１，０．２：．１：１，１．２：．２：２，３：２０，２２：２：８０，８１：９８，９８．２：．２：９９，９９．１：．１：９９．９］；％横坐标的值ｙ＝ｎｏｒｍｉｎｖ（ｘ/１００，０，１）；ｙ＝ｙ－ｙ（１）；
-１２５２-ｆｏｒｉ＝１：ｓｉｚｅ（ｙ，２）
ｌｉｎｅ（［ｙ（ｉ），ｙ（ｉ）］，［０２０００］，′ｃｏｌｏｒ′，′ｋ′）；ｅｎｄ
ｚ＝０：１００：２０００；
河南科学
第31卷第8期
ｆｏｒｉ＝１：２１；％纵坐标的间距设置，可以自己设置ｌｉｎｅ（［０，ｙ（ｅｎｄ）］，［ｚ（ｉ）ｚ（ｉ）］，′ｃｏｌｏｒ′，′ｋ′）；
ｅｎｄ
ｈ＝ｆｉｎｄｏｂｊ（′ｔｙｐｅ′，′ａｘｅｓ′）；ｓｅｔ（ｈ，′ｘｔｉｃｋ′，［］，′ｘｌｉｍ′，［０ｙ（ｅｎｄ）］，′ｙｌｉｍ′，［０２０００］）；％同时去掉ｘ轴和ｙ轴的刻度
ｘｘ＝［０．０１０．０５０．１０．５１．０５１０１５２０３０４０５０６０７０８０８５９０９５９９９９．５９９．９］；％横坐标的设置，可以自己设置ｙｙ＝ｎｏｒｍｉｎｖ（ｘｘ/１００，０，１）；ｙｙ＝ｙｙ－ｙｙ（１）；ｓｅｔ（ｈ，′Ｘｔｉｃｋ′，ｙｙ）；ｓｅｔ（ｈ，′ＸｔｉｃｋｌａｂｅＬ′，｛′０．０１′，′０．０５′，′０．１′，′０．５′，′１．０′，′５′，′１０′，′１５′，′２０′，′３０′，′４０′，′５０′，′６０′，′７０′，′８０′，′８５′，′９０′，′９５′，．．．
′９９′，′９９．５′，′９９．９′｝）；
ｓｅｔ（ｈ，′Ｙｔｉｃｋ′，ｚ）；ｘｌａｂｅｌ（′Ｐ（％）′）；％对ｘ轴进行标注ｙｌａｂｅｌ（′Ｑｍ（ｍ/ｓ）′）；％对ｙ轴进行标注，可根据水文特征值设置２．２求解总体参数的Ｍ文件编写
利用M ａｔｌａｂ平台编制绘制Ｐ－Ⅲ型频率曲线的Ｍ文件，ｃａｎｓｈｕ函数的编程和变量解释如下：％％参数计算ｆｏｒｍａｔｓｈｏｒｔａ＝ｉｎｐｕｔ（′水文学数据′）；Ａ＝ｓｏｒｔ（ａ，′ｄｅｓｃｅｎｄ′）；％对水文学数据进行排序ｎ＝ｌｅｎｇｔｈ（Ａ）；％计算数据的个数ｑ＝１：ｎ；ｑ＝ｑ．′；Ｐ＝ｑ．（/ｎ＋１）；％计算经验频率ｐ＝ｎｏｒｍｉｎｖ（Ｐ，０，１）；％数据Ｐ所对应的正态分布的Ｘ值Ｐ１＝［０．０１；０．１；０．２；０．３３；０．５；１；２；５；１０；２０；５０；７５；９０；９５；９９；］；ｐ１＝ｎｏｒｍｉｎｖ（Ｐ１．/１００，０，１）；％数据ｐ１所对应的正态分布的Ｘ值Ｍ＝ｎｏｒｍｉｎｖ（０．０００１，０，１）；ｐ２＝－Ｍ＋ｐ１；％转换坐标值（理论频率）ｐ３＝－Ｍ＋ｐ；％转换坐标值（经验频率）Ｑ＝ｍｅａｎ（Ａ）；％求平均值Ｋ＝Ａ．/Ｑ；％求模比系数Ｋ１＝Ｋ－１；Ｋ２＝Ｋ１．＾２；Ｃｖ＝ｓｑｒｔ（ｓｕｍ（Ｋ２））/ｓｑｒｔ（２４－１）；％计算变差系数ｄｉｓｐ（［′平均值Ｑ＝′，ｎｕｍ２ｓｔｒ（Ｑ）］）ｄｉｓｐ（［′变差系数Ｃｖ＝′，ｎｕｍ２ｓｔｒ（Ｃｖ）］）２．３
绘制Ｐ－Ⅲ型频率曲线的Ｍ文件编写
利用M ａｔｌａｂ平台编制绘制Ｐ－Ⅲ型频率曲线的Ｍ文件，ｈｕｉｔｕ函数的编程和变量解释如下：％％绘制图像ｐｌｏｔ（ｐ３，Ａ，′＊ｒ′）；％绘制经验频率的散点
2013年8月李清富，等：水文频率曲线绘制的M ａｔｌａｂ程序设计
-１２５３-
ｈｏｌｄｏｎ；ｇｒｉｄｏｎ；
Ｋｐ＝ｉｎｐｕｔ（′输入查表参数′）；Ｘｐ＝Ｑ＊Ｋｐ；％Ｘｐ的理论计算值ｐｌｏｔ（ｐ２，Ｘｐ，′ｋ′）；％第一次配线ｈｏｌｄｏｎ；ｇｒｉｄｏｎ；
Ｋｐ＝ｉｎｐｕｔ（′输入查表参数′）；Ｘｐ＝Ｑ＊Ｋｐ；％Ｘｐ的理论计算ｐｌｏｔ（ｐ２，Ｘｐ，′ｂ′）；％第二次配线ｈｏｌｄｏｎ；ｇｒｉｄｏｎ；
Ｋｐ＝ｉｎｐｕｔ（′输入查表参数′）；
Ｘｐ＝Ｑ＊Ｋｐ；Ｘｐ的理论计算值ｐｌｏｔ（ｐ２，Ｘｐ，′ｇ′）；％第三次配线ｌｅｇｅｎｄ（′原数据散点′，′第一次配线′，′第二次配线′，′第三次配线′）％图例标注
R /m m
P /%
图２
Ｆｉｇ．２
频率计算成果图
Ｔｈｅｆｉｇｕｒｅｏｆｔｈｅｆｒｅｑｕｅｎｃｙｃａｌｃｕｌａｔｉｏｎｒｅｓｕｌｔｓ
３实例
某测站的年径流资料[5]，如表１所示.
表１
Ｔａｂ. １
年份１９５２１９５３１９５４１９５５１９５６１９５７１９５８１９５９
年径流深/ｍｍ
５３８．３６２４．９６６３．２５９１．７５５７．２９９８．０６４１．５３４１．１
某站实测年径流资料
年份１９６０１９６１１９６２１９６３１９６４１９６５１９６６１９６７
年径流深/ｍｍ
９６４．２６８７．３５４６．７５０９．９７６９．２６１５．５４１７．１７８９．３
年份１９６８１９６９１９７０１９７１１９７２１９７３１９７４１９７５
年径流深/ｍｍ
７３２．９１０６４．５６０６．７５８６．７５６７．４５８７．７７０９．０８８３．５
Ｍｅａｓｕｒｅｄａｎｎｕａｌｒｕｎｏｆｆｄａｔａｏｆａｓｔａｔｉｏｎ
首先，将表１中的年径流资料记为变量Ａ，调入M ａｔｌａｂ程序，在ＣｏｍｍａｎｄＷｉｎｄｏｗ运行“ｃａｎｓｈｕ”得变差系数C v ＝０．２６３３１”，取C v ＝０．３，并分别取C s 为２C v 、２．５C v 、３C v ，查“皮尔逊Ⅲ型频率曲线的模比系数K p 值表”得到频率Ｐ（％）０．０１、０．１、０．２、０．３３、０．５、１、２、５、１０、２０、５０、７５、９０、９５、９９时相对应的K p 值，并在M ａｔｌａｂ
中分别
-１２５４-记为变量K p 1、K p ２、K p ３.
河南科学
第31卷第8期
然后，在M ａｔｌａｂ的ＣｏｍｍａｎｄＷｉｎｄｏｗ运行“ｚｕｏｂｉａｏ．ｍ”、“ｈｕｉｔｕ．ｍ”子程序，得到三次配线结果（可根据原散点图与配线情况的配合情况适当修改程序确定配线的次数，直到找到理想曲线）如图２所示.
观察图２中三条适配线与原数据散点图的拟合情况，选择最匹配的适线为最终结果. 若在实际工作中对拟合结果不满意，可重新进行适配，直到结果满意为止.
４结语
通过扩充M ａｔｌａｂ的库函数，编写３个Ｍ文件，绘制了Ｐ－Ⅲ水文频率曲线，通过目估适线法选取一条较理想的理论频率曲线. 通过实例分析，验证了该程序的正确性. 本程序适于水文学中对某些水文特征值的分析计算，具有可重复性操作，方便水文学的各种计算，为水文计算增加了一种新的途径. 参考文献：
［１］唐启义，冯明光．ＤＰＳ数据处理系统：实验设计、统计分析及数据挖掘［Ｍ］．北京：科学出版社，２００７. ［２］耿鸿江．Ｅｘｃｅ在Ｐ－Ⅲ型分布频率计算中的应用研究［Ｊ］．水电能源科学，２００２，２０（３）：４１－４３．
［３］王双银，向友珍，朱晓群，等．基于ＥＸＣＥＬ的水文频率计算软件开发［Ｊ］．西北农林科技大学学报：自然科学版，２００６，３４（４）：
１１３－１１６．
［４］赵培颖，金
冶，张忠孝．ＶｉｓｕａｌＢａｓｉｃ在绘制Ｐ－Ⅲ频率曲线中的应用［Ｊ］．水利规划与设计，２００８（２）：５５－５７．
［５］宋孝玉，马细霞．工程水文学［Ｍ］．郑州：黄河水利出版社，２００９.
［６］高会生，刘童娜，李聪聪．ＭＡＴＬＡＢ实用教程［Ｍ］．2版. 北京：电子工业出版社，２０１１.
［７］黄飞仁，黄汉球，曾英先．ＭＡＴＬＡＢ在水文水能计算中的应用［Ｊ］．红水河，２００６，２５（１）：１０４－１０７. ［８］程银才，范世香，李明华．一种新的水文频率计算方法［Ｊ］．水文，２００８，２８（１）：５９－６０.
［９］许义和，魏晓妹．基于Ｍａｔｌａｂ的Ｐ－Ⅲ型曲线绘制软件的研发与应用［Ｊ］．水电能源科学，２０１０，２８（７）：１５－１７.
（编辑张松林）
水文频率曲线绘制的Matlab程序设计
作者：作者单位：刊名：英文刊名：年，卷(期)：
李清富， 闫鹏飞， 孙静涛， 程婷， Li Qingfu， Yan Pengfei， Sun Jingtao， Cheng Ting
郑州大学水利与环境学院,郑州,450001河南科学
Henan Science2013(8)
1. 唐启义;冯明光 DPS数据处理系统:实验设计、统计分析及数据挖掘 2007
2. 耿鸿江 Exce在P-Ⅲ型分布频率计算中的应用研究[期刊论文]-水电能源科学 2002(03)
3. 王双银;向友珍;朱晓群 基于EXCEL的水文频率计算软件开发[期刊论文]-西北农林科技大学学报(自然科学版)2006(04)
4. 赵培颖;金冶;张忠孝 Visual Basic在绘制P-Ⅲ频率曲线中的应用[期刊论文]-水利规划与设计 2008(02)5. 宋孝玉;马细霞 工程水文学 2009
6. 高会生;刘童娜;李聪聪 MATLAB实用教程 2011
7. 黄飞仁;黄汉球;曾英先 MATLAB在水文水能计算中的应用[期刊论文]-红水河 2006(01)8. 程银才;范世香;李明华 一种新的水文频率计算方法[期刊论文]-水文 2008(01)
9. 许义和;魏晓妹 基于Matlab的P-Ⅲ型曲线绘制软件的研发与应用[期刊论文]-水电能源科学 2010(07)
引用本文格式：李清富. 闫鹏飞. 孙静涛. 程婷. Li Qingfu. Yan Pengfei. Sun Jingtao. Cheng Ting 水文频率曲线绘制的Matlab程序设计[期刊论文]-河南科学 2013(8)
第31卷第8期2013年8月
文章编号：1004-3918（2013）08-1250-05
河南科学
HENAN SCIENCE
Vol.31No.8Aug. 2013
水文频率曲线绘制的M ａｔｌａｂ程序设计
李清富，
摘
闫鹏飞，孙静涛，程婷
（郑州大学水利与环境学院，郑州４５０００１）
要：基于目估适线法，选用皮尔逊Ⅲ型频率曲线作为水文频率分布线型，在M ａｔｌａｂ平台上编写三个Ｍ文件，分
軃、C v 的求解以及水文频率曲线的绘制. 通过多次的适配，找到与实测水别实现海森频率格纸的绘制、频率曲线参数x
文资料配合较好的理想曲线的参数作为设计值并估算指定的水文特征值. 并以某一水文测站２４年的年径流实测资料为基础，检验了方法的正确性. 为水文学中相关计算以及水文频率曲线的绘制提供了一种新的实现方法. 关键词：M ａｔｌａｂ；目估适线法；海森频率格纸；皮尔逊Ⅲ型频率曲线中图分类号：Ｐ３３３．９；ＴＰ３１２
文献标识码：A
ＤｅｖｅｌｏｐｍｅｎｔｏｆM ａｔｌａｂＰｒｏｇｒａｍｆｏｒＨｙｄｒｏｌｏｇｉｃａｌＦｒｅｑｕｅｎｃｙＣｕｒｖｅＤｒａｗｉｎｇ
Li Qingfu ，Ｙan Ｐｅｎｇｆｅｉ，
Ｓun Ｊｉｎｇｔａｏ，
Ｃheng Ｔｉｎｇ
（ＳｃｈｏｏｌｏｆＥｎｖｉｒｏｎｍｅｎｔａｎｄＷａｔｅｒＣｏｎｓｅｒｖａｎｃｙＥｎｇｉｎｅｅｒｉｎｇ，ＺｈｅｎｇｚｈｏｕＵｎｉｖｅｒｓｉｔｙ，Ｚｈｅｎｇｚｈｏｕ４５０００１，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：The Ｐ－ＩＩＩｆｒｅｑｕｅｎｃｙｃｕｒｖｅｉｓｓｅｌｅｃｔｅｄａｓｔｈｅｈｙｄｒｏｌｏｇｉｃａｌｆｒｅｑｕｅｎｃｙｄｉｓｔｒｉｂｕｔｉｏｎｏｆｌｉｎｅａｒ，ｂａｓｅｄｏｎｖｉｓｕａｌ
ａｄａｐｔａｂｉｌｉｔｙｌｉｎｅ．ＴｈｒｅｅＭ－ｆｉｌｅｓｗａｓｗｒｉｔｔｅｎｉｎＭａｔｌａｂｐｌａｔｆｏｒｍｔｏｄｒａｗｔｈｅＨａｚｅｎｆｒｅｑｕｅｎｃｙｇｒｉｄｐａｐｅｒ，ｔｏｓｏｌｖｅ軃、C v ，ｔｏｄｒａｗｔｈｅｈｙｄｒｏｌｏｇｉｃａｌｆｒｅｑｕｅｎｃｙｃｕｒｖｅ，ｒｅｓｐｅｃｔｉｖｅｌｙ．Ｗｉｔｈｔｈｅｔｈｅｆｒｅｑｕｅｎｃｙｃｕｒｖｅｐａｒａｍｅｔｅｒｓｓｕｃｈａｓx
ｍｕｌｔｉｐｌｅａｄａｐｔａｔｉｏｎ，ｔｈｅｉｄｅａｌｃｕｒｖｅｔｈａｔｆｉｔｓｗｅｌｌｗｉｔｈｔｈｅｍｅａｓｕｒｅｄｈｙｄｒｏｌｏｇｉｃａｌｄａｔａｉｓｆｏｕｎｄ．Ｓｏｗｅｃａｎｅｓｔｉｍａｔｅ
ｔｈｅｖａｌｕｅｏｆｔｈｅｓｐｅｃｉｆｉｅｄｈｙｄｒｏｌｏｇｉｃａｌｃｈａｒａｃｔｅｒｉｓｔｉｃｓｂｙｔｈｅｕｓｉｎｇｏｆｔｈｅｉｄｅａｌｃｕｒｖｅ．Ａｎｎｕａｌｒｕｎｏｆｆｍｅａｓｕｒｅｄｄａｔａｏｆａｈｙｄｒｏｌｏｇｉｃａｌｓｔａｔｉｏｎ，２４ｙｅａｒｓｉｓｕｓｅｄｔｏｔｅｓｔｔｈｅｃｏｒｒｅｃｔｎｅｓｓｏｆｔｈｅｍｅｔｈｏｄ．Ｉｔｐｒｏｖｉｄｅｓａｎｅｗｉｍｐｌｅｍｅｎｔａｔｉｏｎｆｏｒｈｙｄｒｏｌｏｇｉｃａｌｃａｌｃｕｌａｔｉｏｎｓａｎｄｔｈｅｄｒａｗｉｎｇｏｆｔｈｅｈｙｄｒｏｌｏｇｉｃａｌｈｙｄｒｏｌｏｇｙｆｒｅｑｕｅｎｃｙｃｕｒｖｅ．
Ｋｅｙｗｏｒｄｓ：M ａｔｌａｂ；ｔａｒｇｅｔａｄａｐｔａｂｉｌｉｔｙｌｉｎｅ；ｔｈｅＨａｚｅｎｆｒｅｑｕｅｎｃｙｇｒｉｄｐａｐｅｒ；Ｐ－ＩＩＩｆｒｅｑｕｅｎｃｙｃｕｒｖｅ我国水文频率计算一直采用皮尔逊型Ⅲ型频率曲线作为水文分布频率曲线. 随着计算机技术的发展，水文频率分析软件应运而生. 如由浙江大学唐启义教授研制的一款综合性的大型统计软件ＤＰＳ（ＤａｔａＰｒｏｃｅｓｓｉｎｇＳｙｓｔｅｍ）统计软件，通用多功能数理统计和数学模型处理软件系统，将数值计算、统计分析、模型模拟以及画线制表等功能融为一体[1]；耿鸿江[2]，王双银[3]等人在研究了基于Ｅｘｃｅｌ的Ｐ－Ⅲ型分布曲线的绘制方法；赵培颖[4]等人研究了ＶｉｓｕａｌＢａｓｉｃ在绘制Ｐ－Ⅲ频率曲线中的应用. 本文在已有学者的研究基础上，对如何应用Ｍａｔｌａｂ绘制水文频率曲线进行了研究，并给出了利用Ｍａｔｌａｂ绘制水文频率曲线的Ｍ文件，具有重复可操作性.
１计算理论
目估适线法
目估适线法[5]是我国估计水文频率曲线统计参数最主要的方法. 它是以经验频率点据为基础，选配一条拟合较好的理论频率曲线，利用频率曲线与经验频率点据相配合来估计水文要素总体的统计规律. 具体步骤如下：
收稿日期：2013-01-10
基金项目：国家自然科学基金项目（５０９７８２３４）
作者简介：李清富（１９６６－），男，河南林州人，教授，博士，从事交通运输工程领域的研究
闫鹏飞（１９９１－），男，河南商丘人，主要从事道路桥梁方面的研究．
１．１
2013年8月李清富，等：水文频率曲线绘制的M ａｔｌａｂ程序设计
-１２５１-
①绘制频率格纸；②将实测资料由大到小排列，计算各项的经验频率，在频率格纸上点绘经验点据. ③选軃、C v 、C s 用皮尔逊Ⅲ型频率曲线作为水文频率分布线型，采用矩法或其他方法估计出频率曲线参数的初估值x 軃、C v 、C s ，查（凭经验初选为C v 的某一倍数）；④根据拟定的x “皮尔逊Ⅲ型频率曲线的模比系数K p 值表”，得到模比系数K p 值，并计算x p 值；⑤以x p 为纵坐标，P 为横坐标，即可得到频率曲线，将此线画在绘有经验点据的軃、C v 、C s 值以及各x 軃、C v 、C s 所对应的曲线，直到找到图上；⑥分析曲线与经验点据的配合情况，通过不断调整x 理想的频率曲线为止；⑦求指定的水文变量变化值.
１．２Ｐ－Ⅲ型曲线密度函数
英国生物学家皮尔逊提出了１３种非正态的分布曲线，其中第Ⅲ型曲线被引入水文计算中，用途比较广泛. Ｐ－Ⅲ型曲线[5，8鄄9]是一条一端有限、一端无限的不对称单峰、正偏曲线，数学上成为伽马分布，其概率密度函数为：
βa -1-β（x-a ）
x-a 0） e ，（ｆx ）=Γ（α）
α
式中：Γ（α）—α的伽马函数；α、β、a 0—ＰⅢ型曲线分布的形状、尺度和位置参数ａ＞０，β＞０.
軃、C v 、C s 具有显然α、β、a 0这３个参数求出来以后该密度函数便可以确定. 3个参数与总体的3个参数x 如下关系：
α=42；
C s
β=2；軃C v C s x
軃a 0=x （1-2C v
. C s
在水文计算中，一般需要求出指定频率Ｐ所相应的随机变量X 的取值x p ，也就是通过对密度曲线进行积分[5]，即
β
P =P （x ≥x p ）=
Γ（α）
α
乙
x p
∞
（x-a 0） e
a -1-β（x-a 0）
d x .
２程序设计
程序设计流程图如图１所示所示，主要包括５个步骤：①求
解总体参数；②查“皮尔逊Ⅲ型频率曲线的模比系数K p 值表”确定K p 值；③绘制海森频率坐标格纸；④绘制Ｐ－Ⅲ型频率曲线；⑤选线.
根据流程图，利用M ａｔｌａｂ[6鄄7]平台编制“ｚｕｏｂｉａｏ．ｍ”、“ｚｕｏｂｉａｏ．ｍ”、
“ｈｕｉｔｕ．ｍ”共3个Ｍ文件，分别用来求解总体参数、绘制海森频率格纸以及绘制Ｐ－Ⅲ型频率曲线，并都保存在Ｍａｔｌａｂｗｏｒｋ子目录中，以便调用. ２．１海森频率格纸的绘制
海森频率格纸是水文计算中常用的一种坐标格纸，其横坐标的分化就是按把标准正态频率曲线拉成一条直线的原理计算出
图１水文频率曲线计算流程图来的. 这种频率格纸的纵坐标仍是普通分格，但横坐标的分格是
Ｆｉｇ．１Ｔｈｅｃａｌｃｕｌａｔｉｏｎｆｌｏｗｃｈａｒｔｏｆｈｙｄｒｏｌｏｇｉｃａｌ
不相等的中间分格较密，越往两端越稀，中间分格在Ｐ＝５０％的两
ｆｒｅｑｕｅｎｃｙｃｕｒｖｅ
[5]
端是对称的.
利用M ａｔｌａｂ平台编制绘制海森频率格纸的Ｍ文件，ｚｕｏｂｉａｏ函数的编程和变量解释如下：ｆｕｎｃｔｉｏｎｚｕｏｂｉａｏ％定义坐标函数
ｘ＝［０．０１：０．０１：．１，０．２：．１：１，１．２：．２：２，３：２０，２２：２：８０，８１：９８，９８．２：．２：９９，９９．１：．１：９９．９］；％横坐标的值ｙ＝ｎｏｒｍｉｎｖ（ｘ/１００，０，１）；ｙ＝ｙ－ｙ（１）；
-１２５２-ｆｏｒｉ＝１：ｓｉｚｅ（ｙ，２）
ｌｉｎｅ（［ｙ（ｉ），ｙ（ｉ）］，［０２０００］，′ｃｏｌｏｒ′，′ｋ′）；ｅｎｄ
ｚ＝０：１００：２０００；
河南科学
第31卷第8期
ｆｏｒｉ＝１：２１；％纵坐标的间距设置，可以自己设置ｌｉｎｅ（［０，ｙ（ｅｎｄ）］，［ｚ（ｉ）ｚ（ｉ）］，′ｃｏｌｏｒ′，′ｋ′）；
ｅｎｄ
ｈ＝ｆｉｎｄｏｂｊ（′ｔｙｐｅ′，′ａｘｅｓ′）；ｓｅｔ（ｈ，′ｘｔｉｃｋ′，［］，′ｘｌｉｍ′，［０ｙ（ｅｎｄ）］，′ｙｌｉｍ′，［０２０００］）；％同时去掉ｘ轴和ｙ轴的刻度
ｘｘ＝［０．０１０．０５０．１０．５１．０５１０１５２０３０４０５０６０７０８０８５９０９５９９９９．５９９．９］；％横坐标的设置，可以自己设置ｙｙ＝ｎｏｒｍｉｎｖ（ｘｘ/１００，０，１）；ｙｙ＝ｙｙ－ｙｙ（１）；ｓｅｔ（ｈ，′Ｘｔｉｃｋ′，ｙｙ）；ｓｅｔ（ｈ，′ＸｔｉｃｋｌａｂｅＬ′，｛′０．０１′，′０．０５′，′０．１′，′０．５′，′１．０′，′５′，′１０′，′１５′，′２０′，′３０′，′４０′，′５０′，′６０′，′７０′，′８０′，′８５′，′９０′，′９５′，．．．
′９９′，′９９．５′，′９９．９′｝）；
ｓｅｔ（ｈ，′Ｙｔｉｃｋ′，ｚ）；ｘｌａｂｅｌ（′Ｐ（％）′）；％对ｘ轴进行标注ｙｌａｂｅｌ（′Ｑｍ（ｍ/ｓ）′）；％对ｙ轴进行标注，可根据水文特征值设置２．２求解总体参数的Ｍ文件编写
利用M ａｔｌａｂ平台编制绘制Ｐ－Ⅲ型频率曲线的Ｍ文件，ｃａｎｓｈｕ函数的编程和变量解释如下：％％参数计算ｆｏｒｍａｔｓｈｏｒｔａ＝ｉｎｐｕｔ（′水文学数据′）；Ａ＝ｓｏｒｔ（ａ，′ｄｅｓｃｅｎｄ′）；％对水文学数据进行排序ｎ＝ｌｅｎｇｔｈ（Ａ）；％计算数据的个数ｑ＝１：ｎ；ｑ＝ｑ．′；Ｐ＝ｑ．（/ｎ＋１）；％计算经验频率ｐ＝ｎｏｒｍｉｎｖ（Ｐ，０，１）；％数据Ｐ所对应的正态分布的Ｘ值Ｐ１＝［０．０１；０．１；０．２；０．３３；０．５；１；２；５；１０；２０；５０；７５；９０；９５；９９；］；ｐ１＝ｎｏｒｍｉｎｖ（Ｐ１．/１００，０，１）；％数据ｐ１所对应的正态分布的Ｘ值Ｍ＝ｎｏｒｍｉｎｖ（０．０００１，０，１）；ｐ２＝－Ｍ＋ｐ１；％转换坐标值（理论频率）ｐ３＝－Ｍ＋ｐ；％转换坐标值（经验频率）Ｑ＝ｍｅａｎ（Ａ）；％求平均值Ｋ＝Ａ．/Ｑ；％求模比系数Ｋ１＝Ｋ－１；Ｋ２＝Ｋ１．＾２；Ｃｖ＝ｓｑｒｔ（ｓｕｍ（Ｋ２））/ｓｑｒｔ（２４－１）；％计算变差系数ｄｉｓｐ（［′平均值Ｑ＝′，ｎｕｍ２ｓｔｒ（Ｑ）］）ｄｉｓｐ（［′变差系数Ｃｖ＝′，ｎｕｍ２ｓｔｒ（Ｃｖ）］）２．３
绘制Ｐ－Ⅲ型频率曲线的Ｍ文件编写
利用M ａｔｌａｂ平台编制绘制Ｐ－Ⅲ型频率曲线的Ｍ文件，ｈｕｉｔｕ函数的编程和变量解释如下：％％绘制图像ｐｌｏｔ（ｐ３，Ａ，′＊ｒ′）；％绘制经验频率的散点
2013年8月李清富，等：水文频率曲线绘制的M ａｔｌａｂ程序设计
-１２５３-
ｈｏｌｄｏｎ；ｇｒｉｄｏｎ；
Ｋｐ＝ｉｎｐｕｔ（′输入查表参数′）；Ｘｐ＝Ｑ＊Ｋｐ；％Ｘｐ的理论计算值ｐｌｏｔ（ｐ２，Ｘｐ，′ｋ′）；％第一次配线ｈｏｌｄｏｎ；ｇｒｉｄｏｎ；
Ｋｐ＝ｉｎｐｕｔ（′输入查表参数′）；Ｘｐ＝Ｑ＊Ｋｐ；％Ｘｐ的理论计算ｐｌｏｔ（ｐ２，Ｘｐ，′ｂ′）；％第二次配线ｈｏｌｄｏｎ；ｇｒｉｄｏｎ；
Ｋｐ＝ｉｎｐｕｔ（′输入查表参数′）；
Ｘｐ＝Ｑ＊Ｋｐ；Ｘｐ的理论计算值ｐｌｏｔ（ｐ２，Ｘｐ，′ｇ′）；％第三次配线ｌｅｇｅｎｄ（′原数据散点′，′第一次配线′，′第二次配线′，′第三次配线′）％图例标注
R /m m
P /%
图２
Ｆｉｇ．２
频率计算成果图
Ｔｈｅｆｉｇｕｒｅｏｆｔｈｅｆｒｅｑｕｅｎｃｙｃａｌｃｕｌａｔｉｏｎｒｅｓｕｌｔｓ
３实例
某测站的年径流资料[5]，如表１所示.
表１
Ｔａｂ. １
年份１９５２１９５３１９５４１９５５１９５６１９５７１９５８１９５９
年径流深/ｍｍ
５３８．３６２４．９６６３．２５９１．７５５７．２９９８．０６４１．５３４１．１
某站实测年径流资料
年份１９６０１９６１１９６２１９６３１９６４１９６５１９６６１９６７
年径流深/ｍｍ
９６４．２６８７．３５４６．７５０９．９７６９．２６１５．５４１７．１７８９．３
年份１９６８１９６９１９７０１９７１１９７２１９７３１９７４１９７５
年径流深/ｍｍ
７３２．９１０６４．５６０６．７５８６．７５６７．４５８７．７７０９．０８８３．５
Ｍｅａｓｕｒｅｄａｎｎｕａｌｒｕｎｏｆｆｄａｔａｏｆａｓｔａｔｉｏｎ
首先，将表１中的年径流资料记为变量Ａ，调入M ａｔｌａｂ程序，在ＣｏｍｍａｎｄＷｉｎｄｏｗ运行“ｃａｎｓｈｕ”得变差系数C v ＝０．２６３３１”，取C v ＝０．３，并分别取C s 为２C v 、２．５C v 、３C v ，查“皮尔逊Ⅲ型频率曲线的模比系数K p 值表”得到频率Ｐ（％）０．０１、０．１、０．２、０．３３、０．５、１、２、５、１０、２０、５０、７５、９０、９５、９９时相对应的K p 值，并在M ａｔｌａｂ
中分别
-１２５４-记为变量K p 1、K p ２、K p ３.
河南科学
第31卷第8期
然后，在M ａｔｌａｂ的ＣｏｍｍａｎｄＷｉｎｄｏｗ运行“ｚｕｏｂｉａｏ．ｍ”、“ｈｕｉｔｕ．ｍ”子程序，得到三次配线结果（可根据原散点图与配线情况的配合情况适当修改程序确定配线的次数，直到找到理想曲线）如图２所示.
观察图２中三条适配线与原数据散点图的拟合情况，选择最匹配的适线为最终结果. 若在实际工作中对拟合结果不满意，可重新进行适配，直到结果满意为止.
４结语
通过扩充M ａｔｌａｂ的库函数，编写３个Ｍ文件，绘制了Ｐ－Ⅲ水文频率曲线，通过目估适线法选取一条较理想的理论频率曲线. 通过实例分析，验证了该程序的正确性. 本程序适于水文学中对某些水文特征值的分析计算，具有可重复性操作，方便水文学的各种计算，为水文计算增加了一种新的途径. 参考文献：
［１］唐启义，冯明光．ＤＰＳ数据处理系统：实验设计、统计分析及数据挖掘［Ｍ］．北京：科学出版社，２００７. ［２］耿鸿江．Ｅｘｃｅ在Ｐ－Ⅲ型分布频率计算中的应用研究［Ｊ］．水电能源科学，２００２，２０（３）：４１－４３．
［３］王双银，向友珍，朱晓群，等．基于ＥＸＣＥＬ的水文频率计算软件开发［Ｊ］．西北农林科技大学学报：自然科学版，２００６，３４（４）：
１１３－１１６．
［４］赵培颖，金
冶，张忠孝．ＶｉｓｕａｌＢａｓｉｃ在绘制Ｐ－Ⅲ频率曲线中的应用［Ｊ］．水利规划与设计，２００８（２）：５５－５７．
［５］宋孝玉，马细霞．工程水文学［Ｍ］．郑州：黄河水利出版社，２００９.
［６］高会生，刘童娜，李聪聪．ＭＡＴＬＡＢ实用教程［Ｍ］．2版. 北京：电子工业出版社，２０１１.
［７］黄飞仁，黄汉球，曾英先．ＭＡＴＬＡＢ在水文水能计算中的应用［Ｊ］．红水河，２００６，２５（１）：１０４－１０７. ［８］程银才，范世香，李明华．一种新的水文频率计算方法［Ｊ］．水文，２００８，２８（１）：５９－６０.
［９］许义和，魏晓妹．基于Ｍａｔｌａｂ的Ｐ－Ⅲ型曲线绘制软件的研发与应用［Ｊ］．水电能源科学，２０１０，２８（７）：１５－１７.
（编辑张松林）
水文频率曲线绘制的Matlab程序设计
作者：作者单位：刊名：英文刊名：年，卷(期)：
李清富， 闫鹏飞， 孙静涛， 程婷， Li Qingfu， Yan Pengfei， Sun Jingtao， Cheng Ting
郑州大学水利与环境学院,郑州,450001河南科学
Henan Science2013(8)
1. 唐启义;冯明光 DPS数据处理系统:实验设计、统计分析及数据挖掘 2007
2. 耿鸿江 Exce在P-Ⅲ型分布频率计算中的应用研究[期刊论文]-水电能源科学 2002(03)
3. 王双银;向友珍;朱晓群 基于EXCEL的水文频率计算软件开发[期刊论文]-西北农林科技大学学报(自然科学版)2006(04)
4. 赵培颖;金冶;张忠孝 Visual Basic在绘制P-Ⅲ频率曲线中的应用[期刊论文]-水利规划与设计 2008(02)5. 宋孝玉;马细霞 工程水文学 2009
6. 高会生;刘童娜;李聪聪 MATLAB实用教程 2011
7. 黄飞仁;黄汉球;曾英先 MATLAB在水文水能计算中的应用[期刊论文]-红水河 2006(01)8. 程银才;范世香;李明华 一种新的水文频率计算方法[期刊论文]-水文 2008(01)
9. 许义和;魏晓妹 基于Matlab的P-Ⅲ型曲线绘制软件的研发与应用[期刊论文]-水电能源科学 2010(07)
引用本文格式：李清富. 闫鹏飞. 孙静涛. 程婷. Li Qingfu. Yan Pengfei. Sun Jingtao. Cheng Ting 水文频率曲线绘制的Matlab程序设计[期刊论文]-河南科学 2013(8)