经验正交函数分解(EOF)方法介绍
(2014-03-31 11:39:02)分类: 数理统计 |
E0F分析方法能够把随时间变化的变量场分解为不随时间变化的空间函数部分以及只依赖时间变化的时间函数部分。空间函数部分概括场的地域分布特点,而时间函数部分则是由场的空间点的变量线性组合所构成,称为主要分量。这些分量的头几个占有原场内空间点所有变量的总方差的很大部分,这就相当于把原来场的主要信息浓缩在几个主要分量上,因而研究主要分量随时间变化的规律就可以代替场的时间变化研究,且可以通过这一分析得出的结果来解释场的物理变化特征。它的优点在于典型场由变量场序列本身的特征来确定,而不是事先人为规定,因而能较好地反映出场的基本结构。这种方法展开收敛速度快,很容易将大量资料信息浓缩集中。它能对有限区域内不规则分布的站点进行分解,且分解的空间结构具有明确的物理意义。Lorenz在20世纪50年代首次将其引入气象和气候研究,现在该方法已在海洋和其他学科中得到了广泛的应用。
注:这是从别人那得到的程序,在此对编此程序的前辈表示感谢!
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),
C ***********************************************
C * N:
C * M:
C * JOB: CONTRAL
PARAMETER;
C
C
*
C
*
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’,
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)
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),
IF(JOB.EQ.0)THEN
GOTO 1
ELSE IF(JOB.EQ.1)THEN
GOTO 2
ELSE
GOTO 3
END IF
DO 10 J=1,N
GOTO 111
GOTO 111
DO 20 I=1,N
DO 20 J=1,N
A(I,J)=0.0
DO 17 K=1,LL
A(J,I)=A(I,J)
WRITE(6,200)
WRITE(6,900)
DO 800 I=1,N
c WRITE(6,801)I,(A(I,J),J=1,N)
c
CALL JACOBI (N,.TRUE.,D,IRT,B,Z,V,A)
DO 400 I=1,N
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
WRITE (6,104)
WRITE(6,500)(IX(I),I=1,N)
WRITE(6,501)(D(I),I=1,N)
DO 502 J=1,N
ILW=IX(J)
DO 503 I=1,N
WRITE (6,105)
DO 600 I=1,N
WRITE(6,601)(V1(I,J),J=1,46)
DO 25 L=1,LL
DO 25 J=1,N
T(L,J)=0.
DO 31 I=1,N
WRITE (6,106)
DO 700 I=1,LL
WRITE(6,701)(T(I,J),J=1,46)
DO 705 I=1,LL
AP=0.0
DO 702 I=1,N
DO 703 I=1,N
AP1=0.0
DO 704 J=1,I
WRITE(6,107)
WRITE(6,108)(H1(I),I=1,N)
RETURN
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
W(J)=W(J)/FLOAT(LL)
DO 30 J=1,N
DO 30 I=1,LL
H(I,J)=F(I,J)-W(J)
RETURN
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
FM(J)=W(J)/FLOAT(LL)
W(J)=0.0
DO 30 I=1,LL
W(J)=SQRT(W(J)/FLOAT(LL))
DO 40 I=1,LL
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
GOTO 10
B(K)=A(K,K)
D(K)=B(K)
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
IF (SM)101,210,101
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)).
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
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
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