首先要說明的是這是up豬曉曉1074和我共同完成的，她負責數(shù)據(jù)收集以及制圖工作，正是她的工作對我有一點啟發(fā)，才寫了相關的代碼。

眾所周知， 去趨勢波動分析(detrended fluctuation analysis，DFA)算法是由 Peng 等提出的，適合分析信號的長程相關性的標度信號分析方法，能系統(tǒng)地去除序列中的各階趨勢，檢測序列的長程冪律相關性，適用于各種非穩(wěn)態(tài)時間序列研究。去趨勢（detrend）處理也可以消除傳感器在獲取數(shù)據(jù)時產(chǎn)生的偏移對后期計算產(chǎn)生的影響。從數(shù)據(jù)中刪除趨勢可以將分析集中在數(shù)據(jù)趨勢本身的波動上。遙感去趨勢的目的主要是消除傳感器在獲取數(shù)據(jù)時產(chǎn)生的偏移對后期計算產(chǎn)生的影響。下圖是采用去趨勢分析最后所做的圖，像元大小為1km*1km。

主要的代碼如下：

clc

clear

[a,R]=geotiffread('G:2000-2022逐月\2000-2022\china_A2000M02_NDVI1km.tif');%先導入某個圖像的投影信息，為后續(xù)圖像輸出做準確

info=geotiffinfo('G:2000-2022逐月\2000-2022\china_A2000M02_NDVI1km.tif');

[m,n]=size(a);

months=11;

k=1;

T=1;

monthsdata=zeros(m*n,months);

t=1:months;

yeardata=zeros(m*n,23);

for month=2:9

? ?file=['G:2000-2022逐月\2000-2022\','china_A2001M0',int2str(month),'_NDVI1km','.tif'];%注意自己的名字形式，這里使用的名字是年china_A2000M02_NDVI1km.tif，根據(jù)這個可修改

? ?bz=importdata(file);

? ?bz=reshape(bz,m*n,1);

? ?monthsdata(:,k)=bz;

? ?k=k+1;

end

for month=10:12

? ?file=['G:2000-2022逐月\2000-2022\','china_A2001M',int2str(month),'_NDVI1km','.tif'];%注意自己的名字形式，這里使用的名字是年china_A2000M10_NDVI1km.tif.tif，根據(jù)這個可修改

? ?bz=importdata(file);

? ?bz=reshape(bz,m*n,1);

? ?monthsdata(:,k)=bz;

? ?k=k+1;

end

monthdata2=zeros(m*n,1);

for i=1:(m*n)

? ?monthdata1=monthsdata;

? ?monthdata2(i,:)=mean(monthdata1(i,:));

end

yeardata(:,T)=monthdata2;

T=T+1;

months=12;

for year=2000:2022 %起始年份

? ?k=1;

? ?monthsdata=zeros(m*n,months);

? ?for month=1:9

? ?file=['G:2000-2022逐月\2000-2022\','china_A',int2str(year),'M0',int2str(month),'_NDVI1km','.tif'];%注意自己的名字形式，這里使用的名字是年china_A2000M02_NDVI1km.tif，根據(jù)這個可修改

? ?bz=importdata(file);

? ?bz=reshape(bz,m*n,1);

? ?monthsdata(:,k)=bz;

? ?k=k+1;

? ?end

? ?for month=10:12

? ?file=['G:\2000-2022逐月\2000-2022\','china_A',int2str(year),'M',int2str(month),'_NDVI1km','.tif'];%注意自己的名字形式，這里使用的名字是年china_A2000M10_NDVI1km.tif.tif，根據(jù)這個可修改

? ?bz=importdata(file);

? ?bz=reshape(bz,m*n,1);

? ?monthsdata(:,k)=bz;

? ?k=k+1;

? ?end

? ?for i=1:(m*n)

? ?monthdata1=monthsdata;

? ?monthdata2(i,:)=mean(monthdata1(i,:));

? ?end

? ?yeardata(:,T)=monthdata2;

? ?T=T+1;

end

realdata=zeros(m*n,23);

for i=1:(m*n)

? ?realdata1=yeardata(i,:);

? ?realdata2=detrend(realdata1,1);

? ?realdata(i,:)=realdata2;

end

xielv=zeros(m,n);

? ?p=zeros(m,n);

? ?years=23;

for i=1:length(realdata)

? ?bz=realdata(i,:);

? ?if max(bz)>-1 %注意這是進行判斷有效值范圍，如果有效范圍是-1到1，則改成max(bz)>-1即可

? ? ? ?bz=bz';

? ? ? ?X=[ones(size(bz)) bz];

? ? ? ?X(:,2)=[1:years]';

? ? ? ?[b,bint,r,rint,stats] = regress(bz,X);

? ? ? ?pz=stats(3);

? ? ? ?p(i)=pz;

? ? ? ?xielv(i)=b(2);

? ?end

end

trend=reshape(xielv,m*n,1);

p=reshape(p,m*n,1);

for f=1:size(trend,1)

? ?if trend(f)>0

? ? ? ?if p(f)<=0.05

? ? ? ? ? ?trend(f)=1;

? ? ? ?elseif p(f)>0.05

? ? ? ? ? ?trend(f)=2;

? ? ? ?end

? ?elseif trend(f)<0

? ? ? ?if p(f)<=0.05

? ? ? ? ? ?trend(f)=3;

? ? ? ?elseif p(f)>0.05

? ? ? ? ? ?trend(f)=4;

? ? ? ?end

? ?else trend(f)=0;

? ? ? ?trend(f)=0;

? ?end

end

trend=reshape(trend,m,n);

geotiffwrite('G:\2000-2022逐月\結果\去趨勢分析結果\2000-2022年預測.tif',trend,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

這段代碼主要分成三個部分，第一個部分即藍色字部分，是數(shù)據(jù)的讀入功能，由于2000年MODIS數(shù)據(jù)沒有1月數(shù)據(jù)，所以有一個單獨讀取2000年數(shù)據(jù)的位置，我在這里還進行了平均值的計算，即最終是對年平均值的計算，不是月數(shù)據(jù)，尤其注意，而且受到計算機硬件限制，我在這里寫了一個循環(huán)，每次計算一個月的數(shù)據(jù)。

第二個部分，即紅色字的部分，屬于去趨勢部分，其中最重要的為這一句：?

realdata2=detrend(realdata1,1);

我這里只是進行了線性的去趨勢，如果需要去除平均值，可以把后面括號里的1改為0；而進行一階差分的話，需要把0改為2。

第三部分，即綠色字體的部分，是使用一元線性回歸的對數(shù)據(jù)進行趨勢分析的部分，這一部分大家可以參考基于Matlab的柵格數(shù)據(jù)一元線性回歸及顯著性檢驗（slope趨勢分析）這篇文章里我講了具體的思路，大家可以看這一個研究。

我自己的水平也不是很高，公開出來呢只是為了方便后面的人在做到和我類似的工作的時候能夠減少一點在程序上面花費的時間，更好的去分析解決自己所遇到的問題，如果幫到你了就給我一個點贊吧，有問題的話大家可以在評論區(qū)里留言，也可以給我發(fā)私信，大家一起交流進步，我也會不定時的回復大家的。 最后再次感謝up豬曉曉1074的幫助。