用Matlab进行MK检验.doc

因而适用性强［12］因此，本文采用非参数的Mann-Kendall检验法对昌马河流域近50年的气用 Matlab 进行 Mann-Kendall 分析、M-K趋势分析法与M-K突变检验的原理介绍1、Mann-Kendall趋势分析法非参数检验，又称任意分布检验， 它不对变量的分布做严格假定， 检验不针对特定的参由于其不对总体分布做严格假定，数，而是模糊地对变量分布的中心位置或分布状态做检验,i候水文要素时间序列显著性检验，定量反映变化趋势的显著性计算公式如下:在公式中,S -1S a0Jvar( S)彳 0,S = 0S +1S c0(Jvar( S)Zc(1.1)n d nSAL sgn(Xk -xji 吕 k 1(1.2)1, 0sgn 但)=* 0, 日=0-1, 9 <0(1.3)var[S] = |n(n—1)(2n 十 5) —E t(t — 1)(2t 十5^/‘18- t J(1.4)式中：Xk , Xi为连续的气候、水文等数据序列,n为数据集合的总长度，t为每个单位的宽度，工表示所有单位的总和衡量趋势大小的指标为:Xi_Xj(1.5)式中：1

Mann-Kendall检验如下所示:零假设H： 3 =0当Zc >Zy2,拒绝H>假设式中：-Z^-'2为标准正态方差，a为显者性检验水平2、Mann-kendall 突变检验气候系统变化是一个不稳定且不连续的变化过程，而检验其变化的常用方法之一就是Mann-kendall突变检验方法［13］，该方法对于变化要素从一个相对稳定状态变化到另一个状态的变化检验非常有效且广泛应用于水文,气候，化学，矿物成分检验等各个方面Mann-ke ndall突变检验方法如下:对于具有n个样本量的时间序列 x，构造一个秩序列:其中Sk 二"ri(k = 2,3「 ，n)(1.61,Xir0,Xi < XjXj(j -1,2, ,i) (1.7)可见，秩序列 2是第i时刻数值大于j时刻数值个数的累计数在时间序列随机独立的假定下，定义统计量:UFk」Sk「E(Sk)]xVar (Sk)(k = 1,2厂,n) (1.8)其中UF1=0, E(Sk)，Var(s<)是累计数2的均值和方差，在为公2，…，人相互独立,且有相同连续分布时，它们可由下式算出:E(sQ二虫日(1.9)Var(s k)二n(n -1)( 2n 5)72(1.10),Xn计算出的统计量序列，给定UFi为标准正态分布，它是按时间序列 X顺序为公2,显著性水平a,查正态分布表，若 UFi Ua，则表明序列存在明显的趋势变化。

把此方法 引用到时间序列的逆序序列中，按 xn , xn-1，…x1，再重复上述过程，同时使 UFk=-UBk,k=n, n-1，…1, UB=0>给定显著性水平 a，将UFk和UBk两个统计量曲线和显著性水平线 绘在同一个图上，若UFk和UBk的值大于0,则表明序列呈上升趋势， 小于0则呈下降趋势当超过临界直线时，表明上升或下降趋势显著，超过临界线的范围确定为突变的时间区域如果UFk和UBk两条曲线出现交点，且交点在临界线之间， 那么交点对应的时刻便是突变开始的时间二、M-K程序介绍1、M-K趋势检验function [slope,zc,za,sign]=MannKendall(x)n J n%计算 s %S = 、sgn(xk - xi)i -1 k ：[]1s=0;len=size(x,2);for m=1:le n-1for n=m+1:le nif x(n) >x(m)s=s+1;elseif x(n)==x(m)s=s+0;elses=s-1;endendendvar[S] = R(n -1)(2 n *5)—送 t(t —1)(2t +5)“18 %计算vars和e% L tvars=le n*(le n-1)*(2*le n+5)/18;Zc%计算zc %ST 7 S h U0zc=(s-1)/sqrt(vars); elsezc=(s+1)/sqrt(vars); end%计算zaza=var(x);sig n=0; zc仁abs(zc); if zc1 >= zasig n=1; elsesig n=0; end%计算倾斜度 n dash = len * ( len - 1 ) / 2; slope仁 zeros( n dash, 1 ); m=1;for k = 1:le n-1.& -Xj=Median u -j_j ： ifor j = k+1:le n.slope1(m) = ( x(j) - x(k) ) / ( j - k ); m = m + 1;en d;en d;slope= media n( slope1 );2、M-K突变检验clc;y=data(1:50,6);% 输入数据 n=len gth(y);Sk=zeros(size(y));UFk=zeros(size(y));s1=0;for i=2: nk% Sk 八 ri (k =2,3, ,n)imfor j=1:iif y(i)>y(j)s仁 s1+1;elses1=s1+0;1,Xj Xjn% 1°,x- xj ( j=1,2，…，i)end;end;Sk(i)=s1;E=i*(i-1)/4;Var=i*(i-1)*(2*i+5)/72;UFk(i)=(Sk(i)-E)/sqrt(Var);en d;E(Sk)n(n _1)4__% UFk[Sk - E (Sk)]XVar (Sk)(k = 1,2 - , n)y2=zeros(size(y));Sk2=zeros(size(y));UBk=zeros(size(y));s2=0;for i=1: ny2(i)=y( n-i+1);en d;for i=2: nfor j=1:iif y2(i)>y2(j)s2=s2+1;elses2=s2+0;en d;en d;Sk2(i)=s2;E=i*(i-1)/4; % Sk2(i)的均值% Var(s k)=n(n -1)(2n 5)72Var=i*(i-1)*(2*i+5)/72; % Sk2(i) 的方差UBk(i)=0-(Sk2(i)-E)/sqrt(Var);en d;UBk2=zeros(size(y));for i=1: nUBk2(i)=UBk (n-i+1);en d;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');。