1. 从“黑箱”到“白箱”：系统辨识的工程价值

在工业控制、信号处理乃至经济建模的日常工作中，我们常常面对一个“黑箱”：你给它一个输入信号，它吐出一个输出信号，但箱子内部的结构、参数，你一无所知。比如，你想预测一个化学反应器的温度变化，或者想设计一个无人机的飞控算法，你首先得知道这个系统“听不听话”——给它一个推力，它会以多快的速度响应？这个响应过程是平滑的还是震荡的？系统辨识，就是把这“黑箱”变成“白箱”的科学与艺术。它不要求你拆开物理设备，而是通过输入输出数据，用数学方法“猜”出系统内部的动态规律，也就是数学模型。

而“最小二乘”，无疑是打开这扇门最经典、最实用的一把钥匙。它不是什么高深莫测的理论，其核心思想朴素得惊人：找一条线（或一个曲面），让所有数据点到这条线的“距离”平方和最小。在系统辨识里，这条“线”就是我们的候选模型，那些“距离”就是模型预测输出与实际观测输出之间的误差。最小二乘法告诉我们，哪个模型参数能让预测和现实最贴合。这门“最小二乘系统辨识课”的上篇，我们就聚焦于最基础、也最核心的部分：如何为你的系统选择一个合适的“辨识模型”，以及最小二乘法是如何在这个模型框架下大显身手的。无论你是自动化专业的学生，还是从事算法开发的工程师，理解这套基础范式，都是后续处理更复杂非线性、时变系统问题的基石。

2. 模型选择：给系统“画像”的第一步

在动手算参数之前，选对模型结构至关重要。这就像画画，你得先决定是用素描还是油画，是写实还是抽象。选错了，后续再精妙的算法也是徒劳。在经典的系统辨识中，主要有两大类模型结构：方程误差模型和输出误差模型。理解它们的区别，是避免后续踩坑的关键。

2.1 方程误差模型：最直接的线性回归视角

方程误差模型，也叫ARX模型，是入门系统辨识最先接触的结构。它的形式非常直观： A(z)y(k) = B(z)u(k) + e(k) 这里， y(k) 是系统输出， u(k) 是系统输入， e(k) 是一个白噪声序列。 A(z) 和 B(z) 是后移算子 z^{-1} 的多项式。展开写，就是一个差分方程： y(k) + a1*y(k-1) + ... + ana*y(k-na) = b1*u(k-1) + ... + bnb*u(k-nb) + e(k) 这个方程告诉我们，当前的输出 y(k) ，是由过去若干时刻的输出 y(k-1)... 、过去若干时刻的输入 u(k-1)... 以及一个当前的噪声 e(k) 共同决定的。

为什么它如此受欢迎？ 核心原因在于，这个模型关于待估参数 a1, a2, ..., b1, b2, ... 是 线性 的。我们把方程稍微变形一下： y(k) = -a1*y(k-1) - ... - ana*y(k-na) + b1*u(k-1) + ... + bnb*u(k-nb) + e(k) 现在，我们把等号右边除了 e(k) 之外的所有项，看作是一组已知的“特征”： -y(k-1) , -y(k-2) , ..., u(k-1) , u(k-2) ...。而 a1, a2, ..., b1, b2, ... 就是这些特征对应的权重系数。看，这完美地契合了多元线性回归 y = θ1*x1 + θ2*x2 + ... + θn*xn + e 的形式！这意味着，我们可以直接套用成熟、高效、解析解明确的最小二乘法来估计参数，计算简单，全局最优。

它的“阿喀琉斯之踵”：噪声假设 方程误差模型有一个很强的假设：噪声 e(k) 是直接加在方程等式两端的。在实际的物理系统中，噪声往往不是这样作用的。更常见的场景是，噪声影响了系统内部状态，然后经过系统自身的动力学特性后才体现在输出上。这就引出了方程误差模型的一个主要缺点：当真实系统的噪声通道与模型假设不符时，即使数据量无限，其参数估计值也会存在偏差，无法收敛到真实值。这种现象在辨识理论中称为“有偏估计”。因此，ARX模型更适合对模型精度要求不是极端苛刻，或者噪声影响较小的初步建模场景。

2.2 输出误差模型：更贴近物理现实的描述

为了克服方程误差模型的偏差问题，输出误差模型（OE模型）被提出。它的结构更贴近我们对许多物理系统的认知： y(k) = [B(z)/F(z)] * u(k) + e(k) 这里， [B(z)/F(z)] 代表一个传递函数，描述了输入 u 到系统“无噪声输出” x(k) 的动态过程。而最终的观测输出 y(k) ，是这个无噪声输出 x(k) 加上一个独立的观测噪声 e(k) 。用差分方程写出来是： x(k) + f1*x(k-1) + ... + fnf*x(k-nf) = b1*u(k-1) + ... + bnb*u(k-nb) y(k) = x(k) + e(k)

模型优势与代价 OE模型清晰地分离了系统的确定性动态（ B/F ）和随机性干扰（ e(k) ）。只要 e(k) 是白噪声，理论上通过一些迭代算法（如预报误差法）可以得到无偏的、一致的参数估计。这听起来很完美，对吧？但代价是：模型关于参数 b1, b2, ..., f1, f2, ... 是 非线性 的。因为无噪声输出 x(k) 本身依赖于参数 f ，这使得模型无法写成关于所有参数的线性回归形式。我们不能再用简单的最小二乘直接求解，而必须借助非线性优化算法（如高斯-牛顿法、梯度下降法）进行迭代搜索，计算复杂，且可能陷入局部最优。

注意 ：在实际工程中，模型结构的选择（ARX还是OE，以及阶次 na, nb, nf 如何确定）本身就是一个重要课题。通常建议从简单的ARX模型开始，利用其计算快的优势进行模型阶次的试探和初步分析，如果残差序列（预测误差）表现出明显的相关性，说明ARX的噪声假设不成立，再考虑使用OE等更复杂的模型。不要一开始就追求复杂模型，简单有效永远是第一原则。

3. 最小二乘原理：误差平方和最小化的几何与统计意义

现在，我们以最经典的ARX模型为例，深入最小二乘法的内核。假设我们通过实验，采集到了一组长度为N的输入输出数据 {u(1), y(1)}, {u(2), y(2)}, ..., {u(N), y(N)} 。我们选定模型阶次 na 和 nb ，想要估计参数向量 θ = [a1, a2, ..., ana, b1, b2, ..., bnb]^T 。

根据ARX模型，对于第 k 个采样时刻（ k > max(na, nb) ），我们可以写出： y(k) = φ(k)^T θ + e(k) 其中， φ(k) 称为 回归向量 或 数据向量 ： φ(k) = [-y(k-1), -y(k-2), ..., -y(k-na), u(k-1), u(k-2), ..., u(k-nb)]^T 你看，这个式子把非线性（关于 y 和 u ）的系统动力学，巧妙地转化成了关于参数 θ 的线性形式。对于所有 k = L+1 到 N （ L = max(na, nb) ），我们可以构建一个庞大的线性方程组：

y(L+1) = φ(L+1)^T θ + e(L+1)
y(L+2) = φ(L+2)^T θ + e(L+2)
...
y(N)   = φ(N)^T θ   + e(N)

写成矩阵形式，简洁有力： Y = Φ θ + E 其中：

• Y = [y(L+1), y(L+2), ..., y(N)]^T 是输出向量。
• Φ = [φ(L+1), φ(L+2), ..., φ(N)]^T 是数据矩阵（也叫信息矩阵或回归矩阵）。
• E = [e(L+1), e(L+2), ..., e(N)]^T 是噪声向量。

我们的目标是找到一组参数 θ ，使得模型预测值 Φθ 尽可能接近真实观测值 Y 。最小二乘准则规定：这个“接近”的程度，用所有误差的平方和来衡量，即代价函数 J(θ) = E^T E = (Y - Φθ)^T (Y - Φθ) 。我们要找到使 J(θ) 最小的那个 θ 。

从几何角度理解 ：向量 Y 存在于一个N维空间中。矩阵 Φ 的列向量张成了一个子空间（列空间）。模型预测值 Φθ 是这个子空间中的一个点。最小二乘解的意义在于，在子空间中寻找一个点 Φθ ，使得它到真实点 Y 的欧几里得距离最短。根据几何知识，这个最短距离是通过 Y 向子空间做 正交投影 得到的。因此，最优的 Φθ 是 Y 在 Φ 列空间上的投影，而误差向量 E = Y - Φθ 垂直于该列空间。

从统计角度理解 ：如果我们假设噪声 e(k) 是零均值、同方差且互不相关的白噪声，那么最小二乘估计量具有一系列优良性质：它是无偏的（估计值的期望等于真值），并且是所有无偏估计中方差最小的（有效估计）。这使得最小二乘在统计意义上也是最优的。

求解这个最优化问题，可以通过对代价函数 J(θ) 求关于 θ 的梯度，并令其为零： ∇J(θ) = -2Φ^T (Y - Φθ) = 0 由此得到著名的 正规方程 ： (Φ^T Φ) θ = Φ^T Y 只要矩阵 Φ^T Φ 是可逆的（这要求数据足够丰富，且输入信号具有一定的激励性，即持续激励条件），我们就可以得到最小二乘参数估计的解析解： θ_LS = (Φ^T Φ)^{-1} Φ^T Y 这个公式干净利落，是所有系统辨识、机器学习线性回归问题的基石。

4. 实操核心：数据矩阵构建与持续激励条件

理论很优美，但落地到代码和实验，有两个细节决定成败：如何正确构建数据矩阵 Φ ，以及如何保证 Φ^T Φ 可逆。很多初学者在这里栽跟头。

4.1 数据矩阵构建的“边界”问题

构建 Φ 矩阵时，第一个实际问题是 数据索引的起始点 。注意我们的回归向量 φ(k) 包含了 y(k-1), ..., y(k-na) 和 u(k-1), ..., u(k-nb) 。这意味着，要构造 φ(L+1) ，我们需要 y(L) 和 u(L) 的数据，而 L = max(na, nb) 。因此， 有效的数据段是从 k = L+1 开始，到 k = N 结束 。你用于拟合的数据长度实际上是 N - L ，而不是 N 。在编程时，一个常见的错误是直接从 k=1 开始循环构造 φ(k) ，导致数组越界。正确的做法是：

import numpy as np
# 假设 y_data, u_data 是长度为 N 的数组
na, nb = 2, 2
L = max(na, nb)
N = len(y_data)
Y = y_data[L:]  # 从索引L到末尾
Phi = []
for k in range(L, N):
    phi_k = np.concatenate([-y_data[k-1:k-na-1:-1], u_data[k-1:k-nb-1:-1]])
    Phi.append(phi_k)
Phi = np.array(Phi)
# 然后求解 theta = np.linalg.inv(Phi.T @ Phi) @ Phi.T @ Y

这里 y_data[k-1:k-na-1:-1] 利用了Python切片来获取 [y(k-1), y(k-2), ..., y(k-na)] ，注意顺序。

4.2 持续激励：让数据“会说话”

第二个，也是更本质的问题是： 什么样的输入数据 u(k) ，才能让我们唯一地、可靠地辨识出所有参数？ 答案就是输入信号必须满足“持续激励”条件。

直观理解：如果你用一个恒定值（比如 u(k)=1 ）去激励系统，你只能得到系统在某个静态工作点附近的信息，无法分辨出系统动态（ a1, a2,... ）中不同模式的影响。这就像你想了解一个弹簧的质量和阻尼系数，只把它压住不动是测不出来的，必须用不同频率的力去推拉它。

从数学上看， Φ^T Φ 的可逆性要求数据矩阵 Φ 是列满秩的。对于ARX模型， Φ 的列由过去的输出和输入数据组成。而过去的输出 y(k-i) 本身又依赖于过去的输入和参数。因此，归根结底，要求输入信号 u(k) 能充分激发系统所有模态。一个经典且实用的选择是 伪随机二进制序列 。它看起来像随机的0/1跳变，但具有周期性和良好的自相关特性，能在一个宽频带内提供近似白噪声的激励，是系统辨识实验中常用的输入信号。

实操心得 ：在仿真中，你可以用PRBS或白噪声作为输入。但在实际物理系统测试中，输入信号必须考虑系统的安全限幅和执行器的物理限制。一个折中的好办法是使用幅值受限的、不同频率的正弦扫频信号，或者幅值随机变化的阶跃信号序列。关键是要让输入有足够丰富的变化，覆盖你关心的频率范围。记录数据时，务必确保输入输出数据是同步采集的，并注意剔除明显的野值。

5. 评估与诊断：你的模型“合格”了吗？

参数 θ 算出来了，模型就有了。但模型质量如何？不能只靠感觉，需要有量化的评估和诊断工具。这里介绍三个最实用的方法。

5.1 拟合优度：量化匹配程度

最直接的指标是 拟合优度 ，通常用归一化的均方误差来表示，例如： FIT = (1 - norm(Y - Y_pred) / norm(Y - mean(Y))) * 100% 其中 Y_pred = Φ θ_LS 是模型的预测输出。 FIT 越接近100%，说明模型对这段训练数据的解释能力越强。但要注意，高拟合优度可能意味着 过拟合 ——模型不仅拟合了系统动态，还拟合了数据中特定的噪声。因此，它通常用于同一组数据上不同模型结构的横向比较，而不是绝对质量的评判。

5.2 残差分析：检验模型假设的“试金石”

残差，就是预测误差序列 ε(k) = y(k) - y_pred(k) 。如果我们的模型（包括其噪声假设）是完美的，那么残差序列应该是一个 白噪声 序列——均值为零，序列自身不同时刻的值互不相关。

如何检验？

1. 绘制残差序列图 ：肉眼观察是否围绕零均值线随机波动，有无明显的趋势或周期性。
2. 计算残差的自相关函数 ：对于一个白噪声，其自相关函数在时滞 τ≠0 时应接近于零。我们可以计算残差的自相关函数，并观察其是否落在95%的置信区间内。如果很多点落在区间外，特别是前几个时滞的相关性显著不为零，则说明残差中存在未建模的动态信息，模型结构（如阶次）可能选择不当。
3. 残差与输入信号的互相关函数 ：一个理想的模型，其残差应与过去的输入信号不相关。如果互相关函数显著不为零，说明模型未能完全捕获输入到输出的动态关系，或者存在非线性未建模。

残差分析是系统辨识中 极其重要 的一步，它比单纯的拟合优度更能揭示模型的本质问题。我个人的经验是，一个 FIT 只有85%但残差接近白噪声的模型，通常比一个 FIT 高达95%但残差自相关严重的模型更可靠、更具泛化能力。

5.3 交叉验证：防范过拟合的黄金准则

最终极的检验，是将模型用在它“没见过”的数据上。这就是 交叉验证 。具体做法：

1. 将你的数据集分为两部分：一部分用于参数估计（训练集），另一部分用于模型验证（测试集）。
2. 只用训练集的数据来构建 Φ_train 和 Y_train ，并计算参数 θ_LS 。
3. 锁定这个 θ_LS ，将其应用到测试集上。用测试集的输入数据和模型公式（注意，此时计算预测输出 y_pred_test 时，使用的过去输出值应是测试集真实的 y 值，而不是模型自己递归预测的值，这称为“仿真验证”模式），得到测试集上的预测输出。
4. 计算模型在测试集上的拟合优度或误差指标。

如果模型在训练集上表现很好，但在测试集上表现大幅下降，这就是典型的过拟合。交叉验证能最真实地反映模型对新数据的预测能力，是评估模型泛化性能的黄金标准。在实际项目中，我强烈建议至少保留30%的数据作为测试集。

6. 一个完整的仿真案例：二阶离散系统辨识

让我们用一个具体的MATLAB/Python仿真例子，把上述所有步骤串起来。假设真实系统是一个二阶离散系统： y(k) - 1.5y(k-1) + 0.7y(k-2) = u(k-1) + 0.5u(k-2) + e(k) 其中 e(k) 是方差为0.1的高斯白噪声。

步骤1：数据生成

% MATLAB
N = 1000;
u = randn(N, 1); % 使用白噪声作为输入，满足持续激励
e = 0.1 * randn(N, 1);
y = zeros(N,1);
y(1:2) = [0; 0]; % 初始条件
for k = 3:N
    y(k) = 1.5*y(k-1) - 0.7*y(k-2) + u(k-1) + 0.5*u(k-2) + e(k);
end
% 划分训练集和测试集
train_ratio = 0.7;
N_train = floor(N * train_ratio);
u_train = u(1:N_train); y_train = y(1:N_train);
u_test = u(N_train+1:end); y_test = y(N_train+1:end);

# Python
import numpy as np
N = 1000
np.random.seed(42)
u = np.random.randn(N)  # 白噪声输入
e = 0.1 * np.random.randn(N)
y = np.zeros(N)
y[:2] = [0, 0]
for k in range(2, N):
    y[k] = 1.5*y[k-1] - 0.7*y[k-2] + u[k-1] + 0.5*u[k-2] + e[k]
# 划分数据集
train_ratio = 0.7
N_train = int(N * train_ratio)
u_train, y_train = u[:N_train], y[:N_train]
u_test, y_test = u[N_train:], y[N_train:]

步骤2：模型阶次选择与数据矩阵构建 我们猜测系统可能是二阶的（ na=2, nb=2 ）。用训练集数据构建ARX模型的数据矩阵。

na, nb = 2, 2
L = max(na, nb)
# 构建训练集数据矩阵
Y_train = y_train[L:]
Phi_train = []
for k in range(L, len(y_train)):
    phi_k = np.concatenate([-y_train[k-1:k-na-1:-1], u_train[k-1:k-nb-1:-1]])
    Phi_train.append(phi_k)
Phi_train = np.array(Phi_train)

步骤3：最小二乘参数估计

# 求解正规方程
theta_LS = np.linalg.inv(Phi_train.T @ Phi_train) @ Phi_train.T @ Y_train
print(f"估计参数: {theta_LS}")
print(f"真实参数: a1=-1.5, a2=0.7, b1=1.0, b2=0.5")
# 注意：我们回归方程写的是 y(k) = -a1*y(k-1) -a2*y(k-2) + b1*u(k-1) + b2*u(k-2) + e(k)
# 所以 theta_LS 的前两个元素对应 -a1, -a2，后两个对应 b1, b2
a1_est, a2_est, b1_est, b2_est = -theta_LS[0], -theta_LS[1], theta_LS[2], theta_LS[3]
print(f"估计的差分方程: y(k) = {a1_est:.4f}*y(k-1) + {a2_est:.4f}*y(k-2) + {b1_est:.4f}*u(k-1) + {b2_est:.4f}*u(k-2)")

运行后，估计参数应该非常接近真实值 [-1.5, 0.7, 1.0, 0.5] ，但由于噪声存在，会有微小偏差。

步骤4：模型评估与诊断

# 1. 训练集拟合
y_pred_train = Phi_train @ theta_LS
fit_train = 100 * (1 - np.linalg.norm(y_train[L:] - y_pred_train) / np.linalg.norm(y_train[L:] - np.mean(y_train[L:])))
print(f"训练集拟合优度: {fit_train:.2f}%")

# 2. 测试集验证 (仿真模式)
# 注意：测试集仿真需要递归计算，因为每一步的预测输出会作为下一步的输入
y_pred_test = np.zeros_like(y_test)
y_pred_test[:L] = y_test[:L] # 用真实值初始化
for k in range(L, len(y_test)):
    # 使用模型和过去的“预测输出”及“真实输入”
    phi_k_test = np.concatenate([-y_pred_test[k-1:k-na-1:-1], u_test[k-1:k-nb-1:-1]])
    y_pred_test[k] = phi_k_test @ theta_LS
fit_test = 100 * (1 - np.linalg.norm(y_test[L:] - y_pred_test[L:]) / np.linalg.norm(y_test[L:] - np.mean(y_test[L:])))
print(f"测试集拟合优度: {fit_test:.2f}%")

# 3. 残差分析 (以训练集残差为例)
residual = y_train[L:] - y_pred_train
# 计算残差自相关函数 (这里简化为计算前20个时滞)
max_lag = 20
acf = np.correlate(residual, residual, mode='full')
acf = acf[len(acf)//2 : len(acf)//2 + max_lag + 1]
acf = acf / acf[0] # 归一化
# 绘制自相关图 (略)，观察是否在置信带内。

通过这个完整流程，你不仅得到了模型参数，还通过测试集验证和残差分析，对模型质量有了定量和定性的评估。如果测试集拟合度尚可且残差接近白噪声，那么这个模型就可以用于后续的控制器设计或系统分析了。

7. 局限与进阶思考：最小二乘辨识的边界

通过上篇的讨论，我们掌握了基于最小二乘和ARX模型进行系统辨识的完整流程。然而，在实际工程中，你会很快遇到它的边界。认识到这些局限，正是迈向更高级辨识方法（如递推最小二乘、广义最小二乘、预报误差法等）的起点。

首先，是噪声模型的局限 。ARX模型假设噪声直接加在方程输出端，这在实际中往往不成立。当噪声通过系统动力学环节（即有色噪声）时，标准最小二乘估计是有偏的。这时就需要引入更复杂的模型，如ARMAX（方程误差模型但噪声为移动平均过程）或BJ（Box-Jenkins）模型，并使用广义最小二乘或极大似然等方法来获得无偏估计。

其次，是时变系统的挑战 。我们假设系统参数是定常的。但对于缓慢变化或突变的系统（如化学反应器催化剂活性衰减、飞机在不同空速下的气动参数），需要能够在线跟踪参数变化的算法，这就是 递推最小二乘 （RLS）及其变种（如带遗忘因子的RLS）的用武之地。

再者，是关于线性与静态的假设 。最小二乘辨识本质上是线性回归，它只能辨识线性系统的参数。对于非线性系统，你需要选择非线性模型结构（如Hammerstein模型、Wiener模型、神经网络等），并使用非线性优化方法进行参数估计，计算复杂度和初始值敏感性会急剧增加。

最后，我想分享一个最容易被忽视的心得： 系统辨识不仅仅是数学和算法，更是对物理对象的理解与实验设计的艺术 。再精巧的算法，如果输入数据质量差（激励不足、测量噪声大、采样不同步），也得不到好模型。在启动辨识算法之前，请务必花时间思考：我的输入信号能激发系统的所有重要模式吗？我的采样频率足够高吗（满足香农定理）？我的传感器数据可靠吗？很多时候，花在改进实验设计和数据预处理上的时间，远比调试算法参数带来的回报大得多。最小二乘给了我们一个强大的工具，但用好这个工具的前提，是深刻理解你的系统和你的数据。