基于尺度不變特征變換篩選穩定特征波長的近紅外光譜模型傳遞方法

洪士軍，黃 雯，張立國，葛 炯*，倪力軍*，欒紹嶸

(1.華東理工大學 化學與分子工程學院，上海 200237;2.上海煙草集團有限責任公司 技術中心理化實驗室，上海 200082)

近紅外光譜(Near infrared spectra,NIRS)是一種簡單、快速的新型分析方法[1]。20世紀80年代后期以來，隨著NIRS儀器的不斷改進、化學計量學方法的應用和計算機的快速發展，NIRS技術被廣泛應用于各個領域[2]。建立一個良好的NIRS模型需要足夠多的樣品數據并進行方法和模型參數優化，模型建立、維護需要人力、物力的持續投入，通常希望所建立的NIRS模型可以在不同儀器之間共享[3]。但由于各儀器的使用環境、壽命及狀態不同，相同測試條件下，同一樣品在同型號的不同儀器上的光譜響應并非完全相同，導致主儀器上所建立的NIRS模型對從儀器樣品的預測誤差可能超出許可范圍，此時需要采用各種模型傳遞方法對從儀器光譜或主機模型進行修正[4]。本課題組提出的根據儀器間光譜差異的方差分析[5-6]、光譜比值分析[7]篩選儀器間信號一致且穩定的波長，基于這些波長的光譜信息在主機建立偏最小二乘(PLS)近紅外光譜模型，并直接用于從機玉米、黃芩樣品中主要成分含量預測的方法，對從機樣品的預測誤差與分段直接校正算法(Piecewise direct standardization，PDS)相當或更低。進一步在一致且穩定波長的基礎上，綜合應用相關系數、無信息變量消除(Uninformative variable elimination，UVE)[8]等波長篩選方法優化波長組合后建立煙葉總植物堿NIRS模型，該模型直接預測從機煙葉樣品總植物堿的誤差滿足企業內控要求[9]。文獻[10]通過分析主、從機光譜在不同波長下的相關性，篩選儀器間光譜響應一致性好的波長，建立了玉米中主要成分含量的NIRS模型，對從機樣品的預測誤差與主機相當。上述研究表明基于儀器間光譜響應一致性好的波長建立NIRS模型，有助于提高模型的穩健性，使模型可在不同儀器間直接共享，而無需在模型傳遞過程進行從機光譜校正或模型校正，比傳統的有標樣模型傳遞方法簡便。但上述文獻報道的波長篩選過程仍需若干主、從機樣品光譜作為波長篩選的信息來源，不能稱之為完全意義的無標樣模型傳遞。

尺度不變特征變換(Scale invariant features transform,SIFT)算法是一種用來提取圖像局部性特征的算法，該方法可在不同尺度空間查找圖像中一些與其周圍區域有高區分度的局部圖像特征，這些特征(又稱關鍵點)具有重復性好、區分度和準確度高的特點，并且在不同尺度空間保持不變，不隨圖像旋轉、亮度變化而改變，對視角變化、仿射變換、噪聲也保持一定程度的穩定性。通過提取這些特征的位置、尺度、旋轉不變量，可實現對圖像的精準識別[11-12]。本文以煙葉總植物堿近紅外光譜模型的建立和轉移為背景，利用SIFT算法提取主機樣品一維近紅外光譜在不同尺度下的關鍵點，這些點即為主、從機光譜中的穩定特征波長。以這些波長的光譜信息建立近紅外光譜模型，可剔除噪聲和無效信息對模型的干擾，獲得穩健的近紅外光譜模型，有望實現模型在不同儀器間的共享或直接轉移。

1 算法與原理

1.1 一維SIFT算法篩選穩定特征波長

一維SIFT算法是二維SIFT算法的特例，其核心是構建不同尺度空間，尋找在不同尺度下均保持穩定不變的關鍵點。主要步驟如下：

(1)構建尺度空間

尺度空間包含兩個尺度概念，一個尺度是對光譜圖像的采樣間隔，不同采樣間隔構成了圖像抽提的尺度；另一個尺度是對光譜信號進行不同尺度變換的尺度。高斯卷積核是實現尺度變換的唯一線性核，一幅一維光譜圖像在不同尺度下的變換由光譜與高斯卷積核卷積獲得：

L(x,σ)=G(x,σ)*I(x)

(1)

上式中I(x)為近紅外光譜數據，x(光譜波長)是空間坐標，x坐標的間隔構成不同尺度空間的光譜I(x)。G(x,σ)是尺度可變高斯函數，又稱為高斯卷積核，其定義如下：

(2)

σ為尺度空間變換因子，是高斯正態分布的方差，其大小決定圖像的平滑尺度。大的σ對應粗糙尺度(低分辨率)圖像的概貌特征，小的σ對應精細尺度(高分辨率)圖像的細節特征。圖1為二維、一維圖像在不同空間尺度(采樣間隔)和不同尺度變換因子σ下得到的高斯卷積示意圖。由該圖可知，對于二維圖像，隨著采樣步長(x的間隔)的增大，尺度空間L(x,σ)變小，不同尺度空間L(x,σ)構成金字塔形狀，稱為高斯金字塔。對于一維圖像，高斯金字塔退化為一個三角形。每個尺度空間L(x,σ)由若干長度相同的層構成不同組Octave1、Octave2…，采樣間隔(x的間隔)不同產生長度不同的組，高斯卷積的尺度變換因子σ大小不同產生同組中的不同層。高斯金字塔的組數O、每組層數S及尺度變換因子σ決定了高斯金字塔的結構。

圖1 二維圖像(A)以及一維圖像(B)的高斯金字塔示意圖Fig.1 Schematic diagram of Gaussian pyramid of two-dimensional image(A) and one-dimensional image(B)

(2)進行尺度空間極值檢測

搜索所有尺度上的光譜點，通過高斯差分(DOG)尺度空間來識別潛在的對尺度和旋轉不變的關鍵點。DOG尺度空間由不同尺度的高斯差分核與光譜圖像卷積生成：

D(x,σ)=(G(x,kσ)-G(x,σ))*I(x)=L(x,kσ)-L(x,σ)

(3)

式(3)中k為相鄰層之間尺度因子大小的比例。第一組第一層的σ稱為初始σ，記為σ0。

(3)高斯金字塔結構參數的優化

采用一維SIFT算法可根據一根光譜確定在光譜特定金字塔結構下的關鍵點集合。基于同一個樣品連續p次測量的精密度光譜數據，可得到p個關鍵點集合。由于儀器噪聲和測試誤差等因素影響，同一樣品連續p次測試的光譜不可能完全重合。因此，這p個集合中的關鍵波長點也不完全重合。本文以重現率(Percent recurrence,PR)和重現度(Reproducibility，RPD)衡量SIFT算法篩選出的波長的穩定性，以PR和RPD最大為目標，采用3因素3水平正交表優化高斯金字塔的結構參數。SIFT篩選波長重現率和重現度的計算公式如下：

PR=(n0×w0+n1×w1+n2×w2)/n

(4)

RPD=(n0×w0+n1×w1+n2×w2)/p

(5)

n0代表p次篩選出來的波長點相同的波長個數，n1代表p次篩選出來的波長相差一個點數的波長個數，n2代表p次篩選出來的波長差異在兩個點數的波長個數，n代表根據p條精密度測試光譜篩選出的波長總數。w0、w1和w2分別為波長點重合度不同情況下的權重因子，對于p次全部重合的波長，其個數n0所乘的權重因子w0應該最大，本文設w0為3；基于該原則，對于波長相差一個點的波長個數n1，其所乘的權重因子w1定義為2，波長相差2個點的波長個數n2所乘的權重因子w2設為1。

1.2 近紅外光譜模型的建立與評價

以所篩選的穩定特征波長下的光譜信息為自變量，采用偏最小二乘(PLS)算法建立待測性質的定量模型。以RMSEC(校正集的均方根殘差)評價模型的擬合性能，RMSEP(驗證集的均方根殘差)和MRE(相對誤差絕對值的平均值)來評估模型的預測性能及傳遞性能。

(6)

(7)

式(6)～(7)中的yi,actual為第i個樣品實測值，yi,predicted為第i個樣品的模型預測值，m為樣品數目，RMSEP與RMSEC計算公式相同，但在計算時將公式(6)中校正集樣本替換成預測集或驗證集樣本。RMSEC越小，則認為模型的擬合性能越好；RMSEP及MRE越接近于0，則認為模型的預測性能越好。煙草企業內控指標要求煙葉總植物堿的MRE小于6%。

本文所有算法在MATLAB平臺完成和編譯。

2 實驗部分

本文共使用3套煙葉樣本數據集，Set A由78個煙葉樣本分別在4臺AntarisⅡ近紅外儀器(賽默飛世爾科技有限公司)上測得的光譜組成，4臺儀器分別為主機M(Master)、3臺從機S1(Slave 1)、S2(Slave 2)和S3(Slave 3)；Set B由279個在主機M上測得的煙葉樣本光譜組成。Set A、Set B中各煙葉樣品的總植物堿采用YC/T 160-2002[14]測定，其含量在0.55%～6.30%之間。某一煙葉樣品在主機M上連續測量3次的精密度近紅外光譜記為Set C。

煙葉經烘箱50 ℃烘干后粉碎，過40～60 目篩。置于可旋轉石英杯中，在(22±4)℃、30%～80%濕度下測試其光譜：分辨率取8 cm-1，波數范圍取4 000～10 000 cm-1，掃描64次，增益取2。

3 結果與討論

3.1 SIFT算法中的參數優化

取組數O的水平為3、4、5，層數S的水平為5、6、7，尺度因子σ0水平取1.5、2.5、3.5。以Set C中的3次精密度測試光譜為信息來源，基于正交表L9(33)篩選出的3個波長集合的重現率與重現度如表1所示。

表1 根據正交表L9(33)篩選出的波長重現率及重現度Table 1 Wavelength reproducibility rate and reproducibility based on orthogonalTable L9(33)

由表1可知，重現度RPD與重現率PR呈正相關趨勢，故只對重現率PR進行方差分析，結果如表2所示。

表2 表1中重現率的方差分析結果Table 2 Variance analysis results of reproducibility rate inTable 1

由表2可知，組數和層數的P值遠大于0.05，說明組數O和層數S的取值對重現率、重現度與結果無顯著影響。而高斯卷積初始尺度因子σ0的P值遠小于0.05，表明σ0值的變化對波長篩選結果有顯著影響。由表1可知，組數O的最優水平為5，層數S的最優水平為5，σ0的最優水平為3.5，最優的各因素水平組合為第7號實驗，對應的波長重現率為97%。

表3 O=4，S=5時，不同σ0下SIFT方法所 篩選波長的重現率與重現度Table 3 Reproducibility rate and reproducibility of wavelengths selected by SIFT with different σ0 when O=4,S=5

由于組數和層數的影響不顯著，為減少運算量，將高斯組數設為4，層數設為5，探索σ0不同取值下的重現率變化情況，如表3所示。

由表3可知，σ0取2.7時重現率可以達到94%，σ0的繼續增大對于重現率的提升貢獻有限，而重現度隨著σ0的增大有較為明顯的提升。根據式(7)可知，重現度反映的是根據3次精密度測試光譜篩選出的穩定特征波長重合的波長個數的均值，其隨σ0增大而增大。根據本課題組前期研究[9]，在篩選的波長個數為80～120范圍內，基于這些波長所建煙葉總植物堿NIRS模型的MRE在5%左右，滿足企業內控要求。σ0取2.7時，可篩選出103個波長。故本文取高斯金字塔組數為4，層數為5，初始σ0取2.7。

3.2 煙葉樣品光譜及SIFT方法所篩選的穩定特征波長

采用Set B數據建立主機模型，以Set A中主機M上測試的78個煙葉樣品為主機外部驗證集，Set A中分別在S1、S2、S3上測試的78個樣品為從機檢驗集。首先將所有的光譜數據經過標準正態變換(SNV)+一階導數預處理，再從Set B的279個建模樣品中利用基于x-y空間距離劃分樣本的SPXY(Sample set partitioning based on jointx-ydistance)方法[15]分別挑選差異最大的前5、10、15、20、30個代表性樣品。采用“3.1”得到的SIFT方法中的優化參數，根據這5(10、15、20、30)個代表性樣品的近紅外光譜分別篩選得到5(10、15、20、30)個穩定特征波長集合，每個集合的波長點數在80～90之間。取這些集合的交集，僅能得到40個左右的波長點，以這幾十個波長的光譜信息建立煙葉總植物堿NIRS模型時，由于模型可調參數太少，其MRE超過6%，不滿足企業內控要求。故本文分別取這5(10、15、20、30)個波長集合的并集作為SIFT算法篩選的穩定特征波長點，以此建立的煙葉總植物堿NIRS模型記為SIFT-PLS模型。除根據5個樣品光譜篩選的波長所建立的SIFT-PLS模型不能保證對3臺從機樣品總植物堿的MRE小于6%之外，取10、15、20及30個樣品光譜篩選波長所建的SIFT-PLS模型預測主、從機樣品總植物堿的MRE均小于6%。限于篇幅，本文只給出10個代表性樣品光譜所篩選的波長集合并集的結果，最終得到304個穩定特征波長點。篩選出的波長及由這些波長繪制的4臺近紅外儀器的平均光譜如圖2B所示。由該圖可見，根據這304個特征波長點的光譜響應所繪的圖雖然不是很光滑，但保留了原樣品光譜(圖2A)的基本圖像特征，說明SIFT方法篩選的波長確實提取出了煙葉樣品的主要光譜特征。

由圖3可知，SIFT方法篩選出的大部分特征波長點在儀器間光譜差異很小的區域，那么以這304個波長的光譜響應為自變量，建立的煙葉總植物堿的SIFT-PLS光譜模型對主、從機樣品預測結果應有良好的一致性。

圖3 主機M與從機S2的SNV+一階導數光譜差譜的 絕對值及SIFT篩選波長Fig.3 The absolute difference spectrum between instrument M and S2 and the wavelengths selected by SIFT algorithm

3.3 不同煙葉總植物堿近紅外光譜模型結果的比較

表4給出了SIFT-PLS模型、全波長PLS(Whole wavelength PLS，簡稱WW-PLS)模型以及WW-PLS模型與SIFT-PLS模型經PDS校正從機光譜后模型(PDS-WW-PLS、PDS-SIFT-PLS)的結果。從表4可以看出，WW-PLS在主機M以及從機S2上的預測誤差滿足企業內控要求(MRE<6%)，其對另外2臺從機樣品的MRE均大于6%。經PDS校正從機光譜后，PDS-WW-PLS模型對S1、S3從機樣品的MRE分別從9.13%、8.18%降為5.30%與6.12%，而對S2的MRE反而從3.85%增大到8.13%；SIFT-PLS模型預測S2從機樣品總植物堿的MRE經PDS校正后從4.24%增大到4.88%。說明NIRS模型直接轉移到從機的誤差滿足要求時，采用PDS校正反而會增大模型傳遞的誤差；在模型傳遞誤差較高時，PDS能起到一定改善作用，但也不能保證WW-PLS對3臺從機的傳遞誤差均滿足企業內控要求。而SIFT-PLS模型對3臺從機樣品的MRE均滿足小于6%的企業內控要求，當采用PDS校正從機光譜后，該模型對3臺從機樣品總植物堿的RMSEP均較SIFT-PLS模型直接轉移有所降低，且PDS-SIFT-PLS模型對從機樣品的MRE均低于PDS-WW-PLS模型。說明SIFT算法篩選波長所建模型較全波長模型更為穩健、精簡，且PDS校正從機光譜后仍能保證SIFT-PLS模型對從機樣品的MRE≤5%。

表4 不同煙葉總植物堿模型預測結果的比較Table 4 Results of total alkaloids contents predicted by different calibration models

4 結 論

利用SIFT算法篩選近紅外光譜的關鍵點，波長篩選過程中無需從機樣品信息，篩選出的波長點保留了光譜的主要特征，其光譜響應不隨光譜尺度變化而改變，是樣品光譜的穩定特征點，以這些波長下的光譜信息為自變量建立的煙葉總植物堿近紅外光譜模型具有穩健、移植性好的特點。該模型直接轉移到3臺從機，無需任何模型傳遞方法進行從機光譜或模型校正，對從機樣品總植物堿的MRE為4.24%～5.81%，經PDS算法校正從機光譜后，PDS-SIFT-PLS模型預測從機樣品總植物堿的MRE為4.48～5.00%，均能滿足企業內控要求；而PDS-WW-PLS模型預測兩臺從機樣品總植物堿的MRE大于6%，表明全波長模型的穩健性和移植性均不及SIFT-PLS模型。

基于SIFT方法篩選近紅外光譜的穩定特征波長建立煙葉總植物堿的近紅外光譜模型，可實現模型的無標樣傳遞，簡便易行。該方法用于其他領域及樣品的模型傳遞有待進一步驗證。