简介面向水文工作者的频率分析工具包围绕读取数据、预处理、分布拟合、频率计算全流程设计适用于洪水、干旱等极端事件概率评估学习或实际项目使用。资源包内共十四个文件压缩包大小约1011KB包含可执行程序、动态链接库、样例数据表、txt说明文档及chm帮助文件可支撑运行和查阅完整的水文频率分析操作流程。已有四千三百七十四人学习或下载该资料。读者可在已有工具基础上结合自身降雨或径流数据完成皮尔逊三型等分布拟合利用交叉验证评估模型效果进而得到洪水重现期对应的特征值资源内置曲线拟合示例便于快速上手有助于理解计算机方法在水文统计中的实际应用对开展水文分析工作有参考价值。1. 水文频率分析中的 cv 与 cs决定设计洪水大小的两个隐形旋钮做水文频率分析的人几乎都跟 cv、cs 这两个参数打过交道。cv 叫变差系数管的是数据序列的离散程度cs 叫偏态系数管的是分布的不对称性。很多人在推求设计洪水时把精力全花在延长序列、考证历史洪水上却忽略了这两个参数才是真正决定“百年一遇”设计值高低的隐形旋钮。我见过某水库除险加固项目矩法初算出来 cv0.42、cs0.30怎么适线都对不上高水端点据最后把 cs 调到 1.20设计洪水一下子就抬高了近两成整个溢洪道复核结论跟着改变。cv、cs 不光是统计参数更直接关系到工程规模、投资和安全值得把它们的估计逻辑、调整规律和坑摸清楚。这篇笔记适合水利、水文、市政排水领域的从业者也适合正在做毕业设计的在校生。2. 先立住概念cv 与 cs 在 P-III 型分布里到底管什么2.1 从矩估计说起cv、cs 的统计定义与计算方法水文频率分析里cv 和 cs 都来自样本矩。设某站年最大洪峰流量序列为 x1, x2, …, xn均值为 x̄标准差为 σ那么 cv 的定义式是Cv σ / x̄它反映的是序列的相对离散程度。Cv 越大说明各年洪水之间的丰枯差异越大频率曲线就越陡同一个重现期下推出来的设计值也越高。我国大多数河流的年最大洪峰流量序列Cv 在 0.2 到 1.5 之间北方干旱半干旱地区往往比南方湿润地区大。cs 的定义式是Cs [ Σ(xi - x̄)³ / n ] / σ³它反映的是分布对均值的偏斜方向和程度。水文变量普遍是正偏的也就是存在少数年份出现极大洪水把分布右侧拖出一条长尾所以 cs 一般大于 0。cs 越大频率曲线在高水端翘得越高尾部设计值增长越快。注意水文上习惯用 n 做分母计算 cs这和统计学里常见的无偏修正版本不同。样本容量不大时这个差异会导致 cs 的系统性偏差后面避坑章节会专门展开。2.2 为什么水文规范普遍用 P-III 型分布承载 cv、cs做了多年频率分析的人都知道可选的理论分布很多比如对数正态分布、广义极值分布GEV国外水文界用 GEV 的越来越多。但我国水利行业的频率分析规范长期以皮尔逊 III 型分布P-III 型为主流。P-III 型分布的概率密度函数由均值、Cv、Cs 三个参数唯一确定它本质上是一种带偏态的伽马分布族能很好地拟合洪水、暴雨这类变量的正偏长尾特征。P-III 型分布之所以在水文界站稳脚跟不是因为它在数学上最优而是因为它的偏态形状参数恰好对应 cs让工程师能直观地把“偏态大小”和“曲线形状”对应起来。而且几十年下来各地区已经积累了海量的 Cv、Cs 经验关系比如很多地方 Cs ≈ (2~4)Cv。这些经验值是拿真金白银的工程实践换来的换一种分布族这些经验就全部作废这是 P-III 型分布的强大惯性。2.3 样本矩的天然缺陷小样本下 cv、cs 为什么不能直接用频率分析的痛点在于实测洪水序列往往只有 20 到 50 年用短序列直接算出来的矩统计误差很大。cv 的矩估计在样本容量小于 30 时系统性偏高cs 更是离谱样本容量小于 30 时cs 的方差极大算出来的值可能从负的到正的翻几倍都不奇怪。这就是为什么水文频率分析不能停留在“矩法算完就交差”而要走“矩法初估 适线修正”的路子。适线法本质上就是用经验频率点据去约束理论频率曲线让理论曲线的形状参数在点据的“引力”下反复调整直到曲线与点据整体贴合。理解了这个逻辑后面调参才不会像无头苍蝇乱撞。3. 用矩法初估 cv、cs最小可复现的计算流程3.1 用 Python 计算样本均值、cv、cs 的完整流程先给出一段可以直接跑通的 Python 代码。数据我用一组模拟的年最大洪峰流量序列代替你自己跑的时候替换成实际序列即可。这段代码计算均值、均方差、cv 和 cs并把关键中间量都打印出来方便核对中间过程。import numpy as np # 模拟一组年最大洪峰流量序列单位 m^3/s # 实际使用时替换为你自己的实测序列 q np.array([1250, 890, 1420, 2100, 680, 1150, 950, 1680, 1300, 780, 2200, 1080, 920, 1450]) n len(q) mean_q np.mean(q) # 样本均值 std_q np.std(q, ddof1) # 样本标准差水文上也有用 ddof0 的见说明 cv std_q / mean_q # 变差系数 # 计算偏态系数这里采用水文上常见的 sigma 全序列求和的写法 sigma np.sqrt(np.sum((q - mean_q) ** 2) / n) cs np.sum((q - mean_q) ** 3) / (n * sigma ** 3) print(样本容量 n , n) print(均值 x_mean %.2f % mean_q) print(标准差 sigma %.2f % std_q) print(Cv %.4f % cv) print(Cs %.4f % cs)这段代码的用意是把矩估计的每一步显式算出来。注意标准差那里我用了 ddof1这是统计学上的无偏估计但后面算 cs 时分子里的离差三次方和除以 n用的又是水文上的常见口径。严格来说cv 与 cs 的口径应该保持一致有些规范推荐 cv 也用全序列 σ 计算这样 cv 和 cs 才能在同一套矩体系下配合 P-III 曲线。实操中我建议先按规范口径统一再来回比较。3.2 从样本矩到频率曲线查出离均系数 Φ 并绘制理论曲线矩法得到 Cv、Cs 初值之后就要把理论频率曲线画出来。P-III 型分布的频率曲线通常用离均系数 Φ 表示。给定一个频率 P理论设计值 xp 为xp x̄ × (1 Cv × Φ(P, Cs))这里的 Φ 可以从规范附表的 P-III 型离均系数表查得也可以用 Python 里的 scipy.stats.pearson3 分布直接算。下面这段代码演示了用 scipy 计算理论频率曲线并与经验频率点据做对照的思路。import numpy as np from scipy import stats # 沿用上一节算出的 x_mean、cv、cs 初值 x_mean mean_q cv_init cv cs_init cs sigma_init x_mean * cv_init # 矩法下标准差与均值、cv 的关系 # 构造 P-III 分布对象skew 是偏态系数loc 是均值scale 是标准差 dist stats.pearson3(skewcs_init, locx_mean, scalesigma_init) # 定义一组频率水文上常用百分比表示 freqs_pct np.array([0.01, 0.1, 0.5, 1, 2, 5, 10, 20, 50, 75, 90, 95]) freqs freqs_pct / 100.0 # 对 exceedance 概率 P频率设计值对应的是 1-P 分位数 design_values dist.ppf(1 - freqs) for p, v in zip(freqs_pct, design_values): print(P%.2f%% design%.2f m^3/s % (p, v))这里有两个容易搞混的参数说明。第一scipy 的 pearson3 用 skew 参数表示偏态系数正偏对应正值与水文习惯一致。第二ppf是分位点函数输入的是累计概率水文上说的“频率 P1%”意味着超过概率为 1%对应累计概率就是 99%所以代码里用的是ppf(1 - freqs)。这一行写反的人特别多算出来的设计值整体错位务必检查。3.3 经验频率点据计算把实测值画到概率格纸上理论曲线有了还得有实测点据去对照。经验频率最常用的公式是数学期望公式P m / (n 1)其中 m 是实测值从大到小排列的序号n 是样本容量。也就是说最大的洪水对应 m1经验频率为 1/(n1)最小的洪水对应 mn经验频率为 n/(n1)。下面这段代码把点据频率与实测值组织成表格数据。import numpy as np import pandas as pd # 仍用前面的序列 q sorted_q np.sort(q)[::-1] # 从大到小排列 n len(sorted_q) m np.arange(1, n 1) # 序号 1,2,...,n empirical_p m / (n 1) # 经验频率用数学期望公式 df pd.DataFrame({ 序号: m, 实测洪峰(m^3/s): sorted_q, 经验频率P: empirical_p, 重现期T(年): 1 / empirical_p }) print(df)经验频率公式的选择对适线结果影响很大。数学期望公式是无偏的但最大值的经验频率等于 1/(n1)意味着实测最大值被赋予了大约 n 年一遇的频率这在样本短时会导致高水端点据“偏矮”。有些地区规范推荐用中位秩公式或 Weibull 公式P 值会略有差异。我一般初算统一用数学期望公式到了复核阶段再按工程等级和规范要求换公式做敏感性分析。3.4 目估适线的六步操作套路矩法初估完、点据画好之后真正的“手工活”是目估适线。常见做法分六步走第一步把经验频率点据绘制在概率格纸上纵坐标为变量值横坐标为概率尺度。第二步把矩法初估的 Cv、Cs 代入 P-III 分布计算理论频率曲线并叠加在图上。第三步检查曲线与点据的偏离情况优先看高水端P0.01%~1%和中部P10%~50%。第四步固定均值 x̄ 不变小幅调整 Cv让曲线陡峭程度贴近点据整体趋势。第五步固定 Cv 后调整 Cs让曲线高水端翘曲程度匹配最大几个点据。第六步反复迭代直到曲线与点据没有系统性的单向偏离再读取设计频率对应的设计值。这套流程听起来简单但实际操作时“小幅调整”到底调多少全凭经验。我的做法是每次调整幅度不超过初值的 10%比如 Cv 初值 0.50一次只调 0.02 到 0.05Cs 同理。调完一轮画一次图而不是在参数空间里大步乱跳否则很容易陷入“怎么调都合不上”的死循环。4. 让 cv、cs“听话”调参规律与设计值合理性校验4.1 cv、cs 联动调整时曲线形状怎么变先单独看 cv。其他参数不变cv 增大频率曲线整体变陡即大频率端低水段数值下降小频率端高水段数值抬升。直观理解就是年际变化越大同样百年一遇的洪水就越夸张。再看 cs。其他参数不变cs 增大频率曲线低水端抬高、高水端也抬高但中间段P20%~70%会被压低整个曲线变得更“歪”。这是因为 cs 增大意味着分布把更多概率质量推向左侧而右侧长尾的极端分位数会显著变大。也就是说cs 对设计值的影响在 P1% 以下时比 cv 还猛这也是前面那个水库案例里 cs 从 0.30 调到 1.20 后设计洪水大涨的原因。实际操作中cv 和 cs 是联动的不能独立乱调。经验规律上是先定均值、再定 Cv最后用 Cs 作微调。Cv 定得过低后面用 Cs 怎么拉都拉不出高水端的翘曲Cv 定得过高曲线会从低水端开始就整体偏离点据。所以调参顺序本身就是一种约束。4.2 设计值合理性校验让曲线外延经得起推敲适线完成不等于工作结束设计值的合理性必须独立检验。常用的校验手段是表格化的对比把相邻几个频率的设计值放在一起看递增规律是否平稳。频率 P(%) 设计值(m^3/s) 相邻倍比 0.01 5980 1.25 0.1 4780 1.21 1 3950 1.18 5 3350 1.13 10 2960 1.11 20 2450 1.09相邻频率设计值的倍比应该随着频率降低而逐渐趋近于 1如果出现 P0.1% 到 P0.01% 之间倍比突然大于 1.5 的情况大概率是 cs 取得太大尾部外延失控了。另一个校验方式是把推求的千年一遇、万年一遇设计值与流域内历史调查洪水对比如果万年一遇设计值还赶不上调查到的历史洪水下限说明参数整体偏低反过来如果百年一遇就远超历史最大洪水几个量级也要怀疑参数是否偏高。4.3 地区综合经验值什么时候能用、什么时候必须放弃各省水文手册和水文图集里基本都给了分区的 Cv 等值线图以及 Cs/Cv 比值的经验范围。对于缺乏资料的工程地址用地区综合值插值是合理且必要的。但地区综合值有两个使用前提一是设计流域与周边参证站的流域特征面积、河长、下垫面不能差太远二是序列中没有经过考证的历史特大洪水。一旦序列里加入了历史特大洪水比如某站 1880 年调查洪水明显比实测期最大值还大这时地区经验值往往会失效。原因很简单地区综合值是基于大量一般站点统计出来的无法解释某个站点独有的特大值。这种情况下要以本站的洪水考证成果为主把特大值作为不连序样本纳入矩估计而不是强行套用 Cs/Cv 比值。5. 水文频率分析避坑cv、cs 相关的五个典型翻车现场5.1 矩法算出的 cs 是负的或接近 0曲线完全不合点据现象直接用矩法算出 Cs-0.05P-III 曲线几乎对称高水端点据明显上翘曲线怎么都压不住那几个大洪水点。原因样本容量小于 30 时Cs 的矩估计方差非常大符号都可能搞反。尤其序列里只有一两个突出大值时三阶矩的分母被极端值放大算出来的 Cs 反而偏低甚至为负。解决放弃矩法 Cs 初值改用地区经验 Cs/Cv 比值作为初值比如先设 Cs3Cv然后再适线调整。同时回头检查特大值是否已经单独处理过如果最大洪水没有作为历史特大值考证而是混在实测序列里矩法 Cs 被污染几乎是必然的。5.2 cs 越调越大尾部设计值涨到离谱现象适线时为了贴住最大那个点据Cs 一路调到 5.0 以上曲线高水端剧烈上翘万年一遇设计值比千年一遇大了近一倍审查时被一眼看出不合理。原因单个特大值对适线产生了“引力”模型用 Cs 的极端增大去拟合这个点牺牲了整个曲线尾部形态。这在点据稀少、高水端只有一两个点时特别容易发生属于典型的过拟合。解决先复核特大值的重现期是否准确如果历史洪水重现期考证错误整个序列的框架就错了。其次给 Cs 设一个工程经验上限我一般控制在 3.0~3.5 之间超过这个值就该怀疑数据或模型选择有问题而不是继续无脑调参。5.3 适线只盯高水端低水端点据整体跑偏现象适线后高水端几个大点贴合很好但 P50%~90% 的低水端点据全部位于曲线上方或下方形成系统性偏离。原因目估适线时人眼天然被大数值吸引潜意识里优先照顾高水端导致曲线中低段被牺牲。低水端对应的是常遇洪水、中小洪水直接影响中小型工程的设计洪水位不能忽视。解决把点据按频率分成高水段P10%、中水段10%~50%、低水段P50%三段适线时分别观察。必要时引入加权最小二乘的量化指标给低水段更高的权重让目估有据可依而不是全凭感觉。5.4 同一条序列不同人适线的结果相差很大cs 能差一倍现象同一组资料A 同学适线得到 Cv0.55、Cs1.50另一位同事得到 Cv0.60、Cs2.80两者算出的百年一遇设计值差 12% 以上。原因目估适线主观性太强。点据少、曲线在 P0.01%~0.1% 区间基本上没有实测点约束外部形态全靠 Cs 在撑不同的初始值会收敛到不同的局部“满意解”。解决固定矩法 Cv 初值不动Cs 只允许在地区经验比值附近做小范围调整调整幅度控制在初值的 ±50% 以内。这样一来不同人的适线结果不至于差出数量级。最终提交成果必须附上参数取值依据和适线图让审查方能看到调整路径。5.5 加入历史洪水前后cv、cs 剧烈变化现象序列原来没有考虑历史洪水时Cv0.48、Cs1.10把 1956 年调查洪水按 150 年重现期加入后Cv 变成 0.65Cs 变成 2.30设计值全面重排。原因历史特大值是一个超高值加入序列后均值、二阶矩、三阶矩全部被拉动尤其对 Cs 的拉动效应是几何级的所以参数变化剧烈是正常现象不是算错了。解决处理历史洪水时一定要用不连序样本的矩公式而不是简单地把特大值混进实测序列里一起算。同时对历史洪水的重现期做敏感性分析比如分别按 100 年、150 年、200 年重算看设计值得变化幅度能否接受。如果某个重现期假设导致结果大幅波动说明资料的考证精度不足以支撑当前结论需要补充调查或多方案比选。6. 把 cv、cs 的估计从手工试错变成可复算的优化适线手工目估适线的结果很难复算换个人又不一样。这里分享一个我常用的收尾技巧把“曲线与点据的离差平方和最小化”变成一个优化问题借助计算机自动搜索 cv、cs 的最优组合再把优化结果作为人工终审的起点。import numpy as np from scipy import stats from scipy.optimize import minimize def fitting_loss(params, x_mean, sorted_q, empirical_p): # params 是 [cv, cs] cv, cs params sigma x_mean * cv dist stats.pearson3(skewcs, locx_mean, scalesigma) theo_p dist.ppf(1 - empirical_p) # 理论设计值 # 用相对离差平方和避免大数值点主导目标函数 return np.mean(((theo_p - sorted_q) / sorted_q) ** 2) # 初值取矩法结果边界约束在合理工程范围内 x0 np.array([cv, cs]) bounds [(0.1, 2.0), (0.1, 4.0)] res minimize(fitting_loss, x0, args(x_mean, sorted_q, empirical_p), methodL-BFGS-B, boundsbounds) print(优化Cv %.4f, 优化Cs %.4f % (res.x[0], res.x[1]))这段代码里我刻意用了相对离差而不是绝对离差原因是洪水数值跨度太大绝对离差会让最大洪水点说了算低水端点据等于被忽略。目标函数括号里除以实测值相当于把每个点的误差百分比均匀看待。优化出来的 Cv、Cs 永远只能作为“初稿”最终要回到概率格纸上人工审视高水端外延是否合理。这已经是我做频率分析的习惯动作了机器搜索定初稿人工判断定终稿两条腿走路。希望这篇整理对你有用也希望你在自己的流域资料上跑一遍拿真实点据去验证 cv、cs 的取值逻辑踩过一轮坑之后你会对这两个参数有完全不一样的感觉。本文还有配套的精品资源点击获取