在前一篇文章中講述了用sen法進(jìn)行長時(shí)間的趨勢分析，但并未對結(jié)果進(jìn)行顯著性檢驗(yàn)彬伦，通常Sen與MK檢驗(yàn)是結(jié)合在一起的模闲，
因此本文主要講述如何進(jìn)行MK檢驗(yàn)。具體代碼如下
% @author yinlichang3064@163.com
clear
[a,R]=geotiffread('D:\GIS\vegetation\output\yearmax\1982.tif'); %先導(dǎo)入投影信息
info=geotiffinfo('D:\GIS\vegetation\output\yearmax\1982.tif');%先導(dǎo)入投影信息
[m,n]=size(a);
cd=34; %34年忘嫉，時(shí)間跨度
datasum=zeros(m*n,cd)+NaN;
p=1;
for year=1982:2015 %起始年份
filename=['D:\qixiang\年全國8kmPET\china',int2str(year),'pet.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 % 有效格點(diǎn)判定荤牍，我這里有效值在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('C:\MATLAB\MK檢驗(yàn)結(jié)果.tif',zc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag); %注意修改路徑
通過上述代碼的運(yùn)行可以得到MK檢驗(yàn)的結(jié)果。上述代碼運(yùn)行時(shí)只需要修改起始年份和年份長度以及文件的名稱庆冕，注意文件名稱
按照規(guī)律來進(jìn)行分布康吵，本文中的名稱是china1982pet.tif，china1983pet.tif...china2015pet.tif,保證能夠按照規(guī)律讀取访递。
假設(shè)讀者已經(jīng)運(yùn)行完了sen代碼和本文中的代碼晦嵌，則可以得到兩個(gè)tif文件，分別是MK檢驗(yàn)結(jié)果和sen的結(jié)果拷姿，進(jìn)而通過以下代碼
來進(jìn)行最終的判斷
[a,R]=geotiffread('D:\GIS\vegetation\output\yearmax\1982.tif'); %先導(dǎo)入投影信息
info=geotiffinfo('D:\GIS\vegetation\output\yearmax\1982.tif');%先導(dǎo)入投影信息
data=importdata('C:\MATLAB\MK檢驗(yàn)結(jié)果.tif');
sen_value=importdata('D:\zhang\基于sen的pet變化趨勢.tif');
sen_value(abs(data)<1.96)=NaN; %MK結(jié)果值高于1.96則認(rèn)為通過了顯著性95%
geotiffwrite('C:\MATLAB\通過顯著性95%的MK+sen趨勢分析結(jié)果.tif',sen_value,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);%注意修改路徑
更多需求惭载，請查看個(gè)人介紹
