基于Matlab计算栅格数据Hurst指数和未来趋势

?各位同学们好，今天分享的基于Matlab计算栅格数据Hurst指数和未来趋势。如果您需要下载或处理遥感数据等方面的帮助，私信或评论。

一、Hurst指数

Hurst指数是一种用于描述未来短时间内变化趋势可持续性的指标，可以在分析年际变化特征方面提供更好的帮助。该指数由英国水文专家Hurst提出，并在地质、遥感和水文等领域得到广泛应用。其中，研究使用了重标极差分析法(R/S)来计算Hurst指数，并应用于分析未来短期变量的变化趋势。这一方法的基本原理包括以下几个部分。

H=0.5时，表示序列是无一致性随机序列，即无法判断未来的趋势。

当0.5<H<1时，表明时间序列存在长期记忆性，即未来的趋势与现在的趋势一致。

当0≤H<0.5时，表示时间序列是不一致的，即未来的趋势与现在的趋势不一致。

二、变化趋势

一般将Slope和Hurst指数叠加得出变量在未来的变化趋势，即

三、Matlab代码

1. Hurst指数计算

clc
clear
[aa,R]=readgeoraster('E:\MOD13A3\qilianshan\2000.tif');
info=geotiffinfo('E:\MOD13A3\qilianshan\2000.tif');
ndvisum=zeros(size(aa,1)*size(aa,2),21);
for year=2000:2020
    filename=strcat('E:\MOD13A3\qilianshan\',int2str(year),'.tif');
    ndvi=importdata(filename);
    ndvi=reshape(ndvi,size(ndvi,1)*size(ndvi,2),1);
    ndvisum(:,year-1999)=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
geotiffwrite('E:\MOD13A3\qilianshan\Hurst\NDVI_Hurst指数.tif',hsum,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
disp('ok!')

2. 未来趋势检验

2.1 Sen趋势检验

%sen趋势
clear;
clc;
[a,R]=geotiffread('E:\MOD13A3\qilianshan\2000.tif');%先导入投影信息
info=geotiffinfo('E:\MOD13A3\qilianshan\2000.tif');
[m,n]=size(a);
cd=2020-2000+1;%根据自己的数据时间跨度修改
datasum=zeros(m*n,cd)+NaN; 
k=1;
for year=2000:2020 %起始年份
    filename=['E:\MOD13A3\qilianshan\',int2str(year),'.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)>-1 %根据自己数据有效值修改，我这里的有效值必须大于-1
        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:\MOD13A3\qilianshan\Hurst\NDVI_sen.tif'];
geotiffwrite(filename,result,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag)

2.2?未来趋势计算

clc
clear
[a,R]=readgeoraster('E:\MOD13A3\qilianshan\Hurst\去Nodate\NDVI_HURST.TIF');  
info=geotiffinfo('E:\MOD13A3\qilianshan\Hurst\去Nodate\NDVI_HURST.TIF');
trend=importdata('E:\MOD13A3\qilianshan\Hurst\去Nodate\NDVI_SEN.TIF');
h=importdata('E:\MOD13A3\qilianshan\Hurst\去Nodate\NDVI_SEN.TIF');
[m,n]=size(trend);
result=zeros(m,n)+6;  % 6是未通过显著性检验
for i=1:m
    for j=1:n
        if trend(i,j)<1&&trend(i,j)>-1   % 有趋势的像元
            if trend(i,j)>0&&h(i,j)<1&&h(i,j)>0.5  
                result(i,j)=1;   % 增-增  
            elseif trend(i,j)>0&&h(i,j)<0.5&&h(i,j)>0
                result(i,j)=2;   % 增-减 
            elseif trend(i,j)<0&&h(i,j)<1&&h(i,j)>0.5 
                result(i,j)=3;   % 减-减 
            elseif trend(i,j)<0&&h(i,j)<0.5&&h(i,j)>0 
                result(i,j)=4;   % 减-增 
            else
                result(i,j)=5;   % 随机
            end
        end   
    end
end
geotiffwrite('E:\MOD13A3\qilianshan\Hurst\NDVI未来趋势.tif',result,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
disp('ok!')