基于多元線性回歸分析的汽油辛烷值損失預(yù)測(cè)建模

徐宗煌

(福州理工學(xué)院 應(yīng)用科學(xué)與工程學(xué)院，福建 福州 350506)

汽油是小型車輛的主要燃料，汽油燃燒產(chǎn)生的尾氣排放對(duì)大氣環(huán)境有重要影響.為此，世界各國(guó)都制定了日益嚴(yán)格的汽油質(zhì)量標(biāo)準(zhǔn).汽油清潔化重點(diǎn)是降低汽油中的硫、烯烴含量，同時(shí)盡量保持其辛烷值.辛烷值是反映汽油燃燒性能的最重要指標(biāo)，并作為汽油的商品牌號(hào).現(xiàn)有技術(shù)在對(duì)催化裂化汽油進(jìn)行脫硫和降烯烴過(guò)程中，普遍降低了汽油辛烷值.辛烷值每降低1個(gè)單位，相當(dāng)于損失約150元/t.以一個(gè)100萬(wàn)t/年催化裂化汽油精制裝置為例，若能將RON損失降低0.3個(gè)單位，其經(jīng)濟(jì)效益將達(dá)到4 500萬(wàn)元.

目前，國(guó)內(nèi)外對(duì)汽油辛烷值的預(yù)測(cè)進(jìn)行了一系列的相關(guān)研究.李煒等[1]針對(duì)汽油樣本數(shù)據(jù)較少的因素，提出了一種基于PSO-LSSVM的汽油辛烷值預(yù)測(cè)模型，并應(yīng)用在汽油實(shí)際調(diào)合系統(tǒng)中.周小偉等[2]分別利用BP神經(jīng)網(wǎng)絡(luò)和多元線性回歸分析，建立了二次反應(yīng)清潔汽油的辛烷值預(yù)測(cè)模型，同時(shí)進(jìn)行了模型驗(yàn)證和對(duì)比分析.韓仲志等[3]基于汽油樣品的近紅外光譜數(shù)據(jù)，采用五折交叉驗(yàn)證法，分析比較了最小均方二乘法(PLS)、神經(jīng)網(wǎng)絡(luò)(ANN)和支持向量機(jī)(SVM)的汽油辛烷值預(yù)測(cè)模型的穩(wěn)定性與精確性.CHEN等[4]為了提高汽油辛烷值的預(yù)測(cè)精度，建立了一種改進(jìn)的基于極限學(xué)習(xí)的模型，并基于汽油辛烷值標(biāo)準(zhǔn)數(shù)據(jù)庫(kù)驗(yàn)證了該模型的有效性.

現(xiàn)階段國(guó)內(nèi)外大多數(shù)學(xué)者只研究了汽油辛烷值的預(yù)測(cè)模型，而未能給出辛烷值損失的相關(guān)預(yù)測(cè)模型.鑒于此，本文基于催化裂化汽油精制裝置所采集的325個(gè)數(shù)據(jù)樣本(每個(gè)數(shù)據(jù)樣本都有354個(gè)操作變量)，運(yùn)用R型聚類法得到建模的25個(gè)主要變量，并利用數(shù)據(jù)挖掘技術(shù)，以辛烷值損失為因變量，采用十折交叉驗(yàn)證法建立基于多元線性回歸分析的汽油辛烷值(RON)損失的預(yù)測(cè)模型.

1 數(shù)據(jù)預(yù)處理

1.1 研究思路

原始數(shù)據(jù)來(lái)自于中石化高橋石化實(shí)時(shí)數(shù)據(jù)庫(kù)及LIMS實(shí)驗(yàn)數(shù)據(jù)庫(kù)，采集時(shí)間為2017年4月至2020年5月，采集操作位點(diǎn)共354個(gè).原始數(shù)據(jù)中，大部分變量數(shù)據(jù)正常，但每套裝置的數(shù)據(jù)均有部分位點(diǎn)存在問(wèn)題：部分變量只含有部分時(shí)間段的數(shù)據(jù)，部分變量的數(shù)據(jù)全部或部分為空值.因此，我們分別針對(duì)只含有部分時(shí)間點(diǎn)、數(shù)據(jù)全部為空值、不在變量操作范圍內(nèi)的樣本利用拉依達(dá)準(zhǔn)則(3σ準(zhǔn)則)[5—6]去除異常值等數(shù)據(jù)預(yù)處理方法，對(duì)285號(hào)和313號(hào)樣本原始數(shù)據(jù)進(jìn)行預(yù)處理.

1.2 預(yù)處理過(guò)程

1.2.1只含有部分時(shí)間點(diǎn)的位點(diǎn) 對(duì)于只含有部分時(shí)間點(diǎn)的位點(diǎn)，若其殘缺數(shù)據(jù)較多，無(wú)法補(bǔ)充，則將此類位點(diǎn)刪除.

通過(guò)對(duì)285號(hào)和313號(hào)樣本原始數(shù)據(jù)分析可知，285號(hào)數(shù)據(jù)樣本和313號(hào)數(shù)據(jù)樣本的位點(diǎn)數(shù)據(jù)均存在全部空值的情況(表1). 由于表1中11個(gè)變量全部時(shí)間點(diǎn)的數(shù)據(jù)均為0，因此將該變量直接剔除.

表1 時(shí)間不連續(xù)數(shù)據(jù)處理過(guò)程記錄

1.2.2部分?jǐn)?shù)據(jù)為空值的位點(diǎn) 對(duì)于部分?jǐn)?shù)據(jù)為空值的位點(diǎn)，空值處用其前后2 h數(shù)據(jù)的平均值代替.

在原始數(shù)據(jù)中，285號(hào)數(shù)據(jù)樣本的位點(diǎn)數(shù)據(jù)不存在部分空值的情況；而313號(hào)數(shù)據(jù)樣本存在部分?jǐn)?shù)據(jù)為空值的位點(diǎn)(表2). 將該變量在數(shù)據(jù)為空值的時(shí)間點(diǎn)直接剔除，最后再取平均值，作為313號(hào)數(shù)據(jù)樣本在該操作變量的數(shù)據(jù).

表2 時(shí)間不連續(xù)數(shù)據(jù)處理過(guò)程記錄

1.2.3不在變量操作范圍內(nèi)的樣本 根據(jù)工藝要求與操作經(jīng)驗(yàn)，總結(jié)出原始數(shù)據(jù)變量的操作范圍，然后采用最大最小的限幅方法剔除一部分不在此范圍內(nèi)的樣本.

首先，根據(jù)354個(gè)操作變量信息的取值范圍，與285號(hào)和313號(hào)樣本原始數(shù)據(jù)進(jìn)行比較，可知

(ⅰ) 變量編號(hào)為22、位號(hào)為S-ZORB.AT_5201.PV的取值范圍為0～5 μg/g.

在313號(hào)樣本原始數(shù)據(jù)中，其變量值處于4.787 2～9.432 3，符合標(biāo)準(zhǔn)取值范圍0～5的只有在2017-05-15 07:48:00時(shí)間點(diǎn)1個(gè)原始數(shù)據(jù)，其余均不在此范圍內(nèi)，因此將該變量直接剔除.

此外，在325個(gè)樣本數(shù)據(jù)中，其變量值處于-0.320 2～32.234 2，符合標(biāo)準(zhǔn)取值范圍0～5的有178個(gè)，約占總樣本的54.77%.

(ⅱ) 變量編號(hào)為43、位號(hào)為S-ZORB.PDC_2502.PV的取值范圍為-1～55 kPa.

在313號(hào)樣本原始數(shù)據(jù)中，其變量值處于-4.755 5～99.710 6，符合標(biāo)準(zhǔn)取值范圍-1～55 Pa的只有20個(gè)原始數(shù)據(jù)，其余均不在此范圍內(nèi)，因此將該變量直接剔除.

(ⅲ) 變量編號(hào)為49、位號(hào)為S-ZORB.SIS_LT_1001.PV的取值范圍為40%～80%.

在285號(hào)樣本原始數(shù)據(jù)中，其變量值處于1 093 484～1 093 667；而在313號(hào)樣本原始數(shù)據(jù)中，其變量值處于1 235 356～1 235 539.該變量全部值均不在標(biāo)準(zhǔn)取值范圍內(nèi)，因此將該變量直接剔除.

(ⅳ)變量編號(hào)為84、位號(hào)為S-ZORB.AI_2903.PV的取值范圍為0.5%～3%.

在285號(hào)樣本原始數(shù)據(jù)中，其變量值處于0.378 548 6～0.380 829；而在313號(hào)樣本原始數(shù)據(jù)中，其變量值均為5.089 727.該變量全部值均不在標(biāo)準(zhǔn)取值范圍內(nèi)，因此亦將該變量直接剔除.

(ⅴ) 變量編號(hào)為111、位號(hào)為S-ZORB.FT_1204.TOTAL的取值范圍為45 000～2 500 000.

在285號(hào)樣本原始數(shù)據(jù)中，其變量值處于7 071 433～7 072 616；而在313號(hào)樣本原始數(shù)據(jù)中，其變量值處于7 989 003～7 990 187.該變量全部值均不在標(biāo)準(zhǔn)取值范圍內(nèi)，因此將該變量直接剔除.

(1)

式中σ為標(biāo)準(zhǔn)誤差.該公式為貝塞爾公式[7].若某個(gè)測(cè)量值xb的剩余誤差vb(b∈[1,n])，滿足

(2)

則認(rèn)為xb是含有粗大誤差值的壞值，應(yīng)予剔除.

這里利用3σ準(zhǔn)則對(duì)325個(gè)樣本數(shù)據(jù)的14個(gè)非操作變量進(jìn)行處理，通過(guò)MATLAB編程，可找出14個(gè)非操作變量中的異常值并剔除(表3).

表3 利用拉依達(dá)準(zhǔn)則(3σ準(zhǔn)則)去除的異常值

1.3 預(yù)處理后的數(shù)據(jù)

經(jīng)過(guò)以上詳細(xì)分析可知：

(ⅰ) 對(duì)于只含有部分時(shí)間點(diǎn)的位點(diǎn)，編號(hào)為25、46、74、96、110、151、205、208、216、299、332的11個(gè)變量全部時(shí)間點(diǎn)的數(shù)據(jù)均為0，將這11個(gè)變量剔除；(ⅱ) 對(duì)于不在變量操作范圍內(nèi)的樣本，編號(hào)為22、43、49、84、111的5個(gè)變量全部值均不在標(biāo)準(zhǔn)取值范圍內(nèi)，將這5個(gè)變量剔除；(ⅲ) 由于產(chǎn)品性質(zhì)中的RON損失不屬于變量，因此也將此變量剔除.

所以最終剩下的數(shù)據(jù)變量為350個(gè).

2 建模主要變量

2.1 研究思路

工程技術(shù)應(yīng)用中經(jīng)常使用先降維后建模的方法，這有利于忽略次要因素，發(fā)現(xiàn)并分析影響模型的主要變量與因素.系統(tǒng)聚類法中的R型聚類法可以在變量眾多的情況下，通過(guò)研究各個(gè)變量之間的相似關(guān)系并按相似關(guān)系進(jìn)行聚類后，根據(jù)相似度的高低,找出影響系統(tǒng)的主要因素.因此，我們根據(jù)預(yù)處理后的數(shù)據(jù)，將剩余的變量利用R型聚類法進(jìn)行聚類分析(為了工程應(yīng)用方便，建議降維后的主要變量在30個(gè)以下).我們選擇聚類數(shù)為25，即將降維后的主要變量控制在25個(gè)，然后從每類中選擇1個(gè)最有代表性的變量，這些變量即為建立降低辛烷值損失模型具有代表性、獨(dú)立性的主要變量.

2.2 R型聚類法

2.2.1基本思想 在實(shí)際工作中，為了避免漏掉某些重要因素，往往在一開始選取指標(biāo)時(shí)盡可能考慮所有的相關(guān)因素，而這樣做的結(jié)果則是變量過(guò)多，變量間的相關(guān)度較高，給統(tǒng)計(jì)分析與建模帶來(lái)極大不便.因此人們希望能夠研究變量間的相似關(guān)系，按照變量的相似關(guān)系將其聚合成若干類，進(jìn)而找出影響系統(tǒng)的主要因素，引入了R型聚類法[8—9].其主要思想[10—11]為：根據(jù)系統(tǒng)的多個(gè)指標(biāo)變量，通過(guò)找出一些能夠具體度量指標(biāo)變量間相似程度的統(tǒng)計(jì)量，然后以這些統(tǒng)計(jì)量為劃分類型的依據(jù)，將各個(gè)相似程度較高的指標(biāo)變量聚合為一類，直到把所有的指標(biāo)變量聚合完畢.

2.2.2變量相似性度量 在對(duì)變量進(jìn)行聚類分析時(shí)，首先要確定變量之間的相似性度量[12]，常用的變量相似性度量方法有相關(guān)系數(shù)法和夾角余弦法，其計(jì)算公式如下.

(ⅰ) 相關(guān)系數(shù)法.

記變量xj的取值為(x1j,x2j,…,xnj)T∈Rn，j=1,2,…,m，則可將變量xj與xk的樣本相關(guān)系數(shù)作為其相似性度量：

(3)

在對(duì)指標(biāo)變量進(jìn)行R型聚類分析時(shí)，最常利用相關(guān)系數(shù)矩陣確定變量之間的相似性.

(ⅱ) 夾角余弦法.

(4)

對(duì)于任意的j,k都有|rjk|≤1，且rjk=rkj.此外，|rjk|越接近于1，表示xj與xk之間的相關(guān)性越大，即越相似；|rjk|越接近于0，則表示兩變量之間的相似性越弱.

2.2.3變量相似性度量 利用R型聚類法，根據(jù)指標(biāo)變量之間的相關(guān)性進(jìn)行變量聚類分析[13]時(shí)，最常用的方法是最大系數(shù)法和最小系數(shù)法，其計(jì)算公式如下.

(ⅰ) 最大系數(shù)法.

在最大系數(shù)法中，定義兩類指標(biāo)變量的距離為

(5)

式中：G在R型聚類分析中通常表示指標(biāo)變量的類；R1(G1,G2)表示兩類指標(biāo)變量中最相似的兩變量xj與xk的相似性度量值.

(ⅱ) 最小系數(shù)法.

在最小系數(shù)法中，定義兩類指標(biāo)變量的距離為

(6)

式中：R2(G1,G2)表示兩類指標(biāo)變量中相似性最小的兩變量xj與xk之間的相似性度量值.

2.3 編程求解

我們選擇聚類數(shù)為25，通過(guò)MATLAB編程，利用相關(guān)系數(shù)法和最大系數(shù)法,分別確定指標(biāo)變量之間的相似性并選取具有代表性的變量，聚類樹形圖如圖1所示.

最終得到25個(gè)最有代表性的指標(biāo)變量(表4)，作為建立降低辛烷值損失模型的主要變量.

3 辛烷值損失預(yù)測(cè)模型

3.1 研究思路

首先，將得到的25個(gè)最有代表性的指標(biāo)變量作為建立降低辛烷值損失模型的主要變量. 建立辛烷值損失預(yù)測(cè)模型需要采用數(shù)據(jù)挖掘技術(shù)，該模型以辛烷值損失為因變量，因此可建立多元線性回歸分析模型進(jìn)行預(yù)測(cè).其次，采用M折交叉驗(yàn)證法對(duì)模型進(jìn)行驗(yàn)證，利用隨機(jī)函數(shù)打亂25個(gè)主要變量的數(shù)據(jù)順序，再將數(shù)據(jù)按均勻分布隨機(jī)分成M份，其中(M-1)份數(shù)據(jù)作為訓(xùn)練集，剩下的1份作為測(cè)試集，依次重復(fù)進(jìn)行M次，在M次的預(yù)測(cè)結(jié)果中挑選預(yù)測(cè)精度最高的1次，作為最終所建立的辛烷值損失預(yù)測(cè)模型，并對(duì)回歸模型和回歸系數(shù)分別進(jìn)行假設(shè)檢驗(yàn).

圖1 指標(biāo)變量聚類樹形圖

表4 模型的25個(gè)主要變量

3.2 模型建立

一元線性回歸分析模型,主要是研究一個(gè)作為自變量的主要影響因素以解釋因變量的變化.而實(shí)際問(wèn)題研究中,一個(gè)變量往往會(huì)受到多個(gè)變量的影響，即辛烷值的損失量與所得25個(gè)主要變量的變化關(guān)系密切.因此，我們建立多元線性回歸模型[14—16](multivariable linear regression model)，研究辛烷值的損失隨25個(gè)主要變量變化的情況，其模型表達(dá)式為

(7)

式中：β0,β1,β2,…,βm,σ2均為與x1,x2,…,xm無(wú)關(guān)的未知參數(shù)；β0,β1,β2,…,βm稱為模型的回歸系數(shù)，用來(lái)描述某一自變量變化引起的因變量變化程度，且β0為常數(shù)項(xiàng).

現(xiàn)得到n個(gè)獨(dú)立觀測(cè)數(shù)據(jù)(yi,xi1,xi2,…,xim)，由 (7) 式可得

(8)

不妨記

(9)

則 (7) 式可表示為

(10)

式中En為n階單位矩陣.

(11)

達(dá)到最小.因此，不妨令

(12)

可得

(13)

則正規(guī)方程組的矩陣解為

(14)

(15)

而這組數(shù)據(jù)的擬合值為

.

(16)

3.2.2M折交叉驗(yàn)證法[17—18]采用M折交叉驗(yàn)證法對(duì)模型進(jìn)行驗(yàn)證.利用隨機(jī)函數(shù)打亂具有25個(gè)主要變量的數(shù)據(jù)順序，再將數(shù)據(jù)按均勻分布隨機(jī)分成M份.不妨取M=10，則挑選其中9份數(shù)據(jù)作為訓(xùn)練集，剩下的1份作為測(cè)試集.依次重復(fù)進(jìn)行10次，根據(jù)訓(xùn)練集訓(xùn)練出模型.將該模型放到測(cè)試集上，得到10次的預(yù)測(cè)結(jié)果.挑選預(yù)測(cè)精度最高的1次，作為最終所建立的辛烷值損失預(yù)測(cè)模型.該方法可以充分利用所有樣本數(shù)據(jù)，其主要步驟如下：

算法1 十折交叉驗(yàn)證法算法

Step 1 利用隨機(jī)函數(shù)打亂具有25個(gè)主要變量的數(shù)據(jù)順序；Step 2 將數(shù)據(jù)按均勻分布隨機(jī)分成M份，取M=10；Step 3 挑選其中9份數(shù)據(jù)作為訓(xùn)練集，剩下的1份作為測(cè)試集；Step 4 依次重復(fù)進(jìn)行10次，根據(jù)訓(xùn)練集訓(xùn)練出模型；Step 5 將模型放到測(cè)試集上，得到10次預(yù)測(cè)結(jié)果；Step 6 挑選預(yù)測(cè)精度最高的1次，作為辛烷值損失預(yù)測(cè)模型.

3.3 模型求解

根據(jù)上述結(jié)果，利用十折交叉驗(yàn)證法確定的訓(xùn)練集和測(cè)試集建立了多元線性回歸模型，通過(guò)MATLAB編程，可得到結(jié)果(見圖2).

同時(shí)將十折交叉驗(yàn)證法隨機(jī)生成的10次訓(xùn)練方案所得到多元線性回歸分析模型的回歸系數(shù)、用于檢驗(yàn)回歸模型的統(tǒng)計(jì)量以及測(cè)試數(shù)據(jù)的平均誤差結(jié)果輸出(見表5).

由表5可以看出，第9次訓(xùn)練方案所得測(cè)試數(shù)據(jù)的平均誤差最小，為0.144 4.因此，我們將第9次訓(xùn)練方案得到的多元線性回歸分析模型作為最終所建立的辛烷值損失預(yù)測(cè)模型，即

y=-6.319 9-0.000 3x1+0.015 5x2+

0.006 9x3-0.014 5x4+0.000 9x5-

0.014 4x6+0.018 3x7-0.003 1x8+

0.008 8x9+0.031 3x10+0.206 4x11-

0.000 1x12+0.000 0x13- 0.032 0x14-

9.993 1x15-0.000 1x16+0.000 0x17+

0.000 9x18-0.046 4x19+0.000 0x20-

0.008 6x21-0.047 0x22+0.000 5x23+

1.498 4x24+0.000 5x25.

(17)

圖2 多元線性回歸模型結(jié)果

表5 10次訓(xùn)練方案的輸出結(jié)果

續(xù)表

3.4 模型檢驗(yàn)

求解出回歸模型的回歸系數(shù)后，需要對(duì)回歸模型和回歸系數(shù)分別進(jìn)行假設(shè)檢驗(yàn)，其檢驗(yàn)過(guò)程如下.

H0:βj=0,j=1,2,…,m，

(18)

其中m=25,n=325.由于

(19)

當(dāng)H0成立時(shí)，有

(20)

在顯著性水平α下，若

F1-α/2(m,n-m-1)

(21)

則接受H0，否則拒絕.

利用MATLAB編程可求得統(tǒng)計(jì)量為F=2.458 3，查表得上α/2分位數(shù)為

2.458 3=F>Fα/2(m,n-m-1)=1.679 9.

因此拒絕 (18) 式的原假設(shè)，模型整體上通過(guò)檢驗(yàn).

3.4.2對(duì)回歸系數(shù)的假設(shè)檢驗(yàn) 當(dāng) (18) 式的原假設(shè)H0被拒絕時(shí)，βj不全為0，但是不排除其中若干個(gè)等于0.因此，應(yīng)進(jìn)一步作(m+1)個(gè)假設(shè)：

(22)

(23)

式中cjj為(XTX)-1中的第j行第j列元素.對(duì)給定的α，若

|tj|

(24)

利用MATLAB編程可求得統(tǒng)計(jì)量t的26個(gè)值分別為-2.429、-1.470、1.326、1.402、-1.414、1.045、-2.405、3.437、-1.890、2.156、3.321、1.369、-0.760、0.975、-0.757、-0.567、-0.256、0.490、1.198、-1.694、0.147、-0.597、0.311、0.564、0.955、0.502.

查表得上α/2分位數(shù)為tα/2(299)=1.968 9.因此，對(duì)于統(tǒng)計(jì)量|tj|

4 結(jié)語(yǔ)

本文給出了基于多元線性回歸分析的汽油辛烷值(RON)損失的預(yù)測(cè)模型.基于催化裂化汽油精制裝置所采集的325個(gè)數(shù)據(jù)樣本，首先，利用拉依達(dá)準(zhǔn)則(3σ準(zhǔn)則)去除異常值等方法對(duì)原始數(shù)據(jù)樣本進(jìn)行預(yù)處理；其次，為了工程應(yīng)用方便，根據(jù)系統(tǒng)聚類法中的R型聚類法，利用相關(guān)系數(shù)法和最大系數(shù)法分別確定指標(biāo)變量之間的相似性并選取具有代表性的變量，得到降維后的25個(gè)主要變量；再次，以辛烷值損失為因變量，建立多元線性回歸分析模型進(jìn)行預(yù)測(cè).最后，采用十折交叉驗(yàn)證法對(duì)模型進(jìn)行驗(yàn)證，得到了辛烷值損失預(yù)測(cè)模型.同時(shí)，對(duì)回歸模型和回歸系數(shù)分別進(jìn)行了假設(shè)檢驗(yàn)，驗(yàn)證了所構(gòu)建的辛烷值損失預(yù)測(cè)模型的合理性.