洪士軍,黃 雯,張立國,葛 炯*,倪力軍*,欒紹嶸
(1.華東理工大學(xué) 化學(xué)與分子工程學(xué)院,上海 200237;2.上海煙草集團有限責(zé)任公司 技術(shù)中心理化實驗室,上海 200082)
近紅外光譜(Near infrared spectra,NIRS)是一種簡單、快速的新型分析方法[1]。20世紀(jì)80年代后期以來,隨著NIRS儀器的不斷改進、化學(xué)計量學(xué)方法的應(yīng)用和計算機的快速發(fā)展,NIRS技術(shù)被廣泛應(yīng)用于各個領(lǐng)域[2]。建立一個良好的NIRS模型需要足夠多的樣品數(shù)據(jù)并進行方法和模型參數(shù)優(yōu)化,模型建立、維護需要人力、物力的持續(xù)投入,通常希望所建立的NIRS模型可以在不同儀器之間共享[3]。但由于各儀器的使用環(huán)境、壽命及狀態(tài)不同,相同測試條件下,同一樣品在同型號的不同儀器上的光譜響應(yīng)并非完全相同,導(dǎo)致主儀器上所建立的NIRS模型對從儀器樣品的預(yù)測誤差可能超出許可范圍,此時需要采用各種模型傳遞方法對從儀器光譜或主機模型進行修正[4]。本課題組提出的根據(jù)儀器間光譜差異的方差分析[5-6]、光譜比值分析[7]篩選儀器間信號一致且穩(wěn)定的波長,基于這些波長的光譜信息在主機建立偏最小二乘(PLS)近紅外光譜模型,并直接用于從機玉米、黃芩樣品中主要成分含量預(yù)測的方法,對從機樣品的預(yù)測誤差與分段直接校正算法(Piecewise direct standardization,PDS)相當(dāng)或更低。進一步在一致且穩(wěn)定波長的基礎(chǔ)上,綜合應(yīng)用相關(guān)系數(shù)、無信息變量消除(Uninformative variable elimination,UVE)[8]等波長篩選方法優(yōu)化波長組合后建立煙葉總植物堿NIRS模型,該模型直接預(yù)測從機煙葉樣品總植物堿的誤差滿足企業(yè)內(nèi)控要求[9]。文獻[10]通過分析主、從機光譜在不同波長下的相關(guān)性,篩選儀器間光譜響應(yīng)一致性好的波長,建立了玉米中主要成分含量的NIRS模型,對從機樣品的預(yù)測誤差與主機相當(dāng)。上述研究表明基于儀器間光譜響應(yīng)一致性好的波長建立NIRS模型,有助于提高模型的穩(wěn)健性,使模型可在不同儀器間直接共享,而無需在模型傳遞過程進行從機光譜校正或模型校正,比傳統(tǒng)的有標(biāo)樣模型傳遞方法簡便。但上述文獻報道的波長篩選過程仍需若干主、從機樣品光譜作為波長篩選的信息來源,不能稱之為完全意義的無標(biāo)樣模型傳遞。
尺度不變特征變換(Scale invariant features transform,SIFT)算法是一種用來提取圖像局部性特征的算法,該方法可在不同尺度空間查找圖像中一些與其周圍區(qū)域有高區(qū)分度的局部圖像特征,這些特征(又稱關(guān)鍵點)具有重復(fù)性好、區(qū)分度和準(zhǔn)確度高的特點,并且在不同尺度空間保持不變,不隨圖像旋轉(zhuǎn)、亮度變化而改變,對視角變化、仿射變換、噪聲也保持一定程度的穩(wěn)定性。通過提取這些特征的位置、尺度、旋轉(zhuǎn)不變量,可實現(xiàn)對圖像的精準(zhǔn)識別[11-12]。本文以煙葉總植物堿近紅外光譜模型的建立和轉(zhuǎn)移為背景,利用SIFT算法提取主機樣品一維近紅外光譜在不同尺度下的關(guān)鍵點,這些點即為主、從機光譜中的穩(wěn)定特征波長。以這些波長的光譜信息建立近紅外光譜模型,可剔除噪聲和無效信息對模型的干擾,獲得穩(wěn)健的近紅外光譜模型,有望實現(xiàn)模型在不同儀器間的共享或直接轉(zhuǎn)移。
一維SIFT算法是二維SIFT算法的特例,其核心是構(gòu)建不同尺度空間,尋找在不同尺度下均保持穩(wěn)定不變的關(guān)鍵點。主要步驟如下:
(1)構(gòu)建尺度空間
尺度空間包含兩個尺度概念,一個尺度是對光譜圖像的采樣間隔,不同采樣間隔構(gòu)成了圖像抽提的尺度;另一個尺度是對光譜信號進行不同尺度變換的尺度。高斯卷積核是實現(xiàn)尺度變換的唯一線性核,一幅一維光譜圖像在不同尺度下的變換由光譜與高斯卷積核卷積獲得:
L(x,σ)=G(x,σ)*I(x)
(1)
上式中I(x)為近紅外光譜數(shù)據(jù),x(光譜波長)是空間坐標(biāo),x坐標(biāo)的間隔構(gòu)成不同尺度空間的光譜I(x)。G(x,σ)是尺度可變高斯函數(shù),又稱為高斯卷積核,其定義如下:
(2)
σ為尺度空間變換因子,是高斯正態(tài)分布的方差,其大小決定圖像的平滑尺度。大的σ對應(yīng)粗糙尺度(低分辨率)圖像的概貌特征,小的σ對應(yīng)精細尺度(高分辨率)圖像的細節(jié)特征。圖1為二維、一維圖像在不同空間尺度(采樣間隔)和不同尺度變換因子σ下得到的高斯卷積示意圖。由該圖可知,對于二維圖像,隨著采樣步長(x的間隔)的增大,尺度空間L(x,σ)變小,不同尺度空間L(x,σ)構(gòu)成金字塔形狀,稱為高斯金字塔。對于一維圖像,高斯金字塔退化為一個三角形。每個尺度空間L(x,σ)由若干長度相同的層構(gòu)成不同組Octave1、Octave2…,采樣間隔(x的間隔)不同產(chǎn)生長度不同的組,高斯卷積的尺度變換因子σ大小不同產(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)進行尺度空間極值檢測
搜索所有尺度上的光譜點,通過高斯差分(DOG)尺度空間來識別潛在的對尺度和旋轉(zhuǎn)不變的關(guā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)鍵點集合?;谕粋€樣品連續(xù)p次測量的精密度光譜數(shù)據(jù),可得到p個關(guān)鍵點集合。由于儀器噪聲和測試誤差等因素影響,同一樣品連續(xù)p次測試的光譜不可能完全重合。因此,這p個集合中的關(guān)鍵波長點也不完全重合。本文以重現(xiàn)率(Percent recurrence,PR)和重現(xiàn)度(Reproducibility,RPD)衡量SIFT算法篩選出的波長的穩(wěn)定性,以PR和RPD最大為目標(biāo),采用3因素3水平正交表優(yōu)化高斯金字塔的結(jié)構(gòu)參數(shù)。SIFT篩選波長重現(xiàn)率和重現(xiàn)度的計算公式如下:
PR=(n0×w0+n1×w1+n2×w2)/n
(4)
RPD=(n0×w0+n1×w1+n2×w2)/p
(5)
n0代表p次篩選出來的波長點相同的波長個數(shù),n1代表p次篩選出來的波長相差一個點數(shù)的波長個數(shù),n2代表p次篩選出來的波長差異在兩個點數(shù)的波長個數(shù),n代表根據(jù)p條精密度測試光譜篩選出的波長總數(shù)。w0、w1和w2分別為波長點重合度不同情況下的權(quán)重因子,對于p次全部重合的波長,其個數(shù)n0所乘的權(quán)重因子w0應(yīng)該最大,本文設(shè)w0為3;基于該原則,對于波長相差一個點的波長個數(shù)n1,其所乘的權(quán)重因子w1定義為2,波長相差2個點的波長個數(shù)n2所乘的權(quán)重因子w2設(shè)為1。
以所篩選的穩(wěn)定特征波長下的光譜信息為自變量,采用偏最小二乘(PLS)算法建立待測性質(zhì)的定量模型。以RMSEC(校正集的均方根殘差)評價模型的擬合性能,RMSEP(驗證集的均方根殘差)和MRE(相對誤差絕對值的平均值)來評估模型的預(yù)測性能及傳遞性能。
(6)
(7)
式(6)~(7)中的yi,actual為第i個樣品實測值,yi,predicted為第i個樣品的模型預(yù)測值,m為樣品數(shù)目,RMSEP與RMSEC計算公式相同,但在計算時將公式(6)中校正集樣本替換成預(yù)測集或驗證集樣本。RMSEC越小,則認為模型的擬合性能越好;RMSEP及MRE越接近于0,則認為模型的預(yù)測性能越好。煙草企業(yè)內(nèi)控指標(biāo)要求煙葉總植物堿的MRE小于6%。
本文所有算法在MATLAB平臺完成和編譯。
本文共使用3套煙葉樣本數(shù)據(jù)集,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上連續(xù)測量3次的精密度近紅外光譜記為Set C。
煙葉經(jīng)烘箱50 ℃烘干后粉碎,過40~60 目篩。置于可旋轉(zhuǎn)石英杯中,在(22±4)℃、30%~80%濕度下測試其光譜:分辨率取8 cm-1,波數(shù)范圍取4 000~10 000 cm-1,掃描64次,增益取2。
取組數(shù)O的水平為3、4、5,層數(shù)S的水平為5、6、7,尺度因子σ0水平取1.5、2.5、3.5。以Set C中的3次精密度測試光譜為信息來源,基于正交表L9(33)篩選出的3個波長集合的重現(xiàn)率與重現(xiàn)度如表1所示。
表1 根據(jù)正交表L9(33)篩選出的波長重現(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)趨勢,故只對重現(xiàn)率PR進行方差分析,結(jié)果如表2所示。
表2 表1中重現(xiàn)率的方差分析結(jié)果Table 2 Variance analysis results of reproducibility rate inTable 1
由表2可知,組數(shù)和層數(shù)的P值遠大于0.05,說明組數(shù)O和層數(shù)S的取值對重現(xiàn)率、重現(xiàn)度與結(jié)果無顯著影響。而高斯卷積初始尺度因子σ0的P值遠小于0.05,表明σ0值的變化對波長篩選結(jié)果有顯著影響。由表1可知,組數(shù)O的最優(yōu)水平為5,層數(shù)S的最優(yōu)水平為5,σ0的最優(yōu)水平為3.5,最優(yōu)的各因素水平組合為第7號實驗,對應(yīng)的波長重現(xiàn)率為97%。
表3 O=4,S=5時,不同σ0下SIFT方法所 篩選波長的重現(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ù)的影響不顯著,為減少運算量,將高斯組數(shù)設(shè)為4,層數(shù)設(shè)為5,探索σ0不同取值下的重現(xiàn)率變化情況,如表3所示。
由表3可知,σ0取2.7時重現(xiàn)率可以達到94%,σ0的繼續(xù)增大對于重現(xiàn)率的提升貢獻有限,而重現(xiàn)度隨著σ0的增大有較為明顯的提升。根據(jù)式(7)可知,重現(xiàn)度反映的是根據(jù)3次精密度測試光譜篩選出的穩(wěn)定特征波長重合的波長個數(shù)的均值,其隨σ0增大而增大。根據(jù)本課題組前期研究[9],在篩選的波長個數(shù)為80~120范圍內(nèi),基于這些波長所建煙葉總植物堿NIRS模型的MRE在5%左右,滿足企業(yè)內(nèi)控要求。σ0取2.7時,可篩選出103個波長。故本文取高斯金字塔組數(shù)為4,層數(shù)為5,初始σ0取2.7。
采用Set B數(shù)據(jù)建立主機模型,以Set A中主機M上測試的78個煙葉樣品為主機外部驗證集,Set A中分別在S1、S2、S3上測試的78個樣品為從機檢驗集。首先將所有的光譜數(shù)據(jù)經(jīng)過標(biāo)準(zhǔn)正態(tài)變換(SNV)+一階導(dǎo)數(shù)預(yù)處理,再從Set B的279個建模樣品中利用基于x-y空間距離劃分樣本的SPXY(Sample set partitioning based on jointx-ydistance)方法[15]分別挑選差異最大的前5、10、15、20、30個代表性樣品。采用“3.1”得到的SIFT方法中的優(yōu)化參數(shù),根據(jù)這5(10、15、20、30)個代表性樣品的近紅外光譜分別篩選得到5(10、15、20、30)個穩(wěn)定特征波長集合,每個集合的波長點數(shù)在80~90之間。取這些集合的交集,僅能得到40個左右的波長點,以這幾十個波長的光譜信息建立煙葉總植物堿NIRS模型時,由于模型可調(diào)參數(shù)太少,其MRE超過6%,不滿足企業(yè)內(nèi)控要求。故本文分別取這5(10、15、20、30)個波長集合的并集作為SIFT算法篩選的穩(wěn)定特征波長點,以此建立的煙葉總植物堿NIRS模型記為SIFT-PLS模型。除根據(jù)5個樣品光譜篩選的波長所建立的SIFT-PLS模型不能保證對3臺從機樣品總植物堿的MRE小于6%之外,取10、15、20及30個樣品光譜篩選波長所建的SIFT-PLS模型預(yù)測主、從機樣品總植物堿的MRE均小于6%。限于篇幅,本文只給出10個代表性樣品光譜所篩選的波長集合并集的結(jié)果,最終得到304個穩(wěn)定特征波長點。篩選出的波長及由這些波長繪制的4臺近紅外儀器的平均光譜如圖2B所示。由該圖可見,根據(jù)這304個特征波長點的光譜響應(yīng)所繪的圖雖然不是很光滑,但保留了原樣品光譜(圖2A)的基本圖像特征,說明SIFT方法篩選的波長確實提取出了煙葉樣品的主要光譜特征。
由圖3可知,SIFT方法篩選出的大部分特征波長點在儀器間光譜差異很小的區(qū)域,那么以這304個波長的光譜響應(yīng)為自變量,建立的煙葉總植物堿的SIFT-PLS光譜模型對主、從機樣品預(yù)測結(jié)果應(yīng)有良好的一致性。
圖3 主機M與從機S2的SNV+一階導(dǎo)數(shù)光譜差譜的 絕對值及SIFT篩選波長Fig.3 The absolute difference spectrum between instrument M and S2 and the wavelengths selected by SIFT algorithm
表4給出了SIFT-PLS模型、全波長PLS(Whole wavelength PLS,簡稱WW-PLS)模型以及WW-PLS模型與SIFT-PLS模型經(jīng)PDS校正從機光譜后模型(PDS-WW-PLS、PDS-SIFT-PLS)的結(jié)果。從表4可以看出,WW-PLS在主機M以及從機S2上的預(yù)測誤差滿足企業(yè)內(nèi)控要求(MRE<6%),其對另外2臺從機樣品的MRE均大于6%。經(jīng)PDS校正從機光譜后,PDS-WW-PLS模型對S1、S3從機樣品的MRE分別從9.13%、8.18%降為5.30%與6.12%,而對S2的MRE反而從3.85%增大到8.13%;SIFT-PLS模型預(yù)測S2從機樣品總植物堿的MRE經(jīng)PDS校正后從4.24%增大到4.88%。說明NIRS模型直接轉(zhuǎn)移到從機的誤差滿足要求時,采用PDS校正反而會增大模型傳遞的誤差;在模型傳遞誤差較高時,PDS能起到一定改善作用,但也不能保證WW-PLS對3臺從機的傳遞誤差均滿足企業(yè)內(nèi)控要求。而SIFT-PLS模型對3臺從機樣品的MRE均滿足小于6%的企業(yè)內(nèi)控要求,當(dāng)采用PDS校正從機光譜后,該模型對3臺從機樣品總植物堿的RMSEP均較SIFT-PLS模型直接轉(zhuǎn)移有所降低,且PDS-SIFT-PLS模型對從機樣品的MRE均低于PDS-WW-PLS模型。說明SIFT算法篩選波長所建模型較全波長模型更為穩(wěn)健、精簡,且PDS校正從機光譜后仍能保證SIFT-PLS模型對從機樣品的MRE≤5%。
表4 不同煙葉總植物堿模型預(yù)測結(jié)果的比較Table 4 Results of total alkaloids contents predicted by different calibration models
利用SIFT算法篩選近紅外光譜的關(guān)鍵點,波長篩選過程中無需從機樣品信息,篩選出的波長點保留了光譜的主要特征,其光譜響應(yīng)不隨光譜尺度變化而改變,是樣品光譜的穩(wěn)定特征點,以這些波長下的光譜信息為自變量建立的煙葉總植物堿近紅外光譜模型具有穩(wěn)健、移植性好的特點。該模型直接轉(zhuǎn)移到3臺從機,無需任何模型傳遞方法進行從機光譜或模型校正,對從機樣品總植物堿的MRE為4.24%~5.81%,經(jīng)PDS算法校正從機光譜后,PDS-SIFT-PLS模型預(yù)測從機樣品總植物堿的MRE為4.48~5.00%,均能滿足企業(yè)內(nèi)控要求;而PDS-WW-PLS模型預(yù)測兩臺從機樣品總植物堿的MRE大于6%,表明全波長模型的穩(wěn)健性和移植性均不及SIFT-PLS模型。
基于SIFT方法篩選近紅外光譜的穩(wěn)定特征波長建立煙葉總植物堿的近紅外光譜模型,可實現(xiàn)模型的無標(biāo)樣傳遞,簡便易行。該方法用于其他領(lǐng)域及樣品的模型傳遞有待進一步驗證。