简介：面向控制系统、信号处理与机器学习方向的学习者与工程人员，这份资源提供最小二乘类参数辨识的MATLAB实现，重点解决递推最小二乘与渐消记忆法在系统参数在线估计中的编程落地问题。压缩包共1个docx文件，约78KB，以文档形式完整收录辨识思路、系统设计说明与MATLAB程序代码，便于对照阅读与复现。资源已有2162人学习，说明其在参数辨识入门与课程实践中具有一定参考价值。内容围绕一个极点位于单位圆内的稳定系统展开，先以500点M序列作为输入、叠加白噪声扰动生成观测数据，再用递推最小二乘完成参数估计，并通过估计方差P、增益矩阵K的更新过程展示算法细节；随后将递推语句改为渐消记忆形式，对比辨识精度变化，并给出参数估计过程的可视化结果。读者可据此理解两类算法的差异、掌握程序结构，并迁移到自身辨识任务中。
1. 递推最小二乘：为什么参数估到最后总是差一口气
跑过系统辨识的人大概都见过这样一张图：四条参数曲线从初值出发，一路收敛，最后稳稳停在真值附近，但就是不肯贴上去，始终留着一道肉眼可见的缝。你以为是噪声太大、数据太少，于是加长采样、换输入信号，缝还在。这个现象不是玄学，而是递推最小二乘（RLS）在有色扰动下的固有偏差——它优化的目标和你想辨识的真参数之间，本来就隔着一层。
这份 MATLAB 程序把这道缝完整地暴露了出来：先用标准 RLS 辨识一个二阶稳定系统，四个参数收敛但带偏；再把增益更新和协方差更新里的标量 1 换成 0.998，偏差立刻被压下去。两段代码只差两个数字，结果却完全不同。搞懂这两个数字为什么能起作用，就搞懂了遗忘因子、协方差矩阵和新息这三者之间的耦合关系，也就明白了参数辨识里最容易被忽略的那层机制。下面从系统构造一路拆到源码级实现。
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)把它写成标准形式，待辨识的参数向量是[-a1, -a2, b1, b2]，真值分别为-1.8、0.9、0.75、0.5。特征方程z^2 - 1.8z + 0.9 = 0的根模长为sqrt(0.9) ≈ 0.949，落在单位圆内，系统稳定。这一点很关键：递推算法每一拍都要用历史输出z(i-1)、z(i-2)构造回归向量，如果系统本身发散，回归量会指数增长，协方差矩阵P的数值条件迅速恶化，辨识根本无从谈起。稳定性是这套方法能跑通的前提，不是可有可无的装饰。
常见的实现思路是，辨识前先对系统做一次开环仿真，确认输出有界、量级合理，再喂给估计算法。程序里z(1)=-1, z(2)=0给了两个非零初值，正是为了让回归向量从第一拍起就非退化。
2.2 M 序列：为什么要用伪随机二进制激励
输入信号由 9 级移位寄存器加异或反馈生成：
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; endM 序列（最大长度伪随机二进制序列）的周期是2^9 - 1 = 511，程序取前 500 点，基本覆盖一个完整周期。选它做激励有两个实打实的理由：一是它的自相关函数接近冲激，频谱足够宽，能持续激励系统的所有模态——如果输入是常数或缓慢变化的信号，回归向量各列近似共线，h'*P*h接近奇异，参数就不可辨识；二是它幅值只有 0 和 1，实现简单，便于在嵌入式或半实物平台上复现。
提示：寄存器初值不能全零，否则异或反馈恒为 0，序列退化为全零，激励失效。这是手写 M 序列最常见的翻车点。
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注意索引对齐：z(i)对应方程里的Z(k+2)，所以输入取M(i-1)、M(i-2)，噪声取v(i-2)。这种错位在移植代码时极易写错，一旦错一拍，回归向量和输出就对不上，辨识结果会系统性偏移，而且从曲线上看不出明显异常，属于典型的黑匣子式故障。
3. 递推最小二乘的源码解剖：一次更新到底算了什么
3.1 回归向量与四个状态量的角色
核心循环只有六行，但每一行都对应最小二乘的一个数学环节：
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)]; endh是回归向量，前两项是取负的历史输出，后两项是历史输入，正好对应差分方程里-a1*z(k+1)-a2*z(k)+b1*u(k+1)+b2*u(k)的结构。Theta是待估参数，P是估计误差的协方差矩阵，K是增益。三者构成一个闭环：P决定这一步该信多少新数据，K把这份信任折算成参数修正量，更新后的P又反映修正之后还剩多少不确定性。
P初值取100*eye(4)而不是零矩阵，含义是“初始参数完全不可信，允许大幅修正”。如果取零，增益K恒为零，参数永远停在初值上，这是新手最容易踩的坑之一。
3.2 增益公式里的那个 1 是什么
K=P*h*inv(h'*P*h+1)中分母的1不是随手写的常数，它对应新息e(k)=z(i)-h'*Theta的方差。标准 RLS 假设扰动是白噪声、方差为 1，所以分母是h'*P*h + σ²，这里σ²=1。如果实际噪声方差是 4，分母却写 1，增益就会被系统性高估，参数在真值附近来回抖，收敛慢且稳态方差大。
inv在这里其实是对标量求逆，写成1/(h'*P*h+1)更高效，MATLAB 的inv对标量也能算，但会多一层函数调用开销。500 点数据下无所谓，换成几万点的在线辨识，这个细节就值得改。
3.3 新息驱动：参数为什么是“一点点挪过去”的
参数更新项K*(z(i)-h'*Theta(:,i-2))里，括号内是预测输出与实测输出之差，即新息。新息大，说明当前模型预测不准，参数修正幅度大；新息小，说明模型已经贴合，修正趋缓。K相当于一个随协方差自适应调整的学习率：初期P大，K大，参数快速向真值靠拢；后期P收缩，K变小，参数趋于稳定。
这种“增益随不确定性衰减”的结构，正是 RLS 相比批处理最小二乘的优势——不需要存全部历史数据，每来一个新样本做一次 4×4 矩阵运算即可，天然适合在线估计。代价是它对初值和噪声统计特性敏感，偏差一旦形成，靠后续数据很难自行纠正。
3.4 协方差更新与 P 的单调收缩
P=(eye(4)-K*h')*P是协方差更新。可以验证，更新后的P是原P在h的正交补空间上的投影再缩放，对角元单调不增。程序把P的四个对角元存进Pstore并画出来，就是为了观察这个收缩过程：如果P对角元很快掉到接近零，说明算法已经“自信”，后续新数据几乎不再改变参数；如果掉不下去，说明激励不足或噪声过大。
这个单调收缩特性也解释了标准 RLS 的偏差为什么改不掉：P一旦收缩，增益K随之变小，算法对持续存在的模型误差越来越“迟钝”，偏差被锁死在收敛值上。要纠正它，必须人为阻止P无限收缩，这正是下一章要动的地方。
4. 渐消记忆：把 1 换成 0.998 之后发生了什么
4.1 遗忘因子的物理含义
把两行改成：
K=P*h*inv(h'*P*h+0.998); % 分母常数从 1 改为 0.998 P=(eye(4)-K*h')*P/0.998; % 协方差更新后整体除以 0.9980.998 就是遗忘因子 λ。它的作用是在每一拍把旧信息按比例“打折”：P除以 λ 相当于把协方差整体放大1/λ，让算法重新变得“不自信”，增益K不会一路衰减到零，新数据始终保有修正能力。λ 越接近 1，记忆越长，越接近标准 RLS；λ 越小，越偏向近期数据，跟踪时变参数的能力越强，但稳态方差也越大。
分母里的 0.998 与P更新里的/0.998必须配对出现，前者调整增益尺度，后者补偿协方差收缩。只改一处，增益和协方差的量纲就对不上，参数会发散或剧烈振荡。这是改这段代码时最容易翻车的地方。
4.2 为什么偏差被压下去了
标准 RLS 的偏差来源是扰动v与回归向量中的历史输出相关——z(i-1)、z(i-2)里含有过去的噪声分量，而回归向量又用它们做输入，形成噪声自相关，最小二乘的正交性条件被破坏，估计不再无偏。遗忘因子并不直接消除这种相关性，但它让算法持续给近期数据更高权重，等效于在一个滑动窗口上做加权最小二乘，窗口内的噪声相关性被稀释，偏差随之减小。
从曲线看，四个参数更快贴近-1.8、0.9、0.75、0.5四条参考线，收敛段更陡，稳态缝隙明显收窄。这不是玄学调参，而是 λ 改变了有效记忆长度，把“无限记忆累积偏差”换成了“有限记忆局部拟合”。
4.3 λ 的取值与跟踪能力的权衡
λ 的选取没有普适公式，常见做法是按有效记忆长度1/(1-λ)估算：λ=0.998 对应约 500 拍记忆，恰好与 500 点数据量匹配；λ=0.99 对应约 100 拍，适合参数缓慢时变；λ=0.95 对应约 20 拍，跟踪快但噪声敏感。程序选 0.998 不是巧合，而是让记忆窗口覆盖整个采样段，既保留足够数据平均噪声，又避免无限记忆锁死偏差。
如果系统参数在运行中跳变，λ 应取小一些，让算法快速跟上；如果参数恒定、只求精度，λ 取 0.999 以上更稳。这个取舍没有标准答案，取决于你对“跟踪速度”和“稳态精度”哪个更在意。
5. 避坑与排查：这段代码跑不对时先看哪里
5.1 参数曲线纹丝不动
现象：四条曲线从初值出发几乎是一条直线，Pstore对角元也不变化。原因通常是P初值给了零矩阵，或者h构造时用了全零数据。解决：确认P=100*eye(4)，并检查z和M的前几个点是否非零、索引是否错位。
5.2 参数发散或数值溢出
现象：曲线指数增长，P对角元迅速变大。原因多为遗忘因子只改了一处，或者 λ 取得过小（如 0.9 以下）导致协方差被反复放大。解决：确保K分母和P更新里的 λ 成对出现，λ 不低于 0.95；必要时在P更新后加对称化P=(P+P')/2抑制数值漂移。
5.3 辨识结果整体偏移一个固定量
现象：曲线收敛但四条都偏同一方向。原因通常是回归向量索引错一拍，或差分方程里输入项符号写反。解决：对照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)逐项核对h的四个分量，特别注意M的索引比z滞后一拍。
5.4 收敛极慢、500 点不够用
现象：曲线到末尾仍在缓慢爬升。原因可能是 M 序列周期与采样长度不匹配，激励不充分，或噪声方差远大于 1 而分母仍写 1。解决：确认 M 序列长度覆盖一个周期以上，按实际噪声方差调整增益分母常数，必要时延长采样。
5.5 换一组随机噪声结果差异巨大
现象：每次运行曲线形状差别明显。这是单次噪声实现的正常波动，不是代码错误。解决：固定随机种子rng(0)做对比实验，或多次蒙特卡洛取平均，再判断算法本身的偏差水平。
6. 从两个数字看参数辨识的设计取舍
这段程序最值得琢磨的地方，是它用两个数字（分母的 1 和 0.998）把“无偏”和“稳定”这对矛盾摆在了台面上。标准 RLS 在数学上是最优的——在扰动为白噪声、模型结构精确的前提下，它给出最小方差无偏估计；但现实里扰动往往与回归量相关，模型阶次也未必精确，此时“最优”反而变成“有偏且不可纠正”。遗忘因子牺牲了一部分统计最优性，换来对模型误差的鲁棒性，本质是用偏差换跟踪能力。
这个取舍在工程上随处可见。批处理最小二乘可以一次算完、精度高，但不适合在线；RLS 省内存、能实时更新，代价是对初值和噪声统计敏感；渐消记忆进一步放弃长期记忆，换取对时变参数的适应。没有哪种方法全面占优，选型取决于系统是恒定还是时变、噪声是白还是有色的、算力是否允许存全部历史。
验证这套逻辑是否真的成立，最直接的办法是做对照实验：固定同一组 M 序列和噪声，分别跑 λ=1、0.998、0.99 三组，把Theta四条曲线叠在一张图上，同时记录Pstore对角元的衰减速度。你会看到 λ 越小，收敛越快、稳态抖动越大，偏差和方差此消彼长，曲线不会骗人。再换几组不同种子重复，观察偏差的均值和方差，就能判断某个 λ 是不是真的适合当前场景，而不是凭感觉调参。
我自己的习惯是，每次改这类递推算法的参数，都先把P的对角元和增益K的范数打出来看一遍，确认它们的变化趋势符合预期，再去看参数曲线——曲线好看不代表机制对，中间状态量才是黑匣子里的真相。从那以后我每次动遗忘因子，都强制走一遍“固定种子、对照三组、看 P 衰减”的流程，省下不少返工。希望帮到你。
本文还有配套的精品资源，点击获取
