递推最小二乘参数辨识MATLAB程序详解与渐消记忆改进

简介：最小二乘类参数辨识MATLAB程序包，面向自动化、控制理论与信号处理方向的学习者与科研人员，解决系统参数在线估计与辨识的实现问题。程序基于递推最小二乘与渐消记忆最小二乘算法，配合500点M序列输入和白噪声扰动，在稳定系统模型上演示了参数初值设定、增益矩阵K与方差P的递推更新及结果展示。程序包为1个docx文档，容量78KB，含完整代码、关键注释和辨识结果分析，并解释了渐消因子0.998下K与P公式的调整方式，便于对比两种算法的精度与收敛差异。文档还介绍白噪声、增益矩阵、估计方差等概念，提供参数初值选取经验，便于理解数学基础。已有2158人学习使用，适合掌握最小二乘辨识原理并开展MATLAB编程实践，也可用于课程设计或论文仿真。

1. 最小二乘参数辨识：这份MATLAB程序能帮你省掉哪些推导

做控制、信号处理的人，迟早会碰到参数辨识：系统模型有了，但里面的系数不知道，要用输入输出数据把它估出来。我最初自己推RLS公式时，被增益矩阵和协方差更新绕得晕头转向，直到拿到一份能直接跑的最小二乘类MATLAB程序，才真正理解递推过程。这份程序用500点M序列作为激励，对二阶系统做递推最小二乘，并演示了渐消记忆法修正稳态偏差。它适合刚刚接触辨识理论、想用MATLAB验证算法、或者在课程设计中需要完整代码的同学。下面我把这段代码拆开讲清楚，包括参数怎么设、哪里容易踩坑。

2. 从系统模型到递推最小二乘：为什么用M序列和500个数据点

2.1 系统方程与辨识目标：先把对象写清楚

程序里设计的系统是： Z(k+2)=1.8Z(k+1)-0.9Z(k)+0.75u(k+1)+0.5u(k)+v(k)

其中v(k)是白噪声。这个方程把未来输出表示成过去两拍输出、过去两拍输入和噪声的组合。实际生成观测数据时，程序用的是z(i)=1.8z(i-1)-0.9z(i-2)+0.75M(i-1)+0.5M(i-2)+v(i-2)，两者只是下标平移了两个周期，本质相同。把方程改写成回归形式：

z(k) = [-z(k-1); -z(k-2); u(k-1); u(k-2)]' * theta + v(k-2)

其中theta = [-1.8; 0.9; 0.75; 0.5]。注意h向量里带负号，待估参数是-1.8而不是1.8。程序最后画图时，拿Theta(1,:)和-1.8比，别搞反。

系统稳定性可以验证：特征方程是λ^2 -1.8λ +0.9=0，判别式为3.24-3.6=-0.36，两个共轭复根模长sqrt(0.9)≈0.9487，位于单位圆内，所以系统稳定。稳定是递推辨识能收敛的前提之一。如果极点靠近单位圆甚至在外，输出会发散，或者需要很长的数据才能保证持续激励。这里极点模长0.9487已经比较接近1，实际上收敛速度不算快，早期瞬态会影响初期的参数估计，所以用500点数据是有道理的。

参数辨识的目标，就是用输入M(i)和输出z(i)序列，把theta里的四个数估出来。程序里没有直接给出真值，而是通过初始Theta=[1;1;1;1]和递推公式逐步逼近。这个“从错误初值逼近真值”的过程，恰好就是RLS最直观的教学演示。如果你把初值改成[0;0;0;0]再跑一遍，会发现收敛过程相似，只是路径不同。

2.2 M序列与白噪声：激励信号怎么选

系统辨识需要输入信号能激发系统的所有动态模态。单位阶跃或单一正弦往往只能激励部分频率，参数估计会不稳定。程序用的是M序列，即伪随机二进制序列。M序列的频谱接近白噪声，且幅值恒定，属于持续激励信号，能保证参数可辨识。生成代码如下：

x=[0 1 0 1 1 0 1 1 1]; %寄存器初值，不能全零 n=503; %总共产生503个点 N=500; %实际使用前500点 for i=1:n temp=xor(x(4),x(9)); M(i)=x(9); %输出M序列 for j=9:-1:2 x(j)=x(j-1); end x(1)=temp; end

x是9位移位寄存器，反馈抽头在第4位和第9位，异或后送回第1位。这个配置对应的周期是2^9-1=511，取前500点。注意循环产生了503个点，比N多3个，这是给后续生成观测序列时留的余量，避免下标越界。

白噪声用randn(1,N)生成，均值为0、方差为1的标准正态随机数，恰好对应标准RLS递推公式分母里的那个“1”。如果实际噪声方差不是1，比如你把randn的结果乘以sqrt(2)，那K公式里的1就要改成2，否则增益计算不匹配。很多初学者直接照搬公式，换了噪声却不改分母，结果收敛慢或偏差大，这是第一个容易踩的坑。

为什么选500点？二阶系统有4个参数，理论上每个参数至少需要几个数据点的信息量，但考虑到噪声方差为1、系统极点模长0.9487接近单位圆，早期瞬态需要一段时间衰减。500点覆盖了M序列的一个完整周期，既保证激励充分，又不会让计算量过大。如果换成100点，你会看到参数曲线还没稳定就结束了。

2.3 递推最小二乘的四个核心变量：Theta、P、K、h

标准RLS代码段如下：

P=100*eye(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

四个核心变量说清楚：

Theta是参数估计向量，这里用4×(N+1)矩阵存每一拍的结果。Theta(:,i-1)就是第i个数据点对应的估计值，初值故意设成[1;1;1;1]，离真值[-1.8;0.9;0.75;0.5]很远，这样才能看出收敛过程。

P是误差协方差矩阵，初值100*eye(4)，表示对初始估值非常不确信。P越大，第一拍增益就越大，修正幅度也大。随着递推进行，P会逐渐变小，增益也随之变小。

K是增益向量，由P和h计算得到，决定当前新息z(i)-h'*Theta的权重。新息越大，修正越多。

h是回归向量，由上一拍输出、上两拍输出、上一拍输入、上两拍输入组成，作用是把当前数据映射到参数空间。

递推逻辑只有一句话：用预测误差乘以增益修正参数，再更新P。注意K公式里的inv，在小维度里没问题，但换成除法K=P*h/(h'Ph+1)会更稳。后面的避坑章节我会细说。

Pstore单独存对角线，是为了观察协方差收敛情况。如果Pstore下降很快并趋于平稳，说明数据信息在积累；如果一直震荡，多半是输入激励不足或模型阶次不对。另外，Pstore定义成4×(N+1)，但循环里实际只填到第N列，最后一列保持0，不影响画图。

3. 跑通标准RLS程序：Matlab代码逐行拆解

3.1 生成M序列的移位寄存器实现

M序列生成代码在上一章已经出现，但有几个细节值得单独说。

寄存器初值x不能全零。如果全部是0，异或结果永远是0，移位后序列也全是0，输入就没有任何激励。这里用[0 1 0 1 1 0 1 1 1]，非零初值，且包含足够的0/1混合。

抽头位置选择第4和第9，这是基于本原多项式。抽头不同，序列周期可能会缩短。如果周期缩短，伪随机性变差，频谱会出现较强的离散线谱，激励不够平均。我一般会在仿真前先检查M序列的周期：把M序列重复两次，用find(diff(M)==1)之类的办法找重复，但最简单的是直接看前100点是否有明显规律。如果看到长串连续0或连续1，说明寄存器配置不好。

M序列的幅值这里是0和1，相当于有直流分量。严格来说，持续激励信号要求频谱丰富，带直流也可以，但直流分量会让输出z也带上偏置，在有色噪声或非零均值环境下会影响参数估计。后面渐消记忆章节你会看到，原程序最终结果仍有一点偏差，输入有直流也是原因之一。把0改成-1可以消除直流，这是常见做法。

3.2 生成观测序列z(k)的递推式

观测数据生成段：

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(1)和z(2)是初始状态，分别设为-1和0。由于系统稳定，初始状态的影响会随时间衰减，但前几个数据点仍会受影响。这相当于真实系统在开始采集前已有未知初始状态，辨识算法需要靠后期数据慢慢纠正。如果你希望避免初始瞬态干扰，可以先让系统空跑50拍，再从第51拍开始记录。不过在这个例子里，500点足够长，末尾的收敛结果基本不受影响。

这里有一个容易被忽视的下标映射：z数组从1开始，M数组从1开始。当i=3时，用M(2)和M(1)，对应u(k-1)和u(k-2)。也就是说，第1个可用输出是z(3)，它依赖前两拍输出和M的前两个点。因此z的长度是N+2，多出两个元素给初值。

3.3 递推循环：增益矩阵与协方差更新

把递推循环单独拿出来，加一个残差变量便于调试：

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); e = z(i) - h'*Theta(:,i-2); Theta(:,i-1)=Theta(:,i-2)+K*e; P=(eye(4)-K*h')*P; Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; end

这段代码是RLS的核心。h'Theta(:,i-2)就是用上一拍参数对当前输出做预测，e是预测误差。Ke就是对参数的修正量。P更新是标准公式，没有遗忘因子，所以P会随着数据增加不断减小。

如果你打印Pstore的最后一列，会发现P对角线已经非常小，比如1e-3以下。这意味着参数估计的置信度很高，但也意味着后续新数据的修正能力很弱。如果此时发现最终估计偏了，想再拉回来几乎不可能。这正是第4章要解决的问题。

调试时建议在循环里加一个计数器，隔50拍打印一次Theta的当前值，观察是否单调靠近真值。如果中途出现跳跃，可能是h向量构造错了，比如下标偏移一位。我在做模拟项目X时，就曾把M(i-2)误写成M(i-1)，结果前200拍参数完全乱跳，后来逐行对公式才找到问题。

3.4 画图与收敛性判断

程序最后画图：

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

Theta是4×(N+1)，因为循环从i=3到N+2，列号从1到N。画图时i取1:N+1，正好覆盖所有列。真值线用ones(1,N+1)乘对应常数，画成水平直线。注意Theta(1,:)对应的是-1.8，Theta(2,:)对应0.9，Theta(3,:)对应0.75，Theta(4,:)对应0.5。

判断收敛不能只看最后一拍，还要看曲线是否在末尾平稳。如果末尾还在明显波动，说明数据长度不够或遗忘因子太大。也可以把Pstore画出来：

figure(2); plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:));

Pstore下降速度反映参数信息积累速度。如果Pstore在几拍内突然掉到接近0，说明P初值可能设太大了；如果Pstore几乎不变，说明输入激励不够。

4. 渐消记忆最小二乘：把遗忘因子0.998塞进哪里

4.1 为什么普通RLS会有稳态偏差

原程序注释说“辨识结果会与其精确值始终存在偏差”。这不是程序bug，而是标准RLS在这个场景下的固有表现。

原因有两层。第一，数据长度有限。500点虽然是完整M序列周期，但实际噪声是randn生成的一次性样本，有限样本噪声均值不一定为0，方差也不完全等于1，这些随机波动会残留在参数估计里。第二，标准RLS的P矩阵随递推收敛到很小，增益K趋近于零。当P到后期特别小时，参数修正量K*e几乎为零，参数曲线变成一条水平线。如果前期的数据有较大的噪声干扰，后期就没机会再纠正了。

你可以做一个小实验：把P初值从100改成1，会发现最终偏差变大，因为P太小导致早期修正不足。反过来把P初值改成10000，会发现初阶段参数振荡剧烈，但最终偏差不一定更小。这说明P初值也是一个需要配合数据长度调节的参数。

4.2 两行语句的改动：遗忘因子的作用位置

渐消记忆最小二乘只改两行，非常优雅：

K=P*h*inv(h'*P*h+0.998); P=(eye(4)-K*h')*P/0.998;

对比标准RLS：

K=P*h*inv(h'*P*h+1); P=(eye(4)-K*h')*P;

改动后的P更新多除了一个0.998。因为0.998小于1，所以P每次更新后不是缩小，而是被放大到约1.002倍。P不再趋近于零，增益K也保持在一个合适范围。这样即使到了后期，新数据到来时仍能对参数进行修正，把早期残留的偏差逐步拉回真值。

为什么K公式里的1也改成0.998？这一项对应观测噪声方差，在渐消记忆版本里，遗忘因子λ以类似的形式参与分母。完整的渐消记忆RLS通常写成：

lambda = 0.998; K = P*h / (lambda + h'*P*h); P = (eye(4)-K*h')*P / lambda;

记住lambda=0.998，所以分母里那个数就是lambda，P更新除以lambda。0.998非常接近1，基本保留了全部历史数据，但又没有完全让P衰减到零，属于“既维持原有记忆，又保留修正能力”。如果lambda等于1，就退化成标准RLS。

4.3 遗忘因子选多大：0.95~0.999的权衡

遗忘因子决定了旧数据被淡忘的速度。lambda越接近1，遗忘越慢，稳态估计方差越小，但对参数变化的跟踪越迟钝；lambda越小，跟踪越快，但参数曲线抖动越大。

我做某跨平台系统仿真时，曾把lambda设成0.98测试一个参数缓慢漂移的系统，效果比0.998好很多，但代价是稳态波动明显增大。在真实应用里，选lambda的经验是：先设0.99，看收敛曲线，再根据实际方差微调。如果滤波后曲线平稳但跟踪太慢，减小lambda；如果曲线毛刺太多，增大lambda。

一个常用技巧是让lambda先取0.98跑前100点，之后切回0.998，这样可以兼顾初期的快速收敛和后期的低波动。把这段切换逻辑写进循环里并不复杂：

if i < 100 lambda = 0.98; else lambda = 0.998; end K = P*h / (lambda + h'*P*h); P = (eye(4)-K*h')*P / lambda;

这样做之后，你会发现前100点参数曲线迅速接近真值，后面平稳微调，整体效果比单一lambda更好。这个变体在一些教科书里叫“变遗忘因子RLS”，很适合在课程设计里展示。

5. 参数辨识避坑指南：初值、噪声、激励和收敛判断

5.1 收敛失败：先查输入激励和初值

现象：Theta曲线发散，或长时间在真值附近振荡不收敛。

原因：最常见的是输入信号不满足持续激励。如果输入是常值或正弦单频，h矩阵的行向量之间线性相关，信息矩阵不满秩，参数不可辨识。另外，P初值取得过小也会导致早期修正不足。比如P=eye(4)而不是100*eye(4)，第一拍增益就小，从远初值靠近真值会很慢。

解决：确保输入用M序列或叠加多个不同频率的信号。M序列要检查周期，如果抽头配置不对，序列周期变短，频谱出现线谱，相当于激励不充分。初值方面，P初值建议至少十倍于参数的量级。比如参数真值是1.8量级，P初值取100就合理。我把P初值降为1试过，参数曲线需要近200拍才靠近真值，明显变慢。

5.2 数值爆炸：inv、NaN和P奇异

现象：递推几个循环后，Theta突然变成NaN或Inf，或者P矩阵对角线急剧增加。

原因：维度不高时inv一般没问题，但如果P初值极大、输入又极小，K可能会非常大，导致P更新出现数值振荡。另外，如果数据重复导致h'Ph+1接近0，inv的值爆炸。还有，如果系统本身不稳定，输出z发散，RLS也会跟着发散。

解决：把inv改成右除，写成K=P*h/(h'Ph+1)，数值稳定性更好。在分母里加一个很小的正数作为正则项，比如h'Ph+1e-6，能防止矩阵奇异。还要检查观测序列是不是真的稳定：画一下z的曲线，如果振幅随时间增长，那就不是辨识算法的问题，而是系统设置或数据生成代码下标错了。我见过有人把z(i-1)误写成z(i+1)，输出立刻发散。

5.3 稳态偏差：直流、有色噪声和渐消记忆

现象：参数曲线很平滑，但不是停在真值上，而是停在真值附近的一个固定偏差点。

原因：标准RLS后期增益太小，早期偏差无法修正；或者M序列是0/1，带直流分量，使得输出包含偏置，估计参数会吸收这个偏置；或者噪声不是纯白噪声，而是有色噪声，最小二乘估计会有偏。

解决：对于固定偏差，用渐消记忆RLS，如第4章。对于直流分量，把M序列的0改成-1，1仍为1，这样输入均值为0。对于有色噪声，可以考虑广义最小二乘或增广最小二乘，但那是另一个主题。有一个快速验证方法：把randn噪声改成零序列v=zeros(1,N)，再跑一遍标准RLS，如果偏差消失，说明偏差主要来自噪声样本的有限均值；如果偏差还在，说明是输入直流或模型结构问题。

5.4 数据长度不足：500点到底够不够

现象：把N从500改成100后，参数曲线不收敛，末尾仍明显波动。

原因：100点不足以覆盖M序列完整周期（511），激励频率不完整，且瞬态影响占比更大。另外，系统极点模长0.9487，时间常数约为1/(1-0.9487)≈19.5拍，100点只包含约5个时间常数，勉强够，但如果有噪声，统计精度就不足。

解决：N至少取一个M序列周期，最好取2个周期。本例500接近511，已经足够。如果你改周期或系统，可以用经验公式：N至少是待估参数个数乘10到20倍，对于4个参数，最少40到80点，但那是理想情况，有噪声时建议翻倍再翻倍。

5.5 收敛判断：别只看最后一点

现象：参数曲线在最后200拍已经平了，但最后一拍值仍与真值有偏差，有人误以为算法失败。

原因：RLS的收敛是一个统计过程，最终估计是在噪声影响下的一个概率结果，不会严格等于真值。曲线平了只代表参数不再大幅更新，不代表就是真值。

解决：把最后50拍的Theta取平均作为最终估计值，能有效降低单点随机波动。也可以看残差序列的均方根值，如果接近噪声标准差1，说明辨识质量可以接受。我一般会这样写：

theta_final = mean(Theta(:,end-49:end),2);

这会稍微改善估计精度，尤其是当稳态偏差不大时。

6. 验证与进阶：用残差分析和改变输入信号确认辨识结果

6.1 残差白噪声检验

参数辨识出来之后，不能只看曲线是否贴近真值。一个更严格的做法是检验残差序列是否接近白噪声。残差是最终模型对每个观测输出的预测误差。写一段代码：

theta_hat = mean(Theta(:,end-49:end),2); % 用最后50拍均值 res = zeros(1,N); for i=3:N+2 h=[-z(i-1); -z(i-2); M(i-1); M(i-2)]; res(i-2)=z(i)-h'*theta_hat; end figure; subplot(2,1,1); plot(res); title('残差序列'); subplot(2,1,2); autocorr(res); grid on;

如果残差自相关函数在非零延迟处基本落在95%置信区间内，说明模型已经捕获了系统的动态，剩下的只是白噪声。如果自相关明显超界，比如延迟为1或2时超出，说明模型阶次偏低，或者噪声本身有色。对于本例，残差应该接近白噪声。如果残差均值明显不为零，多半是输入直流或模型偏差。

6.2 不同遗忘因子对比实验

把渐消记忆RLS封装成一个循环，分别取lambda=1、0.998、0.97、0.95，记录每种情况最后的估计误差。你几乎肯定能看到：lambda=1时误差最大，lambda=0.998时误差减小，lambda=0.95时误差可能略微增大但跟踪能力最强。这个对比实验是理解遗忘因子最好的方式。我建议把四条Theta(1,:)曲线画在同一张图里，你能清楚看到哪个收敛最快、哪个波动最大。实际跑下来，lambda=0.95的曲线在末尾会像毛刺一样跳动，而lambda=0.999的曲线几乎是一条平线，但位置偏高。

6.3 扩展到更高阶系统与变体

这套程序框架非常容易扩展。把h向量从4维改成2n维，Theta初值也改成相同长度，就能辨识一个n阶系统。递推公式本身不关心阶次，只关心h与Theta的维度匹配。我曾经把一个二阶模型扩展到四阶，只改了z生成代码里的系数、h的构造和真值线常量，其他都没动。

扩展时最容易翻车的是下标偏移。高阶差分方程往往涉及z(i-1)到z(i-n)，输入m(i-1)到m(i-n)，写循环时务必画一个时间轴，验证i的起点和z的初值个数。我的习惯是先用零噪声跑一遍，参数能精确回到真值，再上噪声。这样能把代码bug和噪声影响分开。

还有一个进阶技巧：在标准RLS基础上加入输出权值限制，可以防止参数突变。具体做法是每次更新后检查Theta的每个分量是否超出合理范围，超了就按上界方向拉回。这个技巧在工业现场特别管用，因为传感器偶尔会有野值，一旦野值进入递推，参数可能瞬间跑飞。从那以后我每次做辨识，都强制走一遍“零噪声验证、残差白噪声检验、遗忘因子对比”这三步，再也不敢只看估计曲线就收工。希望帮到你。

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