matlab物理场的正交分解

经验正交分解的原理

简介

经验正交函数分析方法(empirical orthogonal function,缩写为EOF)，也称特征向量分析(eigen vector analysis)，或者主成分分析(principal component analysis,缩写PCA)，是一种分析矩阵数据中的结构特征，提取主要数据特征量的一种方法。Lorenz在1950年代首次将其引入气象和气候研究，现在在地学及其他学科中得到了非常广泛的应用。

分析中通常特征向量对应的是空间样本，称空间特征向量或者空间模态；对应的是时间变化称为时间系数。

特点

（1）主要特色：能够有效地体现物理场的主要信息，保留次要信息并排除外来的随机干扰。
（2）无固定的函数形式：这种正交函数展开不像三角函数展开、球函数展开那样有固定的展开形式。
（3）图形是由场本身决定的，不是事先人为地给定典型场函数的。
（4）具有收敛快、能更好地反应出场的基本结构的特征。

主要应用

（1）可在有限的区域中不规则分布的站点进行：空间不同站点/同一点的不同时间、不同高度的多种要素进行综合分析。
（2）可用于要素场分析、垂直结构分析、动力模型垂直分层等。

基本思想

经验正交函数分解的基本思想是把包含m个空间点（或p个变量）的n个时次的观测场的序列（i = 1,2,…,m; j = 1,2,…,n）分解成相互正交的时间函数与相互正交的空间函数的乘积之和，常把空间函数看做典型场，时间函数看做典型场的权重系数，则不同时间的要素场是若干个典型场按不同的权重线性叠加的结果，各个场之间的差别就在于各典型场的系数不同。

$$ nXm = nEm \times m\Phi_{m}^{T} $$

EOF函数简介

EOF函数的原型为：

$$ [V,EOFs,EC,error]=EOF(D,p) $$

本站提供的下载:EOF.m

输入参数

D：输入的物理场数据，为一个二维矩阵，每一行假定为一个样本，每一列为一个变量。因此D的一列表示一个时间序列的变量。
P：可选项，指EOF分析结果的模态数，其取值应小于时间间隔和样本数

输出参数

V：与EOFs有关的特征向量；
EOFs：时间系数；
EC：空间模态，由p个模态的列向量组成的矩阵；
error：重建误差；

经验正交分解实例

数据下载

前往SODA: Simple Ocean Data Assimilation下载数据
选择相关变量，以“TAUY”为例。
合理选择研究区域，时间范围等

选择合适数据下载后，认真核对数据基本信息

在数据下载页选择合适的数据格式，我们选择NetCDF

现在开始提取数据，提取的数据名为data.cdf
建议根据下载的数据变量名、范围、时间等命名，便于查找。如tauy_1958_2008_NWPacific.cdf

数据读取

%% read data
file = 'tauy_1958_2008_NWPacific.cdf';
% ncdisp(file);
time = ncread(file,'time');
lon = ncread(file,'lon');
lat = ncread(file,'lat');
tauy = ncread(file,'tauy');
missing_value = ncreadatt(file,'tauy','missing_value');

数据预处理

（1）显示数据的范围
使用m_map工具箱进行地图的显示：

1	>> figure;m_proj('mercator','lat',[20 60],'lon',[100 160]);m_coast('patch',[0.7 0.7 0.7]);m_grid;

（2）选择海区

1
2
3

% set range
la = 46:59; 
lo = 145:155;

（3）去趋势处理

xx = find(lon>=min(lo) & lon<=max(lo) );  
yy = find(lat>=min(la) & lat<=max(la));
tauy_sub = zeros(length(xx),length(yy),612);
tauy_sub = tauy(xx,yy,:)
%% detrend
tauy_de = zeros(size(tauy_sub));
for ii = 1:length(xx)
    for jj = 1:length(yy)
        tauy_de(ii,jj,:) =  detrend(squeeze(tauy_sub(ii,jj,:)));
    end
end

EOF计算

%% EOF
tauy_de = reshape(tauy_de,length(xx)*length(yy),612);
[V,EOFs,EC,error] = EOF(tauy_de,10);
EC = reshape(EC,length(xx),length(yy),10);

计算方差贡献率

1	V = [2.3544 0.1429 0.0625 0.0410 0.0305 0.0085 0.0077 0.0073 0.0024 0.0021]

计算各模态的方差贡献率：

1 2	>> V/sum(V) ans =[0.8853 0.0537 0.0235 0.0154 0.0115 0.0032 0.0029 0.0028 0.0009 0.0008]

绘图

根据方差贡献率，画出方差贡献率大于90%的前若干个模态的空间模态和时间系数

%% figure
figure;
% EOF2
plot(EOFs(:,2));
ylabel('PC2','FontSize',18);
xlabel('Year','FontSize',18)
axis([1 612 -0.2 0.2]);
set(gca,'xtick',25:12*5:576,'xticklabel',{1960:5:2005},'FontSize',14);
set(gca,'tickdir','out' ,'LineWidth',1.5);
set(gca,'Position',[0.075 0.125 0.875 0.15]);
box off
% EOF1
subplot(323)
plot(EOFs(:,1));
ylabel('PC1','FontSize',18);
axis([1 612 -0.2 0.2]);
set(gca,'xtick',25:12*5:576,'xticklabel',{},'FontSize',14);
set(gca,'tickdir','out' ,'LineWidth',1.5);
set(gca,'Position',[0.075 0.3 0.875 0.15]);
box off
% EC2
subplot(322)
m_proj('Mercator','lat',[45 60],'lon',[130 160]);
m_pcolor(lon(xx),lat(yy),EC(:,:,2)');
caxis([-1.0  2.5])  %设置颜色映射范围
shading interp;
m_coast('patch',[0.7 0.7 0.7]);
m_grid('linestyle','none','box','fancy','tickdir','out');
set(gca,'Position',[0.5 0.5 0.45 0.4]);
title('EOF2','FontSize',18);
 
% EC1
subplot(321)
m_proj('Mercator','lat',[45 60],'lon',[130 160]);
m_pcolor(lon(xx),lat(yy),EC(:,:,1)');
caxis([-1.0  2.5])
shading interp;
m_coast('patch',[0.7 0.7 0.7]);
m_grid('linestyle','none','box','fancy','tickdir','out');
set(gca,'Position',[0.05 0.5 0.45 0.4]);
title('EOF1','FontSize',18);
 
% colorbar
h = colorbar;
caxis([-1.0  2.5])
set(h,'position',[0.9,0.5,0.02,0.4])
set(h,'ytick',-1.0:0.7:2.5,'yticklabel',{-1.0:0.7:2.5})