DMD实战指南：从流场快照到动态模态分解的完整实现

简介：这份资源是面向动力系统数据分析学习者与科研人员的MATLAB版动态模式分解（DMD）实现包，适合具备一定线性代数与MATLAB基础、希望将高维时间序列降维并提取低维动态模式的中高级用户。包内共3个文件，包含1个m脚本、1个pdf说明与1个md文档，压缩包约146KB，脚本可用于数据准备、矩阵构造、SVD谱分解、DMD模式与频率增益计算及重构预测，文档则辅助理解算法流程与关键概念。资源围绕流体动力学、信号处理、图像分析与机械结构健康监测等场景展开，读者可据此掌握从奇异值分解到DMD模式提取的完整链路，并了解扩展DMD、控制DMD等变体的思路，便于迁移到自身研究或工程问题中。目前已有501人学习下载，可作为入门DMD并动手复现的轻量参考。

1. DMD 不是玄学：从一段非定常流场数据说起

手里有一组随时间演化的流场快照，比如圆柱绕流的速度场，几百个时间步，每个时刻几万个网格点。你想知道：这里面到底藏着几种主导模态？它们的频率是多少？增长率是正还是负？传统做法是做 FFT，但 FFT 要求信号平稳，而且对瞬态、间歇性现象几乎无能为力。Dynamic Mode Decomposition（DMD）就是冲着这个问题来的——它把时间序列拆成若干空间模态，每个模态配一个特征值，实部管增长衰减，虚部管振荡频率。我第一次用它处理 PIV 实验数据时，最直观的感受是：终于不用靠肉眼在涡量图里数涡了。这篇笔记面向已经拿到快照矩阵、准备动手算的工程师，从数据排布讲到参数调优，再到踩过的坑，尽量把每一步都落到可复现的命令和代码上。

2. 快照矩阵怎么排：DMD 的数据入口与三种预处理

2.1 从原始数据到 X 和 X‘

DMD 的输入本质上就两个矩阵。假设你有 m 个时间步，每个时间步的场量展平成 n 维列向量，记第 k 个时刻为 x_k。那么：

X = [x_1, x_2, ..., x_{m-1}] X' = [x_2, x_3, ..., x_m]

X 和 X' 各是 n×(m-1) 的矩阵。DMD 要找的是一个最佳线性算子 A，使得 X' ≈ A X。注意这里没有常数项、没有控制输入，纯靠数据本身推演化。如果你的数据来自 CFD，n 可能是 10^5 量级，m 可能只有几百，矩阵极度扁长，所以后面必须做降维。

常见做法是先把每个快照减去时间平均值。不减均值的话，零频模态会混进结果里，特征值落在单位圆上，看起来像个“稳定模态”，实际上只是均值场。我一般会显式做这一步：

import numpy as np # snapshots: shape (n, m), 每一列是一个时刻的场 snapshots = np.load("flow_snapshots.npy") mean_field = snapshots.mean(axis=1, keepdims=True) X_all = snapshots - mean_field X = X_all[:, :-1] # 前 m-1 列 Xp = X_all[:, 1:] # 后 m-1 列 print(X.shape, Xp.shape)

逻辑说明：先减均值再切片，保证 X 和 X‘ 用的是同一个均值场。参数上，mean_field 的 keepdims=True 是为了保持 (n,1) 形状，方便广播。如果你的数据里有明显的瞬态起始段，比如前 50 步还没充分发展，建议直接丢掉，不要用滤波去补，DMD 对非物理瞬态很敏感。

2.2 归一化与加权：容易被忽略的一步

如果快照的各个分量量纲不同，比如速度场和压力场混在一起，或者网格非均匀，直接做 DMD 会让大量纲分量主导结果。常见做法是对每个分量做无量纲化，或者引入网格权重。我一般会构造一个对角权重矩阵 W，把内积改成加权内积。在均匀网格上 W 就是单位阵乘网格面积，非均匀网格上每个点的权重正比于其控制体积。

# 假设 dx, dy 是网格间距，二维场展平后每个点的面积权重 dx, dy = 0.01, 0.01 n_points = X.shape[0] weights = np.full(n_points, dx * dy) W = np.diag(weights) # 加权内积下的“能量”范数 energy = np.sqrt(np.real(np.sum(np.conj(X) * (W @ X), axis=0))) print("各快照加权能量:", energy[:5])

逻辑说明：W 的作用是让 DMD 在加权空间里找最优算子，避免网格密集区域被过度加权。参数上，dx 和 dy 要和你 CFD 或实验的实际分辨率一致。如果懒得算权重，至少在非均匀网格上做一次面积归一化，否则模态空间分布会失真。

2.3 降维：SVD 截断阶数怎么定

直接对 X 做 SVD：X = U Σ V*。U 的列是 POD 模态，Σ 的对角元是奇异值。DMD 的降维就是保留前 r 个奇异值对应的子空间。r 选多少？这是 DMD 里最像玄学的一步。我的血泪经验是：不要只看奇异值衰减曲线拐点，还要看你要提取的模态频率是否落在保留的子空间里。

U, S, Vh = np.linalg.svd(X, full_matrices=False) # 保留能量占比 99.9% energy_ratio = np.cumsum(S**2) / np.sum(S**2) r = np.searchsorted(energy_ratio, 0.999) + 1 print(f"保留阶数 r = {r}, 能量占比 = {energy_ratio[r-1]:.6f}") U_r = U[:, :r] S_r = np.diag(S[:r]) V_r = Vh[:r, :].conj().T

逻辑说明：energy_ratio 是累计能量占比，searchsorted 找到第一个超过 99.9% 的位置。参数 0.999 可以调，噪声大的数据用 0.99，干净数据用 0.9999。注意 SVD 的复杂度是 O(n m^2)，如果 m 很大，可以先对时间维度做降采样，或者用随机 SVD。我一般会同时看 r 从 10 到 50 的结果，如果主导频率对 r 不敏感，说明模态分离得好；如果频率随 r 跳变，说明截断太狠或数据信噪比不够。

3. 从算子到模态：精确 DMD 与投影 DMD 的代码实现

3.1 精确 DMD 的七行核心代码

精确 DMD 的思路是：在 POD 子空间里构造低维算子 Ã = U_r* X‘ V_r S_r^{-1}，然后对 Ã 做特征分解，再把特征向量投影回高维空间。下面是最小可复现版本：

# 精确 DMD Atilde = U_r.conj().T @ Xp @ V_r @ np.linalg.inv(S_r) eigvals, W_tilde = np.linalg.eig(Atilde) # 高维 DMD 模态 Phi = Xp @ V_r @ np.linalg.inv(S_r) @ W_tilde # 连续时间特征值 dt = 0.01 # 采样时间间隔 omega = np.log(eigvals) / dt # 按模态振幅排序 b = np.linalg.lstsq(Phi, X[:, 0], rcond=None)[0] amplitude = np.abs(b) idx = np.argsort(amplitude)[::-1]

逻辑说明：Atilde 是 r×r 小矩阵，特征分解很快。Phi 的每一列是一个 DMD 模态，omega 的实部是增长率，虚部是角频率。b 是初始时刻各模态的振幅，用最小二乘从第一帧快照投影得到。参数 dt 必须和你的快照间隔一致，单位是秒。如果 dt 设错，频率会整体缩放，这是最常见的翻车点之一。

3.2 投影 DMD：当 n 远大于 m 时的省内存写法

精确 DMD 里 Phi 是 n×r 矩阵，如果 n 是 10^6，r 是 100，Phi 就是 10^8 个浮点数，内存吃不消。投影 DMD 不显式构造 Phi，而是把模态表示为 POD 基的线性组合：Phi = U_r W_tilde。这样只需要存 U_r 和 W_tilde。

# 投影 DMD：模态用 U_r @ W_tilde 表示 Phi_proj = U_r @ W_tilde # 重构任意时刻的场：x(t) ≈ Phi_proj @ (b * np.exp(omega * t)) t = 0.5 x_recon = Phi_proj @ (b * np.exp(omega * t)) print("重构场形状:", x_recon.shape)

逻辑说明：Phi_proj 和精确 DMD 的 Phi 在数学上等价，但计算路径不同。投影 DMD 的代价是每次重构都要做 U_r @ W_tilde 的乘法，但省了显式存 Phi 的内存。参数上，如果你只需要频率和增长率，不需要空间模态，那连 Phi_proj 都不用算，直接看 omega 就行。

3.3 参数表：dt、r、窗口长度怎么配

这张表是我自己调参时的速查卡。dt 和 m 是数据采集阶段就定死的，r 和权重是后处理阶段可调的。如果频率结果对 r 敏感，优先检查数据质量和均值扣除，而不是继续调 r。

4. 避坑与排查：DMD 翻车现场记录

4.1 特征值全落在单位圆外，增长率大得离谱

现象：算出来的 omega 实部全是正的，而且数值很大，重构场几秒后就爆炸。原因：X 和 X‘ 的构造顺序反了，或者 dt 单位搞错（比如毫秒当成秒）。解决：检查 Xp 是不是 X 右移一列，dt 是不是实际采样间隔。我一般会先算一个已知频率的正弦信号做自检，比如 10 Hz 信号，dt=0.001，看 DMD 能不能还原出 10 Hz。

4.2 模态空间分布全是噪声，看不出物理结构

现象：Phi 的列看起来像随机噪声，没有涡结构。原因：SVD 截断阶数 r 太大，把噪声子空间也保留进来了。解决：降低 r，或者对 X 做一次低通滤波。另一个可能是快照数量太少，m < 2r，导致 Ã 欠定。经验规则是 m 至少是 r 的 3 倍。

4.3 频率出现混叠，高频模态对不上

现象：DMD 给出的频率和 FFT 对不上，或者出现负频率。原因：采样率不够，Nyquist 频率低于信号最高频。解决：提高采样率，或者对数据做降采样前先低通滤波。DMD 不会自动抗混叠，它只反映你给的数据。

4.4 重构场和原始场差很多，但频率看起来对

现象：omega 合理，但用 Phi 和 b 重构的场和原始快照对不上。原因：b 是用第一帧算的，如果第一帧恰好落在瞬态里，振幅就不准。解决：用最小二乘拟合所有时刻的振幅，或者从稳态段选一帧作为参考。我一般会做一次全时间最小二乘：b = (Phi* Phi)^{-1} Phi* X，这样对噪声更鲁棒。

4.5 内存溢出，SVD 跑不动

现象：n=10^6，m=500，SVD 直接 OOM。原因：full_matrices=False 已经省了一半，但 U 还是 n×m。解决：用随机 SVD 或者增量 SVD，只算前 r 个奇异值和向量。scikit-learn 的 randomized_svd 可以指定 n_components=r，内存占用和 n×r 成正比。

5. 进阶技巧：用 DMD 做频率排序与模态选择

5.1 按振幅和频率联合筛选模态

DMD 算完一堆模态，怎么挑出真正有用的？我一般会画一张振幅-频率散点图，横轴是频率，纵轴是振幅，每个点代表一个模态。主导模态通常振幅大、频率集中。下面这段代码把模态按振幅排序，并输出前 10 个的频率和增长率：

# 按振幅排序，输出前 10 个模态信息 idx = np.argsort(np.abs(b))[::-1] for i in idx[:10]: freq = np.imag(omega[i]) / (2 * np.pi) growth = np.real(omega[i]) print(f"模态 {i}: 频率 = {freq:.2f} Hz, 增长率 = {growth:.4f}, 振幅 = {np.abs(b[i]):.4f}")

逻辑说明：freq 是普通频率，单位 Hz；growth 是增长率，正数表示增长，负数表示衰减。参数上，如果你只关心稳定振荡，可以过滤掉 |growth| > 0.1 的模态，它们要么是瞬态要么是数值噪声。

5.2 用 sparsity-promoting DMD 自动选模态

如果模态太多，手动挑很累。sparsity-promoting DMD 在优化目标里加 L1 正则，让不重要的模态振幅自动归零。常见做法是交替方向乘子法，但实现起来代码量不小。我一般先用振幅阈值做粗筛，比如保留振幅大于最大振幅 1% 的模态，再人工看频率是否合理。这样比直接上稀疏优化快得多，而且可解释性好。

5.3 验证重构精度：三个必看指标

算完 DMD 一定要验证。我固定看三个指标：第一，重构场和原始场的相对误差，一般要求小于 5%；第二，主导频率和 FFT 峰值频率的偏差，小于 1% 算合格；第三，模态空间分布是否满足边界条件，比如壁面处速度是否为零。如果这三个都过了，结果基本可信。

# 重构所有时刻并计算相对误差 t_all = np.arange(X.shape[1]) * dt X_recon = np.zeros_like(X, dtype=complex) for k, t in enumerate(t_all): X_recon[:, k] = Phi_proj @ (b * np.exp(omega * t)) rel_err = np.linalg.norm(X - X_recon) / np.linalg.norm(X) print(f"重构相对误差: {rel_err:.4f}")

逻辑说明：X_recon 是复数矩阵，取实部才是物理场。rel_err 用 Frobenius 范数算，小于 0.05 说明模态选得不错。参数上，如果误差大，先增加 r，再检查 dt 和均值扣除。

5.4 我自己的习惯

每次跑完 DMD，我会把 omega 的实部和虚部画在一张复平面上，单位圆画出来，一眼就能看出哪些模态在增长、哪些在衰减、哪些在振荡。这个图比任何数字都直观。另外，我习惯把 dt、r、m、均值扣除标志写进结果文件头，过一个月再回来看，不用猜当时怎么设的。DMD 本身不复杂，复杂的是数据预处理和结果解读，把这两头做扎实，中间的计算就是几行代码的事。希望帮到你。

本文还有配套的精品资源，点击获取