改進的蒙特卡洛方法用于煙絲近紅外光譜定量分析中奇異樣本的識別

胡蕓，劉劍，白興，阮藝斌，張辭海，姬厚偉，李博巖

1 貴州中煙工業(yè)有限責(zé)任公司技術(shù)中心，貴陽市小河經(jīng)濟技術(shù)開發(fā)區(qū)開發(fā)大道96 號 550009；

2 貴州醫(yī)科大學(xué)公共衛(wèi)生與健康學(xué)院，貴陽市花溪大學(xué)城 550025

在煙草及煙草制品的總植物堿、總糖、總氮等化學(xué)組分的近紅外光譜定量分析中，校正模型往往是在大量樣本基礎(chǔ)上建立的。模型質(zhì)量在很大程度上依賴于建模樣本數(shù)據(jù)的準(zhǔn)確性。異常樣本的存在會嚴重干擾定量模型質(zhì)量，影響模型預(yù)測的準(zhǔn)確性[1-4]。因此，奇異樣本的識別是建立可靠準(zhǔn)確的定量分析模型的前提。

奇異樣本一般是指那些落在主體樣本之外的樣本，其觀測數(shù)據(jù)常偏離主體樣本呈現(xiàn)的模式。造成奇異樣本的因素可能是實驗操作誤差、儀器誤差、測量誤差以及樣本性質(zhì)變化等。奇異樣本的診斷與識別一直都是多元校正研究中的熱點，涉及的方法主要有經(jīng)典識別方法（如殘差法、馬氏距離法[5]、帽子矩陣法等）和穩(wěn)健回歸方法（如最小一乘估計、極大似然估計等）。診斷方法著重于先檢測奇異樣本，剔除后再用經(jīng)典方法建模。穩(wěn)健回歸方法是在解析過程中估計數(shù)據(jù)分布的主體，然后構(gòu)建穩(wěn)健模型，從而降低數(shù)據(jù)中奇異樣本對回歸估計產(chǎn)生的不利影響。

劉智超等[6]基于樣本在蒙特卡洛交叉驗證中的統(tǒng)計規(guī)律，建立一定數(shù)量的模型，然后依據(jù)預(yù)測誤差平方和預(yù)測殘差平方和（PRESS）排序并統(tǒng)計每個樣本在不同模型中出現(xiàn)的頻次對異常樣本進行識別。Cao等[7]通過研究樣本在模型群中預(yù)測誤差的分布，發(fā)現(xiàn)預(yù)測誤差的統(tǒng)計特征可用于識別正常樣本和奇異樣本，由此提出蒙特卡洛方法奇異樣本檢測方法（MonteCarlo outlier detection, MCOD）。與傳統(tǒng)方法相比，MCOD 在一定程度上降低了由掩蔽效應(yīng)帶來的模型風(fēng)險，具有較高的奇異樣本識別能力。在此基礎(chǔ)上，Zhang 等[8-9]提出使用正常樣本建模，對可疑異常樣本進行二次驗證，大大降低了可疑異常樣本的誤判和掩蔽效應(yīng)帶來的建模風(fēng)險。

本文首次把改進的蒙特卡洛奇異樣本（improved MonteCarlo outlier detection, IMCOD）識別方法應(yīng)用于卷煙成品煙絲的化學(xué)成分（總植物堿、總糖、還原糖、總氮、鉀和氯）的定量建模中，提高校正模型的穩(wěn)健性與預(yù)測準(zhǔn)確性，具有較大的應(yīng)用意義。

1 材料與方法

1.1 儀器與材料

Antaris 傅立葉近紅外分析儀（美國Thermo 公司）、Futura 型8 通道連續(xù)流動化學(xué)分析儀（法國Alliance公司）、FED260 型烘箱（德國Binder 公司）、410 型單通道火焰光度計、VELP 消化爐、ML204 型電子分析天平（梅特勒-托利多儀器上海有限公司）、HY-6 型雙層調(diào)速多用振蕩器（江蘇金壇榮華儀器制造有限公司）。

192 個卷煙樣品由貴州中煙工業(yè)有限責(zé)任公司提供。

1.2 實驗方法

1.2.1 樣品處理

每個卷煙樣品取2 盒，然后剝開煙支取煙絲混合均勻后裝入一次性密封袋低溫（0℃～4℃）保存，待測。

1.2.2 化學(xué)值測定方法

按照煙草行業(yè)標(biāo)準(zhǔn)，測定了煙葉樣品中總糖、還原糖、煙堿、總氮、氯和鉀6 種常規(guī)化學(xué)成分的含量。相關(guān)標(biāo)準(zhǔn)如下：《YC/T 31—1996 煙草及煙草制品試樣的制備和水分測定烘箱法》；《YC/T 468—2013 煙草及煙草制品總植物堿的測定連續(xù)流動(硫氰酸鉀)法》；《YC/T 159—2002 煙草及煙草制品水溶性糖的測定連續(xù)流動法》；《YC/T 161—2002 煙草及煙草制品總氮的測定連續(xù)流動法》；《YC/T 217—2007 煙草及煙草制品鉀的測定連續(xù)流動法》；《YC/T 162—2002 煙草及煙草制品氯的測定連續(xù)流動法》。

1.2.3 光譜采集與數(shù)據(jù)處理

實驗室測試環(huán)境溫度控制在（22±2）℃，相對濕度控制在（40±10）%。近紅外儀器開機預(yù)熱不低于1 h，用ValPro 程序校檢合格后使用。低溫保存過的煙絲樣品應(yīng)在實驗室測試環(huán)境下平衡不少于4 h。將裝有煙絲（壓平后，樣品厚度≥20 mm）的樣品杯放置于旋轉(zhuǎn)臺，采用積分球漫反射的方式測量光譜，掃描范圍為4000～10000 cm-1，分辨率為8 cm-1，掃描次數(shù)為64 次。每個樣品重復(fù)掃描3 次，每次采集樣品光譜前均須測量背景光譜。所使用的軟件為儀器自帶的RESULT-Integration 工作流（Workflow）設(shè)計軟件、RESULT-Operation 操作軟件和TQ Analyst 光譜分析軟件。

1.3 算法

1.3.1 蒙特卡洛奇異樣本檢測方法

一般來說，穩(wěn)健診斷方法多是基于評價模型殘差，即擬合殘差之上。但由于奇異點的存在，擬合殘差對奇異點是不靈敏的，只有預(yù)測殘差才對奇異點很敏感，它是擬合誤差與Cook 距離之和。這就是蒙特卡洛采樣方法檢測奇異點的原理[7]。

圖1 為蒙特卡洛采樣方法的示意圖。該方法對m個樣本數(shù)據(jù)進行N 次采樣，N 值可以設(shè)定為樣本數(shù)的平方或者一個較大的數(shù)值。每次采樣時，按照一個設(shè)定的比例（如70%～80%）隨機選擇樣本，如圖中白色標(biāo)記的數(shù)據(jù)樣本點，構(gòu)建校正模型，再預(yù)測剩余的樣本（黑色標(biāo)記）。N 次采樣后，可獲得參數(shù)估計分布的蒙特卡洛近似分布，由此評估和分析樣本預(yù)測誤差。

圖1 蒙特卡洛采樣示意圖Fig.1 Schematic diagram of Monte Carlo Sampling

基于不同性質(zhì)的樣本點具有不同的分布，可以計算每個樣本點預(yù)測誤差分布的均值的絕對值（mean value）和標(biāo)準(zhǔn)偏差（standard deviation）。以預(yù)測誤差的均值的絕對值為X 軸，預(yù)測誤差的標(biāo)準(zhǔn)偏差為Y 軸作圖，得到一個四區(qū)域散點圖（圖2）。其中，左下角區(qū)域涵蓋了大部分的數(shù)據(jù)樣本點，對應(yīng)的均值和偏差值都比較小，表示正常樣本；左上角區(qū)域是X 方向的奇異樣本，具有較大的偏差和較小的均值；右下角區(qū)域是Y 方向的奇異樣本或者模型奇異樣本，對應(yīng)小的偏差和大的均值；右上角區(qū)域是特殊的奇異樣本或者不正常樣本，可能是某些誤差造成的。該區(qū)域散點圖提供了奇異樣本可視化診斷信息，可根據(jù)實際情況選擇多數(shù)主體樣本的預(yù)測誤差和標(biāo)準(zhǔn)偏差均值的2～2.5倍作為區(qū)域劃分的量化依據(jù)[7]。

圖2 基于樣本預(yù)測誤差統(tǒng)計的不同類型樣本的判斷Fig.2 Determinant of different types of samples based on statistics of prediction errors

1.3.2 改進的蒙特卡洛奇異樣本檢測方法

在利用蒙特卡洛奇異樣本檢測方法時，掩蔽效應(yīng)的存在可能使得正常樣本與異常樣本的區(qū)域分界線不是很清晰。所謂的掩蔽效應(yīng)是指多個奇異樣本存在于數(shù)據(jù)時，一些奇異樣本會使得其它奇異樣本不能被檢測方法診斷出來。因此，提出了改進的蒙特卡洛奇異樣本檢測方法。IMCOD 的基本步驟如圖3 所示：（1）循環(huán)N 次MCOD 采樣[7]，計算全部樣本的預(yù)測誤差分布的均值與標(biāo)準(zhǔn)偏差；（2）將具有較小預(yù)測誤差均值和標(biāo)準(zhǔn)偏差的樣本視為正常樣本，剩下的樣本視為可疑樣本；（3）按一定比例把正常樣本數(shù)據(jù)集隨機劃分為訓(xùn)練集和驗證集；（4）用訓(xùn)練集樣本構(gòu)建校正模型，然后預(yù)測驗證集樣本和可疑樣本，計算預(yù)測誤差并保存；（5）重復(fù)步驟（3）和（4）N 次后，得到正常樣本和可疑樣本的預(yù)測誤差分布的均值和標(biāo)準(zhǔn)差，判斷異常樣本。這些樣本預(yù)測誤差是在正常樣本建模的基礎(chǔ)上產(chǎn)生的，故能有效降低掩蔽效應(yīng)的影響，更加清晰地判斷正常和異常樣本的邊界。值得注意的是與MCOD 方法相比，IMCOD 方法采用了已篩選出的正常樣本構(gòu)建模型來預(yù)測可疑樣本，可能增加了模型的過擬合風(fēng)險。針對該問題，使用蒙特卡洛交互驗證方法[9]，合理估計訓(xùn)練集正常樣本校正模型的維度（即主成分數(shù)），以有效降低去除可疑樣本后數(shù)據(jù)集中樣本減少而帶來的過擬合風(fēng)險。一般來說，模型的主成分數(shù)越大，過擬合的可能性越大，模型越不穩(wěn)健，預(yù)測能力越弱。樣本數(shù)據(jù)的最后的回歸模型應(yīng)該是由去除IMCOD 方法檢測的異常樣本后，合并非異常樣本與正常樣本所構(gòu)建的。

圖3 改進的蒙特卡洛奇異樣本檢測方法的流程圖Fig.3 Flowchart of IMCOD method

2 結(jié)果與討論

2.1 卷煙樣品的近紅外光譜與化學(xué)參考值

近紅外光譜的奇異樣本主要涉及光譜異常和化學(xué)異常樣本。光譜異常樣本又稱為自變量X方向的奇異樣本，在高維變量空間中遠離樣本分布主體。一般來說，成品卷煙煙絲的近紅外光譜比較復(fù)雜（圖4），樣本的顆粒大小及樣本背景成分的差異均可以引起光譜差異，可能造成一些樣本的光譜不屬于樣本的總體分布，成為奇異樣本。另外，光譜的相似性高，若不借助化學(xué)計量學(xué)方法，很難直接快速識別光譜異常樣本。近紅外光譜定量分析模型建立的另一個關(guān)鍵要素是借助標(biāo)準(zhǔn)分析方法獲得樣本的化學(xué)參考值。這些數(shù)值的測量可能是有誤的，從而形成化學(xué)異常樣本，也稱為響應(yīng)變量Y方向的奇異樣本。這類奇異樣本經(jīng)常遠離化學(xué)指標(biāo)Y的數(shù)據(jù)主體，在建模過程中經(jīng)常會引起較大的預(yù)測誤差。因此，化學(xué)異常樣本的識別是建立一個穩(wěn)定模型的關(guān)鍵因素。表1 列出了成品卷煙中6 種化學(xué)成分的含量信息。若僅從數(shù)據(jù)的角度看，很難得知樣品中總植物堿、總糖、還原糖、氯、總氮和鉀化學(xué)參考值是否正常。再者就是結(jié)合了以上兩種性質(zhì)的奇異點，即同時在X方向與Y方向的奇異點。在建立校正模型之前剔除奇異樣本，能夠建立穩(wěn)健的校正模型，從而提高模型對樣本的預(yù)測準(zhǔn)確性。

表1 卷煙樣品中6 種化學(xué)成分的含量信息Tab.1 The contents of six components in cigarette cut samples

圖4 卷煙煙絲樣品的原始近紅外光譜Fig.4 Original NIR spectra of cigarette cut samples

2.2 卷煙煙絲奇異樣本識別

采用化學(xué)計量學(xué)多元散射校正（MSC）和SavGol一階導(dǎo)數(shù)方法對卷煙煙絲樣本的近紅外光譜原始數(shù)據(jù)進行預(yù)處理，使用4200～7600 cm–1區(qū)域內(nèi)的光譜數(shù)據(jù)構(gòu)建定量分析模型（圖4）。

對于煙絲中的6 種化學(xué)成分含量，依照上MCOD和IMCOD 方法檢測步驟，分別選取80%的樣本作為訓(xùn)練集，建立偏最小二乘（PLS）定量模型，預(yù)測剩余的樣本的化學(xué)成分含量。隨機采樣2000 次，得到每個樣本的預(yù)測誤差，計算其對應(yīng)的平均值的絕對值和標(biāo)準(zhǔn)偏差，如圖5、6 所示。大部分樣本都處于四分區(qū)圖的左下角（零值附近），也就是說這些正常樣本的預(yù)測誤差的均值和標(biāo)準(zhǔn)偏差都很小。根據(jù)煙絲樣本中化學(xué)成分的不同，選擇預(yù)測誤差的2.5～3.5 倍和標(biāo)準(zhǔn)偏差均值的1.5～2.0 倍作為判別正常樣本與奇異樣本的量化指標(biāo)。

圖5 MCOD(a、b、c)和IMCOD(d、e、f)方法診斷卷煙樣品中總植物堿、總糖和還原糖的異常樣本Fig.5 Outlier detection for total alkaloids, total sugars and reducing sugars in cigarette samples using by MCOD and IMCOD

圖6 MCOD(a、b、c)和IMCOD(d、e、f)方法診斷卷煙樣品中總氮、鉀和氯的異常樣本Fig.6 Outlier detection for total nitrogen, potassium and chlorine in cigarette samples using by MCOD and IMCOD

圖5a 中MCOD 方法檢測總植物堿的結(jié)果表明8個樣本（#148、#153、#157、#158、#159、#160、#179、#192）與左下方區(qū)域的正常樣本差異明顯，診斷為異常樣本。樣本#32 和#35 接近正常樣本，不易識別。在此基礎(chǔ)上，選取預(yù)測誤差的均值小于0.15%和平均偏差小于0.04%的樣本視為正常樣本集，運行IMCOD方法得到圖5d 的預(yù)測結(jié)果。可以看出，正常樣本與異常樣本的界限更加清晰，樣本#32 與#35 和正常樣本的區(qū)別更加顯著。原先被視為異常樣本的#192 樣本的預(yù)測誤差變小，分布在正常樣本區(qū)域范圍之內(nèi)。樣本#11、#12、#15 和#57 被識別為異常樣本。

圖5b 為總糖的MCOD 方法診斷異常樣本結(jié)果：樣本#157 和#159 的預(yù)測誤差的標(biāo)準(zhǔn)偏差較大，可能為光譜異常點。樣本#35 的預(yù)測誤差的均值較大，判斷為化學(xué)異常點。用IMCOD 方法則識別出8 個異常樣本（圖5e），#148 和#160 也為光譜方向的奇異樣本。樣本#31、#98 和#158 也是奇異樣本。類似地，用MCOD和IMCOD 方法分別識別還原糖、總氮、鉀和氯含量樣品集中的奇異樣本，除還原糖外，其他5 種化學(xué)成分相關(guān)的異常樣本中都有樣本#157、#158、#159 和#160。近紅外光譜數(shù)據(jù)的主成分分析結(jié)果表明這些樣本的光譜殘差遠大于主體樣本的光譜殘差（圖7），驗證這四個樣本的光譜異常。實際上，奇異樣本的識別并不容易。若實驗條件允許，應(yīng)通過重復(fù)篩選出的光譜異?；蚧瘜W(xué)成分異常樣本的光譜掃描和化學(xué)成分分析實驗，來確證識別的正確性。

圖7 192 個卷煙煙絲的主成分分析殘差圖Fig.7 Plot of PCA residuals of 192 cigarette cut samples

使用原始的全部樣品建立PLS 定量模型，并與用MCOD 和IMCOD 方法剔除異常樣本后建立的模型進行比較（表2）。MCOD 和IMCOD 方法剔除異常樣本后，化學(xué)成分的定量分析模型的相關(guān)系數(shù)和校正均方根誤差都得到明顯改善和提高。例如，總糖含量模型的相關(guān)系數(shù)從0.936 提高到0.966，校正均方根誤差從0.913%降到0.648%。IMCOD 方法檢測奇異樣本是基于正常樣本建立預(yù)測模型基礎(chǔ)上的，降低了多個奇異點的掩蔽效應(yīng)，更加準(zhǔn)確地識別異常樣本。

表2 全部樣本建模和剔除異常樣本建模的結(jié)果比較Tab.2 Comparison of modeling results between sample with outliers and sample without outliers

2.3 卷煙樣品的近紅外光譜定量模型

采用MCOD 和IMCOD 方法分別剔除異常樣本后，優(yōu)化和建立了成品卷煙煙絲的6 種化學(xué)成分的定量分析模型（表3）。由表可知，用IMCOD 方法剔除異常樣本后，所建模型的相關(guān)系數(shù)得到提高，且校正集均方根誤差（RMSEC）和預(yù)測集均方根誤差（RMSEP）變小。IMCOD 方法將總氮含量的校正模型的相關(guān)系數(shù)從0.874 提高到0.919，預(yù)測均方根誤差降低了14.8%，模型的預(yù)測能力得到很大改善。計算模型穩(wěn)健性評價參數(shù)SEP/SEC 的值均小于1.2[10]，表明模型對預(yù)測樣本具有較好的穩(wěn)定性。定量分析模型的相關(guān)系數(shù)（R）均大于0.91。這是因為樣品粒徑及其分布狀態(tài)會影響近紅外光譜的預(yù)測效果。固體樣品的粒徑大小改變了近紅外光在樣品粒子中的光程，影響近紅外光對樣品的反射性和穿透性，體現(xiàn)在樣品對近紅外光的吸收和散射強度。因此，粒徑較大或粒徑不均勻樣品的預(yù)測準(zhǔn)確性和精密性會差一些。在實際應(yīng)用中，若使用煙絲模型預(yù)測未知樣本的含量，雖然準(zhǔn)確度有所降低，但考慮到煙絲的模型不需要樣品制打磨步驟，可以大大節(jié)約時間和成本。

表3 卷煙煙絲樣品中6 種化學(xué)成分的模型結(jié)果Tab.3 The model results of six components in cigarette cut samples

3 結(jié)論

改進的蒙特卡洛奇異樣本檢測方法主要是在正常樣本的基礎(chǔ)上構(gòu)建N 個定量模型來獲得全部樣本預(yù)測誤差的分布，通過樣本分布間的差異性判斷正常樣本和奇異樣本。方法的應(yīng)用不僅能夠得到正常樣本與異常樣本的區(qū)域散點圖，還提供了奇異樣本可視化的診斷信息，提高了奇異樣本識別的準(zhǔn)確性，有效降低了掩蔽效應(yīng)帶來的建模風(fēng)險。卷煙成品煙絲中6 種化學(xué)成分的定量分析模型的建立和應(yīng)用也證實了該方法的合理性和準(zhǔn)確性。對于邊界附近的異常樣本，可以通過優(yōu)化正常樣本的預(yù)測誤差平均值和標(biāo)準(zhǔn)偏差，增加運行次數(shù)等進一步確定。這也說明了在保證光譜數(shù)據(jù)和基礎(chǔ)數(shù)據(jù)準(zhǔn)確的前提下，才能建立準(zhǔn)確的可靠的近紅外分析模型。