频率响应系统辨识方法：从实测数据建模的完整示例

从实测数据建模：频率响应系统辨识的实战指南

你有没有遇到过这样的场景？手头有一个“黑盒子”系统——可能是电机驱动器、电源环路，也可能是某个声学装置。你知道它有输入和输出，但内部结构复杂甚至完全未知。你想设计控制器，却苦于没有数学模型；想优化性能，又缺乏对带宽、谐振点的准确判断。

这时候，频率响应系统辨识就是你的破局利器。

不同于依赖先验知识的理论建模，这种方法只关心“输入-输出”关系，通过实验测量直接构建系统的动态模型。整个过程就像给系统做一次“频域体检”，不仅能看清它的增益与相位变化，还能提取出可用于仿真的传递函数。

本文将以一个完整的工程实例为线索，带你走完从信号激励到模型拟合的全过程，深入剖析每个环节的技术细节与常见陷阱。无论你是控制工程师、嵌入式开发者，还是从事振动分析或电源设计的技术人员，都能从中获得可复用的方法论。

如何设计有效的激励信号？

一切始于输入信号的选择。别小看这一步——很多辨识失败的根本原因，就藏在一开始施加的那个测试信号里。

为什么不能随便用白噪声？

我见过不少初学者直接用白噪声去“轰”被测系统，结果要么信噪比太低，要么执行机构因峰值过大而饱和。问题在于：理想很丰满，现实很骨感。

虽然白噪声理论上能激发所有频率，但它具有很高的峰值因子（Crest Factor），即瞬时峰值远高于有效值。这对功率级设备是致命的：比如一个平均功率合适的信号，可能在某几个时刻突然冲高，导致MOSFET过压或电机堵转。

所以，真正实用的激励信号必须满足几个硬性条件：

常见激励类型怎么选？

✅ 推荐首选：Chirp 信号（线性/对数扫频）

import numpy as np # 生成对数chirp信号，0.1Hz ~ 5kHz，持续2秒 fs = 20000 # 采样率 t = np.linspace(0, 2, int(fs*2), endpoint=False) u = signal.chirp(t, f0=0.1, f1=5000, t1=2, method='logarithmic')

优势：
- 能量集中在当前扫描频率附近，信噪比高
- 实现简单，适合现场快速测试
- 对非强非线性系统表现良好

⚠️ 注意事项：避免扫得太快！确保每个频率都有足够时间让系统达到稳态响应，否则会低估相位滞后。

✅ 高精度实验室推荐：Multisine 信号

这类信号由多个离散正弦波叠加而成，频率成分可控，且可通过优化相位降低峰值因子。

# 构造multisine：仅包含特定频率点 frequencies = np.logspace(np.log10(0.1), np.log10(5000), 100) U = np.zeros(len(t), dtype=complex) for f in frequencies: omega = 2 * np.pi * f # 添加随机相位并限制峰值 phase = np.random.uniform(0, 2*np.pi) U += np.exp(1j*(omega*t + phase)) u_multisine = np.real(U) # 取实部作为激励

这种信号特别适合需要重复验证的高精度建模任务，比如滤波器原型验证或材料阻抗测量。

❌ 慎用：纯白噪声 / PRBS

除非你在做在线辨识或系统监控，否则不要轻易使用原始白噪声。如果非要使用，务必进行幅值限幅处理或结合窗函数平滑。

数据采集不是“录下来就行”

很多人以为只要把输入输出信号存下来就万事大吉了，其实这才是挑战的开始。

同步采集：差一毫秒都不行

假设你用两块独立的数据采集卡分别记录输入和输出，即使各自时钟都很准，微小的时间漂移也会在高频段引入严重的相位误差。例如，在1kHz处，1ms的延迟就会造成360°的相位偏移！

✅正确做法：
- 使用同一块多通道DAQ设备
- 所有通道共享同一个时钟源和触发信号
- 若必须异步采集，需事后通过互相关对齐时间轴

抗混叠滤波不可省

奈奎斯特准则说“采样率至少是最高频率的两倍”，但这只是理论下限。实际中建议采样率至少为最大分析频率的5~10倍，并在前端加入模拟低通滤波器（抗混叠滤波器），防止高频干扰折叠进有用频段。

加窗处理：双刃剑

FFT要求信号是周期性的，但现实中大多数激励都不是整周期截断的。这就导致了频谱泄漏——能量“洒”到邻近频率上，影响估计精度。

解决办法是加窗，常用Hanning窗或Flat Top窗：

def preprocess(time, u_raw, y_raw): # 去除直流分量（趋势项） u_detrend = signal.detrend(u_raw) y_detrend = signal.detrend(y_raw) # 应用Hanning窗 window = np.hanning(len(u_detrend)) u_win = u_detrend * window y_win = y_detrend * window return time, u_win, y_win

⚠️但注意：加窗虽能减少泄漏，也会降低频率分辨率。如果你关心两个靠得很近的谐振峰，过度加窗会让它们“糊”在一起。因此，是否加窗、用什么窗，要根据具体需求权衡。

怎么算出可靠的频率响应？别再只用 Y/U 了！

最简单的想法是：对输入 $ u(t) $ 和输出 $ y(t) $ 做FFT，然后逐点相除得到 $ H(f) = Y(f)/U(f) $。听起来合理，但在真实世界中几乎一定会翻车。

为什么？因为总有噪声存在。一旦输入或输出中含有噪声，简单的比值估计就会严重偏离真实值。

更稳健的做法：用功率谱估计 FRF

工业界广泛采用的是基于互功率谱与自功率谱的 H1 估计法：

$$
H_1(f) = \frac{G_{yu}(f)}{G_{uu}(f)}
$$

其中：
- $ G_{yu} $ 是输出与输入的互功率谱
- $ G_{uu} $ 是输入的自功率谱

这个方法假设输入无噪声、输出有噪声，这在大多数传感器测量场景中是成立的（毕竟你能精确控制输入信号）。相比直接相除，H1估计对方差更鲁棒，尤其适合信噪比较低的情况。

实现上我们通常用Welch 方法（分段平均）来进一步降低估计方差：

from scipy.signal import csd, welch # 参数设置 fs = 20000 # 采样率 nperseg = 4096 # 每段长度 overlap = 0.75 # 重叠比例（75%） # 计算自功率谱和互功率谱 f, P_uu = welch(u_windowed, fs=fs, nperseg=nperseg, noverlap=int(nperseg*overlap)) f, P_yy = welch(y_windowed, fs=fs, nperseg=nperseg, noverlap=int(nperseg*overlap)) f, P_yu = csd(y_windowed, u_windowed, fs=fs, nperseg=nperseg, noverlap=int(nperseg*overlap)) # H1估计 H1 = P_yu / P_uu

判断结果可信吗？看相干函数！

光有 $ H_1(f) $ 还不够，你还得知道哪些频率的数据是可靠的。这就是相干函数（Coherence）的作用：

$$
\gamma^2(f) = \frac{|G_{yu}(f)|^2}{G_{yy}(f) G_{uu}(f)}
$$

它的取值范围是 [0,1]：
- $ \gamma^2 \approx 1 $：说明该频率下输入与输出高度相关，FRF估计可信
- $ \gamma^2 < 0.8 $：可能存在噪声污染、非线性失真或外部干扰，结果应谨慎对待

coherence = np.abs(P_yu)**2 / (P_yy * P_uu) # 绘图查看 plt.figure() plt.semilogx(f, coherence) plt.xlabel('Frequency (Hz)') plt.ylabel('Coherence') plt.title('Coherence Function γ²(f)') plt.ylim(0, 1.1)

💡经验法则：如果某个频段的相干值偏低，先检查以下几点：
- 是否存在外部干扰源（如开关电源噪声）？
- 激励信号在该频段是否有足够能量？
- 系统是否进入了非线性区（如磁饱和）？

如何从频响数据拟合成可用的传递函数？

现在你有了 $ H_1(f) $，也就是一堆复数形式的频率响应点。下一步是要把它变成一个可以放进Simulink仿真、用于PID调参的传递函数模型。

模型结构怎么定？

常见的选择有三种：
1.传递函数形式：$ G(s) = \frac{b_m s^m + \cdots + b_0}{a_n s^n + \cdots + a_0} $
2.零极点增益模型：$ G(s) = k \frac{(s - z_1)\cdots(s - z_m)}{(s - p_1)\cdots(s - p_n)} $
3.状态空间模型：适用于多输入多输出系统

对于单输入单输出（SISO）系统，推荐从低阶传递函数开始尝试。

拟合策略：别盲目追求高阶！

新手常犯的错误是不断往上加阶数，直到曲线完美贴合。殊不知这很可能陷入了过拟合：模型记住了噪声，却失去了物理意义。

✅ 正确做法：
- 先观察Bode图的大致形状：有几个转折点？是否有明显谐振峰？
- 根据特征初步判断模型阶数。例如：
- 一阶惯性环节 → 2阶传递函数
- 含谐振的机械系统 → 至少4阶（二阶振荡环节）
- 使用AIC/BIC信息准则辅助选择最优阶数

实战代码：最小二乘拟合传递函数

import control as ct from scipy.optimize import least_squares def tf_response(params, s): """计算传递函数在s=jω处的响应""" num = np.polyval(params['num'], s) den = np.polyval(params['den'], s) return num / den def residual(params, w, H_meas): s = 1j * w H_pred = tf_response(params, s) err_real = np.real(H_meas - H_pred) err_imag = np.imag(H_meas - H_pred) return np.hstack([err_real, err_imag]) # 示例：拟合一个3阶分母、2阶分子的传递函数 initial_guess = { 'num': [1.0, 2.0, 1.0], # 分子系数 [b2, b1, b0] 'den': [1.0, 3.0, 3.0, 1.0] # 分母系数 [a3, a2, a1, a0] } result = least_squares( residual, x0=np.concatenate([initial_guess['num'], initial_guess['den']]), args=(2*np.pi*f[1:], H1[1:]), # 忽略DC点 method='lm' # Levenberg-Marquardt算法 ) # 重构传递函数 num_fit = result.x[:3] den_fit = result.x[3:] sys_identified = ct.TransferFunction(num_fit, den_fit, True)

💡 提示：初始值很关键！你可以先用手绘方式估算截止频率和阻尼比，或者用矢量拟合（Vector Fitting）工具自动提供初值。

工程应用中的那些“坑”与对策

坑1：系统明明是非线性的，还适用吗？

频率响应本质上基于线性时不变（LTI）假设。如果你的系统存在严重非线性（如死区、饱和、摩擦），那不同幅值下的FRF会不一样。

✅ 对策：
- 在多个幅值下重复测试，观察FRF一致性
- 选择线性度较好的工作区间进行辨识
- 或改用非线性辨识方法（如Hammerstein模型）

坑2：温度变了，模型就不准了？

确实如此。尤其是电机、电源等功率器件，温升会导致电感、电阻参数漂移。

✅ 对策：
- 在热稳定状态下测试
- 定期重新辨识，用于自适应控制或健康监测

坑3：如何用于控制器设计？

一旦你得到了可靠的传递函数模型，就可以直接用于：
-稳定性分析：计算增益裕度、相位裕度
-PID整定：使用Ziegler-Nichols或回路整形方法
-前馈设计：构造逆模型补偿动态滞后
-故障诊断：对比新旧FRF，检测参数偏移

写在最后：掌握这项技能意味着什么？

当你学会用chirp信号“唤醒”一个沉默的系统，看着屏幕上逐渐浮现的Bode图，那一刻你会感受到一种独特的掌控感。

你不再依赖厂商提供的模糊规格书，也不必盲调PID参数碰运气。你可以：
- 准确说出系统的带宽是多少Hz
- 指出哪个频率存在机械共振
- 设计陷波滤波器精准抑制振荡
- 评估控制方案的鲁棒性边界

这正是数据驱动建模的力量。

随着边缘计算和智能传感器的发展，这类技术正越来越多地集成到自动化测试平台中，成为预测性维护、数字孪生、自愈控制系统的核心组件。

如果你正在做控制、电源、音频或机电一体化相关的工作，强烈建议把这套方法纳入你的工具箱。它不会让你一夜成为专家，但会在每一次调试中默默为你节省数小时的试错时间。

如果你在实践中遇到了具体的难题——比如“为什么我的相干函数总是很低？”、“怎样处理含有延迟的系统？”——欢迎在评论区留言，我们一起探讨解决方案。