基于组合辛普森积分的振动测试时域积分方法

许晓红, 倪春锋, 颜晓青

(沙洲职业工学院, 江苏 张家港 215600)

1 前 言

振动测试是机械量测试的一个重要领域,关系到机器设备的工作性能和使用寿命,强烈的振动噪声还会对人体的生理健康产生危害。振动加速度反映振动的冲击力度,振动速度又叫振动烈度,反映振动能量的大小,振动位移直接反映振动体的位置变化量。理论上对于加速度、速度和位移的测量都有相应的传感器和测试方法,但是振动测量领域往往涉及到大型旋转设备、运动车辆和大型建筑结构或桥梁,这些环境对于位移测量来说,很难找到位移传感器安装所需要的相对静止结构。而加速度传感器不需要相对于待测结构静止的安装位置,直接将加速度传感器与待测结构刚性连接即可[1～5]。所以实际上往往都是通过加速度传感器获得振动加速度,然后通过积分获得振动速度和振动位移,这样通过一个传感器得到3个参量,也具有一定的经济性。

目前常用的积分方法有频域积分和时域积分。频域积分先对原始信号作傅里叶变换转化为频域信号,再对频域信号进行积分运算,最后再进行傅里叶逆变换得到时域信号。由于信号的处理经过了变换域的正变换和逆变换,容易引起截断误差。时域积分对测得的加速度信号直接进行积分,会出现干扰趋势项,需要拟合趋势项以去除信号偏移。在一般的低频振动测试中,考虑到实时分析对资源和效率的要求,采用时域积分更为准确有效[6,7]。本文对时域测量中的数值积分问题进行了研究。

2 测量原理

振动加速度可以通过加速度传感器获得,由于振动加速度初始值往往不为零,再有有些加速度传感器存在着零点偏移,这给通过振动加速度求取振动速度和振动位移的积分过程带来了累积误差:

(1)

式中:a(t)为振动加速度;v(t)为振动速度;s(t)为振动位移;a0,a1,a2为积分过程带来的累计误差。

式(1)所示的累积误差实际上是一种低频信号,表示了速度v′(t)和位移s′(t)的变化趋势,这个趋势项可以通过曲线拟合的方法提取,常用的拟合方法有多项式拟合法和样条函数拟合法[8,9],将积分信号v′(t)和s′(t)分别减去各自的趋势项,就可得到实际的振动速度和振动位移变化曲线。除此之外,加速度测量通道中的高频噪声也对振动速度和振动位移的积分计算产生影响,所以需要采用一个低通滤波器移除高频噪声。整个的数据处理与计算原理方框图如图1所示,由此得到振动速度和振动位移曲线,并由波形曲线进行参数评定。

图1 测量原理框图

3 组合辛普森积分方法

在整个测量过程中,关键在于积分运算。对于被积函数是确定性函数,可以精确得出积分函数,但实际振动加速度都是非确定性函数,需要通过数值积分方法,逼近得出积分函数。常用的数值积分方法有矩形公式法和梯形公式法,它们实际上是用相邻点之间的矩形和梯形逼近积分区域,具有计算简单的特点,但是逼近精度较差。而辛普森积分是基于牛顿-科特斯公式,它是以二次曲线逼近的方式取代矩形或梯形积分公式,具有更高的逼近精度,并且它的计算公式也相对简单,是一种工程计算中的实用方法[10,11]。牛顿-科特斯公式如下:

(2)

(3)

(4)

对于振动速度和振动位移的测量,实际需要根据振动速度和位移的波形曲线,来求取相应的参数,例如说峰峰值、均方根值等,这就需要计算每一点的积分值。而辛普森计算公式只是计算了定积分,并且是3个采样点之间的曲线面积,所以还不能完全直接应用。根据牛顿-莱布尼兹公式,如式(5)所示,一个定积分式的值,就是原函数在上限的值与原函数在下限的值的差。本文在给出振动速度和位移的初始值前提下,通过对辛普森公式的改进,得出递推算法,可求出每一点的振动速度与位移,进而确定振动速度与位移曲线。而对于这两个初始值,工程上可以取零值。

(5)

式中F为f的积分函数。

对于式(4),可以看出辛普森积分计算了相邻2个采样区域的定积分,涉及了3个采样点。由于它不是针对单一采样区域,且只涉及整个积分区域中首尾两点的积分值,所以相邻的辛普森积分序列不可能覆盖所有采样点。为了得到每一点的积分值,提出了一种组合辛普森积分方法,即对采样点中奇数序列采用一套辛普森积分计算公式,对偶数序列采用另一套辛普森积分计算公式,两者组合得到所有点的积分值。

这里设加速度采样值序列为a{(0),a(1),…,a(n-2),a(n-1)},其中n为偶数。由辛普森积分原理,将奇数序列a{(1),a(3),…,a(n-3),a(n-1)}组成一个积分序列,将偶数序列a{(0),a(2),…,a(n-4),a(n-2)}组成一个积分序列。这样可以得到如下的积分计算关系:

(6)

对于振动速度初始值v(0)可以取为0,而对于v(1)可以采用梯形积分公式计算法求得,见式(7)。2部分积分计算组合,就可得到全部采样点的振动速度值,从而可以得到波形曲线,进行参数计算。

(7)

同理,对于振动位移的积分运算,也可在振动速度去趋势项后,使用该方法处理。

4 算法校准验证

为了对上述算法进行验证,可以采用一个确定函数信号,由于确定函数通过理论计算可以求得其积分函数,将这个积分函数作为标准信号与辛普森算法得到的积分函数相比较进行验证。由于周期信号都可以展开傅里叶级数,这里选取确定函数为一正弦函数10 sin(200 π t),其积分函数如式(8),波形如图2所示。

(8)

图2 标准正弦函数和其积分函数

将标准正弦函数10sin(200 π t)增加随时间变化趋势项,合成得到10sin(200 π t)+50t+2,采用第3节所述组合辛普森积分算法对其积分,并对积分结果进行趋势项拟合,获得积分变化趋势项,如图3所示。得到的趋势项p(t)如式(9)所示。将积分结果进行去趋势项处理,如图4所示,可以看出组合辛普森积分去趋势项后的曲线与理论计算的积分曲线高度一致,说明了本方法的处理能力。

p(t)=24.98t2+2.003t+0.0158

(9)

图3 组合辛普森积分及趋势项拟合曲线

图4 去趋势项积分曲线

辛普森积分具有3次代数精度,对于小于等于3次的积分多项式可以精确积分,如果f∈C4(a,b),则存在截断误差项e(f):

(10)

对于组合辛普森积分,需要将截断误差进行累积,设M=max|f(4)(x)|,则总的截断误差E(f)绝对值为:

(11)

式中f(4)(ξ)为函数f的4阶导数。

从式(11)可以看出,总的截断误差E(f)为h4的高阶无穷小,即o(h4),当采样步长h足够小的情况下,总的截断误差也相当小。例如对于上面的标准正弦信号,采样步长为0.000 1,所得总的截断误差小于0.173 2×10-6,可以忽略。

对于振动测试来说,为了保证信号采集的完整性,采样频率通常很高,也即采样周期很短,所以从计算精度和计算效率的角度来考虑,辛普森积分可以满足振动测试的时域积分要求。

5 Matlab振动信号测试

Matlab中包含了一个可视化地震波信号quake,它是一个开放的数据源,记录了1989年美国加利福尼亚州洛马·普雷塔大地震中200 Hz采样频率的地震波信号,它通过加速度传感器记录了东-西、南-北和垂直3个方向的地震波信号[12],实验通过这些数据可以得到3个方向的振动速度和振动位移,图5分别给出了采用本文所述的组合辛普森积分方法得到的3个方向8～15 s间的振动速度和振动位移波形,采样间隔为5 ms,每个方向共 1 401 个采样点。

图5 地震波振动测试

实际上Matlab提供的演示程序quake,也对地震波信号进行了积分处理,得到了3个方向的振动速度和振动位移信号,然而Matlab中采用的积分算法为cumsum,它实际上是矩形积分算法,就是以采样间隔内的矩形面积近似代替这一区域的积分值。图6(a)和图7(a)是Matlab中矩形积分算法得到的振动速度和振动位移波形曲线,图6(b)和图7(b)是本文论述的组合辛普森积分算法得到的振动速度和振动位移波形曲线。从图6(a)和图6(b)、图7(a)和图7(b)的图形曲线对比中可以看出两种曲线形状基本一致,这也说明了本算法的有效性。

图6 2种算法对于振动速度的积分对比

图7 2种算法对于振动位移的积分对比

对于辛普森积分,依据式(11),其截断误差为:

24.31×10-12f(4)(c)

(12)

而矩形积分,其截断误差为:

(13)

从式(12)和式(13)可以看出,组合辛普森积分总的截断误差为o(h4),而矩形积分总的截断误差为o(h2)。所以就算法本身而言,辛普森积分具有更高的计算精度,更小的截止误差。

5 结 论

本文根据振动测试要求,由辛普森积分公式和牛顿-莱布尼兹公式,将原始信号分成奇数序列和偶数序列,推导得出了组合辛普森积分算法。结合对积分初始值的有效处理,得出了相应的递推计算公式。应用此算法可以得到测试数据序列的完整积分波形,对地震波数据的实验,验证了本方法的有效性。并且该算法具有更高的计算精度,非常适于振动测试的时域积分应用。