正确的Mann-Kendall突变检测程序

映山红

魏凤英老师编写的《现代气候统计诊断与预测技术》（第二版）书中介绍了使用Mann-Kendall法来检测时间序列的突变点。大家应该都遇到过这个问题：按照书中所述Mann-Kendall法步骤所编写的突变检测程序，得到的突变点在1930年左右，但书中插图给出的突变点在1920年左右，两者并不一致。

参照网友们编的程序（Mann-Kendall突变分析MATLAB程序（重发修改版）_chenjing_新浪博客 (sina.com.cn)；MATLAB与DPS做Mann-Kendall显著性检验_dreamhottest的博客-CSDN博客_mann-kendall检验 matlab），我们可以得到如下Mann-Kendall突变检测MATLAB程序：
-----------------------------------------------------------------------------------
function [year0,UF,UB2]=MK(year,X,Alpha)
% Mann-Kendall 突变检测
% 输入量
% year:时间序列对应的年份
% X:时间序列,1*N维或N*1维矩阵
% Alpha: 显著水平，默认值为0.05

% 输出量
% year0:突变时间点
% UF: UF统计量，1*N维矩阵
% UB2: UB统计量序列,1*N维矩阵

if (nargin<3)||(isempty(Alpha)), Alpha=0.05;end
%  显著水平Alpha的默认值为0.05

N=length(X);

% 正序列计算
UF=zeros(1,N);  % 时间序列X的UF统计量
Sk=zeros(1,N);  % 时间序列X的秩序列(累计量序列)
r=0;

for i=2:N
   for j=1:i-1
  %原文误写为X(i)>=X(j)，得出突变点在1920年左右，并不正确
      if X(i)>X(j)  
       r=r+1;
      else
       r=r+0;
      end
   end
Sk(i)=r;

% 原文公式错写为 ESk = i*(i+1)/4;
ESk=i*(i-1)/4;
VarSk=i*(i-1)*(2*i+5)/72;
UF(i)=(Sk(i)-ESk)/sqrt(VarSk);
end

% 逆序列计算
Y=X(N:-1:1);     % 取时间序列X的逆序
UB=zeros(1,N);   % 时间序列X的UB统计量
Sk2=zeros(1,N);  % 时间序列Y的秩序
r=0;

for i=2:N
    for j=1:i-1
   %原文误写为Y(i)>=Y(j)，得出突变点在1920年左右，并不正确
       if Y(i)>Y(j)
           r=r+1;
       else
           r=r+0;
       end
    end
    Sk2(i)=r;

% 原文公式错写为 ESk = i*(i+1)/4;
ESk=i*(i-1)/4;
VarSk=i*(i-1)*(2*i+5)/72;
UB(i)=0-(Sk2(i)-ESk)/sqrt(VarSk);
end
% UB是逆序列在逆序时间上的趋势统计量
% 与UF作图寻找突变点时，2条曲线应具有同样的时间轴
% 故再按时间序列逆转UB，得到时间正序的UB2，作图用
UB2=UB(N:-1:1);

% U_Alpha是显著水平Alpha下的正态分位数，alpha=0.05时，U_Alpha=1.96
% 若UF > U_Alpha, 表明序列存在明显的趋势变化
U_Alpha=norminv(1-Alpha/2);


%寻找正序列与逆序列交叉点对应的年份
index=[];
for i=1:N
    if abs(UF(i)-UB2(i))<0.5 %UF与UB2之差小于0.5，即认为二者无差别
       index=[index i];
    end
end
year0=year(index);

figure
h1=plot(year,UF,'r-','linewidth',1);
hold on
h2=plot(year,UB2,'b-','linewidth',1);
YY=get(gca,'Ylim');
axis([min(year) max(year) -max(YY) max(YY)])
xlabel('year')
ylabel('统计量')
hold on
h3=plot(year,U_Alpha*ones(1,N),':','color',[1 0.6 0.3],'linewidth',1.5);
plot(year,-U_Alpha*ones(1,N),':','color',[1 0.6 0.3],'linewidth',1.5)
plot(year,zeros(1,N),'-','linewidth',1)
legend([h1 h2 h3],{'UF统计量','UB统计量','0.05显著水平'},'location','best')
grid on
end

------------------------------------------------------------------------------------------
输入1900-1990年上海年平均气温序列：
X=[15.4  14.6  15.8  14.8  15.0  15.1  15.1  15.0  15.2  15.4 ...
   14.8  15.0  15.1  14.7  16.0  15.7  15.4  14.5  15.1  15.3 ...
   15.5  15.1  15.6  15.1  15.1  14.9  15.5  15.3  15.3  15.4 ...
   15.7  15.2  15.5  15.5  15.6  16.1  15.1  16.0  16.0  15.8 ...
   16.2  16.2  16.0  15.6  15.9  16.2  16.7  15.8  16.2  15.9 ...
   15.8  15.5  15.9  16.8  15.5  15.8  15.0  14.9  15.3  16.0 ...
   16.1  16.5  15.5  15.6  16.1  15.6  16.0  15.4  15.5  15.2 ...
   15.4  15.6  15.1  15.8  15.5  16.0  15.2  15.8  16.2  16.2 ...
   15.2  15.7  16.0  16.0  15.7  15.9  15.7  16.7  15.3  16.1  16.2]';
year=[1900:1990]';


调用函数[year0,UF,UB2]=MK(year,X,Alpha)，检测突变点：
MK(year,X)
结果如下图1，红线与蓝线交叉点即突变点在1930年附近。


上述程序中，计算秩序列Sk时，据书中所述，判断依据是：当xi>xj时，r=1；否则r=0。我们把这个程序记为程序1。如果把判断依据改为：当xi>=xj时，r=1；否则r=0。这个程序记为程序2。程序2中xi=xj也是可以的。把程序1中X(i)>X(j) 和Y(i)>Y(j)改为X(i)>=X(j) 和Y(i)>=Y(j)，就得到程序2。突变检测结果如下图2，这就是书中的插图，红线与蓝线交叉点即突变点在1920年附近。


那么，究竟哪个程序是正确的？兰州干旱气象所的王劲松老师2020年发表在《理论与应用气候》杂志上的论文（Determining the most accurate program for the Mann-Kendall method in detecting climate mutation | SpringerLink）专门研究了这个问题，并认为程序1是正确的。因为使用滑动t检验、Yamamoto法和累积异常法得出的突变点均在1930年代。比如，图3是累积异常法得到的结果，突变点在1932年左右。


Yamamoto法结果如图4，突变点也在1930年代。


李月洪和张正秋1991年发表于《气象》杂志上的文章（百年来上海、北京气候突变的初步分析-A preliminary analysis on abrupt climatic change in Shanghai and Beijing for the last 100 years (nmc.cn)）使用Mann-Kendall检验的结果（下图5）也显示突变点在1932年，气温由1930年代前的负距平转为1930年代后的正距平。


许多文献中介绍Mann-Kendall法的公式书写错误或者不严谨。图6是李月洪等（1991）文中对Mann-Kendall法的介绍：

映山红

对上述Mann-Kendall突变检测程序进行了优化，使得它能在屏幕直接输出准确的突变时间点。建议大家使用下面的程序：
------------------------------------------------------------------------------------------------------------
function [t0,UF,UB2]=MK(t,X,Alpha)
% Mann-Kendall 突变检测
% 输入量
% t:时间序列对应的年份
% X:时间序列,1*N维或N*1维矩阵
% Alpha: 显著水平，默认值为0.05

% 输出量
% t0:突变时间点
% UF: UF统计量，1*N维矩阵
% UB2: UB统计量序列,1*N维矩阵

if (nargin<3)||(isempty(Alpha)), Alpha=0.05;end

N=length(X);

% 正序列计算
UF=zeros(1,N);  % 时间序列X的UF统计量
Sk=zeros(1,N);  % 时间序列X的秩序列(累计量序列)
r=0;

for i=2:N
   for j=1:i-1
  %原文误写为X(i)>=X(j)，得出突变点在1920年左右，并不正确
      if X(i)>X(j)  
       r=r+1;
      else
       r=r+0;
      end
   end
Sk(i)=r;

% 原文公式错写为 ESk = i*(i+1)/4;
ESk=i*(i-1)/4;
VarSk=i*(i-1)*(2*i+5)/72;
UF(i)=(Sk(i)-ESk)/sqrt(VarSk);
end

% 逆序列计算
Y=X(N:-1:1);     % 取时间序列X的逆序
UB=zeros(1,N);   % 时间序列X的UB统计量
Sk2=zeros(1,N);  % 时间序列Y的秩序
r=0;

for i=2:N
    for j=1:i-1
   %原文误写为Y(i)>=Y(j)，得出突变点在1920年左右，并不正确
       if Y(i)>Y(j)
           r=r+1;
       else
           r=r+0;
       end
    end
    Sk2(i)=r;

% 原文公式错写为 ESk = i*(i+1)/4;
ESk=i*(i-1)/4;
VarSk=i*(i-1)*(2*i+5)/72;
UB(i)=0-(Sk2(i)-ESk)/sqrt(VarSk);
end
% UB是逆序列在逆序时间上的趋势统计量
% 与UF作图寻找突变点时，2条曲线应具有同样的时间轴
% 故再按时间序列逆转UB，得到时间正序的UB2，作图用
UB2=UB(N:-1:1);

% U_Alpha是显著水平Alpha下的正态分位数，alpha=0.05时，U_Alpha=1.96
% 若UF > U_Alpha, 表明序列存在明显的趋势变化
U_Alpha=norminv(1-Alpha/2);


%寻找正序列与逆序列交叉点对应的年份
FB=abs(UF-UB2);
index=find(FB==min(FB));
t0=t(index);

figure
h1=plot(t,UF,'r-','linewidth',1);
hold on
h2=plot(t,UB2,'b-','linewidth',1);
xlabel('time')
ylabel('统计量')
hold on
h3=plot(t,U_Alpha*ones(1,N),':','color',[1 0.6 0.3],'linewidth',1.5);
plot(t,-U_Alpha*ones(1,N),':','color',[1 0.6 0.3],'linewidth',1.5)
plot(t,zeros(1,N),'-','linewidth',1)
legend([h1 h2 h3],{'UF统计量','UB统计量','0.05显著水平'},'location','best')
YY=get(gca,'Ylim');
axis([min(t) max(t) -max(abs(YY)) max(abs(YY))])
grid on
end
----------------------------------------------------------------------------------------------------------
在MATLAB中新建脚本，粘贴上述程序，保存为MK.m文件，设置路径，即可直接调用函数MK(t,X)检测突变点

dongyi

总结的很棒

映山红

上述Mann-Kendall突变检测程序仍然存在问题，它只能检测出一个突变点。如果存在多个突变点，怎么办？改进后的程序如下，它能够在屏幕中直接输出所有突变点对应的时间：
-------------------------------------------------------------------------------------------
function [t0,UF,UB2]=MK(t,X,Alpha)
% Mann-Kendall突变检测，输出时间序列的所有突变点
% 输入量
% t:时间序列对应的年份
% X:时间序列,1*N维或N*1维矩阵
% Alpha: 显著水平，默认值为0.05

% 输出量
% t0:突变时间点，可以不止一个
% UF: UF统计量，1*N维矩阵
% UB2: UB统计量序列,1*N维矩阵

if (nargin<3)||(isempty(Alpha)), Alpha=0.05;end
%  显著水平Alpha的默认值为0.05

N=length(X);

% 正序列计算
UF=zeros(1,N);  % 时间序列X的UF统计量
Sk=zeros(1,N);  % 时间序列X的秩序列(累计量序列)
r=0;

for i=2:N
   for j=1:i-1
  %原文误写为X(i)>=X(j)，得出突变点在1920年左右，并不正确
      if X(i)>X(j)  
       r=r+1;
      else
       r=r+0;
      end
   end
Sk(i)=r;

% 原文公式错写为 ESk = i*(i+1)/4;
ESk=i*(i-1)/4;
VarSk=i*(i-1)*(2*i+5)/72;
UF(i)=(Sk(i)-ESk)/sqrt(VarSk);
end

% 逆序列计算
Y=X(N:-1:1);     % 取时间序列X的逆序
UB=zeros(1,N);   % 时间序列X的UB统计量
Sk2=zeros(1,N);  % 时间序列Y的秩序
r=0;

for i=2:N
    for j=1:i-1
   %原文误写为Y(i)>=Y(j)，得出突变点在1920年左右，并不正确
       if Y(i)>Y(j)
           r=r+1;
       else
           r=r+0;
       end
    end
    Sk2(i)=r;

% 原文公式错写为 ESk = i*(i+1)/4;
ESk=i*(i-1)/4;
VarSk=i*(i-1)*(2*i+5)/72;
UB(i)=0-(Sk2(i)-ESk)/sqrt(VarSk);
end
% UB是逆序列在逆序时间上的趋势统计量
% 与UF作图寻找突变点时，2条曲线应具有同样的时间轴
% 故再按时间序列逆转UB，得到时间正序的UB2，作图用
UB2=UB(N:-1:1);

% U_Alpha是显著水平Alpha下的正态分位数，alpha=0.05时，U_Alpha=1.96
% 若UF > U_Alpha, 表明序列存在明显的趋势变化
U_Alpha=norminv(1-Alpha/2);


%寻找正序列与逆序列交叉点对应的年份
FB=UF-UB2;
index=[];
for j=1:N-1
if FB(j)*FB(j+1)<=0
    if abs(FB(j))<abs(FB(j+1))
      index=[index j];
    else  
      index=[index j+1];
    end
end
end
index=unique(index);  %去除矩阵中相同的元素
% 正序列与逆序列交叉点对应的年份
t0=t(index);

figure
h1=plot(t,UF,'r-','linewidth',1);
hold on
h2=plot(t,UB2,'b-','linewidth',1);
xlabel('年份')
ylabel('统计量')
hold on
h3=plot(t,U_Alpha*ones(1,N),':','color',[1 0.6 0.3],'linewidth',1.5);
plot(t,-U_Alpha*ones(1,N),':','color',[1 0.6 0.3],'linewidth',1.5)
plot(t,zeros(1,N),'-','linewidth',1)
legend([h1 h2 h3],{'UF统计量','UB统计量','0.05显著水平'},'location','best')
YY=get(gca,'Ylim');
axis([min(t) max(t) -max(abs(YY)) max(abs(YY))])
grid on
end
----------------------------------------------------------------------------------------------

river

感谢分享，有理有据，分析很透彻

映山红

MK突变检测究竟能否检测出多个突变点？参考M-K突变检验不适合多个突变点-编程作图-气象家园_气象人自己的家园 (06climate.com)，我测试了一下，发现MK法确实不能检测出多个突变点。
-------------------------------------------------------------------------------------------
示例1：利用随机数法，生成阶段式上升序列X，它在时间点30和60处发生了2次均值突变。MK法得到一个交叉点，时间点是45，显然并非真正的突变点。
X=[rand(1,30) rand(1,30)+2 rand(1,30)+4];
figure
plot(X)
grid on
t=1:length(X);
MK(t,X)


示例2：序列X阶段式先升后降，仍在时间点30和60处发生2次均值突变。MK法得到3个交叉点，时间点是78、79、88，仍然不是真正的突变点。
X=[rand(1,30) rand(1,30)+2 rand(1,30)]/2;
figure
plot(X)
grid on
t=1:length(X);
MK(t,X)

吃个肉夹馍

优秀 留个爪

loser.

66666666666

Swallow

映山红 发表于 2021-8-9 11:08
对上述Mann-Kendall突变检测程序进行了优化，使得它能在屏幕直接输出准确的突变时间点。建议大家使用下面的 ...

很需要，不知道有没有python版本的

ironchen

感谢分享!!!!