matlab 编写SPEI计算代码

SPEI（Standardized Precipitation Evapotranspiration Index）是一种基于标准化降水和蒸散发的气候指数，用于评估地区干旱或湿润程度。以下是MATLAB编写SPEI计算代码的示例：

1. 首先，加载数据。假设你已经有了一个包含每月降水和蒸散发数据的矩阵，其中每一行代表一个时间步长，每一列代表一个站点。你可以使用MATLAB的load函数加载数据文件。

data = load('datafile.mat');
precip = data.precip; % 降水数据
evap = data.evap; % 蒸散发数据

2. 计算SPEI。下面是一个计算SPEI的函数，其中alpha是形状参数（通常为1），beta是比例参数，gamma是位置参数，delta是尺度参数，n是时间步长数量。函数返回一个SPEI矩阵，其中每一行代表一个时间步长，每一列代表一个站点。

function spei = calculate_spei(precip, evap, alpha, beta, gamma, delta, n)
% 计算SPEI
% precip: 降水数据矩阵，每行代表一个时间步长，每列代表一个站点
% evap: 蒸散发数据矩阵，每行代表一个时间步长，每列代表一个站点
% alpha, beta, gamma, delta: SPEI参数
% n: 时间步长数量
% spei: SPEI矩阵，每行代表一个时间步长，每列代表一个站点

% 计算PET（潜在蒸散发）
pet = zeros(size(evap));
for i = 1:n
    pet(i,:) = evap(i,:) .* (1 + beta .* (evap(i,:) ./ precip(i,:)));
end

% 计算P-PET
ppet = precip - pet;

% 计算累积P-PET
c = zeros(size(ppet));
for i = 2:n
    c(i,:) = (1 - alpha) .* c(i-1,:) + alpha .* ppet(i-1,:);
end

% 计算SPEI
spei = zeros(size(ppet));
for i = 1:n
    spei(i,:) = (ppet(i,:) - gamma - delta .* c(i,:)) ./ sqrt(delta .* sum((c(i,:) - mean(c(i,:))) .^ 2) ./ (n-1));
end
end

3. 调用函数并输出SPEI结果。

alpha = 1;
beta = 0.5;
gamma = 0;
delta = 1;
n = size(precip, 1);

spei = calculate_spei(precip, evap, alpha, beta, gamma, delta, n);

% 输出结果
disp(spei);