做动力系统数据分析这些年Koopman算子一直是我工具箱里最趁手的“解码器”。它能把混沌、非线性、周期振荡这些看起来分门别类的东西统一放到观测函数空间里去理解。配合动态模态分解DMD我们可以直接从时间序列里估计Koopman算子的谱特性再结合遍历论你还能给这些谱一个统计力学层面的解释。这篇文章是我自己围绕“Koopman算子谱特性的计算研究”做的一次完整记录代码环境是Matlab R2023b核心路径只依赖ode45、svd和eig这些基础函数没有额外工具箱。适合正在做动力系统辨识、非线性信号分解或者准备把DMD写进论文和毕业设计的读者。项目里我并不是简单地把DMD代码跑通就完事而是把遍历论、DMD和谱特性当成一个整体来研究先明确Koopman谱背后的数学含义再在Matlab里实现稳定的计算流程最后用Lorenz系统做验证。这样得到的结论既不是纯理论推演也不是黑箱式的调参结果而是能落到代码里的“计算力学”。下面按我的实际推进顺序来拆解。1. 为什么Koopman、DMD和遍历论要放在一起看1.1 Koopman算子的本质非线性系统的线性“外壳”Koopman理论的核心动作是“换一个视角看动力学”。考虑一个离散映射系统[ x_{k1} 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_{k0}^{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} ]参数取最常见的混沌状态(\sigma10, \rho28, \beta8/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算子和遍历论讨论的都是不变测度下的性质如果轨迹还处在向吸引子收敛的瞬态过程那么你的时间平均结果会受到初始条件影响没法反映平稳状态。这个细节容易被忽略但对谱计算影响很大。采样步长 (dt0.02) 对Lorenz系统来说已经足够细。如果 (dt) 太大DMD会把高频折叠特征值会歪曲如果太小相邻快照几乎线性相关奇异值会非常接近零数值稳定性变差。Lorenz系统的功率谱主峰在低频段0.02已经能覆盖。最后把状态矩阵转置成 (3 \times N)每一列是一个快照。然后构造 (X_1) 和 (X_2)X1 X(:, 1:end-1); X2 X(:, 2:end);注意这里是连续时间系统但DMD的映射是离散一步(dt0.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) 坐标理论上由于系统关于 (x0, y0) 对称长时平均应该趋于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这个拐点大概在 (r20\sim40) 之间。如果选太小模态数量不够重构误差迅速增大选太大后面若干模态幅值很小主要是噪声拟合。还有一个有用的小技巧做两组DMD一组 (r)一组 (r5)比较重叠的主要特征值。如果位置稳定说明这些模态是数据里的真实结构如果到处乱跑那基本就是噪声模态。4.3 伪谱计算的工程细节伪谱最怕算得慢。pseudoSpectrum里对每个网格点做一次svd如果 (r30)300×300网格就是9万次30×30矩阵的SVDMatlab可能要跑几十秒甚至几分钟。我的做法是两个优化第一先粗后细。用 (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看到特征值分布和伪谱形态两者一对比连续谱和离散谱的区别就特别直观。这个对照习惯帮我避开了解读上的很多误判。