task64_timeseries.py - master - exender/coursework-data-analysis

task64_timeseries.py
313 строк14 KB
exender
Create: task49_nlp.py, task48_images.py, task64_timeseries.py, task58_insurance.py
08 июн 2026, 14:44
08 июн 2026, 14:440500eb6
"""
Задача 64: Раннее обнаружение перегрева подшипника по тренду температуры
Датасет: FEMTO Bearing (IEEE PHM 2012)
Источник: https://www.kaggle.com/datasets/vinayak123tyagi/bearing-dataset

Установка зависимостей:
    pip install pandas numpy matplotlib seaborn scipy statsmodels

Структура папки с датасетом:
    bearing-dataset/
        Bearing1_1/
            acc_00001.csv, ..., temp_00001.csv, ...
        Bearing1_2/ ...

Запуск:
    python task64_timeseries.py
"""

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import seaborn as sns
from scipy import stats
from statsmodels.tsa.seasonal import seasonal_decompose
import os
import warnings
warnings.filterwarnings('ignore')

# ─── 0. Загрузка / генерация данных ─────────────────────────────────────────
def load_bearing_data(data_dir='bearing-dataset'):
    """Загружает температурные ряды из папки датасета FEMTO."""
    bearings = {}
    for bearing_name in os.listdir(data_dir):
        path = os.path.join(data_dir, bearing_name)
        if not os.path.isdir(path):
            continue
        temp_files = sorted([f for f in os.listdir(path) if f.startswith('temp')])
        if not temp_files:
            continue
        frames = []
        for f in temp_files:
            df = pd.read_csv(os.path.join(path, f), header=None,
                             names=['hour', 'minute', 'second', 'us', 'temperature'])
            frames.append(df)
        df_all = pd.concat(frames, ignore_index=True)
        # Создаём временной индекс
        df_all['t_sec'] = (df_all['hour'] * 3600 + df_all['minute'] * 60
                           + df_all['second'] + df_all['us'] / 1e6)
        df_all = df_all.sort_values('t_sec').reset_index(drop=True)
        bearings[bearing_name] = df_all['temperature'].values
    return bearings


def generate_synthetic_bearings():
    """Генерирует синтетические температурные ряды, имитирующие FEMTO."""
    print("Генерируем синтетические температурные ряды (датасет не найден)...")
    np.random.seed(42)
    bearings = {}
    configs = [
        ('Bearing1_1', 2803, 28, 52, 0.005, 'normal'),
        ('Bearing1_2', 1714, 30, 68, 0.009, 'fast'),
        ('Bearing1_3', 2375, 27, 79, 0.007, 'moderate'),
    ]
    for name, n, t_start, t_max, rate, mode in configs:
        t = np.linspace(0, n / 10, n)  # секунды
        # Нелинейный рост ближе к отказу
        trend = t_start + (t_max - t_start) * (t / t.max()) ** 2 * rate * 100
        trend = t_start + (t_max - t_start) * np.clip((t / t.max()) ** 2.5, 0, 1)
        noise = np.random.normal(0, 0.4, n)
        seasonal = 0.3 * np.sin(2 * np.pi * t / 60)  # минутный период
        bearings[name] = trend + noise + seasonal
    return bearings


# Попытка загрузить реальные данные
if os.path.isdir('bearing-dataset'):
    bearings = load_bearing_data('bearing-dataset')
    print(f"Загружено {len(bearings)} подшипников из датасета FEMTO.")
else:
    bearings = generate_synthetic_bearings()

# Берём три подшипника для анализа
bearing_names = sorted(list(bearings.keys()))[:3]
print(f"Анализируем: {bearing_names}")

# ─── 1. Визуализация временных рядов ────────────────────────────────────────
def plot_timeseries():
    fig, axes = plt.subplots(len(bearing_names), 1, figsize=(14, 4 * len(bearing_names)), sharex=False)
    if len(bearing_names) == 1:
        axes = [axes]
    
    colors = ['steelblue', 'coral', 'mediumseagreen']
    
    for ax, name, color in zip(axes, bearing_names, colors):
        ts = bearings[name]
        t = np.arange(len(ts)) / 10  # секунды → мин/10
        ax.plot(t, ts, color=color, linewidth=0.8, alpha=0.9)
        ax.axhline(y=70, color='red', linestyle='--', linewidth=1.5,
                   label='Критический порог (70°C)')
        # Зона предотказного предупреждения
        warn_idx = np.where(ts > 60)[0]
        if len(warn_idx) > 0:
            warn_t = warn_idx[0] / 10
            ax.axvspan(warn_t, t[-1], alpha=0.12, color='orange',
                       label='Зона предупреждения (>60°C)')
        ax.set_title(f'Температурный ряд: {name}', fontweight='bold', fontsize=12)
        ax.set_ylabel('Температура (°C)')
        ax.set_xlabel('Время (×0.1 с)')
        ax.legend(loc='upper left', fontsize=9)
        ax.grid(alpha=0.3)
        # Метка максимума
        ax.annotate(f'Макс: {ts.max():.1f}°C', xy=(t[np.argmax(ts)], ts.max()),
                    xytext=(t[len(t)//2], ts.max() - 5),
                    arrowprops=dict(arrowstyle='->', color='red'), fontsize=10)
    
    plt.suptitle('Температурные временные ряды подшипников (FEMTO PHM 2012)',
                 fontsize=14, fontweight='bold', y=1.01)
    plt.tight_layout()
    plt.savefig('fig6_timeseries.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("Рисунок 6 сохранён: fig6_timeseries.png")

plot_timeseries()

# ─── 2. Статистический анализ ───────────────────────────────────────────────
print("\n" + "="*60)
print("ЭТАП 2. СТАТИСТИЧЕСКИЙ АНАЛИЗ")
print("="*60)
stats_data = []
for name in bearing_names:
    ts = bearings[name]
    stats_data.append({
        'Подшипник': name,
        'N': len(ts),
        'Среднее': f"{ts.mean():.2f}",
        'Медиана': f"{np.median(ts):.2f}",
        'СКО': f"{ts.std():.2f}",
        'Мин': f"{ts.min():.2f}",
        'Макс': f"{ts.max():.2f}",
        'Частота (Гц)': '10'
    })
df_stats = pd.DataFrame(stats_data)
print(df_stats.to_string(index=False))

# ─── 3. Диаграммы размаха (выбросы) ─────────────────────────────────────────
def plot_boxplots():
    fig, ax = plt.subplots(figsize=(10, 6))
    data = [bearings[n] for n in bearing_names]
    bp = ax.boxplot(data, labels=bearing_names, patch_artist=True,
                    medianprops=dict(color='red', linewidth=2))
    colors = ['lightsteelblue', 'lightsalmon', 'lightgreen']
    for patch, color in zip(bp['boxes'], colors):
        patch.set_facecolor(color)
    ax.axhline(y=70, color='red', linestyle='--', linewidth=1.5,
               label='Критический порог (70°C)')
    ax.set_title('Диаграммы размаха температурных рядов', fontweight='bold', fontsize=13)
    ax.set_ylabel('Температура (°C)')
    ax.set_xlabel('Подшипник')
    ax.legend()
    ax.grid(alpha=0.3, axis='y')
    plt.tight_layout()
    plt.savefig('fig7_boxplots.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("Рисунок 7 сохранён: fig7_boxplots.png")

plot_boxplots()

# ─── 4. Анализ пропусков ────────────────────────────────────────────────────
print("\n" + "="*60)
print("ЭТАП 3. АНАЛИЗ ПРОПУСКОВ И ВЫБРОСОВ")
print("="*60)
for name in bearing_names:
    ts = pd.Series(bearings[name])
    n_miss = ts.isnull().sum()
    mean, std = ts.mean(), ts.std()
    outliers_3sigma = ((ts < mean - 3*std) | (ts > mean + 3*std)).sum()
    print(f"  {name}: пропуски = {n_miss} ({n_miss/len(ts)*100:.2f}%), "
          f"выбросы (3σ) = {outliers_3sigma}")

# ─── 5. Корреляционный анализ ───────────────────────────────────────────────
def plot_correlation():
    # Выравниваем длины рядов до минимальной
    min_len = min(len(bearings[n]) for n in bearing_names)
    df_corr = pd.DataFrame({n: bearings[n][:min_len] for n in bearing_names})
    corr = df_corr.corr()
    
    fig, ax = plt.subplots(figsize=(7, 5))
    sns.heatmap(corr, annot=True, fmt='.3f', cmap='coolwarm', center=0,
                ax=ax, square=True, linewidths=0.5,
                xticklabels=bearing_names, yticklabels=bearing_names)
    ax.set_title('Матрица корреляций Пирсона между температурными каналами',
                 fontweight='bold', fontsize=11)
    plt.tight_layout()
    plt.savefig('fig8_correlation.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("Рисунок 8 сохранён: fig8_correlation.png")
    print(f"\nМатрица корреляций:\n{corr.round(3)}")

plot_correlation()

# ─── 6. Декомпозиция ряда и SNR ─────────────────────────────────────────────
def plot_decomposition():
    ts_name = bearing_names[0]
    ts = pd.Series(bearings[ts_name])
    
    # Сглаживаем незначительные пропуски
    ts = ts.interpolate()
    
    period = 60  # 60 отсчётов = 6 секунд при 10 Гц
    try:
        decomp = seasonal_decompose(ts, model='additive', period=period, extrapolate_trend='freq')
    except Exception:
        # Если ряд короткий — уменьшаем период
        period = 30
        decomp = seasonal_decompose(ts, model='additive', period=period, extrapolate_trend='freq')
    
    signal = decomp.trend + decomp.seasonal
    noise = decomp.resid.dropna()
    signal_clean = signal.dropna()
    
    snr = 10 * np.log10(np.var(signal_clean) / np.var(noise))
    print(f"\n{'='*60}")
    print(f"ДЕКОМПОЗИЦИЯ ({ts_name}): период сезонности = {period} отсч.")
    print(f"SNR = {snr:.1f} дБ")
    if snr > 20:
        snr_quality = "Отлично"
    elif snr > 10:
        snr_quality = "Хорошо"
    elif snr > 0:
        snr_quality = "Удовлетворительно"
    else:
        snr_quality = "Плохо"
    print(f"Качество сигнала: {snr_quality}")

    fig, axes = plt.subplots(4, 1, figsize=(14, 12), sharex=True)
    t = np.arange(len(ts))
    
    axes[0].plot(t, ts, color='steelblue', linewidth=0.8)
    axes[0].set_title('Исходный ряд', fontweight='bold')
    axes[0].set_ylabel('T (°C)')
    axes[0].grid(alpha=0.3)
    
    axes[1].plot(t, decomp.trend, color='darkorange', linewidth=1.5)
    axes[1].set_title('Тренд', fontweight='bold')
    axes[1].set_ylabel('T (°C)')
    axes[1].grid(alpha=0.3)
    
    axes[2].plot(t, decomp.seasonal, color='mediumseagreen', linewidth=0.8)
    axes[2].set_title('Сезонная компонента', fontweight='bold')
    axes[2].set_ylabel('ΔT (°C)')
    axes[2].grid(alpha=0.3)
    
    axes[3].plot(t, decomp.resid, color='tomato', linewidth=0.7, alpha=0.7)
    axes[3].axhline(0, color='black', linewidth=0.8)
    axes[3].set_title(f'Остатки (шум) | SNR = {snr:.1f} дБ [{snr_quality}]', fontweight='bold')
    axes[3].set_ylabel('Шум (°C)')
    axes[3].set_xlabel('Отсчёт')
    axes[3].grid(alpha=0.3)
    
    plt.suptitle(f'Аддитивная декомпозиция: {ts_name}', fontsize=14, y=1.01)
    plt.tight_layout()
    plt.savefig('fig9_decomposition.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("Рисунок 9 сохранён: fig9_decomposition.png")
    
    # Гистограмма остатков
    fig, ax = plt.subplots(figsize=(8, 5))
    ax.hist(noise.dropna(), bins=40, color='tomato', edgecolor='white', alpha=0.85, density=True)
    # Наложим нормальное распределение
    mu, sigma = noise.mean(), noise.std()
    x = np.linspace(noise.min(), noise.max(), 200)
    ax.plot(x, stats.norm.pdf(x, mu, sigma), 'k-', linewidth=2, label=f'Норм. N({mu:.2f}, {sigma:.2f}²)')
    ax.set_title('Гистограмма остатков (шум)', fontweight='bold')
    ax.set_xlabel('Значение остатка (°C)')
    ax.set_ylabel('Плотность')
    ax.legend()
    ax.grid(alpha=0.3)
    _, p_val = stats.shapiro(noise.dropna().sample(min(500, len(noise.dropna())), random_state=42))
    ax.text(0.02, 0.95, f'Тест Шапиро-Уилка: p={p_val:.4f}',
            transform=ax.transAxes, fontsize=10,
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    plt.tight_layout()
    plt.savefig('fig10_noise_hist.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("Рисунок 10 сохранён: fig10_noise_hist.png")

plot_decomposition()

# ─── 7. Градиент температуры (признак для модели) ───────────────────────────
print("\n" + "="*60)
print("СОЗДАНИЕ ПРИЗНАКОВ ДЛЯ МОДЕЛИ ПРЕДУПРЕЖДЕНИЯ")
print("="*60)
for name in bearing_names:
    ts = pd.Series(bearings[name])
    grad = ts.diff().fillna(0)                 # первая производная
    roll_mean = ts.rolling(window=60).mean()   # скользящее среднее (6 сек)
    roll_std = ts.rolling(window=60).std()     # скользящее СКО
    
    # Момент первого превышения порога warning
    warn_mask = ts > 60
    if warn_mask.any():
        warn_idx = warn_mask.idxmax()
        print(f"  {name}: предупреждение на отсчёте {warn_idx} "
              f"({warn_idx/10:.0f} с до начала опасной зоны)")
    else:
        print(f"  {name}: порог предупреждения не превышен")

print("\n✓ Анализ временных рядов завершён.")
print("  Файлы: fig6_timeseries.png, fig7_boxplots.png, fig8_correlation.png,")
print("         fig9_decomposition.png, fig10_noise_hist.png")