希尔伯特相位熵：从相位同步视角捕捉旋转机械故障的隐性特征（Python）

在做故障诊断时，我们习惯盯着频谱看，哪个频率分量冒出来了，能量分布变没变。

但有一次处理一组早期磨损的振动数据，频谱上几乎看不出异常，轴承却已出现肉眼可见的剥落。振幅没怎么变，频率结构也大致正常。

那一刻我意识到：有些故障信息，并不编码在振幅或频率里，而是藏在相位关系中。

01 振幅和频率之外，还有什么？

旋转机械的振动信号是典型的非线性、非平稳过程，轴承元件在运转中产生周期冲击，冲击之间的时间关系、相位耦合模式，蕴含着故障位置和严重程度的丰富信息。

然而，传统的时域统计量和频域分析方法，本质上捕捉的是信号的振幅分布和频率构成。当故障处于早期阶段，冲击能量弱、信噪比低，振幅和频谱的变化往往不够显著。

相位信息，恰恰是一个被长期低视的诊断维度。

相位描述的是信号振动的节奏和时机，不是多强、多快，而是什么时候发生。不同故障类型，在相空间中会呈现不同的相位同步模式。如果能有效量化这种模式差异，就有望在传统指标失效时提供增量的判别依据。

02 什么是希尔伯特相位熵？

2.1 从相空间重构说起

振动信号的一个样本是时间序列。要让相位关系显现出来，第1步是相空间重构：将1维序列通过时间延迟嵌入，映射为高维空间中的轨迹。

给定嵌入维数和延迟，构造状态向量矩阵。矩阵的每1行，是信号在相空间中的1个位置。于是，振动信号的运动轨迹，就被翻译成了相空间中点的迁移序列。

这一步保留了系统的动态演化结构，是后续相位分析的基础。

2.2 希尔伯特变换提取瞬时相位

对相空间中的每个状态向量进行希尔伯特变换，得到解析信号，解析信号的幅角，就是这个状态点的瞬时相位。

瞬时相位的物理直觉是：在相空间的某一瞬间，系统正处在哪个角度。

2.3 相位锁值与同步量化

有了各状态点的瞬时相位，就可以计算相邻状态之间的相位差。定义相位锁值——1个归一化的同步度量，来刻,2个状态之间相位的一致性：

相位锁值越接近1，表示相邻状态高度同步，系统轨迹在相空间中沿稳定方向推进；
相位锁值越接近0，表示相位发散，系统的运动模式趋于随机或混沌。

不同故障状态下，振动信号的相位同步模式存在系统性差异：正常状态往往呈现较高同步性，而故障引起的冲击扰动会降低局部锁值。

2.4 从锁值分布到希尔伯特相位熵

将相位锁值在[0,1]区间等分，统计锁值落入各子区间的频数，得到经验概率分布。对该分布计算归一化香农熵，即为希尔伯特相位熵。

低熵值：相位锁值集中于少数区间，表明信号存在稳定的相位同步模式；
高熵值：锁值分布均匀，表明信号相位关系紊乱。

换句话说，希尔伯特相位熵用1维标量，编码了整段信号在相空间中的相位组织程度。

03 为什么需要插值补偿？

传统多尺度熵分析的做法是：先对信号分段平均粗粒化，再对缩短后的序列提取熵值。随着尺度增大，粗粒化序列长度递减，导致大尺度下的熵估计方差增大、信息可靠性下降。

针对这一问题，引入线性插值补偿的粗粒化策略：

1. 对第k尺度的粗粒化序列，在相邻点间进行线性插值；
2. 构造一条长度约为原序列两倍减一的细化序列；
3. 对细化序列重复相位熵提取流程，得到当前尺度的熵值。

插值的核心作用： 在保持多尺度分析框架的前提下，显著缓解大尺度下因序列过短造成的信息损失，使各尺度的熵估计更稳定。

将所有尺度下的希尔伯特相位熵拼接，即构成多尺度相位熵特征向量。本实验中，选取20个尺度，最终形成20维特征。

def main():
    # 5.1 加载数据
    print("Loading data...")
    X_all, y_all, samples_per_cls = load_all_data()
    total_samples = X_all.shape[0]
    print(f"Total samples: {total_samples}, Feature length: {X_all.shape[1]}")


    # 5.2 时域波形示例图（每类取第一个样本展示）
    fig, axes = plt.subplots(2, 5, figsize=(15, 6))
    axes = axes.ravel()
    for i in range(10):
        idx = np.where(y_all == i)[0][0]   # 该类第一个样本的全局索引
        ax = axes[i]
        ax.plot(X_all[idx][:500])
        ax.set_title(FAULT_NAMES[i], fontsize=9)
        ax.set_xlabel('Sample points')
        ax.set_ylabel('Amplitude')
    plt.suptitle('Time domain waveforms (first 500 points)')
    plt.tight_layout()
    plt.savefig('time_domain_examples.png', dpi=300)
    plt.show()


    # 5.3 提取所有样本的 IMHPE 特征 (整个数据集)
    print("Extracting IMHPE features... (this may take several minutes)")
    features = extract_features(X_all, SCALES)
    print(f"Feature matrix shape: {features.shape}")  # (n_samples, 20)


    # 5.4 画不同故障类型在各尺度的 IMHPE 均值曲线
    plt.figure(figsize=(10, 6))
    for cls in range(10):
        idx = y_all == cls
        mean_curve = features[idx].mean(axis=0)
        plt.plot(SCALES, mean_curve, marker='o', label=FAULT_NAMES[cls])
    plt.xlabel('Scale factor')
    plt.ylabel('IMHPE')
    plt.title('Mean IMHPE across scales for each fault type')
    plt.legend(fontsize=8)
    plt.grid(True)
    plt.tight_layout()
    plt.savefig('imhpe_curves.png', dpi=300)
    plt.show()


    # 5.5 分类实验 (重复10次随机划分)
    acc_list = []
    n_samples = features.shape[0]
    # 存储最后一次实验的预测结果用于混淆矩阵
    last_test_true = None
    last_test_pred = None
    last_test_feats = None


    print("Starting classification experiments...")
    for r in range(NUM_REPEATS):
        # 随机打乱
        idx = np.random.permutation(n_samples)
        feats_shuffle = features[idx]
        labels_shuffle = y_all[idx]


        split = int(n_samples * (1 - TEST_RATIO))
        train_X = feats_shuffle[:split]
        train_y = labels_shuffle[:split]
        test_X = feats_shuffle[split:]
        test_y = labels_shuffle[split:]


        # 标准化（仅在训练集上计算参数，防止数据泄露）
        scaler = StandardScaler()
        train_X_norm = scaler.fit_transform(train_X)
        test_X_norm = scaler.transform(test_X)


        # 转为 one-hot 标签
        train_T = np.eye(10)[train_y]
        test_T = np.eye(10)[test_y]


        # 训练 ELM
        elm = ELM(n_hidden=HIDDEN_NODES)
        elm.fit(train_X_norm, train_T)


        # 预测
        pred = elm.predict(test_X_norm)
        acc = np.mean(pred == test_y)
        acc_list.append(acc)


        if r == NUM_REPEATS - 1:
            last_test_true = test_y
            last_test_pred = pred
            last_test_feats = test_X_norm


    print(f"Average accuracy over {NUM_REPEATS} runs: {np.mean(acc_list):.4f} ± {np.std(acc_list):.4f}")


    # 5.6 t-SNE 可视化 (使用最后一次测试集特征)
    print("Performing t-SNE visualization...")
    tsne = TSNE(n_components=2, perplexity=30, random_state=42)
    feats_2d = tsne.fit_transform(last_test_feats)
    plt.figure(figsize=(8, 6))
    scatter = plt.scatter(feats_2d[:, 0], feats_2d[:, 1],
                          c=last_test_true, cmap='tab10', alpha=0.7)
    plt.colorbar(scatter, ticks=range(10))
    plt.title('t-SNE visualization of IMHPE features (test set)')
    plt.xlabel('Component 1')
    plt.ylabel('Component 2')
    plt.tight_layout()
    plt.savefig('tsne_imhpe.png', dpi=300)
    plt.show()


    # 5.7 混淆矩阵
    cm = confusion_matrix(last_test_true, last_test_pred)
    plt.figure(figsize=(8, 6))
    plt.imshow(cm, interpolation='nearest', cmap=plt.cm.Blues)
    plt.title('Confusion matrix (last run)')
    plt.colorbar()
    tick_marks = np.arange(10)
    plt.xticks(tick_marks, FAULT_NAMES, rotation=45, fontsize=8)
    plt.yticks(tick_marks, FAULT_NAMES, fontsize=8)
    plt.xlabel('Predicted label')
    plt.ylabel('True label')
    # 在格子里标注数值
    thresh = cm.max() / 2.
    for i in range(cm.shape[0]):
        for j in range(cm.shape[1]):
            plt.text(j, i, format(cm[i, j], 'd'),
                     ha="center", va="center",
                     color="white" if cm[i, j] > thresh else "black")
    plt.tight_layout()
    plt.savefig('confusion_matrix.png', dpi=300)
    plt.show()


    # 5.8 不同尺度下的分类准确率分析 (只用单个尺度特征)
    single_scale_accs = []
    print("Evaluating single-scale performance...")
    for s_idx, s in enumerate(SCALES):
        feat_s = features[:, s_idx].reshape(-1, 1)
        temp_acc = []
        for r in range(NUM_REPEATS):
            idx = np.random.permutation(n_samples)
            f_shuffle = feat_s[idx]
            l_shuffle = y_all[idx]
            split = int(n_samples * (1 - TEST_RATIO))
            tr_X, tr_y = f_shuffle[:split], l_shuffle[:split]
            te_X, te_y = f_shuffle[split:], l_shuffle[split:]


            scaler = StandardScaler()
            tr_X_norm = scaler.fit_transform(tr_X)
            te_X_norm = scaler.transform(te_X)


            tr_T = np.eye(10)[tr_y]
            te_T = np.eye(10)[te_y]
            elm = ELM(n_hidden=HIDDEN_NODES)
            elm.fit(tr_X_norm, tr_T)
            pred = elm.predict(te_X_norm)
            temp_acc.append(np.mean(pred == te_y))
        single_scale_accs.append(np.mean(temp_acc))


    plt.figure(figsize=(10, 5))
    plt.plot(SCALES, single_scale_accs, marker='s', color='crimson')
    plt.xlabel('Scale factor')
    plt.ylabel('Accuracy')
    plt.title('Classification accuracy using single-scale IMHPE feature')
    plt.grid(True)
    plt.tight_layout()
    plt.savefig('single_scale_accuracy.png', dpi=300)
    plt.show()


    # 5.9 总结输出
    print("\n===== Experiment Summary =====")
    print(f"Method: IMHPE + ELM")
    print(f"Number of classes: 10")
    print(f"Feature dimension: {len(SCALES)}")
    print(f"Classification accuracy: {np.mean(acc_list):.2%} ± {np.std(acc_list):.2%}")

参考文章：

希尔伯特相位熵：从相位同步视角捕捉旋转机械故障的隐性特征（Python）

如果你对信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测有疑问，或者需要论文思路上的建议，欢迎学术咨询

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