第３３卷第１期 水文 Ｖ０１．３３ Ｎｏ．１ ２０１３年２月 ＪＯＵＲＮＡＬ ＯＦ ＣＨＩＮＡ ＨＹＤＲ０ＬＯＧＹ Ｆｅｂ．．２０１３ 应用ｍａｔｌａｂ求解小流域推理公式的方法 田景环，梁文涛 （华北水利水电学院，河南郑州４５００１１） 摘 要：针对《水利水电工程设计洪水计算规范》在推求小流域洪水时采用的推理公式，提出以图解法及 牛顿迭代法为基础，用ｍａｔｌａｂ软件在各种汇流类型下求解推理公式的洪峰流量及汇流时间，该方法计 算速度快、精度高。而且简便、实用、易操作。 关键词：推理公式；图解法；牛顿迭代法；ｍａｄａｂ 中图分类号：ＴＶ１２２ 文献标识码：Ａ 文章编号：１０００—０８５２（２０１３）０１—００７９—０３ ｌ 概述 ｈ Ｈｔ—ｍ （３） Ｈ ５。￡ （４） 小流域洪水计算在水利水电工程规划设计中占有 重要地位。１９５６年，原水利科学研究院水文研究所的 ｈＦ ＝ （５） 林平一提出了以推理公式为基础的计算最大流量方 上 ／Ｌ ｌ 法。１９５８年陈家琦等人提出了水利科学院推理公式， 凡 『（１＿凡） ］ — 、 （６） 并很快在实际运用中得到了普及ｌ１】。现行《水利水电工 程设计洪水计算规范》中推荐小流域设计洪水采用水 一 （７） 科院推理公式计算进行计算。求解推理公式洪峰流量 Ｑ 及汇流时间 常采用图解法求解法、迭代法等，但 公式（１）～（７）中： 为河道长度；．，为河道比降；Ｆ为流 图解法比较繁琐，而且手绘图形精确度不高：迭代法的 域面积；Ｑ 为设计洪峰流量；ｍ为汇流参数；／ｘ为平均 收敛速度较慢，虽然引进了牛顿迭代法翻解决了这一 降雨入渗率；ｎ、．ｓ。为暴雨参数；０．２７８为单位换算系数； 问题，但是迭代法求解推理公式所适用的汇流类型还 ｔ 为产流历时；．ｒ为汇流历时； 为单一洪峰的净雨；日 有待进一步研究。本文以图解法和牛顿迭代法为基础， 为ｔ时段最大雨量。 提出利用ｍａｔｌａｂ软件求解小流域推理公式的方法。该 ２．２图解法原理 方法可求解各种同汇流类型下推理公式的洪峰流量 先假设出一组Ｑ 值，根据公式（１）计算出一组相 Ｑ 及汇流时间 ，不仅简便、实用、易操作。且计算速度 应的 值，绘出ｐｍ￣ 关系曲线。假设汇流类型，在ｔ＜ｔ。 快、精度高。 或ｔ ｔ 时，根据公式（２）（３）（４）（５）确定不同汇流类 型时的ｐ ￡关系式。取之前Ｑ ～ｆ计算中求出的 值， ２推理公式及求解原理 令 ＝ ，即可求出一组对应的Ｑ ～ｔ值，将求出的两组 ２．１ 推理公式 ｐｍ￣了、Ｑ ～ｔ值点绘在同一坐标系内，描出两条关系曲 线，两条曲线的交点即为所求的设计洪峰流量值Ｑ 。 仁０Ｉ２７８ （１） 若有交点，求出Ｑ 对应的时间ｔ，并通过公式（７）求出 Ｑｍ＝Ｏ．２７８ｈ ，Ｆ （２） 的ｔ 与之对比，是否符合之前对流域汇流类型的假设， 不符合假设则需重新设定流域汇流类型进行计算。若 ｔ＜ｔ，时． 无交点，则重新假设Ｑ 值，继续上述计算＿３ｌ。 收稿日期：２０１２—０５—１４ 作者简介：田景环（１９６５－），女，河北涿州人，教授，主要从事水力学及河流动力学的科研和教学工作。Ｅ－ｍａｉｌ：ｔｉａｎｊｉｎｇｈｕａｎ＠ｎｃｗｕ．ｅｄｕ．ｅｎ 水文 第３３卷 ２．３迭代法公式 ｆｕｎｃｔｉｏｎ［Ｑｍ，Ａ，Ｂ，Ｑｍ０，ｔ］进行计算。 －Ｑｍ＝Ｏ．２７８（ａ－／ｘ）Ｆ＝０．２７８（巫－ ｒ ｎ／ｚ）Ｆ （）８ 当ｔ＜ｔ ，即全面汇流时，将公式（１）（３）（４）代人公式（２）， 即可得到公式（９），也就是迭代法求解推理公式的方程。 建立Ｍ文件如下，并保存： 式中：Ｏ／为平均暴雨强度。 ２．４迭代法原理 ｆｕｎｃｔｉｏｎ［Ｑｍ，Ａ，Ｂ，Ｑｍ０，ｔ】＝ｔ２（Ｓ，ｍ，Ｊ，Ｌ，Ｕ，ｎ，Ｆ）％ 简单迭代法 原理就是将方程厂（ ）＝Ｏ整理成形如 ＝ （ ）的等价方程，然后给定一个初值‰，将 。代入 公式右端，可得出 。＝ （ 。），再将 代入方程右端得到 ，依次类推，可得到Ｘｋ＋ｌ＝ （钆），ｋ＝０，１，２…如果迭代 方程 ＝ （ ）是收敛的，则必有ｘ＊－￣ｏ（ｘ ）， 即为方程的 解。牛顿迭代法是设Ｘ 为方程．厂（ ）＝０的一个近似根， 把ｆ（ｘ）在瓢处泰勒展开，得到近似方程－厂（钆）－ｆ （钆） ｘ－ｘ ）＝０，从而构造出形如 。 一 的迭代方程， ．，Ｌ 大大提高了迭代方程的收敛速度，减少计算步骤。 ２．５对推理公式的整理 把公式（１）代入（８）中并整理，得到 Ｑｍ＝０．２７８－ （ ）ｎＱ 一０．２７８／￣Ｆ （９） 设Ａ：０．２７８ （ ）ｎ，Ｂ：０．２７８ｐ，Ｆ，则 Ｑｍ＝ＡＱ Ｂ （１０） 可得到迭代公式 Ｑｍ＝Ｑ 。．． （１１） １一 Ｑ 。 设公式（１０）中的Ｂ值为０，则此时Ｑ 取得最大值 Ｑ ｏ＝Ａ ‘ ，将此值作为迭代的初值ｌ ５ｌ。 ３推理公式的ｍａｔｌａｂ求解 求解推理公式即为求解公式（１）与（２）的联立方程。 当ｔ ｔｃ，汇流类型为部分汇流时，ｔ＝Ｔ，将（１）（５） （６）代入（２）中，整理得 Ｑｍ＝Ｉ 】 ：Ｊ］｝ ４／３ 建立Ｍ文件如下，并保存： ｆｕｎｃｔｉｏｎ【Ｑｍ］＝ｆｌ（Ｓ，ｍ，Ｊ，Ｌ，Ｕ，ｎ，Ｆ）％各参数 Ｑｍ＝（（Ｆ ｍ Ｊ （１／３） ｎ Ｓ （（１一ｎ） Ｓ／ｕ） （（１一ｎ）／ｎ））／Ｌ）＾（４／３） ｔ＝０．２７８ Ｌ ｍ Ｊ （１／３） Ｑｍ＾（１／４）） （一１） ｔｃ＝（（１一ｎ） Ｓ／ｕ） （１／ｎ） 在指令窗Ｅｌ输入ｆ１（Ｓ，ｉｎ，Ｊ，Ｌ，Ｕ，ｎ，Ｆ），即可得到 结果ｐ～ｔ、ｔ 。若ｔ＞－ｔ ，则可取，若ｔ＜ｔ ，则运行如下 各参数 Ａ＝（０．２７８ （１一ｎ）） Ｓ （（ｍ Ｊ＾（１／３）） ｎ） （一ｎ） Ｆ Ｂ＝０．２７８＊ｕ＊Ｆ Ｑｍ０＝Ａ （４／（４一ｎ））％迭代初值 ｆｏｒ ｉ—ｌ：１０％迭代次数 Ｑｍ０＝Ｑｍ０一（Ｑｍ０－Ａ Ｑｍ（）＾（ｎ／４）＋Ｂ）／（１一Ａ ｎ／４ Ｑｍ０ （（ｎ一４）／４）） ｅｎｄ Ｑｍ＝Ｑｍ０ ｔ＝０．２７８ Ｌ ｍ Ｊ＾（１／３） Ｑｍ＾（１／４）） （一１） 在ｍａｔｌａｂ的指令窗口中输入ｆ２（Ｓ，ｍ，Ｊ，Ｌ，ｕ，ｎ， Ｆ），即可得到结果Ｑ～ｔ。 计算时先假设ｔ＞－－ｔ ，调用ｆｕｎｃｔｉｏｎ【Ｑｍ】，求出Ｑ～ ｔ、ｔ ，若求得ｔ ，则Ｑｍ为所求设计洪水、ｔ为汇流时 间。若求得ｔ＜ｔ ，调用ｆｕｎｃｔｉｏｎ［Ｑｍ，Ａ，Ｂ，Ｑｍ０，ｔ】，即可 得到结果Ｑ～ｔ。 ４实例 ４．１ 算例１ 已知某流域特征参数如表１，求该流域百年一遇 设计洪峰流量。 表１流域特征参数 Ｔａｂｌｅ １ Ｔｈｅ ｃｈａｒａｃｔｅｒｉｓｔｉｃ ｐａｒａｍｅｔｅｒｓ ｏｆ ｔｈｅ ｗａｔｅｒｓｈｅｄ ４．１．１ 图解法求解 由公式Ｔ＝０．２７８南＝ ，假设 ３０、 ２３５、２４０、２４５ｍ３／ｓ，与之相对应的ｆ分别为２８．９１、 ２８．７５、２８．６０、２８．４６。由此可得Ｑ 一Ｔ关系曲线。假设￡ ￡ ， 则 ： ，于是Ｑ ：０．２７８＿ｈＲ Ｆ＝ ，并取ｔ＝２８．９ｌ、 ２８．７５、２８．６０、２８．４６，可求得与之对应的Ｑ 分别为 ２５１．５、２５２．９、２５４－３、２５５．５ｍ３／ｓ。由此可得Ｑ ～ｆ关系曲 线。在同一坐标纸上，绘制Ｑ ～ 和Ｑ ～ 曲线，两者的 交点就是所求的Ｑ ，即Ｑｍ＝２５９．６ｍ３／ｓ，Ｔ＝２８．Ｏｌｈ。ｔｃ＝２１．２ｈ 第１期 田景环等：应用ｍａｔｌａｂ求解小流域推理公式 ８１ 小于 ，与假设相符。 ４．１．２ ｍａｔｌａｂ求解 大于ｔ。与假设相悖。不可选用。建立并运行ｆｕｎｃｔｉｏｎ ［Ｑｍ，Ａ，Ｂ，Ｑｍ０，ｔ］，在指令窗口中输入： Ｆ２（１７５．４，１，０．００５，５０，２，０．７，１５０） 设ｆ ￡ ，建立并运行ｆｕｎｃｔｉｏｎ［Ｑｍ］，在指令窗口中 输人：ｆ１（６５，０．７，０．００４，４５，２．３，０．７，２３０） 求得Ｑ ＝１０５４．７４９８ｍ％，ｔ＝１４．２６３９ｈ。ｔｃ＝１０６．８４２５ｈ， 大于ｔ。与假设相符，可以选用。 求得Ｑ ＝２５９．１４５０ｍ３／ｓ，ｔ＝２８．０５９９ｈ。ｔｃ＝２１．１９１ｌｈ， 小于ｔ，可以选用。 ４．２算例２ ５结语 针对《水利水电工程设计洪水计算规范》在推求小 流域洪水时采用的推理公式。本文提出了以图解法及 牛顿迭代法为基础．采用ｍａｔｌａｂ软件在各种汇流类型 设某流域特征参数如表２。求该流域百年一遇设 计洪峰流量。 表２流域特征参数 Ｔａｂｌｅ ２ Ｔｈｅ ｃｈａｒａｃｔｅｒｉｓｔｉｃ ｐａｒａｍｅｔｅｒｓ ｏｆ ｔｈｅ ｗａｔｅｒｓｈｅｄ 下求解推理公式的洪峰流量Ｑ 及汇流时间ｆ。实例分 析计算表明，该方法对全面产流和部分产流均适用，不 仅避免了传统图解法计算繁琐、手工绘图不准确的缺 点：克服了迭代法求解部分汇流时计算不适用以及迭代 次数较多、收敛速度慢、计算量大等缺点。而且计算速度 ４．２．１迭代法求解 快、精度高，而且简便、实用、易操作，便于推广使用。 参考文献： 【１］陈家琦，张恭肃．小流域暴雨洪水计算问题【Ｍ］．北京：水利电力出 版社，１９８３．（ＣＨＥＮ Ｊｉａｑｉ，ＺＨＡＮＧ Ｇｏｎｇｓｕ．Ｈｅａｖｙ Ｒａｉｎ ａｎｄ Ｆｌｏｏｄ 由公式Ｑｍｏ＝Ａ ‘ ，选迭代初值为１ １５６．６ｍ３／ｓ，迭 代过程如表３。 表３迭代法计算过程 Ｔａｂｌｅ ３ Ｔｈｅ ｃａｌｃｕｌａｔｉｏｎ ｐｒｏｃｅｓｓ ｏｆ ｔｈｅ ｉｔｅｒａｔｉｏｎ ｍｅｔｈｏｄ Ｃａｌｃｕｌａｔｉｏｎ ｏｆ Ｓｍａｌｌ Ｂａｓｉｎｓ 【Ｍ］．Ｂｅｉｊｉｎｇ：Ｗａｔｅｒ Ｒｅｓｏｕｒｃｅｓ ａｎｄ Ｅｌｅｃｔｉｒｃ Ｐｏｗｅｒ Ｐｒｅｓｓ，１９８３．（ｉｎ Ｃｈｉｎｅｓｅ）） ［２］徐德龙，肖华．小流域设计洪水推理公式计算方法探讨咖．人民长 江，２０００，３１（７）：１３—１４．（ＸＵ Ｄｅｌｏｎｇ，ＸＩＡＯ Ｈｕａ．Ｄｉｓｃｕｓｓｉｏｎ ｏｆ ｄｅ— ｓｉｇｎ ｆｌｏｏｄ ｃａｌｃｕｌａｔｉｏｎ ｏｆ ｓｍａｌｌ ｂａｓｉｎｓ［ｊ］．Ｙａｎｇｔｚｅ Ｒｉｖｅｒ，２０００，３１ （７）：１３－１４．（ｉｎ Ｃｈｉｎｅｓｅ）） 『３１邱林，孙元元，周生通．一种基于ＶＢ求解小流域设计洪峰流量的 图解方法 ［ 水文，２０１２，（１）．（ＱＩＵ Ｌｉｎ，ＳＵＮ Ｙｕａｎｙｕａｎ，ＺＨＯＵ Ｓｈｅｎｇｔｏｎｇ．Ａ ｍｅｔｈｏｄ ｆｏｒ ｄｅｓｉｇｎｅｄ ｐｅａｋ ｆｌｏｏｄ ｃａｌｃｕｌａｔｉｏｎ ｏｆ ｓｍａｌｌ ｂａｓｉｎｓ ｂａｓｅｄ ｏｎ ＶＢ ｇｒａｐｈｉｃａｌ ｍｅｔｈｏｄｓ［Ｊ］．Ｊｏｕｒｎａｌ ｏｆ Ｃｈｉｎａ Ｈｙ— ｄｒｏｌｏｇｙ，２０１２，３２（１）．（ｉｎ Ｃｈｉｎｅｓｅ）） 可求得Ｑ ＝１０５４．７ｍ３／ｓ，ｔ＝１４．３ｈ，ｔｃ＝１０６．８ｈ大于ｔ。 与假设相符，可以选用。 ４．２．２ ｍａｔｌａｂ求解 ［４］李清善，宋士仓．数值方法【Ｍ］．郑州：郑州大学出版社，２００７：１６４— １６６． （ＬＩ Ｑｉｎｇｓｈａｎ，ＳＯＮＧ Ｓｈｉｃａｎｇ． ［Ｍ】．Ｎｕｍｅｒｉｃ Ｍｅｔｈｏｄｓ【Ｍ］． Ｚｈｅｎｇｚｈｏｕ：Ｚｈｅｎｇｚｈｏｕ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ，２００７：１６４－１６６．（ｉｎ Ｃｈｉｎｅｓｅ）） 【５］周学国，滕凯，邹伟．小流域设计洪水计算方法的简化 ．东北水利 设￡ ｔ ，建立并运行ｆｕｎｃｔｉｏｎ［Ｑｍ］，在指令窗口中 输入：ｆ１（１７５．４，１，０．００５，５０，２，０．７，１５０） 水电，２００１，１９（７）．（ＺＨＯＵ Ｘｕｅｇｕｏ，ＴＥＮＧ Ｋａｉ，ＺＯＵ Ｗｅｉ．Ｔｈｅ ｓｉｍ— ｐｌｉｉｆｅｄ ｍｅｔｈｏｄ ｆｏｒ ｓｍａｌｌ ｗａｔｅｒｓｈｅｄ ｄｅｓｉｇｎ ｃａｌｃｕｌａｔｉｏｎ ｍ．Ｗａｔｅｒ Ｒｅ— ｓｏｕｒｃｅｓ＆Ｈｙｄｒｏｐｏｗｅｒ ｏｆ Ｎｏｒｔｈｅａｓｔ Ｃｈｉｎａ，２００１，１９（７）．（ｉｎ Ｃｈｉｎｅｓｅ）） 求得Ｑｍ＝１６２３．６１４１ｍ％，ｔ＝１２．８０５７ｈ。ｔｃ＝１０６．８４２５ｈ， Ａｐｐｌｉｃａｔｉｏｎ ｏｆ Ｍａｔｌａｂ ｉｎ Ｒａｔｉｏｎａｌ Ｆｏｒｍｕｌａ ｆｏｒ Ｗａｔｅｒｓｈｅｄｓ ＴＩＡＮ Ｊｉｎｇｈｕａｎ，ＬＩＡＮＧ Ｗｅｎｔａｏ （ＮｏｒｔｈＣｈｉｎａＩｎｓｔｉｔｕｔｅ ｏｆＷａｔｅｒＣｏｎｓｅｒｖａｎｃｙａｎｄＨｙｄｒｏｅｌｅｃｔｒｉｃＰｏｗｅｒ，Ｚｈｅｎｇｚｈｏｕ４５００１１，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ａｓ ｆｏｒ ｔｈｅ ｒａｔｉｏｎａｌ ｆｏｒｍｕｌａ ｌｉｓｔｅｄ ｉｎ Ｒｅｇｕｌａｔｉｏｎ ｏｒｆ Ｃａｌｃｕｌａｔｉｎｇ Ｄｅｓｉｇｎｅｄ Ｆｌｏｏｄ Ｗａｔｅｒ Ｒｅｓｏｕｒｃｅｓ ａｎｄ Ｈｙｄｒｏｐｏｗｅｒ Ｐｒｏｊｅｃｔｓ ｏｒｆ ｄｅｒｉｖ・ ｉｎｇ ｓｍａｌｌ ｗａｔｅｒｓｈｅｄ ｆｌｏｏｄ，ｗｅ ｐｒｏｐｏｓｅｄ ｔｏ ｄｅｒｉｖｅ ｔｈｅ ｐｅａｋ ｄｉｓｃｈａｒｇｅ ａｎｄ ｆｌｏｗ ｃｏｎｃｅｎｔｒａｔｉｏｎ ｔｉｍｅ ｏｆ ｔｈｅ ｒａｔｉｏｎａｌ ｗｉｔｈ Ｍａｔｌａｂ ｓｏｆｔｗａｒｅ ｉｎ ｔｈｅ ｃｏｎｄｉｔｉｏｎ ｏｆ ｖａｒｉｏｕｓ ｋｉｎｄｓ ｏｆ ｆｌｏｗ ｃｏｎｃｅｎｔｒａｔｉｏｎ．Ｔｈｅ ｍｅｔｈｏｄ ｉｓ ｓｉｍｐｌｅ，ｐｒａｃｔｉｃａｌ，ａｎｄ ｅａｓｙ ｔｏ ｏｐｅｒａｔｅ ｗｉｔｈ ａｄｖａｎｔａｇｅｓ ｏｆ ｈｉｇｈ ｃａｌｃｕｌａｔｉｏｎ ｓｐｅｅｄ ａｎｄ ｐｒｅｃｉｓｉｏｎ． Ｋｅｙ ｗｏｒ￣：ｒａｔｉｏｎａｌ ｆｏｒｍｕｌａ；ｇｒａｐｈｉｃａｌ ｍｅｔｈｏｄ；Ｎｅｗｔｏｎ ｉｔｅｒａｔｉｏｎ ｍｅｔｈｏｄ；Ｍａｔｌａｂ
