★ FEATURED ARTICLE 简介一份围绕位移、速度与加速度信号互转的Matlab配套代码与数据包适合信号处理、振动分析及工程数据预处理场景帮助学习者把微分积分理论落实到可运行脚本中。资源共4个文件包含3个m脚本与1个mat数据文件压缩包约23KB。脚本覆盖角速度与角位移转换、正弦信号生成与测试、加速度与速度之间换算等典型操作数据文件提供示例采样数据便于直接加载验证处理流程。已有6120人学习下载。通过阅读脚本可掌握diff、cumsum等函数在信号求导与积分中的实际运用并了解滤波去噪、数据导入、结果绘图与保存等完整处理思路是巩固物理量与Matlab编程结合能力的高性价比实践案例。1. 位移、速度、加速度信号互转先避开第一个坑搞振动测试、结构健康监测或者工业设备状态检测的同行几乎都会撞上同一个问题手里只有加速度传感器但甲方要的是速度有效值或者位移峰峰值。做模态测试时想从加速度积分出位移却总被趋势项和低频噪声折磨得怀疑人生。反过来手里有位移信号想求速度微分之后发现高频噪声被放大了几个数量级根本没法看。位移、速度、加速度这三个物理量在数学上就是积分和微分的关系方向看似简单实际做起来却是一堆玄学积分漂移、端点失真、相位偏差、频域泄漏每一个都能把结果搞翻车。这篇笔记会把时域积分、频域积分和微分滤波的完整路径拆开讲清楚每一步附能直接跑的代码和参数设置中间章重点写避坑最后给一个用合成信号验证精度的闭环方法方便你验证自己手头代码到底靠不靠谱。2. 位移、速度、加速度三者之间的本质关系积分与微分的两条误差放大路径2.1 从加速度到位移要过两次积分误差按时间指数增长物理关系本身不复杂加速度对时间做一次积分得到速度再积分一次得到位移。反过来说位移对时间求一阶导得到速度二阶导得到加速度。但在离散采样系统里这三条路径的误差特性完全不同。积分操作的误差主要来自两个方面。第一是低频漂移加速度传感器本身存在零偏即使静止状态下输出的也不是零而是一个小幅直流量。这个直流成分经过一次积分会变成一条向上或向下的斜线经过两次积分会变成一条抛物线。很多入门者看到二次积分出来的位移是条抛物线第一反应是程序写错了其实这是直流偏置和低频噪声累积的必然结果。第二是数值积分的累积误差。梯形法、辛普森法都是基于离散点之间的面积近似采样率越高、单次积分步长越小累积误差越小但误差会随着积分点数增加而积累。这里有个常见误区提高采样率能有效减小单步误差却不能直接消除累积漂移真正要处理的是信号里的低频成分和直流偏置。所以在做任何积分之前先做高通滤波或趋势项去除是业内公认的标准前置操作。2.2 反向微分路径的噪声放大比积分更让人头疼微分是积分逆运算但它的噪声特性正好相反。理想情况下一个正弦位移信号求导后是余弦速度信号幅值乘以角频率ω。这意味着每求一次导高频分量的幅值就被放大ω倍。如果信号里混入高频噪声而噪声往往能量集中在高频段微分之后噪声的幅值会被放大到噪声淹没真实信号的程度。实际工程中我一般不会直接对原始信号做差分来求速度因为一阶差分相当于差分滤波器幅频响应不是平坦的它在低频段接近真实导数在高频段会出现严重失真。更常见的做法是用频域微分或平滑差分滤波器也就是先对信号做低通滤波再求导或者在频域里做一次高通积分器的逆运算。这里有个反直觉的事实对位移信号求速度优先考虑的是低通滤波的参数而不是微分算法本身因为算法差异带来的误差远小于噪声放大带来的误差。2.3 三种信号对应的工程应用场景决定了处理策略不同场景对互转结果的精度要求差距很大。像汽车平顺性评价需要的是加速度积分成速度后的均方根值对积分结果的绝对相位不敏感重点是去掉趋势项像建筑结构模态分析关心的是位移的频域幅值这时候频域积分比时域积分更合适像伺服系统位置环调试位移微分得到速度要用于实时控制对延迟极其敏感频域方法根本不能用必须考虑在线式滤波微分。所以完全脱离应用场景去谈互转方法是没法落地的。我会在下一章先把时域积分讲透因为它是理解误差机制的起点。时域积分虽然漂移难缠但胜在理解直观、代码透明、适合实时处理把它的坑摸清楚之后再看频域积分就会顺畅很多。3. 时域积分从加速度到位移漂移处理与最小可运行代码3.1 用 Python 写一个加速度积分的完整流程先给出一段经过验证的时域积分代码它做了三件事去均值、梯形积分、再用高通滤波消除残余趋势。这段代码可以直接用于处理加速度时程数据。import numpy as np from scipy.signal import butter, filtfilt def accel_to_displacement(acc, dt, highpass_cutoff0.5): # acc: 加速度时程数组 (m/s^2) # dt: 采样间隔 (s) # highpass_cutoff: 高通滤波截止频率 (Hz) # 第一步去均值消除直流分量 acc acc - np.mean(acc) # 第二步梯形积分得到速度 vel np.zeros_like(acc) vel[0] 0.0 # 假设初始速度为0 for i in range(1, len(acc)): vel[i] vel[i-1] (acc[i-1] acc[i]) * 0.5 * dt # 第三步速度再积分得到位移 disp np.zeros_like(vel) disp[0] 0.0 for i in range(1, len(vel)): disp[i] disp[i-1] (vel[i-1] vel[i]) * 0.5 * dt # 第四步高通滤波去除积分产生的低频漂移 nyquist 0.5 / dt norm_cutoff highpass_cutoff / nyquist b, a butter(2, norm_cutoff, btypehigh) disp_filtered filtfilt(b, a, disp) return vel, disp_filtered代码逻辑上值得说明的是四步的次序。先去均值这一步很多人会忽略但它直接影响第一次积分的结果加速度里的直流分量哪怕只有 0.01 m/s²积分 100 秒后速度里就会多出 1 m/s 的线性趋势。梯形积分是数值积分里最稳妥的入门选择它用的是相邻两点的平均值比矩形法精度高又不像高阶牛顿-柯特斯公式那样对噪声敏感。最后的高通滤波是整个流程的关键它用零相位滤波的 filtfilt 函数消除了积分造成的抛物线型趋势但代价是引入了边界效应所以两端的数据会有些失真使用时要注意裁剪。3.2 高通滤波截止频率怎么设先看信号的最低关心频率高通滤波的截止频率是最难调的一个参数。设得太低滤不掉积分漂移设得太高会把真实有用的低频位移信号一并滤掉。我一般遵循一个经验原则截止频率设为信号中你需要保留的最低频率分量的三分之一到二分之一。比如被测结构的固有频率是 2 Hz那你至少要把 0.7 Hz 以上的信号保住高通截止频率可以设在 0.3 Hz 到 0.5 Hz 之间。但这个参数需要结合实际采样时长去看。高通滤波本质上是一种对低频成分的抑制滤波器的过渡带大约从截止频率的 0.5 倍处开始如果采样时长只有 10 秒低于 0.1 Hz 的频率分量本来就只有一个半周期滤波结果会在端点处严重扭曲。遇到这种情况要么加长采样时长要么在滤波前对信号两端做延拓。3.3 积分结果里出现抛物线趋势现象、原因、对策如果加速度去均值之后仍然在二次积分里出现明显抛物线最常见的元凶是加速度信号里的低频噪声没有完全去除。去均值只去掉了直流分量但依然会有 0.1 Hz 以下的缓慢变化分量。梯形积分相当于一个理想的数字积分器它对频率的响应幅度与 1/ω 成正比也就是频率越低放大倍数越大。处理这种问题我的习惯是把高通滤波放到第一次积分之前而不是等到二次积分之后再做。虽然上一节代码里我放在了最后那是为了代码演示简单。实际工程中更稳妥的做法是加速度先去均值再经过一次高阶高通滤波再做两次积分。滤波越早在误差放大之前把低频成分压掉后面的积分累积误差就越小。如果两段都做一次在加速度时域、一次在位移后效果最稳健但要注意滤波器叠加带来的额外相位和幅值失真。4. 频域积分用 FFT 做无漂移积分的完整代码与参数设置4.1 频域积分的数学原理除法远好过加法时域积分的最大痛点就是漂移而频域积分通过一个技巧完美绕开了这个问题。对加速度幅值谱做积分相当于在频域里除以 jω即每个频率分量的幅值变为原来的 1/ω相位滞后 90 度。因为直流分量位于 ω0 处除法会出现无穷大所以在频域积分中直接丢掉零频分量这就从根本上消除了直流偏置引起的趋势项。频域积分的第二个优势是能精确控制频带。你可以只对关心的频率范围做积分比如在 1 Hz 到 100 Hz 频带内积分先把噪声频带直接置零再反变换。这种带通积分在模态测试里几乎是标准做法。但它也有明显缺陷要求整段数据一次性处理不能实时输出对长的非周期信号做 FFT 会引入泄漏需要加窗函数频谱分辨率受限于采样时长频率分辨率就是 1/T。4.2 用 Python 实现带通频域积分解决泄漏与边界问题import numpy as np from scipy.fft import rfft, irfft def freq_integral(signal, dt, low_cut, high_cut): # signal: 加速度时程数组 # dt: 采样间隔 # low_cut, high_cut: 积分频带范围 (Hz) n len(signal) # 去掉均值缓解泄漏影响 signal signal - np.mean(signal) # 单边FFT从实数信号获得正频率谱 spectrum rfft(signal) freqs np.fft.rfftfreq(n, ddt) # 生成积分算子幅值乘以 1/(j*2*pi*f)零频置0 integral_op np.zeros_like(spectrum) for i, f in enumerate(freqs): if f 0 and low_cut f high_cut: integral_op[i] 1.0 / (1j * 2 * np.pi * f) # 带通范围内的谱做积分范围外的置零 integrated_spectrum spectrum * integral_op # 逆变换回时域取实部 result irfft(integrated_spectrum, nn) return result这段代码最核心的是循环里构造的积分算子。只有在通带范围内的频率分量才会被保留和乘以 1/jω通带外的分量包括全部噪声频段被直接置零。相比时域积分它不需要事后滤波结果天然没有漂移。代码里使用 rfft 而不是 fft是因为输入信号是实信号用单边谱可以节省一半计算量。逆变换用 irfft 自动处理对称关系不需要手动构造负频率部分。使用这段代码时要注意它的输出是速度信号想继续到位移就再调用一次。两次调用之间不需要任何额外的去趋势处理因为频谱里的零频分量已经被丢弃了。4.3 频域积分的三个隐藏参数频谱分辨率、窗函数和重叠率第一个隐藏参数是频谱分辨率。FFT 的分辨率是 1/T_signal如果采样时长是 10 秒分辨率是 0.1 Hz。这意味着信号里如果有 0.05 Hz 的真实分量它会被平摊到相邻的频率点上积分后幅值会被低估。对于低频积分的场景我建议保证采样时长至少是被关心最低频率周期的 5 到 10 倍。第二个参数是窗函数。直接对截断信号做 FFT如果信号两端不是从零开始会产生频谱泄漏。常见的做法是先乘以汉宁窗再变换但这样会改变信号幅值积分前需要做幅度恢复或者直接用平顶窗。不过在我的经验里如果信号长度内包含了足够多的周期泄漏影响很小不加窗有时比加窗更省事因为避免了窗函数带来的幅值失真。第三个参数是重叠率它主要在有加窗分析时才有意义。如果要对长信号分段做频域积分再拼接段与段之间需要 50% 到 75% 的重叠否则拼接点会出现台阶。对于一次性整段积分的场景这个参数不存在。5. 信号互转中最容易翻车的五个坑现象、原因、对策5.1 加速度积分出的位移是条漂亮抛物线趋势项滤不掉现象数据明明很干净积分两次后位移包络呈抛物线状向上飘调高通截止频率也没有明显改善。原因加速度计零偏没有完全消除或者信号里有低于滤波截止频率的超低频分量。高通滤波器的阶数不够时过渡带太宽低频成分仍有部分残留这些残留在积分中持续累积。对策先检查原始加速度信号做一次功率谱密度分析看峰值是否集中在 0 到 1 Hz 频段。如果是零偏问题用更高阶的高通滤波器阶数从 2 提高到 4如果信号本身存在真实超低频运动就不要强行用高通滤波改用频域带通积分把不关心的低频段直接置零。5.2 对位移求导得到的速度信号全是毛刺真实波形被淹没现象位移信号经过 np.diff 求导后结果布满高频毛刺信噪比极差。原因直接差分是减法两次相邻采样点的噪声是独立的减完之后的噪声标准差变成原始噪声的√2倍。同时高频噪声幅值被放大 ω 倍而位移信号里真实速度的高频成分并没有那么多噪声占比迅速上升。对策不要对原始位移直接差分。先进行低通滤波截止频率选择在你要关心速度的最高频率分量的 1.5 倍位置再中央差分。如果用于实时控制改用一阶惯性滤波加差分的组合先低通后差分差分后用滑动平均降噪。5.3 积分结果两端严重偏离中间看起来完美现象积分出的位移曲线两端出现喇叭口状发散中间一段与参考值高度吻合。原因零相位滤波 filtfilt 是双向滤波两端会预测并补出数据导致边界失真时域积分本身以初始值作为起点初值误差会在整个积分区间内累积。对策对时域积分开头和结尾各丢弃总时长 5% 到 10% 的数据点保留中间稳定段。对频域积分在变换前对信号两端做镜像延拓或多项式延拓变换后把延拓部分裁剪掉。更彻底的方案是改用频域积分因为它的积分结果天然不依赖初始条件。5.4 相同算法在不同采样率下结果差异很大低频段特别明显现象同一段加速度信号用 200 Hz 重采样后积分的结果和 1000 Hz 时完全不同位移幅值差出 20% 以上。原因数字积分器的幅度响应在低频段严重受采样率影响。频率 f 的信号在采样率 fs 下做梯形积分等效幅度增益为 sin(πf/fs)/(πf/fs)当 f 接近 fs 时增益明显下跌。更关键的是不同采样率下信号经过抗混叠滤波器后低频段相位响应也不同。对策统一先对原始模拟信号做抗混叠滤波再以固定采样率采集如果已经是离散数据重采样前必须经过低通滤波。处理完成后用后文第 6 章的合成信号法做一次全流程校验两侧采样率下都应该逼近真实幅值。5.5 加速度、速度、位移的单位和幅值方向反了检查半天排查不出来现象积分结果幅值方向不错但数值大小比理论值小几个数量级或者速度信号是负的。原因加速度传感器输出的单位可能是 g、m/s²、mg 三者混用频域积分过程中没有除以 2π符号约定里加速度方向定义为负时积分后速度符号也会翻转。对策在代码里显式保存一条单位换算常量。先把 g 转换为 9.80665 m/s²再做积分。频域积分中检查算子是 1/jω 还是 1/(j·2π·f)这两个是等价的前者要求角频率 ω后者要求 Hz 频率。最稳妥的办法是在数据预处理阶段强制约定加速度向上为正写入注释后续所有积分和微分结果沿用同一约定。6. 用合成正弦信号验证互转精度一个能救命的闭环测试法无论时域还是频域互转代码写完都别急着上真实数据先用合成信号做闭环验证。这个方法能帮你在一分钟内定位问题是在信号处理还是环境噪声上。构造一组已知的位移正弦信号频率取 5 Hz幅值取 2 mm对它求两次导数得到加速度加入少量高斯白噪后作为算法的输入。正确的结果应该是从加速度积分出的位移在去掉边界段之后能还原出接近 2 mm 幅值的正弦波形。我把这个验证写成一个极简脚本import numpy as np fs 1000 T 20.0 t np.arange(0, T, 1/fs) f0 5.0 disp_true 0.002 * np.sin(2*np.pi*f0*t) vel_true 0.002 * 2*np.pi*f0 * np.cos(2*np.pi*f0*t) acc_true -0.002 * (2*np.pi*f0)**2 * np.sin(2*np.pi*f0*t) acc_meas acc_true np.random.normal(0, 0.01*max(acc_true), len(t))然后把 acc_meas 丢进你写的时域积分或频域积分流程中取中间 80% 的区段计算积分位移与 disp_true 的均方根误差。如果误差在 5% 以内算法正常如果误差超过 20%重点检查滤波截止频率是否压掉了 5 Hz 分量或者频域积分是否把 5 Hz 排在了通带边缘。我习惯把这个验证脚本固定保留下来每次更换传感器或者改变采样率时重新跑一遍。真实数据里往往藏着各种说不清的干扰一旦积分结果异常先用合成信号确认算法本身靠谱再回头排查传感器和采集链路能省下大量无效排障时间。希望这个闭环测试法能帮到你少走几段弯路。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
