MATLAB实战：如何用Hankel矩阵法搞定三阶系统辨识（附完整代码）

从脉冲响应到状态空间：Hankel矩阵法在MATLAB中的工程化系统辨识实践

在控制工程和信号处理的实际项目中，我们常常会遇到一个经典问题：面对一个物理系统或“黑箱”，如何仅通过其输入输出数据，反向推导出描述其内部动态的数学模型？这就是系统辨识的核心任务。对于从事自动化、机器人或振动分析领域的工程师和研究者而言，掌握一种稳健、可编程的辨识方法，意味着能将实验数据快速转化为可仿真、可控制的设计模型，极大地缩短研发周期。

Hankel矩阵法，正是这样一把连接理论与实践的钥匙。它绕开了复杂的频域分析，直接从系统的脉冲响应序列出发，通过构造具有特殊结构的Hankel矩阵，巧妙地提取出系统的阶次与参数。这种方法理论优美，计算清晰，特别适合在MATLAB这类数值计算环境中实现。今天，我们不打算重复教科书上的矩阵推导，而是聚焦于如何将这一方法工程化——从数据采集的坑点、代码实现的细节，到结果分析的技巧，手把手带你走完一个完整的三阶系统辨识流程。你会发现，理论上的“完美辨识”与实际操作中的“误差控制”，完全是两回事。

1. 工程起点：理解脉冲响应与采样周期的博弈

在开始写代码之前，我们必须先理清一个基础但至关重要的问题：从哪里获得可靠的脉冲响应数据？对于很多新手来说，直接对真实物理系统施加一个理想的脉冲激励是不现实的，甚至可能损坏设备。因此，仿真环境下的“数字孪生”成为了学习和验证算法的首选。

我们以一个典型的三阶系统作为辨识对象，其传递函数为： [ G(s) = \frac{s + 15}{s^3 + 5s^2 + 4s + 15} ] 在MATLAB中，我们可以轻松地创建这个系统模型，并获取其离散化后的脉冲响应。这里就引出了第一个工程决策点：采样周期T0的选择。

注意：采样周期并非越小越好。过小的T0会导致数据量剧增、矩阵条件数变差，并可能引入数值计算问题；过大的T0则会因信息丢失（混叠）而导致辨识失败。

为了直观感受采样周期的影响，我们可以先进行一个简单的预实验。下面的代码生成了不同采样周期下，系统脉冲响应的对比图。

% 定义原连续系统 b0 = 1; b1 = 15; a0 = 1; a1 = 5; a2 = 4; a3 = 15; sys_continuous = tf([b0 b1], [a0 a1 a2 a3]); % 设置不同的采样周期 T0_list = [0.05, 0.2, 0.5, 1.0]; figure('Position', [100, 100, 1200, 800]); for i = 1:length(T0_list) T0 = T0_list(i); sys_discrete = c2d(sys_continuous, T0, 'zoh'); % 使用零阶保持器离散化 subplot(2,2,i); impulse(sys_discrete, 10); % 绘制前10个采样点的脉冲响应 grid on; title(sprintf('采样周期 T0 = %.2f 秒', T0)); xlabel('时间 (秒)'); ylabel('幅度'); end sgtitle('不同采样周期下的系统离散脉冲响应');

运行这段代码，你会立刻观察到脉冲响应形态的显著变化。T0=0.05时，响应曲线光滑，细节丰富；而T0=1.0时，曲线变得稀疏，很多动态细节已经丢失。这直观地告诉我们，采样周期直接决定了我们用于构建Hankel矩阵的“原材料”的质量。

在实际工程中，一个经验法则是采样频率至少为系统预期最高频率的10倍。对于这个三阶系统，我们可以通过damp(sys_continuous)命令查看其极点频率，从而估算一个合理的采样范围。

2. Hankel矩阵的构造：从数据序列到系统信息的封装

拿到了脉冲响应序列g(1), g(2), g(3), ...之后，下一步就是构建Hankel矩阵。这是整个算法的核心数据结构。一个p x q的Hankel矩阵H定义为其逆对角线上的所有元素相等，更具体地说，它由序列的前p+q-1项填充而成，其第i行第j列的元素为g(i+j-1)。

对于一个n阶系统，我们通常构造一个n x n的方阵H(n,1)，有时也称为“能观性”或“能控性”Hankel矩阵。其形式如下：

[ H(n,1) = \begin{bmatrix} g(1) & g(2) & \cdots & g(n) \ g(2) & g(3) & \cdots & g(n+1) \ \vdots & \vdots & \ddots & \vdots \ g(n) & g(n+1) & \cdots & g(2n-1) \end{bmatrix} ]

这个矩阵的神奇之处在于，它的秩理论上就等于系统的阶次n，并且其元素中蕴含了足以重构系统状态空间模型（A, B, C矩阵）的全部信息。在我们的辨识场景中，目标是求出传递函数的系数，这可以通过求解一个由Hankel矩阵构成的线性方程组来实现。

在MATLAB中，构造这个矩阵非常简洁。假设我们已经从离散系统的脉冲响应中提取了足够多的数据点，存储于向量g中（g(1)对应第一个采样点的响应值）。

% 假设已知系统阶次 n=3，且已获得脉冲响应序列 impulse_response % impulse_response 是一个列向量，长度至少为 2*n n = 3; % 系统阶次 g = impulse_response(1:2*n); % 取前2n个点，实际可能需要更多 % 构造 n x n 的 Hankel 矩阵 H H = zeros(n, n); for i = 1:n for j = 1:n H(i, j) = g(i + j - 1); end end % 更高效的向量化构造方法： % i = (1:n)'; j = 1:n; % H = g(i + j - 1);

然而，这里隐藏着一个工程上的关键点：我们如何确定系统的阶次n？在实际辨识中，n往往是未知的。一种实用的方法是构建一个较大的Hankel矩阵，然后通过分析其奇异值（Singular Value Decomposition, SVD）来确定有效阶次。

% 探索性阶次判定 max_n = 10; % 假设最大可能阶次 L = length(g); H_large = zeros(max_n, max_n); for i = 1:max_n for j = 1:max_n if (i+j-1) <= L H_large(i, j) = g(i+j-1); else H_large(i, j) = 0; % 数据不足时补零 end end end % 进行奇异值分解 [U, S, V] = svd(H_large); singular_values = diag(S); % 绘制奇异值谱 figure; stem(1:max_n, singular_values(1:max_n), 'filled', 'LineWidth', 2); grid on; xlabel('奇异值序号'); ylabel('奇异值大小'); title('Hankel矩阵的奇异值谱'); set(gca, 'YScale', 'log'); % 对数坐标更容易观察阶跃

在奇异值谱图中，你会看到前几个奇异值显著大于后面的值，形成一个明显的“拐点”。拐点之后奇异值接近零或数量级骤降，这个拐点对应的序号通常就是系统有效阶次n的一个良好估计。这个方法比盲目假设阶次要可靠得多。

3. 算法实现：求解系数与离散到连续的转换

确定了阶次n并构造好Hankel矩阵H后，我们就可以根据算法推导出的公式来求解系统参数。核心是求解以下两个线性方程组：

1. 求解分母多项式系数a_i： [ H \cdot \begin{bmatrix} a_n \ a_{n-1} \ \vdots \ a_1 \end{bmatrix} = - \begin{bmatrix} g(n+1) \ g(n+2) \ \vdots \ g(2n) \end{bmatrix} ] 这里a_i对应离散传递函数分母1 + a_1 z^{-1} + ... + a_n z^{-n}的系数。

2. 求解分子多项式系数b_i： [ \begin{bmatrix} b_1 \ b_2 \ \vdots \ b_n \end{bmatrix} = - T \cdot \begin{bmatrix} g(1) \ g(2) \ \vdots \ g(n) \end{bmatrix} ] 其中矩阵T是一个下三角Toeplitz矩阵，由分母系数a_i构成。这一步将脉冲响应与分子系数联系起来。

下面的代码块完整实现了这一过程，并包含了从离散传递函数（Z域）转换回连续传递函数（S域）的关键步骤。

function [sys_identified_continuous, sys_identified_discrete] = hankel_identify(sys_original, T0, n) % HANKEL_IDENTIFY 使用Hankel矩阵法进行系统辨识 % 输入： % sys_original: 原始连续系统 (tf 或 ss 对象)，用于生成脉冲响应数据 % T0: 采样周期 % n: 预设的系统阶次 % 输出： % sys_identified_continuous: 辨识得到的连续系统 % sys_identified_discrete: 辨识得到的离散系统 % 1. 离散化原系统并获取脉冲响应 sysd_original = c2d(sys_original, T0, 'zoh'); [g, t_imp] = impulse(sysd_original); % 确保有足够的数据点 (至少 2n 个) N_needed = 2*n; if length(g) < N_needed warning('脉冲响应数据点不足，尝试延长仿真时间或减小采样周期。'); % 可以在这里扩展仿真时间重新获取脉冲响应 [g, t_imp] = impulse(sysd_original, (N_needed-1)*T0); end g = g(:); % 确保是列向量 % 2. 构造 Hankel 矩阵 H 和向量 h H = zeros(n, n); for i = 1:n for j = 1:n H(i, j) = g(i + j - 1); end end h = -g(n+1 : 2*n); % 等式右边的向量 % 3. 求解分母系数 a (注意：方程是 H * a = h) % 使用反斜杠运算符求解，它比直接求逆更稳定 a_coeff = H \ h; % a_coeff = [a_n; a_{n-1}; ...; a_1] % 4. 构造 Toeplitz 矩阵 T 并求解分子系数 b T = eye(n); % 初始化单位矩阵 for i = 2:n for j = 1:i-1 T(i, j) = a_coeff(n - (i - j) + 1); % 根据a_coeff的顺序填充下三角 end end b_coeff = -T \ g(1:n); % b_coeff = [b_1; b_2; ...; b_n] % 5. 构建辨识出的离散传递函数 % 注意MATLAB中tf函数的系数顺序：降幂排列 numd_identified = [0; b_coeff]'; % 分子，补上 b_0=0（因为通常g(0)=0） dend_identified = [1; flipud(a_coeff)]'; % 分母，补上 a_0=1 sys_identified_discrete = tf(numd_identified, dend_identified, T0); % 6. 使用双线性变换（Tustin变换）将离散系统转换回连续系统 sys_identified_continuous = d2c(sys_identified_discrete, 'tustin'); % 可选：输出辨识结果 fprintf('辨识完成。\n'); fprintf('离散系统传递函数系数：\n'); fprintf(' 分子: '); disp(numd_identified); fprintf(' 分母: '); disp(dend_identified); end

这段代码有几个值得强调的工程细节：

• 数据长度检查：确保脉冲响应序列长度至少为2n，否则矩阵方程无法构建。
• 稳定求解：使用H \ h（MATLAB的左除运算）代替inv(H)*h，前者基于QR或LU分解，数值稳定性更高，尤其当H接近奇异时。
• 系数顺序：MATLAB的tf函数要求多项式系数按s或z的降幂排列。我们求解出的a_coeff是[a_n; ...; a_1]，需要翻转并补上首项1才能构成分母多项式。
• 双线性变换：d2c(..., 'tustin')是将离散模型转换回连续模型的常用方法，它能保持稳定性，但在高频段存在频率畸变。这也是后续误差的一个来源。

4. 误差分析与工程验证：不止于看图

代码跑通了，得到了辨识模型，但这远远不是终点。如何定量评价辨识效果？除了肉眼对比阶跃响应曲线，我们还需要更严谨的指标和更深层的分析。

4.1 多维度性能评估

一个完整的评估应该包括时域、频域和定量指标。

% 评估辨识效果 function evaluate_identification(sys_original, sys_identified, T0) % 输入：原系统，辨识系统，采样周期 % A. 时域对比：阶跃响应 figure('Position', [50, 50, 1400, 500]); subplot(1,3,1); [y_orig, t_orig] = step(sys_original); [y_iden, t_iden] = step(sys_identified); plot(t_orig, y_orig, 'b-', 'LineWidth', 1.5); hold on; plot(t_iden, y_iden, 'r--', 'LineWidth', 1.5); grid on; xlabel('时间 (秒)'); ylabel('幅度'); title('阶跃响应对比'); legend('原系统', '辨识系统', 'Location', 'best'); % B. 频域对比：伯德图 subplot(1,3,2); bode(sys_original, 'b-', sys_identified, 'r--'); grid on; legend('原系统', '辨识系统', 'Location', 'best'); title('伯德图对比'); % C. 定量误差计算 subplot(1,3,3); % 在相同时间点上比较 t_eval = 0:0.01:20; % 评估时间点 y_orig_eval = step(sys_original, t_eval); y_iden_eval = step(sys_identified, t_eval); error = y_iden_eval - y_orig_eval; % 计算多种误差指标 MSE = mean(error.^2); % 均方误差 RMSE = sqrt(MSE); % 均方根误差 MAE = mean(abs(error));% 平均绝对误差 % 拟合优度 R-squared SS_res = sum(error.^2); SS_tot = sum((y_orig_eval - mean(y_orig_eval)).^2); R2 = 1 - (SS_res / SS_tot); % 绘制误差曲线 plot(t_eval, error, 'k-', 'LineWidth', 1); grid on; xlabel('时间 (秒)'); ylabel('误差'); title(sprintf('误差曲线 (MSE=%.2e, R^2=%.4f)', MSE, R2)); % 在命令行输出详细指标 fprintf('\n========== 辨识性能评估 ==========\n'); fprintf('采样周期 T0 = %.2f 秒\n', T0); fprintf('均方误差 (MSE): %.4e\n', MSE); fprintf('均方根误差 (RMSE): %.4f\n', RMSE); fprintf('平均绝对误差 (MAE): %.4f\n', MAE); fprintf('拟合优度 (R-squared): %.4f\n', R2); if R2 > 0.95 fprintf('评价: 辨识效果优秀。\n'); elseif R2 > 0.85 fprintf('评价: 辨识效果良好。\n'); elseif R2 > 0.70 fprintf('评价: 辨识效果一般。\n'); else fprintf('评价: 辨识效果较差，需检查采样周期或数据质量。\n'); end end

4.2 采样周期影响的系统性研究

现在，让我们系统地研究采样周期T0对辨识精度的影响。我们将测试一组T0值，并记录关键指标。

% 系统化研究采样周期的影响 T0_test = [0.05, 0.1, 0.2, 0.3, 0.5, 0.8, 1.0, 1.5]; n = 3; % 已知系统阶次 results = table('Size', [length(T0_test), 5], ... 'VariableTypes', {'double', 'double', 'double', 'double', 'double'}, ... 'VariableNames', {'T0', 'MSE', 'RMSE', 'MAE', 'R2'}); for idx = 1:length(T0_test) T0 = T0_test(idx); [sys_id_cont, ~] = hankel_identify(sys_continuous, T0, n); % 计算误差指标 t_eval = 0:0.01:15; y_orig = step(sys_continuous, t_eval); y_iden = step(sys_id_cont, t_eval); error = y_iden - y_orig; MSE = mean(error.^2); RMSE = sqrt(MSE); MAE = mean(abs(error)); SS_res = sum(error.^2); SS_tot = sum((y_orig - mean(y_orig)).^2); R2 = 1 - (SS_res / SS_tot); results(idx, :) = {T0, MSE, RMSE, MAE, R2}; end % 绘制误差指标随T0变化曲线 figure('Position', [100, 100, 1000, 400]); subplot(1,2,1); semilogy(results.T0, results.MSE, 'o-', 'LineWidth', 2, 'MarkerSize', 8); grid on; xlabel('采样周期 T0 (秒)'); ylabel('MSE (对数坐标)'); title('均方误差 (MSE) vs 采样周期'); subplot(1,2,2); plot(results.T0, results.R2, 's-', 'LineWidth', 2, 'MarkerSize', 8); grid on; xlabel('采样周期 T0 (秒)'); ylabel('拟合优度 R^2'); title('拟合优度 vs 采样周期'); ylim([0, 1.05]);

运行这段分析代码，你会得到类似下表的量化结果，以及更直观的趋势图：

从数据和图表中可以清晰地看到，随着T0增大，MSE呈指数上升趋势，R²则快速下降。在T0=0.2附近是一个性能拐点，之后辨识质量急剧恶化。这背后的原因主要有两点：

1. 混叠效应：采样频率过低，无法捕获系统的高频动态，导致信息永久丢失。
2. 双线性变换的畸变：d2c使用的Tustin变换在频率较高时（接近奈奎斯特频率）存在非线性畸变，T0越大，这个畸变影响的范围越向低频扩展。

4.3 应对噪声的鲁棒性改进

上述流程基于“干净”的仿真脉冲响应。但真实世界的数据总是伴有噪声。为了提升算法的工程实用性，我们必须考虑噪声的影响。一个常见的改进是在构造Hankel矩阵前，对脉冲响应数据进行预处理，例如使用滑动平均滤波。

% 为脉冲响应添加高斯白噪声，并测试滤波效果 noise_level = 0.02; % 噪声标准差为信号最大幅值的2% g_clean = impulse_response; g_noisy = g_clean + noise_level * max(abs(g_clean)) * randn(size(g_clean)); % 简单的滑动平均滤波 window_size = 5; b = (1/window_size)*ones(1, window_size); a = 1; g_filtered = filter(b, a, g_noisy); % 补偿滤波引入的相位延迟 g_filtered = g_filtered(ceil(window_size/2):end); figure; subplot(3,1,1); plot(g_clean, 'b-', 'LineWidth', 1.5); title('干净的脉冲响应'); grid on; subplot(3,1,2); plot(g_noisy, 'r-'); title('带噪声的脉冲响应'); grid on; subplot(3,1,3); plot(g_filtered, 'g-', 'LineWidth', 1.5); title('滤波后的脉冲响应'); grid on; xlabel('采样点');

使用滤波后的数据g_filtered进行Hankel矩阵构造和参数求解，通常能显著提升在噪声环境下的辨识鲁棒性。此外，也可以考虑使用更高级的总体最小二乘法来求解H * a = h方程，该方法能同时考虑系数矩阵H和观测向量h中的误差，比普通最小二乘法更抗噪。

5. 从仿真到实战：扩展应用与注意事项

掌握了基础的三阶系统辨识流程后，我们可以将其思路扩展到更广泛的场景。

场景一：高阶或未知阶次系统对于阶次未知的系统，前面提到的奇异值分解法是确定n的关键。在实际代码中，可以设定一个阈值（例如，最大奇异值的1%或0.1%），将大于该阈值的奇异值个数作为系统阶次的估计。

场景二：多输入多输出系统MIMO系统的Hankel矩阵辨识要复杂得多，需要构建分块Hankel矩阵，并利用其进行状态空间实现（如Ho-Kalman算法或子空间辨识）。MATLAB的系统辨识工具箱中的n4sid命令就是基于这类子空间方法。

场景三：从实验数据直接辨识如果你拥有的是系统的输入输出时间序列数据，而非脉冲响应，则需要先估计脉冲响应。这可以通过计算输入输出的互相关函数，或使用诸如频域估计等方法来实现。

重要提示：Hankel矩阵法对数据的“持续激励性”有要求。如果输入信号不能充分激发系统的所有模态，构造出的Hankel矩阵可能会秩亏，导致辨识失败。在实践中，使用伪随机二进制序列作为输入激励信号，是获取高质量脉冲响应估计的常用手段。

最后，分享一个我在调试中常遇到的坑：当采样周期非常小时，脉冲响应的初始几个采样点值可能极其接近零，导致Hankel矩阵H的条件数非常大（病态）。此时直接求逆或求解线性方程会引入巨大数值误差。解决方法一是适当增加采样周期，二是在求解时使用正则化技术，例如Tikhonov正则化，在(H'*H + lambda*I) \ (H'*h)中引入一个小的正则化参数lambda，以改善数值稳定性。