经验分享 | SEN+Mk趋势分析(matlab代码分享)

代码分享

方法介绍：Sen 斜率估计用于计算趋势值，通常与MK非参数检验结合使用。即首先计算Sen趋势值，然后使用MK方法判断趋势显著性

示例：1984-2018NDVI年最大值趋势分析

注意：在对NDVI进行趋势分析时，绝对值0.1以下的NDVI值需要去除

代码1：MKTrend(代码2会用)

function MKResult = MKTrend(X,Alpha)
%tic
​
% 时间序列数据的Mann-Kendall趋势分析
%
% 原假设：H0  Beta == 0
% 当 |Z| <= pNorm ,表明在当前显著性水平Alpha下接受原假设（无显著变化）
% 当 Z > pNorm , 表明在当前显著性水平Alpha下认为有上升趋势：此时Beta > 0
% 当 Z < -pNorm ,表明在当前显著性水平Alpha下认为有下降趋势：此时Beta < 0
%
% MKResult = [Z,pNorm];    (统计量S的标准化)Z统计量，上侧分位数
%
% Input
% X                             - 时间序列数据，向量
% Alpha                         - 显著性水平
%
% Output
% Beta                          - 衡量趋势大小的倾斜度或斜率(整体变化趋势的速率)
% MKResult                      - MK趋势检验的Z统计量及Z分位数，[Z,pNorm]
% S_stdS                        - 统计量S及标准差，[S,sqrt(varS)];
% ConInterval                   - 置信区间
%  -----------------------------------------------------
​
if (nargin < 2) || (isempty(Alpha)), Alpha = 0.05; end
​
​
N = length(X);
if N < 8, error('Error:MannKendallTrend: 要求数据样本数目不小于8！'); end
% if N < 41, error('Error:MannKendallTrend: 要求数据样本数目不小于40！'); end
​
S = 0;                                           % 时间序列X的MK趋势检验的统计量S
for ii = 1:(N - 1)
for jj = ii:N
S = S + sign(X(jj) - X(ii));
end
end
​
% 计算方差
varS = ( N * ( N - 1 ) * ( 2 * N + 5 ) ) / 18;
​
% 计算 Z 统计量，该统计量服从标准正态分布
if S == 0
Z = 0;
elseif S > 0
Z = (S - 1) / sqrt(varS);
else
Z = (S + 1) / sqrt(varS);
end
​
% 原假设：H0  Beta == 0
% 当 |Z| <= pNorm ,表明在当前显著性水平Alpha下接受原假设
% 当 Z > pNorm , 表明在当前显著性水平Alpha下认为有上升趋势
% 当 Z < -pNorm ,表明在当前显著性水平Alpha下认为有下降趋势
pNorm = norminv(1 - Alpha / 2);
MKResult = [Z,pNorm];                            % MK趋势检验结果
​
%toc
% end of function MannKendallTrend
%%
end

代码2：调用代码1。MKTrendAnalysis(需要修改的地方已在代码中注释)

clc; clear;imgdir1 = 'D:\homework2\data\data\NDVI year\'; %%修改为所要处理的数据路径addpath(genpath(imgdir1));%% MK趋势分析filenames = dir([imgdir1 '*.tif']);for i = 1:numel(filenames)    data(:,:,i) = single(imread(filenames(i).name));  %% 原始数据 end%%[row,col, N]=size(data);timeslice = N;beg = 1984; %%数据起始年份last = 2018; %%数据结束年份NA = data(1,1,1);​%MK_para=zeros(row,col,2);K=zeros(row,col)*NaN;Z=zeros(row,col)*NaN;X=zeros(1,timeslice)*NaN;t=beg:1:last;Alpha=0.05; %%置信区间for i=1:row    i    for j=1:col        if ismember(data(1,1,1),data(i,j,:))  % 当某位置的时间序列里有无效的数据时, assign NaN to Z and K            Z(i,j)=-9999;            K(i,j)=-9999;        else​                MKResult=MKTrend(data(i,j,:),Alpha);​             X=squeeze(data(i,j,:));             p=polyfit(t',X,1);             K(i,j)=p(1);   %% 变化量                          Z(i,j)=MKResult(1);  %% 显著性                  end            endend%%[~, R0] = readgeoraster('D:\homework2\data\data\NDVI year\ndvi_1984.tif'); %%输入一幅标准的栅格数据来获取属性信息info = geotiffinfo('D:\homework2\data\data\NDVI year\ndvi_1984.tif');  %%输入一幅标准的栅格数据来获取属性信息geoTags = info.GeoTIFFTags.GeoKeyDirectoryTag;outPath = 'D:\homework2\data\data\NDVI year\';  %%输出路径outName1 = [outPath, 'Z.tif'];  %%输出数据名称geotiffwrite(outName1,Z,R0,'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag); %%若输出Slope值，将本行中的Z改为K即可

ARCGIS制图

首先将更改栅格图像Nodata值

在arcgis目录处找到需要处理的栅格数据 - 选中数据 - 单击右键 - 打开属性 - 点击编辑 - 将Nodata值改为-9999 - 计算 - 确定 - 应用

然后对Z值和Slope值分别进行重分类，重分类的结果用栅格计算器相乘，得到最终趋势变化图

其中分类标准如下：

对于Slope：-0.0005 - 0.0005为稳定区域，大于或等于0.0005为植被改善区域，小于-0.0005为植被退化区域

对于Z值：置信水平为0.05，-1.96 - 1.96之间为不显著，绝对值大于1.96为显著。

最终可以分为五类：严重退化、轻微退化、稳定不变、轻微改善、明显改善

结果图如下所示：

END

编辑 | 南波婉帆

审核 | 南波婉琳