用Matlab进行MK检验.doc
《用Matlab进行MK检验.doc》由会员分享，可在线阅读，更多相关《用Matlab进行MK检验.doc（5页珍藏版）》请在冰点文库上搜索。
用Matlab进行Mann-Kendall分析
一、M-K趋势分析法与M-K突变检验的原理介绍
1、Mann-Kendall趋势分析法
非参数检验，又称任意分布检验，它不对变量的分布做严格假定，检验不针对特定的参数，而是模糊地对变量分布的中心位置或分布状态做检验，由于其不对总体分布做严格假定，因而适用性强[12]。
因此，本文采用非参数的Mann-Kendall检验法对昌马河流域近50年的气候水文要素时间序列显著性检验，定量反映变化趋势的显著性。
计算公式如下：
（1.1）
在公式中，
（1.2）
（1.3）
（1.4）
式中：
xk，xi为连续的气候、水文等数据序列，n为数据集合的总长度，t为每个单位的宽度，Σ表示所有单位的总和。
衡量趋势大小的指标为：
（1.5）
式中：
1Mann-Kendall检验如下所示：
零假设H0:
β=0
当，拒绝H0假设。
式中：
为标准正态方差，α为显著性检验水平。
2、Mann-kendall突变检验
气候系统变化是一个不稳定且不连续的变化过程，而检验其变化的常用方法之一就是Mann-kendall突变检验方法[13]，该方法对于变化要素从一个相对稳定状态变化到另一个状态的变化检验非常有效。
且广泛应用于水文，气候，化学，矿物成分检验等各个方面。
Mann-kendall突变检验方法如下：
对于具有n个样本量的时间序列x，构造一个秩序列：
（1.6）
其中
（）（1.7）
可见，秩序列是第时刻数值大于j时刻数值个数的累计数。
在时间序列随机独立的假定下，定义统计量：
（1.8）
其中UF1=0，，是累计数的均值和方差，在相互独立，且有相同连续分布时，它们可由下式算出：
（1.9）
（1.10）
为标准正态分布，它是按时间序列x顺序计算出的统计量序列，给定显著性水平a，查正态分布表，若，则表明序列存在明显的趋势变化。
把此方法引用到时间序列的逆序序列中，按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%
s=0;
len=size（x,2）;
form=1:
len-1
forn=m+1:
len
ifx（n）>x（m）
s=s+1;
elseifx（n）==x（m）
s=s+0;
else
s=s-1;
end
end
end
%计算vars和e%
vars=len*（len-1）*（2*len+5）/18;
%计算zc%
ifs>0
zc=（s-1）/sqrt（vars）;
else
zc=（s+1）/sqrt（vars）;
end
%计算za
za=var（x）;
sign=0;
zc1=abs（zc）;
ifzc1>=za
sign=1;
else
sign=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）;%
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%
forj=1:
i
ify（i）>y（j）
s1=s1+1;
else
s1=s1+0;%（）
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）;%
end;
y2=zeros（size（y））;
Sk2=zeros（size（y））;
UBk=zeros（size（y））;
s2=0;
fori=1:
n
y2（i）=y（n-i+1）;
end;
fori=2:
n
forj=1:
i
ify2（i）>y2（j）
s2=s2+1;
else
s2=s2+0;
end;
end;
Sk2（i）=s2;
E=i*（i-1）/4;%Sk2（i）的均值
Var=i*（i-1）*（2*i+5）/72;%Sk2（i）的方差%
UBk（i）=0-（Sk2（i）-E）/sqrt（Var）;
end;
UBk2=zeros（size（y））;
fori=1:
n
UBk2（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'）;
升级会员 