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循環來實現。通過計算程序能有效提高處理數據的速度，對學生綜合科學研究素質和實際應用水平的提高有較好的幫助。