2. 轨道不平顺样本的异常值识别与剔除

异常类型与识别思路

实测轨道不平顺序列中常见的异常有三类：尖峰（单点或极少数点幅值远超邻域，多由轨面异物、激光测头瞬时失锁引起）、丢帧（采样点缺失或保持为同一值，表现为局部方差骤降或重复值）、传感器跳变（相邻点差分出现阶跃，之后基线整体平移，多由传感器重新锁定或里程计滑移引起）。三类异常在时域上的形态不同，因此需要组合两类判据：

• 统计阈值判据：以滑动窗口中位数与 MAD（中位数绝对偏差）构造稳健 z 分数，识别偏离局部主体的孤立点，对尖峰最敏感。
• 局部差分判据：计算一阶差分 d[i]=x[i]−x[i−1]，用差分的中位数绝对偏差设定门限，识别阶跃型跳变；同时用连续相同值计数识别丢帧保持段。

关键约束是：轨道不平顺本身含有真实的高频成分（短波不平顺），阈值过松会漏检，过紧会把真实短波当成尖峰削掉，从而在后续 PSD 高频段造成能量损失。因此剔除策略采用“只替换被判为异常的点、保留其余原始值”，并在剔除前后对比时域波形与差分序列，确认真实成分未被破坏。

可运行的识别与剔除实现

下面给出一个自洽的 Python 实现。输入为已按里程排序、等空间间隔采样的单通道高低不平顺序列（单位 mm，空间采样间隔 0.25 m）。示例数据由“低频真实不平顺 + 真实短波 + 人工注入的三类异常”合成，便于验证算法是否只命中注入点。

import numpy as np

rng = np.random.default_rng(20240501)

# ---------- 1. 构造示例序列：真实成分 + 注入异常 ----------
dx = 0.25                      # 空间采样间隔, m
n = 4000
s = np.arange(n) * dx

# 真实长波 + 短波不平顺（模拟轨道高低, mm）
truth = (1.8 * np.sin(2 * np.pi * s / 40.0)
         + 0.6 * np.sin(2 * np.pi * s / 6.0)
         + 0.25 * np.sin(2 * np.pi * s / 1.2)
         + rng.normal(0, 0.08, n))

x = truth.copy()
spike_idx = np.array([350, 351, 1200, 2600, 2601, 2602])
x[spike_idx] += np.array([6.5, 5.8, -7.2, 8.0, -6.4, 7.1])

drop_idx = np.arange(1800, 1812)          # 丢帧：保持上一有效值
x[drop_idx] = x[1800 - 1]

jump_at = 3200                            # 传感器跳变：阶跃 + 基线平移
x[jump_at:] += 3.0

# ---------- 2. 稳健统计阈值：滑动中位数 + MAD ----------
def rolling_median(a, win):
    half = win // 2
    pad = np.pad(a, half, mode="reflect")
    return np.array([np.median(pad[i:i + win]) for i in range(len(a))])

def mad(a):
    m = np.median(a)
    return np.median(np.abs(a - m))

win = 41                                  # 约 10 m 窗口
med = rolling_median(x, win)
resid = x - med
sigma_robust = 1.4826 * mad(resid)        # 稳健标准差估计
z = np.abs(resid) / max(sigma_robust, 1e-9)
mask_spike = z > 6.0                      # 统计阈值判据

# ---------- 3. 局部差分判据：阶跃跳变 ----------
d = np.diff(x, prepend=x[0])
d_sigma = 1.4826 * mad(d)
mask_jump = np.abs(d) > 8.0 * max(d_sigma, 1e-9)

# ---------- 4. 丢帧判据：连续相同值 ----------
same = np.zeros(n, dtype=bool)
run = 1
for i in range(1, n):
    run = run + 1 if x[i] == x[i - 1] else 1
    if run >= 3:                          # 连续 3 点及以上完全相同
        same[i - run + 1:i + 1] = True

mask_bad = mask_spike | mask_jump | same

# ---------- 5. 剔除：异常点用局部中位数替换 ----------
x_clean = x.copy()
x_clean[mask_bad] = med[mask_bad]

print("注入尖峰点数:", len(spike_idx))
print("统计阈值命中:", int(mask_spike.sum()))
print("差分判据命中:", int(mask_jump.sum()))
print("丢帧判据命中:", int(same.sum()))
print("总剔除点数:", int(mask_bad.sum()))
print("剔除前后 RMS 变化: %.4f -> %.4f mm"
      % (np.sqrt(np.mean((x - truth) ** 2)),
         np.sqrt(np.mean((x_clean - truth) ** 2))))

代码中三个判据各自独立，最后取并集。替换值使用滑动中位数而非全局均值，避免在长波起伏段引入额外偏差。运行后应观察到：统计阈值判据主要命中注入的尖峰点，差分判据命中跳变起始点，丢帧判据命中保持段；总剔除点数略大于注入点数属于正常现象（跳变点附近的差分也会超限）。

剔除前后的时域对比与判读

下图给出同一段序列在剔除前后的时域波形与一阶差分序列对照，用于判断异常是否被正确清除、真实短波是否保留。

判读要点：

• 尖峰在原始差分中表现为一对方向相反的大幅跳变（先上冲后回落），剔除后差分回到门限内；
• 丢帧保持段在原始序列中表现为一段水平直线，剔除后恢复为与邻域一致的起伏；
• 传感器跳变在差分中只有一个方向的阶跃，剔除后基线平移被消除，但跳变点之后的真实起伏形态保持不变；
• 若剔除后差分序列出现成片超门限的振荡，说明阈值过紧，真实短波被误判，应放宽 z 门限或增大滑动窗口。

参数选择与验证清单