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

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

文章共1,339字 · 阅读需要大约5分钟

一键AI生成摘要，助你高效阅读

问答

GeoLab 219

13992人浏览 · 2022-06-10 17:23:46

GeoLab 219 · 2022-06-10 17:23:46 发布

代码分享

方法介绍：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即可