5・3自然正交函数分析(EOF)程序近年來，自然正交函数(乂称经验正交函数)展开在气象上应用比较广泛。

这种正交函数展开不彖三角函数展开、球函数展开那样有固定的展开形式。

它无固定的函数形式，不是事先人为地给定典型场函数，图形是由场木身来决定的，它具有收敛快又能更好地反映岀场的基木结构的特征。

它可以在有限的区域屮进行，既可以取空间不同站点进行分解，也可以对同一站点的不同吋间、不同高度的多种要素进行综和分析。

因此它在气彖中具有广泛的应用，可用于气象要素场分析、大气垂直结构分析、动力模型垂直分层等。

5. 3.1功能计算要素场的自然正交函数分解。

5. 3. 2方法说明口然止交函数分解是针对气彖要素场进行的，它的基本思想是把包含P个空间点(或P个变量)的n个时次的观测场随时间进行分解，即将某一区域的气象要素场序列Fq (i=l, 2,・・・,p；j=l,2,…,n,即p个空间点的n个时次的观测资料)分解成相互正交的时间函数与相互正交的空间函数的乘积Z和，常把空间函数VW看作典型场，时间函数看作典型场的权重系数，则不同时间的要素场是若干个典型场按不同权重线性叠加的结果，各个场之间的差别就在于各典型场的系数不同。

则气象耍素场可以表示为PEj =》％tkj = Vig+Vj2t2j+・・・+Viptpj (5. 3. 1)k=l英中Fq表示第i个场中的第j个测点的观测值。

可将(5.3.1)是写为矩阵的形式F =VT( 5 . 3 . 2 )式中F为pxn阶的均值为0的资料阵，V为pxp阶的空间函数阵，卩为pxn阶的时间函数阵。

由于V和0是根据场的资料阵F进行分解而得到的，分解的函数没有固定的函数形式，因而称为“经验”的，另外，我们还要求这种分解具有“正交”性，即要求满足下式PVk V, =X v ik v n =0 (kHl)i=1(5. 3. 3 )n兀齐-£，kjtij =0 (k H 1)冃事实上，我们对(5. 3. 2)式右乘厂可得FF =VTTV r( 5 . 3 . 4 )因FF'是pxp阶对称阵，其元素为距平变量的交义积。

根据实对称矩阵的分解定理有FF =VAV f(5. 3. 5 )其小A是FF'矩阵的特征值组成的对角阵，V是对应的特征向量为列向量组成的矩阵。

比较(5.3. 4)和(5. 3. 5 )式可知TT r = A(5 . 3. 6 )乂根据特征向量的性质有VV =W=I式中为I单位矩阵。

显然(5 . 3 . 6 )和(5 . 3. 7 )式满足(5 . 3 . 3 )式的要求。

由此町知空间函数矩阵可从FF'矩阵的特征向量求得，而时间函数则可利用(5. 3. 2)式左乘V'得到，即T =VF( 5. 3. 8)5. 3. 3子程序语句CALL EOF (X, P, N, XF)5. 3. 4哑元说明X——输入变量，二维实型数组，大小为PxN,存放原始观测值。

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

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

XF——输出变量，二维实型数组，人小为PxN,存放恢复值。

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 ⑻,DIMENSIONS,N)::TREAL ⑻,DIMENSION(P)::B,GM,GAREAL(8),DIMENSIONS,LW)::VFREAL(8),DIMENSION(LW,N)::TF! 求X乘以X的转置，即A=XX 'DO 1=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 1= 1 ,PGA(I)=()DO J=1,IGA(I)=GA(I)+B(J)END DOEND DODO 1=1,PGM(I)=GA(I)/GA(P)END DODO 1=1,PDO J=1,PV1(I,J)=V(J,I)END DOEND DOT=MATMUL(V1,X)WRITE(12；(H特征值”))WRITE( 12；(<P>I10)')(I,I= 1 ,P)WRITE( 12；(3X,<P>D 10.4)*)BWRITE(12；(H解释的方差(%)u nWRITE( 12；(<P>I7),)(I,I= 1 ,P)WRIT E( 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")1)WRITE( 12；(<P>F10.4),)((T(I,J),J=1 ,P),I= 1 ,N)DO 1=1,PDO J=1,LWVF(I,J)=V(IJ)END DOEND DODO 1=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),D1MENS1ON(N,N)::A,VREAL(8),DIMENSION(N)::BREAL ⑻::FM,CN,SN,OMEGA,X,YINTEGER:: P,QL=1V=0DO 1= 1 ,NV(I,I)=1END DO 10 FM=0DO 1=I,NDO J=1,NIF(I/=J.AND.ABS(A(I,J))>FM)THENFM=ABS(A(IJ))P=1Q=JEND IFEND DOEND DOIF(FMvEPS)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 1=1,NIF(I/=P.AND.1/=Q)THENFM=A(I,P)A(I,P)二FM*CN+A(I,Q)*SNA(I,Q)=-FM*SN+A(I,Q)*CNEND IFEND DODO 1=1,NFM=V(I,P)V(I,P)=FM*CN+V(I,Q)*SNV(I,Q)=-FM*SN+V(I,Q)*CNEND DODO 1=1,NB(I)=A(I,1)END DOGOTO 10!将特征值按人小排列15 M=N20 IF(M>0)THENJ=M-1M=0DO 1= 1 ,JIF(B(I).LT.B(I+1))THENB1=B(I)B(I)=B(I+1)B(I+1)=BIM=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©),DIMENSIONS,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)J= 1 ,P) ENDDOCLOSE(12) XV=O DO 1=1,P DO J=1,NXV(I)=XV(I)+X(I,J) END DO XV(I)=XV(I)/N END DO DO 1=I,P DO J=1,NX(I,J)=X(I,J)-XV(I) END DO END DOOPEN( 12,FILE-EOF.DAT) CALL EOF(X,P,N,LW,XF) ERROR=X-XF DO 1=1,P DO J=1,NXF(1,J)=XF(I,J)+XV(1) END DO END DOWRITE(12；(H 恢复场？)WRITE( 12；(<LP>F8.1 )*)((XF(I, J), J= 1 ,P),I 二 1 ,N) WRITE(12；(H 误差场WRITE( 12；(<LP>F8.1 )')((ERROR(I,J),J 二 1 ,P),I= 1 ,N) CLOSE(12) END输出结果为: 特征值：.39260+03 . 3108D+03 .2749D+03 . 2655D+03 . 2395D103 .19630)0317 18 19 20.1922D+03 ・ 1148D+03 . 1110D+03 . 9482D+02 解释的方差(%)1 23 4 5 6 7 89 10 78.018& 02 93. 19 95. 7797. 25 9& 2398.7199. 02 99. 1899.3111 1213 14 15 16 17 18 19 20 99. 44 99. 53 99. 62 99. 70 99. 78 99.8499. 90 99. 94 99. 97 100. 00可见第一主分量解释的方差达78.01%,前三个主分呈解释的方差达93.19%。