EDA365电子论坛网

标题: 利用MATLAB实现Mann-Kendall突变检测(mk突变检测) [打印本页]

---

作者: haidaowang    时间: 2020-2-3 09:40
标题: 利用MATLAB实现Mann-Kendall突变检测(mk突变检测)

%最近写论文需要用到MK检验法，网上收集到大量的matlab代码，但是没有一个代码能够
%完全正确运行或者分析信息不全，结合多位网友编写的MK检验法，经过我的改编，顺利得到
%正确的运行结果，谢谢各位网友，希望对有需要的盆友有帮助
% Mann-Kendall突变检测
% 数据序列y
% 结果序列UFk，UBk2
%--------------------------------------------
/ f: c5 S% C6 W4 \7 I%读取excel中的数据，赋给矩阵y
+ J$ w! q- l) j* p% Z%获取y的样本数
6 r+ o$ Y7 F* k0 Y6 S$ t! J8 c%A为时间和径流数据列
A=xlswrite('数据.xls')
x=A(:,1);%时间序列
y=A(:,2);%径流数据列
  p( H) J' ]2 |, m& ?; r; }N=length(y);
n=length(y);
) z/ C6 b$ T2 w0 ?1 t' F% 正序列计算---------------------------------
% 定义累计量序列Sk，长度=y，初始值=0
: }' L, ?) {. f; |) pSk=zeros(size(y));
% 定义统计量UFk，长度=y，初始值=0
UFk=zeros(size(y));
7 D. ^5 M3 Q, Z6 T. J" W# {% 定义Sk序列元素s
s = 0;
7 [) [# A- ^( f( Z1 H6 Z1 |% i从2开始，因为根据统计量UFk公式，i=1时，Sk(1)、E(1)、Var(1)均为0
! ^1 {3 v: S. ~8 h% 此时UFk无意义，因此公式中，令UFk(1)=0
for i=2:n
   for j=1:i
+ c1 J" Q& G1 v) {( q3 Q         if y(i)>y(j)
           s=s+1;
! T! K  P  D8 Q' @- [) }' u         else
* j$ I# M7 e. v3 Y6 t4 i' H! n, F; Y           s=s+0;
" d% q! }3 a1 O         end;
   end;
+ v/ e5 {* J" p" E! }   Sk(i)=s;
   E=i*(i-1)/4; % Sk(i)的均值
  Var=i*(i-1)*(2*i+5)/72; % Sk(i)的方差
& g8 o. {3 b+ P' _# p# [  UFk(i)=(Sk(i)-E)/sqrt(Var);
% X% |* o  b1 m: aend;
% ------------------------------正序列计算end
% 逆序列计算---------------------------------
% 构造逆序列y2，长度=y，初始值=0
( {- _: c% ], }; U$ t' jy2=zeros(size(y));
% 定义逆序累计量序列Sk2，长度=y，初始值=0
9 U+ I8 Q6 t7 s3 r9 xSk2=zeros(size(y));
% 定义逆序统计量UBk，长度=y，初始值=0
; U1 F  p  E, t2 T* u" s' [+ CUBk=zeros(size(y));
5 R% o- q0 X3 @; g' k; p3 n5 G% s归0
s=0;
% 按时间序列逆转样本y
% 也可以使用y2=flipud(y);或者y2=flipdim(y,1);
for i=1:n
    y2(i)=y(n-i+1);
end;
# E2 u% D4 ^, }0 n( g% i从2开始，因为根据统计量UBk公式，i=1时，Sk2(1)、E(1)、Var(1)均为0
3 ?! C  S! n+ V6 A8 z% 此时UBk无意义，因此公式中，令UBk(1)=0
for i=2:n
) Q+ I1 t7 `$ z9 R   for j=1:i
, l% D6 g% ?' [7 j2 O         if y2(i)>y2(j)
           s=s+1;
9 X( H& \( l4 B. B# ]         else
# Q' {9 I5 q  {+ g           s=s+0;
% W/ X! }8 r% e         end;
; P. H, c& F8 Z5 s. ^, M# a   end;
& F/ C$ n% w/ s: C, Q   Sk2(i)=s;
   E=i*(i-1)/4; % Sk2(i)的均值
  Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差
0 j3 `" u8 J, c8 f1 n% 由于对逆序序列的累计量Sk2的构建中，依然用的是累加法，即后者大于前者时s加1，
$ `, m  C5 S% ?+ D6 r% 则s的大小表征了一种上升的趋势的大小，而序列逆序以后，应当表现出与原序列相反
% 的趋势表现，因此，用累加法统计Sk2序列，统计量公式(S(i)-E(i))/sqrt(Var(i))
% 也不应改变，但统计量UBk应取相反数以表征正确的逆序序列的趋势
1 R' j) [: i, w  UBk(i)=0-(Sk2(i)-E)/sqrt(Var);
5 A8 z8 D( y; oend;
% ------------------------------逆序列计算end
# F# a* A0 \6 _! n# X; [7 ~% 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量
" n& ^) b2 S+ t% 与UFk做图寻找突变点时，2条曲线应具有同样的时间轴，因此
' X! W0 N8 s0 \+ c) |8 j1 Q: t) r% D% 再按时间序列逆转结果统计量UBk，得到时间正序的UBk2，做图用
2 F8 q; b& y1 Q4 y) kUBk2=zeros(size(y));
/ o- A( g: Y5 M8 g- k. R3 ~% 也可以使用UBk2=flipud(UBk);或者UBk2=flipdim(UBk,1);
for i=1:n
   UBk2(i)=UBk(n-i+1);
6 C  s# H+ N  G) X  U; tend;
9 @) e5 q4 l4 `. c+ \2 @% K0 T% 做突变检测图时，使用UFk和UBk2
3 M+ @0 n: n/ u% 写入目标xls文件：f:\test2.xls
% 目标表单：Sheet1
% W- d4 x* f" N5 P% 目标区域：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 on
; ], G) z9 H2 t5 Q( Gplot(x,UBk2,'b-.','linewidth',1.5);
: N% [9 M4 k$ Aplot(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);
% q- C+ c7 C2 oylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
$ c+ @# _8 ?% S. n%grid on
1 v  R  X- F) }hold on
1 I- W4 w" |) U7 E7 V8 zplot(x,0*ones(N,1),'-.','linewidth',1);
: V; c+ K7 _% p! s  X, ]3 jplot(x,1.96*ones(N,1),':','linewidth',1);
plot(x,-1.96*ones(N,1),':','linewidth',1);
+ S, v3 [& o% @( w$ _) X

4 M5 x+ e; N: N, g0 m3 J: T

---

作者: ExxNEN    时间: 2020-2-3 19:03
谢谢分享程序

---

欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/)

Powered by Discuz! X3.2