基于Theil-Sen斜率和Mann-Kendall检验的栅格数据趋势分析

Sen’s Slope和Mann-Kendall检验

Sen趋势分析
利用时序数据计算每个像素点的Sen趋势值，判断某种变量（如蒸散、降水、温度等）的趋势（上升、下降或无变化）。
使用中位数斜率方法（Sen’s Slope）来量化趋势，适用于处理非参数时间序列数据。
Mann-Kendall显著性检验
对趋势进行显著性检验，判断趋势是否显著（例如，在95%置信水平下）。
计算Mann-Kendall统计量 SS 和对应的标准化Z值，以衡量趋势的显著性。
输出结果
保存趋势结果（Sen趋势值）、Mann-Kendall统计值（Z值）和显著性趋势结果为GeoTIFF文件，供进一步分析和可视化使用。

代码

clc;
clear;
close all;

% === 输入文件及参数设置 ===
input_folder = 'E:\Annual\';  % 输入数据文件夹
output_sen_trend = 'E:\sen_trend.tif';  % 输出文件路径
output_mk = 'E:\mk_stat.tif';  % MK检验输出路径
output_significant_mk = 'E:\significant_mk.tif';  % 显著性MK结果路径

start_year = 1982;  % 起始年份
end_year = 2021;    % 结束年份

sample_file = fullfile(input_folder, sprintf('year_ET_%d.tif', start_year));
[a, R] = geotiffread(sample_file);
info = geotiffinfo(sample_file);
[m, n] = size(a);

cd = end_year - start_year + 1;

% === 数据加载 ===
datasum = zeros(m * n, cd) + NaN;
k = 1;

for year = start_year:end_year
    filename = fullfile(input_folder, sprintf('year_ET_%d.tif', year));
    data = importdata(filename);
    data = reshape(data, m * n, 1);  % 数据展平
    datasum(:, k) = data;  % 将数据存储到数组中
    k = k + 1;
end

% === Sen 趋势计算 ===
sen_trend = NaN(m, n);  % 初始化 Sen 趋势结果矩阵

for i = 1:size(datasum, 1)  
    data = datasum(i, :);  % 获取每个像素点的时间序列数据
    
    if min(data) >= 0  % 确保有效值大于等于0
        valuesum = [];  % 存储计算的倾斜值
        
        % 计算 Sen 趋势（差分比率）
        for k1 = 2:cd
            for k2 = 1:(k1 - 1)
                cz = data(k1) - data(k2);  % 差值
                jl = k1 - k2;  % 时间差
                value = cz / jl;  % 计算倾斜值
                valuesum = [valuesum; value];  % 存储倾斜值
            end
        end        
        % 使用中位数作为 Sen 趋势的值
        sen_trend(i) = median(valuesum);  
    end
end

% === 输出 Sen 趋势结果 ===
sen_trend = reshape(sen_trend, m, n);  % 重新调整为原始地理大小

% 检查输出文件路径
disp(['保存 Sen 趋势栅格到文件：', output_sen_trend]);

% 保存 Sen 趋势结果
try
    geotiffwrite(output_sen_trend, sen_trend, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
    disp('Sen 趋势栅格已成功保存！');
catch
    disp('保存 Sen 趋势栅格时出错！');
end

% === MK 显著性检验 ===
mk_stat = nan(m * n, 1); % 初始化 MK 统计量数组

for i = 1:size(datasum, 1)
    data = datasum(i, :);
    if all(~isnan(data)) && min(data) >= 0  % 确保数据有效且大于等于0
        s = 0;
        for k1 = 2:cd
            for k2 = 1:(k1 - 1)
                diff = data(k1) - data(k2); % 差值
                s = s + sign(diff); % 累积符号函数
            end
        end
        mk_stat(i) = s; % MK 统计量
    end
end

% 计算方差和 Z 值
variance = cd * (cd - 1) * (2 * cd + 5) / 18; % 方差计算
z_mk = mk_stat ./ sqrt(variance); % Z 值计算
z_mk = reshape(z_mk, m, n); % 重新调整为原始地理大小

% 保存 MK 统计结果
try
    geotiffwrite(output_mk, z_mk, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
    disp('MK统计栅格已成功保存！');
catch
    disp('保存 MK统计栅格时出错！');
end

% === 显著性检验 ===
significance_threshold = 1.96; % 显著性阈值 (95% 信心水平)
significant_mk = sen_trend; % 复制 Sen 趋势值

% 超过显著性阈值的值置为 NaN
significant_mk(abs(z_mk) < significance_threshold) = NaN;

% 保存显著性 MK 结果
try
    geotiffwrite(output_significant_mk, significant_mk, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
    disp('显著性 MK 结果栅格已成功保存！');
catch
    disp('保存显著性 MK 结果栅格时出错！');
end

disp('所有结果已保存！');

出图效果