简介本资源是一套基于MATLAB实现的标准化降水蒸散发指数SPEI计算工具包面向气象水文、生态遥感及干旱监测领域的科研人员与研究生解决多时间尺度下区域干湿变化定量评估问题。压缩包共21个文件含13个核心MATLAB脚本如SPEI_Cal.m、SPEI_Batch.m、AddCell.m等支撑数据预处理、批量计算与结果可视化、4份PDF文献含刘珂团队关于两种潜在蒸散发算法的SPEI对比分析及英文手册、2个说明文本、1个Excel结果模板和1个独立exe可执行程序整体大小为16.25MB。已有1614人学习下载表明其在实际科研中具备较高复用价值。用户可直接按示例数据格式运行主程序灵活修改时间尺度如1/3/12个月结果自动输出至stacking子目录配套PDF文献提供算法原理与验证案例m文件模块化设计便于二次开发特别适合开展中国区域干旱时空演变研究或课程实验实践。 做干旱监测的人大概率都绕不开SPEI这个指数尤其是用MATLAB处理气象站点数据的时候。标题里的“SPEI.zip_SPEI matlab_matlab提取SPEI_mevpp_stacking_蒸散发”一眼看过去像是从某个下载页面复制过来的代码包名加搜索关键词但拆开看其实是一条很完整的技术链路用MATLAB从气象数据里提取SPEI中间涉及潜在蒸散发的估算而这个估算环节用到了stacking集成学习mevpp可以理解为一个蒸散发预测模块或者说脚本。这篇博文就把这条链路完整讲清楚从SPEI的计算原理到MATLAB里怎么实现再到为什么要用stacking来搞蒸散发最后附上我在实际数据上跑通的流程和踩过的坑。这篇内容适合谁正在做干旱监测、水文分析、农业气象相关研究的同学尤其是手里攒了一批站点气象数据、想用MATLAB快速算出SPEI但又被各种工具包和版本问题折磨的人。另外对机器学习感兴趣、想看看stacking怎么落地到气象数据上的朋友也能从中找到一些可以“抄作业”的思路。1. 项目整体设计与思路拆解1.1 先搞清楚这几个关键词到底在说什么先说SPEI全称Standardized Precipitation Evapotranspiration Index标准化降水蒸散发指数。它和SPI标准化降水指数最大的区别就是SPEI在降水的基础上把“蒸散发”这个消耗项也纳入了水分盈亏的计算。说白了一个地方干不干旱不光看下了多少雨还要看这些水够不够蒸发、够不够植物蒸腾用。SPEI这个指数就是基于“降水减去潜在蒸散发PET”的累积值来计算的兼顾了降水对干旱的补给和温度升高带来的蒸发加剧效应。这个指数在两个场景下特别好用。一是气候变化背景下的干旱趋势分析因为温度升高会让PET变大即使降水不变SPEI也会反映出干旱加剧的趋势而SPI是反映不出来的。二是不同时间尺度的干旱监测SPEI可以像SPI一样支持1个月、3个月、6个月、12个月、24个月等不同累积尺度用来对应不同类型的干旱。再说mevpp和stacking。从项目的命名习惯来看mevpp应该是一个和蒸散发预测相关的模块或函数集合可能全称是类似Model-based Evapotranspiration Prediction Package这样的东西。在实际项目里这一部分通常承担的任务是在蒸散发站点观测数据不足、或者需要把单站PET扩展到区域尺度的时候用气象要素作为输入特征通过机器学习模型来估算PET或实际蒸散发ET。而stacking就是这一步用的核心方法——一种集成学习策略把多个基础模型的预测结果作为新的特征再训练一个元模型来综合决策。之所以用stacking而不是随便选一个单模型后面我会详细说。1.2 为什么是MATLAB而不是Python这个可能是很多人纠结的问题。我的看法是如果你做SPEI计算、做气象数据的批处理、做结果的可视化出图MATLAB在某些环节真的更顺手。尤其体现在几个地方数据接口方便。站点气象数据很多以Excel、CSV、NetCDF格式存储MATLAB的readtable、readmatrix、ncread等函数处理起来很直接不需要像Python那样频繁地pandas、numpy来回转换。科研出图质量高。MATLAB的绘图语法虽然老但胜在稳定出图细节可控很多气象、水文学术论文的图还是MATLAB风格居多。SPEI相关的现成代码资源多。标题里的SPEI.zip就是一个典型网上有不少科研团队或个人分享的MATLAB版SPEI计算包下载下来调一调就能用比从零开始写省事得多。当然Python在机器学习的生态上确实更丰富scikit-learn、XGBoost、LightGBM这些库用起来非常方便。所以在实际项目中我推荐的混合方案是数据整理和SPEI计算用MATLAB蒸散发的机器学习预测部分用Python完成如果你更习惯MATLAB后面我也会给出一个MATLAB手写stacking的替代方案最终再把预测的PET结果拿回MATLAB里参与SPEI计算。这种混合工作流在科研中非常常见不要有“非此即彼”的包袱。1.3 整体技术链路把整个项目拆开可以梳理成下面这条链路站点气象数据准备降水、气温最高/最低/平均、湿度、风速、日照时数等。蒸散发估算先用传统方法如Thornthwaite或观测数据得到初步PET然后用stacking模型拟合PET与气象因子之间的关系用来做插补或区域扩展。SPEI计算基于降水和PET数据计算逐月水分盈亏按时间尺度累积拟合概率分布标准化得到SPEI值。结果分析与可视化干旱等级划分、趋势分析、时间序列图输出。这条链路的关键点在于第2步和第3步的衔接。很多人在算SPEI的时候直接把PET用Thornthwaite公式一带而过但对于数据质量差、站点稀疏的研究区这一步做得粗糙会直接影响SPEI的可靠性。用stacking模型来优化蒸散发估算算是这个项目里比较有亮点的地方。2. SPEI的计算原理与MATLAB实现2.1 SPEI的计算流程拆成五步SPEI的计算原理不复杂但每一步都有讲究。标准流程可以拆成五步第一步计算逐月潜在蒸散发PET。常用方法有两种。Thornthwaite方法只需要月平均气温和站点纬度计算简单适合数据资料匮乏的地区也是SPEI原始论文作者推荐的方法。Penman-Monteith方法需要最高温、最低温、湿度、风速、日照时数等数据精度更高但对数据要求也高。实际项目中如果你手头资料齐全建议用Penman-Monteith如果只有气温数据Thornthwaite也是一个能接受的选择。第二步计算逐月水分盈亏。公式很直白D_i P_i - PET_i其中P_i是第i个月的降水量PET_i是第i个月的潜在蒸散发量D_i就是当月的水分盈亏。D为正值说明水分有盈余负值说明水分亏缺。第三步在不同时间尺度上累积。比如计算3个月尺度的SPEI就是把当前月及前两个月的水分盈亏加起来D_k^3 sum(D_i, ik-2 to k)这个累加过程可以平滑短期波动反映中期水分状况。12个月尺度则反映年际尺度的干湿变化。第四步对累积序列拟合概率分布。SPEI默认使用log-logistic分布也称为Fisk分布来拟合累积水分盈亏序列。这一步是整个计算中最核心也最容易出问题的地方数据量太少或者极端值太多都可能导致拟合失败。第五步标准化转换。将拟合得到的累积概率值转换为标准正态分布的Z值也就是最终的SPEI值。转换后SPEI就具有了统一的可比性不同地区、不同时段的SPEI值可以直接对比。2.2 MATLAB里的SPEI.zip到底怎么用网上流传的SPEI.zip包里面通常包含几样东西一个主函数可能是spei.m、一个计算PET的函数pet.m或thornthwaite.m、一个log-logistic分布拟合的子函数以及一份示例数据。我第一次拿到这个包的时候也折腾了一阵子这里把核心用法直接说明白。假设你的数据结构是第一列年份第二列月份第三列降水单位mm第四列气温单位摄氏度那么核心调用逻辑如下% 读取数据 data readmatrix(example_station.csv); year data(:,1); month data(:,2); precip data(:,3); temp data(:,4); % 计算PET这里以Thornthwaite方法为例 % lat是站点纬度度需要根据实际情况填写 lat 36.5; pet thornthwaite(temp, lat); % 计算SPEI时间尺度设为12个月 scale 12; % spei函数的输入一般是逐月降水、逐月PET、时间尺度 % 返回的spi_out就是SPEI值序列 spei_value spei(precip, pet, scale);有些版本的spei函数还需要你传入数据起始年份或者告诉你输出结果中哪些位置是有效的因为前面的数据在累积过程中会被截断前scale-1个月通常没有有效值。用的时候要仔细看一眼代码注释别直接把整个序列拿去用。2.3 分布拟合和标准化的细节log-logistic分布拟合这一步很多MATLAB实现用的是L矩法L-moments来估计分布参数而不是最大似然法。原因是气象水文序列容易出现厚尾分布L矩法对异常值更稳健估计出来的参数更稳定。标准化的过程是对数化后的累积概率P按照下面的规则转换当P大于0.5时取W sqrt(-2 * log(1 - P))然后计算Z值SPEI。当P小于等于0.5时取W sqrt(-2 * log(P))用另一套系数计算Z值。这一步在MATLAB里通常都有现成的函数实现不需要自己从头写但你要理解这个过程不然边界值出现异常时不知道怎么排查。2.4 SPEI的分级标准算出SPEI之后怎么判断干旱等级学术界一般参照下面的划分方法SPEI值等级 2.0极端湿润1.5 ~ 2.0严重湿润1.0 ~ 1.5中等湿润-1.0 ~ 1.0正常-1.5 ~ -1.0中等干旱-2.0 ~ -1.5严重干旱 -2.0极端干旱这个分级和SPI基本一致方便对比使用。3. 蒸散发预测stacking集成学习的落地3.1 为什么蒸散发需要机器学习来预测先回答一个很多人会问的问题蒸散发不是有公式可以算吗为什么还要用机器学习原因有三个。第一Penman-Monteith公式虽然精度高但需要的输入要素太多。很多气象站只有降水和气温记录湿度、风速、辐射这些要素要么缺失要么记录年代不连续。这种情况下直接用简化公式会带来较大误差不如用机器学习把数据中的非线性关系挖掘出来。第二站点不够密的时候需要做空间插值或者区域扩展机器学习可以把地形、纬度、海拔这些因子一起纳入比简单的反距离权重插值靠谱得多。第三实际蒸散发ET和潜在蒸散发PET之间存在复杂的关系受植被、土壤水分的影响这种关系很难用单一物理公式刻画而机器学习可以捕捉到这些隐含规律。标题里的mevpp从我理解的角度来看就是负责这一部分工作的模块。它做的事情本质上是一个回归任务用气象特征预测蒸散发值或者说建立气象因子到蒸散发之间的映射模型。3.2 特征工程哪些因子要进模型我用过的特征组合给各位一个参考基本覆盖了主要的气象驱动因子气温类平均气温、最高气温、最低气温、气温日较差。降水类当月降水、前1~3个月降水滞后效应。湿度类平均相对湿度、最小相对湿度。风类平均风速、最大风速。辐射类日照时数作为辐射的代理变量、太阳总辐射如有。时间类月份、季节编码用于捕捉季节性周期。地理类纬度、海拔、坡度、坡向做空间预测时用。这些特征在使用前要做标准化或归一化处理特别是对SVR这类对特征尺度敏感的模型。树模型随机森林、XGBoost对尺度不敏感但统一标准化是通用做法避免某些特征主导。3.3 模型的stacking框架怎么搭Stacking的核心思想是多模型各显神通让一个元模型学习如何“综合”它们的意见。用一句话概括第一层是基学习器第二层是元学习器。具体到蒸散发预测这个任务我的做法是选出三个差异明显的基学习器。比如随机森林Random Forest擅长捕捉非线性关系、对异常值稳健XGBoost极端梯度提升在中等规模数据上往往精度很高SVR支持向量回归擅长高维小样本。三者差异越大stacking的效果越好。如果三个模型结构类似比如都是深层树模型集成效果会大打折扣这一点要特别注意。用交叉验证生成元特征。这是stacking最关键的一步也是最容易翻车的一步。为了防止数据泄漏不能用全体训练集直接训练基学习器然后预测训练集来生成第二层特征那样元模型会把“记忆”当“学习”。正确做法是用K折交叉验证把训练集分成5折每次用4折训练、1折预测轮转5次后每个训练样本都得到了一个“诚实”的预测值把这些预测值拼起来作为元特征。测试集上的预测则通过5个模型的平均获得。元学习器选择简单的线性模型。第二层模型不需要太复杂一般用岭回归Ridge Regression或者带L2正则的线性回归。因为基学习器的输出已经包含了复杂的非线性信息元模型只需要学习一个最优加权组合即可。用过于复杂的元模型反而容易过拟合。评价指标。蒸散发预测常用的评价指标有决定系数R²、均方根误差RMSE、平均绝对误差MAE。实际项目中R²在0.7以上基本可以接受用于区域估算0.85以上算很好。3.4 MATLAB和Python怎么配合这里说实话Python的scikit-learn提供了现成的StackingRegressor用起来极其方便代码量少不少。写出来大概是这样from sklearn.ensemble import RandomForestRegressor, StackingRegressor from sklearn.svm import SVR from sklearn.linear_model import Ridge from xgboost import XGBRegressor from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler import pandas as pd # 读取特征和目标 data pd.read_csv(et_features.csv) X data.drop([PET], axis1) y data[PET] # 标准化 scaler StandardScaler() X_scaled scaler.fit_transform(X) # 定义基学习器和元学习器 base_models [ (rf, RandomForestRegressor(n_estimators500, random_state42)), (xgb, XGBRegressor(n_estimators500, learning_rate0.05, random_state42)), (svr, SVR(kernelrbf, C10, epsilon0.1)) ] meta_model Ridge(alpha1.0) # 构建stacking模型 stack_model StackingRegressor(estimatorsbase_models, final_estimatormeta_model, cv5) # 训练和评估 X_train, X_test, y_train, y_test train_test_split(X_scaled, y, test_size0.2, random_state42) stack_model.fit(X_train, y_train) print(R^2:, stack_model.score(X_test, y_test))核心就这些scikit-learn的封装已经把交叉验证生成元特征的细节全部处理好了你不用自己写循环。如果你坚持在MATLAB里完成整个流程也不是不行。MATLAB的Statistics and Machine Learning Toolbox里提供了fitrtree回归树、fitrensemble集成学习包括随机森林、fitrsvmSVR等函数。可以手动实现stacking用fitrtree、fitrensemble、fitrsvm训练三个基模型。手动写一个K折交叉验证循环每折都对训练数据预测一次拼出元特征矩阵。用fitrlinear或ridge回归作为元模型输入元特征矩阵进行训练。代码量会大一些但可行。我之前在MATLAB R2021a上试过几百个站点的数据规模跑起来完全没问题。4. 实操过程与核心环节实现4.1 数据准备与清洗决定成败的前置环节做气象数据分析的人都有体会数据清洗花的时间往往比建模还多。降水和气温数据里常见的问题包括缺测值、异常值比如降水出现负数、气温超过物理极限、站点迁址导致的数据不连续、还有时间戳错乱等。我的建议是在进入SPEI计算之前先把数据质量检查做了检查4个方面缺失率检查逐月数据中缺测比例超过10%的站点要么插补要么直接剔除。极值检查降水0且500mm单月极端暴雨除外根据研究区气候调整气温在-50℃~50℃之间。连续缺测检查连续缺测超过3个月的插补难度大谨慎处理。时间连续性检查月份是否连续有没有跳月、重复月份。缺失数据的插补我常用的方法是先用临近站点的数据做线性回归插补回归不了的就用多年同月平均值代替。后者精度差一些但对SPEI计算的影响相对可控因为SPEI本身就是基于累积值的指数对单点误差有一定的平滑作用。下面是一段MATLAB的数据清洗示例% 读取原始数据 data readmatrix(raw_station.csv); year data(:,1); month data(:,2); precip data(:,3); temp data(:,4); % 标记缺测值假设缺测用-9999表示 precip(precip -9000) NaN; temp(temp -9000) NaN; % 逐月插补用多年平均代替缺测 [uy, ~, iy] unique(year); [um, ~, im] unique(month); for i 1:length(precip) if isnan(precip(i)) month_idx month(i); precip(i) nanmean(precip(month month_idx)); end if isnan(temp(i)) month_idx month(i); temp(i) nanmean(temp(month month_idx)); end end % 检查插补结果 fprintf(降水缺测插补后剩余NaN数量%d\n, sum(isnan(precip))); fprintf(气温缺测插补后剩余NaN数量%d\n, sum(isnan(temp)));4.2 完整流程跑一遍从原始数据到SPEI序列以西北某研究区的一个站点为例手头有1970-2020年逐月降水和气温数据我们一步步算。第一步计算PET。使用Thornthwaite方法只需要月均温和纬度。MATLAB里可以这样写function pet thornthwaite(temp, lat) % temp: 逐月平均气温单位摄氏度 % lat: 站点纬度单位度 % 输出pet: 逐月潜在蒸散发单位mm % Thornthwaite方法需要先用月均温计算年热指数I T max(temp, 0); % 低于0度的月份不参与计算 I sum((T / 5).^1.514); % 年热指数 % 计算系数a a 6.75e-7 * I^3 - 7.71e-5 * I^2 1.79e-2 * I 0.49; % 逐月计算未校正的PET pet_unadj 16 * (10 * T / I).^a; pet_unadj(T 0) 0; % 根据纬度和月份计算日长校正系数 month_days [31 28 31 30 31 30 31 31 30 31 30 31]; for m 1:12 % 日均日照时数N用纬度近似 N daily_sunshine_hours(lat, m); % 每月天数修正 days_correction month_days(m) / 30; % 每天12小时标准日照修正 sunshine_correction N / 12; % 校正系数 K(m) days_correction * sunshine_correction; end % 应用到每个月份 months mod((1:length(temp))-1, 12) 1; pet pet_unadj .* K(months); end这里daily_sunshine_hours是一个计算逐月日均日照时数的辅助函数根据纬度和月份用天文公式计算。实际项目中可以直接用经验公式不必追求非常精确。第二步计算水分盈亏序列D然后按12个月尺度累积。D precip - pet; scale 12; n length(D); D_cum zeros(n, 1); for i scale:n D_cum(i) sum(D(i-scale1:i)); end D_cum(1:scale-1) NaN; % 前scale-1个值无有效累积第三步对累积序列拟合log-logistic分布并标准化。这一步建议直接使用现成的spei函数核心逻辑是% 提取有效累积值进行分布拟合 valid_D D_cum(~isnan(D_cum)); % 用L矩法估计log-logistic分布参数均值、尺度、形状 % 具体实现略这里是示意 [mu_s, alpha_s, xi_s] fit_loglogistic_Lmom(valid_D); % 计算累积概率 % 当D_cum不为0时F(x) 1 / (1 (alpha_s / (x - gamma_s))^beta_s) % 再标准化转换得到SPEI spei_out zeros(n, 1); spei_out(:) NaN; for i scale:n x D_cum(i); F 1 / (1 (alpha_s / (x - xi_s))^mu_s); % 概率转标准正态Z值 if F 0.5 W sqrt(-2 * log(F)); spei_out(i) W - (2.515517 0.802853*W 0.010328*W^2) / ... (1 1.432788*W 0.189269*W^2 0.001308*W^3); else W sqrt(-2 * log(1 - F)); spei_out(i) -(W - (2.515517 0.802853*W 0.010328*W^2) / ... (1 1.432788*W 0.189269*W^2 0.001308*W^3)); end end上面的系数是Abramowitz和Stegun给出的标准正态分位数近似公式系数在很多SPEI实现里都能看到。第四步出图把SPEI时间序列画出来。% 构造日期向量 dates datetime(year, month, 1); % 绘图 figure(Color, w); plot(dates, spei_out, b-, LineWidth, 1.2); hold on; yline(0, k--); hold off; xlabel(时间); ylabel(SPEI-12); title(某站点12个月尺度SPEI时间序列); grid on; set(gca, FontName, Helvetica, FontSize, 12);4.3 参数选择与调优经验之谈在实际项目中有几个参数需要你根据研究区特点做调整不要照搬默认值。时间尺度scale的选择。如果你想分析农业干旱与土壤墒情和作物生长关系密切选1~3个月尺度更合适如果你想研究水文干旱与径流、水库蓄水关系密切选6~12个月如果你想看长期气候变化趋势下的干旱演变12~24个月尺度更稳。我在实际项目中一般会同时算SPEI-1、SPEI-3、SPEI-12分别对应不同应用场景。Thornthwaite方法的纬度参数。这个方法对纬度比较敏感纬度输入错误会直接影响PET的计算。如果站点多建议整理一个站点信息表统一读取纬度而不是手动输入。stacking模型方面基学习器的超参数可以先粗调再细调。随机森林主要调n_estimators和max_depthXGBoost主要调learning_rate和n_estimatorsSVR主要调C和epsilon。不要一上来就网格搜索先手动跑一两组感受一下数据规模和时间成本再做精细调参。5. 常见问题与排查技巧实录5.1 数据缺失导致的累积序列间断SPEI计算过程中如果原始数据中间缺了几个月累积序列就会断掉。有些版本的SPEI代码会直接把含缺测的累积值设为NaN导致结果序列上有空洞。如果你发现SPEI结果中间有突变的空值先检查是不是原始降水或气温数据有缺测而不是怀疑代码逻辑。解决办法在计算累积值之前先把缺测值插补完整再进入累积流程。如果你不想插补又想保留原始数据特征可以考虑用“缺测月份不参与累积只对有效月份累积”的改进方案但这样不同时间点的累积窗口长度不一致对比时需要谨慎。5.2 log-logistic分布拟合失败这是最容易出问题的地方。当你发现SPEI结果全部为NaN或者分布拟合时报错通常有三个原因一是数据长度太短。一般建议至少30年以上的逐月数据SPEI计算结果才稳定。如果只有10年数据累积序列样本量不够L矩法估计的参数方差很大。二是数据中存在大量相同值比如干旱区很多月份降水为0导致累积序列的分布出现极端偏态拟合算法不收敛。三是序列中有离谱的异常值。比如某个极端暴雨月份降水是正常月份的10倍这个值会严重拉偏拟合参数。5.3 边界值溢出SPEI出现无穷大标准化转换时如果累积概率F非常接近0或1W值会变得很大SPEI可能溢出为无穷大。很多SPEI代码里会有截断处理把SPEI限制在-3.09到3.09之间。如果没有截断你的结果图里可能会出现突兀的尖峰。我一般会检查一下最大值最小值如果超过了合理范围手动加上截断。5.4 用实测数据检验结果SPEI算完之后怎么看合不合理我用三个方法检验与历史干旱记录对比。把历次明显干旱年份的SPEI值调出来看是否与实际情况吻合。与SPI对比。在降水主导的研究区SPEI和SPI的序列走势应该大体一致但SPEI可能会反映出更多的温度信号。看趋势。多年SPEI序列做线性回归或Mann-Kendall趋势检验看研究区是干化还是湿化这个趋势与已有的文献报道对比一下。6. 几个踩坑后的个人体会最后说几点我自己的经验。不要在数据清洗上省时间。我踩过最大的坑就是拿到一个“看起来完整”的站点数据直接跑了SPEI结果发现中间有两年的降水数据被静默替换成了-9999导致SPEI序列出现一个巨大的假干旱峰值。从那以后所有数据进来先做缺测和极值扫描再进主流程。SPEI.zip这类代码包拿来之后第一件事是打开源码读一遍搞清楚输入数据格式和每个参数的含义。很多论文作者写的代码注释并不完善变量命名也很随意直接跑了报错再回去看代码反而更浪费时间。我现在拿到任何开源MATLAB代码第一件事是建一个demo脚本用包自带的示例数据跑通一遍再替换成自己的数据。stacking模型的收益不是必然的。有时候单个XGBoost已经表现很好加了一个表现较差的SVR进来反而拖后腿。我在项目里的做法是先单独评估每个基学习器的性能只保留那些交叉验证R²还不错的模型进入stacking宁可基学习器少而精也不要多而杂。另外蒸散发预测结果和SPEI计算之间有一个关键细节预测出来的PET是估计值存在不确定性。如果直接把估计值与实测降水一起算SPEI会把模型的误差传导到SPEI结果里。稳妥的做法是对PET预测值做一次不确定性分析比如用多组stacking模型的重复训练结果取区间然后考察SPEI在PET波动下的稳定性。如果SPEI结果对PET非常敏感那么你要么提升模型精度要么在结果解释时留有余地。最后一个建议也是这些年做干旱数据分析最大的体会不要迷信任何一个指数也不要迷信任何一个模型。SPEI再科学它也只是一个概况性的指标只能告诉你某个地区某个时段相对于其历史状况的干湿偏离程度但干旱的成因、影响、应对需要结合当地的水文、农业、社会经济情况综合判断。工具是死的应用是活的搞清楚每个参数背后的物理意义和每个模型假设的适用边界比单纯跑出一个漂亮的结果重要得多。本文还有配套的精品资源点击获取
MATLAB实现SPEI干旱指数：从蒸散发预测到Stacking集成学习完整指南
简介本资源是一套基于MATLAB实现的标准化降水蒸散发指数SPEI计算工具包面向气象水文、生态遥感及干旱监测领域的科研人员与研究生解决多时间尺度下区域干湿变化定量评估问题。压缩包共21个文件含13个核心MATLAB脚本如SPEI_Cal.m、SPEI_Batch.m、AddCell.m等支撑数据预处理、批量计算与结果可视化、4份PDF文献含刘珂团队关于两种潜在蒸散发算法的SPEI对比分析及英文手册、2个说明文本、1个Excel结果模板和1个独立exe可执行程序整体大小为16.25MB。已有1614人学习下载表明其在实际科研中具备较高复用价值。用户可直接按示例数据格式运行主程序灵活修改时间尺度如1/3/12个月结果自动输出至stacking子目录配套PDF文献提供算法原理与验证案例m文件模块化设计便于二次开发特别适合开展中国区域干旱时空演变研究或课程实验实践。 做干旱监测的人大概率都绕不开SPEI这个指数尤其是用MATLAB处理气象站点数据的时候。标题里的“SPEI.zip_SPEI matlab_matlab提取SPEI_mevpp_stacking_蒸散发”一眼看过去像是从某个下载页面复制过来的代码包名加搜索关键词但拆开看其实是一条很完整的技术链路用MATLAB从气象数据里提取SPEI中间涉及潜在蒸散发的估算而这个估算环节用到了stacking集成学习mevpp可以理解为一个蒸散发预测模块或者说脚本。这篇博文就把这条链路完整讲清楚从SPEI的计算原理到MATLAB里怎么实现再到为什么要用stacking来搞蒸散发最后附上我在实际数据上跑通的流程和踩过的坑。这篇内容适合谁正在做干旱监测、水文分析、农业气象相关研究的同学尤其是手里攒了一批站点气象数据、想用MATLAB快速算出SPEI但又被各种工具包和版本问题折磨的人。另外对机器学习感兴趣、想看看stacking怎么落地到气象数据上的朋友也能从中找到一些可以“抄作业”的思路。1. 项目整体设计与思路拆解1.1 先搞清楚这几个关键词到底在说什么先说SPEI全称Standardized Precipitation Evapotranspiration Index标准化降水蒸散发指数。它和SPI标准化降水指数最大的区别就是SPEI在降水的基础上把“蒸散发”这个消耗项也纳入了水分盈亏的计算。说白了一个地方干不干旱不光看下了多少雨还要看这些水够不够蒸发、够不够植物蒸腾用。SPEI这个指数就是基于“降水减去潜在蒸散发PET”的累积值来计算的兼顾了降水对干旱的补给和温度升高带来的蒸发加剧效应。这个指数在两个场景下特别好用。一是气候变化背景下的干旱趋势分析因为温度升高会让PET变大即使降水不变SPEI也会反映出干旱加剧的趋势而SPI是反映不出来的。二是不同时间尺度的干旱监测SPEI可以像SPI一样支持1个月、3个月、6个月、12个月、24个月等不同累积尺度用来对应不同类型的干旱。再说mevpp和stacking。从项目的命名习惯来看mevpp应该是一个和蒸散发预测相关的模块或函数集合可能全称是类似Model-based Evapotranspiration Prediction Package这样的东西。在实际项目里这一部分通常承担的任务是在蒸散发站点观测数据不足、或者需要把单站PET扩展到区域尺度的时候用气象要素作为输入特征通过机器学习模型来估算PET或实际蒸散发ET。而stacking就是这一步用的核心方法——一种集成学习策略把多个基础模型的预测结果作为新的特征再训练一个元模型来综合决策。之所以用stacking而不是随便选一个单模型后面我会详细说。1.2 为什么是MATLAB而不是Python这个可能是很多人纠结的问题。我的看法是如果你做SPEI计算、做气象数据的批处理、做结果的可视化出图MATLAB在某些环节真的更顺手。尤其体现在几个地方数据接口方便。站点气象数据很多以Excel、CSV、NetCDF格式存储MATLAB的readtable、readmatrix、ncread等函数处理起来很直接不需要像Python那样频繁地pandas、numpy来回转换。科研出图质量高。MATLAB的绘图语法虽然老但胜在稳定出图细节可控很多气象、水文学术论文的图还是MATLAB风格居多。SPEI相关的现成代码资源多。标题里的SPEI.zip就是一个典型网上有不少科研团队或个人分享的MATLAB版SPEI计算包下载下来调一调就能用比从零开始写省事得多。当然Python在机器学习的生态上确实更丰富scikit-learn、XGBoost、LightGBM这些库用起来非常方便。所以在实际项目中我推荐的混合方案是数据整理和SPEI计算用MATLAB蒸散发的机器学习预测部分用Python完成如果你更习惯MATLAB后面我也会给出一个MATLAB手写stacking的替代方案最终再把预测的PET结果拿回MATLAB里参与SPEI计算。这种混合工作流在科研中非常常见不要有“非此即彼”的包袱。1.3 整体技术链路把整个项目拆开可以梳理成下面这条链路站点气象数据准备降水、气温最高/最低/平均、湿度、风速、日照时数等。蒸散发估算先用传统方法如Thornthwaite或观测数据得到初步PET然后用stacking模型拟合PET与气象因子之间的关系用来做插补或区域扩展。SPEI计算基于降水和PET数据计算逐月水分盈亏按时间尺度累积拟合概率分布标准化得到SPEI值。结果分析与可视化干旱等级划分、趋势分析、时间序列图输出。这条链路的关键点在于第2步和第3步的衔接。很多人在算SPEI的时候直接把PET用Thornthwaite公式一带而过但对于数据质量差、站点稀疏的研究区这一步做得粗糙会直接影响SPEI的可靠性。用stacking模型来优化蒸散发估算算是这个项目里比较有亮点的地方。2. SPEI的计算原理与MATLAB实现2.1 SPEI的计算流程拆成五步SPEI的计算原理不复杂但每一步都有讲究。标准流程可以拆成五步第一步计算逐月潜在蒸散发PET。常用方法有两种。Thornthwaite方法只需要月平均气温和站点纬度计算简单适合数据资料匮乏的地区也是SPEI原始论文作者推荐的方法。Penman-Monteith方法需要最高温、最低温、湿度、风速、日照时数等数据精度更高但对数据要求也高。实际项目中如果你手头资料齐全建议用Penman-Monteith如果只有气温数据Thornthwaite也是一个能接受的选择。第二步计算逐月水分盈亏。公式很直白D_i P_i - PET_i其中P_i是第i个月的降水量PET_i是第i个月的潜在蒸散发量D_i就是当月的水分盈亏。D为正值说明水分有盈余负值说明水分亏缺。第三步在不同时间尺度上累积。比如计算3个月尺度的SPEI就是把当前月及前两个月的水分盈亏加起来D_k^3 sum(D_i, ik-2 to k)这个累加过程可以平滑短期波动反映中期水分状况。12个月尺度则反映年际尺度的干湿变化。第四步对累积序列拟合概率分布。SPEI默认使用log-logistic分布也称为Fisk分布来拟合累积水分盈亏序列。这一步是整个计算中最核心也最容易出问题的地方数据量太少或者极端值太多都可能导致拟合失败。第五步标准化转换。将拟合得到的累积概率值转换为标准正态分布的Z值也就是最终的SPEI值。转换后SPEI就具有了统一的可比性不同地区、不同时段的SPEI值可以直接对比。2.2 MATLAB里的SPEI.zip到底怎么用网上流传的SPEI.zip包里面通常包含几样东西一个主函数可能是spei.m、一个计算PET的函数pet.m或thornthwaite.m、一个log-logistic分布拟合的子函数以及一份示例数据。我第一次拿到这个包的时候也折腾了一阵子这里把核心用法直接说明白。假设你的数据结构是第一列年份第二列月份第三列降水单位mm第四列气温单位摄氏度那么核心调用逻辑如下% 读取数据 data readmatrix(example_station.csv); year data(:,1); month data(:,2); precip data(:,3); temp data(:,4); % 计算PET这里以Thornthwaite方法为例 % lat是站点纬度度需要根据实际情况填写 lat 36.5; pet thornthwaite(temp, lat); % 计算SPEI时间尺度设为12个月 scale 12; % spei函数的输入一般是逐月降水、逐月PET、时间尺度 % 返回的spi_out就是SPEI值序列 spei_value spei(precip, pet, scale);有些版本的spei函数还需要你传入数据起始年份或者告诉你输出结果中哪些位置是有效的因为前面的数据在累积过程中会被截断前scale-1个月通常没有有效值。用的时候要仔细看一眼代码注释别直接把整个序列拿去用。2.3 分布拟合和标准化的细节log-logistic分布拟合这一步很多MATLAB实现用的是L矩法L-moments来估计分布参数而不是最大似然法。原因是气象水文序列容易出现厚尾分布L矩法对异常值更稳健估计出来的参数更稳定。标准化的过程是对数化后的累积概率P按照下面的规则转换当P大于0.5时取W sqrt(-2 * log(1 - P))然后计算Z值SPEI。当P小于等于0.5时取W sqrt(-2 * log(P))用另一套系数计算Z值。这一步在MATLAB里通常都有现成的函数实现不需要自己从头写但你要理解这个过程不然边界值出现异常时不知道怎么排查。2.4 SPEI的分级标准算出SPEI之后怎么判断干旱等级学术界一般参照下面的划分方法SPEI值等级 2.0极端湿润1.5 ~ 2.0严重湿润1.0 ~ 1.5中等湿润-1.0 ~ 1.0正常-1.5 ~ -1.0中等干旱-2.0 ~ -1.5严重干旱 -2.0极端干旱这个分级和SPI基本一致方便对比使用。3. 蒸散发预测stacking集成学习的落地3.1 为什么蒸散发需要机器学习来预测先回答一个很多人会问的问题蒸散发不是有公式可以算吗为什么还要用机器学习原因有三个。第一Penman-Monteith公式虽然精度高但需要的输入要素太多。很多气象站只有降水和气温记录湿度、风速、辐射这些要素要么缺失要么记录年代不连续。这种情况下直接用简化公式会带来较大误差不如用机器学习把数据中的非线性关系挖掘出来。第二站点不够密的时候需要做空间插值或者区域扩展机器学习可以把地形、纬度、海拔这些因子一起纳入比简单的反距离权重插值靠谱得多。第三实际蒸散发ET和潜在蒸散发PET之间存在复杂的关系受植被、土壤水分的影响这种关系很难用单一物理公式刻画而机器学习可以捕捉到这些隐含规律。标题里的mevpp从我理解的角度来看就是负责这一部分工作的模块。它做的事情本质上是一个回归任务用气象特征预测蒸散发值或者说建立气象因子到蒸散发之间的映射模型。3.2 特征工程哪些因子要进模型我用过的特征组合给各位一个参考基本覆盖了主要的气象驱动因子气温类平均气温、最高气温、最低气温、气温日较差。降水类当月降水、前1~3个月降水滞后效应。湿度类平均相对湿度、最小相对湿度。风类平均风速、最大风速。辐射类日照时数作为辐射的代理变量、太阳总辐射如有。时间类月份、季节编码用于捕捉季节性周期。地理类纬度、海拔、坡度、坡向做空间预测时用。这些特征在使用前要做标准化或归一化处理特别是对SVR这类对特征尺度敏感的模型。树模型随机森林、XGBoost对尺度不敏感但统一标准化是通用做法避免某些特征主导。3.3 模型的stacking框架怎么搭Stacking的核心思想是多模型各显神通让一个元模型学习如何“综合”它们的意见。用一句话概括第一层是基学习器第二层是元学习器。具体到蒸散发预测这个任务我的做法是选出三个差异明显的基学习器。比如随机森林Random Forest擅长捕捉非线性关系、对异常值稳健XGBoost极端梯度提升在中等规模数据上往往精度很高SVR支持向量回归擅长高维小样本。三者差异越大stacking的效果越好。如果三个模型结构类似比如都是深层树模型集成效果会大打折扣这一点要特别注意。用交叉验证生成元特征。这是stacking最关键的一步也是最容易翻车的一步。为了防止数据泄漏不能用全体训练集直接训练基学习器然后预测训练集来生成第二层特征那样元模型会把“记忆”当“学习”。正确做法是用K折交叉验证把训练集分成5折每次用4折训练、1折预测轮转5次后每个训练样本都得到了一个“诚实”的预测值把这些预测值拼起来作为元特征。测试集上的预测则通过5个模型的平均获得。元学习器选择简单的线性模型。第二层模型不需要太复杂一般用岭回归Ridge Regression或者带L2正则的线性回归。因为基学习器的输出已经包含了复杂的非线性信息元模型只需要学习一个最优加权组合即可。用过于复杂的元模型反而容易过拟合。评价指标。蒸散发预测常用的评价指标有决定系数R²、均方根误差RMSE、平均绝对误差MAE。实际项目中R²在0.7以上基本可以接受用于区域估算0.85以上算很好。3.4 MATLAB和Python怎么配合这里说实话Python的scikit-learn提供了现成的StackingRegressor用起来极其方便代码量少不少。写出来大概是这样from sklearn.ensemble import RandomForestRegressor, StackingRegressor from sklearn.svm import SVR from sklearn.linear_model import Ridge from xgboost import XGBRegressor from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler import pandas as pd # 读取特征和目标 data pd.read_csv(et_features.csv) X data.drop([PET], axis1) y data[PET] # 标准化 scaler StandardScaler() X_scaled scaler.fit_transform(X) # 定义基学习器和元学习器 base_models [ (rf, RandomForestRegressor(n_estimators500, random_state42)), (xgb, XGBRegressor(n_estimators500, learning_rate0.05, random_state42)), (svr, SVR(kernelrbf, C10, epsilon0.1)) ] meta_model Ridge(alpha1.0) # 构建stacking模型 stack_model StackingRegressor(estimatorsbase_models, final_estimatormeta_model, cv5) # 训练和评估 X_train, X_test, y_train, y_test train_test_split(X_scaled, y, test_size0.2, random_state42) stack_model.fit(X_train, y_train) print(R^2:, stack_model.score(X_test, y_test))核心就这些scikit-learn的封装已经把交叉验证生成元特征的细节全部处理好了你不用自己写循环。如果你坚持在MATLAB里完成整个流程也不是不行。MATLAB的Statistics and Machine Learning Toolbox里提供了fitrtree回归树、fitrensemble集成学习包括随机森林、fitrsvmSVR等函数。可以手动实现stacking用fitrtree、fitrensemble、fitrsvm训练三个基模型。手动写一个K折交叉验证循环每折都对训练数据预测一次拼出元特征矩阵。用fitrlinear或ridge回归作为元模型输入元特征矩阵进行训练。代码量会大一些但可行。我之前在MATLAB R2021a上试过几百个站点的数据规模跑起来完全没问题。4. 实操过程与核心环节实现4.1 数据准备与清洗决定成败的前置环节做气象数据分析的人都有体会数据清洗花的时间往往比建模还多。降水和气温数据里常见的问题包括缺测值、异常值比如降水出现负数、气温超过物理极限、站点迁址导致的数据不连续、还有时间戳错乱等。我的建议是在进入SPEI计算之前先把数据质量检查做了检查4个方面缺失率检查逐月数据中缺测比例超过10%的站点要么插补要么直接剔除。极值检查降水0且500mm单月极端暴雨除外根据研究区气候调整气温在-50℃~50℃之间。连续缺测检查连续缺测超过3个月的插补难度大谨慎处理。时间连续性检查月份是否连续有没有跳月、重复月份。缺失数据的插补我常用的方法是先用临近站点的数据做线性回归插补回归不了的就用多年同月平均值代替。后者精度差一些但对SPEI计算的影响相对可控因为SPEI本身就是基于累积值的指数对单点误差有一定的平滑作用。下面是一段MATLAB的数据清洗示例% 读取原始数据 data readmatrix(raw_station.csv); year data(:,1); month data(:,2); precip data(:,3); temp data(:,4); % 标记缺测值假设缺测用-9999表示 precip(precip -9000) NaN; temp(temp -9000) NaN; % 逐月插补用多年平均代替缺测 [uy, ~, iy] unique(year); [um, ~, im] unique(month); for i 1:length(precip) if isnan(precip(i)) month_idx month(i); precip(i) nanmean(precip(month month_idx)); end if isnan(temp(i)) month_idx month(i); temp(i) nanmean(temp(month month_idx)); end end % 检查插补结果 fprintf(降水缺测插补后剩余NaN数量%d\n, sum(isnan(precip))); fprintf(气温缺测插补后剩余NaN数量%d\n, sum(isnan(temp)));4.2 完整流程跑一遍从原始数据到SPEI序列以西北某研究区的一个站点为例手头有1970-2020年逐月降水和气温数据我们一步步算。第一步计算PET。使用Thornthwaite方法只需要月均温和纬度。MATLAB里可以这样写function pet thornthwaite(temp, lat) % temp: 逐月平均气温单位摄氏度 % lat: 站点纬度单位度 % 输出pet: 逐月潜在蒸散发单位mm % Thornthwaite方法需要先用月均温计算年热指数I T max(temp, 0); % 低于0度的月份不参与计算 I sum((T / 5).^1.514); % 年热指数 % 计算系数a a 6.75e-7 * I^3 - 7.71e-5 * I^2 1.79e-2 * I 0.49; % 逐月计算未校正的PET pet_unadj 16 * (10 * T / I).^a; pet_unadj(T 0) 0; % 根据纬度和月份计算日长校正系数 month_days [31 28 31 30 31 30 31 31 30 31 30 31]; for m 1:12 % 日均日照时数N用纬度近似 N daily_sunshine_hours(lat, m); % 每月天数修正 days_correction month_days(m) / 30; % 每天12小时标准日照修正 sunshine_correction N / 12; % 校正系数 K(m) days_correction * sunshine_correction; end % 应用到每个月份 months mod((1:length(temp))-1, 12) 1; pet pet_unadj .* K(months); end这里daily_sunshine_hours是一个计算逐月日均日照时数的辅助函数根据纬度和月份用天文公式计算。实际项目中可以直接用经验公式不必追求非常精确。第二步计算水分盈亏序列D然后按12个月尺度累积。D precip - pet; scale 12; n length(D); D_cum zeros(n, 1); for i scale:n D_cum(i) sum(D(i-scale1:i)); end D_cum(1:scale-1) NaN; % 前scale-1个值无有效累积第三步对累积序列拟合log-logistic分布并标准化。这一步建议直接使用现成的spei函数核心逻辑是% 提取有效累积值进行分布拟合 valid_D D_cum(~isnan(D_cum)); % 用L矩法估计log-logistic分布参数均值、尺度、形状 % 具体实现略这里是示意 [mu_s, alpha_s, xi_s] fit_loglogistic_Lmom(valid_D); % 计算累积概率 % 当D_cum不为0时F(x) 1 / (1 (alpha_s / (x - gamma_s))^beta_s) % 再标准化转换得到SPEI spei_out zeros(n, 1); spei_out(:) NaN; for i scale:n x D_cum(i); F 1 / (1 (alpha_s / (x - xi_s))^mu_s); % 概率转标准正态Z值 if F 0.5 W sqrt(-2 * log(F)); spei_out(i) W - (2.515517 0.802853*W 0.010328*W^2) / ... (1 1.432788*W 0.189269*W^2 0.001308*W^3); else W sqrt(-2 * log(1 - F)); spei_out(i) -(W - (2.515517 0.802853*W 0.010328*W^2) / ... (1 1.432788*W 0.189269*W^2 0.001308*W^3)); end end上面的系数是Abramowitz和Stegun给出的标准正态分位数近似公式系数在很多SPEI实现里都能看到。第四步出图把SPEI时间序列画出来。% 构造日期向量 dates datetime(year, month, 1); % 绘图 figure(Color, w); plot(dates, spei_out, b-, LineWidth, 1.2); hold on; yline(0, k--); hold off; xlabel(时间); ylabel(SPEI-12); title(某站点12个月尺度SPEI时间序列); grid on; set(gca, FontName, Helvetica, FontSize, 12);4.3 参数选择与调优经验之谈在实际项目中有几个参数需要你根据研究区特点做调整不要照搬默认值。时间尺度scale的选择。如果你想分析农业干旱与土壤墒情和作物生长关系密切选1~3个月尺度更合适如果你想研究水文干旱与径流、水库蓄水关系密切选6~12个月如果你想看长期气候变化趋势下的干旱演变12~24个月尺度更稳。我在实际项目中一般会同时算SPEI-1、SPEI-3、SPEI-12分别对应不同应用场景。Thornthwaite方法的纬度参数。这个方法对纬度比较敏感纬度输入错误会直接影响PET的计算。如果站点多建议整理一个站点信息表统一读取纬度而不是手动输入。stacking模型方面基学习器的超参数可以先粗调再细调。随机森林主要调n_estimators和max_depthXGBoost主要调learning_rate和n_estimatorsSVR主要调C和epsilon。不要一上来就网格搜索先手动跑一两组感受一下数据规模和时间成本再做精细调参。5. 常见问题与排查技巧实录5.1 数据缺失导致的累积序列间断SPEI计算过程中如果原始数据中间缺了几个月累积序列就会断掉。有些版本的SPEI代码会直接把含缺测的累积值设为NaN导致结果序列上有空洞。如果你发现SPEI结果中间有突变的空值先检查是不是原始降水或气温数据有缺测而不是怀疑代码逻辑。解决办法在计算累积值之前先把缺测值插补完整再进入累积流程。如果你不想插补又想保留原始数据特征可以考虑用“缺测月份不参与累积只对有效月份累积”的改进方案但这样不同时间点的累积窗口长度不一致对比时需要谨慎。5.2 log-logistic分布拟合失败这是最容易出问题的地方。当你发现SPEI结果全部为NaN或者分布拟合时报错通常有三个原因一是数据长度太短。一般建议至少30年以上的逐月数据SPEI计算结果才稳定。如果只有10年数据累积序列样本量不够L矩法估计的参数方差很大。二是数据中存在大量相同值比如干旱区很多月份降水为0导致累积序列的分布出现极端偏态拟合算法不收敛。三是序列中有离谱的异常值。比如某个极端暴雨月份降水是正常月份的10倍这个值会严重拉偏拟合参数。5.3 边界值溢出SPEI出现无穷大标准化转换时如果累积概率F非常接近0或1W值会变得很大SPEI可能溢出为无穷大。很多SPEI代码里会有截断处理把SPEI限制在-3.09到3.09之间。如果没有截断你的结果图里可能会出现突兀的尖峰。我一般会检查一下最大值最小值如果超过了合理范围手动加上截断。5.4 用实测数据检验结果SPEI算完之后怎么看合不合理我用三个方法检验与历史干旱记录对比。把历次明显干旱年份的SPEI值调出来看是否与实际情况吻合。与SPI对比。在降水主导的研究区SPEI和SPI的序列走势应该大体一致但SPEI可能会反映出更多的温度信号。看趋势。多年SPEI序列做线性回归或Mann-Kendall趋势检验看研究区是干化还是湿化这个趋势与已有的文献报道对比一下。6. 几个踩坑后的个人体会最后说几点我自己的经验。不要在数据清洗上省时间。我踩过最大的坑就是拿到一个“看起来完整”的站点数据直接跑了SPEI结果发现中间有两年的降水数据被静默替换成了-9999导致SPEI序列出现一个巨大的假干旱峰值。从那以后所有数据进来先做缺测和极值扫描再进主流程。SPEI.zip这类代码包拿来之后第一件事是打开源码读一遍搞清楚输入数据格式和每个参数的含义。很多论文作者写的代码注释并不完善变量命名也很随意直接跑了报错再回去看代码反而更浪费时间。我现在拿到任何开源MATLAB代码第一件事是建一个demo脚本用包自带的示例数据跑通一遍再替换成自己的数据。stacking模型的收益不是必然的。有时候单个XGBoost已经表现很好加了一个表现较差的SVR进来反而拖后腿。我在项目里的做法是先单独评估每个基学习器的性能只保留那些交叉验证R²还不错的模型进入stacking宁可基学习器少而精也不要多而杂。另外蒸散发预测结果和SPEI计算之间有一个关键细节预测出来的PET是估计值存在不确定性。如果直接把估计值与实测降水一起算SPEI会把模型的误差传导到SPEI结果里。稳妥的做法是对PET预测值做一次不确定性分析比如用多组stacking模型的重复训练结果取区间然后考察SPEI在PET波动下的稳定性。如果SPEI结果对PET非常敏感那么你要么提升模型精度要么在结果解释时留有余地。最后一个建议也是这些年做干旱数据分析最大的体会不要迷信任何一个指数也不要迷信任何一个模型。SPEI再科学它也只是一个概况性的指标只能告诉你某个地区某个时段相对于其历史状况的干湿偏离程度但干旱的成因、影响、应对需要结合当地的水文、农业、社会经济情况综合判断。工具是死的应用是活的搞清楚每个参数背后的物理意义和每个模型假设的适用边界比单纯跑出一个漂亮的结果重要得多。本文还有配套的精品资源点击获取