首先对原始轴承振动信号进行多尺度特征提取，包括时域的均方根值、峰值和峭度特征，时频域的小波能量特征，以及频域的频谱能量分布特征，形成高维特征空间。

接着构建轻量级物理信息神经网络架构，该网络不仅包含常规的深度学习组件，还创新性地引入了黎曼几何正则化机制，通过计算特征间的余弦相似度来保持数据流形结构，同时嵌入物理先验知识约束模型输出的单调性以确保退化过程符合物理规律。

在模型训练阶段，系统采用自适应权重平衡策略，让数据损失、物理约束损失和黎曼正则化损失在学习过程中自动调整比例，配合梯度裁剪和余弦退火学习率调度确保训练稳定性。

完成训练后，系统从当前监测时间点开始，采用逐步前向预测策略模拟轴承退化轨迹，当预测的特征值超过预设阈值时即判定为故障发生时间，从而计算出剩余使用寿命。最后，系统通过丰富的可视化工具展示完整的退化曲线、预测轨迹和RUL对比结果，并生成详细的性能评估报告，为工业决策提供数据支持。

算法步骤说明

数据预处理与特征提取

读取轴承振动信号数据，将连续的振动信号分割成固定长度的数据段，对每个数据段同步提取多种类型的特征指标，包括反映信号整体能量的均方根值、表征冲击成分的峰值幅度、描述分布形状的峭度统计量、通过小波变换获得的时频能量分布特征，以及通过傅里叶变换计算的频谱能量分区特征，形成能够全面表征轴承健康状态的多维特征向量。

模型架构设计与初始化

构建一个轻量级的神经网络结构，采用全连接层配合双曲正切激活函数，在网络中引入黎曼几何正则化模块，通过计算特征向量之间的余弦相似度来保持原始数据的流形结构，同时设计物理约束模块，确保模型预测的退化趋势满足单调非减的基本物理规律，所有约束权重均设为可学习参数以实现自适应平衡。

标准化处理与训练准备

对输入的时间序列数据和输出的特征值分别进行标准化处理，消除量纲影响并加速训练收敛，将故障开始时间之后到当前监测时间点之间的数据作为训练样本，时间点作为输入，对应的特征值作为输出标签。

多目标损失函数优化

定义包含三个核心组件的复合损失函数，包括衡量预测精度的数据拟合损失、保证物理合理性的约束损失和维护数据结构的黎曼正则化损失，采用自适应权重机制自动调整各项损失的重要性比例，使用带权重衰减的优化算法更新网络参数，实施梯度裁剪防止训练不稳定，应用余弦退火策略动态调整学习率。

训练过程监控与早停

在训练过程中持续监控验证集上的损失变化，当损失连续多个周期不再显著下降时提前终止训练，保存表现最佳的模型参数，记录训练过程中的各项损失变化趋势以备分析。

退化轨迹预测与阈值判断

从当前监测时间点开始，将时间坐标输入训练好的网络模型，递归预测未来时间点的特征值变化趋势，形成完整的退化轨迹预测，当预测的特征值首次超过预设的故障阈值时，记录该时间点为预测的故障发生时刻。

剩余使用寿命计算

将预测的故障发生时间减去当前监测时间，得到轴承的预测剩余使用寿命，与实际已知的剩余使用寿命进行对比，计算预测误差和准确率指标。

结果可视化与性能评估

生成多层次可视化图表，包括训练过程损失曲线、轴承完整退化历程、预测轨迹对比、剩余使用寿命误差分析等，计算各类统计性能指标，形成综合评估报告，为维护决策提供直观的数据支持。

# Import necessary modules
import os
import time
import scipy.io
import scipy.stats
import pywt
from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
from sklearn.preprocessing import StandardScaler
from matplotlib.patches import Rectangle

print("PyTorch Version:", torch.__version__)

# Load data files
PHM_path = 'PHM'
PHM_bearing_files = [os.path.join(PHM_path, file) for file in os.listdir(PHM_path)]

# Enhanced feature extraction
def mat_to_arr_enhanced(file):
    """Enhanced feature extraction"""
    h = scipy.io.loadmat(file)['h'].reshape(-1)
    h2 = h.reshape(-1, 2560)

    # Basic features
    kurtosis = np.array([scipy.stats.kurtosis(i) for i in h2])
    rms = np.array([np.mean(i**2)**0.5 for i in h2])
    rms = np.convolve(rms, [0.3, 0.4, 0.3], mode='same')
    ma = np.array([np.max(np.abs(i)) for i in h2])

    # Time-frequency features
    wavelet_features = []
    for segment in h2:
        coeffs = pywt.wavedec(segment, 'db4', level=3)
        energies = [np.sum(c**2) for c in coeffs]
        wavelet_features.append(energies)
    wavelet_features = np.array(wavelet_features)

    # Frequency domain features
    freq_features = []
    for segment in h2:
        fft_vals = np.abs(np.fft.rfft(segment))
        freq_features.append([
            np.sum(fft_vals[:10]),
            np.sum(fft_vals[10:50]),
            np.sum(fft_vals[50:]),
            np.argmax(fft_vals)
        ])
    freq_features = np.array(freq_features)

    # Combine features
    all_features = np.concatenate([
        rms.reshape(-1, 1),
        ma.reshape(-1, 1),
        kurtosis.reshape(-1, 1),
        wavelet_features,
        freq_features
    ], axis=1)

    FPT = int(len(h2)) * 1700 / 2560
    print(f"Fault Progression Time (FPT): {FPT:.2f}, Feature Dimension: {all_features.shape}")

    return h, FPT, all_features

# Simplified Riemannian Regularization
class LightweightRiemannianPINN(nn.Module):
    """Lightweight Riemannian-regularized PINN with stability improvements"""

    def __init__(self, input_dim=1, hidden_dim=32, dropout_rate=0.3):
        super(LightweightRiemannianPINN, self).__init__()

        self.input_dim = input_dim
        self.hidden_dim = hidden_dim

        # Simplified network
        self.network = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.Tanh(),
            nn.Dropout(dropout_rate),
            nn.Linear(hidden_dim, hidden_dim // 2),
            nn.Tanh(),
            nn.Dropout(dropout_rate),
            nn.Linear(hidden_dim // 2, 1)
        )

        # Gentle physics weight
        self.physics_weight = nn.Parameter(torch.tensor(0.05))

        # Lightweight Riemannian regularization parameter
        self.riemannian_weight = nn.Parameter(torch.tensor(0.01))

        # Initialize weights
        self._initialize_weights()

        print(f"Lightweight Riemannian PINN: Hidden Layer={hidden_dim}, Dropout={dropout_rate}")

    def _initialize_weights(self):
        """Stable initialization"""
        for m in self.modules():
            if isinstance(m, nn.Linear):
                nn.init.xavier_uniform_(m.weight, gain=nn.init.calculate_gain('tanh'))
                if m.bias is not None:
                    nn.init.constant_(m.bias, 0)

    def forward(self, x):
        """Forward propagation - maintain reasonable output range"""
        output = self.network(x)
        # Use gentler output limits
        output = torch.sigmoid(output) * 3.0  # Limit to 0-3, more consistent with actual RMS range
        return output.squeeze(-1)

    def riemannian_regularization(self, features):
        """Lightweight Riemannian regularization"""
        if len(features) < 2:
            return torch.tensor(0.0, device=features.device)

        # Calculate similarity between features (using cosine similarity)
        features_norm = F.normalize(features, p=2, dim=1)
        similarity = torch.mm(features_norm, features_norm.t())

        # Riemannian regularization: maintain local neighborhood structure
        riemannian_loss = torch.mean(1.0 - similarity) * 0.01  # Very small weight

        return riemannian_loss * torch.sigmoid(self.riemannian_weight)

    def physics_loss(self, x, y_pred):
        """Gentle physics constraints"""
        x_tensor = x.clone().detach().requires_grad_(True)
        y_pred_tensor = self.forward(x_tensor)

        # First-order derivative
        dy_dx = torch.autograd.grad(
            y_pred_tensor, x_tensor,
            grad_outputs=torch.ones_like(y_pred_tensor),
            create_graph=True,
            retain_graph=True
        )[0]

        # Gentle non-negative constraint (allow small fluctuations)
        loss_non_negative = torch.mean(F.relu(-dy_dx + 0.2)) * 0.1

        return loss_non_negative

担任《Mechanical System and Signal Processing》《中国电机工程学报》《宇航学报》《控制与决策》等期刊审稿专家，擅长领域：信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测

参考文章：

基于黎曼几何正则化物理信息神经网络的滚动轴承退化趋势预测（Pytorch） - 哥廷根数学学派的文章 
https://zhuanlan.zhihu.com/p/1987166729718239698