５．３自然正交函数分析（ＥＯＦ）程序近年来，自然正交函数（又称经验正交函数）展开在气象上应用比较广泛。
这种正交函数展开不象三角函数展开、球函数展开那样有固定的展开形式。
它无固定的函数形式，不是事先人为地给定典型场函数，图形是由场本身来决定的，它具有收敛快又能更好地反映出场的基本结构的特征。
它可以在有限的区域中进行，既可以取空间不同站点进行分解，也可以对同一站点的不同时间、不同高度的多种要素进行综和分析。
因此它在气象中具有广泛的应用，可用于气象要素场分析、大气垂直结构分析、动力模型垂直分层等。
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 vF （５．３．１）其中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'=（５．３．８）VF5.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::PINTEGER::NINTEGER::LWREAL(8),DIMENSION(P,N)::X,XFREAL(8),DIMENSION(P,P)::A,V,V1REAL(8),DIMENSION(P,N)::TREAL(8),DIMENSION(P)::B,GM,GAREAL(8),DIMENSION(P,LW)::VFREAL(8),DIMENSION(LW,N)::TF! 求X乘以X的转置,即A=XXˊDO I=1,PDO J=1,PA(I,J)=0DO K=1,NA(I,J)=A(I,J)+X(I,K)*X(J,K)END DOEND DOEND DO! 用Jacobi法求A的特征值和特征向量! 返回时B存放矩阵的全部特征值，V存放特征向量为列组成的矩阵CALL JCB(A,P,1.0E-6,V,B,L)DO I=1,PGA(I)=0DO J=1,IGA(I)=GA(I)+B(J)END DOEND DODO I=1,PGM(I)=GA(I)/GA(P)END DODO I=1,PV1(I,J)=V(J,I)END DOEND DOT=MATMUL(V1,X)WRITE(12,'(" 特征值")')WRITE(12,'(<P>I10)')(I,I=1,P)WRITE(12,'(3X,<P>D10.4)')BWRITE(12,'(" 解释的方差(%)")')WRITE(12,'(<P>I7)')(I,I=1,P)WRITE(12,'(3X,<P>F7.2)')GM*100WRITE(12,'(" 特征向量为列组成的矩阵，即空间函数V")')WRITE(12,'(<P>F7.4)')((V(I,J),J=1,P),I=1,P)WRITE(12,'(" 时间函数T")')WRITE(12,'(<P>F10.4)')((T(I,J),J=1,P),I=1,N)DO I=1,PDO J=1,LWVF(I,J)=V(I,J)END DOEND DODO I=1,LWDO J=1,NTF(I,J)=T(I,J)END DOEND DOXF=MATMUL(VF,TF)ENDSUBROUTINE JCB(A,N,EPS,V,B,L)! A:调用时存放实对称矩阵! B:返回时存放矩阵的全部特征值! V:存放特征向量,其中第i列为与第i个特征值相对应的特征向量! EPS:存放精度要求REAL(8),DIMENSION(N,N)::A,VREAL(8),DIMENSION(N)::BREAL(8)::FM,CN,SN,OMEGA,X,YINTEGER::P,QL=1V=0DO I=1,NV(I,I)=1END DO10 FM=0DO I=1,NIF(I/=J.AND.ABS(A(I,J))>FM)THENFM=ABS(A(I,J))P=IQ=JEND IFEND DOEND DOIF(FM<EPS)THENL=1GOTO 15END IFIF(L>100)THENL=0WRITE(12,'("L=",I3,2X,"没有达到精度要求")')L GOTO 15END IFL=L+1X=-A(P,Q)Y=(A(Q,Q)-A(P,P))/2OMEGA=X/SQRT(X*X+Y*Y)IF(Y<0)OMEGA=-OMEGASN=1+SQRT(1-OMEGA*OMEGA)SN=OMEGA/SQRT(2*SN)CN=SQRT(1-SN*SN)FM=A(P,P)A(P,P)=FM*CN*CN+A(Q,Q)*SN*SN+A(P,Q)*OMEGA A(Q,Q)=FM*SN*SN+A(Q,Q)*CN*CN-A(P,Q)*OMEGA A(P,Q)=0A(Q,P)=0DO J=1,NIF(J/=P.AND.J/=Q)THENFM=A(P,J)A(P,J)=FM*CN+A(Q,J)*SNA(Q,J)=-FM*SN+A(Q,J)*CNEND IFEND DODO I=1,NIF(I/=P.AND.I/=Q)THENFM=A(I,P)A(I,P)=FM*CN+A(I,Q)*SNA(I,Q)=-FM*SN+A(I,Q)*CNEND IFEND DOFM=V(I,P)V(I,P)=FM*CN+V(I,Q)*SNV(I,Q)=-FM*SN+V(I,Q)*CNEND DODO I=1,NB(I)=A(I,I)END DOGOTO 10! 将特征值按大小排列15 M=N20 IF(M>0)THENJ=M-1M=0DO I=1,JIF(B(I).LT.B(I+1))THENB1=B(I)B(I)=B(I+1)B(I+1)=B1M=IDO K=1,NV1=V(K,I)V(K,I)=V(K,I+1)V(K,I+1)=V1END DOENDIFENDDOGOTO 20ENDIFEND5.3.6例对1991年5月1日至8月31日123天某纬圈上20个格点的850hPa高度场进行自然正交分解。
PROGRAM MAININTEGER,PARAMETER::P=20 !资料的空间点数INTEGER,PARAMETER::N=123 !资料的时间长度REAL(8),DIMENSION(P,N)::X,XF,ERRORREAL(8),DIMENSION(P)::XV !X的平均值INTEGER::LP=P;LW=3OPEN(12,FILE='F58Z850.DAT')DO IT=1,123READ(12,'(<LP>F8.1)')(X(I,IT),I=1,P)END DOCLOSE(12)XV=0DO I=1,PDO J=1,NXV(I)=XV(I)+X(I,J)END DOXV(I)=XV(I)/NEND DODO I=1,PDO J=1,NX(I,J)=X(I,J)-XV(I)END DOEND DOOPEN(12,FILE='EOF.DAT')CALL EOF(X,P,N,LW,XF)ERROR=X-XFDO I=1,PDO J=1,NXF(I,J)=XF(I,J)+XV(I)END DOEND DOWRITE(12,'(" 恢复场")')WRITE(12,'(<LP>F8.1)')((XF(I,J),J=1,P),I=1,N)WRITE(12,'(" 误差场")')WRITE(12,'(<LP>F8.1)')((ERROR(I,J),J=1,P),I=1,N)CLOSE(12)END输出结果为：特征值：1 2 3 4 5 6 7 8.2488D+06 .3194D+05 .1649D+05 .8222D+04 .4707D+04 .3137D+04 .1538D+04 .9665D+039 10 11 12 13 14 15 16.5258D+03 .4200D+03 .3926D+03 .3108D+03 .2749D+03 .2655D+03 .2395D+03 .1963D+0317 18 19 20.1922D+03 .1148D+03 .1110D+03 .9482D+02解释的方差(%)1 2 3 4 5 6 7 8 9 1078.01 88.02 93.19 95.77 97.25 98.23 98.71 99.02 99.18 99.3111 12 13 14 15 16 17 18 19 2099.44 99.53 99.62 99.70 99.78 99.84 99.90 99.94 99.97 100.00可见第一主分量解释的方差达78.01%，前三个主分量解释的方差达93.19%。