频带集成谱本征正交分解与频率归因

《Physical Review Fluids》：Band-ensemble spectral proper orthogonal decomposition with frequency attribution

　　本研究提出了一种频带系综谱本征正交分解（bSPOD）方法。该方法受频率平滑技术的启发——频率平滑是一种用于降低功率谱密度估计中估计方差的方法——并将其推广至SPOD。该算法从全时记录的单次傅里叶变换所获得的连续傅里叶系数中估计SPOD模态，从而避免了时间分段。本研究中，bSPOD被应用于人工测试数据以及宽频-音调空腔流动的粒子图像测速数据集。与更常用的基于Welch方法的SPOD公式相比，bSPOD降低了频谱泄漏，允许更高的频率分辨率，并在相当的计算成本下保留了音调分量的频率信息。这些特性使得估计方差降低，同时保持音调分量的低偏差，使bSPOD对宽频-音调流动尤为有效。

I. 引言

湍流流场本质上是高维且混沌的。在混沌动力学中识别有组织的时空模式，长期以来一直是流体力学的核心问题[1-6]。这些所谓的相干结构通常与高能脉动相关，并且往往是不利现象的根源，例如声产生[7,8]、翼型上的失速胞形成[9]、燃烧不稳定性[10,11]以及结构振动[12]。

处理此类混沌数据的技术在过去一个世纪中不断演变。在此背景下，功率谱密度（PSD）是一个核心工具，它描述了信号功率在时间频率上的分布。在流体力学中，PSD被广泛用于处理速度探针和压力测压孔信号。PSD是自相关函数的傅里叶变换，但通常通过周期图进行估计[13,14]。由于湍流是混沌的，因此使用平均来降低估计方差。对于PSD的估计值，该方差可表示为
（1）
并量化了有限样本谱估计与其理论值之间的偏差。由于后者通常先验未知，估计方差通常不能直接评估。然而，估计方差随着平均过程中使用的统计独立实现数量的增加而减小。在实践中，高估计方差表现为噪声大、抖动或不规则的谱估计。因此，术语“估计方差”常被更宽松地用于描述这些统计涨落，并在本工作中以此含义使用。在平均方法中，Bartlett[15]的方法最为常用：假设遍历性，该方法将记录分成多个段，计算每段的PSD，然后对结果进行平均。这降低了估计方差，但代价是频率分辨率降低，而频率分辨率由段长决定。Welch[16]通过引入重叠段并对每段施加窗函数以降低频谱泄漏，扩展了这一方法。

Thomson[17]引入了多窗估计器，它将全时记录乘以一组正交窗函数以获得近似独立的实现。这些加窗信号经傅里叶变换后在各窗函数间进行平均。这产生与全记录DFT相同频率网格上的PSD，但有效分辨率较低，因为估计值等于真实频谱与由窗函数傅里叶变换形成的谱窗的卷积[18]。在匹配有效分辨率下，多窗方法通常降低频谱泄漏，并在标准高斯假设下提供与分段估计相当或更低的方差[19]。其代价是计算成本，因为每个窗函数都需要计算一次全记录变换。

一种相关的方差降低方法是频率平滑[20-22]。Daniell[20]首先提出对相邻频率的周期图值进行平均，Bartlett[21]指出了其与截断自相关的联系。Blackman和Tukey[22]随后通过在大时延处对自相关进行加窗，正式确立了滞后窗自相关估计器。当平滑滤波器和滞后窗构成傅里叶对时，频率平滑和滞后窗两种视角是等价的。简单来说，该方法通过对相邻频率上的全记录周期图进行平均来降低估计方差，代价是有效频率分辨率降低。

时间分段、多窗法和频率平滑都能降低估计方差，但会增加偏差，即期望谱估计与真实谱之间的偏差，因为无论时域还是频域的有限窗操作都会导致功率扩散到相邻频率单元。在PSD变化剧烈的区域，例如频谱中的音调峰处，更强的平均会使窄谱特征变宽。三种方法之间的一个关键区别是：时间分段首先将记录分成较短的段，计算每段的周期图，然后进行平均；有效频率分辨率由段长决定，并且在整个频谱上基本均匀。相比之下，多窗法允许窗函数的数量和权重随频率变化，这些共同决定有效频率分辨率。类似地，频率平滑先计算单个全记录周期图，然后在频率上平均，因此平滑滤波器及其带宽可以逐频率选择，实现局部的偏差-方差权衡。

在流体力学中，不仅时间分辨数据，空间分辨数据也常常可用，例如来自粒子图像测速（PIV）或同时采样的压力或速度探针阵列。因此，除了时间统计量外，还可以评估空间统计量以提取相干结构。空间结构可以通过本征正交分解（POD）来识别，该方法通过两点空间协方差矩阵的特征值分解，得到能够最优捕获流动能量的正交空间模态。该方法遵循Karhunen-Loève框架，也称为主成分分析[23,24]；由Lumley[1]引入流体力学，并已被广泛使用和扩展。虽然POD不限于空间相干性，但最初主要用于提取空间结构[25,26]；为了将这种最常见的POD变体与其他形式区分开来，通常使用“纯空间POD”这一术语[6]。

Sirovich[2]引入了快照方法，通过求解时间相关（快照）矩阵的特征问题来降低POD的计算成本，该矩阵与空间协方差矩阵共享非零特征值。Borée[27]提出了扩展POD，用于在第二个变量中识别与参考测量最相关的模式。Sieber等人[28]应用了时间滤波POD，产生频谱集中在特定频带内的模态。Mendez等人[29]将多分辨率（小波）分析与POD结合；所得的多尺度POD可跨尺度分离结构，且不限于统计平稳过程。其自然扩展是时空POD，它识别完全空间和时间分辨信号的能量最优基，无需先验尺度限制[1,30,31]。然而，这需要对时空协方差进行特征值分解，计算成本很高。最近还引入了引入时间条件以预测极端事件[32]和生成预报[33]的修正版本。

许多相关流动是统计平稳的，因此时空POD并非必要，能量最优的空间模态可以作为频率的函数来识别[1,34]。这些模态通过对交叉谱密度（CSD）矩阵进行特征分解获得，该方法由Picard和Delville[35]首次称为谱本征正交分解（SPOD）。SPOD在每个频率处识别能量最优基[1]。其应用依赖于大数据集以估计跨频率的CSD。实践中，使用Welch分段：假设记录是遍历的，将其分成块，计算每块的傅里叶模态，并从块实现系综中估计CSD[36]。SPOD被广泛应用[35,37-40]，也因其与可解算子分析[6]的联系而受到关注，后者为数据驱动模态提供了基于物理的模型。SPOD的变体和扩展也被提出。Colanera等人[41]引入了一种鲁棒的POD公式以减轻对异常值的敏感性，改善了噪声数据集上的模态收敛性。Blanco等人[42]应用时间平移来对齐平流结构，从而改善SPOD模态收敛性。Heidt和Colonius[43]将SPOD扩展到循环平稳流动，其中宽带统计量是周期性的而非平稳的。

尽管SPOD在识别时间和空间相干结构方面取得了成功，但采用Welch平均的SPOD继承了有限记录谱估计的估计方差-偏差权衡[6,44]。提高频率分辨率（更长的段）改善了对尖锐谱特征的频率和功率估计精度，但在宽带区域产生高方差和慢收敛。增加块数（更短的段）降低方差并改善收敛，但降低了频率分辨率并使尖锐特征附近的估计产生偏差。为改善这一权衡，Schmidt[45]引入了多窗SPOD，利用了Thomson[17]的多窗估计器。如上所述，多窗法仍然涉及偏差-方差权衡，但在匹配有效分辨率下，其PSD估计性能优于Welch分段[19]。为降低多窗法带来的显著计算成本增加，Schmidt[45]提出了一种混合SPOD，将Welch分段与多窗法结合。在此基础上，Yeung和Schmidt[44]开发了一种受Riedel和Sidorenko[46]启发的自适应SPOD公式，基于正弦窗函数。加窗傅里叶模态估计从信号的单次傅里叶变换获得，使用零填充至两倍记录长度，从而大幅降低成本。遵循随频率自适应估计器的思想[46,47]，窗函数数量通过对齐-收敛准则选择：宽带区域使用更多窗函数以确保收敛，而在频谱峰附近，音调能量导致快速收敛，使用较少窗函数以保持低谱偏差。通过这种方式，自适应SPOD方法改善了谱估计的收敛性，降低了宽带区域的方差，并保持了音调分量附近的低偏差。然而，该方法是迭代的，在收敛前可能需要多次特征值分解，导致高计算成本。

在本研究中，我们提出了一种基于频带内傅里叶模态系综的SPOD模态计算算法，称为频带系综SPOD（bSPOD）。其思想受频率平滑启发，后者通过对相邻频率上的全记录周期图进行平均来降低PSD估计方差。对于独立的PSD估计，Welch分段比频率平滑更常用。然而，当与频带系综策略中的POD结合时，该方法充分发挥了其潜力。对常用SPOD算法和Welch方法的修改很小，但优势显著：bSPOD降低了频谱泄漏，并为每个模态提供数据驱动的频率估计，从而在保持音调特征值附近低偏差的同时实现宽带区域的方差降低，无需依赖自适应过程。计算成本与Welch分段的SPOD相当。由于平滑在计算频谱后应用，bSPOD允许在不重新计算傅里叶变换的情况下修改有效频率分辨率。

接下来，我们首先介绍标准的基于Welch的SPOD，然后描述bSPOD。这两种方法在人工数据集上进行比较以突出bSPOD的优势，并在宽频-音调空腔流动的实验PIV测量中展示其对实际数据的适用性。结果和特征被讨论，研究以主要发现的总结结束。

II. 方法

在本节中，我们概述基于Welch的SPOD算法及其频带系综自适应。我们专注于时间离散公式，并请读者参考文献了解进一步细节[6,34,36,44]。

本研究的一个关键焦点是宽频-音调流动构型。SPOD在平稳性假设下估计CSD。真实谱线对应无限时间相干性，因此相对于宽带分量本质上是奇异的。在有限数据中，线谱、窄带随机过程和循环平稳动力学之间的区别无法通过傅里叶变换识别。因此，本方法应操作性地解释为对有限宽频-音调数据集的分解，而不是作为潜在过程允许严格平稳遍历谱表示的证据。

A. 谱本征正交分解

我们考虑一个离散时空信号，具有个空间自由度和个时间样本。这里，我们使用空间和时间作为坐标，因为它们是流体动力学中最常见的。多个空间方向和物理变量可以堆叠到索引中。时间序列被分成个长度为的块。为了增加块数，相邻段可以重叠。在下文中，我们关注无重叠的核心算法，并在第IV.E节中单独考虑重叠。每个数据块通过时间离散傅里叶变换转换到频率空间， (2)其中 表示虚数单位， 表示数据块索引， 表示频率索引， 为用于减小频谱泄漏的锥形窗（窗函数）。带帽符号表示复数傅里叶振幅。频率间隔由窗长决定，即 ，其中 表示两个连续时间样本之间的时间步长。对于每个频率区间，不同数据块的傅里叶系数可以排列成数据矩阵 (3)其中每一列 是频率 处与数据块 相关联的空间分辨傅里叶系数。 的形状为 。通过对所有数据块取平均，得到频率 处的互谱密度矩阵估计 (4)其中 表示共轭转置。前置因子为 (5)它对窗能量进行归一化，并确保SPOD特征值具有功率谱密度的单位（每单位频率的方差）。 的定义与式(2)中的正向变换归一化一致，并确保满足帕塞瓦尔定理[44]。对 进行特征分解得到SPOD模态 (6)其中 是定义离散内积的权重矩阵。 的列中所含的特征向量即为SPOD模态，它们关于该内积正交归一，即 。模态的谱功率密度包含在对角矩阵 中。如果空间自由度的数量大于数据块数量，通常需要求解另一个具有相同非零特征值的特征值问题，称为快照法[2] (7)其中 是展开系数，在纯空间POD中即时间系数。该特征值问题的大小为 。SPOD模态可由下式得到 (8)该表达式表明，SPOD中的展开系数 量化了 中每个数据块傅里叶系数对某个SPOD模态 的贡献大小。总体而言，这产生 个SPOD模态，可根据特征值中包含的功率密度进行排序。所有模态都对应各自的傅里叶区间频率 ，其中 (9)与式(2)中的定义相比，频率索引发生了平移，使得零频率位于中心。由于窗长 由时间序列的总时长 和数据块数量 决定，SPOD面临频率分辨率与统计收敛性之间的权衡。提高频率分辨率会减少傅里叶系数样本数，从而导致特征值收敛性变差，表现为估计的模态功率谱密度方差增大。相反，使用更短的窗口来增加样本数可改善统计收敛性，但会降低频率分辨率并增加频谱泄漏。这导致频谱锐度损失，通常称为频谱偏差，当不同的音调成分因泄漏而模糊到相邻频率区间时尤为显著。这种频谱方差与频谱偏差之间的权衡是谱估计中的已知问题。在考虑这一权衡时，务必注意SPOD模态构成了每个频率区间动力学的最优能量基。由于有限分段长度，每个区间通过频谱泄漏接收到来自一定频率范围的贡献，其中最强贡献来自邻近的未分辨频谱成分。一种量化该主要贡献频带有效宽度的方法是等效噪声带宽，对于矩形窗，其等于 [48]。因此，考虑更短分段所获得的额外样本仅当在加宽的有效带宽 内贡献的流动动力学保持足够相似时，才能改善谱估计的收敛性。这要求流动动力学不随频率剧烈变化。我们注意到，锥形化通过降低谱窗函数旁瓣的幅度来减小频谱泄漏，但代价是增大了主瓣宽度，进而增加了有效带宽[48]。简言之，锥形化减小了远距离泄漏，但增加了局部频谱模糊。 B. 频带集成谱正交分解频带集成SPOD基于谱估计中使用的频率平滑方法，并扩展到时空数据。与基于Welch的算法相比，bSPOD从一个决定性的差异开始：数据不分割成时间块，而是对整个时间序列进行离散傅里叶变换， (10)为了区分符号与基于Welch的方法，我们使用波浪线。因此，傅里叶变换的频率间隔 相对于基于Welch的SPOD按总长与块长之比增大，即 。在当前表述中不施加锥形化。虽然锥形化并非根本性地与该方法不兼容，但bSPOD依赖于全长度记录的单次傅里叶变换，其中更长的记录长度会为隐含假设的矩形窗产生更窄的谱窗，从而自然地减少频谱泄漏，相比基于Welch的SPOD中使用的更短分段。锥形化可以进一步抑制谱窗的旁瓣，但代价是加宽其主瓣，从而增加相邻傅里叶系数之间的相关性，如前面所述。与Welch方法不同，数据矩阵不是从固定频率的重复实现构建的，而是由 个连续傅里叶模式构建， (11)其中 表示频率 处的全长傅里叶模式。为了说明与基于Welch的表述的联系，我们考虑Welch表述中数据块数量等于频带集成形式中考虑的连续傅里叶模式数量，即 ，且不使用重叠的情况。在这种情况下，Welch-SPOD矩阵 包含在间隔为 的频率网格上定义的傅里叶模式的 个实现，而频带集成矩阵 包含在增强分辨率 下的 个连续傅里叶模式。对应的互谱密度矩阵的模拟可近似为 (12)其中 针对傅里叶变换步骤中提高的频率分辨率进行了调整。按照这种前置因子定义，选择 意味着两种方法都以相同的频率间隔 考虑动力学，并使用相同的每个估计的样本数。这可以通过比较式(4)和(12)看出：在不重叠且 的情况下，Welch-SPOD中的CSD估计可以解释为 个实现的平均，每个定义在有效带宽 下[见式(4)]；假设不施加锥形化（ ），前置因子变为 。在bSPOD中，在相同条件下，CSD定义[式(12)]中的前置因子 简化为 ，使得所得矩阵表示具有相同有效频率带宽的互谱密度。剩余因子 在两种方法之间不同，因为bSPOD不平均 个独立实现，每个定义在频率间隔 下，如同SPOD中采用Welch分割那样。相反，它聚合了频带 内 个连续傅里叶模式的贡献，每个定义在 倍更精细的分辨率 。因此，前置因子确保两个CSD矩阵在有效带宽 内表示相同的总功率。同时，以更高频率分辨率从连续傅里叶模式估计CSD是有利的，因为通过展开系数可以恢复带内信息，如下文所示。由于展开系数的重要性，bSPOD依赖于应用快照法 (13)其大小为 。bSPOD模式如下 (14)使得 。类似于基于Welch的SPOD中的快照法，式(14)中的展开系数量化了 中每个元素对 中每个bSPOD模式的贡献。由于在bSPOD中， 中的每个元素是定义在窄且唯一的频率区间上的傅里叶模式，该信息可用于获得每个模式的带内频率信息。具体而言，我们通过使用展开系数作为频率权重，为每个bSPOD模式分配一个具有代表性的、数据驱动的频率。 的列包含每个离散傅里叶模式对给定bSPOD模式的贡献；因此我们定义频率权重为 (15)其中 是 的元素， 是控制加权类型的指数。排名靠前的 个bSPOD模式 的数据驱动频率估计通过权重 与 中包含的离散傅里叶模式的频率 之间的矩阵向量积获得[见式(11)]， (16)选择 ，权重变为 ，每个模式被归因于窗口中心频率。当 ，平均值按傅里叶模式对bSPOD模式的相对贡献进行加权。当 ，每个bSPOD模式被分配其单个最大贡献傅里叶模式的频率。总体而言，数据驱动的频率归因允许每个bSPOD模式，尽管定义在频率带上，被关联到更精确的带内频率估计。这通过从全时间记录的单次傅里叶变换获得的相邻傅里叶模式形成CSD矩阵成为可能。关键是展开系数，当通过连续傅里叶模式的POD计算最优基时，这些系数允许保留高分辨率傅里叶模式的带内频率信息。我们强调，基于Welch的SPOD和bSPOD都试图近似相同的底层CSD算子。因此，只要CSD近似充分收敛且底层假设成立，所得模态就保留SPOD相关的理论性质，包括与可解分析的既定联系[6]。bSPOD的开源MATLAB实现可在网上获取[49]。 C. 频带集成滤波器长度影响bSPOD的主要参数是频窗长度 ，如上所述，它在Welch-SPOD中的作用类似于数据块数量 。对于小 ，bSPOD模式仅从少数几个连续傅里叶模式计算，这提高了频率分辨率，但降低了统计收敛性，因为特征值问题较小，所考虑频带内的动力学必须由少量模式表示。增加 会加入连续傅里叶模式：模式总体上必须表示更宽频带内的动力学，但更大的问题也提供了更多在模式间分配功率的灵活性。如果频率内容随频率变化缓慢，随着更多傅里叶模式加入，收敛性会改善。特别有趣的是极限 和 。在极限情况 ，数据矩阵 简化为单列， (17)因此相应的特征值问题变得平凡，bSPOD模式与归一化离散傅里叶模式本身重合，其相关特征值为其幅值平方 。在相反的极限 ，所有傅里叶模式都被包含，模式简化为纯空间POD模式。这可以通过定义去均值的时空数据矩阵 (18)其中 为空间自由度， 为时间自由度，以及酉DFT矩阵 (19)来证明。这使我们能够将bSPOD数据矩阵写为 。因此，根据式(12)，bSPOD互谱密度估计变为 (20)其中前置因子与式(12)中的定义略有不同，因为傅里叶变换基于酉归一化。用时间步长 归一化后，式(20)成为时间平均协方差矩阵。因此，bSPOD特征值问题简化为纯空间POD协方差问题。由此我们得出结论，参数 跨越了从 （模式简化为DFT模式，仅捕捉时间相干性）到 （模式成为纯空间POD模式，最优表示整个时间序列的空间二阶统计量）的过渡。 的中间值产生bSPOD模式，这些模式最优表示局部频带内的空间二阶统计量，从而表示时间和空间相干结构。对于基于Welch的SPOD，同样达到这些极限：使用单个数据块（不分割），CSD估计的秩为1，模式等于归一化DFT模式；使用 个长度为1的数据块，CSD矩阵变为协方差矩阵，因为单个时间样本的DFT就是样本本身，SPOD模式变为POD模式。因此，在这些参数极限下，bSPOD收敛到与其他SPOD算法相同的端点，包括基于Welch的SPOD和Sieber等人的时间滤波POD[28]。bSPOD的一个关键优势是频带集成滤波器长度 不必在整个谱上恒定，而是可以随频率变化。实际上，在空间波长预期会发生快速变化的模态区域，可以选择较短的频带；而在模态变化较慢的区域，可以使用较长的频带。这在低频情况下尤其重要，因为相对频率变化自然更大。由于空间波长通常通过色散关系与时间频率耦合，这可能导致相干结构空间范围发生显著变化。例如，对于湍流射流，基于Welch方法的SPOD可能受到这种效应的影响，因此推荐使用可变频率箱宽度[50]。bSPOD的优势在于，可以在不重新计算傅里叶变换的情况下，在频谱中局部调整频带集合的滤波长度，从而使得有效的bSPOD频率分辨率能够随频率变化。相比之下，基于Welch方法的SPOD算法在频率分辨率（块长度）改变时需要重新计算傅里叶系数，因此修改频率分辨率需要完全重新处理数据。D. 差异总结作为方法部分的最后一部分，我们总结了基于Welch方法的SPOD和bSPOD方法之间的差异。两种情况下数据矩阵形成的示意图如图1所示。图1. 基于Welch方法和频带集合SPOD的示意图比较。基于Welch方法的算法（左侧）首先将时间序列分割成块，然后通过离散傅里叶变换将每个块变换到傅里叶空间（标记为DFT的黑色弯曲箭头）。这得到每个频率的傅里叶模态实现。然后，将给定频率的傅里叶模态组装成数据矩阵，从中估计CSD矩阵并计算SPOD模态[参见公式(4)–(8)]。频率箱的数量由总信号长度和块数决定。CSD估计和由此得到的SPOD模态定义在相应的傅里叶箱频率上，并且具有等效噪声带宽，假设使用矩形窗。bSPOD（右侧）反转了顺序：在完整记录上计算一次傅里叶变换，得到在细密频率网格上的傅里叶模态。与无重叠的基于Welch方法相比，间隔小了倍。然后由相邻的傅里叶模态形成块矩阵。CSD矩阵通过近似求和这些狭窄间隔模态的贡献，并求解特征值问题以获得bSPOD模态[参见公式(12)–(14)]。类似于基于Welch的方法，bSPOD频率箱的数量由总时间长度和滤波器长度决定。由连续傅里叶模态集合所跨越的CSD估计和相应SPOD模态的有效带宽为。如果在bSPOD中使用的连续傅里叶模态数量等于基于Welch方法SPOD中的块数，则得到的CSD矩阵具有相同大小，并以相同的有效频率带宽表示交叉谱密度，。在通常的SPOD假设下，即与所包含的频谱内容相关的流动动力学保持足够相似，两种近似应产生相同的结果。关键区别在于，基于Welch方法的SPOD依赖时间分割来估计CSD矩阵，而bSPOD则从定义在更细频率网格上的连续傅里叶系数的贡献来估计模态，该网格是通过对全长记录进行单次傅里叶变换获得的。这一差异产生了若干理论和实际影响。首先，用于估计CSD的不同方法导致对SPOD展开系数的不同解释。在Welch-SPOD中，系数量化了单个时间块对模态的贡献。在bSPOD中，系数则量化了频带内相邻傅里叶系数的贡献。这种频率分辨表示使得上述频率归因过程成为可能，并提供了关于频带内模态功率谱分布的额外信息。其次，基于Welch方法的SPOD采用块分割，并假设块之间统计独立。相比之下，bSPOD不引入这一假设，而是依赖完整时间信号的傅里叶变换。因此，bSPOD依赖于一个单一的长时间观测周期，在此期间统计量保持足够平稳，并且不直接适用于组合同一过程的多个较短记录。在无限观测时间极限下，不同频率的傅里叶系数理论上不相关，从而为CSD矩阵的近似提供独立贡献。虽然有限记录长度通过频谱泄漏在相邻傅里叶系数之间引入相关性，但这些效应随着记录长度的增加而减小。因此，对于足够长的时间序列，相邻傅里叶系数提供近似独立的实现。对于，两种方法都定义了与相同有效带宽相关联的CSD矩阵。然而，在直接比较中，应注意Welch-SPOD将模态估计与傅里叶箱频率相关联，公式(9)，而bSPOD将每个估计定义为由连续傅里叶模态跨越的区间。如果没有重叠且从第一个傅里叶模态开始集合，则导致bSPOD频率网格相对于Welch-SPOD网格偏移约，后者仅隐式表示有效带宽。然而，在实践中，这一区别并不重要，因为两种方法的目标都是估计模态功率谱密度随频率的变化。此外，bSPOD算法设计带来略微不同的计算成本特征。一次长度的傅里叶变换比次长度的傅里叶变换更昂贵。然而，bSPOD避免了与分割、锥形化和块循环相关的开销成本。额外的计算成本可能因为bSPOD依赖快照方法而产生，当空间自由度数量相对于频率窗口长度较小时，该方法更昂贵，并且因为频率归因所需的额外操作。然而，总体计算成本与基于Welch方法的SPOD相当，因为特征值问题通常是最昂贵的步骤，在可比配置下具有相同大小（）。作为回报，bSPOD提供了明显的优势：在完整时间序列上进行傅里叶变换减少了频谱泄漏，并且频率归因步骤使得能够改进对音调成分的频率估计，如第IV节所示。最后，有效的bSPOD频率分辨率可以在不重新计算傅里叶变换的情况下变化，使得该方法适用于自适应频率分辨率策略[44]。III. 数据库 A. 人工宽带-音调信号为了评估bSPOD的性能，我们构造了一个统计平稳且包含宽带和音调特征的人工时空信号。目标是模拟湍流的典型特征，同时保持对频谱内容的完全控制。信号定义在长度为的周期一维域上，在等间距点采样。时间离散为个样本，时间步长为。相应的网格为（21）其中和。信号的宽带部分通过叠加两个独立信号，得到，每个信号代表一个具有随机相位和幅度的对流波包。两个信号的目标时间频谱选为具有自然频率和阻尼的单个振荡器的响应：（22）其中表示角频率。信号在频率空间中构造，然后逆变换到物理域。在每个离散频率处，我们绘制独立的复高斯系数，（23）其中是独立的标准正态随机变量，是频率箱宽度。这种构造确保，使得期望功率谱与目标匹配，而随机相位使实现统计平稳。与的缩放确保信号的方差等于规定功率谱密度的积分。因此，除了随机实现外，两个宽带信号共享相同的频谱分布，而差异主要来自它们的空间结构，接下来考虑这一点。为了引入空间相干性，每个频率分量通过线性色散关系分配一个波数（24）以速度建模对流，偏移量为，这两个宽带分量选择不同。为了避免引入优选的空间原点，两个宽带分支通过随机偏移位移，其中表示均匀分布。这确保集合在空间上统计均匀。得到的信号为（25）为了补充宽带信号，我们叠加五个离散音调：（26）具有振幅、频率、空间波数和随机相位。这些时间序列产生尖锐的频谱峰，为宽带信号添加音调成分。完整的信号通过叠加宽带和音调分量得到（27）带有预定的标量权重用于两个宽带信号，以将它们分为主导和次主导模态。在所示结果中，我们设置，，。对于两个宽带信号，我们使用，，。添加了五个音调，其参数为。 B. 宽带-音调腔流为了在实验数据上展示该方法，bSPOD应用于湍流开放腔流的高速单点PIV数据集[51]。该数据集先前已被用于基准测试谱估计技术，例如基于多锥估计的自适应SPOD算法以改进音调识别[44]，以及Zhang等人[51]利用高速压力测量估计非时间分辨PIV快照的谱POD模态。腔体深度为，轴向范围为，面外宽度为。PIV测量在等间距网格上提供中心平面轴向和垂直速度分量，共有个快照，采样频率为。图2显示了时间平均轴向和垂直速度场、波动垂直速度分量的一个快照，以及波动垂直速度分量的谱功率密度。提取该分量的位置在右上角的图中用标记标出。PSD使用Welch方法计算，窗口长度为1600个样本，重叠率为50%。在940、1540和2130 Hz处观察到三个明显的峰，分别对应于第二、第三和第四罗斯特模态。图2. Zhang等人[51]的马赫0.6腔流单点PIV数据集。左图：时间平均轴向（上）和垂直（下）速度场；粗黑线标记腔体边界。右图：波动（均值减去的）垂直速度的快照（上）和绿色点处波动垂直速度的PSD（下）。红色垂直线表示三个罗斯特模态（940、1540和2130 Hz）。IV. 结果 A. bSPOD的特征为了突出bSPOD的特性，我们首先考虑人工宽带-音调数据集。图3顶面板显示了来自基于Welch方法的SPOD和bSPOD的前三个模态的特征值谱，以及根据信号构造已知的解析真实值。对于宽带分量，显示了空间积分功率谱密度；对于音调分量，显示了积分音调功率，因为它们的能量集中在单一频率上。SPOD和bSPOD结果显示在，这产生了细密的频率分辨率。两种方法最初都在没有重叠或加窗的情况下应用，以比较核心算法而不添加额外特征。在每个频带内，bSPOD模态根据公式(16)分配频率，权重指数为。由于音调分量的谱密度对应于狄拉克脉冲，我们将两种SPOD变体的相关峰特征值乘以频率间隔，以便与音调的总功率进行直接比较。图3中红色框内的相应特征值代表音调功率。图3. 人工数据集的模态功率谱密度和音调功率。将基于Welch方法的SPOD和bSPOD与三种配置的理论值进行比较：其中（上），其中（中），以及使用Hann锥形化的基于Welch方法的SPOD（下）。蓝色和橙色标记分别表示来自基于Welch方法的SPOD和bSPOD的前三个特征值。红色框内的标记代表音调功率（右轴），而其他所有标记和线表示功率谱密度（左轴）。黑线显示解析宽带功率谱密度（），黑色方块表示真实音调功率（）；参见第III.A节。虚线垂直线标记音调频率。两种方法都准确捕捉了音调的频率和功率。两个宽带分量的谱也得到了很好的表示，尽管观测到较高的估计方差，这在特征值估计（标记）围绕期望的ground-truth值（黑色实线）的噪声分布中可见。高方差源于块数较少和短带集成滤波器宽度，这反映了频率分辨率与统计收敛之间的预期权衡。前导和第一子导宽带模式分别由相应的前导和子导SPOD（bSPOD）特征值复现。在包含音调的频带中，主导SPOD（bSPOD）模式代表音调，使模式排序移动一位：前导宽带分量作为第一子导模式出现，子导宽带分量作为第二子导模式出现，依此类推。这些移动后的SPOD和bSPOD特征值在图3的顶部面板中用绿色标记边缘突出显示。总体而言，所示SPOD和bSPOD结果之间的差异很小，除了bSPOD显示更少的谱泄漏。这可以从第二子导模式（方形标记）的幅度中观察到，该幅度纯粹是泄漏的结果，对于bSPOD比SPOD低约一个数量级。这种减少源于在傅里叶变换中使用全时间序列，而不是基于Welch的SPOD中具有时间分割的较短窗口。

为了降低估计方差，块数和带集成滤波器长度增加了四倍，这使两种方法的频率分辨率以相同因子降低。相应的特征值谱如图3中间图所示。对于bSPOD，结果再次针对$N_f=...$（式(15)）绘制。如具有更高频率分辨率的图中所示，bSPOD准确复现了音调和宽带分量，现在宽带信号的方差显著降低。相比之下，基于Welch的SPOD特征值偏离解析参考：频率分辨率不足以解析音调，短窗口导致严重泄漏效应。这不仅影响音调峰值附近的区域，还影响较高频率处的子导模式，其中泄漏掩盖了物理模式与伪模式之间的分离。

在基于Welch的SPOD中，谱泄漏通常通过对每个时间片段施加锥形窗来缓解。图3底部面板中的蓝色标记显示了使用汉宁窗作为锥形时获得的相应特征值。bSPOD特征值相对于中间面板不变，并再次绘制以进行直接比较。锥形化显著降低了远离音调频率处的谱泄漏；然而，由泄漏引起的伪次级子导模式分支对于基于Welch的SPOD仍然比bSPOD更明显。此外，在基于Welch的SPOD中锥形化的副作用是增加了音调附近区域的泄漏。这是因为加窗对应于与窗频率响应的卷积：与矩形窗相比，汉宁窗具有更宽的主瓣，这增加了音调峰值附近的谱涂抹[48]。

在bSPOD中，通过计算全时间序列上的傅里叶变换自然减少泄漏，无需锥形化。在bSPOD内施加锥形会进一步降低伪模式分支的幅度。总体而言，这些结果表明，即使施加锥形化，bSPOD相对于基于Welch的公式也能产生更好的频率选择性和降低的谱泄漏。

B. 应用于实验数据
在将所提出的方法应用于人工测试数据集后，我们转向更相关的情况，并将基于Welch的SPOD和bSPOD应用于均值减去腔流速度快照，$N_t=...$和$N_f=...$。为了公平比较，方法再次配置为$N_f=N_{blk}$；即，基于Welch的SPOD中的块数等于bSPOD中考虑的连续傅里叶模式数。不应用重叠，确保两种方法产生相同数量的模式估计，具有相同的有效频率间距$\Delta f$。对于基于Welch的SPOD，使用汉宁窗来减少谱泄漏，而在bSPOD中，频率归属基于$N_f$；参见式(15)。

图4顶部面板显示了两种方法在频率分辨率为$\Delta f=...$（$N_f=...$）时获得的三个最能量模式的功率谱密度。虚线垂直线突出显示了三个Rossiter模式的频率，这些频率从点态PSD获得，如第III.B节所述。在这种高频率分辨率下，两种方法都很好地捕捉到三个音调Rossiter峰。然而，高频率分辨率以大的估计方差为代价，在远离音调的区域可见为升高的噪声水平，其中前导模式代表宽带动力学。观察到的噪声水平表明谱没有很好地收敛，需要更多模式，但这会降低频率分辨率。这种权衡在图4的中间和底部面板中说明，分别显示了对于$\Delta f=...$（$N_f=...$）和$\Delta f=...$（$N_f=...$）的相同谱。尽管两种方法经历相同的频率分辨率降低，但bSPOD通过数据驱动的频率归属保持准确的音调频率估计。

图4. 基于Welch（蓝色）和带集成（橙色）公式的前三个SPOD特征值随频率的变化，$N_f=N_{blk}$和三种频率分辨率：$\Delta f=...$（顶部）、$\Delta f=...$（中间）和$\Delta f=...$（底部）。垂直红线表示从单点PSD识别的Rossiter模式频率（见图2）。黄色箭头和圆圈标记了其对应模式形状在第IV.C节中显示的特征值。

相反，这种带内频率信息对于Welch-SPOD不可用。此外，Welch方法仍然容易受到谱泄漏的影响：来自相邻音调的能量污染谱估计，导致音调分量的有偏表示。这些结果突出了bSPOD在分析宽带-音调流中的优势。通过减少谱泄漏并利用展开系数中包含的信息，bSPOD实现比基于Welch的SPOD公式更清晰的谱定位和音调与宽带贡献的更清晰分离。

C. 模式收敛
在讨论了特征值谱及其收敛行为之后，我们现在关注模式形状。对应于图4中黄色箭头和圆圈标记的特征值的模式如图5和图6所示。这些图比较了从$N_f$个相邻傅里叶模式计算的bSPOD模式与在相应区间的下界傅里叶频率处评估的Welch-SPOD模式。与上界频率处评估的模式的比较得出相同的结论，因此为简洁起见省略。

图5. 使用带集成SPOD（左）和基于Welch的SPOD（右）在$\Delta f=...$的有效频率分辨率下，在$f=1540$ Hz处的Rossiter模式的垂直速度分量的实部。相应的特征值由图4顶部面板中的黄色十字标记。

图6. 在两种有效频率分辨率下计算的前导宽带模式的垂直速度分量的实部：$\Delta f=...$（顶部）和$\Delta f=...$（底部）。结果分别显示带集成SPOD（左）和基于Welch的SPOD（右）。相应的特征值由图4顶部（$\Delta f=...$）和底部（$\Delta f=...$）面板中的黄色十字标记。

我们首先考虑精细频率分辨率情况$f=1540$ Hz下的音调分量，对应于图4顶部面板中显示的特征值谱。图5显示了垂直速度模式的实部，左侧为bSPOD结果，右侧为Welch-SPOD结果。正如预期，两种模式形状几乎相同，因为它们是从等大小的CSD矩阵估计的，并且与相同的频率范围（1525–1550 Hz）相关联。有趣的是，尽管谱中的方差很高（见图4顶部面板），两种模式似乎都很好地收敛。这是因为音调分量在该频率区间中携带大部分功率，对交叉谱密度产生明确定义的贡献，进而产生清晰可识别的特征向量，即使在高的频率分辨率下也是如此。

对于功率较低的宽带动力学，观察到不同的行为。图6顶行显示了在相同频率分辨率下，在425–450 Hz频带中前导模式为宽带（见图4顶部面板）时SPOD和bSPOD的垂直速度模式的实部。这里，高频率分辨率以及因此较少的傅里叶模式（$N_f=...$）对CSD矩阵的估计有贡献，导致噪声和收敛性差的模式估计。重要的是，SPOD和bSPOD模式表现出相似的噪声水平和整体质量，表明在这种情况下两种方法都没有提供明显优势。图6底行显示了较低频率分辨率$\Delta f=...$下的相应模式形状，使得该模式代表400–500 Hz频带中的主导动力学。对于两种方法，频率间距的增加增加了特征值问题的规模（$N_f=...$），意味着更多正交模式有助于表示频带内的动力学。如SPOD通常观察到的，这同样降低了Welch和带集成方法在图6中看到的噪声水平。

需要注意的是，这种改进是主导流动力学在所考虑的频带中不明显变化的结果；对于有限数据集，增加实现次数，无论是通过添加额外块（SPOD）还是连续傅里叶模式（bSPOD），必然增加对每个谱估计有贡献的有效频率范围。如果模式形状随频率强烈变化，例如由于主导波长的强色散，增加有效间距会混合不相似的动力学，并可能抵消增加实现带来的好处。在这种情况下，只能通过考虑更长的时间序列来实现更好的收敛。

为了证明两种方法对于子导模式也产生非常相似的模式形状，图7显示了在800–900 Hz频带中的第一子导bSPOD模式的实部，并将其与在$f=825$ Hz处计算的相应SPOD模式进行比较。两种模式都表现出一些噪声，这由于子导模式的能量较低而符合预期。然而，它们的整体模式形状非常相似。

图7. 在800–900 Hz频带中第一子导宽带模式的垂直速度分量的实部。结果显示带集成SPOD（左）和基于Welch的SPOD（右）。相应的特征值由图4底部面板（$\Delta f=...$）中的黄色十字标记。

最后，图8比较了SPOD和bSPOD相对于可用数据量的模式收敛行为。因此，数据集从100%逐步缩短到25%，步长为25%，同时保持频率间距恒定在$\Delta f=...$。这意味着表示每个有效频率带内动力学的模式数随数据减少成比例减少。

图8. 通过将快照数减少到完整时间记录的75%、50%和25%来评估前导（顶部）和第一子导（底部）SPOD特征值的模式收敛。结果显示带集成SPOD（左）和基于Welch的SPOD（右）。收敛性通过保持频率分辨率固定并将从缩减时间序列计算的前导SPOD模式形状与从完整时间记录获得的相应模式对齐来量化[参见式(28)]。

图8顶行显示了从缩减数据集$\psi_{\text{red}}$和完整数据集$\psi_{\text{full}}$获得的前导模式（最大特征值）之间的对齐度，定义为
$$ \text{align} = |\langle \psi_{\text{red}}, \psi_{\text{full}} \rangle| \quad (28) $$
这是常用的模态相似性度量[52]。对齐值为1意味着模式完全对齐，而值0表示模式形状正交。底行显示了第一子导模式的相应对齐度。如预期，图8显示两种方法的对齐度都随着考虑的数据量增加而增加。尽管个别频率的确切值不同，SPOD和bSPOD表现出可比的收敛行为。在给定的数据缩减水平下，第一子导模式（底部）的整体对齐度低于前导模式（顶部），表明子导模式的收敛通常需要更多数据。这与其较低的能量含量和因此较低的信号噪声比一致。然而，SPOD和bSPOD的整体收敛趋势对于子导模式也保持非常相似。

我们得出结论，在等效条件$N_f=N_{blk}$下，对应于相同的有效频率带宽，动力学由相同数量的POD模式表示，SPOD和bSPOD产生相似的模式形状和可比的收敛特性。

D.频率归因

到目前为止，我们已经看到，bSPOD在实现与基于Welch的SPOD相当的模态收敛性的同时，通过频率归因步骤减少了频谱泄漏并改进了特征值的频率估计。这一优势源于bSPOD由连续的傅里叶模态构成互谱密度（CSD）估计，因此快照法特征值（展开系数，等式(13)）量化了每个离散傅里叶模态对给定bSPOD模态的贡献（等式(14)）。归一化后，这些系数充当频率权重，bSPOD利用这些权重推断每个模态的带内频率估计（等式(16)）。本节中，我们将详细阐述宽带空腔流动的频率归因步骤。

图9展示了前导空腔流动模态在 和两种分辨率（ 和 ，分别对应上部和下部）下的频率权重 。相应的特征值谱在上面图4的中间和底部面板中给出。权重绘制为对800至1700 Hz之间每个频带贡献的离散傅里叶频率的函数；偶数和奇数频带分别用黑色和灰色表示，以强调按频带的定义。在大多数频带中，权重围绕均值 （绿色水平线）波动，没有明显的模式。然而，在包含音调的频带中（用红色垂直线标记），音调附近的权重升高，并对bSPOD模态的贡献更强，反映了音调更高的功率及其与周围宽带动力学的清晰可分性。

图9. 前导频带集合模态的频率权重 与频率的关系，对应 和 （上部）以及 （下部）；参见等式(15)。交替的灰白背景和交替的黑/灰曲线表示bSPOD中使用的频带。蓝色和黄色垂直线分别标记使用 和 获得的bSPOD频率估计。红色垂直线表示罗西特模态频率（见图2）。

蓝色虚线标记使用 [等式(16)] 的归因步骤得到的bSPOD频率估计，从每个频带中贡献最大的单个傅里叶模态中选择频率。在含音调的频带中，这些估计非常接近罗西特频率（红色）。在 处，两个显示的罗西特频率都靠近频带上限。算法考虑到了这一点，将接近音调（红色）的频率（蓝色）不仅分配给包含该音调的频带，也分配给相邻频带，将归因频率置于相邻频带的下边缘。这产生了比基于Welch的SPOD更尖锐的模态功率密度（图6，中部），后者使用固定频率网格。

图9中的黄色虚线显示了使用所示权重（ ）得到的bSPOD频率估计。这些估计在表示音调频率方面不如 估计准确，但相比固定选择（如频带中心频率）仍然提供了更好的频率估计。在没有音调的宽带频带中， 估计倾向于选择更接近频带中心的频率，这与功率在频带内分布一致。在这些宽带部分，即使贡献接近均匀， 也可能产生接近频带边缘的频率估计。因此，期望对含音调频带依赖 ，对宽带频带依赖 ；这样的切换可以通过特征值功率谱密度随频率的局部变化来触发。在所有情况下，频率权重的归一化（等式(15)）确保估计频率保持在各自频带内，因此 和 都不会产生带外频率估计。因此，与依赖固定频率网格的基于Welch的SPOD相比，频率估计只能得到改善。

E. 频带集合SPOD中的重叠

对于基于Welch的SPOD，通常对时间片段进行重叠以增加有限时间信号的傅里叶估计数量。同样的原理可用于bSPOD，通过将频带集合窗口移动小于 ，使连续频率块重叠。这产生额外的估计，从而细化模态功率谱的频率网格。

然而，应用重叠时需要考虑两点。首先，尽管bSPOD模态以与所选频带长度 相关的频率间隔定义，但功率估计的有效间隔会随重叠而改变，因为傅里叶模态在重叠块中被重用。实际上，在频率上积分功率密度以恢复信号总功率时，必须考虑导致频率分辨率提高的重叠因子。

其次，在包含音调分量的宽带信号中，如果所有模态都简单地分配给频带中心频率（ ），重叠会人为地将音调的能量扩展到多个bSPOD频带。对于使用频率平滑的PSD估计，这会产生频谱中的块状跳跃，这可能是PSD更常通过Welch分段而非频率平滑估计的原因之一。在bSPOD中，这个问题由数据驱动的频率归因步骤避免，该步骤由POD展开系数启用。它确保每个bSPOD模态与其贡献最大的频率相关联，并使音调功率正确集中。

为了说明这种行为，图10展示了空腔流动数据的前三个bSPOD特征值，对应 （ ；与图4底部面板比较）和75%的频带重叠。与无重叠的情况相比，频率分辨率提高了四倍。需要注意的是，相邻模态部分共享相同的傅里叶分量，导致统计相关的估计。结果显示分别使用 和 的频率归因，以及无归因的情况（每个bSPOD模态被分配频带中心频率 ）。

图10. 空腔流动数据的前三个bSPOD特征值与频率的关系，使用 （ ）和75%频带重叠计算。特征值分别显示为固定频率网格（ ，将每个特征值分配给频带中心）和两种数据驱动的频率归因， 和 [见等式(15)]。

不出所料，没有频率归因（ ）时，频谱类似于经典频率平滑：音调峰表现为块状升高。使用 时，音调模态被归因到接近罗西特频率的频率，从而在频带重叠的情况下仍保持频谱中尖锐的音调峰。然而，在频谱的宽带部分， 可能产生不均匀或不连续的频率趋势，包括小跳跃，这在图10中的第一个次主导模态中尤为明显。使用 产生平滑且分辨率良好的宽带频谱，但音调峰的定位不如 尖锐，尽管仍优于 。这些观察结果表明了一种混合策略：在宽带频带中使用 ，在含音调频带中使用 。如上所述，这种切换可以通过检测模态功率谱密度随频率的强局部变化来触发。

我们得出结论，在bSPOD中，相邻频带之间引入重叠可以增加模态谱的有效频率分辨率，而不会将音调分量涂抹到相邻频带，前提是使用适当的频率归因。在基于Welch的SPOD中，重叠通常用于增加固定长度块的数量。或者，基于Welch的SPOD中的重叠可以解释为在固定块数下允许更长的块长度，从而对频率分辨率产生类似的效果。bSPOD的关键优势在于重叠是在傅里叶变换之后应用的，并且可以在整个频谱上变化，而基于时间分段的SPOD依赖于固定频率网格，不提供这种灵活性。

V. 结论

本工作提出了一种计算统计平稳数据SPOD模态的方法。所提出的方法称为频带集合SPOD，灵感来自频率平滑，一种降低功率谱密度估计估计方差的方法。这里，该方法扩展到SPOD，使得能够估计空间模态及其模态功率谱密度。与更常用的基于Welch的SPOD公式不同，后者通过对时间分段获得的同一频率箱上的傅里叶模态估计互谱密度矩阵，bSPOD则从窄间隔的相邻傅里叶系数近似互谱密度矩阵，从而近似SPOD模态。这些系数从全时间序列的单一傅里叶变换获得，无需分段。

与基于Welch的SPOD相比，bSPOD具有几个优势。首先，通过避免时间分段，频谱泄漏自然减少。其次，每个傅里叶系数对每个bSPOD模态的贡献可以通过展开系数跟踪，从而保留带内频率信息。该信息可用于获得各个bSPOD模态的数据驱动频率估计。通过结合这些特性，bSPOD降低了估计方差，同时保持对音调分量的低偏差，相对于基于Welch的公式，在偏差-方差权衡中提供了有利的平衡。因此，bSPOD特别适用于宽带音调流动，其中希望在不损害音调特征值估计精度的情况下降低方差。

bSPOD公式进一步允许通过改变滤波器长度随频率的变化来调整这种权衡。这实现了可变频带上的模态分解，具有随频率变化的模态数量。这种灵活性使bSPOD非常适合基于收敛标准选择模态数量的自适应策略，如Yeung和Schmidt [44]所建议的。鉴于这些特性以及与基于Welch的SPOD相当的计算成本，bSPOD为湍流流动的谱模态分析提供了实用的改进。