在向量空间中,点乘(或内积)可以用来度量两个向量在同一方向上的投影。当两个向量越接近平行时,点乘的结果越大,表示两个向量越相似。我们可以将点乘理解为一个信号a与另一个信号b的投影。
在分析处理信号时,仅仅知道信号内包含哪些频率成分往往是不够的,有时候我们更关心这些成分在时间轴上的分布情况。为了解决这一问题,短时傅里叶变换(Short-Time Fourier Transform, STFT)应运而生。它通过对信号施加一个时间窗(平移因子)来实现局部频率分析。例如,在处理一个包含10000个数据点的信号时,傅里叶变换会一次性处理所有数据点,从而损失了时间的信息。而短时傅里叶变换可以将信号分割成1000份,每部分包含10个数据点,并分别进行傅里叶变换,从而获得每个部分的频率信息。当窗口变得足够小,它就能反映出瞬时的频率信息。STFT的数学表达式定义为:
a是尺度因子,控制小波的伸缩(当a>1 时,小波被拉伸,捕捉低频信息;当 0<a<1 时,小波被压缩,捕捉高频信息);b是平移因子,控制小波在时间轴上的位置。
离散小波变换(DWT)我还没学,这里就先不讲了。对比STFT和CWT的公式,不难看出,CWT其实就是把短时傅里叶变换的正弦波基函数替换成了一个和时间、频率有关的二元函数ψ(t)。所以说,母小波才是小波变换的核心,我们常用的Morl大概是长这样子的:
为了加深理解,下面我们通过一段代码演示一下连续小波变换的应用:
首先,我构造了四个不同频率的正弦波信号:10 Hz、25 Hz、50 Hz和120 Hz,并将4个信号串联起来得到一个总信号signal。
import numpy as np
import matplotlib.pyplot as plt
import pywt
# 采样率和时间向量
fs = 500 # 采样率
duration = 1 # 总时长
time_per_signal = duration / 4 # 每个信号的时长
time = np.linspace(0, time_per_signal, int(fs * time_per_signal), endpoint=False)
# 第一个频率成分
signal1 = np.sin(2 * np.pi * 10 * time)
# 第二个频率成分
signal2 = np.sin(2 * np.pi * 25 * time)
# 第三个频率成分
signal3 = np.sin(2 * np.pi * 50 * time)
# 第四个频率成分
signal4 = np.sin(2 * np.pi * 120 * time)
# 合并四个信号
signal = np.concatenate((signal1, signal2, signal3, signal4))
接着,调用fft函数进行快速傅里叶计算,并绘制出FFT频谱图:
fft = np.fft.fft(signal)
fft_freq = np.fft.fftfreq(len(signal), d=1/fs)
# 只显示正频率部分
n = len(signal)
n_half = n // 2 # 只取前半部分
plt.plot(fft_freq[:n_half], np.abs(fft)[:n_half])
plt.xlabel('频率')
plt.ylabel('幅度')
plt.title('FFT 频谱图')
plt.show()
从频谱图中,我们能看到在横坐标10 Hz、25 Hz、50 Hz和120 Hz处有4个明显的峰值,代表着总信号有4个主要的频率成分。可由于每个频率分量是对应整个时间轴的,所以我们无法知道这几个频率的先后顺序是怎么样的。假如某一时刻发生了瞬时突变、脉冲或震荡,我们也无法定位出来的。
为了继续了解信号在时间和频率上的局部信息,我们先选定一个中心频率的母小波(不指定中心频率也行,若不指定,PyWavelets 会使用母小波的默认中心频率。例如,Morlet 小波的默认中心频率约为 0.8125 Hz),然后再用一系列在不同尺度和位置上的母小波对原始信号进行卷积。即可得到了不同时间位置、不同频率成分的相似性,这个相似性的大小就是小波系数。
在 pywt 库中,我们可以使用 pywt.cwt 函数来实现连续小波变换,通过设置母小波类型、尺度和其他参数来调整结果:
# 设置CWT参数
# 设置采样频率,这里与信号的采样率相同
sampling_rate = fs
# 尺度数量
totalscal = 128
# 小波基函数
wavename = 'morl'
# 小波函数中心频率
fc = pywt.central_frequency(wavename)
# 常数c
cparam = 2 * fc * totalscal
# 尺度序列
scales = cparam / np.arange(1, totalscal + 1)
# 进行CWT连续小波变换
coefficients, frequencies = pywt.cwt(signal, scales, wavename, 1.0 / sampling_rate)
# 小波系数矩阵绝对值
amp = abs(coefficients)
wavename参数:母小波类型,这里我选了morl。
执行 pywt.cwt 函数后,它返回两个输出:
1、小波系数矩阵(coefficients):一个二维数组,表示信号在不同尺度和时间位置上的小波系数。其大小为len(scales), len(signal)
,即行数等于尺度数,列数等于信号的采样点数。其绝对值可以理解为信号在对应尺度和时间点上的“局部能量”,大值意味着该时刻和尺度下的频率成分较强。在时频分析时,通过它我们就可以观察信号的时频特征。
2、频率系列(frequencies):一个数组,表示每个尺度对应的实际频率。长度为 len(scales)
,每个元素对应于 scales
中的一个尺度。提供了将尺度映射到频率的转换关系,这样我们就可以在时频图中标注出实际的频率,而不仅仅是尺度。
# 绘制时频图
plt.figure(figsize=(14,7))
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
# 根据采样频率sampling_period 生成时间轴 t
t = np.linspace(0, 1.0, fs, endpoint=False)
# 使用 np.meshgrid 函数生成时间 t 和频率 frequencies 的网格
T, F = np.meshgrid(t, frequencies)
# 绘制伪彩色图
plt.pcolormesh(T, F, amp, cmap='jet', shading='auto')
# 添加一个颜色条表示小波系数的幅度值
plt.colorbar(label='振幅')
plt.xlabel('时间')
plt.ylabel('频率')
plt.title('CWT时频图')
plt.show()
CWT时频图的颜色表示不同的能量,越接近红色,表示信号在该时间和频率上的强度越高,而蓝色表示强度较低。从下图中,我们能明显看到4个不同颜色在时间轴上的分布情况:最开始出现的是低频信号,越往后频率越高。4个频率成分随时间是等长传播的,没有发现明显的瞬时突变或振荡。总体来说,该信号的频率结构较为稳定。
有了上述理论基础,在航空维修领域,我们应该如何利用它来辅助我们判断故障呢?举个例子,某飞机机组反映升降率波动异常,ECAM无警告,地面自检BMC1/2无故障信息。我们知道,导致升降率波动的可能原因有:引气压力波动、组件活门间歇性卡滞、外流活门间歇性卡阻或CPC输出的增压控制信号异常。把QAR数据译码并可视化,从时域图中,我们是无法直接判断出哪个部件是故障源的。下图只截取了故障发生期间的数据:
虽然客舱升降率VSCB在波动,但其它的参数曲线相对比较平滑,无论是引气压力、组件流量还是外流活门,各参数值都并未呈现来回波动的情况。于是我尝试了使用快速傅里叶变换转换到频域后,观察各参数的频率分量的数量。从下图中,能明显看到PACK_FLOW_R2的频率分量数量有别于PACK_FLOW_R1。但由于缺乏频率分量出现的时间,仅仅靠它来作为故障判断依据还是不够科学的。
针对傅里叶变换无法同时给出时域频域信息的缺陷,我再次尝试了小波变换。下图是小波变换后的伪色彩图,颜色从蓝色到红色代表小波系数的强度,即信号的频率成分的能量或幅值。红色表示能量较高,蓝色区域则表示能量较低。y轴表示频率,从1Hz到4Hz。x轴表示时间,从0到1800秒。
应该怎么解析这几张图呢?下面我们着重从压力信号和流量信号去分析:
1、压力信号 (PRECOOL_PRESS1 和 PRECOOL_PRESS2):
PRECOOL_PRESS2:与 PRECOOL_PRESS1 类似,能量大部分时间较低,接近 1000 到 1750 秒时有较明显的能量增加(颜色变浅,接近绿色),可能意味着该时刻有突发的信号活动或变化。
2、流量信号(PACK_FLOW_R1和PACK_FLOW_R2):
PACK_FLOW_R2:能量变化比较显著。从 1000 秒到接近 1500 秒时,频率成分的能量有所增加,尤其是在低频段,颜色逐渐变亮。这表明在这一时间段内,PACK_FLOW_R2 信号的低频成分更活跃或有较高的能量。
通过上述时频图中得到的信息,我们可知:虽然PACK_FLOW_R2 和 PRECOOL_PRESS2 都在接近 1000 秒和 1750 秒时有异常能量增加(这两个时间点上的异常可能代表了故障的发生时间),但PRECOOL_PRESS2振幅的增加稍微晚了一些(在时频图中,某一时刻的能量突然增加通常表明该信号直接受到故障影响)。PACK_FLOW_R2 虽然高频段的能量增加不如 PRECOOL_PRESS2 明显,但它在 1000 秒后出现了更广泛的整体能量分布,PACK_FLOW_R2 的能量变化更加持久,延续到较长的时间段。而 PRECOOL_PRESS2 主要在高频段(3 Hz 以上)出现了较显著的能量增加,青绿色区域较为集中,但它的变化幅度相对集中在高频和较短时间窗口内。从能量变化幅度的时刻以及持续时间来看,PACK_FLOW_R2
可能是故障的源头,PRECOOL_PRESS2
可能是对流量故障的响应。综上,导致客舱升降率波动的故障源为右组件的FCV。故障最终也在更换右侧FCV后排除,印证了我们的分析。
总结一下,小波变换因具有出色的时频局部性和多尺度分析能力,使其更适合捕获不同时间和频率尺度上的信号局部特征。其适应非线性和非平稳信号、捕捉动态变化、可以和神经网络的结合等特性,使其成为故障检测领域中一个强大的工具。在航空维修领域,我们可以通过将多变量的时序数据进行小波变换,比较各个参数在某个时间段内特定频率成分的增强或减弱或幅值的变化来辅助我们判断故障源。如果想故障隔离得更精准,还需要对数据进行预处理,包括去噪、标准化等。更关键的是一定要选择合适的小波基函数和尺度,才能捕捉到参数之间的微小变化。