跳到主要内容

气象数据分析方法合集

普通相关系数

计算方法

为了度量两个变量之间的线性关联程度,常用的指标是普通相关系数,即 Pearson 相关系数。 任意两个变量 x 和 y ,样本(时间长度)为 n,其标准差分别为 σx\sigma_xσy\sigma_y 则它们之间的 Pearson 相关系数的计算公式是:

r=1ni=1n(xixˉσx)(yiyˉσy)r=\frac{1}{n}\sum_{i=1}^{n}(\frac{x_i-\bar{x}}{\sigma_x})(\frac{y_i-\bar{y}}{\sigma_y})

另外一种计算公式是:

r=i=1n(xixˉ)(yiyˉ)i=1n(xixˉ)2i=1n(yiyˉ)2=SxySxxSyyr=\frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum_{i=1}^n(x_i-\bar{x})^2\sum_{i=1}^{n}(y_i-\bar{y})^2}}=\frac{S_{xy}}{S_{xx}S_{yy}}

r 可以客观的度量两个因子之间的关联程度。| r | \leq 1 ,正相关表示两者之间变化性质相同,负相关表示两个变量变化性质相反。

两个变量之间的相关关系可以大致表示为:

y=a+bxy=a+bx

其中 a 是截距,b 是回归系数, y 是气温, x 是降水量。 利用最小二乘法。

b=i=1n(yiyˉ)(xixˉ)i=1n(xixˉ)(xixˉ)=SxySxxb = \frac{\sum_{i=1}^{n}(y_i-\bar{y})(x_i-\bar{x})}{\sum_{i=1}^{n}(x_i-\bar{x})(x_i-\bar{x})} =\frac{S_{xy}}{S_{xx}}

截距 a 可以由下式计算得到:

a=i=1nyinb(i=1nxi)n=yˉbxˉa=\frac{\sum_{i=1}^{n}y_i}{n}-b\frac{(\sum_{i=1}^{n}x_i)}{n}=\bar{y}-b\bar{x}

SxyS_{xy} 表示 xy 之间的协方差, SxxS_{xx} 表示 x 与自己的协方差,即方差。根据相关系数的定义,可以得到:

r=i=1n(xixˉ)(yiyˉ)i=1n(xixˉ)2i=1n(yiyˉ)2=Sxyσxσyr=\frac{\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y} )}{\sqrt{\sum_{i=1}^n(x_i-\bar{x})^2\sum_{i=1}^{n}(y_i-\bar{y})^2}} = \frac{S_{xy}}{\sigma_x\sigma_y}

结合前几式可以得到:

b=rσyσxb=r\frac{\sigma_y}{\sigma_x}

解释的方差比例 ( 方差解释率或方差贡献率 ) 计算:

计算公式

rs=1回归模型计算出值的方差样本值原有的方差r_s=1-\frac{回归模型计算出值的方差}{样本值原有的方差}

在主成分分析和因子分析中,方差解释率表示提取的主成分/因子对原有变量的解释能力,方差解释率越大,解释能力越强,越能体现原始变量的关键影响因素,提取的主成分或因子越有效。

主成分分析中的方差贡献率则为,该特征值/总特征值之和。

在回归中,决定系数也可以称为方差解释率,他代表了自变量对因变量所解释或决定的比率,决定系数越逼近 1 ,拟合效果或解释效果越好。

显著性检验

做显著性检验要判断 样本与我们所知的总体所做的假设之间的差异是 纯属机会变异 还是 由我们所做的假设与总体的真实情况之间不一致所引起的。 显著性检验 是针对我们对总体所做的假设做检验,其原理就是“小概率事件实际不可能性原理” 来接受或否定假设。

抽样实验会产生抽样误差,对实验资料进行比较分析时,不能凭借两个结果(平均数等等)的不同就作出结论,而是要进行统计学分析,鉴别出两者差异是抽样误差引起的,还是由特定的实验处理引起的。

未检验相关系数有别于 0 的置信度,可以用 Z 统计量检验

Z=n32ln1+r1rZ = \frac{\sqrt{n-3}}{2}ln\frac{1+r}{1-r}

t 统计量检验

t=rn21r2t=\frac{r\sqrt{n-2}}{\sqrt{1-r^2}}

自由度的估计:

[!info] 自由度 自由度是指当以样本的统计量来估计总体参数时,样本中独立或能自由变化的数据的个数,称为该统计量的自由度自由度计算公式:样本个数(k)-样本受约束的个数(k)-1

简单估计:随机样本数减 2 ,即 v = n - 2 实际上气候变量的一个突出特点就是具有红噪声谱,即不同时间的数据之间不是完全独立的(不是随机的)。

气候变量某一时刻的状况对后面的状况是有影响的。因此,时间序列的有效自由度要比要比n-2要小。这会影响对相关系数信度的估计和假设结论的判断。

很多气候变量有很强的持续性或者很高的自相关,例如海温。因此进行相关系数的显著性检验时,需要首先对时间序列的有效自由度进行估计。

估计有效自由度的方法有很多。红噪声时间序列的自相关系数随落后时间步长减少,自相关系数越大则独立样本数(有效自由度)越小。有效自由度(v)和样本数(时间序列长度n)之间的关系:

vn=12ln[r(Δt)]\frac{v}{n}=-\frac{1}{2}ln[r(\Delta t)]

vn=[1r(Δt)2][1+r(Δ)2]\frac{v}{n}=\frac{[1-r(\Delta t)^2]}{[1+r(\Delta )^2]}

滞后相关

如果两个时间序列包含有行波,或者任何随时间传播的信号,需要考虑它们地位相差。

如果正好是同位相,则计算地相关系数接近为+1,如果位相差为180度,则相关系数为-1,当位相差为90度,则相关系数接近0.

为解决上述问题,需要计算滞后相关系数,即考虑信号的时间差。同理,也可以通过计算时间序列之间的不同时间滞后相关,可以发现它们之间的可能的位相差。

相关系数的优缺点

优点:是定量描述变量之间的关系,简单直观,容易解释。 缺点

  • 数据中包含两组或更多组不同性质的数据时,不能混在一起同时处理,需要对不同的组分别对待。
  • 如果用反映线性关系的相关系数来描述非线性关系是不合适的。
  • 相关系数对异常值非常敏感。

相关场的信度检验

气候数据不仅在时间上前后有联系,在空间上也有很大的连续性。即空间上也不是完全独立的,因此,空间场的相关是否显著还需要其他要求。但目前缺乏严格的检验方法

秩相关系数

计算方法

Spearman 秩相关系数,鲁棒性比较好,对异常值不是很敏感。

任意两个变量xix_iyiy_i ,样本(时间长度)为n,分别将xix_iyiy_i按各自大小排序,各自序号分别记为DxiD^i_xDyiD^i_y ,那么序号之间的普通相关系数就是Spearman秩序相关系数(rsr_s 即:

rs=i=1n(DxiDxˉ)(DyiDyˉ)i=1n(DxiDxˉ)2i=1n(DyiDyˉ)2r_s=\frac{\sum_{i=1}^n(D_x^i-\bar{D_x})(D_y^i-\bar{D_y})}{\sqrt{\sum_{i=1}^n(D_x^i-\bar{D_x})^2\sum_{i=1}^n(D_y^i-\bar{D_y})^2}}

可以简化为:

rs=16i=1n(DxiDyi)2n(n21)r_s=1-\frac{6\sum_{i=1}^n(D_x^i-D_y^i)^2}{n(n^2-1)}

信度检验

秩相关系数的显著性检验可以利用Z统计量来检验。

Z=rs1n1Z=\frac{r_s}{\sqrt{\frac{1}{n-1}}}

总结

在实际工作中,如果怀疑资料序列中存在异常数据,没法进行排除和订正的时候,用秩相关系数进行分析更为稳妥。

谐波分析

谐波分析又称调和分析,指用三角函数来拟合数字信号或数字序列。根据拟合函数可以对不同的信号周期,位相及振幅的情况进行了解。

简单谐波分析

气候变量变量序列中包含多种时间尺度的变化,因此经常需要考虑如何反映不同时间尺度的变化。

高阶谐波分析

简单谐波分析都是考虑只取一个周期循环的情况,如果序列中还有其他更短的周期,就需要用高阶谐波进行拟合。数学上任何一个信号或者数据序列,都可以表示一系列的不同周期谐波的和。这种情况,可以理解为将原始数据中的不同周期信号进行了分解。

分解出来的各个谐波只是从数学方法得到的,是否具有物理意义需要根据实际情况来判断,也与具体的研究对象和研究目的有关。

功率谱分析

客观定量地检测时间序列中地显著周期,最常用地方法是功率谱。功率谱地计算常用的算法有两种,一是直接计算,二是落后自相关方法。

直接计算法

就是利用谐波分析方法,计算不同阶数的谐波振幅,即Ck2C_k^2 ,振幅大表示能量强,因此也称功率谱(也称功率谱值和功率谱密度)。

f=wk2π=knf=\frac{w_k}{2\pi}=\frac{k}{n} 周期   p=1f=nk周期~~~p=\frac{1}{f}=\frac{n}{k}

通常可以用功率谱图来直观表示各周期的振幅或方差贡献,横坐标是频率(f)或者周期(p)。由于频率和周期随k不是线性变化的,所以有时也用谐波阶数(k)做横坐标,纵坐标可以是Ck2C_k^2 ,或者是标准化处理后的Ck2C_k^2

落后自相关方法

对一个时间序列先求其不同落后时间步长(τ\tau)的自相关或者自协方差,然后对自相关或者自协方差函数进行谐波分析,以此来检测周期。

如果一个时间序列 有5年周期,那么落后步长为5、10、15时,自相关或者自协方差函数就会周期性地出现极大值,用谐波分析就能很好地检测出来。

落后步长长的话,会提高对低频部分的检测分辨率,但是计算落后相关使用的资料会变少而降低资料的自由度从而影响结果的可靠性。

理论分析指出落后步长可以取n10n3\frac{n}{10}-\frac{n}{3},实际计算中常常使用τ=n3\tau=\frac{n}{3} 左右。计算出来的功率谱值如果有峰值,那么对应的周期可能明显,是否显著还需进行显著性检验。根据所分析的原时间序列的特性判断是红噪声谱还是白噪声谱,然后分别用不同的方法进行检验。

计算方法

  • 决定最大落后步长τ=M\tau=M
  • 计算落后自相关系数
R(τ)=1nτt=1nτ(ytyσy)(yt+τyσy)R(\tau)=\frac{1}{n-\tau}\sum_{t=1}^{n-\tau}(\frac{y_t-\overline{y}}{\sigma_y})(\frac{y_{t+\tau}-\overline{y}}{\sigma_y})

其中τ\tau =0,...,M。

  • 计算功率谱粗谱密度
S^k=BkM[R(0)+2τ=1M1R(τ)cos(πkτM)+R(M)cos(πk)]\hat{S}_k=\frac{B_k}{M}[R(0)+2\sum_{\tau=1}^{M-1}R(\tau)cos(\frac{\pi k\tau}{M})+R(M)cos(\pi k)]

其中系数BkB_k 为 ^ir95i9

Bk={1,k=1,...,M112,k=0,MB_k=\left\{\begin{array}{ll} 1, & k=1,...,M-1 \\ \frac{1}{2},& k=0,M\\ \end{array}\right.
  • 计算平滑功率谱密度值,上式中S^k\hat{S}_k 为功率谱粗谱密度值,通常还有对其进行平滑处理来消除随机噪声的影响,以便得到比较光滑和平稳的谱密度值。常用的平滑方法有Barlett窗方法(矩形或者三角形窗)。
  • 信度检验。上面得到的谱密度值如果有极大值,表示对应的频率和周期信号较强。是否显著需要进行显著性检验。

先要对原时间序列的噪声谱性质进行判断。 给定置信度α=0.05\alpha=0.05 有判据Rα=1+tαn2n1R_{\alpha}=\frac{-1+t_{\alpha}\sqrt{n-2}}{n-1} ,如果原时间序列(yty_t)的一阶自相关系数R(1)【R(1)为落后一年的相关系数】比RαR_{\alpha} 大,则为红噪声谱,反之为白噪声谱。 tαt_{\alpha}α\alpha 信度水平和v=n-2 时的t分布值。

检验不同周期的显著性还需要计算出红色噪音谱的95%置信度的上限值:

Sr0.05=Srχ0.052(2nM2)/M=Sk1R(1)21+R(1)22R(1)cos(πkM)χ0.052(2nM2)/M\begin{equation} \begin{split} S_r^{0.05}&=S_r\frac{\chi^2_{0.05}}{(2n-\frac{M}{2})/M} \\ &=S_k\frac{1-R(1)^2}{1+R(1)^2-2R(1)cos(\frac{\pi k}{M})}\frac{\chi ^2_{0.05}}{(2n-\frac{M}{2})/M} \end{split} \end{equation}

如果原时间序列是白噪声谱,则95%置信度的白噪声上限为:

Sr0.05=Skχ0.052(2nM2)/MS_r^{0.05}=S_k\frac{\chi^2_{0.05}}{(2n-\frac{M}{2})/M}

这里χ0.052\chi^2_{0.05} 表示显著水平α=0.05\alpha=0.05 及自由度v=(2nM2)/Mv=(2n-\frac{M}{2})/M 时的χ2\chi^2 (可以查表得到)。不管是红噪声谱检验还是白噪声谱检验,如果SkS_k 超过SrS_r的95%信度上限阈值(Sr0.05S_r^{0.05}) ,则该周期是显著的。 每个k对应的周期是Pk=2MkP_k=\frac{2M}{k} , 当k等于最大落后步长时,对应的周期是2;当k=1时对应的周期是2倍最大落后步长,即所能表示的最长周期。 如果功率谱值有多个显著极大值,表示有多个显著周期。相邻的多个周期中常取中间的最突出的那个峰值。

注意事项

功率谱是对全局(即整个时间序列)频率或者周期的估计。如果某一时期某一周期比较明显,而其他时期该周期不明显,则功率谱分析可能反映不出来。 对此可以分时间段来研究周期随时间的可能变化。

功率谱分析 中对不同频率段的检测是不连续的,分辨率也是不同的。

落后自相关方法 ,功率谱值对落后步长的变化也是比较敏感的。落后步长变化通常会影响分析结果,特别是低频周期部分的结果。

资料准备,最好先对资料进行标准化处理。然后进行去掉其趋势。如果有趋势的话,会影响谱密度的分布特征,使大部分的能量集中在低频部分,降低对高频变化信号的检测能力。

交叉谱分析

交叉谱又称互谱(Cross spectral analysis)。用途是计算两个时间序列是否有显著的周期,两者之间各种周期的关联程度,以及相位差。

对于两个时间序列yty_txtx_t数据,进行交叉谱分析的步骤是:

  • 计算各自即相互之间的落后相关系数Rxx(τ)R_{xx}(\tau),Ryy(τ)R_{yy}(\tau),Rxy(τ)R_{xy}(\tau)Ryx(τ)R_{yx}(\tau)。其中τ=1,...,M\tau=1,...,M。 最大落后步长M的取法与功率谱一致;
  • 计算x和y各自的功率谱粗谱密度及平滑谱密度 。 计算公式在 [[数据处理/气象数据分析方法合集#^ir95i9]],得到Sx(k)S_x(k)Sy(k)S_y(k)
  • 计算x和y的粗协谱
P^xy(k)=BkM{Rxy(0)+τ=1M1[Rxy(τ)+Ryx(τ)]cos(πkτM)+Rxy(M)cos(πk)}\hat{P}_{xy}(k)=\frac{B_k}{M}\{R_{xy}(0)+\sum_{\tau=1}^{M-1}[R_{xy}(\tau)+R_{yx}(\tau)]cos(\frac{\pi k \tau }{M})+R_{xy}(M)cos(\pi k) \}

k=0-M;其中

Bk={1,k=1,...,M112,k=0,MB_k=\left\{\begin{array}{ll} 1, & k=1,...,M-1 \\ \frac{1}{2},& k=0,M\\ \end{array}\right.

然后对粗协谱进行平滑处理以消除随机噪声。以下是一个常用的二项系数平滑公式:

{Pxy(0)=12P^xy(0)+12P^xy(1)Pxy(k)=14P^xy(k1)+12P^xy(k)+14P^xy(k+1)k=1,...,M1Pxy(M)=12P^xy(M1)+12P^xy(M)\begin{equation} \left\{\begin{array}{l} P_{xy}(0)=\frac{1}{2}\hat{P}_{xy}(0)+\frac{1}{2}\hat{P}_{xy}(1) & \\ P_{xy}(k)=\frac{1}{4}\hat{P}_{xy}(k-1)+\frac{1}{2}\hat{P}_{xy}(k)+\frac{1}{4}\hat{P}_{xy}(k+1) & k=1,...,M-1 \\ P_{xy}(M)=\frac{1}{2}\hat{P}_{xy}(M-1)+\frac{1}{2}\hat{P}_{xy}(M) & \end{array}\right. \end{equation}

其中k对应的周期是Pk=2MkP_k=\frac{2M}{k}

  • 计算x和y的正交谱
Q^xy(k)=1Mτ=1M1[Rxy(τ)Ryx(τ)]sin(πkτM)\hat{Q}_{xy}(k)=\frac{1}{M}\sum_{\tau=1}^{M-1}[R_{xy}(\tau)-R_{yx}(\tau)]sin(\frac{\pi k \tau}{M})

k=1,...,Mk=1,...,M,Q^xy(0)=Q^xy(M)=0\hat{Q}_{xy}(0)=\hat{Q}_{xy}(M)=0 。同样需要对粗正交谱进行平滑以得到平稳的正交谱Q^xy(τ)\hat{Q}_{xy}(\tau) :

{Qxy(0)=12Q^xy(0)+12Q^xy(1)=12Q^xy(1)Qxy(k)=14Q^xy(k1)+12Q^xy(k)+14Q^xy(k+1)k=1,...,M1Qxy(M)=12Q^xy(M1)+12Q^xy(M)=12Q^xy(M1)\left\{\begin{array}{l} Q_{xy}(0)=\frac{1}{2}\hat{Q}_{xy}(0)+\frac{1}{2}\hat{Q}_{xy}(1)=\frac{1}{2}\hat{Q}_{xy}(1) & \\ Q_{xy}(k)=\frac{1}{4}\hat{Q}_{xy}(k-1)+\frac{1}{2}\hat{Q}_{xy}(k)+\frac{1}{4}\hat{Q}_{xy}(k+1) & k=1,...,M-1 \\ Q_{xy}(M)=\frac{1}{2}\hat{Q}_{xy}(M-1)+\frac{1}{2}\hat{Q}_{xy}(M)=\frac{1}{2}\hat{Q}_{xy}(M-1) & \end{array}\right.
  • 计算凝聚谱。x和y的凝聚谱值COxy2(k)CO^2_{xy}(k)由下式计算: COxy2(k)=Pxy2(k)+Qxy2(k)Sx(k)Sy(k)CO^2_{xy}(k)=\frac{P_{xy}^2(k)+Q^2_{xy}(k)}{S_x(k)S_y(k)}

k=1,...,Mk=1,...,M

  • ** 凝聚谱值 ( COxy2(k)CO^2_{xy}(k) )如果有极大值,**
  • 就表示 x 和 y 在该相同的频率 k2M\frac{k}{2M} ,即周期 2Mk\frac{2M}{k}
  • 上都有较强的能量。是否显著需要检验。构造统计量:
F=(v1)COxy2(k)1COxy2(k)F=\frac{(v-1)CO^2_{xy}(k)}{1-CO^2_{xy}(k)}

v=(2nM2)v=(2n-\frac{M}{2}) ,给定显著水平α\alpha 查 F 分布表中第一自由度v1=2,第二自由度v2=2(v-1)的阈值FαF_{\alpha} ,如果凝聚谱值超过FαF_\alpha ,则对应的周期2Mk\frac{2M}{k} 是显著的。同时,x和y中该周期对应的信号之间的位相差也可以求出,x超前y的位相ϕ(k)\phi(k)为:

ϕxy(k)=arctan(Qxy(k)Pxy(k))\phi_{xy}(k)=arctan(\frac{Q_{xy}(k)}{P_{xy}(k)})

式中ϕxy(k)\phi_{xy}(k)以弧度表示,如果换算成时间,则为ϕxy(k)Mkπk\frac{\phi_{xy}(k)M}{k\pi k}

时间序列滤波分析

气候要素的时间序列中包含多种时间尺度变化,而很多时候我们希望能保留某些频率的成分,而去除掉其他的成分。

前面的谐波分析是一种简单的滤波过程。如果只有半年的资料,就不能通过谐波来得到低频年变化信号。

滤波的方法有很多种,按滤波的实现方法可以分为对称权重滤波递归滤波。按用途可以分为低通滤波器高通滤波器带通滤波器以及带阻滤波器

气候研究中最常用的滤波器都是对称权重的数字滤波器,因此滤波工作的一个重要内容是设计合适的数字滤波器,得到相应的滤波权重系数(w)。

频率响应

气候研究中最常用的滤波器都是对称权重的数字滤波器(w),用权重系数对原时间序列x进行滤波可以得到所需要的时间序列即:

yt=k=LLwkxt+ky_t=\sum_{k=-L}^{L}w_kx_t+k

权重系数wk=wkw_k=w_{-k} ;其长度是2L+1。

对于给定的频率f,时间序列在滤波前后可能的变化包括该频率的振幅强度位相。 如果权重系数是对称的,则 位相 是不变的,但其 振幅强度 可能会改变。

如果用谐波来表示y和x的话,对应频率f的谐波振幅的谐波振幅分别Cy(f)C_y(f)Cx(f)C_x(f) ,那么其比值:

H(f)=Cy(f)Cx(f)H(f)=\frac{C_y(f)}{C_x(f)}

就是该频率的响应,H称为频率响应函数。类似上述例子的对称权重系数滤波器,在信号处理中也称有限脉冲响应(FIR)滤波器。其影响在时间上是有限的,如1-2-1滤波器,就影响3个相邻的数据。这类滤波器的频率响应函数可以表达为:

H(f)=w0+2k=1Lwkcos(2kπfΔt)H(f)=w_0+2\sum_{k=1}^{L}w_kcos(2k\pi f\Delta t)

f为频率,wkw_k 是第k个权重系数,w0w_0 是中心权重函数,Δt\Delta t 是资料时间步长,通常资料是等时间步长的,可以当成1。如是月分辨率,其步长基本单位是1月,如果是年分辨率的资料,Δt\Delta t就是1年。这样上式就可以简化为:

H(f)=w0+2k=1Lwkcos(2kπf)H(f)=w_0+2\sum_{k=1}^{L}w_kcos(2k\pi f)

从1-2-1滑动平均的频率响应函数看,在f=0处即线性趋势上,滤波后完全不变,对高频部分f=0.5即周期为2的分量则完全消除掉了。该滤波器可以看成是一个低通滤波器。

几种简单滑动平均

二项式系数滤波器的滤波器的权重系数由二项式系数求得:

Bk=(nk)=n!k!(nk)!;k=0,1,...,nB_k=\begin{pmatrix} n\\k \end{pmatrix} =\frac{n!}{k!(n-k)!};k=0,1,...,n

如果要计算3个二项系数,则n=3-1=2,B0=2!0!×2!=1B_0=\frac{ 2!}{0!\times 2!}=1 ;B1=2!1!×1!=2B_1=\frac{ 2!}{1!\times 1!}=2 ;B2=2!2!×0!=1B_2=\frac{ 2!}{2!\times 0!}=1 这显然就是1-2-1滤波器(相应的权重为14\frac{1}{4},24\frac{2}{4},14\frac{1}{4}),对于n=4时,系数为1,4,6,4,1,权重为:116\frac{1}{16},416\frac{4}{16},\frac{6}{16}$$\frac{4}{16},116\frac{1}{16}

随着权重系数的增加,对高频部分的削弱就更明显。

高斯滤波器权重系数遵从正态分布。

设计滤波器

给定截断频率fsf_s ,希望f>fsf>f_s 的部分完全去掉,而对ffsf \le f_s 的部分完全保留。这个理想的滤波器的频率响应函数为:

H(f)={1,0ffs0,fs<f12H(f)=\left\{\begin{matrix} 1,& 0 \le f \le f_s \\ 0,& f_s < f \le \frac{1}{2} \end{matrix}\right.

EOF分析

经验正交函数分析方法(empircial orthogonal function,EOF),也称为特征向量分析(eigenvector analysis),或者主成分分析(principal component analysis ,[[PCA]]) ,是一种分析矩阵数据中的结构特征,提取主要数据特征量的一种方法。 地学数据分析中通常特征向量对应的是空间样本,所以也称空间特征向量或者空间模态;主成分对应的是时间变化,也称时间系数。

原理与算法

原理: [[EOF]]

显著性检验

可以证明

i=1mXi2ˉ=k=1mλk=k=1mPCk2ˉ\sum_{i=1}^{m} \bar{X_i^2} =\sum_{k=1}^{m}\lambda_k =\sum_{k =1}^{m} \bar{PC_k^2}

这说明矩阵X的方差大小可以简单的用的特征根的大小来表示。λ\lambda 越高说明其对应的模态越重要,对总方差的贡献越大。第k个模态对总的方差解释率为:

λki=1mλk×100%\frac{\lambda_k}{\sum_{i=1}^{m}\lambda_k}\times100\%

即使是随机数或者虚假数据,放在一起进行EOF分析,也可以将其分解成一系列的空间特征向量主成分。因此,实际资料分析中得到的空间模态是否是随机的,需要进行统计检验。 North等的研究指出,在95%置信度水平下特征根的误差:

Δλ=λ2N\Delta \lambda =\lambda \sqrt{\frac{2}{N^*}}

λ\lambda是特征根,NN^* 是数据的有效自由度。将λ\lambda 按顺序依次检查,标上误差范围。如果前后两个λ\lambda之间误差范围有重叠,那么他们之间就没有显著区别。

结果展示

通常情况下,主成分是有单位的,即反映的是矩阵X的单位,而空间特征向量是无量纲的。实际应用中常常对EOF分析得到的主成分特征向量进行标准化处理得到新的PC^*和EOF^*

PC(k)=PC(k)λkPC^*(k)=\frac{PC(k)}{\sqrt{\lambda_k}} EOF(k)=EOF(k)λkEOF^*(k)=EOF(k)\sqrt{\lambda_k}

或者简单的将PC标准化,使得其平均值为0,标准差为1。再将它与原始资料矩阵X进行回归分析,这样就得到PC变化一个单位时,变量X对应的响应的空间特征及其强度。这样得到的回归系数空间分布空间特征向量的分布是相似的,但是回归系数可以看出相应的变化的数量大小。

数据性质与预处理

(1)误差 (2)资料的处理。原始场,距平场与标准化场 (3)空间样本点。大范围的空间数据,特别需要注意资料空间代表性。非均匀场与均匀分布场;空间抽样;面积加权。

时空转换

有时空间样本,m远大于时间序列长度n,计算m×mm\times m 矩阵的特征根很困难,可以考虑对其进行时空转换。矩阵A=1nXXA=\frac{1}{n}XX^{'}B=XXB=XX^{'}的特征根不同,但是特征向量是一样的。而可以证明C=XXC=XX^{'}C=XXC^*=X^{'}X有相同的特征根,但特征向量不同。因此,通过时空转换可以求XXX^{'}X矩阵的特征根,进而计算XXXX^{'} 的矩阵的特征根,进而计算XXXX^{'} 矩阵的特征向量。即有:

C×V=V×C^*\times V^*=V^*\times \wedge

VV^*CC^*的特征向量,\wedge 是特征根对角矩阵。根据VV^*是可以求出C的特征向量的,首先计算Va=X×VV_a=X \times V^* ;对VaV_a 进行处理得到C的前n个特征向量VkV_k

Vk=1λkVa(:,k)V_k=\frac{1}{\sqrt{\lambda_k}}V_a(:,k)

得到特征向量V后,就可以计算相应的主成分。

PC=VT×XPC=V^T\times X

前面计算得到的EOF维数是m×mm\times m ,而通过时空转换得到的EOF维数只有m×nm\times n 。即只能得到前n个特征向量。

SVD分析

EOF分析 中一次只分析了一个变量X,地学中常常涉及多个要素场之间的关系。分析多个要素场关系的方法有很多,有混合EOF(CEOF),奇异值分解(SVD)分析,典型相关(CCA)等。

原理: [[SVD]]