MATLAB栅格数据计算Hurst指数

nzx123

使用NPP数据进行Hurst指数计算。
裁剪的形状是1km像元大小的756, 502的矩形，无数据的像元值为0
数据格式是浮点型

源代码为
[aa,R]=geotiffread('C:\Users\Tim\Desktop\never\2000.tif');
info=geotiffinfo('C:\Users\Tim\Desktop\never\2000.tif');
nppsum=zeros(size(aa,1)*size(aa,2),15);
for year=2000:2014
    filename=strcat('C:\Users\Tim\Desktop\never\',int2str(year),'.tif');
    npp=importdata(filename);
    npp=reshape(npp,size(npp,1)*size(npp,2),1);
    nppsum=reshape(npp,size(npp,1)*size(npp,2),1);
    nppsum(:,year-1999)=npp;
end
hsum=zeros(size(aa,1),size(aa,2))+NaN;
for kk=1:size(nppsum,1)
    npp=nppsum(kk,:);
    if min(npp)>0
        npp_cf=[];
        for i=1:length(npp)-1
            npp_cf1=npp(i+1)-npp(i);
            npp_cf=[npp_cf,npp_cf1];
        end
        M=[];
        for i=1:size(npp_cf,2)
            M1=mean(npp_cf(1:i));
            M=[M,M1];
        end
        S=[];

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

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

        RS=S(2:size(npp_cf,2)).\RR(2:size(npp_cf,2));
        T=[];
        for i=1:size(npp_cf,2)
            T1=i;
            T=[T,T1];
        end
        lag=T(2:size(npp_cf,2));                  
        g=polyfit(log(lag/2),log(RS),1);               
        H=g(1);
        hsum(kk)=H;
        clear der
    end
end
geotiffwrite('2000-2014山东NPP_Hurst指数3.tif',hsum,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);



计算出的结果如图


请问这是什么原因？

nzx123

没有人知道吗