弱形式DMD：基于Galerkin投影与积分平滑的抗噪声模态分解方法

1. 项目概述：从“快照”到“弱形式”的DMD进化

动态模态分解（Dynamic Mode Decomposition, DMD）在过去十几年里，已经从一个流体力学领域的专用工具，演变成了数据驱动动力学建模的“瑞士军刀”。无论是分析实验中的粒子图像测速（PIV）数据，还是从金融时间序列中提取主导模式，DMD都能提供一套优雅的数学框架，将复杂的高维动态数据分解为一系列具有固定频率和衰减/增长率的模态。其核心魅力在于，它完全从数据出发，无需预先知道系统的控制方程，就能挖掘出内在的、线性的动力学结构。

然而，这把“军刀”在实际打磨中遇到了两块硬骨头。第一块是噪声。无论是传感器误差、环境干扰还是数值计算中的舍入误差，噪声都会污染我们采集的“快照”数据矩阵。传统的DMD算法，尤其是基于奇异值分解（SVD）的精确DMD，对噪声异常敏感。噪声会污染我们估计的Koopman算子，导致提取的模态频率失真、衰减率计算错误，甚至产生完全虚假的、由噪声主导的模态。第二块是时间步长限制。DMD的理论基础建立在均匀时间采样的假设上。当数据的时间步长（Δt）过大时，会引发所谓的“混叠”问题，即高频动力学信息被错误地折叠到低频模态中，导致模态识别完全失败；而时间步长过小，则可能捕捉不到系统缓慢演变的特征，且会积累更多的数值误差。

我最初接触DMD时，就被这两个问题困扰了很久。直到深入研究了基于Galerkin方法的弱形式DMD（Weak-form DMD），才找到了一个更具鲁棒性的解决方案。这个方法的精髓，不在于对算法本身进行小修小补，而是从根本上改变了我们“看待”数据的方式——从直接处理离散的快照点，转向处理这些快照在时间上的“积分”或“加权平均”，也就是所谓的“弱形式”。这听起来有点抽象，但你可以把它想象成：与其相信每一个可能被噪声污染的瞬时测量值，不如去相信数据在短时间窗口内表现出来的整体趋势。这种方法不仅能够有效平滑噪声，更重要的是，它放松了对数据采样率的苛刻要求，为处理非均匀采样甚至大数据间隔的数据提供了新的可能。接下来，我将详细拆解这套方法的思路、实现细节以及我在实际应用中的心得体会。

2. 核心思路：Galerkin投影与弱形式积分

要理解弱形式DMD，必须先理解两个核心概念：Galerkin投影和弱形式。这并非DMD的独创，而是计算数学和有限元方法中的经典思想。

2.1 Galerkin投影：在子空间中寻找最佳近似

首先，我们明确DMD的目标：从一个高维状态空间的时间序列数据中，找到一个线性算子A（即Koopman算子的有限维近似），使得 x_{k+1} ≈ A x_k。这里的x_k是第k个时间步的快照（一个列向量）。

Galerkin方法的核心思想是投影。我们并不要求在全局状态空间中都严格满足 x_{k+1} = A x_k，而是将这个等式“投影”到一个精心选择的、维度更低的“试探函数空间”或“权重函数空间”中，要求在这个子空间的意义下等式成立。这就好比，我们无法精确描述一个复杂曲面的每一点，但我们可以用一组简单的基函数（比如多项式）去拟合它，并保证在某种平均意义下（比如最小二乘）误差最小。

在弱形式DMD的语境下，我们构造两个数据矩阵X和Y，其中X的每一列是x_k，Y的每一列是x_{k+1}。传统DMD直接求解最小二乘问题：min ||Y - AX||_F。这等价于要求每个快照对的残差 (y_k - A x_k) 在欧几里得范数下最小。而Galerkin方法则引入一个权重矩阵W，将问题转化为：找到A，使得加权残差W^T(Y - AX)在某种意义下为零或最小。这个W就代表了我们的“试探函数空间”。通过巧妙地设计W，我们可以将噪声的影响平均掉，或者融入我们对系统动力学的先验知识。

2.2 弱形式积分：从点到区间，平滑噪声与混叠

“弱形式”是第二个关键。在微分方程数值解中，弱形式通过对方程两边乘以一个“试验函数”并进行积分，将微分方程转化为积分方程。积分过程天然具有低通滤波和平滑的效果，因为高频振荡在积分区间内会正负相消。

我们将这一思想移植到时间序列数据上。传统的“强形式”数据对是瞬时值：(x(t_k), x(t_{k+1}))。而“弱形式”数据对，则是状态向量在时间窗内的积分或加权平均。例如，我们可以定义： [ \tilde{x}k = \int{t_{k}-\tau}^{t_{k}+\tau} \phi(s) x(s) ds ] 其中φ(s)是一个紧支撑的权重函数（比如一个钟形函数或分段多项式）。这样构造出的新数据对 (\tilde{x}k, \tilde{x}{k+1})，其包含的噪声是原数据在时间窗内噪声的积分平均，方差显著降低。同时，对于时间步长问题，如果我们选择的积分窗口τ与系统的特征时间尺度相匹配，即使原始采样点稀疏，这个积分过程也能在一定程度上捕捉到动力学的演变趋势，缓解因采样不足导致的混叠效应。

两者的结合：弱形式DMD的本质，就是先对原始快照序列进行弱形式处理（时间积分/卷积），得到一组平滑的、抗噪声的“广义快照”，然后再对这组新快照应用基于Galerkin投影的DMD算法。这里的Galerkin投影体现在权重函数φ(s)的选择上，不同的φ(s)定义了不同的“观察”系统动力学的视角。

注意：在实际数值实现中，我们处理的都是离散数据。因此，“积分”操作会退化为离散卷积。假设我们有连续m个时间点的快照数据，我们可以通过一个滑动窗口和一组卷积核权重来生成弱形式数据矩阵。这是整个算法实现中最需要仔细设计的部分。

3. 算法实现：从理论到代码的四个关键步骤

理解了核心理念，我们来看如何一步步实现它。我将结合MATLAB/Python风格的伪代码进行解释，重点说明每个步骤的意图和注意事项。

3.1 第一步：数据预处理与弱形式变换

这是区别于传统DMD的最关键一步。假设我们有一个高维状态的时间序列，排列成矩阵X_raw，其每一列x_raw(:, k)对应时刻t(k)的快照。

1. 选择权重函数（卷积核）： 权重函数φ的选择决定了弱形式的特性。常见选择有：

• 矩形窗：phi = ones(window_length, 1)/window_length。最简单，相当于简单移动平均。对平滑白噪声有效，但频域特性较差（旁瓣高）。
• 高斯窗：phi = gausswin(window_length)。具有优良的时频局部性，平滑效果自然。这是我最常使用的。
• 多项式窗（如基于Legendre多项式）：能与高阶Galerkin方法结合，在保留动力学信息方面更有优势，但更复杂。

2. 执行离散卷积（弱形式变换）： 对X_raw的每一行（即每个状态变量随时间的变化）独立进行与权重函数φ的卷积操作。这生成了弱形式数据矩阵X_weak。

% 假设 X_raw 大小为 [n, m]， n为状态维度， m为时间点数 window_len = 5; % 滑动窗口长度，需为奇数以便对称 phi = gausswin(window_len); % 高斯窗 phi = phi / sum(phi); % 归一化，保证不影响均值 X_weak = zeros(n, m); for i = 1:n X_weak(i, :) = conv(X_raw(i, :), phi, 'same'); % ‘same’ 保持输出长度与输入一致 end

注意事项：

• 卷积操作会使数据序列两端（边界）失真，因为缺少足够的数据进行完全卷积。处理方式有：① 直接舍弃两端数据（最稳妥）；② 进行数据镜像对称填充后再卷积。在DMD中，由于我们通常需要成对的数据矩阵，直接舍弃两端是最简单的，但要确保舍弃后仍有足够的数据量。
• 窗口长度window_len是超参数。太小则平滑效果不足，太大则会过度平滑，抹除真实的快速动力学特征。一个经验法则是将其设置为系统主要频率对应周期的1/4到1/2采样点数。

3.2 第二步：构建数据矩阵与降维

经过弱形式变换后，我们得到平滑的数据X_weak。接下来，像传统DMD一样，我们构建两个数据矩阵：

X = X_weak(:, 1:end-1); % 前m-1个弱形式快照 Y = X_weak(:, 2:end); % 后m-1个弱形式快照，与X一一对应

现在，我们处理的是Y ≈ A X的关系，但这里的X和Y已经是抗噪声能力更强的“广义快照”。

降维（SVD与截断）： 即使经过平滑，为了数值稳定和提取主导模式，我们通常对X进行奇异值分解（SVD）并截断。

[U, S, V] = svd(X, 'econ'); r = 15; % 截断秩，可根据奇异值谱的“拐点”或能量占比（如99%）确定 Ur = U(:, 1:r); Sr = S(1:r, 1:r); Vr = V(:, 1:r);

实操心得： 弱形式变换后，数据矩阵X的奇异值衰减通常会比原始数据更快、更平滑。这意味着所需的截断秩r可能比处理原始数据时更小，因为噪声贡献的奇异值被有效抑制了。你可以通过绘制奇异值（diag(S)）来直观感受，弱形式数据的奇异值曲线在初始快速下降后，会更快地进入一个平坦的“噪声平台”。

3.3 第三步：在降维子空间中求解近似Koopman算子

这是Galerkin投影思想的核心体现。我们不直接在高维空间求A，而是在由Ur张成的r维主导子空间中，求解一个降维的算子Ã。

将A投影到Ur的子空间上：à = Ur' * A * Ur。但我们不知道A。利用Y ≈ A X和X ≈ Ur * Sr * Vr'，我们可以推导出Ã的近似表达式：

% 计算降维的Koopman算子近似 Atilde = Ur' * Y * Vr / Sr; % 这是核心公式之一 % 等价于 Atilde = Ur' * Y * Vr * inv(Sr);

这个Ã是一个r×r的矩阵，它描述了在主导模态坐标下的线性动力学。相比于直接求n×n的A（n可能成千上万），求解Ã不仅计算量小，而且数值稳定性极高。

为什么这是Galerkin投影？我们可以这样理解：Ur的列向量构成了我们的试探函数空间基。要求解A，我们强制残差Y - A X与这个空间正交（即投影为零）。数学上，这等价于Ur' (Y - A X) = 0。结合X的SVD近似，就导出了上面的Atilde公式。因此，我们是在Ur空间的意义下，找到了一个最优的线性动力学近似。

3.4 第四步：特征分解与模态重构

得到Ã后，我们对其进行特征分解：

[W, Lambda] = eig(Atilde); % W是特征向量矩阵，Lambda是对角特征值矩阵

Ã的特征值μ_j包含了系统的动力学信息。连续时间的增长率λ_j和频率ω_j可以通过以下关系得到：

dt = t(2) - t(1); % 原始数据的时间间隔 lambda_j = log(diag(Lambda)) / dt; % 连续时间特征值 growth_rate = real(lambda_j); % 增长率/衰减率 frequency = imag(lambda_j) / (2*pi); % 频率 (Hz)

关键的一步：重构高维DMD模态。Ã的特征向量W是r维子空间中的模态。我们需要将它们映射回原始的高维物理空间。高维的DMD模态Φ（每一列是一个模态）可以通过下式重构：

Phi = Y * Vr / Sr * W; % 这是另一个核心公式

为什么是这个公式？可以从A Φ ≈ Φ Λ和à = Ur' A Ur的关系推导出来。这个重构公式确保了模态与特征值Λ的对应关系。重构出的Phi的每一列，就是我们可以可视化或解释的物理空间中的空间模式。

至此，我们得到了弱形式DMD的全部输出：特征值（包含频率和增长率）和对应的空间模态Phi。

4. 参数选择与调优经验

弱形式DMD引入了几个关键超参数：权重函数类型、窗口长度和SVD截断秩。它们的设置直接影响结果质量。

4.1 权重函数与窗口长度的联合调优

权重函数和窗口长度共同决定了弱形式变换的“滤波器”特性。

• 场景一：抑制高频测量噪声如果你的数据含有明显的高斯白噪声，目标是平滑。建议使用高斯窗。窗口长度的选择可以参考数据的时间自相关函数。一个实用的方法是：逐步增加窗口长度，观察重构数据（将DMD模态和动力学叠加回去）与原始平滑数据的误差。当误差不再显著下降时，即到达了“收益递减”点。通常，窗口长度覆盖2-3个你感兴趣的最短周期对应的采样点即可。

• 场景二：处理非均匀采样或大时间步长当数据点稀疏时，矩形窗或简单的线性权重窗可能更鲁棒。此时的目标不是平滑，而是构造一个合理的“时间局部平均”来代表该时刻的状态。窗口长度应大致等于（或略小于）采样间隔。这相当于承认我们无法分辨该间隔内的细节，转而使用一个“区块”代表符。

• 场景三：分离重叠频率的模式如果你怀疑系统中有频率非常接近的模式，并且被噪声掩盖，可以尝试使用频谱旁瓣更低的窗函数，如汉明窗或布莱克曼窗。但这通常需要更长的窗口长度，并可能牺牲一些时间分辨率。

我的经验法则：在无先验知识时，从高斯窗开始。将窗口长度设置为总时间点数的5%~10%作为初始尝试。然后，通过检查DMD特征值的分布来诊断：如果特征值大量分布在单位圆（离散时间）或虚轴（连续时间）附近一个模糊的带状区域，说明噪声仍占主导，可能需要增大窗口长度或改变窗函数。

4.2 SVD截断秩r的确定

截断秩r的选择永远是在“保留信息”和“剔除噪声”之间权衡。

1. 奇异值谱观察法：绘制X（弱形式变换后）的奇异值，看是否存在明显的“拐点”或“肘部”。拐点之前是信号主导，之后是噪声平台。将r设定在拐点处。弱形式数据的优势在于，这个拐点通常比原始数据更明显。
2. 能量占比法：设定一个阈值，如99%，选择最小的r使得前r个奇异值的平方和占总平方和的比例超过该阈值。即sum(sigma(1:r).^2) / sum(sigma.^2) > 0.99。
3. 后验验证法：尝试不同的r值，计算DMD重构数据与原始数据（或弱形式数据）的误差。绘制误差随r变化的曲线，选择误差曲线开始进入平台期对应的r。同时，观察DMD模态的物理合理性。如果r过大，可能会出现空间结构杂乱无章、无法解释的噪声模态。

一个综合策略：我通常先根据奇异值谱定一个范围，然后在这个范围内选取几个r值分别计算。重点关注那些增长率接近零（对应稳定振荡）或为较小负值（轻微阻尼）的模态，它们的空间结构是否稳定、可解释。选择那个能产生最清晰、最稳定物理模态的r。

4.3 与正则化技术的结合

弱形式本身是一种正则化。但有时为了进一步增强鲁棒性，尤其是当数据量很少时，可以结合其他正则化技术。

• Tikhonov正则化：在求解最小二乘问题min ||Y - AX||_F时，加入一个惩罚项α||A||_F。这可以在算法第三步求解Ã时融入。对于弱形式DMD，由于数据已经平滑，所需的正则化参数α通常非常小，甚至为零。
• 总最小二乘DMD：传统DMD是普通最小二乘，假设误差只在Y中。总最小二乘同时考虑X和Y的误差，理论上更适用于双方都有噪声的情况。弱形式变换后，双方误差相关性增强，此时标准DMD已足够，总最小二乘带来的额外计算开销可能得不偿失。

我的建议是：优先调优弱形式参数（窗函数和长度），这通常是效果最显著的。只有在参数调优后，仍存在少数明显由噪声引起的、增长率异常大的虚假模态时，才考虑引入一个很小的Tikhonov正则化进行压制。

5. 实战案例：含噪流体涡旋数据分解

为了让大家有更直观的感受，我模拟了一个经典的流体力学示例：一个平面内两个旋转频率不同、且有轻微阻尼的涡旋，叠加了高强度噪声。

数据生成：

1. 在二维网格上定义两个空间模态：两个位置不同的高斯涡旋。
2. 每个模态的随时间演化是exp((λ + iω)t)，其中 λ 为衰减率（负值），ω 为角频率。
3. 将两个模态的时空演化线性叠加，得到干净的流场序列X_clean。
4. 加入信噪比（SNR）为10dB的高斯白噪声，得到观测数据X_raw。

对比实验： 我们分别应用：

• 精确DMD：直接对X_raw操作。
• 弱形式DMD：对X_raw使用高斯窗（长度=7）进行弱形式变换后，再执行DMD。

结果分析：

这个案例清晰地展示了，在面对强噪声污染时，弱形式DMD如何通过“时间局部平均”这一思想，稳定地提取出被淹没的真实物理模式。

6. 常见陷阱与排查指南

即使理解了原理和步骤，在实际编码和应用中依然会踩坑。下面是我总结的一些典型问题及解决方法。

6.1 问题一：弱形式变换后，DMD模态全部是实模态（没有振荡）

现象：计算出的特征值没有虚部，或者虚部非常小，对应的模态是静止的空间模式，没有波动特性。

可能原因与排查：

1. 窗口长度过长：这是最常见的原因。过长的平滑窗口相当于一个极强的低通滤波器，把系统固有的振荡频率成分也过滤掉了。解决：逐步减小窗口长度，观察特征值频谱是否开始出现非零的虚部。
2. 权重函数选择不当：例如，使用了对称的矩形窗，且窗口长度与系统的振荡周期成整数倍关系，可能会在某些特定情况下导致相位信息抵消。解决：尝试换用高斯窗等具有更柔和边界的窗函数。
3. 数据本身可能就是非振荡的：首先确认你的原始系统是否真的存在振荡动力学。可以用傅里叶变换（FFT）检查一下主要时间序列的频谱。

6.2 问题二：重构误差很大，即使使用了很多模态

现象：用DMD模态和动力学重构出的数据X_reconstructed，与用于拟合的弱形式数据X_weak相比，误差仍然很大。

排查步骤：

1. 检查数据矩阵构建：确认X和Y是否正确对应。一个常见的低级错误是索引错位。
2. 检查SVD截断：截断秩r可能设得太小，丢失了关键动力学信息。检查奇异值谱，看看是否在设定的r之后还有较大的奇异值。尝试增大r。
3. 检查弱形式变换的边界效应：卷积的边界处理（如‘same’模式）可能在数据两端引入畸变，而这些畸变数据也被用于DMD拟合。解决：在弱形式变换后，主动舍弃数据序列开头和结尾的floor(window_len/2)个数据点，再用剩下的干净数据构建X和Y。
4. 系统非线性强：DMD本质是线性模型。如果系统非线性很强，单一的线性算子A无法很好地近似整个动力学。此时，弱形式DMD可能也无能为力，需要考虑非线性DMD变体（如Kernel DMD, Extended DMD）。

6.3 问题三：特征值物理意义不合理

现象：计算出的连续时间特征值λ的实部（增长率）非常大（正或负），或者频率高得离谱。

诊断与解决：

1. 时间单位：确认dt的单位是否正确。如果dt是毫秒，但你误以为是秒，那么计算出的频率会差1000倍。
2. 数值不稳定：如果Ã的条件数很大，其特征分解会非常敏感。这通常是因为SVD截断秩r选择过大，包含了数值上为零的奇异值。解决：检查Sr矩阵，确保其最小奇异值没有小到接近机器精度。可以设置一个更大的截断阈值。
3. 噪声模态：增长率绝对值很大的模态（无论是正增长还是负衰减）很可能是噪声模态。真正的物理模态通常具有较小的增长率（接近中性稳定或弱阻尼）。解决：根据增长率/衰减率对模态进行筛选，只保留那些|growth_rate| < threshold的模态进行分析。这个阈值可以根据具体物理问题设定。

6.4 问题四：如何处理非均匀采样数据？

弱形式DMD的一个潜在优势是处理非均匀采样数据，但需要调整。

方法：关键在于弱形式变换中的“积分”需要根据非均匀的时间点进行加权。将离散卷积改为非等间距的加权平均。 假设时间点为t(1), t(2), ..., t(m)，状态为x(t(1)), ..., x(t(m))。 对于目标时刻t_k，我们选择一个时间窗口[t_k - τ, t_k + τ]，然后找到所有落在这个窗口内的原始数据点{t_j}及其对应的状态{x(t_j)}。 弱形式状态\tilde{x}_k计算为：\tilde{x}_k = Σ_j w_j * x(t_j) / Σ_j w_j其中权重w_j可以根据距离t_k的远近，用高斯核等函数计算：w_j = exp(-(t_j - t_k)^2 / (2σ^2))。 这样，我们就为每个非均匀的t_k生成了一个“局部加权平均”的快照，从而构造出X_weak。后续的DMD步骤保持不变。这种方法实质上是将非均匀数据通过局部平滑，“映射”到了一个隐含的、更规则的时间基础上。