简介：面向控制系统、信号处理与机器学习方向的学习者与工程人员，这份资源提供一套最小二乘类参数辨识的MATLAB程序，用于对系统参数进行在线估计与辨识。程序以递推最小二乘算法为主线，构造极点位于单位圆内的稳定系统，生成500点M序列输入并叠加白噪声扰动，定义参数估计值Theta、估计方差P与增益矩阵K完成递推更新，并借助MATLAB绘图功能可视化参数收敛过程；随后改用渐消记忆最小二乘算法，通过调整增益与方差更新语句获得更准确的辨识结果。资源包共1个docx文件，约78KB，以文档形式整理程序代码、仿真结果与算法说明，便于对照理解递推最小二乘与渐消记忆法的差异。已有2161人学习，适合希望掌握参数辨识基本流程、复现仿真实验并理解算法实现细节的读者参考。

1. 从一组“收敛了但始终差一点”的曲线说起

做系统辨识的人大多遇到过这种场面：递推最小二乘跑完五百个点，四条参数曲线稳稳地收敛了，可把估计值和真值摆在一起看，每条都差那么一截，不随数据量增加而消失。这不是代码写错了，而是标准递推最小二乘在有色噪声或模型阶次不匹配时的固有偏差。这份最小二乘类的参数辨识 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)

待估参数共四个：a1 = -1.8、a2 = 0.9、b1 = 0.75、b2 = 0.5。写成回归形式就是 z(k) = h'·θ + v(k)，其中回归向量 h = [-z(k-1); -z(k-2); u(k-1); u(k-2)]，参数向量 θ = [a1; a2; b1; b2]。这里有个容易看花眼的地方：程序里 Theta 的第一分量对应的是 -1.8 而不是 1.8，因为回归向量里 z 项带了负号，画图时参考线也画在 -1.8 和 0.9 上，对照时别把符号搞反。

特征多项式 z² - 1.8z + 0.9 的极点模长为 √0.9 ≈ 0.949，落在单位圆内，系统稳定。这一点很关键：只有稳定系统，输出才不会发散，递推过程才有意义。如果换成极点跑到圆外的对象，输出指数增长，协方差矩阵 P 很快数值爆炸，辨识曲线会直接飞出去，那不是算法问题，是对象选错了。

2.2 输入信号为什么选 M 序列

辨识的激励信号必须“足够丰富”，能把系统的动态模态都激发出来。程序用 9 级移位寄存器生成 500 点 M 序列（伪随机二进制序列），初值x=[0 1 0 1 1 0 1 1 1]，反馈抽头取第 4 位和第 9 位异或：

x = [0 1 0 1 1 0 1 1 1]; % 9 级寄存器初值，不能全零 n = 503; % 多跑 3 拍，保证取到 500 个有效点 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

M 序列的周期是 2⁹ - 1 = 511，取 500 点基本覆盖一个完整周期，自相关接近冲激，频谱平坦，对二阶对象来说激励充分。这里有个血泪经验：寄存器初值不能全零，全零状态下异或永远输出零，序列直接死掉；另外抽头位置决定了序列的周期和统计特性，换抽头等于换了一套激励，辨识结果的可比性就没了，做对照实验时别随手改。

2.3 噪声与观测序列的生成

白噪声用v = randn(1,N)生成，方差为 1。观测序列递推生成时要注意索引对齐：z 预分配 N+2 个元素，前两个给初值z(1)=-1, z(2)=0，从 i=3 开始按差分方程算，输入取 M(i-1) 和 M(i-2)，噪声取 v(i-2)。这种“多留两拍”的写法是为了让回归向量 h 在 i=3 时就能凑齐 z(i-1)、z(i-2)、u(i-1)、u(i-2) 四项，避免开头几拍索引越界。很多同学在这里翻车，把 z 的长度设成 N，循环里一取 z(i-1) 就报错，或者错位一拍导致辨识出的 b1、b2 互换。

3. 递推最小二乘的实现：从公式到可跑代码

3.1 三个核心矩阵的初始化

递推最小二乘每拍做三件事：算增益、更新参数、更新协方差。初始化直接决定前几十拍的收敛形态：

P = 100*eye(4); % 协方差初值，大数表示“对初值极不信任” Pstore = zeros(4, N+1); % 记录 P 对角线，观察收敛 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 取 100·I 是常规做法，数值越大表示初始不确定性越大，算法前期步子迈得大、收敛快；取太小则前期几乎不动，看起来像“算法没生效”。Theta 初值故意设成 [1;1;1;1]，和真值 [-1.8;0.9;0.75;0.5] 差得很远，正好用来观察收敛过程。如果初值直接给真值，曲线一条直线，什么信息都看不出来，答辩时也讲不出东西。

3.2 递推主循环与增益更新

主循环从 i=3 跑到 N+2，每拍构造回归向量、算增益、更新参数和协方差：

for i = 3:N+2 h = [-z(i-1); -z(i-2); M(i-1); M(i-2)]; % 回归向量，z 项带负号 K = P*h*inv(h'*P*h + 1); % 增益，分母的 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

分母里的+1对应噪声方差 σ²=1，如果噪声方差改成 0.1，这里要同步改成 0.1，否则增益尺度不对，收敛速度会明显偏慢。inv(h'*P*h+1)在四维情况下是标量求逆，开销可以忽略；但如果把阶次提到十阶以上，建议改成左除K = P*h/(h'*P*h+1)，数值更稳。参数更新项z(i) - h'*Theta(:,i-2)是新息，也就是当前观测与预测之差，新息越大说明模型越不准，增益就把它按比例分配到四个参数上。

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('参数估计过程');

跑完会看到：曲线在前 50 拍快速收敛，之后基本走平，但四条都停在离参考线一小段距离的位置，而且这段距离不随 N 增大而缩小。这就是标准递推最小二乘的偏差问题——当噪声与回归向量相关（有色噪声、或模型阶次偏低）时，最小二乘估计不再是无偏的。读图时重点看两件事：收敛速度（前段斜率）和稳态偏差（后段与参考线的间距），前者反映 P 初值和增益设计，后者反映算法本身的统计性质。

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;

第一处把分母的 1 换成 0.998，第二处在协方差更新后除以 0.998。0.998 就是遗忘因子 λ，取值在 0 到 1 之间，越接近 1 记忆越长、越接近标准递推；越小则越“健忘”，对近期数据加权更重。这里 λ=0.998 对应的有效记忆长度约为 1/(1-λ) = 500 拍，和 500 点数据量刚好匹配，所以稳态偏差被明显压下去。

4.2 遗忘因子为什么能消偏差

标准递推最小二乘等价于对全部历史数据做等权最小二乘，旧数据和新数据权重一样。当噪声有色时，旧数据里的噪声分量会持续“拖住”参数估计，形成系统性偏差。渐消记忆相当于给每拍数据乘上一个指数衰减权重，越老的数据权重越低，参数估计跟着系统当前的统计特性走，偏差自然被削弱。代价是估计方差会略微变大，曲线尾部抖动比标准版明显一点，这是精度和稳定性之间的正常交换，不是程序出错。

4.3 参数怎么调、边界在哪

λ 的选取没有万能值，通常的做法是按数据长度和系统时变程度折中：数据 500 点、系统定常，λ 取 0.995 到 0.999 都合理；系统参数缓慢时变，λ 要往小取，比如 0.98，让算法跟得上变化；λ 取到 0.99 以下时，有效记忆只剩一百拍左右，曲线抖动会明显加剧，答辩时容易被问“你这估计方差是不是太大了”。另外 λ 不能取 1，取 1 就退化成标准递推，改了等于没改；也不能取太小，太小相当于只用最近几拍数据，协方差矩阵 P 会周期性膨胀，数值上容易出问题。我一般会先固定 λ=0.998 跑一遍看稳态偏差，再上下各调一档做对比，把三条曲线放进同一张图里，偏差和抖动的权衡关系就直观了。

5. 避坑与排查：这类辨识程序最容易翻车的五处

5.1 曲线收敛但四条参数整体错位

现象是曲线形状正常、收敛也快，但估出来的 b1、b2 和真值对不上，甚至互换。原因多半是观测序列生成时输入索引错了一拍，M(i-1)和M(i-2)写反，或者 z 的循环起点不是 i=3。解决方法是拿前五拍数据手工代入差分方程验算一遍，确认 z(3) 确实由 z(2)、z(1)、M(2)、M(1) 算出，索引对齐后再跑。

5.2 前几拍参数剧烈跳变甚至出现 NaN

现象是曲线开头几十拍上下乱窜，偶尔冒出 Inf 或 NaN。原因通常是 P 初值给得过大（比如 1e6）叠加回归向量元素量级差异大，h'*P*h数值条件变差，求逆失稳。解决办法是把 P 初值降到 100 量级，或者把inv换成左除，必要时对输入输出做归一化，让回归向量各分量量级接近。

5.3 改了遗忘因子结果反而更差

现象是把 λ 从 0.998 调到 0.95 后，稳态偏差没降、抖动却大幅增加。原因是 λ 太小导致有效记忆长度远小于系统收敛所需的数据量，参数还没稳住旧信息就被丢光了。解决方法是按 1/(1-λ) 估算有效记忆长度，保证它至少是系统时间常数的若干倍，先试 0.995 再逐步下调，每次只动一档并保留对照图。

5.4 换一个被辨识对象就发散

现象是原对象跑得好好的，换成三阶或极点靠近单位圆的系统后输出发散、P 矩阵数值爆炸。原因是对象本身不稳定或接近临界稳定，开环输出指数增长，任何辨识算法都救不回来。解决办法是先检查特征多项式根的模长，确认全部落在单位圆内且留有余量，必要时先做闭环采集数据再辨识。

5.5 噪声方差改了但增益没跟着改

现象是把randn乘上 0.1 把噪声方差改成 0.01 后，收敛速度肉眼可见地变慢。原因是增益公式分母里的+1仍按方差为 1 写，尺度和实际噪声不匹配。解决办法是把分母常数同步改成噪声方差，或者干脆把它写成一个变量sigma2，在程序开头统一定义，改噪声时只改一处。

6. 验收与答辩：用三张图把方案讲圆

验收这类辨识程序，不看代码行数，看三样东西能不能自圆其说。第一张图是标准递推最小二乘的收敛曲线，要能指出前段收敛速度和后段稳态偏差的具体数值，比如“第 50 拍基本收敛，稳态偏差约 0.05”；第二张图是渐消记忆版的对照曲线，把两版画在同一坐标系里，偏差缩小了多少、尾部抖动增加了多少，用数字说话；第三张图是协方差对角线 Pstore 的衰减曲线，它证明算法确实在收敛而不是靠初值撑着。三张图摆出来，从“有问题”到“怎么改”再到“改的代价”，逻辑链就完整了。

答辩时最容易被追问的是参数选择依据，提前把下面这张表准备好，被问到直接报数：

还有一个常被忽略的细节：答辩演示时别只放最终收敛图，把前 20 拍的局部放大图也备一张。评审最爱问“你怎么证明它是在辨识而不是在拟合初值”，局部放大图里参数从初值一路爬向真值的过程，比任何解释都有说服力。我早年第一次讲这个项目，只放了稳态图，被追问“曲线平了是不是因为算法停了”，当场卡壳；从那以后我每次做辨识类演示，都强制把前段收敛过程、稳态偏差数值、协方差衰减三样一起放，再留一张参数对照表兜底。希望这套结构能帮到你，把手里那份最小二乘辨识程序讲成一份站得住的完整方案。

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