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);
转载本文请联系原作者获取授权，同时请注明本文来自李明乾科学网博客。
链接地址：https://m.sciencenet.cn/blog-3459054-1285674.html?mobile=1
收藏
