Matlab处理气象数据（八）数据的逐点变化趋势

2024-07-15 02:43| 来源: 网络整理| 查看: 265

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

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

M-K趋势检验原理

定义检验统计量 S S S： S = ∑ i = 2 n ∑ j = 1 i − 1 s i g n ( X i − X j ) S=\sum_{i=2}^n\sum_{j=1}^{i-1}sign(X_i-X_j) S=i=2∑nj=1∑i−1sign(Xi−Xj) 其中， s i g n ( ) sign() sign()为符号函数。当 X i − X j X_i-X_j Xi−Xj小于、等于或大于0时， s i g n ( X i − X j ) sign(X_i-X_j) sign(Xi−Xj)分别为-1、0或1；M-K统计量公式 S S 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 = ( 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 ) 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 ) ] );

得到的结果：

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 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 % 观测数据逐点变化趋势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

即可得到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

最高、最低温度的与之类似，得到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); % 退耕前的逐点趋势 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 % 退耕后的逐点趋势 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

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

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

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

【本文地址】

公司简介

联系我们