经验正交函数分解matlab,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

数据读取1

2

3

4

5

6

7

8%% read data

file = &＃39;tauy_1958_2008_NWPacific.cdf&＃39;;

% ncdisp(file);

time = ncread(file,&＃39;time&＃39;);

lon = ncread(file,&＃39;lon&＃39;);

lat = ncread(file,&＃39;lat&＃39;);

tauy = ncread(file,&＃39;tauy&＃39;);

missing_value = ncreadatt(file,&＃39;tauy&＃39;,&＃39;missing_value&＃39;);

数据预处理

(1)显示数据的范围

使用m_map工具箱进行地图的显示：

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

(2)选择海区

1

2

3% set range

la = 46:59;

lo = 145:155;

(3)去趋势处理

1

2

3

4

5

6

7

8

9

10

11xx = 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计算1

2

3

4%% 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);

计算方差贡献率1V = [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%的前若干个模态的空间模态和时间系数

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49%% figure

figure;

% EOF2

plot(EOFs(:,2));

ylabel(&＃39;PC2&＃39;,&＃39;FontSize&＃39;,18);

xlabel(&＃39;Year&＃39;,&＃39;FontSize&＃39;,18)

axis([1 612 -0.2 0.2]);

set(gca,&＃39;xtick&＃39;,25:12*5:576,&＃39;xticklabel&＃39;,{1960:5:2005},&＃39;FontSize&＃39;,14);

set(gca,&＃39;tickdir&＃39;,&＃39;out&＃39; ,&＃39;LineWidth&＃39;,1.5);

set(gca,&＃39;Position&＃39;,[0.075 0.125 0.875 0.15]);

box off

% EOF1

subplot(323)

plot(EOFs(:,1));

ylabel(&＃39;PC1&＃39;,&＃39;FontSize&＃39;,18);

axis([1 612 -0.2 0.2]);

set(gca,&＃39;xtick&＃39;,25:12*5:576,&＃39;xticklabel&＃39;,{},&＃39;FontSize&＃39;,14);

set(gca,&＃39;tickdir&＃39;,&＃39;out&＃39; ,&＃39;LineWidth&＃39;,1.5);

set(gca,&＃39;Position&＃39;,[0.075 0.3 0.875 0.15]);

box off

% EC2

subplot(322)

m_proj(&＃39;Mercator&＃39;,&＃39;lat&＃39;,[45 60],&＃39;lon&＃39;,[130 160]);

m_pcolor(lon(xx),lat(yy),EC(:,:,2)&＃39;);

caxis([-1.0 2.5]) %设置颜色映射范围

shading interp;

m_coast(&＃39;patch&＃39;,[0.7 0.7 0.7]);

m_grid(&＃39;linestyle&＃39;,&＃39;none&＃39;,&＃39;box&＃39;,&＃39;fancy&＃39;,&＃39;tickdir&＃39;,&＃39;out&＃39;);

set(gca,&＃39;Position&＃39;,[0.5 0.5 0.45 0.4]);

title(&＃39;EOF2&＃39;,&＃39;FontSize&＃39;,18);

% EC1

subplot(321)

m_proj(&＃39;Mercator&＃39;,&＃39;lat&＃39;,[45 60],&＃39;lon&＃39;,[130 160]);

m_pcolor(lon(xx),lat(yy),EC(:,:,1)&＃39;);

caxis([-1.0 2.5])

shading interp;

m_coast(&＃39;patch&＃39;,[0.7 0.7 0.7]);

m_grid(&＃39;linestyle&＃39;,&＃39;none&＃39;,&＃39;box&＃39;,&＃39;fancy&＃39;,&＃39;tickdir&＃39;,&＃39;out&＃39;);

set(gca,&＃39;Position&＃39;,[0.05 0.5 0.45 0.4]);

title(&＃39;EOF1&＃39;,&＃39;FontSize&＃39;,18);

% colorbar

h = colorbar;

caxis([-1.0 2.5])

set(h,&＃39;position&＃39;,[0.9,0.5,0.02,0.4])

set(h,&＃39;ytick&＃39;,-1.0:0.7:2.5,&＃39;yticklabel&＃39;,{-1.0:0.7:2.5})

结果图