基于局域均值分解的冲击响应谱时域波形合成新方法

师培峰，邱 伟，王振田，吴瑞斌，陈柯勋，王晓斌



基于局域均值分解的冲击响应谱时域波形合成新方法

师培峰，邱 伟，王振田，吴瑞斌，陈柯勋，王晓斌

（北京强度环境研究所，北京 100076）

冲击响应谱时域波形合成技术是目前实验室环境下模拟实际复杂冲击环境的重要手段。文章在研究局域均值分解（LMD）方法和冲击响应谱时域波形合成原理的基础上，对实测冲击时域波形进行分解得到一组具有不同频率分布的PF信号分量，然后对不同频率分布的PF信号分量进行聚类分析并重新构造得到一组新的PF信号分量，使得PF信号分量的频率−幅值特性与设定的冲击响应谱试验条件一致。将重构后的PF信号分量组合即可得到合成后的冲击时域波形。数值仿真分析结果表明，基于LMD的冲击响应谱时域波形合成新方法生成的冲击时域波形能更好模拟真实冲击环境，冲击响应谱具有较好的控制精度。

冲击响应谱；局域均值分解；波形合成；聚类分析；数值仿真

0 引言

冲击试验是用于考核产品的结构和功能在安装、运输和使用过程中对非重复性冲击环境适应性的重要手段[1]。通常是给定具有一定容差范围的冲击响应谱作为试验条件。目前，在实验室条件下，广泛采用振动台对产品进行冲击模拟试验。采取经典波形（半正弦、后峰锯齿波等）试验控制技术虽然实现方式简单，但往往无法模拟实际复杂的冲击环境，因此目前多采用波形综合法。其合成的冲击信号以等效损伤为依据，基于冲击响应谱对应的时域波形合成的方法和原理，通过对衰减正弦波、合成小波等波形进行合成[2-3]来满足冲击响应谱的瞬态时域波形，从而使合成的波形对应的冲击响应谱在一定容差范围内满足冲击试验制定的试验条件。

由于瞬态冲击信号是一种非平稳信号，目前合成时域波形的冲击响应谱与实际冲击环境往往有差异，且合成过程中需要对波形的幅值参数、时间参数和波形数等进行修正，而修正过程对经验值的依赖较大，并且修正迭代次数较多。局域均值分解（local mean decomposition, LMD）可以实现对非平稳信号（包括实测冲击信号）进行自适应分解[4]，得到一组具有特定时频分布特性的PF信号分量。在分析PF信号分量频率分布特征基础上，利用数据聚类方法对PF信号分量进行聚类分析，从而将具有相同频率或者频率分布误差最小的PF信号分量构造成为新的PF信号分量，并对新的PF信号分量进行幅值修正使得PF信号分量的频率−幅值特性与设定的冲击响应谱试验条件一致，最后将修正后新的PF信号分量加以组合，就可以有效快速完成冲击响应谱时域波形的合成与优化，为高效、准确进行冲击试验提供有力的工具。

1 LMD方法及其特性

LMD是一种新的信号自适应时频分析方法。

其基本原理是利用平滑的局域均值函数和局域包络函数将信号自适应分解成一组乘积函数（production function, PF）之和。该乘积函数由一组调频信号和包络信号组成，通过乘积函数可以准确获取信号的时间−频率曲线和时间−幅值曲线。LMD的具体算法流程见文献[5]所述。将一个给定的非平稳振动信号标记为()，对()进行LMD，按照LMD算法将()分解为个PF信号分量和残余信号()，即

根据文献[5]中的LMD原理可知，残余信号R()为一单调函数，其周期可以认为无限长，故其频率接近0Hz；而且该信号的幅值与其他PF信号分量相比，可以忽略不计。因此本文后续工作基于PF信号开展，而不考虑残余信号R()。

信号()分解后的第个PF信号分量通过包络函数a()与纯调频函数s()相乘得到。即

从LMD的分解原理可以看出，对于一个给定的非平稳信号，通过LMD可自适应得到一组频率由高到低分布[6]的PF信号，进而得到非平稳信号的时频分布特性。

对图1所示的实测冲击信号进行LMD，得到13个PF信号。选取其中的第5和第6个PF信号分量（标记为PF-5和PF-6）进行分析，其时间历程曲线和时频分布特性曲线分别如图2和图3所示。

图1 实测冲击信号时间历程曲线

图2 PF信号分量时间历程

图3 PF信号分量时频分布特性

对图2进行分析可以看出，通过LMD不同频率的信号被提取出来，且LMD具有自适应滤波的功能。其中信号分量PF-5和PF-6的中心频率分别为861Hz和450Hz。但从图3可以看出，这两个信号分量的频率分布分别围绕中心频率波动，具有一定的离散性。

2 基于LMD的冲击时域波形合成

对特定环境下的实测冲击信号进行LMD，对分解后的PF信号进行处理后，将一系列的PF信号合成并优化得到满足冲击响应谱条件的冲击响应信号，合成后的冲击时域波形最后通过振动台实现冲击时域波形复现，从而完成对冲击环境的模拟。与传统的冲击时域波形合成方法相比，基于LMD生成的PF信号合成的冲击信号在满足冲击响应谱试验条件的基础上，能更好地模拟产品在实际冲击环境下的冲击响应。

2.1 传统的冲击时域波形合成方法

传统的冲击响应谱时域波形生成[7]通常采用波形综合法，其总体思路是以频域为基础进行时域波形的合成，然后再对波形进行修正。波形综合法合成满足冲击响应谱条件的瞬态波形过程为：通过选取一组频率、幅值、持续时间不同的基本脉冲波形，对这些波形延时求和，综合得到一个合成波形；通过修正后，该波形的冲击响应谱在一定精度范围内满足设定的冲击试验条件。振动台环境下通常采用3种波形作为基本波形进行合成，即衰减正弦波、合成小波和衰减余弦波。其中合成小波和衰减余弦波可以满足上述条件，而衰减正弦波则需要加修正项。以合成小波为基本波形进行冲击时域信号合成的方法如下：

1）选择基波函数。基波函数表达式为

其中：A为波形的幅值；f为波形频率；N为半正弦波的个数；d为基波的波形延迟时间；d为基波的持续时间。这些参数根据冲击试验规定的冲击响应谱试验条件确定[8]。

2）通过将个基波相加求和，得到合成后的冲击时域波形函数，即

3）将各参数代入式(4)和式(5)得到合成后的冲击时域波形。对冲击时域波形进行冲击响应谱计算，并将冲击响应谱与试验条件进行比较；如果不满足试验条件，则需对上述一系列参数进行修正。

4）将最终满足试验条件的合成后的冲击时域波形通过振动台实现时域波形复现，从而完成冲击试验[9]。

2.2 基于LMD的冲击时域波形合成新方法

与传统的冲击信号时域波形合成方法[10]不同，基于LMD的冲击时域波形合成新方法以实测冲击信号为基础进行时域波形的合成，其算法流程为：

1）对信号()进行低通滤波，滤波的截止频率为H，设冲击响应谱试验条件的频率上下限分别为H和L。

2）对低通滤波后的信号进行LMD，得到个PF信号分量和残余信号()。将这个PF信号分量标记为PFs(=1,…,)。

3）根据式(2)计算方法得到第个PF信号分量PFs的时频分布特性曲线。

PFs(=1,…,)的时频分布特性曲线标记为FT_PFs(=1,…,)。设定第个时间−频率曲线FT_PFs的时间序列为Ts,t(=1,…,；=0, …,)对应的频率序列为Fs,t(=1,…,;=0,…,)，即采样点(Ts,t,Fs,t)构成时频分布特性曲线FT_PFs(=1,…,)所有的点。

4）按照1/3倍频程分布规律计算得到频率区间[L,H]的中心频率值及1/3倍频程带宽。则冲击响应谱频率区间[L,H]的频率个数为，=3×Num+1，其中Num为lg(H/L)/lg2的计算结果取整。

将频率区间[L,H]的1/3倍频程中心频率标记为oct(=1,…,)，相邻的两个中心频率比值关系为oct+1/oct=21/3，每个中心频率oct的1/3倍频程带宽为[oct/21/6,oct×21/6]。冲击响应谱试验条件中，每个中心频率oct(=1,…,)对应设定的冲击响应谱幅值标记为A(=1,…,)。

5）将时频曲线Fs,t(=1,…,;=0,…,)的×个(Ts,t,Fs,t)作为训练样本，以个频率值oct(=1,…,)作为初始的聚类中心，采用K-means聚类算法以频率为计算属性对上述所有训练样本进行聚类分析[11-12]，通过迭代运算将所有的(Ts,t,Fs,t)分别聚集到个不同的频率类别中。

6）计算所有训练样本值(Ts,t,Fs,t)与第个聚类中心oct的距离(i,j,t)，(i,j,t)=|Fs,t−oct|。设定判据如下：

如果满足式(6)的条件，则说明训练样本值(Ts,t,Fs,t)被聚集到oct为中心的类别中，且Fs,t在oct的带宽范围内。

7）将第个（聚类中心为oct）类别集合中Ts,t按时间顺序排列，取Ts,t对应的()加速度值为新的幅值，可以得到一个新的时间历程曲线，标记为NPFs，这样一共能得到条时间历程曲线NPFs(=1,…,)。新的时间历程曲线NPFs是一个调频曲线，其中心频率为oct，1/3倍频程带宽为[oct/21/6,oct×21/6]。

8）对曲线NPFs采用改进递归数字滤波器法[13]计算冲击响应谱幅值a，与设定的参考幅值A进行比较，如果|a−A|＞，则直接乘以修正系数A/a对NPFs加速度幅值进行修正。

9）将条时间历程曲线NPFs(=1,…,)叠加得到新的时间历程信号()，为满足冲击试验开始和结束时刻台面加速度、速度、位移为0的必要条件，生成一组修正信号()，()的信号频率在冲击试验频带范围之外，同时满足式(7)的条件，即

建立“一带一路”争端解决机制和机构是中国在“一带一路”倡议落地时期的重要举措。长期以来，国际商事争端解决较多依赖于调解和仲裁机制，随着“一带一路”倡议的全面落地，如何实现诉讼、调解、仲裁的有效衔接，这已成为“一带一路”纠纷解决机制的重要课题。2018年6月29日，中国最高人民法院分别在深圳和西安设立了第一、第二国际商事法庭，7月1日颁布施行了《关于设立国际商事法庭若干问题的规定》，为形成“一站式”国际商事纠纷解决机制迈出了坚实的一步。国际商事法庭是中国全新的尝试，是否有必要设置国际商事法庭？国际商事法庭如何与国内商事审判制度协调？如何设计国际商事法庭？本文将就上述问题展开探讨。

最终合成新的冲击时域波形记为syn()，则有：

通过上述9个步骤，可以合成新的冲击响应时域信号。对冲击时域波形进行信号上采样，然后通过振动台将冲击信号复现，即可按照冲击试验条件快速完成冲击试验。上述算法流程如图4所示。

图4 冲击响应谱时域波形信号合成新方法流程

3 仿真分析

对某产品进行冲击试验，试验条件为：10Hz，0.8；160Hz，20；2000Hz，20。试验持续时间≤1s，冲击响应谱容差为±3dB。该产品冲击环境下实测信号（参见图1所示）的冲击响应谱曲线与试验条件规定的参考谱如图5所示。

图5 信号冲击响应谱及试验参考谱曲线

从图5可以看出，对冲击环境下的实测振动信号进行冲击响应谱分析，其结果无法满足冲击试验条件。具体表现为：在低频段（10～100Hz）的冲击响应谱量级小于设定的幅值，而在高频段（160～2000Hz）的冲击响应谱大于设定的幅值。

根据2.2节中的算法描述，对图1所示的实测冲击信号进行LMD，得到13组PF信号分量和1组残余信号，其中残余信号为单调函数且幅值最大值为0.005。冲击响应谱频率区间为10～2000Hz，根据算法流程中式(5)的计算方法，最终聚类中心个数为23。也就是说，通过对初始分解得到的13个PF信号进行聚类分析，得到23个新的PF信号分量。由于这23个新分量为纯调频函数，其频率值与冲击试验条件中的频率值一一对应，通过与冲击试验条件频率-幅值关系对比，可以方便实现对PF信号分量进行幅值修正。根据冲击试验前后冲击试验台加速度、速度和位移为0的条件，生成修正信号。修正信号的频率成分为高于2000Hz的高频、小量级信号，不会对分析频率范围内有效冲击信号产生影响。

根据冲击试验条件，分别利用传统冲击时域波形合成方法（波形综合法）与基于LMD的冲击时域波形合成新方法合成的时间历程曲线如图6和图7所示，其对应的冲击响应谱曲线如图8所示。

图6 波形综合法合成的时域波形

图7 本文新方法合成的时域波形

图8 冲击响应谱对比图

图8的结果表明，传统方法和基于LMD的新方法合成的冲击时域波形均能够满足设定的冲击试验条件的要求。新方法合成的冲击响应谱波动较小，具有较好平坦度；而传统方法合成的冲击响应谱波动和容差较大。此外，新方法合成的冲击时域波形持续时间与实测冲击信号接近，能更好地模拟实际冲击环境。

4 结束语

本文提出了一种基于LMD的冲击时域波形合成新方法。与传统的冲击时域波形合成方法中以合成小波和衰减余弦波作为基本合成波形不同，本文方法以实测冲击信号为基础，以对其进行LMD后得到的PF信号分量为基本合成波形，通过对PF信号分量进行聚类分析和重新优化组合得到。从2种方法的实际效果来看，它们合成的冲击信号都能够满足设定的试验条件，但应用本文方法生成的冲击响应谱在工作频带内波动更小、曲线更平坦。本文方法可广泛应用于试验室环境的冲击试验，具有一定的推广价值。

[1] 振动与冲击手册编写组. 振动与冲击测试技术[M]. 北京: 国防工业出版社, 1990: 5-8

[2] ESI Group. AutoSEA/SHOCK user’s guide[Z]. New York: Vibro-acoustic Inc,1999: 69-84

[3] 刘洪英, 马爱军. 冲击响应谱时域合成算法研究[J]. 航天医学与医学工程, 2002, 15(6): 437-441

LIU H Y, MA A J. A study for time-history waveform synthesis of algorithm in shock response spectrum (SRS)[J]. Space Medicine & Medical Engineering, 2002, 15(6): 437-441

[4] 张亢, 程军圣, 杨宇. 基于有理样条函数的局部均值分解方法及应用[J]. 振动工程学报, 2011, 24(1): 443-454

ZHANG K, CHENG J S, YANG Y. The local mean decomposition method based on rational spline and its application[J]. Journal of Vibration Engineering, 2011, 24(1): 443-454

[5] SMITH J S. The local mean decomposition and its application to EEG perception data[J]. Journal of the Royal Society Interface, 2005, 2(5): 443-454

[6] 任达千, 杨世锡, 吴昭同, 等. 基于LMD信号瞬时频率求取方法及实验[J]. 浙江大学学报(工学版), 2009, 43(3): 523-528

REN D Q, YANG S X, WU Z T, et al. Instantaneous frequency extraction method and experiment based on LMD[J]. Journal of Zhejiang University (Engineering Science), 2009, 43(3): 523-528

[7] 何志勇, 宋少伟, 邓长华, 等. 瞬态时域数据合成冲击响应谱算法研究[J]. 火箭推进, 2013, 39(5): 56-59

HE Z Y, SONG S W, DENG C H, et al. Research on transient time-domain data synthesis of shock response spectrum[J]. Journal of Rocket Propulsion, 2013, 39(5): 56-59

[8] Spectral dynamics: shock synthesis operating manual[M]. San Jose, CA: Spectral Dynamics Inc, 2005: B1-B3

[9] 金恂叔. 航天器爆炸冲击环境及其模拟试验技术[J]. 航天器环境工程, 1999, 26(3): 1-12

JIN X S. Explosion impact environment of spacecraft and its simulation test technology[J]. Spacecraft Environment Engineering, 1999, 26(3): 1-12

[10] 田芳怡, 田晓涛. 基于演变谱理论和正交化HHT法的冲击响应谱时域信号合成方法[J]. 航天器环境工程, 2014, 31(1): 25-28

TIAN F Y, TIAN X T. A time domain signal synthetic method for shock response spectrum based on evolutionary spectra theory and orthogonal HHT[J]. Spacecraft Environment Engineering, 2014, 31(1): 25-28

[11] HAN J W, KAMBER M. Data mining: concepts and techniques[M]. 2nd ed. San Francisco: Morgan Kaufmann Publishers, 2006: 383-386

[12] 谢娟英, 高红超. 基于统计相关性与K-means的区分基因子集选择算法[J]. 软件学报, 2014, 29(9): 2050-2075

XIE J Y, GAO H C. Statistical correlation and K-means based distinguishable gene subset selection algorithms[J]. Journal of Software, 2014, 29(9): 2050-2075

[13] 王翠荣, 施广富, 郭军. 固体火箭发动机冲击信号响应谱分析[J]. 固体火箭技术, 2003, 26(2): 57-60

WANG C R, SHI G F, GUO J. Analysis of impact signal response spectrum of solid rocket motors[J]. Journal of Solid Rocket Technology, 2003, 26(2): 57-60

（编辑：闫德葵）

A new method for time domain waveform synthesis of shock response spectrum based on local mean decomposition

SHI Peifeng, QIU Wei, WANG Zhentian, WU Ruibin, CHEN Kexun, WANG Xiaobin

(Beijing Institute of Structure and Environment Engineering, Beijing 100076, China)

The technique of time domain waveform synthesis of shock response spectrumis an important approach to simulate the complex shock environment in laboratory.The measured time domain waveform is decomposed by the local mean decomposition (LMD) method, and a set of PF signal components with different frequency distributions are obtained based on the LMD and the time domain waveform synthesis of the shock response spectrum.Then, a new set of PF signal components is reconstructed by using the data clustering analysis method for reconstructing the prior PF signal components.The frequency-amplitude characteristics of the new PF signal components are made to be consistent with the set of the shock response spectrum test conditions.The impulse signal can therefore be reconstructed by properly combining the reconstructed signal components.Simulation results show that the time domain waveform obtained from the shock response spectrum generated by the new method based on the LMD method can well represent the real shock environment, and the shock response spectrum has a good control accuracy.

shock response spectrum (SRS); local mean decomposition (LMD); waveform synthesis; data clustering analysis; numerical simulation

V416.2; TB115.2

A

1673-1379(2018)02-0135-06

10.3969/j.issn.1673-1379.2018.02.006

师培峰（1982—），男，硕士学位，高级工程师，主要研究方向为信号处理、数据挖掘及振动控制系统设计等。E-mail: captain.spf@163.com。

2017-11-29；

2018-04-01