自然正交函数分析(EOF)程序
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
5.3自然正交函数分析(EOF)程序
近年来,自然正交函数(又称经验正交函数)展开在气象上应用比较广泛。这种正交函数展开不象三角函数展开、球函数展开那样有固定的展开形式。它无固定的函数形式,不是事先人为地给定典型场函数,图形是由场本身来决定的,它具有收敛快又能更好地反映出场的基本结构的特征。它可以在有限的区域中进行,既可以取空间不同站点进行分解,也可以对同一站点的不同时间、不同高度的多种要素进行综和分析。因此它在气象中具有广泛的应用,可用于气象要素场分析、大气垂直结构分析、动力模型垂直分层等。
5.3.1功能
计算要素场的自然正交函数分解。
5.3.2方法说明
自然正交函数分解是针对气象要素场进行的,它的基本思想是把包含p个空间点(或p个变量)的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 (5.3.1)
其中F ij 表示第i 个场中的第j 个测点的观测值。
可将(5.3.1)是写为矩阵的形式
VT F = (5.3.2)
式中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.3)
事实上,我们对(5.3.2)式右乘T '可得
V T VT F F ''=' (5.3.4)
因F F '是p p ⨯阶对称阵,其元素为距平变量的交叉积。根据实对称矩阵的分解定理有 V V ΛF F '=' (5.3.5)
其中Λ是F F '矩阵的特征值组成的对角阵,V 是对应的特征向量为列向量组成的矩阵。比较(5.3.4)和(5.3.5)式可知
ΛT T =' (5.3.6)
又根据特征向量的性质有
I V V V V ='=' (5.3.7)
式中为I单位矩阵。显然(5.3.6)和(5.3.7)式满足(5.3.3)式的要求。由
F'矩阵的特征向量求得,而时间函数则可利用(5.3.2)式左此可知空间函数矩阵可从F
乘V'得到,即
T'
=(5.3.8)
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