基于可見_近紅外高光譜圖像的藥品快速鑒別研究

基于可見_近紅外高光譜圖像的藥品快速鑒別研究

江蘇雙利合譜科技有限公司

　　0 前言

　　化橘紅又名化皮、化州橘紅，為蕓香科植物化州柚的未成熟果實的外層果皮。前者習(xí)稱“毛橘紅”，后者習(xí)稱“光七爪”、“光五爪”?；偌t不僅具有治咳化啖、健胃行氣、醒酒功能，而還是人體美容的最佳原料，有廣闊的市場前景。研究表明，揮發(fā)油，黃酮類化合物，多糖以及香豆素類化合物等是化橘紅的主要有效成分。吳宋夏證實了化橘紅中的黃酮苷具有鎮(zhèn)痛、抗炎、解熱的作用；香豆素類化合物具有抗氧化、抗菌抗病毒之功效。不同的品種有效成分的含量不一樣，功效不一樣，并且在價格上也相差較大，以正品皮的效果最佳。因此市場上存在許多用化橘紅的正品果、偽品果、偽品皮冒充正品皮，損害了消費者利益，也沖擊了種植優(yōu)良品種的農(nóng)民們的利益。

　　目前對正品皮常用的鑒別方法主要有性狀鑒定、顯微鑒定、高效液相色譜法。這些方法雖然各有優(yōu)勢，但是存在不同程度上的主觀性強(qiáng)、需要預(yù)處理、實驗過程復(fù)雜等缺點，不能滿足市場快速、可靠檢測的需要。本研究分別利用400-1000nm、1000-2500nm兩款成像高光譜相機(jī)獲取正品皮、偽品皮、正品果、偽品果四種化橘紅成分的高光譜信息，利用光譜指數(shù)(DVI、NDVI)、偏最小二乘判別分析(PLS-DA)和極限學(xué)習(xí)法(ELM)三種方法分別構(gòu)建四種不同成分的鑒別模型，并用獨立樣本數(shù)據(jù)對不同的模型進(jìn)行驗證。

　　1 材料與方法

　　1.1 材料

　　試驗用的化橘紅四種不同成分正品皮、偽品皮、正品果、偽品果由廣東省食品藥品職業(yè)技術(shù)學(xué)校提供，其中正品皮樣本32個、正品果10個、偽品果11個，偽品皮7個。樣本經(jīng)粉碎均勻后，各取5g放置于培養(yǎng)皿上，標(biāo)號，用于高光譜相機(jī)的光譜采集。

　　1.2 高光譜圖像采集

　　高光譜圖像數(shù)據(jù)采集采用江蘇雙利合譜科技有限公司的 GaiaSorter高光譜分選儀系統(tǒng)（V10E、N25E-SWIR）。該系統(tǒng)主要由高光譜成像儀、面陣列相機(jī)、鹵素?zé)艄庠?、暗箱、計算機(jī)組成，如圖1。實驗儀器參數(shù)設(shè)置如表1。

圖1 GaiaSorter 高光譜分選儀

表1   GaiaSorter 高光譜分選儀系統(tǒng)參數(shù)

　　在進(jìn)行高光譜圖像采集時，需要設(shè)置相機(jī)曝光時間，平臺移動速度以及物鏡之間的距離。這 3 個參數(shù)相互影響，圖像調(diào)節(jié)的目的是使采集的圖像大小合適，清晰，不變形失真。經(jīng)過反復(fù)嘗試，物鏡高度設(shè)置為 31 cm，曝光時間設(shè)置為10ms，平臺移動速度分別設(shè)置為 6.0 mm/s（400-1000 nm）、16mm/s（1000-2500 nm）。圖像采集軟件采用江蘇雙利合譜科技有限公司提供的高光譜成像系統(tǒng)采集軟件完成。圖像處理采用 ENVI5.3 軟件進(jìn)行處理。在進(jìn)行圖像處理之前，先要對采集的光譜圖像進(jìn)行圖像校正，圖像校正公式如下：

  (1)

　　式中，R_ref 是校正過的圖像，DN_raw 是原始圖像，DN_white為白板校正圖像，DN_dark 是黑板校正圖像。高光譜圖像的光譜與圖像之間有著對應(yīng)的關(guān)系，據(jù)此，在正品皮、偽品皮、正品果、偽品果四種樣本上選取 25×25 像素點的感興趣區(qū)域，以感興趣區(qū)域所有像素的光譜平均值作為該樣本的平均光譜。

　　1.3  光譜噪聲去除

　　試驗得到光譜含有由儀器和試驗條件等引起的噪聲，對這些噪聲的處理有助于減少噪聲對光譜分析的影響，突出光譜的有效信息。SG 平滑算法可以有效消減光譜數(shù)據(jù)中的隨機(jī)噪聲，消噪效果受平滑點數(shù)的影響，本文中選擇SG 二次多項式 7 點平滑對光譜數(shù)據(jù)進(jìn)行處理（何勇，2013）。

　　1.4  特征波長選擇

　　光譜信息之間存在大量的冗余和共線性信息特征，對光譜有效信息的提取產(chǎn)生了較大的干擾，且大量光譜數(shù)據(jù)造成模型復(fù)雜、計算量大的問題。本文采用連續(xù)投影算法（successive projections algorithm，SPA）進(jìn)行特征波長的選擇，以減少信息冗余和共線性的影響，簡化模型，減少計算量。

　　SPA 是一種特征變量前向選擇算法，在光譜特征波長中取得了廣泛的應(yīng)用。本文采用 SPA 算法對去噪處理后的光譜進(jìn)行特征波長選擇。進(jìn)行SPA 計算時，以建模集的光譜數(shù)據(jù)和類別賦值為輸入，設(shè)置選擇特征波長數(shù)的范圍為 5～30。

　　1.5 光譜指數(shù)

　　光譜指數(shù)的產(chǎn)生來源于植被指數(shù)，植被指數(shù)是指利用衛(wèi)星不同波段探測數(shù)據(jù)組合而成的，能反映植物生長狀況的指數(shù)。植物葉面在可見光紅光波段有很強(qiáng)的吸收特性，在近紅外波段有很強(qiáng)的反射特性，這是植被遙感監(jiān)測利用衛(wèi)星不同波段探測數(shù)據(jù)組合而成的，能反映植物生長狀況的指數(shù)。植物葉面在可見光紅光波段有很強(qiáng)的吸收特性，在近紅外波段有很強(qiáng)的反射特性，這是植被遙感監(jiān)測的物理基礎(chǔ)，通過這兩個波段測值的不同組合可得到不同的植被指數(shù)。光譜指數(shù)是通過任意兩波段組合或三波段組合成各種光譜指數(shù)，如NDVI、DVI等，探尋最佳的波段組合用于各個領(lǐng)域的模型構(gòu)建等。

歸一化植被指數(shù)(Normalized difference vegetable index)

(2)

差值植被指數(shù)(Difference vegetable index)

(3)

其中，λ1和λ2代表任意波長的反射率，波段范圍為400-1000 nm與1000-2500 nm。

　　1.6 判別分析方法

　　偏最小二乘法判別分析(PLS-DA ,Partial least squares discrimination analysis)是一種用于判別分析的多變量統(tǒng)計分析方法。判別分析是一種根據(jù)觀察或測量到的若干變量值，來判斷研究對象如何分類的常用統(tǒng)計分析方法。其原理是對不同處理樣本（如觀測樣本、對照樣本）的特性分別進(jìn)行訓(xùn)練，產(chǎn)生訓(xùn)練集，并檢驗訓(xùn)練集的可信度（Luna et al., 2013）。本文分別基于全光譜、特別波長光譜建立 PLS-DA 判別分析模型，通過建立光譜數(shù)據(jù)與類別特征之間的回歸模型，進(jìn)行判別分析。

　　1.7 極限學(xué)習(xí)機(jī)

　　極限學(xué)習(xí)機(jī)（extreme learning machine）ELM是一種簡單易用、有效的單隱層前饋神經(jīng)網(wǎng)絡(luò)SLFNs學(xué)習(xí)算法。2004年由南洋理工大學(xué)黃廣斌副教授提出。傳統(tǒng)的神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)算法（如BP算法）需要人為設(shè)置大量的網(wǎng)絡(luò)訓(xùn)練參數(shù)，并且很容易產(chǎn)生局部最優(yōu)解。極限學(xué)習(xí)機(jī)只需要設(shè)置網(wǎng)絡(luò)的隱層節(jié)點個數(shù)，在算法執(zhí)行過程中不需要調(diào)整網(wǎng)絡(luò)的輸入權(quán)值以及隱元的偏置，并且產(chǎn)生唯一的最優(yōu)解，因此具有學(xué)習(xí)速度快且泛化性能好的優(yōu)點。本文中隱含層神經(jīng)元個數(shù)從 1 到 520（288）以步長 1 進(jìn)行尋優(yōu)，以最小訓(xùn)練誤差下的神經(jīng)元個數(shù)為 ELM 模型隱含層神經(jīng)元個數(shù)。

　　1.8  評價指標(biāo)

　　回歸模型得到的樣本的預(yù)測值不是整數(shù)，需要設(shè)置閾值以判斷樣本的歸屬。本文中閾值設(shè)置為 0.5 ，預(yù)測值小數(shù)點大于或等于0.5則加1歸整，小于0.5則減1歸整?？傮w識別精度是指正確識別的個數(shù)除以總數(shù)，正品皮識別精度是指正品皮正確識別的個數(shù)除以正品皮的總數(shù)，正品皮識別錯誤率指數(shù)被錯誤分為正品皮的個數(shù)除以正品皮的總數(shù)。

　　2  結(jié)果與分析

　　2.1  化橘紅不同成分的原始光譜曲線

　　本試驗采用的V10E 相機(jī)獲取的是400-1000 nm波長范圍共520個波段的可見/近紅外光譜數(shù)據(jù)，N25E-SWIR相機(jī)獲取的是1000-2500 nm波長范圍共288個波段的近紅外光譜數(shù)據(jù)，全部樣本的原始光譜圖如圖 2所示，正品皮、偽品皮、正品果、偽品果的光譜比較圖如圖3所示。

　　從圖2和圖3可以看到，總體而言，無論是400-1000 nm或1000-2500 nm波長范圍內(nèi)，正品皮的光譜反射率值低于其他三種成分的光譜曲線，從曲線變化趨勢來看四種不同成分并沒有十分明顯的差異。本研究按照Kennard-Stone 算法將樣本分成建模集和預(yù)測集，其中建模集 38 個樣本，預(yù)測集32個樣本。正品皮、偽品皮、正品果、偽品果分別賦值為 1、2、3、4（表 2），不同化橘紅成分建模集和預(yù)測集樣本的劃分如表 2所示。

圖2 全部化橘紅樣本的原始反射光譜圖

圖3 化橘紅不同成分反射光譜曲線圖

表2  化橘紅不同成分類別賦值與建模集合檢驗集樣本劃分

　　2.2  主成分分析定性分析

　　對化橘紅的四種不同成分的光譜數(shù)據(jù)進(jìn)行主成分分析（principal component analysis，PCA） 可知，400-1000 nm范圍內(nèi)，第一主成分（principle component 1, PC1）的貢獻(xiàn)率為88.36%，PC2 的貢獻(xiàn)率為7.24%，PC1 和 PC2 累計貢獻(xiàn)率為95.6%，PC1 和 PC2 能夠解釋絕大部分的變量；在1000-2500 nm范圍內(nèi)，第一主成分（principle component 1, PC1）的貢獻(xiàn)率為93.27%，PC2 的貢獻(xiàn)率為3.72%，PC1 和 PC2 累計貢獻(xiàn)率為97.0%，PC1 和 PC2 也能夠解釋絕大部分的變量；圖4分別為400-1000 nm與1000-2500 nm范圍內(nèi)，第一主成分與第二主成分的散點分布圖。

圖4  400-1000nm（左）、1000-2500nm（右）第一主成分與第二主成分的得分散點分布圖

圖5 主成分組合識別正品皮（左：400-1000 nm，右：1000-2500 nm）

　　由圖 4 可知，在400-1000 nm與1000-2500 nm范圍內(nèi)，各不同成分之間較難直接區(qū)分出來，得分圖中均有重合的地方。從圖5亦可以看出無論是400-1000 nm或1000-2500 nm，主成分的成分組合均未能完全從正品皮、偽品皮、正品果、偽品果四種樣本中識別出正品皮，從圖5可知，識別出的正品皮亦有假的正品皮。因此需要對光譜數(shù)據(jù)進(jìn)行進(jìn)一步的分析和處理，以鑒別化橘紅不同成分。

　　2.3  光譜指數(shù)

　　本研究使用的成像高光譜儀V10E范圍是400-1000nm，N25E-SWIR范圍是1000-2500 nm，通過不同波段的重新組合，形成不同形式的光譜指數(shù)(歸一化光譜指數(shù)，NDVI；差值光譜指數(shù)，DSI)。

　　運(yùn)用Matlab軟件編程，將兩兩組合的所有波段構(gòu)建的光譜指數(shù)與各類別賦值，計算相應(yīng)的決定系數(shù)(R²)，繪制決定系數(shù)圖。圖中顏色從藍(lán)色到紅色變化，圖像顏色越紅，表示決定系數(shù)越大，圖像顏色越藍(lán)，表示決定系數(shù)越小。圖6為化橘紅各成分分類賦值與NDVI、DVI決定系數(shù)（R2）二維圖。從圖6可知，在400-1000 nm范圍內(nèi)， DVI模型決定系數(shù)最高的兩波段在綠光范圍內(nèi)，分別是538.91 nm和543.75nm；NDVI模型絕對系數(shù)最高的兩波段989.76 nm和670.39 nm。在1000-2500 nm范圍內(nèi)，NDVI和DVI兩個光譜指數(shù)模型決定系數(shù)最高的兩波段組合均為是1820.72 nm和1787.12 nm。綜合400-1000 nm與1000-2500 nm

圖 6化橘紅各成分分類賦值與NDVI、DVI決定系數(shù)（R²）二維圖

　　兩組數(shù)據(jù)構(gòu)建NDVI和DVI模型時，NDVI和DVI兩個光譜指數(shù)模型決定系數(shù)最高的兩波段組合與1000-2500 nm范圍內(nèi)相同。比較400-1000 nm和1000-2500 nm的決定系數(shù)圖可知，在1000-2500 nm范圍內(nèi)，NDVI與DVI構(gòu)建的模型決定系數(shù)較高。

　　根據(jù)圖6可知，在400-1000 nm范圍內(nèi)，DVI_{(538.91/543.75)}與NDVI_{(989.76/670.39)}與類別賦值構(gòu)建的模型決定系數(shù)最高，分別為0.643和0.640。在1000-2500 nm范圍內(nèi)，DVI_{(1820.72 /1787.12)}與NDVI_{(1820.72 /1787.12)}與類別賦值構(gòu)建的模型決定系數(shù)最高，分別為0.861和0.834。圖7分別為DVI_{(538.91/543.75)}、NDVI_{(989.76/670.39)} 、DVI_{(1820.72 /1787.12)}、NDVI_{(1820.72 /1787.12)}與類別賦值的線性擬合散點分布圖，從圖中可知，在1000-2500 nm范圍內(nèi)構(gòu)建的NDVI、DVI光譜指數(shù)與類別賦值擬合度最高，且變化曲線顯著。

 

 

圖7  DVI、NDVI光譜指數(shù)與類別賦值的散點分布圖

　　運(yùn)用獨立的數(shù)據(jù)，分別對圖7中DVI_{(538.91/543.75)}、NDVI_{(989.76/670.39)} 、DVI_{(1820.72 /1787.12)}、NDVI_{(1820.72 /1787.12)}與類別賦值構(gòu)建的模型進(jìn)行檢驗，檢驗結(jié)果如圖8和表3所示。根據(jù)查表可知，1-15與28-32為正品皮，15-25號為16-20正品果、21-24為偽品果，25-27為偽品皮。從圖8可知，正品皮、正品果、偽品果、偽品皮之間均存在不同程度的錯誤識別。表3基于DVI、NDVI光譜指數(shù)模型檢驗化橘紅樣本的精度評價表。表中分別統(tǒng)計了總體識別精度、正品皮識別精度和正品皮識別錯誤率。從表3中可知，總體識別精度、正品皮識別精度最高的均為DVI_{(1820.72 /1787.12)}構(gòu)建的模型，分別是66%和75%，正品皮識別錯誤率最低的則為DVI_{(1820.72 /1787.12)}、NDVI_{(1820.72 /1787.12)} 構(gòu)建的模型，均為5%。

  

圖 8 基于NDVI、DVI檢驗集預(yù)測結(jié)果對比圖

表2 基于光譜指數(shù)模型檢驗化橘紅樣本的精度評價

	總體識別精度%	正品皮識別精度%	正品皮識別錯誤率%
400-1000 (DVI)	53	65	15
400-1000 (NDVI)	50	70	20
1000-2500 (DVI)	66	75	5
1000-2500 (NDVI)	59	65	5

 

　　2.4  特征波長選擇

　　本文以建模集樣本的光譜數(shù)據(jù)和類別賦值為輸入，利用連續(xù)投影算法SPA選擇特征波長。選出的特征波長的個數(shù)如表 3所示。從表 3 可知，400-1000 nm范圍內(nèi)所選擇的特征波段為15個，1000-2500 nm范圍內(nèi)所選擇的特征波段較少，為5個。

表 3在400-1000 nm與1000-2500 nm范圍內(nèi)SPA 算法選擇的特征波長個數(shù)

范圍	波段位置/nm
400-1000 nm	395.23, 396.40, 399.90, 401.06, 422.13, 501.59, 670.39, 735.79, 743.39, 872.13, 940.31, 962.16, 974.58, 981.45, 992.73
1000-2500 nm	1461.59, 1714.23, 2038.97, 2329.25, 2574.24

 

　　2.5  偏最小二乘判別分析

　　分別將建模集合檢驗集的全波段光譜及特征波段光譜作為偏最小二乘法判別分析輸入變量，從而獲取預(yù)測樣本類別的賦值。圖9為檢驗集與預(yù)測值的類別賦值圖，并針對圖9進(jìn)行總體識別精度、正品皮識別精度、正品皮識別錯誤率的統(tǒng)計，如表4所示。從表4中可知，總體識別精度、正品皮識別精度最高的均為在1000-2500 nm范圍內(nèi)的PLS-DA構(gòu)建的模型，分別是78%和90%，正品皮識別錯誤率最低的則為1000-2500 nm范圍內(nèi)的PLS-DA、PLS-DA 構(gòu)建的模型，均為5%。從表4可以看出，無論是400-1000 nm或1000-2500 nm，基于全波段的PLS-DA模型總體識別率和正品皮識別率均高于基于特征波段的PLS-DA模型，而正品皮的錯誤識別率，無論是400-1000 nm或1000-2500 nm范圍內(nèi)，基于全波段的PLS-DA模型與基于特征波長的PLS-DA模型的錯誤識別率相同，400- 1000 nm范圍錯誤識別率均為10%，1000-2500范圍錯誤識別率則為5%。

圖 9 基于PLS-DA檢驗集預(yù)測結(jié)果對比圖

表4基于PLS-DA模型檢驗化橘紅樣本的精度評價

	總體識別精度%	正品皮識別精度%	正品皮識別錯誤率%
400-1000 (PLS-DA)	72	65	10
400-1000 (PLS-DA-SPA)	63	55	10
1000-2500 (PLS-DA)	78	90	5
1000-2500 (PLS-DA-SPA)	72	70	5

 

　　2.6  極限學(xué)習(xí)機(jī)

　　分別將建模集合檢驗集的全波段光譜及特征波長光譜作為極限學(xué)習(xí)機(jī)的輸入變量，從而獲取預(yù)測樣本類別的賦值。圖10為檢驗集與預(yù)測值的類別賦值圖，并針對圖10進(jìn)行總體識別精度、正品皮識別精度、  正品皮識別錯誤率的統(tǒng)計，如表5所示。從表5中可知，總體識別精度、正品皮識別精度最高的均為在1000-2500 nm范圍內(nèi)的ELM與ELM-SPA構(gòu)建的模型，分別是84%和95%，正品皮識別錯誤率最低的則為1000-2500 nm范圍內(nèi)的ELM與ELM-SPA 構(gòu)建的模型，均為5%。從表5可知，在400-1000 nm范圍內(nèi)，基于特征波段光譜的ELM模型總體識別率與正品皮識別率均高于基于全波段的ELM模型，對于正品皮的識別錯誤率，基于全波段與基于特征波段的錯誤識別率相同；在1000-2500 nm范圍內(nèi)，無論是基于全波段光譜或基于特征波段的ELM模型，其總體識別率、正品皮識別率、正品皮識別錯誤率均相同，分別為84%、95%和5%。

 

  

圖 10 基于ELM檢驗集預(yù)測結(jié)果對比圖

表4 基于ELM模型檢驗化橘紅樣本的精度評價

	總體識別精度%	正品皮識別精度%	正品皮識別錯誤率%
400-1000 (ELM)	75	75	10
400-1000 (ELM-SPA)	78	85	10
1000-2500 (ELM)	84	95	5
1000-2500 (ELM-SPA)	84	95	5

 

　　2.7  光譜指數(shù)模型，PLS-DA 模型，和 ELM 模型的比較

　　綜合對比光譜指數(shù)模型，PLS-DA 模型，和 ELM 模型的識別效果可知，無論是光譜指數(shù)模型，PLS-DA 模型或ELM 模型，基于1000-2500 nm范圍內(nèi)構(gòu)建的模型，其預(yù)測值的總體識別率、正品皮識別率均高于400-1000 nm范圍內(nèi)的模型，且正品皮的識別錯誤率也低于400-1000nm范圍內(nèi)的模型。在光譜指數(shù)模型、PLS-DA 模型和 ELM 模型的模型中，ELM 模型的識別準(zhǔn)確性最高，其次是PLS-DA模型，最后是光譜指數(shù)模型?；谔卣鞑ǘ喂庾V的PLS-DA模型其識別準(zhǔn)確性低于基于全波段光譜的PLS-DA的模型，但是基于特征波段光譜的ELM模型在400-1000 nm范圍內(nèi)，其識別準(zhǔn)確性高于基于全波段光譜的ELM模型，在1000-25000 nm范圍內(nèi)，其識別準(zhǔn)確性與基于全波段光譜的ELM模型相同。

　　3  結(jié)論與討論

　　基于V10E與N25E-SWIR兩款成像高光譜相機(jī)，分別獲取正品皮、正品果、偽品皮、偽品果四種化橘紅成分400-1000 nm與1000-2500 nm范圍的光譜反射率，采用 SG 平滑算法對提取出的光譜數(shù)據(jù)進(jìn)行去噪處理，同時采用 SPA 算法對去噪后的光譜提取特征波長，并分別基于全波段光譜、特征波段光譜建立 PLS-DA 判別模型和 ELM 模型，同時采用全波段循環(huán)，探尋最佳的NDVI、DVI兩個光譜指數(shù)構(gòu)建判別模型，用于鑒別正品皮、正品果、偽品皮、偽品果，取得了比較好的識別效果。基于特征波段光譜與全波段光譜建立的 ELM 模型取得了最佳效果，總體識別精度、正品皮識別精度、正品皮識別錯誤率分別為84%、95%和5%。在實際運(yùn)用中，考慮到識別時間與成分，基于SPA算法提取的特征波段構(gòu)建的ELM模型效果最佳。本論文研究結(jié)果為高光譜成像技術(shù)在藥品真?zhèn)蔚葯z測中的應(yīng)用提供了可行性。