干旱研究里有个绕不开的坎气象干旱和农业干旱明明是两个不同维度的东西一个看降水亏缺一个看土壤水分和作物响应但它们之间又存在明显的传导和滞后关系。单独算一个干旱指数的频率分布只能回答今年降水偏少到什么程度回答不了这种气象条件下农业干旱发生的可能性有多大。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族因为形式简单、参数少、能刻画非对称依赖在干旱领域占了绝大多数。Copula类型依赖特征适用场景参数个数Clayton下尾依赖强上尾弱气象-农业干旱同时偏枯1Gumbel上尾依赖强下尾弱丰水年洪涝联合1Frank对称依赖无尾部依赖依赖结构较温和1Gaussian对称无尾部依赖线性相关为主1相关系数t对称双尾依赖极端事件双向关联2相关自由度干旱是越干越相关的典型下尾事件所以Clayton通常是首选。但不要凭直觉定后面会讲怎么用拟合优度检验来选。1.3 联合概率、条件概率、重现期分别对应什么实际含义这三个指标是干旱风险评估的核心输出含义完全不同联合概率P(SPI ≤ a 且 SSI ≤ b)气象干旱和农业干旱同时达到某个等级的概率。用于评估复合干旱风险。条件概率P(SSI ≤ b | SPI ≤ a)已知气象干旱达到某等级时农业干旱也达到某等级的概率。这是预警里最有用的指标因为它回答了气象干旱已经发生了农业干旱会不会跟上。联合重现期在给定联合概率下事件平均多少年发生一次。注意重现期有且和或两种定义算出来的数值差别很大论文里必须写清楚用的是哪种。2. 数据准备与边缘分布拟合的实操细节2.1 干旱指数的选择与时间尺度匹配气象干旱用SPI农业干旱用SSI或SMI这是主流做法。但有个容易被忽略的点时间尺度必须匹配。SPI-33个月尺度对应的是短期降水亏缺而土壤湿度对降水的响应通常滞后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); % h0表示不能拒绝原假设拟合可接受 % 对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或1Copula的对数似然函数会变成无穷大。处理方法是在边界处做微小截断比如把小于1e-6的值设为1e-6大于1-1e-6的设为1-1e-6。这个操作对结果影响极小但不做的话程序直接报错。2.3 用Kendall秩相关系数初步判断依赖强度在选Copula之前先算一下Kendalls tau和Spearmans 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和BIClogL -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≤-1u0≈0.159时SSI≤-1v0≈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会低估不确定性。处理方法有两种一是用块Bootstrapblock 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低的区域联合分析带来的信息增量有限不如单独分析。这个判断能帮你把有限的精力集中在真正需要联合分析的区域上。