自然正交函数分析(EOF)程序

５．３自然正交函数分析（ＥＯＦ）程序

近年来，自然正交函数（又称经验正交函数）展开在气象上应用比较广泛。这种正交函数展开不象三角函数展开、球函数展开那样有固定的展开形式。它无固定的函数形式，不是事先人为地给定典型场函数，图形是由场本身来决定的，它具有收敛快又能更好地反映出场的基本结构的特征。它可以在有限的区域中进行，既可以取空间不同站点进行分解，也可以对同一站点的不同时间、不同高度的多种要素进行综和分析。因此它在气象中具有广泛的应用，可用于气象要素场分析、大气垂直结构分析、动力模型垂直分层等。

5.3.1功能

计算要素场的自然正交函数分解。

5.3.2方法说明

自然正交函数分解是针对气象要素场进行的，它的基本思想是把包含ｐ个空间点（或ｐ个变量）的n 个时次的观测场随时间进行分解，即将某一区域的气象要素场序列ij F （i=1, 2, …,p ；j=1,2,…,n ，即p 个空间点的n 个时次的观测资料）分解成相互正交的时间函数与相互正交的空间函数的乘积之和，常把空间函数ik v 看作典型场，时间函数kj t 看作典型场的权重系数，则不同时间的要素场是若干个典型场按不同权重线性叠加的结果，各个场之间的差别就在于各典型场的系数不同。则气象要素场可以表示为

∑=+++==p 1k pj ip j 22i j 11i k j ik ij t v t v t v t v

F （５．３．１）

其中F ij 表示第i 个场中的第j 个测点的观测值。

可将(5.3.1)是写为矩阵的形式

VT F = （５．３．２）

式中F 为n p ⨯阶的均值为0的资料阵，V 为p p ⨯阶的空间函数阵，T 为n p ⨯阶的时间函数阵。由于V 和T 是根据场的资料阵F 进行分解而得到的，分解的函数没有固定的函数形式，因而称为“经验”的，另外，我们还要求这种分解具有“正交”性，即要求满足下式

⎪⎪⎭

⎪⎪⎬⎫≠=='≠=='∑∑==)l k (0t t t t )l k (0v v v v n 1j lj k j l k p 1i il ik l k （5．3．３）

事实上，我们对（５．３．２）式右乘T '可得

V T VT F F ''=' （５．３．４）

因F F '是p p ⨯阶对称阵，其元素为距平变量的交叉积。根据实对称矩阵的分解定理有 V V ΛF F '=' （５．３．５）

其中Λ是F F '矩阵的特征值组成的对角阵，V 是对应的特征向量为列向量组成的矩阵。比较（５．３．４）和（５．３．５）式可知

ΛT T =' （５．３．６）

又根据特征向量的性质有

I V V V V ='=' （５．３．７）

式中为I单位矩阵。显然（５．３．６）和（５．３．７）式满足（５．３．３）式的要求。由

F'矩阵的特征向量求得，而时间函数则可利用（５．３．２）式左此可知空间函数矩阵可从F

乘V'得到，即

T'

=（５．３．８）

V

F

5.3.3子程序语句

CALL EOF(X,P,N,XF)

5.3.4哑元说明

X——输入变量，二维实型数组，大小为P⨯N，存放原始观测值。

P——输入整型变量，空间格点数。

N——输入整型变量，序列的时间长度。

XF——输出变量，二维实型数组，大小为P⨯N，存放恢复值。

5.3.5子程序

SUBROUTINE EOF(X,P,N,LW,XF)

INTEGER::P

INTEGER::N

INTEGER::LW

REAL(8),DIMENSION(P,N)::X,XF

REAL(8),DIMENSION(P,P)::A,V,V1

REAL(8),DIMENSION(P,N)::T

REAL(8),DIMENSION(P)::B,GM,GA

REAL(8),DIMENSION(P,LW)::VF

REAL(8),DIMENSION(LW,N)::TF

! 求X乘以X的转置,即A=XXˊ

DO I=1,P

DO J=1,P

A(I,J)=0

DO K=1,N

A(I,J)=A(I,J)+X(I,K)*X(J,K)

END DO

END DO

END DO

! 用Jacobi法求A的特征值和特征向量

! 返回时B存放矩阵的全部特征值，V存放特征向量为列组成的矩阵

CALL JCB(A,P,1.0E-6,V,B,L)

DO I=1,P

GA(I)=0

DO J=1,I

GA(I)=GA(I)+B(J)

END DO

END DO

DO I=1,P

GM(I)=GA(I)/GA(P)

END DO

DO I=1,P

V1(I,J)=V(J,I)

END DO

END DO

T=MATMUL(V1,X)

WRITE(12,'(" 特征值")')

WRITE(12,'(

I10)')(I,I=1,P)

WRITE(12,'(3X,

D10.4)')B

WRITE(12,'(" 解释的方差(%)")')

WRITE(12,'(

I7)')(I,I=1,P)

WRITE(12,'(3X,

F7.2)')GM*100

WRITE(12,'(" 特征向量为列组成的矩阵，即空间函数V")')

WRITE(12,'(

F7.4)')((V(I,J),J=1,P),I=1,P)

WRITE(12,'(" 时间函数T")')

WRITE(12,'(

F10.4)')((T(I,J),J=1,P),I=1,N)

DO I=1,P

DO J=1,LW

VF(I,J)=V(I,J)

END DO

END DO

DO I=1,LW

DO J=1,N

TF(I,J)=T(I,J)

END DO

END DO

XF=MATMUL(VF,TF)

END

SUBROUTINE JCB(A,N,EPS,V,B,L)

! A:调用时存放实对称矩阵

! B:返回时存放矩阵的全部特征值

! V:存放特征向量,其中第i列为与第i个特征值相对应的特征向量! EPS:存放精度要求

REAL(8),DIMENSION(N,N)::A,V

REAL(8),DIMENSION(N)::B

REAL(8)::FM,CN,SN,OMEGA,X,Y

INTEGER::P,Q

L=1

V=0

DO I=1,N

V(I,I)=1

END DO

10 FM=0

DO I=1,N