简介：皮尔逊III型曲线拟合与绘图资源面向水文、气象及统计领域的科研与工程人员，解决连续偏态随机变量的频率分析问题，尤其适合处理河流流量、降雨量等自然观测数据。资源包含3个MATLAB脚本文件（.m），压缩包整体仅3KB，文件虽精简但分工明确，覆盖曲线计算、参数求解和图形绘制等核心环节，可支撑从原始数据到可视化结果的全流程操作。目前已有1303人学习下载，常被用于水资源评估、防洪设计等工作场景。利用该资源可基于最大似然估计获取形状参数κ、尺度参数θ、位置参数μ，理解数据对称性与偏斜程度，绘制概率密度曲线，并沿曲线插值得到累积概率1%、5%等对应流量，直接为风险决策与工程规划提供量化依据。脚本结构简洁、接口清晰，便于结合具体项目二次修改，也可作为教学演示模块嵌入既有分析框架。
1. 皮尔逊III型曲线：水文频率分析绕不开的那个"偏态分布"
做设计洪水、暴雨强度公式、年径流频率分析的工程师，抽屉里大概率躺着一个叫"皮尔逊III"的老压缩包，标题里那个 .rar 就是这类东西的典型：一套把实测点据拟合成皮尔逊III型理论曲线、再画到概率格纸上的程序，doingni2 是里面做配线迭代的核心函数，p3 曲线的参数就是均值、变差系数 Cv、偏态系数 Cs 这三个数。它能解决什么问题？把零星的历史流量或雨量观测值，外推到百年一遇、千年一遇的设计值——这是水利、市政、桥涵设计里绕不开的环节，也是行业规范推荐的线型。适合谁？做水文分析、水工设计、排水防涝计算的从业者，以及被课程设计里的频率分析折磨的在校生。下面这套流程，我按老程序包的逻辑用 Python 重写了一遍，照着跑就能复现，不用纠结原包装的是 Fortran 还是别的。
2. 从密度函数到三个统计参数：皮尔逊III型凭什么被水文界选中
2.1 密度函数与三参数的关系
皮尔逊III型曲线本质上是带位置参数的 Gamma 分布，密度函数写成：
f(x) = β^α / Γ(α) · (x - a₀)^(α-1) · e^(-β(x - a₀))
自变量 x 必须大于下限 a₀，否则密度无定义。水文上不直接使用 α、β、a₀ 这三个形式参数，而是重新参数化成均值 x̄、变差系数 Cv、偏态系数 Cs，换算关系是：
α = 4 / Cs²
β = 2 / (x̄ · Cv · Cs)
a₀ = x̄ · (1 - 2Cv / Cs)
这套换算关系是所有后续计算的地基。注意看第一个式子：Cs 越大，α 越小，分布形态变化越剧烈；当 Cs = 2Cv 时 a₀ = 0，曲线下限恰好落在原点；Cs 超过 2Cv 后 a₀ 变负，理论下限延伸到负值，这在洪水分析里可以接受，但如果 Cs 逼近甚至超过 4Cv，曲线在小概率端会明显上翘，配线很容易翻车。所以老工程师反复叮嘱：Cs 与 Cv 的倍比关系不要乱调，行业习惯是 Cs ≈ (2~4)Cv，能不动就不动。
水文界选中 P-III 型不是因为它数学上最优，而是因为三点务实理由：一是水文极值系列一般是正偏的，尾部拖得长，正态分布描述不了这种"大值离群"的现象；二是它有明确下限 a₀，贴近"河道流量不会小于某个值"的物理直觉；三是它只有三个参数，实测系列哪怕只有二三十年，也能勉强拟合出相对稳定的结果，换成四参数、五参数的分布，短系列下参数互相抵消，算出来的设计值反而更不靠谱。
2.2 矩法初估三参数的具体计算
配线的第一步从来不是直接画曲线，而是先把三个参数的初始值估出来。最常见的做法是矩法：用一阶原点矩估均值，二阶中心矩估方差，三阶中心矩估偏态，再反推 Cv 和 Cs。公式是：
x̄ = (1/n) Σxᵢ
Cv = sqrt( Σ(Kᵢ - 1)² / (n - 1) )，其中 Kᵢ = xᵢ / x̄
Cs = Σ(Kᵢ - 1)³ / ((n - 1) · Cv³)
注意分母用的是 n - 1 而不是 n，这是水文界沿用统计学无偏修正的习惯。虽然三阶矩在短系列下"无偏"也救不回来，但写法必须和同行保持一致，否则你算出的 Cv 会比别人小一截，适线时理论线整体偏低，最后交出去的设计值对不上审图意见。
用 Python 实现这段计算，核心代码只有十几行：
import numpy as np def p3_moments(x): # x: 实测序列，例如某站年最大洪峰流量 [m3/s] x = np.asarray(x, dtype=float) n = x.size x_mean = x.mean() # 一阶原点矩：均值 k = x / x_mean - 1.0 # 模比系数减 1 cv = np.sqrt(np.sum(k**2) / (n - 1)) cs = np.sum(k**3) / ((n - 1) * cv**3) return x_mean, cv, cs # 示例：某站 20 年最大洪峰流量 q = [1250, 980, 860, 1420, 760, 1130, 920, 1350, 700, 1080, 1550, 890, 1210, 950, 1370, 840, 1180, 1520, 1020, 880] mean_q, cv_q, cs_q = p3_moments(q) print(f"均值={mean_q:.1f}, Cv={cv_q:.3f}, Cs={cs_q:.3f}")这段代码的逻辑分三步：先算均值，再把每个值换算成模比系数并求离差平方和、三阶累积量，最后套公式输出。参数含义上，Cv 描述系列的相对离散程度，洪水系列一般在 0.2~1.5 之间，超过 1.5 说明流域调蓄能力很差、洪水暴涨暴落；Cs 描述偏斜程度，偏态系列算出来通常在 0.5~2 之间。但矩法算出的 Cs 只是粗糙初值，它受样本里最大的那几个值影响非常大，删掉一个 1550 的年份，Cs 可能从 1.2 直接掉到 0.6。所以初值归初值，最终要靠后面配线迭代去修正。
2.3 参数估计路径的选择与边界
除了矩法，还有极大似然法和概率加权矩法两条路。极大似然法对 P-III 型没有闭式解，要迭代求解，而且偏态系数大时似然面很平，迭代经常在几个相近参数之间来回震荡；概率加权矩法对长系列表现不错，但样本小于 30 时三阶概率加权矩的估计方差很大，算出来的 Cs 比矩法还飘。行业里几十年沉淀下来的做法就是"矩法初估 + 适线修正"，不是因为它统计上最优，而是因为稳定、可复现、审图专家看得懂。
提示：如果你用的是短系列（n < 20），矩法给的 Cs 可以直接扔进垃圾桶，别当真。后面避坑章会专门讲这个。
3. doingni2 的真相：配线迭代到底在迭代什么
3.1 经验频率：先把实测点放到概率轴上
配线前必须先给每个实测值一个"概率位置"，这就是经验频率。最常用的是数学期望公式：
P = m / (n + 1)
m 是排序后从大到小的序号，n 是样本长度。也有用中位秩公式 P = (m - 0.3) / (n + 0.4) 的，短系列时结果会整体往右偏一点。水文界习惯用 m/(n+1)，我一般也用它，关键是同一份报告里全用同一个公式。排序方向要小心：频率计算里第 1 位是最大值，对应最小概率 1/(n+1)，排反了整条曲线全部错位。
另外一个容易被忽略的细节：经验点据的 y 值直接用实测值，x 值是对应的频率 P，但在计算机实现里，P 不能直接作为横坐标送进绘图函数，后面第 4 章会专门讲坐标换算。这里先记住一点：经验频率只是把每个实测值得到了一个"坐标位置"，它本身不涉及任何分布假设。
3.2 doingni2 的核心逻辑：输入三参数，输出整条理论曲线
标题里的 doingni2，在压缩包里对应的就是那种"给定 p3 曲线的参数，算出一整条理论频率曲线"的子程序。常见的老程序结构是：传入均值、Cv、Cs 三个参数，内部对一组概率 P 逐个求解离均系数 φ，再换算成设计值 x = x̄·(1 + Cv·φ)。doingni2 这名字多半是作者随手起的简称，不重要，重要的是这个函数的迭代逻辑。
离均系数 φ 没有解析表达式，需要数值求解。P-III 型变量标准化后令 t = (x - x̄) / (x̄·Cv)，其概率密度是 α = 4/Cs² 的 Gamma 分布密度，累积概率通过数值积分得到，再用二分法反解出给定 P 对应的 t，也就是 φ。我第一次写这段时在积分上限上吃过亏：如果上限取无穷大，Cs 较大时尾部衰减很慢，精度差；实际取均值加 15 个标准差足够，精度和速度都兼顾。下面是完整可跑的计算函数：
from scipy.stats import gamma def phi_p3(cs, p, eps=1e-10): # 给定偏态系数Cs与累积概率P，返回P-III型离均系数phi alpha = 4.0 / cs**2 # 形状参数 loc = -2.0 / cs # 标准化后的位置参数 scale = cs / 2.0 # 标准化后的尺度参数 lo, hi = -20.0, 60.0 # 二分法初值范围，必须足够宽 for _ in range(100): mid = 0.5 * (lo + hi) cdf = gamma.cdf(mid, alpha, loc=loc, scale=scale) if cdf < p: lo = mid else: hi = mid return 0.5 * (lo + hi) def p3_curve(mean_q, cv, cs, probs): # 输入三参数与概率数组，返回理论频率曲线对应的设计值数组 phis = [phi_p3(cs, p) for p in probs] return mean_q * (1.0 + cv * np.array(phis))这段代码里，phi_p3 用 SciPy 的 gamma 分布做累积概率计算，省去自己写数值积分。二分法的 lo 和 hi 要覆盖分布支撑的关键区间：Cs 很大时分布右尾拖得很长，hi 取 60 是安全的；Cs 很小时分布接近正态，区间自动落到 -20 到 60 内。p3_curve 则是把离均系数线性换算回变量域，注意这里均值、Cv、Cs 三个参数是并列输入，缺一个都不行——这就是标题里"p3曲线的参数"指代的三个数。
3.3 完整配线流程：目标函数决定曲线贴哪里
有了理论曲线生成函数，配线就变成一个参数寻优问题。老程序里 doingni2 这类函数往往配合双重循环做网格搜索，现代实现也可以直接用优化器。但目标函数的选择比算法重要得多：最小二乘默认对每个经验点一视同仁，对洪水频率分析来说，小概率端 P=0.1%~1% 那两三个点才是设计值的关键，中部点再多也只是背景。常见做法是给尾部加权重，或者干脆只在 P≤20% 的范围内做最小二乘，效果都比全局等权好。
我一般让残差在概率轴上乘一个权重 w = 1/(P + 0.1)，P 越小权重越高，相当于自动化了老配线工"眼睛只盯左上角"的经验。完整拟合流程如下：
from scipy.optimize import minimize def fit_p3_weighted(emp_p, emp_x, mean_q, cv0, cs0): # emp_p: 经验频率数组, emp_x: 对应的实测值 emp_p = np.asarray(emp_p); emp_x = np.asarray(emp_x) weights = 1.0 / (emp_p + 0.1) # 小概率端加权 def obj(params): cv, cs = params theory = p3_curve(mean_q, cv, cs, emp_p) return np.sum(weights * (emp_x - theory) ** 2) res = minimize(obj, [cv0, cs0], bounds=[(0.05, 2.0), (0.1, 6.0)], method='L-BFGS-B') return res.x[0], res.x[1], res.fun # 用前面矩法算出的初值继续 emp_p = np.arange(1, 21) / 21.0 # 数学期望公式 P = m/(n+1) emp_x = np.array(q) # 排序前要先从大到小排列 order = np.argsort(-emp_x) emp_p = emp_p[order] # m=1对应最大值 emp_x = emp_x[order] cv_fit, cs_fit, rss = fit_p3_weighted(emp_p, emp_x, mean_q, cv_q, cs_q) print(f"适线结果: Cv={cv_fit:.3f}, Cs={cs_fit:.3f}, 残差={rss:.1f}")注意代码里有个容易漏的步骤：经验频率的序号 m 必须与从大到小排序后的实测值一一对应，所以先对 x 排序，再让 emp_p 跟着排序。weights 数组里 P 是小数，P=0.01 时权重约 9，P=0.5 时权重约 1.7，差别接近 5 倍，这保证了设计洪水最关心的小概率端点据在拟合中占主导地位。L-BFGS-B 的边界设置也很关键：Cv 下限 0.05 是防止退化到零方差，Cs 上限 6 对应 α=0.11，再大数值积分就不稳定了。
4. 概率格纸坐标换算与绘图：让理论线在图上"贴"住点据
4.1 概率格纸的本质是正态分位刻度
为什么皮尔逊III型曲线在普通直角坐标上是条弯弧，翻到概率格纸上却接近直线？因为概率格纸的横坐标根本不是 P 本身，而是 P 对应的标准正态分位 u = Φ⁻¹(P)。把曲线画在"x 轴是 u、y 轴是变量值"的坐标系里，接近正态的分布会被拉直；P-III 型是偏态分布，拉不直但曲率比普通坐标小得多，人眼更容易判断点据贴合度。
这个原理直接决定了绘图实现方式：先做一次 P 到 u 的换算，把数据点画在 u 坐标上，再把 x 轴刻度标签回写成 P。P=0% 和 P=100% 时 u 趋向正负无穷，程序里要裁剪范围，一般显示区间取 P = 0.01% 到 99.99% 足够满足工程需要，对应 u 大约是 -3.72 到 3.72。
4.2 概率格纸绘制函数与出图样例
用 matplotlib 实现概率格纸的标准做法如下，核心是刻度的换算和回标：
import matplotlib.pyplot as plt from scipy.stats import norm def setup_prob_axis(ax, p_min=0.0001, p_max=0.9999, n_ticks=13): # 把横轴改成 u = Phi^{-1}(P)，并回标概率刻度 ticks = np.linspace(norm.ppf(p_min), norm.ppf(p_max), n_ticks) labels = [f"{norm.cdf(u)*100:.2f}" for u in ticks] ax.set_xticks(ticks) ax.set_xticklabels(labels) ax.set_xlim(norm.ppf(p_min), norm.ppf(p_max)) probs = np.array([0.0001, 0.001, 0.01, 0.1, 1, 5, 20, 50, 80, 95, 99, 99.9, 99.99]) / 100.0 theory_x = p3_curve(mean_q, cv_fit, cs_fit, probs) u_theory = norm.ppf(probs) u_emp = norm.ppf(emp_p) # 经验点据也换算到u坐标 fig, ax = plt.subplots(figsize=(8, 5)) ax.plot(u_theory, theory_x, '-', lw=1.6, label='P-III理论频率曲线') ax.plot(u_emp, emp_x, 'o', ms=5, label='经验点据') setup_prob_axis(ax) ax.set_ylabel('流量 (m³/s)') ax.set_xlabel('频率 P (%)') ax.grid(True, ls='--', alpha=0.5) ax.legend() plt.show()这段代码的关键有两点：一是理论曲线和经验点据必须画在同一个 u 坐标里，如果只把理论线换到 u、点据还用 P，两者就会错开十万八千里；二是刻度标签回写用 norm.cdf，从 u 反算回 P，标签间距在正态分位上均匀、在概率值上不均匀，这是概率格纸的视觉特征。经验点据在图上通常画成空心圆或三角，理论线画实线，两条线同图对比，肉眼判断是否"贴住"。
4.3 Cv 与 Cs 对曲线形态的控制作用
参数对曲线形态的影响在概率格纸上一目了然。Cv 决定理论线的整体斜率：Cv 越大，线越陡，从 P=50% 到 P=1% 的抬升越明显；Cv 偏小，整条线会贴着点据的下包络走，设计值偏低，这在工程上是偏危险的。Cs 决定曲线的弯曲方向和程度：Cs 增大，小概率端抬高、曲线上翘得越厉害；Cs 接近 0 时，曲线退化成接近直线的高斯分布形态。所以在联调时有个实用经验：先定 Cv 让理论线穿越点据主体，再微调 Cs 让左上角的尾部点尽量落在线上，两个参数交替动，每次只调一个。老配线工说的"Cv 管斜率、Cs 管弯头"，在图上就是这个意思。
5. 皮尔逊III型拟合避坑指南：四条高频翻车现场
5.1 短系列下 Cs 是个黑匣子，矩法估值像抽签
现象：同一条河流相邻两站的 20 年系列，矩法算出来的 Cs 一个 1.5、一个 0.4，看起来完全不像同一个地区的洪水特性。换一段年份重新计算，某站的 Cs 又变成 0.8。
原因：三阶矩对样本里的极端值极度敏感。最大那一年流量占系列总量的比重很大，它的离均差三次方几乎决定了 Cs 的值；而短系列里这样的极端值只有一个，样本一变 Cs 就跟着跳。
解决：n < 30 的系列，别把矩法 Cs 当最终值。先用 Cv 和地区经验倍比定 Cs，常见做法是设 Cs = 2.5Cv 作为初值，再放进配线优化里只做小幅修正（比如限定在 2Cv ~ 3.5Cv 之间搜索），比让优化器自由放飞 Cs 可靠得多。我在 fit_p3_weighted 的边界里把 Cs 上限压到 6，实际工程中更常见的是直接固定倍比。
5.2 概率格纸横坐标被当成线性刻度
现象：画出来的图点据全部挤在右侧三分之一处，左侧一大片空白，理论线也堆在一起，完全看不出贴合情况。
原因：直接把 P 值当横坐标线性摆放了。P=0.1% 到 P=50% 在普通坐标轴上只占很小一段，概率格纸要求的正态分位变换没做，点据自然挤成一团。
解决：遵照第 4 章的 setup_prob_axis，所有数据点横坐标一律用 norm.ppf(P)，刻度标签回写概率值。检查是否做对有一个笨办法：P=50% 的 u 是 0，P=0.1% 的 u 是 -3.09，两者在轴上的距离应该比 P=50% 到 P=90%（u=1.28）远一截，这才是概率格纸应有的样子。
5.3 Cs/Cv 倍比失控导致理论线小概率端翘上天
现象：适线时为了让左上角最后一个点贴住线，把 Cs 调到 5 倍甚至 6 倍 Cv，结果理论曲线在 P=0.01% 处突然向上猛翘，设计值比点据外推趋势高出 40% 以上，明显不合理。
原因：Cs 过大时 α = 4/Cs² 很小，分布右尾变得极厚，分位数在小概率端爆炸式增长。这是 P-III 型分布本身的数学性质，不是绘图错误。
解决：给 Cs 设定硬约束，推荐上界 4Cv 或绝对值 3.0~4.0，哪个小用哪个。如果尾部点据实在贴不上，优先怀疑经验频率公式或者系列代表性不足，而不是无脑加 Cs。另外可以对比一下 Cs 取 3Cv 和 4Cv 时 P=0.1% 设计值的差异，若超过 10%，说明这个系列对偏态系数太敏感，设计成果要标注为"对尾部参数敏感"，提醒决策者这份成果的风险。
5.4 二分法求离均系数在极小概率处迭代失败
现象：P=0.01% 时 phi_p3 返回 NaN 或者一个明显偏大的数（比如 8.7），导致理论线在横轴最左端出现尖刺。把 P 改成 0.1% 就正常。
原因：二分法的初值区间不够宽。Cs=1.5 时 P=0.01% 对应的分位数大约在 5~6 之间，但如果 Cs 接近 4，尾部极厚，分位数可能超过 10；而我把 lo 固定从 -20 开始，理论上不会漏，真正的问题在于 Gamma 分布参数很大或很小时，scipy 的累积概率计算在极端分位处精度下降，迭代 100 次还不收敛。
解决：给二分法加收敛判断，循环里同时检查区间宽度，宽度小于 1e-8 就提前退出；再对 P 加一个最小值保护，P < 1e-5 时直接按 P = 1e-5 计算，避免概率值太小导致 cdf 计算下溢。改完的 phi_p3 在 P=0.01% 时依然能稳定返回结果。
5.5 经验频率公式混用，不同报告参数对不上
现象：核对上一版报告的 P-III 参数，发现同样的实测系列，Cv 一致但 Cs 差 0.3，设计洪水差 8%。追溯发现上一版用的是 P = (m - 0.3)/(n + 0.4)，新算的用了 m/(n + 1)。
原因：两个经验频率公式在短系列时的概率位置差异很大，尤其首项：n=20 时 m/(n+1) 首项概率是 4.76%，(m-0.3)/(n+0.4) 首项是 3.43%，点据位置整体移动，理论线为了贴合点据只能跟着移动，参数自然偏差。
解决：固定一套公式走到底，并在成果表里注明"经验频率采用数学期望公式 P=m/(n+1)"。做对比分析时必须用同一公式重算，不能拿两份不同公式的结果比较。项目交接时把公式写进计算书，别只留一张参数表。
6. 一个收尾技巧：用残差盯住适线是否真的到位
配线完成后，别急着把 Cv、Cs 填进报告，先做一个残差检查。把每个经验点据的实测值 xᵢ 与同概率下理论值 x̂ᵢ 的差求出来，除以理论值得到相对残差，然后按概率从大到小排列，观察残差序列是否随机分布在零线两侧。如果小概率端残差全是正的、大概率端全是负的，说明 Cv 系统性偏大，整体斜率需要下调；如果两端残差正、中间负，呈现明显的"碗形"，说明 Cs 偏大导致曲线弯头太急，需要减小倍比。这一步五秒就能看完，但能省掉返工交底的麻烦。
def check_fit(emp_p, emp_x, mean_q, cv, cs): theory = p3_curve(mean_q, cv, cs, emp_p) rel_res = (emp_x - theory) / theory * 100.0 order = np.argsort(emp_p) return emp_p[order], rel_res[order] p_sorted, rel = check_fit(emp_p, emp_x, mean_q, cv_fit, cs_fit) for p_val, r in zip(p_sorted, rel): print(f"P={p_val*100:6.2f}% 相对残差={r:+7.2f}%")残差打印出来后，P=0.1% 到 10% 这几行尤其要看。相对残差控制在 ±10% 以内算合格，超过 20% 的点一定要解释原因，要么是系列里有历史特大洪水需要做历史洪水处理，要么是配线权重没给足尾部。我自己在这上面吃过亏：有一回拟合结果整体漂亮，唯独 P=1% 的点偏了 18%，没深究就出了报告，评审专家一眼看出那个点对应的是流域内一次著名的大洪水年份，数据代表性有问题，当场要求返工。从那以后，残差检查成了我出具任何频率分析成果前雷打不动的习惯——配线不是画完图就结束，最后一道关是看残差，希望帮到你。
本文还有配套的精品资源，点击获取
