#基于matlab、ArcGIS的NDVI的sen mk显著性检验趋势分析

1、sen代码：

% @author yinlichang3064@163.com [a,R]=geotiffread(‘D:\NDVI\7NDVI_Year\2006.tif’);%先导入投影信息 info=geotiffinfo(‘D:\NDVI\7NDVI_Year\2006.tif’); [m,n]=size(a); cd=2020-2006+1;%时间跨度，根据需要自行修改 datasum=zeros(mn,cd)+NaN;
k=1; for year=2006:2020 %起始年份 filename=[‘D:\NDVI\7NDVI_Year\’,int2str(year),’.tif’]; data=importdata(filename); data=reshape(data,mn,1); datasum(:,k)=data;

k=k+1;

end

result=zeros(m,n)+NaN;

for i=1:size(datasum,1)

data=datasum(i,:);
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];
        end
    end
    value=median(valuesum);
    result(i)=value;
end

end

filename=[‘D:\NDVI\7NDVI_Year\趋势\sen.tif’];

geotiffwrite(filename,result,R,’GeoKeyDirectoryTag’,info.GeoTIFFTags.GeoKeyDirectoryTag)

2、mk代码：

% @author yinlichang3064@163.com

clear

[a,R]=geotiffread(‘D:\NDVI\7NDVI_Year\2006.tif’); %先导入投影信息

info=geotiffinfo(‘D:\NDVI\7NDVI_Year\2006.tif’);%先导入投影信息

[m,n]=size(a);

cd=15; %时间跨度

datasum=zeros(m*n,cd)+NaN;

p=1;

for year=2006:2020 %起始年份

 filename=['D:\NDVI\7NDVI_Year\',int2str(year),'.tif'];
data=importdata(filename);
data=reshape(data,m*n,1);
datasum(:,p)=data;         %
p=p+1;

end

sresult=zeros(m,n)+NaN;

for i=1:size(datasum,1) %

data=datasum(i,:);
if min(data)>0       % 有效格点判定，我这里有效值在0以上
    sgnsum=[];  
    for k=2:cd
        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];
        end
    end  
    add=sum(sgnsum);
    sresult(i)=add;  
end

end

vars=cd(cd-1)(2*cd+5)/18;

zc=zeros(m,n)+NaN;

sy=find(sresult==0);

zc(sy)=0;

sy=find(sresult>0);

zc(sy)=(sresult(sy)-1)./sqrt(vars);

sy=find(sresult<0);

zc(sy)=(sresult(sy)+1)./sqrt(vars);

geotiffwrite(‘D:\NDVI\7NDVI_Year\趋势\MK检验结果.tif’,zc,R,’GeoKeyDirectoryTag’,info.GeoTIFFTags.GeoKeyDirectoryTag); %注意修改路径

%注意：以下代码不要用，因为小于1.96的全设置成空白，图像不美观。 [a,R]=geotiffread(‘D:\NDVI\7NDVI_Year\2006.tif’); %先导入投影信息

info=geotiffinfo(‘D:\NDVI\7NDVI_Year\2006.tif’);%先导入投影信息

data=importdata(‘D:\NDVI\7NDVI_Year\趋势\MK检验结果.tif’);

sen_value=importdata(‘D:\NDVI\7NDVI_Year\趋势\sen.tif’);

sen_value(abs(data)<1.96)=NaN; %MK结果值高于1.96则认为通过了显著性95%

geotiffwrite(‘D:\NDVI\7NDVI_Year\趋势\通过显著性95%的MK+sen趋势分析结果.tif’,sen_value,R,’GeoKeyDirectoryTag’,info.GeoTIFFTags.GeoKeyDirectoryTag);%注意修改路径

3、做完前面两步，打开ArcGIS把两幅tif图像导入进去
进行重分类操作，「重分类」工具在「ArcToolbox」-「Spatial Analyst工具」-「重分类」-「重分类」
后续操作参考这里，若链接没了，看如下图
(1)、对sen图进行重分类，设置的值如下(SNDVI就是sen图各像元的NDVI)：

• SNDVI≤−0.0005 植被退化
• −0.0005≤SNDVI≤0.0005 植被生长稳定
• SNDVI≥0.0005 植被改善

• 重分类结果：
  • -1：退化
  • 0：稳定
  • 1：改善

(2)、对mk图进行重分类，设置的值如下：

• |Zs|≤1.96 未通过95%置信度检验，不显著
• |Zs|≥1.96 通过95%置信度检验，显著

• 重分类结果：
  • 1不显著
  • 2显著

4、用「ArcToolbox」-「Spatial Analyst工具」-「数学分析」-「乘」工具将分类后的两幅图像进行(栅格)相乘，得到最终趋势变化图像(参考下图)。后续可以进行改变符号颜色、标注、添加图例等操作。

• -2严重退化
• -1轻微退化
• 0稳定不变
• 1轻微改善
• 2明显改善

参考：

袁丽华,蒋卫国,申文明,刘颖慧,王文杰,陶亮亮,郑华,刘孝富.2000—2010 年黄河流域植被覆盖的时空变化.生态学报,2013,33(24): 7798鄄 7806. NDVI时间序列分析之Sen+MK分析全过程梳理 基于matlab 的长时间栅格数据的Sen+MK显著性检验趋势分析