Theil-Sen Median斜率估计和Mann-Kendall趋势分析：以多年NPP数据为例

一、理论

Theil-Sen Median方法又称为Sen斜率估计，是一种稳健的非参数统计的趋势计算方法。该方法计算效率高，对于测量误差和利群数据不敏感，适用于长时间序列数据的趋势分析。其计算公式为：

Mann-Kendall(MK)检验是一种非参数的时间序列趋势性检验方法，其不需要测量值服从正太分布，不受缺失值和异常值的影响，适用于长时间序列数据的趋势显著检验。其过程如下：对于序列Xt = x1,x2,⋯xn，先确定所有对偶值（ xi , xj , j > i ）中xi与xj的大小关系（设为S）。做如下假设：H0,序列中的数据随机排列，即无显著趋势；H1，序列存在上升或下降趋势。检验统计量S计算公式为（在Z的计算公式中，当S>0时，分子为 S-1）

同样采用双边趋势检验，在给定显著性水平下，在正态分布表中查得临界值 Z1－α/2。当|Z|≤Z1－α/2时，接受原假设，即趋势不显著; 若|Z| ＞Z1－α/2，则拒绝原假设，即认为趋势显著。本文给定显著性水平 α = 0.05，则临界值 Z1－α/2 = ±1.96，当Z的绝对值大于1.65、1.96和2.58时，表示趋势分别通过了信度为90%、95%和99%的显著性检验。趋势显著性的判断方法见表1。

二、代码

clear
[a,R]=geotiffread('E:\Tool\Matlab\MK_Sen\NPP_CASA\NPP_2000.tif'); %先导入投影信息
info=geotiffinfo('E:\Tool\Matlab\MK_Sen\NPP_CASA\NPP_2000.tif');%先导入投影信息
[m,n]=size(a);
cd=2020-2000+1;      %21年，时间跨度
datasum=zeros(m*n,cd)+0;   %生成一个像素个数*年数的矩阵
p=1;
for year=2000:2020      %起止年份filename=['E:\Tool\Matlab\MK_Sen\NPP_CASA\NPP_',int2str(year),'.tif'];     %读入文件名 如D:\qixiang\年全国8kmPET\china2000.pet.tif(china2000pet.tif)data=importdata(filename);  %导入数据data=reshape(data,m*n,1);   %reshape 改变矩阵形式为m*n行、1列datasum(:,p)=data;          %把每年的数据依次放到datasum的每一列p=p+1;
end
sresult=zeros(m,n);
result=zeros(m,n);
for i=1:size(datasum,1)data=datasum(i,:);if min(data)>0       % 有效格点判定，我这里有效值在0以上sgnsum=[];  for k=2:cd       %作用类似于sgn函数    xj-xi>0,sgn=1; xj-xi=0,sgn=0; xj-xi<0,sgn=-1;   (后减前)for j=1:(k-1)sgn=data(k)-data(j);if sgn>0sgn=1;elseif sgn<0sgn=-1;elsesgn=0;endendsgnsum=[sgnsum;sgn];  %在sgnsum后面再加上sgnendend  add=sum(sgnsum);sresult(i)=add;  %检验统计量Send
end
for i=1:size(datasum,1)         data=datasum(i,:);          %第i列赋值给dataif min(data)>0              %判断是否是有效值,我这里的有效值必须大于0  valuesum=[];for k1=2:cdfor k2=1:(k1-1)cz=data(k1)-data(k2);    %(后减前)jl=k1-k2;value=cz./jl;valuesum=[valuesum;value];  %在valuesum后面再加上valueendendvalue=median(valuesum);   result(i)=value;   %Sen趋势度B  B>0上升   B<0下降end
endvars=cd*(cd-1)*(2*cd+5)/18;
zc=zeros(m,n);
sy=find(sresult==0);    %|Z|>1.96变化显著，|Z|<=1.96时变化不显著
zc(sy)=0;                             %S=0时
sy=find(sresult>0);
zc(sy)=(sresult(sy)-1)./sqrt(vars);   %S>0时
sy=find(sresult<0);
zc(sy)=(sresult(sy)+1)./sqrt(vars);   %S<0时result1=reshape(result,m*n,1);
zc1=reshape(zc,m*n,1);
tread=zeros(m,n);
for i=1:size(datasum,1)%result1(i)if result1(i)>0   % Sen趋势B>0  上升     if abs(zc1(i))>=2.58    %极显著上升tread(i)=4;elseif (1.96<=abs(zc1(i)))&&(abs(zc1(i))<2.58)      %显著上升  tread(i)=3;elseif (1.645<=abs(zc1(i)))&&(abs(zc1(i))<1.96)     %微显著上升tread(i)=2;else     %不显著上升tread(i)=1;endelseif result1(i)<0  % Sen趋势B<0  下降 if abs(zc1(i))>=2.58   %极显著下降tread(i)=-4;elseif (1.96<=abs(zc1(i)))&&(abs(zc1(i))<2.58)  %显著下降tread(i)=-3;elseif (1.645<=abs(zc1(i)))&&(abs(zc1(i))<1.96)   %微显著下降tread(i)=-2;else          %不显著下降tread(i)=-1;endelsetread(i)=0;      % 无变化end
end
geotiffwrite('E:\Tool\Matlab\MK_Sen\NPP_CASA\NPPtread.tif',tread,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag)%注意修改路径

三、结果

Theil-Sen Median斜率估计和Mann-Kendall趋势分析：以多年NPP数据为例相关推荐

python实现Theil-Sen Median斜率估计和Mann-Kendall趋势分析
python实现Theil-Sen Median斜率估计和Mann-Kendall趋势分析我的输入数据长这样,直接上代码 # -*- codeing = utf-8 -*-import numpy ...
在模仿中精进数据分析与可视化01——颗粒物浓度时空变化趋势（Mann–Kendall Test）
本文是在模仿中精进数据分析与可视化系列的第一期--颗粒物浓度时空变化趋势(Mann–Kendall Test),主要目的是参考其他作品模仿学习进而提高数据分析与可视化的能力,如果有问题和建议,欢迎在评 ...
【趋势分析方法五】MATLAB实现Sen‘s斜率趋势分析
Sen's斜率趋势分析 Sen's 斜率估计是一种非参数的趋势分析方法,能够有效衡量时间序列的趋势变化量,因而被广泛用于计算水文.气象时间序列趋势的变化程度.该方法基于 Mann-Kendall 趋势 ...
【Matlab】多年NDVI数据的sen趋势分析
具体内容及详细教程请关注微信公众号:夫也的笔记最近老师在课上安排了一个汇报作业,要求对某一地区的多年NDVI数据进行趋势分析,对于小白的我来说,真是!唉!但是功夫不负有心人,在我充分利用百度搜索.谷 ...
基于Python的多时相数据合成
文中的示例代码及数据可关注公众号回复20230105下载,公众号二维码见文末. 1 多时相数据合成由于云覆盖.季节积雪.传感器故障等多种因素的影响,导致从遥感数据中提取的地表参数存在空间分布上的数据 ...
matlab中NDVI趋势分析,【文献阅读笔记】Sen+MK NDVI趋势分析的一些问题
[文献阅读笔记]Sen+MK NDVI趋势分析的一些问题引言 Sen+MK一般用于非参数的趋势估计,最开始用于水文学,计算河流流量的时间变化,后来被引入到多个学科中. Theil–Sen Media ...
【案例教程】地球科学数据（ERA5、雪深、积雪覆盖、海温、植被指数、土地利用）处理实践
[查看原文]地球科学常见数据的处理实践技术应用在地球科学中,不同数据根据具体学科的特色存储为多种数据格式.在科研工作中需要将多种数据进行综合使用分析,因此需要寻找学习 ...
Sen+MK趋势分析
Sen+MK趋势分析结果原理实现非平稳时间序列突变检测 -- Bernaola Galvan分割算法 Sen 斜率估计用于计算趋势值,通常与MK非参数检验结合使用.即首先计算Sen趋势值,然 ...
NDVI时间序列分析之Sen+MK分析全过程梳理
NDVI时间序列分析之Sen+MK分析全过程梳理 Sen斜率估计用于计算趋势值,通常与MK非参数检验结合使用,即先计算Sen趋势值,然后使用MK方法判断趋势显著性. 原理 Theil-Sen Medi ...

Theil-Sen Median斜率估计和Mann-Kendall趋势分析：以多年NPP数据为例

Theil-Sen Median斜率估计和Mann-Kendall趋势分析：以多年NPP数据为例相关推荐

最新文章

热门文章