机器学习系列(12)| 连续小波变换

职场   2024-11-05 12:05   四川  
傅里叶变换和小波变换都是信号处理领域中的重要工具,但它们在时频域分析的方式上存在明显差异。傅里叶变换的核心功能是压缩信号的时间维度,将信号从时域映射到频域,使得信号可以用不同频率的正弦波进行表示。换句话说就是,傅里叶变换能让我们识别出构成信号的频率成分信息。

用小时候背过的古诗《题西林壁》来作类比就很合适,在上图中,一个连续函数在时间域上的曲线看起来像连绵不绝的山脉;而从频率域看,其实是不同幅值、频率的正弦波组合,看起来像一座座高耸的山峰。对应诗句所说的:横看成岭侧成峰,远近高低各不同。发现没有?傅里叶变换其实就是用数学的方式将时域函数在基坐标上做了一个简单的投影(乘以正弦基函数)。

在向量空间中,点乘(或内积)可以用来度量两个向量在同一方向上的投影。当两个向量越接近平行时,点乘的结果越大,表示两个向量越相似。我们可以将点乘理解为一个信号a与另一个信号b的投影。

对于时间域的连续函数f(t),可以看成是一个连续无限维度的向量a,正弦基函数就是向量b。这样就可以通过用点乘的方法将信号从时域转化为频域,来分析信号的频率成分了。可世界并不仅仅是我们眼中的实部,所以为了从更高维空间来思考世界,傅里叶变换就可以看作是时域函f(t)与复指数函数的内积,从而把实数域扩展到带虚数的复数域。这样一来,傅里叶变换的表达公式就可以写成:

其中,ω是角频率,每个可以看作是这个基坐标系中的一个基向量。积分过程实际上是在测量f(t) 在每个频率ω上的投影强度。从公式也能看出,积分区间是从时间轴的负无穷到正无穷,无论每个频率成分出现在什么时间,对积分结果的影响都是一样的。

在分析处理信号时,仅仅知道信号内包含哪些频率成分往往是不够的,有时候我们更关心这些成分在时间轴上的分布情况。为了解决这一问题,短时傅里叶变换(Short-Time Fourier Transform, STFT)应运而生。它通过对信号施加一个时间窗(平移因子)来实现局部频率分析。例如,在处理一个包含10000个数据点的信号时,傅里叶变换会一次性处理所有数据点,从而损失了时间的信息。而短时傅里叶变换可以将信号分割成1000份,每部分包含10个数据点,并分别进行傅里叶变换,从而获得每个部分的频率信息。当窗口变得足够小,它就能反映出瞬时的频率信息。STFT的数学表达式定义为:

其中,f(τ) 是原始信号,g(τ-t) 是一个以t为中心的窗口函数(比如高斯窗),用于在时间t 附近对信号进行截取。
尽管这种方法能够提供瞬时频率信息,但它也存在一个很明显的缺陷:窗口大小是固定的,这限制了其在不同频率下的应用效果。当频率很高时,窗长显得太宽;当频率很低时,窗长又不够。

为了解决STFT窗口大小固定的问题,小波变换引入了一个尺度因子,从而在时间和频率上提供了更灵活的分析窗口(对于高频成分,使用窄的时间窗口来捕捉快速变化的细节;对于低频成分,使用宽的时间窗口来捕捉平缓变化的整体轮廓)。Wavelet 源自于法语“小波浪”。可以从两个维度来理解:一是在时域中能量有限,二是在频域中表现为带通滤波器的特性。
小波变换的核心思想是将信号与一种称为母小波的函数进行卷积,通过缩放和平移得到不同的分辨率。连续小波变换(CWT)的数学表达式为:

其中:
 是对母小波ψ(t)进行缩放和平移得到的小波基函数,除以根号a的绝对值是为了做归一化处理。

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)
scales参数:表示尺度序列,是一个数组,用于调节母小波的宽度。这里我设置了总尺度数量 totalscal 并利用小波的中心频率 fc(在pywt库中,可以通过pywt.central_frequency来获取特定小波基函数的中心频率)计算出 cparam,尺度序列 scales 即可按照频率反比关系生成。这种方法会根据中心频率自动调整尺度范围,保证了尺度和频率的合理匹配,使得时频分析结果能够覆盖整个频率范围。

wavename参数:母小波类型,这里我选了morl。

1.0 / sampling_rate:小波变换的采样周期,它是采样率的倒数,表示每个样本之间的时间间隔。

执行 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个频率成分随时间是等长传播的,没有发现明显的瞬时突变或振荡。总体来说,该信号的频率结构较为稳定。

用一个3D视图来看,可以更直观一些:

下面再来个信号异常的例子,如下图所示,一个平稳信号在某个时刻出现了高频的异常振荡。下面我们尝试用CWT来捕捉这些异常振动发生的具体时间和频率:

经过小波变换后,从3D时频图中我们可以得到这些信息:在整个时间范围内,低频信号是连续且平稳的,在约0.4~0.5秒时间段出现一个短暂的高频振动峰值。


有了上述理论基础,在航空维修领域,我们应该如何利用它来辅助我们判断故障呢?举个例子,某飞机机组反映升降率波动异常,ECAM无警告,地面自检BMC1/2无故障信息。们知道,导致升降率波动的可能原因有:引气压力波动、组件活门间歇性卡滞、外流活门间歇性卡阻或CPC输出的增压控制信号异常。把QAR数据译码并可视化,从时域图中,我们是无法直接判断出哪个部件是故障源的。下图只截取了故障发生期间的数据:

虽然客舱升降率VSCB在波动,但其它的参数曲线相对比较平滑,无论是引气压力、组件流量还是外流活门,各参数值都并未呈现来回波动的情况。于是我尝试了使用快速傅里叶变换转换到频域后,观察各参数的频率分量的数量。从下图中,能明显看到PACK_FLOW_R2的频率分量数量有别于PACK_FLOW_R1。但由于缺乏频率分量出现的时间,仅仅靠它来作为故障判断依据还是不够科学的。

针对傅里叶变换无法同时给出时域频域信息的缺陷,我再次尝试了小波变换。下图是小波变换后的伪色彩图,颜色从蓝色代表小波系数的强度,即信号的频率成分的能量或幅值红色表示能量较高,蓝色区域则表示能量较低。y轴表示频率,从1Hz到4Hz。x轴表示时间,从0到1800秒

应该怎么解析这几张图呢?下面我们着重从压力信号和流量信号去分析:

1、压力信号 (PRECOOL_PRESS1 和 PRECOOL_PRESS2):

PRECOOL_PRESS1频率范围从1 到 4 Hz,整体能量水平较低,主要表现为深蓝色。在接近 1750 秒时有轻微的能量增加(颜色变浅),可能表示在该时刻信号发生了一些变化,但整体影响较小。

PRECOOL_PRESS2:与 PRECOOL_PRESS1 类似,能量大部分时间较低,接近 1000 到 1750 秒时有较明显的能量增加(颜色变浅,接近绿色),可能意味着该时刻有突发的信号活动或变化。

2、流量信号(PACK_FLOW_R1和PACK_FLOW_R2):

PACK_FLOW_R1:信号的能量在整个时间段内都较低,只有在接近 1750 秒时出现微弱的能量增加。

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后排除,印证了我们的分析。


总结一下,小波变换因具有出色的时频局部性和多尺度分析能力,使其更适合捕获不同时间和频率尺度上的信号局部特征。其适应非线性和非平稳信号、捕捉动态变化、可以和神经网络的结合等特性,使其成为故障检测领域中一个强大的工具。在航空维修领域,我们可以通过将多变量的时序数据进行小波变换,比较各个参数在某个时间段内特定频率成分的增强或减弱或幅值的变化来辅助我们判断故障源。如果想故障隔离得更精准,还需要对数据进行预处理,包括去噪、标准化等。更关键的是一定要选择合适的小波基函数和尺度,才能捕捉到参数之间的微小变化。


感谢大家对机务论坛的支持,关注机务论坛,倾听机务心声!航企优秀的方面必定宣传,不足的地方也必须指出,让领导们重视问题,解决问题,营造更好的机务维修环境。

征稿:
所见所闻,个人感悟,个人成长历程及感人故事。
特别征稿:我师傅的故事!
同时,征集劳动仲裁案例,分享案例,让更多的小伙伴能了解劳动纠纷的解决方式,通过劳动仲裁维护自己的合法权益。




评论区留言,同意的点赞
扫码添加小编微信
匿名爆料

民航机务论坛微信公众平台
改名为:机务论坛
发布最新行业动向 深入解读政策法规
开辟维修工程专栏 交流飞机排故经验
分享前沿技术应用 预测职业发展前景
行业大咖讲经布道 业界专家授业解惑
致力打造一流的民航机务朋友圈----机务论坛
关注机务论坛,倾听机务心声!
投稿邮箱:duanwei0615@163.com


机务论坛
民航机务论坛改名为:机务论坛 发布最新行业动向 深入解读政策法规 开辟维修工程专栏 交流飞机排故经验 分享前沿技术应用 预测职业发展前景 行业大咖讲经布道 业界专家授业解惑 致力打造一流的民航机务朋友圈----机务论坛
 最新文章