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>0
                    sgn=1;
                else
                    if sgn<0
                        sgn=-1;
                    else
                        sgn=0;
                    end
                end
                sgnsum=[sgnsum;sgn];  %在sgnsum后面再加上sgn
            end
        end  
        add=sum(sgnsum);
        sresult(i)=add;  %检验统计量S
    end
end
for i=1:size(datasum,1)         
    data=datasum(i,:);          %第i列赋值给data
    if min(data)>0              %判断是否是有效值,我这里的有效值必须大于0  
        valuesum=[];
        for k1=2:cd
            for k2=1:(k1-1)
                cz=data(k1)-data(k2);    %(后减前)
                jl=k1-k2;
                value=cz./jl;
                valuesum=[valuesum;value];  %在valuesum后面再加上value
            end
        end
        value=median(valuesum);   
        result(i)=value;   %Sen趋势度B  B>0上升   B<0下降
     end
end

vars=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;
        end
        elseif 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;
        end
        else
            tread(i)=0;      % 无变化
    end
end
geotiffwrite('E:\Tool\Matlab\MK_Sen\NPP_CASA\NPPtread.tif',tread,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag)%注意修改路径

  三、结果