简介极大似然法是系统辨识中常用的参数估计手段能从输入输出数据中反推模型参数递归极大似然RML则通过数据逐次到达实现参数在线更新非常适合实时系统。这套资料将两者从数学原理到工程实现做了完整串联面向控制类专业学生、科研人员以及对参数估计感兴趣的工程师。压缩包共19个文件6个MATLAB脚本实现了RML递推、极大似然计算和m序列PRBS激励生成9张结果图直观展示白噪声与有色噪声下的参数估计误差、输出误差另有2个.mat数据文件、1份PPT演示文稿和1个Simulink工程文件整体仅865KB轻量易用。目前已有187人参与学习。借助代码、图表和演示文稿读者既能理解RML每次迭代中的似然函数递推细节也能通过不同噪声条件的对比掌握算法鲁棒性还能在Simulink工程中自行验证。这份材料对想动手复现系统辨识算法、理解在线更新机制的读者来说是一份短小精悍的参考资料。1. 系统辨识里的极大似然法这份 RML 资源能替你省掉多少弯路做系统辨识的人迟早会遇到一类怪问题输入输出数据明明很干净阶次、延迟、采样周期全都按经验调好了可最小二乘辨识出来的参数就是和真实值对不上偏差还带着明显规律。原因通常是噪声不是白噪声而是和输入相关或被系统动态“染色”的有色噪声。极大似然法把噪声模型也放进待估参数里从统计上消除了这种偏差而递推极大似然RML是其中工程上最实用的递推实现。这份资源把原理推导、递推算法和配套验证代码放在一起适合做系统辨识课程设计、课题研究或者手里有现场数据需要建模的工程师。2. 从最小二乘的偏差到 RML 的修正噪声模型为什么必须估计2.1 有色噪声怎么把最小二乘拉偏的先说一个我反复见过的现象有人拿一组带测量噪声的数据做最小二乘辨识损失函数显示拟合误差很小可模型参数却和物理常识对不上。问题出在最小二乘的假设上——它默认残差是零均值、方差恒定、互不相关的白噪声序列。真实场景里噪声经过被控对象或传感器后往往带上相关性上一时刻的噪声没消散又混进了下一时刻的输出里。典型模型是带色噪声的方程y(t) a1 y(t-1) a2 y(t-2) b1 u(t-1) b2 u(t-2) e(t) c1 e(t-1)这里的 e(t) 是白噪声c1 e(t-1) 就是有色成分。最小二乘估计时回归矩阵里的 y(t-1) 和当前残差里的 e(t-1) 是相关的。这个相关性直接违反了最小二乘的 Gauss-Markov 假设结果就是参数估计不再无偏。用数学语言说最小二乘的估计公式是 θ̂ (ΦᵀΦ)⁻¹ΦᵀY当噪声项和回归变量相关时E[Φᵀe] ≠ 0于是估计量的期望值里多了一个偏差项。这个偏差不是多采点数据就能消除的数据越多只是让参数“稳定地”偏向错误值。我见过有人在现场纠结了三天反复滤波、去掉坏点最后发现是噪声相关性问题没处理。2.2 RML 的核心思路把噪声多项式 C(q) 也拉进优化目标极大似然法的思路和最小二乘不一样。它不再假设残差必须白而是直接给残差建一个模型 C(q)e(t)把 C 多项式里的系数 c1...cnc 一并当成待估参数。这样一来优化目标就不是简单的“平方损失最小”而是“给定输入和模型出现当前观测序列的概率最大”。把似然函数对 θ 求极值得到的参数估计在理论上具有渐近无偏性和渐近有效性。这是它的统计优越性所在。但直接做极大似然要解一个非线性优化问题每步都要算高维的梯度或 Hessian数据一长就扛不住。RML 就是在递推框架里做近似极大似然每一步用当前的参数估计去更新噪声模型再用滤波后的回归量计算梯度方向递推地朝着似然极大点逼近。它本质上是伪线性回归每一步还是最小二乘的形式但回归量经过了 1/C(q) 的滤波把相关性问题逐渐剔除掉。正因为每一步都在“估计参数—更新噪声模型—再滤波数据”的循环里走RML 才能同时拿到 A、B、C 三组参数。这个“同时”很关键如果先把 C 定死再去估计 A、B等于把噪声模型当成已知条件现场根本没有这个条件。2.3 RML、递推最小二乘和增广最小二乘的分水岭在哪有人会问增广最小二乘ELS也估计噪声模型为什么还要用 RML区别在梯度方向。ELS 回归向量里直接放 y(t-i) 和残差估计 ê(t-i)不做滤波。当噪声模型系数较大时残差估计和回归量之间仍然互相牵扯收敛起来会偏甚至出现参数来回振荡的情况。RML 对回归量做了 1/C(q) 的预滤波得到的梯度方向更接近真实似然梯度。你可以这么理解ELS 是在“一边往前跑一边修正方向盘”RML 是“先看清水流方向再打舵”。看起来只是多了一步滤波但在 C 系数较大、信噪比不高的场景下两者收敛结果差别很大。实际案例里对一个 c10.8 的强色噪声系统ELS 辨识出的 b1 可能偏差 20% 以上RML 一般能收敛到真实值附近。这也是系统辨识教科书里把 RML 归入精度更高一类算法的原因。3. 手写递推极大似然RML主循环代码与三个关键参数3.1 先定 ARMAX 模型结构和符号约定写代码前先把模型结构定死否则后面全是符号混乱。本文统一采用 ARMAX 形式A(q)y(t) B(q)u(t) C(q)e(t)其中A(q) 1 a1 q⁻¹ ... ana q⁻ⁿᵃB(q) b1 q⁻¹ ... bnb q⁻ⁿᵇC(q) 1 c1 q⁻¹ ... cnc q⁻ⁿᶜ参数向量 θ [a1...ana, b1...bnb, c1...cnc]共 nanbnc 个。预测误差写为ε(t) y(t) Σ ai·y(t-i) - Σ bi·u(t-i)注意这里 A(q) 的常数项是 1所以公式里 y(t-i) 前面是正号。很多人的代码在这里翻车把符号写成负号结果整个递推过程振荡发散。后面应对排查章节会专门再说。RML 的梯度向量来自对预测器求导。对 aᵢ 求导得到的是 -yf(t-i)对 bᵢ 求导得到 uf(t-i)对 cᵢ 求导得到 -ef(t-i)。其中 yf、uf、ef 是 y、u、ε 经过 1/C(q) 滤波后的序列这个滤波用的是当前时刻的 C 参数估计值。3.2 递推主循环滤波回归器与增益更新下面这段 Python 代码实现单轮 RML 更新可以直接嵌进你自己的辨识脚本里。它对应的是 ARMAX 结构参数 na、nb、nc 由模型阶次决定。import numpy as np def rml_step(y, u, theta, P, lam, hist): RML 单步递推 y: 当前输出观测 u: 当前输入 theta: 参数向量 [a1..ana, b1..bnb, c1..cnc] P: 协方差矩阵 lam: 遗忘因子0lam1 hist: 字典保存历史 yf/uf/ef 和原始 y/u na hist[na]; nb hist[nb]; nc hist[nc] a theta[:na] b theta[na:nanb] c theta[nanb:] # 1. 预测误差A(q)y - B(q)u eps y sum(a[i] * hist[y][-(i1)] for i in range(na)) \ - sum(b[i] * hist[u][-(i1)] for i in range(nb)) # 2. 用当前 C(q) 做 1/C 滤波 yf y - sum(c[i] * hist[yf][-(i1)] for i in range(nc)) uf u - sum(c[i] * hist[uf][-(i1)] for i in range(nc)) ef eps - sum(c[i] * hist[ef][-(i1)] for i in range(nc)) # 3. 构造滤波后的梯度向量 psi psi np.zeros(na nb nc) for i in range(na): psi[i] -hist[yf][-(i1)] for i in range(nb): psi[na i] hist[uf][-(i1)] for i in range(nc): psi[na nb i] -hist[ef][-(i1)] # 4. RML 增益矩阵与协方差更新 Ppsi P psi denom lam psi Ppsi K Ppsi / denom theta_new theta K * eps P_new (P - np.outer(Ppsi, K)) / lam # 更新历史缓存 hist[y].append(y); hist[u].append(u) hist[yf].append(yf); hist[uf].append(uf); hist[ef].append(ef) for key in [y, u, yf, uf, ef]: if len(hist[key]) max(na, nb, nc) 1: hist[key].pop(0) return theta_new, P_new, eps代码逻辑分四步。第一步算预测误差第二步做 1/C 滤波这是 RML 区别于普通递推最小二乘的灵魂所在。第三步把滤波后的回归量组装成梯度向量 psi。第四步用 RML 增益公式更新参数和协方差矩阵。参数说明lam 是遗忘因子控制对历史数据的“记忆长度”P 是协方差矩阵决定参数修正步长hist 里保存的是最近几拍的历史数据用于构造回归向量。这里 yf、uf、ef 的滤波系数用的都是当前时刻的 c 参数也就是说噪声模型和系统模型在同步进化。实际调用时外层只需要循环访问数据点不断调用 rml_step并把 theta 和 P 传回即可。注意历史缓存要保留至少 max(na,nb,nc) 拍数据否则索引会越界。3.3 遗忘因子、P0、迭代策略怎么选遗忘因子 lam 是 RML 里最敏感的参数。它接近 1 时算法几乎记住全部历史稳态精度高但跟踪慢小于 0.95 时对新数据反应快但协方差矩阵容易膨胀导致参数抖动。我的习惯是先跑几组数据从 0.98 起步往下调到 0.95观察参数轨迹是否平滑。如果参数在真实值附近来回抽动就把 lam 往 1 调如果跟踪滞后明显就往 0.95 调。接收现场数据有缓慢漂移时0.96 到 0.97 是个常用区间。P0 是协方差初值它决定前面几步参数修正的幅度。P0 取得大前几步参数会大幅跳动收敛快但容易过冲取得小收敛慢而且如果数据激励不足参数会在原地“磨蹭”很久。一般取 P0 1000·I 起步数据方差大就降一点数据方差小就升一点。迭代策略上固定用一遍数据从头到尾只扫一次往往不够。常见做法是至少扫两遍第一遍用 RLS 或者大 P0 快速得到一个粗糙初值第二遍带着这个初值把 P0 降到较小值再扫一遍让参数精细收敛。这种“冷启动 热收敛”的做法比单纯加大迭代次数稳定得多。4. 用 RML 辨识一个带噪声的二阶系统仿真数据与参数收敛4.1 数据生成二阶对象 强色噪声理论说得再多不如跑一个能对照真值的例子。我们构造一个二阶对象参数真实值为 a1-1.5a20.7b11.0b20.5噪声系数 c10.8这意味着噪声相关性很强正好是 RML 发挥优势的场景。import numpy as np import matplotlib.pyplot as plt np.random.seed(42) N 2000 u np.random.choice([-1, 1], sizeN) # PRBS 输入激励充分 a_true np.array([-1.5, 0.7]) b_true np.array([1.0, 0.5]) c_true np.array([0.8]) y np.zeros(N) e np.zeros(N) for t in range(2, N): e[t] np.random.randn() * 0.3 # A(q)y B(q)u C(q)e y[t] (b_true[0]*u[t-1] b_true[1]*u[t-2] - a_true[0]*y[t-1] - a_true[1]*y[t-2] e[t] c_true[0]*e[t-1]) plt.figure(figsize(10, 3)) plt.plot(y[:300]) plt.title(有色噪声下的系统输出) plt.show()这段代码里真实参数 a、b、c 都注入了。噪声 e 是白噪声但经过 C(q) 之后进入输出的部分变成了色噪声。u 用伪随机二值序列这样系统才能被充分激励——如果输入是恒定值任何递推估计算法都救不回来因为输出变化里根本没有足够信息。注意 e[t] 的方差调成了 0.3不算高。如果噪声方差再大到了信噪比接近 1 的程度RML 虽然仍稳定但收敛速度会明显变慢需要更长的数据。4.2 跑 RML参数收敛过程与终值对比用上一章实现的 rml_step 函数跑完整段数据。外层循环比较简单我一般把参数轨迹存下来最后画收敛曲线。如果参数轨迹在中间有大的跳变说明 P0 或遗忘因子不合适还需要回调。na, nb, nc 2, 2, 1 n_theta na nb nc theta np.zeros(n_theta) P 2000 * np.eye(n_theta) lam 0.98 hist { na: na, nb: nb, nc: nc, y: list(y[:max(na, nb, nc)]), u: list(u[:max(na, nb, nc)]), yf: [0.0]*max(na, nb, nc), uf: [0.0]*max(na, nb, nc), ef: [0.0]*max(na, nb, nc), } theta_hist [] for t in range(max(na, nb, nc), N): theta, P, _ rml_step(y[t], u[t], theta, P, lam, hist) theta_hist.append(theta.copy()) theta_hist np.array(theta_hist) print(终值参数:, theta) print(真实参数:, np.concatenate([a_true, b_true, c_true]))跑完这段theta 的收敛终值应当非常接近 [-1.5, 0.7, 1.0, 0.5, 0.8]。如果偏差超过 5%优先查历史缓存是不是从同一时刻开始的——hist 里初始的 y 和 u 必须和递推起点一致否则前几百拍梯度方向全是错的。参数轨迹收敛的快慢也说明问题。a 和 b 通常领跑c 收敛慢一些这是正常现象因为 c 只能通过预测误差的“二次传播”被估计。看到 c 慢不用慌把数据扫第二遍就稳了。4.3 残差检验RML 到底把有色成分抽干没有参数对了只是第一步还得看残差序列是不是变成了白噪声。做法是把辨识出来的 A、B 参数代回去重新算预测误差然后检查残差的自相关函数。eps_final np.zeros(N) for t in range(max(na, nb, nc), N): eps_final[t] y[t] sum(theta[i]*y[t-i-1] for i in range(na)) \ - sum(theta[nai]*u[t-i-1] for i in range(nb)) # 残差自相关看 0 之外的值是否落入 ±1.96/sqrt(N) 置信带 acf np.correlate(eps_final - eps_final.mean(), eps_final - eps_final.mean(), modefull) acf acf[len(acf)//2:] / acf[len(acf)//2] plt.bar(range(1, 21), acf[1:21]) plt.hlines([1.96/np.sqrt(N), -1.96/np.sqrt(N)], 1, 20)如果 acf[1] 明显超出置信带说明有色成分没抽干净常见原因是 C 阶次 nc 给低了。直接把 nc 加一重跑很多问题就消失了。这一步是模型验证的核心比盯着参数曲线吹牛可靠得多。5. RML 实操避坑五个常见问题与对应的排查手段5.1 现象迭代几步之后全部变成 NaN程序直接崩原因最常见的是 P 矩阵发散也就是递推里 denom 过小增益 K 爆炸。这往往由两个叠加因素造成遗忘因子太小P 更新时除以 lam 让矩阵持续膨胀数据段里存在一段输入不变的区域回归量几乎为常数协方差矩阵快速病态。解决先把 lam 调到 0.98 以上P 更新公式保持 (P - Ppsi·Kᵀ)/lam 不变但检查输入序列是不是长期恒值。若数据本身激励不足任何算法都白搭重新设计实验或加随机扰动。还有一个细节如果数据有缺失值别用 0 填空缺失段直接跳过不参与更新。5.2 现象参数收敛了但收敛位置明显偏离真值原因RML 是伪线性回归似然函数可能存在多个局值尤其在初始参数选得不好时容易陷入次优点。另外P0 取得太小会让算法“不敢动”还没走到最优位置就停在原地。解决先用递推最小二乘RLS跑一遍把得到的参数当作 RML 初值然后重置 P 再扫第二遍。我习惯第一遍 P05000第二遍 P0200。如果还不够把数据切成长短不同的两段分别从不同初值启动 RML看收敛终值是否一致。多个终值里取似然最大的那个基本不会错。5.3 现象噪声参数 c 一直在变C(q) 的根出了单位圆原因迭代过程中 c 参数被推过了头C 多项式的某个根半径超过 1滤波环节 1/C(q) 就变成了不稳定滤波器yf、uf、ef 全部发散参数轨迹随后失控。解决每步更新后检查 np.roots(c) 的模长如果最大模长超过 0.99就把对应系数径向收缩。实际做法是把所有 c 系数乘一个统一比例因子例如把最大根模长拉到 0.98。这个“投影”操作在文献里叫稳定性投影工程里直接对系数下“重手”也能凑合关键是别让滤波状态先炸掉。5.4 现象A、B 参数看着正常残差检验却始终不通过原因参数没大问题但模型结构不对——最常见是纯延迟 d 没有估计或者噪声阶次 nc 给低了。残差里残留的正是没建模的那一拍延迟或高频噪声。解决先看残差自相关从第几拍开始越界。若 acf[1]、acf[2] 都越界先加 nc若 ar 和 ma 都加了还不行试延迟 d1 或 d2把 b 项改成 b·u(t-1-d)。这类问题在真实数据中出现频率极高不要一上来就怀疑算法。5.5 现象参数稳定但物理意义不对比如增益正负号反过来原因数据没有去均值或去趋势。现场数据几乎总有非零均值回归向量里没有常数项均值就被逼着分摊进了 a、b 系数里导致直流增益出现系统性偏差。解决辨识前先减掉 y 和 u 的均值如果趋势明显做一阶差分或先做 detrend。辨识完之后把均值加回去再计算静态增益用于验证物理方向对不对。这一步虽然简单却是我见过最多翻车的细节没有之一。6. 把 RML 结果拿去交付阶次选择、交叉验证与残差检验6.1 用信息准则把阶次定下来RML 要求给定 na、nb、nc但现场没人告诉你阶次是多少。我一般会用 AIC 或 BIC 扫一遍从 na1、nb1、nc1 起逐步升阶每升一次跑一遍 RML记录对数似然值。AIC 是 -2·lnL 2·n_thetaBIC 是 -2·lnL n_theta·lnN。阶次升高时 lnL 会改善但惩罚项也在变大选最小值对应的阶次。注意 AIC 偏向高阶BIC 更保守。我通常两个都看一致选 BIC 的结果AIC 当作“上限参考”。数据较短时BIC 选出低阶模型反而更稳。6.2 交叉验证训练段辨识验证段检验阶次定了之后还有些不踏实因为辨识那段数据的参数可能过拟合。我的做法是把数据切成两段前 60% 做辨识后 40% 只做验证。把训练段得到的固定参数代入验证段计算验证段残差的均方根和自相关。如果训练段残差小、验证段残差大说明模型过拟合了降阶次或加数据长度。验证段残差的自相关检验是最后一道关方法同 4.3 节置信带取 ±1.96/√N。全部通过这份参数才能拿出去用。这套流程我从第一份系统辨识作业起就在用。后来做现场数据也靠着这几步躲开了不少麻烦——最惨的一次是在现场拿到一批带强均值漂移的数据没去均值直接辨识结果增益符号反了被验收方当众指出来。从那以后我每次跑 RML 前都强制过一遍“去均值、扫阶次、残差白化检验”三板斧再没出过这种低级问题。希望帮到你。本文还有配套的精品资源点击获取