侯秦脈 朱 偉 鄒 象 劉時賢 吳彥農(nóng)

（生態(tài)環(huán)境部核與輻射安全中心 北京 102445）

《中國的核安全》首次發(fā)布于2019年，該書是中國政府發(fā)表的第一部核安全白皮書，白皮書全面介紹了中國核安全事業(yè)的發(fā)展歷程，闡述中國核安全的基本原則和政策主張，分享中國核安全監(jiān)管的理念和實踐，闡明中國推進全球核安全治理進程的決心和行動［1］。截至2019年6月，已安全穩(wěn)定運行累計300余堆年，未發(fā)生過國際核與放射事件分級表（International Nuclear Event Scale，INES）2級及以上的事件或事故，且0級偏差和1級異常事件發(fā)生率呈下降趨勢。該報告提供了1991～2018年我國運行核電廠的運行事件數(shù)量，同時《核安全年報》也提供了2019～2020年的運行事件數(shù)量［2］，此領域研究目前多為對運行事件本身的分析和一般的趨勢性描述。鄭麗馨等［3］指出運行事件對了解機組性能、運行經(jīng)驗趨勢的評價有重要的貢獻。胡攀［4］認為對低級別事件趨勢分析能夠跟蹤管理異常趨勢并提高異常趨勢的識別效率。劉時賢等［5］從誤碰方面分析了運行事件對核電廠的安全性和經(jīng)濟性造成影響。修炳林等［6］指出2015～2016年新機組投運高峰帶來運行事件的增加。但缺乏有效的工具對運行事件數(shù)量變化趨勢及預測進行研究，因此，需要研究一套方法進行核電廠運行事件數(shù)量的趨勢分析和預測。

差分自回歸移動平均模型（Autoregressive Mobile Average Model，ARIMA）是經(jīng)典的時間序列預測方法，武云云等［7］探討了運用ARIMA模型在核電站外圍環(huán)境放射性水平動態(tài)分析和預測，為核電站外圍環(huán)境影響評價和核事故放射性污染提供參考。胡衍坤等［8］指出ARIMA模型能夠較好地體現(xiàn)時間序列數(shù)據(jù)中的線性特征，但是單一的ARIMA模型對非線性變化難以充分有效地處理，結合網(wǎng)絡結構模型能有效提高河流水質(zhì)預測精度7%。在深度學習算法中，長短期記憶神經(jīng)網(wǎng)絡模型（Long and Short-Term Memory，LSTM）由于其特殊的網(wǎng)絡結構，在處理事件序列問題時，比傳統(tǒng)的神經(jīng)網(wǎng)絡更快更容易收斂到最優(yōu)解，非常適合處理核電廠運行事件指標這種時序數(shù)據(jù)。因此，本文建立ARIMA與LSTM組合模型進行核電廠運行事件預測的研究，通過Mann-Kendall檢驗法（簡稱M-K檢驗法）探討運行事件變化特性以及預測的合理性，并結合白皮書中我國運行核電廠運行事件數(shù)據(jù)，分析運行事件變化趨勢，為核安全監(jiān)管提供數(shù)據(jù)支持。

1 研究方法描述

1.1 M-K檢驗法

M-K檢驗法是一種非參數(shù)檢驗方法［9］，對樣本分布無要求，能有效排除個別異常值對時間序列的影響，該方法廣泛應用于大氣和地質(zhì)等行業(yè)時間序列突變點檢測研究中［10］。在突變分析方面，M-K檢驗能夠查找出數(shù)據(jù)異常點，確定突變區(qū)域。

根據(jù)核電廠運行事件數(shù)時間序列的秩序列為：

在時間序列隨機獨立的假定下，定義統(tǒng)計量：

式中：E（Sk）、Var（Sk）分別是秩序列Sk的均值和方差。

根據(jù)UBk=?UFk，同樣方法計算反序列得到統(tǒng)計量UBk。依據(jù)UBk和UFk可得到UB曲線和UF曲線，差值大于零則呈下降趨勢，小于零呈上升趨勢，差值等于零為相交點［11?12］，如果相交點在0.05顯著性水平下，則該點是突變點。

1.2 ARIMA模型

ARIMA是一種廣泛用于時間序列分析和預測的模型，對處理線性時間序列有較好的優(yōu)勢［13］。非平穩(wěn)時間序列通過該方法差分后轉化為平穩(wěn)時間序列，模型對因變量滯后值以及其隨機誤差進行回歸，融合了兩種模型：自回歸模型（Autoregressive，AR）和移動平均模型（Moving Average，MA）。組間隨機變量具有的依存關系，既包括本身的變化規(guī)律同時也包括了外部影響因素，該方法有較高的預測精度，有助于解釋預測變動的規(guī)律。

首選需要確保預測時，ARIMA模型的時間序列是平穩(wěn)的或近似平穩(wěn)的，自相關系數(shù)只有時間間隔一個變量，即時間變化不影響均數(shù)和方差。基于同質(zhì)性是非平穩(wěn)時間序列的普遍特征，可以通過差分處理轉化為平穩(wěn)時間序列，差分次數(shù)就是其次的階數(shù)。那么完整的差分自回歸移動平均過程ARIMA（p，d，q）模型可以用滯后算子B表示為［14］：

式中：p是自回歸滯后階；d是非平穩(wěn)序列轉化為平穩(wěn)序列的差分次數(shù)；q是移動平均階數(shù)，偏自相關函數(shù)值顯著不為零的最高階數(shù)；ei是i時刻的偏差；Фi是自回歸部分的系數(shù)；B是滯后算子；θi是移動平均部分的系數(shù)［15］。

ARIMA模型建立方法如下：

1）模型識別，驗證時間序列數(shù)據(jù)平穩(wěn)性，非平穩(wěn)數(shù)據(jù)需進行差分處理，利用自相關函數(shù)圖（Autocorrelation Function，ACF）和偏自相關函數(shù)圖（Partial Autocorrelation Function，PACF）進行分析，確定自回歸滯后階數(shù)、差分次數(shù)和移動平均階數(shù)。

2）參數(shù)估計和模型檢驗。用特定時間序列的數(shù)據(jù)，估計模型的參數(shù)并進行檢驗［16］，判斷該模型是否恰當，如不恰當，返回步驟1）。

3）趨勢預測應用，并對設定的將來某個時間段給出預測值。

1.3 LSTM模型

LSTM是一種遞歸神經(jīng)網(wǎng)絡（Recursive Neural Network，RNNs）的改進［17］，可以解決長時間周期依賴的預測，并且對于有限數(shù)據(jù)樣本深入分析其長期趨勢，解決RNNs距離感知能力對梯度的依賴，在各行業(yè)時間序列預測上LSTM有著廣泛的應用［18］。

LSTM模型在遞歸神經(jīng)網(wǎng)絡的基礎上添加了遺忘、輸入、輸出等三個門層，解決了有限時間和數(shù)據(jù)樣本的問題［19］，通過梯度下降得到模型的最佳誤差參數(shù)。LSTM的第一步是確定遺忘門層遺忘哪些信息，借助Sigmoid層決定前一時間周期的節(jié)點數(shù)據(jù)。讀取上一層節(jié)點數(shù)據(jù)輸出ht?1和輸入數(shù)據(jù)Xt，并為單元狀態(tài)Ct?1輸出0~1的值ft，0和1分別表示“不允許任何量通過”和“允許任何量通過”的邏輯意義［20］。ft通過Sigmoid函數(shù)的輸出值得到，分配給Ct?1決定前一個狀態(tài)的值是否丟棄，ft為0表示“遺忘”，1表示“保留”，Wf是權重矩陣，bf是偏置項，其表達式為：

下一步是在輸入層中決定輸入哪些信息。Sigmoid函數(shù)決定將更新哪些值，添加到新候選值Ct的tanh函數(shù)輸入值（?1~1）組成一個向量，并乘以Sigmoid函數(shù)值，得到輸出值。其表達式為：

最后需要通過輸出門層輸出函數(shù)的參數(shù)。首先運行一個Sigmoid層，將得到輸出參數(shù)信息ot，包含決定輸出的單元狀態(tài)信息［21］，最終得到將要輸出的部分ht，其表達式為：

1.4 組合預測模型

核電廠運行事件數(shù)據(jù)具有線性特征和非線性特征，將運行事件數(shù)據(jù)記為At，分解為線性和非線性兩個部分，分別用Bt表示線性部分和Ct表示非線性部分，表達式如下：

第一步，時間序列模型建模；對時間數(shù)據(jù)序列At進行建模，然后將訓練數(shù)據(jù)進行擬合并進行預測，通過自相關圖和偏自相關圖建立ARIMA（p，d，q）模型，輸出時間數(shù)據(jù)序列Bt'，其數(shù)據(jù)序列描述線性部分，得到偏差et的時間數(shù)據(jù)序列。

第二步，LSTM殘差預測；使用LSTM模型對偏差的時間數(shù)據(jù)序列et進行擬合。對數(shù)據(jù)序列At中非線性部分數(shù)據(jù)信息進行挖掘，識別Ct中與et的異常關系，通過對et進行LSTM處理，得到預測結果Ct'。

第三步，模型組合；通過ARIMA（p，d，q）模型表達線性部分和LSTM模型表達非線性部分建立混合模型，得到混合預測模型表達式。

因此，將刻畫數(shù)據(jù)序列線性關系的ARIMA（p，d，q）模型和描述數(shù)據(jù)序列非線性關系的LSTM模型的各自優(yōu)點有機結合，能夠顯示混合模型的優(yōu)越性和提高預測的準確度。

2 實驗數(shù)據(jù)和變化特征分析

2.1 實驗數(shù)據(jù)及趨勢分析

實驗運行事件數(shù)量數(shù)據(jù)來源于2019年9月國務院新聞辦公室發(fā)表的《中國的核安全》白皮書，白皮書提供了1991～2018年我國運行核電廠的運行事件數(shù)量。國家核安全年報提供了2019～2020年的運行事件數(shù)量，將用作組合模型預測值的對比。

采用線性回歸的方法對1991～2018年運行事件數(shù)量時間序列進行線性趨勢分析（圖1），運行事件數(shù)量傾向率為線性方程的斜率，在0.95置信限下運行事件數(shù)量趨勢函數(shù)為f（x）=69.252?0.017 5x。由圖1可以看出，1）運行事件數(shù)量與時間呈負線性相關趨勢。2）年度運行事件波動性變化明顯，波動范圍在8~93，1993年、2002年和2016年達到波峰，1991年、2000年和2009年達到波谷，經(jīng)歷了枯-豐-正常-枯-豐-正常-枯-豐階段，最大值與最小值相差10倍。

圖1 1991～2018年運行事件數(shù)量及趨勢Fig.1 Number and trend of operation events from 1991 to 2018

2.2 運行事件數(shù)量時間序列突變分析

采用M-K方法檢驗突變點，對1991～2018年我國運行核電廠的運行事件數(shù)量進行檢驗，檢驗結果見圖2。結果表明：年度運行事件的突變檢驗曲線存在交匯點，圖2中時間序列M-K突變檢驗在0.05顯著性水平下，UF統(tǒng)計量與UB統(tǒng)計量從1995～2013年無交點即無突變點，這一時期我國核電機組數(shù)量以年均不到1臺機組緩慢發(fā)展，運行事件數(shù)量無顯著突變。

圖2 運行事件數(shù)量M-K檢驗結果Fig.2 M-K test results of the number of operation events

而在1991年、1994年、2014年、2016年和2017年存在交點，說明1991～1994年以及2014～2018年為運行事件序列顯著突變年份，查閱核電機組數(shù)量從1991～1994年核電機組從無到有的摸索階段，同時2014～2018年我國核電機組數(shù)量年均增長4.2臺，對于核電行業(yè)來說這兩個階段都是快速發(fā)展期，處在新人員積累經(jīng)驗的階段，可以說明突變點伴隨著行業(yè)的高速發(fā)展，同時也說明機組數(shù)量年均增長1臺不會引起運行事件數(shù)量突變，而機組數(shù)量年均增長4臺會引起國內(nèi)運行事件數(shù)量的突變，應對方法將是在高速發(fā)展階段前提前儲備人才，加強人員職業(yè)能力培養(yǎng)，核安全監(jiān)管隊伍規(guī)模也應該做出相應調(diào)整。

3 運行事件時間序列預測

從文獻［8］中可以看出，客觀事物具有線性和非線性特征，本文§2趨勢分析和突變分析反映了運行事件也具有線性和非線性的雙重特征。對于運行事件的線性關系，通過R軟件調(diào)用Tseries和Forecast程序，實現(xiàn)ARIMA模型對運行事件數(shù)量的預測。對于運行事件的非線性關系，通過R軟件調(diào)用基于Tensorflow框架的keras程序，實現(xiàn)LSTM模型的殘差預測。

采用組合模型按照5個步驟對我國核電廠運行事件數(shù)量建模并預測：

1）平穩(wěn)化處理1991～2018年運行事件數(shù)量時間序列。首先通過ADF檢測檢測未差分時間序列是否平穩(wěn)，未差分的p-value=0.22，大于0.05，可知此時間序列為非平穩(wěn)序列；經(jīng)過一階差分后p-value達到0.05，可知1階差分后時間序列平穩(wěn)化，確定數(shù)值1作為ARIMA模型參數(shù)d的值。

2）ARIMA模型參數(shù)p、q定階。通過一階差分后運行事件數(shù)量時間序列的自相關圖和偏自相關圖，如圖3所示。圖3（a）中說明，有3個自相關系數(shù)顯著不為0，參數(shù)q可以推斷取3或2（考慮3是偶然出現(xiàn)的因素，則取2是顯著不為0的階數(shù)），圖3（b）中說明，有2個偏自相關系數(shù)顯著不為0，參數(shù)p可以推斷取2或1（考慮2是偶然出現(xiàn)的因素，則取1是顯著不為0的階數(shù)）。

圖3 一階差分后的自相關圖ACF(a)和偏自相關圖PACF(b)Fig.3 ACF(a)and PACF(b)maps of first order difference

根據(jù)參數(shù)p、q初步定階結果，ARIMA模型可暫定為ARIMA（3，1，2）、ARIMA（3，1，1）、ARIMA（2，1，2）、ARIMA（2，1，1），等待赤池檢驗確定最終模型參數(shù)。由赤池信息準則（Akaike Information Criterion，AIC）可知，AIC結果大小決定著模型優(yōu)劣，數(shù)值越小表明模型越佳，模型參數(shù)估計見表1。從表1中AIC值可知，ARIMA（2，1，2）為運行事件預測的最佳模型。

表1 ARIMA模型參數(shù)估計Table 1 ARIMA model parameter estimation

3）白噪聲檢驗。通過計算Ljung-Box檢驗統(tǒng)計量，用于檢驗給定時間序列中獨立的零假設。檢驗ARIMA（2，1，2）模型殘差的白噪聲，殘差的白噪聲檢驗統(tǒng)計量結果p-value=0.998 3，大于0.05，說明模型殘差顯著性檢驗沒有通過，確定了殘差無自相關，ARIMA（2，1，2）模型能夠較好地擬合核電廠運行事件數(shù)量時間序列。

4）ARIMA模型預測。本文采用ARIMA（2，1，2）模型對1991～2018年運行事件數(shù)數(shù)量進行預測，并預測了在95%置信限值下的2019～2021年運行事件數(shù)量預測結果，如表2所示。

表2 2019～2021年運行事件數(shù)量預測Table 2 Predicted values of operation events from 2019 to 2021

5）LSTM殘差預測，合并得到組合模型預測值。采用LSTM模型對ARIMA模型預測運行事件數(shù)量殘差進行訓練和預測，訓練集和測試集比例分別占數(shù)據(jù)集的2/3和1/3，對于訓練集控制訓練迭代50次，訓練損失隨迭代次數(shù)快速下降，在20次后趨于平穩(wěn)，說明擬合良好。計算得到運行事件數(shù)量殘差的預測結果，將ARIMA運行事件數(shù)量預測結果與LSTM運行事件數(shù)量殘差預測結果相加得到組合模型預測結果，組合模型預測結果如圖4所示。ARIMA預測值相對真實值滯后的問題在時域分析是普遍存在的，如果運行事件是一個平滑序列的預測，出現(xiàn)滯后的概率就會偏低。如果在新的時刻運行事件出現(xiàn)了一個跳變被認為是客觀規(guī)律，不能判斷為異常，那么預測將會滯后。從圖4看出，2019～2021年運行事件趨勢穩(wěn)中向下，有利于減輕運行板塊的核安全監(jiān)管壓力，實現(xiàn)在不增員的前提下突出監(jiān)管重點，科學協(xié)調(diào)審評和監(jiān)督的資源向在建板塊傾斜，有利于板塊輪動和提升監(jiān)管能力。

圖4 組合模型運行事件數(shù)量預測結果Fig.4 Predicted values of combined model

4 模型評價

根據(jù)ARIMA模型和ARIMA-LSTM組合模型預測的運行事件數(shù)量結果，分別計算模型預測運行事件數(shù)量的均方誤差（Mean-Square Error，MSE）、根誤差（Root Mean Square Error，RMSE）和平均絕對百分比誤差（Mean Absolute Percentage Error，MAPE）（表3），可用于量化模型預測的準確性。

由表3可以看出，ARIMA（2，1，2）模型運行事件預測值的MSE、RMSE、MAPE分別為234.92、15.32和36.5%，ARIMA-LSTM組合模型的MSE、RMSE、MAPE分 別 為223.94、14.96和33.6%。ARIMA（2，1，2）-LSTM組合模型比ARIMA（2，1，2）模型運行事件數(shù)量的預測精度提高約3%。可以看出，ARIMA模型能夠較好地體現(xiàn)時間序列數(shù)據(jù)中的線性特征，但是單一的ARIMA模型對非線性變化難以充分有效的處理，胡衍坤等結合網(wǎng)絡結構模型能有效提高河流水質(zhì)預測精度7%，與本文運行事件數(shù)量的預測精度提高能力相近，在具有更多固定成分的不穩(wěn)定時間序列上表明LSTM做得更好，隨著未來運行數(shù)據(jù)的增加，預測精度將會更加顯著。

表3 模型評價參數(shù)Table 3 Model evaluation parameters

根據(jù)國家核安全局發(fā)布的《國家核安全局2019年報》和《國家核安全局2020年報》，可以得到2019～2020年的核電廠運行事件數(shù)量分別為31起和20起，與組合模型預測值的2019年27.64和2020年21.08非常接近，證明該組合模型能夠精確地預測核電廠運行事件數(shù)量，能夠很好地描述運行事件數(shù)量時間序列的趨勢。有利于提前調(diào)配在建板塊和運行板塊的核安全監(jiān)管審評力量，依據(jù)預測趨勢變化大小，在趨勢下降時調(diào)配資源到在建板塊，在趨勢上升時適當從在建板塊調(diào)配資源到運行板塊，前瞻性地調(diào)配有限監(jiān)管力量，實現(xiàn)板塊間合理輪動。

5 結語

本文研究利用M-K檢驗法對《核安全白皮書》中1991～2018年核電廠運行事件數(shù)量進行了趨勢和突變分析，利用ARIMA（2，1，2）模型對核電廠運行事件數(shù)量進行預測，并利用LSTM模型對非線性部分的殘差進行預測，結合兩模型預測值得到最終組合模型的預測值。通過與2019～2020年《國家核安全局年報》核電廠運行事件數(shù)量對比和驗證，其預測結果能夠準確地描述實測值，表明其預測結果具有良好的可信度。

ARIMA-LSTM組合模型比ARIMA單一模型的預測預測精度提高了3%，充分發(fā)揮了LSTM模型和ARIMA模型的優(yōu)點，前者優(yōu)化非線性部分，后者擬合線性部分，高效地處理了核電廠運行事件數(shù)量存在的非線性問題?；旌夏Ｐ蛯貌罘肿匀换貧w的ARIMA模型的殘差進行優(yōu)化，有效解決了單一模型存在的誤差問題。研究結果表明，通過對原始數(shù)據(jù)的優(yōu)化，組合模型能夠提高核電廠運行事件數(shù)量預測的精度。

作者貢獻聲明侯秦脈：提出研究思路、設計研究方案、進行數(shù)據(jù)計算、數(shù)據(jù)分析及論文的撰寫；朱偉：負責完善研究方案、數(shù)據(jù)收集、參與數(shù)據(jù)計算、稿件的審閱與修訂；鄒象：負責提供技術支持與指導；劉時賢：負責研究進度的監(jiān)督與指導；吳彥農(nóng)：負責收集參考文獻。