基于matlab的长时间栅格序列的残差分析

在量化人类活动对生态参量的影响时，常常引入残差分析法来进行量化，本文在matlab平台下对NDVI进行残差分析，选取的自变量是降水和温度。代码如下所示，通过下面代码，能够获取残差的空间分布及残差趋势以及统计得到的研究区每年的残差值。

%author yinlichang3064@163.com
[aa,R]=geotiffread('D:\年NDVI\2000anveage_ndvi.tif');%先投影信息
info=geotiffinfo('D:\年NDVI\2000anveage_ndvi.tif');
[m,n]=size(aa);

begin_year=2000;%用户修改
end_year=2017;%用户修改
long=end_year-begin_year+1;
temsum=zeros(m*n,long);
presum=zeros(m*n,long);
ndvisum=zeros(m*n,long);
k=1;
for year=begin_year:end_year 
     temp=importdata(['D:\年NDVI\',int2str(year),'TEMP.tif']) ; %根据自己名称修改，本数据名称为2000TEMP.tif
     pre=importdata(['D:\年NDVI\',int2str(year),'PRE.tif']) ; %根据自己名称修改，本数据名称为2000PRE.tif
     ndvi=importdata(['D:\年NDVI\',int2str(year),'anveage_ndvi.tif']) ; %根据自己名称修改，本数据名称为'anveage_ndvi.tif'
     %注意数据的有效范围
     temp(temp<-1000)=NaN;%温度有效范围
     pre(pre<0)=NaN;%有效范围大于0
     ndvi(ndvi<-1)=NaN; %有效范围是-1到1
     temsum(:,k)=reshape(temp,m*n,1);
     presum(:,k)=reshape(pre,m*n,1);
     ndvisum(:,k)=reshape(ndvi,m*n,1);
     k=k+1;
end
%多元回归，ndvi=a*pre+b*tem
cc=zeros(m,n)+NaN;
ccsum=zeros(m*n,long)+NaN;
for i=1:m*n
    pre=presum(i,:)';
    if min(pre)>=0 %进行筛选有效范围
        ndvi=ndvisum(i,:)';
        tem=temsum(i,:)';
        X=[ones(size(ndvi)),tem,pre];
        [b,bint,r,rint,stats] = regress(ndvi,X);
        cc1=ndvi-b(1)-b(2).*tem-b(3).*pre;
        cc1=cc1';
        ccsum(i,:)=cc1;
        cc1=cc1';
        X=[ones(size(ndvi)),[1:long]'];
        [b,bint,r,rint,stats] = regress(cc1,X);
        cc(i)=b(2);
    end
end
filename=['D:\年NDVI\2000-2017年残差的趋势.tif'];
geotiffwrite(filename,cc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag)
data2sum=[];
for i=1:long
    data1=ccsum(:,i);
    data2=mean(~isnan(data1));
    data2sum=[data2sum;data2];
    data1=reshape(data1,m,n);
    filename=['D:\年NDVI\残差',int2str(i-1+begin_year),'.tif'];
    geotiffwrite(filename,data1,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag)
end
data2sum=[[begin_year:end_year]',data2sum];
xlswrite('D:\年NDVI\每年的残差均值.xlsx',data2sum)

更多需求，请查看个人介绍