包络解调在故障诊断中的应用-广义检波解调原理

前言

笔者在做故障诊断的时候，一直在用希尔伯特包络解调，但好像也一直没有深究包络解调技术背后的原理，近来有空正好整理一下相关的内容。

包络解调

包络解调（Envelope Demodulation） 是一种从调幅（AM）信号中提取调制信号的信号处理技术，其核心目的是分离出隐藏在载波信号中的低频调制信息（即包络），以便分析周期性或瞬态冲击特征。
在机械系统中，故障（如轴承损伤、齿轮断齿）通常引发周期性冲击，这些冲击被系统的固有频率（载波频率）调制，形成调幅信号：
$$x(t)=A(1+Bcos(2\pi f_{n}t))cos(2\pi f_{z}t)$$
其中$f_n$是低频调制信号（故障引起的周期性冲击），$f_z$是载波频率（如机械结构的共振频率）。那么包络解调的目的就是要从$x(t)$中提取$f_n$，进而识别机械结构健康状态。
目前常见的包络解调方法也有很多了，如绝对值分析、检波滤波、平方解调、希尔伯特解调、谱峭度解调、平方包络谱等，其中绝对值分析、检波滤波、平方解调三者都属于广义检波解调技术，它们的解调原理都是对信号取绝对值或平方，然后再进行低通滤波，有很大的相似之处，下面就分别展开叙述。

绝对值分析

原理

对上述$x(t)$取绝对值，可得：
$$z(t)=|x(t)|=A|1+Bcos(2\pi f_{n}t)||cos(2\pi f_{z}t)|$$
傅里叶技术的变换公式如下：
$$f(t)=\frac{1}{2}{a}_0+\sum _{n=1}^{\infty }[a_{n}cos(nwt)+b_{n}sin(nwt)]$$
由于$|cos(2\pi f_{z}t)|$是偶函数，因此所有的$b_n$项都是0，只需要计算$a_0$和$a_n$，展开后只有余弦项，而取绝对值则会将周期减半，因此基频由$f_z$变为$2f_z$，将$|cos(2\pi f_{z}t)|$展开为傅里叶级数的时候应为如下形式：
$$|cos(2\pi f_{z}t)|=\frac{a_0}{2}+\sum _{n=1}^{\infty }{a}_{n}cos(2\pi n2f_{z}t)$$

最终可得：
$$|cos(2\pi f_{z}t)|=\frac{2}{\pi }+\frac{4}{3\pi }cos(2\pi 2f_{z}t)-\frac{4}{15\pi }cos(2\pi 4f_{z}t)+...$$
同样地将$|1+Bcos(2\pi f_{n}t)|$也可展开为傅里叶级数，不过此处需要分情况讨论，如果$B<1$，该式恒大于0：
$$|1+Bcos(2\pi f_{n}t)|=1+Bcos(2\pi f_{n}t)$$
如果$B>1$，由于是偶函数，因此只有余弦项：
$$|1+Bcos(2\pi f_{n}t)|=\frac{a_0}{2}+\sum _{n=1}^{\infty }{a}_{n}cos(2\pi n{f}_{n}t)$$
展开式中含有$f_n$的各谐波分量。
依据卷积定理，时域相乘等于频域卷积，因此取绝对值之后的频谱中只含有调制频率$f_n$及其高次谐波，以及以载波频率$f_z$的偶次谐波为中心，调制频率及其高次谐波为边频带的频率成分。对其进行低通滤波，只要选取合适的截止频率，即可得到包含有调制频率及其高次谐波的解调谱。

Python示例

下面我们通过一个例子来演示一下绝对值分析的流程。首先我们创建一个调制频率0.5Hz的$f_n$以及一个载波频率25Hz的$f_z$，将二者相乘得到调幅信号$x_t$：

import numpy as np
import matplotlib.pyplot as plt


plt.rcParams['font.sans-serif'] = ['SimHei']  # 步骤一（替换sans-serif字体）
plt.rcParams['axes.unicode_minus'] = False   # 步骤二（解决坐标轴负数的负号显示问题）


A = 3.0
B = 2.0
f_n = 0.5
f_z = 25
f_s = 8192
theta = np.linspace(0, 8, f_s * 8)
x_n = np.cos(2 * np.pi * f_n * theta)
x_z = np.cos(2 * np.pi * f_z * theta)
x_t = A * (1 + B * x_n) * x_z

# 绘图对比
plt.figure(figsize=(10, 6))
plt.subplot(311)
plt.plot(theta, x_n, label='0.5Hz', linewidth=2)
plt.ylabel('幅值')
plt.legend()
plt.subplot(312)
plt.plot(theta, x_z, '--', label='25Hz')
plt.ylabel('幅值')
plt.legend()
plt.subplot(313)
plt.plot(theta, x_t, '--', label='调幅信号')
plt.xlabel('时间')
plt.ylabel('幅值')
plt.suptitle('测试信号')
plt.grid(True)
plt.show()

对调幅信号使用傅里叶变换，我们可以看到其频域上有载波信号为中心，调制频率为边频带的频率成分。

# 傅里叶变换
spectrum = np.fft.fft(x_t) / len(x_t)
spectrum = np.abs(spectrum)[:len(spectrum) // 2] * 2  # 单边频谱

# 计算频率轴
frequencies = np.fft.fftfreq(len(x_t), d=1 / f_s)[:len(spectrum)]
# 限制频率范围为 0-30 Hz
mask = (frequencies > 0) & (frequencies <= 30)
frequencies = frequencies[mask]
spectrum = spectrum[mask]

plt.figure(figsize=(10, 6))
plt.plot(frequencies, spectrum)
plt.xlabel('频率')
plt.ylabel('幅值')
plt.title('频谱')
plt.grid(True)
plt.show()

我们使用绝对值分析后，再进行傅里叶变换，可以看到频谱中只含有调制频率$f_n$及其高次谐波，以及以载波频率$f_z$的偶次谐波为中心，调制频率及其高次谐波为边频带的频率成分。

# 绝对值解调
z_t = np.abs(x_t)
# 傅里叶变换
spectrum = np.fft.fft(z_t) / len(z_t)
spectrum = np.abs(spectrum)[:len(spectrum) // 2] * 2  # 单边频谱

# 计算频率轴
frequencies = np.fft.fftfreq(len(z_t), d=1 / f_s)[:len(spectrum)]
# 限制频率范围为 0-200 Hz
mask = (frequencies > 0) & (frequencies <= 200)
frequencies = frequencies[mask]
spectrum = spectrum[mask]

plt.figure(figsize=(10, 6))
plt.plot(frequencies, spectrum)
plt.xlabel('频率')
plt.ylabel('幅值')
plt.title('频谱')
plt.grid(True)
plt.show()

如果选择合适的截止频率进行低通滤波，则可以将调制信号及其谐波解调出来。

检波滤波

原理

检波滤波解调的原理就是将调幅信号$x_t$的负半周置为零：
$$\begin{array}{rl}{z}_a(t)& =0.5|x(t)|+0.5x(t)\\ & =0.5A|1+Bcos(2\pi f_{n}t)||cos(2\pi f_{z}t)|+0.5A(1+Bcos(2\pi f_{n}t))cos(2\pi f_{z}t)\\ & =0.5z(t)+0.5A(1+Bcos(2\pi f_{n}t))cos(2\pi f_{z}t)\end{array}$$
对上式进行频谱分析，将会发现频率成分只含有调制频率$f_n$及其高次谐波（幅值为绝对值解调分析方法的一半），以及以载波频率$f_z$和其偶次谐波为中心，调制频率及其高次谐波为边频带的频率成分。

Python示例

我们同样用上面的测试信号来测试。

# 检波解调
z_a = 0.5 * np.abs(x_t) + 0.5 * x_t
# 傅里叶变换
spectrum = np.fft.fft(z_a) / len(z_a)
spectrum = np.abs(spectrum)[:len(spectrum) // 2] * 2  # 单边频谱

# 计算频率轴
frequencies = np.fft.fftfreq(len(z_a), d=1 / 8192)[:len(spectrum)]
# 限制频率范围为 0-30 Hz
mask = (frequencies > 0) & (frequencies <= 200)
frequencies = frequencies[mask]
spectrum = spectrum[mask]

plt.figure(figsize=(10, 6))
plt.plot(frequencies, spectrum)
plt.xlabel('频率')
plt.ylabel('幅值')
plt.title('频谱')
plt.grid(True)
plt.show()

从频谱上可以看到，和我们理论的分析一致，后续选择合适的截止频率进行低通滤波，则可以将调制信号及其谐波解调出来。

平方解调

原理

平方解调的原理就是对$x(t)$取平方：
$$\begin{array}{rl}{z}_s(t)& =A^2(1+Bcos(2\pi f_{n}t){)}^2(cos(2\pi f_{z}t){)}^2\\ & =A^2(1+2Bcos(2\pi f_{n}t)+B^{2}cos(2\pi f_{n}t{)}^2)cos(2\pi f_{z}t{)}^2\\ & =A^2(1+\frac{B^2}{2}+2Bcos(2\pi f_{n}t)+\frac{B^2}{2}cos(2\pi 2f_{n}t))(\frac{1}{2}+\frac{1}{2}cos(2\pi 2f_{z}t))\\ & =(\frac{1}{2}{A}^2+\frac{1}{4}{A}^2{B}^2)+A^{2}Bcos(2\pi f_{n}t)+\frac{1}{4}{A}^2{B}^{2}cos(2\pi 2f_{n}t)+(\frac{1}{2}{A}^2+\frac{1}{4}{A}^2{B}^2)cos(2\pi 2f_{z}t)+A^{2}Bcos(2\pi f_{n}t)cos(2\pi 2f_{z}t)+\frac{1}{4}{A}^2{B}^{2}cos(2\pi 2f_{n}t)cos(2\pi 2f_{z}t)\end{array}$$
依据卷积定理，对$z_s(t)$进行傅里叶变换后得到的频谱中包含直流分量、一倍调制频率$f_n$、二倍调制频率$2f_n$、二倍载波频率成分$2f_z$和以$2f_z$为中心，一倍频和二倍调制频率为边带的频率成分。

Python示例

我们同样用上面的测试信号来测试。

# 平方解调
z_s = x_t ** 2
# 傅里叶变换
spectrum = np.fft.fft(z_s) / len(z_s)
spectrum = np.abs(spectrum)[:len(spectrum) // 2] * 2  # 单边频谱

# 计算频率轴
frequencies = np.fft.fftfreq(len(z_s), d=1 / 8192)[:len(spectrum)]
# 限制频率范围为 0-30 Hz
mask = (frequencies > 0) & (frequencies <= 200)
frequencies = frequencies[mask]
spectrum = spectrum[mask]

plt.figure(figsize=(10, 6))
plt.plot(frequencies, spectrum)
plt.xlabel('频率')
plt.ylabel('幅值')
plt.title('频谱')
plt.grid(True)
plt.show()

从频谱中可以看到，和理论的频率成分一致。

广义检波解调的局限

实际上，对于广义检波解调方法来说，仍然存在两个方面的局限性：
1、将不包括故障调制信息的两时域相加信号，也以其频率之差作为解调信号解出；
2、由于处理的是数字信号，因而在解调过程中可能会产生特有的解调混频效应；
对于第一个问题，我们可以生成两个测试信号进行模拟，一个信号$f_1$为24Hz，一个信号$f_2$为25Hz：

# 两个信号相加
f_n = 24
f_z = 25
f_s = 8192
theta = np.linspace(0, 8, f_s * 8)
x_n = np.cos(2 * np.pi * f_n * theta)
x_z = np.cos(2 * np.pi * f_z * theta)
x_p = x_n + x_z
z_t = np.abs(x_p)
# 傅里叶变换
spectrum = np.fft.fft(z_t) / len(z_t)
spectrum = np.abs(spectrum)[:len(spectrum) // 2] * 2  # 单边频谱

# 计算频率轴
frequencies = np.fft.fftfreq(len(z_t), d=1 / f_s)[:len(spectrum)]
# 限制频率范围为 0-200 Hz
mask = (frequencies > 0) & (frequencies <= 100)
frequencies = frequencies[mask]
spectrum = spectrum[mask]

plt.figure(figsize=(10, 6))
plt.plot(frequencies, spectrum)
plt.xlabel('频率')
plt.ylabel('幅值')
plt.title('频谱')
plt.grid(True)
plt.show()

如果两个信号的频率都比较大，而且接近，那么低通滤波后会在频谱上显示$f_2-f_1$的基频及其高次谐波。

总结

本篇文章介绍了广义检波解调的原理，包括绝对值解调、检波滤波解调、平方解调三种方法，虽然笔者在实际工作中很少用，但也可以适当了解，作为自己的知识储备。