Matlab在气象专业教学中的应用——气象要素的M-K检验突变分析

胡 琦， 马雪晴， 胡莉婷， 王雅婧， 徐 琳， 潘学标

(中国农业大学 资源与环境学院，北京 100193)

0 引 言

气候变化已经成为全球关注的焦点问题，正在产生的影响严重威胁着自然界和人类的安全[1]，气候变化议题也受到各国科学家的持续关注，相关研究不断深化。大量事实表明气候不是缓慢变化的，而是从一种稳定状态跳跃式地转变到另一种稳定状态，称为气候突变现象[2]，其理论研究是近代气候学一个新兴的研究领域[3]。M-K检验法是世界气象组织推荐的用于提取序列变化趋势的有效工具，目前已经被广泛用于气候参数和水文序列的分析中[4-6]。M-K法是中国农业大学应用气象学专业的必修课“应用气候学实习”课程中讲授的一个重要的气候分析方法。M-K方法的原理和计算过程较为繁琐，若利用常规的计算手段和软件(如EXCEL等)效率较低，无法满足时代信息化教学和学生的实际业务应用。因此，将计算机编程的方法应用于气象专业的教学教育之中，在教学活动中应用是现代信息技术，对提高本科教学质量至关重要。

Matlab软件是一门集算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言[7]。目前，Matlab软件已经广泛应用到各个教学领域，如数学[8]、物理[9]、生化[10]、医学[11]、音乐[12]、金融学[13]等。关于Matlab在气象专业教学中应用的研究较为罕见，目前尚未有关于Matlab软件实现M-K检验的教学文献，为此，本文选取了某气象站点多年的温度观测数据，介绍了应用Matlab软件实现M-K法中统计量的自动计算，并如何基于Matlab-GUI实现界面操作和自动制图以提高处理数据的效率。

1 方法原理

M-K检验法最初由曼(H.B.Mann)和肯德尔(M.G.Kendall)提出了原理并发展了这一方法，是世界气象组织推荐的用于提取序列变化趋势的有效工具[14]。M-K检验法不受个别异常值的干扰，能够客观反映时间序列趋势，目前已经被广泛用于气候参数和水文序列的分析中。M-K法可以根据输出的两个序列(UF和UB)明确突变的时段和区域。具体计算方法及参数如下[15]：

对于具有n个样本量的时间序列X，构造一秩序列：

(1)

可见，秩序列sk是第i时刻数值大于j时刻数值个数的累计数,易知k=1时,s1=0。

在时间序列随机独立的假定下，定义统计量:

(2)

式中,UF1=0，E(sk)，Var(sk)是累计数sk的均值和方差，在x1,x2,…,xn相互独立，且有相同连续分布时，它们可由下式算出：

(3)

UF为标准正态分布，它是按时间序列x顺序x1,x2,…,xn计算出的统计量序列。按时间序列x逆序xn,xn-1,…,x1，再重复上述过程，构造逆序列UB。

若UF值大于0，则表明序列呈上升趋势，小于0则表明呈下降趋势。当它们超过临界置信水平直线时(检验置信水平α=0.05时，置信水平线为±1.96)，表明上升或下降趋势显著，超过临界线的范围确定为出现突变的时间区域。如果UF和UB两条曲线出现交点，且交点在临界线之间，那么交点对应的时刻便是突变开始的时间。

2 案例分析

表1为某气象站1961～2015年年平均温度数据序列，试分析该地区的年平均温度是否有气候突变发生。

表1 某气象站1961～2015年年平均温度

2.1 数据输入

将数据存储在txt文档中，存储的txt数据无表头，第1列为年份，第2列为年平均温度，空格分隔符。数据可放在任意路径下(即不提前设定数据存储路径)，可以通过弹出的窗口进行选择文件。

打开Matlab软件，在命令窗口中输入代码：

clear %清空环境变量

clc %清空窗口内容

[filename filepath]=uigetfile(′*.*′,′请选择文件′);

%filename为文件名，filepath为文件路径

if isequal(filename,0)

msgbox(′您没有正确选择文件夹′);

return;

end

Data=textread([filepath, filename]); %读取文件，并赋值给矩阵Data

2.2 计算UF和UB序列值

2.2.1 计算正序列UF值

计算UF正序列值的难点在于秩序列sk的计算，根据式(1)，秩序列sk是第i时刻数值大于j时刻数值个数的累计数，可以用双重for循环计算。

接着在命令窗口中输入代码：

y=Data(:,2);%平均温度序列

Sk=zeros(size(y)); %定义累计量序列Sk，长度=y，初始值=0,Sk(1)=0

UFk=zeros(size(y)); % 定义统计量UFk，长度=y，初始值=0, UFk(1)=0

s = 0; % 定义Sk序列的元素s

for i=2:length(y)

for j=1:i

if y(i)>y(j)

s=s+1;

else

s=s+0;

end

end

Sk(i)=s;

E=i*(i-1)/4; % Sk(i)的均值,见式(3)

Var=i*(i-1)*(2*i+5)/72; % Sk(i)方差,见式(3)

UFk(i)=(Sk(i)-E)/sqrt(Var);%正序列UF值,见式(2)

end

2.2.2 计算逆序列UB值

与计算UF正序列值类似，逆序列UB值的计算代码如下：

Sk2=zeros(size(y)); % 定义逆序累计量序列Sk2，长度=y，初始值=0,Sk(2)=0

UBk=zeros(size(y)); % 定义逆序统计量UBk，长度=y，初始值=0，UBk(1)=0

s=0;

y2=flipud(y); % 按时间序列逆转平均温度序列

for i=2:length(y2)

for j=1:i

if y2(i)>y2(j)

s=s+1;

else

s=s+0;

end

end

Sk2(i)=s;

E=i*(i-1)/4; %均值

Var=i*(i-1)*(2*i+5)/72; %方差

UBk(i)=0-(Sk2(i)-E)/sqrt(Var);

end

UBk2=flipud(UBk); %逆序列UB值

2.3 画图及保存图片

Matlab软件中plot() 函数能够针对向量或矩阵的列来绘制曲线，是绘制二维图形的最基本函数。本例中的M-K检验图须调用函数plot(x,y)来绘制M-K检验值(UF线、UB线)和0.05显著性水平下的置信水平直线，其中以年份为横坐标值，M-K检验值为纵坐标值绘制曲线。

绘制M-K检验图的代码如下：

x=Data(:,1);%年份序列

n=length(x);%年份序列的长度

figure %做图

set(gcf,′unit′,′centimeters′,′position′,[3 5 12 6]) %设置图形的位置及大小

set(gca,′Position′,[.15 .18 .8 .78]);%设置图片比例大小

plot(x,UFk,′r-′,′linewidth′,1.5);%画UF线

hold on

plot(x,UBk2,′b-.′,′linewidth′,1.5);%画UB线

plot(x,1.96*ones(n,1),′k:′,′linewidth′,1);

axis([min(x),max(x),-5,5]);%设置X轴范围和间距

legend(′UF统计量′,′UB统计量′,′0.05显著水平′);%设置图例

legend(′Orientation′,′horizon′,′location′,′south′);%图例横排

legend(′boxoff′) %设置图例为无边框

set(gca, ′FontSize′, 8) %设置图例和坐标轴字体大小

xlabel(′年 Year′,′FontName′,′TimesNewRoman′,′FontSize′,10);%X轴标题

ylabel(′统计量 MK Value′,′FontName′,′TimesNewRoman′,′Fontsize′,10);%Y轴标题

hold on

plot(x,-1.96*ones(n,1),′k:′,′linewidth′,1);

plot(x,0*ones(n,1),′k-.′,′linewidth′,1);%图片绘制完毕

图片绘制完毕后保存图片，图片保存的位置、格式和名称可以通过对话框自行选择，图片可以保存为.jpg、.bmp、.png、.gif等格式。代码如下：

[FileName,PathName] = uiputfile({′*.jpg′,′JPEG(*.jpg)′;…

′*.bmp′,′Bitmap(*.bmp)′;…

′*.gif′,′GIF(*.gif)′;…

′*.png′,′png file(*.png)′;…

′*.*′, ′All Files (*.*)′},…

′Save Picture′,′图片′);%选择图片保存的位置，名称和格式

if FileName==0

msgbox(′No Filename Selected′);

return;

else

saveas(gca,[PathName,FileName]);%保存图片

end

图片结果如图1所示。由图可知，该地区1961～2015年气温呈显著上升趋势，UF和UB统计量有交点且交点在置信直线范围之间，表明气温在1989年前后发生了突变。

图1 某地区平均温度长时间序列M-K统计量变化图

2.4 基于Matlab-GUI的界面操作

基于Matlab-GUI制作了应用气候学实习-突变分析和检验软件v1.0，该软件能够支持.txt、.excel和.mat格式的文件，同时支持单个文件输入和批量文件输入。为与课程相匹配，设置了原理方法介绍模块，方便学生学习和理解M-K检验的原理和计算方法(见图2)。

图2 应用气候学实习-突变分析和检验软件主界面

数据输入成功后，软件“数据信息”面板中能够显示输入数据的名称、格式、存储路径和要素值主要信息等内容。点击“计算”按钮，“结果面板”会显示M-K检验的结果，包括突变点个数、突变年份和显著性(图3(a))。批量文件输入时，“结果面板”会显示每一个文件的名称及相应的M-K检验的结果，如图3(b)所示。

(a)

(b)

图3 应用气候学实习-突变分析和检验软件计算界面

对于单个文件输入时，软件提供了查看“详细参数”和“绘图”的功能，详细参数包括了M-K检验计算时的中间参数sk,E(sk)，Var(sk),UFk,UBk等，如图4(a)所示；绘图功能能够显示M-K检验图和添加突变点信息(图4(b))。

(a)

(b)

图4 应用气候学实习-突变分析和检验软件中查看详细参数和绘图功能

3 结 语

M-K突变检验是世界气象组织推荐提取序列变化趋势的有效工具，也是“应用气候学实习”课程中重要的教学内容。本文以某气象站1961—2015年年平均温度为例，详述了使用Matlab软件计算MK突变检验的UF和UB统计量，并基于Matlab/GUI实现了界面操作和自动制图。MK突变检验计算的难点在于秩序列sk的计算，Matlab编程时可以通过双重for循环来实现。通过计算程序能有效提高处理数据的速度，对学生综合科学研究素质和实际应用水平的提高有较好的帮助。