利用Matlab实现Mann-Kendall（MK）突变检验函数_mann kendall检验方法怎么确定突变时间

利用Matlab实现Mann-Kendall（MK）突变检验函数

一、MK突变检验

        1、一般取显著性水平α=0.05，那么临界值U0.05= ±1.96 。将UFk和UBk两个统计量序列曲线和±1.96 两条直线均绘在一张图上。

         2、若UFk和UBk的值大于0，则表明序列呈上升趋势，小于0 则表明呈下降趋势。 当它们超过临界直线时，表明上升或下降趋势显著，超过临界线的范围确定为出现突变的时间区域。

        3、如果UFk和UBk两条曲线出现交点，且交点在临界线之间，那么交点对应的时刻便是突变开始的时间。

         4、Mann-Kendall突变检测方法的简要计算步骤:

（1）计算顺序时间序列的秩序列，按照上述公式计算UFx；

（2）计算逆序时间序列的秩序列，按照上述公式计算UBx;

（3）给定显著性水平，如a=0.05，对于临界值为U0.05 = ±1.96，将UFx与UB，两个二个统计量序列曲线与U0.05 =±1.96两条直线绘制在一个平面直角坐标α=0.10对应U0.10 = ±1.28，1=0.01对应U0.01= ±2.32。

（4）分析绘制出的UFR与UB,曲线图，若UF,或UBx的值大于0，则表明序列呈上升趋势，小于0则呈下降趋势。当它们超出临界直线时，表明上升或下降趋势显著。超过临界线的范围确定为出现突变的时间区域。若UFx与UB,两条曲线出现交叉点，且交叉点在临界线之间，它们交叉点对应的时刻便是突变开始的数据。

二、在Matlab上的两种实现方法

代码一

data1=load('D:\1\2.txt');
 
x=data1;
 
year=1987:2021; 
 
%% 突变检验
 
for i=2:length(x)
 
    r(i)=0;
 
    for j=1:i
 
        if x(i)>x(j)
 
           r(i)=r(i)+1;
 
        end
 
    end
 
end
 
for k=2:length(x)
 
    S(k)=sum(r(1:k));
 
    E(k)=k*(k-1)/4;
 
    Var(k)=k*(k-1)*(2*k+5)/72;
 
    UF(k)=(S(k)-E(k))./sqrt(Var(k));
 
end
 
x1=x(end:-1:1);
 
for i=2:length(x)
 
    r1(i)=0;
 
    for j=1:i
 
        if x1(i)>x1(j)
 
           r1(i)=r1(i)+1;
 
        end
 
    end
 
end
 
for k=2:length(x)
 
    S1(k)=sum(r1(1:k));
 
    E1(k)=k*(k-1)/4;
 
    Var1(k)=k*(k-1)*(2*k+5)/72;
 
   UB(k)=-(S1(k)-E1(k))./sqrt(Var1(k));
 
end
 
%% 绘图
 
figure(1)
 
plot(year,data1)
 
xlabel('Year','FontSize',12);
 
ylabel('Sunspot','FontSize',12);
 
set(gca,'FontSize',12);
 
figure(2)
 
plot(year,UF,'r-','MarkerSize',2,'linewidth',1.5);
 
hold on
 
plot(year,UB(end:-1:1),'b-','MarkerSize',2,'linewidth',1.5);
 
plot(year,1.96*ones(length(x),1),'k--','linewidth',1);
 
plot(year,-1.96*ones(length(x),1),'k--','linewidth',1);
 
xlabel('Year','FontSize',12);
 
ylabel('UF&UB','FontSize',12);
 
set(gca,'FontSize',12);
 
legend('UF','UB');

代码二：

 
[filename,pathname] = uigetfile('*.txt','请选择打开的数据文件');
file = [pathname, filename];
data = importdata(file);
x=data(:,1);%时间序列
y=data(:,2);%数据列
 
N=length(y);
 
n=length(y);
 
% 正序列计算---------------------------------
 
% 定义累计量序列Sk，长度=y，初始值=0
 
Sk=zeros(size(y));
 
% 定义统计量UFk，长度=y，初始值=0
 
UFk=zeros(size(y));
 
% 定义Sk序列元素s
 
s = 0;
 
% i从2开始，因为根据统计量UFk公式，i=1时，Sk(1)、E(1)、Var(1)均为0
 
% 此时UFk无意义，因此公式中，令UFk(1)=0
 
for 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，初始值=0
 
y2=zeros(size(y));
 
% 定义逆序累计量序列Sk2，长度=y，初始值=0
 
Sk2=zeros(size(y));
 
% 定义逆序统计量UBk，长度=y，初始值=0
 
UBk=zeros(size(y));
 
% s归0
 
s=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)=0
 
for 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('D:\12.xls',UFk,'Sheet1','A1');
 
xlswrite('D:\12.xls',UBk2,'Sheet1','B1');
 
figure(3)%画图
 
plot(x,UFk,'r-','linewidth',1.5);
 
hold on
 
plot(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)','FontName','TimesNewRoman','FontSize',12);
 
ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
 
%grid on
 
hold 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);

三、参考文献：

Mann-Kendall突变检验原理及实现