气候序列的周期提取方法-谱分析

1、 时间序列分析方法时间序列分析方法 存在两种基本的时间序列分析方法：存在两种基本的时间序列分析方法： 时间域分析方法；时间域分析方法； 离散数据离散数据Markov Chains 连续数据连续数据自回归过程自回归过程 频率域分析方法；频率域分析方法； 时间序列分析方法类似于理论分布，即用几个时间序列分析方法类似于理论分布，即用几个 参数作为数据的代表，但理论分布并不考虑数参数作为数据的代表，但理论分布并不考虑数 据的排序特征，而这里的时间序列方法是对数据的排序特征，而这里的时间序列方法是对数 据的排序特征进行推断，从而也可用于对未来据的排序特征进行推断，从而也可用于对未来 数据特征的推断，这要

2、求数据应满足平稳性。数据特征的推断，这要求数据应满足平稳性。 时间序列分析方法时间序列分析方法 存在两种基本的时间序列分析方法：存在两种基本的时间序列分析方法： 时间域分析方法；时间域分析方法； 离散数据离散数据Markov Chains 连续数据连续数据自回归过程自回归过程 频率域分析方法；频率域分析方法； 频率域分析方法频率域分析方法 频率域分析方法根据不同时间尺度（或频率域分析方法根据不同时间尺度（或 频率）的贡献来表证数据；频率）的贡献来表证数据； 每一时间尺度可由一对每一时间尺度可由一对sine和和cosine函数函数 表示；表示； 整个的时间序列就是由不同尺度的整个的时间序列就是由

3、不同尺度的sine 和和cosine函数的叠加构成；函数的叠加构成； 通常我们对单个尺度的波更加的感兴趣。通常我们对单个尺度的波更加的感兴趣。 频率域分析方法频率域分析方法 因此频率域分析方法涉及到将包含因此频率域分析方法涉及到将包含n个点个点 的原始数据转化为一系列的周期函数；的原始数据转化为一系列的周期函数； 虽然，直观上频率域分析方法较难以让虽然，直观上频率域分析方法较难以让 人接受；人接受； 但是，这种方法在大气科学分析中是非但是，这种方法在大气科学分析中是非 常常用的，也是非常重要的方法，常能常常用的，也是非常重要的方法，常能 为我们提供原数据重要的信息。为我们提供原数据重要的信息。

4、 谐波分析谐波分析 （Harmonic analysis） 谐波分析（谐波分析（Harmonic analysis） 谐波分析是将一系列谐波分析是将一系列sine和和cosine函数叠函数叠 加在一起来表征原始数据的振荡或波动；加在一起来表征原始数据的振荡或波动； Cosine和和sine函数特点函数特点 cos(2)cos( ),intkk is anyeger sin(2)sin( ),intkk is anyeger cos()sin( ) 2 sin()cos( ) 2 一个谐波函数表征一个简单的一个谐波函数表征一个简单的 时间序列时间序列存在的问题存在的问题 即便是时间序列具有很好的

5、正弦曲线即便是时间序列具有很好的正弦曲线 的特征，我们用正弦曲线表征该数据的特征，我们用正弦曲线表征该数据 时，仍然存在以下问题：时，仍然存在以下问题： 1. 数据是时间的函数，而三角函数却是角度数据是时间的函数，而三角函数却是角度 的函数；的函数； 2. 余弦和正弦曲线波动的范围在余弦和正弦曲线波动的范围在+1和和-1之间，之间， 而数据的振荡范围通常不能满足此限制；而数据的振荡范围通常不能满足此限制； 3. 余弦曲线极值位于余弦曲线极值位于 处，而处，而 此位置对于正弦曲线则为平均值。此位置对于正弦曲线则为平均值。 02and 一个谐波函数表征一个简单的一个谐波函数表征一个简单的 时间序列

6、时间序列解决方法解决方法 1. 将数据记录长度将数据记录长度n看作为一个周期或基看作为一个周期或基 本周期，则有：本周期，则有： 为基本频率，下标为基本频率，下标1表示整个数表示整个数 据只有一个完整的循环。据只有一个完整的循环。 360 360 / 2 2 / t timeunitst cyclen timeunits cylcen t timeunitst or cyclen timeunits cylcen 1 2 n 一个谐波函数表征一个简单的一个谐波函数表征一个简单的 时间序列时间序列解决方法解决方法 2. 将将cosine或或sine函数向上或向下移动到函数向上或向下移动到 原数据

7、的基本水平处，然后拉伸或压原数据的基本水平处，然后拉伸或压 缩到与原数据相同的振幅范围。缩到与原数据相同的振幅范围。 ： 振幅，则振幅的最大和最小值为振幅，则振幅的最大和最小值为 1 2 cos() t t yyC n 1 C 1 C 一个谐波函数表征一个简单的一个谐波函数表征一个简单的 时间序列时间序列解决方法解决方法 3. 应将谐波函数水平的移动，从而与原应将谐波函数水平的移动，从而与原 序列的脊和槽相匹配。序列的脊和槽相匹配。 为位相角，即将为位相角，即将cosine函数向右移动函数向右移动 角度角度 ，则新的函数在，则新的函数在 处达处达 到最大值。到最大值。 11 2 cos() t

8、 t yyC n 1 1 1 2 n t 举例举例1 1943-1989年年Ithaca平均每月温度；平均每月温度； t=1表示表示1月，月，t=2表示表示2月，月，； 整个序列的年平均温度为整个序列的年平均温度为46.1 ； 原数据近似正弦曲线；原数据近似正弦曲线； 最暖的最暖的7月平均温度为月平均温度为68.8 ，最冷的，最冷的1 月平均温度为月平均温度为22.2 单个谐波的振幅和位相估计单个谐波的振幅和位相估计 由三角函数的特性：由三角函数的特性： 知：知： 111 cos()cos( )cos( )sin( )sin( ) 111111 11 111 111 222 cos()cos(

9、 )cos()sin( )sin() 22 cos()sin() cos( ) sin( ) ttt CCC nnn tt AB nn AC BC 谐波与多元线性回归谐波与多元线性回归 当对上式进行变量代换可转化为一般的当对上式进行变量代换可转化为一般的 多元线性回归方程：多元线性回归方程： 由此，可利用最小二乘法估计参数，且：由此，可利用最小二乘法估计参数，且： 12 1112 22 cos(),sin() , tt xx nn AbBb 22 1/2 111 CAB 位相角位相角 计算公式：计算公式： 且满足且满足 1 111 1 1111 1 tan (/),0 tan (/),0 /2

10、,0 BAA BAA A 1 02 参数参数 在满足时间步长相等，且无缺测的前提在满足时间步长相等，且无缺测的前提 下，通过最小二乘估计得到参数：下，通过最小二乘估计得到参数： 1 1 1 1 22 cos() 22 sin() n t t n t t t Ay nn t By nn 举例举例2 同上例，有下表同上例，有下表 t ycos(2/12)tsin(2/12)t sin(2/12) t ytcos(2/12) t yt -86.417-110.3290.0000.000552.9Sums: 0.00027.4000.0001.00027.412 -19.65034.034-0.500

11、0.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.87 0.000-64.3000.000-1.00064.36 27.400-47.4570.500-0.86654.85 38.450-22.2000.866+0.50044.44 32.2000.0001.0000.00032.23 19.65811.3500.8660.50022.72 11.10019.2250.5000.8

12、6622.21 t 举例举例2 则可计算得到平均值为则可计算得到平均值为46.1 与前述例子的结果非常接近，但位相角与前述例子的结果非常接近，但位相角 向右移动向右移动8度（大约度（大约1周的时间），则结周的时间），则结 果更接近实际数据。果更接近实际数据。 1 1 (2/12)( 110.329)18.39 (2/12)( 86.417)14.40 A B 221/2o 1 18.39( 14.40 )23.36 FC 1 1 tan ( 14.40/ 18.39) 180218 oo 高阶谐波高阶谐波 由于年循环近似于正弦曲线，因此单个由于年循环近似于正弦曲线，因此单个 谐波便可很好的拟合

13、；谐波便可很好的拟合； 但这并不意味着单个谐波可以很好的表但这并不意味着单个谐波可以很好的表 征任何时间序列；征任何时间序列； 类似于多元回归问题，增加更多的余弦类似于多元回归问题，增加更多的余弦 波也会提高谐波分析的拟合结果，但也波也会提高谐波分析的拟合结果，但也 同样面临可能会过度拟合的问题。同样面临可能会过度拟合的问题。 高阶谐波高阶谐波 已证明包含已证明包含n个点的原始数据通过叠加个点的原始数据通过叠加 n/2个谐波函数便可以完全表征原数据，个谐波函数便可以完全表征原数据， 即存在一个可以通过所有点的谐波。即存在一个可以通过所有点的谐波。 频率（圆频率）为：频率（圆频率）为： /2 1

14、 /2 1 2 cos 22 cossin n tkk k n kk k kt yyC n ktkt yAB nn 2 k k n 解释解释 K=2，二阶谐波，其振幅和位相分别为，二阶谐波，其振幅和位相分别为 C2和和 ； 表示时间表示时间t从从0到到n/2有一个循环，而从有一个循环，而从n/2 到到n还有一个完整的循环，即还有一个完整的循环，即2波；波； 同样，同样，k=3表示整个时间过程存在表示整个时间过程存在3个波。个波。 2 参数参数 在时间步长均等，且无缺测点的前提下，在时间步长均等，且无缺测点的前提下， 通过求解多元回归方程可得到：通过求解多元回归方程可得到： 1 1 22 cos

15、() 22 sin() n kt t n kt t kt Ay nn kt By nn 22 1/2 kkk CAB 1 1 tan (/),0 tan (/),0 /2,0 kkk kkkk k BAA BAA A 离散离散Fourier变换变换 Ak和和Bk为为Fourier系数系数 过度拟合过度拟合 有关多元回归分析表明，当拟合线通过有关多元回归分析表明，当拟合线通过 所有数据点时，复相关系数为所有数据点时，复相关系数为100%，为，为 过度拟合。过度拟合。 同样，当谐波方程中包含同样，当谐波方程中包含n/2个谐波时，个谐波时， 也为过度拟合。也为过度拟合。 过度拟合过度拟合 由于每个谐

16、波项都包含由于每个谐波项都包含2个参数，即振幅个参数，即振幅 和位相，同时方程中还包含截距，即样和位相，同时方程中还包含截距，即样 本平均值：本平均值： 当当n为偶数时，可以有为偶数时，可以有n/2个谐波，个谐波， ， n/2个振幅个振幅+ (n/2-1)个位相个位相+样本平均值，方样本平均值，方 程包含有程包含有n个参数；个参数； 当当n为奇数时，可以有为奇数时，可以有(n-1)/2个谐波，参数个谐波，参数 为为(n-1)/2个振幅个振幅+ (n-1)/2个位相个位相+样本平均样本平均 值，则方程仍包含有值，则方程仍包含有n个参数。个参数。 /2 0 n 谐波个数的选取谐波个数的选取 通常并

17、不需要将所有的通常并不需要将所有的n/2个谐波均用于个谐波均用于 拟合原数据，而是仅用几个谐波便可以拟合原数据，而是仅用几个谐波便可以 很好的表征原数据的变化特征；很好的表征原数据的变化特征； 仅当分析的目的是准备表征原数据时，仅当分析的目的是准备表征原数据时， 选用选用n/2个谐波进行分析。个谐波进行分析。 举例举例3 原数据为原数据为El PASO地区地区1948-1983年连续年连续 5天无降水的相对频率的年变化；天无降水的相对频率的年变化； 举例举例3 由原数据可见夏季较为湿润，而春、秋由原数据可见夏季较为湿润，而春、秋 两季较为干燥两季较为干燥非对称的年循环；非对称的年循环； 同时叠

18、加有不规则的、短期波动；同时叠加有不规则的、短期波动； 采用前采用前3个谐波来拟合原数据，得到平均个谐波来拟合原数据，得到平均 值以及前两个谐波参数为：值以及前两个谐波参数为： 11 22 61.4% 13.6%,720.4 13.8%,2721.51 o o y C C 举例举例3图图 谱分析谱分析 （Spectral analysis） 谐波方程与多元回归方程差异谐波方程与多元回归方程差异 谐波方程与多元回归方程不同在于每个谐波方程与多元回归方程不同在于每个 谐波之间是彼此独立的，因此每个谐波谐波之间是彼此独立的，因此每个谐波 参数可以独立的计算出来，当有新的谐参数可以独立的计算出来，当有

19、新的谐 波加入方程时，则已在谐波方程中的谐波加入方程时，则已在谐波方程中的谐 波参数保持不变。波参数保持不变。 谐波函数彼此独立谐波函数彼此独立 谐波函数彼此独立的特性来自于谐波函数彼此独立的特性来自于sine和和 cosine函数彼此正交性；函数彼此正交性； 对于任一整数对于任一整数k和和j有：有： 1 11 22 cos()sin()0 2222 cos()cos()sin()sin()0, n t nn tt ktjt nn ktjtktjt for kj nnnn 谐波函数彼此独立谐波函数彼此独立 例如考虑例如考虑2个预报因子：个预报因子： 它们的相关为：它们的相关为： 13 cos2

20、/ ,cos2 (2 )/ xt nxtn 1 3 1133 1 22 1/2 1133 11 1 221/2 11 ()() ()() 22 2 cos()cos() 0 22 2 cos ()cos () n t x xnn tt n t nn tt xxxx r xxxx tt nn tt nn 单个谐波方差贡献单个谐波方差贡献 由于每个谐波之间彼此独立，因此它们由于每个谐波之间彼此独立，因此它们 对方程的方差贡献并不随谐波方程的变对方程的方差贡献并不随谐波方程的变 化而变化，如对第化而变化，如对第k个谐波而言，其方差个谐波而言，其方差 贡献为：贡献为： 为原序列的样本方差为原序列的样本

21、方差 2 2 2 2 (1) k k y n C R ns 2 y s 谐波总方差贡献谐波总方差贡献 谐波方程中所有谐波的方差贡献为：谐波方程中所有谐波的方差贡献为： 若方程中有若方程中有n/2个可能的谐波，则：个可能的谐波，则： 22 k k in the equation RR 2 1R 离散功率谱离散功率谱 （周期图（周期图/Fourier线谱）线谱） 概述概述 从谱分析的观点看，一个时间序列的谱从谱分析的观点看，一个时间序列的谱 可以显示出不同频率的振荡对该序列变可以显示出不同频率的振荡对该序列变 化的贡献；化的贡献； Panofsky和和Brier(1958)类似的解释为：类似的解释

22、为： 一个光谱显示出不同波长或频率的光对一个光谱显示出不同波长或频率的光对 给定光源能量的贡献。给定光源能量的贡献。 功率谱图功率谱图 由于振幅是频率的函数，因此最简单的由于振幅是频率的函数，因此最简单的 功率谱图由振幅的平方构成功率谱图由振幅的平方构成离散功离散功 率谱；率谱； 也可以为标准谱密度，即也可以为标准谱密度，即 由于功率谱图虽然提供了不同频率谐波由于功率谱图虽然提供了不同频率谐波 对原数据的贡献如何，但并没有提供位对原数据的贡献如何，但并没有提供位 相角的信息，即没有提供不同频率谐波相角的信息，即没有提供不同频率谐波 随时间的变化信息，从而无法重构时间随时间的变化信息，从而无法重

23、构时间 序列。序列。 2 k R 功率谱图功率谱图 通常功率谱图的纵坐标为对数坐标，这通常功率谱图的纵坐标为对数坐标，这 样有利于更好的绘制出对主要由少数几样有利于更好的绘制出对主要由少数几 个频率的谐波主导的时间序列的谱图；个频率的谐波主导的时间序列的谱图； 否则对能量贡献较小的其它谐波将很难否则对能量贡献较小的其它谐波将很难 直接从图中看出其贡献量。直接从图中看出其贡献量。 基本频率和基本频率和Nyquist频率频率 根据频率定义：根据频率定义： 当当k=1时，基本频率，时，基本频率， 当当k=n/2时，时，Nyquist频率，频率， 2 k k n 1 2 n /2n 功率谱图横坐标功率

24、谱图横坐标 功率谱图横坐标可以是频率功率谱图横坐标可以是频率 也可以是频率也可以是频率 则频率可以从基本频则频率可以从基本频1/n到到Nyquist频率频率1/2 可以是周期可以是周期 还可以是波数还可以是波数k k k f 2 k k k f n k 21 k kk n kf 举例举例4离散离散Fourier变换变换 Month19871988kTk(months)AkBkCk 121.420.6124-0.140.440.46 217.922.5212-23.76-2.2023.86 335.932.938-0.990.391.06 447.743.646-0.46-1.251.33 55

25、6.456.554.8-0.02-0.430.43 666.361.964-1.49-2.152.62 770.971.673.43-0.53-0.070.53 865.869.983-0.34-0.210.40 960.157.992.671.560.071.56 1045.445.2102.40.130.220.26 1139.540.5112.180.520.110.53 1231.326.71220.790.79 举例举例4 1987和和1988年年Ithaca地区观测的月平均地区观测的月平均 温度；温度； 在未作谱分析前，我们已知该数据的基在未作谱分析前，我们已知该数据的基 本特征为

26、存在年循环；本特征为存在年循环； 即周期近似为即周期近似为12个月的循环。个月的循环。 离散离散Fourier变化结果变化结果 举例举例4结果分析结果分析 n=24，可以由，可以由12个不同频率的谐波构成；个不同频率的谐波构成； 贡献最大是周期为贡献最大是周期为12个月的谐波，即年个月的谐波，即年 循环；循环； 表中表中B12=0，23个个Fourier系数，外加样系数，外加样 本平均值；本平均值； 举例举例5 分析分析1951-1979年年SOI指数功率谱；指数功率谱； 原序列存在准周期变化；原序列存在准周期变化； 举例举例5离散功率谱图离散功率谱图 (线谱线谱/离散谱离散谱) 举例举例5结

27、果分析结果分析 谱图显示出谱图显示出SOI指数序列存在典型的指数序列存在典型的3-7 年的振荡；年的振荡； 谱图并未画出所有频率的贡献；谱图并未画出所有频率的贡献； 较低频变化反映出我们感兴趣的物理现较低频变化反映出我们感兴趣的物理现 象，即象，即ENSO循环。循环。 位相谱位相谱 在离散功率谱估计中，也可以画出位相在离散功率谱估计中，也可以画出位相 随频率随频率/波数波数/周期的变化，称之为位相谱周期的变化，称之为位相谱 ( ) 位相谱能为我们提供不同频率分量的初位相谱能为我们提供不同频率分量的初 始信息，因此包含了信号的大量信息；始信息，因此包含了信号的大量信息； 用相同的振幅谱和不同的位

28、相谱重构的用相同的振幅谱和不同的位相谱重构的 信号是完全不同的。信号是完全不同的。 k 离散功率谱计算离散功率谱计算 由由 计算计算Fourier系数，进而得系数，进而得 到振幅；到振幅； 但上述方式仅适用于但上述方式仅适用于 ； 通常对通常对k=n/2用下式计算：用下式计算： 1 1 22 cos() 22 sin() n kt t n kt t kt Ay nn kt By nn 1,2,( /2 1)kn 11 /2 /2 122 ( /2)1 coscos, 2 0, 0, nn tt tt n n nt yytn even Annn n odd Bn even or odd 离散功率

29、谱计算离散功率谱计算 通常为了使起始时刻位于通常为了使起始时刻位于0时刻，通常离时刻，通常离 散功率谱计算的不同波数散功率谱计算的不同波数k的的Fourier系系 数估计为：数估计为： 且用原序列或距平序列计算结果是一致且用原序列或距平序列计算结果是一致 的。的。 1 1 22 cos(1) 22 sin(1) n kt t n kt t k Ayt nn k Byt nn 离散功率谱计算离散功率谱计算 实际上，采用上述方式计算离散实际上，采用上述方式计算离散Fourier 变换是非常低效的；变换是非常低效的； 原因是存在多次反复调用完全相同的计原因是存在多次反复调用完全相同的计 算过程，如上

30、例中：算过程，如上例中： t=4，k=1， t=4，k=2， (47.7)sin(2 )(1)(4)/24(47.7)(0.866)41.31 ooo FFF (47.7)sin(2 )(2)(4)/24(47.7)(0.866)41.31 ooo FFF 谱计算谱计算快速快速Fourier变换变换 (Fast Fourier Transforms, FFTs) 采用采用FFTs可以解决上述问题；可以解决上述问题； 很多计算软件或程序中包含有一个或多很多计算软件或程序中包含有一个或多 个个FFT程序；程序； 对于较长序列，对于较长序列，FFT相对于回归方法计相对于回归方法计 算速度提高很多，对

31、于算速度提高很多，对于n=100，FFT的速的速 度是后者的度是后者的100倍，而倍，而n=10000，则提高，则提高 750倍。倍。 FFTs 通常通常FFTs采用的是欧拉复指数计算方程：采用的是欧拉复指数计算方程： 这样处理的方式主要是基于表达方式简这样处理的方式主要是基于表达方式简 便，且易于处理，而其数学含义与前述便，且易于处理，而其数学含义与前述 形式是完全一致的。形式是完全一致的。 cos()sin() i t etit 复指数形式下的谐波方程复指数形式下的谐波方程 则原谐波方程变为：则原谐波方程变为： 是复是复Fourier系数系数 /2 2/ 1 n ik n t tk k y

32、yH e kkk HAiB k H Aliasing(假名假名)问题问题 功率谱分析周期存在的一个问题就是功率谱分析周期存在的一个问题就是 “假名假名”问题；问题； 在功率谱分析中，可能序列存在的重要在功率谱分析中，可能序列存在的重要 的物理过程频率比的物理过程频率比Nyquist频率更高频，频率更高频， 但这些短周期振荡无法由直接分辨出来，但这些短周期振荡无法由直接分辨出来， 则它们的作用会体现在较长周期（则它们的作用会体现在较长周期（ 和和 之间）中；之间）中； 即这些高频变化被即这些高频变化被“冒名顶替冒名顶替”了，就了，就 是是“假名假名”现象。现象。 1 /2n Aliasing形成

33、原因形成原因 f=4/5f=1/5 Aliasing形成原因形成原因 上图说明实际的时间序列的物理过程包上图说明实际的时间序列的物理过程包 含比含比Nyquist频率高的变化部分；频率高的变化部分； 当对数据的取样少，即数据间隔过大而当对数据的取样少，即数据间隔过大而 不能体现这种快变化；不能体现这种快变化； 而这种比而这种比Nyquist频率高的频率并不会因频率高的频率并不会因 此从数据真实变化中消失；此从数据真实变化中消失； 同时这些高频的作用便虚假的体现在可同时这些高频的作用便虚假的体现在可 分辨的频率中。分辨的频率中。 Aliasing形成原因形成原因 在时间序列中可能真实存在比在时间

34、序列中可能真实存在比Nyquist频频 率高的变化部分；率高的变化部分； 由于这部分高频信号无法由谱分析得到，由于这部分高频信号无法由谱分析得到， 从而这些高频信号的能量会体现在频率从而这些高频信号的能量会体现在频率 介于介于0到到Nyquist频率之间的信号中；频率之间的信号中； 也就是说也就是说0到到Nyquist频率顶替了其它高频率顶替了其它高 频信号的作用。频信号的作用。 Aliasing形成原因形成原因 造成假名现象的根本原因是数据离散的造成假名现象的根本原因是数据离散的 观测结果，最短的周期大小直接由数据观测结果，最短的周期大小直接由数据 取样间隔决定，因此，很可能在取样间取样间隔

35、决定，因此，很可能在取样间 隔中存在一定的周期变化，而无法在得隔中存在一定的周期变化，而无法在得 到的数据中体现出来。到的数据中体现出来。 Aliasing如何体现如何体现 无法分辨的高频变化能量会错误的添加无法分辨的高频变化能量会错误的添加 到可分辨的到可分辨的n/2个可分辨的频率变化上；个可分辨的频率变化上； 对于无法分辨的高频频率对于无法分辨的高频频率 ，可体，可体 现在可分辨的频率部分现在可分辨的频率部分 ，有：，有： ,int. 2,int. A A fjfj any postiveeger jj any positveeger AA f or f or Aliasing如何体现如何

36、体现 上式的含义是，高频率的会以上式的含义是，高频率的会以Nyqusit频频 率的整数倍作为每次折叠的率的整数倍作为每次折叠的“折叠折叠”点，点， 对称的叠加到可分辨的频率部分，因此对称的叠加到可分辨的频率部分，因此 Nyqusit频率又称为：频率又称为：“folding” frequency； 如：如： 意味着，稍高于意味着，稍高于 的频率会叠加在接近低频的部分，而稍的频率会叠加在接近低频的部分，而稍 低于低于 的频率会叠加在接近的频率会叠加在接近Nyquist频频 率部分；率部分； 22 23 nAn fff 2 A f 3 A f Aliasing现象举例现象举例 如何避免如何避免Ali

37、asing 一但我们选定数据，则没有办法去除一但我们选定数据，则没有办法去除“Aliasing” 现象；现象； 当然，我们可以通过提高数据分辨率的方法尽当然，我们可以通过提高数据分辨率的方法尽 可能的避免该现象；可能的避免该现象； 或者根据已知的物理过程，来确定资料或者根据已知的物理过程，来确定资料/样本样本 的变化率，从而达到去除假名现象；的变化率，从而达到去除假名现象； 但对于探索性的研究，即不知道物理本质的问但对于探索性的研究，即不知道物理本质的问 题，是没有办法去除假名现象的，但我们期望题，是没有办法去除假名现象的，但我们期望 接近接近Nyquist频率的功率谱能量接近频率的功率谱能量

38、接近0，则可能，则可能 说明高频部分的能量很小，而上图不符合该假说明高频部分的能量很小，而上图不符合该假 设设。 连续功率谱估计连续功率谱估计 上述有关谱计算的方法是基于离散功率上述有关谱计算的方法是基于离散功率 谱得到的；谱得到的； 另一种计算方法为连续功率谱估计方法；另一种计算方法为连续功率谱估计方法； 这种方法基于时间序列的自相关过程作这种方法基于时间序列的自相关过程作 间接的估计。间接的估计。 连续功率谱估计连续功率谱估计 对于以任一周期对于以任一周期T变化的时间函数变化的时间函数 ， 可以展开为傅立叶级数：可以展开为傅立叶级数： ( )x t 0 1 0 1 22 ( )cossin

39、 2 ( )cossin kk k k kkkk k kk x taatbt TT k T x taatbt 连续功率谱估计连续功率谱估计 上式参数为：上式参数为： 2 20 2 2 2 2 1 ( ) 2 ( )cos 2 ( )sin T T T kk T T kk T ax t dt T ax ttdt T bx ttdt T 连续功率谱估计连续功率谱估计 通常在气候序列分析中，我们希望了解通常在气候序列分析中，我们希望了解 在无限时间序列中各种连续频率的波动在无限时间序列中各种连续频率的波动 结构，因此有必要把离散谱推广到连续结构，因此有必要把离散谱推广到连续 谱上面，即把周期谱上面，

40、即把周期T推广至无限大；推广至无限大； 连续功率谱估计连续功率谱估计 将参数带入傅氏级数公式，并利用三角将参数带入傅氏级数公式，并利用三角 函数和差化积公式有：函数和差化积公式有： 2 1 2 12 ( )( )( )cos() T Tk k x txdxtd TT 连续功率谱估计连续功率谱估计 Fourier思想：一个非周期信号可以看成思想：一个非周期信号可以看成 是周期无限长的周期信号，当周期增加是周期无限长的周期信号，当周期增加 时，基本频越小，则成谐波关系的各分时，基本频越小，则成谐波关系的各分 量在频率上就会越来越接近，当周期变量在频率上就会越来越接近，当周期变 的无穷大时，离散的线

41、谱就形成为一个的无穷大时，离散的线谱就形成为一个 连续谱，即从求和变成了积分。连续谱，即从求和变成了积分。 即：任意两个相邻谐波之间的频率差趋即：任意两个相邻谐波之间的频率差趋 近于近于0 1 22 (1)2 kk kk TTT 连续功率谱估计连续功率谱估计 离散功率谱只给出了频谱在离散点离散功率谱只给出了频谱在离散点 （ ）上的值，而无法反映这些点）上的值，而无法反映这些点 之间的频谱内容；之间的频谱内容； 把离散功率谱绘制成连续的谱线，并不把离散功率谱绘制成连续的谱线，并不 表示它就是连续功率谱；表示它就是连续功率谱； 因此必须采用连续功率谱计算这些点的因此必须采用连续功率谱计算这些点的

42、谱值。谱值。 2 k k n 连续功率谱估计连续功率谱估计 在把周期在把周期T推广至无穷大过程中，任意推广至无穷大过程中，任意 两个相邻波数对应的谐波间的频率差无两个相邻波数对应的谐波间的频率差无 限接近于限接近于0，则有：，则有： 22 1 22 0 2 ( )( )( )cos () 2 ,0 1 ( )( )cos () TT TT k T x txdxtd T x tdxtd 连续功率谱估计连续功率谱估计 则进一步令：则进一步令： 并设：并设： ( )( )cos ( )( )sin axd bxd 22 22 ( ) sin ( ) ( )( ) ( ) cos ( ) ( )( )

43、 b ab a ab 连续功率谱估计连续功率谱估计 连续谱中的位相谱为连续谱中的位相谱为 ： 并记振幅谱为：并记振幅谱为： 则有：则有： 222 ( )( )( )Aab 0 0 0 0 0 1 ( )( )cos () 1 ( )coscos)( )sinsin) 1 ( )cos( )sin 1 ( )cos ( )cossin ( )sin 1 ( )cos ( ) 1 ( )cos( x tdxtd dxt dxt d atbt d Att d At d At 0 ) d ( ) 连续功率谱估计连续功率谱估计 利用欧拉公式，上式变为利用欧拉公式，上式变为 则有：则有： 令令 则有：则有

44、： 1 ( )( ) 2 i t x tFed ()() 0 () 1 ( )( ) 2 1 ( ) 2 itit ii t x tAeed Aeed () ( )( ) i FAe 连续功率谱估计连续功率谱估计 称为时间函数称为时间函数 的谱，可的谱，可 以表示为：以表示为： 带入带入 的表达式，则有：的表达式，则有： ( )x t ( )( )cos( )sin ( )(cossin) ( ) i t Fxdixd x ttit dt x t edt ( )F () ( )( )( )(cos ( )sin ( ) ( )( ) i FAeAi aib ( ), ( )ab 连续功率谱估计

45、连续功率谱估计 上述时间函数与它的复谱的关系本来上述时间函数与它的复谱的关系本来 可以作为时间序列分析的工具，但由可以作为时间序列分析的工具，但由 于复谱是复数，使用不方便；于复谱是复数，使用不方便； 而且，我们感兴趣的仅是某个周期波而且，我们感兴趣的仅是某个周期波 动所起的主要作用，即仅对它们的方动所起的主要作用，即仅对它们的方 差贡献感兴趣，因此在时间序列分析差贡献感兴趣，因此在时间序列分析 中常用的是功率谱，也可以称为方差中常用的是功率谱，也可以称为方差 谱或能谱；谱或能谱； 连续功率谱估计连续功率谱估计 由前述分析可知，功率谱分析的是不同由前述分析可知，功率谱分析的是不同 周期波动的方

46、差贡献；周期波动的方差贡献； 来源于物理学中功率的概念；来源于物理学中功率的概念； 若电阻为若电阻为1个单位，瞬时电压为个单位，瞬时电压为 ， 则它的瞬时功率为则它的瞬时功率为 ，它的总能量，它的总能量 为：为： ( )x t 2 ( )x t dt 2 ( )x t 连续功率谱估计连续功率谱估计 对于数学期望为零的时间序列对于数学期望为零的时间序列 ，其连，其连 续功率谱可表示为：续功率谱可表示为： ( )x t 2 1 ( )( )( ) 2 1 ( )( ) 2 1 ( ) () 2 i t i t x tdtx tFeddt Fdx t edt FFd 连续功率谱估计连续功率谱估计 对

47、于式子：对于式子： 其中其中 为共轭复谱为共轭复谱 功率谱密度：功率谱密度： 2 1 ( )( ) () 2 x tdtFFd ( )( ) i t Fx t edt * ()( )FF * ( )F 2 * ( )( )( )( )SFFF 连续功率谱估计连续功率谱估计 连续功率谱估计并不是由上式计算得到，连续功率谱估计并不是由上式计算得到， 而是通过时间函数的自相关函数间接估而是通过时间函数的自相关函数间接估 计得到；计得到； 对于标准化时间函数对于标准化时间函数 ，其数学期望，其数学期望 为为0，方差为，方差为1，则其自相关函数可表示，则其自相关函数可表示 为：为： ( )x t ( )

48、( ) ()x t x tdt 连续功率谱估计连续功率谱估计 对于公式：对于公式： 其中其中 为落后时间步长为落后时间步长 上式可表述为：上式可表述为： 由由Fourier反变换得：反变换得： ( )( ) ()x t x tdt 1 ( )( ) 2 i Sed 0 ( )( )( )cos( )sin ( )cos2( )cos i Seddid dd 谱计算谱计算连续功率谱估计连续功率谱估计 第一步：计算样本自相关函数；第一步：计算样本自相关函数； m为最大落后时间长度为最大落后时间长度 第二步：求粗谱估计；第二步：求粗谱估计； 利用数值积分中的梯形法得到：利用数值积分中的梯形法得到：

49、( ) (0,1,2,)rm 1 1 1 11 2(0)( )cos( )cos 22 m kkk Srrr mm m 谱计算谱计算连续功率谱估计连续功率谱估计 由于数据本身是离散点，则同样满足频由于数据本身是离散点，则同样满足频 率从率从0 到到0.5的变化；的变化； 若取样点数为若取样点数为m，则波数，则波数k最大可取到最大可取到 m/2； 若令若令l=2k，则，则l取值为从取值为从0到到m之间的任之间的任 一整数，有粗谱估计：一整数，有粗谱估计： 1 1 1 (0)2( )cos( )cos m l l Srrr ml mm 谱计算谱计算连续功率谱估计连续功率谱估计 对上式作平滑处理，以

50、消除粗谱估计的对上式作平滑处理，以消除粗谱估计的 小波动；小波动； 采用二项系数平滑得到：采用二项系数平滑得到： 001 11 1 11 22 111 (11) 424 11 22 llll mmm SSS SSSSlm SSS 谱计算谱计算连续功率谱估计连续功率谱估计 上述平滑效果相当于原粗谱公式乘以一上述平滑效果相当于原粗谱公式乘以一 个窗函数个窗函数 则有最终功率谱计算式：则有最终功率谱计算式： 1 1cos 2m 1 1 (0)( )(1 cos)cos (0,1,) 1 (0,) 1 (0,) 2 m l l l Bl Srrlm mmm lm B lm 绘图绘图 以波数以波数l为横

51、坐标，平滑功率谱密度估计为横坐标，平滑功率谱密度估计 值值Sl为纵坐标绘图；为纵坐标绘图； 横坐标也可以为频率和周期：横坐标也可以为频率和周期： 2 l l l l m m T l l f m 注意事项注意事项 连续功率谱估计的方法并不只有上述方法，自连续功率谱估计的方法并不只有上述方法，自 相关函数和平滑方案都可以有一些差别；相关函数和平滑方案都可以有一些差别； 计算过程中人为的选择最大落后长度计算过程中人为的选择最大落后长度m，会改，会改 变谱估计曲线：变谱估计曲线： m取值小，虽然谱曲线比较光滑，但不容易确定峰取值小，虽然谱曲线比较光滑，但不容易确定峰 点，即不好分辨出主要周期；点，即不

52、好分辨出主要周期； m取值大，谱曲线本身的波动较多，但并不意味着取值大，谱曲线本身的波动较多，但并不意味着 这些峰值对应的频率存在确定的周期；这些峰值对应的频率存在确定的周期； 因此，通常因此，通常m的取值范围在序列长度的取值范围在序列长度n的的1/10- 1/3之间。之间。 举例举例7连续功率谱估计连续功率谱估计 00.10.20.30.40.5 F 0.0001 0.001 0.01 0.1 S(l) 举例举例6离散功率谱估计离散功率谱估计 01234567891011 K 0 0.4 0.8 1.2 1.6 C2k 离散功率谱检验离散功率谱检验 原假设：原假设： 有检验统计量：有检验统计

53、量： 遵从分子自由度为遵从分子自由度为2，分母自由度为，分母自由度为 n-2-1的的F分布。分布。 ()()0 kk E aE b 22 222 1 ()/2 2 11 ()/(2 1) 22 kk kk ab F sabn 举例举例6检验检验 01234567891011 K 0 0.4 0.8 1.2 1.6 C2k 连续功率谱检验连续功率谱检验 通常，在连续功率谱图上，功率谱图的通常，在连续功率谱图上，功率谱图的 峰点（极大值）所对应的周期可定为主峰点（极大值）所对应的周期可定为主 要周期；要周期； 但这些周期是否显著，要通过显著性检但这些周期是否显著，要通过显著性检 验方可确定；验方可

54、确定； 显著性检验通常是与非周期性随机过程显著性检验通常是与非周期性随机过程 作比较，即白噪音和红噪音过程。作比较，即白噪音和红噪音过程。 红噪音过程的功率谱红噪音过程的功率谱 在红噪音过程中每一时刻的观测值仅与在红噪音过程中每一时刻的观测值仅与 前一时刻的观测值有关，即一阶自回归前一时刻的观测值有关，即一阶自回归 过程；过程； 红噪音过程功率谱谱密度为：红噪音过程功率谱谱密度为： 2 2 1 ( ) 1 2 cos S 白噪音过程的功率谱白噪音过程的功率谱 白噪音过程是一无周期性的纯随机过程；白噪音过程是一无周期性的纯随机过程； 其自相关函数为：其自相关函数为： 其功率谱密度：其功率谱密度：

55、 ( )( )(0)1 i Sed 10 ( ) 00 功率谱的显著性检验功率谱的显著性检验过程过程1 原假设：所要检验的气象要素的总体谱原假设：所要检验的气象要素的总体谱 为某一非周期的随机过程的谱；为某一非周期的随机过程的谱； 建立检验统计量，某一频率上的样本谱建立检验统计量，某一频率上的样本谱 估计值与假设过程的平均谱估计值之比估计值与假设过程的平均谱估计值之比 遵从被其自由度去除的遵从被其自由度去除的 分布：分布： 2 2 0 k k S S 功率谱的显著性检验功率谱的显著性检验过程过程2 在上式中，自由度的计算和样本容量在上式中，自由度的计算和样本容量n及及 所取的最大落后长度所取的

56、最大落后长度m有关，即：有关，即： 在显著性水平为在显著性水平为0.05下，拒绝原假设的下，拒绝原假设的 要求是：要求是： (2)/ 2 m nm 2 0.05 0kk SS 功率谱的显著性检验功率谱的显著性检验过程过程3 判断用那一种过程来检验，有判据：判断用那一种过程来检验，有判据： 在显著性水平为在显著性水平为0.05下，有临界值：下，有临界值： 若实际计算的落后一个时刻的自相关系数若实际计算的落后一个时刻的自相关系数 有关系：有关系： ，即序列存在明显的持续性，即序列存在明显的持续性， 则用红噪声过程检验，否则用白噪声过程则用红噪声过程检验，否则用白噪声过程 检验。检验。 0 1 1.

57、6452 1 n r n 0 (1)rr 功率谱的显著性检验功率谱的显著性检验过程过程4 对于红噪声过程检验有平均红噪声谱为：对于红噪声过程检验有平均红噪声谱为： 其中：其中： 为估计的样本平均谱值。为估计的样本平均谱值。 2 0 2 1(1) (0,1,) 1(1)2 (1)cos k r SSkm k rr m S 功率谱的显著性检验功率谱的显著性检验过程过程5 对于白噪声过程检验有平均白噪声谱为：对于白噪声过程检验有平均白噪声谱为： 0k SS 举例举例7检验检验 存在存在5年和年和2.5年年2个显著的周期个显著的周期 00.10.20.30.40.5 f 0.0001 0.001 0.

58、01 0.1 S(l) 最大熵谱最大熵谱 (Maximum-Entropy Spectral) 概述概述连续功率谱估计缺点连续功率谱估计缺点 在连续功率谱的谱图中，功率谱值随频率变化在连续功率谱的谱图中，功率谱值随频率变化 是比较平滑的；是比较平滑的； 如果我们分析时间序列的主要目的是寻找序列如果我们分析时间序列的主要目的是寻找序列 的主要周期，则总希望所寻找的主要周期所对的主要周期，则总希望所寻找的主要周期所对 应的功率谱值比较突出，比较大，从而容易分应的功率谱值比较突出，比较大，从而容易分 辨出较精确的周期；辨出较精确的周期； 此外，气象要素的时间序列，通常受观测资料此外，气象要素的时间序

59、列，通常受观测资料 的限制，容量较小，而在连续功率谱估计中，的限制，容量较小，而在连续功率谱估计中， 自相关函数估计与样本容量有关，常会带来真自相关函数估计与样本容量有关，常会带来真 正谱估计的误差。正谱估计的误差。 概述概述最大熵谱最大熵谱 自自1980年以来，最大熵谱分析方法逐渐年以来，最大熵谱分析方法逐渐 被人们所认识，并得到广泛应用；被人们所认识，并得到广泛应用； 该方法的特点是高分辨率和短时性，即该方法的特点是高分辨率和短时性，即 只需较少的资料就可以得到较高的分辨只需较少的资料就可以得到较高的分辨 率；率； 且该方法提取的主次周期更符合实际。且该方法提取的主次周期更符合实际。 概述

60、概述熵熵 “熵熵”（entropy）是德国物理学家）是德国物理学家 Rudolf Clausius在在1850年创造的一个术年创造的一个术 语，语， 他用这个词来表示任何一种能量在他用这个词来表示任何一种能量在 空间中分布的均匀程度，能量分布的越空间中分布的均匀程度，能量分布的越 均匀，熵就越大，如果我们所考虑的系均匀，熵就越大，如果我们所考虑的系 统，能量是完全均匀分布的，则这个系统，能量是完全均匀分布的，则这个系 统的熵就达到最大值。统的熵就达到最大值。 概述概述热力学中的熵热力学中的熵 物理学上熵是指热能除以温度所得的商，标志物理学上熵是指热能除以温度所得的商，标志 热能转化为功的程度；