基于時間序列諧波分析的死亡事故總量波動趨勢分析

顏 峻

（中國勞動關系學院 安全工程系，北京 100048）

0 引言

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

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

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

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

收集我國自2011年1月至2017年11月造成人員死亡（包括下落不明）的較大及以上生產安全事故月度統計數據，繪制時間序列，如圖1[10]。月度事故序列存在長期波動下降趨勢，為非平穩(wěn)時間序列。為了進一步驗證序列具有的非平穩(wěn)特征，利用E-views[11]軟件中的ADF Fisher單位根檢驗方法[12]對月度事故序列進行平穩(wěn)性檢驗，其檢驗模型為：

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

式中：

t—表示序列的時間變化趨勢，月；

α—常數項；

β—時間項系數；

δ—前一期變量系數；模型假設均為H0：δ=0，即存在單位根。從圖1可以看出，事故序列具有明顯的隨時間下降趨勢，因此檢驗模型應包含時間趨勢項βt，平穩(wěn)性檢驗結果，見表1。

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

表1為包含時間趨勢項的ADF單位根檢驗結果，檢驗結果為穩(wěn)定序列。如果把檢驗模型中的時間趨勢項去掉后進行平穩(wěn)性檢驗，檢驗結果不能拒絕原假設，即事故序列為非平穩(wěn)序列。綜合兩種不同形式的檢驗模型所得平穩(wěn)性檢驗結論，月度事故序列為趨勢平穩(wěn)序列，即含有顯著的時間趨勢項。

2 長期穩(wěn)定趨勢分析

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

通過HP濾波可以對短期波動項進行約束，將長期穩(wěn)定性趨勢項YtT分離出來。HP濾波過程是使損失函數（3）達到最小，

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

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

圖2 H P濾波后的趨勢項YtTFig.2 Trend composition after HP filter

3 短期波動周期分析

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

式中：

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°?？梢姡ㄟ^傅里葉變換將時域上的連續(xù)事故波動序列變換為頻域上的非周期離散的函數。圖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種不同的變化趨勢成分，即長期穩(wěn)定變化趨勢和短期波動變化趨勢。長期穩(wěn)定變化趨勢成分研究表明，事故序列屬于趨勢平穩(wěn)序列，并以平均每月約1起的速度穩(wěn)定下降；而短期波動變化趨勢圍繞長期趨勢線波動，同樣屬于平穩(wěn)序列，通過傅里葉變換得到的事故波動序列在頻域上的傅里葉級數展開形式表明，造成死亡的較大及以上生產安全事故波動變化趨勢存在3個較明顯的周期，即4.1個月、11.7個月和27.3～41個月。

5 結論

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

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

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

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