4. 速度信号微分求加速度的滤波参数选择

为什么不能直接对速度做一阶差分

城轨牵引试验中，速度通道通常来自速度传感器或牵引逆变器估算转速，采样率常见为 100 Hz～1 kHz，且本身带有量化误差与随机噪声。若直接做一阶差分 a[k] = (v[k] − v[k−1]) / Δt，噪声会被放大 1/Δt 倍：以 200 Hz 采样、速度噪声 0.01 m/s 为例，单点差分噪声可达 2 m/s²，与城轨牵引加速度量级（约 0.8～1.2 m/s²）同量级，完全不可用。因此必须同时做两件事：加大差分步长（降噪但引入相位滞后与平滑失真）和低通滤波（进一步压制高频噪声，但同样带来相位滞后）。本案例的任务就是在这两个参数之间找到适用于牵引试验的折中组合。

差分步长与截止频率的作用机理

设速度序列 v[n]，采样间隔 Ts。中心差分步长 m（前后各取 m 点）的加速度估计为：

a[n] = (v[n+m] − v[n−m]) / (2·m·Ts)

该差分器等效为一个长度为 2m+1 的滑动平均与微分复合，其幅频特性在低频近似为 ω，但会抑制高频；同时群延迟约为 m·Ts。随后串联一阶低通（或二阶 Butterworth）滤波，截止频率 fc 决定残余噪声与相位滞后的平衡。工程上常用二阶 Butterworth 低通，其相位滞后在 fc 附近约为 90°/级，在远低于 fc 的牵引工况频段（0～2 Hz）滞后较小。

用实测速度序列做参数扫描

下面给出一段可直接运行的 Python 实现：构造一段带噪声的城轨牵引速度序列（0～3 s 加速，加速度约 1.0 m/s²，叠加 0.01 m/s 高斯噪声与 0.005 m/s 量化），然后对 m ∈ {1,3,5,10,20} 与 fc ∈ {2,5,10,20} Hz 做网格扫描，统计稳态段加速度标准差（噪声指标）与阶跃响应上升时间（滞后指标），最后按“噪声 ≤ 0.05 m/s² 且总滞后 ≤ 50 ms”的约束筛选推荐组合。

import numpy as np
from scipy.signal import butter, filtfilt

fs = 200.0
Ts = 1.0 / fs
t = np.arange(0, 6.0, Ts)

# 构造示例速度：0~3 s 匀加速 1.0 m/s²，之后匀速
a_true = np.where(t < 3.0, 1.0, 0.0)
v_true = np.cumsum(a_true) * Ts
rng = np.random.default_rng(42)
v_meas = np.round(v_true + rng.normal(0, 0.01, t.size), 3)  # 量化到 1 mm/s

def diff_acc(v, m):
    a = np.full_like(v, np.nan)
    a[m:-m] = (v[2*m:] - v[:-2*m]) / (2 * m * Ts)
    return a

def lowpass(a, fc, order=2):
    b, aa = butter(order, fc / (fs / 2), btype='low')
    return filtfilt(b, aa, a)

results = []
for m in [1, 3, 5, 10, 20]:
    a_raw = diff_acc(v_meas, m)
    for fc in [2, 5, 10, 20]:
        a_f = lowpass(np.nan_to_num(a_raw), fc)
        steady = a_f[int(4.0/Ts):int(5.5/Ts)]          # 匀速段
        noise_std = np.std(steady)
        # 阶跃响应：加速度从 0 到 1 的 10%~90% 上升时间
        seg = a_f[int(2.8/Ts):int(3.4/Ts)]
        t_seg = t[int(2.8/Ts):int(3.4/Ts)]
        i10 = np.argmax(seg >= 0.1); i90 = np.argmax(seg >= 0.9)
        rise = t_seg[i90] - t_seg[i10]
        results.append((m, fc, noise_std, rise))

print(f"{'m':>3} {'fc/Hz':>6} {'噪声σ/(m/s²)':>14} {'上升时间/ms':>12}")
for m, fc, s, r in results:
    print(f"{m:>3} {fc:>6.0f} {s:>14.4f} {r*1000:>12.1f}")

# 约束筛选：噪声 ≤ 0.05 m/s²，上升时间 ≤ 50 ms
ok = [(m, fc, s, r) for m, fc, s, r in results if s <= 0.05 and r <= 0.05]
print("\n满足约束的组合：", ok)

示例运行会得到类似下表的结果（数值为示例数据下的预期量级，实际以现场数据为准）：

选择依据与落地检查

推荐组合的判据不是“噪声最小”或“滞后最小”，而是在滞后容差内取噪声最小。城轨牵引试验中，工况切分与冲击率计算的时刻误差通常要求不超过 50 ms，因此把总滞后（差分群延迟 + 滤波器相位滞后）控制在 30～40 ms 以内是合理的工程目标。对 200 Hz 采样，m = 5 对应 25 ms 群延迟，配合 fc = 5 Hz 二阶 Butterworth 在 2 Hz 牵引频段附加约 15～20 ms 滞后，总滞后约 40 ms，同时把噪声压到 0.05 m/s² 以下，是较稳妥的起点。

落地时还需做三项检查：

• 边界处理：中心差分在序列首尾各 m 点无定义，需用前向/后向差分补齐或直接截断，并在结果中标注有效区间，避免把边界 NaN 当作异常值。
• 零相位滤波的取舍：filtfilt 是零相位滤波，不引入相位滞后，但会使用未来数据，不能用于在线判定；离线分析可用 filtfilt，在线或实时判定必须用 lfilter 并显式补偿滞后。
• 参数随采样率缩放：若采样率从 200 Hz 变为 1 kHz，同样的物理滞后对应 m 放大 5 倍，fc 保持不变（物理频率），不能直接沿用 m 的数值。