皮尔逊III型曲线拟合与绘图：水文频率计算实战指南

简介：皮尔逊III型曲线是水文频率分析中常用的概率分布模型，可用于河流流量、降雨量等连续变量的拟合与设计值推求。这份RAR压缩包面向水文学、统计学及相关专业的科研人员与工程师，帮助解决P-III型分布参数估计、理论曲线绘制和不同频率流量插值计算等问题。包内共3个文件，均为MATLAB脚本（.m），分别对应参数拟合、频率插值与绘图环节，整体仅3KB，轻量易用。目前已有1302人学习浏览。借助其中的脚本，使用者可以直接读取实测数据进行参数求解，绘制皮尔逊III型概率密度曲线，并通过插值得到指定累积概率下的流量值，为防洪规划、水资源评估和风险决策提供量化支撑。脚本结构清晰，适合需要快速开展频率分析或教学演示的读者参考。

1. 皮尔逊III型曲线：水文频率计算里的那根“隐形拐杖”

做水文、市政排水和防洪影响评价的同行，大概都见过“皮尔逊III型曲线拟合与绘图.rar”这种压缩包工具：解压后就是一两个Excel或小软件，把年最大流量、年最大24小时暴雨量往里一贴，点个按钮，输出一条频率曲线和几个参数。皮尔逊III型曲线其实就是频率分析里最常用的理论线型，用几十年的观测值外推重现期10年、50年、100年对应的设计值。这套方法听着古老，却是工程审查和报告中少不了的“硬通货”。但多数人把工具当黑匣子用，填数据、看曲线、抄参数，很少有人较真P3曲线的参数到底合理不合理，坑也多半埋在这里。下面按我自己的实操顺序，把原理、拟合、绘图和踩过坑的地方拆开讲。

2. 从样本矩到P3曲线参数：均值、Cv、Cs的估计与选择

2.1 皮尔逊III型的密度函数与水文的三个参数

皮尔逊III型分布本质上是从伽马分布平移出来的偏态分布族，支持域有下界，右边拖尾，正好匹配洪峰流量、暴雨量这类“大部分年份小、偶尔出现极端大值”的样本。它的密度函数写成：

f(x) = β^α / Γ(α) · (x − δ)^(α−1) · e^(−β(x−δ))，x ≥ δ

其中α是形状参数，β是尺度参数，δ是位置参数。统计软件里用scipy.stats.gamma就能直接描述它，α和β决定形态和胖瘦，δ决定起点。

但水文计算不这么用。规范里和水文手册上，P3曲线的参数统一写成三个：样本均值 m̄、离势系数 Cv、偏态系数 Cs。这样做的好处是参数本身带有物理含义：m̄反映量级，Cv反映相对离散程度，Cs反映分布不对称程度。做设计时查《水文频率计算手册》里的离均系数 Φ 表，用的也是这三件套。

从统计参数转到水文参数的换算关系不复杂，Cs大于0时：

α = 4 / Cs² β = 2 / (m̄ · Cv · Cs) δ = m̄ · (1 − 2·Cv/Cs)

这三个式子很关键，后面写代码拟合分位数全得靠它们。洪水系列的Cs绝大多数为正偏，少数负偏场景（比如枯水流量分析）可以在代码里做镜像处理，这点到第3章再展开。

2.2 先算经验频率：样本排序是第一步

拟合曲线之前，每个观测值必须先有一个“频率位置”。以某水文站1980—2009年共30年的年最大洪峰流量为例，把流量从大到小排好，第m位的经验频率用数学期望公式算：

P_m = m / (n + 1)，用百分数表示就是 m/(n+1)×100%

这是国内用得最多的公式，也叫Weibull公式。为什么用n+1而不是n？直观上说，最大的那个样本不应该被当成“100%发生”，否则它就没有重现期概念了；分母取n+1能让最大值的经验频率落在n/(n+1)，给外推留出一点余量。这个细节看似小，但对尾部拟合影响很大。

如果样本里包含历史洪水调查值或特大值，经验频率公式要改成不连续系列计算方法，也就是把特大值与实测系列分段排频。实际操作时，我一般先直接把所有数据当作连续系列走一遍矩法和适线法，等确认趋势了再单独处理特大值。切忌一开始就把特大值硬塞进排序里，会让Cv被拉得非常高。

2.3 矩法初值：均值、Cv、Cs怎么算

矩法能给出P3曲线参数的初始估计，算法很直观。设样本为 x₁, x₂, …, xₙ，令 K_i = xᵢ / m̄，则：

m̄ = Σxᵢ / n Cv = sqrt( Σ(Kᵢ − 1)² / (n − 1) ) Cs = Σ(Kᵢ − 1)³ / (n · Cv³)

注意Cv用的是n−1，而Cs用的是n。这是因为Cv对应标准差的无偏修正，而三阶矩在样本量不够时很难谈无偏，工程上更习惯这个简化写法。常规的洪水样本算出来，Cv通常在0.3~1.5之间，Cs一般在Cv的2~4倍。如果你算出来的Cs超过了4倍Cv，或者出现负值，先别急着进入下一步，回去查原始数据里是不是有录入错误或异常极大值。

矩法算出的参数只能作为初值。特别是Cs，样本量少于30时它天生不稳定，一个特大值就能让它从1.5跳到3.0。所以后面必须接一道“适线”工序。

2.4 适线法：为什么还要人眼微调

适线法分目估适线和机助适线两种。目估适线是老工程师的做法：在频率格纸上点出经验点，先按矩法参数画一条理论频率曲线，再通过微调Cv、Cs（有时也动均值）让曲线尽量从点群中间穿过。机助适线则是把“曲线靠近点群”变成一个最小二乘问题，让计算机迭代三个参数。

这里有一个容易被忽略的原则：适当牺牲大频率段（比如50%、70%）的贴合度，优先保证小频率段（0.1%、1%）和中等频率段（5%、20%）的拟合精度。因为设计频率越稀遇，对曲线越敏感。机助适线如果不加权重，常常会把10年一遇以下的部分拟合得极好，百年一遇部分却偏得离谱。

所以我的习惯是：矩法出初值，机助最小二乘迭代，迭代完再看概率纸上的点群关系做最后手调。下面一章就把这套流程落到可运行的Python代码上。

3. 用Python跑通皮尔逊III拟合：最小脚本与参数调优

3.1 样本准备和经验频率计算

先准备一组典型样本。这里用某水文站15年的年最大洪峰流量序列，单位m³/s，排好序并计算经验频率：

import numpy as np import pandas as pd # 年最大洪峰流量，单位：m³/s data = np.array([2850, 1980, 1520, 1310, 1100, 860, 740, 620, 480, 350, 1120, 950, 680, 540, 410]) xs = np.sort(data)[::-1] # 从大到小排序 n = len(xs) emp_p = np.arange(1, n + 1) / (n + 1) # 数学期望公式 P = m/(n+1) print("排序后流量:", xs) print("经验频率:", emp_p)

经验频率数组emp_p表示每个排序值的超过概率，也就是“大于等于该值的年份比例”。排序方向和概率一一对应，后面拟合残差也要用同一套顺序。

pandas在这一步只用来预览数据，实际计算用numpy即可。数据量少时直接写numpy更直观，不容易引入DataFrame索引层面的错位。

3.2 P3分位数函数：正负偏态统一封装

算理论频率曲线需要一个能反推分位数的函数：给一个超过概率p，返回对应的设计值x_p。用SciPy的伽马分布分位数实现，并处理Cs为零或负值的情况：

from scipy import stats def p3_quantile(p, mean, cv, cs): """皮尔逊III型分位数函数，p 为超过概率（0~1）。""" if abs(cs) < 1e-6: # 偏态接近0时用正态近似，避免伽马形状参数爆炸 return stats.norm.isf(p, loc=mean, scale=mean * cv) if cs > 0: alpha = 4.0 / (cs ** 2) beta = 2.0 / (mean * cv * cs) delta = mean * (1.0 - 2.0 * cv / cs) return delta + stats.gamma.isf(p, alpha, scale=1.0 / beta) # 负偏：构造镜像变量求分位数 cs_abs = abs(cs) alpha = 4.0 / (cs_abs ** 2) beta = 2.0 / (mean * cv * cs_abs) delta = mean * (1.0 + 2.0 * cv / cs_abs) return 2.0 * mean - (delta + stats.gamma.ppf(p, alpha, scale=1.0 / beta))

正偏分支先按水文重参数化公式算出 α、β、δ，再调用gamma.isf(p, ...)，它表示“右侧尾概率等于p时的分位数”，正好对应“超过概率p的设计值”。负偏分支用镜像变量解决：让一个正偏变量 X′ 与原序列同均值，实际变量取 X = 2·m̄ − X′，偏度自动反转。

这段代码是后面所有拟合和绘图的基础。手动在Excel里查表时，查的其实也是同一套离均系数关系，只是这里由计算机直接算出来。

3.3 最小二乘拟合三个参数：least_squares实现

矩法给出初值后，用scipy.optimize.least_squares做机助适线。目标函数比较“理论分位数与实测排序值”的差异，残差采用相对误差：

from scipy.optimize import least_squares def fit_residual(params, xs): mean, cv, cs = params if cv <= 0: return np.full(len(xs), 1e6) emp_p = np.arange(1, len(xs) + 1) / (len(xs) + 1) x_fit = np.array([p3_quantile(p, mean, cv, cs) for p in emp_p]) return (x_fit - xs) / xs # 相对残差 mean0 = xs.mean() cv0 = xs.std(ddof=1) / mean0 cs0 = np.mean(((xs - mean0) / (xs.std(ddof=1))) ** 3) # 有偏样本偏度 bounds = ([mean0 * 0.5, 0.02, 0.0], [mean0 * 1.5, 2.0, 5.0]) res = least_squares(fit_residual, x0=[mean0, cv0, cs0], bounds=bounds) mean_f, cv_f, cs_f = res.x print(f"拟合结果: m={mean_f:.1f}, Cv={cv_f:.3f}, Cs={cs_f:.3f}")

残差除以实测值，等价于对数尺度上的近似，能避免大流量观测值在目标函数里权重过高。边界约束把Cs限制在0~5之间，Cv限制为0.02~2.0，均值允许在初值0.5~1.5倍之间浮动。做设计洪水时，均值通常不会偏离样本均值太远，这个bounds设置合理。

注意CS的下界设为0.0，因为洪峰流量理论上不应负偏。如果样本算出负偏度，要么是数据里混合了不同类型洪水，要么是历史特大值权重不够，直接约束成正值反而能暴露问题。

3.4 用拟合参数反推设计值

有了三个拟合参数，百年一遇、五十年一遇就能直接算：

qs = [0.01, 0.02, 0.05, 0.10] # 超过概率 design = {f"p={p:.2f}": p3_quantile(p, mean_f, cv_f, cs_f) for p in qs} for k, v in design.items(): print(k, f"{v:.0f} m³/s")

p=0.01表示超过概率1%，即百年一遇；p=0.02是五十年一遇。算出的设计值如果比样本最大值高出一大截，先核对cs参数是不是过大，再回头检查样本长度。

最小二乘结果只能说明“拟合线离点群近”，不能保证外推合理性。水文规范里常说“外推不超过样本长度的2~3倍”，这是经验法则，也是审查专家盯的重点。你的样本只有15年，报一个500年一遇设计值给下游工程用，从统计上没人能反驳，但从风险控制角度极不负责。

4. 皮尔逊III曲线绘图：概率格纸与四条出图路径

4.1 概率格纸：为什么横轴要“扭一下”

皮尔逊III频率曲线通常不画在线性概率横轴上。因为频率从50%到1%跨越太大，线性轴上曲线尾部会被压成一条竖线，根本看不出拟合好坏。频率格纸的做法是把横轴按正态分位数重排：让常用频率刻度在图上等距展开，P-III曲线在格纸上接近一条直线，工程上来回直线滑动标定非常方便。

常见的换算关系是：横轴坐标 =norm.isf(P)，其中P为超过概率。几个锚点如下：

图上概率值从右往左递减，最左边是0.01%，最右边是99%。很多初看P-III绘图的人会对着横轴方向犯迷糊，记住“坐标值越大代表频率越小”就不容易乱。

4.2 用Matplotlib绘制规范频率图

把理论曲线和经验点画在同一张概率格纸上，这是科研绘图场景下最常用的做法。代码可以直接接上一章的拟合结果：

import matplotlib.pyplot as plt # 生成理论频率曲线在概率格纸上的坐标 freq_p = np.array([0.01, 0.05, 0.1, 0.5, 1, 2, 5, 10, 20, 50, 75, 90, 99]) freq_p = freq_p / 100.0 x_theory = np.array([p3_quantile(p, mean_f, cv_f, cs_f) for p in freq_p]) x_pos = stats.norm.isf(freq_p) # 横轴坐标 fig, ax = plt.subplots(figsize=(8, 5)) ax.plot(x_pos, x_theory, 'r-', lw=2, label='P-III理论曲线') # 经验频率点 emp_p = np.arange(1, n + 1) / (n + 1) x_emp = np.sort(data)[::-1] ax.scatter(stats.norm.isf(emp_p), x_emp, s=30, facecolors='none', edgecolors='k', label='经验点') # 标出常用频率刻度 ticks_p = np.array([0.05, 0.1, 0.5, 1, 2, 5, 10, 20, 50, 80, 90, 95, 99]) ax.set_xticks(stats.norm.isf(ticks_p / 100.0)) ax.set_xticklabels([f"{t:.2g}%" for t in ticks_p]) ax.set_xlabel("频率 P(%)") ax.set_ylabel("洪峰流量(m³/s)") ax.legend() plt.show()

scipy.stats.norm.isf把超过概率P转成标准正态分布右侧分位数，直接用这个值当横轴坐标，刻度标签仍写回百分比。这样图上曲线被“拉直”，经验点偏离曲线多远一目了然。

图中的红线是理论频率曲线，黑色空心点是经验点。如果一个点在尾部明显凸出来，往往是样本里有特大值或Cs参数偏低；如果所有点都在曲线同一侧，说明均值或Cv整体偏离。绘图不是为了好看，是为了诊断。

4.3 用Origin复现并微调曲线：审查汇报更顺手

工程量大的单位里，报告插图还是以Origin为主。把Python算出来的设计值表复制进Origin工作表，A列填频率百分比，B列填设计值，然后用Plot菜单的Symbol画经验点、Line画理论线。关键一步是双击横轴，在Scale选项卡把Type改成Probability，Origin会自动按概率刻度展示。

如果不想在Origin里重算分布，就用笨办法：把第3章算出的多条频率-设计值序列全部粘进去，画成多站点对比图。做防洪影响评价时，不同频率下的设计值常常要放进同一张报告图，Origin的多图层组合功能比Matplotlib省事。

唯一的坑是Origin的概率坐标轴默认可能从左向右递增，而水文频率格纸习惯从右向左。在Axis选项卡里勾上Reverse即可。

4.4 其他可替代方案：Matlab、Excel与网页端绘图

Matlab同样可以完成整套拟合，用fitdist做分布拟合或直接用nlinfit做自定义最小二乘，画图用probplot。如果以前写过Matlab脚本，半路上车并不难。不过Matlab的统计工具箱在一些旧版本里没有直接支持P-III分布，需要自己写分位数函数，本质还是本章的代码逻辑。

Excel本身适合快速填数和存档，但它的图表类型没有天然概率轴，需要自己用NORM.S.INV(1-P)生成一列虚拟横轴坐标，再配合散点图。这个方法适合给外单位看过程数据，不适合当正式成果图。

如果要把多站点P-III曲线做到内网系统或Web页面上，Qt绘图和Canvas绘图都是值得考虑的方向。桌面端用Qt的QCustomPlot控件可以自定义轴刻度，浏览器端可以直接用ECharts这类开源绘图库，把经过norm.isf变换后的横轴坐标传给折线图引擎就行。前端绘图的重点是不要在浏览器里重新实现伽马分位数计算，后台算好坐标再传JSON最省事。

5. 皮尔逊III拟合避坑指南：4个翻车场景与参数保护

5.1 样本太短，Cs虚高：15个点拟合的百年一遇别急着报出去

现象：某小型水库入库流量只有12年实测资料，矩法算出Cv=1.1，Cs=3.2。机助适线后曲线看起来“完美”穿过经验点，但外推百年一遇洪水达到样本最大值的2.8倍，专家质疑后项目停工补资料。

原因：样本量不足时，三阶矩估计的方差极大。一个特大值年份就能把Cs从1.5推到3以上，而这恰巧是P-III曲线尾部抬升最敏感的参数。机助适线只保证“拟合残差小”，不保证外推稳定。

解决：优先延长系列长度，把相邻水文站的同期资料和相关分析结果用上；如果实在无法延长，就用地区综合法，参考《水文手册》中同类气候区的Cv、Cs比值，把初值锁定在合理范围再做适线。机助拟合的边界也不能放开，Cs上限先按4倍Cv设，宁可损失一点尾部贴合度，也不能让曲线末端失去控制。

5.2 优化器把Cs迭代成怪异值：约束必须写进目标函数

现象：用least_squares拟合时没加bounds，某次迭代后Cs变成7.2，曲线在0.01%频率处飙到几万，直接击穿物理常识。代码却显示“优化成功”。

原因：优化器只认残差，不认物理。尾部几个极稀疏的经验点对残差的拉动作用，在某些样本分布下会被放大，参数因此被推到数学可行但物理荒谬的区域。

解决：把参数约束写进优化调用，而不是等结果出来再判断。实际我会在目标函数里再加一道“软约束”：Cs > 4*Cv时在残差上增加惩罚项。硬边界加软惩罚双保险，曲线不会突然失控。每次拟合后还要打印出Cv/Cs的比值，如果超过5，默认这次拟合不可用，回头检查数据和初值。

5.3 概率轴横轴方向不一致：两张图看起来完全不像同一条线

现象：同一个样本，A用Python按线性频率画图，得到一条向上凹的曲线；B用Origin概率轴画图，曲线却接近直线。两人把图放在一起比对，外推设计值明明一致，但视觉差异让甲方误以为有人算错了。

原因：P-III曲线在线性频率轴上本来就呈急弯，这是分布形态决定的，不代表拟合错误。不同人用不同轴类型，图形习惯差异掩盖了结果一致性。

解决：内部技术交流明确规定，正式成果图一律使用概率格纸，横轴方向为“频率小值在右、大值在左”；给外部单位提交图件时在图注里写清坐标轴含义。画图代码里要统一封装坐标转换函数，不要这个脚本写“线性”，另一个脚本写“概率”，回头对比时才知道坑多深。

5.4 老式rar工具在Windows上的经典问题：路径、宏与杀毒误报

现象：从网上下载的“皮尔逊III型曲线拟合与绘图.rar”解压后，运行一个名为“计算绘图.xlsm”的文件，提示宏被禁用或按钮无反应；有时解压工具还会隔离“不知名可执行文件”。

原因：这类老工具很多是用Excel VBA写的，引用的ActiveX控件在新版Office里被默认禁用；部分压缩包被第三方改过结构，杀毒软件会对宏文件误报。解压到中文带空格路径下，VBA内部读取相对路径也会失效。

解决：先把压缩包解压到纯英文路径，比如D:\P3_tool，不要直接放到桌面或中文目录；再在Excel里打开“文件-选项-信任中心宏设置”，启用宏后重新打开。被杀毒隔离的文件，先恢复到隔离区并在沙箱环境单独跑一遍，确认行为正常再使用。标题这类rar工具包往往版本老旧，指望它处理现代大数据量样本并不现实，更适合用来做参数对照和快速试算。

6. 验证拟合结果与批量出图：KS检验和多站自动处理

6.1 用KS检验判断拟合优度

机助适线之后，光看曲线不够，做一个非参数检验能确认参数是否可接受。Kolmogorov-Smirnov检验比较理论分布和经验分布的累计概率差异，指标不依赖于绘图习惯：

def p3_cdf(x, mean, cv, cs): if abs(cs) < 1e-6: return stats.norm.cdf(x, loc=mean, scale=mean * cv) if cs > 0: alpha = 4.0 / (cs ** 2) beta = 2.0 / (mean * cv * cs) delta = mean * (1.0 - 2.0 * cv / cs) return stats.gamma.cdf(x, alpha, loc=delta, scale=1.0 / beta) cs_abs = abs(cs) alpha = 4.0 / (cs_abs ** 2) beta = 2.0 / (mean * cv * cs_abs) delta = mean * (1.0 + 2.0 * cv / cs_abs) return stats.gamma.sf(2.0 * mean - x, alpha, loc=delta, scale=1.0 / beta) ks_stat, ks_p = stats.kstest(xs, lambda t: p3_cdf(t, mean_f, cv_f, cs_f)) print(f"KS统计量={ks_stat:.3f}, p值={ks_p:.3f}")

p值大于0.05说明样本与P-III分布没有显著偏离，可以作为拟合可用的依据之一。但KS检验对尾部不敏感，所以最终判定还是要结合0.1%、1%这些设计频率处的相对偏差是否在工程可接受范围内。

6.2 批量拟合多个站点并自动出图

实际项目里往往是几十个断面一起做频率分析。逐一手动拟合再导出参数，既慢又容易出错。整理一个统一入口批量处理CSV文件是最划算的投入：

import glob from pathlib import Path results = [] for csv_path in glob.glob("data/*.csv"): xs = np.genfromtxt(csv_path, delimiter=",", skip_header=1) xs = np.sort(xs)[::-1] mean0 = xs.mean() cv0 = xs.std(ddof=1) / mean0 cs0 = np.mean(((xs - mean0) / (xs.std(ddof=1))) ** 3) res = least_squares(fit_residual, x0=[mean0, cv0, cs0], bounds=bounds) results.append([Path(csv_path).stem, *res.x]) # 按 4.2 节代码生成 PNG 图片，文件名用站名 print("站点,均值,Cv,Cs") for row in results: print(",".join(f"{v:.3f}" for v in row))

批量跑完先看参数矩阵，再抽查概率图。如果某几个站点的Cv或Cs明显偏离区域规律，优先补资料或复核原始数据，而不是直接改参数“把曲线掰好看”。

我自己的习惯是每批拟合结果都留一个JSON参数文件，里面记录样本年份、矩法初值、机助适线终值、KS p值。日后复核或做地区综合时，回查起来快得多。这条习惯也救过我一次：前期把Cs参数写错了一位小数，幸好留有中间结果文件，对比时直接定位到源头，否则整组设计值都得重算。这套皮尔逊III型曲线拟合与绘图的流程，看起来是行老手艺，但细节处仍然需要较真。希望帮到你。

本文还有配套的精品资源，点击获取