递推最小二乘参数辨识：从RLS偏差到渐消记忆法改进

简介：面向控制系统、信号处理与机器学习方向的学习者与工程人员，这份资源提供了一套最小二乘类参数辨识的MATLAB实现，用于对系统参数进行在线估计与辨识。程序以极点位于单位圆内的稳定系统为例，先产生500点M序列作为输入，叠加方差为1的白噪声扰动生成观测序列，再分别用递推最小二乘与渐消记忆最小二乘完成参数估计，并通过增益矩阵与估计方差递推更新，最后绘制参数估计过程曲线，便于对比两种算法的收敛效果与偏差差异。压缩包共1个docx文件，约78KB，以文档形式整理程序代码、设计思路与辨识结果说明，适合边看边在MATLAB中复现。该资源已有2161人学习，可作为参数辨识入门练习与算法对比的参考素材。

1. 从一组“收敛了但不对”的曲线说起

跑过系统辨识的人大概都遇到过这种场面：递推最小二乘的曲线明明收敛了，四条参数轨迹稳稳停住，可拿真值一对比，稳态误差就摆在那儿，怎么调初值都消不掉。这份最小二乘类的参数辨识 MATLAB 程序，讲的就是这么一件事——先用递推最小二乘（RLS）把一组二阶系统的参数估出来，发现结果有偏，再换成渐消记忆（遗忘因子）最小二乘，把偏差压下去。它适合正在学系统辨识、需要一份能直接跑通并对照修改的 MATLAB 代码的开发者，也适合做控制、信号处理方向、想快速验证辨识算法效果的人。整份程序不依赖任何工具箱，纯脚本，复制进 MATLAB 就能出图。

2. 系统模型与数据生成：先把“真值”造出来

辨识类程序有个绕不开的前提：你得先知道真值，才能判断估得准不准。这份代码的做法是先人为设定一个二阶差分方程，用已知系数生成观测序列，再让算法去“猜”这些系数。猜得对不对，一眼就能看出来。

2.1 被辨识系统的差分方程

程序里设定的系统是：

Z(k+2) = 1.8*Z(k+1) - 0.9*Z(k) + 0.75*u(k+1) + 0.5*u(k) + v(k)

对应到标准形式A(z⁻¹)Z(k) = B(z⁻¹)u(k) + v(k)，待辨识的参数向量是[a1, a2, b1, b2] = [-1.8, 0.9, 0.75, 0.5]。注意符号：差分方程里是+1.8*Z(k+1)，移到等式左边变成-1.8，所以a1 = -1.8。这个符号约定是后面构造回归向量h的关键，写反了曲线会朝反方向跑。

特征方程z² - 1.8z + 0.9 = 0的根模长为√0.9 ≈ 0.949 < 1，极点落在单位圆内，系统稳定。这一点很重要：如果系统本身发散，观测序列z会指数增长，递推过程数值上直接崩掉，辨识根本无从谈起。所以设计辨识实验时，先确认对象稳定，是血泪经验换来的第一步。

2.2 M 序列激励信号的生成

辨识要“激励”系统，输入信号必须持续激励，否则某些参数根本不可观。程序用 9 级线性反馈移位寄存器产生 M 序列（伪随机二进制序列），周期2⁹ - 1 = 511，取前 500 点：

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); % 输出取末级 for j = 9:-1:2 x(j) = x(j-1); % 寄存器整体右移 end x(1) = temp; % 反馈值填入首位 end

逻辑说明：xor(x(4), x(9))是反馈多项式，抽头位置决定了序列的周期和随机性；内层循环做移位，把每一位往后推一格；M(i)取移位前的末级作为输出。参数上，N=500是样本数，够 RLS 收敛；n=503比N多留几点，是因为后面生成z时要用到M(i-1)、M(i-2)，索引会往前错位，多算几点防止越界。

提示：寄存器初值不能全零，全零状态下异或永远是零，序列退化成常数，激励就失效了。这是新手最容易翻车的地方之一。

2.3 白噪声扰动与观测序列合成

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

randn产生标准正态分布，方差为 1，正好对应正文说的“方差为 1 的白噪声”。z开N+2长度，是因为差分方程是二阶的，前两个点作为初始条件，从第 3 点开始递推。v(i-2)的索引错位是为了让噪声和z(i)对齐——z从 3 开始，v从 1 开始，差 2。初始条件z(1)=-1, z(2)=0是任意给的，只要系统稳定，初始条件的影响会随时间衰减，对稳态辨识结果没影响。

3. 递推最小二乘核心：四行代码里的估计逻辑

数据造好之后，重头戏是 RLS 的递推循环。这段代码不长，但每一行都对应一个明确的数学含义，改错一个符号或者一个转置，结果就完全变样。

3.1 回归向量与参数初值

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]; % 增益矩阵初值

P是 4×4 的估计误差协方差矩阵，初值100*eye(4)表示对初始参数“完全没把握”，让算法快速信任新数据。Theta初值给[1;1;1;1]，和真值[-1.8; 0.9; 0.75; 0.5]差得很远，正好用来观察收敛过程。Pstore只存P的对角线，因为对角元素才代表各参数的估计方差，画图时看它收敛到零就说明估计越来越确定。

3.2 递推更新循环

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

逐行拆解：

• h是回归向量，由过去两个输出和两个输入组成。注意z前面带负号，因为差分方程移项后a1*Z(k+1)变成-a1*Z(k+1)，回归向量里放-z才能让Theta直接对应[a1, a2, b1, b2]。
• K是增益矩阵，h'*P*h + 1里的1是噪声方差（这里假设为 1）。如果噪声方差不是 1，这个位置要换成实际方差，否则增益尺度不对，收敛速度会偏。
• Theta更新用的是“新息”z(i) - h'*Theta(:,i-2)，即当前观测减去预测值，再乘增益修正参数。
• P更新是标准的 Riccati 递推，(eye(4) - K*h')*P保证协方差单调递减。

索引上有个细节：Theta(:,i-1)对应的是用第i个观测更新后的结果，所以画图时i=1:N+1正好对上Theta的列数。这种错位在 MATLAB 里很常见，调试时打印几个点核对一下最稳妥。

3.3 结果可视化与偏差现象

i = 1:N+1; figure(1) plot(i, Theta(1,:), i, Theta(2,:), i, Theta(3,:), i, Theta(4,:), ... i, -1.8, i, 0.9, i, 0.75, i, 0.5) grid on; title('参数估计过程');

四条实线是估计值，四条水平线是真值参考。跑出来会看到：曲线确实收敛了，但稳态值和参考线之间始终差一截。这不是代码写错了，而是标准 RLS 在有色噪声或者模型阶次不匹配时的固有偏差——这里噪声虽然是白的，但递推过程中旧数据权重不衰减，早期误差累积下来就形成了稳态偏置。正文里那句“辨识结果会与其精确值始终着偏差”，说的就是这个现象。

4. 渐消记忆最小二乘：两行改动带来的精度提升

发现偏差之后，正文给的解法非常直接：把遗忘因子引进来，只改两行。这个改动小到容易被忽略，但背后的思路值得说清楚。

4.1 遗忘因子的引入方式

% 原版 K = P*h*inv(h'*P*h + 1); P = (eye(4) - K*h')*P; % 渐消记忆版 K = P*h*inv(h'*P*h + 0.998); P = (eye(4) - K*h')*P / 0.998;

0.998就是遗忘因子 λ。它的作用是在每一步把旧协方差“放大”一点（除以小于 1 的数等于放大），让算法逐渐“忘掉”老数据，更重视新观测。λ 越接近 1，记忆越长，越接近标准 RLS；λ 越小，跟踪时变参数的能力越强，但稳态噪声也越大。0.998 对应的时间常数大约是1/(1-0.998) = 500步，和样本数N=500量级匹配，所以在这个例子里效果刚好。

4.2 参数选择与效果对比

改完之后重跑，四条估计曲线会更贴近真值参考线。但要注意，渐消记忆不是万能的：如果 λ 取得太小（比如 0.9），曲线会在真值附近抖动，稳态方差反而变大；如果系统是时变的，λ 小一点合理，但如果是定常系统，λ 太小学要付出噪声代价。我一般会先跑标准 RLS 看偏差量级，再拿 0.99、0.995、0.998 各试一遍，挑收敛和抖动折中的那个。

注意：改P更新那行时，除号的位置别写错。写成P = (eye(4) - K*h')*P*0.998是反的，协方差会越缩越小，算法提前“僵死”，参数停在半路不动。

5. 避坑与排查：跑不出正确曲线时先看这几处

这份代码短，但短代码的坑往往更隐蔽。下面几条是实际调试时最常撞上的，按“现象 → 原因 → 解决”整理。

5.1 曲线发散或出现 NaN

现象：Theta某几列突然冲到1e10量级，或者图上是NaN。

原因：inv(h'*P*h + 1)里的矩阵接近奇异，通常是P初值给得过大（比如1e6*eye(4)）加上早期h接近零向量，数值上求逆不稳定。

解决：P初值控制在100到1000量级；如果还不行，把inv换成pinv或者用左除K = P*h/(h'*P*h + 1)，MATLAB 的左除对病态矩阵更稳。

5.2 参数收敛到错误符号

现象：估计曲线收敛了，但a1停在+1.8附近而不是-1.8。

原因：回归向量h里z的符号写反了，或者差分方程移项时符号约定没统一。

解决：回到 2.1 节，确认h = [-z(i-1); -z(i-2); ...]的负号在，且Theta初值和真值参考线的符号一致。符号问题没有玄学，就是逐项对。

5.3 遗忘因子改了但结果没变

现象：把1改成0.998之后，曲线和原来几乎一样。

原因：只改了K那行，忘了改P更新那行的/0.998；或者改完没清工作区，跑的还是旧变量。

解决：两行必须成对改；每次跑之前clear一下，或者把脚本包成函数，避免工作区残留变量干扰。

5.4 M 序列周期不够导致激励不足

现象：某些参数收敛特别慢，或者P对角线降不下去。

原因：N=500接近 M 序列周期 511，尾部数据相关性变强，激励的“ richness ”下降。

解决：把寄存器级数加到 10 级（周期 1023），或者换用正弦叠加、随机二进制序列做输入，保证在整个样本区间内持续激励。

5.5 噪声方差与增益分母不匹配

现象：收敛速度明显偏慢，或者稳态抖动比预期大。

原因：h'*P*h + 1里的1假设噪声方差为 1，但实际randn的方差确实是 1，如果换成0.01*randn而分母没跟着改，增益尺度就错了。

解决：噪声方差改多少，分母常数就跟着改多少；渐消记忆版则是λ和噪声方差共同决定分母，调的时候一起看。

6. 二次开发：从固定阶次到可配置辨识框架

这份代码拿来直接跑没问题，但要复用到别的系统上，硬编码的地方得先拆掉。我的习惯是把脚本改成一个带参数的函数，把阶次、遗忘因子、样本数都提出来，这样换对象时只改配置不改逻辑。

6.1 参数化改造思路

function [Theta, Pstore] = rls_identify(z, u, na, nb, lambda, P0) % z: 观测序列, u: 输入序列 % na, nb: 输出/输入阶次 % lambda: 遗忘因子，1 表示标准 RLS % P0: 协方差初值 n = na + nb; N = length(z) - max(na, nb); Theta = zeros(n, N+1); Theta(:,1) = ones(n,1); P = P0 * eye(n); Pstore = zeros(n, N+1); for k = max(na,nb)+1 : length(z) h = [-z(k-1:-1:k-na); u(k-1:-1:k-nb)]; K = P*h / (h'*P*h + lambda); Theta(:,k-max(na,nb)) = Theta(:,k-max(na,nb)-1) + ... K*(z(k) - h'*Theta(:,k-max(na,nb)-1)); P = (eye(n) - K*h')*P / lambda; Pstore(:,k-max(na,nb)) = diag(P); end end

改造点有三个：阶次na、nb变成入参，回归向量用切片z(k-1:-1:k-na)自动拼，不再手写四行；lambda统一控制标准/渐消两种模式，传 1 就是标准 RLS；P0暴露出来方便调初值。这样换一个三阶系统，只改调用时的na=3, nb=2就行。

6.2 验证改造是否正确的办法

改完之后别急着信结果，拿原例子的数据回代：na=2, nb=2, lambda=1跑一遍，和原脚本的Theta逐列对比，误差应该在1e-10量级；再传lambda=0.998，看稳态偏差是否比标准版小。两步都对上，说明参数化没引入逻辑错误。我一般还会故意传一个lambda=0.9，看曲线是否如预期那样抖动变大——如果没变化，说明lambda根本没进到更新式里，八成是某处还写着硬编码的1。

6.3 一个容易忽略的索引习惯

原脚本里Theta(:,i-1)和Theta(:,i-2)的错位，在参数化之后变成了k-max(na,nb)和k-max(na,nb)-1，看着绕，但本质是“当前列存更新后的值，前一列是更新前的值”。我踩过的坑是：改造时图省事把索引统一成k，结果第一列永远是初值、最后一列越界，画出来曲线整体右移一格。从那以后我每次改索引，都先打印size(Theta)和前三个点的值核对，确认列数和样本对齐再往下跑。希望这份拆解能帮你少走一遍我走过的弯路，把这份最小二乘参数辨识程序跑通、改顺。

本文还有配套的精品资源，点击获取