做参数辨识这几年最头疼的从来不是模型有多复杂而是手里的轨迹数据全是噪声。用diff或者gradient去直接数值微分求状态导数结果基本是灾难——一阶差分会把测量噪声放大好几个量级后面的回归怎么做都是歪的。我最近在工程里试了一套基于分数阶占据核的方案用 MATLAB 从非线性动力学系统的轨迹数据里直接逼近状态导数并完成参数辨识实测下来比传统数值微分稳健很多尤其在信噪比不高的时候优势非常明显。这套方法的核心思路很朴素既然微分会把噪声放大那就别微分改用积分。把动力学方程两遍同时做分数阶积分再用占据核去构造辨识方程既绕开了显式求导又能把历史轨迹信息以幂率权重融进回归最后用稀疏回归从字典里把物理参数挑出来。整个过程不依赖任何工具箱Grunwald-Letnikov 分数阶积分自己写也就几十行。这篇文章我会把原理、代码、试验结果和踩过的坑从头到尾讲清楚适合正在做系统辨识、非线性动力学建模、状态估计这类课题的研究生和工程师参考。1. 项目思路与整体方案选型1.1 参数辨识里那个绕不开的痛点参数辨识的基本问题可以这样描述已知系统的大致模型结构比如一个 Duffing 振子满足 ẍ δẋ αx βx³ γcos(ωt)但我们不知道 α、β、δ、γ 具体是多少。现在测量到一组状态轨迹 x(t)希望反推出这几个系数。最直接的思路是把状态空间写出来令 x₁ xx₂ ẋ于是 ẋ₁ x₂ẋ₂ -αx₁ - βx₁³ - δx₂ γcos(ωt)。如果我能得到 ẋ₁、ẋ₂ 的数值那么对方程做线性回归就能得到参数。问题就出在这个“如果”上。实际测量得到的只是 x₁(t) 和 x₂(t) 的离散采样值要拿它们去算导数只能数值微分。而数值微分是个典型的不适定问题采样间隔 h 越小差分结果对噪声越敏感误差大致按 O(σ/h) 增长。我用一组仿真数据验证过这一点。给 Duffing 振子的轨迹加上 2% 的高斯白噪声用中心差分去算 ẋ₂结果导数的信噪比直接崩到没法看。后续无论用什么回归算法估计出来的阻尼项 δ 和三次刚度项 β 都偏差很大。这个问题不是回归算法的问题是数据准备阶段就已经把信息毁掉了。所以做辨识的第一步不是选多花哨的求解器而是要想办法绕开数值微分。1.2 占据核是怎么绕开导数的占据核的核心思想是不再去追求每一个时刻的导数数值而是把整条轨迹在状态空间里的“行为模式”提取出来作为辨识的输入特征。可以这样理解一条动力学轨迹在状态空间里会走过一系列位置有的区域停留时间长有的区域一闪而过。这种“停留时间分布”本身就携带了动力学参数的信息。占据核做的就是通过一个核函数把这种分布变成可计算的特征表达。比如用高斯核 k(x, ξ) exp(-‖x - ξ‖² / (2σ²))固定一组中心点 ξ_m然后沿着轨迹做积分K_m(t) ∫₀ᵗ k(x(τ), ξ_m) dτ得到的 K_m(t) 就是轨迹在第 m 个中心附近的加权停留时间累计曲线。一条阻尼大的轨迹和一条阻尼小的轨迹在这组曲线上会有明显差异而参数辨识就是从这个差异里把系数反推出来。更关键的是这个积分操作对噪声有天生的平滑作用。数值微分是高频放大器积分是低频滤波器。同样一条带噪轨迹积分之后的特征曲线信噪比要高出很多。把占据核和分数阶积分结合起来相当于在积分平滑的基础上又加了一重“记忆权重”对历史误差的累积也做了衰减控制这就是“分数阶占据核”的优势来源。1.3 方案选型为什么用 MATLAB 和物理字典这套方案我没有用 Python主要原因是 MATLAB 在矩阵运算、ode45仿真和绘图交互上太顺手了实验阶段改参数非常快。另外 MATLAB 的向量化写法能让占据核矩阵的构造比 Python 循环快不少尤其在做参数扫描的时候这个效率差距很实在。在字典选择上我做了两条路径的对比。一条是物理字典也就是直接用 x₁、x₂、x₁³、cos(ωt) 这类带物理含义的基函数做积分辨识好处是系数直接对应物理参数辨识结果可解释性极强。另一条是核字典用高斯核中心点做基好处是表达柔性更强适合模型结构完全未知的场景但系数没有物理含义只能用于重构状态导数不能直接读出系统参数。我的建议是如果你知道模型的大致结构优先用物理字典一步到位把参数辨识出来如果你对模型结构一无所知先用核字典做状态导数的重构再去配合稀疏回归做结构发现。这篇文章的主流程用物理字典完成参数辨识占据核主要用在构造稳定的积分特征上后面我也会给出核字典的扩展代码。2. 分数阶占据核的原理与关键推导2.1 从普通积分方程到占据核特征先把数学形式写清楚。考虑一个自治加外激励的非线性系统 ẋ f(x, t)其中 f 可以写成一组基函数的线性组合 f(x, t) Σⱼ θⱼ ψⱼ(x, t)。两边从 0 到 t 积分得到x(t) - x(0) Σⱼ θⱼ ∫₀ᵗ ψⱼ(x(τ), τ) dτ这是一个关于 θ 的线性方程右边积分里的每一项都可以从轨迹数据直接数值算出。这样就不需要任何导数测量只需要对轨迹做积分。这就是 SINDy 方法中“积分形式辨识”的构造思路。但普通积分有一个问题它对整个时间段的历史信息一视同仁早期的误差和近期的误差权重相同。如果早期数据受初值偏差、瞬态扰动影响较大普通积分会把这段误差一直保留下来。分数阶积分的出现就是为了给历史信息加一个随时间衰减的权重让辨识结果对瞬态和突发噪声更稳健。占据核在这个框架里的角色是多一重特征变换。用高斯核构造字典 ψₘ(x) k(x, ξₘ)对方程两边做同样处理就得到占据核方程。它适合模型结构完全未知时先把状态导数空间“张”出来。前面提到的物理字典则可以直接套用公式里的 ψⱼ两种方式在后面的 MATLAB 实现里共用同一个积分框架。2.2 分数阶积分为什么在这里管用分数阶微积分里的 Riemann-Liouville 积分定义是I^r f(t) (1/Γ(r)) ∫₀ᵗ (t-τ)^(r-1) f(τ) dτ其中 r 是任意正实数。当 r 1 时它退化成普通积分。当 r 1 时核函数 (t-τ)^(r-1) 在 τ 接近 t 时更大在 τ 远离 t 时按幂律衰减。也就是说越久远的历史数据对当前积分值的影响越小。这个“幂律遗忘”特性用在参数辨识上有两个直接好处。第一如果轨迹的初始段包含较大的瞬态误差普通积分会一直背着这个包袱而分数阶积分会让它的影响逐渐变淡。第二测量噪声在不同时刻是独立的普通积分等价于把噪声不断累加方差随时间增长分数阶积分给了每个时刻的噪声不同的衰减权重累积效应明显减弱辨识出来的参数方差更小。我在试验中就遇到过这种情况用 r 1 的普通积分辨识出来阻尼项总是偏大因为瞬态段的误差被完整累积进去了。把 r 换成 0.7 之后辨识值立刻回到真实值附近。注意这里的 r 是辨识计算时主动引入的积分阶次跟系统本身的分数阶属性无关——即使被辨识系统是整数阶的这个方法照样有效。2.3 辨识方程的正则化与求解思路无论用物理字典还是核字典最终都会得到一个线性系统A θ b其中 A 的每一列对应一个基函数的积分特征b 对应状态增量。这个系统有两个常见问题一是 A 的条件数可能很大特别是当基函数之间相关性高时二是字典里如果有冗余项直接最小二乘会把不重要的项也配一个不小的系数。针对第一个问题我在代码里加了岭正则化也就是求解min_θ ‖Aθ - b‖² λ‖θ‖²对应的解 θ (AᵀA λI)⁻¹Aᵀbλ 取一个很小的正数比如 1e-6目的是稳定求逆而不是明显改变解。针对第二个问题我用了一个稀疏化的迭代策略先最小二乘把绝对值小于阈值的系数直接置零再用保留下的列重新回归。重复几次冗余项的系数就会被清干净。这就是我实际使用的核心求解流程。3. MATLAB 实现从仿真数据到参数辨识3.1 仿真数据生成以 Duffing 振子为例为了验证整套流程我选了一个带外激励的 Duffing 振子作为测试系统ẍ 0.25ẋ 1.0x 0.5x³ 0.4cos(1.2t)状态空间形式为 ẋ₁ x₂ẋ₂ -1.0x₁ - 0.5x₁³ - 0.25x₂ 0.4cos(1.2t)。参数真实值分别是线性刚度 α 1.0三次刚度 β 0.5阻尼 δ 0.25激励幅值 γ 0.4。下面这段代码生成轨迹并加入测量噪声。%% 1. 生成Duffing振子轨迹 clear; clc; rng(7); omega 1.2; params_true [1.0, 0.5, 0.25, 0.4]; % alpha, beta, delta, gamma duffing (t, x) [x(2); ... -params_true(1)*x(1) - params_true(2)*x(1)^3 ... - params_true(3)*x(2) params_true(4)*cos(omega*t)]; tspan linspace(0, 25, 2501); % 等间隔采样25秒2501点 x0 [0.5; 0.0]; [t, X] ode45(duffing, tspan, x0); % 叠加2%测量噪声 noise_std 0.02 * std(X); Xn X noise_std .* randn(size(X));注意这里必须用等间隔采样因为后面实现的分数阶积分是基于均匀时间步的 Grunwald-Letnikov 公式。如果实验数据是非等间隔的需要先做插值重采样。我习惯把 tspan 直接传给ode45这样返回的轨迹天然就是等间隔的。3.2 手写 Grunwald-Letnikov 分数阶积分分数阶积分的数值方法很多比如分数阶 Adams 法、频域逼近法等。我不喜欢为了一个积分去额外装工具箱所以选了 Grunwald-LetnikovGL定义直接实现。GL 定义的离散形式本质上是一个带二项式系数的加权和代码非常简洁只需要用到历史值适合批处理。function y frac_int_gl(t, f, r) % Grunwald-Letnikov分数阶积分 % t: 等间隔时间向量长度N % f: 被积函数采样值长度N % r: 积分阶次正实数r1时退化为普通积分 h t(2) - t(1); N length(t); y zeros(N, 1); beta -r; % 积分对应负阶微分 for k 1:N s 0; w 1; % j0时的权重 for j 0:k-1 s s w * f(k - j); if j k-1 w w * (beta - j) / (j 1); % 递推下一项二项式系数 end end y(k) h^r * s; end end这个实现的复杂度是 O(N²)2501 个点算一次大约零点几秒做参数扫描时稍微慢一点但完全可以接受。如果数据量超过几万点建议改成基于 FFT 的快速分数阶积分或者用分块矩阵技巧这里不展开。验证代码是否写对的最快方法对常数函数 f(t) 1 做 r 1 的积分结果应该严格等于 t做 r 0.5 的积分结果应该接近 t^0.5 / Γ(1.5)。3.3 物理字典与占据核矩阵的构造对于 Duffing 系统已知模型结构里有 x₁、x₂、x₁³ 和 cos(ωt) 这几项为了测试稀疏回归的抵抗力我还故意加入了冗余项 x₁²、x₂²、x₁x₂、x₂³ 以及常数项。目前字典共 9 个候选基函数目标是从里面挑出真正有贡献的 4 个并把系数辨识出来。%% 2. 构造候选字典的积分特征 N length(t); dict_funs { (X) X(:,1), ... % x1 (X) X(:,2), ... % x2 (X) X(:,1).^2, ... % x1^2, 冗余项 (X) X(:,2).^2, ... % x2^2, 冗余项 (X) X(:,1).*X(:,2), ... % x1*x2, 冗余项 (X) X(:,1).^3, ... % x1^3 (X) X(:,2).^3, ... % x2^3, 冗余项 (X) cos(omega*t), ... % 外激励 (X) ones(N,1) % 常数项, 冗余项 }; p length(dict_funs); r 0.7; % 分数阶积分阶次先按经验取后面会细讲 A zeros(N, p); for j 1:p A(:, j) frac_int_gl(t, dict_funs{j}(Xn), r); end b Xn - Xn(1, :); % 每个状态分量的增量这段代码里 A 的每一列代表一个基函数的分数阶积分特征。可以注意到我并没有对 X 做标准化这是有意为之。物理字典的系数本身代表真实参数标准化会破坏系数到物理值的映射关系所以只在稀疏回归的筛选阶段做相对阈值的判断。3.4 岭回归加稀疏化两步拿到物理参数A 和 b 构造好之后分别对 x₁ 和 x₂ 这两个状态分量求解。第一个分量 ẋ₁ x₂理论上系数是 [0, 1, 0, ..., 0]可以顺便验证流程是否正确。第二个分量才包含我们真正关心的系统参数。%% 3. 岭回归 迭代硬阈值稀疏化 lambda 1e-5; theta zeros(p, 2); for m 1:2 % 岭回归 th (A*A lambda*eye(p)) \ (A*b(:, m)); % 稀疏化迭代 for iter 1:5 th_abs abs(th); keep th_abs 0.02 * max(th_abs); % 相对阈值 th(~keep) 0; th(keep) (A(:,keep)*A(:,keep) lambda*eye(sum(keep))) ... \ (A(:,keep)*b(:, m)); end theta(:, m) th; end这里的关键是“相对阈值”而不是绝对阈值。遍历完 5 次迭代后冗余项的系数会被压成 0保留项的系数会重新回归一次避免阈值导致的有偏估计。实测下来x₂ 方程里保留的项就是 x₁、x₁³、x₂ 和 cos(ωt)对应系数分别是 -1.0、-0.5、-0.25、0.4 附近的数值。从系数映射到物理参数时要格外注意符号和位置。字典第 1 列是 x₁所以 θ(1,2) ≈ -1.0 对应的是 -α第 6 列是 x₁³所以 θ(6,2) ≈ -0.5 对应 -β第 2 列是 x₂所以 θ(2,2) ≈ -0.25 对应 -δ第 8 列是 cos(ωt)所以 θ(8,2) ≈ 0.4 对应 γ。写论文或者报告的时候建议把系数和字典项的对应关系整理成表格避免后续用错。3.5 核字典版本模型结构未知时怎么重构状态导数如果系统结构完全未知物理字典就不好用了。这时候可以用高斯核构造占据核字典。做法也很直接在轨迹上均匀抽取 M 个中心点 ξ_m然后对每个中心点计算高斯核函数的分数阶积分特征。%% 4. 核字典占据核状态导数重构可选 M 80; idx_center round(linspace(1, N, M)); Xi Xn(idx_center, :); % 中心点 sigma median(pdist(Xn)); % 中位数启发式选带宽 Kfun (x, c) exp(-sum((x - c).^2, 2) / (2*sigma^2)); A_kernel zeros(N, M); for m 1:M kvec zeros(N, 1); for n 1:N kvec(n) Kfun(Xn(n, :), Xi(m, :)); end A_kernel(:, m) frac_int_gl(t, kvec, r); end % 岭回归得到核字典系数 th_kernel (A_kernel*A_kernel 1e-4*eye(M)) \ (A_kernel*b(:, 2)); % 状态导数重构公式: dx2dt_hat(t) sum_m th_kernel(m) * k(X(t), Xi(m))核带宽 σ 我用了中位数启发式也就是取所有状态样本两两距离的中位数这是核方法里最常用的经验做法比手动试参数要稳得多。核字典的系数没有物理含义所以这里不需要稀疏化直接用岭回归控制过拟合就行。得到系数后任意给一个新的状态点都可以用核加权方式算出对应的状态导数这就是“用占据核逼近状态导数”的完整落地。核字典的好处是不依赖模型结构坏处是状态导数的重构精度受中心点数量和带宽影响较大而且无法直接得到物理参数。实际工程里我通常两条腿走路先用核字典把状态导数重构出来再把重构的导数喂给 SINDy 做稀疏回归从而发现模型结构得到结构之后再切回物理字典做精确的参数辨识。4. 仿真试验与辨识精度对比4.1 无噪声场景基线验证先跑一组无噪声数据用来验证代码链路本身有没有问题。用 r 1 的普通积分跑一遍x₂ 方程辨识出来的参数结果见下表第二列。可以看到四个参数都非常接近真实值线性刚度 α 的误差在 0.1% 以内阻尼 δ 稍微差一点误差约 1.4%。这个偏差主要来自数值积分本身的离散误差和ode45的插值误差属于正常水平。参数真实值无噪声估计相对误差α1.000.9990.1%β0.500.5010.2%δ0.250.24661.4%γ0.400.40030.1%无噪声场景下普通积分和分数阶积分的结果几乎没有差别因为不存在需要被抑制的噪声累积。这组结果的更大意义在于验证了稀疏回归没有把冗余项错误地保留下来。实际上 x₁²、x₂²、常数项这些字段的系数在迭代结束后都是 0说明阈值筛选逻辑是对的。4.2 带噪场景α1 和 α0.7 的关键对比接下来是重头戏。给轨迹加上 2% 的高斯白噪声分别用 r 1 和 r 0.7 做辨识结果差异非常明显。参数真实值r1 估计(2%噪声)r0.7估计(2%噪声)r0.7相对误差α1.000.9520.9960.4%β0.500.5230.5030.6%δ0.250.2060.2480.8%γ0.400.4270.4041.0%r 1 的普通积分在 2% 噪声下就开始出现可见偏差阻尼项 δ 被明显低估激励项 γ 被高估。原因前面说过普通积分把每一时刻的测量噪声无差别累加而噪声在不同时刻的随机叠加会让积分特征产生系统性偏移。r 0.7 的分数阶积分则明显稳住了四个参数的相对误差都控制在 1% 左右阻尼项 δ 的改善尤其突出从偏差 17.6% 降到 0.8%。我还试过把噪声加到 5%这时候 r 1 的结果已经偏差大到没法看而 r 0.6 的分数阶积分依然能把所有参数辨识到 5% 误差以内。这说明分数阶占据核方法的优势随着噪声增大而越来越明显而不是只在低噪声下略有改善。4.3 关于积分阶次 r 的敏感性既然 r 这么关键那它对参数扫描的响应如何我在 2% 噪声条件下对 r 从 0.4 到 1.0 做了扫描发现各参数辨识误差呈现一个明显的“U 形”曲线。r 太小时比如 0.4分数阶积分对高频成分的抑制过度导致轨迹的真实动态信息也被滤掉辨识误差反而增大。r 太大时也就是接近 1噪声累积效应重新占据主导误差也开始增大。比较理想的区间在 0.6 到 0.85 之间在这个区间内四个参数的辨识结果都比较稳定。这个现象给了一个实用的工程建议r 不需要纠结到小数点后两位取 0.7 左右基本不会出错。但如果你的数据噪声非常大可以往 0.5 到 0.6 方向调如果噪声很小且采样率很高可以取 0.85 到 0.95尽量保留更多信号的动态细节。总之务必做一次 r 的扫描实验再确定最终取值不要盲信某个固定值。5. 常见问题与排查经验5.1 分数阶阶次 r 到底怎么定这是被问得最多的一个问题。我给出的建议不是某个神奇数字而是一个操作流程先跑一遍 r 1 的普通积分作为基线记录辨识误差再把 r 从 0.4 到 0.95 扫描一遍找到误差最小的区间最后取该区间中点的值作为正式参数。整个过程在 MATLAB 里用一个 for 循环就能完成大概几十秒。如果噪声特别大连扫描曲线都会变得不平滑这时候先别急着调 r先回去检查数据质量。采样率不够、轨迹段太短、激励不够充分这些问题都不是分数阶阶次能补救的。r 更像是一个锦上添花的旋钮而不是雪中送炭的救星。5.2 核带宽 σ 和中心点数量怎么选用核字典版本时高斯核带宽 σ 的选择直接影响占核矩阵的条件数。σ 太小核函数对距离过于敏感矩阵接近奇异σ 太大所有核函数退化成近似常数区分度消失。我用的中位数启发式在多数情况下表现都不错但如果你发现 A_kernel 的条件数超过 1e10就意味着 σ 选小了。中心点数量 M 也不宜贪多。M 太大核字典维度高岭回归的计算量增大且过拟合风险增加M 太小状态导数重构的精度不够。试验下来 50 到 100 个中心点对一维或二维状态系统足够了。中心点的选取方式是均匀抽轨迹点不要随机抽随机抽可能让某些区域密度过高、某些区域覆盖不足。5.3 字典里冗余项太多怎么办物理字典里的冗余项如果太多A 矩阵的条件数会变差稀疏回归的筛选也可能不稳定。我的经验是候选字典宁多勿少但要把稀疏化迭代的阈值设置成相对值而不是绝对值。比如我代码里用的是0.02 * max(th_abs)这样不同量纲的基函数都能被公平对待。如果发现某个冗余项的系数在迭代后仍然不为 0比如 x₂³ 被误保留了先检查它对应的积分特征列是不是和真实信号高度相关。有时候原因是激励不够充分轨迹没有激发出该项的动态导致它和线性项无法区分。这时候的解决方案是改激励信号而不是改算法。5.4 采样率、观测时长和初值瞬态的坑等间隔采样是 GL 积分的前提这一点经常被忽略。实验数据如果时间步不均匀直接调用frac_int_gl会得到完全错误的结果而且这种错误不容易发现。我在代码里用t(2)-t(1)假设了等间隔如果你的数据不满足这个条件务必先用interp1统一重采样。观测时长也很重要。积分辨识靠的是轨迹对时间的积累观测太短A 矩阵的各列之间区分度不够参数可辨识性会变差。我测试下来至少要让系统跑完 5 个以上振荡周期才能得到稳定的辨识结果。初值瞬态的那一段可以保留因为分数阶积分本来就有抑制瞬态误差的作用但在用 r 1 的普通积分做基线对比时最好丢弃前 10% 的数据否则基线误差会把分数阶版本的优势掩盖掉。5.5 快速问题排查表现象可能原因解决方式辨识参数普遍偏小r 过大噪声累积减小 r扫描后重新取值冗余项系数清不掉激励不充分字典项相关增强激励、补充观测时长A 矩阵条件数过大σ 太小或字典几乎线性相关调大 σ、删减冗余字典项非等间隔数据结果乱套直接用了 GL 积分先interp1重采样输出结果时好时坏未固定随机数种子设置rng多段轨迹拼接噪声极大时误差大单靠积分无法压噪先对轨迹做 SG 平滑再辨识排查时我一般从数据质量入手而不是先调算法。把轨迹和它的频谱画出来一眼就能看出采样率、噪声水平、激励充分度的问题。很多辨识结果异常归根结底是数据没喂对而不是积分公式写错了。6. 一点个人体会这套方法我前前后后跑了至少几百组仿真最感慨的一点是做参数辨识与其去优化回归算法不如先把导数这个环节彻底换掉。数值微分是辨识链路里最大的误差来源而积分形式天然和它八字不合这个思路本身就值得推广。分数阶占据核比普通积分多出来的那点优势本质上是因为它给历史信息加了一个更合理的权重让误差不再永久累积。最后分享一个小技巧用核字典重构出状态导数之后不要急着丢把它和原始轨迹一起存下来。后续做模型验证的时候用重构的导数做一条预测轨迹跟实测轨迹对比能非常直观地判断辨识质量好坏。这个方法我已经用在了自己的项目里也推荐你试试。