3. 塔架振动数据预处理与异常段截取

为什么原始振动时序不能直接做频谱

塔架振动在线监测系统通常以 50～200 Hz 连续采样，原始数据里同时混着四类成分：塔架结构响应、风轮转频及其谐波、传感器温漂与安装基线漂移、以及启停/偏航/变桨等瞬态过程。直接对整段数据做 FFT，会因趋势项造成低频能量泄漏、因瞬态冲击抬高宽带底噪、因转速变化使 1P/3P 谱峰展宽甚至分裂。因此预处理的目标很明确：保留结构响应与稳态转频成分，压制趋势、噪声与瞬态段，并按转速工况切出稳定运行段。

处理链与关键参数

本例采用一条可复现的处理链，参数按 100 Hz 采样、塔架一阶频率约 0.28 Hz 的典型海上风机设定（示例参数，需按实际机组调整）：

• 去趋势：对每段数据做一阶多项式最小二乘拟合后相减，消除温漂与基线漂移；对长记录可分段去趋势，段长取 60 s，避免把真实的低频结构响应一并减掉。
• 低通滤波：4 阶 Butterworth 低通，截止 20 Hz，零相位（前后向）滤波，避免相位延迟影响后续峰值时刻判读；20 Hz 已覆盖塔架前几阶与 3P 上限。
• 重采样：统一降到 50 Hz（满足 Nyquist 对 20 Hz 的 2.5 倍余量），便于多测点、多设备数据对齐。
• 工况截取：以 SCADA 10 min 平均转速为工况标签，取转速在目标工况 ±3% 内、且持续 ≥120 s 的区段；同时剔除发电机并网/脱网前后 60 s、偏航动作期间、变桨速率超限时段。

可复现的处理实现

下面给出一个自洽的最小实现，输入为一段示例振动序列与同步转速序列，输出为按工况切好的稳定段列表。示例数据用合成信号构造（含趋势、1P/3P、塔架一阶与启停瞬态），便于验证处理链是否正确，不代表任何实测结果。

import numpy as np
from scipy import signal

FS_RAW = 100.0          # 原始采样率 Hz
FS_OUT = 50.0           # 目标采样率 Hz
FC     = 20.0           # 低通截止 Hz
SEG_S  = 60.0           # 去趋势分段长度 s
MIN_LEN_S = 120.0       # 稳定段最短长度 s
RPM_TOL   = 0.03        # 转速容差 ±3%

def make_demo(n_sec=1800, seed=0):
    """合成示例：趋势 + 1P/3P + 塔架一阶 + 启停/偏航瞬态"""
    rng = np.random.default_rng(seed)
    t = np.arange(0, n_sec, 1/FS_RAW)
    # 转速：0-120s 启机，120-1500s 稳定 9.6rpm，1500-1560 偏航，1560-1800 稳定
    rpm = np.full_like(t, 9.6)
    rpm[t < 120] = np.linspace(0, 9.6, (t < 120).sum())
    rpm[(t >= 1500) & (t < 1560)] = 9.6 + 0.8*np.sin(2*np.pi*(t[(t>=1500)&(t<1560)]-1500)/30)
    f1p = rpm / 60.0
    # 塔架一阶 0.28 Hz，1P≈0.16 Hz，3P≈0.48 Hz
    x = (0.02*t                                  # 温漂趋势
         + 0.8*np.sin(2*np.pi*0.28*t)            # 塔架一阶
         + 0.5*np.sin(2*np.pi*f1p*t)             # 1P
         + 0.3*np.sin(2*np.pi*3*f1p*t)           # 3P
         + 0.05*rng.standard_normal(t.size))     # 宽带噪声
    # 启停瞬态冲击
    x[t < 120] += 1.5*np.exp(-((t-60)/8)**2)
    return t, x, rpm

def detrend_seg(x, fs, seg_s=SEG_S):
    """分段一阶多项式去趋势"""
    n = int(seg_s*fs)
    y = x.copy()
    for i in range(0, len(x), n):
        seg = x[i:i+n]
        if len(seg) < 10:
            continue
        k = np.arange(len(seg))
        p = np.polyfit(k, seg, 1)
        y[i:i+n] = seg - np.polyval(p, k)
    return y

def lowpass_zero_phase(x, fs, fc=FC, order=4):
    b, a = signal.butter(order, fc/(fs/2), btype='low')
    return signal.filtfilt(b, a, x)

def resample(x, fs_in, fs_out):
    g = np.gcd(int(fs_in), int(fs_out))
    return signal.resample_poly(x, int(fs_out)//g, int(fs_in)//g)

def extract_stable(t, x, rpm, fs, min_len_s=MIN_LEN_S, tol=RPM_TOL):
    """按转速工况截取稳定段，剔除启停/偏航"""
    rpm_mean = np.mean(rpm[rpm > 0.5]) if np.any(rpm > 0.5) else 0
    if rpm_mean == 0:
        return []
    in_band = np.abs(rpm - rpm_mean) / rpm_mean <= tol
    # 剔除转速变化率过大的点（启停/偏航）
    drpm = np.abs(np.gradient(rpm, t))
    smooth = in_band & (drpm < 0.02)
    segs, start = [], None
    for i, ok in enumerate(smooth):
        if ok and start is None:
            start = i
        elif not ok and start is not None:
            if (i - start)/fs >= min_len_s:
                segs.append((start, i))
            start = None
    if start is not None and (len(smooth)-start)/fs >= min_len_s:
        segs.append((start, len(smooth)))
    return [(t[a], t[b-1], x[a:b]) for a, b in segs]

# ---- 执行 ----
t, x, rpm = make_demo()
x_dt = detrend_seg(x, FS_RAW)
x_lp = lowpass_zero_phase(x_dt, FS_RAW)
x_rs = resample(x_lp, FS_RAW, FS_OUT)
t_rs = np.arange(len(x_rs)) / FS_OUT
rpm_rs = resample(rpm, FS_RAW, FS_OUT)

segs = extract_stable(t_rs, x_rs, rpm_rs, FS_OUT)
print(f"原始时长 {len(x)/FS_RAW:.0f}s，输出采样率 {FS_OUT}Hz，稳定段数 {len(segs)}")
for i, (a, b, s) in enumerate(segs):
    print(f"  段{i+1}: {a:.1f}s ~ {b:.1f}s, 长度 {b-a:.1f}s, 均值 {s.mean():.4f}, 标准差 {s.std():.4f}")

预期输出形态（示例数据下的合理结果，非实测）：稳定段数 2，分别对应启机后与偏航后的两段稳态，段长均超过 120 s；去趋势后段内均值接近 0，标准差主要由 0.28 Hz 塔架一阶与 1P/3P 贡献。若稳定段数为 0 或段长不足，优先检查转速容差与最短段长是否与机组实际转速波动匹配。

结果判读与常见坑