加载中…
个人资料
  • 博客等级:
  • 博客积分:
  • 博客访问:
  • 关注人气:
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
正文 字体大小:

经验正交函数分解(EOF)方法介绍

(2014-03-31 11:39:02)
分类: 数理统计

E0F分析方法能够把随时间变化的变量场分解为不随时间变化的空间函数部分以及只依赖时间变化的时间函数部分。空间函数部分概括场的地域分布特点,而时间函数部分则是由场的空间点的变量线性组合所构成,称为主要分量。这些分量的头几个占有原场内空间点所有变量的总方差的很大部分,这就相当于把原来场的主要信息浓缩在几个主要分量上,因而研究主要分量随时间变化的规律就可以代替场的时间变化研究,且可以通过这一分析得出的结果来解释场的物理变化特征。它的优点在于典型场由变量场序列本身的特征来确定,而不是事先人为规定,因而能较好地反映出场的基本结构。这种方法展开收敛速度快,很容易将大量资料信息浓缩集中。它能对有限区域内不规则分布的站点进行分解,且分解的空间结构具有明确的物理意义。Lorenz在20世纪50年代首次将其引入气象和气候研究,现在该方法已在海洋和其他学科中得到了广泛的应用。

        经验正交函数方法是一种使用特征技术的统计方法,能够将变量场时间和空间变化分离,并用尽可能少的模态表达出主要的时间和空间变化。它是变量场空间相关结构或多变量综合信息的有用方法。但是在有些情况下,应用EOF方法最大可能地反映了所有m个原变量的变化和信息,接下来的每个EOF也都尽其所能地表达所有m个原变量的尚未被其他EOF表示出的部分。这样,当变量个数很多(m很大),而且变量间的相关只在局部变量之间时,经验正交函数方法就会过分强调m个变量的整体相关机构,而使重要的局部相关结构被掩盖。

    这里贴不了数学公式,等能贴了,我再把具体的原理贴出来。我把用fortran写的程序代码贴上来吧。不过代码需要做些修改才能使用的。

注:这是从别人那得到的程序,在此对编此程序的前辈表示感谢!

THIS IS A PROGRAM FOR CACULATING EMPIRICAL ORTHOGONAL FUNCTION 
PROGRAM EOF 
PARAMETER(N=48,M=160,JOB=1) 
DIMENSION F(N,M),T(N,M),IX(M),V(M,M),A(M,M), 
     & D(M),W(M),H(N,M),B(M),Z(M),H1(M),V1(M,M),FM(M) 
C *********************************************** 
C * N:   SAMPLE SIZE                            * 
C * M:   NUMBER OF STATION OR GRID POINT        * 
C * JOB: CONTRAL PARAMETER;                     * 
          WHEN JOB=0 ORIGINAL DATA;              * 
C *      WHEN JOB=1 DEPATURE DATA;              * 
C *      WHEN JOB=2 STANDARDIZE DATA.           * 
C * F(N,M): ORIGINAL DATA FIELD                 * 
C * D(M): EIGENVALUE                            * 
C * V(M,M): CHARACTERISTIC VECTOR               * 
C * T(N,M): TIME COE FFICIENT                   * 
C *********************************************** 
OPEN(11,FILE=’F:\paperprogram\monthlylh\summerrain.dat’, 
                  status=’old’) 
READ(11,*)((F(I,J),J=1,M),I=1,N) 
OPEN(6,FILE=’rain2.dat’,STATUS=’NEW’) 
CALL SEOF(M,N,JOB,F,H,T,IX,V,A,D,W,FM,B,Z,H1,V1) 
close(11) 
close(6)

STOP 
END

SUBROUTINE SEOF(N,LL,JOB,F,H,T,IX,V,A,D,W,FM,B,Z,H1,V1) 
DIMENSION F(LL,N),T(LL,N),IX(N),V(N,N),A(N,N),D(N), 
     & W(N),H(LL,N),B(N),Z(N),H1(N),V1(N,N),FM(N) 
IF(JOB.EQ.0)THEN 
GOTO 1 
ELSE IF(JOB.EQ.1)THEN 
GOTO 2 
ELSE 
GOTO 3 
END IF 
  1 DO 10 I=1,LL 
DO 10 J=1,N 
  10 H(I,J)=F(I,J) 
GOTO 111 
  2 CALL DEP(N,LL,F,H,W) 
GOTO 111 
  3 CALL NOR(N,LL,F,H,FM,W) 
  111 CONTINUE 
DO 20 I=1,N 
DO 20 J=1,N 
A(I,J)=0.0 
DO 17 K=1,LL 
  17 A(I,J)=A(I,J)+H(K,I)*H(K,J) 
A(J,I)=A(I,J) 
  20 CONTINUE 
WRITE(6,200) 
  200 FORMAT(30X,’—————EOF—————’)  
WRITE(6,900) 
  900 FORMAT(//20X,’****** CORVARINCE ******’) 
DO 800 I=1,N 
c WRITE(6,801)I,(A(I,J),J=1,N) 
801 FORMAT(2X,I3/(2X,8F15.4)) 
  800 CONTINUE 
CALL JACOBI (N,.TRUE.,D,IRT,B,Z,V,A) 
DO 400 I=1,N 
  400 IX(I)=I 
DO 401 I=1,N-1 
DO 402 J=I+1,N 
IF(ABS(D(I)).GE.ABS(D(J))) GOTO 402 
W1=D(I) 
D(I)=D(J) 
D(J)=W1 
K=IX(I) 
IX(I)=IX(J) 
IX(J)=K 
  402 CONTINUE 
  401 CONTINUE 
WRITE (6,104) 
  104  FORMAT (//20X,’******* EIGENVALUE VALUE ******’) 
WRITE(6,500)(IX(I),I=1,N) 
  500 FORMAT(/20X,’MAX-MIN’/(12X,20I5)) 
WRITE(6,501)(D(I),I=1,N) 
  501 FORMAT(2X,8F15.5) 
DO 502 J=1,N 
ILW=IX(J) 
DO 503 I=1,N 
  503 V1(I,J)=V(I,ILW) 
  502 CONTINUE 
WRITE (6,105) 
  105  FORMAT (//20X, ‘****** CHARACTERISTIC VECTOR ******’) 
DO 600 I=1,N 
WRITE(6,601)(V1(I,J),J=1,46) 
  601 FORMAT(10X,46F10.2) 
  600 CONTINUE 
DO 25 L=1,LL 
DO 25 J=1,N 
T(L,J)=0. 
DO 31 I=1,N 
  31 T(L,J)=T(L,J)+F(L,I)*V1(I,J) 
  25 CONTINUE 
WRITE (6,106) 
  106 FORMAT (///20X, ‘****** TIME COEFFICIENT ******’) 
DO 700 I=1,LL 
WRITE(6,701)(T(I,J),J=1,46) 
  701 FORMAT(10X,46F10.2) 
  700 CONTINUE 
DO 705 I=1,LL 
  705 CONTINUE 
AP=0.0 
DO 702 I=1,N 
  702 AP=AP+D(I) 
DO 703 I=1,N 
AP1=0.0 
DO 704 J=1,I 
  704 AP1=AP1+D(J) 
  703 H1(I)=1.0*AP1/AP 
WRITE(6,107) 
  107 FORMAT(//20X,’****** ACCUMULATE PROPORTION ****** ‘) 
WRITE(6,108)(H1(I),I=1,N) 
  108 FORMAT(/10X,10F10.5) 
RETURN 
END

SUBROUTINE DEP(N,LL,F,H,W) 
DIMENSION F(LL,N),H(LL,N),W(N) 
DO 10 J=1,N 
W(J)=0.0 
DO 20 I=1,LL 
  20 W(J)=W(J)+F(I,J) 
W(J)=W(J)/FLOAT(LL) 
  10 CONTINUE 
DO 30 J=1,N 
DO 30 I=1,LL 
H(I,J)=F(I,J)-W(J) 
  30 CONTINUE 
RETURN 
END

SUBROUTINE NOR(N,LL,F,H,FM,W) 
DIMENSION F(LL,N),H(LL,N),FM(N),W(N) 
DO 10 J=1,N 
W(J)=0.0 
DO 20 I=1,LL 
  20 W(J)=W(J)+F(I,J) 
FM(J)=W(J)/FLOAT(LL) 
W(J)=0.0 
DO 30 I=1,LL 
  30 W(J)=W(J)+(F(I,J)-FM(J))**2 
W(J)=SQRT(W(J)/FLOAT(LL)) 
DO 40 I=1,LL 
  40 H(I,J)=(F(I,J)-FM(J))/W(J) 
  10 CONTINUE 
RETURN 
END 
SUBROUTINE JACOBI(N,EV,D,IRT,B,Z,V,A) 
DIMENSION D(N),B(N),Z(N),V(N,N),A(N,N) 
LOGICAL EV 
IF(. NOT. EV) GOTO 20 
DO 10 K=1,N 
DO 10 L=1,N 
IF (K-L) 5,6,5 
   6 V (K,K)=1.0 
GOTO 10 
   5 V(K,L)=0.0 
  10 CONTINUE 
  20 DO 30 K=1,N 
B(K)=A(K,K) 
D(K)=B(K) 
  30 Z(K)=0.0 
IRT=0 
DO 200 I=1,50 
SM=0.0 
N1=N-1 
DO 40 K=1,N1 
K1=K+1 
DO 40 L=K1,N 
  40 SM=SM+ABS(A(K,L)) 
IF (SM)101,210,101 
  101 IF(I-4)50,60,60 
  50 TRESH=0.2*SM/(FLOAT(N)*FLOAT(N)) 
  GOTO 70 
  60 TRESH=0.0 
  70  DO 180 K=1,N1 
K1=K+1 
DO 180 L=K1,N 
G=100.0*ABS(A(K,L)) 
IF (I. GT. 4. AND. ABS(D(K))+G. EQ. ABS(D(K)). 
     AND. ABS(D(L))+G.EQ. ABS(D(L))) GOTO 170 
IF (ABS(A(K,L)). LE . TRESH) GOTO 180 
H=D(L)-D(K) 
IF (ABS(H)+G.EQ.ABS(H)) GOTO 80 
THETA=0.5*H/A(K,L) 
T=1.0/(ABS(THETA)+SQRT(1.0+THETA*THETA)) 
IF (THETA. LT. 0.0) T=-T 
GOTO 90 
  80 T=A(K,L)/H 
  90    C=1.0/SQRT(1.0+T*T) 
S=T*C 
H=T*A(K,L) 
Z(K)=Z(K)-H 
Z(L)=Z(L)+H 
D(K)=D(K)-H 
D(L)=D(L)+H 
A(K,L)=0.0 
KM1=K-1 
IF (KM1) 120,120,100 
100 DO 110 J=1,KM1 
G=A(J,K) 
H=A(J,L) 
A(J,K)=C*G-S*H 
110 A(J,L)=S*G+C*H 
120 L1=L-1 
IF (L1-K1)140,300,300 
300 DO 130 J=K1,L1 
G=A(K,J) 
H=A(J,L) 
A(K,J)=C*G-S*H 
130  A(J,L)=S*G+C*H 
140 L1=L+1 
IF (L1. GT. N) GOTO 150 
DO 145 J=L1,N 
G=A(K,J) 
H=A(L,J) 
A(K,J)=C*G-S*H 
145 A(L,J)=S*G+C*H 
150 IF (.NOT. EV) GOTO 160 
DO 155 J=1,N 
G=V(J,K) 
H=V(J,L) 
V(J,K)=C*G-S*H 
155 V(J,L)=S*G+C*H 
160 IRT=IRT+1 
GOTO 180 
170 A(K,L)=0.0 
180 CONTINUE 
DO 200 K=1,N 
D(K)=B(K)+Z(K) 
B(K)=D(K) 
200 Z(K)=0.0 
210 RETURN 
END

0

阅读 收藏 喜欢 打印举报/Report
  

新浪BLOG意见反馈留言板 欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

新浪公司 版权所有