企业官网建设流程全解析

简介：这份资源面向气象水文、生态地理及农业干旱研究方向的科研人员与研究生，提供标准化降水蒸散发指数（SPEI）的MATLAB计算实现，用于量化区域干湿变化与干旱程度。压缩包共21个文件，约16.25MB，以13个m脚本为核心，辅以4份pdf文献、2个txt示例数据、1个xlsx结果表及1个exe程序，覆盖计算、绘图与批量处理等环节。使用者按示例数据格式运行SPEI_Cal.m即可得到结果，并在stacking中查询，时间尺度可自行调整，亲测可用。包内还包含基于两种潜在蒸散发算法的SPEI分析文献，以及AddData、SPEI_Batch、mkfig等辅助脚本，便于理解蒸散发输入、批量计算与图形输出流程。目前已有1622人学习下载，适合希望快速复现SPEI计算、对照文献验证结果并开展干旱分析的研究者参考。

1. 从 SPEI.zip 说起：一份 MATLAB 蒸散发与 SPEI 计算链路到底解决了什么

如果你手头有一个叫SPEI.zip的压缩包，里面是 MATLAB 脚本，还夹着mevpp、stacking、蒸散发这些关键词，那它大概率不是单一功能的玩具脚本，而是一条从气象栅格数据到干旱指数的完整处理链。SPEI（标准化降水蒸散指数）本身是干旱监测里最常用的指标之一，但真正让一线从业者头疼的从来不是公式，而是蒸散发怎么算、时间尺度怎么滑、栅格怎么对齐、缺测怎么补。这套东西要解决的就是把降水、气温、蒸散发三类输入，经过潜在蒸散发估算、水分盈亏累积、概率分布拟合、标准化转换，最终输出一张能直接进 GIS 或论文的 SPEI 栅格。适合谁？做区域干旱评估的研究生、做农业气象业务的技术员、以及手里有站点或栅格气象数据但卡在蒸散发环节的工程师。下面我按自己复现这类项目的顺序，把每一步拆开讲。

2. 蒸散发与 SPEI 的底层逻辑：为什么 mevpp 和 stacking 会出现在同一个包里

2.1 SPEI 计算链里蒸散发为什么不能随便替代

SPEI 的核心是水分盈亏D = P - PET，其中 PET 是潜在蒸散发。很多人图省事直接用 Thornthwaite 法，因为它只需要月均温和纬度，计算极快。但 Thornthwaite 在干旱区和高海拔区偏差很大，因为它对温度的非线性响应过于敏感，忽略了辐射和风速。如果你手里有mevpp这个关键词，它通常指向 MODIS 蒸散发产品（MOD16）或者基于 Penman-Monteith 的改进算法。用mevpp类数据替代 Thornthwaite，SPEI 的干湿趋势会发生系统性偏移，尤其在夏季峰值月份。我一般会先确认输入蒸散发的时空分辨率是否和降水一致，再决定是直接重采样还是重新计算 PET。这一步选错，后面所有标准化都是白做。

2.2 stacking 在气象栅格处理里到底 stack 什么

stacking在这个语境下不是机器学习里的模型堆叠，而是把多时相栅格按时间维堆成一个三维数组或一个长表。MATLAB 里常见做法是用cat(3, ...)或者把每个月的栅格读成矩阵后存进 cell，再转成double数组。为什么要 stack？因为 SPEI 需要至少 12 个月以上的滑动窗口，你不可能每个月单独算完再拼，必须在内存里保持一个lat × lon × time的结构。stacking 的另一个作用是统一 NaN 掩膜：如果某个月某个像元缺测，整个时间序列在该像元上都要标记，否则分布拟合会报错。我见过太多人在这里翻车，后面拟合出来全是 NaN 还找不到原因。

2.3 从 D 序列到 SPEI：分布拟合的选型与参数

算出 D 序列后，SPEI 的标准做法是用 log-logistic 分布拟合，然后做概率转换再取标准正态分位数。MATLAB 里可以用fitdist配合自定义 log-logistic，或者直接用ksdensity做非参数估计。参数上最关键的是滑动窗口长度：3 个月适合农业干旱，12 个月适合水文干旱，48 个月适合地下水位响应。窗口越长，序列两端缺失越多，通常要牺牲前 11 个月。我一般会先跑 12 个月窗口，检查拟合的 p 值是否大于 0.05，如果大量像元不通过，就换非参数方法或者对 D 做 Box-Cox 变换。这一步没有后悔药，拟合不通过硬转标准正态，出来的 SPEI 会失真。

3. 用 MATLAB 跑通 SPEI 最小链路：从栅格读取到标准化输出

3.1 环境准备与数据目录约定

MATLAB 版本建议 R2020b 以上，因为fitdist和geotiffread的稳定性在旧版上差异很大。如果你用的是 R2023b 或更新，注意中文注释乱码问题，可以在预设里把编码改成 UTF-8，否则脚本里的中文变量名会变成乱码导致报错。数据目录我习惯这样组织：data/precip/放降水栅格，data/temp/放气温，data/pet/放蒸散发，output/放结果。每个文件夹里按YYYYMM.tif命名。下面是最小读取和 stacking 的代码。

% 读取指定年份范围的降水、气温、PET 栅格并 stack 成三维数组 years = 2000:2020; months = 1:12; [lat, lon] = deal([]); P = []; T = []; PET = []; for y = years for m = months fn_p = sprintf('data/precip/%04d%02d.tif', y, m); fn_t = sprintf('data/temp/%04d%02d.tif', y, m); fn_e = sprintf('data/pet/%04d%02d.tif', y, m); if ~exist(fn_p,'file') || ~exist(fn_t,'file') || ~exist(fn_e,'file') warning('缺失 %04d-%02d，跳过', y, m); continue; end [p, R] = readgeoraster(fn_p); [t, ~] = readgeoraster(fn_t); [e, ~] = readgeoraster(fn_e); if isempty(P) [lat, lon] = deal(R.LatitudeLimits, R.LongitudeLimits); end P = cat(3, P, double(p)); T = cat(3, T, double(t)); PET = cat(3, PET, double(e)); end end fprintf('stack 完成：%d 个时相\n', size(P,3));

这段代码的逻辑很直接：逐月读栅格，用readgeoraster拿到数据和地理参考，然后沿第三维拼接。参数上注意double转换必须在cat之前做，否则整型拼接会溢出。R.LatitudeLimits只取一次，用来后面写输出文件。如果某个文件缺失，这里选择跳过并警告，但实际业务中更稳妥的做法是记录缺失月份，后面在时间序列里插入 NaN，而不是直接跳过导致时间轴错位。

3.2 计算水分盈亏 D 并处理缺测像元

拿到 P 和 PET 后，D = P - PET。但直接相减之前要检查两者的空间范围是否完全一致。如果 PET 来自 MODIS，分辨率可能是 500m，而降水是 0.05 度，必须先重采样到同一网格。我一般用imresize对 PET 做双线性插值到降水网格，或者用interp2做经纬度插值。下面假设已经对齐。

% 计算水分盈亏 D，并标记全时间序列缺测的像元 D = P - PET; % 检查每个像元的时间序列有效长度 valid_count = sum(~isnan(D), 3); min_valid = 12; % 至少 12 个月有效才参与拟合 mask_bad = valid_count < min_valid; % 对有效像元，用线性插值补内部短缺口 for i = 1:size(D,1) for j = 1:size(D,2) if mask_bad(i,j), continue; end ts = squeeze(D(i,j,:)); nan_idx = isnan(ts); if any(nan_idx) && sum(~nan_idx) >= min_valid ts(nan_idx) = interp1(find(~nan_idx), ts(~nan_idx), find(nan_idx), 'linear', 'extrap'); D(i,j,:) = reshape(ts, 1, 1, []); end end end fprintf('有效像元比例：%.2f%%\n', 100 * (1 - sum(mask_bad(:))/numel(mask_bad)));

这里的关键参数是min_valid，我设 12 是因为后续 12 个月滑动窗口至少需要 12 个点。interp1的extrap选项要慎用，如果缺口在序列两端，外推可能产生虚假趋势，更安全的做法是两端缺口直接保留 NaN，只补中间。循环写法在 MATLAB 里效率不高，如果栅格很大，建议用arrayfun或把时间维移到第一维做向量化，但为了可读性，小区域这样写没问题。

3.3 滑动累积与 log-logistic 拟合

SPEI 的滑动窗口通常用movsum对 D 做累积，窗口长度k取 3、12、48。累积后再对每个像元做分布拟合。MATLAB 没有内置 log-logistic，需要自定义概率密度和累积分布函数。下面给出一个可用的拟合与转换流程。

% 对每个像元做 12 个月滑动累积并拟合 log-logistic，输出 SPEI k = 12; [nlat, nlon, nt] = size(D); SPEI = nan(nlat, nlon, nt - k + 1); % 自定义 log-logistic 的 CDF loglogcdf = @(x, alpha, beta) 1 ./ (1 + (alpha ./ x).^beta); for i = 1:nlat for j = 1:nlon if mask_bad(i,j), continue; end ts = squeeze(D(i,j,:)); if sum(~isnan(ts)) < k, continue; end acc = movsum(ts, [k-1 0], 'omitnan'); acc = acc(k:end); % 去掉前 k-1 个不完整窗口 if numel(acc) < 30, continue; end % 样本太少不拟合 % 用矩估计初值，再用 mle 拟合 mu = mean(acc, 'omitnan'); sd = std(acc, 'omitnan'); beta0 = pi / (sd * sqrt(3)); alpha0 = mu; try params = mle(acc, 'pdf', @(x,a,b) (b/a).*(x/a).^(b-1)./(1+(x/a).^b).^2, ... 'start', [alpha0, beta0], 'LowerBound', [eps, eps]); cdf_vals = loglogcdf(acc, params(1), params(2)); cdf_vals = min(max(cdf_vals, 1e-6), 1-1e-6); spei_ts = norminv(cdf_vals, 0, 1); SPEI(i, j, :) = reshape(spei_ts, 1, 1, []); catch continue; end end end fprintf('SPEI 计算完成，有效时间步：%d\n', size(SPEI,3));

这段代码里movsum的[k-1 0]表示向后累积，即当前月加上前 11 个月。mle拟合 log-logistic 时，start给矩估计初值能显著提高收敛率。norminv把概率转成标准正态分位数，就是 SPEI。注意cdf_vals的截断，避免 0 和 1 导致无穷大。如果某个像元拟合失败，这里直接跳过，输出保持 NaN，后面出图时用掩膜盖掉。参数k可以根据研究尺度改，但改完要重新检查有效样本数，30 是我个人经验下限，低于这个数拟合结果不稳定。

4. 避坑与排查：SPEI 计算里最容易翻车的五个地方

4.1 现象：SPEI 全图都是 NaN，但 D 序列明明有值

原因通常是mask_bad把大部分像元标记为无效，而mask_bad的判断依据是valid_count < min_valid。如果降水或 PET 栅格里有大量 0 值被误当成有效，或者栅格范围没对齐导致D里全是 NaN，就会触发。解决方法是先单独检查P和PET的sum(~isnan(...),3)，确认两者有效像元重叠区域。另一个常见原因是movsum的omitnan在窗口内全 NaN 时返回 NaN，导致acc长度不够。

4.2 现象：拟合报错 “The data must be non-negative” 或类似

log-logistic 分布定义在正数域，但 D = P - PET 在干旱区经常为负。直接对负值拟合会报错。常见做法是对 D 做平移，比如D_shift = D - min(D) + eps，拟合后再把 CDF 对应回去。但平移会改变分布形状，更严谨的做法是用三参数 log-logistic，增加位置参数。MATLAB 的mle支持自定义三参数 PDF，但收敛更慢。我一般先试平移，如果拟合 p 值可接受就用，否则换非参数ksdensity做经验 CDF。

4.3 现象：SPEI 时间序列在窗口边界出现剧烈跳变

这是因为滑动累积的前 k-1 个时间步被丢弃后，剩下的序列起点和终点样本数不对称。尤其是 48 个月窗口，丢掉前 47 个月后，如果总时间长度只有 60 个月，拟合样本只有 13 个，分布参数极不稳定。解决方法是保证总时间长度至少是窗口的 3 倍以上，或者改用更短的窗口。另一个技巧是对累积序列做去趋势，避免气候变化信号被当成干旱信号。

4.4 现象：输出栅格在 GIS 里打开后上下颠倒

MATLAB 的栅格读取默认第一行对应最北，但某些 GeoTIFF 的R参考对象里LatitudeLimits是递增的。写输出时如果直接用geotiffwrite而不翻转，就会颠倒。解决方法是写之前用flipud翻转矩阵，或者用maprefcells明确指定参考。我习惯在readgeoraster后立刻检查R.LatitudeLimits(1) > R.LatitudeLimits(2)，如果不是就翻转。

4.5 现象：MATLAB 版本升级后mle对自定义 PDF 的接口变了

R2023a 之后mle对函数句柄的参数检查更严格，如果 PDF 里用了未定义的变量或者参数顺序不对，会直接报错而不是警告。解决方法是把自定义 PDF 写成独立函数文件，而不是匿名函数，并且在函数里加assert检查输入。另外，R2023b 的中文注释乱码问题也会导致函数名解析失败，建议脚本里变量名和函数名全用英文，注释可以用中文但保存为 UTF-8。

5. 进阶技巧：用 stacking 思路做多模型 SPEI 集合与验证

如果你已经不满足于单一 PET 输入的 SPEI，可以走集合路线：用 Thornthwaite、Penman-Monteith、mevpp 三种 PET 各算一套 SPEI，然后按像元做中位数或加权平均。这就是stacking在更高层的用法——不是堆栅格，而是堆模型输出。具体做法是把三套 SPEI 沿第四维拼接，得到lat × lon × time × model，然后用median(..., 4)得到集合 SPEI。验证时可以用站点 SPEI 做相关分析，或者用历史干旱事件做命中率检验。我一般会保留三套原始结果，因为不同 PET 在湿润区和干旱区的表现差异很大，集合值能平滑掉部分偏差，但也会掩盖极端事件。最后说一个我自己的习惯：每次跑完 SPEI，先画一张 12 个月窗口的 Hovmöller 图，看时间维有没有断裂，再画空间分布图看有没有条带伪影。这两个图能拦住八成以上的低级错误。希望帮到你。

本文还有配套的精品资源，点击获取