基于非負(fù)矩陣分解的卷煙配方維護方法

王 楠，雒興剛*，張忠良，于 洋，唐加福，湯建國

1. 杭州電子科技大學(xué)管理學(xué)院，杭州市杭州經(jīng)濟開發(fā)區(qū)白楊街道2 號大街1158 號 310018

2. 東北大學(xué)信息科學(xué)與工程學(xué)院，沈陽市和平區(qū)文化路3 號巷11 號 110819

3. 云南中煙工業(yè)有限責(zé)任公司技術(shù)中心，云南省昆明市五華區(qū)紅錦路367 號 650231

卷煙配方是由多種不同的單等級片煙按照一定比例配伍而成，在卷煙生產(chǎn)過程中，卷煙配方維護是維持卷煙品牌質(zhì)量穩(wěn)定的重要內(nèi)容[1]。傳統(tǒng)的卷煙配方維護是一項復(fù)雜的工作[2]，需要感官評吸和測量儀器綜合評定替換后的卷煙配方的感官和煙氣指標(biāo)是否滿足要求。這種依靠專家評吸和分析檢測的方式效率低、耗時長、成本高，因此，開發(fā)卷煙配方智能化維護技術(shù)和手段顯得必要且迫切[3]。

近年來，隨著機器學(xué)習(xí)、數(shù)據(jù)挖掘及人工智能在各領(lǐng)域的廣泛應(yīng)用與發(fā)展，數(shù)據(jù)挖掘等智能技術(shù)在卷煙產(chǎn)品開發(fā)與生產(chǎn)過程中的應(yīng)用也有一些報道[4-5]。王曉輝等[6]為了克服基于煙葉感官分組方法易受人為主觀影響的缺點，提出了根據(jù)煙葉的化學(xué)指標(biāo)進行聚類分析的方法對煙葉進行分類。賴燕華等[7]采用相似度分析、主成分分析和聚類分析對卷煙質(zhì)量進行綜合分析，為卷煙產(chǎn)品的質(zhì)量控制提供有效的評價方法。此外，還有研究采用數(shù)據(jù)挖掘技術(shù)構(gòu)建煙葉化學(xué)成分與感官和煙氣質(zhì)量之間的關(guān)系模型。Feng 等[8]將遺傳算法、神經(jīng)網(wǎng)絡(luò)、支持向量機（SVM）和模糊集方法與專家知識結(jié)合起來建立了產(chǎn)品配方設(shè)計智能化系統(tǒng)。郭東鋒等[9]分析了化學(xué)成分與主流煙氣指標(biāo)之間的相關(guān)關(guān)系。張圣男[10]采用遺傳算法優(yōu)化支持向量機的參數(shù)建立了基于支持向量機的卷煙煙氣指標(biāo)預(yù)測模型。過偉民等[11]分析了煙葉不同理化指標(biāo)對感官指標(biāo)的影響情況。王濤[12]提出了用回歸函數(shù)估計SVM 方法進行建模來解決感官評估多類劃分問題。張忠良等[13]將基于代價敏感的反饋神經(jīng)網(wǎng)絡(luò)應(yīng)用于卷煙感官評估中，可有效識別卷煙感官評估中的高代價樣本。雖然已有研究對卷煙配方智能化維護進行了一些探索，但這些研究大多根據(jù)單等級片煙特征進行聚類分析或者根據(jù)卷煙的物理化學(xué)指標(biāo)與感官質(zhì)量及煙氣指標(biāo)之間的映射關(guān)系建模，缺乏考慮單等級片煙配伍性問題。最近，雒興剛等[14]運用基于關(guān)聯(lián)規(guī)則挖掘技術(shù)試圖獲得單等級片煙之間的配伍信息，然而該方法局限于配方中部分非主干單等級片煙之間的替換，無法維護配方中任意單等級片煙缺失的情況。

基于此，本文中提出了一種基于非負(fù)矩陣分解（Non-negative Matrix Factorization，NMF）技術(shù)[15]的卷煙配方維護方法。與其他矩陣分解技術(shù)（如奇異值分解SVD，主成分分析PCA）相比，NMF 的分解結(jié)果有直觀的語義解釋，且配方可以看作是單等級片煙的線性組合。而且NMF 具有基于部分表示的特點，可以結(jié)合不同的煙草數(shù)據(jù)特征來表示數(shù)據(jù)，發(fā)現(xiàn)配方數(shù)據(jù)中單等級片煙間的配伍關(guān)系。另外，NMF 方法可以對大量的數(shù)據(jù)進行壓縮和概括，節(jié)省存儲空間，可提高數(shù)據(jù)處理速度。因此，本文中所提出的方法不僅可以與卷煙配方數(shù)據(jù)相結(jié)合，挖掘出單等級片煙之間的配伍規(guī)則，而且對于配方中任何單等級片煙缺失的情況都能進行有效處理，可為卷煙配方設(shè)計人員推薦合理的卷煙配方方案。

1 材料與方法

1.1 數(shù)據(jù)準(zhǔn)備

選擇某煙草企業(yè)3 年（2010 年1 月—2012 年12 月）的數(shù)據(jù)進行研究。其中包括24 個卷煙品牌，涉及共1 534 個配方及其包含的1 220 種單等級片煙。單等級片煙包括國外102 份（6 個不同國家），國內(nèi)1 018 份（14 個不同省份），13 個品種，497個不同等級，6 種不同顏色。其中，包括24 種再造煙葉，5 種出口片煙及1 191 種國內(nèi)片煙。

單等級片煙數(shù)據(jù)文件包含1 220 種單等級片煙及其28 個屬性信息。配方數(shù)據(jù)文件包含1534 種配方（編號FBI0001～FBI1534），列出了每種配方包含的 單 等 級 片 煙（AD0001～AD0716，DLC0001～DLC2072，不連續(xù)）共1220 個，及其所屬的品牌。其中任意配方基本保持由10 至37 份單等級片煙經(jīng)過一定的比例配制而成。

將配方數(shù)據(jù)轉(zhuǎn)換成二進制配方矩陣Y（n×m），其中，n 表示配方的數(shù)量，m 表示單等級片煙的數(shù)量。如果配方i 含有單等級片煙j，則Yij為1，否則為0。

1.2 方法

1.2.1 非負(fù)矩陣分解

NMF 方法是Lee 和Seung 在1999 年提出的一種新的矩陣分解技術(shù)[15]，該方法已在圖像處理與識別、文本挖掘、語音處理、故障檢測與診斷等領(lǐng)域得到了廣泛的應(yīng)用[16-26]。

非負(fù)矩陣分解（NMF）是所有數(shù)據(jù)都是非負(fù)數(shù)的矩陣分解方法。對于任意給定的一個非負(fù)矩陣Y，NMF 算法能夠?qū)ふ业揭粋€非負(fù)矩陣W 和一個非負(fù)矩陣H，使得滿足：

其中：Y 是一個（n×m）卷煙配方矩陣；W 是一個包含k 個潛在特征的（n×k）基矩陣；H 也是一個包含k 個潛在特征的（k×m）矩陣。參數(shù)k 是分解的等級并確定矩陣分解的潛在特征數(shù)量。W 和H 的最優(yōu)值可以通過最小化Y 和WH 之間的平方誤差來確定[27-28]。

1.2.2 基于NMF 的卷煙配方維護模型

圖1 基于NMF 的卷煙配方維護模型框架Fig.1 Framework of NMF-based cigarette blend maintenance model

為了確定模型系數(shù)矩陣M，應(yīng)用NMF 方法分解卷煙配方矩陣，如圖2 所示。NMF 方法將識別給定配方數(shù)據(jù)的信息，并用k 個潛在特征總結(jié)這些信息，所以k 的值對于數(shù)據(jù)的表示至關(guān)重要。當(dāng)k的值太小時，數(shù)據(jù)的表達就會不足；當(dāng)k 的值太大時，就會過度擬合數(shù)據(jù)。k 的最優(yōu)值將在下文驗證期間確定。

k 值確定后，配方數(shù)據(jù)矩陣的分解為：

其中Yold_formulas和Xold_formulas是配方的二進制表示，采用最小二乘方法估計M[29]：

1.2.3 模型參數(shù)優(yōu)化及測試方法

圖2 NMF 分解卷煙配方矩陣Fig.2 NMF factorized cigarette blend matrix

將配方數(shù)據(jù)分為訓(xùn)練集+驗證集和測試數(shù)據(jù)集，其中，訓(xùn)練集+驗證集共占配方數(shù)據(jù)的90%，測試集占10%。在訓(xùn)練和驗證模型時采用基于10 折交叉驗證的方法，將訓(xùn)練集+驗證集隨機平均分為10 份，在10 次循環(huán)中，每次選取9 份作為訓(xùn)練集用于訓(xùn)練模型，1 份作為驗證集用于調(diào)整模型參數(shù)，最后測試集用于測試模型性能[30]。

模型完成配方的能力是通過驗證模型是否能夠檢索從配方中消除的單等級片煙來完成的。對于驗證集中的每個配方i，隨機選擇一個單等級片煙j 并從配方中消除，這個單等級片煙在兩種配方中是可以相同的，也可以是不同的。在實際驗證集矩陣中，將消除的單等級片煙j 對應(yīng)的“1”替換為“0”，即，配方完成驗證后將恢復(fù)到它們的初始狀態(tài)

在驗證模型時，模型為每個配方返回1 220 個單等級片煙預(yù)測值的有序列表。接著，確定被消除的單等級片煙在輸出單等級片煙有序列表中的排列等級。如果消除的單等級片煙在列表中具有最高預(yù)測值，則其排名是1（最高等級），驗證過程如圖3 所示。

圖3 驗證過程Fig.3 Validation process

通過反復(fù)試驗，確定矩陣分解參數(shù)k 的范圍取為30～40，根據(jù)式（2）訓(xùn)練集配方的二進制矩陣Ytrain被分解成W 和H：

Ytrain和Xtrain是訓(xùn)練集中的二進制配方。當(dāng)Xold_formulas=Xtrain時，由式（3）估計出模型系數(shù)矩陣，所以驗證集配方的預(yù)測值為：

其中Xvalidate是驗證集中的二進制配方矩陣。Yvalidate，pred根據(jù)圖3 驗證過程確定的Xvalidate中每個消除單等級片煙在推薦序列表中的等級。針對每個驗證集，分別計算所選范圍中的每個k 值對應(yīng)的平均數(shù)等級，選擇平均等級最小的k 值作為當(dāng)前驗證集的最優(yōu)k 值。在所有驗證集的10 個最優(yōu)值k 中，選擇中位數(shù)值作為模型的最優(yōu)k 值。

為了評估模型挖掘單等級片煙間配伍規(guī)則的性能，對模型進行測試，并修改測試集配方y(tǒng)ij=1→yij=0，其中i 為測試的配方之一，j 為配方i 中隨機選擇的單等級片煙。此時的實際訓(xùn)練集為訓(xùn)練集+驗證集，即Yvalidate+train，共1 380 個配方。t 設(shè)定為10，使用驗證期間得到的最優(yōu)k 值，通過NMF 方法分解矩陣Yvalidate+train。使用式（3）計算模型系數(shù)矩陣M，其中Xold_formulas=（Xvalidate+train）。

測試集包含154 個配方，對于測試集中的每個配方，確定消除單等級片煙在推薦列表中的等級，將用于計算如下3 種模型性能度量指標(biāo)[30]：①在推薦單等級片煙的有序列表中消除單等級片煙的平均數(shù)等級。②在推薦單等級片煙的有序列表中消除單等級片煙的中位數(shù)等級。③消除的單等級片煙位于推薦單等級片煙列表的前10 位的百分比。

1.2.4 基于相似度的單等級片煙替換方案

如上所述，在推薦的候選列表S 中選擇最合適的1 種（或1 組）單等級片煙替換缺失的單等級片煙，考慮單等級片煙缺失情況的多樣性，提出了兩種替換方案，即一對一替換和多對多替換兩種情況。本研究中采用歐幾里得距離計算缺失單等級片煙和候選單等級片煙的化學(xué)成分信息（總糖、還原糖、總煙堿、總氮、鉀、氯、蛋白質(zhì)、施木克值）的相似度[7]，距離值越小則意味著相似度越高。具體方案如下：

（1）一對一替換（One Versus One，OVO）。當(dāng)某配方缺失一種單等級片煙，計算缺失的一個單等級片煙與推薦的t 個單等級片煙的相似度并排序，從中選擇相似度最大的一個進行替換。

（2）多 對 多 替 換（Multiple Versus Multiple，MVM）。當(dāng)某配方缺失x（1＜x＜t）種單等級片煙，在推薦的t 個單等級片煙中選擇x 個進行替換。其中多對多替換有3 種方法：

a. 基于一對一的多對多替換（MVM based on OVO，MVM_OO）。每次假設(shè)這x 個缺失配方中只缺少其中一個單等級片煙，在推薦的t 個單等級片煙中選擇等級最高的一個單等級片煙保存，x次替換后，保存的x 個單等級片煙即為替換單等級片煙。

b. 基于最大相似度的多對多替換（MVM based on maximum similarity，MVM_MS）。分別計算每個缺失單等級片煙與t 個推薦的單等級片煙的相似度，在t 個推薦單等級片煙中找出與每個缺失單等級片煙相似度最高的單等級片煙，所組成的單等級片煙組即為缺失單等級片煙的替換單等級片煙。

c. 基于平均相似度的多對多替換（MVM based on average similarity，MVM_AS）。舉例說明：當(dāng)缺失單等級片煙不大于10 個時，假設(shè)某品牌的配方中，缺失DLC01、DLC02、DLC03 3 個單等級片煙時，在推薦的10（t=10）個單等級片煙中進行排列組合C（10，3），選出3 個單等級片煙為1 組，共120 組。首先，分別計算組合的3 個單等級片煙和缺失單等級片煙的平均相似度距離ave_d1，ave_d2，ave_d3。對3 個平均相似度距離求和，產(chǎn)生120 個和值sum_d120-n（n=1,2,…,120），選出和值最小的1組，所對應(yīng)的單等級片煙ZDLCn1，ZDLCn2，ZDLCn3即為替換單等級片煙，如圖4 所示。

圖4 MVM_AS 單等級片煙替換選擇步驟Fig.4 Procedure for MVM_AS substitution of single grade strips

本研究中選用某品牌卷煙的一個配方FBI0282（包含23 種單等級片煙）進行實驗。假設(shè)配方FBI0282 中分別缺失1 種，3 種單等級片煙，產(chǎn)生2 個單等級片煙組：①單等級片煙組1：缺失單等級片煙AD0280。②單等級片煙組2：缺失單等級片煙AD0280、AD0603、DLC0277。使用圖1 框架建立模型，t 設(shè)定為30，這時訓(xùn)練使用的數(shù)據(jù)集包含所有配方（1 534 個）。用以上兩組單等級片煙作為模型實例輸入，分別進行OVO 配方維護和MVM 配方維護實驗，并選擇替換方案。

1.2.5 卷煙配方維護效果評估方法

為了評估基于NMF 的卷煙配方維護模型的效果，基于文獻[10]和[12]的研究成果，采用基于支持向量機的預(yù)測模型建立成品卷煙8 種化學(xué)指標(biāo)（總糖、還原糖、總煙堿、總氮、鉀、氯、蛋白質(zhì)、施木克值）與6 種感官指標(biāo)（光澤、香氣、諧調(diào)、雜氣、刺激性、余味）及3 種煙氣指標(biāo)（一氧化碳、焦油、煙氣煙堿）之間的映射關(guān)系，預(yù)測經(jīng)過維護的卷煙配方的感官和煙氣指標(biāo)，再計算與原配方的感官和煙氣指標(biāo)的差值，差值越小，維護效果越好[9]。采用libsvm 庫的SVR 模型進行預(yù)測，即s=3；核函數(shù)選取徑向基（RBF）函數(shù)，其中懲罰參數(shù)c=2，核參數(shù)g=l。

2 結(jié)果與討論

2.1 模型測試結(jié)果及分析

10 次實驗對154 個配方的測試結(jié)果如表1 所示。由表1 中k 的均值可知，k 的最優(yōu)值約等于37，這意味著需要大約37 個潛在特征來表達配方數(shù)據(jù)。觀察表1 的等級部分，第1 個度量（平均等級）說明平均每一個被消除的單等級片煙可以在推薦的單等級片煙有序列表中的第29 位置（共1 220 個單等級片煙）找到；第2 個度量（中位等級）說明50%的消除單等級片煙的等級小于6。第3 個度量表明被消除的單等級片煙被發(fā)現(xiàn)在前10 等級的配方占測試配方（共154 個單等級片煙）的63.99%，由第5 個度量可知被消除的單等級片煙排在前30的配方占測試配方的百分比達到了82.48%；以上說明模型表現(xiàn)很好。

表1 10 次實驗測試結(jié)果Tab.1 Results of 10 experimental tests

應(yīng)用實驗2 的數(shù)據(jù)，繪制消除單等級片煙的等級統(tǒng)計圖像，見圖5。從圖中可以看到，累計百分比曲線很快地達到一個較高值，然后上升緩慢，最后接近于水平，說明NMF 模型可以更好地了解單等級片煙信息矩陣，推薦的單等級片煙與輸入單等級片煙組具有良好配伍性。從圖5 還可以看出，一些消除的單等級片煙具有非常低的等級，表明該模型不能夠檢索這些被消除的單等級片煙。通過查閱單等級片煙信息表，這些低級等級的單等級片煙往往是相當(dāng)罕見的單等級片煙，通常不與某測試集配方之外的單等級片煙結(jié)合，而只與該測試集配方中的單等級片煙組合。由于這個配方未被用來訓(xùn)練模型，在測試過程中模型將不識別這些單等級片煙組合。當(dāng)配方中很少使用某種單等級片煙時，上述情況出現(xiàn)的概率很高。因此這不是NMF 模型的問題，僅與被消除的單等級片煙和所用測試集配方的特征密切相關(guān)。

圖5 測試數(shù)據(jù)等級分布和累計等級百分比分布Fig.5 Level distribution and cumulative level percentage distribution of test data

2.2 模型系數(shù)矩陣分析

通過仔細(xì)觀察模型系數(shù)矩陣M，M 是一個方陣，其行和列都代表1 220 個單等級片煙，該矩陣給出了每一對單等級片煙的相互影響權(quán)值，這些權(quán)值反映了每對單等級片煙組合的好壞程度。其中一些權(quán)值是負(fù)值，表示一種單等級片煙的存在限制了另一種單等級片煙在推薦單等級片煙列表中的出現(xiàn)，這意味著這兩種單等級片煙不能形成良好的組合。而正值則表示這兩種單等級片煙能很好組合，并且這些單等級片煙之一的存在將增加另一種單等級片煙在推薦單等級片煙列表中排在前十的機會，值越大，推薦該單等級片煙的機會就越高。一些單等級片煙對于其余單等級片煙具有接近零的值，不會對二進制配方向量與矩陣M的乘積造成影響。這說明這些單等級片煙的存在不會影響推薦單等級片煙的列表，并且不會被推薦添加到配方中。

2.3 替換方案分析及選擇

為了選擇最佳的替換方案，預(yù)測新配方的感官和煙氣指標(biāo)，并計算與原配方的感官和煙氣指標(biāo)的平均絕對偏差，針對不同的方案分別用相同的單等級片煙組1 和單等級片煙組2 作為輸入，進行20 次替換的實驗結(jié)果見圖6。從圖6 可以看出，OVO 比3V3 的維護效果好，原因是缺失單等級片煙數(shù)量越多，配方的不穩(wěn)定性越高，替換的效果就會受到影響。3V3 的3 種方法的平均絕對偏差圍繞一個穩(wěn)定值上下波動。進一步計算OVO 和3V3的3 種方法各項指標(biāo)的平均相對偏差，如表2 所示。從表2 可以看出：①OVO 在感官和煙氣指標(biāo)的維護中，相比于3V3 方法，各項指標(biāo)平均相對偏差均較低，各項指標(biāo)綜合平均相對偏差（0.28%）在4 種方法中也達到了最小，達到了非常好的效果，與上述結(jié)論相同。②在3 種3V3 方法中，3V3_AS在煙氣指標(biāo)的維護中達到了非常好的效果，維護的煙氣指標(biāo)的平均相對偏差均達到了最?。欢?V3_MS 和3V3_OO 在感官指標(biāo)的維護中達到了較好的效果，兩者相差不大。③3V3_AS 的綜合平均相對偏差在3 種方法中最小（0.80%）。因為3V3_AS 方法能考慮到各替代單等級片煙間的相互影響，從整體上維護各項指標(biāo)的穩(wěn)定，而3V3_MS 和3V3_OO 僅考慮單一候選單等級片煙與缺失單等級片煙相似性。與其余兩種3V3 方法相比，3V3_AS 方法能更有效維護原品牌綜合指標(biāo)的穩(wěn)定。因此，多對多替換應(yīng)該考慮選擇基于平均相似度的方法（MVM_AS）。

圖6 新配方與原配方的平均絕對偏差Fig.6 Average absolute deviation between new blend and original blend

表2 新配方與原配方的平均相對偏差Tab.2 Average relative deviation between new blend and original blend （%）

2.4 實例應(yīng)用結(jié)果及分析

基于上文結(jié)論，選擇相同的單等級片煙組1 和單等級片煙組2 作為輸入，分別采用OVO 和3V3_AS 方法進行一次卷煙配方維護實驗，推薦出替換單等級片煙及缺失單等級片煙信息如表3 所示。在表3 中，基于OVO 方法推薦的單等級片煙為DLC0568，因此當(dāng)單等級片煙AD0280 缺失時，可用推薦的單等級片煙DLC0568 替換缺失單等級片煙AD0280。3V3_AS 是基于平均相似度的三對三替換，因此當(dāng)配方中缺失單等級片煙AD0280、AD0603、DLC0277 時，用 單 等 級 片 煙DLC0977、DLC0245、AD0495 進行替換。

對按照表3 結(jié)果進行配方維護前后的卷煙樣品進行感官評吸和煙氣指標(biāo)測定，并計算平均絕對偏差，結(jié)果如表4 所示。表4 中，配方1和配方2 分別是OVO 和3V3_AS 配方維護實驗形成的新配方。由表4 可以看出，配方1 和配方2與原配方的感官和煙氣指標(biāo)的平均絕對偏差都較?。慌浞? 的平均絕對偏差遠(yuǎn)小于配方2 的平均絕對偏差，與上文結(jié)論一致；配方2 與原配方的感官質(zhì)量和煙氣指標(biāo)的平均絕對偏差相對略大一些，但是在可接受范圍之內(nèi)，可以進行維護。

表3 缺失的單等級片煙和模型推薦的單等級片煙的特征屬性及化學(xué)成分分析①Tab.3 Characteristics and chemical component analysis of missed and recommended single grade strips

表4 配方維護前后卷煙配方的感官和煙氣指標(biāo)①Tab.4 Sensory and smoke indexes of cigarette blend before and after blend maintenance

因此，實驗結(jié)果證明了此模型在進行一對一替換和多對多替換時，均能有效維護原品牌卷煙的感官質(zhì)量和煙氣指標(biāo)的穩(wěn)定。

3 結(jié)論

①提出了采用非負(fù)矩陣分解（Non-negative Matrix Factorization，NMF）方法構(gòu)建模型來完成卷煙配方的維護。②NMF 模型能夠檢索配方中被消除的單等級片煙，可以很好地維護配方的完整性。③NMF 卷煙配方維護模型可以較好地推薦與單等級片煙組形成良好組合的單等級片煙，使新卷煙配方中的單等級片煙具有良好的配伍性。在候選單等級片煙的選擇中，將單等級片煙的化學(xué)信息加入模型，提高了替換單等級片煙與缺失單等級片煙的相似性，增強了卷煙配方維護的可靠性。④一對一（OVO）替換和基于平均相似度方法的多對多（MVM_AS）替換都能有效維護原品牌卷煙的感官質(zhì)量和煙氣指標(biāo)的穩(wěn)定。綜上，本研究中提出的基于NMF 的卷煙配方維護模型，可為卷煙配方設(shè)計人員提供配方替代推薦，有助于企業(yè)維持卷煙感官質(zhì)量的穩(wěn)定性，提高工作效率。