简介：这是一篇面向控制系统、信号处理与机器学习等领域学习者的最小二乘参数辨识程序说明文档。资源共包含一个docx格式文档，压缩包大小七十八KB，内容围绕递推最小二乘算法和渐消记忆最小二乘算法展开：先构造一个极点位于单位圆内的稳定系统，生成五百点M序列作为输入，并叠加白噪声扰动；随后通过增益矩阵与估计方差矩阵的递推更新，实现系统参数的在线估计与辨识，同时借助绘图功能将参数估计过程可视化。文档还针对渐消记忆法给出关键改进，通过引入遗忘因子削弱旧数据影响，使估计结果更接近真实值；其中包含完整的递推公式、参数初值设定、循环迭代流程以及仿真图形输出，便于读者直接对照学习或改造成其他辨识对象。已有两千一百五十八人学习下载，适合需要掌握最小二乘类参数辨识原理、递推算法实现及仿真验证方法的初学者和相关技术开发人员。
1. 最小二乘参数辨识：一套能直接改着用的 MATLAB 递推程序
做控制系统或信号处理的人，大概率都遇到过这种事：系统结构大致清楚，但里面的系数拿不准。传递函数写不出来，状态空间更是无从下手。这时候最直接的办法就是给系统一个已知的激励信号，记录输出，再用参数辨识把模型里的未知系数一点一点“逼”出来。这篇笔记要拆的是一份操作型很强的 MATLAB 程序——它用最小二乘递推算法（RLS）辨识一个二阶系统的四个参数，然后给出渐消记忆最小二乘的改进版本，顺带把 M 序列激励、白噪声扰动、增益矩阵更新这些细节全部摊开来讲。适合刚接触辨识理论但想尽快跑通一个完整例子的从业者，也适合那些已经在用工具箱、但想搞明白内部迭代逻辑的人。
这程序里辨识的对象是：Z(k+2)=1.8Z(k+1)-0.9Z(k)+0.75u(k+1)+0.5u(k)+v(k)。四个待辨识参数分别是 -1.8、0.9、0.75、0.5（注意程序里 h 向量带负号，所以实际迭代目标是 -1.8 和 0.9 这两个带符号的量）。极点都在单位圆内，系统稳定，所以数据是收敛可用的。M 序列做输入、白噪声做扰动，这是辨识里最常见的测试组合。下面从算法原理到每一段代码怎么改，一步步来过。
2. 递推最小二乘的数学内核：为什么是 P、K、Theta 三个量在滚
2.1 最小二乘的核心思想与递推化的必要性
最小二乘参数辨识的基本思路很简单：把系统写成 y(k) = h'(k) * θ + v(k) 的形式，其中 h(k) 是回归向量（由过去的输入输出组成），θ 是待辨识参数向量，v(k) 是噪声。我们的目标是找一组 θ，让残差平方和最小。批量最小二乘可以直接用矩阵求逆算出 θ=(H'H)^(-1)H'Z，但问题是：如果数据是一批一批到的，或者系统参数本身可能缓慢漂移，每次都重新做一次矩阵求逆，计算量会随数据量增长而爆炸，而且老数据会把新数据的信息“稀释”掉。
递推最小二乘就是为了解决这个矛盾。它把批处理公式拆成一个迭代式子：每来一个新数据，就用增益矩阵 K(k) 对当前的参数估计 Theta(k-1) 做一个修正。修正量等于 K(k) 乘以预测误差（实际输出减去当前参数预测的输出）。这个过程中有三个量在同步滚：参数估计值 Theta、误差协方差矩阵 P、增益矩阵 K。P 矩阵本质上记录了参数估计的不确定度，K 则是把 P 和当前回归向量 h 组合起来算出的修正步长。三者配合，理论上可以从任意初值出发，逐步逼近真实参数。
这个程序用的系统方程是 Z(k+2)=1.8Z(k+1)-0.9Z(k)+0.75u(k+1)+0.5u(k)+v(k)。把它改写成最小二乘的标准形式，回归向量 h(k)=[-z(k-1); -z(k-2); u(k-1); u(k-2)]（程序中循环从 i=3 开始，所以 h 里取的是 z(i-1)、z(i-2)、M(i-1)、M(i-2)，对应前面的 z(k+1)、z(k)、u(k+1)、u(k)），待辨识参数 θ=[a1; a2; b1; b2]，其中 a1=-1.8、a2=0.9、b1=0.75、b2=0.5。注意这里回归向量里的 z 和 u 都带了负号或正号，符号问题在后面看代码时要特别留神——很多人自己写的时候就是在这里把符号搞反，导致辨识出来的参数符号全错。
2.2 RLS 迭代公式与程序中的变量一一对应
RLS 的标准迭代公式有三个：
增益矩阵：K(k) = P(k-1) * h(k) / (h'(k) * P(k-1) * h(k) + λ)。λ 是遗忘因子，标准 RLS 取 λ=1，渐消记忆法取 λ<1（如 0.998）。
协方差更新：P(k) = (I - K(k) * h'(k)) * P(k-1) / λ。
参数更新：Theta(k) = Theta(k-1) + K(k) * (z(k) - h'(k) * Theta(k-1))。
程序里的P=100*eye(4)对应 P 的初始值，Theta(:,1)=[1;1;1;1]对应参数初值，K=[1;1;1;1]只是初始化一个列向量，在循环里会被真正计算的值覆盖。h=[-z(i-1);-z(i-2);M(i-1);M(i-2)]就是当前时刻的回归向量。循环体里三行代码和上面三个公式一一对应。P 的初值取得比较大（100 倍单位阵），表示一开始对参数估计非常不确定，这样 K 的初始修正步长比较大，参数能快速从初值 [1;1;1;1]（这是一个故意给偏的初值）收敛到真值附近。
这里有个值得注意的细节：程序里 h'Ph + 1 中的 1 是噪声方差假设。如果实际系统噪声方差不是 1，这个 1 应该相应调整。但在这个例子里，v=randn(1,N) 产生的白噪声方差恰好是 1，所以不用改。很多人直接把工具箱里的 RLS 拿过来用，根本没意识到这个数字和噪声方差挂钩，等换了系统换了噪声，辨识结果变差了却找不到原因——这就是其中一条隐蔽的坑。
3. 程序逐段拆解：M 序列、白噪声、观测数据和递推主循环
3.1 M 序列产生：移位寄存器与异或反馈的配合
程序开头这段是产生 M 序列（伪随机二进制序列）作为输入激励：
x=[0 1 0 1 1 0 1 1 1]; % 寄存器初值（9级） n=503; % 脉冲总数 N=500; % 实际使用的点数 for i=1:n temp=xor(x(4),x(9)); % 第4级与第9级异或，作为反馈 M(i)=x(9); % 输出取第9级 for j=9:-1:2 x(j)=x(j-1); % 寄存器右移一位 end x(1)=temp; % 反馈值写入第1级 end逻辑说明：这是一个 9 级线性反馈移位寄存器（LFSR），反馈多项式取第 4 级和第 9 级的异或。每来一个时钟脉冲，寄存器内容右移一位，最高位丢弃，最低位补上反馈值。输出的 M 序列取的是第 9 级的值。这个序列是周期性的，周期为 2^9 - 1 = 511。程序产生 503 个点，取前 500 个使用，所以实际用的是一整个周期内的连续片段。
参数说明：寄存器初值 x 不能全零，否则序列永远输出 0，这一点要特别注意。n 设为 503 而 N 设为 500，多出来的 3 个点是留给后面生成观测数据时用的——观测数据 z 需要从 i=3 开始循环，所以输入序列必须多准备几个点才能让下标对齐。反馈抽头选择第 4 和第 9 级不是随意的，9 级寄存器用这对抽头能保证序列周期达到最大长度 511。如果你要改寄存器级数，抽头位置也要跟着变，否则周期会缩短，激励信号的频谱覆盖就不够均匀。
3.2 白噪声与观测数据生成：验证模型时的关键路径
v=randn(1,N); % 方差为1的白噪声 z=zeros(N+2,1); z(1)=-1; z(2)=0; for i=3:N+2 z(i)=1.8*z(i-1)-0.9*z(i-2)+0.75*M(i-1)+0.5*M(i-2)+v(i-2); end逻辑说明：z 是系统的观测输出序列。z(1)=-1、z(2)=0 是两个初始状态，不是为了凑数，而是为了让循环从 i=3 开始时有前两个输出值可用。第 i 个观测值由前两个输出、前两个输入和当前噪声共同决定。公式里的系数 1.8、-0.9、0.75、0.5 就是我们要辨识的真值。注意噪声项是 v(i-2)，和输出 z(i) 对齐——因为回归向量里的输入输出都是滞后项，噪声也要对应滞后两个时刻才合理。
参数说明：randn(1,N) 产生标准正态分布随机数，均值 0 方差 1，正好匹配 RLS 公式里分母上的“1”。如果你把噪声方差调大，比如改成 2*randn(1,N)，那么 RLS 公式里的那个 1 也要相应改成 4（方差值），否则增益矩阵的步长和噪声水平不匹配，参数估计的波动会明显变大。程序里 z 的长度是 N+2，比 M 序列多两个点，也是为了让循环下标从 3 开始后能恰好用到 M 序列的全部 500 个点。
3.3 递推主循环：P、K、Theta 三个量的更新顺序是成败关键
P=100*eye(4); % 估计方差初值，给大 Pstore=zeros(4,N+1); Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; Theta=zeros(4,N+1); % 参数估计值存储 Theta(:,1)=[1;1;1;1]; % 故意给偏的参数初值 K=[1;1;1;1]; % 增益矩阵初值（占位） for i=3:N+2 h=[-z(i-1);-z(i-2);M(i-1);M(i-2)]; % 回归向量 K=P*h*inv(h'*P*h+1); % 增益矩阵 Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2)); % 参数更新 P=(eye(4)-K*h')*P; % 协方差更新 Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; % 记录方差对角元 end逻辑说明：循环里三步依次是算 K、更新 Theta、更新 P，顺序是固定的——必须先用旧的 P 算 K，再用 K 修正 Theta，最后才更新 P。换顺序会导致 P 和 K 不匹配，参数收敛轨迹会乱掉。存储数组 Pstore 和 Theta 的列数都是 N+1（从 1 到 501），但循环体里 Theta(:,i-1) 的下标范围是 2 到 501，所以初值放在第 1 列，最后一个估计值放在第 501 列。画图时 i=1:N+1 正好覆盖全部估计点。
参数说明：P 的初值取 100eye(4)，表示初始协方差为 100，即参数不确定度很高。这个值如果取太小（比如 1），参数收敛会明显变慢；取太大（比如 10^6）则前几步的参数估计会剧烈震荡。Theta 初值取 [1;1;1;1] 离真值 [-1.8;0.9;0.75;0.5] 有不小距离，这样能直观看到收敛过程。回归向量 h 里的符号很关键：因为差分方程是 z(k+2)=1.8z(k+1)-0.9*z(k)+...，移到等式左边后，z 项的系数分别变成 -1.8 和 +0.9，所以 h 里对应位置是 -z(i-1) 和 -z(i-2)，这样辨识出来的 θ(1) 和 θ(2) 会收敛到 -1.8 和 0.9。
4. 渐消记忆最小二乘：遗忘因子 0.998 背后的原理和改法
4.1 为什么标准 RLS 会有稳态偏差
标准 RLS 的推导里，P 矩阵会随着迭代次数增加而不断缩小。这意味着后期增益矩阵 K 的模越来越小，参数的修正步长也随之变小。理论上如果系统是线性定常的、噪声是零均值的，参数会收敛到真值附近，但实际仿真中你往往会发现估计值和真值之间始终存在一个小的偏差，而且这个偏差不随迭代次数增加而消失。原因在于：RLS 的代价函数对所有历史数据等权重看待，越老的数据积累越多，在总代价里占比越大，新数据对参数修正的贡献相对被稀释。如果系统本身有轻微时变，或者噪声不是纯白噪声（比如有色噪声），这个偏差会更明显。
把系统输出和参数估计代入程序画出的图上，你会发现四条参数曲线都收敛了，但每条曲线和真值之间都有一条肉眼可见的“缝”。这个缝的大小取决于噪声水平、P 初值和数据长度。P 初值取 100、数据取 500 点时，偏差通常较小，但不为零——这正是程序注释里写的“辨识结果会与其精确值始终存在着偏差”。
4.2 改成渐消记忆法的关键：两个语句而已
渐消记忆法的思路是给历史数据加一个指数衰减的权重，越老的数据权重越小，这样新数据能在代价函数中占据更大比例，从而持续修正参数。实现上只需要改两处：
% 原语句 K=P*h*inv(h'*P*h+1); P=(eye(4)-K*h')*P; % 改为渐消记忆法（遗忘因子 λ=0.998） K=P*h*inv(h'*P*h+0.998); P=(eye(4)-K*h')*P/0.998;逻辑说明：分母上的 1 改成 0.998，对应遗忘因子 λ。P 更新时除以 0.998（即乘以 1/0.998 ≈ 1.002），效果是让 P 矩阵不会像标准 RLS 那样快速收缩到接近零。这样每步结束后 P 维持在一个“有弹性”的量级，下一个数据点到来时 K 仍然有足够大的模，参数修正步长不会衰减到微乎其微。综合效果就是：估计值能持续跟随真值，稳态偏差明显缩小。
参数说明：遗忘因子 λ 的取值范围通常是 0.95～1.0。λ 越接近 1，记忆越长，参数越稳定但跟踪越慢；λ 越小（比如 0.95），记忆越短，跟踪时变参数的能力越强，但稳态波动也越大。0.998 是一个比较保守的值，适合参数基本不变但希望消除稳态偏差的场景。如果你的系统参数是缓变的，可以考虑取 0.98～0.99；如果参数变化较快，可能需要 0.95～0.97，但这时要考虑噪声放大的副作用。注意在这里噪声方差假设从 1 变成了 λ=0.998，所以分母的含义不只是噪声方差，而是噪声方差与遗忘因子的组合——这个细节很多人会忽略。
4.3 两种算法的结果对比怎么看
对比标准 RLS 和渐消记忆法的参数估计曲线，最直观的差异在稳态段。标准 RLS 的曲线虽然收敛，但每条曲线和真值之间始终有偏差；渐消记忆法的曲线几乎贴近真值线。P 矩阵的对角元变化也解释了这一点：标准 RLS 的 P 对角元持续衰减到接近零，渐消记忆法的 P 对角元稳定在一个较小的非零值附近。这个差异在你需要高精度参数值时非常关键——比如你辨识出来的参数接下来要用于设计控制器，偏差可能导致闭环极点偏出期望位置。
5. 避坑指南：从初值设定到噪声匹配，四个最容易翻车的地方
5.1 现象：P 初值给太小导致收敛极慢
现象：把 P=100eye(4) 改成 P=eye(4) 或 P=0.1eye(4)，参数估计曲线需要很长时间才能逼近真值，500 个数据点可能不够用。原因：P 初值代表对参数初始不确定度的估计，给小了等于默认“初始参数已经很准”，增益矩阵 K 的初始步长就小，参数只能慢慢爬过去。解决：保持相对较大的初值（如 100 倍单位阵），或者根据参数可能的变化范围手动估计一个合适的数量级。我一般的做法是先跑一次标准 RLS 看 P 对角元收敛到多少，再反推一个合理的初值，一般取收敛值的 10～100 倍。
5.2 现象：遗忘因子太小导致参数估计波动大
现象：把遗忘因子从 0.998 改成 0.95，参数曲线虽然在真值附近波动，但波动的幅度明显变大，看起来“毛刺”很多。原因：遗忘因子越小，有效数据窗口越短，等效于用更少的数据做估计，方差自然增大。解决：先试 0.99～0.998 这个区间，确认稳态偏差确实被消除后再尝试适当减小。如果偏差还在，先检查是不是噪声方差没匹配，而不是一味调小遗忘因子。我自己的习惯是：只在对参数有明确时变需求时才把遗忘因子调到 0.98 以下，否则尽量保持在 0.99 附近。
5.3 现象：回归向量的符号写反，辨识出的参数符号全错
现象：把 h=[-z(i-1);-z(i-2);M(i-1);M(i-2)] 写成 h=[z(i-1);z(i-2);M(i-1);M(i-2)]，结果 θ(1) 收敛到 1.8 而不是 -1.8，θ(2) 收敛到 -0.9 而不是 0.9。原因：差分方程移项后系数符号变了，但回归向量没有同步调整。解决：先把差分方程写成标准回归形式 y(k)=h'(k)*θ+v(k)，逐项确认符号，再写代码。这一步最基础但最容易错，而且错的隐蔽——最终参数的绝对值是对的，正负号却是反的。
5.4 现象：M 序列初值全零导致激励信号为零，辨识完全失败
现象：寄存器初值改成 x=[0 0 0 0 0 0 0 0 0]，程序运行后参数曲线完全不动，或者发散到异常值。原因：全零初值的 LFSR 永远输出 0，系统完全没有被激励，观测数据里不含任何参数信息。解决：确保寄存器初值至少有一个 1，最省事的就是保持程序原来的初值 [0 1 0 1 1 0 1 1 1]。如果你自己重新生成 M 序列，记得在代码里加一个判断：如果 sum(x)==0，就报错提示重新初始化。
5.5 现象：噪声方差与 RLS 公式中的常数不匹配，估计值偏离真值
现象：把噪声改成 v=2*randn(1,N)，但 RLS 公式里的常数保持 1，最后参数估计偏离真值，且曲线波动变大。原因：公式里的常数同时承担了噪声方差设定和遗忘因子两个角色，噪声实际方差和假设值不一致时，增益矩阵的步长不是最优的。解决：如果你改了噪声方差，记得同步修改inv(h'*P*h+常数)里的常数，改成噪声方差的数值。这是最容易被忽略的一个参数，因为它藏在矩阵求逆的表达式里，不是单独一句话。
6. 把辨识结果落地：残差校验 + 利用带遗忘因子的 RLS 做在线跟踪
当你把渐消记忆法的程序跑通、参数曲线也收敛到真值附近之后，还有一个值得做的验证步骤：检查模型输出和实际观测之间的残差。做法很简单，把辨识出的参数代回系统方程，用同一段输入重新仿真一遍，然后算残差：
% 辨识完成后，用估计参数计算模型输出 z_est = zeros(N+2,1); z_est(1) = -1; z_est(2) = 0; for i=3:N+2 z_est(i) = -Theta(1,end)*z_est(i-1) - Theta(2,end)*z_est(i-2) ... + Theta(3,end)*M(i-1) + Theta(4,end)*M(i-2); end residual = z(3:N+2) - z_est(3:N+2); % 残差序列 figure(2) subplot(2,1,1) plot(1:N, residual) title('模型残差序列') subplot(2,1,2) histogram(residual, 30) title('残差直方图')逻辑说明：残差等于实际输出减模型输出。如果辨识效果好，残差应该接近白噪声——均值接近零、直方图呈钟形。如果残差里还有明显的周期性或趋势，说明模型结构没对（比如阶次不够高），或者输入信号没有充分激励出系统的所有模态。
参数说明：注意这里模型输出用 Theta(1,end)、Theta(2,end) 取的是最后一列的估计值。Theta(1,end) 是辨识出的 -a1（因为回归向量里 h(1)=-z(i-1)），所以计算模型输出时要写成-Theta(1,end)*z_est(i-1)。这个符号细节和前面回归向量里的符号是呼应的，别在这里搞反。
一个我自己常用的习惯是：在线辨识场景里不只要画参数曲线，还要同时画 P 矩阵的迹（或对角元和），因为它能直观反映“当前参数可信度”。如果迹在迭代后期还在持续下降，说明系统还没进入稳态；如果迹已经稳住但有波动，说明遗忘因子在起作用。从那以后我每次跑辨识实验都会强制走一遍“观察 P 迹 → 检查残差 → 确认符号”这三步，再把结果用于后续设计。这一步看着多花几分钟，实际能挡掉大量后期的低级返工。希望这份程序和拆解思路能帮到你，尤其是那些刚开始接触递推辨识、想尽快把理论落到代码上的朋友——按上面的步骤改参数、看曲线、检查残差，你很快就能摸清这套算法的手感。
本文还有配套的精品资源，点击获取
