基于matlab实现Sen+MK趋势分析_grace数据做senmk趋势分析在matlab有效值为多少

一、理论
Theil-Sen方法，也就是Sen的斜率估计，是一个用于趋势检测的稳健非参数统计技术，具有较高的计算效率。这个方法能够很好地处理测量误差和异常值，适合分析长期序列数据的变化趋势。其运算公式为：

\beta =mean\left( \frac{X_j-X_i}{j-i} \right) ,2000\leq i<j\leq 202

（以上是AxMath代码）式中，Xj和Xi分别为时间序号为j和i对应的变量；β > 0表示序列在时间上呈增长趋势；β< 0表示呈降低趋势；β = 0，则表明序列趋势不变。

使用M-K（Mann-Kendall）非参数统计方法进行变量的趋势检验，检验方法为:

S=\sum_{i=1}^{n-1}{\sum_{j=i+1}^n{\sgn}}\left( X_j-X_i \right) 
​
​
sgn\left( X_j-X_i \right) =\left\{ \begin{array}{l}
  \quad 1,\quad X_j>X_i\\
  \quad 0,\quad X_j=X_i\\
  \quad -1,\quad X_j<X_i\\
\end{array} \right. 
​
​
Var\left( S \right) =\frac{n\left( n-1 \right) \left( 2n+5 \right)}{18}
​
​
Z=\left\{ \begin{array}{l}
  \frac{S-1}{\sqrt{Var\left( S \right)}}\text{，}S>0\\
  0,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,S=0\\
  \frac{s+1}{\sqrt{Var\left( S \right)}}\text{，}S<0\\
\end{array} \right. 

式中，Z为标准化统计量；var（S）为方差；n为序列中的数据个数。

二、matlab代码

%mk趋势
clear
[a,R]=geotiffread('F:\SIF\CSIF\resample\year\2000.tif'); %先导入投影信息
info=geotiffinfo('F:\SIF\CSIF\resample\year\2000.tif');%先导入投影信息
[m,n]=size(a);
cd=2020-2000+1;       %34年，时间跨度  
datasum=zeros(m*n,cd)+NaN; 
p=1;
for year=2000:2020     %起始年份
     filename=['F:\SIF\CSIF\resample\year\',int2str(year),'.tif'];
    data=importdata(filename);
    data=reshape(data,m*n,1);
    datasum(:,p)=data;         %
    p=p+1;
end
sresult=zeros(m,n)+NaN;
​
for i=1:size(datasum,1)        %
    data=datasum(i,:);
    if min(data)>0       % 有效格点判定，我这里有效值在0以上
        sgnsum=[];  
        for k=2:cd
            for j=1:(k-1)
                sgn=data(k)-data(j);
                if sgn>0
                    sgn=1;
                else
                    if sgn<0
                        sgn=-1;
                    else
                        sgn=0;
                    end
                end
                sgnsum=[sgnsum;sgn];
            end
        end  
        add=sum(sgnsum);
        sresult(i)=add; 
    end
end
vars=cd*(cd-1)*(2*cd+5)/18;
zc=zeros(m,n)+NaN;
sy=find(sresult==0);
zc(sy)=0;
sy=find(sresult>0);
zc(sy)=(sresult(sy)-1)./sqrt(vars);
sy=find(sresult<0);
zc(sy)=(sresult(sy)+1)./sqrt(vars);
geotiffwrite('F:\SIF\CSIF\resample\trend\sif_MK检验结果.tif',zc,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag); %注意修改路径
​
%sen趋势
clear;
clc;
[a,R]=geotiffread('E:\GLEAM_DAY+MOU+YEAR\v3.7a\china_year\E\ET2000_china.tif');%先导入投影信息
info=geotiffinfo('E:\GLEAM_DAY+MOU+YEAR\v3.7a\china_year\E\ET2000_china.tif');
[m,n]=size(a);
cd=2020-2000+1;%根据自己的数据时间跨度修改
datasum=zeros(m*n,cd)+NaN; 
k=1;
for year=2000:2020 %起始年份
    filename=['E:\GLEAM_DAY+MOU+YEAR\v3.7a\china_year\E\ET',int2str(year),'_china.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:\GLEAM_DAY+MOU+YEAR\v3.7a\china_year\trend\E\E_sen.tif'];
geotiffwrite(filename,result,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag)

欢迎各位同学们交流~