Sen斜率计算的matlab代码（可批量计算）

   代码是从网络下载，具体位置已记不清楚，对其进行了部分修改与解释，可进行批量数据的处理，得到多系列的Sen斜率。[注：如果数据系列长度不一致，计算结果会错误，但是最长的数据会有结果]

   欢迎讨论。

   原始与修改的matlab代码如下：

（1）网络下载原始代码【可进行单站斜率估计】

% where   z   = Mann-Kendall Statistic
%         sl  = Sen's Slope Estimate
%         lcl = Lower Confidence Limit of sl
%         ucl = Upper Confidence Limit of sl
%         y   = Time Series of Data
%         dt  = Time Interval of Data

A=xlsread('1.xls')
x=A(:,1);%时间序列
y=A(:,2);%径流数据列
n=length(y);
disp( 'Sen''s Nonparametric Estimator:' );
% calculate slopes
ndash = n * ( n - 1 ) / 2;
s = zeros( ndash, 1 );
i = 1;
for k = 1:n-1,
  for j = k+1:n,
    s(i) = ( y(j) - y(k) ) / ( j - k ) ;
    i = i + 1;
  end;
end;
% the estimate
sl = median( s );
disp( [ '     Slope Estimate = ' num2str( sl ) ] );

%----------------------------------------------------
% the end

（2）修正后的代码【可进行批量计算】

%         sl  = Sen's Slope Estimate
A=xlsread('data.xlsx')%读取样本数据
[a,b]=size(A);%确定数据的行列数为a和b
Senslope=zeros(1,b)%构建预输出结果矩阵，1行，b列
for c=1:b%令c从1列循环至b列
x=A(:,c);%径流数据列
n=length(x);%计算数据系列的长度
ndash = n * ( n - 1 ) / 2;
s = zeros( ndash, 1 );
i = 1;
for k = 1:n-1,
  for j = k+1:n,
    s(i) = ( x(j) - x(k) ) / ( j - k ) ;
    i = i + 1;
  end;
end;
% the estimate
sl = median( s );
Senslope(1,c)=sl;
end
disp( '计算的Sen斜率估计如下：  ' );
disp(Senslope);