Matlab处理气象数据（八）数据的逐点变化趋势_mann-kendall趋势检验matlab代码

前面所做的都是从时间尺度上研究的变化趋势，而从空间尺度上分析，能够更加直观地看出温度变化的地理位置。

M-K（Mann-Kendall）是世界气象组织推荐并被广泛用于实际研究的非参数检验方法，是时间序列趋势分析方法之一。它不要求被分析样本遵从一定分布，同时也不受其它异常值的干扰，对于非正统分布的气象数据，M-K秩次相关检验具有更加突出的适用性。

M-K趋势检验原理

定义检验统计量$S$：
$$S=\sum _{i=2}^n\sum _{j=1}^{i-1}sign(X_i-X_j)$$
其中， $sign()$为符号函数。当$X_i-X_j$小于、等于或大于0时，$sign(X_i-X_j)$分别为-1、0或1；M-K统计量公式$S$大于、等于或小于0时分别为：
Z = { ( S − 1 ) / n ( n − 1 ) ( 2 n + 5 ) / 18 S > 0 0 S = 0 ( S + 1 ) / n ( n − 1 ) ( 2 n + 5 ) / 18 S < 0 Z=

$$\begin{cases} {(S-1)}/{\sqrt{n(n-1)(2n+5)/18}} \qquad & S \gt 0 \\ 0 \qquad & S = 0 \\ {(S+1)}/{\sqrt{n(n-1)(2n+5)/18}} \qquad & S \lt 0 \\ \end{cases}$$ Z=⎩⎪⎨⎪⎧​(S−1)/n(n−1)(2n+5)/18 ​0(S+1)/n(n−1)(2n+5)/18 ​​S>0S=0S<0​
$Z$为正值表示增加趋势，负值表示减少趋势。$Z$的绝对值在大于等于1.28、1.64和2.32时表示分别通过了信度90%、95%和99%的显著性检验。

M-K趋势检验Matlab代码

% function [ z, sl, lcl, ucl ] = mk( y )
y=DT;
n = length( y );
dt=1;

% 计算统计量
s = 0;
for k = 1:n-1,
    for j = k+1:n,
        s = s + sign( y(j) - y(k) );
    end;
end;

% 方差
v = ( n * ( n - 1 ) * ( 2 * n + 5 ) ) / 18;

% 检验统计量
if s == 0,
	z = 0;
	elseif s > 0,
	z = ( s - 1 ) / sqrt( v );
else
	z = ( s + 1 ) / sqrt( v );
end; 

% 标准值
nor = 1.64;
% 结果
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:' );

% 计算斜率
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;

% 估计
sl = median( s );
disp( [ ' Slope Estimate = ' num2str( sl ) ] );
% 方差 ( 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 ) ] );
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

得到的结果：

DT:
 n = 35
 Mean Value = 1.3474
 Z statistic = 5.4249
 Upward trend detected
Sens Nonparametric Estimator:
 Slope Estimate = 0.018411
 Lower Confidence Limit = 0.014871
 Upper Confidence Limit = 0.021697



Tem1:
 n = 35
 Mean Value = 5.2637
 Z statistic = 2.8687
 Upward trend detected
Sens Nonparametric Estimator:
 Slope Estimate = 0.020431
 Lower Confidence Limit = 0.0088048
 Upper Confidence Limit = 0.030395

Tem2:
 n = 35
 Mean Value = 6.6111
 Z statistic = 4.8569
 Upward trend detected
Sens Nonparametric Estimator:
 Slope Estimate = 0.036645
 Lower Confidence Limit = 0.027685
 Upper Confidence Limit = 0.046547
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

8.1 分别对NCEP和观测数据进行逐点趋势分析

%NCEP数据逐点变化趋势s1
% 导入数据T1
Y1=reshape(T1,[72,128,12,35]);
Y1=nanmean(Y1,3);
Y1=squeeze(Y1);

for j=1:128
	for i=1:72
		if isnan(Y1)==1
			continue
		else
			y=Y1(i,j,:);
			n = 35;
			dt=1;
		
			% 计算斜率
			ndash = n * ( n - 1 ) / 2;
			s = zeros( ndash,1 );
			s=nan;
			r = 1;
			for p = 1:n-1
				for q = p+1:n
					s(r) = ( y(q) - y(p) ) / ( q - p ) / dt;
					r = r + 1;
				end
			end 
			s1(i,j) = median( s );
		end
	end
end
pcolor(s1);
shading flat
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

---

% 观测数据逐点变化趋势s2
% 导入数据T2
Y2=reshape(T2,[72,128,12,35]);
Y2=nanmean(Y2,3);
Y2=squeeze(Y2);

for j=1:128
	for i=1:72
		if isnan(Y2)==1
			continue
		else
			y=Y2(i,j,:);
			n = 35;
			dt=1;
		
			% 计算斜率
			ndash = n * ( n - 1 ) / 2;
			s = zeros( ndash,1 );
			s=nan;
			r = 1;
			for p = 1:n-1
				for q = p+1:n
					s(r) = ( y(q) - y(p) ) / ( q - p ) / dt;
					r = r + 1;
				end
			end
			s2(i,j) = median( s );
		end
	end
end
pcolor(s2);
shading flat
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

即可得到NCEP数据和观测数据的平均温度逐点变化趋势，用同样的方法可以得到最高温度、最低温度的变化趋势。（图略）

8.2 平均、最高、最低温度差值的趋势

% 两套数据年平均值逐点变化趋势sy
% 导入数据T1、T2
Y1=reshape(T1,[72,128,12,35]);
Y1=nanmean(Y1,3);
Y1=squeeze(Y1);
Y2=reshape(T2,[72,128,12,35]);
Y2=nanmean(Y2,3);
Y2=squeeze(Y2);

for j=1:128
	for i=1:72
		Y=Y2-Y1;
		y=Y(i,j,:);
		n = 35;
		dt=1;
		% 计算斜率
		ndash = n * ( n - 1 ) / 2;
		s = zeros( ndash,1 );
		s=nan;
		r = 1;
		for p = 1:n-1
			for q = p+1:n
				s(r) = ( y(q) - y(p) ) / ( q - p ) / dt;
				r = r + 1;
			end
		end
		sy(i,j) = median( s );
	end
end
pcolor(sy);
shading flat
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

最高、最低温度的与之类似，得到NCEP数据和观测数据平均温度、最高温度、最低温度差值的逐点变化趋势。（图略）

8.3 退耕还林前后的温度差值对比

仍然以1999年为界限，1981-1999年为退耕前，2000-2013年为退耕后。
保存之前由T1、T2计算得到的Y1、Y2两套数据，以便后面直接调用，其格式为72*128*35。

% 导入数据Y1、Y2
Y=Y2-Y1;
Y=Y(:,:,3:35);
YA=Y(:,:,1:19);
YB=Y(:,:,20:33);
1
2
3
4
5

% 退耕前的逐点趋势
for j=1:128
	for i=1:72
		y=YA(i,j,:);
		n = 19;
		dt=1;
		% 计算斜率
		ndash = n * ( n - 1 ) / 2;
		s = zeros( ndash,1 );
		s=nan;
		r = 1;
		for p = 1:n-1
			for q = p+1:n
				s(r) = ( y(q) - y(p) ) / ( q - p ) / dt;
				r = r + 1;
			end
		end
		sya(i,j) = median( s );
	end
end
pcolor(sya);
shading flat
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

% 退耕后的逐点趋势
for j=1:128
	for i=1:72
		y=YB(i,j,:);
		n = 14;
		dt=1;
		% 计算斜率
		ndash = n * ( n - 1 ) / 2;
		s = zeros( ndash,1 );
		s=nan;
		r = 1;
		for p = 1:n-1
			for q = p+1:n
				s(r) = ( y(q) - y(p) ) / ( q - p ) / dt;
				r = r + 1;
			end
		end
		syb(i,j) = median( s );
	end
end
pcolor(syb);
shading flat
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

得到NCEP数据和观测数据平均温度差值1981-1999年和2000-2013年的逐点变化趋势。（图略）

最高、最低温度差值退耕前后的计算与之类似。

---

相关链接：
Matlab处理气象数据——目录