MATLAB最小二乘参数辨识：批量与递推算法实现及避坑指南

简介：这是一份面向控制、信号处理等领域学习者与科研人员的MATLAB参数辨识程序文档，聚焦最小二乘类算法在系统参数在线估计与辨识中的应用。资源以单个docx文档形式提供，体积仅78KB，内含完整可运行的递推最小二乘与渐消记忆最小二乘代码实现，并配有参数估计过程可视化与结果对比说明。文档基于一个极点位于单位圆内的稳定系统，生成500点M序列作为输入并叠加白噪声扰动，逐步展示增益矩阵、估计方差及参数更新公式的推导与编程细节，同时对比了两种算法在相同条件下的辨识精度差异。该资源目前已有2158人学习下载，适合正在学习系统辨识、自适应控制或MATLAB数值仿真的读者，可将其作为理解最小二乘参数辨识原理、快速上手编程实践的直接参考资料。

1. 最小二乘类的参数辨识MATLAB程序：从“凑模型”到“算参数”的第一步

做过系统辨识的人都有过这种经历：手里握着几十秒的输入输出数据，面对一个“大概是二阶、带纯延迟”的工业过程，想把传递函数里的四个系数定下来。有人用频域法，画Bode图再手动拟合，费眼费神，遇到噪声大一点的采样序列，曲线毛刺能让人怀疑人生。我自己的做法是直接上最小二乘类的参数辨识MATLAB程序——把问题变成“构造回归矩阵，解一次超定方程”，几十行代码就能拿到参数估计值。这套方法不挑对象，线性模型基本都能用，既能离线批量算，也能改成递推形式做在线辨识。适合正在做课程设计、写论文需要辨识模型的学生，也适合在现场调试时想快速拿到一个可用的线性模型、后续再做控制的工程师。读完这篇，你能跑通批量最小二乘和递推最小二乘两套程序，知道噪声环境下怎么补救，也清楚那些“同一份数据，不同人辨识结果不一样”的坑到底出在哪。

2. 为什么参数辨识首选最小二乘：模型假设、解析解与适用边界

2.1 先把要辨识的模型写出来：ARX结构与回归形式

参数辨识的第一步不是写代码，而是把“要辨识的对象”用数学结构固定下来。工业场景里最常见的做法是把对象写成ARX模型，也就是自回归各态外部输入模型，离散形式如下：

y(k) + a1·y(k-1) + a2·y(k-2) + … + ana·y(k-na) = b1·u(k-1) + b2·u(k-2) + … + bnb·u(k-nb) + e(k)

这里y(k)是当前时刻输出，u是输入，e(k)是噪声项，na和nb分别是输出和输入的阶次。这个模型为什么好用？因为它把“过去时刻的输出”和“过去时刻的输入”当成已知量，把未知参数a和b当成待求系数，于是整个问题就变成了一个线性回归问题。把上式改写成矩阵形式会更直观：定义回归向量φ(k) = [-y(k-1), -y(k-2), …, u(k-1), u(k-2), …]ᵀ，参数向量θ = [a1, a2, …, b1, b2, …]ᵀ，那么每个时刻的观测都能写成一个线性方程y(k) = φ(k)ᵀ·θ + e(k)。

这个形式是后面所有代码的地基。无论是批量最小二乘、递推最小二乘还是辅助变量法，本质上都是在处理“从一堆这样的线性方程里把θ解出来”这件事。所以动手写程序前，先想清楚你的对象能不能用ARX结构描述。大多数线性定常系统、带稳定工作点的工业过程，这层近似是合理的。如果你的对象是强非线性、时变或闭环运行的系统，ARX只是起点，后面需要更复杂的结构，但最小二乘类方法仍然是入口。

2.2 解析解与数值实现：为什么不要用inv()

对于批量处理的数据，假设我们有N组观测，把所有方程堆起来就得到矩阵形式：Y = Φ·θ + E。其中Φ的第k行就是φ(k)ᵀ，是一个N×(na+nb)的矩阵，Y是N×1的输出向量。最小二乘的准则是最小化残差平方和J(θ) = ‖Y - Φ·θ‖²。对θ求导并令梯度为零，能得到著名的正规方程：

θ_hat = (Φᵀ·Φ)⁻¹·Φᵀ·Y

理论形式很简洁，但落到MATLAB里我有一个很明确的习惯：永远不要直接用inv(Phi'*Phi)去算。原因有两个，一是Φᵀ·Φ的条件数是原矩阵条件数的平方，数据稍微有些共线，这个矩阵就可能病态到求逆结果完全失真；二是MATLAB的反斜杠运算符有更聪明的求解路径。对超定方程，直接写theta_hat = Phi \ Y，MATLAB会自动选择QR分解路径来求解最小二乘问题，数值稳定性比凑出正规方程再求逆好得多。

有一个经验值供参考：当cond(Phi'*Phi)超过1e6时，inv的解法大概率已经开始丢失有效数字。尤其现场采集的数据，传感器噪声、量化误差、相邻采样点的相关性都会抬高条件数。你可以在程序里用cond()函数确认一下，但顺手写成Phi\Y就能绕开这个雷区。这一点不夸张地说，能救回不少“明明公式对但结果不对”的翻车现场。

2.3 适用边界：白噪声、充分激励与数据长度

最小二乘不是万能药，它的理论保证建立在几个前提上。第一，噪声项e(k)必须是白噪声，也就是与回归向量φ(k)不相关。这个条件在实际中很难严格满足——比如输出测量噪声经过传感器滤波后天然带相关性，又比如闭环运行的系统里，噪声会通过反馈回路污染输入信号。一旦噪声有色，批量最小二乘的参数估计就是有偏的，偏差大小跟噪声强度和相关性有关。第二，输入信号必须持续激励。阶跃信号虽然工程上最容易做，但它的频谱能量集中在低频，对二阶以上的模型往往激不出足够的动态信息，导致参数不可辨识或方差极大。第三，数据长度要够。理论上N大于参数个数就能算，但工程上我一般要求N至少是参数个数的20倍以上，否则估计方差会大得没法看。

这三个边界条件决定了你选什么工具：数据离线、噪声白、激励充分，用普通批量最小二乘即可；数据在线到达、参数慢变，用递推最小二乘加遗忘因子；噪声明显有色，就要上辅助变量法或增广最小二乘。下面两张表可以帮你快速选型。

3. 批量最小二乘MATLAB程序：从仿真数据到参数估计的一条龙实现

3.1 先造一份能验证结果的仿真数据

学辨识最怕的就是拿真机数据练手，因为真实参数未知，算出来的结果没法判断对错。所以我每次教学和自测，都是先造一份“答案已知”的仿真数据，跑通流程后再上真数据。这个习惯能帮你快速定位是程序写错了还是数据本身的问题。下面这段代码生成一个二阶ARX系统的输入输出序列，真实参数设为a1=-0.8、a2=0.15、b1=0.8、b2=0.2，噪声取方差0.01的高斯白噪声。

% 生成仿真数据：y(k)+a1*y(k-1)+a2*y(k-2)=b1*u(k-1)+b2*u(k-2)+e(k) rng(42); % 固定随机种子，保证可复现 N = 1000; % 数据长度 u = randn(N, 1); % 白噪声激励信号 e = 0.01 * randn(N, 1); % 过程噪声，方差为 0.01^2 a1 = -0.8; a2 = 0.15; % 真实参数：输出侧 b1 = 0.8; b2 = 0.2; % 真实参数：输入侧 y = zeros(N, 1); for k = 3:N y(k) = -a1*y(k-1) - a2*y(k-2) + b1*u(k-1) + b2*u(k-2) + e(k); end

这段代码里有一个细节值得说明：差分方程用的系数是“模型里带正号”的a1、a2，而递推计算时负号被显式写了出来。这样做的意图是让真实参数和后面辨识结果直接对照，不用来回变符号。随机种子rng(42)的作用是让每次运行生成相同的数据序列，方便复现和调试。如果你希望更接近工程环境，可以把输入u从randn换成PRBS（伪随机二进制序列），它的频谱更平直、幅值可控，在后文避坑部分我还会专门提。

3.2 批量辨识主程序：构造回归矩阵、求解与残差分析

数据有了，接下来就是核心的辨识程序。这一步最关键的活是把数据矩阵Φ构造对——行对应时刻，列对应回归变量，顺序必须和参数向量θ一一对应。我习惯的写法是先用阶次na、nb声明维度，再循环逐行填充，这样逻辑清楚，阶次改动时不容易漏改。

% 批量最小二乘辨识 na = 2; nb = 2; % 模型阶次 L = max(na, nb); % 构造回归矩阵所需的最小起始点 Phi = zeros(N - L, na + nb); % 回归矩阵，每行对应一个时刻 Y = y(L+1:N); % 输出向量，与Phi行数对齐 for k = L+1:N Phi(k-L, :) = [-y(k-1) -y(k-2) u(k-1) u(k-2)]; % 注意列顺序与theta一致 end theta_hat = Phi \ Y; % 核心一步：反斜杠求解最小二乘问题 a1_hat = theta_hat(1); a2_hat = theta_hat(2); b1_hat = theta_hat(3); b2_hat = theta_hat(4); % 残差分析 resid = Y - Phi * theta_hat; % 预测残差 V = resid' * resid / length(Y); % 残差方差 fprintf('估计参数: a1=%.4f, a2=%.4f, b1=%.4f, b2=%.4f\n', ... a1_hat, a2_hat, b1_hat, b2_hat); fprintf('真实参数: a1=%.4f, a2=%.4f, b1=%.4f, b2=%.4f\n', a1, a2, b1, b2); fprintf('残差方差: %.6f\n', V);

这里有几个关键点必须说透。第一，回归矩阵从第L+1行开始构造，是因为前L个时刻缺少足够的y延迟项，硬凑只会引入无效行。第二，Phi的列顺序必须和参数向量的位置严格一致，我先把y的延迟排前面、u的延迟排后面，这样做的好处是参数解释直观。第三，theta_hat = Phi \ Y这一行是整个程序的核心，反斜杠在超定系统下自动走QR分解路径，比(Phi'*Phi)(Phi'*Y)多一层数值保护。跑完这段程序，你会看到估计参数与真实参数在小数点后两到三位吻合，残差方差非常接近0.01²，这就说明程序和模型结构都对。

3.3 结果验证：拟合优度与残差自相关

拿到参数不等于完事。我的习惯是立刻做两个检验：一是拟合优度，二是残差自相关。拟合优度的指标看R²，定义如下：

% 拟合优度检验 SST = sum((Y - mean(Y)).^2); % 总平方和 SSE = sum(resid.^2); % 残差平方和 R2 = 1 - SSE / SST; % R² 越接近1越好 fprintf('R² = %.6f\n', R2); % 残差自相关检验（粗略版） maxLag = 20; R_auto = zeros(maxLag, 1); for lag = 1:maxLag R_auto(lag) = resid(lag+1:end)' * resid(1:end-lag) / (resid' * resid); end figure; stem(R_auto); title('残差自相关函数'); xlabel('滞后');

R²接近1说明模型解释了绝大部分输出变化，但在辨识这件事上，R²好不代表模型好，因为ARX模型本身就包含过去的输出项，哪怕参数全错，只要数据平滑，R²也可能很高。所以我更看重残差自相关：如果残差在滞后1到20阶都接近0，说明残差是白噪声，模型已经把动态信息榨干了；如果某个滞后处明显超出±1.96/sqrt(N)的置信带，说明还有动态没建模进去，或者噪声本身有色。上面这段程序没有调用任何工具箱函数，全部用矩阵运算手写，方便你移植到自己的工程里。

4. 递推最小二乘：在线辨识程序与遗忘因子调参实战

4.1 从批量到递推：RLS更新公式的工程含义

批量最小二乘要求数据全部到位才能算，但很多场景数据是一点一点来的：在线监测、自适应控制、设备老化跟踪。这时批量法有两个硬伤：一是每来一个新数据就要重新求解一次矩阵方程，计算量随N线性增长；二是老数据和新数据权重相同，参数如果已经缓慢漂移，老数据会拖住估计值不放。递推最小二乘（RLS）就是为了解决这两个问题而存在的。

RLS的核心是一条递推更新链，每来一个新观测(k时刻的φ(k)和y(k))，就按下面的公式更新一次参数：

增益K(k) = P(k-1)·φ(k) / (λ + φ(k)ᵀ·P(k-1)·φ(k)) 参数θ(k) = θ(k-1) + K(k)·(y(k) - φ(k)ᵀ·θ(k-1)) 协方差P(k) = (I - K(k)·φ(k)ᵀ)·P(k-1) / λ

这里的λ是遗忘因子，取值范围通常在0.95到1之间。λ=1表示不遗忘，所有历史数据权重相同，RLS退化为普通最小二乘的递推形式；λ越小，旧数据被遗忘得越快，参数跟踪能力越强，但噪声也会跟着放大。工程直觉可以这样记：λ决定“记忆长度”，0.98大致相当于记得最近约50个点的信息，0.95大约20个点。这个对应关系不是精确的，但它能帮你把抽象的λ映射到实际需要的时间尺度上。

4.2 带遗忘因子的RLS的MATLAB实现

下面这段代码是完整的RLS辨识实现，直接接着第3章的仿真数据运行，每来一个数据点就更新一次参数，并把参数轨迹记录下来方便绘制收敛曲线。

% 递推最小二乘辨识，带遗忘因子 n = na + nb; % 待辨识参数个数 theta_rls = zeros(n, 1); % 参数初值，先给零向量 P = 1e6 * eye(n); % 协方差初值，大数代表初始不确定度高 lambda = 0.98; % 遗忘因子，0.95~0.99之间调 theta_hist = zeros(n, N); % 记录每步的参数，用于画收敛曲线 for k = L+1:N phi_k = [-y(k-1) -y(k-2) u(k-1) u(k-2)]'; % 当前回归向量，列向量 y_k = y(k); % 当前输出 K = P * phi_k / (lambda + phi_k' * P * phi_k); % 增益向量 e_k = y_k - phi_k' * theta_rls; % 预测误差 theta_rls = theta_rls + K * e_k; % 参数更新 P = (eye(n) - K * phi_k') * P / lambda; % 协方差更新 theta_hist(:, k) = theta_rls; % 记录 end % 绘制参数收敛曲线 figure; plot(theta_hist(1,:), 'b', 'LineWidth', 1.5); hold on; plot(theta_hist(2,:), 'r', 'LineWidth', 1.5); plot(theta_hist(3,:), 'g', 'LineWidth', 1.5); plot(theta_hist(4,:), 'm', 'LineWidth', 1.5); yline(a1, '--b'); yline(a2, '--r'); yline(b1, '--g'); yline(b2, '--m'); legend('a1','a2','b1','b2','真实值','Location','best'); xlabel('采样时刻 k'); ylabel('参数估计值'); title('RLS参数收敛曲线');

这段代码里有几个必须解释清楚的点。第一，P的初值取1e6·I，物理含义是“初始时对参数一无所知，置信度极低”，这样第一次更新时增益K会很大，参数迅速从零跳到大概范围。如果初值取1或更小，参数收敛会慢很多，甚至会一路缓爬。第二，分母里lambda加上一个标量φᵏᵀ·P·φₖ，这个标量本质是“当前回归向量在新息方向上的不确定性”，数值上永远是正数，所以不存在除零风险。第三，P的更新里除以λ，是为了放大协方差、保留对新数据的开放度，这也是遗忘因子的数学作用。

运行后你会看到参数曲线在前200个点快速收敛，之后在真实值附近小幅波动。波动幅度取决于λ和噪声方差：λ越接近1波动越小，但参数慢变时跟不上；λ越小波动越大，但能跟上缓变。这个矛盾就是调参的实质。

4.3 遗忘因子和初值的调参经验

实际调试RLS的时候，我一般按下面的顺序做。先固定λ=1跑一遍，看参数能不能收敛到稳定值，这一步验证模型结构和数据质量；确认没问题后，再把λ往小调，观察参数轨迹的跟踪速度和波动幅度。一个常见做法是让λ=0.98起步，如果参数曲线抖得厉害就提到0.99，如果发现参数明显滞后于实际变化就降到0.95。这个“步进搜索”的方式比直接猜一个值高效得多。

还有一个容易忽略的细节：RLS的收敛轨迹和输入信号的激励强度强相关。如果某段时间输入几乎不变（比如阀门停在某个开度），φ(k)的信息量不足，P会慢慢膨胀或停滞，此时参数会出现漂移。工程上我见过有人把λ调成0.9想加快跟踪，结果输入一停参数立刻乱跳。正确做法是在线监测预测误差e_k：如果e_k长期很小，说明系统处于稳定工况，此时应该暂停更新或把λ临时调回0.999，防止参数被噪声牵着走。

5. 最小二乘参数辨识避坑专题：5个我反复踩的坑

5.1 数据没去均值：直流分量把辨识结果整体带偏

现象：辨识出来的模型稳态增益和真实值差一大截，但残差看起来不大，R²还挺高。

原因：采集的y和u往往带有直流分量，也就是工作点不为零。ARX模型里如果没有显式的常数项c，回归方程Y = Φ·θ根本无法描述直流偏移，最小二乘只能拿“倾斜”的参数去硬凑这个偏移，结果就是动态参数都被污染。我在某温控系统辨识时就翻过这个车，气体阀门开度信号有2V偏置，温度信号还有1.5V偏置，直接辨识出的增益比真值低了将近40%。

解决：辨识前先做去均值处理，或者干脆用detrend函数去掉线性趋势。最稳妥的做法是分别计算u和y的均值并在构造回归矩阵前减掉。注意减均值要在同一时间段内计算，别用整段数据的全局均值去减前面的数据，那样会引入额外的斜坡失真。

5.2 用inv()解正规方程：病态矩阵下的“翻车”现场

现象：程序没报错，但算出来的参数数值大得离谱，正负号都在乱跳，比如a1估计出-3.7，而真实值是-0.8。

原因：直接用inv(Phi'*Phi)时，ΦᵀΦ的条件数是Φ的条件数的平方。现场数据各列之间天然有相关性，比如y(k-1)和y(k-2)在采样周期短时几乎是线性相关的，条件数轻易超过1e8，矩阵求逆的结果已经丢失了所有有效数字。这不是算法错了，是数值方法选错了。

解决：第一选择是phi = Phi \ Y，反斜杠自动处理病态问题。第二选择是pinv(Y)或对Phi做QR分解再回代。诊断手段是算一下cond(Phi'*Phi)，如果发现超过1e6，首先要怀疑的不是算法而是数据——要么激励不够，要么阶次选高了，要么数据段里有大段重复模式。修好数据问题后，条件数自然会降下来。

5.3 激励信号不丰富：阶跃信号辨识不出二阶参数

现象：用阶跃响应数据辨识二阶模型，b2估计值接近0，而且换一段数据结果差异极大，参数方差大得没法看。

原因：阶跃信号只有两个幅值状态，切换一次之后系统就进入稳态，之后的数据几乎不携带动态信息。从频域看，阶跃的能量集中在低频，二阶模型的极点信息大多在中高频段，激不出来。从数学上看，输入序列接近常值向量，回归矩阵的相关性极高，信息矩阵接近奇异，参数自然不可辨识。这是最小二乘最经典的失败模式之一。

解决：换激励信号。标准做法是用PRBS（伪随机二进制序列），它的频谱近似白噪声、幅值可控，对系统冲击小。MATLAB里如果没有系统辨识工具箱，可以手写一个简单的m序列生成器，或者退一步用randn生成随机幅值序列直接激励。注意激励信号的幅值要根据对象的线性范围来定，曾经有一个案例，某同事把PRBS幅值取得过大，系统进入了非线性区，辨识出的模型在后续控制器设计时完全不匹配。

% 简单的PRBS生成（备选方案） nBits = 10; N_prbs = 2^nBits - 1; % m序列长度 reg = ones(1, nBits); % 移位寄存器初值 prbs = zeros(N_prbs, 1); for p = 1:N_prbs prbs(p) = reg(1); fb = xor(reg(10), reg(1)); % 反馈抽头，不同抽头对应不同序列 reg = [fb, reg(1:end-1)]; end u = (prbs - 0.5) * 2; % 转成±1幅值序列

5.4 有色噪声使估计有偏：残差检验与补救方向

现象：参数估计和真值明显不一致，残差方差很小但自相关函数在低滞后处超出置信带。

原因：最小二乘的无偏性依赖“回归向量与噪声不相关”这个白噪声假设。如果测量噪声经过传感器低通滤波，或者系统有未建模的高频动态，残差就会带上相关性，此时θ̂的期望值不再是θ真值。这个问题在开环数据上也可能出现，闭环系统里则几乎是必现，因为反馈回路会把噪声耦合回输入。

解决：先做残差自相关检验来确认问题存在，然后换用辅助变量法（IV）或者增广最小二乘。IV的基本思路是找到一组与噪声不相关、但与信号强相关的工具变量来替代回归矩阵。后文第6章会给出具体代码。这里想强调一点：残差小不是模型好的充分条件，一定要做残差自相关检验，这是区分“模型解释得不错”和“模型统计上正确”的分水岭。

5.5 遗忘因子调太小：递推参数震荡的排查

现象：RLS刚开始收敛得挺好，跑到某一时刻参数突然剧烈震荡，甚至发散，把控制器输出都带抖了。

原因：遗忘因子设太小（比如0.9以下），模型只记住最近十几个点的信息，此时如果输入信号恰好处于低激励时段，信息矩阵接近奇异，参数估计的方差会急剧放大。这是个放大器效应：遗忘因子越大，对噪声越不敏感，但跟踪能力越差；太小则反过来。我见过某装置的在线辨识迫降旁路阀抖动，排查出来就是λ从0.98被运维人员改成了0.92。

解决：不要单独调λ，要同时观察预测误差e_k和输入信号的活跃程度。工程上我通常的做法是做一个“条件更新”：只有当输入变化超过某个阈值时才执行递推更新，否则只保持参数不变。这样既保留了遗忘因子的跟踪能力，又避开了低激励时段的参数漂移。此外，可以把λ设成0.97到0.995之间，别低于0.95。

6. 再进一步：辅助变量法与阶次辨识，让最小二乘类工具更耐打

6.1 用独立数据段验证辨识结果：习惯比算法更重要

做完辨识之后，我建议你养成一个固定习惯：把数据从中间切开，用前60%做辨识，后40%做验证。验证的做法很简单，把辨识出的参数代入模型，用后40%的输入重建输出序列，再与真实输出对比，计算验证集上的预测误差。这个方法成本极低，却能暴露两类问题：一是过拟合，训练集上R²很高、验证集上一塌糊涂，说明阶次偏高或数据量不够；二是时变，如果前段辨识的参数在后段明显失配，说明系统参数随时间漂移了，这时就该考虑RLS加遗忘因子，而不是抱着批量解不放。

% 用后40%数据做验证 N_val = round(N * 0.4); Y_val = y(N - N_val + 1 : N); U_val = u(N - N_val + 1 : N); y_pred = zeros(N_val, 1); for k = 3:N_val y_pred(k) = -a1_hat*y_pred(k-1) - a2_hat*y_pred(k-2) + ... b1_hat*U_val(k-1) + b2_hat*U_val(k-2); end val_err = Y_val - y_pred; fprintf('验证集残差标准差: %.4f\n', std(val_err));

注意这里验证时用的是预测输出y_pred自身的过去值，而不是验证集里真实的y过去值，这样才能真正检验模型的动态预测能力。如果只用真实y代入，那验证的就是一步预测能力，两者差距很大。这一步的代码很短，但我认为它是整套辨识流程里最能拦下废模型的关口。

6.2 辅助变量法：有色噪声下的5行代码补救

当残差检验确认噪声有色后，最快的补救手段是辅助变量法（IV）。核心思路是用一组与噪声不相关的“工具变量”Z替代回归矩阵Φ做投影。常见的工具变量选择是输入的延迟序列，因为输入u与外生白噪声不相关，而u的延迟与Φ中的信号成分强相关。实现代码非常短，构造一个Z矩阵再算一次投影即可。

% 辅助变量法：使用输入延迟作为工具变量 Z = zeros(N - L, na + nb); for k = L+1:N Z(k-L, :) = [u(k-1) u(k-2) u(k-3) u(k-4)]; % 工具变量 end theta_iv = (Z' * Phi) \ (Z' * Y); % IV估计

这段代码能工作的前提是Z与Φ列空间有足够的重合度，同时Z与噪声不相关。用u的延迟做工具变量，在开环数据下这两个条件都容易满足。替换掉原来的Phi \ Y之后，有偏问题会明显缓解。代价是方差比普通最小二乘略大，这是用偏差换方差的经典权衡。如果你有系统辨识工具箱，直接用辅助变量法命令会更方便，但手写这段代码能让你真正理解它在干什么。

6.3 阶次选择：在最小二乘框架里用信息准则定阶

模型阶次na、nb选错是辨识失败的重要原因，但选多少其实可以由数据告诉你。信息准则的做法是：对一组候选阶次分别做批量最小二乘，计算残差方差V，再按下面的公式打分，取最小值对应的阶次。

% 阶次辨识：AIC / BIC 信息准则 maxOrder = 4; AIC_mat = zeros(maxOrder, maxOrder); BIC_mat = zeros(maxOrder, maxOrder); for na_t = 1:maxOrder for nb_t = 1:maxOrder L_t = max(na_t, nb_t); Phi_t = zeros(N - L_t, na_t + nb_t); for k = L_t+1:N Phi_t(k-L_t, :) = [-y(k-1:-1:k-na_t) u(k-1:-1:k-nb_t)]; end Y_t = y(L_t+1:N); theta_t = Phi_t \ Y_t; resid_t = Y_t - Phi_t * theta_t; V_t = resid_t' * resid_t / length(Y_t); npar = na_t + nb_t; AIC_mat(na_t, nb_t) = length(Y_t) * log(V_t) + 2 * npar; BIC_mat(na_t, nb_t) = length(Y_t) * log(V_t) + npar * log(length(Y_t)); end end [minAIC, idxAIC] = min(AIC_mat(:)); [minBIC, idxBIC] = min(BIC_mat(:)); [na_AIC, nb_AIC] = ind2sub(size(AIC_mat), idxAIC); [na_BIC, nb_BIC] = ind2sub(size(BIC_mat), idxBIC);

这里有个使用前提：噪声必须是白噪声，否则V的估计本身有偏，信息准则会偏好高阶模型。所以在定阶之前，先确认你已经处理好噪声相关性问题。另一个值得注意的坑是，仿真数据里噪声设成0时，V会小到让log(V)变成负无穷，信息准则直接失效。我自己的调试经历里不止一次被这个问题坑过，所以仿真时噪声方差不要低于1e-6。把这一段跑通后，你手头的最小二乘类参数辨识工具链就齐了：批量LS打底、RLS应对在线、IV处理有色噪声、信息准则定阶，再加上交叉验证习惯，这套组合能覆盖绝大多数线性系统辨识需求。希望你也能少走几步弯路，一次就把参数辨识到位。

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