基于模式識別和譜圖映射的通用油品調(diào)合模型

王 瑩，杜中元，李澤飛

(中國科學院 自動化研究所，北京 100190)

傳統(tǒng)油品調(diào)合優(yōu)化模型一般直接從油品的各個性質(zhì)出發(fā)建立質(zhì)量約束，無論是汽油、柴油還是原油，需要分析的油品性質(zhì)指標都很多，而且某些性質(zhì)在成品油與組分油之間是較為復雜的非線性關(guān)系，并沒有一個統(tǒng)一的模型能用來描述，不同煉油廠必須根據(jù)自己的生產(chǎn)情況去尋找適用的模型[1]。相關(guān)文獻中有很多關(guān)于油品性質(zhì)調(diào)合關(guān)系的模型，如柴油凝點性質(zhì)，有指數(shù)模型、換算因子法、換算指數(shù)法等[1-2]；汽油辛烷值性質(zhì)，有Stewart法、調(diào)合因素法、指數(shù)對數(shù)模型、ETHYL RT-70法等[1-3]。在實際應用中，當生產(chǎn)條件變化導致組分油波動時，原有的模型很可能不再具有適用性，需要再采集新的樣本去擬合新的模型，這就造成了油品調(diào)合任務的繁雜和高成本。

針對以上問題，筆者提出了一種通用的基于模式識別和譜圖映射的油品調(diào)合優(yōu)化模型，將不統(tǒng)一、繁雜的非線性性質(zhì)卡邊約束轉(zhuǎn)化為訓練集特征空間上成品油與組分油得分向量間的線性調(diào)合約束，這種建模思想適用于包括汽油、柴油、原油等不同油品的調(diào)合，適用于不同的煉油廠，也適用于同一個煉油廠生產(chǎn)條件的波動。

油品的性質(zhì)是由各類烴化合物的各種X—H結(jié)構(gòu)基團的特征和含量決定的，而油品的近紅外譜圖是一種能直接反映油品烴類化合物的組成、結(jié)構(gòu)信息的數(shù)據(jù)[4-5]。油品調(diào)合是組分油的分子擴散達到完全均勻的混合為止[1]，沒有化學反應的產(chǎn)生。成品油中各種C—H、N—H和O—H鍵結(jié)構(gòu)基團的含量與組分油中各種C—H、N—H和O—H鍵結(jié)構(gòu)基團含量的線性加和基本一致，所以成品油的近紅外光譜圖的吸收度也是組分油譜圖吸收度的線性加和。用主成分分析法對油品的譜圖提取出特征向量——得分向量，經(jīng)過推導也可以得出成品油與組分油的得分向量間存在線性加和關(guān)系。

基于上述理論依據(jù)，筆者提出基于模式識別和譜圖映射的油品調(diào)合優(yōu)化模型，訓練集最大限度的覆蓋煉油廠正常生產(chǎn)條件下的所有工況，包含各種代表性組分油和成品油樣本的近紅外譜圖，用主成分分析法運算后在特征空間形成分類模型；將參與調(diào)合的組分油譜圖投影到特征空間，并用樸素貝葉斯分類器進行分類判別到訓練集的某類樣本中；在成品油類中用核密度估計方法量化樣本的分布并確定一個約束區(qū)域；最后利用組分油和成品油得分向量之間的線性加和關(guān)系建立調(diào)合油品的質(zhì)量約束。

1 成品油與組分油得分向量之間的線性加和關(guān)系

成品油的近紅外光譜圖與組分油近紅外譜圖的線性加和關(guān)系用數(shù)學表達式表示，如式(1)所示。

(1)

對訓練集樣本以及上述成品油和組分油樣本譜圖進行主成分分析法運算，可以得到式(2)。

(2)

將式(2)代入式(1)，兩邊同時乘以L并化簡得到式(3)。

(3)

選擇前f個主成分上的得分，用式(4)表示。

(4)

從式(4)可以看出，成品油的得分向量與組分油的得分向量是線性加和關(guān)系。

2 訓練集樣本在特征空間的分類和調(diào)合組分油樣本的分類判別

對訓練集樣本作主成分分析，確定能涵蓋原始樣本信息98%左右的主成分個數(shù)f，這樣不僅能得到反映每個油品特征的得分向量，而且在主成分構(gòu)成的特征平面或特征空間上不同油品落入不同的區(qū)域，不同油品各自分開，屬于不同的類別[4-8]。

采集參與調(diào)合的組分油近紅外譜圖，并投影到訓練集樣本構(gòu)成的特征空間上，采用樸素貝葉斯分類器[9-10]對組分油樣本進行分類判別到訓練集的某類樣本中。生產(chǎn)條件的波動會導致組分油落在所屬類區(qū)域的不同位置，得到不同的得分向量，將對應的得分向量帶入到下文中的式(8)用于建立調(diào)合油品質(zhì)量約束。

用D=(D1,D2,…,Df)表示樣本屬性集(得分向量)，C表示類變量(油品類別)，樸素貝葉斯分類器用先驗概率P(D)、P(C)和類條件概率P(D|C)通過式(5)表示后驗概率，并將樣本判別到后驗概率最大的類別中，對于確定的屬性集分母P(D)是相同的，只需判別到分子取最大值的類別中。

(5)

樸素貝葉斯分類器假設樣本屬性條件獨立，所以對于給定的類別c，特定的屬性值d=(d1,d2,…,df)，類條件概率可以用式(6)表示。

(6)

類的先驗概率可以通過訓練集的各類樣本出現(xiàn)的次數(shù)來估計(c類先驗概率=c類樣本的數(shù)量/樣本總數(shù))。為了估計屬性的分布參數(shù)，可以假設訓練集數(shù)據(jù)滿足高斯分布或核密度分布，用樣本值估計出屬性分布函數(shù)后計算某個屬性值上的概率。

3 核密度估計方法確定調(diào)合油品的落入?yún)^(qū)域

要使得某幾種組分油按一定比例調(diào)合出滿足質(zhì)量指標要求的成品油，調(diào)合油品的得分向量也必須落入對應的成品油所屬區(qū)域的類中。因此，在成品油類中確定一個子區(qū)域作為調(diào)合油品得分向量必須落入的區(qū)域，以此建立成品油的質(zhì)量約束。

核密度估計[9-12](Kernel density estimation)又名Parzen窗(Parzen window)，用于所研究的樣本概率密度函數(shù)形式未知的場合，直接利用訓練數(shù)據(jù)對概率密度進行估計。設s1,s2…sn為一元獨立同分布F的n個樣本點，F(xiàn)的概率密度函數(shù)用核密度估計的公式如下：

(7)

K(x)為核函數(shù)，h>0為一個平滑參數(shù)，稱作帶寬(Bandwidth)，也叫窗口。由于高斯內(nèi)核方便的數(shù)學性質(zhì)，經(jīng)常使用K(x)=φ(x)，φ(x)為高斯核函數(shù)。

筆者選用的核函數(shù)也是高斯核函數(shù)，窗寬采用可使平均積分平方誤差最小的最優(yōu)窗寬，在各主成分上分別進行一維的核密度函數(shù)估計。確定調(diào)合油品落入?yún)^(qū)域的方法如下：

分別在第1到第f主成分軸上按2個步驟確定區(qū)間范圍：(1)根據(jù)樣本的核密度函數(shù)確定樣本最集中的點，即核密度函數(shù)最大值對應的點a；(2)設定概率水平β，在a兩側(cè)確定上、下分位數(shù)λ1、λ2，λ1

4 譜圖調(diào)合優(yōu)化模型

假設調(diào)合一種成品油(調(diào)合多種成品油時對每種成品油用相同的方法分別建立質(zhì)量約束)，在組分油和成品油樣本確定的特征空間上，根據(jù)組分油和成品油的得分向量以及按第3節(jié)中介紹的方法確定的約束區(qū)域，建立質(zhì)量約束表達式如式(8)。

(8)

同時構(gòu)建模型中物料平衡、產(chǎn)量需求等其他約束條件，并且與目標函數(shù)一起構(gòu)成了整個優(yōu)化模型。借助于優(yōu)化軟件即可求解。

4.1 主成分個數(shù)對求解結(jié)果的影響

主成分個數(shù)f的大小影響得分向量代表原始譜圖信息的能力，進而影響不同油品的分類效果，如在二維特征平面有交疊的不同油品可能在三維特征空間就是完全分開的。因此，應該選擇較大的f使得得分向量包含原始數(shù)據(jù)的信息多，不同油品在特征空間中的區(qū)分效果好；但f也不能太大，以免包含噪音信息。f的大小同時影響著解的可行域的大小。f偏小時，解的自由度較大，得出的配方可能使得調(diào)合油品質(zhì)量不合格或質(zhì)量過剩太多；當f偏大時，解的自由度較小，可能導致所求問題無解。

4.2 約束區(qū)域?qū)η蠼饨Y(jié)果的影響

當主成分個數(shù)確定后，約束區(qū)域的大小和形狀都會對求解結(jié)果產(chǎn)生影響。概率水平β取值的大小直接影響約束區(qū)域的大小。當β取值太小時，可行解的范圍就會縮小，優(yōu)化問題可能會無解；β取值太大時，約束區(qū)域較大，邊界部分對應的樣本點可能較稀疏，所求出的優(yōu)化配方調(diào)合出的成品油可能是不合格的或質(zhì)量過剩較多。

5 譜圖調(diào)合優(yōu)化模型應用實例

采集的某煉油廠的柴油數(shù)據(jù)，訓練集數(shù)據(jù)是包含0#柴油(D-0#，國Ⅳ標準)、-35#柴油(D-35#，國Ⅳ標準)、一套常壓蒸餾裝置側(cè)一線(1CDUC1)、一套常壓蒸餾裝置側(cè)二線(1CDUC2)、二套重油FCC裝置輕柴油(2RFCCD)、加氫裂化柴油(HCD)、加氫精制柴油(HTD)、加氫裂化煤油(HCK)、加氫精制煤油(HTK)9種油品312個近紅外數(shù)據(jù)。

有2個-35#成品油樣本是由5個組分油：一套常壓蒸餾裝置側(cè)一線、二套重油FCC裝置輕柴油、加氫精制煤油、加氫裂化柴油、加氫裂化煤油，分別按φ(1CDUC1)=25%、φ(2RFCCD)=10%、φ(HTK)=20%、φ(HCD)=15%、φ(HCK)=30%，以及φ(1CDUC1)=25%、φ(2RFCCD)=5%、φ(HTK)=15%、φ(HCD)=10%、φ(HCK)=45%兩組體積分數(shù)調(diào)合出來的。采集這5個組分油和2個調(diào)合成品油樣品的近紅外譜圖并投影到訓練集樣本確定的特征空間中，得到各自的得分向量。將5個組分油的譜圖按照這兩組比例線性加和得到2個調(diào)合譜圖，也投影到特征空間中。圖1顯示了訓練集樣本在前2個主成分軸確定的特征平面上的分類，其中7個空心圓圈()為5個組分油和2個調(diào)合成品油的譜圖對應的投影，2個星()為2個線性加和的調(diào)合譜圖對應的投影。由圖1可以看出，調(diào)合成品油的譜圖與線性加和的調(diào)合譜圖的得分向量偏差很小。

5.1 調(diào)合問題

用上述5個組分油調(diào)合-35#柴油，組分油得分向量、產(chǎn)量、價格、密度和調(diào)合成品油的價格列于表1。當主成分個數(shù)為3時，能涵蓋訓練集97.9%的信息，因此在問題求解中選擇的主成分個數(shù)為3，得分向量只列出了前三維的數(shù)據(jù)。

調(diào)合要求為：調(diào)合出滿足國Ⅳ指標要求的-35#柴油(訓練集樣本的-35#柴油滿足的是國Ⅳ指標)，質(zhì)量過剩盡量小，組分油用量不超過產(chǎn)量，在滿足這些條件的同時，要求收益最大。

調(diào)合問題中，二套重油FCC裝置輕柴油(2RFCCD)、加氫裂化柴油(HCD)是凝固點較高的組分油，產(chǎn)量較大，價格較便宜，追求收益大就會多消耗這2種組分油，但用量太多會使得調(diào)合出的成品油不合格，因此用量應適當。

圖1 訓練集樣本在特征平面上的分類以及成品油譜圖和調(diào)合譜圖在上面的投影Fig.1 Classification of training samples on feature plane and projection of product oil spectra and blended spectra on itPC1—The first principal components;PC2—The second principal components; Projection of the components and blended product oils spectra; Projection of blended spectra;(1) Projection of blended product oil spectrum and blended spectrum with φ(1CDUC1)=25%, φ(2RFCCD)=10%, φ(HTK)=20%, φ(HCD)=15%, φ(HCK)=30%; (2) Projection of blended product oil spectrum and blended spectrum with φ(1CDUC1)=25%, φ(2RFCCD)=5%, φ(HTK)=15%, φ(HCD)=10%, φ(HCK)=45%

Sampletm/tPrice/(RMB·t-1)ρ/(kg·m-3)1CDUC1(0.210,0.147,0.089)700063827882RFCCD(-1.815,-0.072,0.075)87104807883HTK(-0.149,0.059,0.093)78806078808HCD(0.659,0.047,0.128)159805134792HCK(0.169,0.438,0.124)66006412778D-35#6249

5.2 問題求解

選擇的主成分個數(shù)f=3。采用核密度法估計-35#柴油樣本的分布情況來確定約束區(qū)域。由于高凝點組分油的產(chǎn)量較充裕而且價格便宜，所以求解時令概率水平稍微小一些，以免可行解范圍較大時得出的配方不合格，取概率水平β=40%。圖2為用核密度估計方法在第一主成分軸上確定的約束區(qū)間。圖3為在特征空間-35#柴油樣本類中確定的約束區(qū)域。

表2列出了在前3個主成分軸上根據(jù)樣本的分布確定的約束區(qū)間[λ1,λ2]。

按照式(8)建立質(zhì)量約束不等式，同時列出組分油用量約束表達式和目標函數(shù)。模型建立好之后，調(diào)用MATLAB的linprog函數(shù)進行優(yōu)化求解，得出5種組分油的用量配方為m(1CDUC1)=7000 t，m(2RFCCD)=3248 t，m(HTK)=7880 t，m(HCD)=7560 t，m(HCK)=2692 t，對應的5種組分油的體積分數(shù)為φ(1CDUC1)=25%，φ(2RFCCD)=10%，φ(HTK)=28%，φ(HCD)=27%，φ(HCK)=10%。按照該用量配方調(diào)合出的成品油性質(zhì)：凝固點-36℃，冷濾點-31℃，閃點54℃，十六烷值49，十六烷值指數(shù)48、餾程50%為218℃、餾程90%為273℃，密度804 kg/m3，達到-35#柴油國Ⅳ標準指標要求。

圖2 采用核密度估計法在第一主成分軸上確定約束區(qū)間Fig.2 Constraint interval determined by kernel density estimation method on first principal components(a) Kernel density function;(b) Inverse cumulative probability function;a—Point corresponding to the maximum of kernel density function;λ1, λ2—Lower and upper limits of constraint interval;pa, p1, p2—Probability corresponding to a, λ1, λ2;the difference between p1 and p2 is β

圖3 采用核密度估計法在特征空間成品油樣本類中確定的約束區(qū)域Fig.3 Constraint region determined by kernel density estimation method in the product oil samples group in the feature space Points correspond to the product oil samples;Cross corresponds to the point which has maximum value of kernel density function on each of the first three principal components; Cuboid corresponds to the constraint region.

Constraint intervalλ1λ2PC1-0.00570.0646PC20.09870.1494PC3-0.03850.1025

5.3 主成分個數(shù)和約束區(qū)域的選擇對調(diào)合結(jié)果的影響

進一步考察了主成分個數(shù)和約束區(qū)域的選擇對調(diào)合結(jié)果的影響。當選擇主成分個數(shù)分別為2、3、4時，涵蓋的訓練集信息依次為94.5%、97.9%、98.9%，其他條件不變，仍取β=40%。經(jīng)過計算，當主成分個數(shù)f=4時求解出的結(jié)果與f=3時相同。所以，表3列出的是主成分個數(shù)f分別為2和3時的模型優(yōu)化求解的調(diào)合配方。

表3 主成分個數(shù)(f)分別為2和3時用模型進行優(yōu)化求解得到的配方Table 3 Formula optimized by the model choosing two or three principal components

根據(jù)該優(yōu)化配方調(diào)合出的成品-35#柴油性質(zhì)：f=2時，凝固點為-33℃，冷濾點為-28℃；f=3時，凝固點為-36℃，冷濾點為-31℃。當f=2時，得分向量包含的原始樣本信息不夠充分，信息有損失，而且可行解的范圍較大，高凝點組分油產(chǎn)量足夠多而且價錢便宜，所以優(yōu)化求解得到的結(jié)果中用的高凝點組分油量較多，根據(jù)各組分油用量調(diào)合出的成品柴油的凝固點和冷濾點較高，不能滿足-35#柴油國Ⅳ標準指標要求。

仍取f=3，建立質(zhì)量約束區(qū)域時，依次取概率水平β為30%、35%、40%、50%調(diào)用模型算法進行求解。當β為30%時約束區(qū)域太小，模型沒有可行解。表4列出了β為35%、40%、50%的優(yōu)化用量。表5為根據(jù)優(yōu)化配方調(diào)合出的成品柴油的凝固點和冷濾點。

表4 概率水平(β)分別為35%、40%、50%時模型求解的優(yōu)化配方Table 4 Formula optimized by the model choosing the probability (β) equals 35%, 40% or 50%

表5 概率水平(β)分別為35%、40%、50%時根據(jù)模型優(yōu)化求解的配方調(diào)合出的成品柴油的凝固點和冷濾點Table 5 Solidifying point and cold filter plugging point of the blended diesel based on the formula optimized by the model choosing the probability (β) equals 35%, 40% or 50%

由實驗可以看出，約束區(qū)域較小時，模型可能無解；約束區(qū)域較大時，邊界部分對應的樣本點較稀疏，根據(jù)所求出的優(yōu)化配方調(diào)合出的成品油可能不合格，對于其他調(diào)合問題也可能質(zhì)量過剩較多。

6 結(jié) 論

(1)基于模式識別和譜圖映射的油品調(diào)合優(yōu)化模型，以油品的近紅外譜圖為基礎(chǔ)，在特征空間上根據(jù)成品油和組分油得分向量間的線性調(diào)合關(guān)系來建立模型，解決了傳統(tǒng)上從油品性質(zhì)出發(fā)建立調(diào)合模型時面臨的非線性、不統(tǒng)一、復雜等問題。筆者所提出的建模思想能適應生產(chǎn)條件波動導致的組分油性質(zhì)波動，適合不同的煉油廠，也適合柴油、汽油、原油等不同油品的調(diào)合。

(2)在模型中影響求解結(jié)果的因素和參數(shù)有主成分個數(shù)、譜圖投影的特征空間、確定約束區(qū)域的概率水平、每個主成分軸上的約束區(qū)間；這些影響因素反映了煉油廠正常生產(chǎn)工況下組分油和成品油性質(zhì)波動的特征，反映了組分油和成品油的組成結(jié)構(gòu)，針對具體調(diào)合問題應該合理確定。

符號說明：

A=(a1,a2,…,af)——在特征空間的成品油類中確定的一個中心點；

ai——第i個主成分軸上核密度函數(shù)最大值對應的點；

a——某個主成分軸上核密度函數(shù)最大值對應的點；

b*——成品油的得分向量，k維向量；

b——成品油取f個主成分的得分向量，f維向量；

C——類變量(油品類別)；

c——特定的某個類別；

D——樣本屬性集(得分向量)；

d——特定的某個屬性值；

ex——殘差光譜；

et*,et——得分向量偏差；

f——主成分個數(shù)；

f*(x)——密度估計函數(shù)；

h——核密度函數(shù)的窗寬；

K(x)——核函數(shù)；

k——光譜數(shù)據(jù)的維數(shù)；

L=(l1,l2,lf,…,lk)——訓練集樣本進行主成分運算后在特征空間的基；

LT——L的轉(zhuǎn)置；

M——組分油個數(shù)；

m——用量(質(zhì)量)，t；

P(C)——類別的先驗概率；

P(D)——屬性的先驗概率；

P(C|D)——類的后驗概率；

P(D|C)——類條件概率；

s1,s2…sn——用于估計核密度函數(shù)的樣本點，n為樣本點個數(shù)；

ti——組分油i取f個主成分的得分向量，f維向量；

x——成品樣光譜或中心化后的光譜，k維向量；

yi——組分油i的光譜或中心化后的光譜(減去平均光譜后的光譜)，k維向量；

β——設定的概率水平，樣本在約束區(qū)間[λ1,λ2]分布的概率為β；

φ(x)——高斯核函數(shù)；

φ——體積分數(shù)，%；

φi，φj——組分油i、j調(diào)入的體積分數(shù)，%；

λ1,λ2——某個主成分軸上確定的約束區(qū)間的上、下限；

ρ——密度，kg/m3；

[1] 蔡智,黃維秋,李偉民,等.油品調(diào)合技術(shù)[M].北京:中國石化出版社,2005.

[2] 謝可堃, 王志剛, 張曉光, 等. 汽柴油調(diào)合常用質(zhì)量指標的設計計算方法[J].煉油技術(shù)與工程, 2009, 39(9): 57-60.(XIE Kekun, WANG Zhigang, ZHANG Xiaoguang, et al. Design and calculation method of common quality index for gasoline and diesel blending[J].Petroleum Refinery Engineering, 2009, 39(9): 57-60.)

[3] 王偉, 李澤飛, 黃燕. 基于油品性質(zhì)的汽油調(diào)合辛烷值模型的選取[J].石油學報(石油加工), 2006, 22(6): 39-44．(WANG Wei, LI Zefei, HUANG Yan. Choice of octane number models for gasoline based on feedstock quality[J].Acta Petrolei Sinica (Petroleum Processing Section), 2006, 22(6): 39-44.)

[4] 陸婉珍, 袁洪福, 徐廣通, 等. 現(xiàn)代近紅外光譜分析技術(shù)[M].北京: 中國石化出版社, 2000.

[5] 馮新滬, 史永剛. 近紅外光譜及其在石油產(chǎn)品分析中的應用[M].北京: 中國石化出版社, 2002.

[6] 禇小立, 袁洪福, 陸婉珍. 光譜結(jié)合主成分分析和模糊聚類方法的樣品聚類與識別[J].分析化學研究報告, 2000, 28(4): 421-427.(CHU Xiaoli, YUAN Hongfu, LU Wanzhen. Sample clustering and identification based on spectral principal component analysis and fuzzy clustering method[J].Chinese Journal of Analytical Chemistry,2000, 28(4): 421-427.)

[7] 劉倩, 孫培艷, 高振會, 等. 衰減全反射傅里葉變換紅外光譜技術(shù)結(jié)合模式識別進行油品鑒別[J].光譜學與光譜分析, 2010, 30(3): 663-665.(LIU Qian, SUN Peiyan, GAO Zhenhui, et al. Attenuated total reflectance FTIR spectroscopy and pattern recognition for oil identification [J].Spectroscopy and Spectral Analysis, 2010, 30(3): 663-665.)

[8] 劉穎榮, 許育鵬, 楊海鷹, 等. 汽油樣品類型的模式識別研究與應用[J].色譜，2004, 22(5): 482-485.(LIU Yingrong, XU Yupeng, YANG Haiying, et al. Investigation and application of gasoline sample identity technique[J].Chinese Journal of Chromatography, 2004, 22(5): 482-485.)

[9] SERGIOS T H, KONSTANTINOS K. 模式識別[M].李晶皎譯. 第4版. 北京：電子工業(yè)出版社, 2016.

[10] 謝中華. MATLAB統(tǒng)計分析與應用: 40個案例分析[M].北京: 北京航空航天大學出版社, 2010.

[11] 徐東斌, 黃磊, 劉昌平. 自適應核密度估計運動檢測方法[J].自動化學報，2009, 35(4): 379-385.(XU Dongbin, HUANG Lei, LIU Changping. Adaptive kernel density estimation for motion detection[J].Acta Automatica Sinica, 2009, 35(4): 379-385.)

[12] 淦文燕, 李德毅. 基于核密度估計的層次聚類算法[J].系統(tǒng)仿真學報, 2004, 16(2): 302-309.(GAN Wenyan, LI Deyi. Hierarchical clustering based on kernel density estimation[J].Journal of System Simulation, 2004, 16(2): 302-309.)