用Matlab进行MK检验

用Matlab进行Mann-Kendall分析一、M-K趋势分析法与M-K突变检验的原理介绍1、Mann-Kendall趋势分析法非参数检验，又称任意分布检验，它不对变量的分布做严格假定，检验不针对特定的参数，而是模糊地对变量分布的中心位置或分布状态做检验，由于其不对总体分布做严格假定，因而适用性强[12]。因此，本文采用非参数的Mann-Kendall检验法对昌马河流域近50年的气候水文要素时间序列显著性检验，定量反映变化趋势的显著性。计算公式如下：0,)var(10,00,)var(1SSSSSSSZc(1.1)在公式中，111)sgn(ninikikxxS(1.2)0,10,00,1)sgn((1.3)18)52)(1()52)(1(]var[ttttnnnS(1.4)式中：xk，xi为连续的气候、水文等数据序列，n为数据集合的总长度，t为每个单位的宽度，Σ表示所有单位的总和。衡量趋势大小的指标为：ijjixxMedianji,(1.5)式中：1jin，β所代表的是斜率，正的则是上升，负的代表下降，值的大大小代表趋势明显与否。Mann-Kendall检验如下所示：零假设H0:β=0当21ZZc，拒绝H0假设。式中：21Z为标准正态方差，α为显著性检验水平。2、Mann-kendall突变检验气候系统变化是一个不稳定且不连续的变化过程，而检验其变化的常用方法之一就是Mann-kendall突变检验方法[13]，该方法对于变化要素从一个相对稳定状态变化到另一个状态的变化检验非常有效。且广泛应用于水文，气候，化学，矿物成分检验等各个方面。Mann-kendall突变检验方法如下：对于具有n个样本量的时间序列x，构造一个秩序列：),,3,2(S1knkriki(1.6）其中1,0,ijiijxxrxx（1,2,,ji）(1.7)可见，秩序列kS是第i时刻数值大于j时刻数值个数的累计数。在时间序列随机独立的假定下，定义统计量：),,2,1()()]([nksVarsEsUFkkkk(1.8)其中UF1=0，)(ksE，)(ksVar是累计数ks的均值和方差，在nxxx,,,21相互独立，且有相同连续分布时，它们可由下式算出：4)1()E(sknn(1.9)72)52)(1()Var(sknnn(1.10)iUF为标准正态分布，它是按时间序列x顺序nxxx,,,21计算出的统计量序列，给定显著性水平a，查正态分布表，若aiUUF，则表明序列存在明显的趋势变化。把此方法引用到时间序列的逆序序列中，按xn，xn-1，…x1，再重复上述过程，同时使UFk=-UBk，k=n，n-1，…1，UB=0。给定显著性水平α，将UFk和UBk两个统计量曲线和显著性水平线绘在同一个图上，若UFk和UBk的值大于0，则表明序列呈上升趋势，小于0则呈下降趋势。当超过临界直线时，表明上升或下降趋势显著，超过临界线的范围确定为突变的时间区域。如果UFk和UBk两条曲线出现交点，且交点在临界线之间，那么交点对应的时刻便是突变开始的时间。二、M-K程序介绍1、M-K趋势检验function[slope,zc,za,sign]=MannKendall(x)%计算S%111)sgn(ninikikxxSs=0;len=size(x,2);form=1:len-1forn=m+1:lenifx(n)x(m)s=s+1;elseifx(n)==x(m)s=s+0;elses=s-1;endendend%计算vars和e%18)52)(1()52)(1(]var[ttttnnnSvars=len*(len-1)*(2*len+5)/18;%计算zc%0,)var(10,00,)var(1SSSSSSSZcifs0zc=(s-1)/sqrt(vars);elsezc=(s+1)/sqrt(vars);end%计算zaza=var(x);sign=0;zc1=abs(zc);ifzc1=zasign=1;elsesign=0;end%计算倾斜度ndash=len*(len-1)/2;slope1=zeros(ndash,1);m=1;fork=1:len-1,forj=k+1:len,slope1(m)=(x(j)-x(k))/(j-k);%ijjixxMedianji,m=m+1;end;end;slope=median(slope1);2、M-K突变检验clc;y=data(1:50,6);%输入数据n=length(y);Sk=zeros(size(y));UFk=zeros(size(y));s1=0;fori=2:n%),,3,2(S1knkrikiforj=1:iify(i)y(j)s1=s1+1;elses1=s1+0;%1,0,ijiijxxrxx（1,2,,ji）end;end;Sk(i)=s1;%4)1()E(sknnE=i*(i-1)/4;Var=i*(i-1)*(2*i+5)/72;UFk(i)=(Sk(i)-E)/sqrt(Var);%),,2,1()()]([nksVarsEsUFkkkkend;y2=zeros(size(y));Sk2=zeros(size(y));UBk=zeros(size(y));s2=0;fori=1:ny2(i)=y(n-i+1);end;fori=2:nforj=1:iify2(i)y2(j)s2=s2+1;elses2=s2+0;end;end;Sk2(i)=s2;E=i*(i-1)/4;%Sk2(i)的均值Var=i*(i-1)*(2*i+5)/72;%Sk2(i)的方差%72)52)(1()Var(sknnnUBk(i)=0-(Sk2(i)-E)/sqrt(Var);end;UBk2=zeros(size(y));fori=1:nUBk2(i)=UBk(n-i+1);end;year=1961:2010;year=year';xlswrite('D:\test.xls',year,'Sheet1','A1');xlswrite('D:\test.xls',UFk,'Sheet1','B1');xlswrite('D:\test.xls',UBk2,'Sheet1','C1');