简介：这份资源聚焦皮尔逊III型曲线的拟合与绘图，面向水文学、统计学及水资源管理领域的学习者与工程技术人员，帮助解决流量、降雨等连续随机变量的概率分布建模与频率分析问题。压缩包共3个文件，均为m脚本文件，整体约3KB，分别对应曲线拟合、参数计算与绘图等核心环节，便于在MATLAB环境中直接调用与二次修改。内容围绕形状参数κ、尺度参数θ和位置参数μ的估计展开，涵盖最大似然估计思路、概率密度函数绘制，以及通过插值获取不同累积概率下流量数据的方法，可用于风险评估与防洪规划中的极端流量推算。目前已有1303人学习下载，适合希望掌握P3曲线参数求解、图形展示与频率插值完整流程的读者参考实践。

1. 皮尔逊III型曲线拟合与绘图：从参数估计到设计频率插值的完整链路

做水文频率分析的人，迟早会撞上“适线法”这三个字。手里攥着三十年甚至五十年的年最大流量序列，要在几张概率格纸上把经验点据排好，再徒手调一条皮尔逊III型曲线去逼近它们——这个过程在过去靠的是透明尺和肉眼，现在则完全可以交给代码。我拿到的这份皮尔逊III型曲线拟合与绘图.rar就是干这件事的：它把参数估计、曲线绘制、频率插值三件事拆成了三个 MATLAB 脚本，pl_p3.m、p3_cz.m、p3_nnd.m，分别对应绘图、参数初值计算和拟合迭代。如果你正在做防洪评价、水库调洪或者水资源论证，需要从实测流量序列里推求 P=1%、P=2% 这类设计频率对应的洪峰流量，这套东西能直接省掉你重新推导公式的时间。它适合两类人：一类是水文水资源专业的学生，想搞明白适线法背后的数学到底怎么落地；另一类是一线工程师，手头有数据、有 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是整个包里最核心的拟合脚本。它的思路是：给定一组参数，计算理论频率曲线与经验点据之间的离差平方和，然后用优化算法最小化这个离差。目标函数通常有两种形式：

1. 纵向离差：对每个经验频率点，计算理论分位数与实测值之差。
2. 横向离差：对每个实测值，计算理论累积概率与经验概率之差。

水文上更常用纵向离差，因为流量值的量纲更直观。下面是我常用的目标函数写法：

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并手动提供梯度。不过对水文序列来说，样本量通常不超过 100，fminsearch完全够用。

3.2 经验频率的计算与绘点位置

拟合之前必须先算经验频率。国内水文规范常用的公式是数学期望公式：

$$P_m = \frac{m}{n+1} \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负责绘图。它的核心是生成一条从 P=0.01% 到 P=99.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 插值获取设计频率对应的流量

有了拟合参数，插值就是一句话的事。比如你要 P=1% 的设计洪峰流量：

P_design = [0.01, 0.02, 0.05, 0.10]; % 设计频率 x_design = mu_fit + theta_fit * gaminv(P_design, kappa_fit, 1); fprintf('P=1%% 设计流量：%.2f m^3/s\n', x_design(1)); fprintf('P=2%% 设计流量：%.2f m^3/s\n', x_design(2)); fprintf('P=5%% 设计流量：%.2f m^3/s\n', x_design(3)); fprintf('P=10%% 设计流量：%.2f m^3/s\n', x_design(4));

gaminv的第二个参数是形状参数 κ，第三个参数是尺度参数，这里固定为 1，因为尺度已经通过 θ 体现了。注意P_design用的是小数形式，不是百分数。如果你输入 1 而不是 0.01，gaminv会返回无穷大，因为累积概率不能超过 1。

注意：插值得到的设计流量是理论值，实际工程中还要考虑样本代表性、历史特大洪水处理等问题。这套代码只负责数学拟合，不负责水文合理性判断。

4.3 拟合优度检验的补充手段

拟合完之后，怎么判断拟合得好不好？除了看图形，还可以算几个统计量：

这些指标可以在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 小的一端）的拟合精度。如果目标函数对所有点一视同仁，小频率端的偏差会被大频率端的点“平均”掉。

解决：给目标函数加权。常用的做法是给 P<10% 的点赋更大的权重，比如权重取 $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 曲线、算出 P=1% 的设计流量，这只是点估计。评审专家经常会问：这个设计流量的置信区间是多少？这时候就需要 Bootstrap 方法。思路很简单：从原始序列里有放回地重采样 1000 次，每次重新拟合 P3 分布，算出 P=1% 的设计流量，然后取 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')); % 计算 P=1% 设计流量 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('P=1%% 设计流量 95%% 置信区间：[%.2f, %.2f] m^3/s\n', ci_lower, ci_upper);

这段代码跑 1000 次大概需要几分钟，取决于你的机器性能。如果嫌慢，可以降到 500 次，精度也够用。Bootstrap 的好处是不需要假设分布形式，完全由数据驱动。缺点是计算量大，而且如果原始序列太短（比如少于 20 年），重采样出来的样本可能仍然缺乏代表性。

我一般会在最终报告里同时给出点估计和区间估计。如果区间宽度超过点估计的 30%，说明数据的不确定性太大，设计流量要谨慎使用。这时候我会建议补充历史洪水调查或者用区域综合法来佐证。

从那以后我每次做频率分析，都会强制走一遍 Bootstrap，哪怕评审没要求。因为水文计算里最怕的就是“看起来精确”的假象，一个点估计值背后如果没有不确定性量化，心里总是不踏实。希望这套脚本和上面的补充代码能帮到你，少走一些我当年走过的弯路。

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