最近看到一些博主分享关于 Sen+MK 检验的代码，对于新手来说可能有点复杂。我们编写了一段 MATLAB 代码，能够一次性解决这些问题，简化操作流程。我们还准备了几个关于趋势检验的空间分布图，供大家参考。
一、Sen's Slope和Mann-Kendall检验
在研究中，我们常结合 Theil-Sen 斜率估计（Sen's Slope）和 Mann-Kendall (MK) 检验对栅格数据进行趋势分析。以下是这些方法的具体公式：
二、MATLAB代码
只需要更改栅格数据的输入输出路径、年份和栅格的有效值。
clc;clear;close all;% === 输入文件及参数设置 ===input_folder = 'E:\1981_2021_E\ET\ET\Annual\'; % 输入数据文件夹output_sen_trend = 'E:\1981_2021_E\ET\ET\Annual\result\sen_trend.tif'; % 输出文件路径output_mk = 'E:\1981_2021_E\ET\ET\Annual\result\mk_stat.tif'; % MK检验输出路径output_significant_mk = 'E:\1981_2021_E\ET\ET\Annual\result\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_yearfilename = 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 % 确保有效值大于等于0valuesum = []; % 存储计算的倾斜值% 计算 Sen 趋势（差分比率）for k1 = 2:cdfor k2 = 1:(k1 - 1)cz = data(k1) - data(k2); % 差值jl = k1 - k2; % 时间差value = cz / jl; % 计算倾斜值valuesum = [valuesum; value]; % 存储倾斜值endend% 使用中位数作为 Sen 趋势的值sen_trend(i) = median(valuesum);endend% === 输出 Sen 趋势结果 ===sen_trend = reshape(sen_trend, m, n); % 重新调整为原始地理大小% 检查输出文件路径disp(['保存 Sen 趋势栅格到文件：', output_sen_trend]);% 保存 Sen 趋势结果trygeotiffwrite(output_sen_trend, sen_trend, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);disp('Sen 趋势栅格已成功保存！');catchdisp('保存 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 % 确保数据有效且大于等于0s = 0;for k1 = 2:cdfor k2 = 1:(k1 - 1)diff = data(k1) - data(k2); % 差值s = s + sign(diff); % 累积符号函数endendmk_stat(i) = s; % MK 统计量endend% 计算方差和 Z 值variance = cd * (cd - 1) * (2 * cd + 5) / 18; % 方差计算z_mk = mk_stat ./ sqrt(variance); % Z 值计算z_mk = reshape(z_mk, m, n); % 重新调整为原始地理大小% 保存 MK 统计结果trygeotiffwrite(output_mk, z_mk, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);disp('MK统计栅格已成功保存！');catchdisp('保存 MK统计栅格时出错！');end% === 显著性检验 ===significance_threshold = 1.96; % 显著性阈值 (95% 信心水平)significant_mk = sen_trend; % 复制 Sen 趋势值% 超过显著性阈值的值置为 NaNsignificant_mk(abs(z_mk) < significance_threshold) = NaN;% 保存显著性 MK 结果trygeotiffwrite(output_significant_mk, significant_mk, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);disp('显著性 MK 结果栅格已成功保存！');catchdisp('保存显著性 MK 结果栅格时出错！');enddisp('所有结果已保存！');
三、ArcGIS制图
ABC是MATLAB输出的栅格，我们需要把D做出来。以下是具体的步骤：
1. 将得到三个栅格数据，图中分别是MK检验、SEN趋势和MK显著区域。
2. 对significant_mk栅格转面。
3. 更改矢量的样式。
四、成为会员
