气候突变检测

所有变量的变化方式不外乎两种基本形式：一种是连续性变化；另一种是不连续性变化。

不连续变化现象的特点是突发性，所以人们称不连续现象为“突变”。

对突变有不同的理解和定义。

其实，突变可以理解为一种质变，当量变达到一定的限度时发生的质变。

形形色色的突变现象向传统的分析方法提出了挑战。

20世纪60年代末期，法国数学家托姆（Rene Thom）创立了突变理论，很快突变理论风靡一时。

经过十几年从理论到实际应用方面的改进与完善，在科学界产生了很大影响。

随后，突变理论在数学、生物、天文、地震、气象、社会科学等领域得到了广泛应用。

突变理论是以常微分方程为数学基础的（谷松林1993），其精髓是关于奇点的理论。

其要点在于考察某种系统或过程从一种稳定状态到另一种稳定状态的飞跃。

从统计学的角度，可以把突变现象定义为从一个统计特性到另一个统计特性的急剧变化，即从考察统计特征值的变化来定义突变，如考察均值、方差状态的急剧变化。

目前，突变统计分析还相当不成熟。

针对常见的突变问题，人们借助统计检验、最小二乘法、概率论等发展了一些行之有效的检验方法，主要涉及检验均值和方差有无突然漂移、回归系数有无突然改变，以及事件的概率有无突然变化等方面。

这里仅介绍几种在检测气候突变现象中最常用的方法。

顺便指出，突变理论研究中最为活跃同时争议最大的就是有关应用问题。对一些物理机制目前还不甚明确的突变现象，人们很难给予解释。有时使用的检测方法不当，可能会得出错误的结论。因此，建议在确定某气候系统或过程发生突变现象时，最好使用多种方法进行比较。另外，要指定严格的显著性水平进行检验。运用气候学的专业知识对突变现象进行判断也十分重要。

5.1 滑动t检验

5.1.1 方法概述

滑动t检验是通过考察两组样本平均值的差异是否显著来检验突变。其基本思想是把一气候序列中两段子序列均值有无显著差异看做来自两个总体均值有无显著差异的问题来检验。如果两段子序列的均值差异超过了一定的显著性水平，可以认为均值发生了质变，有突变发生。

对于具有n个样本量的时间序列t，人为设置某一时刻为基准点。基准点前后两段子序列和的样本分别为和，两段子序列平均值分别为 和，方差分别为 和。定义统计量：

其中：

 

方程（5.1）遵从自由度的t分布。
 

这一方法的缺点是子序列时段的选择带有人为性。为避免任意选择子序列长度造成突变点的漂移，具体使用这一方法时，可以反复变动子序列长度进行试验比较，以提高计算结果的可靠性。

5.1.2 计算步骤：

(1) 确定基准点前后两子序列的长度，一般取相同长度，即和。

(2) 采取滑动的方法连续设置基准点，依次按方程计算统计量。由于进行滑动的连续计算，可得到统计量序列。

给定显著性水平 ，查 t 分布表（附表2）得到临界值。若 ，则认为基准点前后的两子序列均值无显著差异；否则认为在基准点时刻出现了突变。

在编制程序计算时，滑动计算两子序列平均值和，相当于执行两子序列的滑动平均过程。设子序列长度为，以前个数据之和为基数，依次减去前一个数并加上后一个数求平均，这是第一个序列的滑动平均过程。第二个滑动平均是以第个至个数据之和为基数，再依次减去前一个数并加上后一个数求平均。再用滑动的方式依次计算两子序列各自的方差。

5.1.3 计算结果分析

根据t统计量曲线上的点是否超过t值来判断序列是否出现过突变。如果出现过突变，可以确定大致的时间。另外，根据诊断出的突变点分析突变前后序列的变化趋势。

应用实例[5.1]使用滑动t检验检测1911—1995年中国年平均气温等级序列的突变。这里n=85，两子序列长度n=10。给定显著性水平α=0.01，按照t分布自由度，。这里为了方便编制程序，给定，实际上给出了更严格的显著性水平。将计算出的t统计量序列和 绘成图5.1。从图5.1中看出，在1920年以后，有两处t统计量超过了0.01显著性水平，一处是正值（出现在1920年），另一处是负值（出现在1950年）。说明在近85年来，中国年平均气温出现过两次明显的突变。20世纪20年代经历了一次由冷到暖的转变，50年代经历了一次由增暖转为变冷的明显突变。尽管70年代末80年代初，中国气温与全球气温同步在回升，但没有达到显著性水平。

5.2 克拉默（Cramer）法
5.2.1 方法概述
克拉默（Cramer）方法的原理与t检验类似，区别仅在于它是用比较一个子序列与总序列的平均值的显著差异来检测突变。

设总序列 和子序列 的均值分别为 和 ，总序列方差为，定义统计量：

5.3 山本（Yamamoto）法

山本方法是从气候信息与气候噪声两部分来讨论突变问题的。由于是由山本（R. Yamamoto）最先将信噪比用于确定日本地面气温、降水、日照时数等序列的突变，故称其为山本法（Yamamoto, 1986）。

5.4 曼-肯德尔（Mann-Kendall）法

曼-肯德尔法是一种非参数统计检验方法。在第3章统计检验中介绍的检验方法都是参数方法，即假定了随机变量的分布。非参数检验方法，亦称无分布检验，其优点是不需要样本遵从一定的分布，也不受少数异常值的干扰，更适用于类型变量和顺序变量，计算也比较简便。

由于最初由曼（H.B. Mann）和肯德尔（M.G. Kendall）提出了原理并发展了这一方法，故称其为曼-肯德尔法。然而，当时这一方法仅用于检测序列的变化趋势。后来，经其他人进一步完善和改进，才形成目前的计算格式。

5.5 佩蒂特（Pettitt）方法

佩蒂特（Pettitt）方法是一种与曼-肯德尔法相似的非参数检验方法。由于是由佩蒂特（A.N.Pettitt）最先用于检测突变点的，故将其称为佩蒂特法（Pettitt 1979）。与曼-肯德尔法一样，构造形如方程(5.4)的一秩序列。不同的是，r;是分3种情况定义的，即：

5.6 勒帕热（Le Page）法

勒帕热（Le Page）法是一种无分布双样本的非参数检验方法。它的统计量是由标准的威氏（Wilcoxon）检验和安氏-布氏（Ansari-Bradley）检验之和构成的。由于将两个检验联合在一起的原理最早是由勒帕热（Platt Le Page）提出的，因此将其称为勒帕热检验（Yonetani 1992）。已有研究证明，与其他检验相比，它是一种十分有效的检验方法。但是，迄今为止，它还没有像上述介绍的方法那样广泛地应用到气候研究领域。

5.7 BG(Bernaola-Galvan)分割算法

针对非平稳时间序列的特点，美国波士顿大学的Bernaola-Galvan等人（2001）提出了一种突变检测的分割算法，简称BG分割算法，并将其应用到人的心率序列的突变检测中。随后，封国林等人（2005）通过理想时间序列验证了这一算法的有效性，并利用这一方法检测和分析了北半球树木年轮距平宽度序列不同尺度的突变特征。

Bernaola-Galvan等人（2001）提出的检测方法的主要思想是将非平稳时间序列的突变检测问题视为一个分割问题，即将非平稳时间序列看作由多个具有不同平均值的子序列构成。此方法的目的是找出各子序列之间最大差值的平均值的位置。

6气候序列周期提取方法

近年来，提取时间序列振荡周期的统计方法发展十分迅速。从离散的周期图、方差分析过渡到连续谱分析。然而，周期图不能处理周期的相位突变和周期振幅的变化。方差分析在具体实施时，对原序列寻找一个隐含的显著周期的统计推断是十分巧妙的，但用剩余序列推断第二和第三个周期时，从假设检验意义上讲，就很牵强。就其结果而言，上述两种方法及经典的谐波分析均是从时间域上研究气候序列中周期振荡的方法，它们将气候序列中的周期性视为正弦波，有其固有的局限性，这里不作介绍。

1807年法国数学家傅里叶(J.B.J. Fourier)提出了在有限时间间隔内定义的任何函数均可以用正弦分量的无限谐波的叠加来表示，这样就出现了与时域相对应的频域。特别是1965年出现快速傅里叶变换以来，使频域分析走向实用并迅速拓展。这里将重点介绍以傅里叶变换概念为基础的功率谱、交叉谱以及以自回归模型为基础的最大峡谱。

近年来，又出现了研究周期现象的新技术（奇异谱分析）和时频结构分析的新方法（小波分析），使得提取气候序列周期的技术有了新的飞跃，这些内容将在本章进行介绍。

6. 1功率谱

功率谱分析是以傅里叶变换为基础的频域分析方法。其意义在于将时间序列的总能量分解到不同频率上的分量。根据不同频率波的方差贡献，可以诊断出序列的主要周期，从而确定周期的主要频率，即序列隐含的显著周期。功率谱是一种应用非常广泛的分析周期的方法。有关功率谱的概念、谱分解及傅里叶变换的算法，许多书籍中都有详尽的阐述（黄嘉佑等1984）。这里仅给出有关提取显著周期的具体方法、计算流程及结果分析要点。

6.2 最大熵谱

如上所述，连续功率谱估计需要借助谱窗函数对粗谱进行平滑以求得结果。因此，其统计稳定性和分辨率都与选择的窗函数相关。例如，方程(6.7)就是一种对应于汉宁窗函数的平滑公式。由于使用了与分析系统毫无关系的窗函数，有时可能得出虚假的结论。另外，在连续功率谱估计中，自相关函数估计与样本量大小相关，这也会导致谱估计误差，影响分辨率。可见，功率谱存在分辨率不高和有可能产生虚假频率分量等缺点。

由于功率谱不需要由时间序列本身提供某种参数模式，因此是一种非参数谱估计。Burg (1967) 提出了一种称之为“最大似然”谱估计的方法，从而将谱估计推进到一个新的阶段。最大似然谱的基本思想是以信息论中似然的概念为基础，选择这样一种谱估计，即在外推已知时间序列的自相关函数时，其外推原则是使相应序列在未知点上取值的可能性具有最大的不确定性，也就是不对结果进行人为的主观干预，因而所得信息最多。

最大熵谱估计与确定时间序列的参数模式（如自回归模型）有关，是一种参数谱估计。最大似然谱具有高分辨率等优点，尤其适用于短序列，因此受到广泛关注。

6.3 交叉谱

在实际问题中，我们不仅要研究单个气候序列的频域结构和周期特性，还要分析不同序列在频域变化上的相互关系。因此，需要讨论多个序列（这里仅限两个序列）之间的交叉谱。

6.4 多维最大功率谱

在上一节介绍了使用交叉谱研究两个气候序列之间的凝聚和相位关系。在提取单个气候序列的周期时，最大功率谱表现出相对于普通功率谱的优越性。这一节将最大功率谱推广到多变量形式——多维最大功率谱，它描述的是多个不同气候时间序列之间的交叉能量关系，是一种估计复合谱（海金1986）。多维最大功率谱最早主要应用于地质信号处理中，后来推广到雷达、通讯的信号处理领域，各种计算方法在应用中相继应运而生。在雷达、通讯领域通常将其称为多信道最大功率谱。因此，在方法的描述上有关术语通常使用通讯信号处理的用法。

6.5奇异谱分析

奇异谱分析（Singular Spectrum Analysis, SSA）是一种基于时间序列的动态重构，与经验正交函数（Empirical Orthogonal Function, EOF）相关联的统计技术。它已广泛应用于时间序列信号处理中。SSA的具体操作过程是，将样本量为n的时间序列 \(r(t)\) 按给定的嵌套空间维数m（称为窗口长度）构造成一个资料矩阵。当这个资料矩阵计算出明显成对的特征值，且相应的EOF几乎呈周期性或正交关系时，通常对应着信号中的振荡行为。可以看出，数学上，SSA相当于在延滞坐标上表达的EOF，也可以看作是EOF的一种特殊应用。分解的空间结构与时间尺度有关，可以有效地从一个有限的、含有噪声的时间序列中提取信息。

Broomhead等（1986）最先将SSA引入非线性动力学研究，后来由Vautard（1992）和Ghil等（1991）进行了一系列改进，并应用于研究气候序列的周期振荡现象中。SSA的优点主要表现在两方面：首先，其滤波器不像通常的谱分析需要预先给定，而是根据资料自身最优确定。因此，它适用于确定和寻找噪声系统中的弱信号。特别是，它不需要假定时间序列由不同频率的正弦波叠加而成，也无需将一个本质上是非线性振荡的信号分解为大量正弦波之叠加来讨论。其次，对于嵌套空间维数m的限定，可以使得对振荡的转换进行时间定位。SSA是一种特别适合于识别隐含在气候序列中的弱信号的新型统计技术，也是研究周期振荡现象的有力工具。

6. 6小波分析

小波分析（Wavelet Analysis），又称多分辨率分析（Multiresolution Analysis），是近几年国际上十分热门的一个前沿领域，被认为是傅里叶分析方法的突破性进展。1982年，法国地质学家J. Morlet在分析地震波的局部性质时引入了小波概念（Grossman等，1985）。之后，Grossman等（1985）和Meyer（1990）又对小波进行了一系列深入研究，使小波理论有了坚实的数学基础。进入20世纪90年代，小波分析成为众多学科共同关注的热点。在信号处理、图像处理、地震勘探、数字电路、物理学、应用数学、力学、光学等诸多科技领域得以广泛应用（崔锦秦，1994；秦前清等，1994）。小波分析因其对信号处理具有特殊优势而很快得到气象学家们的重视，并应用于气象和气候序列的时频结构分析中，取得不少引人注目的研究成果（Weng等，1994）。在气候诊断中，广泛使用的傅里叶变换可以显示出气候序列不同尺度的相对贡献，而小波变换不仅可以给出气候序列变化的尺度，还可以显现出变化的时间位置。后者对于气候预测是十分有用的（Arnedo等，1988；Meyer等，1992）。需要指出的是，小波分析是一种基本数学手段，它可以应用在多种领域，可以从统计学角度研究，也可以应用在动力学乃至人工智能中。这里仅介绍用小波分析进行气候序列小波分解的具体方法及主要分析内容。

7气候变量场时空结构的分离
 

某一区域的气候变量场通常由许多个观测站点或网格点构成，这给直接研究其时空变化特征带来困难。如果能用个数较少的几个空间分布模态来描述原变量场，且又能基本涵盖原变量场的信息，是一项很有实用价值的工作。也就是来寻找某种数学表达式，将变量场的主要空间分布结构有效地分离出来。气候统计诊断中应用最为普遍的办法是把原变量场分解为正交函数的线性组合，构成为数很少的互不相关典型模态，代替原始变量场，每个典型模态都含有尽量多的原始场的信息。其中经验正交函数 (Empirical Orthogonal Function, EOF) 分解技术就是这样一种方法。

EOF最早是由 Pearson (1902) 提出来的。20世纪50年代中期，Lorenz将其引入大气科学研究中。由于计算条件的限制，直至20世纪70年代初才在我国的气候研究领域中使用。20世纪70年代中期以后，随着计算机技术的迅速发展，EOF分解技术在气候诊断研究中得以充分应用。之所以被广泛使用，还由于它具有一系列突出的优点：第一，它没有固定的函数，不像有些分解需要以某种特殊函数为基函数，如球谐函数等；第二，它能在有限区域对不规则分布的站点进行分解；第三，它的展开收敛速度快，很容易将变量场的信息集中在几个模态上；第四，分离出的空间结构具有一定的物理意义。正因为如此，EOF已成为气候科学研究中分析变量场特征的主要工具。以EOF为气候特征分析手段的研究成果颇丰，揭示出许多有价值的气候变化事实。

近10年来，气候统计诊断方法有了很大的进展，其中以 EOF 为基础的变量场分解方法的飞跃发展格外引人注目。针对气候变量场特征分析的需要，发展了揭示气象场空间结构和时间相关特征的扩展经验正交函数 (Extended Empirical Orthogonal Function, EEOF)、着重表现空间的相关性分布结构的旋转经验正交函数 (Rotated Empirical Orthogonal Function, REOF)、可以揭示空间行波结构的复经验正交函数 (Complex Empirical Orthogonal Function, CEOF) 和描述动力系统非线性变化特征的主振荡型 (Principal Oscillation Patterns, POPs)。这些方法使气候统计诊断研究开阔了视野，并进入了一个更高水平。本章就以 EOF 为基础，介绍上述几种方法的特点、计算步骤及作者对计算结果分析的一些认识。

当然，EOF的应用范围远不止这一章所包含的内容。第6章叙述的奇异谱分析就是与EOF有联系的统计技术。尤其是近年来，EOF分析方法在应用方面发展十分迅速。利用EOF是正交函数这一基本事实，发展了以EOF为基函数对强迫气候信号进行检测和估计 (North等1995)、对循环稳态型气候时间序列信号进行检测和估计 (Kim 1996) 等技术。张邦林等 (1991) 还提出了基于EOF的气候数值模拟及模式设计的新构思。另外，EOF还被用来作为气候变量缺测资料插补的工具 (孙照渤等1991)。