尧图网站开发 
1. 从非线性的痛点聊起为什么Koopman算子让人上头搞动力学分析的人十有八九都被非线性系统折磨过。以前做线性系统分析特征值一求稳定性、收敛速率、模态形状全出来了一套组合拳干净利落。但换到非线性系统局部线性化只能保一小块区域的近似精度想刻画全局演化行为就非常吃力。这也是动态模态分解DMD和Koopman算子这些年在流体力学、电力系统、神经科学、金融时序分析里这么火的原因——它们给了一套把非线性系统“强行”拉回线性框架的操作方法。先说我个人的结论Koopman算子本身是一个无穷维线性算子它把有限维非线性动力系统的演化操作提升到了一个函数空间上。在这个空间里系统的状态转移变成了线性作用于是特征值、特征函数、模态分解这些经典武器又能用了。代价也很直白无穷维嘛工程上全是近似遍历论则提供了一个观察长期行为和空间平均之间关系的理论透镜而DMD就是从离散数据快照里构造有限维近似矩阵的最常用手段。这篇内容适合的人很清楚正在做非线性动力学仿真的研究生、工程师想用Matlab快速验证Koopman/DMD方法效果的研究者以及对谱分析感兴趣但对无穷维算子有点发怵的入门选手。下面我按自己实际做这套计算的顺序来写先讲理论到底在干嘛再给Matlab实现框架最后列几个我踩过的坑。2. 理论层面必须想清楚的三个问题2.1 遍历论到底在说什么遍历论研究的是保测变换的长期统计行为。通俗点讲一个动力系统如果具有遍历性那么沿着几乎任何一条轨道的长时间时间平均会趋近于整个状态空间上的空间平均。这个性质在物理上非常重要因为很多宏观可观测量是从时间序列估计的如果系统不遍历你用一段轨迹去推断系统全局行为是会出大偏差的。从Koopman算子谱的视角看遍历性和谱密切相关。一个保测系统的Koopman算子如果是遍历的那么其谱在常数特征值1上具有某种单重性更强的混合性质则对应连续谱成分。这意味着计算Koopman算子的谱特性不仅能提取动力学模式还能从理论上判断系统是否具有遍历性、混合性甚至周期结构。简单说谱里藏着系统的“长期行为签名”。之前做过一个简单例子用了一个准周期系统两个不可约的频率比是无理数。理论上它的Koopman算子谱应该包含在单位圆上稠密的点谱数值上如果你用有限维DMD逼近特征值会密密麻麻地落在单位圆附近构成近似连续谱的形态。这就是遍历论和谱计算结合的一个直观体现。2.2 Koopman算子的定义与谱特性考虑离散时间动力系统[ x_{k1} F(x_k), \quad x \in \mathcal{M} ]观测函数是 ( g: \mathcal{M} \rightarrow \mathbb{C} )Koopman算子 (\mathcal{K}) 定义在观测函数空间上[ \mathcal{K}g g \circ F ]也就是说(\mathcal{K}g(x) g(F(x)))。这看起来就是把“取观测”和“系统演化一步”交换了顺序先在当前状态计算观测函数再对演化后的状态求值。(\mathcal{K}) 的线性性直接从函数加法和数乘的线性得出这一点和底层系统F是否非线性无关。谱特性中最核心的对象是特征值和特征函数。若存在非零函数 (\phi) 和非零复数 (\lambda) 满足[ \mathcal{K}\phi \lambda \phi ]则沿轨道有[ \phi(x_k) \lambda^k \phi(x_0) ]于是任何可以展开成特征函数线性组合的观测函数其演化都能表示为特征值幂次的线性叠加——整个长期预测就被谱分解拿下了。DMD做的就是从数据里估出这些 (\lambda) 以及对应的模态特征函数在采样点处的值。2.3 DMD与Koopman算子之间的桥DMD的基本输入是一组快照对 ((x_k, x_{k1}))也就是实测或仿真得到的离散状态序列。把快照排成矩阵X和X目标是找一个最优线性算子A使得[ X \approx A X ]最常用的解法是SVD截断[U, S, V] svd(X, econ); % 截断到rank r Ur U(:, 1:r); Sr S(1:r, 1:r); Vr V(:, 1:r); A_tilde Ur * X * Vr / Sr; [W, D] eig(A_tilde); Phi X * Vr / Sr * W; % 精确DMD模态 lambda diag(D);很多人问DMD和Koopman算子之间到底什么关系。严谨地说DMD在快照数量趋于无穷、且观测空间足够丰富时可以看作Koopman算子在某个有限维子空间上的投影近似。它得到的特征值对应Koopman谱的有限截断模态对应特征函数在该子空间上的表示。所以DMD不是Koopman算子的精确计算但它是目前工程上最实用的估计手段。还有个常见误解以为先有了Koopman理论再发明了DMD。实际历史顺序反过来的。DMD最早源于流体领域的模态分解需求后来研究者才意识到DMD其实是在有限维空间上近似Koopman算子。理解这一点很有用因为当你发现DMD结果不理想时很多理论工具能帮你分析原因——比如谱污染、连续谱近似、观测函数选取不充分等。3. Matlab代码实现从快照生成到谱特性输出3.1 生成动力学数据快照先拿一个经典非线性系统开刀Duffing振子% 参数 delta 0.3; alpha -1; beta 1; Gamma 0.5; omega 1.2; f (t, x) [x(2); ... -delta*x(2) - alpha*x(1) - beta*x(1).^3 Gamma*cos(omega*t)]; N 10000; dt 0.01; tspan linspace(0, N*dt, N1); x0 [1; 0]; [t, x] ode45(f, tspan, x0); X_raw x.;这里采样了10001个状态点每个状态是二维的位移和速度。对单条轨迹做DMD快照矩阵就是X X_raw(:, 1:end-1); Xp X_raw(:, 2:end);快照矩阵的大小是 (n_states, N-1)每列是一个状态向量。需要注意的是对于混沌系统或周期激励系统采样时长要足够长最好覆盖多个特征时间尺度否则DMD谱会有严重偏差。3.2 标准DMD完整实现下面这个函数是我平时常用的完整版包含了中心化处理、SVD截断、特征值和DMD模态计算function [lambda, Phi, A_tilde] dmd(X, Xp, r) % 中心化去掉时间平均分量 X_mean mean(X, 2); X_centered X - X_mean; Xp_centered Xp - X_mean; [U, S, V] svd(X_centered, econ); if nargin 3 || isempty(r) % 自动截断保留99%的能量 energy cumsum(diag(S).^2) / sum(diag(S).^2); r find(energy 0.99, 1); end Ur U(:, 1:r); Sr S(1:r, 1:r); Vr V(:, 1:r); A_tilde Ur * Xp_centered * Vr / Sr; [W, D] eig(A_tilde); lambda diag(D); Phi Xp_centered * Vr / Sr * W; end这段代码的关键步骤有三处。第一中心化。很多初学者忘了这一步直接把原始快照丢进去。结果是第一个特征值会接近1对应的模态几乎是“常数偏移量”这会污染其他模态的能量占比。实测下来中心化之后再算DMD模态更干净特别是系统的均值不为零时效果差距非常明显。第二SVD截断秩r的选择。固定r容易踩坑。r太小会丢失弱模态r太大会带进噪声。我用过两种自适应方法一是按能量占比上面代码里就是二是按奇异值的gap。能量占比法简单设置0.99基本够用奇异值gap法适合谱结构分明、奇异值有明显跳变的信号。第三精确DMD与标准DMD的区别。上面的写法算的是精确DMD的模态通过X和Vr、Sr构造Phi标准的投影DMD则是 ( \Phi U_r W )。两种方法在数据无噪声时几乎一致但有噪声时精确DMD对噪声更敏感投影DMD更稳定。如果数据质量不高建议用投影DMDPhi_proj Ur * W;3.3 谱特性与误差分析算完特征值后先看它们在复平面上的位置。单位圆对应中性稳定振荡/准周期圆内是衰减圆外是增长。对耗散系统物理上合理的特征值应该在单位圆内或恰好圆上跑出圆外且能量占比高的模态基本可以判定是数值伪影。% 特征值谱绘制 figure; plot(real(lambda), imag(lambda), bx, MarkerSize, 10); hold on; theta linspace(0, 2*pi, 200); plot(cos(theta), sin(theta), k--); axis equal; xlabel(Re(\lambda)); ylabel(Im(\lambda)); title(Koopman/DMD特征值谱);再算每个模态的能量占比这对应Koopman算子谱的强度分布% 模态幅度 alpha Phi \ X_centered(:, 1); energy abs(alpha).^2; energy energy / sum(energy); % 绘制谱强度分布 stem(1:length(energy), energy, filled); xlabel(模态序号); ylabel(归一化能量);误差分析这个环节容易被忽略但非常重要。我们可以用DMD模型预测未来若干步与真实演化比较% 预测 x_pred zeros(size(X_centered)); x_pred(:, 1) X_centered(:, 1); for k 1:N-1 x_pred(:, k1) Phi * diag(lambda.^k) * alpha; end err sqrt(sum((X_centered - x_pred).^2, 1)) ./ sqrt(sum(X_centered.^2, 1)); figure; semilogy(err); xlabel(时间索引); ylabel(相对误差); grid on;预测误差的走向非常直观地反映了DMD模型的好坏。如果误差在短时间内呈指数增长说明谱的截断或噪声抑制出了问题如果误差缓慢增长说明DMD抓住了主要动力学残差来自截断的高阶Koopman模式。3.4 遍历性与谱连续性的数值判断理论部分讲了遍历性对应Koopman算子谱的某些特征。数值上怎么判断一个系统是否遍历常见做法是比较时间平均和空间平均。空间平均需要对大量不同初始条件进行采样计算它们的长期平均然后与单个轨迹的时间平均对比% 时间平均 Ntraj 20; t_obs 2000; x_avg_time mean(X_raw(1, 1:t_obs), 2); % 空间平均多轨迹 x_avg_space 0; for i 1:Ntraj x0_i [randn; randn]; [~, x_i] ode45(f, tspan(1:t_obs), x0_i); x_avg_space x_avg_space mean(x_i(:, 1)); end x_avg_space x_avg_space / Ntraj; fprintf(时间平均: %f\n空间平均: %f\n, x_avg_time, x_avg_space);两者接近说明在该观测函数下系统表现出遍历性差距大则要怀疑系统是否有多个遍历分支多个吸引子或准周期环。Koopman谱在这时也会有对应的提示——如果特征值谱中除了0频率外还有很多集中分布的低频成分往往对应长时间尺度的慢变过程时间平均收敛就会非常慢。3.5 一个端到端的Matlab算例周期强迫Duffing系统把上面所有环节拼起来给出一个可以直接跑的完整脚本。这里我用了Gamma0.8的强强迫情形系统会呈现丰富的次谐波和准周期行为DMD谱会有多个显著模态% --- 系统定义 --- delta 0.3; alpha -1; beta 1; Gamma 0.8; omega 1.2; f (t, x) [x(2); ... -delta*x(2) - alpha*x(1) - beta*x(1).^3 Gamma*cos(omega*t)]; N 20000; dt 0.005; tspan linspace(0, N*dt, N1); [t, x] ode45(f, tspan, [1; 0]); % --- 准备快照 --- X_raw x.; X X_raw(:, 1:end-1); Xp X_raw(:, 2:end); X_mean mean(X, 2); X_c X - X_mean; Xp_c Xp - X_mean; % --- DMD --- r 30; [U, S, V] svd(X_c, econ); Ur U(:, 1:r); Sr S(1:r, 1:r); Vr V(:, 1:r); A_tilde Ur * Xp_c * Vr / Sr; [W, D] eig(A_tilde); lambda diag(D); Phi Ur * W; % 投影DMD模态更稳 % --- 模态排序 --- alpha Phi \ X_c(:, 1); [~, idx] sort(abs(alpha), descend); lambda_sorted lambda(idx); Phi_sorted Phi(:, idx); alpha_sorted alpha(idx); % --- 绘图频谱 --- freq abs(angle(lambda_sorted)) / (2*pi*dt); figure; stem(freq, abs(alpha_sorted).^2 / sum(abs(alpha_sorted).^2), filled); xlabel(频率 (Hz)); ylabel(归一化能量); title(DMD频谱);运行这段代码你会看到频谱在强迫频率0.19Hz左右处有主峰还有谐波分量。如果强迫幅度够大还会在低频段看到组合频率成分这些对应的特征值都在单位圆附近。数一下显著模态的数量基本就得到了系统约化模型的维度——这对后续做模型降阶很有用。4. 核心难点的实操拆解4.1 快照数量与采样间隔怎么选先给结论公式采样间隔要满足奈奎斯特条件至少每周期采两个点快照数量要能覆盖系统最慢的时间尺度一般至少几十个周期。但在Koopman/DMD框架下还有个容易被忽略的点快照矩阵的列数即样本数直接决定了SVD截断能保留多少模。列数太少连系统的基本模态数都凑不齐再牛的算法也白搭。我一般先做一次快速测试% 观察奇异值谱 [~, S, ~] svd(X_c, econ); semilogy(diag(S), o);如果奇异值没有明显衰减到底说明快照数量或者状态维数不够需要增加采样长度。如果奇异值衰减得很快几十个奇异值后就是噪声平台说明系统内在维度低可以用很小的r做模型。采样间隔也很讲究。间隔过大快照之间的关联弱DMD会误判高频模态间隔过小相邻快照几乎一样数值上X和X过于接近SVD会退化模态对噪声极其敏感。经验做法先用功率谱或自相关分析估计系统主频率再取主周期的大约20~50个采样点也就是每个周期采20到50次。4.2 截断秩r怎么定能量法、gap法和交叉验证前面提过能量法。实际操作中我推荐“能量gpare检验”双保险energy cumsum(diag(S).^2) / sum(diag(S).^2); r1 find(energy 0.999, 1); gap abs(diag(S(1:end-1, 1:end-1))) - abs(diag(S(2:end, 2:end))); r2 find(gap mean(gap) 2*std(gap), 1); r max(r1, r2); % 取两者较大值保证不丢主要模态gap法的逻辑是奇异值谱如果存在“悬崖”比如从第15个到第16个突然掉了一个量级那说明前15个是信号后面全是噪声取r15就足够了。能量法的问题是遇到噪声水平较高时噪声能量会推高总能量导致需要更多模态才能达到99%的能量占比反而引入了噪声模态。两者结合取较大值兼顾能量覆盖和噪声抑制。还有更严谨的方案训练/测试快照交叉验证。前70%快照用于DMD学习后30%用于预测比较选择预测误差最小的r。这种方式可靠但计算量稍大。对于离线分析完全没问题如果做实时估计就得小心了。4.3 系统在快照窗口内不平稳怎么办DMD隐含的假设是底层系统是时不变的自治系统或者至少窗口内近似时不变。如果系统参数慢变比如Duffing的delta随时间逐渐增大直接用DMD会把慢变过程当作一个模态提取出来产生误导。应对策略有几种。一是加窗DMD把长快照序列切分为滑动窗口每个窗口算一次DMD观察特征值轨迹的变化二是用扩展DMDEDMD把观测函数扩充为包含延迟或非线性基函数让Koopman近似在更大的函数空间上更准确三是如果知道参数的变化规律可以尝试做参数化Koopman算子但Matlab实现复杂度高通常只在研究阶段使用。工程上我最常用的是滑动窗口。代码不复杂window 1000; step 200; lambda_all []; for i 1:step:N-window Xw X(:, i:iwindow-1); Xwp X(:, i1:iwindow); [~, lambda_i, ~] dmd(Xw, Xwp, r); lambda_all [lambda_all; lambda_i.]; end然后绘制特征值随窗口变化的轨迹图可以直观看到哪些模态在漂移、哪些模态保持稳定。实时监测系统退化或参数漂移时这个方法非常好用。5. 常见问题与排查技巧实录5.1 DMD特征值和精确特征值对不上这是因为DMD的特征值实际上是Koopman算子特征值在有限维子空间上的投影而不是精确的谱。特别是当观测函数不能span整个特征函数空间时DMD特征值会发生偏移。解决办法是增加状态变量比如加入延迟嵌入或扩充观测函数集。延迟嵌入是实际中提升效果最明显的手段。将观测从 ( x_k ) 扩展成 ( [x_k; x_{k-1}; ...; x_{k-d}] )相当于用历史信息增强观测空间的表达能力d 5; X_delay zeros(d*size(X, 1), size(X, 2)-d); for i 1:d X_delay((i-1)*size(X, 1)1:i*size(X, 1), :) X(:, i:end-di-1); end延迟嵌入之后重新算DMD特征值通常会更接近真实谱尤其是对混沌系统效果立竿见影。5.2 所有特征值都接近1或0特征值全接近1通常是没做中心化处理。原始数据中包含常数偏移或慢变量DMD会把均值部分作为零频模态提取出来其他模态被压缩到频率0附近。处理方法是减去时间平均或者对数据进行高通滤波。特征值全接近0是另一个问题多见于快照矩阵中相邻列高度线性相关SVD后有效秩很小系统信息全被数值噪声掩盖。此时加大采样间隔减少列间相关性或者增加采样点数让信噪比提升。5.3 DMD预测发散但真实系统不发散预测发散最常见原因是截断秩r过大把数值噪声模态也保留了。这些模态的特征值模长大于1在预测步数增加后指数爆炸。解决方法是检查特征值模长分布剔除模长大于1的伪模态。还有一种隐蔽的情况真实系统确实收敛到吸引子但DMD预测在中间时刻超过物理边界。比如Lorenz系统状态变量不可能无限大但DMD模态线性叠加时可能产生正反馈。这种时候需要在DMD模型后接一个饱和非线性层或者在损失函数中加入物理约束。简单起见可以先尝试只保留最显著的几个模态降低模型复杂度过拟合。5.4 谱强度和预期不符谱强度分布受初始快照的影响很大。如果初始状态在吸引子附近但没充分进入吸引子前几个快照的瞬态成分会主导DMD模态幅度掩盖了长期行为的主导模态。我的做法是先抛弃前10%的快照让系统演化充分后再开始收集DMD数据。如果数据来自仿真直接延长烧录spin-up时间如果来自实验也要尽量等待系统进入稳态后再采样。此外模态幅度计算用的是伪逆 ( \alpha \Phi^ x_1 )当模态矩阵病态时这个计算会把噪声放大。可以改用正则化最小二乘alpha (Phi * Phi 1e-6 * eye(size(Phi, 2))) \ (Phi * X_c(:, 1));这个小改动往往能让谱强度分布更稳定。5.5 遍历性判断的数值陷阱遍历性判断很容易被有限时间窗口坑到。遍历理论中的极限是时间趋于无穷数值仿真永远是有限时间必然有误差。特别是对次谐波或准周期系统时间平均收敛极慢你可能需要非常长的轨迹才能让时间平均和空间平均接近。一种实用技巧是分窗口考察时间平均的收敛曲线% 移动平均 cum_mean cumsum(X_raw(1, :)) ./ (1:length(X_raw)); figure; semilogx(1:length(X_raw), abs(cum_mean - x_avg_space));如果对数-对数图上曲线斜率接近-1说明时间平均以1/T速率收敛这是遍历性的典型表现如果曲线出现平台或回升说明系统存在多个时间尺度的慢过程或者多吸引子共存这时需要更长的采样长度或者在多个初始条件下平均。6. 扩展方向从单轨迹DMD到Koopman计算的进阶工具记住一点标题里的“谱特性计算”不是算完特征值就结束了。真正能体现Koopman算子威力的是进一步构造可用的数值工具。6.1 扩展DMD变分自编码器与基函数网格EDMD的核心思路是把观测函数从状态向量提升到一组基函数的张成空间。例如对二维系统用高斯基函数% 网格中心 cx -2:0.4:2; cy -2:0.4:2; [CX, CY] meshgrid(cx, cy); centers [CX(:), CY(:)]; sigma 0.5; % 对每个快照计算基函数值 Psi zeros(size(centers, 2), size(X, 2)); for j 1:size(X, 2) d2 sum((X_c(:, j) - centers).^2, 1); Psi(:, j) exp(-d2 / (2*sigma^2)); end % 然后用 Psi 替换 X 做DMDEDMD对非线性系统的逼近能力远强于常规DMD代价是基函数数量和布局需要调参。高斯基函数中心选取可以通过K-means聚类来自动生成避免手动布点。6.2 库普曼特征函数的数值近似算出Koopman特征值后下一步自然是特征函数。DMD模态提供的是特征函数在采样点上的值但想要在状态空间任意位置计算特征函数需要插值或回归。一个简单有效的方法是用径向基函数回归% 用RBF拟合特征函数 rbf_phi (x, c, s) exp(-sum((x - c).^2, 1) / (2*s^2)); A_rbf zeros(size(centers, 2), size(centers, 2)); b zeros(size(centers, 2), 1); for i 1:size(centers, 2) for j 1:size(centers, 2) A_rbf(i, j) rbf_phi(centers(:, i), centers(:, j), sigma); end end b Phi(:, 1); % 取第一阶特征函数 coef A_rbf \ b; % 任意状态x处特征函数值 phi_val (x) sum(coef .* rbf_phi(x, centers, sigma));有了特征函数的解析代理函数就可以在相空间中画等值线直观看到特征函数的结构分析Koopman算子谱对应的动力学几何。6.3 应用数据驱动模型预测控制算完Koopman谱后最实在的应用是做模型预测控制MPC。把非线性系统提升到Koopman坐标后动力学近似线性MPC变成标准的线性二次规划。Matlab中可以用内置的fmincon或者第三方MPC工具箱实现。大致流程是先用DMD/EDMD训练一个Koopman模型在线测量得到当前提升状态然后在提升空间里求解有限时域最优控制问题把第一个控制量施加到系统。我跑过一个简单的Duffing镇定任务Koopman-MPC的效果比直接在非线性模型上用非线性MPC更快而且控制律更平滑——因为线性模型没有局部线性化导致的抖动。6.4 如果数据是时间序列而非完整状态很多时候你只有单个传感器的时间序列没有完整的系统状态。这时要用Hankel矩阵构造延迟坐标function H hankel_matrix(x, d) n length(x) - d 1; H zeros(d, n); for i 1:d H(i, :) x(i:ni-1); end endHankelDMD就是在Hankel矩阵上做DMD相当于把Takens嵌入思想和Koopman理论结合非常适合仅含单通道测量的场景。这个方法在故障诊断和生物信号分析中用得很广。在Matlab中如果不想自己写全套可以用系统辨识工具箱中的某些函数做辅助验证。不过说实话DMD和Koopman这套东西自己写完全够用核心代码量其实不大理解透彻比套工具箱重要得多。7. 最后分享几个实战心得第一永远先做中心化再谈DMD。很多人问为什么我的特征值和别人的对不上十有八九是均值处理的问题。Koopman算子作用于函数空间均值分量对应零频模态不处理会让整个谱偏移。第二跑Koopman相关计算不要迷信高阶截断。无穷维算子的数值近似永远有截断误差刻意追求特征值无限逼近真值往往事倍功半。关键是确定你关心的频率或时间尺度然后选择能捕获这些尺度的最小模态集。第三矩阵数据一定要检查数值条件数。Matlab的svd和eig很鲁棒但A_tilde这个矩阵的条件数如果太大超过1e8整个特征分解就不可信了。出现这种情况降低r或者改投影DMD都能缓解。第四把延迟嵌入当作默认选项。即使你认为状态已完整测量加几个延迟坐标往往能提升谱的稳定性对噪声时的DMD结果尤其明显。这套流程我在多个系统上跑过从简单的受迫摆到混沌Lorenz再到流体仿真数据核心步骤都是一样的快照准备、SVD截断、特征分解、谱可视化、模态分析。上手之后你会发现非线性系统在Koopman视角下其实没那么可怕关键是选对观测空间和截断策略。如果你刚接触这块建议先从Duffing振子或Van der Pol振荡器入手把上面脚本跑通再换系统参数观察特征值怎么变化。这个“换参数——看谱——调截断”的循环比读十篇理论文章都更能建立直觉。 
1. 从非线性的痛点聊起为什么Koopman算子让人上头搞动力学分析的人十有八九都被非线性系统折磨过。以前做线性系统分析特征值一求稳定性、收敛速率、模态形状全出来了一套组合拳干净利落。但换到非线性系统局部线性化只能保一小块区域的近似精度想刻画全局演化行为就非常吃力。这也是动态模态分解DMD和Koopman算子这些年在流体力学、电力系统、神经科学、金融时序分析里这么火的原因——它们给了一套把非线性系统“强行”拉回线性框架的操作方法。先说我个人的结论Koopman算子本身是一个无穷维线性算子它把有限维非线性动力系统的演化操作提升到了一个函数空间上。在这个空间里系统的状态转移变成了线性作用于是特征值、特征函数、模态分解这些经典武器又能用了。代价也很直白无穷维嘛工程上全是近似遍历论则提供了一个观察长期行为和空间平均之间关系的理论透镜而DMD就是从离散数据快照里构造有限维近似矩阵的最常用手段。这篇内容适合的人很清楚正在做非线性动力学仿真的研究生、工程师想用Matlab快速验证Koopman/DMD方法效果的研究者以及对谱分析感兴趣但对无穷维算子有点发怵的入门选手。下面我按自己实际做这套计算的顺序来写先讲理论到底在干嘛再给Matlab实现框架最后列几个我踩过的坑。2. 理论层面必须想清楚的三个问题2.1 遍历论到底在说什么遍历论研究的是保测变换的长期统计行为。通俗点讲一个动力系统如果具有遍历性那么沿着几乎任何一条轨道的长时间时间平均会趋近于整个状态空间上的空间平均。这个性质在物理上非常重要因为很多宏观可观测量是从时间序列估计的如果系统不遍历你用一段轨迹去推断系统全局行为是会出大偏差的。从Koopman算子谱的视角看遍历性和谱密切相关。一个保测系统的Koopman算子如果是遍历的那么其谱在常数特征值1上具有某种单重性更强的混合性质则对应连续谱成分。这意味着计算Koopman算子的谱特性不仅能提取动力学模式还能从理论上判断系统是否具有遍历性、混合性甚至周期结构。简单说谱里藏着系统的“长期行为签名”。之前做过一个简单例子用了一个准周期系统两个不可约的频率比是无理数。理论上它的Koopman算子谱应该包含在单位圆上稠密的点谱数值上如果你用有限维DMD逼近特征值会密密麻麻地落在单位圆附近构成近似连续谱的形态。这就是遍历论和谱计算结合的一个直观体现。2.2 Koopman算子的定义与谱特性考虑离散时间动力系统[ x_{k1} F(x_k), \quad x \in \mathcal{M} ]观测函数是 ( g: \mathcal{M} \rightarrow \mathbb{C} )Koopman算子 (\mathcal{K}) 定义在观测函数空间上[ \mathcal{K}g g \circ F ]也就是说(\mathcal{K}g(x) g(F(x)))。这看起来就是把“取观测”和“系统演化一步”交换了顺序先在当前状态计算观测函数再对演化后的状态求值。(\mathcal{K}) 的线性性直接从函数加法和数乘的线性得出这一点和底层系统F是否非线性无关。谱特性中最核心的对象是特征值和特征函数。若存在非零函数 (\phi) 和非零复数 (\lambda) 满足[ \mathcal{K}\phi \lambda \phi ]则沿轨道有[ \phi(x_k) \lambda^k \phi(x_0) ]于是任何可以展开成特征函数线性组合的观测函数其演化都能表示为特征值幂次的线性叠加——整个长期预测就被谱分解拿下了。DMD做的就是从数据里估出这些 (\lambda) 以及对应的模态特征函数在采样点处的值。2.3 DMD与Koopman算子之间的桥DMD的基本输入是一组快照对 ((x_k, x_{k1}))也就是实测或仿真得到的离散状态序列。把快照排成矩阵X和X目标是找一个最优线性算子A使得[ X \approx A X ]最常用的解法是SVD截断[U, S, V] svd(X, econ); % 截断到rank r Ur U(:, 1:r); Sr S(1:r, 1:r); Vr V(:, 1:r); A_tilde Ur * X * Vr / Sr; [W, D] eig(A_tilde); Phi X * Vr / Sr * W; % 精确DMD模态 lambda diag(D);很多人问DMD和Koopman算子之间到底什么关系。严谨地说DMD在快照数量趋于无穷、且观测空间足够丰富时可以看作Koopman算子在某个有限维子空间上的投影近似。它得到的特征值对应Koopman谱的有限截断模态对应特征函数在该子空间上的表示。所以DMD不是Koopman算子的精确计算但它是目前工程上最实用的估计手段。还有个常见误解以为先有了Koopman理论再发明了DMD。实际历史顺序反过来的。DMD最早源于流体领域的模态分解需求后来研究者才意识到DMD其实是在有限维空间上近似Koopman算子。理解这一点很有用因为当你发现DMD结果不理想时很多理论工具能帮你分析原因——比如谱污染、连续谱近似、观测函数选取不充分等。3. Matlab代码实现从快照生成到谱特性输出3.1 生成动力学数据快照先拿一个经典非线性系统开刀Duffing振子% 参数 delta 0.3; alpha -1; beta 1; Gamma 0.5; omega 1.2; f (t, x) [x(2); ... -delta*x(2) - alpha*x(1) - beta*x(1).^3 Gamma*cos(omega*t)]; N 10000; dt 0.01; tspan linspace(0, N*dt, N1); x0 [1; 0]; [t, x] ode45(f, tspan, x0); X_raw x.;这里采样了10001个状态点每个状态是二维的位移和速度。对单条轨迹做DMD快照矩阵就是X X_raw(:, 1:end-1); Xp X_raw(:, 2:end);快照矩阵的大小是 (n_states, N-1)每列是一个状态向量。需要注意的是对于混沌系统或周期激励系统采样时长要足够长最好覆盖多个特征时间尺度否则DMD谱会有严重偏差。3.2 标准DMD完整实现下面这个函数是我平时常用的完整版包含了中心化处理、SVD截断、特征值和DMD模态计算function [lambda, Phi, A_tilde] dmd(X, Xp, r) % 中心化去掉时间平均分量 X_mean mean(X, 2); X_centered X - X_mean; Xp_centered Xp - X_mean; [U, S, V] svd(X_centered, econ); if nargin 3 || isempty(r) % 自动截断保留99%的能量 energy cumsum(diag(S).^2) / sum(diag(S).^2); r find(energy 0.99, 1); end Ur U(:, 1:r); Sr S(1:r, 1:r); Vr V(:, 1:r); A_tilde Ur * Xp_centered * Vr / Sr; [W, D] eig(A_tilde); lambda diag(D); Phi Xp_centered * Vr / Sr * W; end这段代码的关键步骤有三处。第一中心化。很多初学者忘了这一步直接把原始快照丢进去。结果是第一个特征值会接近1对应的模态几乎是“常数偏移量”这会污染其他模态的能量占比。实测下来中心化之后再算DMD模态更干净特别是系统的均值不为零时效果差距非常明显。第二SVD截断秩r的选择。固定r容易踩坑。r太小会丢失弱模态r太大会带进噪声。我用过两种自适应方法一是按能量占比上面代码里就是二是按奇异值的gap。能量占比法简单设置0.99基本够用奇异值gap法适合谱结构分明、奇异值有明显跳变的信号。第三精确DMD与标准DMD的区别。上面的写法算的是精确DMD的模态通过X和Vr、Sr构造Phi标准的投影DMD则是 ( \Phi U_r W )。两种方法在数据无噪声时几乎一致但有噪声时精确DMD对噪声更敏感投影DMD更稳定。如果数据质量不高建议用投影DMDPhi_proj Ur * W;3.3 谱特性与误差分析算完特征值后先看它们在复平面上的位置。单位圆对应中性稳定振荡/准周期圆内是衰减圆外是增长。对耗散系统物理上合理的特征值应该在单位圆内或恰好圆上跑出圆外且能量占比高的模态基本可以判定是数值伪影。% 特征值谱绘制 figure; plot(real(lambda), imag(lambda), bx, MarkerSize, 10); hold on; theta linspace(0, 2*pi, 200); plot(cos(theta), sin(theta), k--); axis equal; xlabel(Re(\lambda)); ylabel(Im(\lambda)); title(Koopman/DMD特征值谱);再算每个模态的能量占比这对应Koopman算子谱的强度分布% 模态幅度 alpha Phi \ X_centered(:, 1); energy abs(alpha).^2; energy energy / sum(energy); % 绘制谱强度分布 stem(1:length(energy), energy, filled); xlabel(模态序号); ylabel(归一化能量);误差分析这个环节容易被忽略但非常重要。我们可以用DMD模型预测未来若干步与真实演化比较% 预测 x_pred zeros(size(X_centered)); x_pred(:, 1) X_centered(:, 1); for k 1:N-1 x_pred(:, k1) Phi * diag(lambda.^k) * alpha; end err sqrt(sum((X_centered - x_pred).^2, 1)) ./ sqrt(sum(X_centered.^2, 1)); figure; semilogy(err); xlabel(时间索引); ylabel(相对误差); grid on;预测误差的走向非常直观地反映了DMD模型的好坏。如果误差在短时间内呈指数增长说明谱的截断或噪声抑制出了问题如果误差缓慢增长说明DMD抓住了主要动力学残差来自截断的高阶Koopman模式。3.4 遍历性与谱连续性的数值判断理论部分讲了遍历性对应Koopman算子谱的某些特征。数值上怎么判断一个系统是否遍历常见做法是比较时间平均和空间平均。空间平均需要对大量不同初始条件进行采样计算它们的长期平均然后与单个轨迹的时间平均对比% 时间平均 Ntraj 20; t_obs 2000; x_avg_time mean(X_raw(1, 1:t_obs), 2); % 空间平均多轨迹 x_avg_space 0; for i 1:Ntraj x0_i [randn; randn]; [~, x_i] ode45(f, tspan(1:t_obs), x0_i); x_avg_space x_avg_space mean(x_i(:, 1)); end x_avg_space x_avg_space / Ntraj; fprintf(时间平均: %f\n空间平均: %f\n, x_avg_time, x_avg_space);两者接近说明在该观测函数下系统表现出遍历性差距大则要怀疑系统是否有多个遍历分支多个吸引子或准周期环。Koopman谱在这时也会有对应的提示——如果特征值谱中除了0频率外还有很多集中分布的低频成分往往对应长时间尺度的慢变过程时间平均收敛就会非常慢。3.5 一个端到端的Matlab算例周期强迫Duffing系统把上面所有环节拼起来给出一个可以直接跑的完整脚本。这里我用了Gamma0.8的强强迫情形系统会呈现丰富的次谐波和准周期行为DMD谱会有多个显著模态% --- 系统定义 --- delta 0.3; alpha -1; beta 1; Gamma 0.8; omega 1.2; f (t, x) [x(2); ... -delta*x(2) - alpha*x(1) - beta*x(1).^3 Gamma*cos(omega*t)]; N 20000; dt 0.005; tspan linspace(0, N*dt, N1); [t, x] ode45(f, tspan, [1; 0]); % --- 准备快照 --- X_raw x.; X X_raw(:, 1:end-1); Xp X_raw(:, 2:end); X_mean mean(X, 2); X_c X - X_mean; Xp_c Xp - X_mean; % --- DMD --- r 30; [U, S, V] svd(X_c, econ); Ur U(:, 1:r); Sr S(1:r, 1:r); Vr V(:, 1:r); A_tilde Ur * Xp_c * Vr / Sr; [W, D] eig(A_tilde); lambda diag(D); Phi Ur * W; % 投影DMD模态更稳 % --- 模态排序 --- alpha Phi \ X_c(:, 1); [~, idx] sort(abs(alpha), descend); lambda_sorted lambda(idx); Phi_sorted Phi(:, idx); alpha_sorted alpha(idx); % --- 绘图频谱 --- freq abs(angle(lambda_sorted)) / (2*pi*dt); figure; stem(freq, abs(alpha_sorted).^2 / sum(abs(alpha_sorted).^2), filled); xlabel(频率 (Hz)); ylabel(归一化能量); title(DMD频谱);运行这段代码你会看到频谱在强迫频率0.19Hz左右处有主峰还有谐波分量。如果强迫幅度够大还会在低频段看到组合频率成分这些对应的特征值都在单位圆附近。数一下显著模态的数量基本就得到了系统约化模型的维度——这对后续做模型降阶很有用。4. 核心难点的实操拆解4.1 快照数量与采样间隔怎么选先给结论公式采样间隔要满足奈奎斯特条件至少每周期采两个点快照数量要能覆盖系统最慢的时间尺度一般至少几十个周期。但在Koopman/DMD框架下还有个容易被忽略的点快照矩阵的列数即样本数直接决定了SVD截断能保留多少模。列数太少连系统的基本模态数都凑不齐再牛的算法也白搭。我一般先做一次快速测试% 观察奇异值谱 [~, S, ~] svd(X_c, econ); semilogy(diag(S), o);如果奇异值没有明显衰减到底说明快照数量或者状态维数不够需要增加采样长度。如果奇异值衰减得很快几十个奇异值后就是噪声平台说明系统内在维度低可以用很小的r做模型。采样间隔也很讲究。间隔过大快照之间的关联弱DMD会误判高频模态间隔过小相邻快照几乎一样数值上X和X过于接近SVD会退化模态对噪声极其敏感。经验做法先用功率谱或自相关分析估计系统主频率再取主周期的大约20~50个采样点也就是每个周期采20到50次。4.2 截断秩r怎么定能量法、gap法和交叉验证前面提过能量法。实际操作中我推荐“能量gpare检验”双保险energy cumsum(diag(S).^2) / sum(diag(S).^2); r1 find(energy 0.999, 1); gap abs(diag(S(1:end-1, 1:end-1))) - abs(diag(S(2:end, 2:end))); r2 find(gap mean(gap) 2*std(gap), 1); r max(r1, r2); % 取两者较大值保证不丢主要模态gap法的逻辑是奇异值谱如果存在“悬崖”比如从第15个到第16个突然掉了一个量级那说明前15个是信号后面全是噪声取r15就足够了。能量法的问题是遇到噪声水平较高时噪声能量会推高总能量导致需要更多模态才能达到99%的能量占比反而引入了噪声模态。两者结合取较大值兼顾能量覆盖和噪声抑制。还有更严谨的方案训练/测试快照交叉验证。前70%快照用于DMD学习后30%用于预测比较选择预测误差最小的r。这种方式可靠但计算量稍大。对于离线分析完全没问题如果做实时估计就得小心了。4.3 系统在快照窗口内不平稳怎么办DMD隐含的假设是底层系统是时不变的自治系统或者至少窗口内近似时不变。如果系统参数慢变比如Duffing的delta随时间逐渐增大直接用DMD会把慢变过程当作一个模态提取出来产生误导。应对策略有几种。一是加窗DMD把长快照序列切分为滑动窗口每个窗口算一次DMD观察特征值轨迹的变化二是用扩展DMDEDMD把观测函数扩充为包含延迟或非线性基函数让Koopman近似在更大的函数空间上更准确三是如果知道参数的变化规律可以尝试做参数化Koopman算子但Matlab实现复杂度高通常只在研究阶段使用。工程上我最常用的是滑动窗口。代码不复杂window 1000; step 200; lambda_all []; for i 1:step:N-window Xw X(:, i:iwindow-1); Xwp X(:, i1:iwindow); [~, lambda_i, ~] dmd(Xw, Xwp, r); lambda_all [lambda_all; lambda_i.]; end然后绘制特征值随窗口变化的轨迹图可以直观看到哪些模态在漂移、哪些模态保持稳定。实时监测系统退化或参数漂移时这个方法非常好用。5. 常见问题与排查技巧实录5.1 DMD特征值和精确特征值对不上这是因为DMD的特征值实际上是Koopman算子特征值在有限维子空间上的投影而不是精确的谱。特别是当观测函数不能span整个特征函数空间时DMD特征值会发生偏移。解决办法是增加状态变量比如加入延迟嵌入或扩充观测函数集。延迟嵌入是实际中提升效果最明显的手段。将观测从 ( x_k ) 扩展成 ( [x_k; x_{k-1}; ...; x_{k-d}] )相当于用历史信息增强观测空间的表达能力d 5; X_delay zeros(d*size(X, 1), size(X, 2)-d); for i 1:d X_delay((i-1)*size(X, 1)1:i*size(X, 1), :) X(:, i:end-di-1); end延迟嵌入之后重新算DMD特征值通常会更接近真实谱尤其是对混沌系统效果立竿见影。5.2 所有特征值都接近1或0特征值全接近1通常是没做中心化处理。原始数据中包含常数偏移或慢变量DMD会把均值部分作为零频模态提取出来其他模态被压缩到频率0附近。处理方法是减去时间平均或者对数据进行高通滤波。特征值全接近0是另一个问题多见于快照矩阵中相邻列高度线性相关SVD后有效秩很小系统信息全被数值噪声掩盖。此时加大采样间隔减少列间相关性或者增加采样点数让信噪比提升。5.3 DMD预测发散但真实系统不发散预测发散最常见原因是截断秩r过大把数值噪声模态也保留了。这些模态的特征值模长大于1在预测步数增加后指数爆炸。解决方法是检查特征值模长分布剔除模长大于1的伪模态。还有一种隐蔽的情况真实系统确实收敛到吸引子但DMD预测在中间时刻超过物理边界。比如Lorenz系统状态变量不可能无限大但DMD模态线性叠加时可能产生正反馈。这种时候需要在DMD模型后接一个饱和非线性层或者在损失函数中加入物理约束。简单起见可以先尝试只保留最显著的几个模态降低模型复杂度过拟合。5.4 谱强度和预期不符谱强度分布受初始快照的影响很大。如果初始状态在吸引子附近但没充分进入吸引子前几个快照的瞬态成分会主导DMD模态幅度掩盖了长期行为的主导模态。我的做法是先抛弃前10%的快照让系统演化充分后再开始收集DMD数据。如果数据来自仿真直接延长烧录spin-up时间如果来自实验也要尽量等待系统进入稳态后再采样。此外模态幅度计算用的是伪逆 ( \alpha \Phi^ x_1 )当模态矩阵病态时这个计算会把噪声放大。可以改用正则化最小二乘alpha (Phi * Phi 1e-6 * eye(size(Phi, 2))) \ (Phi * X_c(:, 1));这个小改动往往能让谱强度分布更稳定。5.5 遍历性判断的数值陷阱遍历性判断很容易被有限时间窗口坑到。遍历理论中的极限是时间趋于无穷数值仿真永远是有限时间必然有误差。特别是对次谐波或准周期系统时间平均收敛极慢你可能需要非常长的轨迹才能让时间平均和空间平均接近。一种实用技巧是分窗口考察时间平均的收敛曲线% 移动平均 cum_mean cumsum(X_raw(1, :)) ./ (1:length(X_raw)); figure; semilogx(1:length(X_raw), abs(cum_mean - x_avg_space));如果对数-对数图上曲线斜率接近-1说明时间平均以1/T速率收敛这是遍历性的典型表现如果曲线出现平台或回升说明系统存在多个时间尺度的慢过程或者多吸引子共存这时需要更长的采样长度或者在多个初始条件下平均。6. 扩展方向从单轨迹DMD到Koopman计算的进阶工具记住一点标题里的“谱特性计算”不是算完特征值就结束了。真正能体现Koopman算子威力的是进一步构造可用的数值工具。6.1 扩展DMD变分自编码器与基函数网格EDMD的核心思路是把观测函数从状态向量提升到一组基函数的张成空间。例如对二维系统用高斯基函数% 网格中心 cx -2:0.4:2; cy -2:0.4:2; [CX, CY] meshgrid(cx, cy); centers [CX(:), CY(:)]; sigma 0.5; % 对每个快照计算基函数值 Psi zeros(size(centers, 2), size(X, 2)); for j 1:size(X, 2) d2 sum((X_c(:, j) - centers).^2, 1); Psi(:, j) exp(-d2 / (2*sigma^2)); end % 然后用 Psi 替换 X 做DMDEDMD对非线性系统的逼近能力远强于常规DMD代价是基函数数量和布局需要调参。高斯基函数中心选取可以通过K-means聚类来自动生成避免手动布点。6.2 库普曼特征函数的数值近似算出Koopman特征值后下一步自然是特征函数。DMD模态提供的是特征函数在采样点上的值但想要在状态空间任意位置计算特征函数需要插值或回归。一个简单有效的方法是用径向基函数回归% 用RBF拟合特征函数 rbf_phi (x, c, s) exp(-sum((x - c).^2, 1) / (2*s^2)); A_rbf zeros(size(centers, 2), size(centers, 2)); b zeros(size(centers, 2), 1); for i 1:size(centers, 2) for j 1:size(centers, 2) A_rbf(i, j) rbf_phi(centers(:, i), centers(:, j), sigma); end end b Phi(:, 1); % 取第一阶特征函数 coef A_rbf \ b; % 任意状态x处特征函数值 phi_val (x) sum(coef .* rbf_phi(x, centers, sigma));有了特征函数的解析代理函数就可以在相空间中画等值线直观看到特征函数的结构分析Koopman算子谱对应的动力学几何。6.3 应用数据驱动模型预测控制算完Koopman谱后最实在的应用是做模型预测控制MPC。把非线性系统提升到Koopman坐标后动力学近似线性MPC变成标准的线性二次规划。Matlab中可以用内置的fmincon或者第三方MPC工具箱实现。大致流程是先用DMD/EDMD训练一个Koopman模型在线测量得到当前提升状态然后在提升空间里求解有限时域最优控制问题把第一个控制量施加到系统。我跑过一个简单的Duffing镇定任务Koopman-MPC的效果比直接在非线性模型上用非线性MPC更快而且控制律更平滑——因为线性模型没有局部线性化导致的抖动。6.4 如果数据是时间序列而非完整状态很多时候你只有单个传感器的时间序列没有完整的系统状态。这时要用Hankel矩阵构造延迟坐标function H hankel_matrix(x, d) n length(x) - d 1; H zeros(d, n); for i 1:d H(i, :) x(i:ni-1); end endHankelDMD就是在Hankel矩阵上做DMD相当于把Takens嵌入思想和Koopman理论结合非常适合仅含单通道测量的场景。这个方法在故障诊断和生物信号分析中用得很广。在Matlab中如果不想自己写全套可以用系统辨识工具箱中的某些函数做辅助验证。不过说实话DMD和Koopman这套东西自己写完全够用核心代码量其实不大理解透彻比套工具箱重要得多。7. 最后分享几个实战心得第一永远先做中心化再谈DMD。很多人问为什么我的特征值和别人的对不上十有八九是均值处理的问题。Koopman算子作用于函数空间均值分量对应零频模态不处理会让整个谱偏移。第二跑Koopman相关计算不要迷信高阶截断。无穷维算子的数值近似永远有截断误差刻意追求特征值无限逼近真值往往事倍功半。关键是确定你关心的频率或时间尺度然后选择能捕获这些尺度的最小模态集。第三矩阵数据一定要检查数值条件数。Matlab的svd和eig很鲁棒但A_tilde这个矩阵的条件数如果太大超过1e8整个特征分解就不可信了。出现这种情况降低r或者改投影DMD都能缓解。第四把延迟嵌入当作默认选项。即使你认为状态已完整测量加几个延迟坐标往往能提升谱的稳定性对噪声时的DMD结果尤其明显。这套流程我在多个系统上跑过从简单的受迫摆到混沌Lorenz再到流体仿真数据核心步骤都是一样的快照准备、SVD截断、特征分解、谱可视化、模态分析。上手之后你会发现非线性系统在Koopman视角下其实没那么可怕关键是选对观测空间和截断策略。如果你刚接触这块建议先从Duffing振子或Van der Pol振荡器入手把上面脚本跑通再换系统参数观察特征值怎么变化。这个“换参数——看谱——调截断”的循环比读十篇理论文章都更能建立直觉。