基于Copula的气象-农业干旱联合概率与重现期计算（MATLAB实操）

干旱研究里有个绕不开的坎：气象干旱和农业干旱明明是两个不同维度的东西，一个看降水亏缺，一个看土壤水分和作物响应，但它们之间又存在明显的传导和滞后关系。单独算一个干旱指数的频率分布，只能回答"今年降水偏少到什么程度"，回答不了"这种气象条件下，农业干旱发生的可能性有多大"。Copula函数就是干这件事的——它能把两个边缘分布不同的变量"粘"在一起，构造出联合分布，进而算出联合概率、条件概率、重现期这些真正能支撑风险评估的指标。这篇内容面向的是已经会用MATLAB做基本数据处理、但对Copula还停留在"听说过"阶段的读者，从边缘分布拟合一路讲到联合概率计算和重现期绘制，把中间那些论文里通常一笔带过的细节补全。

1. 先搞清楚Copula到底在解决什么问题

1.1 为什么不能直接对两个干旱指数做联合分布拟合

假设你手上有某站点30年的SPI（标准化降水指数，代表气象干旱）和SSI（标准化土壤湿度指数，代表农业干旱）序列。最直觉的做法是把这两个变量当成二维正态分布来处理，直接估计均值向量和协方差矩阵。问题在于，SPI和SSI的边缘分布未必是正态的，而且它们之间的依赖结构也未必是线性的。用二维正态去拟合，等于同时假设了两件事：边缘是正态的，依赖结构是高斯的。这两个假设在干旱数据上经常不成立。

Copula的思路是把这个问题拆成两步：第一步，分别对SPI和SSI拟合各自的边缘分布，把它们转换成[0,1]区间上的均匀分布；第二步，用一个Copula函数来描述这两个均匀变量之间的依赖结构。这样边缘分布和依赖结构就解耦了，你可以给SPI配一个Gamma分布，给SSI配一个正态分布，然后用Clayton Copula去刻画它们之间的下尾依赖——干旱这种极端事件恰恰是下尾相关的，一个变量取极小值时另一个也倾向取极小值。

1.2 干旱研究里最常用的三类Copula

实际做干旱联合概率，用得最多的就三类：Archimedean族（Clayton、Gumbel、Frank）、椭圆族（Gaussian、t）和极值族。Archimedean族因为形式简单、参数少、能刻画非对称依赖，在干旱领域占了绝大多数。

干旱是"越干越相关"的典型下尾事件，所以Clayton通常是首选。但不要凭直觉定，后面会讲怎么用拟合优度检验来选。

1.3 联合概率、条件概率、重现期分别对应什么实际含义

这三个指标是干旱风险评估的核心输出，含义完全不同：

• 联合概率P(SPI ≤ a 且 SSI ≤ b)：气象干旱和农业干旱同时达到某个等级的概率。用于评估"复合干旱"风险。
• 条件概率P(SSI ≤ b | SPI ≤ a)：已知气象干旱达到某等级时，农业干旱也达到某等级的概率。这是预警里最有用的指标，因为它回答了"气象干旱已经发生了，农业干旱会不会跟上"。
• 联合重现期：在给定联合概率下，事件平均多少年发生一次。注意重现期有"且"和"或"两种定义，算出来的数值差别很大，论文里必须写清楚用的是哪种。

2. 数据准备与边缘分布拟合的实操细节

2.1 干旱指数的选择与时间尺度匹配

气象干旱用SPI，农业干旱用SSI或SMI，这是主流做法。但有个容易被忽略的点：时间尺度必须匹配。SPI-3（3个月尺度）对应的是短期降水亏缺，而土壤湿度对降水的响应通常滞后1到2个月，所以SSI也取3个月尺度比较合理。如果你用SPI-1去对SSI-6，物理意义就对不上了，算出来的依赖结构会很弱，Copula参数估计出来接近独立。

我一般建议先做互相关分析，看SPI和SSI在哪个滞后阶数上相关性最强，再决定时间尺度的搭配。MATLAB里用xcorr就能快速看：

% 假设spi和ssi是等长的列向量 [c, lags] = xcorr(spi, ssi, 12, 'coeff'); [max_c, idx] = max(abs(c)); best_lag = lags(idx); fprintf('最大相关滞后阶数: %d, 相关系数: %.3f\n', best_lag, c(idx));

如果best_lag是正的，说明SSI滞后于SPI，符合物理预期。如果best_lag是0附近且相关性很高，说明两个指数几乎同步，那时间尺度可以直接对齐。

2.2 边缘分布拟合：参数估计与拟合优度检验

边缘分布的选择直接决定Copula的输入。常用的候选分布有：Gamma、Weibull、Log-normal、Normal、Generalized Extreme Value。SPI本身在计算时已经做了正态化处理，所以SPI序列理论上接近标准正态，但实际样本里还是会有偏态，建议还是做一次拟合检验。

MATLAB里拟合边缘分布用fitdist，检验用kstest或adtest：

% 对SPI拟合正态分布 pd_spi = fitdist(spi, 'Normal'); [h, p] = kstest((spi - pd_spi.mu) / pd_spi.sigma); % h=0表示不能拒绝原假设，拟合可接受 % 对SSI拟合Gamma分布 pd_ssi = fitdist(ssi, 'Gamma'); [h2, p2] = kstest(ssi, 'CDF', pd_ssi);

注意：kstest对参数估计后的分布检验偏保守，样本量小于50时p值不太可靠。干旱研究里30到60年的序列很常见，建议同时看AD检验和PPCC（概率点相关系数），三个指标综合判断。

拟合完之后，把原始序列通过概率积分变换转成均匀分布：

u = cdf(pd_spi, spi); v = cdf(pd_ssi, ssi); % u和v现在应该在[0,1]上近似均匀

这里有个坑：如果某个观测值的CDF算出来正好是0或1，Copula的对数似然函数会变成无穷大。处理方法是在边界处做微小截断，比如把小于1e-6的值设为1e-6，大于1-1e-6的设为1-1e-6。这个操作对结果影响极小，但不做的话程序直接报错。

2.3 用Kendall秩相关系数初步判断依赖强度

在选Copula之前，先算一下Kendall's tau和Spearman's rho，这两个秩相关系数不依赖边缘分布，能直接反映依赖结构的强弱和方向。

tau = corr(spi, ssi, 'Type', 'Kendall'); rho = corr(spi, ssi, 'Type', 'Spearman'); fprintf('Kendall tau: %.3f, Spearman rho: %.3f\n', tau, rho);

对于Clayton Copula，参数theta和tau的关系是 theta = 2*tau/(1-tau)。如果tau只有0.1左右，算出来的theta很小，Copula接近独立，这时候做联合概率意义不大。一般tau在0.3以上，联合分析才有实际价值。如果tau偏低，先回头检查时间尺度是否匹配、数据是否有趋势需要去趋势。

3. Copula参数估计与拟合优度检验

3.1 用最大似然估计Copula参数

MATLAB没有内置的Copula工具箱（Statistics and Machine Learning Toolbox里有copulafit和copulacdf，但只支持椭圆族），Archimedean族需要自己写。核心是构造Copula的密度函数，然后对参数做最大似然估计。

以Clayton为例，其密度函数为：

c(u,v;θ) = (1+θ) * (u*v)^(-θ-1) * (u^(-θ) + v^(-θ) - 1)^(-2-1/θ)

对数似然函数就是对所有样本点的c取对数再求和。用fminbnd做一维搜索：

function negLL = clayton_negLL(theta, u, v) if theta <= 0 negLL = 1e10; return; end n = length(u); term1 = n * log(1 + theta); term2 = -(theta + 1) * sum(log(u) + log(v)); term3 = -(2 + 1/theta) * sum(log(u.^(-theta) + v.^(-theta) - 1)); negLL = -(term1 + term2 + term3); end % 估计 tau = corr(spi, ssi, 'Type', 'Kendall'); theta0 = 2*tau/(1-tau); % 用矩估计作为初值 theta_hat = fminbnd(@(t) clayton_negLL(t, u, v), 0.01, 20);

用矩估计值作为初值是个好习惯，因为极大似然在theta接近0时容易收敛到边界。fminbnd的搜索区间上界设20足够了，Clayton参数超过20意味着tau超过0.95，干旱数据里几乎不可能出现。

3.2 拟合优度检验：AIC、BIC与Rosblatt变换

估计完参数不能直接用，得检验拟合好不好。最常用的是AIC和BIC：

logL = -clayton_negLL(theta_hat, u, v); AIC = -2*logL + 2*1; % 1个参数 BIC = -2*logL + log(n)*1;

把Clayton、Gumbel、Frank都估一遍，选AIC最小的。但AIC只能比较相对好坏，不能告诉你"这个Copula拟合得可以接受"。更严格的检验是Rosblatt变换：把(u,v)通过Copula转换成新的变量，如果Copula拟合正确，转换后的变量应该独立且均匀。然后用Cramer-von Mises统计量检验。

实操中我一般两步走：先看AIC选最优，再用Rosblatt变换的p值确认拟合可接受。如果所有候选Copula的p值都小于0.05，说明数据里有Copula刻画不了的复杂依赖结构，这时候要么考虑混合Copula，要么回头检查数据质量。

3.3 尾部依赖系数：Clayton为什么适合干旱

尾部依赖系数衡量的是一个变量取极端值时另一个变量也取极端值的倾向。Clayton的下尾依赖系数是 2^(-1/θ)，上尾依赖系数为0。这意味着Clayton只能刻画"同时偏枯"，不能刻画"同时偏丰"。

干旱研究关心的是偏枯端，所以Clayton天然合适。但如果你研究的是干旱和洪涝的联合风险（比如同一个流域既有干旱又有暴雨），那就需要能刻画双尾依赖的t-Copula。选Copula之前先想清楚你关心的是哪个尾部。

4. 联合概率、条件概率与重现期的计算

4.1 联合概率的两种形式：且与或

这是最容易出错的地方。联合概率有两种定义：

• P(U≤u 且 V≤v)= C(u,v)，两个变量同时小于等于某阈值
• P(U≤u 或 V≤v)= u + v - C(u,v)，至少有一个小于等于某阈值

在干旱风险评估里，"且"对应的是复合干旱事件——气象干旱和农业干旱同时发生；"或"对应的是任一类干旱发生就算。论文里如果不写清楚，审稿人一定会问。

% 计算联合概率 u0 = 0.2; % SPI对应的分位数，比如SPI=-0.84对应0.2 v0 = 0.3; % SSI对应的分位数 P_and = copulacdf('Clayton', [u0, v0], theta_hat); P_or = u0 + v0 - P_and;

4.2 条件概率：预警里最有用的指标

条件概率 P(V≤v | U≤u) = C(u,v)/u。这个指标回答的是"气象干旱已经达到某等级时，农业干旱达到某等级的概率"。

P_cond = copulacdf('Clayton', [u0, v0], theta_hat) / u0;

举个例子，如果SPI≤-1（u0≈0.159）时，SSI≤-1（v0≈0.159）的条件概率是0.45，意味着气象干旱达到中度时，有45%的概率农业干旱也达到中度。这个数字比无条件概率（约0.159）高出一大截，说明两者确实存在明显的依赖关系。

实际预警里，我会把条件概率做成一张表，横轴是SPI的等级，纵轴是SSI的等级，每个格子填条件概率。这样决策者一眼就能看出"当前气象干旱等级下，农业干旱升级的风险有多大"。

4.3 重现期计算：单变量与联合的区别

单变量重现期 T = 1/(1-F(x))，这个大家都熟。联合重现期有两种：

• 联合重现期（且）：T_and = 1 / (1 - u - v + C(u,v))，对应"或"事件的补集
• 联合重现期（或）：T_or = 1 / (1 - C(u,v))

等等，这里要特别小心。不同的文献对重现期的定义有差异，有的用P(且)的倒数，有的用P(或)的倒数。我建议在论文里直接写清楚公式，不要只写"联合重现期"四个字。

T_single_spi = 1 / (1 - u0); T_single_ssi = 1 / (1 - v0); T_and = 1 / (1 - u0 - v0 + P_and); % 注意这里用的是"或"事件的补 T_or = 1 / (1 - P_and);

实际算的时候，通常固定一个变量的重现期（比如SPI的10年一遇），然后算另一个变量在不同重现期下的联合重现期，画成等值线图。这张图是干旱风险评估报告里的标配。

4.4 用MATLAB绘制联合概率等值线图

等值线图能直观展示两个变量在不同组合下的联合概率分布。核心是构造网格，逐点算Copula值，然后contour。

u_grid = linspace(0.01, 0.99, 100); v_grid = linspace(0.01, 0.99, 100); [U, V] = meshgrid(u_grid, v_grid); P_joint = zeros(size(U)); for i = 1:size(U,1) for j = 1:size(U,2) P_joint(i,j) = copulacdf('Clayton', [U(i,j), V(i,j)], theta_hat); end end contour(U, V, P_joint, [0.05 0.1 0.2 0.3 0.5], 'LineWidth', 1.5); xlabel('SPI分位数'); ylabel('SSI分位数'); colorbar;

提示：copulacdf在循环里逐点调用效率很低，100x100的网格要算10000次。如果网格再密一点，建议自己写Clayton的CDF向量化计算：C(u,v) = (u^(-θ) + v^(-θ) - 1)^(-1/θ)。一行就能算完整个矩阵。

5. 实操中踩过的坑与经验总结

5.1 边缘分布拟合不好会直接毁掉Copula结果

这是我踩过最深的坑。早期做实验时，SPI序列直接用正态分布拟合，没做检验，结果Copula参数估计出来偏大，联合概率被高估。后来发现SPI序列在干旱年份有明显的负偏，正态分布拟合的CDF在左尾偏小，导致转换后的u值在低分位处偏离均匀分布，Copula误把这种偏离当成了依赖结构。

解决办法很简单：边缘分布一定要做拟合优度检验，而且要在尾部重点检查。可以画PP图或QQ图，看尾部点是否贴合。如果尾部拟合不好，考虑用非参数方法——直接用经验CDF转换，虽然损失了一点光滑性，但避免了分布假设错误。

5.2 样本量不足时参数估计的不确定性

30年的年尺度数据只有30个点，用来估计Copula参数其实偏少。参数估计的标准误可能达到估计值的20%以上。这种情况下，建议用Bootstrap方法给出参数的置信区间：

n_boot = 1000; theta_boot = zeros(n_boot, 1); for b = 1:n_boot idx = randsample(length(u), length(u), true); tau_b = corr(u(idx), v(idx), 'Type', 'Kendall'); theta_boot(b) = 2*tau_b/(1-tau_b); end ci = prctile(theta_boot, [2.5, 97.5]); fprintf('theta的95%%置信区间: [%.3f, %.3f]\n', ci(1), ci(2));

如果置信区间宽得离谱，说明数据量不够支撑联合分析，这时候要么延长序列，要么用月尺度数据增加样本量（但要注意月尺度数据的自相关性会影响显著性检验）。

5.3 月尺度数据的自相关处理

用月尺度SPI和SSI做Copula，样本量能到几百，但月与月之间的自相关很强，直接做Bootstrap会低估不确定性。处理方法有两种：一是用块Bootstrap（block bootstrap），块长取自相关衰减到不显著时的滞后阶数；二是先对序列做去趋势和去季节处理，再用残差做Copula。

我一般用块Bootstrap，块长取12个月（覆盖一个完整年周期），这样能保留季节内的依赖结构。

5.4 不同Copula选出来的结果差异有多大

实测下来，Clayton和Gumbel在干旱数据上的联合概率差异可以到10%到20%，尤其是当tau在0.4到0.6之间时。所以Copula选择不是走过场，AIC差个2到3就要认真比较。如果两个Copula的AIC很接近，建议把两个结果都报出来，说明结论对Copula选择不敏感，这样更稳妥。

5.5 重现期等值线图的解读陷阱

联合重现期等值线图上，同一条线上的点对应的联合重现期相同，但单变量重现期可能差别很大。比如10年一遇的联合重现期线上，可能有一个点是SPI的5年一遇配SSI的20年一遇，另一个点是SPI的20年一遇配SSI的5年一遇。解读的时候必须结合单变量重现期一起看，不能只看联合重现期。

另外，等值线图的外推要谨慎。Copula在数据边缘区域的拟合效果通常较差，如果等值线延伸到样本范围之外，那部分结果只能作为参考，不能作为定量依据。

6. 从联合概率到实际风险评估的落地思路

算完联合概率和重现期，最终要落到风险评估上。我通常的做法是：先根据历史灾损数据确定一个临界阈值（比如SSI≤-1.5对应减产10%以上），然后算这个阈值下不同SPI等级的条件概率，做成风险矩阵。再结合预报的SPI，就能给出农业干旱的风险等级。

这套流程在MATLAB里可以封装成函数，输入SPI和SSI序列，输出Copula参数、联合概率矩阵、条件概率矩阵和重现期等值线图。封装好之后，换一个站点只需要改数据路径，几分钟就能出结果。

有个细节值得注意：不同站点的最优Copula可能不同。我在北方某站点做的时候Clayton最优，换到南方湿润区，Frank反而更好。所以不要写死Copula类型，每次都要重新做拟合优度检验。这也是为什么我坚持把Copula选择做成自动化流程的一部分，而不是手动指定。

最后分享一个实用技巧：如果要做多个站点的对比分析，建议把每个站点的tau值和最优Copula类型画在一张空间分布图上。tau高的区域说明气象干旱和农业干旱的耦合关系强，这些区域更适合做联合风险评估；tau低的区域，联合分析带来的信息增量有限，不如单独分析。这个判断能帮你把有限的精力集中在真正需要联合分析的区域上。