基于模式識別和譜圖映射的通用油品調合模型

王 瑩，杜中元，李澤飛

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

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

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

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

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

1 成品油與組分油得分向量之間的線性加和關系

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

(1)

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

(2)

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

(3)

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

(4)

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

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

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

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

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

(5)

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

(6)

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

3 核密度估計方法確定調合油品的落入區域

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

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

(7)

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

筆者選用的核函數也是高斯核函數，窗寬采用可使平均積分平方誤差最小的最優窗寬，在各主成分上分別進行一維的核密度函數估計。確定調合油品落入區域的方法如下：

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

4 譜圖調合優化模型

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

(8)

同時構建模型中物料平衡、產量需求等其他約束條件，并且與目標函數一起構成了整個優化模型。借助于優化軟件即可求解。

4.1 主成分個數對求解結果的影響

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

4.2 約束區域對求解結果的影響

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

5 譜圖調合優化模型應用實例

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

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

5.1 調合問題

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

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

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

圖1 訓練集樣本在特征平面上的分類以及成品油譜圖和調合譜圖在上面的投影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 問題求解

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

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

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

圖2 采用核密度估計法在第一主成分軸上確定約束區間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 采用核密度估計法在特征空間成品油樣本類中確定的約束區域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 主成分個數和約束區域的選擇對調合結果的影響

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

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

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

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

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

表5 概率水平(β)分別為35%、40%、50%時根據模型優化求解的配方調合出的成品柴油的凝固點和冷濾點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%

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

6 結 論

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

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

符號說明：

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

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

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

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

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

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

c——特定的某個類別；

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

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

ex——殘差光譜；

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

f——主成分個數；

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

h——核密度函數的窗寬；

K(x)——核函數；

k——光譜數據的維數；

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

LT——L的轉置；

M——組分油個數；

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

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

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

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

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

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

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

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

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

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

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

φ——體積分數，%；

φi，φj——組分油i、j調入的體積分數，%；

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

ρ——密度，kg/m3；

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

[2] 謝可堃, 王志剛, 張曉光, 等. 汽柴油調合常用質量指標的設計計算方法[J].煉油技術與工程, 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] 王偉, 李澤飛, 黃燕. 基于油品性質的汽油調合辛烷值模型的選取[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] 陸婉珍, 袁洪福, 徐廣通, 等. 現代近紅外光譜分析技術[M].北京: 中國石化出版社, 2000.

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

[6] 禇小立, 袁洪福, 陸婉珍. 光譜結合主成分分析和模糊聚類方法的樣品聚類與識別[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] 劉倩, 孫培艷, 高振會, 等. 衰減全反射傅里葉變換紅外光譜技術結合模式識別進行油品鑒別[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版. 北京：電子工業出版社, 2016.

[10] 謝中華. MATLAB統計分析與應用: 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].系統仿真學報, 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.)