气候统计气候序列的周期提取方法课件

气候序列的周期提取方法气候序列的周期提取方法气候序列的周期提取方法气候序列的周期提取方法谱分析谱分析谱分析谱分析时间序列分析方法时间序列分析方法存在两种基本的时间序列分析方法：存在两种基本的时间序列分析方法：时间域分析方法；时间域分析方法；离散数据离散数据Markov Chains连续数据连续数据自回归过程自回归过程频率域分析方法；频率域分析方法；时间序列分析方法类似于理论分布，即用几个时间序列分析方法类似于理论分布，即用几个参数作为数据的代表，但理论分布并不考虑数参数作为数据的代表，但理论分布并不考虑数据的排序特征，而这里的时间序列方法是对数据的排序特征，而这里的时间序列方法是对数据的排序特征进行推断，从而也可用于对未来据的排序特征进行推断，从而也可用于对未来数据特征的推断，这要求数据应满足平稳性。数据特征的推断，这要求数据应满足平稳性。平稳时间序列平稳时间序列将某种随机变量按出现时间的将某种随机变量按出现时间的顺序排列顺序排列起来称为起来称为时间序列时间序列.平稳时间序列是指其平稳时间序列是指其中随机变量的时间序列中随机变量的时间序列,它的前期演变过它的前期演变过程的统计相关规律在未来的一段时间内程的统计相关规律在未来的一段时间内是不变的是不变的,也就是说它的也就是说它的数学期望数学期望值与方值与方差是不变，即变量的分布特征不随时间差是不变，即变量的分布特征不随时间变化。变化。时间序列的平稳性时间序列的平稳性（stationarity）相对而言，如果一个时间序列的均值和相对而言，如果一个时间序列的均值和方差在统计意义上变化很弱，则可以视方差在统计意义上变化很弱，则可以视其为平稳时间序列；其为平稳时间序列；但通常大气变量的时间变化序列是不平但通常大气变量的时间变化序列是不平稳的稳的如温度，存在年循环，对于月际尺度变化，如温度，存在年循环，对于月际尺度变化，不平稳（不平稳（1月月7月月）；）；风速，存在日变化，对于尺度为几个小时的风速，存在日变化，对于尺度为几个小时的分析是非平稳的。分析是非平稳的。时间序列的平稳化处理时间序列的平稳化处理必要性必要性对于绝大多数时间序列的分析方法而言，对于绝大多数时间序列的分析方法而言，其前提条件是时间序列的平稳性；其前提条件是时间序列的平稳性；因此必须将非平稳序列进行平稳性转化。因此必须将非平稳序列进行平稳性转化。时间序列的平稳化处理时间序列的平稳化处理方法方法1数学变换：数学变换：剔除年循环，即体现为数据在不同月具有不剔除年循环，即体现为数据在不同月具有不同的均值和方差，处理为月标准化数据；同的均值和方差，处理为月标准化数据；将数据标准化，即使数据具有定常的均值和将数据标准化，即使数据具有定常的均值和方差（均值为方差（均值为0，方差为，方差为1）；）；时间序列的平稳化处理时间序列的平稳化处理方法方法2数据分级：数据分级：对于分析足够短的数据而言，能够视其为近对于分析足够短的数据而言，能够视其为近平稳序列；平稳序列；如我们可以假定某站点如我们可以假定某站点1月气温来自于相同月气温来自于相同的物理过程，则可以用长度为的物理过程，则可以用长度为31天的气温分天的气温分析其日变化，但这种假设不能扩展至析其日变化，但这种假设不能扩展至7月，月，甚至是甚至是2月。月。频率域分析方法频率域分析方法频率域分析方法根据不同时间尺度（或频率域分析方法根据不同时间尺度（或频率）的贡献来表证数据；频率）的贡献来表证数据；每一时间尺度可由一对每一时间尺度可由一对sine和和cosine函数函数表示；表示；整个的时间序列就是由不同尺度的整个的时间序列就是由不同尺度的sine和和cosine函数的叠加构成；函数的叠加构成；通常我们对单个尺度的波更加的感兴趣。通常我们对单个尺度的波更加的感兴趣。频率域分析方法频率域分析方法因此频率域分析方法涉及到将包含因此频率域分析方法涉及到将包含n个点个点的原始数据转化为一系列的周期函数；的原始数据转化为一系列的周期函数；虽然，直观上频率域分析方法较难以让虽然，直观上频率域分析方法较难以让人接受；人接受；但是，这种方法在大气科学分析中是非但是，这种方法在大气科学分析中是非常常用的，也是非常重要的方法，常能常常用的，也是非常重要的方法，常能为我们提供原数据重要的信息。为我们提供原数据重要的信息。谐波分析谐波分析（Harmonic analysis）谐波分析（谐波分析（Harmonic analysis）谐波分析是将一系列谐波分析是将一系列sine和和cosine函数叠函数叠加在一起来表征原始数据的振荡或波动；加在一起来表征原始数据的振荡或波动；Cosine和和sine函数特点函数特点一个谐波函数表征一个简单的一个谐波函数表征一个简单的时间序列时间序列存在的问题存在的问题即便是时间序列具有很好的正弦曲线的即便是时间序列具有很好的正弦曲线的特征，我们用正弦曲线表征该数据时，特征，我们用正弦曲线表征该数据时，仍然存在以下问题：仍然存在以下问题：1.数据是时间的函数，而三角函数却是角度数据是时间的函数，而三角函数却是角度的函数；的函数；2.余弦和正弦曲线波动的范围在余弦和正弦曲线波动的范围在+1和和-1之间，之间，而数据的振荡范围通常不能满足此限制；而数据的振荡范围通常不能满足此限制；3.余弦曲线极值位于余弦曲线极值位于 处，而处，而此位置对于正弦曲线则为平均值。此位置对于正弦曲线则为平均值。一个谐波函数表征一个简单的一个谐波函数表征一个简单的时间序列时间序列解决方法解决方法1.将数据记录长度将数据记录长度n看作为一个周期或基看作为一个周期或基本周期，则有：本周期，则有：为基本频率，下标为基本频率，下标1表示整个数表示整个数据只有一个完整的循环。据只有一个完整的循环。一个谐波函数表征一个简单的一个谐波函数表征一个简单的时间序列时间序列解决方法解决方法2.将将cosine或或sine函数向上或向下移动到函数向上或向下移动到原数据的基本水平处，然后拉伸或压原数据的基本水平处，然后拉伸或压缩到与原数据相同的振幅范围。缩到与原数据相同的振幅范围。：振幅，则振幅的最大和最小值为振幅，则振幅的最大和最小值为一个谐波函数表征一个简单的一个谐波函数表征一个简单的时间序列时间序列解决方法解决方法3.应将谐波函数水平的移动，从而与原应将谐波函数水平的移动，从而与原序列的脊和槽相匹配。序列的脊和槽相匹配。为位相角，即将为位相角，即将cosine函数向右移动函数向右移动角度角度 ，则新的函数在，则新的函数在 处达处达到最大值。到最大值。举例举例11943-1989年年Ithaca平均每月温度；平均每月温度；t=1表示表示1月，月，t=2表示表示2月，月，；整个序列的年平均温度为整个序列的年平均温度为46.1；原数据近似正弦曲线；原数据近似正弦曲线；最暖的最暖的7月平均温度为月平均温度为68.8，最冷的，最冷的1月平均温度为月平均温度为22.2 单个谐波的振幅和位相估计单个谐波的振幅和位相估计由三角函数的特性：由三角函数的特性：知：知：谐波与多元线性回归谐波与多元线性回归当对上式进行变量代换可转化为一般的当对上式进行变量代换可转化为一般的多元线性回归方程：多元线性回归方程：由此，可利用最小二乘法估计参数，且：由此，可利用最小二乘法估计参数，且：位相角位相角计算公式：计算公式：且满足且满足参数参数在满足时间步长相等，且无缺测的前提在满足时间步长相等，且无缺测的前提下，通过最小二乘估计得到参数：下，通过最小二乘估计得到参数：举例举例2同上例，有下表同上例，有下表-86.417-110.3290.0000.000552.9Sums:0.00027.4000.0001.00027.412-19.65034.034-0.5000.86639.311-42.86724.750-0.8660.50049.510-60.2090.000-1.0000.00060.29-58.109-33.550-0.866-0.50067.18-34.400-59.581-0.500-0.86668.870.000-64.3000.000-1.00064.3627.400-47.4570.500-0.86654.8538.450-22.2000.866+0.50044.4432.2000.0001.0000.00032.2319.65811.3500.8660.50022.7211.10019.2250.5000.86622.21t举例举例2则可计算得到平均值为则可计算得到平均值为46.1 与前述例子的结果非常接近，但位相角与前述例子的结果非常接近，但位相角向右移动向右移动8度（大约度（大约1周的时间），则结周的时间），则结果更接近实际数据。果更接近实际数据。高阶谐波高阶谐波由于年循环近似于正弦曲线，因此单个由于年循环近似于正弦曲线，因此单个谐波便可很好的拟合；谐波便可很好的拟合；但这并不意味着单个谐波可以很好的表但这并不意味着单个谐波可以很好的表征任何时间序列；征任何时间序列；类似于多元回归问题，增加更多的余弦类似于多元回归问题，增加更多的余弦波也会提高谐波分析的拟合结果，但也波也会提高谐波分析的拟合结果，但也同样面临可能会过度拟合的问题。同样面临可能会过度拟合的问题。高阶谐波高阶谐波已证明包含已证明包含n个点的原始数据通过叠加个点的原始数据通过叠加n/2个谐波函数便可以完全表征原数据，个谐波函数便可以完全表征原数据，即存在一个可以通过所有点的谐波。即存在一个可以通过所有点的谐波。频率（圆频率）为：频率（圆频率）为：解释解释K=2，二阶谐波，其振幅和位相分别为，二阶谐波，其振幅和位相分别为C2和和 ；表示时间表示时间t从从0到到n/2有一个循环，而从有一个循环，而从n/2到到n还有一个完整的循环，即还有一个完整的循环，即2波；波；同样，同样，k=3表示整个时间过程存在表示整个时间过程存在3个波。个波。参数参数在时间步长均等，且无缺测点的前提下，在时间步长均等，且无缺测点的前提下，通过求解多元回归方程可得到：通过求解多元回归方程可得到：离散离散Fourier变换变换Ak和和Bk为为Fourier系数系数过度拟合过度拟合有关多元回归分析表明，当拟合线通过有关多元回归分析表明，当拟合线通过所有数据点时，复相关系数为所有数据点时，复相关系数为100%，为，为过度拟合。过度拟合。同样，当谐波方程中包含同样，当谐波方程中包含n/2个谐波时，个谐波时，也为过度拟合。也为过度拟合。过度拟合过度拟合由于每个谐波项都包含由于每个谐波项都包含2个参数，即振幅个参数，即振幅和位相，同时方程中还包含截距，即样和位相，同时方程中还包含截距，即样本平均值：本平均值：当当n为偶数时，可以有为偶数时，可以有n/2个谐波，个谐波，n/2个振幅个振幅+(n/2-1)个位相个位相+样本平均值，方样本平均值，方程包含有程包含有n个参数；个参数；当当n为奇数时，可以有为奇数时，可以有(n-1)/2个谐波，参数个谐波，参数为为(n-1)/2个振幅个振幅+(n-1)/2个位相个位相+样本平均样本平均值，则方程仍包含有值，则方程仍包含有n个参数。个参数。谐波个数的选取谐波个数的选取通常并不需要将所有的通常并不需要将所有的n/2个谐波均用于个谐波均用于拟合原数据，而是仅用几个谐波便可以拟合原数据，而是仅用几个谐波便可以很好的表征原数据的变化特征；很好的表征原数据的变化特征；仅当分析的目的是准备表征原数据时，仅当分析的目的是准备表征原数据时，选用选用n/2个谐波进行分析。个谐波进行分析。举例举例3原数据为原数据为El PASO地区地区1948-1983年连续年连续5天无降水的相对频率的年变化；天无降水的相对频率的年变化；举例举例3由原数据可见夏季较为湿润，而春、秋由原数据可见夏季较为湿润，而春、秋两季较为干燥两季较为干燥非对称的年循环；非对称的年循环；同时叠加有不规则的、短期波动；同时叠加有不规则的、短期波动；采用前采用前3个谐波来拟合原数据，得到平均个谐波来拟合原数据，得到平均值以及前两个谐波参数为：值以及前两个谐波参数为：举例举例3图图谱分析谱分析（Spectral analysis）谐波方程与多元回归方程差异谐波方程与多元回归方程差异谐波方程与多元回归方程不同在于每个谐波方程与多元回归方程不同在于每个谐波之间是彼此独立的，因此每个谐波谐波之间是彼此独立的，因此每个谐波参数可以独立的计算出来，当有新的谐参数可以独立的计算出来，当有新的谐波加入方程时，则已在谐波方程中的谐波加入方程时，则已在谐波方程中的谐波参数保持不变。波参数保持不变。谐波函数彼此独立谐波函数彼此独立谐波函数彼此独立的特性来自于谐波函数彼此独立的特性来自于sine和和cosine函数彼此正交性；函数彼此正交性；对于任一整数对于任一整数k和和j有：有：谐波函数彼此独立谐波函数彼此独立例如考虑例如考虑2个预报因子：个预报因子：它们的相关为：它们的相关为：单个谐波方差贡献单个谐波方差贡献由于每个谐波之间彼此独立，因此它们由于每个谐波之间彼此独立，因此它们对方程的方差贡献并不随谐波方程的变对方程的方差贡献并不随谐波方程的变化而变化，如对第化而变化，如对第k个谐波而言，其方差个谐波而言，其方差贡献为：贡献为：为原序列的样本方差为原序列的样本方差谐波总方差贡献谐波总方差贡献谐波方程中所有谐波的方差贡献为：谐波方程中所有谐波的方差贡献为：若方程中有若方程中有n/2个可能的谐波，则：个可能的谐波，则：离散功率谱离散功率谱（周期图（周期图/Fourier线谱）线谱）概述概述从谱分析的观点看，一个时间序列的谱从谱分析的观点看，一个时间序列的谱可以显示出不同频率的振荡对该序列变可以显示出不同频率的振荡对该序列变化的贡献；化的贡献；Panofsky和和Brier(1958)类似的解释为：类似的解释为：一个光谱显示出不同波长或频率的光对一个光谱显示出不同波长或频率的光对给定光源能量的贡献。给定光源能量的贡献。功率谱图功率谱图由于振幅是频率的函数，因此最简单的由于振幅是频率的函数，因此最简单的功率谱图由振幅的平方构成功率谱图由振幅的平方构成离散功离散功率谱；率谱；也可以为标准谱密度，即也可以为标准谱密度，即由于功率谱图虽然提供了不同频率谐波由于功率谱图虽然提供了不同频率谐波对原数据的贡献如何，但并没有提供位对原数据的贡献如何，但并没有提供位相角的信息，即没有提供不同频率谐波相角的信息，即没有提供不同频率谐波随时间的变化信息，从而无法重构时间随时间的变化信息，从而无法重构时间序列。序列。功率谱图功率谱图通常功率谱图的纵坐标为对数坐标，这通常功率谱图的纵坐标为对数坐标，这样有利于更好的绘制出对主要由少数几样有利于更好的绘制出对主要由少数几个频率的谐波主导的时间序列的谱图；个频率的谐波主导的时间序列的谱图；否则对能量贡献较小的其它谐波将很难否则对能量贡献较小的其它谐波将很难直接从图中看出其贡献量。直接从图中看出其贡献量。基本频率和基本频率和Nyquist频率频率根据频率定义：根据频率定义：当当k=1时，基本频率，时，基本频率，当当k=n/2时，时，Nyquist频率，频率，功率谱图横坐标功率谱图横坐标功率谱图横坐标可以是频率功率谱图横坐标可以是频率也可以是频率也可以是频率 则频率可以从基本频则频率可以从基本频1/n到到Nyquist频率频率1/2可以是周期可以是周期还可以是波数还可以是波数k举例举例4离散离散Fourier变换变换Month19871988kTk(months)AkBkCk121.420.6124-0.140.440.46217.922.5212-23.76-2.2023.86335.932.938-0.990.391.06447.743.646-0.46-1.251.33556.456.554.8-0.02-0.430.43666.361.964-1.49-2.152.62770.971.673.43-0.53-0.070.53865.869.983-0.34-0.210.40960.157.992.671.560.071.561045.445.2102.40.130.220.261139.540.5112.180.520.110.531231.326.71220.790.79举例举例41987和和1988年年Ithaca地区观测的月平均地区观测的月平均温度；温度；在未作谱分析前，我们已知该数据的基在未作谱分析前，我们已知该数据的基本特征为存在年循环；本特征为存在年循环；即周期近似为即周期近似为12个月的循环。个月的循环。离散离散Fourier变化结果变化结果举例举例4结果分析结果分析n=24，可以由，可以由12个不同频率的谐波构成；个不同频率的谐波构成；贡献最大是周期为贡献最大是周期为12个月的谐波，即年个月的谐波，即年循环；循环；表中表中B12=0，23个个Fourier系数，外加样系数，外加样本平均值；本平均值；举例举例5分析分析1951-1979年年SOI指数功率谱；指数功率谱；原序列存在准周期变化；原序列存在准周期变化；举例举例5离散功率谱图离散功率谱图(线谱线谱/离散谱离散谱)举例举例5结果分析结果分析谱图显示出谱图显示出SOI指数序列存在典型的指数序列存在典型的3-7年的振荡；年的振荡；谱图并未画出所有频率的贡献；谱图并未画出所有频率的贡献；较低频变化反映出我们感兴趣的物理现较低频变化反映出我们感兴趣的物理现象，即象，即ENSO循环。循环。位相谱位相谱在离散功率谱估计中，也可以画出位相在离散功率谱估计中，也可以画出位相随频率随频率/波数波数/周期的变化，称之为位相谱周期的变化，称之为位相谱()位相谱能为我们提供不同频率分量的初位相谱能为我们提供不同频率分量的初始信息，因此包含了信号的大量信息；始信息，因此包含了信号的大量信息；用相同的振幅谱和不同的位相谱重构的用相同的振幅谱和不同的位相谱重构的信号是完全不同的。信号是完全不同的。离散功率谱计算离散功率谱计算由由 计算计算Fourier系数，进而得系数，进而得 到振幅；到振幅；但上述方式仅适用于但上述方式仅适用于 ；通常对通常对k=n/2用下式计算：用下式计算：离散功率谱计算离散功率谱计算通常为了使起始时刻位于通常为了使起始时刻位于0时刻，通常离时刻，通常离散功率谱计算的不同波数散功率谱计算的不同波数k的的Fourier系系数估计为：数估计为：且用原序列或距平序列计算结果是一致且用原序列或距平序列计算结果是一致的。的。离散功率谱计算离散功率谱计算实际上，采用上述方式计算离散实际上，采用上述方式计算离散Fourier变换是非常低效的；变换是非常低效的；原因是存在多次反复调用完全相同的计原因是存在多次反复调用完全相同的计算过程，如上例中：算过程，如上例中：t=4，k=1，t=4，k=2，谱计算谱计算快速快速Fourier变换变换(Fast Fourier Transforms,FFTs)采用采用FFTs可以解决上述问题；可以解决上述问题；很多计算软件或程序中包含有一个或多很多计算软件或程序中包含有一个或多个个FFT程序；程序；对于较长序列，对于较长序列，FFT相对于回归方法计相对于回归方法计算速度提高很多，对于算速度提高很多，对于n=100，FFT的速的速度是后者的度是后者的100倍，而倍，而n=10000，则提高，则提高750倍。倍。FFTs通常通常FFTs采用的是欧拉复指数计算方程：采用的是欧拉复指数计算方程：这样处理的方式主要是基于表达方式简这样处理的方式主要是基于表达方式简便，且易于处理，而其数学含义与前述便，且易于处理，而其数学含义与前述形式是完全一致的。形式是完全一致的。复指数形式下的谐波方程复指数形式下的谐波方程则原谐波方程变为：则原谐波方程变为：是复是复Fourier系数系数Aliasing(假名假名)问题问题功率谱分析周期存在的一个问题就是功率谱分析周期存在的一个问题就是“假假名名”问题；问题；在功率谱分析中，可能序列存在的重要的在功率谱分析中，可能序列存在的重要的物理过程频率比物理过程频率比Nyquist频率更高频，但频率更高频，但这些短周期振荡无法由直接分辨出来，则这些短周期振荡无法由直接分辨出来，则它们的作用会体现在较长周期（它们的作用会体现在较长周期（和和 之间）中；之间）中；即这些高频变化被即这些高频变化被“冒名顶替冒名顶替”了，就是了，就是“假名假名”现象。现象。Aliasing形成原因形成原因f=4/5f=1/5Aliasing形成原因形成原因上图说明实际的时间序列的物理过程包上图说明实际的时间序列的物理过程包含比含比Nyquist频率高的变化部分；频率高的变化部分；当对数据的取样少，即数据间隔过大而当对数据的取样少，即数据间隔过大而不能体现这种快变化；不能体现这种快变化；而这种比而这种比Nyquist频率高的频率并不会因频率高的频率并不会因此从数据真实变化中消失；此从数据真实变化中消失；同时这些高频的作用便虚假的体现在可同时这些高频的作用便虚假的体现在可分辨的频率中。分辨的频率中。Aliasing形成原因形成原因在时间序列中可能真实存在比在时间序列中可能真实存在比Nyquist频频率高的变化部分；率高的变化部分；由于这部分高频信号无法由谱分析得到，由于这部分高频信号无法由谱分析得到，从而这些高频信号的能量会体现在频率从而这些高频信号的能量会体现在频率介于介于0到到Nyquist频率之间的信号中；频率之间的信号中；也就是说也就是说0到到Nyquist频率顶替了其它高频率顶替了其它高频信号的作用。频信号的作用。Aliasing形成原因形成原因造成假名现象的根本原因是数据离散的造成假名现象的根本原因是数据离散的观测结果，最短的周期大小直接由数据观测结果，最短的周期大小直接由数据取样间隔决定，因此，很可能在取样间取样间隔决定，因此，很可能在取样间隔中存在一定的周期变化，而无法在得隔中存在一定的周期变化，而无法在得到的数据中体现出来。到的数据中体现出来。Aliasing如何体现如何体现无法分辨的高频变化能量会错误的添加无法分辨的高频变化能量会错误的添加到可分辨的到可分辨的n/2个可分辨的频率变化上；个可分辨的频率变化上；对于无法分辨的高频频率对于无法分辨的高频频率 ，可体，可体现在可分辨的频率部分现在可分辨的频率部分 ，有：，有：Aliasing如何体现如何体现上式的含义是，高频率的会以上式的含义是，高频率的会以Nyqusit频频率的整数倍作为每次折叠的率的整数倍作为每次折叠的“折叠折叠”点，点，对称的叠加到可分辨的频率部分，因此对称的叠加到可分辨的频率部分，因此Nyqusit频率又称为：频率又称为：“folding”frequency；如：如：意味着，稍高于意味着，稍高于 的频率会叠加在接近低频的部分，而稍的频率会叠加在接近低频的部分，而稍低于低于 的频率会叠加在接近的频率会叠加在接近Nyquist频频率部分；率部分；Aliasing现象举例现象举例如何避免如何避免Aliasing一但我们选定数据，则没有办法去除一但我们选定数据，则没有办法去除“Aliasing”现象；现象；当然，我们可以通过提高数据分辨率的方法尽当然，我们可以通过提高数据分辨率的方法尽可能的避免该现象；可能的避免该现象；或者根据已知的物理过程，来确定资料或者根据已知的物理过程，来确定资料/样本的样本的变化率，从而达到去除假名现象；变化率，从而达到去除假名现象；但对于探索性的研究，即不知道物理本质的问但对于探索性的研究，即不知道物理本质的问题，是没有办法去除假名现象的，但我们期望题，是没有办法去除假名现象的，但我们期望接近接近Nyquist频率的功率谱能量接近频率的功率谱能量接近0，则可能，则可能说明高频部分的能量很小，而上图不符合该假说明高频部分的能量很小，而上图不符合该假设设。连续功率谱估计连续功率谱估计上述有关谱计算的方法是基于离散功率上述有关谱计算的方法是基于离散功率谱得到的；谱得到的；另一种计算方法为连续功率谱估计方法；另一种计算方法为连续功率谱估计方法；这种方法基于时间序列的自相关过程作这种方法基于时间序列的自相关过程作间接的估计。间接的估计。连续功率谱估计连续功率谱估计对于以任一周期对于以任一周期T变化的时间函数变化的时间函数 ，可以展开为傅立叶级数：可以展开为傅立叶级数：连续功率谱估计连续功率谱估计上式参数为：上式参数为：连续功率谱估计连续功率谱估计通常在气候序列分析中，我们希望了解通常在气候序列分析中，我们希望了解在无限时间序列中各种连续频率的波动在无限时间序列中各种连续频率的波动结构，因此有必要把离散谱推广到连续结构，因此有必要把离散谱推广到连续谱上面，即把周期谱上面，即把周期T推广至无限大；推广至无限大；连续功率谱估计连续功率谱估计将参数带入傅氏级数公式，并利用三角将参数带入傅氏级数公式，并利用三角函数和差化积公式有：函数和差化积公式有：连续功率谱估计连续功率谱估计Fourier思想：一个非周期信号可以看成思想：一个非周期信号可以看成是周期无限长的周期信号，当周期增加是周期无限长的周期信号，当周期增加时，基本频越小，则成谐波关系的各分时，基本频越小，则成谐波关系的各分量在频率上就会越来越接近，当周期变量在频率上就会越来越接近，当周期变的无穷大时，离散的线谱就形成为一个的无穷大时，离散的线谱就形成为一个连续谱，即从求和变成了积分。连续谱，即从求和变成了积分。即：任意两个相邻谐波之间的频率差趋即：任意两个相邻谐波之间的频率差趋近于近于0连续功率谱估计连续功率谱估计离散功率谱只给出了频谱在离散点（离散功率谱只给出了频谱在离散点（）上的值，而无法反映这些点之间的频）上的值，而无法反映这些点之间的频谱内容；谱内容；把离散功率谱绘制成连续的谱线，并不把离散功率谱绘制成连续的谱线，并不表示它就是连续功率谱；表示它就是连续功率谱；因此必须采用连续功率谱计算这些点的因此必须采用连续功率谱计算这些点的谱值。谱值。连续功率谱估计连续功率谱估计在把周期在把周期T推广至无穷大过程中，任意推广至无穷大过程中，任意两个相邻波数对应的谐波间的频率差无两个相邻波数对应的谐波间的频率差无限接近于限接近于0，则有：，则有：连续功率谱估计连续功率谱估计则进一步令：则进一步令：并设：并设：连续功率谱估计连续功率谱估计连续谱中的位相谱为连续谱中的位相谱为 ：并记振幅谱为：并记振幅谱为：则有：则有：连续功率谱估计连续功率谱估计利用欧拉公式，上式变为利用欧拉公式，上式变为则有：则有：令令则有：则有：连续功率谱估计连续功率谱估计 称为时间函数称为时间函数 的谱，可的谱，可以表示为：以表示为：带入带入 的表达式，则有：的表达式，则有：连续功率谱估计连续功率谱估计上述时间函数与它的复谱的关系本来上述时间函数与它的复谱的关系本来可以作为时间序列分析的工具，但由可以作为时间序列分析的工具，但由于复谱是复数，使用不方便；于复谱是复数，使用不方便；而且，我们感兴趣的仅是某个周期波而且，我们感兴趣的仅是某个周期波动所起的主要作用，即仅对它们的方动所起的主要作用，即仅对它们的方差贡献感兴趣，因此在时间序列分析差贡献感兴趣，因此在时间序列分析中常用的是功率谱，也可以称为方差中常用的是功率谱，也可以称为方差谱或能谱；谱或能谱；连续功率谱估计连续功率谱估计由前述分析可知，功率谱分析的是不同周期波由前述分析可知，功率谱分析的是不同周期波动的方差贡献；动的方差贡献；来源于物理学中功率的概念；来源于物理学中功率的概念；若电阻为若电阻为1个单位，瞬时电压为个单位，瞬时电压为 ，则它的，则它的瞬时功率为瞬时功率为 ，它的总能量为（体现了数，它的总能量为（体现了数学期望为学期望为0的的 的方差）：的方差）：功率谱就是将总能量或总方差分解为各波动的功率谱就是将总能量或总方差分解为各波动的部分能量或方差之和。部分能量或方差之和。连续功率谱估计连续功率谱估计对于数学期望为零的时间序列对于数学期望为零的时间序列 ，其连，其连续功率谱可表示为：续功率谱可表示为：连续功率谱估计连续功率谱估计对于式子：对于式子：其中其中 为共轭复谱为共轭复谱 功率谱密度：功率谱密度：连续功率谱估计连续功率谱估计连续功率谱估计并不是由上式计算得到，连续功率谱估计并不是由上式计算得到，而是通过时间函数的自相关函数间接估而是通过时间函数的自相关函数间接估计得到；计得到；对于标准化时间函数对于标准化时间函数 ，其数学期望，其数学期望为为0，方差为，方差为1，则其自相关函数可表示，则其自相关函数可表示为：为：连续功率谱估计连续功率谱估计对于公式：对于公式：其中其中 为落后时间步长为落后时间步长上式可表述为：上式可表述为：由由Fourier反变换得：反变换得：谱计算谱计算连续功率谱估计连续功率谱估计第一步：计算样本自相关函数；第一步：计算样本自相关函数；m为最大落后时间长度为最大落后时间长度第二步：求粗谱估计；第二步：求粗谱估计；利用数值积分中的梯形法得到：利用数值积分中的梯形法得到：谱计算谱计算连续功率谱估计连续功率谱估计由于数据本身是离散点，则同样满足频由于数据本身是离散点，则同样满足频率从率从0 到到0.5的变化；的变化；若取样点数为若取样点数为m，则波数，则波数k最大可取到最大可取到m/2；若令若令l=2k，则，则l取值为从取值为从0到到m之间的任之间的任一整数，有粗谱估计：一整数，有粗谱估计：谱计算谱计算连续功率谱估计连续功率谱估计对上式作平滑处理，以消除粗谱估计的对上式作平滑处理，以消除粗谱估计的小波动；小波动；采用二项系数平滑得到：采用二项系数平滑得到：谱计算谱计算连续功率谱估计连续功率谱估计上述平滑效果相当于原粗谱公式乘以一上述平滑效果相当于原粗谱公式乘以一个窗函数个窗函数则有最终功率谱计算式：则有最终功率谱计算式：绘图绘图以波数以波数l为横坐标，平滑功率谱密度估计为横坐标，平滑功率谱密度估计值值Sl为纵坐标绘图；为纵坐标绘图；横坐标也可以为频率和周期：横坐标也可以为频率和周期：注意事项注意事项连续功率谱估计的方法并不只有上述方法，自连续功率谱估计的方法并不只有上述方法，自相关函数和平滑方案都可以有一些差别；相关函数和平滑方案都可以有一些差别；计算过程中人为的选择最大落后长度计算过程中人为的选择最大落后长度m，会改，会改变谱估计曲线：变谱估计曲线：m取值小，所反映的波数少，谱的分辨率低；取值小，所反映的波数少，谱的分辨率低；m取值大，虽然估计的谱分量较多，但取值大，虽然估计的谱分量较多，但m太大时，太大时，参与自相关系数计算的样本容量小（参与自相关系数计算的样本容量小（n-m）估计的）估计的自相关函数差；自相关函数差；因此，通常因此，通常m的取值范围在序列长度的取值范围在序列长度n的的1/10-1/3之间。之间。举例举例6离散功率谱估计离散功率谱估计举例举例7连续功率谱估计连续功率谱估计离散功率谱检验离散功率谱检验原假设：原假设：有检验统计量：有检验统计量：遵从分子自由度为遵从分子自由度为2，分母自由度为，分母自由度为 n-2-1的的F分布。分布。举例举例6检验检验连续功率谱检验连续功率谱检验通常，在连续功率谱图上，功率谱图的通常，在连续功率谱图上，功率谱图的峰点（极大值）所对应的周期可定为主峰点（极大值）所对应的周期可定为主要周期；要周期；但这些周期是否显著，要通过显著性检但这些周期是否显著，要通过显著性检验方可确定；验方可确定；显著性检验通常是与非周期性随机过程显著性检验通常是与非周期性随机过程作比较，即白噪音和红噪音过程。作比较，即白噪音和红噪音过程。红噪音过程的功率谱红噪音过程的功率谱在红噪音过程中每一时刻的观测值仅与在红噪音过程中每一时刻的观测值仅与前一时刻的观测值有关，即一阶自回归前一时刻的观测值有关，即一阶自回归过程；过程；红噪音过程功率谱谱密度为：红噪音过程功率谱谱密度为：白噪音过程的功率谱白噪音过程的功率谱白噪音过程是一无周期性的纯随机过程；白噪音过程是一无周期性的纯随机过程；其自相关函数为：其自相关函数为：其功率谱密度：其功率谱密度：功率谱的显著性检验功率谱的显著性检验过程过程1原假设：所要检验的气象要素的总体谱原假设：所要检验的气象要素的总体谱为某一非周期的随机过程的谱；为某一非周期的随机过程的谱；建立检验统计量，某一频率上的样本谱建立检验统计量，某一频率上的样本谱估计值与假设过程的平均谱估计值之比估计值与假设过程的平均谱估计值之比遵从被其自由度去除的遵从被其自由度去除的 分布：分布：功率谱的显著性检验功率谱的显著性检验过程过程2在上式中，自由度的计算和样本容量在上式中，自由度的计算和样本容量n及及所取的最大落后长度所取的最大落后长度m有关，即：有关，即：在显著性水平为在显著性水平为0.05下，拒绝原假设的下，拒绝原假设的要求是：要求是：功率谱的显著性检验功率谱的显著性检验过程过程3判断用那一种过程来检验，有判据：判断用那一种过程来检验，有判据：在显著性水平为在显著性水平为0.05下，有临界值：下，有临界值：若实际计算的落后一个时刻的自相关系数有若实际计算的落后一个时刻的自相关系数有关系：关系：，即序列存在明显的持续性，即序列存在明显的持续性，则用红噪声过程检验，否则用白噪声过程检则用红噪声过程检验，否则用白噪声过程检验。验。功率谱的显著性检验功率谱的显著性检验过程过程4对于红噪声过程检验有平均红噪声谱为：对于红噪声过程检验有平均红噪声谱为：其中：其中：为估计的样本平均谱值。为估计的样本平均谱值。功率谱的显著性检验功率谱的显著性检验过程过程5对于白噪声过程检验有平均白噪声谱为：对于白噪声过程检验有平均白噪声谱为：举例举例7检验检验存在存在5年和年和2.5年年2个显著的周期个显著的周期最大熵谱最大熵谱(Maximum-Entropy Spectral)概述概述连续功率谱估计缺点连续功率谱估计缺点在连续功率谱的谱图中，功率谱值随频率变化在连续功率谱的谱图中，功率谱值随频率变化是比较平滑的；是比较平滑的；如果我们分析时间序列的主要目的是寻找序列如果我们分析时间序列的主要目的是寻找序列的主要周期，则总希望所寻找的主要周期所对的主要周期，则总希望所寻找的主要周期所对应的功率谱值比较突出，比较大，从而容易分应的功率谱值比较突出，比较大，从而容易分辨出较精确的周期；辨出较精确的周期；此外，气象要素的时间序列，通常受观测资料此外，气象要素的时间序列，通常受观测资料的限制，容量较小，而在连续功率谱估计中，的限制，容量较小，而在连续功率谱估计中，自相关函数估计与样本容量有关，常会带来真自相关函数估计与样本容量有关，常会带来真正谱估计的误差。正谱估计的误差。概述概述最大熵谱最大熵谱自自1980年以来，最大熵谱分析方法逐渐年以来，最大熵谱分析方法逐渐被人们所认识，并得到广泛应用；被人们所认识，并得到广泛应用；该方法的特点是高分辨率和短时性，即该方法的特点是高分辨率和短时性，即只需较少的资料就可以得到较高的分辨只需较少的资料就可以得到较高的分辨率；率；且该方法提取的主次周期更符合实际。且该方法提取的主次周期更符合实际。概述概述熵熵“熵熵”（entropy）是德国物理学家）是德国物理学家Rudolf Clausius在在1850年创造的一个术年创造的一个术语，语，他用这个词来表示任何一种能量在他用这个词来表示任何一种能量在空间中分布的均匀程度，能量分布的越空间中分布的均匀程度，能量分布的越均匀，熵就越大，如果我们所考虑的系均匀，熵就越大，如果我们所考虑的系统，能量是完全均匀分布的，则这个系统，能量是完全均匀分布的，则这个系统的熵就达到最大值。统的熵就达到最大值。概述概述热力学中的熵热力学中的熵物理学上熵是指热能除以温度所得的商，标志物理学上熵是指热能除以温度所得的商，标志热能转化为功的程度；热能转化为功的程度；在所分析的特定系统中如果能量密度分配不均，在所分析的特定系统中如果能量密度分配不均，则能量倾向于从密度较高的地方流向密度较低则能量倾向于从密度较高的地方流向密度较低的地方，直到一切达到均匀为止；的地方，直到一切达到均匀为止；熵原文的字意是转变，随着转换的进行，系统熵原文的字意是转变，随着转换的进行，系统趋于平衡态，熵值也越来越大；趋于平衡态，熵值也越来越大；从物质微观热运动看，熵是混乱程度的标志，从物质微观热运动看，熵是混乱程度的标志，系统越无序、越混乱、熵就越大。系统越无序、越混乱、熵就越大。概述概述信息论中的熵信息论中的熵1在信息论中，熵是信息的度量单位；在信息论中，熵是信息的度量单位；信息论的创始人信息论的创始人Shannon在其著作的在其著作的通信的通信的数学理论中数学理论中提出了建立在概率统计模型上的提出了建立在概率统计模型上的信息度量，他把信息定义为信息度量，他把信息定义为“用来消除不确定用来消除不确定性的东西性的东西”；Shannon的公式为：的公式为：l(A)=-logP(A)l(A)为度量事件为度量事件A发生所提供的信息量，发生所提供的信息量，P(A)为事为事件件A发生的概率；发生的概率；和和H=-sum(P*log(P)称为熵。称为熵。在大多数情况下，我们不了解随机事件的概率在大多数情况下，我们不了解随机事件的概率分布，可能仅知道它的平均值。分布，可能仅知道它的平均值。概述概述信息论中的熵信息论中的熵2例如，我们知道班上同学的考试成绩由三个分例如，我们知道班上同学的考试成绩由三个分数档：数档：80，90和和100分，且已知平均成绩为分，且已知平均成绩为90分，则有：分，则有：且有约束条件：且有约束条件：对于有无限多组解的方程组，如何选出最合理对于有无限多组解的方程组，如何选出最合理的解，即如何从可能存在所有分布中选出的解，即如何从可能存在所有分布中选出“最最合理合理”的分布，这种挑选的标准就是最大信息的分布，这种挑选的标准就是最大信息熵原理。熵原理。概述概述信息论中的熵信息论中的熵3最大熵原理的问题可表述为：设有某一最大熵原理的问题可表述为：设有某一离散随机变量离散随机变量X，其概率分布，其概率分布p(x)未知，未知，已知其与若干函数的数学期望：已知其与若干函数的数学期望：求最佳估计求最佳估计由最大熵原理求解该问题，可表述为具由最大熵原理求解该问题，可表述为具有约束条件的优化问题。有约束条件的优化问题。概述概述信息论中的熵信息论中的熵4则具有约束条件的优化问题表述为：则具有约束条件的优化问题表述为：目标函数：目标函数：约束条件为：约束条件为：求解求解熵与最大熵谱熵与最大熵谱对于连续型随机变量对于连续型随机变量x，其熵可写为：，其熵可写为：为为x的概率密度函数，若的概率密度函数，若x遵从数学期遵从数学期望为望为0，方差为，方差为 的的Gaussian分布，则：分布，则：熵与最大熵谱熵与最大熵谱则可导出：则可导出：对于连续型变量，呈正态分布时，其熵对于连续型变量，呈正态分布时，其熵值最大，将熵与功率谱联系得到：值最大，将熵与功率谱联系得到：详细描述见详细描述见Berg(1975)博士论文中博士论文中p16-18页。页。熵与最大熵谱熵与最大熵谱可见将熵与功率谱联系，是求解功率谱可见将熵与功率谱联系，是求解功率谱的另一种方法，当熵值最大时，则功率的另一种方法，当熵值最大时，则功率谱也最大，这样的功率谱就是最大熵谱。谱也最大，这样的功率谱就是最大熵谱。接下来的问题就是如何求解最大熵？接下来的问题就是如何求解最大熵？求解最大熵求解最大熵Berg的最大熵理论的最大熵理论Berg在他在他1975年的博士论文中详细讲述年的博士论文中详细讲述过有关他求解最大熵的方法；过有关他求解最大熵的方法；即给定约束条件：即给定约束条件：对应于对应于p阶自回归过程：阶自回归过程：求解最大熵求解最大熵Berg的最大熵理论的最大熵理论上式中上式中 满足均值为满足均值为0，方差为，方差为则问题化为如何求解则问题化为如何求解可由可由Yule-Walker方程求解上述参数；方程求解上述参数；则最后得到最大熵谱为：则最后得到最大熵谱为：最大熵与连续功率谱比较最大熵与连续功率谱比较在连续功率谱估计中平滑过程会受到临在连续功率谱估计中平滑过程会受到临近频率振动方差的影响；近频率振动方差的影响；最大熵谱是在最小二乘意义下对资料最最大熵谱是在最小二乘意义下对资料最佳拟合的自回归过程的理论谱；佳拟合的自回归过程的理论谱；该方法常用来分析不易分辨的周期。该方法常用来分析不易分辨的周期。最大熵谱的使用最大熵谱的使用对于平稳的时间序列是一种非常有效的对于平稳的时间序列是一种非常有效的提取周期的方法；提取周期的方法；但对于非平稳序列则可能会失效；但对于非平稳序列则可能会失效；MEM最大的问题在于自回归阶数最大的问题在于自回归阶数p的选的选择；择；对于给定的时间序列，随着对于给定的时间序列，随着p的增加，峰的增加，峰值也相应增加；值也相应增加；自回归过程阶数选择问题自回归过程阶数选择问题因此，存在一个矛盾，即提高分辨率因此，存在一个矛盾，即提高分辨率（较高的（较高的p）与较少的虚假峰值（较低的）与较少的虚假峰值（较低的p）；）；SSA-MTM软件提供了默认值以供选择；软件提供了默认值以供选择；p不应超过时间序列的一半；不应超过时间序列的一半；自回归过程阶数选择问题自回归过程阶数选择问题为了克服阶数选择的矛盾，可采取下列为了克服阶数选择的矛盾，可采取下列方式：方式：通过减小通过减小p值来确定有效峰值；值来确定有效峰值；与其它功率谱方法计算结果进行比较；与其它功率谱方法计算结果进行比较；首先使用首先使用SSA将原数据分解为几个分量，每将原数据分解为几个分量，每个分量仅包含几个谐波，因此可以选择较小个分量仅包含几个谐波，因此可以选择较小的的p值分析。值分析。SSA-MTM软件介绍软件介绍软件全称：软件全称：Singular Spectrum Analysis-Multitaper Method；SSA-MTM软件是由软件是由SSA-MTM工作组研工作组研发，研发人员主要来自于发，研发人员主要来自于UCLA(University of California Los Angeles),包包括括Myles Allen,Mike Dettinger,Kayo Ide,Dmitri Kondrashov,Michael Ghil,Mike Mann等。等。SSA-MTM软件介绍软件介绍软件是免费的，可以自由下载；软件是免费的，可以自由下载；但软件有版权，使用该软件应在所撰写但软件有版权，使用该软件应在所撰写的论文或书籍中注明；的论文或书籍中注明；如果认为该软件有助于自己的科研工作，如果认为该软件有助于自己的科研工作，可以通过引用工作组成员的文章回报他可以通过引用工作组成员的文章回报他们的工作们的工作如如Ghil M.,Allen M.R.,Mann M.E.,Vautard R.等撰写的文章。等撰写的文章。SSA-MTM软件介绍软件介绍该软件由一系列的程序构成，主要用于该软件由一系列的程序构成，主要用于数据的谱分析和数据分解（可用于单变数据的谱分析和数据分解（可用于单变量或多变量）；量或多变量）；主要包括以下部分：主要包括以下部分：时间序列的谱估计；时间序列的谱估计；时间序列分解为趋势项、准周期振荡成分以时间序列分解为趋势项、准周期振荡成分以及噪音；及噪音；对分解后的时间成分进行重构。对分解后的时间成分进行重构。SSA-MTM软件介绍软件介绍该软件包括以下几种方法：该软件包括以下几种方法：“Traditional”Correlogram Analysis;Maximum-Entropy Spectral Estimates;Multi-Taper-Method Estimates;Singular-Spectrum Analysis;Multi-Channel Singular-Spectrum Analysis。SSA-MTM软件中的软件中的MEM方法方法数据介绍数据介绍为为Southern Oscillation Index(SOI)月资月资料料长度：长度：1942年年1月月-1999年年6月，共月，共690个月；个月；计算：资料由计算：资料由Tahiti和和Darwin站的月平均海站的月平均海平面气压计算得到：平面气压计算得到：分别将两站月平均海平面气压资料作月标准化分别将两站月平均海平面气压资料作月标准化处理；处理；在求在求Tahiti和和Darwin月标准化数据的差值。月标准化数据的差值。SSA-MTM软件中的软件中的MEM方法方法SSA-MTM软件中的软件中的MEM方法方法举例举例8存在存在5.3年和年和2.51年年2个显著的周期个显著的周期