【AI引导共振定位+物理校验】 基于无监督CAE和谱峭度自动寻优的滚动轴承故障诊断方法

def signal_to_spectrogram(signal):    # 短时傅里叶变换，仅取幅值    f, t, Zxx = stft(signal, fs=FS, nperseg=128)    spec = np.abs(Zxx)                # 丢弃相位    spec = np.log1p(spec)             # 对数压缩动态范围    spec = (spec - spec.min()) / (spec.max() - spec.min() + 1e-8)    spec = cv2.resize(spec, (IMG_SIZE, IMG_SIZE))  # 缩放到固定尺寸    return spec

class UltraTightAutoencoder(nn.Module):    def __init__(self, latent_dim=128):        super().__init__()        self.encoder_cnn = nn.Sequential(            nn.Conv2d(3, 16, 3, 2, 1), nn.ReLU(),            nn.Conv2d(16, 32, 3, 2, 1), nn.ReLU(),            nn.Conv2d(32, 64, 3, 2, 1), nn.ReLU(),            nn.Conv2d(64, 128, 3, 2, 1), nn.ReLU(),        )        self.flatten = nn.Flatten()        self.encoder_linear = nn.Linear(128*14*14, latent_dim)        self.decoder_linear = nn.Linear(latent_dim, 128*14*14)        self.unflatten = nn.Unflatten(1, (128, 14, 14))        self.decoder_cnn = ...  # 转置卷积还原至原尺寸    def forward(self, x):        ...

def kurtogram_find_frequency_band(residual_image, img_size, max_freq_of_image, bandwidth, min_lowcut):    # 若为彩色图，先转为灰度（此处取均值）    gray_residual = np.mean(residual_image, axis=2) if len(residual_image.shape)==3 else residual_image    # 沿频率轴（行）计算峭度    kurtosis_profile = kurtosis(gray_residual, axis=1, fisher=True)    peak_y_pixel = np.argmax(kurtosis_profile)          # 峭度最大行索引    freq_per_pixel = max_freq_of_image / img_size    peak_frequency = peak_y_pixel * freq_per_pixel    lowcut = max(min_lowcut, peak_frequency - bandwidth)    highcut = min(max_freq_of_image, peak_frequency + bandwidth)    return lowcut, highcut

def butter_bandpass_filter(data, lowcut, highcut, fs, order=4):    nyq = 0.5*fs    b, a = signal.butter(order, [low/nyq, high/nyq], btype='band')    return signal.filtfilt(b, a, data)   # 零相位滤波def perform_envelope_analysis(filtered_signal, fs):    analytic = signal.hilbert(filtered_signal)    envelope = np.abs(analytic)    envelope -= np.mean(envelope)    freqs = np.fft.rfftfreq(len(envelope), d=1/fs)    fft_vals = np.abs(np.fft.rfft(envelope)) * 2 / len(envelope)    return freqs, fft_vals, envelope

def auto_diagnose(freqs, fft_amps, bpfo, bpfi, threshold, tolerance):    mask = (freqs > 40) & (freqs <= 500)    valid_freqs, valid_amps = freqs[mask], fft_amps[mask]    f_max = valid_freqs[np.argmax(valid_amps)]    amp_max = valid_amps.max()    if amp_max < threshold:        return "NORMAL", "幅值低于阈值，无故障"    # 计算与BPFO/BPFI的相对误差    err_bpfo = abs(f_max - bpfo)/bpfo    err_bpfi = abs(f_max - bpfi)/bpfi    if min(err_bpfo, err_bpfi) <= tolerance:        # 匹配，输出对应缺陷类型        return "FAULT", "匹配成功"    else:        return "NORMAL", "幅值虽高但物理频率不匹配，归为正常"

# ============================ 2. 定义AI模型结构（卷积自编码器） ============================class UltraTightAutoencoder(nn.Module):    """    轻量级沙漏型卷积自编码器 (CAE)    作用：仅用健康状态谱图训练，作为"正常基线"的压缩与重建模型。         对故障样本会产生较大的重建残差，从而凸显异常能量。    """    def __init__(self, latent_dim=128):        super().__init__()        # ---------- 编码器部分：4层下采样，提取深层特征 ----------        self.encoder_cnn = nn.Sequential(            nn.Conv2d(3, 16, kernel_size=3, stride=2, padding=1),  # 224->112            nn.ReLU(),            nn.Conv2d(16, 32, kernel_size=3, stride=2, padding=1), # 112->56            nn.ReLU(),            nn.Conv2d(32, 64, kernel_size=3, stride=2, padding=1), # 56->28            nn.ReLU(),            nn.Conv2d(64, 128, kernel_size=3, stride=2, padding=1),# 28->14            nn.ReLU(),        )        self.flatten = nn.Flatten()        self.encoder_linear = nn.Linear(128 * 14 * 14, latent_dim) # 压缩到隐向量        # ---------- 解码器部分：4层上采样，还原图像 ----------        self.decoder_linear = nn.Linear(latent_dim, 128 * 14 * 14)        self.relu = nn.ReLU()        self.unflatten = nn.Unflatten(1, (128, 14, 14))        self.decoder_cnn = nn.Sequential(            nn.ConvTranspose2d(128, 64, kernel_size=3, stride=2, padding=1, output_padding=1), # 14->28            nn.ReLU(),            nn.ConvTranspose2d(64, 32, kernel_size=3, stride=2, padding=1, output_padding=1),  # 28->56            nn.ReLU(),            nn.ConvTranspose2d(32, 16, kernel_size=3, stride=2, padding=1, output_padding=1),  # 56->112            nn.ReLU(),            nn.ConvTranspose2d(16, 3, kernel_size=3, stride=2, padding=1, output_padding=1),   # 112->224            nn.Sigmoid(),  # 输出归一化到 [0,1]        )    def forward(self, x):        x = self.encoder_cnn(x)        x = self.flatten(x)        x = self.encoder_linear(x)          # 得到压缩后的潜在表征        x = self.decoder_linear(x)        x = self.relu(x)        x = self.unflatten(x)        x = self.decoder_cnn(x)             # 输出重建谱图        return x# ============================ 3. 信号预处理与AI推理工具 ============================def load_and_preprocess_image(image_path, image_size, device):    """    从磁盘加载谱图图片，并转为模型输入的Tensor格式    """    import cv2    img = cv2.imread(image_path)    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)          # 转为RGB通道    img_resized = cv2.resize(img, (image_size, image_size))    img_normalized = img_resized.astype(np.float32) / 255.0    # 转为 PyTorch 格式: [C, H, W] 并增加 batch 维度    tensor_input = torch.tensor(np.transpose(img_normalized, (2, 0, 1))).unsqueeze(0).to(device)    return tensor_inputdef run_autoencoder_inference(model, tensor_input):    """    执行CAE前向推理，返回重建图与正残差图（ReLU处理）    残差 R = max(0, 输入 - 重建) ，放大模型无法解释的异常区域    """    with torch.no_grad():        tensor_output = model(tensor_input)        tensor_residual = torch.relu(tensor_input - tensor_output)  # 仅保留正向误差    return tensor_output, tensor_residual# ============================ 4. 核心：AI引导共振频带自动定位（谱峭度法） ============================def kurtogram_find_frequency_band(residual_image, img_size, max_freq_of_image, bandwidth, min_lowcut):    """    基于残差图的谱峭度自动寻找最优共振解调频带。    原理：残差图中能量越集中的区域，峭度越高，该频带越可能包含故障冲击成分。    """    # 若残差图为彩色（3通道），转为灰度图（取各通道均值）    if len(residual_image.shape) == 3:        gray_residual = np.mean(residual_image, axis=2)    else:        gray_residual = residual_image    # 沿频率轴（纵轴/行方向）计算峭度，反映各频率行上能量分布的尖锐程度    kurtosis_profile = kurtosis(gray_residual, axis=1, fisher=True)    kurtosis_profile = np.nan_to_num(kurtosis_profile)  # 防止出现NaN    # 找到峭度最大的行索引（像素坐标）    peak_y_pixel = int(np.argmax(kurtosis_profile))    # 将像素坐标映射到真实物理频率（Hz）    freq_per_pixel = max_freq_of_image / img_size    peak_frequency = peak_y_pixel * freq_per_pixel    # 以峰值频率为中心，结合设定带宽生成带通范围，并做边界截断    lowcut = max(min_lowcut, peak_frequency - bandwidth)    highcut = min(max_freq_of_image, peak_frequency + bandwidth)    print(f"\n[AI引导定位] 谱峭度峰值中心频率: {peak_frequency:.1f} Hz")    print(f"[AI引导定位] 推荐带通滤波范围: {lowcut:.1f} - {highcut:.1f} Hz\n")    return lowcut, highcut# ============================ 5. 原始一维信号加载与物理滤波解调 ============================def load_cwru_signal(mat_path, fs):    """    从CWRU数据集格式的.mat文件中提取驱动端加速度信号（DE_time）    并截取前 fs 个点（即1秒数据）用于分析    """    data = sio.loadmat(mat_path)    # 自动查找包含 "DE_time" 的键名（不同版本命名略有差异）    key = next((k for k in data.keys() if "DE_time" in k), None)    if key is None:        raise ValueError(f"在 {mat_path} 中未找到 DE_time 信号")    sig = data[key].flatten()[:fs]  # 取前 fs 点（等效1秒时长）    sig = sig - np.mean(sig)        # 去除直流分量    return sigdef butter_bandpass_filter(data, lowcut, highcut, fs, order=4):    """    零相位数字巴特沃斯带通滤波器    用于提取AI指定的最佳共振频带内的冲击成分    """    nyq = 0.5 * fs    low = lowcut / nyq    high = highcut / nyq    b, a = signal.butter(order, [low, high], btype='band')    return signal.filtfilt(b, a, data)  # filtfilt实现零相位延迟def perform_envelope_analysis(filtered_signal, fs):    """    希尔伯特包络解调 + 包络频谱分析    从高频共振信号中解调出低频故障脉冲序列    """    analytic_signal = signal.hilbert(filtered_signal)    amplitude_envelope = np.abs(analytic_signal)    amplitude_envelope = amplitude_envelope - np.mean(amplitude_envelope)  # 去除均值    n = len(amplitude_envelope)    frequencies = np.fft.rfftfreq(n, d=1/fs)    fft_values = np.abs(np.fft.rfft(amplitude_envelope)) * 2 / n    return frequencies, fft_values, amplitude_envelope# ============================ 6. 物理确定性校验与故障决策 ============================def auto_diagnose(freqs, fft_amps, bpfo, bpfi, threshold, tolerance):    """    物理力学校验：将包络谱主峰频率与理论BPFO/BPFI进行对比，    只有既满足幅值阈值、又满足频率容差（±2%）时才判定为故障，否则视为正常波动。    """    # 限制分析频段（通常为40~500Hz，避开低频转频干扰）    mask = (freqs > 40.0) & (freqs <= 500.0)    valid_freqs = freqs[mask]    valid_amps = fft_amps[mask]    if len(valid_amps) == 0:        return "NORMAL", "有效频段无数据"    max_idx = valid_amps.argmax()    f_max = valid_freqs[max_idx]    amp_max = valid_amps[max_idx]    print("\n========== 物理力学校验 (Physics Verification) ==========")    # 第一层判断：峰值幅值是否达到预警门槛    if amp_max < threshold:        status = "NORMAL"        report = (f"峰值频率 {f_max:.2f} Hz 幅值 ({amp_max:.4f}) 低于阈值，"                  f"结论：无机械故障迹象（正常）")        print(report)        return status, report    # 第二层判断：计算主峰与理论频率的相对误差    err_bpfo = abs(f_max - bpfo) / bpfo    err_bpfi = abs(f_max - bpfi) / bpfi    # 选取匹配度更高（误差更小）的那一个作为最佳匹配    if err_bpfo <= err_bpfi:        best_error, theo_freq, defect_name = err_bpfo, bpfo, "BPFO (外圈故障)"        match_key = "BPFO"    else:        best_error, theo_freq, defect_name = err_bpfi, bpfi, "BPFI (内圈故障)"        match_key = "BPFI"    # 若误差在容忍范围内（如2%），确认故障；否则归为正常（拒绝虚假报警）    if best_error <= tolerance:        status = "FAULT_CONFIRMED"        report = (f"峰值 {f_max:.2f} Hz 与理论 {match_key} ({theo_freq:.2f} Hz) "                  f"匹配度达 {(1-best_error)*100:.1f}%，结论：{defect_name}")    else:        status = "NORMAL"        report = (f"峰值 {f_max:.2f} Hz 虽然幅值高，但与理论频率误差 "                  f"{best_error*100:.1f}% 超出容差 {tolerance*100:.0f}%，"                  f"结论：视为正常波动，拒绝虚警")    print(report)    return status, report# ============================ 7. 端到端一键运行主流程 ============================def run_end_to_end_pipeline(config_path):    """    全流程执行函数：AI引导定位 -> 物理滤波解调 -> 力学公式校验    输入：配置文件路径（包含模型路径、信号参数、诊断阈值等）    """    # ---------- 加载配置 ----------    with open(config_path, "r", encoding="utf-8") as f:        cfg = yaml.safe_load(f)    # 核心超参数提取    MODEL_PATH = cfg["paths"]["model_path"]    TEST_IMAGE_PATH = cfg["paths"]["test_image_path"]   # 由测试信号生成的谱图    TEST_MAT_PATH = cfg["paths"]["test_file_mat"]       # 对应的原始一维信号(.mat)    FS = cfg["signal"]["fs"]                 # 采样率 48kHz    Fr = cfg["signal"]["fr"]                 # 转频（Hz）    BPFO = cfg["signal"]["bpfo_multiplier"] * Fr    BPFI = cfg["signal"]["bpfi_multiplier"] * Fr    IMG_SIZE = cfg["image"]["image_size"]              # 谱图尺寸（如224）    MAX_FREQ = cfg["image"]["max_freq_of_image"]       # 谱图对应的最大频率（如3000Hz）    BANDWIDTH = cfg["filter"]["bandwidth"]             # 共振带半宽度（如150Hz）    MIN_LOWCUT = cfg["filter"]["min_lowcut"]           # 最小截止频率（如500Hz）    THRESHOLD = cfg["diagnosis"]["threshold"]          # 幅值门限    TOLERANCE = cfg["diagnosis"]["tolerance"]          # 频率容差（如0.02）    LATENT_DIM = cfg["model"]["latent_dim"]            # 隐向量维度    DEVICE = "cuda" if torch.cuda.is_available() else "cpu"    # ---------- 步骤1：加载CAE模型 ----------    print("[1] 加载预训练卷积自编码器...")    model = UltraTightAutoencoder(latent_dim=LATENT_DIM).to(DEVICE)    state_dict = torch.load(MODEL_PATH, map_location=DEVICE, weights_only=True)    model.load_state_dict(state_dict)    model.eval()  # 切换为推理模式    # ---------- 步骤2：谱图推理，计算残差并自动定位共振频带 ----------    print("[2] 执行AI谱图推理，计算残差...")    tensor_in = load_and_preprocess_image(TEST_IMAGE_PATH, IMG_SIZE, DEVICE)    _, tensor_res = run_autoencoder_inference(model, tensor_in)    # 将残差张量转为numpy图像格式（HWC），供谱峭度分析使用    residual_img = tensor_res.squeeze().cpu().numpy().transpose(1, 2, 0)    lowcut, highcut = kurtogram_find_frequency_band(        residual_img, IMG_SIZE, MAX_FREQ, BANDWIDTH, MIN_LOWCUT    )    # ---------- 步骤3：回溯原始一维信号 -> 带通滤波 -> 包络解调 ----------    print("[3] 回溯原始信号，执行物理带通滤波与希尔伯特解调...")    raw_signal = load_cwru_signal(TEST_MAT_PATH, FS)    filtered_signal = butter_bandpass_filter(raw_signal, lowcut, highcut, FS)    freqs, fft_amps, envelope = perform_envelope_analysis(filtered_signal, FS)    # ---------- 步骤4：物理力学校验，输出最终诊断 ----------    print("[4] 执行物理频率匹配与决策...")    final_status, final_report = auto_diagnose(        freqs, fft_amps, BPFO, BPFI, THRESHOLD, TOLERANCE    )    # ---------- 可视化展示 ----------    print("\n========== 最终诊断结果 ==========")    print(f"诊断状态: {final_status}")    print(f"详情报告: {final_report}")    # 绘制包络频谱（展示峰值与理论频率线）    plt.figure(figsize=(10, 4))    mask_plot = freqs <= 500    plt.plot(freqs[mask_plot], fft_amps[mask_plot], color='blue', label='包络频谱')    plt.axvline(BPFO, color='orange', linestyle='--', label=f'BPFO={BPFO:.1f}Hz')    plt.axvline(BPFI, color='red', linestyle='--', label=f'BPFI={BPFI:.1f}Hz')    plt.axvline(Fr, color='green', linestyle=':', label=f'Fr={Fr:.1f}Hz')    plt.xlabel('频率 (Hz)')    plt.ylabel('幅值')    plt.title('AI引导共振定位后包络解调频谱与物理频率比对')    plt.legend()    plt.grid(alpha=0.3)    plt.tight_layout()    plt.show()    return final_status, final_report