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

haidaowang

%最近写论文需要用到MK检验法，网上收集到大量的matlab代码，但是没有一个代码能够
! Q* F3 x. {! \8 j5 t5 j+ s%完全正确运行或者分析信息不全，结合多位网友编写的MK检验法，经过我的改编，顺利得到
& g) `6 I3 y" y%正确的运行结果，谢谢各位网友，希望对有需要的盆友有帮助
% Mann-Kendall突变检测
% 数据序列y
% 结果序列UFk，UBk2
%--------------------------------------------
%读取excel中的数据，赋给矩阵y
%获取y的样本数
3 i" S+ ^1 L) V: S) O! ^1 f2 |5 t%A为时间和径流数据列
A=xlswrite('数据.xls')
; h6 ^% Y0 F; q8 e: I* U6 m3 \x=A(:,1);%时间序列
$ F' C  G' k# N, j3 M+ l  hy=A(:,2);%径流数据列
5 G0 K" c! e8 T& U) q4 O; Y6 G7 R, ON=length(y);
3 t. G  @" Q: Z" T) \  Mn=length(y);
3 h, c( e, l- F4 ?5 H/ c! g% 正序列计算---------------------------------
$ V1 g& ]  K# e% 定义累计量序列Sk，长度=y，初始值=0
* l! B& u; ?1 _$ u' dSk=zeros(size(y));
; i! m1 R# c, {: R% N6 r& ^% 定义统计量UFk，长度=y，初始值=0
UFk=zeros(size(y));
4 F: c1 }: \9 \7 t2 {% 定义Sk序列元素s
s = 0;
+ k1 {" ?  j3 x7 y% Q: o% i从2开始，因为根据统计量UFk公式，i=1时，Sk(1)、E(1)、Var(1)均为0
% 此时UFk无意义，因此公式中，令UFk(1)=0
. O0 g3 Z, `+ Y6 M7 w2 Tfor i=2:n
   for j=1:i
         if y(i)>y(j)
           s=s+1;
         else
           s=s+0;
         end;
   end;
) S7 ?  q$ X4 K( ]   Sk(i)=s;
1 R" o' T) Q3 p- z7 J" ~   E=i*(i-1)/4; % Sk(i)的均值
. v" `9 g; N! ~, {( Z$ f1 I2 T5 q% R  Var=i*(i-1)*(2*i+5)/72; % Sk(i)的方差
  UFk(i)=(Sk(i)-E)/sqrt(Var);
# G1 j4 _8 o. A2 D1 z3 c! Y& vend;
, O8 N% [  m: D: r/ `6 E% ------------------------------正序列计算end
% 逆序列计算---------------------------------
% 构造逆序列y2，长度=y，初始值=0
y2=zeros(size(y));
% 定义逆序累计量序列Sk2，长度=y，初始值=0
- @4 D, t# c2 W1 d/ t( tSk2=zeros(size(y));
% 定义逆序统计量UBk，长度=y，初始值=0
UBk=zeros(size(y));
: e6 V- w) B% P4 t; G. t% s归0
0 W5 W4 N# r" E! Q5 W" c8 y0 ~s=0;
; \6 {8 R  r1 E, X, c% Q3 a% 按时间序列逆转样本y
% 也可以使用y2=flipud(y);或者y2=flipdim(y,1);
for i=1:n
8 n1 D+ X4 ]8 g! j5 v    y2(i)=y(n-i+1);
end;
% i从2开始，因为根据统计量UBk公式，i=1时，Sk2(1)、E(1)、Var(1)均为0
4 s3 c# s, s* \% 此时UBk无意义，因此公式中，令UBk(1)=0
9 K$ i* |2 M0 A4 y9 I8 W/ {0 `for i=2:n
# K- c6 h. v4 Y- N   for j=1:i
         if y2(i)>y2(j)
           s=s+1;
8 _) F2 ]: i2 I0 y! Z3 q/ |         else
0 R" z* Z( e( ~% r. `! e$ F6 r, ~           s=s+0;
$ S) ]* c& [# o& T6 j9 w7 }6 x# t         end;
   end;
   Sk2(i)=s;
& e$ t' c5 d% e( \+ C2 l+ W7 G* a4 O   E=i*(i-1)/4; % Sk2(i)的均值
2 J* S' r  b9 b& N8 a1 S  Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差
$ q, W8 u: n! D% 由于对逆序序列的累计量Sk2的构建中，依然用的是累加法，即后者大于前者时s加1，
% 则s的大小表征了一种上升的趋势的大小，而序列逆序以后，应当表现出与原序列相反
% h2 X3 b# y0 _. \; Z9 I% 的趋势表现，因此，用累加法统计Sk2序列，统计量公式(S(i)-E(i))/sqrt(Var(i))
% 也不应改变，但统计量UBk应取相反数以表征正确的逆序序列的趋势
2 S9 B9 }2 S; x; I. t' e% s- F1 ?  UBk(i)=0-(Sk2(i)-E)/sqrt(Var);
  P& P# x1 H1 k8 B  T5 V5 B* P. o; \end;
% ------------------------------逆序列计算end
% 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量
% 与UFk做图寻找突变点时，2条曲线应具有同样的时间轴，因此
4 p2 U4 m3 X6 ?! X) C2 w: V; y% 再按时间序列逆转结果统计量UBk，得到时间正序的UBk2，做图用
2 c. v1 t0 y6 D, UUBk2=zeros(size(y));
4 X8 f7 F, w1 `% L2 I% 也可以使用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
- Y$ Y: L, s. N9 W8 S7 }% 目标区域：UFk从A1开始，UBk2从B1开始
xlswrite('f:\test2.xls',UFk,'Sheet1','A1');
  Q2 g" ~. N$ V. [/ i* Sxlswrite('f:\test2.xls',UBk2,'Sheet1','B1');
figure(3)%画图
- E7 j, [+ q5 P: \plot(x,UFk,'r-','linewidth',1.5);
hold on
' g0 {# ~% A8 a* F2 Fplot(x,UBk2,'b-.','linewidth',1.5);
2 D7 A6 F  h% ~) B* T+ t  _plot(x,1.96*ones(N,1),':','linewidth',1);
axis([min(x),max(x),-5,5]);
legend('UF统计量','UB统计量','0.05显著水平');
xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
5 @5 I1 T, i- w! w8 w7 B$ v6 _  Bylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
%grid on
: q+ _8 R, i6 M1 [4 whold on
5 s: [2 _% j6 U! B8 bplot(x,0*ones(N,1),'-.','linewidth',1);
4 ?( o  H! ]7 J+ \1 g  M6 Qplot(x,1.96*ones(N,1),':','linewidth',1);
plot(x,-1.96*ones(N,1),':','linewidth',1);
$ \3 }% L1 ?4 A2 V, s8 \+ P
* j" L: {% n/ D, Z

ExxNEN

谢谢分享程序