hurst指数_hurst栅格_matalhurst指数基于栅格数据_

[mm,R]=geotiffread('E:\Extract_img19821.tif');%先导入投影信息 info=geotiffinfo('E:\Extract_img19821.tif'); Extractsum=zeros(size(mm,1)*size(mm,2),34); for year=1982:2015 filename=['E:\NDVITEST\Extract_img',int2str(year),'1.tif']; ndvi=importdata(filename); ndvi=reshape(ndvi,size(ndvi,1)*size(ndvi,2),1); Extractsum(:,year-1981)=ndvi; end hsum=zeros(size(mm,1),size(mm,2))+NaN; for kk=1:size(Extractsum,1); ndvi=Extractsum(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:\hurst.tif',hsum,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);