arrioH"air«rainrtructir«mhara«rrrn«d«mcrahampgHnMDui^bree«—fiHna-ili-qi^iikLarmbfw-bii・・dk«sancmLtnAora-fcrihiirr^z«rLtri■urfidrrwhrtVic+hmintrnvidiin4ihArrljeandl*=4-rilgcrfthni
defitiwiciM-iniheMM伸占Mbjrwk«ng如h■领r>«almstouwnceab1部口心草质5M«Ri«虫u.SSTs供占rot占此mtinu姑日wig泌iiigH网用尽时小申十SST占rd皿占由的fKMirMBK
rar■ditiiitaddaia-pb=4iodth<日做■■■〔:
.ndillpi-ndidiEnpraraii,■»■*tt>iEndp^arCibrifirinmladjc-n*!
.WbracarnrrMndadjiaura«dthiiNliiniuingih«drta.
n.Mir・■rH^Jahl'afrom4i・dcTnlnMmwMlh-zuldiM-gifar1h«hif^4£«£tf^rn*M£-■»』!
■Hidirkc4iXftrbdp4ti^«ir4dffrt丫
2)http:
//www.metoffice.gov.uk/hadobs/hadisst/data/download.html
选择SST数据下载:
图2
数据处理:
本文选取下载了1870年1月至2011年3月的全球SST数据,数据格式为*.netcdf,数据名为
HadISST_sst.nc,截取北太平洋(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个主成分的累积方差贡献率占总方差的
0.810567965759345,对前13个主模态的方差贡献率和累积方差贡献率列表格:
方差贝献率
累积方差贝献率
第1模态
0.2190372311380337
0.2190372311380337
第2模态
0.12913861798280676
0.34817584912084043
第3模态
0.1000386758391279
0.44821452495996833
第4模态
0.08507766885727319
0.5332921938172415
第5模态
.0592********
0.5925626300087424
第6模态
.0537********
0.6463602159422434
第7模态
0.0424686060052352
0.688828821947479
第8模态
.028*********
0.717529769716775
第9模态
.022*********
0.740483642439213
第10模态
.021*********
0.761598282240791
第11模态
0.0192806477668055
0.780878930007597
第12模态
0.0163162166248416
0.797195146632438
第13模态
0.0133728191269071
0.810567965759345
图3
现仅列出北太平洋前5个主模态的空间分布填色图及时间序列,并对第1和第3
模态进行分析:
北太平洋第1模态填色图及时间序列图4
北太平洋第2模态填色图及时间序列图5
北太平居尊〕痕蒲菱②倒
北太平洋第3模态填色图及时间序列图6
北太平洋第4模态填色图及时间序列图7
北太平洋第5模态填色图及时间序列图8
对第1主模态进行分析:
北太平洋洋流图9
北太平洋第1模态填色图
60PNH
-aa?
003
•CI.Q4
150PE
lECPW
150如V
1泌轩
图10(a)
北太平洋第1模态填曳图
图10(b)
图10(a)是第1模态空间分布型,它解释海温场总方差的22.9%,此型在北太平洋西、中部
被一片强负值控制,负中心约在170°E,40°N和150°W,40°N附近,而北太平洋东部和北美沿岸为较弱的正值区,说明北太平洋西、中部海温与东部海温是反相关关系,负区与
北太平洋西风漂流区(如图9)吻合。
由美国海洋学家斯蒂文?
黑尔于1996年发现的太平洋年代际振荡(PDO)被科学研究的初步
结果表明其与厄尔尼诺(ElNi?
o)和拉尼娜(LaNina)现象有着极其密切的关系。
该型可以反映和PDO有关的大尺度分布特征,因此这种分布型是全球海洋与大气相互作用的一个重要组成部分,它是北太平洋海温非季节变化的最重要的型式。
monthlyvaluesforlhePDOindex:
1900-January2008
itoo19201940I96019802000
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)
10.0%,此型在北太平洋阿留申群
N附近,而其西南部和日本海海域
图11(a)是第3模态空间分布型,它解释海温场总方差的
岛南部被一片强正值控制,正值中心约在150°W,40
为一片较强的负区与之相互补偿,其东南部北美洲沿岸为较弱的负值区,说明北太平洋中部海表面温度与东、西两侧海表面温度呈反相关。
正、负中心基本上与北太平洋大气活动中心
(阿留申低压与西太平洋高压)对应。
表明正、负区域是海气相互作用最活跃的区域。
查阅资料可知,这种模态的变化,对我国的天气与气候有明显的影响。
图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\HadISST_sst.nc';
fid=netcdf.open(address,'NC_NOWRITE');
sstid=netcdf.inqVarid(fid,'sst');
sst=netcdf.getVar(fid,sstid);%
%**************************************************************************
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_area1(i,:
)=reshape(squ,1,8400);%
end
%**************************************************************************
读取nc格式数据
选取所需要区域的数据
全零数组
执行该指令后sst数据转
将数据转变为二维
%剔除与其它站点相关系数小的站点的数据~简单的认为剔除陆地和冬季结冰点的数据
sst_area1(sst_area1<-10000)=NaN;%陆地和冰点的填充值为
-1.00000001504747e+30~将此值定义为NaN
%i=1;
%forj=1:
8400
%ifsst_area1(i,j)==-1.79999995231628
%sst_area1(i,j)=NaN;%
%i=i+1;
%end
%end
sst_nan=isnan(sst_area1);
冰点的填充值为-1.79999995231628
i=0;
forj=1:
8400
ifsum(sst_nan(:
j))==0;
i=i+1;
sst_region(:
i)=sst_area1(:
j);
end
end
%**************************************************************************
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=[130.5:
269.5];
%lat=[