📚 课程目录第3章船舶旋转机械声振信号预处理与同步对齐

3. 船舶旋转机械声振信号预处理与同步对齐

预处理链路与同步对齐的核心问题

振动通道与声压通道往往来自不同传感器、不同调理链路，甚至不同采样率。直接做互相关或联合谱分析时，会同时叠加三类误差：趋势项（温漂、直流偏置、积分漂移）、频带差异（振动关注高频结构共振，声压关注中低频辐射噪声）、时间基准差异（采样率不同、触发延迟不同）。本节以一台船用齿轮箱试验台为例，给出可复现的预处理与对齐流程。

示例条件（教学假设）：振动通道 fs=51200 Hz，声压通道 fs=48000 Hz，采集时长 10 s，两通道由同一触发源启动但存在约 1.7 ms 的固定延迟；齿轮箱输入轴转频约 12.5 Hz（750 r/min）。目标是把两路信号统一到 51200 Hz、零均值、带限一致，并把通道间残余延迟压到 1 个采样点以内。

关键实现：从原始双通道到对齐矩阵

下面给出可直接运行的 Python 实现，使用 NumPy 与 SciPy，覆盖去趋势、零相位带限滤波、有理数重采样、互相关时延估计与分数延迟对齐。示例数据由已知延迟的合成信号构造，便于验证对齐结果是否收敛到设定值。

import numpy as np
from scipy import signal

fs_v, fs_a = 51200, 48000
dur = 10.0
true_delay = 1.7e-3          # 声压通道相对振动的真实延迟（教学假设）

rng = np.random.default_rng(7)
n_v = int(fs_v * dur)
n_a = int(fs_a * dur)
t_v = np.arange(n_v) / fs_v
t_a = np.arange(n_a) / fs_a

# 合成：齿轮啮合分量 + 宽带随机 + 慢变趋势
f_mesh = 12.5 * 23           # 输入轴转频 12.5 Hz，齿数 23
vib = (0.8*np.sin(2*np.pi*f_mesh*t_v)
       + 0.3*rng.standard_normal(n_v)
       + 0.05*t_v)            # 人为加入线性趋势
ac_base = (0.5*np.sin(2*np.pi*f_mesh*(t_a - true_delay))
           + 0.25*rng.standard_normal(n_a)
           + 0.02*t_a)
ac = ac_base

def detrend_poly(x, order=2):
    t = np.arange(len(x))
    coef = np.polyfit(t, x, order)
    return x - np.polyval(coef, t)

def bandpass_zero_phase(x, fs, lo=20.0, hi=8000.0, numtaps=1025):
    hi = min(hi, 0.45 * fs)
    taps = signal.firwin(numtaps, [lo, hi], pass_zero=False, fs=fs)
    return signal.filtfilt(taps, [1.0], x)

vib_d = detrend_poly(vib)
ac_d  = detrend_poly(ac)

vib_f = bandpass_zero_phase(vib_d, fs_v)
ac_f  = bandpass_zero_phase(ac_d,  fs_a)

# 声压 48k -> 51.2k，有理数比 16/15
ac_rs = signal.resample_poly(ac_f, up=16, down=15)

# 截取等长公共段
n = min(len(vib_f), len(ac_rs))
vib_f, ac_rs = vib_f[:n], ac_rs[:n]

# 互相关估计时延（以振动为参考）
max_lag = int(0.02 * fs_v)   # 搜索 ±20 ms
corr = signal.correlate(ac_rs, vib_f, mode="full")
lags = signal.correlation_lags(len(ac_rs), len(vib_f), mode="full")
mask = np.abs(lags) <= max_lag
lag = lags[mask][np.argmax(corr[mask])]
peak = corr[mask].max() / (np.linalg.norm(ac_rs)*np.linalg.norm(vib_f))
delay_est = lag / fs_v

# 整数移位 + 分数延迟插值对齐
shift = -lag
if shift >= 0:
    ac_al = np.concatenate([np.zeros(shift), ac_rs])[:n]
else:
    ac_al = np.concatenate([ac_rs[-shift:], np.zeros(-shift)])[:n]

frac = delay_est - np.round(delay_est)
if abs(frac) > 1e-6:
    ac_al = signal.resample_poly(ac_al, 1000, 1000)  # 占位，实际用分数延迟滤波器
    # 更稳妥：用 FFT 相位斜坡实现分数延迟
    N = len(ac_al)
    f = np.fft.rfftfreq(N, d=1/fs_v)
    spec = np.fft.rfft(ac_al)
    ac_al = np.fft.irfft(spec * np.exp(-2j*np.pi*f*frac), n=N)

print(f"估计延迟 = {delay_est*1e3:.3f} ms, 真实延迟 = {true_delay*1e3:.3f} ms")
print(f"互相关归一化峰值 = {peak:.3f}")

运行后应观察到：估计延迟与 1.7 ms 接近（误差在 1 个采样点即约 0.0195 ms 量级内），归一化互相关峰值明显高于 0.6。若峰值偏低，通常说明带限不一致或趋势未去除干净，应先回到第②步检查通带设置。

参数选择与常见陷阱

环节	关键参数	取值依据	典型陷阱
去趋势	多项式阶次	2 阶适合温漂与缓慢基线漂移	阶次过高会吞掉低频故障分量
带限滤波	通带 20–8000 Hz	覆盖齿轮啮合与结构共振，避开电源与极低频	用单向滤波引入相位偏移，破坏对齐
重采样	16/15 有理数比	48000 与 51200 的最小整数比	先重采样后滤波导致混叠
时延估计	搜索窗 ±20 ms	覆盖触发链路与传播延迟	搜索窗过窄导致峰值落在边界
分数延迟	FFT 相位斜坡	亚采样精度对齐	直接四舍五入到整数样本，残余延迟超标

重采样必须在抗混叠滤波之后进行。若声压通道原始带宽接近 24 kHz，直接 48k→51.2k 上采样虽不产生混叠，但若后续还要降采样，则必须先低通到新奈奎斯特频率以下。此外，两通道截取公共段时应以触发基准为锚点，而不是各自从数组第 0 点开始，否则会引入整段偏移。

对齐结果的验证方法

对齐是否成功不能只看代码是否跑通，应做三项检查：其一，对齐后两通道互相关的归一化峰值应显著高于对齐前；其二，用已知延迟的合成信号做闭环验证，估计值与设定值之差应小于 1 个采样周期；其三，检查对齐后两通道的相干函数在齿轮啮合频率处是否接近 1，若明显偏低则说明残余延迟或频带不一致仍存在。三项检查全部通过，才可把数据交给后续的联合特征提取环节。