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

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

0 引 言

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

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

1 方法原理

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

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

(1)

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

在時間序列隨機獨立的假定下，定義統(tǒng)計量:

(2)

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

(3)

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

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

2 案例分析

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

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

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

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

打開Matlab軟件，在命令窗口中輸入代碼：

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

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

[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的計算，根據(jù)式(1)，秩序列sk是第i時刻數(shù)值大于j時刻數(shù)值個數(shù)的累計數(shù)，可以用雙重for循環(huán)計算。

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

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

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

UFk=zeros(size(y)); % 定義統(tǒng)計量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)); % 定義逆序統(tǒng)計量UBk，長度=y，初始值=0，UBk(1)=0

s=0;

y2=flipud(y); % 按時間序列逆轉(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 畫圖及保存圖片

Matlab軟件中plot() 函數(shù)能夠針對向量或矩陣的列來繪制曲線，是繪制二維圖形的最基本函數(shù)。本例中的M-K檢驗圖須調(diào)用函數(shù)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]) %設(shè)置圖形的位置及大小

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

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]);%設(shè)置X軸范圍和間距

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

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

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

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

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

ylabel(′統(tǒng)計量 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

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

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

2.4 基于Matlab-GUI的界面操作

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

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

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

(a)

(b)

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

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

(a)

(b)

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

3 結(jié) 語

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