Theil-Sen median趋势分析是一种稳健的非参数计算方法,今天介绍的是Theil-Sen median趋势分析,一般还要结合Mann-Kendall建议,这个MK检验后续文章会写。
计算式如下
公式中SET指计算n(n-1)/2 个数据组合的斜率的中位数,ETi和ETj代表i和j年的ET值。如果SET>0,则ET呈上升趋势,否则,ET呈下降趋势。
Matlab具体代码如下:
[a,R]=geotiffread('E:\GIS\NDVI\2000.tif');%先导入投影信息
info=geotiffinfo('E:\GIS\NDVI\2000.tif');
[m,n]=size(a);
cd=2020-2000+1;%时间跨度,根据需要自行修改
datasum=zeros(m*n,cd)+NaN;
k=1;
for year=2000:2020 %起始年份
filename=['E:\GIS\ENVI\xinjiang',int2str(year),'2000.tif'];
data=importdata(filename);
data=reshape(data,m*n,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=['E:\GIS\基于SEN的NDVI变化趋势.tif'];
geotiffwrite(filename,result,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag)