蔡 爽,胡瑾秋,張來斌
(中國石油大學(北京) 機械與儲運工程學院,油氣資源與探測國家重點實驗室,北京 102249)
煉油化工生產過程復雜,工藝條件苛刻,常伴有高溫、高壓、低溫、真空、大流量、高轉速等極端條件,在一個煉油廠區中,一旦某單一設備或工藝過程出現故障,易借助系統單元之間的相互依存、相互制約關系,觸發連鎖效應,由一種擾動或故障經過一段時間傳播,從一個空間傳播至另一個更廣闊的空間,引起衍生事故,帶來難以估量的社會與環境風險,引發出一系列故障鏈直至引起事故或災難。為了有效抑制煉油化工過程風險傳播,保證生產安全及產品質量,有必要對風險擾動傳播路徑進行分析研究。現有傳播路徑分析方法主要是對于故障根原因的分析,以期降低事故后果。而根原因分析的基礎在于因果圖的建立,主要分為基于過程知識的建模和基于過程數據的建模2大類。基于過程知識的因果性建模方法包括多層流模型(Multilevel Flow Model,MFM)[1]、因果模型(Cause-Effect Model,CE)[2]和符號有向圖(Signed Directed Graph,SDG)[3]等。與多層流模型與CE模型相比,符號有向圖在復雜的化工過程仿真及實際過程中得到了更為廣泛的應用[4-8]。知識建模方法對過程變量間的因果關系進行簡單和圖形化的描述,易于理解和應用;但該類方法過于依賴專家經驗,無法定量描述變量間關聯性的大小。基于過程數據的建模方法主要有互相關函數(Ross-correlation function,CCF)、傳遞熵和貝葉斯網絡模型(Bayesian Network,BN)等。互相關函數[9]通過計算2過程變量間的時滯和關聯系數,確定2變量的關聯程度及傳播方向[9-10]。但互相關函數會產生一些冗余時滯干擾分析結果且僅能分析變量間的線性相關關系。文獻[11]提出基于解析結構模型的信息融合方法,通過構建連接矩陣和可達矩陣建立相應的因果網絡。貝葉斯網絡模型[12]可以反映真實過程的隨機性和不確定性,但其概率的設置常依據主觀經驗,且概率的物理意義并不直接。傳遞熵方法是在信息熵的基礎上于2000年被提出的,這是1種基于概率分布、信息熵及統計方法得出時間序列間因果性的方法,可用于分析變量間的非線性相關關系。文獻[10,11,13]通過計算2變量間的傳遞熵,分析過程變量間由于信息傳遞所帶來的因果性。但有時由該數據驅動方法建立的擾動傳播路徑可能與實際過程不符,有必要結合過程知識檢驗傳播路徑的合理性。
因現有方法多是針對故障根原因的診斷,忽視了對風險傳播過程的預測及預防,本文提出一種基于傳遞熵與核極限學習機(Kernel Extreme Learning Machine,KELM)的煉油化工過程風險傳播路徑分析方法,針對某一工藝擾動分析風險傳播過程,基于傳遞熵分析不同過程變量間的非線性關聯關系及傳播方向,建立煉油化工過程風險傳播推繹模型,并提出一種基于KELM的風險傳播搜索方法,預測未來一段時間內風險的可能傳播路徑,以便操作人員及時切斷風險傳播路徑,保證煉油化工過程生產安全及產品質量。
傳遞熵方法是一種基于概率分布、信息熵及統計方法得出時間序列間因果性的方法,可用于分析煉油化工過程監控變量間的非線性相關關系及兩變量間信息傳遞的方向。核極限學習機方法可用于實現煉油化工過程監控變量的時間序列預測,是一種具有較強泛化能力及穩定性的預測方法。下面分別對2種方法的基本原理作一簡單介紹。
信息熵是信息論中用于度量信息量的一個概念,其定義為
(1)
式中:p(x)表示變量x的概率分布;Hx表示變量x的信息熵。
2個變量x和y的信息熵大小可用聯合熵Hxy表示,其定義為
(2)
式中:p(x,y)為x和y的聯合概率分布。
互信息[14]是信息論中用于表示信息之間相關性的一個概念,反映一個隨機變量中包含的關于另一個隨機變量的信息量,2個變量x和y的互信息Ixy如式(3):

(3)
式中:p(x)和p(y)分別表示變量x和y的概率分布,p(x,y)為x和y的聯合概率分布。
但互信息僅能表示兩變量間關聯性的大小,而無法體現兩變量間信息傳遞的方向。為準確度量動態過程中隨機變量之間的關聯性及傳播方向,T. Schreiber于2000年在信息熵的基礎上提出了傳遞熵分析法[15],這是一種基于概率分布、信息熵及統計方法得出時間序列間因果性的方法。因傳遞熵的計算需要長度較大的時間序列,對于過去普遍數據量較小的時代,其應用領域十分受限。隨著自動化水平的不斷提高,各種傳感器被大量應用于煉油化工過程,可被挖掘的數據量和信息也隨之增大,使得傳遞熵在煉油化工過程中的應用成為可能。
變量y到變量x的傳遞熵定義如式(4):
(4)

y到x的傳遞熵實質為y的信息對于x不確定性大小的改變,即y傳遞給x的信息量的大小。傳遞熵可作為衡量變量間因果性的指標。由于傳遞熵考慮的是變量間的信息量傳遞,而不需要假定變量間具有特定形式的關系,因此具有比Wiener-Granger因果性更好的適用性,尤其是對于具有非線性特征的煉油化工過程變量。
極限學習機(Extreme Learning Machine,ELM)[16]是一種單隱層前饋神經網絡學習算法。傳統的神經網絡學習算法(如BP算法)需要人為設置大量的網絡訓練參數且容易產生局部最優解,而極限學習機只需設置網絡的隱層節點個數,在算法執行過程中不需要調整網絡的輸入權值及隱層神經元的偏置,并且產生唯一的最優解,因此具有學習速度快且泛化性能好的優點。
極限學習機網絡訓練模型如圖1所示:

圖1 極限學習機網絡模型Fig.1 ELM model

(5)
式中:ωi=[ωi1,ωi2,,ωiz]T為第i個隱層節點的輸入權重;βi為第i個隱層節點的輸出權重;bi為第i個隱層節點的偏置;oj為ELM網絡模型的實際輸出值。
(6)
上式可由矩陣表示為
Hβ=U
(7)
β=[β1,β2,…,βL]T,U=[u1,u2,…,uW]T
式中:H為ELM的隱層輸出矩陣;β為輸出權重矩陣;U為期望輸出向量。
可采用最小二乘法求取最優權重向量β*使得實際輸出與期望輸出差值的平方和最小,其解為
β*=HΨU
(8)
式中:HΨ是矩陣H的Moore-Penrose廣義逆。

Ω=HHT,
Ωij=h(xi)·h(xj)=K(xi,xj),i,j=1,2,…,W
(9)
核矩陣Ω替代ELM中的隨機矩陣HHT,利用核函數將z維輸入樣本映射到高維隱層特征空間。核函數K(xi,xj)是核矩陣Ω中第i行第j列的元素,包括RBF核函數、線性核函數和多項式核函數等,通常選擇RBF核函數,其表達式為:
K(xi,xj)=exp(-γ(xi-xj)2),γ>0
(10)
式中:γ為核參數。
將參數I/C添加到HHT中的主對角線上,使其特征根不為零,并由此求權重向量β*為
β*=HT(I/C+HHT)-1U
(11)
式中:I為單位對角矩陣;C為懲罰系數,用于權衡結構風險和經驗風險之間的比例。
可得KELM模型的實際輸出為
(12)
式中:α=(I/C+HHT)-1U,為KELM網絡的輸出權值。
當過程出現異常工況產生擾動時,報警系統將通過聲光形式的報警信號向操作者傳遞異常信息。報警通常由連續過程中存在的擾動引起。由于過程設備之間的關聯性,擾動可能傳播并影響大量過程變量,如果不能及時得到有效控制,可能導致風險的出現甚至災難性的后果。因此,有必要分析過程變量間的關聯關系,當報警發生時及時辨識并有效遏制風險傳播途徑。本文提出基于傳遞熵與KELM的煉油化工過程風險傳播路徑分析方法,如圖2所示。針對某一工藝擾動分析風險傳播過程,建立煉油化工過程風險傳播推繹模型,提出一種風險傳播搜索方法,預測風險傳播路徑,以便操作人員及時采取預防措施,保證煉油化工過程生產安全及產品質量。
煉油化工過程工藝復雜,設備眾多,同一過程設備及相鄰設備間的監控變量往往具有關聯性。采用傳遞熵分析方法可分析不同風險過程變量間的關聯關系,推斷他們固有的因果關系,從而建立煉油化工過程風險傳播推繹模型,方法如下:
對于過程監控變量X和Y,變量Y到X的傳遞熵TY→X的計算公式如式(4),Y到X的傳遞熵實質為Y的信息對于X不確定性大小的改變,即Y傳遞給X的信息量的大小。因此,傳遞熵可以作為衡量變量間因果性的指標。

圖2 基于傳遞熵與KELM的煉油化工過程風險傳播路徑分析方法流程Fig.2 Flowchart of the transfer entropy and KELM based risk propagation analysis method
同樣,可求出變量X到Y的傳遞熵,2變量間的關聯系數由式(13)確定:
(13)
如果ρX,Y=TY→X,表示2變量間傳播方向為Y→X;如果ρX,Y=TX→Y,表示2變量間傳播方向為X→Y;如果ρX,Y=0,表示2變量無因果關系。
因關聯系數由統計方法計算得到,每兩個時間序列可得到一確定值,但若TY→X與TX→Y的差值過小,考慮2變量間的因果性將沒有意義。因此有必要設置合適的閾值對2變量間因果關系的顯著性水平進行檢驗。本文采用如下的假設檢驗方法。
選取73個長度為200的隨機生成序列,令tX→Y=TX→Y-TY→X,通過式(4)計算每兩個序列間的tX→Y值,并求出所有|tX→Y|(tX→Y的絕對值)的均值μtX→Y和標準差σtX→Y,通過式(14)計算閾值[17],判斷2變量間因果關系的顯著性。若關聯系數沒有通過顯著性檢驗,表明2變量間不具備顯著的因果性。
|tX→Y|-μtX→Y≥6σtX→Y
(14)
通過計算可得μtX→Y=0.053 5,σtX→Y=0.353。
對于N個風險過程變量X1,X2,…,XN,可通過計算傳遞熵確定其中每兩變量的關聯系數及傳播方向,并據此建立風險傳播推繹模型,如圖3。模型由有向弧和代表風險過程變量的節點組成。對于任意2過程變量Xi和Xj,若tXi→Xj>0,有向弧由Xi指向Xj,即由上級原因變量指向下級影響變量;反之,則傳播方向相反。

圖3 風險傳播推繹模型示意Fig. 3 Diagram of risk propagation reasoning model
最后結合過程知識檢驗傳播路徑的合理性,對模型進行適當修正。
過程變量發生報警通常由連續過程中存在的風險擾動引起。由于過程設備之間的關聯性,擾動可能傳播并影響大量過程變量,若不能及時得到有效控制,可能導致風險的出現甚至災難性的后果。為此,提出一種風險傳播搜索方法,以便當報警發生時及時辨識風險傳播途徑,提醒操作人員及時采取有效措施,避免異常風險的進一步發展。
風險傳播搜索方法主要包括如下4個步驟:
步驟1:相同的報警可能產生不同的風險傳播路徑。為了準確辨識風險傳播路徑,當某一設備的過程變量Xj發生報警時,將其作為上級原因變量,根據2.1節所提方法建立的風險因素傳播模型搜索與其直接相連的同設備及其相鄰設備中所有下級影響變量節點,例如,圖3中的X4報警,可搜索到其2個下級變量節點X3和X5。
步驟2:若變量Xj在tκ時刻發生擾動,對于Xj的各相關下級變量Xi(i=1,2,…,I;I為下級變量個數),可通過式(15)計算變量Xi的擾動變化率,其值大小可近似反映各相關下級變量受上級變量擾動的影響大小。
擾動變化率(Disturbance Rate, DR)定義如下:
定義1:對于變量Xi,考慮以時刻tκ為中心,選擇時間間隔為[tκ-m,tκ+m] (時間序列長度為2m+1)的變量Xi的時間序列進行最小二乘線性擬合如式(15),所求斜率ai絕對值的大小作為變量Xi的擾動變化率。
其中,變量Xi在tκ+1至tκ+m時刻的值通過KELM方法預測得到,以變量Xi在tκ時刻前一段時間內(ts至tκ時刻)的歷史數據構造KELM的訓練樣本,輸入樣本Xi*和輸出樣本U分別為

z為輸入樣本維數,樣本個數為W=κ-z-s+1,xil(l=s,s+1,…,κ)為變量Xi在第tl時刻的值。
令xik=aitk+bi,其中,k=κ-m,…,κ,…,κ+m;tk=1,2,…,2m+1;i=1,2…,N;N為過程變量個數;xik為標準化后的變量Xi在第tk個時刻的變量值,最小二乘線性擬合公式如式(15):
(15)
變量Xi的標準化公式如式(16):
(16)

步驟3:通過式(13)計算兩過程變量Xi和Xj間的關聯系數,通過式(15)計算變量Xi的擾動變化率,基于所求關聯系數和擾動變化率,根據式(17)計算上級變量Xj對各下級變量Xi的影響因數Ri,比較影響因數Ri的大小,將影響因數最大值對應的下級變量作為其下級影響變量。
下級變量的影響因數Ri越大,其受上級變量的影響越大,因此,通過Ri的大小可以比較上級變量對各相關下級變量產生擾動的影響大小,將Ri的最大值對應的下級變量作為其下級影響變量。若所求影響因數過小,說明該變量的時間序列數據趨于平穩,并受到上級變量影響,因此,根據專家經驗及歷史數據統計,如果Ri值小于閾值Rth(這里依據專家經驗與歷史統計,設為1.75),將不考慮該下級變量。
(17)

步驟4:重復步驟2和步驟3,依次確定可能受到風險影響的各設備過程變量并確定最可能的風險傳播路徑。
催化裂化裝置是煉油化工過程的關鍵裝置。隨著煉油廠生產規模的不斷擴大,反應再生系統、分餾系統結焦、分餾塔沖塔、油漿泵故障等常導致催化裂化裝置發生故障、非計劃停機甚至災害性事故。分餾單元是催化裂化生產裝置的一部分,由反應器來的反應產物(油氣)從底部進入分餾塔,油氣經分餾后得到的富氣、粗汽油、輕柴油、重柴油、回煉油及油漿。沖塔故障是分餾過程中經常發生的現象之一。當氣液相負荷過大時,氣體通過塔板的壓降增大,會使降液管中液面高度增加,液相負荷增加時,出口堰上液面高度增加,當液體充滿整個降液管時,上下塔板液體連成一片,分餾完全破壞,導致沖塔的發生。沖塔故障一旦發生,產品質量將受到嚴重影響。造成分餾塔沖塔的原因有很多,如塔盤掀翻或損壞、塔盤結鹽、換熱器故障、各種機泵故障等。為了保證產品質量,避免沖塔故障的發生,需分析沖塔風險,以下以某石化企業催化裂化分餾單元的沖塔風險過程為例,進行風險傳播推繹分析。
就分餾塔沖塔風險而言,主要對分餾單元進行分析建模。催化裂化生產裝置中的分餾單元主要包括分餾塔、分餾塔頂油氣分離器、回煉油罐、原料油緩沖罐和輕柴油汽提塔。某石化企業催化裂化分餾單元的沖塔風險過程變量如表1所示。

表1 分餾單元沖塔風險過程變量
根據沖塔風險過程變量及某石化企業的歷史數據(采樣間隔為5 s,數據長度為1 000),通過傳遞熵計算(見1.2節式(13))可得到分餾單元的沖塔風險過程變量間的關聯系數,并依據式(14)判斷2變量間因果關系的顯著性;最后結合過程知識修正,建立沖塔過程風險傳播推繹模型如圖4。圖4中各變量之間的關聯系數如表2。

圖4 沖塔過程風險傳播推繹模型Fig.4 Risk propagation reasoning model in the flooding process
以一個發生在某石化公司的分餾塔沖塔事件為例,采用所提出的風險傳播搜索方法進行分析。2016年8月23日,某石化公司催化裂化分餾塔發生沖塔,經現場人員調查分析,沖塔發生原因為一中回流泵故障,引起一中回流量過低發生低報警,使得分餾塔塔頂溫度升高,重柴油側線餾出口溫度升高,分餾塔底液位過低,造成回煉油下返塔流量過低,操作人員發現回煉油下返塔流量過低后及時排查出了回流泵的故障原因,并采取措施使過程恢復到了正常運行狀態。

表2 各過程變量間關聯系數
基于前文所建立的沖塔過程風險傳播推繹模型,進行風險傳播路徑搜索分析,如圖4所示:
步驟1:因一中回流量過低(節點9)發生報警,故將其作為上級變量,根據因果圖搜索與其直接相連的下級變量分餾塔塔頂溫度(節點16)。
步驟2:通過極限學習機預測方法預測節點9發生報警時刻后3分鐘(共36個時刻)節點16的變量值,由式(15)可得其擾動變化率如表3所示,各變量時間序列的標準化后數據擬合曲線如圖5,各變量單位見表1。

表3 各上級變量的擾動變化率及影響因數

圖5 各下級變量時間序列數據擬合曲線Fig. 5 Time sequence fitting curve of each subsequent variable
步驟3:根據式(17)計算節點16的影響因數如表3所示,因節點16的影響因數值Ri=1.905 3,大于閾值Rth(這里依據專家經驗及歷史統計設為1.75),故將分餾塔塔頂溫度(節點16)作為下級影響變量。
步驟4:繼續搜索節點16的下級變量包括重柴油側線餾出口溫度(節點5)、分餾塔塔底液位(節點7)、輕柴油側線餾出口溫度(節點17)、塔頂油氣分離罐液位(節點8)和分餾塔頂循回流溫度(節點2),通過極限學習機預測方法預測節點9發生報警時刻后3 min(共36個時刻)各節點的變量值,由式(15)及式(17)可得其擾動變化率及影響因數如表3所示,通過比較各相關下級變量的影響因數,其最大值為1.771 4,對應節點7,大于閾值Rth,因此將該最大值對應的下級變量分餾塔塔底液位(節點7)作為下級影響變量。繼續搜索節點7的下級變量包括一中返塔溫度(節點10)和回煉油下返塔流量(節點18),由式(15)及式(17)可得其擾動變化率及影響因數如表3所示,通過比較各相關下級變量的影響因數,其最大值為2.335 7,對應節點18,大于閾值Rth,因此將該最大值對應的下級變量回煉油下返塔流量(節點18)作為下級影響變量。繼續搜索節點18的下級變量輕柴油汽提塔液位(節點1),由式(15)及式(17)可得其擾動變化率及影響因數如表3所示,因節點1的影響因數值Ri=1.654 1,小于閾值Rth,搜索結束。
1)通過傳遞熵方法分析不同尺度下沖塔風險過程變量間的非線性關聯關系及風險傳播方向,采用所提風險傳播搜索方法辨識出沖塔風險在整個分餾單元內的最可能傳播路徑為:一中回流量過低(節點9)→分餾塔塔頂溫度(節點16)→分餾塔塔底液位(節點7)→回煉油下返塔流量(節點18),該結論與前文所述的現場異常工況發生過程一致,從而驗證了所提方法的有效性。
2)根據沖塔風險過程中2個上下級變量間的關聯系數及下級變量的擾動變化率計算上級變量對下級變量的影響因數,影響因數的大小近似反映了上級變量對下級變量產生擾動的影響大小,有助于確定沖塔風險過程的傳播路徑,例如,節點7的下級變量包括一中返塔溫度(節點10)和回煉油下返塔流量(節點18),通過比較各相關下級變量的影響因數,其影響因數最大值對應的下級變量回煉油下返塔流量(節點18)即為節點7的下級影響變量。
3)采用所提方法對分餾塔沖塔過程風險傳播路徑進行分析,不僅考慮了存在沖塔風險設備分餾塔本身,還分析了分餾單元內與分餾塔存在關聯的其他設備,避免遺漏某些可能受沖塔影響的風險過程變量,從而通過KELM預測風險傳播路徑,操作人員及時采取預防措施,合理規避風險。
1)針對分餾塔沖塔風險過程,采用傳遞熵方法分析分析不同過程變量間的非線性關聯關系及傳播方向,建立了沖塔風險過程風險傳播推繹模型,并基于KELM預測,提出一種風險傳播搜索方法,可預測未來一段時間內風險的可能傳播路徑。
2)所提方法可辨識風險擾動傳播給某些設備甚至整個生產單元帶來的風險影響,有助于操作人員準確把握風險動態并及時采取合理措施抑制風險的進一步發展。
[1]DAHLSTRAND F. Consequence analysis theory for alarm analysis[J]. Knowledge-Based Systems, 2002, 15(1): 27-36.
[2]WAN Yiming, YANG Fan, LYU Ning, et al. Statistical root cause analysis of novel faults based on digraph models[J]. Chemical Engineering Research & Design, 2013, 91(1): 87-99.
[3]YANG Fan, SHAH S L, XIAO De-yun. Signed directed graph modeling of industrial processes and their validation by Data-Based methods[C]//2010 CONFERENCE ON CONTROL AND FAULT-TOLERANT SYSTEMS (SYSTOL'10).Nice,France, 2010: 387-392.
[4]MAURYA M R, RENGASWAMY R, VENKATASUBRAMANIAN V. A systematic framework for the development and analysis of signed digraphs for chemical processes. 1. Algorithms and analysis[J]. Industrial & Engineering Chemistry Research, 2003, 42(20): 4789-4810.
[5]MAURYA M R, RENGASWAMY R, VENKATASUBRAMANIAN V. A systematic framework for the development and analysis of signed digraphs for chemical processes.2.Control loops and flowsheet analysis[J]. Industrial & Engineering Chemistry Research, 2003, 42(20): 4811-4827.
[6]HE Bo, CHEN Tao, YANG Xianhui. Root cause analysis in multivariate statistical process monitoring: Integrating reconstruction-based multivariate contribution analysis with fuzzy-signed directed graphs[J]. Computers & Chemical Engineering, 2014, 64: 167-177.
[7]MAURYA M R, RENGASWAMY R, VENKATASUBRAMANIAN V. A signed directed graph-based systematic framework for steady-state malfunction diagnosis inside control loops[J]. Chemical Engineering Science, 2006, 61(6): 1790-1810.
[8]CHANG C T, CHEN C Y. Fault diagnosis with automata generated languages[J]. Computers & Chemical Engineering, 2011, 35(2): 329-341.
[9]BAUER M, THORNHILL N F. A practical method for identifying the propagation path of plant-wide disturbances[J]. Journal of Process Control, 2008, 18(7): 707-719.
[10]HAN Liu, GAO Huihui, XU Yuan, et al. Combining FAP, MAP and correlation analysis for multivariate alarm thresholds optimization in industrial process[J]. Journal of Loss Prevention in the Process Industries, 2016, 40: 471-478.
[11]GAO Huihui, XU Yuan, GU Xiangbai, et al. Systematic rationalization approach for multivariate correlated alarms based on interpretive structural modeling and Likert scale[J]. Chinese Journal of Chemical Engineering, 2015, 23(12): 1987-1996.
[12]ABELE L, ANIC M, GUTMANN T, et al. Combining Knowledge Mo ̄deling and Machine Learning for Alarm Root Cause Analysis[J]. IFAC Proceedings Volumes, 2013, 46(9):1843-1848.
[13]GAO J, TULSYAN A, YANG F, et al. A transfer entropy method to quantify causality in stochastic nonlinear systems[J]. IFAC-PapersOnLine, 2016, 49(7): 454-459.
[14]ZHAO Xiaojun, SHANG Pengjian, LIN Aijing. Transfer mutual information: A new method for measuring information transfer to the interactions of time series[J]. Physica a: Statistical Mechanics and Its Applications, 2017, 467: 517-526.
[15]SCHREIBER T. Measuring information transfer[J]. Physical Review Letters, 2000, 85(2): 461-464.
[16]HUANG Guangbin, ZHU Qinyu, SIEW C K. Extreme learning machine: Theory and applications[J]. Neurocomputing, 2006, 70(1/3): 489-501.
[17] HAJIHOSSEINI P, SALAHSHOOR K, MOSHIRI B. Process fault isolation based on transfer entropy algorithm[J]. ISA Transactions, 2014, 53(2): 230-240.