基于回歸分析的燃?xì)夤艿朗鹿事暑A(yù)測

摘 要

摘要：定量化和精細(xì)化管理是企業(yè)發(fā)展的必由之路，燃?xì)馄髽I(yè)老的管理模式是從事故到事故的被動搶修模式，新的管理模式是以風(fēng)險預(yù)評判為基礎(chǔ)的無事故預(yù)防型管理。風(fēng)險預(yù)評判只有建立
摘要：定量化和精細(xì)化管理是企業(yè)發(fā)展的必由之路，燃?xì)馄髽I(yè)老的管理模式是從事故到事故的被動搶修模式，新的管理模式是以風(fēng)險預(yù)評判為基礎(chǔ)的無事故預(yù)防型管理。風(fēng)險預(yù)評判只有建立在定量分析的水平上，才能真正具有系統(tǒng)性、精確性和預(yù)見性，所以管道事故數(shù)據(jù)的獲取和分析是定量風(fēng)險評估方法建立、改進(jìn)和完善的基礎(chǔ)?；诨貧w分析的方法預(yù)測管道事故率，事故率為因變量，自變量可以有很多種(管道內(nèi)徑、管道埋深、管壁厚度等)，本文以管徑為自變量從而分析出遞變性的規(guī)律，從最具普遍性的損害原因(第三方損害)著手，借鑒國外(歐洲)燃?xì)夤艿朗鹿式y(tǒng)計資料和分析方法，得出管道事故率回歸方程，體現(xiàn)影響因素與管道事故率之間的函數(shù)關(guān)系。
關(guān)鍵詞：城市燃?xì)?；天然氣；管道事故；定量風(fēng)險評估；回歸分析；第三方損害；事故統(tǒng)計；事故率預(yù)測
　1 引言
　定量化和精細(xì)化管理是企業(yè)發(fā)展的必由之路，管道定量風(fēng)險評估是提高管線運(yùn)行經(jīng)濟(jì)和安全水平的重要途徑[1]。風(fēng)險的識別和評價基于管道事故方面的數(shù)據(jù)信息，定量風(fēng)險評價技術(shù)正處在發(fā)展階段，歷史數(shù)據(jù)的發(fā)掘能夠揭示事故率和事故原因之間的定量關(guān)系。
　2 管道事故率統(tǒng)計的分類方法
　管道事故率的確定是一個復(fù)雜的問題，并且有很大的不確定性，造成管道事故的原因多種多樣，有的是與時間相關(guān)的(例：腐蝕、疲勞失效)，有的是與時間無關(guān)的(例：第三方損壞、超壓)。所以，管道事故率是設(shè)計參數(shù)、施工質(zhì)量、維護(hù)措施和運(yùn)行環(huán)境的函數(shù)[2]。
　管道事故率的預(yù)測問題可以根據(jù)導(dǎo)致事故的原因分類考慮，首先從最主要事故原因著手分析無疑是行之有效的方法。國內(nèi)外的天然氣管道事故統(tǒng)計表明第三方損害是泄漏事故發(fā)生的主要原因，以往的研究結(jié)果證明第三方損害的危害程度主要受以下參數(shù)影響：埋深、壁厚、位置、保護(hù)措施。歐洲燃?xì)夤艿朗鹿蕯?shù)據(jù)庫(European Gas pipeline Incident data Group，EGIG)報告中有較為詳細(xì)的分類統(tǒng)計。
　EGIG對于第三方損害的統(tǒng)計是基于以下參數(shù)分類的[3]：①管徑；②管壁厚；③管道地理位置；④管道埋深；⑤管道保護(hù)措施。
　需要強(qiáng)調(diào)的是這些管道事故的統(tǒng)計，管道壓力≥1.5MPa，包括管徑DN 100mm以下的管道。在管道事故統(tǒng)計的分類辦法中，按管徑分類可以方便的把統(tǒng)計結(jié)果應(yīng)用于其它地區(qū)各管徑系列管道的風(fēng)險分析。
　北京市中壓以上(包括中壓)天然氣管道事故分類可以采用的方法有[1]：
　① 按管徑分類：主要包括DN 200、300、400、500、700、1000mm；
?、?按埋深分類：≥1.5m(市政路)、≥1.2m：
?、?按管壁厚分類(次高壓A以上)：9.5mm、11.9mm：
?、?按保護(hù)措施分類：有電極保護(hù)、無電極保護(hù)、防腐層3層PE，其它防腐層。
　從以上的分類來看，按埋深、壁厚和保護(hù)措施的分類形式較少，以此為基礎(chǔ)都難以分析出遞變性的規(guī)律，所以按管徑分類的辦法，一方面分類比較詳細(xì)，另一方面可以為其它地域的統(tǒng)計數(shù)據(jù)引入或比較提供方便。
　3 歷史數(shù)據(jù)
　歷史數(shù)據(jù)的有效收集是進(jìn)行定量風(fēng)險評價的基礎(chǔ)，也是政府和企業(yè)制定安全規(guī)范和制定風(fēng)險決策的依據(jù)。現(xiàn)在世界上各個國家和專業(yè)技術(shù)公司都已認(rèn)識到基礎(chǔ)數(shù)據(jù)的極端重要性，已經(jīng)建立起了數(shù)百種由政府和企業(yè)資助的基礎(chǔ)數(shù)據(jù)庫[4]。
　歐洲燃?xì)夤艿朗鹿蕯?shù)據(jù)庫(EGIG)，包括丹麥、西班牙、比利時、芬蘭、挪威、法國、德國、意大利、瑞士、英國、捷克、葡萄牙、瑞典、愛爾蘭、奧地利，總計15個國家，總管長約13×104km(管道壓力≥1.5MPa.包括DN 100mm以下的管道)㈣。這個數(shù)據(jù)庫已經(jīng)在世界各地的燃?xì)夤艿腊踩治鲋袕V泛應(yīng)用，對提高管道安全發(fā)揮了作用。
　歐洲燃?xì)夤艿?970年至2007年間的統(tǒng)計數(shù)據(jù)根據(jù)以下泄漏的尺寸將管道破損進(jìn)行分類記錄[3]。
?、?/span>微孔/裂紋：缺陷直徑＜2cm；②孔洞：2cm≤缺陷直徑≤盾管道直徑；③斷裂：缺陷直徑＞管道直徑。
　M¼所示為EGIG第三方(外部干擾)造成的管道破損面積(事故情形)與管道內(nèi)徑之間的統(tǒng)計關(guān)系[3]，橫坐標(biāo)為管徑，縱高坐標(biāo)為事故率(1000km管道在1年內(nèi)發(fā)生事故的次數(shù)：次/(1000km·a))。對原始圖表[3](圖1)的觀測發(fā)現(xiàn)小管徑的受損率更高，這是因為較小的管徑的壁厚也較小，它被刮破的可能性也更大。
 
　4 回歸分析
　為了能夠利用歐洲燃?xì)夤艿朗鹿蕡蟾娴臄?shù)據(jù)，首先要對其報告中展示的柱狀圖進(jìn)行數(shù)據(jù)的二次生成，原始數(shù)據(jù)顯示對應(yīng)相應(yīng)的管徑范圍事故情形按破損程度分為三類：微孔、孔洞、斷裂。對原始圖表的數(shù)據(jù)二次生成如表1所示。利用回歸分析的方法進(jìn)行預(yù)測時，還需注意歷史數(shù)據(jù)的適用范圍，并根據(jù)現(xiàn)在考慮的具體情況作出相應(yīng)修整，特別是有些風(fēng)險事件尚未發(fā)生，事前很難對其作出準(zhǔn)確事故發(fā)生概率預(yù)測[4]。下面所述是針對尚未發(fā)生的小概率事件而采取的彌補(bǔ)辦法。圖1僅顯示了管徑小于1000mm的數(shù)據(jù)，因為更大管徑的事故尚未出現(xiàn)(事故率處在較低的水平，但并不等于零)，這些缺少的數(shù)據(jù)不能直接計算得到，需要用數(shù)理統(tǒng)計的方法確定這些缺少的數(shù)據(jù)。欠缺數(shù)據(jù)的取值范圍介于零和最后可獲得數(shù)值之間，最后可獲得數(shù)值是指按管徑分級的數(shù)據(jù)中的最后一個管徑系列的事故率觀測值，大于此管徑系列的其它管徑事故率觀測值為零，現(xiàn)在就是用最后可獲得值代替未獲得的數(shù)據(jù)[5]。
表1 第三方造成的損害——管徑與事故率

管徑范圍/英寸
管徑范圍(單位轉(zhuǎn)換mm)
事故情形分類/[次/1000km·a)
總計/[次/1000km·a)]
下限
上限
孔洞
微孔
斷裂
≤5
0+
127
O.181
0.257
0.124
0.562
5～10
127
254
0.067
O.190
0.057
0.314
11～17
279
432
0.038
0.071
O.024
0.133
17～23
432
584
0.019
0.024
O.010
0.052
23～29
584
737
0.006
0.008
0.006
0.020
29～35
737
889
—
—
0.005
0.005
35～41
889
1041
—
—
0.005
0.005
41～47
1041
1194
—
—
—
—
≥47
1194
1194+
—
—
—
—
注：符號“—”代表未能從原始數(shù)據(jù)中讀取的數(shù)據(jù)。原始數(shù)據(jù)的管徑是按英制單位分類的，這里為了便于分析，用mm對原有數(shù)據(jù)進(jìn)行了轉(zhuǎn)換。

　表2顯示了這種數(shù)據(jù)處理方法的結(jié)果(表中填充部分)。這種做法主要考慮欠缺數(shù)據(jù)的缺少是因為這些數(shù)值較小，在現(xiàn)有的運(yùn)營期限內(nèi)案例尚未出現(xiàn)，它的值介于零和最小觀測值之間，取最小觀測值相當(dāng)于取上限值，屬于按最不利情況的取值，防止統(tǒng)計分析的結(jié)果未能體現(xiàn)出可能出現(xiàn)的最危險狀況。

表2 第三方造成的損害——管徑與事故率欠缺數(shù)據(jù)的彌補(bǔ)

管徑范圍/英寸
管徑范圍(單位轉(zhuǎn)換mm)
事故情形分類/[次/(1000km·a)]
下限
上限
孔洞
微孔
斷裂
≤5
0+
127
0.257
0.181
0.124
5～10
127
254
0.19
0.067
0.057
11～17
279
432
0.071
0.038
0.024
17～23
432
584
O.024
0.019
0.010
23～29
584
737
O.008
0.006
0.006
29～35
737
889
0.008
0.006
0.005
35～41
889
1041
0.008
0.006
0.005
41～47
1041
1194
0.008
0.006
0.005
≥47
1194
”94+
0.008
O.006
0.005
注：表中斜體是經(jīng)過彌補(bǔ)的數(shù)據(jù)，用字框顯著標(biāo)出的數(shù)值是彌補(bǔ)時采用的數(shù)值。

　根據(jù)表2數(shù)據(jù)以管徑為自變量(橫坐標(biāo))，事故率為因變量(縱坐標(biāo))繪制散點(diǎn)圖，觀察散點(diǎn)圖發(fā)現(xiàn)事故率與管徑存在遞減的函數(shù)關(guān)系，可以利用曲線線性回歸分析的辦法進(jìn)行試算，可供選擇和函數(shù)關(guān)系有：冪函數(shù)y=xa、指數(shù)函數(shù)y=ax、對數(shù)函數(shù)y=logax。經(jīng)過試算發(fā)現(xiàn)冪函數(shù)的擬合效果最好，確定為事故率回歸分析的函數(shù)類型。

　回歸分析就是從某一組數(shù)據(jù)出發(fā)，判斷變量之間是否存在相關(guān)關(guān)系，確定它們之間的數(shù)學(xué)表達(dá)式即回歸方程，并對其可信程度作統(tǒng)計檢驗。然后利用確定的回歸方程進(jìn)行預(yù)測和控制[6～8]。

表3 管徑與第三方損害事故率

序號
1
2
3
4
5
6
7
8
9
管徑范圍/英寸
≤5
5～10
11～17
17～23
23～29
29～35
35～41
41～47
≥47
x管徑中間值/mm
127
191
356
508
660
813
965
1118
1194
y事故率
/[次/(1000km·a)]
0.562
0.314
O.133
0.052
O 02
0.019
O 019
0.019
0.019
注：表3中第三方損害事故率是指由表2中第三方損害“微孔”、“孔洞”、“斷裂”事故率的加和。

管徑作為自變量，事故率則作為因變量，通過回歸分析預(yù)測某一管徑的管道事故率范圍。管徑與事故率的關(guān)系基本呈倒冪函數(shù)關(guān)系，可用下式表示：
 
式中yi——管徑系列i的事故率，次/(1000km·a)
        a、b——系數(shù)
        xi——第i系列管道的管徑中間值，mm
 
　建立經(jīng)驗回歸方程的目的在于揭示兩個相關(guān)變量x與y之間的內(nèi)在規(guī)律，所建立的回歸方程是否有意義還要經(jīng)過回歸顯著性檢驗。用F檢驗對回歸方程作顯著性檢驗的方法稱為方差分析，其檢驗過程可由一張方差分析表來進(jìn)行[6～8]，見表4，計算案例如表5所示。
表4 方差分析表

方差來源
偏差平方和
自由度
方差
F值
Fα
顯著性
回歸
SR
1
VR=SR/1
F=VR/Ve
F0.05(1，n-2)
**
剩余
Se
n-2
Ve=Se/(n～2)
—
F0.05(1，n-2)
—
總和
ST
n-1
—

　若F＞F0.04(1，n-2)，則為高度顯著，在上表顯著性列標(biāo)“**”；F0.05(1，n-2)≤F≤F0.01(1，n-2)，則為顯著，在顯著性欄標(biāo)“*”；若F＜F0.05(1，n-2)則為不顯著，在顯著性欄不標(biāo)記號[6～8]。
　ST：總的偏差平方和，它反映y1，y2，…yx只的離散程度，ST=SR+Se；
　SR：回歸平方和，反映了由于x的變化而引起y的波動的大小，通也對y的相關(guān)性而引起的：
　Se：剩余平方和，反映出除了x對y的線性影響之外的剩余因素對y所引起的波動大小。
　若回歸方程有意義，即引起y波動主要是由x變化而引起的，其它一切因素是次要的，即要求SR盡可能大，而Se盡可能小?；貧w直線對樣本數(shù)據(jù)點(diǎn)的擬合程度，通常用可決系數(shù)R2來表示R2=SR/ST，R2越接近于1，表明回歸直線對樣本數(shù)據(jù)點(diǎn)的擬合程度越高[6～8]。
表5 管徑與第三方損害事故率方差分析表

方差來源
偏差平方和
自由度
方差
F值
Fα
顯著性
回歸
0.2841
1
0.2841
517.4457
5.591
**
剩余
0.0038
7
0.00054903
—
12.2464
—
總和
0.2879
8
—

　根據(jù)表5數(shù)據(jù)，通過回歸分析可以得到回歸方程：
　y=-0.074199+77.9809/x
　可決系數(shù)：R2=O.9933
　F檢驗[F＞F0.01(1，n-2)]：517.4457＞5.5914；顯著性：高度顯著。
　管徑與事故率的回歸方程通過了F檢驗與R檢驗，即可用來預(yù)測，所謂預(yù)測問題，就是在確定自變量的某一個X0值時求相應(yīng)的變勤的估計值。第三方損害的計算結(jié)果列于表6中，回歸圖如圖3～5所示。
表6 第三方損害一元回歸分析結(jié)果

y事故率/[次/(1000km·a)]
回歸方程
可決系數(shù)
F檢驗
F＞F0.01(1，n-2)
顯著性
F＞F0.01(1，n-2)為高度顯著
“微孔”破損
y=-0.023006+23.1641/x
R2=0.96777
103.359＞5.591
高度顯著
“孔洞”破損
y=-O.035285+38.4419/x
R2=O.98994
342.620＞5.591
高度顯著
“斷裂”破損
y=-0.015763+16.3614/x
R2=O.98145
183.412＞5.591
高度顯著

 
　5 應(yīng)用舉例
　某一地區(qū)的事故原因構(gòu)成如圖6所示，單位管長泄漏事故的發(fā)生率用ф2表示，簡稱“事故率”，單位：[次/(1000km·a)]。
 
　式中i——管道破損情形序號，i=1、2、3分別對應(yīng)管道破損微孔、孔洞、斷裂情形
　　j——事故原因序號，對于天然氣管道泄漏事故可以按照管道破損原因進(jìn)行劃分；j=1，2，3……，6分別對應(yīng)第三方損害、地表移動、操作失誤、腐蝕、材料缺陷、以及其它未知原因。
　該地區(qū)有燃?xì)夤艿劳度脒\(yùn)營的年數(shù)有限，所積累的管道事故數(shù)據(jù)不足，近年來其管道運(yùn)營管理水平與歐洲相似，可以利用EGIG確定出管道事故率的數(shù)量級。例如DN·500·mm的天然氣管道因第三方損害發(fā)生泄漏的事故率為(利用表6中的回歸方程計算，取x=500)：
 
　6 結(jié)語
　第三方損害管道事故率的分析表明，一方面，相對來說管徑越大事故率越低，這也與管道運(yùn)營管理的經(jīng)驗相符，大管徑因為機(jī)械強(qiáng)度更高而受損泄漏的概率更低。另一方面，管道事故數(shù)據(jù)的收集整理能夠為管道的定量風(fēng)險評價提供依據(jù)。
參考文獻(xiàn)：
[1] 杜學(xué)平.北京市天然氣輸配管網(wǎng)風(fēng)險評估的應(yīng)用研究(碩士學(xué)位論文)[D].北京：北京建筑工程學(xué)院.2010
[2] JO Y D，BUM J A A method of quantitative risk assessment for transmission pipeline carrying natural gas[J].Journal of Materials，2005，(8)1-12
[3] European Gas pipeline Incident Data Group(EGIG).The 7th EGIG report[EB/OL]. http：/SNWW egig.nl/nav/publication.htm.2008-12
[4] 周偉國.馬國彬能源工程管理[M].上海：同濟(jì)大學(xué)出版社，2007
[5] MATHER J，BLACKMOREC，PETRIEA，et al An assessment of measures in use for gas pipeline to mitigate against damage third party activity[M].Norwich(United Kingdom)：Health and Safety Executive(HSE)，2001
[6] 劉燕，張應(yīng)輝，杜學(xué)平.基于回歸分析的北京市天然氣供應(yīng)規(guī)模預(yù)測[J].煤氣與熱力，2010，30(6)：A26-A29
[7] 何燦芝，羅漢應(yīng)用數(shù)理統(tǒng)計[M].長沙：湖南大學(xué)出版社，2004
[8] 王巖，隋思連，王愛青數(shù)理統(tǒng)計與MATLAB工程數(shù)據(jù)分析[M].北京：清華大大學(xué)出版社，2006
 
(本文作者：杜學(xué)平1 劉燕1 劉蓉2 1.北京市燃?xì)饧瘓F(tuán)有限責(zé)任公司 北京 100035；2.北京建筑工程學(xué)院 北京 100044)