matlab绘制加速度曲线,加速度积分得到速度和位移的问题小结

在振动信号采集分析中，我们常常使用加速度传感器，是因为加速度传感器具有可靠性高，可测频带宽，结构小巧，抗干扰能力强等优点。但是有些时候，必需要得到速度和位移信号，这个时候就出现了如何通过对加速度积分得到速度和位移的基础问题，这个问题难倒了不少相关人员。在国内的各种教科书上也很少提及，也不知道是为什么。

不管是线性系统还是非线性系统，所产生的信号，在时域满足如下微分方程

其解可以写成

其中v0是积分常数，积分下限t=0表示采样起始时刻，是人为指定的时间轴原点。如果振动没有交流分量，则加速度是常数，可以得到v(t)=v0+at 。对速度积分得到位移的处理类似，所以有

其中d0是积分常数，由于v0和d0通常是未知的，加速度积分得到速度和位移，就存在积分常数问题。

积分常数问题是非常麻烦的问题。首先可以看到，对于任意δv，得到

都对应同样的加速度a(t)，也就是说v0和d0取任意值都对应于同一个测量的加速度和测量的速度，所以要想仅仅凭目前的信息通过积分获取速度和位移，是不可能的，不管是数值积分还是模拟电路积分都不行。本文主要讨论数值积分的问题，关于模拟电路积分存在的问题，推荐大家阅读《压电加速度计和振动前置放大电器》一书，结论就是在高频，放大器增益会下降，不是理想放大器；在甚低频，会受到时间常数或者回路电容限制，也不够精确。这个回复大家可能不满意，坊间不是有各种FFT，各种去趋势项，各种滤波的操作嘛？有人发表在《力学学报》上的论文不也是用模拟积分电路积分了么？我的建议是，谨慎看待那些论文中的说法，你应该有自己的思考和判断。

实际工作中，在无穷多个v0的选择中，可以选择一个最靠谱的数值凑合用一下，因为在这些选项中，可确定最靠谱的那个值一定在某个有限幅值的区间[va，vb]以内。所以，需要增加对被测试物体的状态进行判断，基于这种判断，做出假设，基于这种假设列出约束条件，根据约束条件找到一个方程，然后求出那个看起来好像最靠谱的积分常数。如果假设成立，就得到了一个很Nice的估计；如果假设不成立，误差很大，那就是自娱自乐，你高兴就好。

从可观测性角度看这个问题

下面的讨论限于线性系统，是最近关注某个问题的时候顺便做的笔记。

设振动方程为

则，对应的状态空间模型为

可以得到离散以后的状态方程

可观性的定义：假定A，B，C和D已知，从已经观测的加速度和载荷向量，能否计算得到位移和速度？从式的迭代解可以知道，只需要设法测试得到x(1)即可知道所有后续的状态x(k),k=2,3,4,…。

现在，考虑如下序列

写为矩阵形式

注意到对SDOF模型，上面的矩阵方程始终是一个p×(p+2)矩阵，属于欠估计问题，不可能通过测试加速度同时将载荷和初始状态x(1)求出来。现将该矩阵分割为

其中Op就是可观测矩阵，可以求得

也就是说要联合测试加速度响应，知道可观测矩阵Op、系统矩阵E和载荷才能识别得到所有状态变量，而且可观测矩阵不能奇异，必需是满秩的。2015年Computers &Structures上的一篇关于载荷识别的论文，就是依据对式(8)的简化，假定初始状态x(1)=0，得到

这种零初始条件的假定，在实际问题的应用中，是一厢情愿的。在数据采集器采样了一个数据块处理以后，不能再假定这个数据块采集终了时候的状态变量还是零，所以，至少在后续的计算中，需要考虑 “过去的激励引起的未来的响应的问题”或者x(1)不为0的问题。

同样的道理，实际情况下，我们并没有测试激励力ḟ，还没有做实验获取系统参数，更不能保证传感器布置点所对应的观测矩阵是Regular的，所以单独想要用加速度响应，就完全构造出速度和位移是不可能实现的。就算你得到了，误差可能在80%左右，正确的是量纲。

现在的问题是，我们得到了加速度数据，到底应该考虑哪些实际问题，有没有办法得到一个不是很坏的结果？答案是肯定的，只是这些要求太高了或者假设成立的可能性微乎其乎，所以得到的数据的误差通常较大。

可能的假设

J.S Bendat和A.G. Piersol对信号划分为确定性的和随机的，一个信号是否是确定性的还是随机的，可能存在争议，短时间可能是确定性的，而长时间就变得随机了，或者长时间以后就出现其他“开关条件”触发了不可预测的因素。对确定的信号和随机信号，又分为了如下几种：

在确定性信号中，周期信号是稀有的；在随机信号分类中，各态历经的平稳信号是稀有的。最常见的信号是非周期信号和非平稳信号，是实际中遇到的高概率事件。

所以，下面讨论一些“低概率”出现的假设和“高概率误操作”。

3.1 假设1：位移是零均值的信号

被测试结构没有发生明显的位移现象，所以假定是零均值的。然而，我们知道，有限时间采样以后，得到的是离散数据，这些离散数据的均值并不能严格等于采样区间连续信号的真实均值。而且采样区间真实的均值，极有可能是非零均值的，只是均值很小而已。对速度信号，进行数值积分获取零均值位移信号的方法，同下面的对加速度数值积分获取速度的方法。

3.2 假设2：速度是零均值的信号

设当加速度不为常数时，采样间隔时间dt足够小，一种离散的数值积分为梯形公式

v0任意假定。零均值的获取就是去掉趋势项，得到

上面的式(10)是一种程序语言，代表将有偏估计的速度信号(组成的矩阵)，每个元素去除整体上的已知均值，这样修正得到的v就是零均值的。去掉趋势项的操作，可能存在各种优雅的处理，比如Matlab里面的detrend命令等等，这里不再涉及。

值得一提的是，如果去趋势项的运算是含有非线性因素，那么可能就存在频率混叠的困难，因为采样频率是一定的，这些非线性运算，比如带入平方，开平方根，矩阵求逆，或者取符号运算等等非线性的操作，会引入带宽以上的频率分量，而这可能不是希望的结果，需要后续修正，最好的处理是先保证足够的带宽，最后做无相移的低通滤波。K. Worden教授对去趋势项的研究结论认为，那些“优雅的”去趋势项操作引入的误差很难处理。有趣的是，他早年在MSSP上发表的论文认为，数值积分方案是很Nice的选择；几年以后在他自己的书里对这个问题又做了更细致的论证和分析，他认为数值微分方案才是更好的。

3.3 假设3：速度是周期信号

这种情况，对加速度FFT，在频域积分，再IFFT回来就可以了。

但是我们要知道两点潜在的问题：加速度信号，从时域图上看，一定要是周期的或者Almost周期的，所以分辨率一定要选择得合理才能做到整周期采样；如果频率分辨率选择不合理，可能存在很严重的泄露，加窗号称可以减小泄露，但是，加窗通常对积分以后的速度信号是毁灭性的，加窗的效应保留在时域速度信号里面，无法修正回来；如果有本事做到不加窗，对数据的分辨率选择就提出了很高的要求。如果不是周期信号，怎么选择分辨率都没用啊。

图1 非周期信号+加窗的FFT积分再IFFT结果，以及和精确解(蓝色)对比

3.4 假设4：速度为零时，加速度最大

这个假设，是本人偶然想到的，是在研究某算法时，看到某非线性微分方程强迫响应，然后突然想到，假定速度为零时，对应于动能最小，然后惯性力最大，产生的加速度也是最大的，就这样得到一个结果，在那个例子中，需要在众多的加速度峰值中选择一个，然后确定积分常数，确实产生了较好的效果，但是这个实在太巧合了，再次按照这个方法选择其他峰值做一遍，失效了，现实中是没有这么理想的假设的。

图2 速度(红色)和加速度(蓝色)时域曲线对比

当时的做法是，在某一个明显的加速度峰值附近，选取几个离散加速度点；用多项式进行拟合，得到加速度关于时间的表达式，这一步可以表达为

其中下标k=1,2,3...s只是代表选择的邻近的有限个数据点。

假定系统的阻尼很小，在最大惯性力和最大弹性恢复力平衡时粘性阻尼力为零，这个假设在某些时候是成立的，那么加速度最大的时候，速度可以为零。令上面的表达式为零，选择离散数据中加速度最大的时刻tmax，解方程找到加速度真正的峰值时刻tp，tp非常靠近tmax，即

对加速度的表达式积分得到速度的表达式，积分常数v0待求

tp时刻的速度为零，从而解出积分常数v0

在这个区间的速度表达式就已知了，回到tmax时刻，可以直接得到这个时刻下的速度，应该是一个非常小的量

这样，我们终于找到了“估计的速度小量”，可以根据上面的离散公式，继续完成剩下的时间积分。

因为这个方法并不总是成立，失效以后的某个结果如下所示。

图3 梯形数值积分结果和参考结果(红色)对比，存在一个整体偏移量δv

一个数值积分案例

在网上荡到了一个96年Steven D. Glaser等人在(Preliminary Processing of the Lotung LSST Data报告中，对一个测试到的加速度数据进行数值积分的案例，附录中有源代码，这个数据是地震加速度信号，频率特别低，把源代码贴在这里，给感兴趣的人：

function[a,v,d]=vd(f,Lcut,dT,pre,name)

% function[a,v,d]=vd(f,Lcut,dT,pre,name)

% detrends andfilters the input acc and integrates twice to give velo and disp

% f is the inputacc vector

% Lcut is thelow-end cutoff frequence in Hz

% dT is the timestep

% n is the orderof the Butterworth filter

% pre is thepre-event segment length to be zeroed

zip=zeros(size(zip));

[b,c]=butter(n,Lcut*dT*2,'high');

a=filtfilt(b,c,f);

a=dtrend(a(:,1),1,pre);%%mxl:这个函数已经过时了.

a=detrend(a);

a=a-a(pre+1);

a(1:pre)=zip;

v=dtrend(v(1,:)',1,pre);%%mxl:这个函数已经过时了

v=v-v(pre+1);%整体平移量

v(1:pre)=zip;%前pre个数值设置为0

d=inttrap(v,dT);%%mxl:这个是数值积分，这个函数过时了

[b,c]=butter(n,Lcut*dT*2,'high');

d=filtfilt(b,c,d);

d=dtrend(d(:,1),1,pre);

d=detrend(d); %去趋势项

d=d-d(pre+1); %整体平移量

d(1:pre)=zip; %前pre个数值设置为0

可以看到，他们首先使用了低通滤波，做了去趋势项，然后整体平移了一个量，再将前pre个数值人为置为零。

终极建议

如果同时测试了加速度和位移，那么就特别好办了。所以，如果有条件，尽可能测试加速度和位移。

设在某个区间 [t1,t2] 里面，加速度和位移都可以拟合一个表达式(expression)

在两次积分的情况下可以得到时间常数的表达式，从而解出未知系数

可见速度也同时得到了。

测试位移+测试加速度，可以通过数值微分+低通滤波得到速度；当然也可以直接通过Kalman滤波得到速度，方便、快捷、高效又稳健；或者带有微分算法的低通滤波得到速度，等等。方法总是有很多种。但并不建议测试位移，两次数值微分得到加速度，因为数值微分算法是高通滤波器，会引入高频噪声，采样频率越高，高频噪声问题越严重，两次数值微分得到的加速度误差较大，但也并不是完全不能使用，加上低通滤波可以降低这个误差，总比从加速度积分得到位移要靠谱得多。使用低通滤波算法要注意不要引入相移，平滑处理(Smooth)也是一种低通滤波，效果还行。下面是把不同点数差分算法看成是高通滤波器，然后计算其频率响应得到的曲线，注意是对采样频率fs做了归一化以后的结果。结论是：用5点差分算法是可行的，向前差分或者Matlab的diff函数，最糟糕。

图4 不同差分格式的频响特性

5.1 部分数值实验结果

下面贴出一个位移传感器在低频范围内振动，并做了相对位移和相对加速度采集的结果，从而依据这些测试结果进行数据后处理。

图5 diff会引入高频噪声；加低通滤波能消除这个噪声；smooth是移动平均，也能消除这个噪声

图6 最小二乘法可能计算不准初速度(m/s)，微小的差别最终导致位移无法复现，如图7所示

图7 两种方法计算得到速度再梯形积分得到的位移(mm)对比，红色为测试位移结果黑色和红色几乎完美重合，蓝色为上面least+低通滤波的结果

图8 filter函数会引入相移，而filtfilt函数可以消除相移

图9 最小二乘法初始速度计算不为零，造成了位移恢复结果累积误差较大，产生overshoot现象，精度不及diff+滤波结果

图10 恢复位移和测试结果比较，黑色为diff恢复结果，蓝色为最小二乘加低通滤波结果，最小二乘法可能产生underestimate现象

图11用速度恢复位移和测试结果比较，黑色为diff+fir低通滤波结果，蓝色为diff+移动平均结果，说明diff+移动平均结果也不错

图12 最小二乘法与diff法得到的加速度对比，低通截止0.05

图13 最小二乘法与diff法得到的加速度对比，低通截止0.5

提高低通截止，最小二乘法得到的加速度结果明显优于diff，说明这个时候差分的噪声很难消除。

上面的实验部分实验结果，主要集中在低频，大多数位移测试的也是在低频，少部分如激光测振才能达到较高频率。

有人可能说会有更好的方案，所以，可以求解一些线性或者非线性ODE，得到精确的位移、速度和加速度数据，可以试试对这些加速度进行积分，你可以试试数值积分或者模拟电路积分，并和标准速度和标准位移对比一下，考验一下你的算法的精度。

作者简介

牟小龙，北京理工大学机械与车辆学院博士研究生，擅长线性有限元分析和系统阻抗匹配。

觉得不错，请点赞！

扩展阅读