当前位置：首页 > news >正文

基于MATLAB的标准化降水蒸散指数（SPEI）实现

news 2025/10/21 12:41:40

一、架构

%% 主程序框架
[prec, pet, time] = load_data('input.nc'); % 加载降水与PET数据
prec_acc = accumulate(prec, 3); % 3个月时间尺度累积
pet_acc = accumulate(pet, 3);
d = prec_acc - pet_acc; % 水分盈亏量
spei = calculate_spei(d, 'loglogistic'); % 计算SPEI
plot_spei(spei, time); % 可视化
save_tif(spei, 'spei.tif'); % 输出栅格

二、模块实现

1. 数据读取与预处理

function [prec, pet, time] = load_data(filename)% 支持NetCDF和Excel格式if contains(filename, '.nc')nc = netcdf.open(filename);prec = ncread(nc, 'precipitation');pet = ncread(nc, 'pet');time = ncread(nc, 'time');netcdf.close(nc);elsedata = readtable(filename);prec = data.Precipitation;pet = data.PET;time = datetime(data.Year, data.Month, 1);end% 数据清洗prec = fillmissing(prec, 'linear');pet = fillmissing(pet, 'linear');invalid = (prec < 0) | (pet < 0);prec(invalid) = NaN;pet(invalid) = NaN;
end

2. 潜在蒸散量（PET）计算

function pet = calculate_pet(temp, method)% 支持Thornthwaite和Penman-Monteith方法switch methodcase 'thornthwaite'I = 0.0123 * mean(temp)^3 - 0.245 * mean(temp)^2 + 1.166 * mean(temp) + 0.063;pet = 16 * (10*temp/I).^3;case 'penman_monteith'% 输入参数：temp(温度), rh(湿度), wind(风速), rad(辐射)gamma = 0.000665; % 湿度常数delta = 4098 * exp(17.27*temp./(237.3+temp)) ./ (237.3+temp).^2; % 饱和水汽压斜率es = delta .* (17.27*temp./(237.3+temp)); % 饱和水汽压ea = 0.6108 * exp(17.27*temp./(237.3+temp)); % 实际水汽压Rn = 0.75 * rad; % 净辐射G = 0.1 * temp; % 土壤热通量pet = (0.408*delta*(Rn-G) + gamma*(900/temp+273)*wind*(es-ea)).../ (delta + gamma*(1+0.34*wind));otherwiseerror('未知计算方法');end
end

3. 多时间尺度累积

function acc = accumulate(data, window)% 滑动窗口累积acc = movmean(data, [window-1, 0], 'omitnan', 'Endpoints','discard');% 异常值处理mu = nanmean(acc);sigma = nanstd(acc);acc = (acc - mu) ./ sigma; % 标准化
end

4. SPEI计算核心

function spei = calculate_spei(d, dist)% 参数估计（log-logistic分布）n = length(d);W0 = sum((1./(1+exp(-0.1*(d-(-1.5))))));W1 = sum((1./(1+exp(-0.1*(d-(-0.5))))));W2 = sum((1./(1+exp(-0.1*(d-0.1)))));beta = 2*(W1-W0)/(6*(W1-W0) - 6*W2 + sqrt(36*(W1-W0)^2 - 24*(W1-W0)*W2));alpha = (W0-2*W1)*beta/(gamma(1+1/beta)*gamma(1-1/beta));gamma = W0 - alpha*gamma(1+1/beta)*gamma(1-1/beta);% 标准化计算P = 1./(1+exp(-(d - gamma)/alpha));spei = zeros(size(d));for i = 1:nif P(i) <= 0.5w = -2*log(P(i));spei(i) = w - (2.515517 + 0.802853*w + 0.010328*w^2).../(1 + 1.432788*w + 0.189269*w^2 + 0.001308*w^3);elsew = sqrt(-2*log(1-P(i)));spei(i) = -((2.515517 + 0.802853*w + 0.010328*w^2).../(1 + 1.432788*w + 0.189269*w^2 + 0.001308*w^3) - w);endend
end

三、高级功能

1. 空间可视化（使用m_map工具箱）

function plot_spei(spei, time)% 地理空间绘制m_proj('mercator','long',[73 135],'lat',[18 54]);m_coast('Color','k');m_grid('FontSize',10,'LineStyle',':');% 颜色映射cmap = parula(256);[X,Y] = meshgrid(linspace(73,135,144), linspace(18,54,73));Z = griddata(lon,lat,spei,X,Y,'cubic');% 绘制等值线contourf(X,Y,Z,20,'LineColor','none');hold on;[cs,h] = m_contour(X,Y,Z,[ -2 -1 0 1 2 ],'k');clabel(cs,h);title(sprintf('SPEI指数 (%s)',datestr(time,'yyyy-mm')));colorbar;
end

2. 干旱等级分类

function classes = classify_drought(spei)% 干旱等级划分标准classes = zeros(size(spei));classes(spei < -1.5) = 4; % 极端干旱classes(spei < -1)  = 3; % 严重干旱classes(spei < -0.5) = 2; % 中度干旱classes(spei < 0)    = 1; % 轻度干旱
end

推荐代码 matlab标准化降水蒸散发指数 www.youwenfan.com/contentcsj/53357.html

四、完整工作流

%% 数据准备
filename = 'China_Meteorological_Data.nc';
[prec, pet, time] = load_data(filename);%% PET计算（Penman-Monteith方法）
temp = ncread('temperature.nc','temp');
rh = ncread('humidity.nc','rh');
wind = ncread('wind.nc','wind');
rad = ncread('radiation.nc','rad');
pet = calculate_pet(temp, rh, wind, rad, 'penman_monteith');%% 多时间尺度计算
scales = ;
spei_cell = cell(size(scales));
for i = 1:numel(scales)acc_prec = accumulate(prec, scales(i));acc_pet = accumulate(pet, scales(i));d = acc_prec - acc_pet;spei_cell{i} = calculate_spei(d);
end%% 结果分析
figure;
subplot(2,2,1);
plot_spei(spei_cell{1}, time);
title('1个月尺度SPEI');
subplot(2,2,2);
plot_spei(spei_cell{2}, time);
title('3个月尺度SPEI');
subplot(2,2,3);
plot_spei(spei_cell{3}, time);
title('6个月尺度SPEI');
subplot(2,2,4);
plot_spei(spei_cell{4}, time);
title('12个月尺度SPEI');%% 干旱统计
drought_class = classify_drought(spei_cell{4});
stats = regionprops(drought_class, 'Area');
disp(['极端干旱面积: ', num2str(stats(4).Area), ' km²']);