一種面向實時風險預測的工業報警優先級評估方法*

王 佳，陳 勛

（1.中國軟件評測中心，北京 100048；2.中國電子科技網絡信息安全有限公司，四川 成都 610000）

0 引 言

報警優先級評估是報警管理系統的一項重要內容，主要用來表達不同過程變量超越報警閾值的危害性，能夠幫助工廠操作人員有針對性地處理實時產生的報警，采取更加有效的措施[1-6]。報警變量優先級的確定應該同時考慮對工廠和系統的實時危害程度和響應的緊急程度[7]。

風險危害程度一般有定性和定量兩種分析方法[8]。定性的分析方法有危害和可操作分析（Hazards and operability Analysis）、過程危害分析（Process Hazards Analysis）、傳統的失效模式與影響分析（Failure Mode and Effects Analysis）以及復合有向圖等[9-10]；定量方法包括建立故障樹分析模型[11]、馬爾可夫模型以及因果關系分析圖等。此外，隨著安全分析理論的進一步發展，涌現出貝葉斯（Bayesia）網絡、復雜網絡模糊映射等與多學科相結合的分析方法。上述方法能夠對變量的風險進行描述和建模，但是主要依靠專家經驗，并且由于變量不確定性，不能從實時狀態變化的角度給出變量的風險轉移過程。Chang[12]等通過建立危險發生概率、潛在危險影響和過程安全時間的綜合模型，對整個過程的安全報警進行優先級評估來減少報警的數量，但是需要依靠專家經驗，不能有效運行過程中檢測的數據，“實時性”明顯不足。Foong[13]等結合整個過程中的數據和專家知識，利用模糊推理對產生的報警變量進行優先級排序，但是欠缺對報警發生時間的考慮。由于隨機不確定性因素的影響，各個時刻之間的狀態轉移可用Markov過程描述，但是轉移概率不能直接測量獲取具有時變的特征，每一時刻的轉移概率值不是固定的。目前，粒子濾波對非齊次轉移概率的估計不受過程噪聲和測量噪聲限制。相比Kalman濾波，主要利用一系列離散采樣值近似實現轉移概率的最大后驗估計，得到狀態的估計值[14-15]。

因此，考慮到變量的不確定性因素，本文針對報警風險狀態轉移過程中的狀態和時間兩個方面進行分析，建立幅值裕度和時間裕度定量性分析指標。由于受到外界的干擾轉移概率的隨機時變性，建立動態Markov狀態轉移模型描述報警狀態轉移和時間變化過程，通過粒子濾波估計狀態轉移矩陣，得到系統狀態變化的估計值，由此提出了一個新的報警優先級評價的系統優化方法。

1 報警風險評價指標

1.1 幅值裕度指標

一般情況下，工業過程報警狀態是從低報警到高報警狀態進行演化的，因此在對報警危害的嚴重程度進行評價時，可以選用計算簡單、反映經濟指標或者物理變化過程的幅值裕度作為評價指標。

實時報警系統主要監視過程變量的參數值，一旦超過給定的報警閾值，立即觸發報警。由于工藝的限制，不同工況下工藝參數的報警閾值不一樣，可以設置成定值或者是與工況相聯系的預先設定的變量。圖1為報警設定值的參數關系，其中符號A代表在參數正常運行過程中擾動值和報警整定值之間的差值，B代表為了克服擾動帶來的干擾，擾動值與采取保護動作定值之間的差值。

圖1 報警整定值的設定

由于生產過程中報警的嚴重程度與產生的后果、變量的初始狀態和工況條件等因素有關，因此報警系統應該充分考慮每一個報警變量的嚴重程度來對報警進行優先級劃分。它的危害度函數的定義如下：

其中：Sev(Tc)表示每個報警的嚴重度，Tc為實時狀態值，ρ為報警狀態轉移概率，Tn和Ts分別為報警的警戒閾值線和報警閾值線。幅值裕度指標不僅用概率的動態變化對報警進行分析，還考慮了狀態的實時信息，指出了狀態變量與最嚴重程度的偏離程度，對不同變量所處同一概率狀態下的危害程度。

1.2 時間裕度指標

對于報警系統來說，超過報警之后有高報、高高報的情況下，為了更加明確報警狀態的轉移情況，在原有狀態下增加幾個安全報警狀態，在變量的容許范圍定義報警狀態轉移時間，分別為：正常狀態（0）到低報警狀態（1）的轉移時T0→1=T1-T0；低報警狀態（1）到中報警狀態（2）之間的轉移時T1→2=T2-T1；中報警狀態（2）到高報警狀態（3）之間的轉移時 T2→3=T3-T2[16]。

根據定義的報警轉移時間裕度，時間風險指標為：

圖2為風險指標曲線，報警的風險值與時間有直接關系。報警狀態轉移時間越小，處理報警的時間越短，產生的后果也就越嚴重；報警狀態轉移時間越大，操作者處理報警的過程時間充裕，產生的風險也越小。

圖2 風險指標曲線

2 基于Markov時變模型報警狀態轉移概率預測

2.1 建立馬爾科夫時變模型

根據變量的歷史數據和運行條件預測狀態的發展趨勢，結合幅值裕度指標中的報警狀態轉移概率，了解變量未來的某一時刻運行狀態。

將Markov模型中，報警的狀態空間分成（0，1，2，3）四種元素狀態，其報警狀態轉移過程如圖3所示。通過分析報警狀態的轉移過程的狀態，進一步分析報警變量的危害程度，降低誤報率和漏報率。報警狀態轉移的集合分別為：正常狀態（0）的狀態集合{0|x1,x2,…,xn0}到低報警狀態（1）{1|xn0+1,xn0+2,…,xn1}、低報警狀態（1）{1|xn0+1,xn0+2,…,xn1}到中報警狀態（2）{2|xn1+1,xn1+2,…,xn2}、中報警狀態（2）{2|xn1+1,xn1+2,…,xn2}到 高 報 警 狀 態（3）{3|xn2+1,xn2+2,…,xn3}。由Markov過程確定的狀態轉移概率模型包括過程狀態、轉移概率和初始概率分布情況[17-18]。狀態轉移概率為：

表示時刻t處于狀態i，經過時刻t+1處于狀態j 的概率，初始概率用 π=[π1,π2,…,πn]來表示。

圖3 報警狀態網絡空間轉移過程

報警狀態的轉移規律可采用統計近似方法：

式中nij為從狀態i轉移到狀態j的樣本數，因此通過計算得到時刻t={0,1,2,3,…,m}，正常狀態（0）、低報警狀態（1）、中報警狀態（2）和高報警狀態（3）的狀態概率轉移矩陣pij(0:m)為：

根據馬爾科夫無后性和貝葉斯條件概率公式，有：

式（6）狀態轉移矩陣是不隨時間變化的，而在實際工業過程中，狀態轉移矩陣是時刻發生變化的，因此需要根據系統實時變化的狀態更新狀態轉移概率矩陣，從而更準確地估計風險狀態：

式中，第k個時刻的狀態為π(k)，由矩陣中最大的概率值確定當前的報警狀態，其中狀態轉移矩陣P′的準確程度是最關鍵的因素。

2.2 粒子濾波估計動態狀態轉移矩陣

狀態轉移矩陣不能直接測量，通常是實驗得到，因此存在很多不確定性。粒子濾波通過系統狀態條件分布產生一組隨機樣本，然后不斷調整修正粒子最初的條件分布。

假設P(xk(0))作為先驗知識，那么P(xk(t)|yk(1:t))通過預測方程和概率的實時更新遞推得到，然后利用最新得到的值對后驗概率密度進行估計修正。

預測方程為：

更新方程為：

P(xk(t)|xk(1:t-1))由狀態轉移矩陣獲得，P(yk(t)|xk(t))由預測方程獲取。通過粒子濾波對預測方程和更新方程進行計算，計算得到狀態轉移過程的狀態轉移序列。它的狀態轉移序列的概率越大，代表下一時刻將要發生的情況。

假設變量的采樣時間間隔為Δt，對于報警變量從時刻0運行到時刻m的時間轉移過程的報警轉移時間的計算可分別得到如下內容：

（1）r(tm)=0且r(tm+1)=1，系統從正常狀態轉移到低風險狀態的報警轉移時間為T0→1=m×Δt；

（2）r(tm)=1且r(tm+1)=2，系統從低風險狀態轉移到中風險狀態的報警狀態轉移時間為T1→2=m×Δt-T0→1；

（3）r(tm)=2且r(tm+1)=3，系統從中風險狀態轉移到高風險狀態的報警狀態轉移時間為T2→3=m×Δt-T1→2。

因此，通過已知的初始時刻和初始報警狀態，利用式（7）得到第k時刻處于哪種報警狀態空間（0，1，2，3）的概率情況，且得到處于不同報警狀態空間的統計時間，綜合發生的概率和時間對報警等級進行排序。

3 報警優先級評估計算

依據文中的報警幅值裕度指標和時間裕度指標，綜合分析對報警優先級進行劃分。流程如圖4所示，分為離線訓練和在線監測評價兩個部分。

圖4 報警優先級評估流程

離線部分：

（1）設變量的采樣時間間隔為Δt，初始時刻為k=0，m=0；Δt=0所有的訓練數據均正常，未發生報警事件；

（2）將歷史數據作為訓練數據，得到初始的狀態轉移矩陣；

（3）根據馬爾科夫模型（4），得到各個狀態發生的概率，并且估計報警狀態的轉移時間，由此得到報警狀態的轉移時間；

（4）根據式（1）和式（2），得到訓練數據的報警等級指標；

（5）根據報警安全指標，得到報警等級的分類條件。

在線評價：

（1）當新樣本到來時，計算新樣本屬于的報警風險狀態，結合上一個報警狀態，求取當前采樣點的報警等級指標；

（2）判斷所處的報警等級，計算當前的幅值和時間裕度，根據得到的裕度，預測下一時刻運行狀態報警等級。

4 算例分析

以化工模型田納西-伊斯曼（Tennessee Eastman，TE）過程為例。整個過程有6種操作模式和反應器、冷凝器、汽提塔、氣液分離器及壓縮機5個組成單元。化學反應方程式如下：

A(g)+C(g)+D(g)→G(liq)

A(g)+C(g)+E(g)→H(liq)

A(g)+E(g)→F(liq)

3D(g)+→2F(liq)

其中，A、C、D、E為4種所需的反應氣體。當惰性氣體B分別加入TE反應器中進行反應時，產物分別為液體G、H和副產品F。

選取數據庫中5 500個采樣數據進行分類劃分，得到正常狀態（0）、低報警狀態（1）、中報警狀態（2）和高報警狀態（3），統計采集的各個時刻的狀態數據轉移次數，得到反應器進料流量結果如表2所示。

表2 反應器流量狀態轉移次數統計

根據統計結果計算狀態的轉移概率，得到的初始報警轉移概率矩陣如下：

假設起始時刻都處在正常狀態，初始狀態的概率為π(0)={1,0,0,0}，將π(1)=π(0)P作為當前的狀態概率。表3給出了傳統的Markov狀態轉移變化結果。

表3 傳統的Markov預測結果

由于工業現場環境非常復雜，外界干擾元素影響多，具有不確定性，其每一個變量的狀態轉移概率很難保持不變，是實時發生變化的，因此固定的狀態轉移概率不能反映真實情況，可以將轉移矩陣更改為：

其中，φ=rand(0,1)為干擾變量，代表在區間[0，1]中產生的一系列隨機數。為了更加真實地反映變量的狀態情況，利用粒子濾波對變量的狀態轉移概率矩陣進行估計預測，可以得到動態的狀態轉移矩陣。圖5和圖6為采用3種濾波估計的仿真曲線圖，通過對比不敏感卡爾曼濾波（Unscented Kalman Filter，UKF）、粒子濾波（Particle Filter，PF）和擴展卡爾曼濾波（Extended Kalmar Filter，EKF）3種方式對真實數值的估計曲線和誤差曲線，得出粒子濾波具有更小的誤差。

圖5 報警狀態轉移矩陣估計

由于傳統馬爾可夫預測未能真實反映狀態信息，導致預測結果和實際值有較大的誤差。但是，通過粒子濾波估計狀態轉移矩陣的狀態概率，使不確定性的轉移概率具有較強的魯棒性。其中，利用粒子濾波估計狀態轉移矩陣的狀態概率如表4所示。

圖6 估計誤差

根據馬爾科夫狀態轉移矩陣各個報警狀態的概率值大小來確定危害程度，對得到的報警進行分類，確定當前時刻實際的報警狀態是處在正常狀態（0）、低報警狀態（1）、中報警狀態（2）和高報警狀態（3）中的哪一類。

依據報警風險估計的預測結果，得到系統的報警風險狀態轉移時間分別為T0→1=m×Δt、T1→2=m×Δt-T0→1和T2→3=m×Δt-T1→2。對比各個報警狀態轉移時間的大小，可以得出處理變量的緊急程度。

利用上述方法，選擇化工TE仿真模型中模式一中的反應器進料量進行評估，得到不同時刻的報警分類圖，如圖7所示。

表4 Markov時變模型預測結果

圖7 反應器進料量報警狀態

同理，圖8、圖9得到TE過程的反應器壓力XMEAS8和反應器液位XMEAS9得到的變量報警圖。如果變量所處的報警狀態危害程度發生的概率相同，通過比較報警狀態轉移時間的大小來確定報警狀態的處理的優先級。t=20時刻下報警優先級對比報警狀態和轉移時間得到同一時刻的報警排序結果，如表5所示。

圖8 反應器壓力報警狀態

圖9 反應器液位報警狀態

表5 t=20報警狀態等級劃分

5 結 語

報警分級評價模型能夠幫助工廠操作者引起對重要報警的注意，而不需要在誤報警和等級低的報警上花費更多時間。本文定義了量化報警風險的幅值裕度和時間裕度指標，考慮報警變量的危害程度和報警的轉移時間建立報警優先級評估模型，綜合評價處于不同狀態的報警，以降低同時產生多個報警導致的重大損失。同時，考慮實際工況情況，降低外界不確定因素對變量的影響程度，利用粒子濾波估計Markov狀態轉移過程，預測報警狀態變量的狀態轉移概率，實時預測反應狀態變量的變化。通過化工TE過程算例的計算結果表明，提出的風險量化指標對報警狀態的排序和對報警狀態安全風險評估有一定的指導意義。由于報警系統的處理主要依靠專家經驗，具有一定的主觀性，在下一步的工作中，在實時報警等級管理系統中采用混合技術，通過自學習、自診斷修改模型來提高報警分類和預警的準確性。