MK突变检测程序.doc

MannKendall突变检测算法MATLAB源码题目：MannKendall突变检测算法MATLAB源码% 以下是单边MannKendall突变检测算法% GreenSim团队原创作品，转载请注明% Email:% GreenSim团队主页：/greensim% color=red欢迎访问GreenSim算法仿真团队url=/greensim/greensim/url/colorclearload dataX=MJL;N=length(X);U=zeros(N-1,1);for t=2:N x=X(1:t); S=0; n=length(x); for k=1:(n-1) for j=(k+1):n S=S+sign(x(j)-x(k); end end VarS=n*(n-1)*(2*n+5)/18; if S0 Z=(S+1)/sqrt(VarS); elseif S=0 Z=0; else Z=(S-1)/sqrt(VarS); end U(t-1)=Z;endfigure(1)plot(1:(N-1),U,linewidth,1.5);hold onplot(1:(N-1),1.96*ones(N-1,1),:,linewidth,1);legend(统计量,0.05显著水平);hold onplot(1:(N-1),0*ones(N-1,1),-.,linewidth,1);plot(1:(N-1),1.96*ones(N-1,1),:,linewidth,1);plot(1:(N-1),-1.96*ones(N-1,1),:,linewidth,1);axis(1,N-1,-5,5);xlabel(t (month),FontName,TimesNewRoman,FontSize,12);ylabel(U统计量,FontName,TimesNewRoman,Fontsize,12);%grid on%计算趋势显著水平Alpha=1-normcdf(U(end),0,1);disp(趋势的显著水平为);disp(Alpha);%计算整体趋势的变化速率Qi=zeros(N*(N-1)/2,1);counter=1;for k=1:(N-1) for j=(k+1):N Qi(counter)=(X(j)-X(k)/(j-k); counter=counter+1; endendQ=median(Qi);disp(整体趋势的变化速率为);disp(Q);% 以下是双边MannKendall突变检测算法% GreenSim团队原创作品，转载请注明% Email:% GreenSim团队主页：/greensim% color=red欢迎访问GreenSim算法仿真团队url=/greensim/greensim/url/colorclearload dataX=YJL;%计算UF统计量N=length(X);UF=zeros(N-1,1);for t=2:N x=X(1:t); S=0; n=length(x); for k=1:(n-1) for j=(k+1):n if x(j)x(k) S=S+1; else S=S+0; end end end ES=n*(n+1)/4; VarS=n*(n-1)*(2*n+5)/72; Z=(S-ES)/sqrt(VarS); UF(t-1)=Z;end%计算UB统计量Y=flipud(X);UB=zeros(N-1,1);for t=2:N x=Y(1:t); S=0; n=length(x); for k=1:(n-1) for j=(k+1):n if x(j)x(k) S=S+1; else S=S+0; end end end ES=n*(n+1)/4; VarS=n*(n-1)*(2*n+5)/72; Z=(S-ES)/sqrt(VarS); UB(t-1)=-Z;end%绘图figure(2)plot(1:(N-1),UF,r-,linewidth,1.5);hold onplot(1:(N-1),UB,b-.,linewidth,1.5);plot(1:(N-1),1.96*ones(N-1,1),:,linewidth,1);axis(1,N-1,-4,4);legend(UF统计量,UB统计量,0.05显著水平);xlabel(t (year),FontName,TimesNewRoman,FontSize,12);ylabel(统计量,FontName,TimesNewRoman,Fontsize,12);%grid onhold onplot(1:(N-1),0*ones(N-1,1),-.,linewidth,1);plot(1:(N-1),1.96*ones(N-1,1),:,linewidth,1);plot(1:(N-1),-1.96*ones(N-1,1),:,linewidth,1);%最近写论文需要用到MK检验法，网上收集到大量的matlab代码，但是没有一个代码能够%完全正确运行或者分析信息不全，结合多位网友编写的MK检验法，经过我的改编，顺利得到%正确的运行结果，谢谢各位网友，希望对有需要的盆友有帮助% Mann-Kendall突变检测 % 数据序列y% 结果序列UFk，UBk2%-%读取excel中的数据，赋给矩阵y%获取y的样本数%A为时间和径流数据列A=xlswrite(数据.xls)x=A(:,1);%时间序列y=A(:,2);%径流数据列N=length(y);n=length(y);% 正序列计算-% 定义累计量序列Sk，长度=y，初始值=0Sk=zeros(size(y);% 定义统计量UFk，长度=y，初始值=0UFk=zeros(size(y);% 定义Sk序列元素ss = 0;% i从2开始，因为根据统计量UFk公式，i=1时，Sk(1)、E(1)、Var(1)均为0% 此时UFk无意义，因此公式中，令UFk(1)=0for 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，初始值=0y2=zeros(size(y);% 定义逆序累计量序列Sk2，长度=y，初始值=0Sk2=zeros(size(y);% 定义逆序统计量UBk，长度=y，初始值=0UBk=zeros(size(y);% s归0s=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)=0for 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% 写入目标xls文件：f:test2.xls% 目标表单：Sheet1% 目标区域：UFk从A1开始，UBk2从B1开始xlswrite(f:test2.xls,UFk,Sheet1,A1);xlswrite(f:test2.xls,UBk2,Sheet1,B1);figure(3)%画图plot(x,UFk,r-,linewidth,1.5);hold onplot(x,UBk2,b-.,linewidth,1.5);plot(x,1.96*ones(N,1),:,linewidth,1);axis(min(x),max(x),-5,5);legend(UF统计量,UB统计量,0.05显著水平);xlabel(t (year