【Matlab】水文气象数据分析——MK趋势检验及突变检验

---

前言

• 在时间序列趋势分析中，Mann-Kendall检验是使用广泛的非参数检验方法，是一种定量的方式，被广泛应用于非正态分布的数据趋势分析中，而且该方法可以对数据整体趋势做分析，计算方便。

一、MK趋势检验

1. 定义

• Mann-Kendall单调检验用于检测水文气象时间序列假设检验的趋势，但未指定趋势是线性的还是非线性的。原假设（H0）假设时间序列中不存在显着的增加或减少趋势，而时间序列在备择假设下显示出显著的趋势变化。
• Mann-Kendall单调检验方法中，假定要素的时间序列为X = (x1，x2，… …，xn)其中，n为样本序列的个数，则检测统计量S的计算公式：
  
• Mann(1954)和Kendall(1948)指出，当n≥8时,统计量S近似服从正态分布且均值为E(S)=0，方差Var(S)计算公式如下式：
  

2.代码

% Time Series Trend Detection Tests
% [ z, sl, lcl, ucl ] = trend( y, dt )
% where z = Mann-Kendall Statistic
% sl = Sen's Slope Estimate
% lcl = Lower Confidence Limit of sl
% ucl = Upper Confidence Limit of sl
% y = Time Series of Data
% dt = Time Interval of Data
%y是待检测数据序列
function [ z, sl, lcl, ucl ] = mk( y )
n = length( y );
dt=1;

% Mann-Kendall Test for N > 40
    % disp( 'Mann-Kendall Test;' );
    % if n < 41,
    % disp( 'WARNING - sould be more than 40 points' );
    % end;

% calculate statistic
s = 0;
for k = 1:n-1,
    for j = k+1:n,
        s = s + sign( y(j) - y(k) );
    end;
end;

% variance ( assuming no tied groups )
v = ( n * ( n - 1 ) * ( 2 * n + 5 ) ) / 18;

% test statistic
if s == 0,
z = 0;
elseif s > 0,
z = ( s - 1 ) / sqrt( v );
else
z = ( s + 1 ) / sqrt( v );
end; 

% should calculate Normal value here
nor = 1.96;
% results
disp( [ ' n = ' num2str( n ) ] );
disp( [ ' Mean Value = ' num2str( mean( y ) ) ] );
disp( [ ' Z statistic = ' num2str( z ) ] );
if abs( z ) < nor,
disp( ' No significant trend' );
z = 0;
elseif z > 0,
disp( ' Upward trend detected' );
else
disp( ' Downward trend detected' );
end;
disp( 'Sens Nonparametric Estimator:' );

% calculate slopes
ndash = n * ( n - 1 ) / 2;
s = zeros( ndash, 1 );
i = 1;
for k = 1:n-1,
for j = k+1:n,
s(i) = ( y(j) - y(k) ) / ( j - k ) / dt;
i = i + 1;
end;
end;

% the estimate
sl = median( s );
disp( [ ' Slope Estimate = ' num2str( sl ) ] );
% variance ( assuming no tied groups )
v = ( n * ( n - 1 ) * ( 2 * n + 5 ) ) / 18;
m1 = fix( ( ndash - nor * sqrt( v ) ) / 2 );
m2 = fix( ( ndash + nor * sqrt( v ) ) / 2 );
s = sort( s );
lcl = s( m1 );
ucl = s( m2 + 1 );
disp( [ ' Lower Confidence Limit = ' ...
num2str( lcl ) ] );
disp( [ ' Upper Confidence Limit = ' ...
num2str( ucl ) ] );

3.结果

>> MK
 n = 70
 Mean Value = 1964.5
 Z statistic = 12.2382
 Upward trend detected
Sens Nonparametric Estimator:
 Slope Estimate = 1
 intercept Estimate = 0
 Lower Confidence Limit = 1
 Upper Confidence Limit = 1

二、MK突变检验

1. 定义

• Mann-Kendall突变检验是检验时间序列突变最有效的方法之一，该检验可以明确蚀变的开始时间。本研究采用该方法分析气候因子和流量的变化点。
• 对于具有n个样本量的时间序列x1，x2，… …，xn，构造—秩序列：
  

2.代码

% Mann-Kendall突变检测 
% 数据序列y
 % 结果序列UFk，UBk2
 %--------------------------------------------
 %读取txt中的数据，赋给矩阵y
 %获取y的样本数
%A为时间和数据列
%A=load('d:\matlab2008a32\84data.txt');
%T=load('d:\matlab2008a32\年序列.txt');
load 时间序列.mat
load 暴雨量.mat
A=s


x=T(:,1);%时间序列
y=A(:,1);%气温数据列
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
 % 写入目标xlsx文件：C:\Users\chenfei\Desktop\test2.xlsx
 % 目标表单：Sheet1
 % 目标区域：UFk从A1开始，UBk2从B1开始
%  xlswrite('C:\Users\administrator\Desktop\test2.xlsx',UFk,'Sheet1','A1');
%  xlswrite('C:\Users\administrator\Desktop\test2.xlsx',UBk2,'Sheet1','B1');
figure(3)%画图
set(gcf,'position',[200,300,800,400]);
set(gcf,'color','white')


plot(x,UFk,'r-','linewidth',1.5);
 hold on
 plot(x,UBk2,'b-','linewidth',1.5);
 plot(x,1.96*ones(N,1),'k:','linewidth',2);
%  plot(x,2.56*ones(N,1),'k:.','linewidth',2);
 axis([min(x),max(x),-4,6]);
 legend('UF','UB','\alpha_9_5=\pm1.96');
 xlabel('Year','FontName','Times new roman','FontSize',16,'Fontweight','bold');
 ylabel('Statistic','FontName','Times new roman','FontSize',16,'Fontweight','bold');
 grid on
 hold on
 plot(x,0*ones(N,1),'m-.','linewidth',2);
 plot(x,-1.96*ones(N,1),'k:','linewidth',2);
%  plot(x,-2.56*ones(N,1),'k:.','linewidth',2);
 set(gca,'FontSize',15,'FontName','Times new roman','Fontweight','bold');
 set(gca,'XTicklabel',1965:5:2015);
 

3.结果