顏 峻

（中國勞動(dòng)關(guān)系學(xué)院 安全工程系，北京 100048）

0 引言

我國安全生產(chǎn)形勢(shì)呈現(xiàn)出總體穩(wěn)定、趨向好轉(zhuǎn)的發(fā)展態(tài)勢(shì)，但在事故指標(biāo)不斷下降的同時(shí)，仍存在著反復(fù)波動(dòng)的現(xiàn)象。這一現(xiàn)象在事故時(shí)間序列上表現(xiàn)為長(zhǎng)周期穩(wěn)定變化和短期波動(dòng)變化2種趨勢(shì)。若不將2種趨勢(shì)加以分解，將會(huì)影響對(duì)安全生產(chǎn)事故趨勢(shì)的精準(zhǔn)判斷。因此，本文的主要目的是將原本耦合在一起的2種趨勢(shì)加以分解，分別在時(shí)域和頻域上構(gòu)建描述我國生產(chǎn)安全事故變化趨勢(shì)的特征模型。

事故時(shí)序分析是指基于歷史數(shù)據(jù)統(tǒng)計(jì)的事故時(shí)間序列變化規(guī)律的研究，主要方法包括概率分布函數(shù)、自回歸移動(dòng)平均模型、支持向量機(jī)、相關(guān)向量機(jī)、人工神經(jīng)元網(wǎng)絡(luò)、灰色模型等。到目前為止，研究人員已經(jīng)初步發(fā)現(xiàn)事故時(shí)間序列存在著不同的波動(dòng)特征。文獻(xiàn)[1]采用季節(jié)指數(shù)法分析了月度事故序列中存在的季節(jié)因素，通過建立不同時(shí)段事故數(shù)的概率分布模型，研究事故發(fā)生的時(shí)間段所具有的周期特征。文獻(xiàn)[2]建立了礦山安全事故季節(jié)周期回歸模型，得到事故發(fā)生的趨勢(shì)線（即未來事故的平均數(shù)），事故數(shù)圍繞這一趨勢(shì)線上下波動(dòng)。在趨勢(shì)預(yù)測(cè)方面，由于ARⅠMA模型可以有效的處理非平穩(wěn)時(shí)間序列，因此被廣泛用于事故趨勢(shì)分析和預(yù)測(cè)。文獻(xiàn)[3]針對(duì)航空裝備事故時(shí)序數(shù)據(jù)具有的非平穩(wěn)性、自相關(guān)性特征，構(gòu)建了事故時(shí)序數(shù)列的ARⅠMA（4,1,4）模型，通過一次差分將非平穩(wěn)數(shù)列轉(zhuǎn)換為平穩(wěn)數(shù)列對(duì)事故數(shù)進(jìn)行了預(yù)測(cè)。文獻(xiàn)[4]建立了道路交通事故ARⅠMA（0,1,1）預(yù)測(cè)模型，發(fā)現(xiàn)時(shí)序數(shù)列殘差具有較為明顯的季節(jié)性波動(dòng)趨勢(shì)。為了消除季節(jié)性因素的影響，文獻(xiàn)[5]采用季節(jié)調(diào)整X-12方法去消除事故月度數(shù)據(jù)中的季節(jié)性影響因素后，建立了針對(duì)新序列的灰色GM（1,1）模型和隨機(jī)波動(dòng)項(xiàng)的ARⅠMA模型。文獻(xiàn)[6]采用相關(guān)向量機(jī)（RVM）方法建立了交通事故時(shí)間序列預(yù)測(cè)模型，研究了數(shù)據(jù)中具有的線性關(guān)系和非線性關(guān)系。文獻(xiàn)[7]采用人工神經(jīng)元網(wǎng)絡(luò)方法對(duì)事故趨勢(shì)進(jìn)行了分析，發(fā)現(xiàn)事故變化趨勢(shì)具有長(zhǎng)期穩(wěn)定和短期波動(dòng)2種特征，尤其是月度數(shù)據(jù)序列具有明顯的周期性特征。此外，灰色預(yù)測(cè)模型也被廣泛用于事故趨勢(shì)預(yù)測(cè)，但該方法并沒有區(qū)分時(shí)序數(shù)列包含的長(zhǎng)期趨勢(shì)項(xiàng)和短期波動(dòng)項(xiàng)[8-9]。

綜上所述，生產(chǎn)安全事故時(shí)間序列具有明顯的非平穩(wěn)性特征，同時(shí)受多種因素影響而呈現(xiàn)出2種不同特征項(xiàng)，即長(zhǎng)期穩(wěn)定趨勢(shì)和短期波動(dòng)趨勢(shì)，兩者交織在一起增加了變化規(guī)律研究和預(yù)測(cè)的難度。已有研究表明，事故序列具有一定的周期性，但大多難以分解2種特征趨勢(shì)。雖然X-12等季節(jié)調(diào)整方法可研究事故序列具有的季節(jié)波動(dòng)特性，但仍將長(zhǎng)期平穩(wěn)趨勢(shì)和短期波動(dòng)趨勢(shì)視為一體不能分開。本文將從生產(chǎn)安全事故月度時(shí)間序列變化趨勢(shì)研究著手，對(duì)序列的長(zhǎng)期趨勢(shì)項(xiàng)和短期波動(dòng)項(xiàng)加以分解研究。將采用Hodrick-Prescott濾波方法提取事故序列中的長(zhǎng)期趨勢(shì)項(xiàng)，采用線性回歸模型確定事故變化的長(zhǎng)期趨勢(shì)特征；對(duì)分解后的短期波動(dòng)項(xiàng)，采用諧波分析中的快速傅里葉變換將時(shí)域序列轉(zhuǎn)換到頻域上，以研究事故變化趨勢(shì)在短期上具有的周期性變化特征。

1 事故序列平穩(wěn)性檢驗(yàn)

收集我國自2011年1月至2017年11月造成人員死亡（包括下落不明）的較大及以上生產(chǎn)安全事故月度統(tǒng)計(jì)數(shù)據(jù)，繪制時(shí)間序列，如圖1[10]。月度事故序列存在長(zhǎng)期波動(dòng)下降趨勢(shì)，為非平穩(wěn)時(shí)間序列。為了進(jìn)一步驗(yàn)證序列具有的非平穩(wěn)特征，利用E-views[11]軟件中的ADF Fisher單位根檢驗(yàn)方法[12]對(duì)月度事故序列進(jìn)行平穩(wěn)性檢驗(yàn)，其檢驗(yàn)?zāi)Ｐ蜑椋?/p>

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

式中：

t—表示序列的時(shí)間變化趨勢(shì)，月；

α—常數(shù)項(xiàng)；

β—時(shí)間項(xiàng)系數(shù)；

δ—前一期變量系數(shù)；模型假設(shè)均為H0：δ=0，即存在單位根。從圖1可以看出，事故序列具有明顯的隨時(shí)間下降趨勢(shì)，因此檢驗(yàn)?zāi)Ｐ蛻?yīng)包含時(shí)間趨勢(shì)項(xiàng)βt，平穩(wěn)性檢驗(yàn)結(jié)果，見表1。

表1 月度事故序列平穩(wěn)性檢驗(yàn)（水平值）Tab.1 Monthly accident sequence stationarity test（level）

表1為包含時(shí)間趨勢(shì)項(xiàng)的ADF單位根檢驗(yàn)結(jié)果，檢驗(yàn)結(jié)果為穩(wěn)定序列。如果把檢驗(yàn)?zāi)Ｐ椭械臅r(shí)間趨勢(shì)項(xiàng)去掉后進(jìn)行平穩(wěn)性檢驗(yàn)，檢驗(yàn)結(jié)果不能拒絕原假設(shè)，即事故序列為非平穩(wěn)序列。綜合兩種不同形式的檢驗(yàn)?zāi)Ｐ退闷椒€(wěn)性檢驗(yàn)結(jié)論，月度事故序列為趨勢(shì)平穩(wěn)序列，即含有顯著的時(shí)間趨勢(shì)項(xiàng)。

2 長(zhǎng)期穩(wěn)定趨勢(shì)分析

上述“事故時(shí)間序列為趨勢(shì)平穩(wěn)序列”的檢驗(yàn)結(jié)果證明月度事故序列具有長(zhǎng)期穩(wěn)定趨勢(shì)項(xiàng)。采用Hodrick-Prescott濾波方法[13]對(duì)事故序列中具有的長(zhǎng)期趨勢(shì)項(xiàng)和短期周期性波動(dòng)項(xiàng)進(jìn)行分解。設(shè)Yt為包含長(zhǎng)期穩(wěn)定性趨勢(shì)項(xiàng)和短期周期性波動(dòng)項(xiàng)的月度事故時(shí)間序列；YtT是原事故序列中的長(zhǎng)期穩(wěn)定性趨勢(shì)項(xiàng)；YtC是其中短期周期性波動(dòng)項(xiàng)；原事故序列可由公式（2）表示，

通過HP濾波可以對(duì)短期波動(dòng)項(xiàng)進(jìn)行約束，將長(zhǎng)期穩(wěn)定性趨勢(shì)項(xiàng)YtT分離出來。HP濾波過程是使損失函數(shù)（3）達(dá)到最小，

其中：λ為調(diào)節(jié)參數(shù)，λ越大，估計(jì)趨勢(shì)越光滑；λ趨于無窮大時(shí)，估計(jì)趨勢(shì)將接近線性函數(shù)，通常對(duì)于月度數(shù)據(jù)λ取14400。采用HP濾波方法對(duì)事故序列進(jìn)行了趨勢(shì)周期分解，分解結(jié)果，如圖2?？芍?，原始序列Yt圍繞長(zhǎng)期穩(wěn)定趨勢(shì)序列YtT波動(dòng)，采用最小二乘法對(duì)長(zhǎng)期穩(wěn)定趨勢(shì)序列進(jìn)行擬合，得到長(zhǎng)期穩(wěn)定趨勢(shì)序列方程：

上述模型估計(jì)結(jié)果的調(diào)整R2值為0.97，說明模型的擬合度較好；常數(shù)項(xiàng)和時(shí)間趨勢(shì)項(xiàng)系數(shù)估計(jì)值均達(dá)到1%顯著性水平。至此，完成了對(duì)月度生產(chǎn)安全事故數(shù)序列的長(zhǎng)期趨勢(shì)項(xiàng)分解和擬合。分解結(jié)果表明，我國生產(chǎn)安全事故月度序列具有穩(wěn)定的長(zhǎng)期變化趨勢(shì)，該趨勢(shì)序列符合時(shí)間的一元線性回歸模型，其中時(shí)間項(xiàng)系數(shù)為-0.851表明造成人員死亡（包括下落不明）的較大及以上事故起數(shù)隨著月份變化具有一個(gè)顯著的下降趨勢(shì)特征。該系數(shù)表明，不考慮短期波動(dòng)成分的影響下，造成人員死亡的較大及以上生產(chǎn)安全事故數(shù)以每月接近1起的速度下降。同時(shí)，由圖2可以直觀的看出，該下降速度逐漸變緩，說明從長(zhǎng)期來看事故起數(shù)將圍繞一個(gè)固定值波動(dòng)變化（說明我國總體安全生產(chǎn)狀況趨于穩(wěn)定），該值將受國內(nèi)安全生產(chǎn)技術(shù)水平所控制。

圖2 H P濾波后的趨勢(shì)項(xiàng)YtTFig.2 Trend composition after HP filter

3 短期波動(dòng)周期分析

將月度事故原始序列去除長(zhǎng)期穩(wěn)定下降趨勢(shì)項(xiàng)后得到短周期波動(dòng)趨勢(shì)，如圖3。初步判斷該序列為平穩(wěn)序列，ADF單位根檢驗(yàn)結(jié)果表明，可以顯著拒絕原假設(shè)，即YtC序列為平穩(wěn)時(shí)間序列。將信號(hào)處理中的諧波分析方法用于短周期項(xiàng)波動(dòng)特征的研究。信號(hào)處理原理認(rèn)為，有些信號(hào)在時(shí)域上是很難看出什么特征的，但是如果變換到頻域之后，就很容易看出特征，如頻率、幅值、初相位。傅立葉原理表明[14-15]，任何連續(xù)測(cè)量的時(shí)序或信號(hào)，都可以表示為不同頻率的余弦（或正弦）波信號(hào)的無限疊加，即：

式中：

t—時(shí)間；

Fn—頻率；

圖3 H P濾波后的波動(dòng)項(xiàng)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 月度事故序列傅里葉變換結(jié)果Tab.2 Fourier transform sequence of monthly accidents

采用Matlab編程[16]，對(duì)時(shí)域上的事故波動(dòng)序列進(jìn)行傅里葉變換可得到對(duì)應(yīng)的幅值譜和相位譜。時(shí)域上的造成死亡的較大及以上生產(chǎn)安全事故序列經(jīng)變換后得到的41個(gè)余弦波信號(hào)所對(duì)應(yīng)的頻率和幅值，其中每一個(gè)點(diǎn)對(duì)應(yīng)著一個(gè)頻率和幅度特性。根據(jù)傅里葉變換理論，疊加的余弦波數(shù)量越多越能接近原始曲線，但疊加波數(shù)是無窮多的，是無法實(shí)現(xiàn)的。變換得到的幅值譜，如圖4，其中直流分量和幅值較大的8個(gè)余弦波信號(hào)已在圖中標(biāo)明（每一點(diǎn)括號(hào)中的數(shù)據(jù)分別為特征頻率和幅值，例如頻率為0.0319的余弦波的幅值為6.2759），直流分量和余弦波分量對(duì)月度事故序列的短周期波動(dòng)影響較大。直流分量和余弦波信號(hào)的相位譜，如圖5，圖中標(biāo)明的9組數(shù)據(jù)分別為上述直流分量和幅值較大的8個(gè)余弦波信號(hào)對(duì)應(yīng)的特征頻率和相位，例如上述頻率為0.0319的余弦波的初相位為114.0149°?？梢姡ㄟ^傅里葉變換將時(shí)域上的連續(xù)事故波動(dòng)序列變換為頻域上的非周期離散的函數(shù)。圖4、5中每一點(diǎn)代表了一個(gè)余弦波的頻率、幅值和相位，將這些數(shù)據(jù)代入公式（4）便可累加得到事故序列的短周期波動(dòng)項(xiàng)。

將傅里葉變換結(jié)果匯總于表2。表2中所示為各余弦波和直流分量按照幅值從大向小順序排列結(jié)果，前1～8行為8個(gè)余弦波波形參數(shù)，最后一行為直流分量參數(shù)。在幅值排列前8個(gè)余弦波中，有2個(gè)幅值較大，分別為：頻率為0.0319×10-6Hz，幅值為6.2759，初相位為114.01°，對(duì)應(yīng)的事故周期為11.7個(gè)月；另外1個(gè)的頻率為0.0911×10-6Hz，幅值為4.7811，初相位為163.51°的余弦波，對(duì)應(yīng)的事故周期為4.1個(gè)月；該2個(gè)余弦波構(gòu)成了對(duì)事故波動(dòng)項(xiàng)影響較大的一部分。雖然從第5至8個(gè)余弦波的波形可以看出，在一年中造成死亡的較大及以上生產(chǎn)安全事故波動(dòng)序列存在一個(gè)相對(duì)減弱的周期，即2.1～4.8個(gè)月，但其與上述4.1個(gè)月的強(qiáng)周期基本重合，因此忽略其對(duì)波動(dòng)的影響。另外從長(zhǎng)期來看，事故序列還存在2個(gè)時(shí)間較長(zhǎng)但影響較弱的周期即27.3～41個(gè)月，約3年時(shí)間。表2中最后一行為直流分量，可以看出該直流分量頻率為0、初相位角為0°反映在函數(shù)中直流分量，幅值較小對(duì)事故變化波動(dòng)影響較小，可忽略。

利用直流分量和影響較大的前3個(gè)余弦分量據(jù)將月度生產(chǎn)安全事故的波動(dòng)項(xiàng)表示為傅里葉級(jí)數(shù)形式，即：

4 趨勢(shì)分析

通過研究得到，造成死亡的較大及以上事故月度事故序列包含2種不同的變化趨勢(shì)成分，即長(zhǎng)期穩(wěn)定變化趨勢(shì)和短期波動(dòng)變化趨勢(shì)。長(zhǎng)期穩(wěn)定變化趨勢(shì)成分研究表明，事故序列屬于趨勢(shì)平穩(wěn)序列，并以平均每月約1起的速度穩(wěn)定下降；而短期波動(dòng)變化趨勢(shì)圍繞長(zhǎng)期趨勢(shì)線波動(dòng)，同樣屬于平穩(wěn)序列，通過傅里葉變換得到的事故波動(dòng)序列在頻域上的傅里葉級(jí)數(shù)展開形式表明，造成死亡的較大及以上生產(chǎn)安全事故波動(dòng)變化趨勢(shì)存在3個(gè)較明顯的周期，即4.1個(gè)月、11.7個(gè)月和27.3～41個(gè)月。

5 結(jié)論

基于造成死亡的較大及以上生產(chǎn)安全事故月度數(shù)據(jù)，采用Hodrick-Prescott濾波、線性回歸和諧波分析中的快速傅里葉變換等3種方法，將事故時(shí)間序列分解為長(zhǎng)周期穩(wěn)定項(xiàng)和短周期波動(dòng)項(xiàng)并分別建立了特征模型，得到以下結(jié)論：

（1）月度生產(chǎn)安全事故統(tǒng)計(jì)數(shù)據(jù)在時(shí)間域上表現(xiàn)為趨勢(shì)平穩(wěn)序列，由于受到多重因素影響其變化過程隱含著長(zhǎng)周期穩(wěn)定趨勢(shì)和短周期性波動(dòng)變化等2種趨勢(shì)。

（2）通過HP濾波分解出長(zhǎng)周期項(xiàng)后，通過擬合發(fā)現(xiàn)造成死亡的較大及以上的月度事故序列符合時(shí)間（月份）的一元線性回歸模型，回歸模型能夠較好的解釋事故穩(wěn)定下降的趨勢(shì)。

（3）通過諧波分析方法中的傅里葉變換，將事故序列中的短期波動(dòng)項(xiàng)展開成傅里葉級(jí)數(shù)形式，得到事故序列的頻率—幅值譜和頻率—相位譜，發(fā)現(xiàn)事故短期波動(dòng)具有3種主要周期成份，其一為影響較強(qiáng)的4.1、11.7個(gè)月，其二為影響相對(duì)較弱的27.3～41個(gè)月。