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

(1.華東理工大學(xué) 化學(xué)與分子工程學(xué)院，上海 200237;2.上海煙草集團(tuán)有限責(zé)任公司 技術(shù)中心理化實(shí)驗(yàn)室，上海 200082)

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

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

1 算法與原理

1.1 一維SIFT算法篩選穩(wěn)定特征波長(zhǎng)

一維SIFT算法是二維SIFT算法的特例，其核心是構(gòu)建不同尺度空間，尋找在不同尺度下均保持穩(wěn)定不變的關(guān)鍵點(diǎn)。主要步驟如下：

(1)構(gòu)建尺度空間

尺度空間包含兩個(gè)尺度概念，一個(gè)尺度是對(duì)光譜圖像的采樣間隔，不同采樣間隔構(gòu)成了圖像抽提的尺度；另一個(gè)尺度是對(duì)光譜信號(hào)進(jìn)行不同尺度變換的尺度。高斯卷積核是實(shí)現(xiàn)尺度變換的唯一線性核，一幅一維光譜圖像在不同尺度下的變換由光譜與高斯卷積核卷積獲得：

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

(1)

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

(2)

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

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

(2)進(jìn)行尺度空間極值檢測(cè)

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

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

(3)

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

(3)高斯金字塔結(jié)構(gòu)參數(shù)的優(yōu)化

采用一維SIFT算法可根據(jù)一根光譜確定在光譜特定金字塔結(jié)構(gòu)下的關(guān)鍵點(diǎn)集合?；谕粋€(gè)樣品連續(xù)p次測(cè)量的精密度光譜數(shù)據(jù)，可得到p個(gè)關(guān)鍵點(diǎn)集合。由于儀器噪聲和測(cè)試誤差等因素影響，同一樣品連續(xù)p次測(cè)試的光譜不可能完全重合。因此，這p個(gè)集合中的關(guān)鍵波長(zhǎng)點(diǎn)也不完全重合。本文以重現(xiàn)率(Percent recurrence,PR)和重現(xiàn)度(Reproducibility，RPD)衡量SIFT算法篩選出的波長(zhǎng)的穩(wěn)定性，以PR和RPD最大為目標(biāo)，采用3因素3水平正交表優(yōu)化高斯金字塔的結(jié)構(gòu)參數(shù)。SIFT篩選波長(zhǎng)重現(xiàn)率和重現(xiàn)度的計(jì)算公式如下：

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

(4)

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

(5)

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

1.2 近紅外光譜模型的建立與評(píng)價(jià)

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

(6)

(7)

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

本文所有算法在MATLAB平臺(tái)完成和編譯。

2 實(shí)驗(yàn)部分

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

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

3 結(jié)果與討論

3.1 SIFT算法中的參數(shù)優(yōu)化

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

表1 根據(jù)正交表L9(33)篩選出的波長(zhǎng)重現(xiàn)率及重現(xiàn)度Table 1 Wavelength reproducibility rate and reproducibility based on orthogonalTable L9(33)

由表1可知，重現(xiàn)度RPD與重現(xiàn)率PR呈正相關(guān)趨勢(shì)，故只對(duì)重現(xiàn)率PR進(jìn)行方差分析，結(jié)果如表2所示。

表2 表1中重現(xiàn)率的方差分析結(jié)果Table 2 Variance analysis results of reproducibility rate inTable 1

由表2可知，組數(shù)和層數(shù)的P值遠(yuǎn)大于0.05，說明組數(shù)O和層數(shù)S的取值對(duì)重現(xiàn)率、重現(xiàn)度與結(jié)果無顯著影響。而高斯卷積初始尺度因子σ0的P值遠(yuǎn)小于0.05，表明σ0值的變化對(duì)波長(zhǎng)篩選結(jié)果有顯著影響。由表1可知，組數(shù)O的最優(yōu)水平為5，層數(shù)S的最優(yōu)水平為5，σ0的最優(yōu)水平為3.5，最優(yōu)的各因素水平組合為第7號(hào)實(shí)驗(yàn)，對(duì)應(yīng)的波長(zhǎng)重現(xiàn)率為97%。

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

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

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

3.2 煙葉樣品光譜及SIFT方法所篩選的穩(wěn)定特征波長(zhǎng)

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

由圖3可知，SIFT方法篩選出的大部分特征波長(zhǎng)點(diǎn)在儀器間光譜差異很小的區(qū)域，那么以這304個(gè)波長(zhǎng)的光譜響應(yīng)為自變量，建立的煙葉總植物堿的SIFT-PLS光譜模型對(duì)主、從機(jī)樣品預(yù)測(cè)結(jié)果應(yīng)有良好的一致性。

圖3 主機(jī)M與從機(jī)S2的SNV+一階導(dǎo)數(shù)光譜差譜的 絕對(duì)值及SIFT篩選波長(zhǎng)Fig.3 The absolute difference spectrum between instrument M and S2 and the wavelengths selected by SIFT algorithm

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

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

表4 不同煙葉總植物堿模型預(yù)測(cè)結(jié)果的比較Table 4 Results of total alkaloids contents predicted by different calibration models

4 結(jié) 論

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

基于SIFT方法篩選近紅外光譜的穩(wěn)定特征波長(zhǎng)建立煙葉總植物堿的近紅外光譜模型，可實(shí)現(xiàn)模型的無標(biāo)樣傳遞，簡(jiǎn)便易行。該方法用于其他領(lǐng)域及樣品的模型傳遞有待進(jìn)一步驗(yàn)證。