import numpy as np
import matplotlib.pyplot as plt
import os
import pandas as pd
from tqdm import tqdm
import warnings
from scipy.signal import convolve
import pywt
from scipy.fft import fft, fftfreq
from matplotlib.colors import LinearSegmentedColormap
# 全局设置
warnings.filterwarnings("ignore")
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
class MODWPTProcessor:
def __init__(self, wavelet='db10', level=5):
self.wavelet = wavelet
self.level = level
self.filters = self._get_wavelet_filters()
def _get_wavelet_filters(self):
wavelet_obj = pywt.Wavelet(self.wavelet)
return wavelet_obj.dec_lo, wavelet_obj.dec_hi
def _insert_zeros(self, filter_coeffs, j):
zeros_to_insert = 2 ** (j - 1) - 1
if zeros_to_insert <= 0:
return filter_coeffs
new_length = len(filter_coeffs) * (zeros_to_insert + 1)
new_filter = np.zeros(new_length)
for i, coeff in enumerate(filter_coeffs):
new_filter[i * (zeros_to_insert + 1)] = coeff
return new_filter
def modwpt_decomposition(self, signal):
M = len(signal)
scaling_filter, wavelet_filter = self.filters
H_current = {0: signal.copy()}
for j in range(1, self.level + 1):
H_next = {}
s_tilde = self._insert_zeros(scaling_filter, j)
w_tilde = self._insert_zeros(wavelet_filter, j)
for n in range(2 ** (j - 1)):
parent_coeff = H_current[n]
H_low = convolve(parent_coeff, s_tilde, mode='same', method='auto')
H_high = convolve(parent_coeff, w_tilde, mode='same', method='auto')
if n % 4 == 0 or n % 4 == 3:
H_next[2 * n] = H_low
H_next[2 * n + 1] = H_high
else:
H_next[2 * n] = H_high
H_next[2 * n + 1] = H_high
H_current = H_next
energies = {}
for n in range(2 ** self.level):
coeff = H_current[n]
energies[n] = np.sum(coeff ** 2)
return energies
def read_single_csv_file(file_path):
"""读取单个CSV文件中的水平方向振动数据（第5列）"""
try:
data = np.loadtxt(file_path, delimiter=',', skiprows=0)
except Exception as e:
raise ValueError(f"CSV文件读取失败（{file_path}）：{str(e)}")
if data.ndim < 2:
raise ValueError(f"CSV文件数据格式错误（{file_path}），未检测到有效列")
if data.shape[1] < 5:
raise ValueError(f"CSV文件列数不足（{data.shape[1]}列），需至少5列以提取水平方向数据")
horizontal_data = data[:, 4]
horizontal_data = horizontal_data[~np.isnan(horizontal_data)]
if len(horizontal_data) == 0:
raise ValueError(f"水平方向数据为空（{file_path}），请检查数据完整性")
return horizontal_data
def process_all_files(directory_path):
"""处理目录下的所有CSV文件"""
csv_files = [f for f in os.listdir(directory_path) if f.endswith('.csv')]
csv_files.sort()
print(f"找到 {len(csv_files)} 个CSV文件，开始处理...")
processor = MODWPTProcessor(wavelet='db10', level=5)
num_bands = 2 ** processor.level
energy_matrix = []
file_indices = []
for i, file_name in enumerate(tqdm(csv_files, desc="处理CSV文件")):
file_path = os.path.join(directory_path, file_name)
try:
signal = read_single_csv_file(file_path)
energies = processor.modwpt_decomposition(signal)
energy_vector = [energies[band_idx] for band_idx in range(num_bands)]
energy_matrix.append(energy_vector)
file_indices.append(i)
except Exception as e:
print(f"处理文件 {file_name} 时出错: {e}")
continue
return np.array(energy_matrix), file_indices, csv_files
def normalize_energy_matrix(energy_matrix):
"""将能量矩阵进行归一化：95%正常范围映射到0-1，5%异常映射到1-2"""
# 计算整个能量矩阵的95%分位数
threshold_95 = np.percentile(energy_matrix, 95)
# 创建归一化后的矩阵
normalized_matrix = np.zeros_like(energy_matrix)
# 对每个元素进行归一化
for i in range(energy_matrix.shape[0]):
for j in range(energy_matrix.shape[1]):
if energy_matrix[i, j] <= threshold_95:
# 正常范围：线性映射到 [0, 1]
normalized_matrix[i, j] = energy_matrix[i, j] / threshold_95
else:
# 异常范围：线性映射到 (1, 2]
max_val = np.max(energy_matrix)
if max_val > threshold_95:
normalized_matrix[i, j] = 1 + (energy_matrix[i, j] - threshold_95) / (max_val - threshold_95)
else:
normalized_matrix[i, j] = 1.0
return normalized_matrix, threshold_95
def create_custom_colormap():
"""创建自定义颜色映射：正常范围用蓝色系，异常范围用红色系"""
colors = [
(0, 0, 0.5), # 深蓝 - 最低能量
(0, 0.5, 1), # 天蓝
(0.5, 0.8, 1), # 浅蓝
(1, 1, 0.5), # 黄色 - 正常/异常分界
(1, 0.8, 0), # 橙色
(1, 0.4, 0), # 红色
(0.8, 0, 0) # 深红 - 最高能量
]
return LinearSegmentedColormap.from_list('custom_energy', colors, N=256)
def plot_representative_results(energy_matrix, file_indices, sample_rate=25600):
"""绘制最具代表性的图表（使用调整后的幅值映射）"""
# 归一化能量矩阵
normalized_energy, threshold_95 = normalize_energy_matrix(energy_matrix)
num_files, num_bands = normalized_energy.shape
print(f"能量归一化参数：95%分位数 = {threshold_95:.4f}")
print(f"归一化后能量范围: [{np.min(normalized_energy):.4f}, {np.max(normalized_energy):.4f}]")
bandwidth = sample_rate / (2 ** (5 + 1)) # L=5
frequencies = [i * bandwidth for i in range(num_bands)]
# 创建自定义颜色映射
custom_cmap = create_custom_colormap()
# 1. 窄带能量瀑布图（使用归一化能量）
plt.figure(figsize=(14, 8))
X, Y = np.meshgrid(file_indices, frequencies)
Z = normalized_energy.T
# 设置颜色映射范围为0-2
im = plt.pcolormesh(X, Y, Z, shading='auto', cmap=custom_cmap, vmin=0, vmax=2)
# 添加颜色条和分界线标注
cbar = plt.colorbar(im, label='归一化能量 (0-1:正常, 1-2:异常)', shrink=0.8)
cbar.ax.plot([0, 1], [1, 1], 'w--', linewidth=2) # 在颜色条上标记分界线
plt.xlabel('时间序列（文件序号）', fontsize=12)
plt.ylabel('频率 (Hz)', fontsize=12)
plt.title('滚动轴承归一化窄带能量瀑布图\n(95%正常范围→0-1, 5%异常范围→1-2)',
fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
plt.ylim(0, frequencies[-1] + bandwidth)
# 添加分界线说明
plt.text(0.02, 0.98, f'95%分界线: {threshold_95:.2e}',
transform=plt.gca().transAxes, fontsize=10,
bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.8))
plt.tight_layout()
plt.show()
# 2. 总能量趋势和关键频带能量变化
plt.figure(figsize=(15, 10))
# 2.1 总能量趋势（使用原始能量，但标注分界线）
plt.subplot(2, 2, 1)
total_energy = np.sum(energy_matrix, axis=1)
threshold_total = np.percentile(total_energy, 95)
plt.plot(file_indices, total_energy, 'k-', linewidth=2, label='总能量')
plt.axhline(y=threshold_total, color='r', linestyle='--',
linewidth=2, label=f'95%分界线 ({threshold_total:.2e})')
# 标记异常点
abnormal_indices = np.where(total_energy > threshold_total)[0]
if len(abnormal_indices) > 0:
plt.scatter([file_indices[i] for i in abnormal_indices],
[total_energy[i] for i in abnormal_indices],
color='red', s=30, zorder=5, label='异常点')
plt.xlabel('时间序列（文件序号）')
plt.ylabel('总振动能量')
plt.title('轴承总振动能量趋势\n(红色虚线为95%分界线)')
plt.legend()
plt.grid(True, alpha=0.3)
# 2.2 关键频带归一化能量变化
plt.subplot(2, 2, 2)
# 选择能量最大的几个频带
mean_energy_per_band = np.mean(normalized_energy, axis=0)
top_bands = np.argsort(mean_energy_per_band)[-5:] # 选择能量最大的5个频带
colors = ['red', 'blue', 'green', 'orange', 'purple']
for i, band in enumerate(top_bands):
plt.plot(file_indices, normalized_energy[:, band],
color=colors[i], linewidth=1.5,
label=f'频带{band}({frequencies[band]:.0f}Hz)')
# 添加异常分界线
plt.axhline(y=1, color='black', linestyle='--', linewidth=2, alpha=0.7,
label='正常/异常分界线')
plt.xlabel('时间序列（文件序号）')
plt.ylabel('归一化窄带能量')
plt.title('关键频带归一化能量变化\n(>1表示异常)')
plt.legend()
plt.grid(True, alpha=0.3)
# 2.3 归一化能量分布对比（初期 vs 后期）
plt.subplot(2, 2, 3)
if num_files > 10:
# 取前10%作为初期，后10%作为后期
early_idx = max(1, num_files // 10)
late_idx = num_files - early_idx
early_energy = np.mean(normalized_energy[:early_idx], axis=0)
late_energy = np.mean(normalized_energy[late_idx:], axis=0)
plt.plot(frequencies, early_energy, 'b-', linewidth=2, label='运行初期')
plt.plot(frequencies, late_energy, 'r-', linewidth=2, label='运行后期')
plt.axhline(y=1, color='black', linestyle='--', linewidth=2, alpha=0.7,
label='正常/异常分界线')
plt.xlabel('频率 (Hz)')
plt.ylabel('平均归一化能量')
plt.title('运行初期vs后期归一化能量分布对比')
plt.legend()
plt.grid(True, alpha=0.3)
# 2.4 故障发展指标（基于归一化能量）
plt.subplot(2, 2, 4)
# 计算异常能量占比
abnormal_ratio = np.sum(normalized_energy > 1, axis=1) / num_bands * 100
plt.plot(file_indices, abnormal_ratio, 'g-', linewidth=2, label='异常频带占比')
# 计算移动平均以显示趋势
window_size = min(20, len(abnormal_ratio) // 10)
if window_size > 1:
moving_avg = np.convolve(abnormal_ratio, np.ones(window_size) / window_size, mode='valid')
plt.plot(file_indices[window_size - 1:], moving_avg, 'r--', linewidth=2,
label=f'{window_size}点移动平均')
plt.xlabel('时间序列（文件序号）')
plt.ylabel('异常频带占比 (%)')
plt.title('故障发展指标\n(归一化能量>1的频带占比)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 3. 3D归一化能量瀑布图
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(file_indices, frequencies)
Z = normalized_energy.T
# 使用surface plot
surf = ax.plot_surface(X, Y, Z, cmap=custom_cmap, alpha=0.8,
linewidth=0, antialiased=True, vmin=0, vmax=2)
ax.set_xlabel('时间序列')
ax.set_ylabel('频率 (Hz)')
ax.set_zlabel('归一化能量')
ax.set_title('3D归一化能量瀑布图\n(0-1:正常范围, 1-2:异常范围)')
ax.set_zlim(0, 2)
plt.colorbar(surf, ax=ax, shrink=0.5, aspect=20, label='归一化能量')
plt.show()
return total_energy, abnormal_ratio, threshold_95
def print_health_assessment(total_energy, abnormal_ratio, threshold_95):
"""打印健康状态评估结果"""
print("\n" + "=" * 60)
print("轴承健康状态评估报告（基于归一化能量分析）")
print("=" * 60)
# 计算基于归一化能量的指标
total_threshold = np.percentile(total_energy, 95)
is_abnormal_total = total_energy[-1] > total_threshold
abnormal_band_ratio = abnormal_ratio[-1]
print(f"当前总能量: {total_energy[-1]:.2e} (95%分界线: {total_threshold:.2e})")
print(f"异常频带占比: {abnormal_band_ratio:.1f}%")
print(f"能量归一化95%分界线: {threshold_95:.2e}")
# 健康状态判断（基于多个指标）
if not is_abnormal_total and abnormal_band_ratio < 10:
status = "健康"
color = "绿色"
confidence = "高"
elif abnormal_band_ratio < 30:
status = "早期故障"
color = "黄色"
confidence = "中"
else:
status = "严重故障"
color = "红色"
confidence = "高"
print(f"轴承状态评估: {status} ({color}) - 置信度: {confidence}")
# 趋势分析
early_abnormal = np.mean(abnormal_ratio[:len(abnormal_ratio) // 3])
late_abnormal = np.mean(abnormal_ratio[2 * len(abnormal_ratio) // 3:])
trend = "上升" if late_abnormal > early_abnormal + 5 else "稳定" if abs(
late_abnormal - early_abnormal) < 5 else "下降"
print(f"异常趋势: {trend} (早期{early_abnormal:.1f}% → 近期{late_abnormal:.1f}%)")
if status != "健康":
print("\n建议措施:")
if status == "早期故障":
print("- 加强监测频率（建议每周一次）")
print("- 准备维护计划")
print("- 检查润滑状况")
else:
print("- ⚠️ 立即停机检查")
print("- 更换轴承")
print("- 检查相关传动部件")
else:
print("\n建议措施:")
print("- 继续常规监测")
print("- 保持正常维护计划")
# 主程序
if __name__ == "__main__":
# 数据路径
data_directory = r"D:\成电——研究生\基于数据驱动的故障诊断研究\数据集汇总\phm-ieee-2012-data-challenge-dataset-master\Learning_Set\Bearing1_1"
print("开始处理Bearing1_1数据集...")
# 处理所有文件
energy_matrix, file_indices, csv_files = process_all_files(data_directory)
if len(energy_matrix) == 0:
print("未找到有效数据，请检查文件路径和数据格式")
exit()
print(f"成功处理 {len(energy_matrix)} 个文件")
print(f"窄带能量矩阵形状: {energy_matrix.shape} (时间点 × 频带)")
# 绘制代表性图表（使用新的归一化方法）
total_energy, abnormal_ratio, threshold_95 = plot_representative_results(energy_matrix, file_indices)
# 输出健康评估
print_health_assessment(total_energy, abnormal_ratio, threshold_95)分析代码的合理性
最新发布
振动色在界面设计中可能导致元素边缘模糊，影响易读性。本文探讨了振动色的产生、对文本可读性的危害，以及在设计时应考虑的可访问性问题。结论指出，虽然振动色能吸引眼球，但在UI设计中应谨慎使用。
扫一扫
860
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
