2000—2023 年逐月网格数据nc、tif 混在一起要算 SPEI 的 1、3、6、12 个月结果这是我在各种干旱评估、水资源分析项目里最常接到的需求。标题看着简单真正跑起来全是细节时间维和空间维怎么对齐潜在蒸散怎么估算每个格点的分布参数怎么拟合多个时间尺度怎么循环才不把机器跑死任何一步没想清楚结果拿出来都不敢写结论。这篇文章把整条链路重新捋一遍也把我在实际项目里踩过的坑放在里面适合在气象、生态、农业遥感背景里做数据处理的人。SPEI 的核心是“水分子亏缺”但它又不只是把降水标准化。下面先从原理开始讲理解清楚后再上手比直接抄代码要稳得多。1. 先想清楚SPEI 到底在算什么1.1 从 SPI 到 SPEI缺的不只是水的输入SPI标准化降水指数只拿降水序列做分布拟合然后标准正态化。它有一个天然问题后半段气温上升导致蒸发需求增大但降水不变时SPI 看不出任何异常真正到了田间作物已经感受到水分胁迫了。SPEI 把“潜在蒸散”PET引进来用月降水减去月 PET得到月水分差额D P - PET之后把 D 按时间尺度累加再做分布拟合和标准正态变换。这样做的直接好处是同样的降水量在高温年份会得到更低的 SPEI 值干旱程度会被更早暴露出来。对于 2000 到 2023 年这种气温明显抬升的序列SPEI 比 SPI 更有参考价值。水分子亏缺单位通常用毫米mm。P 和 PET 必须同单位、同空间范围否则后面算出来的指数没有任何可比性。这是我反复强调的第一个前置条件。1.2 数据形态与计算框架我的输入数据一般有两类一类是 netCDF 文件纬度、经度、时间都在变量里另一类是每个月一张的 GeoTIFF时间维靠文件名隐式表达。很多时候项目交付物里两者混着来比如降水是 nc气温是 tif就需要先做“统一化”。统一化包括三件事所有影像重采样到同一网格范围和分辨率时间序列按相同的年份、月份对齐缺失月份必须显式补 NaN单位统一降水如果给的是 mm/月就直接用气温给的是 K 就减 273.15PET 单位必须是 mm/月。我习惯的数据组织方式是一张[time, lat, lon]的三维矩阵time 从 2000-01 依次排到 2023-12共 288 个月。所有格点计算都在这套矩阵上进行输出时再按需要写成 nc 或 tif。2. 核心算法与 MATLAB 实现细节2.1 潜在蒸散 PETThornthwaite 是工程上的合格选择PET 的方法很多Penman-Monteith 精度好但要求风速、湿度、辐射等一堆输入网格尺度上往往拿不到连续 24 年的完整资料。Hargreaves 只需要最高最低温也不错但最常用、和 SPEI 配套最方便的还是 Thornthwaite。Thornthwaite 的月尺度公式是PET 16 * (10 * T / I)^a其中 T 是月均温摄氏I 是年热量指数I Σ ((T_m / 5)^1.514)逐月累加T 小于等于 0 的月份按 0 处理指数 a 是 I 的三次多项式a 6.75e-7 * I^3 - 7.71e-6 * I^2 0.01792 * I 0.49239原始方法还需要根据纬度和月份做日照时数校正。如果只做区域相对比较不做绝对蒸发量标定很多实际项目会省略这一步但如果要发表结果我建议至少保留日长校正系数否则算出来的 PET 在冬季高纬度会被严重低估。在 MATLAB 里我一般不写逐格点循环而是把温度矩阵整体向量化。下面这个函数简化成“按年求热量指数、再算 PET”function PET calc_pet_thornthwaite(T, lat, month) % T: 单月温度矩阵单位℃ [lat, lon] 或 [lon, lat] % lat、month 用于日长校正这里仅保留骨架 T(T 0) 0; if size(T, 3) 12 % 如果是整年 12 个月温度先算热量指数 I sum((T / 5) .^ 1.514, 3); else error(请传入整年温度后逐月计算); end I(I 1) 1; a 6.75e-7 * I.^3 - 7.71e-6 * I.^2 0.01792 * I 0.49239; PET 16 .* (10 .* T ./ I) .^ a; end注意这里的T是当前月份温度但I必须用全年温度算不能只用当月。否则冬季热量指数偏小PET 会被杠杆放大整个序列都会变形。2.2 多时间尺度累加用移动求和而不是循环滚动SPEI-3 表示的是“截至当月往前推 3 个月”的水分子累计值。比如 2023 年 6 月的 SPEI-3对应 2023 年 4 月到 6 月的 D 之和。在 MATLAB 里最直观的写法是for t k:nTime Dk(t, :, :) sum(D(t-k1:t, :, :), 1); end这个循环在小规模数据上没问题但几千个格点、多个尺度会非常慢。我更推荐直接按矩阵操作先生成一个临时列然后用movsum做尾部累加Dk movsum(D, [k-1, 0], Endpoints, discard);这里[k-1, 0]的意思是“往前看 k-1 步往后看 0 步”也就是尾部窗口。注意不能写成movsum(D, k)那是中心窗口会把未来数据混进来结果等于“用了未来信息”在地学指标里是绝对不允许的。movsum丢弃前 k-1 个不完整窗口输出序列会比原序列短。为了后续按月份做分布拟合我更喜欢保留原始时间轴把不足 k 个月的位置写成 NaNDk NaN(size(D)); Dk(k:end, :, :) movsum(D, [k-1, 0], Endpoints, discard);这样时间索引始终和原始月份对齐后面对“年份、月份”做匹配时不会错位。2.3 分布拟合和标准正态变换这是 SPEI 最容易翻车的一步SPEI 的标准化步骤是对 D 序列选择三参数 log-logistic 分布做拟合得到累积概率后再用标准正态分布的逆函数转成 SPEI 值。为什么要用 log-logistic因为 D 的分布经常是偏态的存在比较严重的左尾或右尾三参数 log-logistic 比伽马分布更灵活对负值也能处理。很多人直接拿gamfit拟合 D遇到 D 有负值就会报错或得到 NaN因为这个分布的支撑集是正数。还有一个常见做法是给 D 整体加一个常数再拟合 gamma这是“骗自己”拟合结果在负数段完全是外推算出来的极端干旱根本不可信。我不建议在 MATLAB 里追求完全复刻某个现成统计包但至少要实现“三参数分布 连续概率水位转换”。下面是我在批量脚本里使用的核心子函数先做 log-logistic 的参数估计再做 CDF 和逆正态变换function p loglogistic_cdf_fit(x) % x: 某固定月份、某固定格点、所有年份的 D_k 序列 x x(~isnan(x) isfinite(x)); x x(:); if numel(x) 5 p NaN(size(x)); return; end theta fit_loglogistic_mle(x); g theta(1); % 位置参数 a theta(2); % 尺度参数 b theta(3); % 形状参数 z (x - g) / a; p 1 ./ (1 z .^ (-b)); % 避免 0 和 1 两端饱和否则 norminv 会返回 Inf p max(min(p, 1 - 1e-10), 1e-10); end function theta fit_loglogistic_mle(x) % 负对数似然函数直接写成一个子函数 theta0 [min(x) - 0.1 * range(x), 0.5 * range(x), 1.0]; opts optimset(MaxIter, 300, Display, off); theta fminsearch((th) negll_loglogistic(th, x), theta0, opts); end function nll negll_loglogistic(th, x) g th(1); a th(2); b th(3); if a 0 || b 0 || g min(x) nll 1e15; return; end z (x - g) / a; if any(z 0) nll 1e15; return; end nll -sum(log(b / a) (b - 1) * log(z) - 2 * log(1 z .^ b)); end得到累积概率 p 后SPEI 就用norminv转换SPEI norminv(p);p 如果等于 0 或 1norminv会给出 -Inf 或 Inf通常会在插值、制图时把极端值截断到 -3 和 3 之间。不要小看这个截断很多出图“满屏红色”的鬼故事都是 Inf 没处理。如果格点太多fminsearch在每一个格点、每一个月份、每一个尺度上跑非常耗时。实际我通常把格点展成一维用parfor按格点并行Dmat reshape(D, nTime, []); % 时间 × 格点 nCell size(Dmat, 2); SPEI_mat NaN(size(Dmat)); parfor c 1:nCell for mi 1:12 rows find(miMonth mi ~isnan(Dmat(:, c))); if numel(rows) 5 xs Dmat(rows, c); p loglogistic_cdf_fit(xs); SPEI_mat(rows, c) norminv(max(min(p, 1-1e-10), 1e-10)); end end end这里miMonth是每个时间点所属的月份1 到 12。分月拟合非常关键因为不同月份的 D 分布均值、方差都不一样混在一起拟合等于把冬夏水分差异全部搅浑。3. 完整实操记录从读取数据到写出 TIFF3.1 nc 和 tif 的读取与统一网格先读 nc。我最常用的是ncdisp先看变量结构不要上来就ncreadncdisp(pre_2000_2023.nc);常见的维度顺序有(time, lat, lon)也有(lon, lat, time)。如果是后者读取后要permuteP_raw ncread(pre_2000_2023.nc, pre); % [lon, lat, time] P permute(P_raw, [3, 2, 1]); % [time, lat, lon]同时把经纬度读出来lon ncread(pre_2000_2023.nc, lon); lat ncread(pre_2000_2023.nc, lat);再看时间变量。很多 nc 的时间单位是“days since 1900-01-01”需要换算到年月。我这里假设文件里已经按自然月排列如果不放心可以用datetime转换后检查首尾t0 datetime(2000, 1, 1); tVec t0 calmonths(0:287);再读 tif。如果每个月一个文件文件名形如T_200001.tif用readgeoraster逐月读入fileList dir(T_*.tif); for i 1:numel(fileList) [A, R] readgeoraster(fullfile(fileList(i).folder, fileList(i).name)); % A 是矩阵R 是地理参考对象 endreadgeoraster在老版本 MATLAB 里对应geotiffread。不要先写代码再检查一定要先看一眼数据大小和形状。我遇到过 tif 里第一列是 Y 方向、第二列是 X 方向和 nc 的经纬度方向完全反着的案例不检查后面全白算。如果 nc 和 tif 网格不一致最省事的办法是把 tif 插值到 nc 的经纬网格上。用griddedInterpolant或者interp2都行关键是先确认插值方向不会把经纬度弄反[lonG, latG] meshgrid(lon, lat); [T_grid, R_grid] ... % 对每个 tif 月份做二维插值 T_interp interp2(lonTif, latTif, T_tif, lonG, latG, linear, NaN);这里的lonTif,latTif可以从R对象里用latlim、lonlim和矩阵大小生成。插值后把 288 个月叠成[time, lat, lon]的 T 矩阵和 P 对齐。3.2 按格点批量计算的循环结构处理完数据后我习惯把空间维拍平用二维矩阵做主体计算省掉三层循环里的两层。nTime size(P, 1); nLat size(P, 2); nLon size(P, 3); % 转成 [time, lat*lon] Pmat reshape(P, nTime, nLat * nLon); Tmat reshape(T, nTime, nLat * nLon);然后逐格点算 PET。这一步可以按照整年算因为 Thornthwaite 的热量指数是全年累积的。我先把 Tmat 整成年份块Tyear reshape(Tmat, 12, 24, nLat * nLon); % 12月 × 24年 × 格点 Tyear(Tyear 0) 0; I squeeze(sum((Tyear / 5) .^ 1.514, 1)); % 24 × 格点 a 6.75e-7 * I.^3 - 7.71e-6 * I.^2 0.01792 * I 0.49239; I(I 1) 1;再按月份生成 PET 矩阵PETmat zeros(size(Tmat)); for yr 1:24 for mo 1:12 idx (yr - 1) * 12 mo; Iyr I(yr, :); ayr a(yr, :); Tm Tmat(idx, :); Tm(Tm 0) 0; PETmat(idx, :) 16 .* (10 .* Tm ./ Iyr) .^ ayr; end end这里我没有加日长校正严格写论文时要补上但如果只是做趋势分析和干旱区划这个量级通常够用。之后 D 就是Dmat Pmat - PETmat;Dmat 中凡是原始数据缺测的地方都要保持 NaN不能硬填 0。填 0 会把“数据缺失”伪装成“水分子平衡为 0”后面拟合出的分布完全错掉。核心尺度循环写成这样scales [1, 3, 6, 12]; SPEI_stack cell(1, numel(scales)); for s 1:numel(scales) k scales(s); Dk NaN(nTime, nLat * nLon); Dk(k:end, :) movsum(Dmat, [k-1, 0], Endpoints, discard); % 对每个格点和月份做 log-logistic 拟合 SPEI_stack{s} calc_spei_by_month(Dk); endcalc_spei_by_month就是把前面讲的“分月拟合 norminv”包装成函数。每一个尺度都会生成和原始时间对齐的 288 个月 SPEI 序列月份空间维不变。3.3 结果输出为 GeoTIFF / netCDF算完后把矩阵 reshape 回原来的网格再写文件。写 GeoTIFF 可以用SPEI_grid reshape(SPEI_stack{s}, nTime, nLat, nLon); % 输出 2023 年 12 月尺度 outFile SPEI_12_202312.tif; geotiffwrite(outFile, SPEI_grid(end, :, :), R, CoordRefSysCode, 4326);如果geotiffwrite提示R格式不对检查一下是不是readgeoraster返回的R是mapref对象而geotiffwrite需要georefcells或者世界文件。遇到这种兼容问题我一般直接用geotiffwrite(outFile, data, R, GeoKeyDirectoryTag, ...)绕开。写 nc 也同样简单nccreate(SPEI_6_2000_2023.nc, spei, Dimensions, ... {time, nTime, lat, nLat, lon, nLon}); ncwrite(SPEI_6_2000_2023.nc, spei, SPEI_grid);写好之后一定要读回来做一次“反演对比”随机取一个格点把 nc 里读出的序列和内存里的序列画在一起完全重合才算过。这一步花不了两分钟但能挡住 90% 的维度错位。4. 常见问题与排查实录4.1 时间尺度不对齐SPEI-3 的 1 月为什么总是缺值这是我最常被问到的问题。SPEI-3 的 2000 年 1 月需要 1999 年 11 月到 2000 年 1 月的 D如果数据从 2000 年 1 月开始那 2000 年 1 月的 SPEI-3 本来就无法计算只能给 NaN而不是用 2000 年 1 月三个月的数据凑。解决办法也很粗暴提前多准备 12 个月数据或者把结果序列截断到 2000-03 以后。实际项目中官方常常用“2000-2023”表示计算区间但内行人会写“有效结果区间为 2000-03 到 2023-12尺度为 3”。这个细节在成果说明里写清楚能避免很多不必要的返工。4.2 PET 单位没统一导致全区域极端干旱有次排查某区域 SPEI 结果发现连年都是 -3 以下的极端干旱地图红成一片。检查后发现 PET 数据来自一个日值物理过程模拟单位是 mm/天我把日值直接当成月值用了。降水是 mm/月PET 是 mm/天乘 30水分亏缺瞬间爆表。凡是跨数据源做 D P - PET我强烈建议第一件事就是打印三个数的统计量P 的年累计范围、PET 的年累计范围、D 的年累计范围。P 和 PET 的年累计量应该在同一个数量级如果 PET 年累计只有几十毫米那肯定单位或公式有问题。4.3 分月拟合时样本太少参数估计爆掉2000—2023 只有 24 个样本分到每个月后某个月只有 24 个点。三参数 log-logistic 的fminsearch在小样本下可能不收敛或者收敛到极端参数。我的处理是样本数少于 5 直接输出 NaN对 CDF 做截断防止norminv出 Inf如果某个格点在所有年份都是零降水或者无变化人工剔除不参与拟合可以用parfor配合并行池但要注意随机数种子否则每次结果可能会因初值不同而轻微变化。还有一个更隐蔽的问题当 D 序列里有大量重复值比如长时间完全干旱导致 D 恒为一个常数log-logistic 的尺度参数会趋近于 0CDF 变成阶跃函数。此时我会退一步用经验分布tiedrank代替参数分布p (tiedrank(x) - 0.5) / numel(x);这个“经验分布兜底”不是偷懒工程上能避免一堆统计病态结果。4.4 投影坐标系不一致导致 tif 写不进nc 的坐标通常是经纬度tif 可能是 UTM 投影坐标。如果把经纬度结果直接geotiffwrite给一个 UTM 参考会得到一张位置完全错掉的文件。我的习惯是先在全球范围全部转成 WGS84 经纬度输出再让后续 GIS 工具自行投影。如果项目必须输出 UTM用一个统一的目标 EPSG比如 32650一次性重投影。还有一个常见但不明显的坑geotiffwrite的矩阵行列顺序是“行号对应北纬/上方向”而某些 nc 文件的 lat 是倒序从北往南。写入 tif 前要先用flipud把纬度方向翻正不然出图上下颠倒。4.5 多时间尺度结果之间要自洽1、3、6、12 个月尺度的 SPEI 反映的是不同时间窗口的干旱特征它们不应该互相矛盾但会有滞后和幅度差异。我每次批量计算完之后会随机选 5 个格点把四个尺度画在同一张图上检查相位关系。SPEI-1 波动最剧烈SPEI-12 相对平滑如果 SPEI-12 比 SPEI-1 还毛刺说明时间尺度累加的窗口方向写反了大概率又是movsum用成了中心窗口。这个检查非常廉价却能一眼暴露累加逻辑问题。做大规模网格数据时宁可多画几张诊断图也不要等到交付之后等着被技术审核打回来。最后再分享一个小技巧批量计算 SPEI 最好不要临时在 MATLAB 里手敲整套代码而是把“读取→统一网格→计算 PET→多尺度 D→分月拟合→写 tif”封装成一个脚本只通过 config 参数切换输入路径、变量名和时间范围。我自己的脚本就是一套通用模板换一个区域、换一套数据改三行配置就能跑。这样做的好处不只是省时间更重要的是过半年后想复现当时的结果你还能清楚记得每个环节用的是什么数据、什么边界条件。对做气候指数分析的人来说可复现可能比算得准更重要。