1、 研究两个场相互关系的研究两个场相互关系的SVDSVD方法方法数学上的SVD分解及其性质SVD在气象学中的应用前面学习的EOF(或PCA)方法，主要用于分析单个单个气象要素场的时空变化结构。大气科学研究中，还经常需要分析两个两个气象要素场的时空变化之间的联系，可采用什么方法来实现？方法一：多变量方法一：多变量EOF(MV-EOF)仍然使用EOF方法，只是将不同的气象要素场在空间维空间维依序排列在同一个资料阵X内；方法二：方法二：EOF 结合相关或回归分析结合相关或回归分析对某一要素场X作EOF，取出某一模态的时间序列，与另一要素场Y作相关或回归分析，可以考察X的该模态与Y场的关系。111222
2、111212111211222212221212=,.=.nnnmnmmm nmmnyyyyyxxxxxxxxyyyyxXYMMMMMMMM设任意两个气象要素场，记为X与Y，分别有m1和m2个空间点。方法三：典型相关分析方法三：典型相关分析将x=x1,x2,xm1T的m1个空间变量通过线性组合构成一个新变量u，同时，也将y=y1,y2,ym2T的m2个空间变量通过线性组合构成一个新变量v，使得新变量u与v之间的相关系数达到极大。方法四：奇异值分解方法四：奇异值分解(SVD)方法。方法。能够找到成对出现的空间分布型，代表两场的相互联系着的空间结构。SVD其实是一个数学名词，一种矩阵运算。气象上研
3、究两个场相互关系的方法之所以习惯上称之为“SVD”，是因为：该方法的核心运算是矩阵的SVD分解。SVD(Singular Value Decomposition)的数学定义的数学定义对于任意一般实矩阵(行数和列数可以不等)，例如m1行m2列的矩阵C，可以分解为：12112212T()()()()mmmmmmmmC=UV000其中，120000=00r 111111111121212221212.=,.,.mmmmmm muuuuuuuuuUu uuMMMM222222211121212221212.=,.,.mmmmmm mvvvvvvvvvVvvvMMMMr min(m1,m2)列向量uj和
4、vj都是单位化单位化的；uj称为”左奇左奇异向量异向量”，vj称为“右奇异向量右奇异向量”，U和V都是正交矩阵；j称为“奇异值奇异值”，1 2 r 0T0,1,ijijiju u当时当时T0,1,ijijijv v当时当时12112212T()()()()mmmmmmmmCC=UV矩阵 的奇异值分解也可写为：000TTTT11 12221=.rrrrkkkkCuvu vu vu v12()mmTjju vC其中，每个都是和 形状相同行列 的矩阵。上式使我们联想到多变量资料阵X(m行n列)的EOF：1122()()()=+.+mmm nm mm nXLYl yl yl y12T212T211()
5、EOF;(1)/()=jnjjjttnnjjjtjjkjjjtkmmmnmnyyyjTjjCXulXXvyvyyXX2阵(行列)的SVD分解阵行 列 的分解就相当于 （的单位化的特征向量）(1行列)类似于行 列，只是将进行了单位化，相当于就相当于 的模，若对计算特征值（这时），那么比较就相当于VS因此，对C进行SVD分解也就同时得到了C的EOF结果。通过以上对比分析了解到，C的SVD分解奇异值j的平方就是CCT的特征值，C的左奇异向量uj就是CCT的单位化的特征向量。实矩阵SVD的性质性质性质1：TkkkkkkC uvC vuTTT11TT1T=(=rrjjjjjjjjkrkjjjkkkjkk
6、kkkjkCu vCv uuC uv u uvUu uCvu性质1证明：由得上式等号两遍都同时右乘，得上式用到了 的正交性质：只有的那一项不为0,等于1)同理可证性质性质2：2T2T(1),kkkkkkCCuC Cv是的第 个特征值，对应的特征向量为左奇异向量(2)同时，也是的第 个特征值，对应的特征向量为右奇异向量关于第(2)条：虽然CTCCCT，但他们有相同的特征值，不同的特征向量。(曾在第五章的时空转换内容中讨论过)。关于第(1)条:在比较SVD与EOF时已经理解。TTTTTTTT,()(),=()kkkkkkkkkkkkkuCCCCuuCC C C uC uC CC uvv如果和是的特
7、征值和特征向量，则有()等式两边左乘，得：因此，也是的特征值，对应的特征向量为性质一，单位化后即为性质性质3：如果矩阵C是实对称阵实对称阵(m阶)，那么，左、右奇异向量都等于C的特征向量，奇异值就是C的特征值。T()()()()TT()()()()TTTTTTT0SVD00000,000=000000mmm mm mm mm mm mm mm mkCCC=UVC=VUCCU=VCVVV CV V UV V V VV VV证明：的分解的表达式：将上式转置得：因为所以，把 矩阵左乘，右乘，得：上式是使实对称阵对角化的关系式，因此，奇异向量也是 的kC特征向量，奇异值是对应的 的特征值。即：这时，C
8、的特征值的平方就是CCT(或CTC)特征值，并且C与CCT(或CTC)拥有相同的特征向量。性质性质4：C矩阵的所有元素的平方和，等于所有奇异值的平方和：12222111=mmrijkijkcC如何来理解？根据性质2：C的奇异值k的平方就等于S=CCT的特征值k，即：2kk根据主成分的性质2：C中m1个变量的方差之和等于C的所有主成分(共有r个)的方差之和，即 S=CCT的对角线元素之和（CCT的迹）等于CCT的r个特征值之和：因此得证：C矩阵的所有元素的平方和等于所有奇异值的平方和。11221111()mmmriiijkiijksc 12221111()=mmrrijkkijkkc 112T2
9、21111211()kmmmriiijkiijkmscmCC倘若是的特征值，那么就会有：SVDSVD方法在气象学中的应用方法在气象学中的应用通常在气象学中,SVD的用法是：两个气象要素场，x和y，空间点数分别为m1和m2,时间点数都为n，计算它们的标准化资料阵标准化资料阵如下：12122211111212122111212122221)1)(22(=,.=.nnmmmnmnnmmmnnnmxxxxyyyyyyyyyxxxxxXYMMMMMMMM计算x和y的交叉协方差阵协方差阵 C=(XYT)/n，(或者C=XYT)，C的形状为m1行m2列，对协方差阵对协方差阵C进行进行SVD分解。分解。121
10、2112212TT()()()()()()1=mnn mmmmmmmmmnCXY=UV000X和Y的协方差阵也可用C=YXT（或者(YXT)/n）来计算，这时C的形状为：m2行m1列，这时称Y为左场，X为右场。这时，X被称为左场，Y为右场单位不同的两个场进行协方差运算，通常都要先标准化。以下的讨论中，令X为左场，Y为右场，它们的标准化资料阵如所示：1212112212TT()()()()()()01=00mnn mmmmmmmmmnCXY=UV120000=00rLLMM OML111111111121212221212.=,.,.mmmmmm muuuuuuuuuUu uuMMMM22222
11、2211121212221212.=,.,.mmmmmm mvvvvvvvvvVvvvMMMM空间分布：空间分布：U和V的每一列都是一个空间分布函数；每个场内部的空间型相互正交。左奇异向量左奇异向量uj表示第j模态左场的空间分布；右奇异向量右奇异向量vj表示第j模态右场的空间分布。111222111212111211222212221212=,.=.nnnmnmmm nmmnyyyyyxxxxxxxxyyyyxXYMMMMMMMM时间序列：时间序列：把原观测场X投影到uj，把Y场投影到vj，即可得到第j模态左右场的时间序列时间序列，记为aj和bj ：TTjjjjau Xbv Y左场的时间序列：
12、右场的时间序列：aj与bj的形状都是1行n列。下面考察第j模态两个时间序列的协方差TTTTT11Cov(,)()()jjjjjjjjjxyjnnaba ba bu Xv Yu C v和 的协方差TT000 xyC=UVUV由奇异值分解公式,等号两边同时左乘，同时右乘，得：T=jxyjjju C v取出该式的第 列，即：T0=00 xyU C VCov(,)=jjjjjjaba b因此，和 的协方差等于第 模态的奇异值：第1模态左右场的时间系数的协方差协方差达到最大，第2模态次之,因此，SVD分解的各模态是根据每对时间序列的协方差的大小来排列的。所以，SVD也被称为“最大协方差分析最大协方差分析
13、(MCA)”。以上分析的是同一对SVD模态的两个时间序列之间的关系，同理易证，不属于同一对的时间序列之间的协方差为0.这一性质可用式子表示为：TT,Cov(,)=0kjkjxykkjkjkjka bu C vu u当时，当时kkkCvu也可利用性质：证明。同一场内部的不同模态时间系数之间的协方差(如 Cov(a1,a2)也等于0，请自行证明。22221111Cov(,)(,)1111jjjjjjjnnnnjtjtjtjtttttrababnnnnababab要衡量两个奇异向量场的密切程度，通常还要对和计算相关系数：第j模态两个时间序列的相关系数q 相关系数r(aj,bj)总是大于0的，因为j0
14、q两时间序列的协方差j是按模态从大到小来排列的，它们的相关相关系数系数r(aj,bj)是否也按从大到小排列？按从大到小排列？不是!第一模态的左右场时间序列的相关系数不一定是最大的。q 此相关系数（一般都比较高）不宜使用简单相关系数的检验方法来进行显著性检验，因为它本质上是一种复相关系数。空间型的表示左右两场协同变化的空间型由左、右奇异向量表示，但是，由于uj和vj是单位化的，因此，其数值的大小与格点数有关，缺少物理意义。在分析SVD结果时，通常，将SVD空间型的结果uj和vj都变换成相关系数相关系数分布图。同类相关系数对第j模态来说：左场X的每个格点的序列可以与左场的时间系数aj做相关系数，得
15、到的相关系数空间分布图称为“左场左场的同类的同类(homogeneous,或称同性、齐次)相关系数图相关系数图”。右场Y的每个格点的序列可以与右场的时间系数bj做相关系数，得到的相关系数空间分布图称为“右场右场的同类相关系数图的同类相关系数图”。TTTT22221111111()(,)11111jjjxxjjnnnnjtjtjtjtttttnnnaaaannnnXaX u XXX uC ur X a左场X的同类相关系数（m1行1列）为：这里，设X场的各格点上的序列都进行了标准化,所以标准差为1。同理，右场Y的同类相关系数（m2行1列）可表示为：21(,)1yyjjnjttbnC vr Y b异
16、类异类(heterogeneous，或或“异性异性”、“非齐次非齐次”)相关系数相关系数q左场左场的异类相关系数：的异类相关系数：左场每个格点观测的时间序列(X的某一行)与右场的时间系数bj的相关系数。q右场右场的异类相关系数：的异类相关系数：右场每个格点观测的时间序列(Y的某一行)与左场的时间系数aj的相关系数。TTTT2222211111111()(,)111111jjjxyjjjjnnnnnjtjtjtjtjttttttnnnbbbbbnnnnnXbX v YXY vC vr X bu左场X的异类相关系数（m1行1列）为：同理，右场Y的异类相关系数异类相关系数r(Y,aj)(m2行1列)
17、也与右奇异向量右奇异向量呈常数倍关系：TkkkkkkC uvC vu性质：21(,)1jjjnjttanr Y av所以，左场的异类相关系数左场的异类相关系数与左奇异向量左奇异向量呈简单的常数倍关系。各模态对总协方差平方和的贡献率SVD的性质的性质4：“C的所有元素平方和等于所有奇异值的平方和”，协方差阵C的所有元素平方和|Cxy|2就是总的协方差平方和总的协方差平方和(Squared Covariance)，可表示两个场整体的相关程度：12222111=mmrxyijkijkcC221SCFjjrkkSCF是Squared Covariance Fraction的缩写前K个模态的累积(Cum
18、ulated)协方差平方和的贡献率为：2121CSCFKkkKrkk注意 这里的SCF(或CSCF)是指的对两场总的协方差平方和两场总的协方差平方和的贡献，而不是对两场两场总方差总方差的贡献。可定义第j模态对总的协方差平方和总的协方差平方和的贡献率为：第j模态对左右场各自的方差贡献率前面已知：第1模态左场(X)的空间型(m1行1列)为：u1，时间序列(1行n列)为a1=u1T X；利用X第1模态这对空间和时间函数可对原X场作出估计（类似EOF的估计，空间向量乘以时间序列），记为X(1)(m1行n列)：X(1)=u1a1其中的第t列表示第t时刻的估计场：u1a1t；观测到的实际X场记为：xt,所
19、以，第t时刻的误差向量(E的第t列)可写为t：1121 1=tttttm tauxM误差矩阵可写为：111111212122212.=.nnmmm nxxxxxxxxxXMMMM现在要考察第1模态对左场X的方差贡献率，即：X(1)的总方差占X的总方差的比例。要想直接计算推导X(1)的总方差不容易，我们从另一侧面来分析：考察误差的总方差。11111121212221211.=.=nnmmm n uaX，MMMM11122T11111T1 11 112TTTT11111111222T1111T2221111111111()()11111=mmnnnititttittitntttttnttttttt
20、tnttttttmnntttittttitnnnaanaaanaaanaxannn uxuxu uu xx ux xx xx x122111111()=nmnnittittxannTT111T11,1tttau xx uu u利用:12112211111()1=mnititntttxnaanXa就是左场 的的总方差，是时间序列 的方差。同理可证，取前两模态对X进行估计时，对X场估计的总误差方差为：X的总方差减去第一模时间序列a1的方差再减去第二模时间序列a2的方差。因此，第j模态对左场总方差的贡献即为：第j模态左场时间序列aj的方差。12111mnititn 总的误差方差为 阵各行方差之和：上
21、式表明：第一模态对X场的估计的总误差等于：X场的总方差减去时间序列a1的方差，因此可认为，时间序列 a1的方差就是第一模态对左场的解释方差，即X(1)的总方差；a2的方差就是第二模态对左场的解释方差，即X(2)的总方差；因此，第j模态对左场X的总方差的贡献率VF(Variance Fraction)可用以下两者的比值来计算：1212111VF1njttxjmnititanxn 同理，第j模态对右场Y的总方差的贡献率可用以下两者的比值来计算：2212111VF=1njttyjmnititbnyn 注意注意:“对左右各自场的方差方差贡献率(VFxj或VFyj)”与“总的协方差协方差平方和的贡献率(
22、SCFj)”是不同的概念:qSCF反映时间系数的协方差，总是按照从大到小来排列，第一模态的SCF最大；q而VFxj或VFyj反映的是“SVD模态在各自场各自场的变率中所起作用的大小”，不一定按照从大到小的顺序来排列。SVD的显著性检验通过SVD分解，从两个样本资料中识别出了时间系数之间协方差最大、次大的一对对空间型。但是，它是否说明两个场之间确确实实存在联系？任意两个毫无关系的随机场任意两个毫无关系的随机场，用它们的样本资料做SVD，也能计算出一对对空间型和时间系数。所以，SVD的结果是否可信，还需进一步检验。SVD的显著性检验通常使用“蒙特卡罗（Monte Carlo）方法”来实现。“蒙特卡
23、罗”方法又被称为“随机模拟”或“统计模拟”方法。当我们对某一统计量的概率密度分布无法推导出解析解时，就借助大量的计算机随机试验结果对其进行估计。（1）把其中的一个气象资料场（例如Y）随机打乱时间次序，然后与另一个场（X）作协方差并对协方差阵进行SVD分解，记录下每个模态的|Cxy|2和SCF；（2）重复以上步骤100次，即：重新把Y在时间维打乱次序，对协方差阵进行SVD分解，记录下每个模态的|Cxy|2和SCF；（3）以上得到了100个|Cxy|2；对于每个模态（例如第k模态），有100个SCFk，将他们按从大到小从大到小的顺序排列。（4）检验规则：左右场整体相关性的检验：用真实的X和Y的协方
24、差作SVD分解得到的|Cxy|2与随机试验得到的100个|Cxy|2进行比对，如果大于第5个|Cxy|2，则认为在0.05的显著性水平下，XY两场整体上存在相关性。SVD第k模态的检验：对某一模态(如第k模态)来说，也是将该模态的SCFk与100个SCFk比对，如果大于第5个，表明SCFk通过了信度为0.05的显著性检验。12222111=mmrxyijkijkcC221SCFjjrkk检验对象可以有两个：两场整体上有无相关性(|Cxy|2);第k模态是否显著(SCFk)。蒙特卡罗检验的重要依据：如果要检验的对象(|Cxy|2或SCFk)是显著的，那么，它应该大于绝大多数随机场所得到的结果。S
25、VD的实际应用举例Wallace等(1992)研究了冬季太平洋地区SST（左场）与500hPa高度场（右场）的SVD分解。时间长度为从1946年12月到1985年12月，n=39。两个资料场均做了标准化处理。所采用的海温资料在太平洋地区共有m1=157个格点，高度场资料共m2=125个格点。右图为SVD分解的第一对奇异向量第一对奇异向量，异类相关系数分布图。该图反映出：当海温分布形态为45 N，160 W处有一冷(暖)异常中心、热带地区为暖(冷)异常中心时，相应的北太平洋高度场异常型为“太平洋-北美”遥相关型(PNA)。各统计量：第一模态对协方差平方和的贡献率为：SCF1=52%；第一模态左右场两个时间系数(图略)的相关系数为0.81；第一模态解释SST场（左场）方差的12%；解释高度场（右场）方差的17%；
