项目_‘华为杯’数模研赛复盘_fault_diagnosis_of_bearing_for_high_speed_第一问

34 阅读11分钟

目录


前言

碎碎念:知知为不知 不知为不知。

这是复盘华为杯E题系列,由于每一问都有很多细节值得再思考记忆,所以写的比较长,便一问写一篇。


提示:以下是本篇文章正文内容

赛题

赛题link:

​​​​​​https://www.kdocs.cn/l/cu3nOrmzuECG

一、第一问

1、问题背景与数据理解

分析数据树

核心任务就是从源域数据中选择部分数据进行特征提取

即从原始振动信号中提取能够有效表征轴承故障状态的多维特征向量。这个任务看似简单,但实际上涉及信号处理、特征工程、数据组织等多个环节。在开始编码之前,我首先对数据集的结构进行了分析,发现数据组织遵循了轴承故障诊断领域的标准命名规范(这种命名规范来源于凯斯西储大学(CWRU)轴承数据集),这为后续的数据加载和标签提取提供了一个文件处理的代码思路。  

编辑

这个数据组织超级麻烦。我要先收集所有.mat文件。

然后用scipy解析.mat文件。再依次对每个文件进行预处理:动态查找以_DE_time结尾的key,展平.mat文件里的多维数据作为signal,然后找rpm这个key的值。 然后对signal做 Z-score归一化。

简单说就是找到每个.mat文件里DE_time,signal键值对和rpm键值对

def load_and_preprocess(file_path):
    """
    加载和预处理.mat文件
    
    Parameters:
    -----------
    file_path : str
        .mat文件路径
    
    Returns:
    --------
    normalized_signal : np.array
        归一化后的信号(Z-score标准化)
    rpm : float
        转速(转/分钟)
    """
    # 步骤1:加载.mat文件(MATLAB格式)
    mat_data = scipy.io.loadmat(file_path)
    
    # 步骤2:动态查找数据键(核心逻辑)
    # 【为什么动态查找?】
    # 不同文件的数据键名可能不同,例如:
    # - B007_0_DE_time(故障数据)
    # - X097_DE_time(正常数据)
    # - 其他可能的命名格式
    
    data_key = None
    
    # 策略1:优先查找以_DE_time结尾的键(最常见格式)
    for key in mat_data.keys():
        if key.endswith('_DE_time'):
            data_key = key
            break

    # 如果找不到,抛出异常
    if data_key is None:
        raise ValueError(f"无法在文件 {file_path} 中找到数据键")
    
    # 步骤3:提取信号数据并展平为1D数组
    # .mat文件中的数据可能是多维数组,需要展平为一维
    signal = mat_data[data_key].flatten()

    # 步骤4:尝试从.mat文件中读取RPM
    if 'RPM' in mat_data:
        rpm_value = mat_data['RPM']
        if isinstance(rpm_value, np.ndarray):
            rpm = float(rpm_value.flatten()[0])  # 如果是数组,取第一个值
        else:
            rpm = float(rpm_value)

    
    # 步骤5:Z-score归一化(标准化)
    # 公式:(x - mean) / std
    # 目的:消除不同样本间的幅值差异,使特征具有可比性
    normalized_signal = (signal - np.mean(signal)) / np.std(signal)
    
    return normalized_signal, rpm

#不同的信号幅值存在差异,如果不归一化,幅值大的信号会在特征提取中占主导。
#避免某些特征因为幅值过大而主导模型训练过程。

数据筛选与平衡

首先,12kHz采样率的选择是基于奈奎斯特定理的考虑。对于轴承故障诊断,我们主要关心的频率范围通常在0-6kHz之间(这是12kHz采样率的奈奎斯特频率)。轴承的故障特征频率通常在几百Hz到几千Hz之间,12kHz的采样率已经足够捕捉这些频率成分。更高的采样率(如48kHz)虽然能够提供更多的频率信息,但也会带来数据量的急剧增加和计算成本的上升,经过严肃调研发现, 对于大多数故障诊断任务来说,这些额外的频率信息可能并不必要。

其次,DE(驱动端)数据的选择是基于信号质量的考虑。因为,驱动端的信号通常包含更丰富的故障信息,信号强度也更强。在实际的轴承故障诊断系统中,驱动端往往是首选的监测位置。使用DE数据能够保证特征提取的质量,为后续的分类任务提供更好的基础。

只使用12kHz_DE_data还有一个重要原因:与目标域数据的一致性,方便后面搞迁移学习。目标域数据的采样率是32kHz,虽然与源域的12kHz不同,但通过特征提取,我们可以将这种采样率的差异转化为特征空间的差异,符合迁移学习的场景。如果源域使用了多种采样率的数据,反而会增加数据分布的复杂性,不利于迁移学习。

Normal数据从48kHz_Normal_data目录中读取,主要是是因为正常状态的样本相对较少,需要从所有可用的数据源中收集,以保证数据集的平衡性,数据平衡就是为了避免模型偏向某些类别。

2、计算故障频率

这一步操作主要是服务于后续频域特征提取。因为在提取频域特征之前,需要计算轴承的理论故障特征频率。这些频率是基于轴承的几何参数和运行转速,通过理论公式计算得出的,它们代表了不同类型故障在频域中的特征位置。对于SKF 6205深沟球轴承,其关键参数包括:滚动体数量N_BALLS=9,滚动体直径D_BALL=0.3126英寸,节圆直径D_PITCH=1.537英寸,接触角ALPHA_CONTACT=0度(深沟球轴承的接触角通常为0度)。

​编辑

​编辑

​编辑

外圈故障频率(BPFO, Ball Pass Frequency Outer)的计算公式为:f_BPFO = (N_BALLS/2) × f_r × (1 - (D_BALL/D_PITCH) × cos(α)),其中f_r是转频(rpm/60)。这个公式的物理意义是:当外圈存在故障点时,每个滚动体经过故障点时会产生一次冲击,因此故障频率等于滚动体经过外圈故障点的频率。内圈故障频率(BPFI, Ball Pass Frequency Inner)的计算公式为:f_BPFI = (N_BALLS/2) × f_r × (1 + (D_BALL/D_PITCH) × cos(α)),其物理意义类似,但内圈是旋转的,所以公式中有一个加号。滚动体故障频率(BSF, Ball Spin Frequency)的计算公式为:f_BSF = (D_PITCH/(2×D_BALL)) × f_r × (1 - ((D_BALL/D_PITCH) × cos(α))²),这表示滚动体自身的旋转频率。

def calculate_fault_frequencies(rpm):
    """
    计算理论故障特征频率
    
    【故障频率的物理意义】
    轴承故障会产生周期性的冲击信号,这些冲击的频率与故障类型相关:
    - BPFO (Ball Pass Frequency Outer): 外圈故障频率
    - BPFI (Ball Pass Frequency Inner): 内圈故障频率
    - BSF (Ball Spin Frequency): 滚动体故障频率
    
    【计算公式来源】
    这些公式基于轴承几何学和运动学理论,是轴承故障诊断的经典公式
    
    Parameters:
    -----------
    rpm : float
        转速(转/分钟)
    
    Returns:
    --------
    fault_freqs : dict
        包含各种故障特征频率的字典
    """
    # 转频(旋转频率):rpm转换为Hz
    fr = rpm / 60.0
    
    # 接触角转换为弧度(用于计算)
    alpha_rad = np.deg2rad(ALPHA_CONTACT)
    
    # 外圈故障频率 (BPFO - Ball Pass Frequency Outer)
    # 公式:f_bpfo = (N/2) * fr * (1 - (d/D) * cos(α))
    # 物理意义:当外圈有故障时,滚动体每次经过故障点会产生冲击
    f_bpfo = (N_BALLS / 2) * fr * (1 - (D_BALL / D_PITCH) * np.cos(alpha_rad))
    
    # 内圈故障频率 (BPFI - Ball Pass Frequency Inner)
    # 公式:f_bpfi = (N/2) * fr * (1 + (d/D) * cos(α))
    # 物理意义:当内圈有故障时,滚动体每次经过故障点会产生冲击
    f_bpfi = (N_BALLS / 2) * fr * (1 + (D_BALL / D_PITCH) * np.cos(alpha_rad))
    
    # 滚动体故障频率 (BSF - Ball Spin Frequency)
    # 公式:f_bsf = (D/(2*d)) * fr * (1 - ((d/D) * cos(α))^2)
    # 物理意义:当滚动体有故障时,故障点周期性接触内外圈
    f_bsf = (D_PITCH / (2 * D_BALL)) * fr * (1 - ((D_BALL / D_PITCH) * np.cos(alpha_rad)) ** 2)
    
    return {
        'fr': fr,      # 转频
        'BPFO': f_bpfo,  # 外圈故障频率
        'BPFI': f_bpfi,  # 内圈故障频率
        'BSF': f_bsf     # 滚动体故障频率
    }

这一步操作主要是服务于后续频域特征提取。 由于实际信号中可能存在频率漂移(由于转速波动、轴承磨损等因素)

3、特征工程

时域

经过严肃调研发现,在很多文献中,内圈故障通常导致较高的峰度和脉冲因子,而外圈故障可能更多地影响RMS值。大多数轴承故障诊断系统都采用以下几种时域特征:

RMS(均方根值)是最基本的能量特征,它反映了信号的平均能量水平。在轴承故障诊断中,RMS值会随着故障的严重程度而增加,因为故障会产生额外的振动能量。峰度(Kurtosis)是衡量信号分布尖锐程度的统计量,正常轴承的振动信号通常接近高斯分布(峰度约为3),而故障信号由于存在周期性冲击,其峰度值会显著增大,有时甚至超过10。偏度(Skewness)衡量信号分布的对称性,故障信号由于冲击的存在,往往表现出非对称性,偏度值会偏离0。峰值(Peak)是信号的最大绝对值,它直接反映了冲击的强度。

峭度因子(Clearance Factor)的计算公式为 CF = Peak / (mean(√|x|))²,这是一个无量纲的特征,能够反映信号的冲击特性,对早期故障特别敏感。脉冲因子(Impulse Factor)的计算公式为 IF = Peak / mean(|x|),它衡量峰值相对于平均幅值的倍数,故障信号的这个值通常较大。形状因子(Shape Factor)的计算公式为 SF = RMS / mean(|x|),它反映了信号的形状特征,正常信号和故障信号在这个特征上通常有明显差异。

频域

经过严肃调研发现, 频域分析是轴承故障诊断的核心方法,因为故障会在特定的频率位置产生特征信号。然而,直接对原始信号进行FFT分析往往效果不佳,因为故障信号通常是调制信号,故障频率信息被调制在高频载波上。因此,我采用了包络谱分析(Envelope Spectrum Analysis)的方法,这在轴承故障诊断中算比较经典的了。

包络谱分析的第一步是计算信号的解析信号(Analytical Signal),这是通过Hilbert变换实现的。Hilbert变换能够将实信号转换为复信号,其实部是原信号,虚部是原信号的Hilbert变换。解析信号的模长(即包络)能够提取出信号的幅值调制信息,这正是我们需要的故障特征。计算包络后,我对包络信号进行FFT变换,得到包络谱。在包络谱中,故障频率及其谐波会以明显的峰值形式出现,这使得故障诊断变得直观和可靠。在提取故障频率幅值时,我不仅提取了基频(1倍频)的幅值,还提取了2倍频和3倍频的幅值,并取其中的最大值作为该故障频率的特征值。这样做的原因是,在实际信号中,故障频率的基频可能被噪声淹没,但其谐波成分往往更加明显。通过提取谐波的最大值,能够提高特征提取的鲁棒性。对于BPFO、BPFI、BSF和转频(fr),我都采用了这种处理方式。

除了故障频率幅值,我还计算了频谱质心(Spectral Centroid),这是频谱能量分布的加权平均频率,能够反映信号的主要频率成分位置。频谱质心是一个全局特征,它不依赖于特定的故障频率,而是从整体上描述信号的频率特性,这对于区分不同类型的故障也有一定帮助。

时频域

经过严肃调研,选择了小波包变换,短时傅里叶变换(STFT),经验模态分解(EMD)

4、特征融合合并数据

在for循环里得到文件名filename,对应的故障类,rpm字典再用字典的update函数合并两个字典再用列表list.append字典。 最后将列表转化为np.DataFrame().to_csv 导出特征工程的表格

简单说就是把这一行(一个文件的数据)塞进大列表里。循环结束后,这个列表通常会被转换成 Pandas DataFrame(表格),就像 Excel 一样,每一行是一个样本,每一列是一个特征。

5、总结

第一问涉及了信号处理、特征工程、数据组织等多个方面的知识。通过这个任务,我深入理解了轴承故障诊断的特征提取流程,掌握了时域、频域、时频域特征提取的方法,也学会了如何根据工程实际需求选择合适的数据和特征(选小而精的 不然复杂的数据对后面的模型学习不利)。


放个热乎乎的奖状证 感觉往年的好看。希望能给我补
没拿国一的原因竟是因为一个学校只有1个国一,我们被和谐了 。。。

编辑