胡瑾秋， 張立強， 張來斌

(中國石油大學(北京) 油氣資源與探測國家重點實驗室 機械與儲運工程學院， 北京 102249)

安全預警技術早已被廣泛應用到各石油化工企業(yè)之中。安全預警通過分布式控制系統(tǒng)(DCS)采集的數(shù)據(jù)，對關鍵參數(shù)進行多方位的綜合分析，實時捕捉裝置運行過程中存在的異?，F(xiàn)象[1]。但這類報警方法會因噪聲及采集誤差導致大量錯誤報警，且不能提前發(fā)出警告[2]?，F(xiàn)有的預警方法主要包括主元分析的方法和參數(shù)預測的方法，齊詠生等[3-5]基于傳統(tǒng)的主元分析方法，改進了主元選取方法及統(tǒng)計量的計算方法，實現(xiàn)故障的預測。但主元分析方法面對數(shù)據(jù)變化不太明顯的異常時，會有很高的漏報率。Wang等[6]運用經(jīng)驗模態(tài)分解(EMD)方法對非穩(wěn)態(tài)的時間序列進行分解，并進行多步前向預測，但這類方法仍屬于短期預測。Chang等[7]提出了化工過程裝置的風險預警系統(tǒng)，運用定性方法對裝置的風險值進行打分，但該方法無法準確表達異常工況產(chǎn)生的風險值。

針對以上問題，筆者開展了針對化工過程預警技術的研究。提出了運用故障鏈預測煉化裝置故障發(fā)生后狀態(tài)；應用灰色關聯(lián)分析法和卡爾曼濾波法實現(xiàn)多參數(shù)關聯(lián)預警；以及運用粒子群優(yōu)化最小二乘向量機預測過程參數(shù)實現(xiàn)預警等有效的預警方法[8-10]。

在此基礎上，筆者針對石油化工裝置長周期運行風險超早期精確預警方法開展了相關研究，提出基于損失函數(shù)的石油化工裝置風險預警方法。損失函數(shù)通常用于精確描述目標偏移與預期損失之間的關系，運用損失函數(shù)，計算參數(shù)偏差可能造成的預期損失，并運用剩余時間模型計算異常工況發(fā)生的可能概率，從而得到該參數(shù)偏差所產(chǎn)生的風險，通過監(jiān)測風險值的變化實現(xiàn)裝置的風險預警。經(jīng)過案例分析，相比DCS系統(tǒng)，該方法可提前10 min左右監(jiān)測到異常工況，實現(xiàn)了異常工況的有效預警，為控制石油化工生產(chǎn)過程中異常工況的產(chǎn)生提供了指導。

1 風險預警方法的基本原理

1.1 損失函數(shù)構造

損失函數(shù)根據(jù)其應用范圍不同分為許多種類，常見的如平方損失函數(shù)、加權損失函數(shù)、轉置正態(tài)損失函數(shù)等。Hashemi等[11-12]總結了不同損失函數(shù)的特性及對損失的刻畫能力?？紤]石油化工裝置參數(shù)與實際值的隨機偏離，筆者采用修正型轉置正態(tài)損失函數(shù)。該損失函數(shù)可以通過調節(jié)形狀因子修正損失函數(shù)的形狀，從而更好地刻畫實際損失。損失函數(shù)L如式(1)所示：

(1)

式(1)中，y為狀態(tài)參數(shù)值；Δ為能夠造成最大損失的參數(shù)值到參數(shù)目標值的距離；γ為形狀因子；T為狀態(tài)參數(shù)目標值。損失函數(shù)的變化范圍0≤L≤1。

根據(jù)實際情況的不同，可以得到不同的形狀因子，從而構造出與待測參數(shù)相匹配的損失函數(shù)，進而準確描述參數(shù)偏差帶來的預計損失。由式(1)可以看出，當參數(shù)處于目標值時，損失值為0，參數(shù)偏離目標值越遠，損失值越大。

1.2 異常工況發(fā)生概率計算

運用剩余時間理論對異常工況發(fā)生概率進行計算[13]，該理論通過計算狀態(tài)參數(shù)到達報警閾值所需要的時間，對裝置的安全性進行量化。剩余時間越多，則能夠對異常狀態(tài)采取的措施越充分，裝置的安全性越高；反之，剩余時間越少，則更難以對異常狀態(tài)采取完善措施，異常工況的發(fā)生概率就會越高。

剩余時間t的表達式如下所示：

(2)

式(2)中，ylim為狀態(tài)參數(shù)閾值；ΔV為狀態(tài)參數(shù)的瞬時變化率，該參數(shù)可用最小二乘法對數(shù)據(jù)曲線進行擬合，并對擬合多項式進行一階微分計算求得[14]，最小二乘法的計算步驟如下。

如式(3)所示，x為狀態(tài)參數(shù)對應的時間，S;a為系數(shù)。通過最小二乘法確定最優(yōu)ai，使得在各個點的偏差δ平方和最小，從而得到精確的擬合多項式。

f(x)=anxn+an-1xn-1+…+a1x+a0

(3)

將k對數(shù)據(jù)(xi，yi)帶入式(3)中，得到如下方程：

(4)

將式(4)轉換成矩陣的形式，求得矩陣的唯一1組最優(yōu)近似解，使得偏差δ的平方和最小，從而得到最小二乘擬合多項式。對得到的最小二乘多項式進行一階微分計算，得到ΔV。

異常工況發(fā)生的概率密度函數(shù)f(t)如式(5)所示：

(5)

式(5)中，t為剩余時間；e為自然常數(shù)；λ為狀態(tài)參數(shù)變化過程中允許剩余時間倒數(shù)，其表達式如式(6)所示：

(6)

式(6)中，ΔVt為參數(shù)最大允許變化率，取正常工況下參數(shù)變化速率最大值，即max(ΔV)， 對一個確定的狀態(tài)參數(shù)而言，ΔVt是常數(shù)。

根據(jù)式(2)、式(5)、式(6)，異常工況發(fā)生概率P如式(7)所示：

(7)

1.3 參數(shù)異常的風險計算

異常工況產(chǎn)生的風險包括固有風險和趨勢風險。固有風險是狀態(tài)參數(shù)偏離過程中已經(jīng)產(chǎn)生的風險，其后果嚴重程度為損失函數(shù)的值，即后果嚴重程度為L，異常工況發(fā)生概率P=1，故根據(jù)風險計算公式，固有風險R1為：

(8)

趨勢風險用來表達由于參數(shù)仍具有向異常情況偏離的趨勢所帶來的風險，因損失函數(shù)的最大值為1，故異常狀況的最大后果值為1，則參數(shù)偏離的最大后果嚴重程度為I=1-L，根據(jù)式(7)可以得到異常狀況的發(fā)生概率P，故根據(jù)風險計算公式，參數(shù)的趨勢風險R2為：

(9)

根據(jù)參數(shù)固有風險和趨勢風險的計算公式，可以得到發(fā)生異常狀況的風險：

R=R1+R2

(10)

因此，根據(jù)以上方法，可以將裝置監(jiān)控系統(tǒng)實時狀態(tài)信息轉換為風險信息，實現(xiàn)實時風險評估。

2 風險預警方法實施步驟

風險預警方法的具體實施步驟如下：

(1)采集一個工作時段(約8 h)內的狀態(tài)參數(shù)值共5000組，存為數(shù)值矩陣B1×5000，該組數(shù)據(jù)的平均值記為M，參數(shù)的報警閾值為ylim，聯(lián)鎖報警閾值為Q。

(2)將Q作為造成最大損失的狀態(tài)參數(shù)值，M作為狀態(tài)參數(shù)目標值，則根據(jù)式(1)，令T=M， Δ=Q-M， 公式變形如下：

(11)

(3)設置風險閾值：根據(jù)ALARP原則和帕累托分布的分級策略，將風險不可接受值到最大風險值之間的區(qū)域認定為不可接受的高風險區(qū)域，占總風險范圍的20%[15]。因此，將風險不可接受值設置為風險報警閾值，得到風險閾值S為0.80。

(4)當狀態(tài)參數(shù)y達到聯(lián)鎖報警閾值Q， 即y=Q時，損失函數(shù)達到最大值，即L=1； 當y與目標值M相同，即y=M時，裝置處在最穩(wěn)定狀態(tài)，損失值為0，故L=0。 已知發(fā)生異常狀況的風險參數(shù)包括固有風險R1和趨勢風險R2， 當y達到報警閾值ylim， 且趨勢風險R2=0時，為防止漏報，風險值應大于等于風險閾值，即R1≥0.8； 當y未達到報警閾值ylim， 且趨勢風險R2=0時，為防止誤報，風險值應小于等于風險閾值，即R1≤0.8； 推理可得，當y=ylim時，R1=0.8， 又根據(jù)式(8)，R1=L， 故此時L=0.8， 即當y=ylim時，損失函數(shù)L=0.8。

根據(jù)以上分析，對于某些特定狀態(tài)參數(shù)值yi， 其理想損失值Li已知(例：當yi=Q時，根據(jù)前文，L=1， 事實上，在實際運用過程中，可以根據(jù)實際情況增加損失函數(shù)的特征點，提高損失函數(shù)的匹配度，筆者只針對通用的3個特征點進行計算)；又對于給定狀態(tài)參數(shù)值y=yi和形狀因子γ=γr， 運用式(11)可以計算出其實際損失值Lf。 因此，根據(jù)式(12)構造最小二乘方程組，求出使得偏差δ平方和最小的形狀因子γ， 從而得出最接近理想狀態(tài)的損失函數(shù)。

(12)

(5)計算參數(shù)所能引起的最大變化值：由于石油化工裝置參數(shù)波動明顯，運用最小二乘方法對數(shù)據(jù)集B1×5000進行擬合，并對擬合曲線進行一階微分計算，得到參數(shù)歷史數(shù)據(jù)的瞬時變化速率ΔVj。則得到最大允許變化速率ΔVt=max(ΔVj)。

(6)實時記錄狀態(tài)參數(shù)值y，將y帶入式(11)，得到參數(shù)的損失值為L，則根據(jù)式(8)，可以得到固有風險值R1=L×P， 其中P=1。

(8) 根據(jù)式(10)，可以計算裝置的實時風險值R=R1+R2：

因狀態(tài)參數(shù)變化趨勢具有雙向性，僅當狀態(tài)參數(shù)變化方向與偏離目標值方向相同，即ΔVi與(y-M)同號時，才將趨勢風險R2納入整體風險的計算。

根據(jù)計算出的風險值，結合風險閾值，就可以實現(xiàn)實時風險評估。

3 案例分析

筆者將超早期精確預警方法分別應用到石油化工企業(yè)聚丙烯裝置與氣分裝置的風險預警中，對該方法的有效性進行驗證。

3.1 氣分裝置風險預警

3.1.1 工藝介紹

氣分裝置的主要作用是利用氣體中各組分揮發(fā)度不同，對混合氣體進行分離，該裝置主要包含脫硫部分、脫硫醇部分和氣體分餾部分，其主要任務是通過精餾使混合液態(tài)烴得到分離，從而生產(chǎn)出高純度的丙烯以及其余部分，其中丙烯主要作為聚丙烯生產(chǎn)的原料，其余部分則成為液化氣或其他化工生產(chǎn)原料。

現(xiàn)場軟件通過DCS系統(tǒng)采集數(shù)據(jù)，并與預先設置的報警閾值進行比較，當參數(shù)超過報警閾值時則產(chǎn)生報警，圖1為現(xiàn)場DCS系統(tǒng)監(jiān)測圖。當報警發(fā)生時，留給現(xiàn)場操作人員的相應時間較短，需要采用合適的預警方法實現(xiàn)裝置的早期預警。

圖1 DCS系統(tǒng)監(jiān)測圖Fig.1 Monitoring chart of DCS system

3.1.2 脫乙烷塔塔底液位風險預警

圖2為氣分裝置脫乙烷塔塔底液位變化曲線，其報警閾值ylim=54.5%、報警聯(lián)鎖閾值Q=57%，狀態(tài)參數(shù)歷史數(shù)據(jù)平均值M=50%，實時狀態(tài)參數(shù)值為y，因此，根據(jù)式(11)，構造損失函數(shù)如下：

通過最小二乘法確定γ=3.024，將y帶入構造出的損失函數(shù)，計算出實時的損失值L，則根據(jù)式(8)，可以求得參數(shù)固有風險值R1=L×P，其中P=1。

狀態(tài)參數(shù)實時風險值R為：

R=R1+R2

圖2 脫乙烷塔塔底液位變化曲線Fig.2 Profile of deethane tower bottom liquid level

繪制裝置的實時風險曲線如圖3所示。

圖3 脫乙烷塔塔底液位風險曲線Fig.3 Risk profile of of deethane tower bottom liquid level

可以看出，風險值在695 s開始超出報警閾值，發(fā)出預警，于900 s左右快速上升并反復超出報警限。由圖3可知，該參數(shù)在1250 s時達到報警閾值觸發(fā)DCS報警。結果表明，該方法較DCS報警提前10 min左右。

3.1.3 脫乙烷塔塔頂壓力風險預警

選擇脫乙烷塔塔頂壓力作為目標參數(shù)，其參數(shù)變化曲線如圖4所示，報警閾值ylim=2.773 MPa，報警聯(lián)鎖閾值Q=2.785 MPa，歷史數(shù)據(jù)平均值M=2.75 MPa，實時狀態(tài)參數(shù)值為y，根據(jù)式(11)，構造損失函數(shù)如下：

通過最小二乘法確定γ=0.0133，將y帶入構造出的損失函數(shù)，計算出實時的損失值L，則根據(jù)式(8)，固有風險值R1=L×P， 其中P=1。

狀態(tài)參數(shù)實時風險值R為：

R=R1+R2

圖4 脫乙烷塔塔頂壓力變化曲線Fig.4 Profile of deethane tower top pressure

繪制裝置的風險曲線如圖5所示。

圖5 脫乙烷塔塔頂壓力風險曲線Fig.5 Profile of deethane tower top pressure

由圖5可得，風險曲線1470 s左右達到風險閾值，發(fā)出預警。由圖4知，狀態(tài)參數(shù)于1995 s產(chǎn)生DCS報警。結果表明，該方法較DCS報警提前8 min 左右。

3.2 聚丙烯裝置風險預警

3.2.1 聚丙烯裝置介紹

聚丙烯裝置是利用氣分裝置產(chǎn)出的丙烯氣生產(chǎn)聚丙烯產(chǎn)品的裝置系統(tǒng)，該裝置主要包含丙烯精制系統(tǒng)和聚合閃蒸系統(tǒng)，其中丙烯精制系統(tǒng)令原料丙烯通過各精制塔，通過各物理或化學作用，脫去其中的硫、水、氧、砷等雜質，為下一道聚合工序儲備原料。聚合閃蒸系統(tǒng)通過向聚合釜中投入催化劑、活化劑、二苯基甲基二甲氧基硅烷(DDS)、氫氣等助劑，用熱水升溫激活反應，從而生成聚丙烯，并回收未反應的丙烯至回收計量罐。聚丙烯裝置通過將DCS系統(tǒng)對現(xiàn)場異常狀況進行監(jiān)控，如圖6所示，筆者通過風險預警方法對聚丙烯裝置進行預警。

圖6 聚丙烯裝置DCS監(jiān)測圖Fig.6 Parameter alarms of polypropylene devices

3.2.2 熱水罐液位風險預警

選擇聚丙烯裝置熱水罐液位作為待測參數(shù)，其參數(shù)變化曲線如圖7所示，其報警閾值ylim=82%、報警聯(lián)鎖閾值Q=95%，歷史數(shù)據(jù)平均值M=60%，實時狀態(tài)參數(shù)值為y，因此，根據(jù)式(11)，構造損失函數(shù)如下：

通過最小二乘法確定γ=12.58，將y帶入構造出的損失函數(shù)，計算出實時的損失值L，則根據(jù)式(8)，固有風險值R1=L×P，其中P=1。

狀態(tài)參數(shù)實時風險值R為：

R=R1+R2

繪制裝置的風險圖如圖8所示。

由圖8可知，在1185 s時狀態(tài)參數(shù)風險值超標，此時DCS系統(tǒng)尚未達到報警閾值。這是由于參數(shù)上升速率過快，導致趨勢風險增大所引起的。1225 s 時參數(shù)風險值再次超標并在閾值附近波動。由圖7得，參數(shù)于 2170 s 時超過DCS報警閾值并發(fā)出報警。結果表明，該方法較DCS報警提前16 min 左右。

圖7 熱水罐液位變化曲線Fig.7 Profile of hot-water cylinder liquid level

圖8 熱水罐液位風險曲線Fig.8 Risk profile of hot-water cylinder liquid level

根據(jù)該方法應用于不同狀態(tài)參數(shù)的風險預警情況，將預警結果整理如下：

表1 不同狀態(tài)參數(shù)預警時間Table 1 Early warning time of different parameters

事實證明，該方法可以放大參數(shù)的早期偏差，實現(xiàn)參數(shù)的早期預警，根據(jù)案例分析，該方法對于不同的異常均能很好地檢出，并在參數(shù)到達報警閾值前10 min左右發(fā)出預警信息，對異常工況進行有效預警。

4 結 論

(1)針對石油化工裝置現(xiàn)有的預警技術誤報率較高，以及風險預警方法無法定量分析等不足，將損失函數(shù)應用到石油化工裝置的風險預警中，通過構造合適的損失函數(shù)，定量描述參數(shù)偏離所產(chǎn)生的預期損失，并根據(jù)剩余時間理論計算化工過程異常工況發(fā)生的概率，定量計算出參數(shù)偏離所產(chǎn)生的風險，從而實現(xiàn)異常工況的風險預警。

(2)將該方法分別應用到石油化工企業(yè)的聚丙烯裝置與氣分裝置的風險預警，通過計算相關參數(shù)的風險值，并與實際情況進行比對，從而對方法效果進行驗證。

(3)結果表明，隨著監(jiān)測參數(shù)接近報警限，裝置風險值明顯上升，且參數(shù)變化幅度越大，該風險值上升越快。這說明該風險預警模型能有效地將狀態(tài)參數(shù)的偏離和參數(shù)的變化趨勢轉化為實時的風險信息，對非合理狀態(tài)參數(shù)變化趨勢進行風險預警。該方法用于狀態(tài)參數(shù)發(fā)生偏離的初期辨識，可以提前10 min左右監(jiān)測到異常，對異常工況進行有效預警。