2021, Vol. 37
2. 中国科学院大气物理研究所LASG，北京 100029;
3. 中国民航大学，天津 300300;
4. 东营市气象局，山东 东营 257091;
5. 国家气象中心，北京 100081
2. State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics, Institute of Atmospheric Physic, Chinese Academy of Sciences, Beijing 100029, China;
3. Civil Aviation University of China, Tianjin 300300, China;
4. Dongying Meteorological Bureau, Shandong 257100, China;
5. National Meteorological Center, Beijing 100081, China
在气象科研和业务工作中, 经常需要分析气象变量变化的原因, 这时可以把该气象变量作为研究对象, 分析该气象变量与其它气象变量之间的同期或前期关系。比如, 在其它气象变量中选择一个变量, 分析其对研究对象的影响, 即分析一个变量与另一个变量之间的相关关系。如果二者相关关系显著, 则表明它们之间关系密切, 可能存在因果关系或相互影响的物理过程, 可以继续深入分析其是否存在物理机理上的关联; 如果二者相关关系不显著, 则表明它们之间不存在密切关系。本文把气象中变量相互关系或关联性分析的方法称为相关分析。相关分析在天气[1-5]、气候[6-10]、气候变化[11-12]、农业气象[13]、气象服务[14]等领域广泛应用。
变量之间的关系从统计学的角度大致可以分为两类: 函数关系和相关关系。函数关系要求非常严格, 一般的数据很难满足函数关系。而相关关系要求相对宽松, 所以被人们广泛接受, 这也是相关分析在气象变量关系研究中广泛应用的原因。
虽然相关分析在气象科研和业务中应用广泛, 但很多人对相关分析依然存在某些困惑, 甚至存在误用现象。因此, 本文对气象相关分析的现有成果进行梳理和总结, 可以为关注气象相关分析的科研和业务工作者提供借鉴和启示。此外, 气象数据正在进入“大数据时代”[9, 14], 因此本文还简要综述了相关分析在“气象大数据”中的应用价值和面临的新挑战。
2 气象相关分析方法的分类气象科研与业务中经常使用的相关有: 点(站点或变量)点相关、点场相关、场场相关。点点相关可以看作两个变量之间的相关, 点场相关可以看作一个变量与场中的每一个变量分别作点点相关, 可以归结为点点相关, 因此也属于两个变量之间的相关。场场相关属于多个变量之间的相关。因此, 气象中的相关大致可以分为两类: 两个变量之间的相关和多个变量之间的相关。下面将分别介绍两个变量之间的相关系数和多个变量之间的相关系数。
2.1 两个变量之间的相关系数 2.1.1 两个定距变量之间的相关系数定距变量是指数值变量, 不同的定距变量之间具有数量上的差别, 可以对定距变量进行加减乘除运算。相关图(或散点图)和相关表可以反映两个数值变量之间的相互关系及其相关变化的方向, 但无法确切地表示它们之间相关的程度。度量定距变量相关性最常用的是皮尔逊相关系数[1-2, 10, 15-16]。
相关系数在科学研究中的应用具有悠久的历史。1886年, 英国科学家高尔顿给出了关于遗传的相关和回归概念, 并且思考了它们在生物遗传研究中的可用性和价值[17]。著名统计学家卡尔·皮尔逊对高尔顿的“相关”概念十分着迷, 经过深入研究, 他在前人的研究基础上, 采用极大似然法, 把一个二元正态分布的相关系数最佳值p用样本积矩相关系数r表示, 这被后人称为“皮尔逊相关系数”或“皮尔逊积矩相关系数”(Pearson product-moment correlation coefficient, 简写为PPMCC或PCCs), 文章中常用r或Pearson's r表示。皮尔逊相关系数是最常用的普通相关系数, 常简称为“相关系数”或“简单相关系数”, 两个一维随机变量x =(x1, x2, ……, xn) 和y =(y1, y2, ……, yn)之间相关系数的计算公式为:
| $ r=\frac{\sum\limits_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)}{\sqrt{\sum\limits_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2}} \sqrt{\sum\limits_{i=1}^{n}\left(y_{i}-\bar{y}\right)^{2}}} $ | (1) |
皮尔逊相关系数r可以表示两个随机变量之间线性关系的强弱, 其取值范围为-1≤r≤1。r越趋近于1, 这表示这两个变量之间正线性相关关系越强; 反之, r越趋近于-1, 这时表示这两个变量之间负线性相关关系越强; 而当r等于0或接近于0时, 表示这两个变量之间不存在线性关系或线性关系很弱。r对于不同的相关现象, 名称有所差异, 一般将反映两变量间直线线性相关关系的统计量称为相关系数(相关系数的平方称为判决系数); 将反映两变量间曲线相关关系的统计量称为非线性相关系数、非线性判决系数; 将反映多个变量之间的多元线性相关关系的统计量称为复相关系数、复判决系数。
皮尔逊相关系数是两个变量之间关系的简单单值度量, 并且其形式适合数学运算, 因此, 应用非常广泛。但需要指出的是, 不能不加辨别地机械计算相关系数, 因为皮尔逊相关系数无法识别非线性关系, 并且皮尔逊相关系数对一个或几个离群(异常)点极为敏感。如图 1中的4个子图, 皮尔逊相关系数全都相同, 但显然这4张图的关系存在较大差异。此外, 通常需要使用t检验对皮尔逊相关系数进行检验, 而t检验是基于数据呈正态分布假设的, 当变量数据不服从正态分布时, 即使对大样本, 皮尔逊相关系数的显著性检验也可能存在较大偏差。
|
图 1 子图a，b，c，d的每个水平变量（x）都有相同的平均值9.0和标准差11.0，每个垂直变量（y）也都有相同的平均值7.5和标准差为4.12，对这4个子图来说，皮尔逊（普通）相关系数是相同的，都为rxy=0.816[2] |
对于点场相关, 虽然常用皮尔逊相关系数进行研究, 但也有很多研究者采用一元线性回归研究点场相关。比如许立言等[18]研究欧亚大陆春季融雪与东亚夏季风之间的关系, 就采用了一元线性回归方法。他们选取春季融雪EOF第2模态的标准化时间序列, 对850 hPa风场、500 hPa高度场和中国夏季降水场分别进行线性回归分析。需要指出的是, 如果用两个原始变量(或距平变量)做一元线性回归, 这时由回归系数得到的回归场是带单位的。而相关场是两个标准化变量之间的相关系数得到的, 是没有单位的。回归场与相关场的两个场分布形势略有差别, 例如高度场的距平场, 往往高纬地区距平变化大, 标准差大, 所以回归场对高纬地区之间的关系反映得更大。而相关场则消除了高低纬之间标准差的不同, 因此两个场有相似的地方, 也有不同的地方。对于变量x和y来说, 回归系数b与相关系数rxy之间的关系为
皮尔逊相关系数在气象科研中广泛应用。比如, Wallace等[19]对500 hPa高度场和海平面气压场, 计算了其中一个格点与场的其它格点之间的皮尔逊相关系数, 发现在北半球冬季存在5个冬季遥相关型; 丁一汇等[20]计算了点场之间的皮尔逊相关系数, 研究了亚洲-太平洋季风区中的遥相关关系, 清楚地揭示了东亚夏季风、印度夏季风和西北太平洋夏季风之间的相互作用。
2.1.2 两个定序变量之间的相关系数定序变量也称为等级变量, 其取值具有等级或次序之分。下面介绍3种常用的定序变量相关系数。
(1) γ系数[21]。
γ系数的公式为:
| $ \gamma=\frac{n_{s}-n_{d}}{n_{s}+n_{d}} $ | (2) |
式中, ns、nd分别为同序对、异序对的数目。同序对、异序对的定义为: 某对样本(x1, y1)和(x2, y2), 如果在等级或次序上x1优于x2, 并且y1优于y2, 则称为同序对; 否则, 称为异序对。
这个公式的直观意义是, 相关系数定义为同序对和异序对数目的差与同序对和异序对数目的和之比值。比如: 若ns = 0, 则γ =-1, 也就是说对于这两个变量, 所有的样本对都呈现出异序, 这时可以认为这两个变量之间完全呈负相关关系; 若nd = 0, 则γ = 1, 也就是对于这两个变量来说, 所有的样本对都呈现出同序, 这时可以认为这两个变量之间完全呈正相关关系。
(2) Spearman相关系数[22]。
针对皮尔逊相关系数无法识别非线性关系, 并且对一个或几个离群(异常)点极为敏感的局限性, Spearman相关系数可以作为皮尔逊相关系数的替代方法。
Spearman相关系数有时也被称为级别(顺序)相关系数或秩相关系数, 该相关系数是根据两个变量的秩(排序后的等级或顺序值)进行相关分析, Spearman相关系数可以用来衡量这两个变量间是否存在单调相关关系。两个一维随机变量x =(x1, x2, ……, xn)和y =(y1, y2, ……, yn)的秩之间的Pearman相关系数定义为:
| $ \rho=\frac{\sum\limits_{i=1}^{n}\left(\mathrm{r}_{i}-\bar{r}\right)\left(s_{i}-\bar{s}\right)}{\sqrt{\sum\limits_{i=1}^{n}\left(\mathrm{r}_{i}-\bar{r}\right)^{2}} \sqrt{\sum\limits_{i=1}^{n}\left(s_{i}-\bar{s}\right)^{2}}} $ | (3) |
其中ri和si分别为xi和yi的秩, i = 1, 2, ……, n。当变量里出现相等值(秩结)的时候, 该值对应的秩为这几个值相对应的秩的平均值。Spearman相关系数可以简化为:
ρ的取值范围为-1 ≤ρ ≤1。当一个变量随另一个变量单调递减时, ρ =-1;反之, 当一个变量随另一个变量单调递增时, ρ = 1。
只要两个变量的值是成对的等级数据, 或者是经由连续变量转化得到的等级数据, 就可以用上述Spearman相关系数的公式进行计算, 分析这两个变量之间的关系。Spearman相关系数与变量的分布和样本容量都没有关系, 并且具有鲁棒性和抗干扰性, 即计算结果对个别异常值不敏感。图 2显示了Spearman相关系数与皮尔逊相关系数之间的关联和差异。图中两个一维变量X和Y之间的皮尔逊相关系数为0.88, 而它们之间的Spearman相关系数为1, 表示它们之间的单调相关程度很强, 等于1。
|
图 2 两个变量X和Y的散点图[23] |
(3) Kendall相关系数[24]。
Kendall相关系数是衡量等级变量相关程度的一个统计量, 其主要思想是根据两个变量间序对的一致性来判断其相关性。
设x、y分别是两个一维随机变量, x =(x1, x2, ……, xn)和y =(y1, y2, ……, yn)。把(xi, yi)记为一个序对, 序对之间的关系为下列三种情形: (1)当xi > xj且yi < yj, 或者xi < xj且yi > yj时, 称这个序对是不一致的; (2)当xi > xj且yi > yj, 或者xi < xj且yi < yj时, 称这个序对是一致的; (3)当xi =xj或者yi =yj时, 则称这个序对既不是一致的也不是不一致的。
那么, 两个一维随机变量X和Y之间的Kendall相关系数τ定义为:
| $ \tau=\frac{2 S}{\frac{1}{2} n(n-1)}-1=\frac{4 S}{n(n-1)}-1 $ | (4) |
式中, S表示一致的序对个数。Kendall相关系数τ的取值范围为-1≤τ≤1。当τ =-1时, 表示这两个随机变量具有完全相反的等级相关性; 当τ = 1时, 表示这两个随机变量具有完全一致的等级相关性; 当τ = 0时, 表示这两个随机变量之间相互独立。具体的气象计算实例可参考文献[1-2]。
2.1.3 两个定类变量之间的相关系数定类变量是指变量的值是研究对象的符号或名称, 每个值代表一个类别, 这些值之间相互平等, 没有次序、大小的区别。
计算定类变量之间的相关关系可以借助列联表。列联表是数据按照两个或更多个属性进行分类后所列出的频数表。假设有两个属性X和Y, 属性X有k类, 用Xi表示第i类, i = 1, 2, ……, k; 属性Y有s类, 用Yj表示第j类, j = 1, 2, ……, s。对于n个样本, 用nij代表既属于特征X的第i类又属于特征Y的第j类的样本频数。这样可以得到一个k行s列的列联表(表 1)。
|
表 1 列联表的一般形式 |
下面, 我们虚构一个简单的2×2列联表, 用来示例说明列联表相关系数的计算。假设我们要研究性别(X)与天气预报评分高低(Y)之间的关系, 在调查的200名预报员中, 预报评分偏高的预报员中男性为70名, 女性为25名, 预报评分偏低的预报员中男性为30名, 女性为75名(表 2)。
|
|
表 2 性别与预报评分高低相关分析的2 × 2列联表 |
在定类变量相关系数的计算方法中, Q系数是最简单的一种方法, Q系数只适用于2×2列联表, 公式如下[25-26]:
| $ Q=\frac{n_{11} n_{22}-n_{12} n_{21}}{n_{11} n_{22}+n_{12} n_{21}} $ | (5) |
为了理解Q系数的构造, 我们取两种极端情况。若n12 =n21 = 0, 则Q系数为1;若n11 =n22 = 0, 则Q系数为-1。显然, 这两种情况都表明性别与预报评分高低呈现出完全相关特征, 而正负号表示相关关系方向的不同。在这个假设的例子中, Q系数为1代表男预报员预报评分高, 女预报员评分低; 而Q系数为-1代表女预报员预报评分高, 男预报员评分低。
除了Q系数之外, 还有λ系数、χ2检验、φ系数、C系数、V系数。由于篇幅关系, 不再一一介绍, 具体细节可参见文献[25-27]。
2.2 多个变量之间的相关系数 2.2.1 偏相关系数在二元或者多元回归分析中, 对于变量之间的相关关系, 可用偏相关系数来表示。在研究多个自变量x1, x2, x3, ……, xk与因变量y之间的线性相关程度时, 如果其它自变量保持不变, 只考虑y与其中某一个xi(i = 1, 2, 3, ……, k)之间的关系, 这种相关叫做偏相关。衡量偏相关程度的指标, 就是偏相关系数。例如在二元线性回归模型中, r01, 2表示x2保持不变时y与x1的偏相关系数, r02, 1表示x1保持不变时y与x2的偏相关系数, r12, 0表示y保持不变时x1与x2的偏相关系数。在偏相关系数中, 还可以根据固定自变量数目的多少, 区分为零阶偏相关系数、一阶偏相关系数、K-1阶偏相关系数等。例如, r0i(i = 1, 2, 3, ……, k)表示零阶偏相关系数(即简单相关系数), r02, 1(x1保持不变)称为一阶偏相关系数, r01, 23(x2和x3保持不变)称为2阶偏相关系数, r01, 234(x2、x3和x4保持不变)称为三阶偏相关系数, 依次类推。偏相关系数的具体计算公式和实例可以参考相关文献[28-29]。
2.2.2 复相关系数衡量一个变量(因变量y)与多个变量(自变量x1, x2, ……, xp)之间线性关系的量称为复相关系数。复相关系数的定义要涉及到多元线性回归和剩余残差的概念, 具体计算公式和实例可以参考相关文献[15]。
2.2.3 典型相关系数典型相关系数是先对所研究的两组变量进行主分量分析, 得到新的线性无关的综合指标, 然后再计算两组综合指标之间的线性相关系数, 进而研究这两组变量之间的相关关系[1-2, 10, 15]。典型相关分析(CCA)可以有效地分离两气象场的最大线性相关模态。
奇异值分解(SVD)的出发点与典型相关相同, 也可以用来分析两个气象场的相关模态, 但计算要简便得多[30-32]。从统计学角度讲, CCA推理更加严瑾, 而SVD需要一定的使用条件[10]。
3 相关分析在气象中应用的新进展相关分析在气象中的应用具有悠久的历史。早在20世纪初, Walker在研究世界各地海平面气压变化之间的关系, 提出全球“三大涛动”的概念时, 就采用了相关分析的研究方法。后来, 在大气遥相关的研究中, 也有不少文献采用了相关分析方法[33-34]。下面介绍最近几年相关分析中出现的一些新进展。
3.1 全窗口滑动相关当计算两个时间序列的相关系数时, 由于相关关系往往并不稳定, 所以当子序列的长度变化时, 相关系数经常会发生变化, 有时甚至会发生很大甚至完全相反的变化。那么如何量化评估时间序列之间相关的稳定性呢?Zhao等[35]提出的一种全窗口相关系数方法, 可以有效地解决这个问题。
下面以他们文献中的计算为例介绍该方法的计算过程。
根据概率统计知识, 求相关系数的序列不能太短, 因此规定滑动窗口最小为11(当然也可以根据所研究的具体问题相应调整)。
时间段为1872—2010年, 因此一共有139年。所以滑动窗口为: 11~129年(确保滑动窗口有足够的自由度和时间段)当滑动窗口为11(12, ……, 129)年时, 所有可能时间段的总数为129(128, ……, 11), 开始年从1872, 1873, ……, 到2000(1999, ……, 1882)。因此, 得到129(128, ……, 11)个相关系数。
这样, 可以得到一个两维的滑动相关场(图 3), 该场由
|
图 3 1872—2010年期间冬季北半球雪盖指数（SCFN）与北极涛动指数（IAO）（a）和冬季SCFN与西伯利亚高压指数（ISH）之间（b）的全窗口相关系数 通过0.05显著性检验的相关系数用填色等值线图表示[35]。 |
最后, 根据滑动窗口的大小, 对每个相关系数做显著性检验后, 计算稳定度百分比(PS), PS定义为显著相关的数目相对于总相关数目N的百分比, 这可以作为度量两个序列之间相关稳定程度的指数。图 3中PS为梯形区域内彩色面积的百分比, 对于图 3a来说, PS=70.6%, 意味着在不同的滑动窗口下大部分相关系数是显著的, 因此可以认为SCFN和AO之间的相关整体看是稳定的。而图 3b中PS=47.3%, 意味着多于一半的滑动相关系数是不显著的, 因此SCFN和SH之间的相关整体看不是十分稳定。而且, 通过全窗口滑动相关能够辨识出信号稳定的时段和相应滑动窗口这二维信息, SCFN和AO之间的相关在1930年之后更稳定, 从年际到80 a尺度的相关都是稳定的; 而SCFN与SH在1930年之后, 相关变得不稳定。
3.2 小波相关小波分析可以对局部的时间或空间频率进行分析, 该方法通过平移伸缩运算对信号或函数进行多尺度细化, 最终达到低频处频率细分, 高频处时间细分, 该方法可聚焦到信号的任意细节, 能自动适应时频信号分析的精细化要求, 因此小波分析被很多专家称为“数学显微镜”。有学者用小波分析研究两个时间序列的关系, Grinsted等[36]仿照Torrence等[37]定义了两个时间序列的小波相关为,
| $ R_{n}^{2}(s)=\frac{\left|\mathrm{M}\left(s^{-1} W_{n}^{X Y}(s)\right)\right|^{2}}{\mathrm{M}\left(s^{-1}\left|W_{n}^{X}(s)\right|^{2}\right) \cdot \mathrm{M}\left(s^{-1}\left|W_{n}^{Y}(s)\right|^{2}\right)} $ | (6) |
式中s是尺度, M是平滑算子, 具体计算公式可参考文献[36]。小波相关的这个定义类似于传统的相关系数, 但分子中有平方, 所以不同于传统的相关系数有负值, 小波相关系数全都为正值。小波相关系数通过位相表示传统相关系数负值的意义。
一些研究应用这种方法检测出了两个因子间时间频率的局部化信息和镜像关系[27]。图 4是AO与BMI(波罗的海最大年海冰面积)小波相干谱。从图 4中可以看出, 1860—1900年, 2~6年的时间尺度上存在大块显著相关区域, 但中间的1865—1882年, 2~4年的时间尺度上相关不显著。在1960年附近, 2~3年的时间尺度存在相关区域。在1900—1950年期间, 显著相关的时间尺度从8~11年变化为3~10年。1915—1990年期间, 显著相关的时间尺度从12~16年变化为5~30年。图 4中所有这些显著相关区域, 平均位相角是174 °, 显示了AO和BMI之间的反位相关系, 并且AO略领先于BMI。
|
图 4 标准化的AO和BMI（波罗的海最大年海冰面积）时间序列之间的小波相关 相对于红噪音的5%显著性用粗等值线显示。所有的显著性区域显示反位相行为。相对的位相关系用箭头显示（同位相指向右，反位相指向左，BMI超前AO 90 o指向正下，BMI滞后AO 90 o指向正上）[36]。 |
气象变量的相关关系存在突变现象, 对相关关系的突变检测具有重要的理论和实用价值。此前国内常用的M-K秩和检验与Yamamoto判据, 一次计算都只能检测序列中某一时间尺度的某次突变, 而以前的小波分析, 虽然可以进行多尺度分析, 但没有进行相关关系突变点的显著性检验。
江剑民[38]提出了可以检验两个序列间多尺度平均值突变现象的相干性(同步或反位相变化)的算法, 并将该算法用于尼罗河年最高与最低水位历史序列, 可以较为客观、自动地检测出两序列在不同时间尺度上突变的相干性。Jiang [39]融合小波分析的多尺度分辨功能和统计学4种参数, 即平均值、方差、趋势以及分段子样本的相关系数, 分别简要概述了多尺度突变点的扫描式检测算法。Zhu等[40]进一步发展和完善了其中的相关系数和方差分析2种算法, 并且利用经过正态化处理后的珠江流域中的西江下游马口水文站月平均流量与西江水域月降水量资料, 给出了具体的应用实例, 主要结果见图 5。
|
图 5 a. 马口站月平均流量正态化指数NSI和西江流域月降水量正态化指数NPI之间相关系数的多尺度突变扫描式U检验结果的等值线；b. NSI（粉色虚线）和NPI（绿色虚线）13点高斯滤波低通曲线 两序列样本中分段子样本相关系数突变点用黑色垂直粗线表示和分时段子样本的相关系数用黑色水平粗线表示[40]。 |
图 5a中横坐标为月份时间(j, 或称为参考点), 横坐标表示尺度(n, 即子样本容量), 等值线表示Zhu等[40]中等式(4)定义的统计量Ur(n, j)。式(4)中Uα(n, j)是判断相应尺度相关关系突变是否显著的临界值, α为显著性水平, 比如0.05、0.01等。当样本容量n > 50时, U0.05 = 1.96, U0.01 =2.58;当样本容量n ≤50时, 可采用学生氏t检验。但还需要对序列的非独立性进行订正[39]。理论上, Ur(n, j) < -1.0的局部最小值中心, 表示相关系数显著变小(弱); 而Ur(n, j) > 1.0的局部最大值中心, 表示相关系数显著变大(强)。由于本例总体数值较低, 如放宽到0.4, 可从图 5a中找出14个突变点。比如, 1965年4月与时间尺度42个月交汇处有一正极大值中心; 1968年9月与时间尺度32个月交汇处有负极小值中心。接下来的极值点还有1971年7月与42个月尺度上的正极大值, 1974年9月与42个月尺度上的负极小值, 1978年4月与37个月尺度上的正极大值, 1986年2月37个月尺度上的负极小值。1989年和1998年之间的三个局部正极大值表示相关关系的接连增强。1999—2010年之间的3个负极小值表示相关关系的接连下降, 但期间包含有一个2001年10月于37个月时间尺度上的小正值中心。整个图面上表现出最强的突变点发生在2012年12月37个月的时间尺度上, 超过了0.1的显著性检验。由图 5b可见, 两条经过平滑滤波后的曲线可以直观地验证粗黑实线所标示的突变点及其相对应的子样本时段的相关系数变化情况。这种算法, 在数学界也获得较好的认可与评价。
4 气象大数据中的相关分析 4.1 气象大数据的特点最近几年, 学术界就大数据在体积、类型、速度和价值这四个方面的特征基本达成了共识, 即大数据的4V特征: Volume(大体积), 数据量为ZB级别; Variety(多样性), 除传统的结构化数据之外, 还有大量的非结构化数据; Velocity(高速度), 秒级的数据分析处理以及输出; Value(价值稀薄), 有效信息占信息总量的比例小。气象数据与气象服务对象的数据也基本具备了大数据的上述特征[41-44], 特别是随着气象及气象服务对象数据的进一步丰富, 未来必将进入气象大数据时代[9, 14], 而基于大数据的人工智能技术也已经在天气预报和气象服务中得到应用[45]。
4.2 相关分析在大数据中的作用现有的知识体系是建立在数据稀缺背景下的, 在大数据时代, 人们的思维和工作方式必须发生变革。相关分析已经成为大数据挖掘与分析的关键应用技术和核心科学问题[46], 在大数据分析中相关关系比因果关系更重要[47], 这指明了相关分析在大数据应用和分析中的重要地位。在大数据研究和机器学习中, 相关分析常被称为关联分析, 就是在关系数据、市场交易数据, 或其它可用的信息载体中, 查找可能存在的关联、相关性或因果结构。关联分析是一种简单实用的分析技术, 可以发现存在于大数据集中的相关性, 从而可以描述事物中某些属性可能出现的规律和模式。关联分析常用的方法有FP-G算法、Apriori算法、FreeSpan算法、Prefixspan算法等。在大数据相关分析的应用领域中, 引人关注的是推荐系统, 即基于相关分析度量出物品相似性、用户相似性等特征, 进而对不同的顾客进行精准的产品推荐[48], 从而可以避免传统广告中“普遍撒网”, 推荐目标不准确的缺陷。对于气象服务网站来说, 基于相关分析的产品推荐也有应用前景。在灾害应急管理中大数据相关分析也有着广泛应用[49]。
Wang等[9]指出, 根据大数据的思想和分析方法, 聚类分析、相似分析、相关分析和机器学习技术在大数据时代的气候预测中将被广泛应用。其实, 在气象大数据时代, 这些技术在天气预报和气象服务中也广泛使用。需要指出的是, 聚类分析、相似分析和机器学习, 都离不开相关分析, 因此相关分析是气象大数据分析中十分重要的技术方法。
4.3 气象中的降维方法与相关分析的关系对于气象大数据来说, 降低数据的维度或者简化数据是一项重要工作, 经验正交函数(EOF)分解或主分量分析(PCA)是一种重要的降维方法[50-51], EOF及其变种方法在气象的多个领域都具有重要应用。EOF分解得到的是空间分型模态, 它的某一分量的空间荷载值就是该分量的时间系数与该格点的时间序列样本之间的相关系数。其实也可以看作为相关场, 与用EOF空间模态绝对值最大的点为定点的点场相关得到的皮尔逊相关系数相关场类似[52], EOF通常也称为分型模态。
4.4 大数据相关分析的新方法对大数据来说, 数据关系往往呈现非线性、高维度等复杂特征, 传统的相关分析方法往往难以有效地探测数据的内在结构与规律, 迫切需要新的相关分析计算方法。下面介绍大数据研究领域出现的两种具有代表性的新算法。
(1) 距离相关。Szekely等[53]从特征函数的距离视角考察了两个随机向量之间的非线性相关系数, 为高维数据的非线性分析提供了有效的度量准则。距离相关具有如下两个优点: ①所度量的相关也包含非线性关系, 而不仅仅局限于线性相关关系; ②可以度量任意两个不同维数的随机向量的相关性。但距离相关系数要进行高维向量间的距离计算和矩阵乘积运算, 所以计算耗时很大, 对计算能力要求较高。气候变化研究[54]中使用了该方法。
(2) MIC相关性。2011年, 《Science》上发表了一篇题为“Detecting novel associations in large data sets”的论文[55], 该论文提出了衡量两个变量之间相关关系的一种新方法——最大信息系数(Maximal Information Coefficient, MIC)。MIC具有两个重要性质: 通用性和均等性。通用性是指, 传统的相关系数往往只能度量某种特定类型的函数类型(如线性、指数、对数或周期性函数)变量之间的相关程度, 而MIC可以度量任意函数形式变量之间的相关性, 因此具有通用性。均等性是指, 对于具有相等MIC值, 但函数形式不同的数据, 外加同等程度的噪音, 然后重新计算MIC, 这时MIC的值仍能保持相等, 而传统的相关系数计算方法很难做到这一点。该方法比皮尔逊相关系数、Spearman相关系数等方法更细致地描述和刻画了两个变量之间的相关关系, 但该方法目前在气象领域的应用还较少。
4.5 大数据相关分析面临的挑战和困难由于大数据具有数据规模大、数据类型复杂、价值密度低等特点, 这为相关分析带来了很多挑战和困难, 表现最突出的就是计算效率问题。杨静等[56]针对传统大数据典型相关分析(CCA)方法的高复杂度在面临大数据PB级数据规模时不再适用的现状, 提出了一种基于云模型的大数据CCA方法, 该方法在云计算架构的基础上, 通过云运算将各端点云合并为中心云, 并据此产生中心云滴, 以中心云滴作为大数据的不确定性复原小样本, 在其上施以CCA运算, 进而提出了具有较高计算效率的大数据CCA分析的云模型。Nguyen等(2014) [57]提出了一种高维相关子空间的搜索方法, 该方法基于相关图的极大团进行挖掘分析, 克服了传统Apriori算法采用的逐层搜索模式计算耗时高的弊端, 为大数据中进行多变量高维相关分析提供了一种快速便捷的计算方法。并行计算也是加快计算能力的有效手段, 以MapReduce为代表的非关系数据型库管理技术为大数据分析与处理提供了一种并行处理架构[46], 可以为大数据相关分析的并行计算提供理论支持。云计算和并行计算领域对相关分析的这些研究, 对于开展气象大数据相关分析研究具有重要的借鉴意义。
5 气象相关分析应用中存在的问题气象领域在应用相关分析的时候, 有几个问题是需要注意的。
(1) 需选择有物理意义的气象变量做相关分析。在气象研究中做相关分析时, 不能滥用。只有具有物理意义的变量之间做出的相关, 才能更加稳定, 在预报预测中也更有价值。
(2) 当计算相关的两个变量都包含明显的趋势变化成分时, 原变量之间的相关特征可能被歪曲(夸大或者缩小)。施能等[58]的数值试验结果表明: 两个变量带有相同性质的趋势时, 这两个变量之间的相关系数会增加, 具体表现为正相关系数值被夸大, 负相关系数值被减小; 而当这两个变量带有性质相反的趋势变化时, 则会使这两个变量之间的相关系数减小, 具体表现为正相关的数值被减小, 而负相关的数值被夸大。
(3) 需要做严格的显著性检验。根据概率统计的术语, 我们得到的气象资料应称为样本。无论气象资料的时间有多长, 都是有限的, 而气象变量的总体是无限的。当根据有限的样本推断无限总体的性质时, 必须进行显著性检验。具体检验方法可参考相关文献[59-60]。相关系数的检验不应该称为“信度检验”, 其检验水平也不应该称为“95%置信度”, 而应该称为“显著性检验”, 显著性水平为5%(10%, 1%)[60]。如果把点场相关中相关系数高于或低于某一数值绘图为阴影区, 那么这些区域应该称为“超过或高于显著性水平5%(10%, 1%)的区域”, 相关系数的检验过程应该称为“显著性检验”。显著性水平α需要取小值, 气象中常取5%(10%, 1%)。现在的很多统计软件(比如SAS、SPSS、Matlab、R等)中, 计算相关系数时, 常自动给出P值, 当P值小于5%(10%, 1%)时, 即可以说这两个序列的相关超过显著性水平5%(10%, 1%)。
当多个统计检验结果必须被同时评估的时候出现了特殊问题, 这被称为“检验的多重性问题”[2]。比如求点场相关时, 在相关场中逐点进行相关系数检验完成后, 如果在若干个格点上表现是显著的, 那么是否可以认为这个相关场是显著的呢?这个问题已经由Taleb[61]根据所谓的“无穷猴子理论”进行了有趣的说明。如果我们能以某种方式, 把无穷数量的猴子放在键盘前, 并且允许它们随机的打字, 事实上, 肯定有一只猴子最终能打出Iliad。但是据此推断, 这只猴子与众不同是不合理的。例如, 推断这个猴子接下来比其它的猴子有更高的概率打出Odyssey。假定无限数量的猴子打字, 一只猴子复制出可辨识内容的事实, 并没有提供反对原假设的充分证据, 即这只是一只普通的猴子, 其将来的文字输出, 与其它任何猴子没有什么不同。Livezey等[62]用概率论的观点把含N个格点的相关场的检验看成N次掷硬币试验, 即每个格点只有两种检验结果: 通过和不通过检验的两个互斥事件, 实际上这是一个二项分布检验问题。设检验的显著水平为0.05, 那么对某一格点来说, 成功通过检验的概率为p = 5%, 不成功的概率为q = 95%。根据二项分布, 即可以计算出N次试验中, M次事件成功发生(该问题中则为总格点数N个中有M个格点通过检验)的概率。例如, 当场的总格点数为N=30时, 4个点通过检验的概率为0.045, 5个点通过检验的概率为0.016, ……。这样, 可以计算“至少4个点通过检验”事件的概率为0.045 + 0.016 +……≈0.062。从而可以确定出相当于0.05水平下的临界点数。该例的计算表明, 在显著性水平0.05下, 该相关场至少需要超过4.24个点通过检验时才能认为该相关场是显著的, 即其显著区域面积为总场格点数的14.1%(4.24/30)时, 该场是显著的。类似可以计算, 当总格点数为N=80时, 显著临界区域面积为10%, N=500时, 显著临界区域面积为7%。
两个气象变量的相关系数是否显著, 可以采用t检验, 但当气象变量本身具有强的持续性或高的自相关时, t检验的自由度不能用n-2, 而应该用有效自由度, 有效自由度的计算方法[63-64]可以参考文献[27]。当然, 这时也可以采用蒙特卡洛检验[17, 48]。当对时间序列做了滑动平均后再求相关时, 这时的有效自由度就发生了变化, 就可以采用蒙特卡洛检验, 具体的例子可参考宋燕等[65]的例子。
6 相关分析在气象中应用的未来发展趋势天气预报和气候预测正在进入大数据和智能预报时代, 这对相关分析提出了新要求。大数据的核心是预测, 大数据之所以能够预测未来, 是基于对相关关系的准确把握。传统的统计方法对样本数据的正态性、变量的独立性、变量个数、假设检验等都有较高要求。大数据技术是对传统数量统计学方法的拓展和延伸, 大数据分析侧重于高维建模、复杂网络建模、非参数模型等技术方法从种类繁多、数量庞大的数据中快速获取有价值的信息[57, 66-67]。大数据时代气象科研和业务工作对相关分析提出了新要求, 根据我们的理解, 可能主要表现在以下方面。
(1) 对于不符合正态分布的变量如何做相关, 并且检验其相关的显著性。
(2) 如何对变量之间不同尺度的信号求相关, 以及把这种多尺度的相关应用在天气预报和气候预测中。
(3) 如何求非线性相关?可能是一个值得深入研究的问题。因为传统的皮尔逊相关系数所求得的相关关系只是线性关系。而很多气象变量之间存在复杂的非线性关系, 如何客观定量地描述这种非线性关系, 需要进一步深入研究。
(4) 如何解决相关分析的高效计算问题?未来气象大数据也将具有数据规模大、数据类型复杂、价值密度低等特点, 如何在平衡计算能力和业务时效之间提高计算效率, 是需要深入研究的问题。可能需要在并行计算和云计算方面对相关分析算法做深入研究和改进。
本文系统综述了不同计算形式的相关分析在气象中的应用, 特别是最近几年相关分析的新进展, 尽管仍然会挂一漏万。希望本文能够为关注气象数据分析理论与应用的专家提供借鉴相关分析作为探寻与发现气象变量内在规律的重要工具, 在气象大数据分析与挖掘中具有重要应用, 但在大数据时代面临新的挑战, 尚存在不少问题值得深入研究。
| [1] |
WILKS D S. Statistical Methods in the Atmospheric Sciences(3nd edn)[M]. Amsterdam: Academic Press, 2011.
|
| [2] |
朱玉祥, 李宏毅, 吕行, 译著. 大气科学中的统计方法(第三版)[M]. 北京: 气象出版社, 2017.
|
| [3] |
韦青, 李伟, 彭颂, 等. 国家级天气预报检验分析系统建设与应用[J]. 应用气象学报, 2019, 30(2): 245-256. |
| [4] |
张自银, 赵秀娟, 熊亚军, 等. 基于动态统计预报方法的京津冀雾霾中期预报试验[J]. 应用气象学报, 2018, 29(1): 57-69. |
| [5] |
胡波, 俞燎霓, 滕代高. 高斯过程回归方法在浙江沿海海盗冬季阵风预报中的应用试验[J]. 热带气象学报, 2019, 35(6): 767-779. |
| [6] |
HAM Y G, KIM J H, LUO J J. Deep learning for multi-year ENSO Forecasts[J]. Nature, 2019, 573: 568-572. DOI:10.1038/s41586-019-1559-7 |
| [7] |
王欢, 李栋梁. 气候变暖背景下全球海温对中国东部夏季降水年代际转折的影响[J]. 热带气象学报, 2019, 35(3): 398-408. |
| [8] |
VON STORCH H, ZWIERS F W. Statistical analysis in climate research[M]. Cambridge: Cambridge University Press, 1999.
|
| [9] |
WANG X, SONG L C, WANG G F, et al. Operational climate prediction in the era of big data in China: Reviews and prospects[J]. J Meteor Res, 2016, 30(3): 444-456. DOI:10.1007/s13351-016-6081-3 |
| [10] |
魏凤英. 现代气候统计诊断与预测技术[M]. 北京: 气象出版社, 2007.
|
| [11] |
孙颖, 尹红, 田沁花, 等. 全球和中国区域近50年气候变化检测归因研究进展[J]. 气候变化研究进展, 2013, 9(4): 235-245. DOI:10.3969/j.issn.1673-1719.2013.04.001 |
| [12] |
SUN Y, ZHANG X B, ZWIERS F W, et al. Rapid increase in the risk of extreme summer heat in eastern China[J]. Nature Climate Change, 2014, 4: 1082-1085. DOI:10.1038/nclimate2410 |
| [13] |
吴门新, 庄立伟, 侯英雨, 等. 中国农业气象业务系统(CAgMSS)设计与实现[J]. 应用气象学报, 2019, 30(5): 513-527. |
| [14] |
杨承睿, 任芳, 马楠. 试论大数据在气象服务中的应用[J]. 农业网络信息, 2016(8): 53-55. DOI:10.3969/j.issn.1672-6251.2016.08.013 |
| [15] |
黄嘉佑. 气象统计分析与预报方法(第四版)[M]. 北京: 气象出版杜, 2016.
|
| [16] |
丁一汇, 李霄, 李巧萍. 气候变暖背景下中国地面风速变化研究进展[J]. 应用气象学报, 2020, 31(1): 1-12. |
| [17] |
GALTO N F. Regressio n towards mediocrity in hereditary stature[J]. J Anthropol Inst, 1886, 15: 246-263. |
| [18] |
许立言, 武炳义. 欧亚大陆春季融雪量与东亚夏季风的可能联系[J]. 大气科学, 2012, 36(6): 1 180-1 190. |
| [19] |
WALLACE J M, GUTZLER D S. Teleconnections in the geopotential height field during the northern hemisphere winter[J]. Mon Wea Rev, 1981, 109: 784-812. DOI:10.1175/1520-0493(1981)109<0784:TITGHF>2.0.CO;2 |
| [20] |
丁一汇, 刘芸芸. 亚洲-太平洋季风区的遥相关研究[J]. 气象学报, 2008, 66(5): 670-682. DOI:10.3321/j.issn:0577-6619.2008.05.002 |
| [21] |
GOODMAM L A, KRUSKAL W H. Measure of association for cross classifications, Ⅳ: Simplification of asymptotic variances[J]. Journal of the American Statistical Association, 1972, 67(338): 415-421. DOI:10.1080/01621459.1972.10482401 |
| [22] |
SPEARMAN C. The proof and measurement of association between two things[J]. The American Journal of Psychology, 1904, 15(1): 72-101. DOI:10.2307/1412159 |
| [23] |
樊嵘, 孟大志, 徐大舜. 统计相关分析研究进展[J]. 数学建模及其应用, 2014, 3(1): 1-12. DOI:10.3969/j.issn.2095-3070.2014.01.002 |
| [24] |
KENDALL M G. A new measure of rank correlation[J]. Biometrika, 1938, 30(1/2): 81-39. DOI:10.2307/2332226 |
| [25] |
GOODMAM L A, KRUSKAL W H. Measure of association for cross classifications[J]. Journal of the American Statistical Association, 1954, 49(268): 732-764. |
| [26] |
GOODMAM L A, KRUSKAL W H. Measure of association for cross classifications, Ⅱ: Further discussion and references[J]. Journal of the American Statistical Association, 1959, 54(285): 123-163. DOI:10.1080/01621459.1959.10501503 |
| [27] |
黄嘉佑, 李庆祥. 气象数据统计分析方法[M]. 北京: 气象出版社, 2015.
|
| [28] |
施能. 气象科研与预报中的多元分析方法[M]. 北京: 气象出版社, 2002.
|
| [29] |
施能. 气象统计预报[M]. 北京: 气象出版社, 2009.
|
| [30] |
BRETHERTON C S, SMITH C, WALLACE J M. An intercomparison of methods for finding coupled patterns in climate data[J]. J Climate, 1992, 5: 541-560. DOI:10.1175/1520-0442(1992)005<0541:AIOMFF>2.0.CO;2 |
| [31] |
WALLACE J M, SMITH C, BRETHERTON C S. Singular value decomposition of sea-surface temperature and 500 mb heights anormalies[J]. J Climate, 1992, 5: 561-576. DOI:10.1175/1520-0442(1992)005<0561:SVDOWS>2.0.CO;2 |
| [32] |
史加荣, 杨柳. 基于奇异值分解的气象数据推测[J]. 气象学报, 2020, 78(1): 128-142. |
| [33] |
DING Q H, WANG B. Circumglobal teleconnection in the Northern Hemisphere summer[J]. J Climate, 2005, 18: 3 483-3 505. DOI:10.1175/JCLI3473.1 |
| [34] |
DING Q H, WANG B. Intraseasonal teleconnection between the summer Eurasian wave train and the Indian monsoon[J]. J Climate, 2007, 20: 3 751-3 767. DOI:10.1175/JCLI4221.1 |
| [35] |
ZHAO L, ZHU Y X, LIU H W, et al. A stable snow-atmosphere coupled mode[J]. Climate Dyn, 2016, 47: 2 085-2 104. DOI:10.1007/s00382-015-2952-z |
| [36] |
Grinsted A, Moore J, Jevrejeva S. Application of the cross wavelet transform and wavelet coherence to geophysical time series[J]. Nonlin Process Geophys, 2004, 11: 561-566. DOI:10.5194/npg-11-561-2004 |
| [37] |
TORRENCE C, COMPO G P. A practical guide to wavelet analysis[J]. Bull Am Meteorol Soc, 1998, 79: 61-78. DOI:10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2 |
| [38] |
江剑民. 多尺度突变现象的扫描式t检验方法及其相干性分析[J]. 地球物理学报, 2001, 44(1): 31-39. DOI:10.3321/j.issn:0001-5733.2001.01.004 |
| [39] |
JIANG J M. Scanning detections of multi-scale significant change-points in subseries means, variances, trends and correlations[C]. Sixth International Conference on Fuzzy Systems and Knowledge Discovery, 2009, 285: 609-613.
|
| [40] |
ZHU Y X, JIANG J M, HUANG C X, et al. Applications of multiscale change-point detections to monthly streamflow and rainfall in Xijiang River in Southern China, Part I: correlation and variance[J]. Theor and Appl Climatal, 2019, 136: 237-248. DOI:10.1007/s00704-018-2480-y |
| [41] |
张巍. 大数据在气象服务中的研究与应用[J]. 计算机光盘软件与应用, 2014(18): 27-28. |
| [42] |
赵蓓. 大数据时代对气象服务的推动[J]. 考试周刊, 2014(33): 195-196. DOI:10.3969/j.issn.1673-8918.2014.33.175 |
| [43] |
FAGHMOUS J, KUMAR V. A big data guide to understanding climate change: The case for theory-guided data science[J]. Big Data, 2014, 2: 155-163. DOI:10.1089/big.2014.0026 |
| [44] |
崔巍. 大数据在气象服务中的应用与分析[J]. 低碳技术, 2016(9): 121-122. |
| [45] |
嵇磊, 王在文, 陈敏, 等. 人工智能技术能否提高地面气温预报的精度——记AI Challenger 2018全球天气预报挑战赛[J]. 气象学报, 2019, 77(5): 960-964. |
| [46] |
梁吉业, 冯晨娇, 宋鹏. 大数据相关分析综述[J]. 计算机学报, 2016, 39(1): 1-18. |
| [47] |
盛阳燕, 周涛译. 大数据时代[M]. 杭州: 浙江人民出版社, 2013.
|
| [48] |
LU L Y, MEDO M, YEUNG C H, et al. Recommender systems[J]. Physics Reports, 2012, 1-49. |
| [49] |
LU X, BENGTSSON L, HOLME P. Predictability of population displacement after the 2010 Haiti earthquake[C]. Proceedings of the National Academy of Sciences of the United States of America, 2012, 109(29): 11 576-11 581.
|
| [50] |
HANNACHI A, JOLLIFFE I T, STEPHENSON D B. Empirical orthogonal functions and related techniques in atmospheric science: A review[J]. Int J Climatol, 2007, 27(9): 1 119-1 152. DOI:10.1002/joc.1499 |
| [51] |
李颖, 陈怀亮. 机器学习技术在现代农业气象中的应用[J]. 应用气象学报, 2020, 31(3): 257-266. |
| [52] |
吴洪宝, 吴蕾. 气候变率诊断和预测方法[M]. 北京: 气象出版社, 2005.
|
| [53] |
SZEKELY G J, RIZZO M LK, BAKIROV N K. Measuring and testing dependence by correlation of distances[J]. The Annals of Statistics, 2007, 35(6): 2 769-2 794. |
| [54] |
RACHERLA P N, SHINDELL D T, FALUVEGI G S. The added value to global model projections of climate change by dynamical downscaling: A case study over the continental U. S. using the GISS-ModelE2 and WRF models[J]. J Geophy Res, 2012, 1117(D20118): 1-8. |
| [55] |
RESHEF D N, RESHEF Y A, FINUCANE H K, et al. Detecting novel associations in large data sets[J]. Science, 2011, 334: 1 518-1 524. DOI:10.1126/science.1205438 |
| [56] |
杨静, 李文平, 张建沛. 大数据典型相关分析的云模型方法[J]. 通信学报, 2013, 34(10): 121-134. DOI:10.3969/j.issn.1000-436x.2013.10.015 |
| [57] |
NGUYEN H V, MULLER E, BOM K. A near-linear time subspace search scheme for unsupervised selection of correlation features[J]. Big Data Research, 2014, 1: 37-51. DOI:10.1016/j.bdr.2014.07.004 |
| [58] |
施能, 魏凤英, 封国林. 气象场相关分析及合成分析中的梦特卡洛检验[J]. 南京气象学院学报, 1997, 20(3): 355-359. |
| [59] |
ZWIERS F W, VON STORCH H. On the role of statistics in climate research[J]. Int J Climatol, 2004, 24: 665-680. DOI:10.1002/joc.1027 |
| [60] |
黄嘉佑. 气象中使用统计检验的几个问题[J]. 气象, 2005, 31(7): 3-5. DOI:10.3969/j.issn.1000-0526.2005.07.001 |
| [61] |
TALEB N N. Fooled by Randomness[M]. Texere: New York, 2001.
|
| [62] |
LIVEZEY R E, CHEN W Y. Statistical field significance and its determination by Monte Carlo techniques[J]. Mon Wea Rev, 1983, 111(1): 46-59. DOI:10.1175/1520-0493(1983)111<0046:SFSAID>2.0.CO;2 |
| [63] |
BAYLEY G V, HAMMERSLEY J M. The effective number of independent observations in an auto-correlated time series[J]. J Roy Statist Soc, 1946, 8(1B): 184-197. |
| [64] |
AMBAUM M H. Significance tests in climate science[J]. J Climate, 2010, 23(22): 5 927-5 932. DOI:10.1175/2010JCLI3746.1 |
| [65] |
宋燕, 李智才, 肖子牛, 等. 太阳活动与高原积雪及东亚环流的年代际相关分析[J]. 高原气象, 2016, 35(5): 1 135-1 147. |
| [66] |
郭华东, 王力哲, 陈方, 等. 科学大数据与数字地球[J]. 科学通报, 2014, 59: 1 047-1 054. |
| [67] |
唐亘. 精通数据科学: 从线性回归到深度学习[M]. 北京: 人民邮电出版社, 2018.
|