煙葉總還原糖近紅外光譜穩(wěn)健模型的建立及其在多臺儀器的長期應用

束茹欣，居 雷，吳圣超，倪力軍，王文俊，欒紹嶸*，張立國

（1.上海煙草集團有限責任公司 卷煙煙氣重點實驗室，上海 201315；2.華東理工大學化學與分子工程學院，上海 200237）

近紅外光譜（NIRS）技術作為一種綠色、無損、快速的檢測方法，在各行各業(yè)得到了廣泛應用［1］，基于NIRS技術對產(chǎn)品進行質量檢測的行業(yè)標準、國家標準也不斷完善和推出［2］。該技術以樣品的近紅外光譜信號為自變量，樣品待測組分的定量及定性性質為因變量，利用多元統(tǒng)計或化學計量學方法建立定量或定性的NIRS模型，應用模型對未知樣品進行預測［3-4］。測試環(huán)境、樣品形態(tài)、儀器的工作狀態(tài)及樣品采收期等因素均會對光譜質量產(chǎn)生影響，并可能導致模型在后續(xù)時段及其他儀器上應用時的預測誤差增大［5-7］，需要更新模型或采用傳遞方法使模型能適應新的儀器或樣品［5，8-9］。而對于煙葉這類按年份或季節(jié)采收的農(nóng)產(chǎn)品，通常期望建立的NIRS模型不僅能在不同儀器間直接共享，還能有較長的使用壽命。

研究表明［10-12］，建立一個性能良好的NIRS模型未必需要所有的光譜信息。本課題組［11］前期基于圖像處理的尺度不變特征變換（SIFT）算法，根據(jù)若干主機代表性樣品的光譜信息篩選特征波長建立了煙葉總植物堿NIRS模型，該模型對主機及3臺從機樣品的預測誤差均能滿足企業(yè)內控要求。并進一步在SIFT 篩選波長的基礎上，通過剔除光譜信號間相關性高的波長及信號不穩(wěn)定波長，建立了黃芩苷的NIRS 模型［12］，該模型直接傳遞到從機后的預測誤差滿足測定要求。Ni 等［13］的工作表明，基于SIFT 篩選的重要且獨立性強的波長點建立的NIRS 模型，可直接傳遞到多臺從機實現(xiàn)對從機玉米樣品中蛋白質、油脂含量及煙葉總植物堿含量的準確預測。上述工作表明基于主機樣品光譜，采用合適的多步波長篩選可建立穩(wěn)健的近紅外光譜模型，從而實現(xiàn)模型的無標樣傳遞。

迄今有關NIRS 模型傳遞的報道多為模型在多臺儀器上轉移的研究，而NIRS 模型在同一儀器不同年份下的長期應用以及在不同儀器上的長期應用屬于NIRS模型傳遞的更廣泛的應用場景。這些場景在農(nóng)產(chǎn)品NIRS模型的應用中十分常見，但關于這類場景的研究和應用探索尚不多見［6-7］。

本文以煙葉總還原糖（TRS）NIRS 光譜模型在多臺儀器傳遞及在多臺儀器的長期應用為背景，探究了基于SIFT的多步波長篩選方法實現(xiàn)該目標的可行性，以改進模型的傳遞性能及使用壽命，減少模型傳遞及維護的工作量。

1 實驗部分

1.1 樣品及光譜采集

本文采用3 個煙葉樣本集建立模型并進行驗證。樣本集 A 為2011～2013 年間采收的292 個煙葉樣本，在傅里葉變換近紅外光譜儀A2（Thermo Fisher 公司，型號：ANTARIS Ⅱ NIR）上采集其光譜信息；樣本集 B由2011～2013年另外收集的77個煙葉樣本構成，分別在主機及6臺代號為C1、C2、J、N、P2和S 的從機（Thermo Fisher 公司，型號：ANTARIS Ⅱ NIR）上采集其光譜，用以檢驗模型的預測和傳遞性能。樣本集C 由2014～2020 年積累的180 個按采收年份分成7 個組別的煙葉樣本組成，分別在主機A2及上述6臺從機上采集其光譜，用以檢驗模型在多臺儀器的長期使用壽命。

所有煙葉樣品粉碎后，過40～60 目篩。置于可旋轉石英杯中，在（22±4） ℃、30%～80%濕度下采集其光譜：分辨率為8 cm-1，波長范圍為10 000～4 000 cm-1，掃描64 次，增益取2。各樣品的TRS 含量采用連續(xù)流動法［14］測定。3個樣本集中TRS的分布范圍、均值及方差等信息如表1所示。

表1 煙葉樣本集A、B及C中總還原糖含量信息Table 1 The information of TRS in tobacco samples of Set A，Set B and Set C

Set A中2011～2013年的292個煙葉樣品在主機A2上的近紅外光譜如圖1A所示，Set B中2011～2013年的77個樣品在主機A2及各從機上所測光譜的平均譜圖如圖1B所示。其中，Set A中的樣本用于波長篩選、建立和驗證煙葉總還原糖的NIRS 模型，Set B 中的樣本用于驗證模型的傳遞性能，Set C 用于檢驗模型在后續(xù)年度下在多臺儀器上的長期應用性能。

圖1 主機（A2）上采集的Set A中2011～2013年292個煙葉樣品的原始光譜（A）和Set B中77個樣本在7臺儀器上的平均紅外光譜（B）Fig.1 Spectra of the 292 samples in Set A collected in 2011-2013 tested on the primary（A2）（A） and average spectra of the 77 samples in Set B tested on the primary spectrometer（A2） and six secondary spectrometers（C1，C2，J，N，P2 and S）（B）

1.2 研究思路

從圖像分析角度而言，SIFT 算法［15-16］所篩選的特征波長集合Uc中各波長下的光譜信號與其相鄰信號存在明顯區(qū)別，被視為具有高度區(qū)分特征的波長點。SIFT 算法用于光譜特征點抽提及其參數(shù)的優(yōu)化在文獻［11-12］中已有介紹，本文不再贅述。

SIFT 算法中采用2～10 個代表性主機樣品光譜即可滿足特征波長篩選的需要，故本文首先采用SPXY（Sample portioning method based on X and Y）方法［17］從Set A 中篩選10 個代表性樣品，根據(jù)這10 個樣品的光譜篩選特征波長集合Uc?？紤]到SIFT算法的端部效應，需要結合其他波長篩選方法剔除冗余的、信息量少且不穩(wěn)定的波長點。眾所周知，不同樣品在同一波長下的光譜響應不同且這些光譜響應的差異程度隨波長的不同而改變。若不同樣品在某一波長下的光譜響應差異很大，那么該波長對于建立一個好的NIRS 模型非常重要；反之，如果各樣品在某波長下光譜響應的變化很小，則建立NIRS 模型時可以不必采用該波長的光譜信息。

基于上述考量，本課題組基于樣品光譜標準方差（SDSS）的大小來剔除Uc中信息量少、冗余的波長［13］，SDSS的計算公式如下：

其中，Mij是第i個樣本在波長j下的光譜響應值，m是Set A 中的樣本個數(shù)，是第j個波長下各樣本光譜響應的均值。SDSS（j）越大，則波長j下各樣品光譜響應的差異越大。由此可見，高SDSS值的光譜點可以很好地體現(xiàn)不同樣品間的光譜差異，是多元統(tǒng)計模型中的重要變量。

因此，本文在第二步波長篩選時，將SDSS值小于某一閾值b的波長從Uc中剔除，得到重要特征波長集合Uic，該法簡稱為SIFT-SDSS。眾所周知，水分及濕度對樣品近紅外光譜的穩(wěn)定性有很大影響，即使樣品的待測組分濃度相同，水分含量不同也將導致近紅外光譜形狀的改變［18］，如Gaines 與Windham 發(fā)現(xiàn)小麥的漫反射吸光度全光譜會隨著水分含量的增大而增大［19-21］。為提高所篩選波長下光譜響應的穩(wěn)定性，本文在SIFT-SDSS 篩選的重要特征波長集合Uic的基礎上進一步去除水分吸收系數(shù)大于8 000 的波段（分別位于4 000～4 058 cm-1與5 012～5 278 cm-1）［22］，挑選出重要、穩(wěn)定的特征波長集合Uisc，該波長篩選方法簡稱為SIFT-SDSS-MUS（moisture-unsensitive）。

分別基于Uc、Uic、Uisc與全波長光譜信息，選取Set A 中部分代表性樣本作為建模集，其余樣本作為驗證集，建立煙葉總還原糖的偏最小二乘（PLS）校正模型M-SIFT、M-SIFT-SDSS、M-SIFT-SDSSMUS與M-WW。采用Set B、Set C對各模型傳遞性能及長期應用性能進行檢驗。

1.3 煙葉總還原糖NIRS模型的建立與評價

根據(jù)課題組前期研究結果［23］，采用一階導數(shù)與二次多項式S-G 卷積平滑（光滑點數(shù)為9～15）對煙葉近紅外光譜數(shù)據(jù)進行預處理時，煙葉總植物堿NIRS模型的傳遞性能優(yōu)于采用其他預處理光譜所建的模型。故本文采用該方法對光譜進行預處理，光滑點數(shù)取15。在Set A 中的292個樣本中按照濃度區(qū)間挑選80%的樣本（234 個）作為校正集，剩余20%樣品（58 個）作為內部驗證集，以保障建模集和驗證集樣品個數(shù)沿濃度的分布相似。采用偏最小二乘回歸（PLSR）方法建立煙葉總還原糖的NIRS 模型。預測均方根誤差（RMSEP）及交叉驗證均方根誤差（RMSECV）通常用來評估近紅外光譜模型的預測性能。但是，均方根誤差的水平會隨樣本數(shù)據(jù)而改變，其閾值很難確定。煙草企業(yè)通常以NIRS模型預測值與實測值相對誤差絕對值的平均值（MARE）是否小于6%來判斷模型預測結果的準確度，當MARE小于6%時，則認為該模型可以接受。MARE的計算公式如下：

式中，m表示樣品數(shù)，yi，actual和yi，predicted分別表示第i個樣品的實際含量和預測含量。本文根據(jù)MARE小于6%的內控要求來評價煙葉總還原糖NIRS模型的傳遞及長期應用效果。

PLSR算法中的潛變量個數(shù)（nLVs）對模型結果有重要影響，通常采用留一交叉驗證確定。但本課題組在研究黃芩、玉米與煙葉中主要成分NIRS 模型傳遞過程中［12-13，22］發(fā)現(xiàn)，采用留一交叉驗證確定的nLVs個數(shù)偏大，易造成過擬合，導致模型傳遞性能變差。根據(jù)前若干個潛變量累積貢獻率大于或等于99.9%或者驗證集預測誤差最小確定nLVs時，模型的穩(wěn)健性及傳遞性能較好。故本文采用主機內部驗證集樣品的均方根殘差（RMSEV）最小來確定nLVs。

本文所有算法在Matlab2020a上完成。

2 結果與討論

2.1 煙葉光譜中重要且穩(wěn)定特征波長的篩選

SIFT算法所篩選的光譜點與算法中的參數(shù)O（組數(shù)）、S（層數(shù)）及初始高斯變換參數(shù)σ0密切相關。文獻［11-13］的研究表明，初始高斯變換參數(shù)σ0對篩選結果的影響最大，而O、S則無顯著性影響，建議S取值4～5，O取值3～5。本文參照課題組對煙葉總植物堿NIRS 模型的研究結果，選定SIFT 算法篩選煙葉特征波長點的參數(shù)O=S=5，σ0=3.0［13］。

基于“1.2”所述波長篩選流程及上述SIFT參數(shù)，本文在1.0×10-4～1.5×10-4范圍、間隔1.0×10-4選取不同的SDSS 閾值b，可得到不同的重要且穩(wěn)定特征波長集Uisc，基于各Uisc所建立的M-SIFT-SDSS-MUS 模型在不同nLVs 下的主機驗證集的RMSEV如圖2所示。

圖2 不同閾值b下模型M-SIFT-SDSS-MUS的RMSEV隨nLVs的變化Fig.2 The RMSEV predicted by different PLS models，which were built on Uisc obtained at different threshold of b along with nLVs

由圖2 可知，b=1.4×10-4時，基于波長集合Uisc所建的TRS模型在nLVs為9時，RMSEV最小。故SDSS 的閾值b設為1.4×10-4，此時Uisc中有56個波長，該閾值下SIFT-SDSS 方法所篩選的波長集合Uic中有72個波長。

由圖3可知，SIFT 所選的367個特征波長（以○表示）中很多位于10 000～7 200 cm-1區(qū)間，而6 800～5 300 cm-1、10 000～7 200 cm-1區(qū)間的SDSS 小于閾值1.4×10-4，說明這兩個區(qū)間的光譜信號波動非常小、其對NIRS 模型的影響可以忽略不計。SIFT-SDSS 方法剔除這些波長后得到的重要特征波長集合Uic有72 個點（用▲表示）。使用SIFT-SDSS-MUS 方法進一步從Uic中去除區(qū)間4 058～4 000 cm-1與 5 278～5 012 cm-1的水吸收系數(shù)大于8 000 的波長，最終獲得位于4 800～4 066 cm-1、7 200～6 800 cm-1及5 300～5 278 cm-1區(qū)間的56個重要、穩(wěn)定的煙葉近紅外光譜特征波長點（用*表示）。

圖3 校正集樣品的平均光譜（AS）、樣本標準方差譜（SDSS）和水的吸收系數(shù)譜、主機與從機樣品光譜絕對值差譜的平均譜（MASD）以及Uisc、Uic、Uc三個波長集合點Fig.3 Average spectrum（AS） and standard deviation of samples spectra（SDSS） of calibration samples，absorption coefficient spectrum of water，mean absolute difference spectrum（MASD） between primary and secondary spectra，and the wavelengths of Uisc，Uic，Uc

2.2 煙葉總還原糖NIRS 定量模型預測2011～2013年樣品的結果

表2 給出了基于不同波長集合所建各煙葉總還原糖NIRS 模型在不同儀器上測得的Set B 中77個樣品光譜預測TRS 的MARE。由表可知，全波長模型M-WW 和M-SIFT-SDSS-MUS 直接預測C1、C2、N、J、P2 5 臺從機樣品TRS 含量的MARE 均小于6%，僅對從機S 上樣品TRS 的預測誤差（分別為10.46%與6.07%）不滿足企業(yè)內控要求；而M-SIFT模型預測兩臺從機（C1、N）TRS的MARE大于6%，只有M-SIFT-SDSS模型直接傳遞到6臺從機的MARE 均小于6%，且其在7 臺儀器上的MARE 均值最?。?.02%），M-WW、M-SIFT 與M-SIFTSDSS-MUS 預測7 臺儀器測試樣品TRS 的MARE 均值分別為5.06%、4.56%與4.26%?？傮w而言，當各模型預測與建模集樣品同一時間段采集的Set B 中的77個樣品時，M-SIFT-SDSS 的傳遞性能最佳，MSIFT-SDSS-MUS次之，M-SIFT最差。

表2 不同模型預測Set B中77個樣品的MARE（%）Table 2 MARE（%） of TRS content in the 77 tobacco samples of Set B predicted by different models

表3 全波長模型（M-WW）預測Set C中各年度樣品中總還原糖的MARE（%）Table 3 MARE（%） of TRS content in the samples of Set C predicted by the M-WW model according to the sample collection year

表4 M-SIFT模型預測Set C中各年度樣品中總還原糖的MARE（%）Table 4 MARE（%）of TRS content in the samples of Set C predicted by the M-SIFT model according to the sample collection year

表5 M-SIFT-SDSS模型預測Set C中各年度樣品中總還原糖的MARE（%）Table 5 MARE（%）of TRS content in the samples of Set C predicted by the M-SIFT-SDSS model according to the sample collection year

表6 M-SIFT-SDSS-MUS模型預測Set C中各年度樣品中總還原糖的MARE（%）Table 6 MARE（%）of TRS content in the samples of Set C predicted by the M-SIFT-SDSS-MUS model according to the sample collection year

2.3 煙葉總還原糖NIRS定量模型預測2014～2020年Set C樣品的結果

表3～表6 分別給出了各模型根據(jù)2014～2020 年間樣品在7 臺儀器上測定的近紅外光譜所預測TRS含量的MARE以及按照儀器、年度得到的各MARE的均值。下面分幾個方面對結果進行討論。

2.3.1 煙葉總還原糖NIRS 定量模型在7臺儀器長期應用能力的考察由表3可知，根據(jù)NIRS模型預測TRS 的MARE 應小于6%的企業(yè)內控標準，全波長M-WW 模型在C1 和J 儀器上可以應用到2016 年；在A2、C2儀器可以應用到2015年；而在N、P2儀器該模型2015年就必須更新維護。M-WW 預測儀器S 上各年份樣品TRS 的MARE 為6.02%～16.42%，均大于6%。因此，對于在2014～2020 年間分年度收集的Set C中的樣品而言，其樣品采收期和測試時間不同于建模集樣品，M-WW 在7臺儀器上多年份應用的性能不佳，遠不及表2中該模型對Set B中2011～2013年樣品的預測性能。

由表4 可知，M-SIFT 模型在儀器S、N、C1、A2、P2、C2 上的使用壽命分別為7 年（2014～2020）、5 年（2014～2018）、3 年（2014～2016）、2 年（2014～2015）、2 年（2014～2015）和1 年（2014 年），模型在儀器J 上的長期應用性能最差：對2014 年J 儀器測試樣品TRS 的預測MARE 為8.93%、超過了6%。總體而言，M-SIFT 模型在7 臺儀器上的長期應用性能優(yōu)于M-WW 模型，表明基于SIFT 算法篩選特征波長點建立的模型較全波長模型更穩(wěn)健。

表5中，僅J儀器2014年樣品及P2儀器2020年樣品的MARE 大于6%（分別為7.43%和6.79%），其余MARE 均小于6%。說明M-SIFT-SDSS 模型在A2、C1、C2、N 和S 5 臺儀器上可以從2014 應用到2020 年，在P2 儀器可從2014 應用到2019 年。該模型在7 臺儀器的長期應用性能大大優(yōu)于M-SIFT 模型。表明在SIFT算法篩選的特征波長基礎上進一步挑選可良好區(qū)分樣品差異的重要波長點集合Uis，所建立的M-SIFT-SDSS模型的穩(wěn)健性及長期應用能力較M-SIFT模型有明顯改進。

表6 結果表明，M-SIFT-SDSS-MUS 模型可在C1、C2 與N 儀器上應用7 年（2014～2020）、在儀器S上應用4 年（2014～2017）、儀器A2 及P2 應用3 年（2014～2016），而在儀器J 上應用時，2014 年就必須更新模型。該模型在7臺儀器上的長期應用能力不及M-SIDT-SDSS模型，但優(yōu)于M-SIFT和M-WW 模型。結合表2 結果可知，在Uic重要特征波長集合基礎上，剔除對水分敏感的波長點，并不能改進煙葉總還原糖近紅外光譜模型M-SIFT-SDSS 的傳遞性能及其在多臺儀器的長期應用性能。造成這一結果的可能原因是，總還原糖中含有多個羥基，其在近紅外光譜區(qū)的特征吸收與水的強吸收區(qū)域重合度較高，剔除這些波長點反而會降低TRS 模型的準確度，并削弱模型的傳遞性能和長期應用能力。一般在一個模型的預測性能及傳遞性能優(yōu)良的情況下，進一步的波長篩選或者是采用標樣校正從機預測結果（光譜），得到的結果往往不如原模型。

2.3.2 驗證集樣本數(shù)及分布對煙葉總還原糖NIRS 定量模型預測結果的影響若不按年度區(qū)分，以Set C 中2014～2020 年所有180 個樣本為驗證集，則表3～表6 的最后一行Mean 即為各模型根據(jù)這180 個樣本在不同儀器上測試的光譜所預測TRS 的MARE。由表3 可知，應用全波長模型M-WW 預測這180個樣品TRS 時，A2、C1、C2 和J 4 臺儀器上的MARE（即Mean）值均小于6%，這意味著不區(qū)分2014～2020 年間180 個樣本的采收年份而將其視為整體時，采用2011～2013 年樣品建立的全波長TRS 模型可在主機及C1、C2 和J 3 臺從機上應用，但儀器N、P2 和S 應用該模型預測Set C 中180 個樣品TRS 含量的MARE在6.22%～9.46%之間，不滿足企業(yè)內控要求。

由表4最后一行可知，M-SIFT 模型除在A2和J儀器上預測180個樣品中TRS含量的MARE 大于6%外，在其他5 臺儀器的MARE 均小于6%。表5 和表6 的最后一行表明，7 臺儀器應用M-SIFT-SDSS 和M-SIFT-SDSS-MUS 兩個模型預測Set C 中不按年度劃分的180 個樣品的TRS 含量時，各儀器的MARE均小于6%。

經(jīng)考察建模集、Set B、Set C 樣品數(shù)目隨濃度區(qū)間的分布，發(fā)現(xiàn)建模集樣品、Set B 的77 個樣品及不按年度劃分時Set C 中所有180個樣品的分布圖接近正態(tài)分布，這3個圖非常相似。而Set C 中各年度樣品數(shù)目在20～30 之間，每個年度樣品隨濃度區(qū)間的分布與建模集差異很大。從統(tǒng)計學角度而言，對于20～30范圍的樣本空間，如果有個別樣本的相對誤差絕對值較大，會導致MARE顯著增大。而在77～180 個樣本中有若干相對誤差大的樣本出現(xiàn)時，這一風險會小很多。由此不難理解為何表3、表4 中，M-WW 和M-SIFT 模型預測各儀器多年份樣品的TRS 誤差均很大，但模型用于總體預測180 個樣品時，在多臺儀器的MARE 均小于6%。故增加2014～2020 各年度煙葉樣品數(shù)目并使各年度樣品數(shù)隨TRS 濃度分布圖與建模集盡可能相似，有望降低各模型預測2014～2020 年各年度樣品TRS 的誤差，提高模型在多臺儀器上長期應用的能力。

3 結 論

本文采用多步波長篩選方法，綜合光譜特征、對樣本的區(qū)分度以及對水分的敏感度等因素進行建模，發(fā)現(xiàn)基于從若干代表性主機樣品光譜中挑選重要（或重要且穩(wěn)定）特征波長的兩步與三步波長篩選方法所建的M-SIFT-SDSS 及M-SIFT-SDSS-MUS 模型的傳遞能力及服役年限均顯著優(yōu)于全波長模型，且優(yōu)于一步波長篩選法所建的M-SIFT 模型。其中M-SIFT-SDSS 模型的傳遞能力及在多臺儀器上的長期應用能力最佳，其模型變量個數(shù)僅為全波長模型的4.9%。該模型用于預測2011～2013 年樣品時，可在6 臺從機上直接共享；用于預測2014～2020 年各年度樣品時，可在5 臺儀器的考察年度內長期應用7年、在一臺儀器上應用6年（2014～2019）。本文提出的基于SIFT算法的多步波長篩選方法可實現(xiàn)TRS近紅外光譜模型在6 臺從機上的無標樣傳遞，并延長模型在各儀器上的使用壽命。本方法對其他數(shù)據(jù)集的適應性有待進一步研究。