北太平洋SST19801992的EOF分解及初步分析.docx
《北太平洋SST19801992的EOF分解及初步分析.docx》由会员分享，可在线阅读，更多相关《北太平洋SST19801992的EOF分解及初步分析.docx（13页珍藏版）》请在冰豆网上搜索。
北太平洋SST19801992的EOF分解及初步分析
北太平洋SST（1980-1992）的EOF分解及初步分析
数据来源：
本文利用来自英国气象局Hadley气候预测和研究中心的HadISST海温资料,对北太平洋（20N-80N130E-90W）1980至1992年156个月的SST进行研究,。
本文选取下载了1870年1月至2011年3月的全球SST数据，数据格式为*cdf，分辨率为1°。
以下给出了数据下载的途径和链接：
①点击页面中心处蓝色“maindatapage”：
图1
②选择SST数据下载:
图2
数据处理：
本文选取下载了1870年1月至2011年3月的全球SST数据，数据格式为*cdf,数据名为，截取北太平洋（20N-80N130E-90W）1980年1月至1992年12月156个月份的SST数据进行经验正交函数（EmpiricalOthorgnalFunction）分解,简记为EOF分解，得到该区域该时段的海温时空特征。
在编写Matlab程序过程中，应特别注意：
⑴剔除与其它站点相关系数小的站点的数据~简单的认为剔除陆地和冬季结冰点的数据；
⑵求距平值的协方差矩阵时，要进行逐月平均求距平，而不能是156个月的平均值，否则会导致第1模态的方差贡献率很大；
⑶当变量数m远大于观测样本数n时，导致协方差矩阵mRm=（nXm’）*（nXm）的阶数较大，可先求（nXm）*（nXm’）矩阵的特征值和特征向量，再求（nXm’）*（nXm）的特征值和特征向量，这叫做时空转换；
⑷M文件编写时要尽量减少循环量，提高运算速度；
⑸EOF分析能够有效地体现物理场主要信息，保留次要信息，并排除外来的随机干扰。
数据分析：
用主成分分析（PrincipalComponentAnalysis）的方法，即PCA对结果进行分析:
主成分分析是多元统计分析中一个非常重要的内容，它是一种从多个变量化为少数变量的统计方法。
由于多个变量之间是相互影响的，它们之间的关系是非常复杂的，为简化分析又不损失信息，并提取它们之间相互关系的主要特征，主成分分析利用多个变量之间的相互关系构造一些新变量，这些新变量不仅能综合反映原来多个变量的信息，而且彼此之间是相互独立的，同时是按方差贡献大小排列的。
方差贡献率小的变量通常规律性很差，其实际物理意义也不清晰，因此在实际分析过程中常常视为误差量或噪声而忽略，只取方差贡献率大的变量来研究，从而达到降维分析的目的。
通过对相应数据处理分析，前13个主成分的累积方差贡献率占总方差的，对前13个主模态的方差贡献率和累积方差贡献率列表格：
方差贡献率
累积方差贡献率
第1模态
第2模态
第3模态
第4模态
0.
第5模态
第6模态
0.
第7模态
第8模态
第9模态
第10模态
第11模态
第12模态
第13模态
图3
现仅列出北太平洋前5个主模态的空间分布填色图及时间序列，并对第1和第3模态进行分析：
北太平洋第1模态填色图及时间序列图4
北太平洋第2模态填色图及时间序列图5
北太平洋第3模态填色图及时间序列图6
北太平洋第4模态填色图及时间序列图7
北太平洋第5模态填色图及时间序列图8
对第1主模态进行分析:
北太平洋洋流图9
图10（a）
图10（b）
图10（a）是第1模态空间分布型，它解释海温场总方差的%，此型在北太平洋西、中部被一片强负值控制，负中心约在170°E，40°N和150°W，40°N附近，而北太平洋东部和北美沿岸为较弱的正值区，说明北太平洋西、中部海温与东部海温是反相关关系，负区与北太平洋西风漂流区（如图9）吻合。
由美国海洋学家斯蒂文•黑尔于1996年发现的太平洋年代际振荡（PDO）被科学研究的初步结果表明其与厄尔尼诺（ElNiño）和拉尼娜（LaNina）现象有着极其密切的关系。
该型可以反映和PDO有关的大尺度分布特征，因此这种分布型是全球海洋与大气相互作用的一个重要组成部分，它是北太平洋海温非季节变化的最重要的型式。
MonthlyvaluesforthePDOindex:
1900—January2008图11
图10（c）
观察发现图11（MonthlyvaluesforthePDOindex:
1900—January2008）1980年至1992年时间段的指数和第1模态的时间序列图10（c）有很好的对应关系，可以验证北太平洋海表面温度第1模态空间分布型确实与PDO有很强的相关性。
资料显示，近100多年来，PDO已出现了两个完整的周期：
第一周期的“冷位相”发生于1890年至1924年，而1925年至1946年为“暖位相”；第二周期的“冷位相”出现于1947年至1976年，1977年至90年代后期为“暖位相”。
当PDO现象以“暖位相”形式出现时，北美大陆附近海面的水温就会异常升高，而北太平洋洋面温度却异常下降。
并且，在20-30年的冷、暖位相中，会存在短期的反向指数。
由时间序列可知：
1980年至1988年底，时间序列指数基本为正值，说明图10（a）中蓝色区域海表面温度低于红色区域海表面温度，即北太平洋西、中部海温低于东部海温；1989年初至1992年初，时间序列指数为负值，说明10（a）中蓝色区域海表面温度高于红色区域海表面温度，即北太平洋西、中部海温高于东部海温。
此分析和历史资料相吻合。
对第3主模态进行分析:
图12（a）
图11（a）是第3模态空间分布型，它解释海温场总方差的%，此型在北太平洋阿留申群岛南部被一片强正值控制，正值中心约在150°W，40°N附近，而其西南部和日本海海域为一片较强的负区与之相互补偿，其东南部北美洲沿岸为较弱的负值区，说明北太平洋中部海表面温度与东、西两侧海表面温度呈反相关。
正、负中心基本上与北太平洋大气活动中心（阿留申低压与西太平洋高压）对应。
表明正、负区域是海气相互作用最活跃的区域。
查阅资料可知，这种模态的变化，对我国的天气与气候有明显的影响。
图12（b）
结合第3模态的时间序列可知:
1980年至1983年6月、1987年6月至1988年底及1990年，时间序列指数基本为负值，说明图10（a）中蓝色区域海表面温度高于红色区域海表面温度，即北太平洋西、中部海温高于东部海温；1983年7月至1987年初、1989年及1991年初至1992年底，时间序列指数为正值，说明10（a）中蓝色区域海表面温度低于红色区域海表面温度，即北太平洋西、中部海温高于东部海温。
总结：
本文通过对北太平洋1980年至1992年SST数据处理、EOF分解和初步分析，掌握了主成分分析、时空转换的原理和方法，提高了Matlab的编程和绘图技巧并对主模态和时间序列的分析进行了联系。
在此过程中克服了很多困难，受益匪浅。
但和很多同学相比，仍然有很大差距，将会更加虚心请教，刻苦钻研，以取得不断进步。
参考文献：
[1]左军成.海洋水文环境要素的分析方法和预报
[2]胡基福.气象统计原理与方法
[3]黄嘉佑.气象统计分析与预报方法
[4]杜凌.海洋要素计算（2011）PPT
[5]姜霞.气象统计原理与方法（2011）PPT
M文件：
clear;clc;closeall
address='E:
\oceanelement\';
fid=（address,'NC_NOWRITE'）;
sstid=（fid,'sst'）;
sst=（fid,sstid）;%读取nc格式数据
%**************************************************************************
sst1=sst（1:
90,11:
70,1320:
1475）;%选取所需要区域的数据
sst2=sst（311:
360,11:
70,1320:
1475）;
sst3=zeros（140,60,156）;
sst3（90:
-1:
1,1:
60,1:
156）=sst1;
sst3（140:
-1:
91,1:
60,1:
156）=sst2;
sst=sst3;
%**************************************************************************
sst_area1=zeros（156,8400）;%zeros全零数组
fori=1:
156;
squ=squeeze（sst（:
:
i））;%执行该指令后sst数据转换为二维数组
sst_area1（i,:
）=reshape（squ,1,8400）;%将数据转变为二维
end
%**************************************************************************
%剔除与其它站点相关系数小的站点的数据~简单的认为剔除陆地和冬季结冰点的数据
sst_area1（sst_area1<-10000）=NaN;%陆地和冰点的填充值为+30~将此值定义为NaN
%i=1;
%forj=1:
8400
%ifsst_area1（i,j）==
%sst_area1（i,j）=NaN;%冰点的填充值为
%i=i+1;
%end
%end
sst_nan=isnan（sst_area1）;
i=0;
forj=1:
8400
ifsum（sst_nan（:
j））==0;
i=i+1;
sst_region（:
i）=sst_area1（:
j）;
end
end
%**************************************************************************
%求距平~注意季节的变换
X=zeros（size（sst_region））;%学者给的程序
fori=1:
12
X（i:
12:
144+i,:
）=sst_region（i:
12:
end,:
）-repmat（mean（sst_region（i:
12:
end,:
）,1）,size（sst_region（i:
12:
end,:
）,1）,1）;
end
%s=0;%自己的~很多循环计算慢
%fork=1:
12;
%forj=1:
4344;
%fori=k:
12:
144+k
%s=s+sst_region（i,j）;
%end
%s_ave=s/13;
%fori=k:
12:
144+k
%ave_region（i,j）=s_ave;
%end
%s=0;
%end
%end
%X=sst_region-ave_region;%X1是156*8400的矩阵
%**************************************************************************
%由于变量数（8400）>观测样本数~协方差矩阵的阶数较大~因此可用时空转换的方法提高计算速度
R=X*X';%协方差矩阵R=X*X'是8400*8400的方阵~现定义矩阵R=X'*X是156*156的矩阵
[v,d]=eig（R）;%进行EOF分解~因为X'*X与X*X'的秩相同所以特征值相同~d为x的特征值组成的对角阵~v为X*X'的特征向量~
v=fliplr（v）;%矩阵作左右翻转
d=rot90（d,2）;%矩阵上下翻转后再左右翻转（查看生成的对角阵是由小到大排列的~此指令可使其由大到小排列~fliplr（flipud（d））=rot90（d,2））
diagonal=diag（d）;
spacef=X'*v;
fori=1:
156;
spacef（:
i）=spacef（:
i）/sqrt（diagonal（i））;%空间本征函数
end
timef=X*spacef;%时间本征函数
sum_d=sum（diagonal）;
count=0;
fori=1:
156;
count=count+diagonal（i）;
G1（i）=count/sum_d;%G1（i）是累积方差贡献率
end
fori=1:
156;
G2（i）=diagonal（i）/sum_d;%G2（i）是方差贡献率
end
%**************************************************************************
%将删去的陆地与冰点的填充值补回
sst_area2=zeros（156,8400）;
sst_area2（:
:
）=NaN;
i=1;
spacef2=spacef';
forj=1:
8400
ifsum（sst_nan（:
j））==0;
sst_area2（:
j）=spacef2（:
i）;
i=i+1;
end
end
sst_area3=sst_area2';
%**************************************************************************
%画图
%subplot（2,1,2）%将绘图窗口划分为2*1个子窗口，并在第2个子窗口中绘图
x=1:
156;
plot（x,timef（:
6）,'g'）;
ylim（[-8080]）;
%xlabel（'1980-1992年156个月','fontsize',15,'fontname','隶书'）
ylabel（'INDEX','fontsize',12,'fontname','黑体'）
set（gca,'xtick',（1:
6:
162））
set（gca,'xticklabel',{'1980','','1981','','1982','','1983','','1984','','1985','','1986','','1987','','1988','','1989','','1990','','1991','','1992','','1993'}）
title（'北太平洋第6模态1980至1992年SST时间序列','color','k','fontsize',15,'fontname','幼圆'）
gridon
holdoff
%%subplot（2,1,1）
%lon=[:
];
%lat=[:
];
%m_proj（'EquidistantCylindrical','lat',[],'lon',[]）;
%m_contourf（lon,lat,rot90（reshape（sst_area3（:
6）,140,60）',2）,30,'linestyle','none'）;
%colorbar
%m_coast（'patch',[.95.95.95]）;
%m_coast（'color','k'）;
%m_grid（'linestyle','none','tickdir','out','linewidth',;
%%xlabel（'longitude','fontsize',15,'fontname','comicsansms'）;
%%ylabel（'latitude','fontsize',15,'fontname','comicsansms'）;
%title（['北太平洋第6模态填色图'],'fontsize',15,'fontname','幼圆'）;
lon=[:
];
lat=[:
];
m_proj（'lambert','long',[],'lat',[]）;
m_contourf（lon,lat,rot90（reshape（sst_area3（:
1）,140,60）',2））;
%colorbar;
m_coast（'patch',[1.85.7]）;
m_elev（'contourf',[500:
500:
6000]）;
m_grid（'box','fancy','tickdir','in'）;
colormap（flipud（copper））;
%xlabel（'longitude','fontsize',15,'fontname','幼圆'）;
%ylabel（'latitude','fontsize',15,'fontname','幼圆'）;
title（['北太平洋第1模态填色图'],'fontsize',15,'fontname','幼圆'）;
升级会员 