舰船冲击谱若干计算方法比较研究

谢 浩，冯麟涵，吴静波，李晓文，郭 君

(1.哈尔滨工程大学 船舶工程学院，哈尔滨 150001；2.中国人民解放军92537部队，北京 100007)

舰船冲击谱若干计算方法比较研究

谢 浩1，冯麟涵2，吴静波2，李晓文1，郭 君1

(1.哈尔滨工程大学 船舶工程学院，哈尔滨 150001；2.中国人民解放军92537部队，北京 100007)

在水下爆炸冲击载荷作用下，舰船通常采用冲击谱来考核设备抗冲击性能。但由于没有统一计算标准，不同工程人员采用冲击谱的算法不同，得到的计算结果差异很大。针对现在已有冲击谱的不同计算方法，分别采用直接积分法、递归法、三次样条法、递归滤波法等计算方法对滤波后的信号求解冲击谱，通过改变采样频率和采样精度分析不同计算方法的使用条件、精度、效率。结果表明在采样精度一定时，与其他方法相比，递归法、递归滤波法计算效率高、稳定性好，适合实际工程使用。

振动与波；冲击谱；直接积分法；递归法；三次样条法；递归滤波法

冲击谱是潜艇和水面舰船研究中广泛采用的描述冲击环境的方法，其计算方法正确与否关系着舰艇的生命力，一旦计算错误，会造成不可估量的损失。自从1943年冲击谱被M.A.Biot提出后就被广泛应用到航天、地震以及船舶工程等各个领域[1]。在水下爆炸领域，冲击谱作为考核设备抗冲击性能的依据，其计算方法引起人们广泛关注。1981年David Smallwood在文献[2]中提出改进的递归滤波法，该方法克服传统递归滤波法的缺陷，采用斜坡不变模型；1983年陆圣才在文献[3]中提出一种新的计算方法—三次样条法，该方法原理简单，计算效率高；1984年王树棠在文献[4]中，对现有的冲击谱计算方法进行比较，得到改进的递归滤波法有很大的优越性的结论。1996年S.Smith提出了计算冲击谱的若干准则[5]；2009 年 Allan G.Piersol在文献[6]中介绍了不同激励下的冲击谱曲线，并指明其应用领域；2012年Justin N.Martin对Kelly的计算方法进行修正[7]。虽然冲击谱的计算方法很多，但是在水下爆炸领域中，实际工程人员对冲击谱还没有统一的计算标准。文中基于常见的几种冲击谱计算方法，通过对比分析，找到一种可靠的计算方法，为舰船的抗冲击评估提供使用依据。

1 冲击谱的计算方法

本小节简要介绍根据冲击谱导出的几种主要计算方法。根据其计算原理进行分类，如图1所示。

图1 冲击谱计算方法分类图

其计算公式如下

式中δ——相对位移

ξ——系统阻尼系数

ω——系统无阻尼固有频率

直接采用Duhamel积分法，得到系统的响应

式中ωd——系统有阻尼固有频率

(1)直接积分法[8]

其优点是编程简单，但每求一个δi都要调用SIN和EXP等子程序，因而在i较大时，求和项数多，故计算时间长。

(2)递归直接积分法[8]

(3)递归法[8]

上面递归直接积分法是将冲击输入采用一次多项式去逼近，为了提高精度，下面是另一种采用多项式去逼近的递归方程组

的1阶前差分为

递归法的优点是计算循环进行，简化了程序，其计算结果的精度依赖于所取多项式。

(4)递归数字滤波法[8]

该方法是基于脉冲不变假设的模拟滤波器，数字递归滤波器的公式为

P0、P1、Q1、Q2为滤波系数；

该法具有计算速度快、占用内存小的特点，使用广泛。但对采样频率和采样精度要求高。

(5)改进的递归数字滤波法[2]

该法是用广义斜台函数取代脉冲不变模型的ü(t)=δ(t)。广义斜台函数为

式中u(t-KΔt)是单位阶跃函数，A是在t=KΔt时斜台的斜率。则有如下的斜台不变模型的递归公式

对于计算低频响应时采样率过高而采样精度不够、计算结果误差大的情况，可用如下递归公式解决

(6)三次样条函数法[3]

假定在区间[ti，ti+1]内，函数为

代入方程，求出基本递推公式

该法特点是原理简单，占用内存小，计算时间短。

（7）微分方程法[9]

利用Matlab中的Simulink库去求解微分方程，默认初始位移和速度为零，采用建立如图2所示的流程图。

该方法由于采用Matlab自带的模块，几乎不需要手动编写代码，实现方法简单。采用多个Scope模块，便于查看结构相对基础的位移、速度、加速度随时间变换的曲线，并将计算的结果以矩阵的形式输出，便于直接观察。该方法对采样频率有很高的要求。

图2 微分方程法计算流程图

（8）拉普拉斯变换法[9]

拉普拉斯方法仅适用单输入单输出的线性定常系统，是线性系统的时域表达式。

对式（1）两边同时做拉普拉斯变换，得到

整理得到传递函数为

在Matlab中Continuous模块的Transfer Fcn模块输入上面的传递函数，其程序框如图3所示。

图3 拉普拉斯变换法计算流程

该方法编程简单，直接采用已有的模块计算。

(9)状态空间法[9]

状态空间表达式是由状态方程和输出方程组成，它允许非零值的初始条件。其状态空间方程如式（13）所示

将式(13)代入到图4中的Matlab状态空间计算流程中。

图4 状态空间计算流程图

该方法计算精度高，其最大特点为允许设置初始位移和初始速度，对于某些有初始条件的问题，可以得到很好的计算结果。

2 冲击谱各计算方法的精度及效率比较

2.1 各冲击谱求解方法与采样率的关系

在数据处理及应用方面，抽样定理具有广泛的用处，又称采样定理、取样定理[10]。为了给设备提供有效的冲击环境数据，需选取对设备或船体局部结构有效频段内的数据。因此选取低通滤波截止频率fm为300 Hz[11]（也称为分析频率），并保留200 Hz防卫带，则理论上采样率为1 kHz时，即可保证离散信号的完整。文中以1 kHz为基准采样频率，实际采样频率为基准频率的λ倍，即λ为采样系数

式中fs为采样频率；fmax为最高频率；n至少大于2。

图5所示为圆柱壳舱段受水下爆炸的某工况下典型位置的原始加速度响应时历曲线（采样率均为10 kHz，采样精度为10-5，计算自然频率为1 Hz～10 kHz），计算时长为1.2 s，包含首次气泡脉动。

图5 原始加速度响应时历曲线

通过傅里叶变换能清晰反映滤波前后的差异。

从图6中两曲线对比可以发现，滤波之后，两条曲线300 Hz以下完全一致，滤波之后300 Hz～450 Hz之间有较大降低，450 Hz频率以上的信号基本全都被滤掉。根据以上研究，对加速度数据信号进行300 Hz滤波，阻尼系数为0.01，研究各冲击谱计算方法在同一采样精度、不同采样率下的求解效果。

2.1.1 λ=1时各方法求解结果

从图7中看出，所有冲击谱曲线在150 Hz之前变化趋势一致；微分方程法和状态空间法会在低频部分出现毛刺，但冲击谱曲线变化趋势大体一致；但在300 Hz以上时，直接积分法、递归直接积分法、三次样条法会出现发散，前两者表现为在500 Hz以后交替出现波峰和波谷，后者则是在550 Hz处出现无穷大；递归法、递归滤波法、改进递归滤波法三者曲线基本重合；拉普拉斯法与改进递归滤波法在300 Hz以后出现微小差异。

图6 原始加速度信号及300 Hz滤波后傅里叶谱曲线

图7 λ=1时各求解方法对应的冲击谱曲线

说明在λ=1时，递归法、递归数字滤波法、改进的递归数字滤波法、拉普拉斯法适用。

2.1.2 λ=5时各方法求解结果

从图8中看出，所有曲线在300 Hz之前变化趋势大体一致。

图8 λ=5时各求解方法对应的冲击谱曲线

直接积分法和递归直接积分法的冲击谱曲线在高频区域与其他方法出现不同，在1 kHz以上发散并交替出现波峰和波谷；三次样条法在2.7 kHz以内与其他曲线一致，在2.7 Hz以上发散，值无穷大；微分方程法变化趋势与其方法一致，但会在低频部分出现毛刺，在高频部分与其他方法有微小差异；递推法、递推数字滤波法、改进的递归数字滤波法、拉普拉斯算法、状态空间法曲线完全一致。说明在λ=5时，递归法法、递归数字滤波法、改进的递归数字滤波法、拉普拉斯法适用。

2.1.3 λ=10时各方法求解结果

从图9中看出，直接积分法和递归直接积分法在1 kHz以上与其他方法出现不同，在1 kHz以上发散，交替出现波峰和波谷；三次样条法在5 kHz以内与其他曲线一致，5 kHz以上发散，值无穷大；递推法、递推数字滤波法、改进的递归数字滤波法、拉普拉斯算法、微分方程法、状态空间法曲线完全一致，说明在λ=10时，递归法、递归数字滤波法、改进的递归数字滤波法、拉普拉斯法、微分方程法、状态空间法适用。

图9 λ=10时各求解方法对应的冲击谱曲线图

2.1.4 λ=50时各方法求解结果

直接积分法在λ=50时单次计算时间超过1小时，失去了实际意义，因此不予讨论。

从图10中看出，所有曲线变化趋势完全一致，只有递归直接积分法在4 kHz以后出现微小差异；说明在λ=50时，递归法、递归数字滤波法、改进的递归数字滤波法、三次样条函数法、递归直接积分法、状态空间法、微分方程法、拉普拉斯算法都适用。

图10 λ=50时各求解方法对应的冲击谱曲线图

比较图7和图8可以看到，当采样率较低时，低频段曲线几乎没有差别，高频段的曲线则容易失稳，出现局部的毛刺，即使能够实现计算，但也不能作为正确的冲击谱曲线。采样率的选取直接关系到冲击谱曲线的正确与否。同样通过各曲线失稳点也可以看出冲击谱曲线与计算的自然频率有着直接关系，自然频率越大，所要求的采样频率也越高。

以下考察各计算方法所用时间。计算时长的考核使用计算机硬件为8核CPU，型号为Intel(R)Celeron(R)i7-4790 CPU，主频为 3.6 GHz，内存(RAM)为4.00 GB，软件使用Matlab R2012b，在计算机只运行该软件时，分别统计单次计算时间，列入表1（－号表示计算时间超过一小时）。

表1 各求解方法在不同采样率情况下的求解时间对比/s

从表1可以看出直接积分法占用机时最久，在采样率高时，单次计算时间超过了1小时；递归直接积分法、递归法、递归滤波法、改进递归滤波法占机内存小、效率高，随着采样率的提高，计算时间稳定增加，在容许时间内；微分方程法占机时间仅次于直接积分法，单次超过20 s，效率低；状态空间与拉普拉斯法单次时间超过12 s，随采样率的提高，计算时间增加较大，不利于大数据计算；三次样条计算时间居中，但随着采样率的增加，计算时间成倍增加，耗时长。

2.2 各冲击谱求解方法与采样精度的关系

由上一节可知，在采样率低的情况下，直接积分法、递归直接积分法和三次样条法会出现失稳情况，然而冲击谱曲线与采样的精度也有着密不可分的关系，同样以上节的加速度响应为输入信号，滤波截止频率为300 Hz，阻尼系数为0.01，考虑到冲击加速度信号的最大值超过160 g；分别计算输入信号精确到百位、十位、个位、小数点一位、两位下的冲击谱，以此寻求各种冲击谱计算方法的适应精度。下面以递归法为例说明其采样精度对计算结果的影响。

2.2.1 λ=100时不同精度条件下的结果

图11 λ=100时不同采样精度条件下的冲击谱曲线

从图11中可以看到，当采样精度不同时，各种精度条件下的冲击谱曲线在10 Hz以后基本重合；只有精确到百位时，冲击谱曲线会在低频区出现不一样，表现为明显偏大。说明对于递归法，当加速度响应信号精确到十位，就能满足计算要求。同时也可以说明在采样频率足够高的情况下，采样精度决定低频部分的冲击谱曲线。

2.2.2λ=10时不同精度条件下的结果

从图12中可以看到，当采样精度不同时，各冲击谱曲线在20 Hz以后基本重合；只有精确到百位、十位时，冲击谱曲线会在低频区出现不一样，精确到百位时与其余曲线相比明显偏大，而十位时只在2 Hz之前才有微小区别。说明在λ=10时，加速度响应信号精确到个位，就能满足计算要求。同时也可以看出采样精度决定着冲击谱的低频部分曲线，采样精度越高，出现失稳点的频率越低。

图12 不同精度条件下的冲击谱曲线

2.2.3不同采样率条件下根据同一精度得到的结果

从图13中可以看到，当采样精度为十位时，随着采样率的增加，各冲击谱曲线在100 Hz之前几乎完全重合；当采样频率为1 kHz时，与其他采样率相比，冲击谱曲线在100 Hz之后出现差异，曲线不平滑，有间断的峰值。

说明采样率对高频部分的冲击谱有明显影响。采样率越高，冲击谱曲线出现失稳点的频率越高。

图13 不同采样率下的冲击谱曲线

3 结语

（1）从总体来看，滤波之后数据的采样率要达到1 kHz才能确保各计算方法的正确性，而在未滤波之前，由于结构的高频振动会造成数据产生虚假的峰值，滤波之前原始数据的采样率要达到10 kHz才能保证滤波之后数据的正确性，因此在仿真输出及试验采集信号时时间间隔应小于1×10-4s。

（2）各冲击谱求解方法中，直接积分法占用机时最久，当采样率较低时在高频段曲线发散,不适合实际计算；递归直接积分法计算效率高，但在采样率较低时，曲线会出现交替的波峰和波谷，当采样率超过50 kHz时，能得到稳定的结果，适用于采样率较高的情况；三次样条函数法在采样率较低时出现极值并且占用机时较长，在采样率较高时可以使用；微分方程法、状态空间法编程简单，求解容易，对数据的要求性较低，但在采样率较低时，会出现低频段曲线的毛刺，适合于较高采样率下的计算；递归法、递归数字滤波法、改进的递归数字滤波法、拉普拉斯变换法在采样率超过5 kHz时，具有计算精度好、时间短的特点，适合实际的工程计算。其中递归数字滤波法在效率和精度方面都能得到很好的结果。

（3）采样频率的高低决定着高频段曲线失稳点的频率，采样率越高，失稳点的频率越高；采样精度的高低决定着低频段失稳点的频率，采样精度越高，失稳点的频率越低。冲击谱曲线与计算的自然频率有着直接关系，自然频率越高，所要求的采样频率也越高。

[1]卢来洁，马爱军，冯雪梅.冲击响应谱试验规范评述[J].振动与冲击，2002(21)：18-20.

[2]DAVID SMALLWOOD.An improved recursive formula for calculating shock response spectra[C]//Shock and vibration Symposium.San Diego,CA,1980.

[3]陆圣才.计算冲击谱的一种新方法-样条函数法[J].振动与冲击，1986(3)：71-76.

[4]王树棠.计算冲击响应谱方法的比较[J].强度与环境，1985(3)：11-15.

[5]S SMITH,R MELANDER.Why shock measurements performed at different facilities don’t agree[C].Proceedings of the 66 th Shock and Vibration Symposium:Biloxi,Mississippi,1995.

[6]ALLAN G PIERSOL,THOMAS L PAEZ.Harris’shock and vibration handbook[M].NewYork:McGraw-Hill,2009.

[7]JUSTIN N MARTIN.On the shock-response-spectrum recursive algorithm of Kelly and Richman[J].Shock and Vibration,2015,19(1):19-24.

[8]唐照千，黄文虎.振动与冲击手册[M].北京：国防工业出版社，1990.

[9]李颖.Simulink动态系统建模与仿真基础[M].西安：西安电子科技大学出版社，2004.

[10]郑君里，应启珩，杨为理.信号与系统[M].北京：教育出版社，2000.

[11]姬秀滨.小水线面双体船水下爆炸冲击环境特性分析[D].哈尔滨：哈尔滨工程大学，2014：23-25.

Comparative Study on Several Calculation Methods for Ship’s Shock Spectra

XIE Hao1,FENG Ling-han2,WU Jing-bo2,LI Xiao-wen1,GUO Jun1
（1.College of Shipbuilding Engineering,Harbin Engineering University,Harbin 150001,China;2.92537 Unit of PLA,Beijing 100007,China）

The shock spectrum is usually used to assess impact resistant performance of the equipment under underwater explosion.But,there are several computation methods currently and different algorithms for the shock spectra will lead to quite different calculation results.In this paper,based on different calculation methods,such as direct integration,recursive method,cubic spline method and recursive filtering method,the shock spectra of the filtered signals are solved.The conditions,precision and efficiency of the different calculation methods are analyzed by changing the sampling frequency and sampling precision.The results show that the recursive method and recursive filtering method are more efficient and stable than the other methods in the case of the same sampling precision and are suitable for practical engineering application.

vibration and wave;shock spectrum;direct integration;recursive method;cubic spline;recursive filtering method

O175.1

：A

：10.3969/j.issn.1006-1355.2017.04.023

1006-1355(2017)04-0115-06

2017-03-15

工信部高技术船舶资助项目；黑龙江省博士后科研启动基金资助项目（LBH-Q16054）

谢浩（1993－），男，湖北省洪湖市人，硕士生，主要研究方向为舰船结构振动与冲击、舰船结构动力学。

E-mail:xieaho@hrbeu.edu.cn