利用MATLAB实现Mann-Kendall突变检测(mk突变检测)

haidaowang

%最近写论文需要用到MK检验法，网上收集到大量的matlab代码，但是没有一个代码能够
%完全正确运行或者分析信息不全，结合多位网友编写的MK检验法，经过我的改编，顺利得到
* L4 J; j& j9 R& S; A, r% q2 C%正确的运行结果，谢谢各位网友，希望对有需要的盆友有帮助
% Mann-Kendall突变检测
% 数据序列y
% 结果序列UFk，UBk2
%--------------------------------------------
%读取excel中的数据，赋给矩阵y
* y& ?9 X) j! `9 w2 ]%获取y的样本数
%A为时间和径流数据列
3 v3 F' c: V* V2 p0 W0 OA=xlswrite('数据.xls')
x=A(:,1);%时间序列
% k( T5 X  l4 [1 |5 U- Oy=A(:,2);%径流数据列
N=length(y);
# z7 a8 G! Z9 h3 U4 a" a- on=length(y);
% 正序列计算---------------------------------
. w0 o* Y4 W8 g5 @' B( \% 定义累计量序列Sk，长度=y，初始值=0
Sk=zeros(size(y));
& y( A# h" j6 E' B0 ?! |  R6 F% 定义统计量UFk，长度=y，初始值=0
UFk=zeros(size(y));
% 定义Sk序列元素s
s = 0;
5 e! B3 A% F5 U" F* U7 A% P) K% i从2开始，因为根据统计量UFk公式，i=1时，Sk(1)、E(1)、Var(1)均为0
# Z* X6 @1 b! p% 此时UFk无意义，因此公式中，令UFk(1)=0
7 I6 q1 Y9 \* N- c' X0 M! N* [, ?for i=2:n
   for j=1:i
; d4 g/ r" B  b* n) A! Z# C4 y3 ^         if y(i)>y(j)
7 o) U9 |! ^/ B5 n           s=s+1;
         else
           s=s+0;
% N2 [0 ~5 v. h; ^) [         end;
  D9 _8 s/ x0 |) V) I+ H6 Z4 P   end;
' }) a- m  `1 F( }: P- K8 ]   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);
$ E* |6 u/ v6 Z" i. I  k( Nend;
% ------------------------------正序列计算end
" M9 @6 R  Q- _: S2 [/ A9 w% 逆序列计算---------------------------------
: j# P: t  u" c0 f; s, n0 n% 构造逆序列y2，长度=y，初始值=0
y2=zeros(size(y));
2 ?5 D- T" F. a. D1 Z% 定义逆序累计量序列Sk2，长度=y，初始值=0
Sk2=zeros(size(y));
% 定义逆序统计量UBk，长度=y，初始值=0
UBk=zeros(size(y));
: _! \0 W4 ?& T% F% s归0
* g3 K6 t0 O" p+ X0 @  A% ]0 G1 Ps=0;
% 按时间序列逆转样本y
6 g9 v, k  |) r) p% m$ u% 也可以使用y2=flipud(y);或者y2=flipdim(y,1);
4 T) l6 v, I2 L% w& xfor i=1:n
5 Z7 k8 U4 y4 W1 H. ~$ g    y2(i)=y(n-i+1);
end;
$ }7 k8 o$ Z3 S% i从2开始，因为根据统计量UBk公式，i=1时，Sk2(1)、E(1)、Var(1)均为0
% 此时UBk无意义，因此公式中，令UBk(1)=0
for i=2:n
2 S3 O9 \8 w: q6 B# ]   for j=1:i
4 ]" ?4 |6 s! z8 L         if y2(i)>y2(j)
           s=s+1;
1 J8 |) q, M3 C! i         else
           s=s+0;
         end;
   end;
   Sk2(i)=s;
3 {- ^1 [- Y! |5 V   E=i*(i-1)/4; % Sk2(i)的均值
  Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差
% 由于对逆序序列的累计量Sk2的构建中，依然用的是累加法，即后者大于前者时s加1，
% 则s的大小表征了一种上升的趋势的大小，而序列逆序以后，应当表现出与原序列相反
5 ~' o& t' ?2 l! ~% 的趋势表现，因此，用累加法统计Sk2序列，统计量公式(S(i)-E(i))/sqrt(Var(i))
* u# h8 U, t& {1 {7 V1 ?% 也不应改变，但统计量UBk应取相反数以表征正确的逆序序列的趋势
0 Q5 h( k# P8 ^" f8 N) |% r  y2 L  UBk(i)=0-(Sk2(i)-E)/sqrt(Var);
: Y, B3 Z" ^7 w& T$ [. R# Fend;
% ------------------------------逆序列计算end
% 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量
! y% h, e6 f+ F5 N- h' L; O3 A% 与UFk做图寻找突变点时，2条曲线应具有同样的时间轴，因此
+ I+ H; X4 Q$ W' n- o% 再按时间序列逆转结果统计量UBk，得到时间正序的UBk2，做图用
UBk2=zeros(size(y));
+ S/ u# C7 c/ L% f8 M" L% 也可以使用UBk2=flipud(UBk);或者UBk2=flipdim(UBk,1);
2 m( \+ f% y7 a: L) f3 [4 ^) q# Dfor i=1:n
, p! r/ ?9 o1 Q% ~: F! F* j3 [   UBk2(i)=UBk(n-i+1);
end;
5 Q$ H/ X( [. c  _4 j( X% 做突变检测图时，使用UFk和UBk2
5 s9 c/ C# G+ i0 p1 C- o3 ]" K9 B* D% 写入目标xls文件：f:\test2.xls
% 目标表单：Sheet1
7 H7 i+ S- ^2 g. Y0 U% 目标区域：UFk从A1开始，UBk2从B1开始
xlswrite('f:\test2.xls',UFk,'Sheet1','A1');
xlswrite('f:\test2.xls',UBk2,'Sheet1','B1');
+ A) v1 \+ q$ \4 R% M- hfigure(3)%画图
plot(x,UFk,'r-','linewidth',1.5);
3 P9 T+ Q4 d# N, T+ x$ {hold on
) j* N9 X- {7 {# z) i+ _plot(x,UBk2,'b-.','linewidth',1.5);
plot(x,1.96*ones(N,1),':','linewidth',1);
3 f+ {  f; O- zaxis([min(x),max(x),-5,5]);
legend('UF统计量','UB统计量','0.05显著水平');
, i) W4 E* e; exlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
! v4 q7 S& [# @' R7 t0 i2 P1 n$ k- Q%grid on
/ ^3 i+ g. E4 P# @1 X6 Dhold 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);

* o& }+ \4 @+ Z  v0 s$ V% ]2 [9 q
- X! K, i/ b% U* O
8 Y; U, h9 U1 }

ExxNEN

谢谢分享程序