做动力系统数据分析这些年，Koopman算子一直是我工具箱里最趁手的“解码器”。它能把混沌、非线性、周期振荡这些看起来分门别类的东西，统一放到观测函数空间里去理解。配合动态模态分解（DMD），我们可以直接从时间序列里估计Koopman算子的谱特性；再结合遍历论，你还能给这些谱一个统计力学层面的解释。这篇文章是我自己围绕“Koopman算子谱特性的计算研究”做的一次完整记录，代码环境是Matlab R2023b，核心路径只依赖ode45、svd和eig这些基础函数，没有额外工具箱。适合正在做动力系统辨识、非线性信号分解，或者准备把DMD写进论文和毕业设计的读者。
项目里我并不是简单地把DMD代码跑通就完事，而是把遍历论、DMD和谱特性当成一个整体来研究：先明确Koopman谱背后的数学含义，再在Matlab里实现稳定的计算流程，最后用Lorenz系统做验证。这样得到的结论既不是纯理论推演，也不是黑箱式的调参结果，而是能落到代码里的“计算力学”。下面按我的实际推进顺序来拆解。
1. 为什么Koopman、DMD和遍历论要放在一起看
1.1 Koopman算子的本质：非线性系统的线性“外壳”
Koopman理论的核心动作是“换一个视角看动力学”。考虑一个离散映射系统
[ x_{k+1} = F(x_k) ]
其中状态 (x) 一般是非线性的，可能是连续系统采样而来。Koopman算子 (K) 不直接作用在状态上，而是作用在观测函数 (g) 上：
[ (Kg)(x) = g(F(x)) ]
也就是说，给定一个观测函数 (g)，你先看它在当前状态的值，再通过系统演化一步，这个复合结果就是新函数在 (x) 处的值。这等于把“状态怎么演化”这个非线性问题，转换为“观测函数怎么演化”的线性问题。因为 (K) 对函数是线性的：
[ K(a g_1 + b g_2) = aK g_1 + bK g_2 ]
这个性质非常关键。非线性系统本身可能很难处理，但它的Koopman算子始终是线性算子。代价是这个算子生活在无穷维函数空间里，不能直接放进矩阵里算，必须截断和逼近。
如果Koopman算子有离散特征值和特征函数：
[ K\varphi = \lambda \varphi ]
那么沿着轨迹就有：
[ \varphi(x_k) = \lambda^k \varphi(x_0) ]
等价地说，存在一组“坐标变换”，能把原来的非线性动力学在这些坐标下变成完全解耦的线性系统。这就是很多人口中的“Koopman线性化”。我用一个比较生活化的类比：你没法直接改变一个暴躁用户的行为，但如果你只观察他每天发的消息摘要，这些摘要的变化模式反而可能是一条线性规律。
这个思想在应用上最直接的价值是：线性系统有成熟的特征值、模态、频谱分析工具，只要能把非线性系统“映射”过去，就能用线性工具研究非线性动力学。
1.2 DMD：从快照序列中拼出谱的近似
动态模态分解（DMD）本质上是数据驱动的Koopman算子近似。我们只有观测快照，不知道真实的映射 (F)，但我们可以构造两个快照矩阵：
[ X_1 = [x_0, x_1, \dots, x_{N-1}] ] [ X_2 = [x_1, x_2, \dots, x_N] ]
然后寻找一个矩阵 (A)，使得
[ X_2 \approx A X_1 ]
最优解是最小二乘意义下的：
[ A = X_2 X_1^{+} ]
这里的 (X_1^{+}) 是伪逆。(A) 的特征向量和特征值，就给了DMD模态和特征值。
DMD与Koopman算子的关系要说清楚：如果状态坐标函数本身落在Koopman算子不变子空间里，那么DMD得到的 (A) 就是Koopman算子在这个有限维子空间上的精确矩阵表示。对线性系统，这个条件自动满足；对非线性系统，它只是近似。更一般的做法是把观测函数也扩展进去，称为Extended DMD。但在实际工程里，标准DMD已经能给出非常有价值的定性判断，特别是谱特性的大致分布。
我自己的理解是：DMD是在有限维矩阵空间里“拼”一个Koopman算子的影子。它不一定能精确恢复所有谱点，但可以恢复那些在实际轨迹中能量占主导的谱成分。这正好是我们做谱特性计算时需要的。
1.3 遍历论：给谱一个统计解释
遍历论在这个框架里的角色有点像是“谱类型”的裁判。经典的Birkhoff遍历定理说，如果一个动力系统在不变测度下是遍历的，那么时间平均等于空间平均：
[ \lim_{N\to\infty} \frac{1}{N}\sum_{k=0}^{N-1} g(x_k) = \int g , d\mu ]
这个性质把“长时间观测”和“状态空间上的统计”联系了起来。
在Koopman算子谱理论里，谱的结构能反映动力学的混合程度。粗略地说：
- 纯点谱对应系统运动基本上是周期的或拟周期的；
- 连续谱往往对应混沌和混合，系统在相空间内“充分搅动”。
举个例子，Lorenz系统在经典参数下是混沌的、遍历的。它的Koopman谱就不只是离散的几条谱线，而是存在连续谱成分。这意味着你不管怎么增加DMD截断阶数，谱点都不会稳定地收敛到几个固定位置，而是会逐渐填充一条带状区域。这个经验和周期系统完全不一样。
所以，遍历论不是空泛的数学背景，它直接决定了我们对DMD特征值结果的解释方式。是几个离散点，还是一片连续带？这背后反映的是动力学的整体性质。
2. Matlab代码架构与关键实现
2.1 程序模块划分：别把所有代码塞进一个脚本里
这个项目我一开始也图省事，把所有步骤写在一个大脚本里，结果改参数的时候非常痛苦。后来拆成几个独立模块，清晰很多：
runAll.m：主脚本，负责调用所有步骤，设置全局参数；dataGen.m：生成连续/离散动力系统的快照数据，返回 (X_1, X_2)；dmdCore.m：实现核心DMD算法，输出特征值、模态、连续时间特征值；pseudoSpectrum.m：计算伪谱，用于观察连续谱；traversalTest.m：做遍历性检验，输出时间平均收敛曲线；plotSpectra.m：把所有谱相关的图统一画出来。
模块化带来的好处是：想换一个系统验证，只需要改dataGen.m；想换截断方法，只动dmdCore.m的输入。我在后面调Lorenz参数的时候，几乎只改主脚本里的sigma, rho, beta，其他文件都没碰。
2.2 核心DMD函数：每一行都要知道为什么
下面是我在项目里实际用的DMD核心函数，已经去掉了一堆调试输出，只保留主干：
function [Phi, lambda, omega, b, Atilde, U, S, V] = dmdCore(X1, X2, r, dt) % DMD核心函数 % 输入： % X1, X2 ：快照矩阵，维度为 nState × nSnap，相邻两列相差一个采样周期 % r ：SVD截断秩，传 [] 表示自动选择非零奇异值 % dt ：采样时间步长，用于把离散特征值转换为连续时间特征值 % 输出： % Phi ：DMD模态，维度为 nState × r % lambda ：DMD离散特征值 % omega ：连续时间特征值，omega = log(lambda) / dt % b ：初始条件在模态上的振幅 % Atilde ：降维后的近似矩阵 % 1. 对 X1 做经济SVD [U, S, V] = svd(X1, 'econ'); % 2. 截断秩 if nargin < 3 || isempty(r) s = diag(S); r = length(find(s > 1e-10)); end U_r = U(:, 1:r); S_r = S(1:r, 1:r); V_r = V(:, 1:r); % 3. 低维投影矩阵 Atilde = U_r' * X2 * V_r / S_r; % 4. 特征分解 [W, D] = eig(Atilde); lambda = diag(D); % 5. 连续时间特征值 if nargin >= 4 && ~isempty(dt) omega = log(lambda) / dt; else omega = []; end % 6. 精确DMD模态 Phi = X2 * V_r / S_r * W; % 7. 初始振幅 b = Phi \ X1(:, 1); end这段代码的重点在于第6步。前面特征分解得到的是低维矩阵 ( \tilde{A} ) 的特征向量 (W)，但如果直接把 (W) 当作DMD模态，它是在PCA降维后的坐标空间里的，不是原始状态空间。必须用：
[ \Phi = X_2 V_r S_r^{-1} W ]
映射回原始空间，这是exact DMD的标准做法。我自己第一次写的时候直接用了U_r * W，结果画出来的模态空间位置怎么都不对，后来查文献才发现这里有个坑。
另一个容易忽略的点是b = Phi \ X1(:,1)，这用左除解线性方程得到初始振幅。后续重构轨迹用：
[ x(t) \approx \Phi \cdot \mathrm{diag}(b) \cdot \lambda^{t/dt} ]
这个重构步骤特别适合验证DMD模型的拟合能力。
2.3 谱特性的数值计算要点：特征值、连续时间频率和伪谱
谱特性在Matlab里主要输出三类东西。
第一类是离散特征值lambda，直接画在复平面上。以单位圆为参考：
- 模小于1：模态衰减；
- 模等于1：模态等幅振荡；
- 模大于1：模态增长。
第二类是连续时间特征值omega。对采样步长 (dt)，有：
[ \omega = \frac{\log(\lambda)}{dt} ]
实部是增长率，虚部是角频率。画出来更接近我们熟悉的频率谱，尤其是在分析振荡系统时非常直观。注意取对数时要小心 (dt) 太大导致的相位混叠，后面一章会专门说。
第三类是伪谱（pseudospectrum）。为什么要算伪谱？因为在Lorenz这样的混沌系统里，Koopman谱带有连续谱成分，而特征值只是连续谱上的“采样点”。如果我们只画eig出来的几十个点，很难看出谱是离散还是连续。伪谱能告诉我们，在复平面上哪些区域的谱对扰动敏感，从而间接识别连续谱带。
我写了一个简单的伪谱计算函数：
function rho = pseudoSpectrum(Atilde, zGrid) % 计算 Atilde 的伪谱奇异值分布 % Atilde：低维DMD矩阵，维度为 r×r % zGrid ：复平面网格 r = size(Atilde, 1); rho = zeros(size(zGrid)); for i = 1:numel(zGrid) rho(i) = min(svd(zGrid(i) * eye(r) - Atilde)); end end对每个复数 (z)，计算 ((zI - \tilde{A})) 的最小奇异值。这个值越小，说明 (z) 越接近“谱”。如果特征值周围最小奇异值有很宽的下降带，说明谱对这个区域非常敏感，连续谱通常表现为带状区域。
伪谱的计算量不小，网格点多了会卡。我通常先粗略算 (200 \times 200) 的网格，找到主要亮暗区域后再局部加密。
3. 算例研究：Lorenz系统的谱特性与遍历性
3.1 算例设置与数据生成
我选的算例是经典Lorenz系统：
[ \begin{cases} \dot{x} = \sigma (y - x) \ \dot{y} = x (\rho - z) - y \ \dot{z} = x y - \beta z \end{cases} ]
参数取最常见的混沌状态：(\sigma=10, \rho=28, \beta=8/3)。
数据生成代码：
sigma = 10; rho = 28; beta = 8/3; dt = 0.02; tspan = 0:dt:120; f = @(t,x) [sigma*(x(2)-x(1)); x(1)*(rho-x(3)) - x(2); x(1)*x(2) - beta*x(3)]; [~, X] = ode45(f, tspan, [1; 1; 1]); % 丢弃前500个点，确保轨迹落在吸引子上 X = X(501:end, :)';为什么丢弃瞬态？Koopman算子和遍历论讨论的都是不变测度下的性质，如果轨迹还处在向吸引子收敛的瞬态过程，那么你的时间平均结果会受到初始条件影响，没法反映平稳状态。这个细节容易被忽略，但对谱计算影响很大。
采样步长 (dt=0.02) 对Lorenz系统来说已经足够细。如果 (dt) 太大，DMD会把高频折叠，特征值会歪曲；如果太小，相邻快照几乎线性相关，奇异值会非常接近零，数值稳定性变差。Lorenz系统的功率谱主峰在低频段，0.02已经能覆盖。
最后把状态矩阵转置成 (3 \times N)，每一列是一个快照。然后构造 (X_1) 和 (X_2)：
X1 = X(:, 1:end-1); X2 = X(:, 2:end);注意这里是连续时间系统，但DMD的映射是离散一步，(dt=0.02) 就是这个离散映射的步长。
3.2 DMD谱特性结果
调用核心函数：
r = 30; % 先固定一个较大的截断阶数观察趋势 [Phi, lambda, omega, b, Atilde] = dmdCore(X1, X2, r, dt);先说奇异值分布。对Lorenz数据做SVD后，奇异值衰减不快，不会像线性系统那样有一两个特别大的峰。我见过的结果里，前面20个左右奇异值占主导，之后慢慢下降。这说明系统本质上是无穷维的，DMD的有限维近似需要比较多模态。
画出特征值图：
figure; plot(real(lambda), imag(lambda), 'o'); hold on; theta = linspace(0, 2*pi, 200); plot(cos(theta), sin(theta), 'k--'); axis equal; grid on; xlabel('Re(\lambda)'); ylabel('Im(\lambda)');实际结果中，绝大多数特征值落在单位圆内，但离单位圆并不远。这说明DMD把Lorenz系统的模态解释为一组衰减振荡模态。如果你只看这些特征值，可能误判系统是稳定的，但Lorenz系统明明是有界的混沌吸引子。这是因为DMD的 ( \tilde{A} ) 是一个低维近似，有限维线性系统无法完全复现无穷维混沌的长期行为，它给出的只是局部短期的近似。
再看连续时间特征值omega，画在复平面实部-虚部坐标系：
figure; plot(real(omega), imag(omega), '.'); grid on; xlabel('Re(\omega)'); ylabel('Im(\omega)');会看到虚部对应不同的振荡频率，实部有负有正。正实部的模态在理论上会增长，但它们的振幅很小，对重构轨迹影响不大。这个细节很重要：不要看到正实部就觉得数值发散了，还要结合模态振幅b来看。
我还做了重构对比：用DMD模态重构一段轨迹，和原始轨迹对比。短时间窗口内，DMD能很好地跟踪Lorenz系统的瞬态；长时间窗口误差不断放大。这是因为混沌系统对初始条件敏感，任何模型预测都是有限时间有效。
3.3 遍历性检验与连续谱解释
做遍历性检验时，最简单的做法是看时间平均是否收敛。对Lorenz系统的 (x) 坐标，理论上由于系统关于 (x=0, y=0) 对称，长时平均应该趋于0。实际数值：
N = size(X, 2); M = floor(logspace(1, log10(N), 40)); for k = 1:length(M) avg(k) = mean(X(1, 1:M(k))); end semilogx(M, avg); xlabel('样本数'); ylabel('x 的时间平均');运行后能看到平均在振荡中逐渐衰减并趋近0。这就是Birkhoff遍历定理在数值上的体现。再有名的例子是 (z) 坐标的平均值，它应该收敛到一个与初始条件基本无关的常数。用不同初始条件跑多组，取长时间平均，会发现结果非常接近，说明系统在吸引子上的物理测度是唯一的。这就是遍历性的实际表现。
遍历性对谱计算的启示是什么？回到伪谱上。如果系统只有离散谱，那么特征值位置很稳定，增加或减少DMD截断阶数 (r)，主要特征值位置基本不变。但对Lorenz系统，我反复改变 (r) 值后发现，特征值在复平面上并非固定在几个点，而是随着 (r) 增大逐渐填充某个区域。这正是连续谱的特征。
我还用伪谱函数在单位圆附近算了一下：
z = exp(1i * linspace(0, 2*pi, 300)); rho = pseudoSpectrum(Atilde, z); plot(z, rho, '.');结果里，单位圆附近最小奇异值并不是只有一个尖峰，而是有一段连续的“低值区”。这说明该区域的谱对扰动敏感，典型的连续谱或近似连续谱特征。如果系统是周期振荡，伪谱会出现非常锐利的谷点，对应单一频率。这个对比是做谱特性研究时非常有力的工具。
4. 常见问题、坑点与调参经验
4.1 采样率与快照长度怎么定
很多第一次跑DMD的人问：采样率是不是越高越好？不是。采样率太高会带来两个问题：
- 相邻快照几乎线性相关，(X_1) 的奇异值有几个会非常接近0，伪逆计算不稳定；
- 内存和计算量大幅增加，SVD慢。
采样率太低的后果是频率混叠。连续特征值 (\omega = \log(\lambda) / dt) 的虚部范围受 Nyquist 限制，如果系统真实频率高于 (1/(2dt))，就会被折叠到低频区。
我个人的经验是：先做一个快速傅里叶分析，找到主能量所在的频率范围，采样率设为最高感兴趣频率的10到20倍。比如Lorenz系统主要能量集中在低频，取 (dt) 在0.01到0.05之间就够。如果做混沌系统，还要考虑Lyapunov时间尺度，总观测时间至少覆盖几个Lyapunov时间，否则DMD得到的模态统计意义不强。
快照长度方面，初始的瞬态必须丢弃。我通常在生成数据后去掉前10%-20%的点。
4.2 SVD截断值r：不要死记硬背公式
SVD截断 (r) 的选择没有万能答案。常见的经验方法是看能量占比：
s = diag(S); energy = cumsum(s.^2) / sum(s.^2); r = find(energy > 0.99, 1, 'first');这个方法对线性系统非常好用，但Lorenz这样的连续谱系统会给出很大的 (r)，因为能量衰减慢，0.99可能要到几百阶才满足。这时候你要结合问题背景去选。
我更常用的策略是画奇异值谱：
semilogy(s, 'o');找到曲线下降由陡变缓的拐点位置，取拐点附近的 (r)。对Lorenz，这个拐点大概在 (r=20\sim40) 之间。如果选太小，模态数量不够，重构误差迅速增大；选太大，后面若干模态幅值很小，主要是噪声拟合。
还有一个有用的小技巧：做两组DMD，一组 (r)，一组 (r+5)，比较重叠的主要特征值。如果位置稳定，说明这些模态是数据里的真实结构；如果到处乱跑，那基本就是噪声模态。
4.3 伪谱计算的工程细节
伪谱最怕算得慢。pseudoSpectrum里对每个网格点做一次svd，如果 (r=30)，300×300网格就是9万次30×30矩阵的SVD，Matlab可能要跑几十秒甚至几分钟。我的做法是两个优化：
第一，先粗后细。用 (80\times80) 网格找到感兴趣区域，再局部细化。 第二，不要把伪谱扩大到整个复平面，只算单位圆附近或者eig(lambda)结果周围的带状区域。
另外要注意伪谱里的奇异值变化范围很大，画图时建议用对数色标：
imagesc(real(zGrid), imag(zGrid), 20*log10(rho)); set(gca, 'YDir', 'normal');这样能同时看到深谷和缓坡，比线性色标清楚得多。
4.4 谱特性解读的三个常见误区
第一个误区是把DMD模态直接说成“Koopman特征函数”。DMD模态是通过状态坐标观测函数得到的近似结果，它是否逼近真实特征函数取决于状态坐标是否落在不变子空间里。对非线性系统，这个条件很难满足。我一般在报告里说“DMD模态是Koopman特征函数在某组观测下的投影”，这样更严谨。
第二个误区是把所有落在单位圆内的特征值都解释为“系统稳定”。单位圆内确实对应衰减，但DMD作为一个有限维近似，连续谱上的采样也可能会有人为的衰减。特别是在混沌系统里，很多特征值虚部对应真实的振荡频率，实部则被连续谱误差拉成负值。如果你要判断真实稳定性，必须结合原始系统能量约束，不能只看DMD特征值。
第三个误区是忽略观测函数选择的影响。标准DMD用状态坐标本身做观测，所以谱结果依赖坐标系。同一个物理系统，换一组坐标，DMD特征值会有变化。这也是Koopman框架里为什么会有Extended DMD，增加合适的观测函数能得到更接近真实谱的结果。对我们的标准DMD来说，结论要限定在“给定观测函数”的前提下。
最后分享一个小经验：谱特性计算最怕的不是代码写不出来，而是不知道怎么解读结果。我踩过几次坑之后养成一个习惯，每个算例都同时跑一个已知周期系统做对照。比如先跑一个线性振荡器，确认DMD能给出清晰的两个共轭特征值；再跑Lorenz，看到特征值分布和伪谱形态，两者一对比，连续谱和离散谱的区别就特别直观。这个对照习惯帮我避开了解读上的很多误判。
