- 积分
- 445
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-2-13
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
数据:
1957 | 7.4 (℃) | 1958 | 8.1 | 1959 | 7.9 | 1960 | 8.0 | 1961 | 7.9 | 1962 | 7.6 | 1963 | 7.7 | 1964 | 7.5 | 1965 | 7.7 | 1966 | 8.0 | 1967 | 6.9 | 1968 | 7.6 | 1969 | 7.8 | 1970 | 7.1 | 1971 | 8.0 | 1972 | 7.6 | 1973 | 8.2 | 1974 | 7.6 | 1975 | 7.5 | 1976 | 6.9 | 1977 | 7.5 | 1978 | 7.9 | 1979 | 7.7 | 1980 | 7.6 | 1981 | 7.5 | 1982 | 7.7 | 1983 | 7.1 | 1984 | 6.8 | 1985 | 7.5 | 1986 | 7.7 | 1987 | 8.3 | 1988 | 7.6 | 1989 | 7.5 | 1990 | 7.9 | 1991 | 7.9 | 1992 | 7.5 | 1993 | 7.6 | 1994 | 8.1 | 1995 | 8.1 | 1996 | 7.7 | 1997 | 8.6 | 1998 | 9.1 | 1999 | 8.6 | 2000 | 8.5 | 2001 | 8.6 | 2002 | 8.9 | 2003 | 8.5 | 2004 | 8.5 | 2005 | 8.3 | 2006 | 9.4 | 2007 | 8.9 | 2008 | 8.3 | 2009 | 8.7 | 2010 | 8.8 | 2011 | 8.4 |
1.MK检验
程序(显著性=0.05)
% Mann-Kendall突变检测
% 数据序列y
% 结果序列UFk,UBk2
%--------------------------------------------
%读取excel中的数据,赋给矩阵y
%获取y的样本数
%A为时间和径流数据列
A=xlsread('F:\DATA.xls')
x=A(:,1);%时间序列
y=A(:,2);%径流数据列
N=length(y);
n=length(y);
% 正序列计算---------------------------------
% 定义累计量序列Sk,长度=y,初始值=0
Sk=zeros(size(y));
% 定义统计量UFk,长度=y,初始值=0
UFk=zeros(size(y));
% 定义Sk序列元素s
s = 0;
% i从2开始,因为根据统计量UFk公式,i=1时,Sk(1)、E(1)、Var(1)均为0
% 此时UFk无意义,因此公式中,令UFk(1)=0
for i=2:n
for j=1:i
if y(i)>y(j)
s=s+1;
else
s=s+0;
end;
end;
Sk(i)=s;
E=i*(i-1)/4; % Sk(i)的均值
Var=i*(i-1)*(2*i+5)/72; % Sk(i)的方差
UFk(i)=(Sk(i)-E)/sqrt(Var);
end;
% ------------------------------正序列计算end
% 逆序列计算---------------------------------
% 构造逆序列y2,长度=y,初始值=0
y2=zeros(size(y));
% 定义逆序累计量序列Sk2,长度=y,初始值=0
Sk2=zeros(size(y));
% 定义逆序统计量UBk,长度=y,初始值=0
UBk=zeros(size(y));
% s归0
s=0;
% 按时间序列逆转样本y
% 也可以使用y2=flipud(y);或者y2=flipdim(y,1);
for i=1:n
y2(i)=y(n-i+1);
end;
% i从2开始,因为根据统计量UBk公式,i=1时,Sk2(1)、E(1)、Var(1)均为0
% 此时UBk无意义,因此公式中,令UBk(1)=0
for i=2:n
for j=1:i
if y2(i)>y2(j)
s=s+1;
else
s=s+0;
end;
end;
Sk2(i)=s;
E=i*(i-1)/4; % Sk2(i)的均值
Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差
% 由于对逆序序列的累计量Sk2的构建中,依然用的是累加法,即后者大于前者时s加1,
% 则s的大小表征了一种上升的趋势的大小,而序列逆序以后,应当表现出与原序列相反
% 的趋势表现,因此,用累加法统计Sk2序列,统计量公式(S(i)-E(i))/sqrt(Var(i))
% 也不应改变,但统计量UBk应取相反数以表征正确的逆序序列的趋势
UBk(i)=0-(Sk2(i)-E)/sqrt(Var);
end;
% ------------------------------逆序列计算end
% 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量
% 与UFk做图寻找突变点时,2条曲线应具有同样的时间轴,因此
% 再按时间序列逆转结果统计量UBk,得到时间正序的UBk2,做图用
UBk2=zeros(size(y));
% 也可以使用UBk2=flipud(UBk);或者UBk2=flipdim(UBk,1);
for i=1:n
UBk2(i)=UBk(n-i+1);
end;
% 做突变检测图时,使用UFk和UBk2
% 写入目标xlsx文件:C:\Users\chenfei\Desktop\test2.xlsx
% 目标表单:Sheet1
% 目标区域:UFk从A1开始,UBk2从B1开始
xlswrite('F:\test.xls',UFk,'Sheet1','A1');
xlswrite('F:\test.xls',UBk2,'Sheet1','B1');
figure(3)%画图
figure1 = figure('Color',[1 1 1]);
plot(x,UFk,'r-','linewidth',1.5);
hold on
plot(x,UBk2,'b-.','linewidth',1.5);
plot(x,1.96*ones(N,1),':','linewidth',1);
axis([min(x),max(x),-6,6]);
legend('UF统计量','UB统计量','0.05显著水平');
xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
%grid on
hold on
plot(x,0*ones(N,1),'-.','linewidth',1);
plot(x,1.96*ones(N,1),':','linewidth',1);
plot(x,-1.96*ones(N,1),':','linewidth',1);
图1
2.滑动T检验
(1)程序(步长为5,0.05)
clc,clear;
%% 读取数据
excelFile='F:\DATA.XLS';% 文件路径
myData=xlsread(excelFile);% 读取数据
timeSeries=myData(:,1); % 时间序列数据
dataSeries=myData(:,2); %
dataCount=length(dataSeries);
%% 设置步长与检验值
step=5; % 步长
v=step+step-2; % 计算自由度
ttest=2.306; % 查表得t检验值,修改
len1=step;
len2=step;
x=timeSeries(step:dataCount-step);
for i=step:dataCount-step
n1=dataSeries(i-step+1:i);
n2=dataSeries(i+1:i+step);
mean1=mean(n1);
mean2=mean(n2);
c=(len1+len2)/(len1*len2);
var1=1/len1*sum((n1-mean1).^2);
var2=1/len2*sum((n2-mean2).^2);
delta1=len1*var1+len2*var2;
delta=delta1/(len1+len2-2);
t(i-step+1)=(mean1-mean2)/sqrt(delta*c);
end
%% 制图
figure(1);
figure1 = figure('Color',[1 1 1]);
plot(x,t,'r-','linewidth',1.5);
xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
axis([min(x),max(x),-10,10]);
hold on
plot(x,0*ones(i-step+1,1),'-.','linewidth',1);
plot(x,ttest*ones(i-step+1,1),':','linewidth',1);% 更改数字
plot(x,-ttest*ones(i-step+1,1),':','linewidth',1);% 更改数字
legend('t统计量','0.05显著水平');
图2
(2)程序(步长为10,0.05)
ttest=2.101
图3
3.分析
(1)MK有一个交叉点,初步判断1995年前后发生了突变;
(2)T检验步长为5的时候在1994-1995年后发生突变;
(3)T检验步长为10的时候有两处超过显著水平,其中1987年之后的突变明显。
4.结论
1995年发生过一次明显的突变。
|
|