4、速度与加速度噪声处理：从速度轨迹提取与滤波，及加速度计算中的噪声抑制

速度与加速度是分子动力学轨迹分析里的高频信息。位置轨迹本身就有噪声，一阶差分等于把噪声放大，二阶差分更是灾难。这一章就聊聊怎么从位置轨迹里合理地提取速度和加速度，以及怎样在计算过程中把噪声压下去。

先把核心逻辑理清楚，我画了张图帮你建立整体感觉——咱们后面一步步展开。

速度与加速度噪声处理流程 原始位置轨迹 (含噪声) 中心差分 速度轨迹 (噪声放大) 滤波后速度 (SG/滑动平均) 加速度轨迹 (二阶差分) 干净速度 滤波后 关键原则 • 先滤波再差分 • 保留原始轨迹备份 • 注意时间步一致

4.1 速度轨迹的提取与噪声来源

最直接的提取方法就是中心差分。假设位置为 r(t)，速度为：

v(t) = [r(t+Δt) - r(t-Δt)] / (2Δt)

为什么用中心差分？因为它的截断误差是 O(Δt²)，比向前差分 O(Δt) 好很多。但这里有个麻烦——位置数据的噪声被差分操作放大了。

举个例子：如果位置有 0.01 Å0.01 / Δt。假设 Δt = 0.001 ps，速度误差就是 10 Å/ps。这个量级对于真实速度来说根本不是噪声，而是直接把信号淹没了。我在处理一个蛋白质体系时就遇到这种情况，算出来的速度曲线看起来完全随机，根本没法用。

简单列一下噪声来源：

  • 位置采样噪声：浮点舍入、积分器误差、输出精度截断
  • 时间步长不均匀：有些MD输出帧的时间间隔不一致（尤其并行写入时）
  • 边界跨跃：原子穿过周期边界导致位置跳变，差分时会出离谱值

处理跨周期边界是我的老经验了。记得有一次我忘了unwrap轨迹就直接算速度，结果某个原子在边界处速度达到了 几百 Å/ps，笑死。后来我在轨迹加载时就先做 unwrap，确保位置连续。

4.2 速度滤波策略

滤除速度噪声的方法很多，我个人最常用三个：滑动平均、Savitzky–Golay（SG）滤波，以及基于FFT的低通滤波。下面这张表帮你快速对比。

滤波方法 优势 劣势 适用场景
滑动平均 实现简单，计算快 降低峰值幅度，信号变钝 定性观察趋势，不要求精确峰值
Savitzky–Golay 保留峰值形状，噪声抑制好 参数（窗口/阶数）需调试 定量分析需要保留速度峰高
FFT低通滤波 可精确控制截止频率 可能产生振铃伪影 已知噪声频段，体系周期长

滑动平均。窗口大小为 w，输出为周围 w 个点的平均值。窗口越大越平滑，但也会抹掉真实的速度变化。我对窗宽的经验：w ≈ 5～15，视轨迹总帧数和原子运动快慢而定。

Savitzky–Golay 滤波。这个我特别推荐。它用局部多项式拟合来平滑，对峰值形状保留得非常好。举个例子，速度谱中尖锐的振动模式用滑动平均就压平了，但SG滤波能恢复出干净的峰。

参数怎么调？窗口长度选奇数，多项式阶数一般用2或3。我习惯先试 window=7, order=2，效果不好再调。阶数太高反而会过拟合噪声。

FFT低通滤波。如果你的体系存在明显的特征频率（比如键伸缩 > 3000 cm⁻¹），用FFT直接把高频部分砍掉最直接。但要注意——截断过陡会在时域产生振铃。我建议用巴特沃斯或高斯滤波器过渡。

4.3 加速度计算的挑战

加速度 = 速度的一阶差分，也就是位置的二阶差分：

a(t) = [v(t+Δt) - v(t-Δt)] / (2Δt)
     = [r(t+Δt) - 2r(t) + r(t-Δt)] / (Δt²)

看看分母，Δt 是平方项——噪声放大效应比速度还猛一个数量级。我见过有人直接对原始位置算二阶差分，出来的加速度图全是噪声，根本看不出任何物理意义。

所以我的核心原则是：先滤波，再差分。

具体做法有两种路径：

  1. 路径A：原始位置轨迹 → 中心差分得速度 → 对速度做SG滤波 → 再用滤波后的速度差分得加速度
  2. 路径B：原始位置轨迹 → 先对位置做SG滤波 → 再从平滑位置差分速度和加速度

两种我都试过。路径A更保险，因为速度噪声相对可控，而且滤波参数可以直接用速度数据来调试。路径B如果位置滤波窗口太大，会把真实的高频运动也平滑掉。

哪种更好？取决于你的体系和目标。我做溶剂扩散系数分析时通常用路径A；分析分子内振动模态时则用路径B但严格控制窗口。

4.4 实战：Python 代码实现

下面是一套完整的实现，从轨迹读取到速度提取、SG滤波、加速度计算，一步到位。

import numpy as np
from scipy.signal import savgol_filter

def extract_vel_acc(pos, dt, sg_window=7, sg_order=2):
    """
    pos: (N_frame, N_atom, 3) 位置轨迹，已unwrap
    dt: 时间步长 (ps)
    返回: vel, acc, vel_filtered
    """
    # 1. 中心差分得速度
    vel = np.zeros_like(pos)
    vel[1:-1] = (pos[2:] - pos[:-2]) / (2 * dt)
    # 边界用向前/向后差分
    vel[0]  = (pos[1] - pos[0]) / dt
    vel[-1] = (pos[-1] - pos[-2]) / dt

    # 2. 对每个原子每个维度做SG滤波
    vel_filt = np.zeros_like(vel)
    for i in range(pos.shape[1]):
        for d in range(3):
            vel_filt[:, i, d] = savgol_filter(
                vel[:, i, d], sg_window, sg_order
            )

    # 3. 从滤波后速度算加速度
    acc = np.zeros_like(vel_filt)
    acc[1:-1] = (vel_filt[2:] - vel_filt[:-2]) / (2 * dt)
    acc[0]  = (vel_filt[1] - vel_filt[0]) / dt
    acc[-1] = (vel_filt[-1] - vel_filt[-2]) / dt

    return vel, acc, vel_filt

# 使用示例
dt = 0.001  # ps
pos = np.random.randn(1000, 500, 3)  # 模拟轨迹
vel, acc, vel_filt = extract_vel_acc(pos, dt, sg_window=9, sg_order=2)

注意代码中边界帧的处理。我提一句——边界帧的精度天然较低，我通常会在后续分析中剔除头尾各2～3帧。

4.5 避坑指南

坑1：忘记统一时间步长

MD输出有时不是严格等间隔的。特别是NPT系综下，脚本自己dump帧，可能因为重采样或写入延迟导致Δt不固定。如果直接当成等步长用，差分值会出错。我在处理GROMACS轨迹时遇到过这类问题——解决办法是在读入帧时核验时间数组的时间间隔，有偏差就插值重采样。

坑2：先差分再滤波 vs 先滤波再差分

一定要记住先滤波再差分。我曾经为了省事直接对位置做一次SG滤波后算加速度，但窗口开太大把真实信号滤掉了，分析结果全偏。现在我做任何处理都会保留一份原始轨迹的副本。

坑3：跨周期边界造成的数据异常

没有unwrap的轨迹，差分时会出现巨大的速度尖峰。解决方案：在读取轨迹时就调用 mdtraj 或 MDAnalysis 的unwrap函数，确保所有原子位置在周期边界内是连续的。

个人小技巧：先画出速度曲线看一眼

在输入滤波参数之前，先随机挑几个原子的速度曲线画出来。看看整体量级、噪声形状，大致确定SG的窗口宽度。这比直接套默认参数靠谱很多。

嗯，速度与加速度的噪声处理核心就这些。下一个关键点是怎么把处理后的速度场用来分析扩散、流体力学等宏观性质。不过那是后面章节的事了，这里先把基础打牢。

← 上一章 📖 返回目录 下一章 →