第九章 波谱分析

第九章波谱分析基本原理功率谱估计滤波方法第九章1谱分析的基本原理大气中的运动，无论从空间上，还是时间上，都存在不同尺度的波动现象。时间上的不同尺度的波动现象，体现为不同的“周期”（或“频率”）。任意一个气象要素的时间变化曲线，都可以看成是由各种频率的规则的波动叠加而成，如果我们对不同频率的波动的方差贡献进行比较，即可分析出时间序列的主要周期，这种从频率域上对气象要素时间序列进行分析的方法称为：“谱分析”。波谱的概念最早在物理学中，用光谱的概念来描述不同频率的各种单色光的性质；光谱：复色光经过色散系统（如棱镜、光栅）分光后，被色散开的单色光按波长（或频率）大小而依次排列的图案。傅立叶(Fourier)级数对于周期为T的任意一个时间函数f(t)，在满足狄利克雷条件的情况下，可以展开成多种频率的谐波相叠加的形式：ak，bk称为傅里叶系数，它们的计算方案是：这里的频率是离散的波数k越大，则该谐波的周期越短，频率(fk)或圆频率(ωk)也就越大。理解傅立叶变换幸福的生活，离不开傅里叶变换；广泛应用于音视频信号处理等领域。右图展示：傅里叶级数(蓝线=各灰线相加)逐渐逼近周期函数f(t)(红线)傅里叶变换本质上是把时域的信号转换为不同频率的各分量(谐波)的叠加各分量具有不同的振幅；通过傅氏变换可以分析各频率的振幅！这就是：振幅谱A(ω)ω因此，傅立叶变换通过对时域信号的分析，得到了频域的信息。各谐波分量除了具有各自的振幅以外，还具有各自的位相！如右图所示：位相谱综上，把“时域图像”、“频域的振幅”和“频域的位相”三者放在一起，如下立体图像：为了考察各个谐波的振幅和位相，利用复变函数中的欧拉公式，将以上傅里叶级数改写成：综上，我们得到了“频率为离散形式的”时间函数和复谱的相互关系式：所以，得到了Ck的表达式：非周期信号的傅里叶变换以上分析的是周期为T的时间序列，但大多数实际资料并没有周期，这就需要考虑“非周期信号”的傅里叶变换。“非周期信号”，可以看成是当周期T∞时的“周期信号”。当T∞时，相邻两个频率的间隔是无穷小，因此，可以认为是连续频率，这时得到的是连续谱。把傅里叶系数公式代入傅里叶级数的表达式：利用三角函数的“积化和差”公式，可导出：因此，有：利用欧拉公式，将上式中的余弦用复数表示：复谱的实部相当于离散谱中的傅里叶系数ak,虚部相当于bk。振幅谱A(ω)相当于离散谱的振幅谱Ak。位相谱θ(ω)相当于离散谱的位相谱θk。综上，我们最终得到了“频率为连续形式”的傅氏变换对：第九章2功率谱分析前面通过傅里叶变换，我们把任意时间序列分解为了各种频率(周期)的谐波，并得到每个谐波的振幅和位相；实际气象资料分析时，我们往往对“哪个周期(频率)的波动起主要作用”感兴趣，这就要分析比较各种谐波的方差贡献或振幅，一般称为“功率谱”，这来自物理概念：按照一般的物理概念，若电阻为1个单位，瞬时电压用f(t)表示，则它的瞬时功率为f2(t),它的总能量为：从统计观点看，上式表示数学期望为零的f(t)的总方差。可把这种总能量(或方差)分解为各种频率的振动分量的能量(或方差)之和。某一时间段内的总能量除以时间间隔得到的就是平均功率，把各频率振动分量的平均功率作成谱图，就是功率谱。离散功率谱设时间函数f(t)及g(t)具有相同的周期T，可分别表示为傅里叶级数的形式：将以上两等式两边的项相乘，取积分，并取T时段内的平均，可得：当g(t)=f(t)时，就可以得到时间函数f(t)的方差或平均功率：因此，时间函数f(t)的平均功率(或方差)可以用各谐波振幅的平方之和来表示。因此，不妨将谱值乘以2，只画出离散功率谱图的正半轴一部分(k>0),故令于是，可以通过功率谱图上各谐波振幅的相对大小，来考察时间序列的主要周期。各谐波的周期可以用：Tk=T/k

得出，T为时间函数f(t)所定义的基本周期：n。离散功率谱的估算及检验设一实测的气象时间序列为x1,x2,…,xn,根据前面的介绍，离散功率谱的计算步骤如下：根据奈奎斯特(Nyquist)采样定理，对于离散时间序列，可以分辨出的最小周期是“2倍的采样间隔”，因此最大波数是int(n/2)注意：在谱图上，为了方便分析，通常在波数坐标下面标上对应的周期或频率值，周期值与波数间的关系为：[例1](P237)根据青岛40年(1951-1990)平均气温资料(如下图所示)作离散功率谱估计。先将原始资料距平化，然后计算各波数k的傅氏系数，最大波数k=int(40/2)=20；例如，k=1时，算得的傅氏系数为：则第1谐波的功率谱估计为：以此类推，陆续算出波数k=2,

3,

…,int(n/2)时的功率谱值，画出离散功率谱线。根据各谐波的谱值的大小，分析出时间序列的主要周期。功率谱值Sk2越大，对应的周期Tk越显著。图中显示，谱值最大的是：k=7，其次是k=3，它们对应的周期Tk分别是40/7年，40/3年。但是问题来了：功率谱值大到何种程度时，可以认为该周期是显著的？需要对功率谱值进行显著性检验！课本P237、P238图9.1勘误离散功率谱的检验原假设：H0:E(ak)=E(bk)=0,对于该例，显著性水平为0.05时，临界F值Fc=3.25,根据各个谱值计算F值，大于Fc的谱值所对应的周期可认为是显著的。例如上图中功率谱的最大值为当波数为7时的S72=0.0356，

若要检验该周期是否显著，根据F统计量的公式算得F7=3.51，大于临界Fc=3.25，因此，该周期是显著的；可依此继续逐个判断其他周期是否显著。对于该例，最终的结果是:谱值最大的两个周期40/7和40/3年是显著的，其他周期不显著。(勘误教材P240)要检验某一频率(或波数k)的谱值Sk2=(ak2+bk2)/2所对应的周期Tk=n/k是否显著,换个角度，采用“临界谱值”可以更方便地检验可知，当谱值Sk2=(ak2+bk2)/2大于临界谱值Sc2时，F将大于临界Fc值，这时可认为Sk2所对应的周期Tk是显著的。在0.05的显著性水平下，把Fc=3.25，代入得Sc2=0.0333,因此，离散功率谱图上，谱值大于0.0333的周期都可认为是显著的。左图红色虚线表示Sc2，谱值高于该值的周期可认为是显著的。时间序列离散功率谱例2：时间序列及其离散功率谱离散功率谱离散功率谱例3：太平洋Niño3.4区(5°S-5°N,160°W-120°W)1950-2013年的月平均海温年际异常序列（64年共768个月）及其功率谱离散功率谱时间序列注意上图横坐标从左到右是按波数k从小到大来排列的，于是周期Tk就按从大到小排列。为方便观察主要周期，习惯上，常把横坐标按周期从小到大的顺序线性排列：根据上图的情况，若要更细致地观察主要周期，可以重新设置横轴的范围：可见，ElNiño具有2-7年的周期连续功率谱的估计及检验从连续频率的角度推导或估算各波动的方差大小，就是下面要介绍的“连续功率谱”。时间函数的总能量(或总方差)可以表示为：P236勘误习惯上，把|F(ω)|2称为能量谱密度或功率谱密度，上式表明，时间函数f(t)的总能量等于功率谱密度对所有频率的积分。功率谱密度简称“功率谱”对某一气象要素序列xt(t=1,2,…,n),连续功率谱的估计并非直接从定义式出发计算，而是通过时间序列x的自相关函数作间接估计。

设时间函数xt为标准化时间序列，它的数学期望为0，方差为1，则它的相关函数为：S(ω)为“功率谱密度”(PowerSpectrumDensity;PSD),表征时间序列的功率随频率的分布；上式表明对平稳过程，功率谱密度S(ω)和自相关函数ρ(τ)是一傅氏变换对(维纳-辛钦定理)。于是，“连续功率谱”是通过自相关函数ρ(ω)来估计的。只需计算不同时滞的自相关，然后对时滞τ积分即可。因此这种周期检测方法又称为“自相关函数”法。下面以太平洋Niño3.4区(5°S-5°N,160°W-120°W)1950-2013年的月平均海温年际异常序列为例，介绍计算步骤。（1），由以上公式知，需首先计算样本的自相关系数r(τ)，(τ=0,1,2,…,m),m为最大落后长度。关于“最大落后长度”m的选取：m太小时，用于谱估计的采样点太少，影响功率谱的估算精度；m太大时，用于计算自相关系数的样本长度就太小，影响自相关系数的估算精度。一般m取在n/3~n/10之间。本例，n=12*64=768，可取m=250;右图所示为Nino3.4序列的自相关系数r(τ)（2），利用“梯形法”对τ求积分，计算粗谱。“梯形法”求定积分示意图根据以上求积分的近似公式知，利用自相关系数r(τ)计算S(ωk)时，可采用以下有限求和公式（这里Δx=1,单位为采样间隔）：由于我们拿到的资料数据是离散点，可检测出的最短周期为2（单位为采样时间间隔），周期Tk的变化范围是从2到无穷大，频率fk变化范围是从0到1/2，因此

圆频率ωk的范围是从0到π。取样点数为m时，ωk=2πk/m，波数k的范围为：0~m/2，频率的间隔为1/m。习惯上，令l=2k,（这时,l的范围是0~m）则功率谱可写为：此处1/m表示作了个平均（3），把上一步计算的粗谱进行平滑，以减小误差。即：在两个端点处采取两点平滑（权重为[0.50.5]），在其余点采用权重为[0.250.50.25]的三点平滑.另外，为了保持权重的一致性，端点处需再乘以0.5（即相当于端点处以及邻近端点处的权重均为0.25）。最终，平滑功率谱密度的估算公式为：三点加权平滑(4)作谱图以波数l（或频率、周期）为横坐标轴，以平滑功率谱密度估计值Sl为纵坐标，作谱图。波数l与周期的关系为：通常把横轴按周期Tk从小到大的顺序进行线性排列：Nino3.4区海温距平的功率谱图大的功率谱值所对应的周期是否显著，需要显著性检验。连续功率谱的检验原假设：假设要检验的序列为非周期性随机过程。功率谱的检验过程就是，把要检验的时间序列的功率谱与非周期过程（红噪声或白噪声）的功率谱进行对比，对于某一周期，如果时间序列的谱值大于噪声的谱值，就表明该周期是显著的。红噪声的功率谱第八章已经学过，红噪声即一阶自回归模型，落后τ时刻和落后1时刻的自相关系数的关系为：非周期性随机过程，主要有两种：红噪声和白噪声。可见，红色噪音谱在ω=0时达到最大值，即只有线性趋势而无其他周期震荡。红噪声的功率谱白噪声的功率谱在介绍白噪声的功率谱之前，先回顾“白噪声”在数学上的定义：一个时间连续的随机过程X(t)为白噪声过程当且仅当它的均值函数μ(t)与协方差函数满足以下条件：检验假设所要检验的某一气象要素的总体谱为某一非周期性的随机过程谱，在这一假设成立的条件下，某一频率的谱估计值与所假设的过程的平均谱估计之比，遵从自由度为f的χ2分布，即：下面所剩的问题就是S0l的计算。这取决于我们要假设哪种非周期随机过程。关于S0l的计算，首先要确定非周期随机过程，通常是“红噪声”与“白噪声”二者中的一种，至于要选哪一个，需要考察样本的时滞为1的自相关系数：r1，如果r1接近于0或者是负值，则假设过程为“白噪声”过程，否则设为“红噪声”过程。如何衡量r1是否接近于0？可用‘相关系数的显著性检验’来决定。（1）如果我们最终确定要假设的随机噪声为红噪声，结合红噪声的功率谱，

红噪声的谱值S0l由下式算出：（2）如果我们最终确定要假设随机噪声为白噪声，白噪声的功率谱对于Nino3.4指数功率谱的检验，由于r(1)远大于0，接近于1，所以需要假定随机噪声为红噪声，作出临界谱值如下图红线所示：另例：阅读教材P244[例如]，分别假设随机噪声为红噪声和白噪声，对青岛年均温度序列的功率谱进行检验。第九章3滤波方法前面已经学过，任意时间序列信号可看成是由不同频率(周期)的信号叠加而成的，我们在对时间序列进行分析时，往往根据需要，先滤掉一些不想要的周期信号，以突出我们需要的周期信号，这一过程称为：“滤波”。设一时间函数f(t),经过一个“滤波器”输出为新的时间函数g(t)；“滤波器”是指的一类数学运算的总称。“滤波”就是对原始数据序列进行频率过滤，得到具有目标频率范围的另外一条新序列。设用一脉冲函数δ(t)作为输入，经过滤波器系统后，它的输出记为h(t),

称为脉冲响应。若用L表示线性系统由δ(t)到h(t)的转换，可记作：h(t)=L[δ(t)]输入f(t)输出g(t)滤波器对任意时间函数f(t),根据脉冲函数的筛选性质，可表示成无穷个脉冲的积分形式：类似地，系统的输出g(t)也可以表示为无穷个脉冲响应的积分形式，即：可见，该线性系统由“脉冲响应函数”h(t)来决定。把脉冲响应函数h(t)的谱记为H(ω)：若记F(ω)和G(ω)分别为输入函数f(t)和输出函数g(t)的谱，对于输出谱G(ω)，则有：对于输出函数g(t)的功率谱，功率谱密度为：Sg(ω)=|G(ω)|2=|H(ω)|2|F(ω)|2=|H(ω)|2Sf(ω)上式表明：对于某一频率ω的振动，经过滤波后，

振动方差变为原来的|H(ω)|2倍；若|H(ω)|<1，表明该频率的振动方差有所削弱。若|H(ω)|=0，表明该频率的振动完全被削除。一般地，根据滤波器所要滤掉的频段，滤波可分为：“低通滤波”、“高通滤波”和“带通滤波”三种。下图所示的三种响应函数曲线可以表征这三种滤波。顾名思义，低通(low-pass)滤波：使过滤后的序列主要包含低频振动分量，把高频的分量滤掉，响应函数如图中实线(‘—’)所示；高通(high-pass)滤波：使过滤后的序列主要包含高频振动分量，把低频的分量滤掉，响应函数如图中点线(‘…’)所示；带通(band-pass)滤波：把高频和低频的信号都滤掉，只保留某一段频带的信号，响应函数如图中虚线(‘---’)所示。1.低通滤波器(low-passfilter)“滑动平均”滤波法。对时间序列选取某一尺寸为m（m为奇数，m=2k+1）的窗口，将窗口内所有的异常值做算数平均，将所求出的平均值作为窗口中心点处的新值。从时间序列的始端逐点向末端移动窗口，不断重复以上求平均的过程，直到到达时间序列的末端为止。例如，m=5,要对长度为n的序列x进行“5点滑动”滤波，即窗口的长度为5，得到的滤波后的新序列为y,y(3)=[x(1)+x(2)+x(3)+x(4)+x(5)]/5y(4)=[x(2)+x(3)+x(4)+x(5)+x(6)]/5……y(t)=[x(t-2)+x(t-1)+x(t)+x(t+1)+x(t+2)]/5……y(n-2)=[x(n-4)+x(n-3)+x(n-2)+x(n-1)+x(n)]/5可见，滑动平均的优点是计算简单，但会造成时间序列首尾的数据缺失，使得平滑后的数据长度小于原数据长度。下面导出滑动平均的响应函数。以便在频域中考察该滤波器使不同频率的振动分量的振幅产生了怎样的变化。正弦函数是奇函数，积分为0对于“等权重”

的m点滑动，滑动窗口的长度m=2k+1,权重为：欲对其计算，首先需要知道权重系数hj（j=0,1,2,..,k）的值。讨论：对于周期等于滑动间隔(m)的振动，由于频率为f=1/m,因此，H(f)=0,表示这种振动分量达到完全削除。同样地，对于周期等于m/i的高频振动（频率为i/m，i=1,2,3,…），也有H(f)=0，振动达到完全削除。对于周期大于m的低频振动，频率f<1/m，因此响应H(f)<1,这类周期也有不同程度的削弱，周期越大，削弱程度越小；对于无限长的周期(f—>0)，频率响应H(f)—>1，因此，过滤后无任何削弱。H(f)|H(f)|2例：对月平均的北太平洋（NP）海温异常（1950-2013）EOF第一模态的时间序列(PC1)作低通滤波红线表示13点滑动，因此主要把周期小于1年的信号过滤。绿线表示121点互动，主要把周期小于10年的信号过滤掉，突出年代际信号。前面所讨论的滑动平均都是“等权”滑动，例如对于3点滑动，每点的权重hi(i=-1,0,1)都等于1/3，对于5点滑动，每点的权重都是1/5。。。除了“等权滑动”之外，还可以使用“不等权滑动”，把权重系数设为与二项式系数成正比，称为“二项系数滑动”。因此，m点(m=2k+1)滑动,对应于(a+b)m-1（即(a+b)2k）的二项式系数，滑动系数为：教材262勘误例如，当m=3时，k=1,二项式(a+b)2的系数为：1,2,1，故权重系数为:1/4,1/2,1/4;当m=5时，k=2,二项式(a+b)4的系数为：1,4,6,4,1，故权重系数为:1/16,4/16,6/16,4/16,1/16。二项式系数滑动的频率响应函数为：等权滑动H(f)右图为二项式系数的频率响应函数(上图)，并与等权滑动的响应函数(下图)作对比。可见，二项式系数滑动与等权滑动相比，是更有效的低通滤波器：能够通过更多的低频信号，同时也能滤掉更多的高频信号。二项式系数滑动H(f)2.高通滤波器(high-passfilter)高通滤波是指：使过滤后的振动分量主要保留高频振动分量的滤波。常用的高通滤波方法是：差分过滤“一阶差分过滤”可以表示为：但是由于系数2q的存在，对于很高频的振动，H(f)将大于1，使得系统不稳定。对于一阶(q=1)差分滤波：f=1/6时（周期为6倍的采样间隔），H(f)=1，该频率的振动分量没有任何削弱;对于频率f>1/6的振动，由于H(f)>1，振幅将得到增强；对于二阶差分滤波，这一分界点为f=0.08，对应的周期为12.4倍的采样间隔。3.带通过滤(band-passfilter)如果我们并非只保留低频振动，也非只想保留高频振动，而是想滤除某一段感兴趣的“频率带”的振动时，可采用“带通滤波器”。（1）方案一，采用两个简单的低通滤波器相减，实现带通滤波；例如，对原序列x先做5年滑动平均，削去周期为5年以下的波，剩下周期为5年以上的波，得到新序列记为x1，再把原序列x做9年滑动平均，除去周期为9年以下的波，剩下周期为9年以上的波，得到的新序列记为x2，令x3=x1-x2，得到的新序列x3包含了5-9年周期的振动，相当于实现了带通滤波。这种滤波器的响应函数为：H3(f)=H1(f)–H2(f)对于3点和21点二项式系数滑动，他们各自的响应函数以及二者之差如右图所示：绿色实线表征带通滤波的响应函数。这种办法存在的问题是：对于我们想保留的频段，也存在一定程度的削弱。（2）方案二：利用傅里叶逆变换进行带通滤波；例如，我们想要对某个序列进行滤波，只保留频率在[f1,f2]之间的振动，这时，可把响应函数设计为：还可以根据其他需求进行任意设计，如教材P264，如果想滤出频率f位于[f0/2，2f0]之内的序列，并且希望在频率f0上无任何削弱，那么就将滤波器的响应函数设为：脉冲响应h(t)可以通过对频率响应函数H(f)实施傅里叶逆变换来得到：可将上式离散化进行求解，f的取值范围为：[0,1/2]。f的离散采