当前位置：文档之家 › 正交函数分解(EOF)源代码(Visual Basic 6.0)

正交函数分解(EOF)源代码(Visual Basic 6.0)

  • 格式：docx
  • 大小：18.45 KB
  • 文档页数：8

下载文档原格式

  / 8

合集下载

EOF应用 从数据预处理到详细分析

EOF应用 从数据预处理到详细分析

EOF分析By lqouc 1.什么是EOF，它的作用是什么。

1.1什么是EOF关于EOF 要先从主成分分析说起，主成分分析是多元统计分析中重要的一部分，是一种从多个变量化为少数变量的统计方法，利用多个变量之间相互关系构造一些新的变量，这些新的变量不仅能综合反映原来多个变量的信息，而且彼此之间是相互独立的，同时是按方差贡献大小排列的，这种统计处理方法称为主成分分析。

主成分分析在气象应用中称为经验正交函数（EOF）分解。

1.2E OF的用途对于一个气象要素，我们通常有m个空间点或者台站，有n次观测，这样组成的矩阵中的任意元素就表示了某一空间某一时刻的函数，我们希望能将这样的时空函数分解成空间函数与时间函数两部分的线性组合。

根据主成分的性质，主成分是按其方差贡献大小排列的，而且是相互独立的，那么可以用前几个时间函数与对应的空间函数的线性组合，对原始场做出估计和解释，这就是经验正交函数分解的主要目的。

2.EOF的数据预处理EOF只是个统计学的方法，本身不带有任何物理意义，更不会揣摩作者的意图，所以在数据导入之前需要对数据进行分析和预处理。

以免得到错误的或者不理想的结果。

在此处所说的预处理不是指一般EOF程序中自带的距平或者标准化的处理，虽然这确实有一定的区别。

总之，在做EOF之前，对数据需要有基本的了解，也要对自己的研究目的十分明确。

2.1 数据预处理的必要性例如：想利用EOF 研究极地海平面气压场的年际变化，数据是六十年的月平均的海平面气压格点资料。

首先对手中的资料有基本的判断，月分辨率的资料包含的时间信号的尺度可能有季节内变化、季节变化、年变化、年际变化、年代际变化以及线性趋势。

而我们需要的只是其中的年际变化的信号，所以为了排除干扰必须对数据进行滤波。

这一步是非常有必要的，因为一般来讲，气温、气压、SST这种受太阳辐射影响巨大的要素都具有很强的季节变化，这样的信号远远强于年际变化。

2.2 滤波的方法对于滤波的方法，我们熟悉的有很多，最简单的是做年平均，还有滑动平均、带通滤波、谐波滤波、线性去趋势。

EOF分解程序

EOF分解程序

fid=fopen('HadISST1_SST_1961-1990.txt','r');Num=360;data=zeros(360,180,Num);for i=1:Numaaa=fscanf(fid,'%s',7);data(:,:,i)=fscanf(fid,'%f',[360,180]);endsst1=data(1:90,11:70,1:Num); % 选取所需要区域的数据sst2=data(311:360,11:70,1:Num);sst3=zeros(140,60,Num);sst3(90:-1:1,1:60,1:Num)=sst1;sst3(140:-1:91,1:60,1:Num)=sst2;sst=sst3;for i=1:140for j=1:60for k=1:Numif(sst(i,j,k)==-1000)||((sst(i,j,k)==-32768))页脚内容1sst(i,j,k)=NaN;endendendendsst_area1=zeros(Num,8400); % zeros全零数组for i=1:Num;squ=squeeze(sst(:,:,i)); % 执行该指令后sst数据转换为二维数组sst_area1(i,:)=reshape(squ,1,8400); % 将数据转变为二维endsst_nan=isnan(sst_area1);i=0;for j=1:8400if sum(sst_nan(:,j))==0;i=i+1;sst_region(:,i)=sst_area1(:,j);end页脚内容2end% 求距平~注意季节的变换X=zeros(size(sst_region)); % 学者给的程序for i=1:12X(i:12:Num-12+i,:)=sst_region(i:12:end,:) - repmat( mean(sst_region(i:12:end,:),1) , size(sst_region(i:12:end,:),1), 1); endR=X*X'; % 协方差矩阵R=X*X'是8400*8400的方阵~现定义矩阵R=X'*X是156*156的矩阵[v,d]=eig(R); % 进行EOF分解~因为X'*X与X*X'的秩相同所以特征值相同~d为x的特征值组成的对角阵~v为X*X'的特征向量~[diagonal,I]=sort(diag(d),'descend');v=v(:,I);G=diagonal/sum(diagonal);Gs=0;for i=1:length(G)Gs(i)=sum(G(1:i));if Gs(i)>0.8 break;页脚内容3endendN=length(Gs)%v=fliplr(v); % 矩阵作左右翻转%d=rot90(d,2); % 矩阵上下翻转后再左右翻转(查看生成的对角阵是由小到大排列的~此指令可使其由大到小排列~fliplr(flipud(d))=rot90(d,2))%diagonal=diag(d);spacef=X'*v;for i=1:Num;spacef(:,i)=spacef(:,i)/sqrt(diagonal(i)); % 空间本征函数endtimef=X*spacef; % 时间本征函数sum_d=sum(diagonal);count=0;for i=1:Num;count=count+diagonal(i);G1(i)=count/sum_d; % G1(i)是累积方差贡献率页脚内容4endfor i=1:Num;G2(i)=diagonal(i)/sum_d; % G2(i)是方差贡献率end%**************************************************************************% 将删去的陆地与冰点的填充值补回sst_area2=zeros(Num,8400);sst_area2(:,:)=NaN;i=1;spacef2=spacef';for j=1:8400if sum(sst_nan(:,j))==0;sst_area2(:,j)=spacef2(:,i);i=i+1;endendsst_area3=sst_area2';页脚内容5%**************************************************************************% 画图% subplot(2,1,2) % 将绘图窗口划分为2*1个子窗口，并在第2个子窗口中绘图figure(1)x=1:Num;plot(x,timef(:,1),'g');%ylim([ -80 80 ]);% xlabel('1980-1992年156个月','fontsize',15,'fontname','隶书')ylabel('INDEX','fontsize',12,'fontname','黑体')set(gca,'xtick',(1:6:Num))set(gca,'xticklabel',{'1980','','1981','','1982','','1983','','1984','','1985','','1986','','1987','','1988','','1989','','1990','' ,'1991','','1992','','1993'})title('北太平洋第1模态1980至1992年SST时间序列', 'color', 'k','fontsize',15,'fontname','幼圆')grid onhold off% %subplot(2,1,1)页脚内容6% lon=[130.5:269.5];% lat=[20.5:79.5];% m_proj('Equidistant Cylindrical','lat',[20.5 79.5],'lon',[130.5 269.5]);% m_contourf(lon,lat,rot90(reshape(sst_area3(:,6),140,60)',2),30,'linestyle','none'); % colorbar% m_coast('patch',[.95 .95 .95]);% m_coast('color','k');% m_grid('linestyle','none','tickdir','out','linewidth',1.5);% % xlabel('longitude','fontsize',15,'fontname','comic sans ms');% % ylabel('latitude','fontsize',15,'fontname','comic sans ms');% title(['北太平洋第6模态填色图'],'fontsize',15,'fontname','幼圆');lon=[130.5:269.5];lat=[20.5:79.5];figure(2)m_proj('lambert','lat',[20.5 79.5],'lon',[130.5 269.5]);m_contourf(lon,lat,rot90(reshape(sst_area3(:,1),140,60)',2));页脚内容7% colorbar;m_coast('patch',[1 .85 .7]);%m_elev('contourf',[500:500:6000]);m_grid('box','fancy','tickdir','in');%colormap(flipud(copper));% xlabel('longitude','fontsize',15,'fontname','幼圆');% ylabel('latitude','fontsize',15,'fontname','幼圆');title(['北太平洋第1模态填色图'],'fontsize',15,'fontname','幼圆'); colorbar页脚内容8。

正交函数

正交函数

2.1 用完备正交函数集表示信号将信号分解为正交函数分量的问题与将矢量分解为正交矢量的方法是类似的。

下面我们首先从正交矢量开始讨论，进而引入正交函数集的概念。

为了给出正交函数的概念，并研究正交函数的分解方法，下面我们首先来回顾一下矢量的正交矢量分解。

2.1.1 正交矢量设矢量与矢量的关系如下图所示其中，是它们的夹角，是在上的投影，是用投影来表示时的误差。

由于直角边长小于斜边长，所以根据投影误差可知：如果要用上的矢量来表示，则应选择在上的垂直投影，这时的投影误差最小。

此时(选择垂直投影来表示)，于是，可以求出系数为式中的算子表示求两矢量的内积，定义如下系数C 12表示的是与的近似程度。

当与重合时，=0， =1；随着 增大，C12减小；当 时，C 12，此时与成为相互垂直的矢量，称为正交矢量。

这样，我们就提出这样一个结论：两个矢量是正交的充要条件是它们的内积为0，即<=> 与垂直或正交下面我们来看看，有了正交矢量后，对我们到底有什么好处？对于二维平面上的矢量V 在直角坐标中可分解为x 方向的分量和y 方向上的分量，其中Vx 、Vy 表示x 和y 方向上的正交单位矢量，即为了便于研究矢量分解，把相互正交的两个矢量组成一个二维正交矢量集，在此平面上的任意矢量均可用二维正交矢量集的分量组合来表示。

同样，对于三维矢量V，也可以用一个三维正交矢量集{Vx,Vy,Vz}的分量组合来表示，即V = C1Vx+C2Vy+C3Vz根据此原理，可以把K维空间中的任一矢量分解为K个互相正交的矢量的和。

2.1.2 正交函数下面我们来考虑在区间(t1,t2)内用函数来近似表示，即此时，所选择的C12应使得C12ƒ2(t)与ƒ1(t)之间的均方误差最小。

为求使均方误差最小时的C12，须使。

即交换微分与积分次序，可得上式中的第一项为零，因此（2-2）仿照矢量内积定义，定义两个函数和在区间上的内积为则式(2-2)可以改写为(2-3)可以看出，它与矢量正交分解的系数公式(2-1)很相似。

EOF应用：从数据预处理到详细分析

EOF应用：从数据预处理到详细分析

响因子， 进行简单相关、 复相关和偏相关分析， 确定可能的影响因子。 确定了影响因子之后可以尝试用多元回归分析， 探讨这些因子与研究 要素之间的可预报性。 除了以上提到的分析，还可以根据自己的目的增加分析的内容。 5. 不同类型的 EOF 5.1EOF 本身的变化 对于 EOF 的介绍很多的参考书籍都将其用于时空分离，也就是用在 了空间和时间构成的三维场。但是实际上，我们回归最前面的 EOF 的出处， 可以看出最原本的主成分分析并没有限定要素是时空的函数。 这种方法只是通过引入新变量来达到数组降维的效果。 所以我们可以 在应用中进行多种尝试，只要能在物理上找到合理的解释就没问题。 因为，这终究只是一种数学工具。 举个例子，我们将一个 30 年长度月分辨率的时间序列，写成一个 30*12 的数组，第一维 30 年，第二维是 12 个月，这样以 30 年为我 们通常认为的时间，12 个月为‘空间’ ，进行 EOF 分析，得到的结果 可以揭示不同模态下 12 个月分别在这 30 年中的变化。 除此之外还有很多种用法，在此不再赘述，仅作抛砖引玉。 5.2 多变量 EOF（MV-EOF） EOF 分析时， 不仅会研究某一要素的时空特征， 有时也会研究某现象 的时空特征，而这些现象往往不能用单一的要素来表征，这时候就需 要用到了多变量的 EOF。 例如，研究海洋大陆的季风系统时空变化特征，很可能要考虑到
小时间序列的自由度。3.带通滤波也是常用的方法（本人没用过） ， 其优点是可以选定一定的频率范围，缺点是边界处处理不是很清晰。 4.谐波滤波，以傅里叶函数为基函数对时间序列进行逼近，其优点是 可以较准确的得到选取的频段信号，缺点是选的基函数有局限性，而 且结果和时间序列的长度有关。5.线性去趋势可以去除时间序列的线 性趋势信号，但是需要这一线性趋势通过显著性检验。 2.3 如何合理选定分析对象 上面谈到的是滤波的方法， 但是如果我们的数据是一些大家不熟悉的 数据，我们并不知道它都主要包含何种尺度的信号，也不知道各个主 要尺度信号的强弱，那就需要先对时间序列进行分析。对于时间序列 的分析，我们可以采用 1.谐波滤波，看各个频率的数值大小。2.功率 谱分析，得到显著周期。3.小波分析，同样可以得到时间序列的多尺 度变化特征。 在此，我推荐的方法是结合空间利用方差分析，因为以上的分析我们 都是忽略了空间的影响， 一种要素的时间变化特征是会随着空间变化 的。例如，对中国地区做某一要素的 EOF 分析，得到的结果不能通 过检验（检验的方法，后面再说） ，这个时候我们就需要考虑是否一 些地区的目标信号不强，而另外一些地区目标信号很强，这样的话就 只需要分析目标信号很强的地区，即只对特定区域进行 EOF 分析。 结合空间的方差分析， 首先需要对要素每一个空间点的时间序列进行 滤波，得到各个不同频率的信号（从季节内到线性趋势） 。对每个平 率的信号求方差，得到了各个频率的方差的空间分布。在分析的过程

经验正交函数分解(EOF)

经验正交函数分解(EOF)

图1b 1976 年是个明显的转折点, 在这之前累积曲线基本上呈 上升趋势, 海温以正距平主, 这之后累积曲线呈下降趋势, 海温以负距平为主。19511975 年平均海温距平为12℃, 而1977- 1993 年平均海温距 平为-128℃。这就是说西风 漂流区年平均海温从19511975 年至1977- 1993 年下降 了0148℃。 图1c 1981 年是个明显的转折点, 在这之前累积曲线呈下降趋 势, 海温以负距平为主, 这之 后累积曲线呈上升趋势, 海温 以正距平为主。1951- 1980 年平均海温距平为0.08℃, 1982- 1993 年为0.21℃, 赤道 太平洋年平均海温1981 年后 比1981 年前增加了 b 西风漂流区的年平均海温距平(实线) 和累积曲线(虚线) ;
PDF 文件使用 "pdfFactory Pro" 试用版本创建
0 v v j = ∑ v ki v kj = k =1 1
T i p
i≠ j i= j
i≠ j i= j
性 质
ZiZ
T j
=
∑z
t =1
n
it
z jt
0 = λ i
i , j = 1,2 ,L , m
PDF 文件使用 "pdfFactory Pro" 试用版本创建
三、分解方法
XX
T
= VZZ V
T
T
A = XX
V AV = Λ
T
PDF 文件使用 "pdfFactory Pro" 试用版本创建
T
A为实对称矩阵，根据实对称矩阵分解原 理，一定有 或者
k =1 p
i = 1, 2 , L , m t = 1 , 2 , L , n

正交函数展开

正交函数展开
there is a uniquie sequence of scalars (an) such that
x (1e1 ....... nen ) 0
as
Then (en) is called basis for X. series
k ek
k 1
Which has the sum x is then called the expansion of x
For b.c’s :
c1u1(a) c2u2(a) 0
c1u1(b) c2u2(b) 0
The condition of nontrivial sol. of c1,c2 to be existed if :
u1(a)u1(a) 0 u2 (b)u2 (b)
The Euler Column
x y x y
Here x and y are arbitrary vector in X and is any scalar
Definition of Banach Space A Banach space is a complete normed space.
Lemma (Translation invariance) A metric d induced by a norm on a normed space X satisfies
Normed Space
Definition of Normed Space Here a norm on a vector space X is a real-value function on X whose value
at an x X is denoted by x
x 0 x 0x0
x y
S(x0;r) B~(x0;r) B(x0;r)

正交函数

正交函数
• 正交函數 在高等數學中，函數可視為廣義向量(Generalized Vector) 。 3-Dimension 中，向量的內積性質： (1) (u, v) = (v, u) (2) (ku, v) = k(u, v) ，k 為純量 (3) 若 u = 0 ⇒ (u, u) = 0 若 u ≠ 0 ⇒ (u, u) > 0 (4) (u + v, w) = (u, w) + (v, w)
b0 = 0,
b1 = ±
3, 2
c0 = ±
2 = m3
5 8
Ex. 試證集合{1, cos x, cos 2x, …}在區間[−π , π ] 內正交。
Ex. 試求上例之正交集合的範[Norm] 。
{ } < Ans > 2π , π , π , ......... 。
Ex. 正交集合{1, cos x , cos 2x, ………}的正範集合為何？
。
** 若{φn (x)}在 [a,b]上對權函數 w(x) 成正交，則[a,b]區間內函數 f (x)的廣義傅立葉級數為
f (x) = ∑∞ cnφn (x)
n=0
但係數 cn 為
cn
=
∫b a
f
(
x) w( x)φn
φn
(
x)
2 w
(
x ) dx
其中
φn
(x)
2 w
=
∫
b a
w(
x )φn 2
(
sin 2π x , L
cos 2π x , L
sin 3π x , L
cos π x , L
L 於區間
(−L, L) 內構成一正交集合。

matlab EOF

matlab EOF

一列特征向量值，也称EOF。如λ1 对应的特征向量值称第一个EOF模态， 也就是V 的第一列即EOF1 = V (:, 1)；第λk 对应的特征向量是V 的第k 列， 即EOFk = V (:, k )。 • 计算主成分。将EOF投影到原始资料矩阵X 上，就得到所有空间特征向量对 应的时间系数(即主成分)，即
数据性质与预处理
(1)误差 (2)资料的处理。原始场，距平场，与标准化场 例子：我国160站夏季降水量的EOF分析(图A.17) (3)空间样本点。大范围的空间数据，特别需要注意资料空间代表性。非均匀 场与均匀分布场；空间抽样；面积加权。 北半球1月SLP例子
时空转换
有时空间样本m远大于时间序列长度n，计算m × m矩阵的特征根很困难，可以
1.80 0.60
-1.20 -0.40
46
EOF1 26.1% 4000 3000 Eigenvalue 2000 1000 0 0.04 0.02 0 −0.02 −0.04 −0.06 0 2 4 6 Number 8 10 −0.08
100 50 PC#1 0 −50 −100 1950 1960 1970 1980 1990 2000
-2.20 -4.40
1.80 0.60
-1.20 -0.40
[U,S,V]=svd(X); 得到 U= 0.19 0.98 0.98 -0.19 S= 6.49 0 0 0 0 0 4.23 0 0 0 V= 0.66 -0.49 0.56 0.02 0.67 0.63 -0.73 -0.31 0.53 0.14 0.39 0.03 -0.10 -0.26 -0.02 EOF=U; PC=S*V’; 得到PC= 4.28 -2.07
0.09 -0.32 0.25 0.91 0.06

matlab经验正交函数EOF(转载)

matlab经验正交函数EOF(转载)

练习：利用[E,V]=eig(C)计算矩阵X 的特征向量和主成分%
X=[2 6 1 5 2; 9 4 0 5 4]; X(1,:)=X(1,:)-mean(X(1,:)); X(2,:)=X(2,:)-mean(X(2,:)); 得到X的距平值：X= -1.20 2.80 -2.20 1.80 -1.20 4.60 -0.40 -4.40 0.60 -0.40 %%% co-variance matrix C=X*X’/5; 协方差阵C= 3.76 0.92 0.92 8.24 [EOF,E]=eig(C); % V: eigenvectors; E: eigenvalues PC=EOF’*X; %% reverse the order E=fliplr(flipud(E)) lambda=diag(E); % retain eigenvalues only EOF=fliplr(EOF) PC=flipud(PC) 得到EOF= 0.19 0.98
对角线上的元素为奇异值)，奇异值与特征根成倍数关
1 • 如果矩阵C = n XX T ，C 的特征根为λ，则有
• 如果矩阵C = XX T ，C 的特征根为λ，则有
√ = nλ； √ = λ；
由于该方法是直接对矩阵X 进行分解，所以对内存的要求远小于方法1。计算速度 很快。 两种方法对比练习。
显著性检验
-0.98 0.19
得到特征根E= 8.42 0 0 3.58 得到主成分PC= 4.28 0.15 2.07 -2.82
-4.74 1.31
0.94 -1.65
-0.62 1.10
%%check EOF*EOF’ % = I 检查EOF的正交性得到： 1.00 0 0 1.00 PC*PC’/5 % = lambda 检查PC的正交性得到： 8.42 0.00 0.00 3.58 EOF*PC % =X 可以完全恢复X的距平值： -1.20 2.80 -2.20 4.60 -0.40 -4.40

基于SPEI的华东地区近30年干旱时空特

基于SPEI的华东地区近30年干旱时空特

基于SPEI的华东地区近30年干旱时空特征分布研究袁小军1 徐亚军1 王柳松2 于晓光3（1.南通市测绘院有限公司，江苏 南通 226006；2.河南省军区数据信息室，河南 郑州 450003；3.32023部队，辽宁 大连 116023 ）作者简介：袁小军（1987—），男，汉族，本科，工程师，从事测绘技术质量管理、测绘地理信息应用等工作。

E-mail：****************摘　要：干旱作为全球十大自然灾害之首，对人类社会生产造成的危害不容小觑，对干旱开展有效研究预测具有重要意义。

基于1990～2020年实测气象资料，利用标准化降水蒸散指数，采用经验正交函数分析法等分析检验方法，对比分析华东地区不同时间尺度旱涝趋势在时间和空间两个维度上的变化特征。

结果表明：（1）研究区干旱事件的发生频率呈南高北低的空间分布，但干旱强度在空间上呈北高南低分布；（2）不同时间尺度下研究区南北旱涝变化差异明显，北部地区呈干旱化趋势，而南部相反。

关键词：标准化降水蒸散指数；华东地区；时空分布；干旱；SPEI1 引言干旱是指在一定时间尺度内，由水分收支或水资源供需失衡而造成的缺水现象，通常伴有降水不足和气温升高现象。

当前，干旱已位列全球十大自然灾害之首，而未来的干旱在时长、规模、影响程度及频率上均会远超以往。

尽管通过节能减排降低室温气体排放可降低风险，但低水平的变暖仍能加剧全球大部分地区的干旱程度[1]，因此，干旱研究尤为重要。

干旱研究目前运用较多的指数主要有标准化降水指数（Standardized Precipitation Index，SPI）、帕默尔干旱指数（Palmer Drought Severity Index，PDSI）、标准化降水蒸散指数（Standardized Precipitation Evapotranspiration Index，SPEI）。

降水和温度均为干旱研究的重要因素。

SPI 具有多时间尺度特性，却未考虑温度带来的地表蒸散影响；PDSI 考虑了温度和降水的影响，却无法满足多尺度和重大干旱预测时效性要求；SPEI 指数在延续SPI 多时间尺度优点的前提下，综合考虑气温所带来的水分盈亏影响，满足干旱指数准确性研究，因此成为研究重点，被广泛应用于地区干旱检测、气温要素利用、气候变暖环境下的干旱变化趋势预测等领域。

EOF分析及其应用(教学课件)

EOF分析及其应用(教学课件)
分析表明，南亚夏季风的爆发主要体现在降 水的突然增加和季风雨带的快速推进上，雨 带的时空分布有突变的特点。
第1 模态——降水量的突然增加。 第2 模态——从南向北的快速推进过程。 第3模态——东西分布型态，及在季风爆发
后印度半岛降水快速增加的过程。 第4模态——印度次大陆东海岸降水的准双
周振荡型态。
EOF分析及其应用
中国气象科学研究院
1
一、引 言 经验正交函数（EOF）方法：最早由统计学家
pearson（1902）提出，由Lorenz（1956）引入 气象问题分析中。该方法以场的时间序列为分析对 象，由于对计算条件要求甚高，直到20世纪60年代 后期才在实际工作中得到广泛应用。 近30年来，出现了适合于各种分析目的的EOF分析 方法，如扩展EOF（EEOF）方法，旋转EOF （REOF）方法，风场EOF（EOFW）方法，复变量 EOF（CEOF）方法。
z11 z12 z1n
V
v21
vm1
v22 vm2
v2m
vmm
Z
z 21
zm1
z 22 zm2
z2n
zmn
v j (v1j ,v2 j ,,vmj )T
是第j个典型场，只是空间的函数。
6
第t个空间场可表示为
x1t v11
v12
v1m
x2t
19
第1 模态——降水量的突然增加
20
降水量的第2 模态 -从南向北的快速推进过程
21
降水量的第3 模态
东西分布型态，及在季风爆发后印度半岛降水快速增加 的过程
22
降水量的第4 模态 印度次大陆东海岸降水的准双周振荡型态
23
我国盛夏500 hPa 风场的EOF 分析及 其与大尺度气候异常的关系

经验正交函数分解(EOF)

经验正交函数分解(EOF)
1n
m是空间点，n是时间序列长度。
PDF 文件使用 "pdfFactory Pro" 试用版本创建
气象场的自然正交展开，是将X分解为时间函 数Z和空间函数V两部分,即
X = VZ
或者
xit = ∑ vik z kt = vi1 z 1t +vi 2 z 2t + L + vip z pt
i
特征向量为
vi
PDF 文件使用 "pdfFactory Pro" 试用版本创建
vi = Xu i
λi
V = v1 , v 2 ,L , v m
Z =V X
T
PDF 文件使用 "pdfFactory Pro" 试用版本创建
EOF分析实例
PDF 文件使用 "pdfFactory Pro" 试用版本创建
三、分解方法
XX
T
= VZZ V
T
T
A = XX
V AV = Λ
T
PDF 文件使用 "pdfFactory Pro" 试用版本创建
T
A为实对称矩阵，根据实对称矩阵分解原 理，一定有 或者
PDF 文件使用 "pdfFactory Pro" 试用版本创建
八、时空转换问题 当 m >> n 时，先求出 X T X 的特征值， 然后求
XX
T
的特征向量，这种方法叫时空转换。
令 X T X 的特征值为 λ i ，其特征向量 为 u ， XX T 的特征值也为 λ i ，其
k =1 p
i = 1, 2 , L , m t = 1 , 2 , L , n

气象统计分析与预报经验正交函数分解

气象统计分析与预报经验正交函数分解

实验二经验正交函数分解一、目的和要求：经验正交函数分解（EOF）是统计天气分析中气象要素场最基础的研究模型，是必须理解和掌握的方法之一，是后续课程中许多气象要素场的计算结果的理解的基础理论，也是毕业设计和论文中的基本分析方法。

该方法用个数较少的几个空间分布模态来描述环流形势，而且基本涵盖环流场的信息，既能作为天气分析模型，其方法的延拓又能作为天气预报模型，在实际工作中也有极强的实用意义。

通过该实验，深刻理解气象要素场的统计模型的意义，掌握气象要素场分析的基本方法，为实际预报业务和科研工作打下一定的基础。

二、实验的主要内容：对（0N-90N，60E-120W）850hPa高度场进行经验正交展开（EOF.FOR），输出分析主要参数指标；绘制环流型图和相应的时间系数序列图，并加以分析。

三、步骤：3.1 熟悉资料方法3.1.1 资料提供的资料为NCEP/NCAR 60年（1948年-2007年）逐年1～12月的850hPa高度场资料，资料范围为（90N-90S，0E-360E），网格距为2.5*2.5，纬向格点数为144，经向格点数为73。

资料为NC格式，资料从南到北、自西向东排列，每月为一个记录，按年逐月排放，注意读取方式以及记录长度。

本次实验应用NCEP/NCAR（0N-90N，60E-120W） 58年（1948年-2005年）逐年7月的850hPa高度场资料，纬向格点数为73，经向格点数为37。

3.1.2 方法（经验正交函数分解EOF）EOF（经验正交函数分解）是针对气象要素场进行的，其基本原理是把包含p个空间点（变量）的场随时间变化进行分解。

设抽取样本容量为n的资料.则场中任一空间点i和任一时间点j 的距平观测值ij x 可看成由p 个空间函数ik v 和时间函数kj y (k=1,2,…,p)的线性组合,表示成11221p ij ik kj i j i j ip pj k x v y v y v y v y ===+++∑EOF 功能是从一个气象场多次观测资料中识别出主要空间型及其时间演变规律。

正交函数分解

正交函数分解

正交函数分解
正交函数分解是一种将任意函数表示为一组正交函数的线性组合的方法。

正交函数是指在某个内积空间中，不同函数之间的内积为零，而同一函数的内积为非零的函数。

利用正交函数分解，我们可以将复杂的函数分解为简单的正交函数的和，从而更容易进行数学分析和计算。

这种方法在信号处理、图像处理、量子力学、泛函分析等领域都有着广泛的应用。

一组常见的正交函数包括勒让德多项式、拉盖尔多项式、傅里叶级数等。

在实际应用中，我们可以根据需要选择不同的正交函数来分解所需的函数。

正交函数分解的基本思想是将一个函数表示为一组正交函数的线性组合。

对于任意给定的函数f(x)，我们可以先将其表示为一组基函数φn(x)的线性组合：
f(x)=∑anφn(x)
其中，an为待求系数，φn(x)为一组基函数。

为了满足正交性，我们要求基函数满足以下条件：
∫φn(x)φm(x)dx=0（n≠m）
∫φn(x)φn(x)dx≠0
通过以上条件，我们可以得到an的表达式：
an=∫f(x)φn(x)dx/∫φn(x)φn(x)dx
最终，我们可以将f(x)表示为以下形式：
f(x)=∑∫f(x)φn(x)dx/∫φn(x)φn(x)dxφn(x)
正交函数分解的优点在于，通过寻找合适的基函数，我们可以将复杂的函数变得简单，并且可以方便地进行数学分析和计算。

因此，在实际应用中，正交函数分解也是一种非常重要的数学工具。

EEOF和CEOF

EEOF和CEOF

4.2 扩展经验正交函数分解（EEOF ）经验正交函数的优点在于其能用相对少的函数及时间权重来描述复杂的气象场的变化。

然而，气象场不仅在空间上有高度的相关性，而且在时间上有自相关及交叉相关。

同时利用空间、时间上的相关性压缩资料可能比经典的EOF 分析更有效，也更利于物理解释。

Wear （1982）提出的扩展EOF 分析，其基本方法类似EOF 分析，但资料矩阵包含了几个时间上连续的值。

分解方法：设m 为空间点，n 为时间点。

EEOF(Extended Empirical Orthogonal Function)考虑了时间上的自相关和交叉相关，取1112121221213113(2)2311314134......................................................n m m mn n m n m m mn n m m mn x x x x x x x x x X x x x x x x x x x ---⨯--⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦上面矩阵实际上是由t ，t ＋1，t ＋2三个时刻的资料组成。

第一行到第m 行对应t 时刻的资料，第m ＋1行到第2m 行对应t ＋1时刻，第2m ＋1行到第3m 行对应t ＋2时刻。

T A X X =A 是一个后延相关矩阵，因为它既包含同时刻，又包含后延1个时刻和2个时刻的元素值相乘的和。

A 矩阵为3m 阶，可以求出3m 个特征值123,,...,m λλλ，以及3m 个特征向量123,,...,m v v v 。

其中，每一个特征向量都是3m 维列向量。

由此123(,,...,)m V v v v =T T V X =分解式为X V Z = 分析及应用：EEOF 分析方法与EOF 的分析略有不同。

其特征向量是3m 维，故有1212221223123(.........)t t t i i i im i m i m i m i m i m i m v v v v v v v v v v ++++=对应对应对应所以1v 所对应的m 个空间点上有3张特征向量图，分别对应123,,t t t 时刻。

  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

'*************************************' 全局变量，便于主函数调用。

' VB 6.0 的函数返回的参数偏少，' 使用全局变量在一定程度可以解决这个问题。

'****************************************Public A() As Single ' 协方差/相关系数矩阵APublic V() As Single '特征向量为列组成的矩阵，即空间函数V （EOF）Public T() As Single '时间系数矩阵T（PC）Public B() As Single '特征值λ（E）,按从大到小排列Public GM() As Single '解释的方差(%)（特征向量对X场的累积贡献率）P Public GA() As SinglePublic GB() As Single '个体i特征向量对X场的贡献率ρPublic XF() As Single '模拟结果'********************************************************' 函数名：CovarMat' 函数用途: 计算协方差（相关系数）矩阵' 参数说明：矩阵下标为1：N,从1开始；' X,存放原始观测值，二维实型数组，X(P,P)。

' 返回：计算协方差（相关系数）矩阵。

'*******************************************************Function CovarMat(X() As Single) As Single()Dim XX() As SingleDim P As Integer, N As IntegerDim px As SingleP = UBound(X, 1)N = UBound(X, 2)px = IIf(N > 0, 1 / N, 1)ReDim Preserve XX(1 To P, 1 To P)Dim iAs Integer, j As Integer, k As Integer' 求X乘以X的转置,即A=XXˊFor i = 1 To PFor j = 1 To PXX(i, j) = 0For k = 1 To NXX(i, j) = XX(i, j) + X(i, k) * X(j, k)Next kXX(i, j) = XX(i, j) * pxNext jNext iCovarMat = XXEnd Function'********************************************************' 函数名：EOF' 函数用途: EOF,经验正交分解法(EOF)' 参数说明：矩阵下标为1：N,从1开始；' X,存放原始观测值，二维实型数组，X(P,N)。

' P,整型变量，空间格点数。

' N,整型变量，序列的时间长度。

' XF,返回时存放恢复值，二维实型数组，XF(P,N)。

'*******************************************************Sub EOF(X() As Single, ByRef S() As String)Dim V1() As SingleDim VF() As Single, TF() As SingleDim P As Integer, N As IntegerP = UBound(X, 1)N = UBound(X, 2)ReDimXF(1 To P, 1 To N)ReDimT(1 To P, 1 To N)ReDimA(1 To P, 1 To P)ReDimV(1 To P, 1 To P)ReDimV1(1 To P, 1 To P)ReDimB(1 To P)ReDimGM(1 To P)ReDimGA(1 To P)ReDimGB(1 To P)ReDimVF(1 To P, 1 To P)ReDimTF(1 To P, 1 To N)' 求X的协方差（相关系数）矩阵A = CovarMat(X)Dim iAs Integer, j As Integer, k As Integer, L As Integer' 用Jacobi法求A的特征值和特征向量' 返回时B存放矩阵的全部特征值，V存放特征向量为列组成的矩阵Call Jacobi(A, P, 0.000001, V, B, L)For i = 1 To PGA(i) = 0For j = 1 ToiGA(i) = GA(i) + B(j)Next jNext iFor i = 1 To PGM(i) = GA(i) / GA(P)GB(i) = B(i) / GA(P)Next iFor i = 1 To PFor j = 1 To PV1(i, j) = V(j, i)Next jNext iT = MATMUL(V1, X) '时间系数'输出结果字符串Dim lsAs Integerls = UBound(S)ReDim Preserve S(ls + 2)S(ls) = MA TtoString(B, 1, , "特征值λ（E）")S(ls + 1) = MA TtoString(GB, 1, 100, "个体i特征向量对X场的贡献率ρ")S(ls + 2) = MA TtoString(GM, 1, 100, "解释的方差(%)（特征向量对X场的累积贡献率）P")For i = 1 To PFor j = 1 TolwVF(i, j) = V(i, j)Next jNext iFor i = 1 TolwFor j = 1 To NTF(i, j) = T(i, j)Next jNext iXF = MATMUL(VF, TF) '模拟值End Sub'*******************'函数名：MATtoString'函数作用：矩阵转成字符串输出'参数说明：' MAT，用以存储矩阵数值，最多为2维数组' retS，返回时存储字符串，以换行符（vbcrlf）结尾' nDim，矩阵维数，最大为2' strMatNote,矩阵备注信息，默认为空字符串'********************Function MATtoString(MAT() As Single, _nDim As Integer, _Optional ZoomCoefAs Single = 1, _Optional MatNoteStringAs String = "") As StringDim retString As StringDim N As Integer, i As Integer'如果数组说明为非空串，则添加换行符retString = IIf(Len(MatNoteString) > 0, MatNoteString&vbCrLf, MatNoteString)Select Case nDimCase 1N = UBound(MAT)For i = 1 To NretString = retString& Format(MAT(i) * ZoomCoef, "#0.##") &vbTabNext iretString = retString&vbCrLfCase 2Dim m As Integer, j As IntegerN = UBound(MA T, 1)m = UBound(MA T, 2)For i = 1 To NFor j = 1 To mretString = retString& Format(MAT(i, j) * ZoomCoef, "#0.##") &vbTabNextretString = retString&vbCrLfNext iEnd SelectMATtoString = retStringEnd Function'*************************'函数名：MATMUL'函数作用：矩阵相乘'参数说明：' Mat1，用以存储矩阵1数值' Mat2，用以存储矩阵1数据' 返回：矩阵相乘结果'**************************Function MATMUL(MAT1() As Single, MA T2() As Single) As Single()Dim MatXX() As SingleDim N As Integer, m As Integer, L As Integer, l2 As IntegerDim iAs Integer, j As Integer, k As IntegerN = UBound(MAT1, 1)m = UBound(MAT2, 2)L = UBound(MAT1, 2)l2 = UBound(MAT2, 1)If L <> l2 Then EndReDimMatXX(1 To N, 1 To m)For i = 1 To NFor j = 1 To mMatXX(i, j) = 0For k = 1 To LMatXX(i, j) = MatXX(i, j) + MAT1(i, k) * MAT2(k, j)Next kNext jNext iMATMUL = MatXXEnd Function'********************************************************' 函数名：Jacobi' 函数用途: 用Jacobi法求A的特征值和特征向量' 参数说明：矩阵下标为1：N,从1开始；' A:调用时存放实对称矩阵,A(N,N)' N:矩阵长度' Bx:返回时存放矩阵的全部特征值,B(N)' Vx:特征向量为列组成的矩阵,V(N,N),其中第i列为与第i个特征值相对应的特征向量' EPS:存放精度要求'*********************************************************Private Sub Jacobi(A() As Single, _N As Integer, _EPS As Single, _ByRefVx() As Single, _ByRefBx() As Single, _ByRef L As Integer)' A:调用时存放实对称矩阵' Bx:返回时存放矩阵的全部特征值' Vx:存放特征向量,其中第i列为与第i个特征值相对应的特征向量' EPS:存放精度要求ReDimVx(1 To N, 1 To N)ReDimBx(1 To N)Dim FM As Single, CN As Single, SN As SingleDim Omega As Single, X As Single, Y AsSingleDim P As Integer, Q As IntegerDim iAs Integer, j As Integer, k As IntegerL = 1'初始化特征向量矩阵For i = 1 To NFor j = 1 To NVx(i, j) = 0Next jVx(i, i) = 1Next i'如果计算次数大于给定次数，则跳出循环Do While (L < 100)FM = 0'查找矩阵中非对角元素的最大值，并记录其位置（P,Q)For i = 1 To NFor j = 1 To NIf (i<> j And Abs(A(i, j)) > FM) ThenFM = Abs(A(i, j))P = iQ = jEnd IfNext jNext i'如果最大值小于给定最小值，则跳出循环If (FM < EPS) ThenL = -1Exit DoEnd IfL = L + 1X = -A(P, Q)Y = (A(Q, Q) - A(P, P)) / 2Omega = X / VBA.Sqr(X * X + Y * Y)If (Y < 0) Then Omega = -OmegaSN = 1 + VBA.Sqr(1 - Omega * Omega)SN = Omega / VBA.Sqr(2 * SN)CN = VBA.Sqr(1 - SN * SN)FM = A(P, P)A(P, P) = FM * CN * CN + A(Q, Q) * SN * SN + A(P, Q) * OmegaA(Q, Q) = FM * SN * SN + A(Q, Q) * CN * CN - A(P, Q) * OmegaA(P, Q) = 0 '处理完当前最大非对角元素值，赋值0，下一循环不再计算A(Q, P) = 0 '处理完当前最大非对角元素值，赋值0，下一循环不再计算For j = 1 To 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 IfNext jFor i = 1 To NIf (i<> P Andi<> Q) ThenFM = A(i, P)A(i, P) = FM * CN + A(i, Q) * SNA(i, Q) = -FM * SN + A(i, Q) * CNEnd IfNext iFor i = 1 To NFM = Vx(i, P)Vx(i, P) = FM * CN + Vx(i, Q) * SNVx(i, Q) = -FM * SN + Vx(i, Q) * CNNext iFor i = 1 To NBx(i) = A(i, i)Next iLoop' 将特征值按大小排列m = NDo While (m > 0)j = m - 1m = 0For i = 1 To jIf (Bx(i) <Bx(i + 1)) ThenB1 = Bx(i)Bx(i) = Bx(i + 1)Bx(i + 1) = B1m = iFor k = 1 To NV1 = Vx(k, i) Vx(k, i) = Vx(k, i + 1)Vx(k, i + 1) = V1Next kEnd IfNext iLoopEnd Sub。