matlab经验正交函数EOF(转载)

A.7EOF分析

经验正交函数分析方法(empiricalorthogonalfunction,缩写为EOF)，也称特征向量分析(eigenvectoranalysis)，或者主成分分析(principalcomponentanalysis,缩写PCA)，是一种分析矩阵数据中的结构特征，提取主要数据特征量的一种方法。Lorenz在1950年代首次将其引入气象和气候研究，现在在地学及其他学科中得到了非常广泛的应用。地学数据分析中通常特征向量对应的是空间样本，所以也称空间特征向量或者空间模态；主成分对应的是时间变化，也称时间系数。因此地学中也将EOF分析称为时空分解。

原理与算法

选定要分析的数据，进行数据预处理，通常处理成距平的形式。得到一个数据矩阵Xm×n

计算X与其转置矩阵XT的交叉积，得到方阵

Cm×m=

1

X×XTn

如果X是已经处理成了距平的话，则C称为协方差阵；如果X已经标准化(即C中每行数据的平均值为0，标准差为1)，则C称为相关系数阵 计算方阵C的特征根(λ1,...,m)和特征向量Vm×m，二者满足

Cm×m×Vm×m=Vm×m×∧m×m

其中∧是m×m维对角阵，即

λ1

...

0

0λ2...0

∧= ............

00...λm

一般将特征根λ按从大到小顺序排列，即λ1>λ2>...>λm。因为数据X是真实的观测值，所以λ应该大于或者等于0。每个非0的特征根对应

42

一列特征向量值，也称EOF。如λ1对应的特征向量值称第一个EOF模态，也就是V的第一列即EOF1=V(:,1)；第λk对应的特征向量是V的第k列，即EOFk=V(:,k)。

计算主成分。将EOF投影到原始资料矩阵X上，就得到所有空间特征向量对应的时间系数(即主成分)，即

T

PCm×n=Vm×m×Xm×n

其中PC中每行数据就是对应每个特征向量的时间系数。第一行PC(1,:)就是第一个EOF的时间系数，其他类推。

上面是对数据矩阵X进行计算得到的EOF和主成分(PC)，因此利用EOF和PC也可以完全恢复原来的数据矩阵X，即

X=EOF×PC

有时可以用前面最突出的几个EOF模态就可以拟合出矩阵X的主要特征。此

1

外，EOF和PC都具有正交性的特点，可以证明nPC×PCT=∧；即不同的PC之

间相关为0。E×ET=I。I为对角单位矩阵，即对角线上值为1，其他元素都为0。这表明各个模态之间相关为0，是独立的。

由上面的计算过程可以看出，EOF分析的核心是计算矩阵C的特征根和特征向量。计算矩阵特征根和特征向量的方法很多，下面具体给出Matlab中进行EOF分析的两种不同的方法。具体步骤可参考下面两个框图中的实例。

方法1：调用[EOF,E]=eig(C)，其中EOF为计算得到的空间特征向量，E为特征根。然后计算主成分PC=EOFT×X。需要指出的时，当数据量很大时，例如分析高分辨率的资料(如1km分辨率的NDVI资料)，空间范围很大维数m很容易超过数万个点，则矩阵C的维数是个巨大量，需要占用大量内存，也会导致计算速度异常缓慢。而且很可能超出计算机的计算极限而死机。

方法2：直接对矩阵X进行奇异值分解

X=UVT

其中为奇异值对交阵(对角线上的元素为奇异值)，奇异值与特征根成倍数关系。

43

如果矩阵C=

=； √T

如果矩阵C=XX，C的特征根为λ，则有=由于该方法是直接对矩阵X进行分解，所以对内存的要求远小于方法1。计算速度很快。

两种方法对比练习。

1

XXT，C的特征根为λ，则有

√

显著性检验

可以证明

m i=1

m k=1

m k=1

Xi=

λk=

PCk

这说明矩阵X的方差大小可以简单的用特征根的大小来表示。λ越高说明其对应的模态越重要，对总方差的贡献越大。第k个模态对总的方差解释率为

λk

i=1λi

×100%

即使是随机数或者虚假数据，放在一起进行EOF分析，也可以将其分解成一系列的空间特征向量和主成分。因此，实际资料分析中得到的空间模态是否是随机的，需要进行统计检验。North等(1982)的研究指出，在95%置信度水平下的特征根的误差

λ=λ

2N

λ是特征根，N 是数据的有效自由度，这在前面相关系数分析中已经有介绍(见4页相关内容)。将λ按顺序依次检查，标上误差范围。如果前后两个λ之间误差范围有重叠，那么他们之间没有显著差别。

图A.16是对1949 2002年北半球1月平均海平面气压，做距平处理处理及面积EOF分析的结果。从特征根误差范围看，第一和第二模态存在显著差别，第二和第三模态之间也存在显著差别。但是第三特征根和第四及以后的特征根之间没有显著的差别。如果要分析主要的模态的话，最好只选择前三个进行分析。

44

练习：利用[E,V]=eig(C)计算矩阵X的特征向量和主成分%

X=[26152;

94054];

X(1,:)=X(1,:)-mean(X(1,:));X(2,:)=X(2,:)-mean(X(2,:));得到X的距平值：X=

-1.202.80-2.201.80-1.204.60-0.40-4.400.60-0.40%%%co-variancematrixC=X*X’/5;协方差阵C=

3.760.920.928.24[EOF,E]=eig(C);%V:eigenvectors;E:eigenvaluesPC=EOF’*X;

%%reversetheorderE=fliplr(flipud(E))

lambda=diag(E);%retaineigenvaluesonlyEOF=fliplr(EOF)PC=flipud(PC)得到EOF=

0.190.98

-0.980.19

得到特征根E=

8.420

03.58

得到主成分PC=

4.280.152.07-2.82

-4.74

1.310.94-1.65-0.621.10

%%check

EOF*EOF’%=I

检查EOF的正交性得到：

1.000

01.00

PC*PC’/5%=lambda检查PC的正交性得到：

8.420.000.003.58EOF*PC%=X

可以完全恢复X的距平值：

-1.202.80-2.204.60-0.40-4.40

1.80

0.60-1.20-0.40

45

练习：利用[U,S,V]=svd(X)计算矩阵X的特征向量和主成分

X=[26152;

94054];

X(1,:)=X(1,:)-mean(X(1,:));X(2,:)=X(2,:)-mean(X(2,:));X的距平是：

-1.202.804.60-0.40

-2.20

-4.401.800.60-1.20-0.40

[U,S,V]=svd(X);得到U=

0.190.980.98-0.19S=

6.49000004.23000V=

0.66-0.490.560.020.670.63-0.73-0.310.530.140.390.03-0.10-0.26-0.02EOF=U;PC=S*V’;得到PC=

4.28-2.07

0.09-0.320.250.910.06-0.060.22-0.160.060.96

0.152.82-4.74-1.310.941.65-0.62-1.10

E=S.^2/5;%=lambda

E的数值与上面得到的特征根完全一样即E=：

8.420000

03.58000EOF*PC%=X

可以完全恢复X的距平值：-1.202.80-2.204.60-0.40-4.40

1.80

0.60-1.20-0.40

46

图A.16:北半球1月海平面气压EOF分析的第一特征向量.(a)为特征根及95%信度误差，(b)第一特征向量，(c)第一主成分，(d)第一主成分偏强+σ时海平面气压的变化量(hPa).1949 2002，NCEP/NCAR再分析资料

47

结果展示

通常情况下，主成分是有单位的，即反映的是矩阵X的单位，而空间特征向量是无量纲的。不过实际应用中常常对EOF分析得到的主成分和特征向量进行标准化处理得到新的PC 和EOF

PC(k)

PC (k)=k

EOF(k)=EOF(k)λk

或者是简单地将PC标准化，使得其平均值为0，标 …… 此处隐藏：5329字，全部文档内容请下载后查看。喜欢就下载吧 ……