Matlab在氣象專(zhuān)業(yè)教學(xué)中的應(yīng)用
——?dú)庀笠氐腗-K檢驗(yàn)突變分析

胡 琦， 馬雪晴， 胡莉婷， 王雅婧， 徐 琳， 潘學(xué)標(biāo)

(中國(guó)農(nóng)業(yè)大學(xué) 資源與環(huán)境學(xué)院，北京 100193)

0 引 言

氣候變化已經(jīng)成為全球關(guān)注的焦點(diǎn)問(wèn)題，正在產(chǎn)生的影響嚴(yán)重威脅著自然界和人類(lèi)的安全[1]，氣候變化議題也受到各國(guó)科學(xué)家的持續(xù)關(guān)注，相關(guān)研究不斷深化。大量事實(shí)表明氣候不是緩慢變化的，而是從一種穩(wěn)定狀態(tài)跳躍式地轉(zhuǎn)變到另一種穩(wěn)定狀態(tài)，稱(chēng)為氣候突變現(xiàn)象[2]，其理論研究是近代氣候?qū)W一個(gè)新興的研究領(lǐng)域[3]。M-K檢驗(yàn)法是世界氣象組織推薦的用于提取序列變化趨勢(shì)的有效工具，目前已經(jīng)被廣泛用于氣候參數(shù)和水文序列的分析中[4-6]。M-K法是中國(guó)農(nóng)業(yè)大學(xué)應(yīng)用氣象學(xué)專(zhuān)業(yè)的必修課“應(yīng)用氣候?qū)W實(shí)習(xí)”課程中講授的一個(gè)重要的氣候分析方法。M-K方法的原理和計(jì)算過(guò)程較為繁瑣，若利用常規(guī)的計(jì)算手段和軟件(如EXCEL等)效率較低，無(wú)法滿足時(shí)代信息化教學(xué)和學(xué)生的實(shí)際業(yè)務(wù)應(yīng)用。因此，將計(jì)算機(jī)編程的方法應(yīng)用于氣象專(zhuān)業(yè)的教學(xué)教育之中，在教學(xué)活動(dòng)中應(yīng)用是現(xiàn)代信息技術(shù)，對(duì)提高本科教學(xué)質(zhì)量至關(guān)重要。

Matlab軟件是一門(mén)集算法開(kāi)發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計(jì)算的高級(jí)技術(shù)計(jì)算語(yǔ)言[7]。目前，Matlab軟件已經(jīng)廣泛應(yīng)用到各個(gè)教學(xué)領(lǐng)域，如數(shù)學(xué)[8]、物理[9]、生化[10]、醫(yī)學(xué)[11]、音樂(lè)[12]、金融學(xué)[13]等。關(guān)于Matlab在氣象專(zhuān)業(yè)教學(xué)中應(yīng)用的研究較為罕見(jiàn)，目前尚未有關(guān)于Matlab軟件實(shí)現(xiàn)M-K檢驗(yàn)的教學(xué)文獻(xiàn)，為此，本文選取了某氣象站點(diǎn)多年的溫度觀測(cè)數(shù)據(jù)，介紹了應(yīng)用Matlab軟件實(shí)現(xiàn)M-K法中統(tǒng)計(jì)量的自動(dòng)計(jì)算，并如何基于Matlab-GUI實(shí)現(xiàn)界面操作和自動(dòng)制圖以提高處理數(shù)據(jù)的效率。

1 方法原理

M-K檢驗(yàn)法最初由曼(H.B.Mann)和肯德?tīng)?M.G.Kendall)提出了原理并發(fā)展了這一方法，是世界氣象組織推薦的用于提取序列變化趨勢(shì)的有效工具[14]。M-K檢驗(yàn)法不受個(gè)別異常值的干擾，能夠客觀反映時(shí)間序列趨勢(shì)，目前已經(jīng)被廣泛用于氣候參數(shù)和水文序列的分析中。M-K法可以根據(jù)輸出的兩個(gè)序列(UF和UB)明確突變的時(shí)段和區(qū)域。具體計(jì)算方法及參數(shù)如下[15]：

對(duì)于具有n個(gè)樣本量的時(shí)間序列X，構(gòu)造一秩序列：

(1)

可見(jiàn)，秩序列sk是第i時(shí)刻數(shù)值大于j時(shí)刻數(shù)值個(gè)數(shù)的累計(jì)數(shù),易知k=1時(shí),s1=0。

在時(shí)間序列隨機(jī)獨(dú)立的假定下，定義統(tǒng)計(jì)量:

(2)

式中,UF1=0，E(sk)，Var(sk)是累計(jì)數(shù)sk的均值和方差，在x1,x2,…,xn相互獨(dú)立，且有相同連續(xù)分布時(shí)，它們可由下式算出：

(3)

UF為標(biāo)準(zhǔn)正態(tài)分布，它是按時(shí)間序列x順序x1,x2,…,xn計(jì)算出的統(tǒng)計(jì)量序列。按時(shí)間序列x逆序xn,xn-1,…,x1，再重復(fù)上述過(guò)程，構(gòu)造逆序列UB。

若UF值大于0，則表明序列呈上升趨勢(shì)，小于0則表明呈下降趨勢(shì)。當(dāng)它們超過(guò)臨界置信水平直線時(shí)(檢驗(yàn)置信水平α=0.05時(shí)，置信水平線為±1.96)，表明上升或下降趨勢(shì)顯著，超過(guò)臨界線的范圍確定為出現(xiàn)突變的時(shí)間區(qū)域。如果UF和UB兩條曲線出現(xiàn)交點(diǎn)，且交點(diǎn)在臨界線之間，那么交點(diǎn)對(duì)應(yīng)的時(shí)刻便是突變開(kāi)始的時(shí)間。

2 案例分析

表1為某氣象站1961～2015年年平均溫度數(shù)據(jù)序列，試分析該地區(qū)的年平均溫度是否有氣候突變發(fā)生。

表1 某氣象站1961～2015年年平均溫度

2.1 數(shù)據(jù)輸入

將數(shù)據(jù)存儲(chǔ)在txt文檔中，存儲(chǔ)的txt數(shù)據(jù)無(wú)表頭，第1列為年份，第2列為年平均溫度，空格分隔符。數(shù)據(jù)可放在任意路徑下(即不提前設(shè)定數(shù)據(jù)存儲(chǔ)路徑)，可以通過(guò)彈出的窗口進(jìn)行選擇文件。

打開(kāi)Matlab軟件，在命令窗口中輸入代碼：

clear %清空環(huán)境變量

clc %清空窗口內(nèi)容

[filename filepath]=uigetfile(′*.*′,′請(qǐng)選擇文件′);

%filename為文件名，filepath為文件路徑

if isequal(filename,0)

msgbox(′您沒(méi)有正確選擇文件夾′);

return;

end

Data=textread([filepath, filename]); %讀取文件，并賦值給矩陣Data

2.2 計(jì)算UF和UB序列值

2.2.1 計(jì)算正序列UF值

計(jì)算UF正序列值的難點(diǎn)在于秩序列sk的計(jì)算，根據(jù)式(1)，秩序列sk是第i時(shí)刻數(shù)值大于j時(shí)刻數(shù)值個(gè)數(shù)的累計(jì)數(shù)，可以用雙重for循環(huán)計(jì)算。

接著在命令窗口中輸入代碼：

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

Sk=zeros(size(y)); %定義累計(jì)量序列Sk，長(zhǎng)度=y，初始值=0,Sk(1)=0

UFk=zeros(size(y)); % 定義統(tǒng)計(jì)量UFk，長(zhǎng)度=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)的均值,見(jiàn)式(3)

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

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

end

2.2.2 計(jì)算逆序列UB值

與計(jì)算UF正序列值類(lèi)似，逆序列UB值的計(jì)算代碼如下：

Sk2=zeros(size(y)); % 定義逆序累計(jì)量序列Sk2，長(zhǎng)度=y，初始值=0,Sk(2)=0

UBk=zeros(size(y)); % 定義逆序統(tǒng)計(jì)量UBk，長(zhǎng)度=y，初始值=0，UBk(1)=0

s=0;

y2=flipud(y); % 按時(shí)間序列逆轉(zhuǎn)平均溫度序列

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 畫(huà)圖及保存圖片

Matlab軟件中plot() 函數(shù)能夠針對(duì)向量或矩陣的列來(lái)繪制曲線，是繪制二維圖形的最基本函數(shù)。本例中的M-K檢驗(yàn)圖須調(diào)用函數(shù)plot(x,y)來(lái)繪制M-K檢驗(yàn)值(UF線、UB線)和0.05顯著性水平下的置信水平直線，其中以年份為橫坐標(biāo)值，M-K檢驗(yàn)值為縱坐標(biāo)值繪制曲線。

繪制M-K檢驗(yàn)圖的代碼如下：

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

n=length(x);%年份序列的長(zhǎng)度

figure %做圖

set(gcf,′unit′,′centimeters′,′position′,[3 5 12 6]) %設(shè)置圖形的位置及大小

set(gca,′Position′,[.15 .18 .8 .78]);%設(shè)置圖片比例大小

plot(x,UFk,′r-′,′linewidth′,1.5);%畫(huà)UF線

hold on

plot(x,UBk2,′b-.′,′linewidth′,1.5);%畫(huà)UB線

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

axis([min(x),max(x),-5,5]);%設(shè)置X軸范圍和間距

legend(′UF統(tǒng)計(jì)量′,′UB統(tǒng)計(jì)量′,′0.05顯著水平′);%設(shè)置圖例

legend(′Orientation′,′horizon′,′location′,′south′);%圖例橫排

legend(′boxoff′) %設(shè)置圖例為無(wú)邊框

set(gca, ′FontSize′, 8) %設(shè)置圖例和坐標(biāo)軸字體大小

xlabel(′年 Year′,′FontName′,′TimesNewRoman′,′FontSize′,10);%X軸標(biāo)題

ylabel(′統(tǒng)計(jì)量 MK Value′,′FontName′,′TimesNewRoman′,′Fontsize′,10);%Y軸標(biāo)題

hold on

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

plot(x,0*ones(n,1),′k-.′,′linewidth′,1);%圖片繪制完畢

圖片繪制完畢后保存圖片，圖片保存的位置、格式和名稱(chēng)可以通過(guò)對(duì)話框自行選擇，圖片可以保存為.jpg、.bmp、.png、.gif等格式。代碼如下：

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

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

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

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

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

′Save Picture′,′圖片′);%選擇圖片保存的位置，名稱(chēng)和格式

if FileName==0

msgbox(′N(xiāo)o Filename Selected′);

return;

else

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

end

圖片結(jié)果如圖1所示。由圖可知，該地區(qū)1961～2015年氣溫呈顯著上升趨勢(shì)，UF和UB統(tǒng)計(jì)量有交點(diǎn)且交點(diǎn)在置信直線范圍之間，表明氣溫在1989年前后發(fā)生了突變。

圖1 某地區(qū)平均溫度長(zhǎng)時(shí)間序列M-K統(tǒng)計(jì)量變化圖

2.4 基于Matlab-GUI的界面操作

基于Matlab-GUI制作了應(yīng)用氣候?qū)W實(shí)習(xí)-突變分析和檢驗(yàn)軟件v1.0，該軟件能夠支持.txt、.excel和.mat格式的文件，同時(shí)支持單個(gè)文件輸入和批量文件輸入。為與課程相匹配，設(shè)置了原理方法介紹模塊，方便學(xué)生學(xué)習(xí)和理解M-K檢驗(yàn)的原理和計(jì)算方法(見(jiàn)圖2)。

圖2 應(yīng)用氣候?qū)W實(shí)習(xí)-突變分析和檢驗(yàn)軟件主界面

數(shù)據(jù)輸入成功后，軟件“數(shù)據(jù)信息”面板中能夠顯示輸入數(shù)據(jù)的名稱(chēng)、格式、存儲(chǔ)路徑和要素值主要信息等內(nèi)容。點(diǎn)擊“計(jì)算”按鈕，“結(jié)果面板”會(huì)顯示M-K檢驗(yàn)的結(jié)果，包括突變點(diǎn)個(gè)數(shù)、突變年份和顯著性(圖3(a))。批量文件輸入時(shí)，“結(jié)果面板”會(huì)顯示每一個(gè)文件的名稱(chēng)及相應(yīng)的M-K檢驗(yàn)的結(jié)果，如圖3(b)所示。

(a)

(b)

圖3 應(yīng)用氣候?qū)W實(shí)習(xí)-突變分析和檢驗(yàn)軟件計(jì)算界面

對(duì)于單個(gè)文件輸入時(shí)，軟件提供了查看“詳細(xì)參數(shù)”和“繪圖”的功能，詳細(xì)參數(shù)包括了M-K檢驗(yàn)計(jì)算時(shí)的中間參數(shù)sk,E(sk)，Var(sk),UFk,UBk等，如圖4(a)所示；繪圖功能能夠顯示M-K檢驗(yàn)圖和添加突變點(diǎn)信息(圖4(b))。

(a)

(b)

圖4 應(yīng)用氣候?qū)W實(shí)習(xí)-突變分析和檢驗(yàn)軟件中查看詳細(xì)參數(shù)和繪圖功能

3 結(jié) 語(yǔ)

M-K突變檢驗(yàn)是世界氣象組織推薦提取序列變化趨勢(shì)的有效工具，也是“應(yīng)用氣候?qū)W實(shí)習(xí)”課程中重要的教學(xué)內(nèi)容。本文以某氣象站1961—2015年年平均溫度為例，詳述了使用Matlab軟件計(jì)算MK突變檢驗(yàn)的UF和UB統(tǒng)計(jì)量，并基于Matlab/GUI實(shí)現(xiàn)了界面操作和自動(dòng)制圖。MK突變檢驗(yàn)計(jì)算的難點(diǎn)在于秩序列sk的計(jì)算，Matlab編程時(shí)可以通過(guò)雙重for循環(huán)來(lái)實(shí)現(xiàn)。通過(guò)計(jì)算程序能有效提高處理數(shù)據(jù)的速度，對(duì)學(xué)生綜合科學(xué)研究素質(zhì)和實(shí)際應(yīng)用水平的提高有較好的幫助。