去年接手某个轴承退化趋势预测任务时，我被一个现象卡了很久：振动特征在健康阶段波动极小，各种深度学习模型的预测都贴着实测走；可一旦设备进入早期退化段，模型的预测曲线就开始“犹豫”，要么滞后几个采样周期，要么把退化趋势平滑成一条近似直线，等到真正发生严重磨损时，RUL还停留在“基本健康”的状态。数据量并不少，标签也不至于完全缺失，但趋势就是追不上去。后来我把黎曼几何和物理信息神经网络的思路结合起来，做成了Riemannian-PINN，在同一批数据上用Pytorch重跑，退化中后期的趋势才算真正稳了下来。
这篇文章会拆开讲三件事：为什么传统数据驱动模型在退化趋势预测上会翻车；为什么要把特征放到SPD流形上用黎曼距离度量；以及一个完整的、能跑的Riemannian-PINN训练管线长什么样。适合正在做设备健康管理、故障预测与健康管理、或者对PINN变体感兴趣的朋友参考。这里面的实现细节和踩坑经验，我会尽量写透，方便你直接迁移到自己的数据集上。
1. 退化趋势预测难在哪：三个看似平常却致命的假设
1.1 训练集与测试集的分布漂移
滚动轴承退化预测本质上是一个时间序列回归问题，但和股票预测、气象预报不一样的是，轴承的退化过程几乎不存在“重复样本”。同一型号轴承在相同工况下，由于材料微观缺陷、装配误差、润滑状态不同，退化轨迹的起始时间、加速点、寿命长度都差异很大。更麻烦的是，训练集里往往是某一台设备从健康到失效的全周期数据，测试集却是另一台设备从健康到失效的数据，两者的特征分布天然不一致。
传统LSTM、GRU这类模型，默认训练数据和测试数据同分布。它们会把训练集里“第300天、RMS=2.1”这种具体数值记住，一旦测试设备的早期退化特征和训练设备不同，模型就会拿健康阶段的模式去解释，预测结果自然偏向保守或滞后。分布漂移是所有数据驱动方法在做健康管理时绕不过去的墙，单纯堆模型容量并不能解决它。
1.2 标签稀疏且定义模糊
退化趋势预测的标签通常有两种形式：一种是剩余使用寿命（RUL）百分比，从1递减到0；一种是健康指标（HI），从0递增到1。问题在于，标签在哪个时间点开始“算退化”并没有客观标准。有人用工况阈值切分，有人用3sigma准则检测突变点，但不同方法切出来的退化起点可能相差几十个周期。
更棘手的是，大部分退化数据集只有末端附近有相对可信的标签，前中期几乎都是模糊的“健康”或“轻微退化”。数据驱动模型在这种标签分布下很容易学到一种偷懒的解：把所有时刻都预测成健康值，因为这样在大部分标签点上误差都很小。这也是很多模型验证集指标不错、实际部署却完全失灵的原因之一。
1.3 高维特征在欧氏空间里的距离失真
做退化预测通常不会直接用原始振动波形，而是提取时域、频域、时频域统计量，比如RMS、峭度、峰值因子、谱峭度、包络谱特征等。这些特征在高维空间里构成一个可能存在非线性流形结构的点云。大多数模型默认用欧氏距离度量样本间差异，但真实退化轨迹往往不是直线变化的，欧氏距离会低估或高估某些阶段的真实退化程度。
举个直观的例子：两个时间窗的振动特征，一个来自健康阶段，一个来自早期退化阶段，它们的欧氏距离可能很小；而两个都来自退化末期的特征，欧氏距离反而很大。如果模型完全依赖欧氏距离判断状态，就容易把早期退化当成健康抖动，把晚期状态差异当成不同故障模式。这就是为什么需要在特征表示层面引入更贴合数据几何结构的度量方式。
2. 为什么选择SPD流形：把退化特征放进正确的几何空间
2.1 协方差矩阵为什么是SPD矩阵
在滚动轴承振动监测中，我们通常采集多个通道的振动信号，常见的有加速度计的x、y、z三向，或者同一测点的多组传感器。对每个滑动时间窗，计算各通道之间的协方差矩阵，会得到一个对称正定（SPD）矩阵。协方差矩阵描述的是通道间联合波动特性，天然携带了振动能量分布和相位关系的信息，比单一统计量更完整。
SPD矩阵构成的集合不是普通的平坦空间，而是一个黎曼流形。意思是在这个空间里，两点之间最短路径不是直线，而是测地线；直接拿矩阵元素做减法和欧氏距离，会得到和流形真实结构不一致的结果。比如协方差矩阵A和B，它们的欧氏均值可能不再是正定矩阵，这在物理上根本解释不通。
2.2 欧氏距离在SPD流形上的失效案例
我用一个实际算例说明这个问题。假设某轴承健康窗口的协方差矩阵近似为对角线矩阵，数值都集中在1附近；退化窗口的协方差矩阵则出现明显的能量集中，某个方向的特征值增大到10。此时用欧氏距离计算两者差异，会受矩阵元素绝对值的影响，可能把“健康与早期退化”的距离算得比“早期退化与严重退化”还大，因为前者在多个元素上都有小幅变化，后者只在一个元素上大幅变化。
Log-Euclidean距离是SPD流形上一种常用的黎曼度量，做法是先对矩阵取矩阵对数，再计算欧氏距离。这个操作的物理含义是：它把矩阵间的差异从“绝对数值差异”变成“对数域的相对差异”，避免了大数值元素主导距离计算的问题。实测中，同一对SPD矩阵，欧氏距离和Log-Euclidean距离给出的大小排序经常完全不同，这直接影响到后续聚类、分类和趋势预测的效果。
2.3 工程上构造SPD特征的三个细节
第一，协方差矩阵要保证数值正定。实际计算中，由于窗口长度有限或数值精度问题，协方差矩阵可能出现极小的负特征值，解决办法是加上一个很小的单位阵系数，比如1e-5或1e-6。第二，时间窗长度的选择要兼顾稳定性和敏感性。窗口太短，协方差估计噪声大；窗口太长，退化细节被平滑掉。经验值是取设备旋转周期的10到20倍，具体需要结合采样率实验。
第三，矩阵对数可以作为一种输入特征变换。由于SPD流形不是线性的，直接把协方差矩阵的元素拼成向量喂给神经网络，网络很难学到流形结构。先把矩阵映射到对数域，再取上三角元素展开成向量，相当于把流形上的点投影到其切空间，这样既保留了黎曼距离的结构，又让特征向量进入网络时是“友好”的欧氏表示。这段代码我在后面第4章会详细给出。
3. 把退化动力学写进损失函数：PINN约束设计
3.1 物理信息神经网络为什么能补数据不足
物理信息神经网络（PINN）的核心思想很朴素：如果某个系统遵循已知的动力学规律，那就在损失函数里加一项“违背规律”的惩罚项，强迫网络的输出不仅拟合有限标签，还要满足微分方程约束。对于滚动轴承退化问题，观测数据难免稀疏、带噪声，但退化趋势总体上是有规律可循的，比如早期缓慢、中后期加速，不会完全随机。
有了物理约束，网络就相当于被限制在一个合理的解空间里，不再可能把所有时刻都预测成健康值。因为退化曲线一旦保持水平，物理残差会立刻放大，模型必须让输出随时间呈现符合退化方程的变化。这在标签稀疏的前中期尤其重要，物理约束起到了伪标签的作用。
3.2 退化控制方程怎么选：受限增长模型
轴承退化曲线大致呈S形：初期退化非常缓慢，中期逐渐加速，后期接近寿命终点时变化速率再次放缓（受限于最大退化程度）。这种形态用受限增长模型描述最自然：
dy/dt = r * y * (1 - y / K)
其中y是归一化退化指标，范围在0到1；r是退化速率系数；K是退化上限，通常设为1。这个方程的直观含义是：退化指标越小，增长越慢；越接近上限，增长也被抑制；中间阶段增长最快。相比纯指数模型，受限增长模型不会出现退化指标无限膨胀的离谱预测，更适合做寿命末期校准。
实际使用时，r需要根据设备历史数据大致估计。有一个快速标定的办法：取退化中段两个标签点，用差分近似dy/dt，再除以对应时刻的y*(1-y)，就能得到初始r值。代码里我会把它设成一个可训练的小参数，或者直接作为常数传入。
3.3 自动微分求物理残差
PINN实现的关键一步是用自动微分计算网络输出对时间的导数。Pytorch的torch.autograd.grad可以方便地完成这一点。假设网络输入包含连续时间t，输出退化指标y，那么对t求梯度就能得到dy/dt，然后和方程右侧的ry(1-y/K)计算MSE，就是物理残差损失。
这里有个细节：t在输入网络前要做归一化处理，让数值落在0到1之间，否则梯度会被时间量纲拉偏。同时，求导时要设置create_graph=True，因为物理损失本身还要参与反向传播，需要保留计算图。这个点新手很容易漏，漏了之后梯度链断裂，loss能算出来但模型参数完全不更新。
3.4 黎曼几何正则：让预测速率对齐特征流形变化
只加物理约束还不够，因为物理方程里的r是全局常数，无法刻画不同设备、不同退化阶段之间的个体差异。这时候就需要几何信息来“校准”退化速率。我的做法是：对相邻时间窗的SPD特征矩阵，计算它们的Log-Euclidean距离，作为该时段特征流形上的实际变化幅度；同时计算网络预测的退化指标在该时段的增量，让两者成比例对齐。
这个约束的直观解释是：特征空间里走得越快，退化指标就应该涨得越快。它把“物理规律”和“特征几何”两套信息绑在同一个约束里，既避免了纯物理方程过于僵化，又弥补了纯数据模型缺乏机制约束的问题。实际数据里，只要SPD特征构造合理，这个几何正则能明显改善早期退化阶段的趋势捕捉能力。
4. Pytorch实现Riemannian-PINN：可运行的训练管线
4.1 数据管线：滑动窗口、SPD特征提取与标签生成
先把数据结构设计清楚。假设原始数据是多通道振动信号，形状为(num_samples, n_channels)，采样率固定。先用滑动窗口切段，每个窗口计算一个SPD协方差矩阵，再通过矩阵对数映射成特征向量。标签按照设备全寿命周期归一化，把健康起点标为0，失效终点标为1，中间线性插值。下面是SPD特征提取的核心代码：
import torch import torch.nn as nn import numpy as np def sym_logm(X, eps=1e-6): # X: (batch, n, n) SPD矩阵 evals, evecs = torch.linalg.eigh(X) evals = torch.clamp(evals, min=eps) log_evals = torch.log(evals) # 利用 X = V diag(λ) V^T，logX = V diag(logλ) V^T return (evecs * log_evals.unsqueeze(-2)) @ evecs.transpose(-2, -1) def extract_spd_features(windows): # windows: (num_windows, n_channels, window_len) cov_mats = [] for X in windows: Xc = X - X.mean(dim=-1, keepdim=True) C = (Xc @ Xc.transpose(-1, -2)) / (X.shape[-1] - 1) C = C + 1e-5 * torch.eye(X.shape[-2]) cov_mats.append(C) cov_mats = torch.stack(cov_mats) # (N, n, n) log_mats = sym_logm(cov_mats) # (N, n, n) n = log_mats.shape[-1] idx = torch.triu_indices(n, n) feat_vec = log_mats[:, idx[0], idx[1]] # 取上三角，展平 return cov_mats, feat_vecsym_logm里用特征分解实现矩阵对数，比直接调用SciPy稳定得多，而且支持反向传播。取上三角是为了避免重复输入冗余的下三角元素，同时保留了SPD矩阵对数域的完整信息。如果你的数据是多通道振动，这一步基本可以直接复用。
4.2 模型定义：时间与Log-Euclidean特征的双输入结构
网络结构采用简单的多层感知机，输入拼接连续时间t和Log-Euclidean特征向量z，输出归一化退化指标y。为什么不让网络直接吃cov_mats？因为SPD流形的非线性几何结构会让普通全连接层学得很吃力，而经过矩阵对数投影后，特征向量已经落在线性切空间里，网络只需要拟合一个相对平滑的映射。
class RiemPINN(nn.Module): def __init__(self, feat_dim, hidden=128, num_layers=4): super().__init__() dims = [1 + feat_dim] + [hidden] * num_layers + [1] layers = [] for i in range(len(dims) - 1): layers.append(nn.Linear(dims[i], dims[i + 1])) if i < len(dims) - 2: layers.append(nn.Tanh()) else: layers.append(nn.Sigmoid()) # 退化指标限制在[0,1] self.net = nn.Sequential(*layers) def forward(self, t, z): # t: (batch, 1)，z: (batch, feat_dim) x = torch.cat([t, z], dim=-1) return self.net(x)激活函数选择Tanh而不是ReLU，是因为退化曲线是连续平滑的，Tanh的平滑性更适合做自动微分求导。最后一层用Sigmoid把输出限制在0到1，和归一化标签对齐。隐藏层宽度128、层数4是经验值，数据量小的场景可以缩到64和3，数据量大的话适度加宽。
4.3 组合损失函数：数据拟合、物理残差、黎曼正则
这是整个模型的核心。损失包含三部分：有监督数据点的MSE、退化方程残差、SPD流形变化速率对齐。物理残差依赖自动微分，黎曼正则依赖Log-Euclidean距离。三个损失各自解决一个问题，最后的组合权重需要细心调度。
def log_euclidean_dist(A, B): # A, B: (batch, n, n) SPD矩阵 diff = sym_logm(A) - sym_logm(B) return torch.norm(diff, dim=(-2, -1)) # (batch,) def pinn_loss(model, batch_data, batch_phys, batch_geom, r=0.02, K=1.0): # batch_data: (t_data, z_data, y_data) t_data, z_data, y_data = batch_data y_pred = model(t_data, z_data) loss_data = nn.functional.mse_loss(y_pred, y_data) # batch_phys: (t_phys, z_phys) t_phys, z_phys = batch_phys t_phys.requires_grad_(True) y_phys = model(t_phys, z_phys) dy_dt = torch.autograd.grad( y_phys.sum(), t_phys, grad_outputs=torch.ones_like(y_phys), create_graph=True )[0] f_phys = r * y_phys * (1 - y_phys / K) loss_phys = nn.functional.mse_loss(dy_dt, f_phys) # batch_geom: (A_spd, B_spd, z_a, z_b, t_a, t_b, scale) A_spd, B_spd, z_a, z_b, t_a, t_b, scale = batch_geom d_riemann = log_euclidean_dist(A_spd, B_spd) y_a = model(t_a, z_a) y_b = model(t_b, z_b) delta_y = (y_b - y_a).squeeze(-1) target_delta = scale * d_riemann loss_geom = nn.functional.mse_loss(delta_y, target_delta) return loss_data, loss_phys, loss_geomscale是退化指标增量与SPD特征距离之间的比例常数。实践中可以用带标签的样本做一个最小二乘预标定：先算历史相邻窗口的特征距离，再算对应的标签增量，拟合一个线性系数。代码注释里建议先跑一遍带标签段，标定好scale再进入正式训练。
关于物理损失，我特别提醒一点：r的取值会影响退化曲线的整体陡峭度。如果r设得过大，物理残差会强行把曲线拉成陡峭的S形，即使数据标签不支持；如果r过小，退化趋势又会被过度平滑。稳妥的做法是把r定在0.01到0.05之间，观察物理残差的数量级，再做微调。
4.4 训练主循环与损失权重调度
训练主循环要注意两个问题：一是类别不平衡，三个损失的量纲不同，需要分别用权重配平；二是早期训练阶段，网络还没把数据拟合好，物理和几何约束加入太早会干扰收敛。我的做法是：数据损失权重固定为1，物理损失权重随训练进程从0.01线性增加到0.5，几何损失权重在第50个epoch之后才启用。
def train_pinn(model, dl_data, dl_phys, dl_geom, epochs=300): opt = torch.optim.Adam(model.parameters(), lr=1e-3) sched = torch.optim.lr_scheduler.CosineAnnealingLR(opt, T_max=epochs) for epoch in range(epochs): w_phys = 0.01 + 0.49 * min(epoch / 200, 1.0) w_geom = 0.1 if epoch >= 50 else 0.0 for batch_data, batch_phys, batch_geom in zip(dl_data, dl_phys, dl_geom): loss_data, loss_phys, loss_geom = pinn_loss( model, batch_data, batch_phys, batch_geom ) loss = loss_data + w_phys * loss_phys + w_geom * loss_geom opt.zero_grad() loss.backward() nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0) opt.step() sched.step() if epoch % 50 == 0: print(f"epoch {epoch}: " f"data={loss_data.item():.4f}, " f"phys={loss_phys.item():.4f}, " f"geom={loss_geom.item():.4f}")梯度裁剪非常重要。物理残差对时间求导之后，梯度数值放大了不少，如果不裁剪，训练很容易在初期就发散。Clip到1.0是我在这个任务上试过比较稳的阈值，如果你的数据尺度不同，可以从0.5到2.0之间调。学习率用余弦退火，可以让后期的损失更低，避免在平台期来回震荡。
评估时不要只看整体RMSE，要分阶段看：健康段、早期退化段、中后期退化段分别计算误差。退化趋势预测最容易出问题的是中后期预测偏保守，如果分段误差显示中后期MAE明显偏大，大概率是物理权重不够或者r参数偏小，对应调整即可。
5. 训练过程踩过的坑：权重失衡、梯度爆炸与平台期
5.1 物理残差权重过大导致拟合完全失效
我第一次跑完整管线时，把物理损失权重设成和数据损失相同，结果训练到100个epoch时，预测曲线虽然形态很“标准”的S形，但数值整体偏低，末期退化指标只有0.6左右，和真实标签0.9完全对不上。原因是物理方程是一个全局约束，它要求所有时刻的y都满足dy/dt = ry(1-y/K)，如果这个约束权重过大，模型就会牺牲局部的数据拟合来换取整体曲线形态的“正确”。
解决方式是物理权重从小到大退火。前20个epoch让数据损失主导，网络先大致拟合标签；然后逐渐增加物理约束，让曲线形态从“噪声轨迹”慢慢被拉成平滑的S形。这条经验几乎可以迁移到所有PINN类任务：物理约束是指导性的，不是主导性的，它应该帮助网络在数据稀疏处保持合理外推，而不是取代数据拟合。
5.2 矩阵对数特征值过小时数值不稳定
SPD矩阵的特征值如果非常接近0，sym_logm虽然做了clamp，但梯度在log(λ)的导数位置会变得异常大，导致反向传播不稳定。我遇到过训练十几个epoch后loss突然变成NaN的情况，排查到最终原因是某个时间窗的振动信号接近静默，协方差矩阵的所有元素都趋向0，特征值数量级掉到1e-8，log之后数值直接发散。
两个补救手段：一是把clamp的下限从1e-6提高到1e-4，虽然损失一点精度，但换来稳定；二是在特征提取前，对原始振动信号做标准化，减去均值后除以标准差，让协方差矩阵的元素量级更稳定。如果你用的是多传感器采集的原始波形，这条尤其重要，不同测点的量纲差异会让SPD特征矩阵奇异性明显增加。
5.3 退化起点不一致导致几何正则失控
黎曼几何正则的设计前提是：“特征流形距离越大，退化指标增量越大”。但不同设备进入退化的起点不一致，有的设备在第100天就开始退化，有的第300天才开始。如果几何约束直接套在所有时间窗对上，健康阶段特征的自然漂移会被几何正则误判成退化增量，预测曲线会出现早期虚高的虚假退化。
我的处理方法是给几何正则加一个门控：只有当SPD特征距离超过某个健康参考阈值时，才激活几何对齐损失。这个阈值可以用全寿命数据里前10%窗口的Log-Euclidean距离分布来标定，取95分位数作为“健康漂移上限”。低于阈值的样本对，几何损失直接置0。加了门控之后，早期退化段的虚高预测明显减少。
5.4 训练后期进入平台期：不是过拟合，是约束打架
训练到250个epoch后，三个损失各自看起来都平稳了，但预测曲线始终和真实标签隔着一层薄膜，末期趋势有，就是不够陡。后来我把三个损失的梯度方向单独拿出来分析，发现数据损失在推着曲线上调，物理损失却因为r参数固定而限制曲线增速，几何损失又试图保持特征距离对齐，三个梯度方向在后期相互抵消，形成了一个隐性死锁。
解法是分阶段解耦：前150个epoch把r当成可训练参数，让它在反向传播中学习；150个epoch后冻结r，只更新网络权重。实测下来，末期曲线的陡峭度明显改善，验证集分段误差比之前降低了大概三成。这种“先学参数、再固定约束”的思路，在处理多约束联合训练时非常实用。
5.5 关于评价指标的提醒：RMSE不是一切
最后说一个评估层面的经验。退化趋势预测的最终目标是给维护决策提供依据，所以真正重要的不是所有点的平均误差，而是退化加速点的定位准不准、末期RUL偏不偏。我曾经用RMSE作为唯一指标调参，调到后期RMSE很漂亮，但退化加速点比实际晚了将近50个周期，这类模型在真实场景里根本不敢用。
建议在验证阶段同时看三个指标：整体RMSE、退化加速点定位误差、末期30%阶段的平均绝对误差。如果后两个指标不达标，即使整体RMSE再好看，也要回头调物理权重或几何门控。这套Riemannian-PINN框架的价值，不在于它一定比LSTM的RMSE低多少，而在于它用物理规律和特征几何把预测约束在“像真实退化”的轨迹上，减少那些看似合理实则危险的误判。做工业健康管理的朋友应该都懂，能让模型在关键时刻不犯错，比让它平时更精确重要得多。
