从一堆只有时间戳的传感器读数里把背后那个决定系统演化规律的微分方程直接挖出来——这事放在五年前我觉得是科幻现在做下来发现方法对了就是几行 Matlab 的事。这篇想聊聊我最近跑通的数据驱动项目基于时域数据的分数混沌系统的稀疏识别全程用的Matlab 代码把从生成数据到识别方程再到验证结果的一套流程完整走了一遍。先说清楚这东西是干啥的。所谓稀疏识别本质是从海量候选函数里挑出“刚好够用”的那几项还原出系统的真实动力学方程。传统的系统辨识依赖大量先验假设神经网络拟合又是黑箱给不出方程结构而稀疏识别走的是白箱路线直接给你一个数学表达式。它特别适合这些场景手头有实验或仿真产生的时域数据知道系统大概有哪些变量但不确定具体方程长什么样或者想验证某个分数阶混沌模型是否真的能解释观测数据。如果你是做非线性动力学、故障诊断、气象数据分析、生物系统建模这类工作的这篇文章的操作可以直接参考。1. 项目到底在做什么一句话拆解核心逻辑1.1 数据驱动方程识别的定位数据驱动这个词这两年很热但很多人把它等同于深度学习。实际上在动力学系统这圈子里数据驱动还有一个更朴素、更工程化的含义直接从观测数据中重构控制方程不依赖完整的机理建模。SINDySparse Identification of Nonlinear Dynamics是这类方法的代表第一次看到这个思路的时候我是真觉得妙把系统状态及其函数组合排列成一个大型候选函数库再把系统导数表示成这个库的线性组合用稀疏回归去求解系数——绝大多数学者几行代码就能把这个系数向量算出来。我之前用 SINDy 做过几个整数阶系统效果都不错。但这回面对的是分数阶混沌系统难度直接上了一个台阶。1.2 为什么选分数阶混沌系统当靶子先解释下“分数阶”三个字。这里的分数不是数学考试里的分数而是微积分的一种推广——导数阶次可以是 0.99、1.2 这样的实数描述具有记忆效应和遗传特性的物理过程特别合适。很多实际系统比如粘弹性材料、电极化过程、以及某些生物组织的动力学整数阶模型往往给不出满意的拟合优度换成分数阶模型后却能精确匹配。当这种分数阶系统进入混沌状态也就是对初始条件极端敏感模型识别就变成一道很有挑战性的考题。选择分数混沌系统作为识别对象还有一个实际原因这种系统的时域信号形态极其复杂包含很宽的频率成分和强烈的非线性耦合常规的参数辨识方法基本失效。如果能在这个高难度场景下把方程还原出来那这套方法迁移到普通系统上就是降维打击。1.3 这个方法能带到哪些领域我总结了这套识别流程可以在哪些方向复用加工产线轴承故障诊断采集振动时域信号识别出系统方程后方程特征的变化能更早地暴露故障迹象。脑电/心电生物信号建模生理信号往往呈现混沌特性稀疏识别可以给出可解释的数学模型辅助诊断。化工过程异常状态监测从流场或温度场时序数据中还原反应动力学模型。新能源系统的充放电建模电池和电容器都有明显的分数阶特性识别出精确模型对 SOC 估计很有帮助。这里多说一句很多人会问为什么不用神经网络神经网络拟合精度高但给不出显式方程而且数据量不够时容易过拟合稀疏识别在小样本、有物理约束的场景下是性价比高得多的方案。2. 稀疏识别的原理与工程适配2.1 SINDy 框架回顾先花点篇幅把 SINDy 的基础捋一遍。设想一个连续时间动力学系统dx/dt f(x(t))给定时间序列观测矩阵 X每一行对应一个时刻的状态先估计出导数矩阵 dX/dt每一行对应同一时刻的导数值。接下来构造候选函数库 Theta(X)它的每一列是状态变量的某一个候选函数比如常数项、各个状态分量、二次项、三次项、三角函数等Theta [1, x, y, z, x*y, x*z, y*z, x^2, y^2, z^2, ...]如果系统的真实方程恰好可以稀疏地表示为这些候选函数的线性组合那么就存在一个系数矩阵 Xi使得dX/dt Theta * Xi关键在于 Xi 是稀疏的——每行只有少数非零系数。求解这个欠定方程组时直接最小二乘会得到一堆非零的小数完全没有解释性。SINDy 的做法是用“最小二乘硬阈值”迭代先最小二乘出一个初始解把绝对值小于阈值的系数直接清零再用保留下来的几十项重新做最小二乘反复迭代到稳定。最终得到的 Xi 包含大且显著的系数其余都是零方程结构一目了然。这里可以拿“整理房间”打比方所有东西都放在地上Theta 矩阵你要找出几件真正重要的东西其余全部扔掉。每次只凭“东西的大小”系数绝对值做判断一次不够多扔几轮房间就干净了。2.2 从整数阶到分数阶核心改动在哪整数阶系统的导数好算差分即可。但分数阶系统里的状态变量满足的是分数阶微分方程D^q x(t) f(x(t)), q 为非整数这里 D^q 是分数阶导数算子。常用的定义有 Caputo、Riemann-Liouville 和 Grünwald-LetnikovGL三种。数值上最容易实现的是 GL 定义D^q x(t) lim(h→0) h^(-q) Σ_{k0}^{∞} w_k x(t - kh)其中 w_k 是记忆系数按二项式展开系数递推得到w_0 1, w_k w_{k-1} * (q - (k-1)) / k注意这个求和牵扯到全部历史时刻——这就是分数阶系统“记忆性”的来源也正是它区别于整数阶系统的本质。分数阶导数具有记忆效应这就意味着预测未来时需要知道过去较长一段时间的全部信息。识别流程对应的改动是在构造候选库之前把整数阶导数 dX/dt 换成 GL 公式算出来的分数阶导数 D^q X。候选库的构造方法不变稀疏回归流程也不变。看起来只是换了个导数计算模块但这正是很多复现的人失败的点——直接用整阶差分去逼近分数阶导数识别出来的方程结构再漂亮也是错的。我仿真验证时还特别注意了一个原则生成数据用的是预估校正法FDE12工具箱识别导数用的是 GL 公式这两个数值方案完全独立。如果识别结果能还原出真实的分数阶混沌方程就说明算法本身真的 work而非“一个数值公式出自己的影子”。2.3 候选库怎么搭才不踩坑候选库的设计直接决定识别成败。我见过很多人一股脑堆到五阶多项式结果矩阵条件数爆炸回归出来的系数全在抖。候选库不是越厚越好而是要匹配系统的物理复杂度。一个典型的候选库设计思路Theta [1, x, y, z, xy, xz, yz, x^2, y^2, z^2, x^3, y^3, z^3]对 Lorenz 这种经典混沌系统通常最高三次方就够用。如果系统里可能存在正弦、余弦项再把 trig 项加进去但每一项都会显著增加矩阵的维度和数值敏感度。这里有个工程细节各候选列的数值尺度差异非常悬殊。x、y 本身是 O(10)xy 是 O(100)x^3 是 O(1000)。直接拿原始数据构造候选库再回归数值上很容易出问题。我的做法是先对候选库做标准化每列减去均值除以标准差求出稀疏系数后再映射回原始坐标这样既保持了回归的数值稳定性又能直接读出方程的物理系数。2.4 稀疏回归的收敛性问题STLSQ 这个迭代过程虽然简单但也不是随便跑跑都能收敛。核心在于阈值 lambda 的选择。lambda 太小保留一堆机小数方程看起来就是“假大空”lambda 太大把真实项也砍掉了系统直接变成错误的退化模型。实际调试中我不能做到一步到位。正确姿势是从比较大的 lambda 开始比如 0.1 的量级看识别出的候选函数集确认某些项被错误裁剪后再逐渐调小。如果每次从零开始搜整个参数空间既费时间又不容易找到全局最优。项目里我一般是先扫描几个量级的 lambda确定合理区间再在这个区间里细化搜索。3. Matlab 代码实现全流程3.1 准备环境和数据生成运行环境很简单R2022b 或以上版本即可不需要额外工具箱FDE12 需要在网上搜索下载免费。第一步生成作为“观测数据”的时域序列。我用的靶系统是经典的分数阶 Lorenz 系统也是验证稀疏识别时最常用的基准模型D^q1 x σ(y - x) D^q2 y ρx - y - xz D^q3 z xy - βz参数取 σ10、ρ28、β8/3阶次 q1q2q30.995这是常见的混沌参数组合。初始条件取 [0.1, 0.1, 0.1]时间步长 h0.01时长 T40s丢掉前面 5s 的瞬态每隔 10 步采一个点。生成数据的核心代码% 读取 FDE12 函数后生成时域数据 T 40; h 0.01; tspan 0:h:T; y0 [0.1; 0.1; 0.1]; q [0.995; 0.995; 0.995]; % FDE12 为分数阶预估校正求解器自行搜索下载 [t, y] FDE12(q, fracLorenzFunc, y0, tspan); % 丢弃前 5s 瞬态数据避免初始过渡段干扰 startIdx find(t 5, 1); t t(startIdx:end); y y(startIdx:end, :); % 降采样降低候选库规模 sampleEvery 10; t t(1:sampleEvery:end); y y(1:sampleEvery:end, :); save(fracLorenzData.mat, t, y);其中fracLorenzFunc是系统的分数阶右端函数function f fracLorenzFunc(t, x) sigma 10; rho 28; beta 8/3; f [sigma*(x(2) - x(1)); rho*x(1) - x(2) - x(1)*x(3); x(1)*x(2) - beta*x(3)]; end生成后我先画了时域曲线和三维相图确认系统确实在混沌状态——Lorenz 的蝴蝶形状还在只是形态比整数阶的略“胖”一些说明分数阶记忆项引入了额外复杂性。3.2 分数阶导数计算识别过程中最关键的一步用 GL 公式计算 D^q x、D^q y、D^q z。我的 MATLAB 实现直接用向量化卷积加速核心代码很短function dq glFD(x, q, dt) % 基于 Grunwald-Letnikov 定义的分数阶导数 % x : 输入时域信号列向量 % q : 分数阶阶次比如 0.995 % dt: 采样时间间隔 % 返回 D^q x(t)长度与 x 相同 N length(x); % 递推生成二项式记忆系数 w_k w zeros(1, N); w(1) 1; for k 2:N w(k) w(k-1) * (q - (k-2)) / (k-1); end % 用卷积实现历史加权求和然后除以 dt^q dq conv(x, w, same) / (dt^q); dq dq(:); end这里面的细节值得多说一句GL 定义要求对所有历史时刻加权数据长度 N 越长记忆系数 w 的累计效应越完整所得导数越接近真实值。但conv的 same 写法在序列两端会有截断效应我实测下来开头和结尾各 2%~5% 的数据点误差会偏大。处理办法是丢弃两端各 50 个点再进行后续回归省心且稳。这个函数的循环乍一看是逐项递推其实 O(N) 复杂度Matlab 跑 3000 个点完全没压力。我就是这样算出了三个状态变量的分数阶导数load fracLorenzData.mat; dt t(2) - t(1); q 0.995; dxq glFD(y(:,1), q, dt); dyq glFD(y(:,2), q, dt); dzq glFD(y(:,3), q, dt); % 丢弃边界效应区域 cut 50; xdata y(cut:end-cut, :); dxdata [dxq, dyq, dzq]; dxdata dxdata(cut:end-cut, :);3.3 候选库与 STLSQ 稀疏回归有了状态矩阵 X 和导数矩阵 dXdt接下来进入核心部分构造候选库、调用 STLSQ。候选库我用的是 0~3 阶多项式及交叉项共 20 列。构造函数和稀疏回归函数都给你放在这里function Theta poolData(x, p) % 构造多项式候选库 % x : N x m 状态矩阵N 个采样点m 个状态变量 % p : 最高多项式阶次 % 返回 Theta : N x C 候选函数矩阵 [mRow, mCol] size(x); nTerms nchoosek(mCol p, p) * p; % 粗略上界下面用动态增加 Theta ones(mRow, 1); % 常数项 % 逐阶构造 for ord 1:p % 获取所有 mCol 个变量的 ord 次组合含交叉项 idxList myNchooseK(mCol, ord); % 见下方说明 for i 1:size(idxList, 1) term ones(mRow, 1); for j 1:ord term term .* x(:, idxList(i, j)); end Theta [Theta, term]; end end end % 简洁的替代方案直接用多项式符号生成组合 % 工程上更稳妥的做法是使用字符串数组维护每一项的名称myNchooseK在工程中可以借助nmultichoosek等函数实现这里不再展开。为了不让代码读起来像天书调试时更好用的是用字符串记录每个候选函数的名字识别出来后可以直接拼出表达式打印出来function names poolNames(mCol, p) names {1}; for ord 1:p idxList myNchooseK(mCol, ord); for i 1:size(idxList, 1) s ; for j 1:ord s [s, sprintf(x%d*, idxList(i,j))]; end s(end) []; % 去掉末尾乘号 names{end1} s; %#okAGROW end end endSTLSQ 核心迭代代码function Xi sparsifyDynamics(Theta, dXdt, lambda, iters) % 稀疏回归最小二乘 硬阈值迭代 % Theta : N x C 候选库 % dXdt : N x m 导数矩阵 % lambda: 稀疏阈值 % iters : 最大迭代轮数 % 初始最小二乘解 Xi Theta \ dXdt; for k 1:iters smallidx abs(Xi) lambda; % 小于阈值的位置置零 Xi(smallidx) 0; for ind 1:size(Xi, 2) bigidx ~smallidx(:, ind); % 只用保留下来的候选列重新求解 Xi(bigidx, ind) Theta(:, bigidx) \ dXdt(:, ind); end end end这个函数看起来简单但它背后的哲学是每次只保留下绝对系数足够大的准则项再用它们重新拟合给那些被误伤的真实小系数项复活的机会。我实测在 Lorenz 这种系统上一般 5~10 轮迭代后系数就稳定了。3.4 模型输出与拟合效果对比跑完回归后把稀疏系数组织成方程与真实方程做对比。结果长这样% 假设状态列顺序为 [x, y, z]候选库前几列对应 [1, x, y, z, ...] Theta poolData(xdata, 3); names poolNames(3, 3); lambda 0.05; Xi sparsifyDynamics(Theta, dxdata, lambda, 10); % 打印识别出的方程 for i 1:3 fprintf(D^q x%d , i); first true; for j 1:length(names) if abs(Xi(j, i)) 1e-6 if ~first fprintf( ); end fprintf(%.3f * %s, Xi(j, i), names{j}); first false; end end fprintf(\n); end在干净数据下识别结果和真实方程几乎一模一样σ、ρ、β 的估计值跟真值贴合到小数点后两三位非物理候选列的系数都在 1e-6 量级以下直接滤掉了。这套结果说明“观测数据-分数阶导数-候选库-稀疏系数-方程结构”的链条是通的。4. 参数怎么调、坑怎么填实测经验汇总4.1 关键参数速查表参数推荐取值影响与备注候选库最高阶 p3~4超过 4 阶矩阵条件数急剧上升谨慎使用稀疏阈值 lambda0.01~0.1从大到小扫描确认保留项后再细化迭代次数 iters5~2010 轮内基本收敛太多轮无益采样时间 dt0.005~0.02太小导数误差大太大丢失高频细节记忆项长度 N尽量用全序列截断会引入导数误差必要时两端丢弃丢弃边界点序列两端各 50 点消除 GL 卷积边界效应4.2 噪声现实从干净数据到实测数据前面跑的是干净仿真数据但真实传感器信号不可能这么完美。我给生成的时域数据叠加上信噪比 20dB 的白噪声再把识别流程跑一遍方程结构基本还能保留但系数误差上升到 10% 左右。如果信噪比降到 10dB稀疏回归就很容易把噪声项当成真实候选函数方程开始“注水”。应对噪声的经验有三条导数做平滑分数阶导数本身对噪声的放大比整数阶更严重因为记忆系数对历史噪声是长期累积的。我建议先用 Savitzky-Golay 滤波器对原始时序做平滑再去算 GL 导数实测噪声敏感度能降低一个量级。增大记忆长度分数阶系统的分数阶导数计算依赖足够长的历史数据数据太短时导数本身就是错的噪声再一叠加稀疏回归直接失效。至少保证上千个采样点。lambda 适当上调有噪声时不能照搬干净场景的阈值适当调大可以滤掉噪声项但注意别砍掉真项。4.3 常见问题速查Q1识别出来的方程全是零怎么回事先查导数是否计算正确。我对这种问题最常做的就是直接打印dyq(1:10)看数值量级常常发现分数阶导数比状态本身小了 2~3 个数量级导致候选库系数约等于零。这时需要先确认 q 与系统定义一致再检查 dt 是否跟积分步长一致——这两个都对了导数量级才会正常。Q2识别结果包含大量微小非零系数稀疏性不够。这是 lambda 取太小或者候选库列之间存在强共线性。前者调大 lambda后者做标准化的同时剔除一些相似项比如 x^2 与 x*y 在部分区间高度相关。我一般用列之间相关系数矩阵看一眼超过 0.95 的候选对只保留大概率物理相关的那一个。Q3使用 FDE12 生成的数据和 GL 导数之间实测有系统偏差。这是正常的。FDE12 是预估校正格式GL 是直接离散化定义两者在数值上是近似方法收敛速度不同。不要把偏差理解为代码 bug采用 3.2 里的两端截断策略可以显著减小这部分影响。Q4Matlab 提示矩阵接近奇异或条件数过大。候选库阶数太高或者候选库包含常数项与某个状态列几乎线性相关比如某个状态在数据段内均值接近常数列。建议降低最高阶次或者先对状态变量做中心化后再构造库。Q5跑 FDE12 时速度很慢。分数阶数值积分本身就带有记忆性复杂度高于整数阶求解器。时间步长从 0.01 放到 0.02通常可接受同时把数据量控制在 3000 点内单次生成几秒就能跑完。这套方法调试时我最大的体会是稀疏识别不是无脑黑盒工具它对中间步骤的数值处理极其敏感。很多人复现 SINDy 失败问题往往不在稀疏回归本身而在导数的计算——分数阶系统更是如此。只要你把 GL 导数、候选库尺度、稀疏阈值这三关把好了从时域数据到方程结构的全过程是很稳定的。最后分享一个小技巧识别出方程后别急着部署拿一组全新的初始条件重新用识别方程做数值仿真再与真实观测时域波形叠图画对比。如果这条曲线和观测曲线在足够长时间内仍然纠缠在一起说明识别方程不是“拟合得好”而是真的掌握了系统的动力学。这一步的验证价值比任何评价指标都更能说明问题。