从经典谱估计到现代谱估计

调和分析可以追溯到古代对时间的研究。18世纪，伯努利、欧拉、拉格朗日研究了波动方程及其正弦解。19世纪早期，傅立叶证明了在有限时间段内定义的任何函数都可以用正弦和余弦分量的无穷次谐波之和来表示。1898年，Schuster在傅里叶分析的基础上拟合了待分析的信号，研究了太阳黑子数的周期变化，提出了周期图的概念。1930年，维纳发表了经典论文《广义调和分析》，对平稳随机过程的自相关函数和功率谱密度给出了精确的定义，并证明了它们之间的关系，从而为功率谱分析奠定了坚实的统计学基础。在1934中，辛钦还独立证明了自相关函数与功率谱之间的傅里叶变换关系，即维纳-辛钦定理。根据这个定理(详见第一章)，平稳离散随机信号x(n)的自相关函数rxx(m)为

rxx(m)= E[x(m+n)x *(n)](4-1)

和功率谱Pxx(ejω)形成傅立叶变换对，即

地球物理信息处理基础

如果x(n)是遍历的，它的自相关函数可以通过时间平均法从它的一个采样时间序列中得到，即

地球物理信息处理基础

在大多数应用中，X(n)是实信号，因此上述公式可以写成

地球物理信息处理基础

事实上，通常在时域中只能观察到随机信号的有限数量的采样值(例如，n个值),这可以表示为

xN(n)={x(0)，x(1)，…，x(N-1)}={x(n)，n=0，1，…，N-1}

它的自相关函数只能由这N个采样数据来估计，并且经常使用有偏估计。

地球物理信息处理基础

这是一种渐近一致的估计，称为抽样自相关函数。

采用采样自相关函数的傅里叶变换作为功率谱的估计。该方法由R.Blackman和J.Tukey在1958中提出，被称为功率谱估计的自相关法(简称BT法)。这种方法需要先得到有限观测数据的估计自相关函数，然后根据方程(4-2)计算功率谱。在快速傅立叶变换(FFT)算法被提出之前，它是最流行的功率谱估计方法。

1965年，Cooley和Tukey完善了著名的FFT算法，使傅里叶变换的计算速度提高了两个数量级，计算量显著减少，使DFT变换迅速广泛应用于各个领域，特别是工程实践中。根据公式(4-5)，是x(n)和x(-n)的卷积运算，因为

地球物理信息处理基础

如果x(n)的傅里叶变换是x(EJω)，那么x(-n)的傅里叶变换等于x *(EJω)。对公式(4-5)的两端进行傅里叶变换，得到

地球物理信息处理基础

这说明，通过直接对随机数据进行离散傅里叶变换，然后取其幅度的平方，再对多个样本进行这种运算，取其平均值作为功率谱的估计，即舒斯特周期图，这种谱估计由于不需要计算自相关函数，而是直接计算功率谱，所以引起了广泛的关注。

周期图和自相关法及其改进方法被称为功率谱估计的经典方法，周期图和自相关法是经典功率谱估计的两种基本方法。由于FFT的出现，周期图和自相关常常结合在一起，步骤如下:

(1)为xN(n)补全N个零，求；

(2)通过傅里叶变换，得到| m |≤m = n-1；

(3)对于加窗函数v(m)，此时| m |≤ m < < n-1，我们得到；

(4)傅立叶变换的使用，即

地球物理信息处理基础

周期图法得到的功率谱估计方差不随样本长度的增加而趋于零。令人惊讶的是，无论数据记录有多长，用周期图法和自相关法得到的估计都不是好的功率谱估计。其实随着录制长度的增加，这两个估计的随机波动会更加严重！此外，它们还有以下两个难以克服的固有缺点。

(1)频率分辨率(区分两个相邻频率分量的能力)不高。由于它们的频率分辨率(赫兹)与数据记录长度(秒)成反比(即δf = k/Tp = k/NT，k为常数，Tp=NT为数据记录长度，t为采样周期)，在实际应用中一般无法获得长的数据记录，即只能有限地观测数据，未观测的数据视为0。这样，如果只有n个观测数据，而信号仍然与n以外的数据强相关，那么估计的功率谱就会有很大偏差。

(2)对于有限的观测数据，相当于在时域将信号乘以矩形窗函数，因此相当于在频域将实功率谱与sinc函数进行卷积。因为sinc函数不同于δ函数，它有主瓣和旁瓣，这使得卷积功率谱不同于真实功率谱。sinc函数的主瓣不是无限窄的，导致功率谱扩展到附近的频域，造成频谱模糊，降低了频谱的分辨率。同时，由于sinc函数的旁瓣，能量漏入旁瓣(称为旁瓣泄漏)，即造成频谱间的干扰。强信号功率谱的旁瓣影响弱信号功率谱的检测。严重时会造成主旁瓣的极大失真，使微弱信号无法被检测到，或者将旁瓣误认为是信号，造成假信号。为了改进经典的功率谱估计，可以使用各种窗函数，但结果是以增加主瓣的宽度来换取旁瓣的降低，所以功率谱分辨率低是经典功率谱估计的致命缺点。

为了克服上述缺点，人们进行了长期的努力，提出了平均、加窗平滑等方法，在一定程度上改善了经典功率谱估计的性能。实践证明，基于傅立叶变换的经典功率谱估计方法对于长数据记录确实是实用的。然而，频率分辨率和功率谱估计稳定性之间的矛盾无法用经典方法从根本上解决，尤其是在数据记录较短的情况下。这促进了现代功率谱估计方法的发展。

现代功率谱估计方法主要基于随机过程的参数模型，称为参数模型法。虽然现代功率谱估计技术的研究和应用主要始于20世纪60年代，但实际上时间序列模型在非工程领域早已被采用。例如，Yule在1927使用自回归模型预测和描述经济时间序列的发展趋势，Walker在1931使用自回归模型，而Prony早在1795就使用指数模型进行拟合。在统计学和数值分析领域，人们也采用了模型方法。

现代功率谱估计主要是针对经典功率谱估计(周期图和自相关法)分辨率和方差性能差的问题。受地震学研究中线性预测滤波的启发，1967 Burg提出了最大熵谱估计方法，在提高分辨率方面做了最有意义的探索。Parzen在1968正式提出了自回归谱估计方法。1971年，Van der Bos证明了一维最大熵谱估计与自回归谱估计等价。1972中出现的谱估计的Prony方法类似于数学中的自回归方法。目前，基于自回归滑动平均模型的谱估计比自回归模型具有更高的频率分辨率和更好的性能。Pisarenko在1973中提出的谐波分解法提供了一种可靠的频率估计方法。1981年，施密特提出了谱估计的MUSIC算法。因此，现代功率谱分析主要包括四种方法:ARMA谱分析、最大似然、熵谱估计和特征分解。ARMA谱分析是一种建模方法，即通过平稳的线性信号过程建立模型来估计功率谱密度；熵谱估计包括最大熵谱和最小交叉法；特征分解又称为特征构造法和子空间法，包括Pisarenko调和分解法、Prony法、MUSIC法和ESPRIT法(利用旋转不变技术估计参数)。

现代功率谱估计研究仍以一维功率谱分析为主，大多基于二阶矩(相关函数、方差、功率谱密度)。然而，由于功率谱密度是频率的实函数，缺乏相位信息，基于高阶谱的谱估计方法正在引起人们的关注，尤其是双谱估计和三谱估计的研究受到了高度重视。其他研究如多维谱估计和多通道谱估计也在发展中。这些新方法有望在信息提取、相位估计和非线性描述方面得到更多的应用。