2021-10-24

Hurst指数——长时间序列NDVI预测（Matlab）

    基于重标极差(Ｒ/S)分析方法的 Hurst 指数最早是由英国水文学家 Hurst 在研究尼罗河水库流量和储存能力的关系时提出 。本文采用常用的Ｒ/S 分析法，以1982-2015年NDVI变化状况为例，进行逐像元分析，最后形成一幅栅格图层。
   本文利用Matlab来进行实现。主要代码如下：

下面展示一些 内联代码片。

//处理范围主要依据读取的tif数据文件，主体代码基本满足需要不需要过多更改。

clear
[aa,R]=geotiffread('D:\Global\NDVI3g\GIMMSraster\raster\最大合成\GIMMS_NDVI2015.tif');**%获取投影信息（顺便添加一幅tiff即可）**
info=geotiffinfo('D:\Global\NDVI3g\GIMMSraster\raster\最大合成\GIMMS_NDVI2015.tif');
ndvisum=zeros(size(aa,1)*size(aa,2),34); **%34为所使用的tiff数据的数目，本文使用的是1982-2015的年NDVI数据34期**
for year=1982:2015
    filename=strcat('D:\Global\NDVI3g\GIMMSraster\raster\最大合成\GIMMS_NDVI',int2str(year),'.tif');**%NDVI数据存储位置，由于本文的tiff命名格式：GIMMS_NDVI年份.tif，例如，GIMMS_NDVI1982.tif，所以strcat('**文件地址**\**前缀**（如，GIMMS_NDVI）',int2str()year,'**保存格式**（.tif）')**
    ndvi=importdata(filename);
    ndvi=reshape(ndvi,size(ndvi,1)*size(ndvi,2),1);
    ndvisum(:,year-1981)=ndvi;
end
hsum=zeros(size(aa,1),size(aa,2))+NaN;
for kk=1:size(ndvisum,1);
    ndvi=ndvisum(kk,:);
    if min(ndvi)>0
        ndvi_cf=[];
        for i=1:length(ndvi)-1
            ndvi_cf1=ndvi(i+1)-ndvi(i);
            ndvi_cf=[ndvi_cf,ndvi_cf1];
        end
        M=[];
        for i=1:size(ndvi_cf,2)
            M1=mean(ndvi_cf(1:i));
            M=[M,M1];
        end
        S=[];

        for i=1:size(ndvi_cf,2)
            S1=std(ndvi_cf(1:i))*sqrt((i-1)/i);
            S=[S,S1];
        end

        for i=1:size(ndvi_cf,2)
            for j=1:i
                der(j)=ndvi_cf(1,j)-M(1,i);
                cum=cumsum(der);
                RR(i)=max(cum)-min(cum);
            end
        end

        RS=S(2:size(ndvi_cf,2)).\RR(2:size(ndvi_cf,2));
        T=[];
        for i=1:size(ndvi_cf,2)
            T1=i;
            T=[T,T1];
        end
        lag=T(2:size(ndvi_cf,2));                  
        g=polyfit(log(lag/2),log(RS),1);                
        H=g(1);
        hsum(kk)=H;
        clear der
    end
end
outname='D:\result\Hurst.tif';**%保存位置（自行设定）**
geotiffwrite(outname,hsum,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);