多时间序列数据MK突变检验突变点提取

之前推送了一篇MK突变分析的文章【Manner-Kendall(M-K)---突变检验】，针对的是如何获得单个时间序列数据的突变点，这个时候突变点可以直接从图中看出。后台有同学询问如何自动提取多个时间序列的突变点，即假如数据有很多列，如何自动判断提取。这个的场景比如我们需要知道区域的降雨突变时间，我们获取了区域内n个气象站点每个站点的逐年降雨观测数据，现在要找到各个站点降雨突变点，从而进行统计，确定区域整个的一个突变点或者突变的时间范围。今天使用Matlab尝试做一下。

1数据31年，100列数据（kk.xlsx的2到101列，代表100个站点）。逐站点提取突变点。2实现

A=xlsread('C:\Users\Desktop\mk突变检验\MKtest\kk.xlsx','Sheet1');x=A(:,1); %时间序列for m=2:101  %第m列 第2列到第101列 代表100个站点 y=A(:,m); %数据列    % MK突变检验代码 ......    ......     % UFK和UBK2交点即为突变点 f1=UFk;    f2=UBk2;    j=1;    for i=1:30 x=i:(i+1); k1=f1(i+1)-f1(i); b1=(i+1)*f1(i)-i*f1(i+1); y1=k1*x+b1; k2=f2(i+1)-f2(i); b2=(i+1)*f2(i)-i*f2(i+1); y2=k2*x+b2; % plot(y1); % hold on % plot(y2); syms x y [xx,yy]=solve(y==k1*x+b1,y==k2*x+b2,x,y); if ~isnan(xx)&xx>=i&xx<(i+1)&yy<1.96&yy>-1.96 pp(j)=double(xx)+1984; j=j+1; end end c{m}=pp;end% 练习数据及完整代码后台回复“多数据突变”获取。3结果站点的突变检验结果储存在c中（1行100 列的元胞），每个站点可以有0到多个突变点。实际中我们可能只需要一个突变点，就需要对检验的突变点进行统计，取频次较大的那个突变时间作为整个区域的突变点。或者划分时间段，统计各个时间段内突变点，取某个时间段作为突变的一个时间范围。这里我们逐年统计突变年份，可以发现2009年发生突变的站点数明显比其年份要多，所以可以认为2009年是整个区域的降雨的突变年份。如需源码，请后台回复：多数据突变有疑问可后台留言或加微信wxid-3ccc交流▼更多精彩推荐，敬请关注我们▼