王 迪，馮偉華，郭軍偉，王 銳，劉惠民，宗國浩，劉紹鋒，王永勝，趙 樂

中國煙草總公司鄭州煙草研究院，鄭州高新技術(shù)產(chǎn)業(yè)開發(fā)區(qū)楓楊街2號 450001

近紅外光譜分析技術(shù)是一種結(jié)合儀器科學(xué)、化學(xué)計量學(xué)等學(xué)科對樣品定性或定量的二次分析技術(shù)，以無損、快速、無污染等優(yōu)點而受到越來越多研究人員的青睞，已經(jīng)在食品、藥物、石油、農(nóng)業(yè)及煙草等領(lǐng)域得到廣泛應(yīng)用［1-3］。為了提高近紅外預(yù)測模型的適用性，往往通過模型轉(zhuǎn)移方法將主機模型應(yīng)用到子機儀器上以避免重復(fù)建模［4］，常用的模型轉(zhuǎn)移算法如截距/斜率校正（S/B）［5］、光譜空間變換（SST）［6］、直接標準化（DS）［7］等均需要主機光譜與子機光譜的波數(shù)一致才能夠計算。在日常分析過程中發(fā)現(xiàn)，即使是相同品牌的傅里葉近紅外光譜儀，在更新?lián)Q代之后，由于激光器的變化，導(dǎo)致其光譜成像波數(shù)不一致，造成模型轉(zhuǎn)移算法失效，而插值方法作為有效解決這一類似問題的方法之一，常被用在模型轉(zhuǎn)移前處理過程中，使子機光譜與主機光譜變換一致后再進行模型轉(zhuǎn)移。

插值［8］是離散函數(shù)逼近的重要方法，根據(jù)函數(shù)在已知有限點處的取值，利用插值方法能夠估算出函數(shù)在其他未知點處的近似值。通過研究Zero、Slinear、Quadratic、Cubic、Nearest［9-11］5種插值方法，分別對光譜數(shù)據(jù)進行插值及Savitzky-Golay［12］（SG）平滑處理。而SG作為近紅外光譜經(jīng)典平滑預(yù)處理方法被廣泛采用，本項目組在前期的研究過程中曾使用該方法對光譜預(yù)處理后實現(xiàn)了建模應(yīng)用。在此基礎(chǔ)上，本研究中提出了一種新的基于SG平滑的插值方法（SG-Inter）對光譜進行變換，從而達到將主機、子機兩種儀器光譜補齊的目的；通過利用Zero、Slinear、Quadratic、Cubic、Nearest以及SG-Inter方法對光譜進行插值處理得到新的子機光譜數(shù)據(jù)后，再經(jīng)過光譜空間變換（SST）將光譜轉(zhuǎn)移后預(yù)測煙草總植物堿、總糖、總氮、還原糖、氯、鉀、淀粉、新植二烯和綠原酸共9種具有代表性的煙草化學(xué)指標的質(zhì)量分數(shù)，根據(jù)預(yù)測結(jié)果統(tǒng)計、分析、評價不同插值方法對模型轉(zhuǎn)移效果的影響。

1 材料與方法

1.1 材料、試劑和儀器

收集2020年全國各省級中煙工業(yè)有限責(zé)任公司的初烤煙葉樣品，共800個；采用Kennard-Stone（KS）方法［13］篩選出200個代表性樣品，分別采集其近紅外光譜。

主機MPA（1代）和子機TANGO（2代）光譜儀（德國Bruker公司）；ZM200型粉碎機（德國Retsch公司）；BSA124S型電子天平（感量0.000 1 g，德國Satorious公司）；AA3連續(xù)流動分析儀（德國BRAN+RUBBE公司）。

1.2 方法

1.2.1 樣品處理與化學(xué)指標檢測

參考行業(yè)標準［14］干燥處理樣品，直至可用手指捻碎。將樣品通過粉碎機粉碎研磨，并過0.250 mm（60目）分樣篩，混勻后裝入密封袋。分別采用行業(yè)標準，無相關(guān)標準方法的采用文獻方法，測定樣品中的總植物堿［15］、總糖和還原糖［16］、總氮［17］、氯［18］、鉀［19］、淀粉［20］、新植二烯［21］以及綠原酸［22］的質(zhì)量分數(shù)。

1.2.2 光譜采集與預(yù)處理

將儀器的激光能量范圍設(shè)置為4 000～10 000 cm-1，分辨率設(shè)置為8 cm-1，掃描次數(shù)設(shè)置為64次，分別用主機MPA（1代）和子機TANGO（2代）在相同的條件下采集近紅外光譜，采用Savitzky-Golay平滑（窗口17，2次多項式，1階導(dǎo)數(shù)）進行光譜預(yù)處理。

1.2.3 主機模型

選擇項目組在前期研究中所構(gòu)建的模型作為煙草近紅外分析主機模型，用于本研究中模型轉(zhuǎn)移效果評價。

1.2.4 插值方法的研究與應(yīng)用

采用Zero、Slinear、Quadratic、Cubic、Nearest 5種插值方法分別對原始光譜以及SG平滑后的光譜數(shù)據(jù)進行插值處理，將子機數(shù)據(jù)與主機數(shù)據(jù)補齊，隨后利用光譜空間變換（SST）方法對插值后的光譜進行處理，并預(yù)測模型轉(zhuǎn)移后煙葉中總植物堿、總糖、總氮、還原糖、氯、鉀、淀粉、新植二烯和綠原酸的質(zhì)量分數(shù)。

1.2.5 基于Savitzky-Golay的平滑插值方法

基于Savitzky-Golay平滑插值（SG-Inter）的模型轉(zhuǎn)移，即采用平滑處理與插值同步計算來進行模型轉(zhuǎn)移。本研究中，SG-Inter（窗口大小為17，1階導(dǎo)數(shù)，2次多項式）隨滑動窗口在子機光譜上的移動，求得窗口內(nèi)子機光譜數(shù)據(jù)的2階多項式并1階求導(dǎo)，把子機光譜SG平滑點的最鄰近主機光譜波長點值代入到函數(shù)中求得插值。把插值變換融合到數(shù)據(jù)預(yù)處理過程中，保留原始光譜的變化趨勢，將波數(shù)為1 456的子機光譜平滑插值到與主機光譜相同的1 555個波長點，子機與主機光譜波數(shù)調(diào)整一致后進行模型轉(zhuǎn)移并預(yù)測總植物堿、總糖、總氮、還原糖、氯、鉀、淀粉、新植二烯和綠原酸的質(zhì)量分數(shù)。具體計算原理如下：

Savitzky-Golay平滑是一種卷積滑動窗口的加權(quán)平均算法，設(shè)濾波窗口的寬度w=2i+1，i是半窗寬度，i=1，2，3，…，n；x代表數(shù)據(jù)點在窗口內(nèi)的相對位置，其取值為［-i，…，0，…，i］，數(shù)據(jù)點所在位置對應(yīng)的函數(shù)值為P（x）。根據(jù)窗口內(nèi)的數(shù)據(jù)點，構(gòu)造一個n階多項式進行擬合得到f（x）表達式，見公式（1）［12］：

式中：cn0，cn1，…，cnn分別代表擬合函數(shù)f（x）中擬合數(shù)據(jù)點的系數(shù)，經(jīng)過最小二乘擬合，得到殘差E的表達式，見公式（2）［12］：

目標是使殘差E趨于最小，將公式（2）中各項系數(shù)導(dǎo)數(shù)εz分別設(shè)置為0，z=（0，1，2，3，…，n），得到公式（3）［12］：

將公式（3）化簡得到公式（4）：

當滑動窗口大小與平滑階數(shù)固定時，將待擬合窗口［P（-i），…，P（0），…，P（i）］內(nèi)數(shù)據(jù)帶入公式（4），可求得多項式系數(shù)列表［cn0，cn1，…，cnn］T。如圖1所示，設(shè)平滑窗口w=5，每次求解窗口內(nèi)第w/2=3個位置的平滑多項式（實心圓點為原始信號點，空心圓點為待插值信號點，實心方點為待插值信號點的最鄰近原始信號點），隨著卷積核窗口以步長為1向前平移，Savitzky-Golay平滑方法將重新擬合卷積核窗口平移后的平滑多項式，并求得窗口中心點所在位置的平滑后數(shù)值。本研究中將已知光譜波長點的對應(yīng)數(shù)值xi與最鄰近待預(yù)測波長點對應(yīng)數(shù)值xj進行等比例縮放得到比例系數(shù)m=xi/xj，根據(jù)已知波長點值的位置yi，求得待預(yù)測波長點的相對位置yj=yi·m，將yj代入到平滑多項式f（x）中求解，得到待預(yù)測波長點的數(shù)據(jù)值。經(jīng)過卷積核窗口的滑動，最終得到平滑在線插值后的光譜數(shù)據(jù)列表［f（y0），f（y1），…，f（yn）］，再經(jīng)過一系列數(shù)據(jù)預(yù)處理操作后，將轉(zhuǎn)移后光譜數(shù)據(jù)列表應(yīng)用于主機模型進行化學(xué)成分預(yù)測分析。

圖1 平滑插值（SG-Inter）原理圖（w=5）Fig.1 Smoothing interpolation theory（Window=5）

1.2.6 光譜轉(zhuǎn)移及化學(xué)成分預(yù)測

以MPA（1代）掃描的標準樣品光譜為主機光譜，MPA（2代）掃描的標準樣品光譜為子機光譜，將子機光譜平滑求導(dǎo)后插值處理，使子機光譜與主機光譜的橫坐標光譜波數(shù)調(diào)整一致，并按照光譜空間變換方法（SST）［6］將調(diào)整后的子機光譜轉(zhuǎn)移至主機光譜，使用主機模型對轉(zhuǎn)移后的子機光譜進行化學(xué)指標分析預(yù)測，待得到指標預(yù)測結(jié)果后再統(tǒng)計分析。

1.2.7 數(shù)據(jù)處理方法

使用計算機編程語言Python實現(xiàn)本研究中的各種模型計算。

2 結(jié)果與分析

2.1 插值方法對光譜的直接影響

2.1.1 原始光譜直接插值后的譜圖分析

使用Zero、Slinear、Quadratic、Cubic、Nearest 5種插值方法對子機原始光譜進行插值處理，并與原始光譜對比得到圖2。可知，原始光譜能夠與這幾種插值方法處理后的光譜數(shù)據(jù)基本重合。說明插值方法能夠解決子機光譜與主機光譜波數(shù)不一致的問題。

圖2 光譜直接插值譜圖對比Fig.2 Comparison of spectra after direct interpolation

模型轉(zhuǎn)移過程中通常需要對光譜進行平滑求導(dǎo)預(yù)處理以消除儀器噪音的影響，將光譜平滑處理對比后得到圖3。從圖3中能明顯看出不同方法插值后的光譜與原始光譜平滑處理后差異較大，通過對幾處差異明顯的出峰點局部放大后，可以清晰、直觀地看出各插值方法處理后的光譜與原始光譜的差異。由于原始光譜直接插值時只考慮了插值區(qū)間內(nèi)的數(shù)據(jù)處理，沒有將數(shù)據(jù)點左右兩側(cè)數(shù)據(jù)值變化趨勢考慮在內(nèi)，但光譜SG平滑預(yù)處理時，需要對每個平滑窗口的表達式進行1階求導(dǎo)，即光譜直接插值的結(jié)果與原始光譜在譜圖出峰點位置的凸函數(shù)表達不一致，導(dǎo)致平滑后的數(shù)據(jù)存在明顯差異。

圖3 插值光譜平滑求導(dǎo)（1階）預(yù)處理后的譜圖對比Fig.3 Comparison of spectra after interpolation and 1st derivative smoothing

2.1.2 光譜平滑后再插值處理的譜圖分析

利用Zero、Slinear、Quadratic、Cubic、Nearest 5種插值方法對經(jīng)過SG平滑后的子機光譜進行插值，同時利用本研究中提出的基于Savitzky-Golay平滑的插值方法（SG-Inter）對光譜進行處理后與原始光譜平滑后對比得到圖4?？芍交蠊庾V數(shù)據(jù)插值與原始光譜平滑后的數(shù)據(jù)能夠基本重合。說明平滑后再插值可以解決子機光譜與主機光譜波數(shù)不一致而導(dǎo)致的模型轉(zhuǎn)移問題。進一步將出峰點光譜局部放大后，發(fā)現(xiàn)經(jīng)過Cubic、Quadratic及SG-Inter插值方法處理后的數(shù)據(jù)能夠較好地與原始光譜相重合，而其他幾類方法如Zero、Nearest等均存在較為明顯的波動。

圖4 平滑求導(dǎo)（1階）后再插值的光譜對比Fig.4 Comparison of spectra after 1st derivative smoothing and interpolation

由于光譜的平滑過程已經(jīng)對原始數(shù)據(jù)進行降噪處理，結(jié)合譜圖變化，從光譜形態(tài)變化的角度分析后不難發(fā)現(xiàn)，插值方法能夠在此基礎(chǔ)上充分獲取數(shù)據(jù)點前后的變化趨勢，更易于尋找插值函數(shù)多項式使子機光譜與主機光譜的波數(shù)調(diào)整一致。因此，先平滑再插值的效果要優(yōu)于先插值再平滑的效果，并有效支撐模型轉(zhuǎn)移和化學(xué)指標定量預(yù)測。

2.2 插值方法對模型轉(zhuǎn)移預(yù)測結(jié)果的影響

使用不同方法插值后的光譜進行模型轉(zhuǎn)移及指標預(yù)測，研究插值方法對指標預(yù)測的定量影響。為了進一步對子機光譜的預(yù)測結(jié)果進行統(tǒng)計分析，選擇各指標的均方根誤差（Root mean square error，RMSE）、模型決定系數(shù)（R2）、相對分析誤差［23］（Residual predictive deviation，RPD）作為評價指標，分別對模型轉(zhuǎn)移效果進行分析比較。

2.2.1 基于原始光譜直接插值的模型轉(zhuǎn)移預(yù)測結(jié)果分析

利用Zero、Slinear、Quadratic、Cubic、Nearest 5種插值方法對子機原始光譜直接插值，通過對插值后的光譜進行空間變換方法（SST）［6］實現(xiàn)光譜轉(zhuǎn)移，將各指標預(yù)測結(jié)果與主機模型預(yù)測結(jié)果統(tǒng)計分析后得到表1?？芍?，基于原始光譜直接插值的模型轉(zhuǎn)移及指標預(yù)測結(jié)果的RMSE 遠大于主機模型交叉驗證結(jié)果。說明模型的預(yù)測精度距離主機模型還有一定差距。同樣，基于光譜直接插值的轉(zhuǎn)移模型的R2普遍較小，說明模型與預(yù)測指標之間的關(guān)聯(lián)關(guān)系還不夠緊密。此外，通常認為，若RPD＜1.4，表明所建模型預(yù)測結(jié)果不可靠；若1.4≤RPD≤2.0，表明模型的預(yù)測結(jié)果可以接受；若RPD＞2.0，則表明模型的預(yù)測準確性很高，能夠用于模型分析［23］。從表1 可知，整體上，基于Cubic、Quadratic 等非線性插值方法處理后的各指標預(yù)測結(jié)果比其他插值后的預(yù)測效果好，模型的RPD 均大于1.4，在可接受范圍內(nèi)，但相對于主機模型的預(yù)測性能還有較大的差距。

表1 原始光譜直接插值下9種煙草化學(xué)指標的RMSE、R2、RPD值Tab.1 RMSE，R2 and RPD values of nine tobacco chemical indexes after interpolation

就總植物堿、總糖、總氮、還原糖、氯、鉀、淀粉、新植二烯和綠原酸的質(zhì)量分數(shù)而言，綠原酸的模型預(yù)測能力最差，且預(yù)測誤差最大。對于模型轉(zhuǎn)移方法，基于子機原始光譜直接插值的不同方法整體預(yù)測效果均比較差，轉(zhuǎn)移模型決定系數(shù)均較小，預(yù)測誤差相對較大，尤其是基于Slinear、Quadratic、Cubic和Nearest插值方法處理的光譜，模型轉(zhuǎn)移預(yù)測能力均處于同一水平；經(jīng)過Zero 插值處理的光譜模型轉(zhuǎn)移預(yù)測效果最差。

2.2.2 基于光譜平滑數(shù)據(jù)插值的模型轉(zhuǎn)移預(yù)測結(jié)果

將子機原始光譜SG平滑預(yù)處理后，利用Zero、Slinear、Cubic、Nearest和Quadratic方法以及本研究中提出的基于Savitzky-Golay的平滑插值（SG-Inter）方法對數(shù)據(jù)進行變換，各化學(xué)指標定量預(yù)測的統(tǒng)計分析結(jié)果見表2。

表2 光譜平滑后再插值下的RMSE、R2、RPD值Tab.2 RMSE，R2 and RPD values after smoothing and interpolation

可知，對于總糖和氯來說，SG-Inter插值處理的光譜預(yù)測誤差明顯優(yōu)于其他插值方法；對于總植物堿、總氮、鉀、還原糖、淀粉和新植二烯來說，無論是基于平滑后再Cubic、Quadratic等插值方法還是平滑插值（SG-Inter）的模型轉(zhuǎn)移預(yù)測誤差RMSE均較小且差別不大，但基于Zero插值方法的上述9種化學(xué)指標模型轉(zhuǎn)移預(yù)測結(jié)果相對不穩(wěn)定且誤差較大；對于多酚類物質(zhì)綠原酸，所建主機模型決定系數(shù)不是很高，對不同插值方法處理后的光譜進行模型轉(zhuǎn)移，化學(xué)指標預(yù)測結(jié)果誤差較大，且模型不夠穩(wěn)健，后續(xù)仍需要對主機模型進行優(yōu)化，對子機模型轉(zhuǎn)移進行更深入的研究。綜合來看，盡管與主機模型相比，這9種化學(xué)指標的子機模型預(yù)測RMSE均大于主機模型交叉驗證的RMSE，且外部驗證誤差稍大于主機模型內(nèi)部交叉驗證誤差還可以被理解，但仍然有提升的空間；而對于總植物堿、總糖、總氮、還原糖和新植二烯來說，SG-Inter插值和Nearest、Cubic、Quadratic等插值方法處理后的光譜整體上模型預(yù)測能力與主機模型相當且差距不大，且均較為良好；但對于氯、鉀、綠原酸和淀粉來說，與主機模型預(yù)測結(jié)果的R2相比仍有一定差距。

綜上，針對本研究中所實測的200個代表性樣品，基于所提出的SG-Inter插值方法處理子機光譜數(shù)據(jù)，在經(jīng)過光譜空間變換（SST）實現(xiàn)子機光譜轉(zhuǎn)移后，對9種指標進行主機模型預(yù)測，結(jié)果表明：除綠原酸外，其余8種指標的模型轉(zhuǎn)移預(yù)測結(jié)果的RMSE均相對最小，轉(zhuǎn)移模型的R2與RPD在多種指標中的評價性能均良好、略優(yōu)或與其他插值方法效果相當。表明SG-Inter方法可以順利支撐模型轉(zhuǎn)移及預(yù)測工作的開展，而且該方法將光譜平滑預(yù)處理過程與插值計算結(jié)合在一起，在保證模型預(yù)測能力和精度的前提下，能夠大大提升模型轉(zhuǎn)移的工作效率。而對于綠原酸等其他指標的模型轉(zhuǎn)移工作需要進一步從近紅外光譜形成機理著手分析數(shù)據(jù)，減少SG-Inter方法處理過程中造成的光譜精度損失，從而提升指標預(yù)測的準確度。

3 結(jié)論

①使用插值的方法能夠有效解決子機光譜與主機光譜波數(shù)及波長點不一致而導(dǎo)致模型轉(zhuǎn)移中主機模型失效的問題。②子機光譜直接插值后與原始光譜能相對重合，但平滑處理后的光譜差異顯著，且模型轉(zhuǎn)移預(yù)測結(jié)果誤差較大；子機光譜SG平滑處理后再插值的光譜與原始光譜能相對重合，且模型轉(zhuǎn)移預(yù)測效果良好。③基于Savitzky-Golay平滑插值（SG-Inter）的模型轉(zhuǎn)移對多種指標質(zhì)量分數(shù)的預(yù)測效果略優(yōu)于Cubic、Quadratic、Slinear等5種插值方法，指標預(yù)測誤差相對較小，模型預(yù)測能力較好，且計算流程清晰，原理簡單，能夠?qū)⑵交幚砼c插值分析集為一體并同步高效運算，可減少預(yù)處理環(huán)節(jié)的分析誤差，提高模型轉(zhuǎn)移的工作效率。