3. 轨道不平顺去趋势与基线漂移修正

问题定位：长波基线漂移如何污染PSD低频段

轨道检测车输出的高低不平顺序列中，除真实轨道不平顺外，还叠加了由车体姿态、传感器零漂、轨道整体坡度与竖曲线引起的长波基线分量。这类分量在时域上表现为缓慢起伏，在频域上则集中到极低空间频率（通常对应波长数十米至数百米），会显著抬高PSD低频段幅值，甚至掩盖真实的长波不平顺。本例对一段含基线漂移的高低序列分别做1阶、2阶、3阶多项式去趋势，并叠加Butterworth高通滤波，比较各方案对低频段PSD的影响，最终确定预处理方案。

数据与参数设定

取一段空间采样间隔 Δx = 0.25 m、长度 2000 m 的高低不平顺序列（示例数据，由真实长波趋势叠加随机不平顺合成，用于演示方法）。空间采样频率 fs = 1/Δx = 4 1/m，对应 Nyquist 空间频率 2 1/m。高通滤波采用零相位（filtfilt）Butterworth，截止波长 λc = 100 m，即截止空间频率 fc = 1/λc = 0.01 1/m，滤波器阶数取 2 阶（双向滤波后等效 4 阶）。

去趋势阶次选择依据：轨道整体坡度对应1阶（线性）趋势；竖曲线与长波沉降对应2阶（抛物线）趋势；3阶及以上容易把真实的长波不平顺一并拟合掉，造成低频段“过度修正”。

关键实现（Python）

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

# 示例：合成含基线漂移的高低序列
dx = 0.25                      # 空间采样间隔 (m)
x = np.arange(0, 2000, dx)     # 里程 (m)
rng = np.random.default_rng(42)
# 真实不平顺：多尺度随机成分
true_irr = 0.8*np.sin(2*np.pi*x/60) + 0.3*rng.standard_normal(x.size)
# 基线漂移：线性坡度 + 抛物线竖曲线
drift = 0.004*x + 2e-6*(x-1000)**2
signal = true_irr + drift

def detrend_poly(y, x, order):
    coef = np.polyfit(x, y, order)
    return y - np.polyval(coef, x)

def highpass_spatial(y, dx, lam_c=100.0, order=2):
    fc = 1.0/lam_c              # 截止空间频率 (1/m)
    fs = 1.0/dx                 # 空间采样频率 (1/m)
    b, a = butter(order, fc/(fs/2), btype='high')
    return filtfilt(b, a, y)

results = {}
for n in (1, 2, 3):
    y_dt = detrend_poly(signal, x, n)
    y_hp = highpass_spatial(y_dt, dx, lam_c=100.0)
    results[n] = y_hp

# 低频段PSD对比（Welch，空间频率）
from scipy.signal import welch
for n, y in results.items():
    f, Pxx = welch(y, fs=1/dx, nperseg=1024)
    low = f < 0.02
    print(f"n={n}, 低频段(<0.02 1/m)平均PSD={Pxx[low].mean():.3e}")

结果判读与方案确定

运行上述代码后，比较三种阶次在低频段（空间频率 < 0.02 1/m，对应波长 > 50 m）的平均PSD：1阶去趋势后低频段仍保留明显翘起，说明线性趋势不足以描述竖曲线引起的抛物线漂移；2阶去趋势后低频段幅值显著下降且曲线趋于平坦；3阶去趋势后低频段进一步下降，但需检查是否出现负值或振荡伪峰——若出现，说明拟合过度，削掉了真实长波成分。

综合判断：本例数据适合采用 2阶多项式去趋势 + λc=100 m 零相位Butterworth高通 的预处理方案。验证方法为：对比处理前后时域波形，确认基线漂移被移除而真实不平顺形态保留；对比低频段PSD，确认幅值下降且无异常伪峰；检查处理序列无负值截断或端点振荡。