在進行柵格數(shù)據(jù)長時間序列分析的時候，常常涉及到轉(zhuǎn)折點部分，前后的趨勢發(fā)生了很大的變化，如何進行提取呢？本文采用分段線性回歸模型來進行提取。基本思路是對前后兩個序列進行擬合，以殘差平方和最小的點作為潛在的轉(zhuǎn)折點，同時評估轉(zhuǎn)折點前后的趨勢是否通過顯著性檢驗，如果兩邊都通過檢驗，則該潛在轉(zhuǎn)折點為實際轉(zhuǎn)折點，否則不是實際轉(zhuǎn)折點。具體代碼如下：

@author yinlichang3064@163.com
[a,R]=geotiffread('F:\校級課題項目\data\屏障帶\2002降水.tif');%先導(dǎo)入投影信息
info=geotiffinfo('F:\校級課題項目\data\屏障帶\2002降水.tif');
[m,n]=size(a);
begin_year=2000;
end_year=2018;
cd=end_year-begin_year+1;
presum=zeros(m*n,cd);
kk=1;
for year=begin_year:end_year
    filename=strcat('F:\校級課題項目\data\屏障帶\',int2str(year),'降水.tif');
    data=importdata(filename);
    data=reshape(data,m*n,1);
    presum(:,kk)=data;
    kk=kk+1;
end
xlcd=5;% 假設(shè)最短的趨勢年份為5年
BPsum=zeros(m,n);
slope1=zeros(m,n);
slope2=zeros(m,n);
for i=1:m*n
    pre=presum(i,:);
    if min(pre)>0
        res_sum=[];
        p=[];
        slopesum=[];
        for k=xlcd:cd-xlcd
            pre1=pre(1:k)';
            pre2=pre(k:cd)';
            % 
            x1=[ones(size(pre1,1),1),[1:k]'];
            [b1,~,r1,~,stats1] = regress(pre1,x1);
            
            x2=[ones(size(pre2,1),1),[k:cd]'];
            [b2,~,r2,~,stats2] = regress(pre2,x2);
            res=sumsqr(r1)+sumsqr(r2);
            res_sum=[res_sum;res];
            pz=[stats1(3),stats2(3)];
            p=[p;pz];
            slopesum=[slopesum;[b1(2),b2(2)]];
        end
       BP=find(res_sum==min(res_sum));% 找到擬合殘差平方和最小的點
       % 進一步檢驗前后是否通過了顯著性檢驗
       pz=p(BP,:);
       slope=slopesum(BP,:);
       if pz(1)<0.1 && pz(2)<0.1 %以前后顯著性p0.1來判斷
           BPsum(i)=BP+begin_year+xlcd-2;
           slope1(i)=slope(1);
           slope2(i)=slope(2);
       end
    end
end
result_mulu='D:\';
filename1=[result_mulu,'轉(zhuǎn)折點年份.tif'];
geotiffwrite(filename1,BPsum,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

filename2=[result_mulu,'第一階段的趨勢值.tif'];
geotiffwrite(filename2,slope1,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

filename3=[result_mulu,'第二階段的趨勢值.tif'];
geotiffwrite(filename3,slope2,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

在matlab中運行上述代碼即可得到最終的結(jié)果,更多需求，請查看個人介紹，后期有許多內(nèi)容分享在公眾號，歡迎大家捧場關(guān)注一波