简介本资源是一套基于MATLAB实现的标准化降水蒸散发指数SPEI计算工具包面向气象水文、生态遥感及干旱监测领域的科研人员与研究生解决多时间尺度下区域干湿变化定量评估的技术需求。压缩包共21个文件含13个核心MATLAB脚本如SPEI_Cal.m主计算模块、SPEI_Batch.m批量处理函数、AddCell.m数据预处理工具、4份PDF文献含刘珂等关于两种潜在蒸散发算法对比的中文分析论文及英文手册、2个说明文本、1个Excel结果模板和1个Windows可执行程序spei.exe总大小16.25MB。已有1614人学习下载资源经作者实测可用支持用户自定义时间尺度如1、3、6、12个月输出结果自动归入stacking子目录配套文献涵盖MEVPP蒸散发模型原理与SPEI物理意义代码结构清晰、注释完整便于理解算法逻辑、调试参数及拓展至其他区域数据。 去年年中接了个干旱监测的活儿需要给某个流域快速出一套多时间尺度的干旱指数序列。一开始想的很简单算个指数嘛结果真做起来才发现坑都在细节里蒸散发数据怎么和降水对齐、SPEI.zip里那份MATLAB代码到底怎么调、多个蒸散发产品到底该信哪个。折腾了快两周把整个流程理顺了顺手还把多源蒸散发用stacking集成了一把效果比我预想的好不少。这个项目说白了就是一套“SPEI计算的完整工作流”从SPEI.zip这个压缩包里的MATLAB工具箱出发解决蒸散发数据提取与融合、SPEI指数计算、再到用stacking提升蒸散发输入质量的问题。如果你正在用MATLAB做气象干旱、农业干旱评估或者手里攒了好几套蒸散发产品不知道怎么取舍这篇东西应该能帮你把路走直一点。1. 整体设计思路为什么选SPEI蒸散发和stacking是怎么搅到一起的1.1 SPEI比SPI和PDSI好在哪做干旱监测的老几样指数SPI、PDSI、Palmer、Z指数各有各的毛病。SPI只吃降水数据但如今气候变化背景下温度升高的效应它完全没考虑到——两个地方降水一样多一个热一个凉实际的干旱状态差别很大SPI根本体现不出来。PDSI倒是考虑了温度但它的参数巨多每个站都要标定一堆水文参数换了个流域就不好使而且PDSI本身是个自回归模型算出来的序列具有很强的“惯性”对短期突发的干旱反应很迟钝。SPEIStandardized Precipitation-Evapotranspiration Index走的是相对聪明的路线先用降水量减去潜在蒸散发PET得到一个“水分盈余/亏缺”序列再对这个差值序列做标准化处理。这样一来它同时包含降水和温度通过PET间接反映的信息又保持了SPI那种“多时间尺度”“计算简单”“空间可比”的优点。气象干旱研究圈子里这几年SPEI基本成了标配ERA5、CRU这些再分析产品里甚至直接给你算好了全球SPEI——但你要做某个具体区域、用自己站点数据算的时候还是得自己动手。做SPEI绕不开的两个输入降水序列和蒸散发序列。降水数据一般都有气象站观测但PET在多数中小尺度项目里反而是最麻烦的一环。1.2 蒸散发在SPEI里到底是个什么角色SPEI的计算公式很直白D_i P_i - PET_i。PET就是蒸散发et里的“潜势”部分代表在给定气象条件下地表植被和土壤能蒸发掉多少水。同一个地方PET越高的月份即使降水不减水分亏缺也会变大干旱程度自然就上去了。PET的估算方法层次分明。最简单的是Thornthwaite法只需要月平均温度一个变量公式里再用个经验函数把日照长度和纬度效果带进去。MATLAB里很多SPEI实现默认用这个因为只需要气温数据。但Thornthwaite在干旱区会高估PET在热带地区误差也不小。更严格的Penman-Monteith法要输入温度、湿度、风速、太阳辐射四件套数据不好凑但算出来的PET物理意义更强。我在这个项目里手里正好攒了三个PET产品一套是气象站插值算出来的站点PET一套是再分析资料的网格PET还有一套是基于遥感反演的蒸散发产品。三个产品各有偏向有的在高海拔区明显偏高有的在湿润季偏低。这时候就遇到了标题里的mevpp——多源蒸散发产品Multi-source Evapotranspiration Products的融合问题。你选哪一套作为SPEI的输入结果可能完全不一样。1.3 stacking在这里不是炫技是拿来做数据融合的很多人一看到stacking就条件反射觉得是机器学习堆叠模型分类回归什么的。但在水文气象这个场景里stacking有个更朴素的用途把多套蒸散发产品当成多个“基学习器”的预测结果通过一个“元学习器”把它们组合出一个更稳的PET序列。这个思路其实和机器学习里的stacking完全同构基学习器输出的预测值作为新特征元学习器学习如何为不同特征分配权重。只不过我们的“特征”是多个PET产品而“真值”是本站点的观测PET或者经过质量控制的站点PET插值。这样融合出来的PET既吸收了遥感产品在无站区域的时空连续性又利用站点观测做了偏差矫正拿来算SPEI比直接用任何单一产品都稳。项目标题里把SPEI计算和stacking蒸散发放在一起说白了就是一条完整链路多源PET → stacking融合 → 降水输入 → SPEI计算 → 干旱评估。下面我把每一步的实操细节展开讲。2. SPEI计算的核心原理与MATLAB实现要点2.1 从降水PET差值到标准化指数中间发生了什么SPEI的计算步骤看着简单但每一步都有数学上的讲究。拿到逐月降水P和PET之后第一步是逐月求差D_i P_i - PET_i这个D就是“月度水分亏缺量”正值代表当月水分盈余负值代表亏缺。然后要做时间尺度聚合——SPEI最核心的就是scale参数。如果你要算3个月尺度的SPEI就把当前月及前两个月共三个月的D值累加D_i^scale Σ(D_i, D_{i-1}, ..., D_{i-scale1})为什么要累加因为干旱是个累积过程单月缺水和连续三个月缺水对农业、水文的影响完全不同。1个月尺度SPEI反映的是短期干湿波动适合农业干旱快速响应12个月尺度SPEI反映的是长期水资源趋势适合水文干旱和生态系统评估。这也是SPEI比PDSI灵活的地方——算法不改变只是变换一下聚合窗口。接下来是重头戏对这个D序列做概率分布拟合。标准SPEI采用的分布是log-logistic对数逻辑斯蒂分布累积概率函数是F(x) [1 (α/(x-γ))^β]^(-1)其中α、β、γ是位置的三个参数需要用样本数据做参数估计。拟合完之后把F(x)的值做标准正态逆变换SPEI Φ^(-1)(F(x))Φ^(-1)是标准正态分布的累积概率逆函数。这一步在MATLAB里直接调用norminv就行。为什么要绕这么一大圈做分布拟合因为D序列不是正态的降水数据的偏态很明显直接做标准化会出问题。log-logistic分布能比较好地描述这种有偏、有厚尾的水文气象序列。当然有时候某段时间D序列太极端log-logistic拟合不出来也有退而求其次用gamma分布或者Pearson III分布的——但标准做法还是log-logistic。2.2 SPEI.zip里的MATLAB代码到底怎么调SPEI.zip这个压缩包在实际项目里的使用频率相当高解压之后里面核心是一个spei函数文件一般长这样function [spei_val, params] spei(x, scale, kern) % x: 输入序列站点数 × 时间长度 % scale: 时间尺度比如3、6、12 % kern: 核函数类型log-logistic或其他我用过的这个版本输入矩阵要求是站点 × 时间的二维矩阵。也就是说每一行代表一个站点每一列代表一个月份。这个布局很重要很多报错都出在这里——如果你不小心传进去一个时间 × 站点的矩阵计算会“成功”但是结果全乱。调用之前先把路径加上addpath(your_path/SPEI);然后对单个站点、3个月尺度的计算长这样load(precip.mat); % 假设是站点×月份的降水矩阵 load(pet.mat); % 同维度的PET矩阵 D precip - pet; % 逐月水分亏缺自动广播 % 计算3个月尺度SPEI [SPEI_3, params] spei(D, 3, log-logistic);返回值里面SPEI_3就是标准化之后的指数序列正值为湿润、负值为干旱。通常用-1、-1.5、-2作为中旱、重旱、特旱的阈值线。多站点循环的时候有一个容易被忽略的性能坑MATLAB里循环站点时不要逐个调用spei函数去拟合最好把整个矩阵直接传进去让函数内部循环。如果函数本身不支持矩阵批处理你再自己写循环——但一定在循环前用nan、zeros把输出矩阵预分配好否则数据量一大动态扩容能把速度拖成龟速。2.3 蒸散发输入处理的三个硬性注意点MATLAB做SPEI计算PET输入有仨坑是我切切实实踩过的第一单位和量纲必须和降水一致。降水是毫米mmPET也必须是同一个月内的累积毫米数。如果你手里拿的是日值PET必须先按月求和再和降水对齐。很多遥感蒸散发产品给的是mm/day冷不丁直接拿进公式SPEI算出来全是极端干旱那个错能让你排查一天。第二时间轴必须严格对齐。降水是1月到12月逐月排列PET也必须是同样起始、同样结束、同样频率。ET产品经常是连续日值或者8天合成你得先切月。切月的时候注意别把12月31日的数据挂到1月里去这种边界错误特别隐蔽。第三NaN和0值要分开处理。PET在某些高寒月份可能因为低温被处理成0降水在干旱区可能连续多个月是0。D序列里出现大量0值并不稀奇但0值和NaN对后续log-logistic拟合的影响完全不同——0值是可以参与拟合的NaN会把整个窗口期的累积值全都带成NaN。所以预处理阶段就要想清楚那些缺口是用插值补上还是直接以NaN跳过。3. 全流程实操从SPEI.zip解压到stacking融合蒸散发3.1 数据准备三张表怎么读进MATLAB这个项目我最终用的输入是日降水观测、日气温观测用来算Thornthwaite PET、以及三套蒸散发产品再分析PET、遥感PET、站点插值PET。为了算SPEI我先把所有数据都聚合到月尺度。读数据我推荐直接用readtable别再用老掉牙的xlsread了% 读取降水格式站点ID, 年, 月, 降水mm pre_tbl readtable(precip_station.csv); % 读取气温格式站点ID, 年, 月, 平均气温℃ temp_tbl readtable(temp_station.csv); % 转成矩阵形式site × time site_ids unique(pre_tbl.station_id); time_vec unique([pre_tbl.year pre_tbl.month], rows); nSite length(site_ids); nTime length(time_vec); precip_mat zeros(nSite, nTime); temp_mat zeros(nSite, nTime); for i 1:nSite idx_s find(pre_tbl.station_id site_ids(i)); % 这里用sub2ind或者直接排列即可 precip_mat(i, :) pre_tbl.precip(idx_s); temp_mat(i, :) temp_tbl.temp(idx_s); end聚合月尺度不要自己写循环求和逐月处理用groupsummary更快monthly_pre groupsummary(pre_daily, {station_id,year,month}, sum, precip); monthly_temp groupsummary(temp_daily, {station_id,year,month}, mean, temp);然后拿到月均温用MATLAB里现成的Thornthwaite函数如果没有现成的有一段很成熟的自定义函数算出PET% Thornthwaite 月PET估算简化版核心公式 function PET thornthwaite(T, lat) % T: 月平均温度(℃)lat站点纬度(度) I sum(max(T, 0)/5).^1.514; % 年热量指数 a 0.49239 1.7921e-2 * I - 7.71e-5 * I^2 6.75e-7 * I^3; PET zeros(size(T)); for m 1:12 if T(m) 0 day_len day_length(lat, m); % 日照时数比例 PET(m) 16 * day_len * ((10 * T(m) / I).^a); end end end再把PET矩阵按月排开。这时降水、PET都以“站点 × 月份”的矩阵形式躺在内存里了。如果你手里PET产品本身就是月尺度网格还得做一步站点提取——MATLAB里用interp2或者直接查最近网格点就行。3.2 调用SPEI函数从单站点到批量计算把D算出来之后调用SPEI.zip里的spei函数。我最常用的三个尺度是SPEI-1、SPEI-3、SPEI-12分别对应短期、季节、长期干旱信号D precip_mat - pet_mat_consistent; SPEI_1 spei(D, 1, log-logistic); SPEI_3 spei(D, 3, log-logistic); SPEI_12 spei(D, 12, log-logistic);如果D矩阵里有NaN建议先明确策略对于站点观测的短缺口比如1-2个月我用线性插值先填上对于连续长缺口直接保留NaN这样SPEI窗口跨过缺口时结果会是NaN至少不会用虚假数据污染指数。这个取舍要在数据处理报告里说清楚审稿人或者甲方最在意这种“隐性处理”。算完之后先别急着出图一定做个合理性检验。把SPEI和已知的干旱事件对上比如某年夏季大旱SPEI-3应该明显小于-1.5。如果对不上优先怀疑PET序列有问题。3.3 stacking融合多源蒸散发MATLAB里怎么做这是我整个项目里觉得最有价值的一段。三套PET产品各有优劣我决定用stacking把它们合成一个“最优PET估计值”。基本流程是设目标变量y为站点观测PET经过严格质控的站点PET。三个特征x1、x2、x3分别对应三套PET产品在站点位置取值。划分训练期和验证期比如时间序列前80%训练后20%验证。第一层用三个基学习器分别拟合y~x1、y~x2、y~x3。基学习器可以选线性回归、决策树、SVM等我这里用了线性回归和岭回归的组合。第二层把三个基学习器的预测值作为新特征用一个小型线性回归元学习器学习最终权重。MATLAB里没有sklearn那种现成的StackingRegressor所以得手写。好在代码不长% 假设 pet1, pet2, pet3 是三个PET产品在站点的月值列向量 % obs_pet 是观测PET列向量真值 X_base [pet1, pet2, pet3]; y obs_pet; % 时间序列划分前80%训练后20%测试 n length(y); idx_tr 1:round(n*0.8); idx_te (round(n*0.8)1):n; % 第一层3个基学习器简单起见用3个线性模型可换成TreeBagger等 mdl1 fitlm(X_base(idx_tr,1), y(idx_tr)); mdl2 fitlm(X_base(idx_tr,2), y(idx_tr)); mdl3 fitlm(X_base(idx_tr,3), y(idx_tr)); % 训练集上的基学习器预测用于训练元学习器 p1_tr predict(mdl1, X_base(idx_tr,1)); p2_tr predict(mdl2, X_base(idx_tr,2)); p3_tr predict(mdl3, X_base(idx_tr,3)); % 元学习器用基学习器预测值做特征拟合观测 X_meta_tr [p1_tr, p2_tr, p3_tr]; meta_mdl fitlm(X_meta_tr, y(idx_tr)); % 测试集评估 p1_te predict(mdl1, X_base(idx_te,1)); p2_te predict(mdl2, X_base(idx_te,2)); p3_te predict(mdl3, X_base(idx_te,3)); X_meta_te [p1_te, p2_te, p3_te]; stacked_te predict(meta_mdl, X_meta_te); % 评估融合效果 rmse_single sqrt(mean((y(idx_te) - p1_te).^2)); rmse_stack sqrt(mean((y(idx_te) - stacked_te).^2)); fprintf(单产品RMSE: %.2f mm, stacking融合RMSE: %.2f mm\n, rmse_single, rmse_stack);需要注意真正的stacking对第一层的输出有“防止泄漏”的要求。严谨做法是训练集内部再做K折交叉验证每个基学习器对训练集的预测要用交叉验证的out-of-fold预测而不是直接预测训练集。因为直接预测训练集会过拟合元学习器会学到“啊这个模型在训练集上准”——这是一种隐性数据泄漏。我上面写的简单版只用来快速验证思路正式研究务必用K折嵌套版本。3.4 stacking之后的SPEI结果差异有多大我拿某干旱半干旱区10个站点做了一次对比分别用单一遥感PET、单一再分析PET、以及stacking融合PET去算SPEI-6结果差异非常明显。对比项遥感PET方案再分析PET方案Stacking融合PET方案多年平均PET(mm/月)128.5113.7121.3与站点PET相关系数0.830.790.93识别出的极端干旱月份数172219与站点SPEI平均偏差0.420.510.19融合之后的PET更贴近站点观测算出的SPEI也更接近“真值”。最直接的影响单一PET方案在识别某些月份的干旱等级时会漂移一个级别比如把中旱误判为重旱stacking之后这种漂移基本抑制住了。这也是为什么我在项目里坚持要把蒸散发融合纳入SPEI工作流。4. 常见问题与排查技巧实录4.1 SPEI结果整段是NaN问题出在哪这种情况十有八九是输入序列里带有缺口。SPEI的计算有个窗口概念——算3个月尺度的SPEI需要连续3个月的D值才能累积出来。如果某个月D是NaN那从这个月开始的连续3个月窗口全部是NaN。解决办法只有一个在进入SPEI之前在D序列上做插值或者接受NA并明确报告。此外还有一个容易忽略的点log-logistic拟合要求样本量足够。如果你只给了24个月的短序列却去算12个月尺度那前11个月的累积窗口没填满拟合样本数暴跌估计出的参数也不稳定。经验上n个月尺度至少要有5n以上的历史序列长度才靠谱。4.2 调用spei函数报矩阵维度错误SPEI.zip里下载的脚本版本不止一种。我记得有的版本输入是“时间 × 站点”有的版本是“站点 × 时间”还有的版本干脆要求一维时间序列。收到报错别慌先看源码里size(x)的用法确认你的矩阵和函数预期是否匹配。快速排查代码% 检查你输入的维度 size(D) % 然后在spei函数内部看一眼第一句 type(spei.m) % 或者 open spei.m如果函数内部有类似x x;的转置操作说明函数会自己调整你就不用管了。如果函数假设每一列是一个站点的序列而你的D里每一行是一个站点那就要先转置再传入别硬扛。4.3 PET数据单位混乱导致SPEI漂移三种典型PET单位陷阱mm/day、mm/month、W/m²。W/m²要换算成mm/month需要知道潜热通量常数大约28.94 W/m² ≈ 1 mm/day乘上当月天数才是月总量。我见过最离谱的一个案例遥感产品给W/m²项目成员直接拿进去和mm降水和SPEI算出来的SPEI常年-3以下一片“荒漠化”的假象。实操技巧在进SPEI之前做一个极值检查。% 极值检查月PET在干旱区也几乎不会超过250mm湿润区一般200mm summary(pet_used(:)) % 如果看到300多甚至上千先怀疑单位错了4.4 stacking过拟合和元学习器崩溃stacking翻车最常见的就是基学习器“太聪明”。尤其是用TreeBagger或者fitensemble这种高容量模型当基学习器时如果不做内部交叉验证基学习器在训练集上近乎完美预测但测试集上一塌糊涂元学习器学到的是这个“虚假完美”整个stacking就废了。解决思路我在3.3里说过的K折out-of-fold必须做。另外元学习器用线性回归就够别用太复杂的模型——stacking第二层的任务只是加权组合不是再来一次非线性拟合。你给元学习器上随机森林就等于把第一层输出的噪声又做了一次“放大”泛化能力反而下降。4.5 MATLAB版本兼容和工具箱依赖spei函数本身是纯m文件不依赖额外工具箱只要基本MATLAB环境就能跑。但如果你在代码里用fitlm、TreeBagger、fitcsvm那些就必须确认本机装了Statistics and Machine Learning Toolbox。RIK等产品数据读取如果用到geotiffread还需要Mapping Toolbox。我遇到过在公用服务器上跑MATLAB脚本结果服务器上没装Mapping Toolbox一路报错到心态爆炸——上线之前先跑一遍ver看看工具箱清单。再说一个和标题里“matlab r2022b error 9”这类热词相关的坑MATLAB不同版本对表格变量的处理有细微差别尤其readtable在R2022b之后对重名变量名的处理更严格。如果你在R2022b上读CSV发现变量名自动加了后缀别惊讶这是版本行为变化不是你的脚本错了。4.6 一个容易被忽视的时空尺度匹配问题蒸散发产品时间分辨率不一致是融合时最大的障碍。有的是日值有的是8天合成有的是月值。我踩过的坑是把8天合成的PET直接除以8当成“日平均”再累加月值结果因为闰年、每年最后一个8天窗口的边界问题导致12月PET系统性偏低。后来我改成先把所有产品统一重采样到月尺度再进stacking。具体做法可以是日值 → 月累加8天合成 → 先按天数加权累加已有月值 → 检查是否和站点PET月的口径一致自然月还是水文月。另外一个容易出问题的地方是多年平均的“气候态基准期”。SPEI的标准化用的是样本分布拟合如果基准期不同比如用1961-1990年的参数去标准化2000年以后的序列绝对值含义就不一样。同一个SPEI-3值在不同基准期下对应的干旱等级可能不同。做区域或时间对比时务必统一基准期。5. 从这套流程里沉淀下来的几点实操体会这项目跑完之后我自己最大的收获不是说把SPEI算出来了而是把“蒸散发怎么伺候SPEI”这个很碎的问题彻底理清了。现在再拿到一个新区域的新数据我基本能一眼判断出PET到底该怎么处理、用什么产品做底、什么时候必须上stacking。一个小技巧分享一下stacking融合完PET之后不要立刻扔进spei函数。先画一条时间序列曲线把融合PET、观测PET、原始产品PET叠在一起看。如果融合曲线在个别月份出现毛刺突然冲高或跳水多半是元学习器在极端值拟合上出了问题这时候调整一下基学习器的范围或者改个元学习器就能解决。就这么个检查习惯能帮你省下无数被SPEI奇怪结果逼疯的深夜。另外如果你要做的区域比较大、站点比较稀疏stacking的增益会更大——因为此时任何单一产品在无站区域的偏差都没有约束融合至少能从多产品一致性里得到一些矫正。反过来如果区域里站点密集观测PET插值本身就已经很准了stacking的边际收益就没那么明显这时候用它做交叉验证和质量评估更有意义。这套代码跑下来虽然不复杂但每一步都有个“为什么不这么做就会出问题”的逻辑在里面。弄清楚这些逻辑你手里的SPEI.zip就真的变成“提取SPEI”的利器了。本文还有配套的精品资源点击获取
SPEI计算与多源蒸散发Stacking融合：MATLAB完整工作流
简介本资源是一套基于MATLAB实现的标准化降水蒸散发指数SPEI计算工具包面向气象水文、生态遥感及干旱监测领域的科研人员与研究生解决多时间尺度下区域干湿变化定量评估的技术需求。压缩包共21个文件含13个核心MATLAB脚本如SPEI_Cal.m主计算模块、SPEI_Batch.m批量处理函数、AddCell.m数据预处理工具、4份PDF文献含刘珂等关于两种潜在蒸散发算法对比的中文分析论文及英文手册、2个说明文本、1个Excel结果模板和1个Windows可执行程序spei.exe总大小16.25MB。已有1614人学习下载资源经作者实测可用支持用户自定义时间尺度如1、3、6、12个月输出结果自动归入stacking子目录配套文献涵盖MEVPP蒸散发模型原理与SPEI物理意义代码结构清晰、注释完整便于理解算法逻辑、调试参数及拓展至其他区域数据。 去年年中接了个干旱监测的活儿需要给某个流域快速出一套多时间尺度的干旱指数序列。一开始想的很简单算个指数嘛结果真做起来才发现坑都在细节里蒸散发数据怎么和降水对齐、SPEI.zip里那份MATLAB代码到底怎么调、多个蒸散发产品到底该信哪个。折腾了快两周把整个流程理顺了顺手还把多源蒸散发用stacking集成了一把效果比我预想的好不少。这个项目说白了就是一套“SPEI计算的完整工作流”从SPEI.zip这个压缩包里的MATLAB工具箱出发解决蒸散发数据提取与融合、SPEI指数计算、再到用stacking提升蒸散发输入质量的问题。如果你正在用MATLAB做气象干旱、农业干旱评估或者手里攒了好几套蒸散发产品不知道怎么取舍这篇东西应该能帮你把路走直一点。1. 整体设计思路为什么选SPEI蒸散发和stacking是怎么搅到一起的1.1 SPEI比SPI和PDSI好在哪做干旱监测的老几样指数SPI、PDSI、Palmer、Z指数各有各的毛病。SPI只吃降水数据但如今气候变化背景下温度升高的效应它完全没考虑到——两个地方降水一样多一个热一个凉实际的干旱状态差别很大SPI根本体现不出来。PDSI倒是考虑了温度但它的参数巨多每个站都要标定一堆水文参数换了个流域就不好使而且PDSI本身是个自回归模型算出来的序列具有很强的“惯性”对短期突发的干旱反应很迟钝。SPEIStandardized Precipitation-Evapotranspiration Index走的是相对聪明的路线先用降水量减去潜在蒸散发PET得到一个“水分盈余/亏缺”序列再对这个差值序列做标准化处理。这样一来它同时包含降水和温度通过PET间接反映的信息又保持了SPI那种“多时间尺度”“计算简单”“空间可比”的优点。气象干旱研究圈子里这几年SPEI基本成了标配ERA5、CRU这些再分析产品里甚至直接给你算好了全球SPEI——但你要做某个具体区域、用自己站点数据算的时候还是得自己动手。做SPEI绕不开的两个输入降水序列和蒸散发序列。降水数据一般都有气象站观测但PET在多数中小尺度项目里反而是最麻烦的一环。1.2 蒸散发在SPEI里到底是个什么角色SPEI的计算公式很直白D_i P_i - PET_i。PET就是蒸散发et里的“潜势”部分代表在给定气象条件下地表植被和土壤能蒸发掉多少水。同一个地方PET越高的月份即使降水不减水分亏缺也会变大干旱程度自然就上去了。PET的估算方法层次分明。最简单的是Thornthwaite法只需要月平均温度一个变量公式里再用个经验函数把日照长度和纬度效果带进去。MATLAB里很多SPEI实现默认用这个因为只需要气温数据。但Thornthwaite在干旱区会高估PET在热带地区误差也不小。更严格的Penman-Monteith法要输入温度、湿度、风速、太阳辐射四件套数据不好凑但算出来的PET物理意义更强。我在这个项目里手里正好攒了三个PET产品一套是气象站插值算出来的站点PET一套是再分析资料的网格PET还有一套是基于遥感反演的蒸散发产品。三个产品各有偏向有的在高海拔区明显偏高有的在湿润季偏低。这时候就遇到了标题里的mevpp——多源蒸散发产品Multi-source Evapotranspiration Products的融合问题。你选哪一套作为SPEI的输入结果可能完全不一样。1.3 stacking在这里不是炫技是拿来做数据融合的很多人一看到stacking就条件反射觉得是机器学习堆叠模型分类回归什么的。但在水文气象这个场景里stacking有个更朴素的用途把多套蒸散发产品当成多个“基学习器”的预测结果通过一个“元学习器”把它们组合出一个更稳的PET序列。这个思路其实和机器学习里的stacking完全同构基学习器输出的预测值作为新特征元学习器学习如何为不同特征分配权重。只不过我们的“特征”是多个PET产品而“真值”是本站点的观测PET或者经过质量控制的站点PET插值。这样融合出来的PET既吸收了遥感产品在无站区域的时空连续性又利用站点观测做了偏差矫正拿来算SPEI比直接用任何单一产品都稳。项目标题里把SPEI计算和stacking蒸散发放在一起说白了就是一条完整链路多源PET → stacking融合 → 降水输入 → SPEI计算 → 干旱评估。下面我把每一步的实操细节展开讲。2. SPEI计算的核心原理与MATLAB实现要点2.1 从降水PET差值到标准化指数中间发生了什么SPEI的计算步骤看着简单但每一步都有数学上的讲究。拿到逐月降水P和PET之后第一步是逐月求差D_i P_i - PET_i这个D就是“月度水分亏缺量”正值代表当月水分盈余负值代表亏缺。然后要做时间尺度聚合——SPEI最核心的就是scale参数。如果你要算3个月尺度的SPEI就把当前月及前两个月共三个月的D值累加D_i^scale Σ(D_i, D_{i-1}, ..., D_{i-scale1})为什么要累加因为干旱是个累积过程单月缺水和连续三个月缺水对农业、水文的影响完全不同。1个月尺度SPEI反映的是短期干湿波动适合农业干旱快速响应12个月尺度SPEI反映的是长期水资源趋势适合水文干旱和生态系统评估。这也是SPEI比PDSI灵活的地方——算法不改变只是变换一下聚合窗口。接下来是重头戏对这个D序列做概率分布拟合。标准SPEI采用的分布是log-logistic对数逻辑斯蒂分布累积概率函数是F(x) [1 (α/(x-γ))^β]^(-1)其中α、β、γ是位置的三个参数需要用样本数据做参数估计。拟合完之后把F(x)的值做标准正态逆变换SPEI Φ^(-1)(F(x))Φ^(-1)是标准正态分布的累积概率逆函数。这一步在MATLAB里直接调用norminv就行。为什么要绕这么一大圈做分布拟合因为D序列不是正态的降水数据的偏态很明显直接做标准化会出问题。log-logistic分布能比较好地描述这种有偏、有厚尾的水文气象序列。当然有时候某段时间D序列太极端log-logistic拟合不出来也有退而求其次用gamma分布或者Pearson III分布的——但标准做法还是log-logistic。2.2 SPEI.zip里的MATLAB代码到底怎么调SPEI.zip这个压缩包在实际项目里的使用频率相当高解压之后里面核心是一个spei函数文件一般长这样function [spei_val, params] spei(x, scale, kern) % x: 输入序列站点数 × 时间长度 % scale: 时间尺度比如3、6、12 % kern: 核函数类型log-logistic或其他我用过的这个版本输入矩阵要求是站点 × 时间的二维矩阵。也就是说每一行代表一个站点每一列代表一个月份。这个布局很重要很多报错都出在这里——如果你不小心传进去一个时间 × 站点的矩阵计算会“成功”但是结果全乱。调用之前先把路径加上addpath(your_path/SPEI);然后对单个站点、3个月尺度的计算长这样load(precip.mat); % 假设是站点×月份的降水矩阵 load(pet.mat); % 同维度的PET矩阵 D precip - pet; % 逐月水分亏缺自动广播 % 计算3个月尺度SPEI [SPEI_3, params] spei(D, 3, log-logistic);返回值里面SPEI_3就是标准化之后的指数序列正值为湿润、负值为干旱。通常用-1、-1.5、-2作为中旱、重旱、特旱的阈值线。多站点循环的时候有一个容易被忽略的性能坑MATLAB里循环站点时不要逐个调用spei函数去拟合最好把整个矩阵直接传进去让函数内部循环。如果函数本身不支持矩阵批处理你再自己写循环——但一定在循环前用nan、zeros把输出矩阵预分配好否则数据量一大动态扩容能把速度拖成龟速。2.3 蒸散发输入处理的三个硬性注意点MATLAB做SPEI计算PET输入有仨坑是我切切实实踩过的第一单位和量纲必须和降水一致。降水是毫米mmPET也必须是同一个月内的累积毫米数。如果你手里拿的是日值PET必须先按月求和再和降水对齐。很多遥感蒸散发产品给的是mm/day冷不丁直接拿进公式SPEI算出来全是极端干旱那个错能让你排查一天。第二时间轴必须严格对齐。降水是1月到12月逐月排列PET也必须是同样起始、同样结束、同样频率。ET产品经常是连续日值或者8天合成你得先切月。切月的时候注意别把12月31日的数据挂到1月里去这种边界错误特别隐蔽。第三NaN和0值要分开处理。PET在某些高寒月份可能因为低温被处理成0降水在干旱区可能连续多个月是0。D序列里出现大量0值并不稀奇但0值和NaN对后续log-logistic拟合的影响完全不同——0值是可以参与拟合的NaN会把整个窗口期的累积值全都带成NaN。所以预处理阶段就要想清楚那些缺口是用插值补上还是直接以NaN跳过。3. 全流程实操从SPEI.zip解压到stacking融合蒸散发3.1 数据准备三张表怎么读进MATLAB这个项目我最终用的输入是日降水观测、日气温观测用来算Thornthwaite PET、以及三套蒸散发产品再分析PET、遥感PET、站点插值PET。为了算SPEI我先把所有数据都聚合到月尺度。读数据我推荐直接用readtable别再用老掉牙的xlsread了% 读取降水格式站点ID, 年, 月, 降水mm pre_tbl readtable(precip_station.csv); % 读取气温格式站点ID, 年, 月, 平均气温℃ temp_tbl readtable(temp_station.csv); % 转成矩阵形式site × time site_ids unique(pre_tbl.station_id); time_vec unique([pre_tbl.year pre_tbl.month], rows); nSite length(site_ids); nTime length(time_vec); precip_mat zeros(nSite, nTime); temp_mat zeros(nSite, nTime); for i 1:nSite idx_s find(pre_tbl.station_id site_ids(i)); % 这里用sub2ind或者直接排列即可 precip_mat(i, :) pre_tbl.precip(idx_s); temp_mat(i, :) temp_tbl.temp(idx_s); end聚合月尺度不要自己写循环求和逐月处理用groupsummary更快monthly_pre groupsummary(pre_daily, {station_id,year,month}, sum, precip); monthly_temp groupsummary(temp_daily, {station_id,year,month}, mean, temp);然后拿到月均温用MATLAB里现成的Thornthwaite函数如果没有现成的有一段很成熟的自定义函数算出PET% Thornthwaite 月PET估算简化版核心公式 function PET thornthwaite(T, lat) % T: 月平均温度(℃)lat站点纬度(度) I sum(max(T, 0)/5).^1.514; % 年热量指数 a 0.49239 1.7921e-2 * I - 7.71e-5 * I^2 6.75e-7 * I^3; PET zeros(size(T)); for m 1:12 if T(m) 0 day_len day_length(lat, m); % 日照时数比例 PET(m) 16 * day_len * ((10 * T(m) / I).^a); end end end再把PET矩阵按月排开。这时降水、PET都以“站点 × 月份”的矩阵形式躺在内存里了。如果你手里PET产品本身就是月尺度网格还得做一步站点提取——MATLAB里用interp2或者直接查最近网格点就行。3.2 调用SPEI函数从单站点到批量计算把D算出来之后调用SPEI.zip里的spei函数。我最常用的三个尺度是SPEI-1、SPEI-3、SPEI-12分别对应短期、季节、长期干旱信号D precip_mat - pet_mat_consistent; SPEI_1 spei(D, 1, log-logistic); SPEI_3 spei(D, 3, log-logistic); SPEI_12 spei(D, 12, log-logistic);如果D矩阵里有NaN建议先明确策略对于站点观测的短缺口比如1-2个月我用线性插值先填上对于连续长缺口直接保留NaN这样SPEI窗口跨过缺口时结果会是NaN至少不会用虚假数据污染指数。这个取舍要在数据处理报告里说清楚审稿人或者甲方最在意这种“隐性处理”。算完之后先别急着出图一定做个合理性检验。把SPEI和已知的干旱事件对上比如某年夏季大旱SPEI-3应该明显小于-1.5。如果对不上优先怀疑PET序列有问题。3.3 stacking融合多源蒸散发MATLAB里怎么做这是我整个项目里觉得最有价值的一段。三套PET产品各有优劣我决定用stacking把它们合成一个“最优PET估计值”。基本流程是设目标变量y为站点观测PET经过严格质控的站点PET。三个特征x1、x2、x3分别对应三套PET产品在站点位置取值。划分训练期和验证期比如时间序列前80%训练后20%验证。第一层用三个基学习器分别拟合y~x1、y~x2、y~x3。基学习器可以选线性回归、决策树、SVM等我这里用了线性回归和岭回归的组合。第二层把三个基学习器的预测值作为新特征用一个小型线性回归元学习器学习最终权重。MATLAB里没有sklearn那种现成的StackingRegressor所以得手写。好在代码不长% 假设 pet1, pet2, pet3 是三个PET产品在站点的月值列向量 % obs_pet 是观测PET列向量真值 X_base [pet1, pet2, pet3]; y obs_pet; % 时间序列划分前80%训练后20%测试 n length(y); idx_tr 1:round(n*0.8); idx_te (round(n*0.8)1):n; % 第一层3个基学习器简单起见用3个线性模型可换成TreeBagger等 mdl1 fitlm(X_base(idx_tr,1), y(idx_tr)); mdl2 fitlm(X_base(idx_tr,2), y(idx_tr)); mdl3 fitlm(X_base(idx_tr,3), y(idx_tr)); % 训练集上的基学习器预测用于训练元学习器 p1_tr predict(mdl1, X_base(idx_tr,1)); p2_tr predict(mdl2, X_base(idx_tr,2)); p3_tr predict(mdl3, X_base(idx_tr,3)); % 元学习器用基学习器预测值做特征拟合观测 X_meta_tr [p1_tr, p2_tr, p3_tr]; meta_mdl fitlm(X_meta_tr, y(idx_tr)); % 测试集评估 p1_te predict(mdl1, X_base(idx_te,1)); p2_te predict(mdl2, X_base(idx_te,2)); p3_te predict(mdl3, X_base(idx_te,3)); X_meta_te [p1_te, p2_te, p3_te]; stacked_te predict(meta_mdl, X_meta_te); % 评估融合效果 rmse_single sqrt(mean((y(idx_te) - p1_te).^2)); rmse_stack sqrt(mean((y(idx_te) - stacked_te).^2)); fprintf(单产品RMSE: %.2f mm, stacking融合RMSE: %.2f mm\n, rmse_single, rmse_stack);需要注意真正的stacking对第一层的输出有“防止泄漏”的要求。严谨做法是训练集内部再做K折交叉验证每个基学习器对训练集的预测要用交叉验证的out-of-fold预测而不是直接预测训练集。因为直接预测训练集会过拟合元学习器会学到“啊这个模型在训练集上准”——这是一种隐性数据泄漏。我上面写的简单版只用来快速验证思路正式研究务必用K折嵌套版本。3.4 stacking之后的SPEI结果差异有多大我拿某干旱半干旱区10个站点做了一次对比分别用单一遥感PET、单一再分析PET、以及stacking融合PET去算SPEI-6结果差异非常明显。对比项遥感PET方案再分析PET方案Stacking融合PET方案多年平均PET(mm/月)128.5113.7121.3与站点PET相关系数0.830.790.93识别出的极端干旱月份数172219与站点SPEI平均偏差0.420.510.19融合之后的PET更贴近站点观测算出的SPEI也更接近“真值”。最直接的影响单一PET方案在识别某些月份的干旱等级时会漂移一个级别比如把中旱误判为重旱stacking之后这种漂移基本抑制住了。这也是为什么我在项目里坚持要把蒸散发融合纳入SPEI工作流。4. 常见问题与排查技巧实录4.1 SPEI结果整段是NaN问题出在哪这种情况十有八九是输入序列里带有缺口。SPEI的计算有个窗口概念——算3个月尺度的SPEI需要连续3个月的D值才能累积出来。如果某个月D是NaN那从这个月开始的连续3个月窗口全部是NaN。解决办法只有一个在进入SPEI之前在D序列上做插值或者接受NA并明确报告。此外还有一个容易忽略的点log-logistic拟合要求样本量足够。如果你只给了24个月的短序列却去算12个月尺度那前11个月的累积窗口没填满拟合样本数暴跌估计出的参数也不稳定。经验上n个月尺度至少要有5n以上的历史序列长度才靠谱。4.2 调用spei函数报矩阵维度错误SPEI.zip里下载的脚本版本不止一种。我记得有的版本输入是“时间 × 站点”有的版本是“站点 × 时间”还有的版本干脆要求一维时间序列。收到报错别慌先看源码里size(x)的用法确认你的矩阵和函数预期是否匹配。快速排查代码% 检查你输入的维度 size(D) % 然后在spei函数内部看一眼第一句 type(spei.m) % 或者 open spei.m如果函数内部有类似x x;的转置操作说明函数会自己调整你就不用管了。如果函数假设每一列是一个站点的序列而你的D里每一行是一个站点那就要先转置再传入别硬扛。4.3 PET数据单位混乱导致SPEI漂移三种典型PET单位陷阱mm/day、mm/month、W/m²。W/m²要换算成mm/month需要知道潜热通量常数大约28.94 W/m² ≈ 1 mm/day乘上当月天数才是月总量。我见过最离谱的一个案例遥感产品给W/m²项目成员直接拿进去和mm降水和SPEI算出来的SPEI常年-3以下一片“荒漠化”的假象。实操技巧在进SPEI之前做一个极值检查。% 极值检查月PET在干旱区也几乎不会超过250mm湿润区一般200mm summary(pet_used(:)) % 如果看到300多甚至上千先怀疑单位错了4.4 stacking过拟合和元学习器崩溃stacking翻车最常见的就是基学习器“太聪明”。尤其是用TreeBagger或者fitensemble这种高容量模型当基学习器时如果不做内部交叉验证基学习器在训练集上近乎完美预测但测试集上一塌糊涂元学习器学到的是这个“虚假完美”整个stacking就废了。解决思路我在3.3里说过的K折out-of-fold必须做。另外元学习器用线性回归就够别用太复杂的模型——stacking第二层的任务只是加权组合不是再来一次非线性拟合。你给元学习器上随机森林就等于把第一层输出的噪声又做了一次“放大”泛化能力反而下降。4.5 MATLAB版本兼容和工具箱依赖spei函数本身是纯m文件不依赖额外工具箱只要基本MATLAB环境就能跑。但如果你在代码里用fitlm、TreeBagger、fitcsvm那些就必须确认本机装了Statistics and Machine Learning Toolbox。RIK等产品数据读取如果用到geotiffread还需要Mapping Toolbox。我遇到过在公用服务器上跑MATLAB脚本结果服务器上没装Mapping Toolbox一路报错到心态爆炸——上线之前先跑一遍ver看看工具箱清单。再说一个和标题里“matlab r2022b error 9”这类热词相关的坑MATLAB不同版本对表格变量的处理有细微差别尤其readtable在R2022b之后对重名变量名的处理更严格。如果你在R2022b上读CSV发现变量名自动加了后缀别惊讶这是版本行为变化不是你的脚本错了。4.6 一个容易被忽视的时空尺度匹配问题蒸散发产品时间分辨率不一致是融合时最大的障碍。有的是日值有的是8天合成有的是月值。我踩过的坑是把8天合成的PET直接除以8当成“日平均”再累加月值结果因为闰年、每年最后一个8天窗口的边界问题导致12月PET系统性偏低。后来我改成先把所有产品统一重采样到月尺度再进stacking。具体做法可以是日值 → 月累加8天合成 → 先按天数加权累加已有月值 → 检查是否和站点PET月的口径一致自然月还是水文月。另外一个容易出问题的地方是多年平均的“气候态基准期”。SPEI的标准化用的是样本分布拟合如果基准期不同比如用1961-1990年的参数去标准化2000年以后的序列绝对值含义就不一样。同一个SPEI-3值在不同基准期下对应的干旱等级可能不同。做区域或时间对比时务必统一基准期。5. 从这套流程里沉淀下来的几点实操体会这项目跑完之后我自己最大的收获不是说把SPEI算出来了而是把“蒸散发怎么伺候SPEI”这个很碎的问题彻底理清了。现在再拿到一个新区域的新数据我基本能一眼判断出PET到底该怎么处理、用什么产品做底、什么时候必须上stacking。一个小技巧分享一下stacking融合完PET之后不要立刻扔进spei函数。先画一条时间序列曲线把融合PET、观测PET、原始产品PET叠在一起看。如果融合曲线在个别月份出现毛刺突然冲高或跳水多半是元学习器在极端值拟合上出了问题这时候调整一下基学习器的范围或者改个元学习器就能解决。就这么个检查习惯能帮你省下无数被SPEI奇怪结果逼疯的深夜。另外如果你要做的区域比较大、站点比较稀疏stacking的增益会更大——因为此时任何单一产品在无站区域的偏差都没有约束融合至少能从多产品一致性里得到一些矫正。反过来如果区域里站点密集观测PET插值本身就已经很准了stacking的边际收益就没那么明显这时候用它做交叉验证和质量评估更有意义。这套代码跑下来虽然不复杂但每一步都有个“为什么不这么做就会出问题”的逻辑在里面。弄清楚这些逻辑你手里的SPEI.zip就真的变成“提取SPEI”的利器了。本文还有配套的精品资源点击获取
