李敏睿，王海清，姚竣瀚，毛 奇

(中國石油大學(xué)(華東)安全科學(xué)與工程系，山東 青島 266580)

0 引言

國內(nèi)化工企業(yè)尤其是涉及“兩重點(diǎn)一重大”的大型企業(yè)按照中華人民共和國應(yīng)急管理部〔2014〕116號文《關(guān)于加強(qiáng)化工安全儀表系統(tǒng)管理的指導(dǎo)意見》的要求廣泛開展安全聯(lián)鎖系統(tǒng)評估與改造,以適應(yīng)當(dāng)前嚴(yán)峻的危化品安全生產(chǎn)形勢要求。安全完整性等級(SIL)，是衡量安全聯(lián)鎖系統(tǒng)安全水平的關(guān)鍵指標(biāo)之一，即在一定時(shí)間和一定條件下，安全相關(guān)系統(tǒng)執(zhí)行所規(guī)定的安全功能的概率[1]。根據(jù)功能安全標(biāo)準(zhǔn)IEC61508，SIL以聯(lián)鎖回路的平均需求時(shí)失效概率PFDavg作為驗(yàn)證依據(jù)。

功能安全標(biāo)準(zhǔn)IEC61508(下文統(tǒng)一簡稱為IEC61508)計(jì)算PFDavg時(shí)均假設(shè)相關(guān)設(shè)備的失效率為恒定值[2](即壽命分布服從指數(shù)分布)，以化工安全設(shè)備為背景的PFDavg算法和SIL驗(yàn)證中，安全管理人員均是采用恒定失效率計(jì)算。但實(shí)際中，設(shè)備由于老化[3]、疲勞磨損[4]等因素，絕大多數(shù)機(jī)械、電子設(shè)備[5]的失效率隨使用時(shí)間緩慢上升，在計(jì)算時(shí)導(dǎo)致PFDavg與實(shí)際值存在較大誤差。不僅如此，IEC61508針對不同應(yīng)用范圍，推薦不同參數(shù)估計(jì)置信度[6]，影響SIL驗(yàn)證的準(zhǔn)確性，增加驗(yàn)證的不確定性[7]，使得化工企業(yè)的生產(chǎn)裝置處于欠保護(hù)或者過保護(hù)的狀態(tài)，進(jìn)而可能誤導(dǎo)風(fēng)險(xiǎn)評判結(jié)果。

鑒于此，本文定量比較研究安全聯(lián)鎖回路在恒定和非恒定失效率情況下的2種算法，并對存在的偏差進(jìn)行定量評估分析，以更好落實(shí)《關(guān)于加強(qiáng)化工安全儀表系統(tǒng)管理的指導(dǎo)意見》在國內(nèi)化工企業(yè)的應(yīng)用。首先使用蒙特卡洛仿真，獲取符合工程實(shí)際的失效數(shù)據(jù)，進(jìn)而使用極大似然估計(jì)求解出威布爾參數(shù)和指數(shù)參數(shù)，并運(yùn)用威布爾分布PFDavg算法、指數(shù)分布PFDavg算法及IEC61508近似算法分別在70%和90%的2種典型置信度下得到3者PFDavg值隨預(yù)防性維修周期的變化趨勢，對SIL驗(yàn)證結(jié)果的影響進(jìn)行對比分析。

1 失效截尾數(shù)據(jù)的蒙特卡洛仿真

本文主要研究目的是對比安全聯(lián)鎖系統(tǒng)設(shè)備的失效數(shù)據(jù)組分別在IEC61508近似算法，指數(shù)分布算法與威布爾分布算法下的可靠性評估差異，在結(jié)合與參考了安全聯(lián)鎖設(shè)備失效數(shù)據(jù)庫、研究對象的大量歷史失效數(shù)據(jù)后，得出1組可信度較高的失效參數(shù)，后使用這組參數(shù)在蒙特卡洛仿真中模擬生成一系列失效數(shù)據(jù)。這種方法能夠滿足本文同時(shí)針對指數(shù)分布參數(shù)估計(jì)和威布爾分布參數(shù)估計(jì)的需求。根據(jù)蒙特卡洛仿真的原理[8]，基于MATLAB設(shè)計(jì)蒙特卡洛仿真算法流程，如圖1所示。

圖1 基于MATLAB蒙特卡洛仿真算法流程

蒙特卡洛仿真使用的仿真方法為概率積分變換方法，這種方法是將任意隨機(jī)變量轉(zhuǎn)換為均勻分布的隨機(jī)變量，代入概率密度分布函數(shù)的逆函數(shù)，循環(huán)多次計(jì)算得到一系列失效時(shí)間樣本值。

2 參數(shù)的極大似然估計(jì)

2.1 威布爾分布參數(shù)估計(jì)

假設(shè)失效數(shù)據(jù)服從威布爾分布，運(yùn)用極大似然法和試錯(cuò)法求出威布爾分布的形狀參數(shù)α和尺度參數(shù)β的估計(jì)值。使用極大似然法[9]求解威布爾參數(shù)估計(jì)值，如式(1)所示：

(1)

(2)

(3)

將式(3)帶入式(2)中，如式(4)～(5)所示：

(4)

(5)

2.2 指數(shù)分布的參數(shù)估計(jì)

假設(shè)失效數(shù)據(jù)服從指數(shù)分布，可由指數(shù)分布I型截尾數(shù)據(jù)下的極大似然法求出失效率λ，指數(shù)分布的極大似然估計(jì)函數(shù)，如式(6)所示：

(6)

(7)

2.3 參數(shù)的置信區(qū)間估計(jì)

根據(jù)IEC61508，參數(shù)估計(jì)的置信度η不低于70%，同時(shí)在用于判斷結(jié)構(gòu)約束時(shí)要求置信度η不低于90%[10],在極大似然法估計(jì)的基礎(chǔ)上，對參數(shù)進(jìn)行不同置信區(qū)間估計(jì)以得出區(qū)間估計(jì)值。分別以70%和90%置信度為條件，研究和對比可靠性參數(shù)的區(qū)間估計(jì)值及其影響。

依據(jù)指數(shù)分布參數(shù)的區(qū)間估計(jì)算法[9]，可計(jì)算出指數(shù)參數(shù)估計(jì)區(qū)間[λL,λU]，如式(8)～(9)所示：

(8)

(9)

式中：η為區(qū)間估計(jì)的置信度；λL為指數(shù)參數(shù)區(qū)間估計(jì)下限值；λU為指數(shù)參數(shù)區(qū)間估計(jì)上限值；T為所有樣本的總失效時(shí)長，h；χ2為卡方分布。

依據(jù)威布爾分布參數(shù)的區(qū)間估計(jì)算法[9,11]，可計(jì)算出不同置信度下形狀參數(shù)和尺度參數(shù)的估計(jì)區(qū)間，對于α的估計(jì)區(qū)間[αL,αU]，如式(10)～(11)所示：

(10)

(11)

式中：αL為形狀參數(shù)的區(qū)間估計(jì)下限值；αU為形狀參數(shù)的區(qū)間估計(jì)上限值。

對于β估計(jì)區(qū)間[βL,βU]，如式(12)～(13)所示：

(12)

(13)

式中：βL為尺度參數(shù)區(qū)間估計(jì)下限值；βU為尺度參數(shù)區(qū)間估計(jì)上限值；d1和d2分別為中間系數(shù)。

對于d1和d2計(jì)算，如式(14)～(18)所示：

(14)

(15)

(16)

(17)

(18)

式中：A3,A4,A5和A6分別為區(qū)間估計(jì)的中間系數(shù)；x為標(biāo)準(zhǔn)正態(tài)分布的1-η分位點(diǎn)。

3 化工安全聯(lián)鎖設(shè)備的PFDavg算法和SIL驗(yàn)證

工程中為保證安全聯(lián)鎖系統(tǒng)在發(fā)生緊急情況時(shí)正常運(yùn)行，定期執(zhí)行設(shè)備的預(yù)防性維修計(jì)劃(PM)，在這個(gè)周期里進(jìn)行1次檢修和維護(hù)。記T0為設(shè)備的預(yù)防性維修周期，給出求解緊急放空執(zhí)行單元PFDavg的3種算法[12]。

(19)

(20)

(21)

4 案例分析

以某蠟油加氫裝置的典型安全聯(lián)鎖系統(tǒng)回路[13]的緊急放空閥(EDV)為例(注：視為1oo1結(jié)構(gòu)，現(xiàn)場實(shí)際的EDV為1用1備)。該回路中，催化加氫反應(yīng)器的頂部氣相產(chǎn)物，經(jīng)過熱高壓分離器分離的熱高分氣體進(jìn)入熱高分氣換熱至184 ℃，再經(jīng)過熱高分器空冷器A01001冷卻至50 ℃后進(jìn)入冷高壓分離器。自冷高壓分離器頂部出來的循環(huán)氫經(jīng)脫硫塔入口分液罐分液后，進(jìn)入D03001循環(huán)氫脫硫塔底部，緊急放空閥的聯(lián)鎖回路，如圖2所示。

圖2 緊急放空閥的聯(lián)鎖回路

該回路中，由壓力變送器PT實(shí)時(shí)監(jiān)測冷高壓分離器內(nèi)壓力，并將數(shù)據(jù)傳至邏輯解算器LS單元，當(dāng)冷高壓分離器內(nèi)壓強(qiáng)超過2.1 MPa時(shí)，由邏輯解算器LS提供觸發(fā)信號給緊急放空閥，打開閥門以0.7 MPa/min的泄壓速率將超壓氣體排出。

結(jié)合SIS失效數(shù)據(jù)庫中緊急放空閥的失效統(tǒng)計(jì)數(shù)據(jù)[14]，以及參考實(shí)際工程中服從威布爾分布的緊急放空閥失效參數(shù)，設(shè)置形狀參數(shù)為1.49，尺度參數(shù)為6 025 003，通過蒙特卡洛仿真模擬生成200組隨機(jī)失效數(shù)據(jù)。

根據(jù)加氫裝置實(shí)際運(yùn)行的情況，在得到模擬樣本后，進(jìn)行截尾時(shí)長為175 200 h的I型截尾處理，得到150組完整失效數(shù)據(jù)和50組截尾失效數(shù)據(jù)。需要強(qiáng)調(diào)的是，采用蒙特卡洛模擬的截尾失效數(shù)據(jù)作為測試數(shù)據(jù)，并不影響比較統(tǒng)計(jì)分析的有效性，且可以比實(shí)際的工程采集數(shù)據(jù)提供更客觀的誤差分析結(jié)果[15]。

使用試錯(cuò)法，先是試錯(cuò)α估計(jì)值區(qū)間的選取，考慮到設(shè)備運(yùn)行周期內(nèi)失效率是1個(gè)緩慢上升的過程，并查閱相關(guān)閥的失效參數(shù)庫，得知α值在1.0至1.7范圍內(nèi)，因此試錯(cuò)值假設(shè)在該區(qū)間，依次帶入式(4)，求得α估計(jì)值偏差量，并求得似然估計(jì)求解的α和β估計(jì)值，如表1所示。

表1 似然估計(jì)求解α和β估計(jì)值

由表1計(jì)算得到的數(shù)據(jù)，繪制α估計(jì)值偏差量與β估計(jì)值雙坐標(biāo)曲線，如圖3所示。

圖3 α估計(jì)值偏差量與β估計(jì)值曲線

根據(jù)式(14)～(18)，威布爾區(qū)間估計(jì)可算出70%置信度時(shí)d1=0.124，d2=-0.116，90%置信度時(shí)d1=0.198，d2=-0.198，代入式(12)～(13)中，算出β估計(jì)區(qū)間的上限值和下限值。

表2 威布爾分布和指數(shù)分布的參數(shù)估計(jì)區(qū)間

表3 70%和90%置信度時(shí)和估計(jì)值

分別在70%和90%這2種不同的參數(shù)估計(jì)置信度下，研究不同T0時(shí)，威布爾分布，指數(shù)分布和IEC61508的PFDavg計(jì)算結(jié)果的變化曲線，并分析這3種算法結(jié)果的差異。結(jié)合工程實(shí)際，為便于比較將安全聯(lián)鎖設(shè)備的T0分別取1，2，3，4，5 a。得到70%和90%置信度時(shí)3 種算法的PFDavg變化趨勢。70%置信度下3種算法的PFDavg與T0關(guān)系，如圖4所示；90%置信度下3種算法的PFDavg與T0關(guān)系，如圖5所示。

圖4 70%置信度下3種算法的PFDavg與T0的關(guān)系

圖5 90%置信度下3種算法的PFDavg與T0的關(guān)系

由圖4～5可知，70%和90%置信度下指數(shù)分布和IEC61508的PFDavg曲線差值均較小，證明IEC61508簡化算法替代指數(shù)分布算法的合理性。

鑒于指數(shù)分布和IEC61508的2種算法的計(jì)算結(jié)果相差不大，因此將僅在70%和90%置信度時(shí)分析指數(shù)分布算法與威布爾分布算法的PFDavg計(jì)算值之間的相對誤差和絕對誤差，并進(jìn)行對比。指數(shù)分布與威布爾分布PFDavg的相對誤差分析，如圖6所示；指數(shù)分布與威布爾分布PFDavg的絕對誤差分析，如圖7所示。

圖6 指數(shù)分布與威布爾分布PFDavg的相對誤差分析

圖7 指數(shù)分布與威布爾分布PFDavg的絕對誤差分析

IEC61508算法的PFDavg值隨著T0的增大，絕對誤差(即IEC61508和指數(shù)分布PFDavg相比于威布爾分布PFDavg的數(shù)值偏差)增大，而相對誤差值減小。絕對誤差作為SIL驗(yàn)證中數(shù)值誤差的直接表現(xiàn)，表明在指定預(yù)防性維修周期內(nèi)存在較大誤差，證明工程上壽命服從威布爾分布的設(shè)備，簡化為采用指數(shù)分布算法計(jì)算其PFDavg的不合理性，并將導(dǎo)致顯著的風(fēng)險(xiǎn)研判偏差。

指數(shù)分布和IEC61508算法的SIL驗(yàn)證結(jié)果基本一致，2者的驗(yàn)證結(jié)果不再區(qū)分考慮。根據(jù)IEC61508中對SIL驗(yàn)證的等級劃分規(guī)則，得出70%置信度時(shí)2種PFDavg算法的SIL驗(yàn)證結(jié)果，如表4所示；90%置信度時(shí)2種PFDavg算法的SIL驗(yàn)證結(jié)果，如表5所示。

表4 70%置信度時(shí)2種PFDavg算法的SIL驗(yàn)證結(jié)果

表5 90%置信度時(shí)2種PFDavg算法的SIL驗(yàn)證結(jié)果

用模擬的方法得出，對于通過SIL驗(yàn)證等級為SIL2的設(shè)備，如果實(shí)際壽命服從威布爾分布，則實(shí)際SIL驗(yàn)證是有可能達(dá)到SIL3的。在表4～5中，預(yù)防性維修周期為1～3 a時(shí)，由威布爾分布算法得到的SIL驗(yàn)證結(jié)果普遍比指數(shù)分布和IEC61508算法的SIL驗(yàn)證結(jié)果高1個(gè)等級，可見服從威布爾分布的設(shè)備在使用IEC61508算法進(jìn)行計(jì)算時(shí)將存在相當(dāng)大的偏差，導(dǎo)致化工企業(yè)SIS升級改造項(xiàng)目投資過大。70%與90%置信度下的威布爾分布算法的SIL驗(yàn)證結(jié)果存在一定差異，考慮到高置信度下計(jì)算得到的PFDavg值高于低置信度計(jì)算PFDavg值，在PFDavg數(shù)值計(jì)算使用70%置信度，結(jié)果較為折中，而使用90%置信度結(jié)果較為保守，驗(yàn)證了IEC61508中，PFDavg算法推薦置信度為70%，數(shù)據(jù)采集和處理推薦置信度為90%的這一建議基本合理。

5 結(jié)論

1)對于聯(lián)鎖回路的執(zhí)行單元等容易發(fā)生機(jī)械疲勞磨損的設(shè)備，其失效率在使用過程中緩慢上升，因此宜使用威布爾分布對設(shè)備失效時(shí)間數(shù)據(jù)進(jìn)行分析和PFDavg計(jì)算；此外本文案例雖考慮的是1oo1結(jié)構(gòu)，整個(gè)計(jì)算流程亦適用于其他表決結(jié)構(gòu)。

2)如果對于壽命服從威布爾分布的設(shè)備，使用IEC61508算法計(jì)算的PFDavg可能會造成SIL驗(yàn)證的誤差，低估了設(shè)備可以達(dá)到的SIL級別，造成化工廠的SIS改造投資增大。

3)區(qū)間估計(jì)的置信度越高，IEC61508算法計(jì)算的PFDavg與實(shí)際PFDavg偏差越小。隨著預(yù)防性維修周期的延長，在相同置信度下IEC61508的PFDavg與實(shí)際PFDavg間的相對誤差減小，絕對誤差增大。因此化工企業(yè)的允許短周期維護(hù)的聯(lián)鎖回路策略，采用IEC61508以及高置信度參數(shù)估計(jì)的驗(yàn)證結(jié)果基本可信。