简介这份资源聚焦皮尔逊III型曲线的拟合与绘图面向水文学、统计学及水资源管理领域的学习者与工程技术人员帮助解决流量、降雨等连续随机变量的概率分布建模与频率分析问题。压缩包共3个文件均为m脚本文件整体约3KB分别对应曲线拟合、参数计算与绘图等核心环节便于在MATLAB环境中直接调用与二次修改。内容围绕形状参数κ、尺度参数θ和位置参数μ的估计展开涵盖最大似然估计思路、概率密度函数绘制以及通过插值获取不同累积概率下流量数据的方法可用于风险评估与防洪规划中的极端流量推算。目前已有1303人学习下载适合希望掌握P3曲线参数求解、图形展示与频率插值完整流程的读者参考实践。1. 皮尔逊III型曲线拟合与绘图从参数估计到设计频率插值的完整链路做水文频率分析的人迟早会撞上“适线法”这三个字。手里攥着三十年甚至五十年的年最大流量序列要在几张概率格纸上把经验点据排好再徒手调一条皮尔逊III型曲线去逼近它们——这个过程在过去靠的是透明尺和肉眼现在则完全可以交给代码。我拿到的这份皮尔逊III型曲线拟合与绘图.rar就是干这件事的它把参数估计、曲线绘制、频率插值三件事拆成了三个 MATLAB 脚本pl_p3.m、p3_cz.m、p3_nnd.m分别对应绘图、参数初值计算和拟合迭代。如果你正在做防洪评价、水库调洪或者水资源论证需要从实测流量序列里推求 P1%、P2% 这类设计频率对应的洪峰流量这套东西能直接省掉你重新推导公式的时间。它适合两类人一类是水文水资源专业的学生想搞明白适线法背后的数学到底怎么落地另一类是一线工程师手头有数据、有 MATLAB但不想从零写优化器。2. 皮尔逊III型曲线的参数体系κ、θ、μ 到底怎么定2.1 三参数与水文统计量的对应关系皮尔逊III型分布的概率密度函数写作$$f(x) \frac{1}{\theta \Gamma(\kappa)} \left(\frac{x-\mu}{\theta}\right)^{\kappa-1} e^{-\frac{x-\mu}{\theta}}$$其中 κ 是形状参数θ 是尺度参数μ 是位置参数。但水文上很少直接用这三个参数更常见的是用均值 $\bar{x}$、变差系数 $C_v$、偏态系数 $C_s$ 来刻画。两者之间的换算关系是$\mu \bar{x}(1 - 2C_v / C_s)$$\theta \bar{x} \cdot C_v \cdot C_s / 2$$\kappa 4 / C_s^2$这个换算关系是整套代码的基石。p3_cz.m里做的第一件事就是根据样本的 $\bar{x}$、$C_v$、$C_s$ 算出 κ、θ、μ 的初值。为什么要用初值因为后续的拟合是一个非线性优化问题直接扔给fminsearch或者lsqnonlin很容易陷进局部极小值。我一般会先用矩法估一组初值再让优化器去微调。注意当 $C_s 2C_v$ 时μ 会变成负数这意味着分布的下界在负半轴。对流量数据来说这显然不合理但数学上没问题拟合时不用刻意回避。2.2 矩法初值估计的代码实现p3_cz.m的核心逻辑就是矩法估计。下面是我根据脚本逻辑重写的等价代码方便你理解每一步在算什么function [mu, theta, kappa] p3_cz(x) % P3_CZ 根据样本序列用矩法估计皮尔逊III型参数初值 % 输入 x: 年最大流量序列列向量 % 输出 mu, theta, kappa: 位置、尺度、形状参数 x x(:); % 强制列向量 n length(x); xbar mean(x); % 样本均值 sigma std(x); % 样本标准差 Cv sigma / xbar; % 变差系数 % 偏态系数用三阶中心矩计算 m3 sum((x - xbar).^3) / n; Cs m3 / sigma^3; % 参数换算 mu xbar * (1 - 2*Cv/Cs); theta xbar * Cv * Cs / 2; kappa 4 / Cs^2; end这段代码里最关键的是Cs的计算。很多初学者会直接用skewness(x)但 MATLAB 自带的skewness默认做了偏差校正算出来的值和手算三阶矩有细微差别。在水文频率计算中通常要求用无偏估计或者矩法估计具体用哪个取决于你的规范要求。我一般会两种都算一遍对比一下差异如果Cs差别超过 0.05就要检查数据里有没有异常值。参数说明xbar是均值sigma是标准差Cv是变差系数Cs是偏态系数。mu、theta、kappa就是最终要传给拟合函数的初值。注意Cs不能为零否则kappa会无穷大代码会报错。实际数据里Cs一般不会恰好为零但如果你的序列特别对称Cs可能接近 0.01这时候kappa会到 40000数值上仍然可算但拟合会非常慢。3. 适线法拟合从初值到最优参数的迭代过程3.1 目标函数与优化器的选择p3_nnd.m是整个包里最核心的拟合脚本。它的思路是给定一组参数计算理论频率曲线与经验点据之间的离差平方和然后用优化算法最小化这个离差。目标函数通常有两种形式纵向离差对每个经验频率点计算理论分位数与实测值之差。横向离差对每个实测值计算理论累积概率与经验概率之差。水文上更常用纵向离差因为流量值的量纲更直观。下面是我常用的目标函数写法function sse p3_obj(params, x, p_emp) % P3_OBJ 计算理论曲线与经验点据的纵向离差平方和 % params: [mu, theta, kappa] % x: 实测流量序列 % p_emp: 经验频率按大小排位后计算 mu params(1); theta params(2); kappa params(3); % 参数合法性检查 if theta 0 || kappa 0 sse 1e10; % 返回一个大数迫使优化器远离非法区域 return; end % 计算理论分位数 x_theo mu theta * gaminv(p_emp, kappa, 1); % 纵向离差平方和 sse sum((x - x_theo).^2); end这里用到了gaminv它是 Gamma 分布的逆累积分布函数。皮尔逊III型分布和 Gamma 分布的关系是如果 $Y \sim \text{Gamma}(\kappa, 1)$那么 $X \mu \theta Y$ 就服从皮尔逊III型分布。所以gaminv(p, kappa, 1)算出来的是标准 Gamma 分位数乘以 θ 再加 μ 就得到 P3 分位数。优化器我一般用fminsearch因为它不需要计算梯度对初值不敏感。但fminsearch的收敛速度慢如果数据量大可以换成fminunc并手动提供梯度。不过对水文序列来说样本量通常不超过 100fminsearch完全够用。3.2 经验频率的计算与绘点位置拟合之前必须先算经验频率。国内水文规范常用的公式是数学期望公式$$P_m \frac{m}{n1} \times 100%$$其中 m 是降序排列的序号n 是样本容量。pl_p3.m里用的就是这个公式。但有些文献会用中值公式或者海森公式区别在于对极端值的处理。我建议你固定用数学期望公式因为它在国内规范里最通用评审专家不会挑毛病。% 经验频率计算 x_sorted sort(x, descend); n length(x_sorted); m (1:n); p_emp m / (n 1); % 数学期望公式这段代码把实测序列降序排列然后给每个点算一个经验频率。注意p_emp是累积概率从大到小对应流量从大到小。绘图时横轴是频率通常用对数刻度或者概率格纸刻度纵轴是流量。3.3 完整拟合流程的代码串联把上面的片段串起来就是一个完整的拟合脚本% 读取数据假设数据存为 data.txt一列流量值 x load(data.txt); % 步骤1矩法估计初值 [mu0, theta0, kappa0] p3_cz(x); % 步骤2计算经验频率 x_sorted sort(x, descend); n length(x_sorted); m (1:n); p_emp m / (n 1); % 步骤3优化拟合 params0 [mu0, theta0, kappa0]; options optimset(Display, iter, MaxIter, 2000); params_opt fminsearch((p) p3_obj(p, x_sorted, p_emp), params0, options); % 步骤4输出最优参数 mu_fit params_opt(1); theta_fit params_opt(2); kappa_fit params_opt(3); fprintf(拟合结果mu%.4f, theta%.4f, kappa%.4f\n, mu_fit, theta_fit, kappa_fit);这段代码跑完你就得到了最优的 μ、θ、κ。Display设为iter可以看到每一步的迭代信息方便判断是否收敛。如果迭代超过 2000 次还没收敛多半是初值太离谱或者数据里有异常值把目标函数带偏了。提示fminsearch对参数的量纲很敏感。如果 μ 是几百、θ 是几十、κ 是个位数优化器会在不同方向上以不同步长搜索效率很低。我一般会先把数据标准化拟合完再还原。不过对这套脚本来说数据量不大不标准化也能跑。4. 绘图与频率插值从参数到设计流量4.1 概率格纸上的理论曲线绘制pl_p3.m负责绘图。它的核心是生成一条从 P0.01% 到 P99.99% 的理论曲线然后叠加经验点据。绘图的关键在于横轴的刻度变换——概率格纸不是线性的而是把正态分布的分位数作为刻度。MATLAB 里可以用norminv来做这个变换% 生成理论曲线 p_theo logspace(-4, -0.0001, 500); % 从0.01%到99.99% x_theo mu_fit theta_fit * gaminv(p_theo, kappa_fit, 1); % 绘图 figure; semilogx(p_theo*100, x_theo, b-, LineWidth, 1.5); hold on; semilogx(p_emp*100, x_sorted, ro, MarkerFaceColor, r); xlabel(频率 P (%)); ylabel(流量 (m^3/s)); grid on;这里用semilogx是因为概率格纸的横轴在对数尺度下近似线性。如果你想要更精确的概率格纸效果可以用norminv把概率转成正态分位数然后在线性坐标下绘图% 正态分位数变换 z_theo norminv(1 - p_theo); z_emp norminv(1 - p_emp); figure; plot(z_theo, x_theo, b-, LineWidth, 1.5); hold on; plot(z_emp, x_sorted, ro, MarkerFaceColor, r); xlabel(正态分位数); ylabel(流量 (m^3/s)); grid on;两种画法效果差不多第一种更直观第二种更符合传统概率格纸的习惯。我一般用第一种因为出图快评审也能看懂。4.2 插值获取设计频率对应的流量有了拟合参数插值就是一句话的事。比如你要 P1% 的设计洪峰流量P_design [0.01, 0.02, 0.05, 0.10]; % 设计频率 x_design mu_fit theta_fit * gaminv(P_design, kappa_fit, 1); fprintf(P1%% 设计流量%.2f m^3/s\n, x_design(1)); fprintf(P2%% 设计流量%.2f m^3/s\n, x_design(2)); fprintf(P5%% 设计流量%.2f m^3/s\n, x_design(3)); fprintf(P10%% 设计流量%.2f m^3/s\n, x_design(4));gaminv的第二个参数是形状参数 κ第三个参数是尺度参数这里固定为 1因为尺度已经通过 θ 体现了。注意P_design用的是小数形式不是百分数。如果你输入 1 而不是 0.01gaminv会返回无穷大因为累积概率不能超过 1。注意插值得到的设计流量是理论值实际工程中还要考虑样本代表性、历史特大洪水处理等问题。这套代码只负责数学拟合不负责水文合理性判断。4.3 拟合优度检验的补充手段拟合完之后怎么判断拟合得好不好除了看图形还可以算几个统计量检验指标计算公式判断标准均方根误差 RMSE$\sqrt{\frac{1}{n}\sum(x_i - \hat{x}_i)^2}$越小越好一般要求小于均值的 10%确定性系数 R²$1 - \frac{\sum(x_i - \hat{x}_i)^2}{\sum(x_i - \bar{x})^2}$越接近 1 越好一般要求大于 0.9相对偏差$\frac{1}{n}\sum\frac{x_i - \hat{x}_i这些指标可以在p3_nnd.m里加几行代码算出来。我一般会把它们打印在图上方便一眼判断拟合质量。5. 避坑与排查那些年我踩过的皮尔逊III型曲线坑5.1 现象拟合结果对初值极度敏感换一组初值结果就变了原因目标函数存在多个局部极小值fminsearch从不同初值出发会收敛到不同解。尤其是当Cs较大时目标函数曲面非常崎岖。解决不要只跑一次拟合。我一般会用三组初值分别拟合矩法初值、Cs加 0.1 的初值、Cs减 0.1 的初值然后取目标函数最小的那组。如果三组结果差异很大说明数据本身不适合用 P3 分布要考虑换线型。5.2 现象gaminv返回 NaN 或 Inf原因kappa或theta在优化过程中变成了负数或零。gaminv要求形状参数和尺度参数都为正。解决在目标函数里加合法性检查非法参数直接返回一个大数。上面p3_obj里的if theta 0 || kappa 0就是干这个的。另外fminsearch没有边界约束如果初值离合法区域太远优化器可能一开始就跑到非法区域。这时候可以用fminbnd或者手动加惩罚项。5.3 现象经验点据在两端偏离理论曲线很远原因两端点据的权重在最小二乘中是一样的但水文上更关心大频率端P 小的一端的拟合精度。如果目标函数对所有点一视同仁小频率端的偏差会被大频率端的点“平均”掉。解决给目标函数加权。常用的做法是给 P10% 的点赋更大的权重比如权重取 $1/P$ 或者 $1/\sqrt{P}$。修改p3_obj里的sse计算w 1 ./ p_emp; % 频率越小权重越大 sse sum(w .* (x - x_theo).^2);这样拟合出来的曲线在小频率端会更贴近经验点据更符合工程需求。5.4 现象数据里有零值或负值gaminv报错原因皮尔逊III型分布的定义域是 $x \mu$如果 μ 是负数定义域可以覆盖零和负值。但如果数据里真的有零或负值而拟合出的 μ 是正数就会出现理论分位数小于 μ 的情况gaminv会返回 NaN。解决先检查数据。流量数据一般不会有负值但零值可能出现比如断流。如果零值很少可以直接剔除如果零值很多说明 P3 分布不适合要考虑用对数正态分布或者 Weibull 分布。5.5 现象绘图时横轴刻度不对概率格纸看起来是歪的原因semilogx只是对数刻度不是真正的概率格纸刻度。概率格纸的刻度是正态分位数不是对数。解决用norminv做变换然后在线性坐标下绘图。或者用 MATLAB 的probplot函数它自带概率格纸刻度。不过probplot默认是正态分布要改成 P3 分布需要自定义。我一般直接用norminv手动变换灵活且可控。6. 进阶技巧用 Bootstrap 量化设计流量的不确定性拟合出一条 P3 曲线、算出 P1% 的设计流量这只是点估计。评审专家经常会问这个设计流量的置信区间是多少这时候就需要 Bootstrap 方法。思路很简单从原始序列里有放回地重采样 1000 次每次重新拟合 P3 分布算出 P1% 的设计流量然后取 2.5% 和 97.5% 分位数作为 95% 置信区间。n_boot 1000; x_design_boot zeros(n_boot, 1); for i 1:n_boot % 有放回重采样 x_boot randsample(x, length(x), true); % 矩法初值 [mu0, theta0, kappa0] p3_cz(x_boot); % 经验频率 x_boot_sorted sort(x_boot, descend); m_boot (1:length(x_boot)); p_boot m_boot / (length(x_boot) 1); % 拟合 params_boot fminsearch((p) p3_obj(p, x_boot_sorted, p_boot), ... [mu0, theta0, kappa0], ... optimset(Display, off)); % 计算 P1% 设计流量 x_design_boot(i) params_boot(1) params_boot(2) * ... gaminv(0.01, params_boot(3), 1); end % 置信区间 ci_lower prctile(x_design_boot, 2.5); ci_upper prctile(x_design_boot, 97.5); fprintf(P1%% 设计流量 95%% 置信区间[%.2f, %.2f] m^3/s\n, ci_lower, ci_upper);这段代码跑 1000 次大概需要几分钟取决于你的机器性能。如果嫌慢可以降到 500 次精度也够用。Bootstrap 的好处是不需要假设分布形式完全由数据驱动。缺点是计算量大而且如果原始序列太短比如少于 20 年重采样出来的样本可能仍然缺乏代表性。我一般会在最终报告里同时给出点估计和区间估计。如果区间宽度超过点估计的 30%说明数据的不确定性太大设计流量要谨慎使用。这时候我会建议补充历史洪水调查或者用区域综合法来佐证。从那以后我每次做频率分析都会强制走一遍 Bootstrap哪怕评审没要求。因为水文计算里最怕的就是“看起来精确”的假象一个点估计值背后如果没有不确定性量化心里总是不踏实。希望这套脚本和上面的补充代码能帮到你少走一些我当年走过的弯路。本文还有配套的精品资源点击获取