MATLAB对NDVI去趋势分析代码及思路阐述
首先要说明的是这是up猪晓晓1074和我共同完成的,她负责数据收集以及制图工作,正是她的工作对我有一点启发,才写了相关的代码。
众所周知, 去趋势波动分析(detrended fluctuation analysis,DFA)算法是由 Peng 等提出的,适合分析信号的长程相关性的标度信号分析方法,能系统地去除序列中的各阶趋势,检测序列的长程幂律相关性,适用于各种非稳态时间序列研究。去趋势(detrend)处理也可以消除传感器在获取数据时产生的偏移对后期计算产生的影响。从数据中删除趋势可以将分析集中在数据趋势本身的波动上。遥感去趋势的目的主要是消除传感器在获取数据时产生的偏移对后期计算产生的影响。下图是采用去趋势分析最后所做的图,像元大小为1km*1km。

主要的代码如下:
clc
clear
[a,R]=geotiffread('G:2000-2022逐月\2000-2022\china_A2000M02_NDVI1km.tif');%先导入某个图像的投影信息,为后续图像输出做准确
info=geotiffinfo('G:2000-2022逐月\2000-2022\china_A2000M02_NDVI1km.tif');
[m,n]=size(a);
months=11;
k=1;
T=1;
monthsdata=zeros(m*n,months);
t=1:months;
yeardata=zeros(m*n,23);
for month=2:9
file=['G:2000-2022逐月\2000-2022\','china_A2001M0',int2str(month),'_NDVI1km','.tif'];%注意自己的名字形式,这里使用的名字是年china_A2000M02_NDVI1km.tif,根据这个可修改
bz=importdata(file);
bz=reshape(bz,m*n,1);
monthsdata(:,k)=bz;
k=k+1;
end
for month=10:12
file=['G:2000-2022逐月\2000-2022\','china_A2001M',int2str(month),'_NDVI1km','.tif'];%注意自己的名字形式,这里使用的名字是年china_A2000M10_NDVI1km.tif.tif,根据这个可修改
bz=importdata(file);
bz=reshape(bz,m*n,1);
monthsdata(:,k)=bz;
k=k+1;
end
monthdata2=zeros(m*n,1);
for i=1:(m*n)
monthdata1=monthsdata;
monthdata2(i,:)=mean(monthdata1(i,:));
end
yeardata(:,T)=monthdata2;
T=T+1;
months=12;
for year=2000:2022 %起始年份
k=1;
monthsdata=zeros(m*n,months);
for month=1:9
file=['G:2000-2022逐月\2000-2022\','china_A',int2str(year),'M0',int2str(month),'_NDVI1km','.tif'];%注意自己的名字形式,这里使用的名字是年china_A2000M02_NDVI1km.tif,根据这个可修改
bz=importdata(file);
bz=reshape(bz,m*n,1);
monthsdata(:,k)=bz;
k=k+1;
end
for month=10:12
file=['G:\2000-2022逐月\2000-2022\','china_A',int2str(year),'M',int2str(month),'_NDVI1km','.tif'];%注意自己的名字形式,这里使用的名字是年china_A2000M10_NDVI1km.tif.tif,根据这个可修改
bz=importdata(file);
bz=reshape(bz,m*n,1);
monthsdata(:,k)=bz;
k=k+1;
end
for i=1:(m*n)
monthdata1=monthsdata;
monthdata2(i,:)=mean(monthdata1(i,:));
end
yeardata(:,T)=monthdata2;
T=T+1;
end
realdata=zeros(m*n,23);
for i=1:(m*n)
realdata1=yeardata(i,:);
realdata2=detrend(realdata1,1);
realdata(i,:)=realdata2;
end
xielv=zeros(m,n);
p=zeros(m,n);
years=23;
for i=1:length(realdata)
bz=realdata(i,:);
if max(bz)>-1 %注意这是进行判断有效值范围,如果有效范围是-1到1,则改成max(bz)>-1即可
bz=bz';
X=[ones(size(bz)) bz];
X(:,2)=[1:years]';
[b,bint,r,rint,stats] = regress(bz,X);
pz=stats(3);
p(i)=pz;
xielv(i)=b(2);
end
end
trend=reshape(xielv,m*n,1);
p=reshape(p,m*n,1);
for f=1:size(trend,1)
if trend(f)>0
if p(f)<=0.05
trend(f)=1;
elseif p(f)>0.05
trend(f)=2;
end
elseif trend(f)<0
if p(f)<=0.05
trend(f)=3;
elseif p(f)>0.05
trend(f)=4;
end
else trend(f)=0;
trend(f)=0;
end
end
trend=reshape(trend,m,n);
geotiffwrite('G:\2000-2022逐月\结果\去趋势分析结果\2000-2022年预测.tif',trend,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
这段代码主要分成三个部分,第一个部分即蓝色字部分,是数据的读入功能,由于2000年MODIS数据没有1月数据,所以有一个单独读取2000年数据的位置,我在这里还进行了平均值的计算,即最终是对年平均值的计算,不是月数据,尤其注意,而且受到计算机硬件限制,我在这里写了一个循环,每次计算一个月的数据。
第二个部分,即红色字的部分,属于去趋势部分,其中最重要的为这一句:
realdata2=detrend(realdata1,1);
我这里只是进行了线性的去趋势,如果需要去除平均值,可以把后面括号里的1改为0;而进行一阶差分的话,需要把0改为2。
第三部分,即绿色字体的部分,是使用一元线性回归的对数据进行趋势分析的部分,这一部分大家可以参考基于Matlab的栅格数据一元线性回归及显著性检验(slope趋势分析)这篇文章里我讲了具体的思路,大家可以看这一个研究。
我自己的水平也不是很高,公开出来呢只是为了方便后面的人在做到和我类似的工作的时候能够减少一点在程序上面花费的时间,更好的去分析解决自己所遇到的问题,如果帮到你了就给我一个点赞吧,有问题的话大家可以在评论区里留言,也可以给我发私信,大家一起交流进步,我也会不定时的回复大家的。 最后再次感谢up猪晓晓1074的帮助。