1、Mann-Kendall趋势分析法
非参数检验，又称任意分布检验，它不对变量的分布做严格假定，检验不针对特定的参
数，而是模糊地对变量分布的中心位置或分布状态做检验，由于其不对总体分布做严格假定，
因而适用性强
[12]
。因此，本文采用非参数的Mann-Kendall检验法对昌马河流域近50年的气
候水文要素时间序列显著性检验，定量反映变化趋势的显著性。世界气象组织建议使用M-K趋势检验法
进行趋势分析。
Storch等建议在进行M-K检验之前对时间序列进行预白化处理，以消除序列中的自相关性，减少自相关性对M-K
检验结果的影响。预白化方法设有时间序列y1,y2,y3,...,yn，n为样本量，首先计算滞后1（一阶自相关系数）的序列
自相关系数c，若c<0.1，则此时间序列可以直接应用于M-K分析；若c>0.1，则序列须经预白化处理得到新的时间序列，
即y2-cy1，y3-cy2,y4-cy3，...,yn-cyn-1，然后将M-K用于处理以后的时间序列趋势分析。
有一种笨方法~把该时间序列（假设为x）的最后一个数去掉，再在序列最前面加个0，再命名为y
然后cor（x，y）
所得的结果就是一阶自相关系数
function [slope,zc,za,sign]=MannKendall(x)
%S的计算
s=0;
len=size(x,2);
for m=1:len-1
for n=m+1:len
if x(n)>x(m)
s=s+1;
elseif x(n)==x(m)
s=s+0;
else
s=s-1;
end
end
end
%计算vars和e%
vars=len*(len-1)*(2*len+5)/18;
%计算zc
if s>0
zc=(s-1)/sqrt(vars);
else
zc=(s+1)/sqrt(vars);
end
%计算za
za=var(x);
sign=0;
zc1=abs(zc);
if zc1 >= za
sign=1;
else
sign=0;
end
%计算倾斜度
ndash = len * ( len - 1 ) / 2;
slope1= zeros( ndash, 1 );
m=1;
for k = 1:len-1,
for j = k+1:len,
slope1(m) = ( x(j) - x(k) ) / ( j - k )
m = m + 1;
end;
end;
slope= median( slope1 );
function [z,s1,lc1,uc1]=trendMK(y)
A=txtread('气温年序列.txt');
y=A(:,2);
n = length( y );
dt=1;
% 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.64;
% 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 ) ] );
