算法特点
1. 多频带谱幅融合创新
改进传统单频带健康指数的局限性，实现全频段谱幅信息的自适应加权融合
通过优化权重分配，实现关键故障频带的自动识别与增强
2. 多目标协同优化
将可分性、单调性、稀疏性三个关键性能指标统一到优化框架中
提出信噪比优化模型，平衡健康指数在不同退化阶段的动态特性
3. 凸优化理论应用
将健康指数构建转化为凸优化问题，保证全局最优解
引入松弛变量处理单调性约束，提高模型可行性和鲁棒性
4. 可解释性设计
优化权重可视化，揭示各频带对健康状态的贡献度
故障特征与物理频带对应，实现数据驱动→物理理解的闭环
5. 自适应退化阶段划分
结合K-means聚类与能量特征，动态识别不同退化阶段
避免传统等分方法的假设偏差，更符合实际退化过程
6. 稀疏性增强
针对故障特征的稀疏特性，引入L1正则化提取关键故障频率
降低冗余信息干扰，提高健康指数对早期微弱故障的敏感性
算法步骤
1. 数据预处理
加载振动信号数据，进行分段加窗处理
计算频谱和包络谱，获取频域特征矩阵
2. 退化阶段识别
动态划分正常、早期退化、严重退化等阶段
建立阶段标签用于可分性约束
3. 可分性优化建模
构建正常与退化阶段的最大分离度目标函数
引入类间距最大化约束，增强健康指数的判别能力
4. 单调性约束嵌入
添加相邻时间点健康指数的单调递增约束
使用松弛变量处理非严格单调情况，提高模型可行性
5. 稀疏性权重学习
针对故障特征的稀疏分布，采用L1正则化提取关键频率
通过交叉验证优化正则化参数，平衡稀疏性与表达能力
6. 多目标协同求解
融合可分性、单调性、稀疏性构建统一优化问题
使用凸优化求解器（ECOS/SCS）获取全局最优权重
7. 健康指数计算与评估
应用优化权重计算广义健康指数
评估单调性、可分性、趋势性等关键性能指标
8. 结果可视化与解释
绘制健康指数退化曲线，标注各阶段边界
可视化权重分布，解释关键故障频带的贡献
import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.fft import fft, fftfreq from scipy.signal import hilbert import cvxpy as cp from sklearn.cluster import KMeans from sklearn.linear_model import LassoCV, Ridge import warnings warnings.filterwarnings('ignore') # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False class GeneralizedHealthIndex: """ 基于谱幅融合的广义健康指数构建与优化模型 """ def __init__(self, fs=12000, segment_length=1024, overlap=0.5): """ 初始化参数 参数: ---------- fs : int 采样频率 (Hz) segment_length : int 分段长度 overlap : float 重叠率 (0-1) """ self.fs = fs self.segment_length = segment_length self.overlap = overlap self.weights = None self.health_indices = None def load_data(self, data_dir='data'): """ 加载西储大学轴承数据 参数: ---------- data_dir : str 数据目录路径 返回: ---------- data_dict : dict 加载的数据字典 labels : dict 数据标签 """ import os data_files = ['98raw.txt', '106raw.txt', '119raw.txt', '131raw.txt'] labels = { '98raw.txt': 'Normal', '106raw.txt': 'Ball Fault', '119raw.txt': 'Inner Race Fault', '131raw.txt': 'Outer Race Fault' } data_dict = {} for file in data_files: file_path = os.path.join(data_dir, file) if os.path.exists(file_path): # 读取数据 data = pd.read_csv(file_path, header=None, names=['vibration']) data_dict[file] = data['vibration'].values print(f"Loaded {file}: {labels[file]}, length: {len(data_dict[file])}") else: print(f"Warning: {file_path} not found!") return data_dict, labels def segment_signal(self, signal, max_segments=None): """ 信号分段 参数: ---------- signal : numpy array 原始振动信号 max_segments : int, optional 最大分段数 返回: ---------- segments : list 分段后的信号列表 """ segments = [] step = int(self.segment_length * (1 - self.overlap)) max_start = len(signal) - self.segment_length + 1 if max_segments is not None: max_start = min(max_start, max_segments * step) for i in range(0, max_start, step): segment = signal[i:i + self.segment_length] segments.append(segment) return np.array(segments) def compute_spectrum(self, signal_segments): """ 计算信号的频谱 参数: ---------- signal_segments : numpy array 分段后的信号，形状为 (n_segments, segment_length) 返回: ---------- spectrum_magnitude : numpy array 频谱幅值，形状为 (n_segments, n_frequencies) frequencies : numpy array 频率数组 """ n_segments = signal_segments.shape[0] n_frequencies = self.segment_length // 2 + 1 spectrum_magnitude = np.zeros((n_segments, n_frequencies)) for i in range(n_segments): # 应用汉宁窗 windowed_signal = signal_segments[i] * np.hanning(self.segment_length) # 计算FFT fft_result = fft(windowed_signal) # 取幅值（正频率部分） magnitude = np.abs(fft_result[:n_frequencies]) spectrum_magnitude[i] = magnitude # 频率数组 frequencies = fftfreq(self.segment_length, 1/self.fs)[:n_frequencies] return spectrum_magnitude, frequencies def compute_envelope_spectrum(self, signal_segments): """ 计算信号的包络谱 参数: ---------- signal_segments : numpy array 分段后的信号 返回: ---------- envelope_spectrum : numpy array 包络谱幅值 envelope_frequencies : numpy array 包络谱频率 """ n_segments = signal_segments.shape[0] n_frequencies = self.segment_length // 2 + 1 envelope_spectrum = np.zeros((n_segments, n_frequencies)) for i in range(n_segments): # 计算解析信号（希尔伯特变换） analytic_signal = hilbert(signal_segments[i]) # 计算包络 envelope = np.abs(analytic_signal) # 对包络信号进行FFT windowed_envelope = envelope * np.hanning(self.segment_length) fft_result = fft(windowed_envelope) magnitude = np.abs(fft_result[:n_frequencies]) envelope_spectrum[i] = magnitude # 频率数组 envelope_frequencies = fftfreq(self.segment_length, 1/self.fs)[:n_frequencies] return envelope_spectrum, envelope_frequencies担任《Mechanical System and Signal Processing》《中国电机工程学报》《宇航学报》《控制与决策》等期刊审稿专家，擅长领域：信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测。基于多频带谱幅自适应加权融合与凸优化理论的装备全生命周期可解释退化评估广义健康指数智能构建算法（Python）
参考文章：
智瞳·工业装备健康透视—基于谱幅自适应融合的退化评估智能算法（Python） - 哥廷根数学学派的文章
https://zhuanlan.zhihu.com/p/2000114112613599005
