基于Savitzky-Golay平滑插值的煙草近紅外光譜模型轉移

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

中國煙草總公司鄭州煙草研究院，鄭州高新技術產業開發區楓楊街2號 450001

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

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

1 材料與方法

1.1 材料、試劑和儀器

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

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

1.2 方法

1.2.1 樣品處理與化學指標檢測

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

1.2.2 光譜采集與預處理

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

1.2.3 主機模型

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

1.2.4 插值方法的研究與應用

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

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

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

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

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

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

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

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

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

1.2.6 光譜轉移及化學成分預測

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

1.2.7 數據處理方法

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

2 結果與分析

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

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

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

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

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

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

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

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

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

由于光譜的平滑過程已經對原始數據進行降噪處理，結合譜圖變化，從光譜形態變化的角度分析后不難發現，插值方法能夠在此基礎上充分獲取數據點前后的變化趨勢，更易于尋找插值函數多項式使子機光譜與主機光譜的波數調整一致。因此，先平滑再插值的效果要優于先插值再平滑的效果，并有效支撐模型轉移和化學指標定量預測。

2.2 插值方法對模型轉移預測結果的影響

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

2.2.1 基于原始光譜直接插值的模型轉移預測結果分析

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

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

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

2.2.2 基于光譜平滑數據插值的模型轉移預測結果

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

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

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

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

3 結論

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