基于时间序列谐波分析的死亡事故总量波动趋势分析

颜 峻

（中国劳动关系学院 安全工程系，北京 100048）

0 引言

我国安全生产形势呈现出总体稳定、趋向好转的发展态势，但在事故指标不断下降的同时，仍存在着反复波动的现象。这一现象在事故时间序列上表现为长周期稳定变化和短期波动变化2种趋势。若不将2种趋势加以分解，将会影响对安全生产事故趋势的精准判断。因此，本文的主要目的是将原本耦合在一起的2种趋势加以分解，分别在时域和频域上构建描述我国生产安全事故变化趋势的特征模型。

事故时序分析是指基于历史数据统计的事故时间序列变化规律的研究，主要方法包括概率分布函数、自回归移动平均模型、支持向量机、相关向量机、人工神经元网络、灰色模型等。到目前为止，研究人员已经初步发现事故时间序列存在着不同的波动特征。文献[1]采用季节指数法分析了月度事故序列中存在的季节因素，通过建立不同时段事故数的概率分布模型，研究事故发生的时间段所具有的周期特征。文献[2]建立了矿山安全事故季节周期回归模型，得到事故发生的趋势线（即未来事故的平均数），事故数围绕这一趋势线上下波动。在趋势预测方面，由于ARⅠMA模型可以有效的处理非平稳时间序列，因此被广泛用于事故趋势分析和预测。文献[3]针对航空装备事故时序数据具有的非平稳性、自相关性特征，构建了事故时序数列的ARⅠMA（4,1,4）模型，通过一次差分将非平稳数列转换为平稳数列对事故数进行了预测。文献[4]建立了道路交通事故ARⅠMA（0,1,1）预测模型，发现时序数列残差具有较为明显的季节性波动趋势。为了消除季节性因素的影响，文献[5]采用季节调整X-12方法去消除事故月度数据中的季节性影响因素后，建立了针对新序列的灰色GM（1,1）模型和随机波动项的ARⅠMA模型。文献[6]采用相关向量机（RVM）方法建立了交通事故时间序列预测模型，研究了数据中具有的线性关系和非线性关系。文献[7]采用人工神经元网络方法对事故趋势进行了分析，发现事故变化趋势具有长期稳定和短期波动2种特征，尤其是月度数据序列具有明显的周期性特征。此外，灰色预测模型也被广泛用于事故趋势预测，但该方法并没有区分时序数列包含的长期趋势项和短期波动项[8-9]。

综上所述，生产安全事故时间序列具有明显的非平稳性特征，同时受多种因素影响而呈现出2种不同特征项，即长期稳定趋势和短期波动趋势，两者交织在一起增加了变化规律研究和预测的难度。已有研究表明，事故序列具有一定的周期性，但大多难以分解2种特征趋势。虽然X-12等季节调整方法可研究事故序列具有的季节波动特性，但仍将长期平稳趋势和短期波动趋势视为一体不能分开。本文将从生产安全事故月度时间序列变化趋势研究着手，对序列的长期趋势项和短期波动项加以分解研究。将采用Hodrick-Prescott滤波方法提取事故序列中的长期趋势项，采用线性回归模型确定事故变化的长期趋势特征；对分解后的短期波动项，采用谐波分析中的快速傅里叶变换将时域序列转换到频域上，以研究事故变化趋势在短期上具有的周期性变化特征。

1 事故序列平稳性检验

收集我国自2011年1月至2017年11月造成人员死亡（包括下落不明）的较大及以上生产安全事故月度统计数据，绘制时间序列，如图1[10]。月度事故序列存在长期波动下降趋势，为非平稳时间序列。为了进一步验证序列具有的非平稳特征，利用E-views[11]软件中的ADF Fisher单位根检验方法[12]对月度事故序列进行平稳性检验，其检验模型为：

图1 月度较大及以上事故序列Fig.1 Monthly larger and above accident sequence

式中：

t—表示序列的时间变化趋势，月；

α—常数项；

β—时间项系数；

δ—前一期变量系数；模型假设均为H0：δ=0，即存在单位根。从图1可以看出，事故序列具有明显的随时间下降趋势，因此检验模型应包含时间趋势项βt，平稳性检验结果，见表1。

表1 月度事故序列平稳性检验（水平值）Tab.1 Monthly accident sequence stationarity test（level）

表1为包含时间趋势项的ADF单位根检验结果，检验结果为稳定序列。如果把检验模型中的时间趋势项去掉后进行平稳性检验，检验结果不能拒绝原假设，即事故序列为非平稳序列。综合两种不同形式的检验模型所得平稳性检验结论，月度事故序列为趋势平稳序列，即含有显著的时间趋势项。

2 长期稳定趋势分析

上述“事故时间序列为趋势平稳序列”的检验结果证明月度事故序列具有长期稳定趋势项。采用Hodrick-Prescott滤波方法[13]对事故序列中具有的长期趋势项和短期周期性波动项进行分解。设Yt为包含长期稳定性趋势项和短期周期性波动项的月度事故时间序列；YtT是原事故序列中的长期稳定性趋势项；YtC是其中短期周期性波动项；原事故序列可由公式（2）表示，

通过HP滤波可以对短期波动项进行约束，将长期稳定性趋势项YtT分离出来。HP滤波过程是使损失函数（3）达到最小，

其中：λ为调节参数，λ越大，估计趋势越光滑；λ趋于无穷大时，估计趋势将接近线性函数，通常对于月度数据λ取14400。采用HP滤波方法对事故序列进行了趋势周期分解，分解结果，如图2。可知，原始序列Yt围绕长期稳定趋势序列YtT波动，采用最小二乘法对长期稳定趋势序列进行拟合，得到长期稳定趋势序列方程：

上述模型估计结果的调整R2值为0.97，说明模型的拟合度较好；常数项和时间趋势项系数估计值均达到1%显著性水平。至此，完成了对月度生产安全事故数序列的长期趋势项分解和拟合。分解结果表明，我国生产安全事故月度序列具有稳定的长期变化趋势，该趋势序列符合时间的一元线性回归模型，其中时间项系数为-0.851表明造成人员死亡（包括下落不明）的较大及以上事故起数随着月份变化具有一个显著的下降趋势特征。该系数表明，不考虑短期波动成分的影响下，造成人员死亡的较大及以上生产安全事故数以每月接近1起的速度下降。同时，由图2可以直观的看出，该下降速度逐渐变缓，说明从长期来看事故起数将围绕一个固定值波动变化（说明我国总体安全生产状况趋于稳定），该值将受国内安全生产技术水平所控制。

图2 H P滤波后的趋势项YtTFig.2 Trend composition after HP filter

3 短期波动周期分析

将月度事故原始序列去除长期稳定下降趋势项后得到短周期波动趋势，如图3。初步判断该序列为平稳序列，ADF单位根检验结果表明，可以显著拒绝原假设，即YtC序列为平稳时间序列。将信号处理中的谐波分析方法用于短周期项波动特征的研究。信号处理原理认为，有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征，如频率、幅值、初相位。傅立叶原理表明[14-15]，任何连续测量的时序或信号，都可以表示为不同频率的余弦（或正弦）波信号的无限叠加，即：

式中：

t—时间；

Fn—频率；

图3 H P滤波后的波动项YtCFig.3 Volatile composition after HP filter

a0—直流分量；

Pn—相位角；

An—幅值，起。

图4 月度事故幅值频谱Fig.4 Monthly accident amplitude spectrum

图5 月度事故相位频谱Fig.5 Monthly accident phase spectrum

表2 月度事故序列傅里叶变换结果Tab.2 Fourier transform sequence of monthly accidents

采用Matlab编程[16]，对时域上的事故波动序列进行傅里叶变换可得到对应的幅值谱和相位谱。时域上的造成死亡的较大及以上生产安全事故序列经变换后得到的41个余弦波信号所对应的频率和幅值，其中每一个点对应着一个频率和幅度特性。根据傅里叶变换理论，叠加的余弦波数量越多越能接近原始曲线，但叠加波数是无穷多的，是无法实现的。变换得到的幅值谱，如图4，其中直流分量和幅值较大的8个余弦波信号已在图中标明（每一点括号中的数据分别为特征频率和幅值，例如频率为0.0319的余弦波的幅值为6.2759），直流分量和余弦波分量对月度事故序列的短周期波动影响较大。直流分量和余弦波信号的相位谱，如图5，图中标明的9组数据分别为上述直流分量和幅值较大的8个余弦波信号对应的特征频率和相位，例如上述频率为0.0319的余弦波的初相位为114.0149°。可见，通过傅里叶变换将时域上的连续事故波动序列变换为频域上的非周期离散的函数。图4、5中每一点代表了一个余弦波的频率、幅值和相位，将这些数据代入公式（4）便可累加得到事故序列的短周期波动项。

将傅里叶变换结果汇总于表2。表2中所示为各余弦波和直流分量按照幅值从大向小顺序排列结果，前1～8行为8个余弦波波形参数，最后一行为直流分量参数。在幅值排列前8个余弦波中，有2个幅值较大，分别为：频率为0.0319×10-6Hz，幅值为6.2759，初相位为114.01°，对应的事故周期为11.7个月；另外1个的频率为0.0911×10-6Hz，幅值为4.7811，初相位为163.51°的余弦波，对应的事故周期为4.1个月；该2个余弦波构成了对事故波动项影响较大的一部分。虽然从第5至8个余弦波的波形可以看出，在一年中造成死亡的较大及以上生产安全事故波动序列存在一个相对减弱的周期，即2.1～4.8个月，但其与上述4.1个月的强周期基本重合，因此忽略其对波动的影响。另外从长期来看，事故序列还存在2个时间较长但影响较弱的周期即27.3～41个月，约3年时间。表2中最后一行为直流分量，可以看出该直流分量频率为0、初相位角为0°反映在函数中直流分量，幅值较小对事故变化波动影响较小，可忽略。

利用直流分量和影响较大的前3个余弦分量据将月度生产安全事故的波动项表示为傅里叶级数形式，即：

4 趋势分析

通过研究得到，造成死亡的较大及以上事故月度事故序列包含2种不同的变化趋势成分，即长期稳定变化趋势和短期波动变化趋势。长期稳定变化趋势成分研究表明，事故序列属于趋势平稳序列，并以平均每月约1起的速度稳定下降；而短期波动变化趋势围绕长期趋势线波动，同样属于平稳序列，通过傅里叶变换得到的事故波动序列在频域上的傅里叶级数展开形式表明，造成死亡的较大及以上生产安全事故波动变化趋势存在3个较明显的周期，即4.1个月、11.7个月和27.3～41个月。

5 结论

基于造成死亡的较大及以上生产安全事故月度数据，采用Hodrick-Prescott滤波、线性回归和谐波分析中的快速傅里叶变换等3种方法，将事故时间序列分解为长周期稳定项和短周期波动项并分别建立了特征模型，得到以下结论：

（1）月度生产安全事故统计数据在时间域上表现为趋势平稳序列，由于受到多重因素影响其变化过程隐含着长周期稳定趋势和短周期性波动变化等2种趋势。

（2）通过HP滤波分解出长周期项后，通过拟合发现造成死亡的较大及以上的月度事故序列符合时间（月份）的一元线性回归模型，回归模型能够较好的解释事故稳定下降的趋势。

（3）通过谐波分析方法中的傅里叶变换，将事故序列中的短期波动项展开成傅里叶级数形式，得到事故序列的频率—幅值谱和频率—相位谱，发现事故短期波动具有3种主要周期成份，其一为影响较强的4.1、11.7个月，其二为影响相对较弱的27.3～41个月。