||
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);
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-27 07:08
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社