基于伏安電子舌的枸杞產(chǎn)地快速辨識

殷廷家 - 楊正偉 - 國婷婷 - 王志強(qiáng) -孫 霞 李彩虹 - 袁文浩 -

(1. 山東理工大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，山東 淄博 255049；2. 山東理工大學(xué)農(nóng)業(yè)工程與食品科學(xué)學(xué)院，山東 淄博 255049)

枸杞(LyciumbarbarumL)屬茄科枸杞屬植物，是一種名貴的保健品，具有抗衰老、抗氧化、滋腎潤肺、益精明目的功效[1]。中國枸杞主要分布在新疆、陜西、甘肅、寧夏等省份，其中以寧夏枸杞最為著名。受利益驅(qū)動，市場上常會出現(xiàn)枸杞產(chǎn)地冒充現(xiàn)象，嚴(yán)重?fù)p害了消費(fèi)者的經(jīng)濟(jì)利益[2]。傳統(tǒng)枸杞產(chǎn)地辨識方法主要根據(jù)果實(shí)的質(zhì)量、大小、色澤等外觀特征進(jìn)行辨別[3]，但該方法對專業(yè)技能要求較高。張波等[4]采用理化分析方法對枸杞內(nèi)的多糖、甜菜堿含量進(jìn)行分析，并根據(jù)其差異性對不同產(chǎn)地枸杞進(jìn)行區(qū)分，但此類方法檢測流程復(fù)雜、分析周期較長。湯麗華等[5]采用近紅外光譜技術(shù)對枸杞進(jìn)行檢測，依據(jù)其光譜特性對不同產(chǎn)地枸杞進(jìn)行鑒別，但此類檢測儀器價(jià)格昂貴、不易攜帶，無法滿足現(xiàn)場檢測需求。

電子舌是一種利用多傳感器陣列結(jié)合模式識別方法對液態(tài)樣本進(jìn)行分析的新型儀器，具有操作簡單、檢測迅速、方便攜帶、再現(xiàn)性好、客觀性強(qiáng)等優(yōu)點(diǎn)。近年來已成功應(yīng)用于食品溯源[6]、摻假檢測[7]、貨期預(yù)測[8]等多個領(lǐng)域。伏安電子舌的模式識別技術(shù)是影響其性能的關(guān)鍵因素，其過程通常分為特征提取和分類決策兩個步驟。特征提取的主要作用是對傳感器陣列產(chǎn)生的高維稀疏信號進(jìn)行壓縮和信息抽取。韓劍眾等[9]通過提取電子舌信號頂點(diǎn)和拐點(diǎn)值作為特征點(diǎn)，對魚肉的品質(zhì)和新鮮度進(jìn)行了評價(jià)。劉晶晶等[10]通過提取峰電流及極小值電流作為特征點(diǎn)，對雞蛋不同等級以及儲存時(shí)間進(jìn)行研究。任奇鋒等[11]利用響應(yīng)信號與時(shí)間軸包圍面積絕對值的和作為特征值，對不同種類的牛奶、啤酒以及不同等級的茶葉進(jìn)行了檢測。此類特征提取方法主要依據(jù)個人經(jīng)驗(yàn)，容易造成有效信息的丟失且很難進(jìn)行統(tǒng)一的推廣。另一類方法是基于空間變換思想，如陳茂晴[12]采用主成分分析(Principal component analysis，PCA)對電子舌信號特征提取，對金耳發(fā)酵過程進(jìn)行監(jiān)測，但PCA方法在將原始信號向低維空間進(jìn)行投影的過程中容易造成分類特征的缺失[13]。Lu等[14]采用快速傅里葉變換(Fast Fourier transform，F(xiàn)FT)對電子舌信號進(jìn)行特征提取，對不同地區(qū)的大米進(jìn)行了區(qū)分，但傅里葉變換無法對信號的時(shí)頻域局部信息進(jìn)行有效表達(dá)，影響了特征提取效果。此外離散小波變換 (Discrete wavelet transform，DWT)[15-16]也被應(yīng)用于伏安電子舌的特征提取領(lǐng)域，DWT雖然能夠在不同尺度上提取信號的時(shí)頻域細(xì)節(jié)信息，但其使用時(shí)需要結(jié)合實(shí)測信號預(yù)先選擇小波基和分解尺度，增加了使用的復(fù)雜性和難度。

希爾伯特—黃變換(Hilbert-Huang transform，HHT)[17]是一種以瞬時(shí)頻率為基礎(chǔ)的非線性、非平穩(wěn)信號時(shí)頻域聯(lián)合分析方法。其過程包含經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical mode decomposition，EMD)和Hilbert譜分析。與FFT、DWT相比，HHT徹底擺脫了Heisenberg測不準(zhǔn)原理制約，在時(shí)間和頻率上均能達(dá)到很高的精度。HHT同時(shí)具有較強(qiáng)的自適應(yīng)性，可在無需預(yù)設(shè)任何的基函數(shù)的條件下分析并提取信號的局部化特征，目前已經(jīng)在腦電波分析[18]、地震波形檢測[19]、故障診斷[20]等領(lǐng)域得到廣泛應(yīng)用，但將HHT應(yīng)用于伏安電子舌信號特征提取目前尚未見報(bào)道。

本研究擬采用伏安型電子舌(Voltammetry electronic tongue，VE-tongue)結(jié)合HHT-LDA模式識別方法對不同地區(qū)枸杞進(jìn)行快速鑒別。針對電子舌響應(yīng)信號的特點(diǎn)，采用EMD的改進(jìn)算法集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble empirical mode decomposition，EEMD)對電子舌原始信號進(jìn)行多尺度分解，得到一組IMF分量，分別求取IMF的奇異譜熵和Hilbert邊際譜作為電子舌信號的特征向量，并利用LDA建立組合枸杞產(chǎn)地非線性組合預(yù)測模型，以期為枸杞的快速、低成本產(chǎn)地鑒別和溯源提供新的方法和思路。

1 材料與方法

1.1 材料與樣本制備

枸杞：分別取自寧夏回族自治區(qū)中寧縣、新疆維吾爾自治區(qū)精河縣、青海省諾木洪和甘肅省景泰縣，每個產(chǎn)地20組樣本，采摘日期均為2017年6～9月份。將5 g干枸杞樣本放入200 mL純凈沸水中浸泡10 min，隨后用濾紙將濾液濾出，冷卻至室溫等待檢測。

1.2 伏安型電子舌

本試驗(yàn)采用實(shí)驗(yàn)室自行研制的伏安型電子舌系統(tǒng)[21]，其系統(tǒng)結(jié)構(gòu)如圖1所示。該系統(tǒng)由傳感器陣列、信號調(diào)理電路、數(shù)據(jù)采集卡(NI-6002)和LabVIEW上位機(jī)軟件組成。傳感器陣列由標(biāo)準(zhǔn)的三電極系統(tǒng)組成，包括8個工作電極(鉑、金、鈦、鈀、銀、鎢、鎳和玻碳)，1個Ag/AgCl 電極作參比電極和1個鉑輔助電極。檢測時(shí)，LabVIEW軟件控制數(shù)據(jù)采集卡在輔助電極與工作電極之間施加大幅方波脈沖伏安信號(Large amplitude pulse voltammetry，LAPV)，每個工作電極的激勵信號范圍為-1～1 V，電位階躍為0.2 V。在LAPV信號激勵下，工作電極表面發(fā)生電化學(xué)反應(yīng)產(chǎn)生微弱電流信號，該信號經(jīng)信號調(diào)理電路處理并經(jīng)數(shù)據(jù)采集卡進(jìn)行數(shù)據(jù)采集后，送至LabVIEW軟件進(jìn)行模式識別處理。

1.3 數(shù)據(jù)處理方法

1.3.1 集合經(jīng)驗(yàn)?zāi)B(tài)分解 EMD方法可以將非平穩(wěn)信號分解成一組不同尺度、互不耦合的IMF分量，但實(shí)際應(yīng)用中發(fā)現(xiàn)EMD會產(chǎn)生模態(tài)混疊現(xiàn)象，導(dǎo)致IMF分量中存在多種的信號震蕩模態(tài)，從而造成IMF部分物理意義的缺失。EEMD通過在原始信號添加小幅白噪聲來均衡信號，以提高信號在頻率范圍內(nèi)的均勻性和時(shí)間尺度上的連續(xù)性，從而有效解決了EMD的模態(tài)混疊現(xiàn)象，極大提高了經(jīng)驗(yàn)?zāi)B(tài)分解的準(zhǔn)確度[22]。對于一給定的非平穩(wěn)信號X(t)，EEMD分解過程如下：

1. 參比電極 2. 輔助電極 3. 工作電極 4. 信號調(diào)理電路 5. 數(shù)據(jù)采集卡 6. LabVIEW上位機(jī)

(1) 在原始信號中添加滿足(0,(αε)2)正態(tài)分布的白噪聲nj(t)，α和ε分別為添加的噪聲幅值和信號的標(biāo)準(zhǔn)差，獲得加噪后的信號：

Xj(t)=X(t)+nj(t)，

(1)

式中：

Xj(t)——第j次加噪后的信號；

X(t)——原始非平穩(wěn)信號；

nj(t)——第j次添加的高斯白噪聲。

(2) 對第j次添加白噪聲后的信號Xj(t)進(jìn)行EMD分解，得到一組IMF分量cij(t)(j=1,2,3,……,J)，其中cij(t)為第j次分解后獲得的第i階IMF分量。

(3) 對cij(t)取平均值得到ci(t)，以消減因多次添加白噪聲信號對實(shí)際IMF的干擾：

(2)

(4) 整個分解過程結(jié)束后，將IMF信號疊加則可以重構(gòu)原信號：

(3)

式中：

X(t)——原始非平穩(wěn)信號；

ci(t)——第i階IMF分量；

rn(t)——多次分解后的殘余量，通常為常數(shù)或者單調(diào)函數(shù)，僅含有信號變化趨勢信息，信息量較少，因此后續(xù)分析中將此分量忽略。

1.3.2 IMF奇異譜熵 奇異譜分析(Singular spectrum analysis，SSA)通過對原始離散的時(shí)間序列信號構(gòu)建軌跡矩陣，然后通過分解、重構(gòu)得到代表原始序列的特征信息，目前已應(yīng)用于多種時(shí)間序列的分析[23-24]中。IMF奇異譜熵[25]則是融合了奇異譜分析以及信息熵的概念，其實(shí)現(xiàn)步驟如下：

(1) 將各階IMF分量ci(t)構(gòu)建成一個新的軌跡矩陣：

Ai=

(4)

式中：

Ai——第i階IMF分量的軌跡矩陣；

M——信號的數(shù)據(jù)點(diǎn)數(shù)；

L——分段長度。

(2) 對軌跡矩陣Ai進(jìn)行奇異值分解

Ai=UXVT，

(5)

式中：

U∈R(M-L+1)×(M-L+1)；

V∈RL×L——正交矩陣；

VT——V的轉(zhuǎn)置矩陣；

X∈R(M-L+1)×(L)——對角矩陣。

對角元素X1,X2,……,Xm為矩陣Ai的奇異值且m=min(M-L+1,L)。

(3) 計(jì)算奇異譜熵，定義信號的奇異譜熵為：

(6)

式中：

E——奇異譜熵；

1.3.3 Hilbert邊際譜 Hilbert邊際譜為每個頻率值所對應(yīng)的能量總和，能夠準(zhǔn)確地將信號中包含的實(shí)際頻率成分反映出來并從統(tǒng)計(jì)意義上表征各個頻率點(diǎn)的累計(jì)幅值分布。其求解過程為：

(1) 對EEMD分解獲得的有效IMF分量進(jìn)行Hilbert變換，將實(shí)信號轉(zhuǎn)變成解析信號，然后將所有IMF分量的解析信號進(jìn)行疊加，并在重構(gòu)信號時(shí)去掉殘余量，這樣得到的時(shí)頻分布即為Hilbert譜：

(7)

式中：

P——柯西主值；

ci(t)——第i階IMF分量；

H[ci(t)]——對第i階IMF分量作Hilbert變換。

(2) 構(gòu)造其解析信號：

Zi(t)=ci(t)+jH[ci(t)]=ai(t)ejθi(t)，

(8)

(9)

(10)

式中：

Zi(t)——解析信號；

ai(t)——ci(t)的解析信號幅值，V；

θi(t)——解析信號相位，rad。

(3) 對θi(t)求導(dǎo)則可得到瞬時(shí)頻率ω(t)：

(11)

(4) 去除殘余分量rn(t)后，可以得到時(shí)頻分布的Hilbert譜：

(12)

式中：

Re[]——對信號取實(shí)部。

(5)H(ω,t)在時(shí)間上的積分可以得到Hilbert邊際譜：

(13)

1.3.4 LDA分類識別算法 LDA是一種有監(jiān)督的機(jī)器學(xué)習(xí)算法，可以通過已知的樣本類別來計(jì)算并尋找能夠?qū)崿F(xiàn)最好分類效果的線性判別函數(shù)。其基本思想是將高維的數(shù)據(jù)樣本映射到最佳鑒別矢量空間，以達(dá)到提取分類信息和壓縮特征維數(shù)的目的，映射后數(shù)據(jù)具有類內(nèi)間距小且類間間距大的特點(diǎn)，從而使得不同類別之間差異更加明顯。采用累計(jì)貢獻(xiàn)率衡量其效果，一般當(dāng)累計(jì)貢獻(xiàn)率≥85%時(shí)認(rèn)為該模型能夠較好的表示原始數(shù)據(jù)信息，使不同類別樣本能夠更好地區(qū)分[26]。

2 結(jié)果與分析

2.1 電子舌響應(yīng)信號

圖2為枸杞溶液響應(yīng)信號，在大幅脈沖信號作用下，電子舌可以獲取到豐富的樣本信息。從圖2中可以看出，8個工作電極的響應(yīng)信號存在差異，依次對8個工作電極進(jìn)行數(shù)據(jù)采集，每次采集共可獲取8 000個數(shù)據(jù)點(diǎn)，則4種枸杞樣本最終可獲得80×8 000的數(shù)據(jù)矩陣。

圖2 枸杞溶液信號

2.2 基于HHT的信號特征提取結(jié)果

2.2.1 EEMD結(jié)果 采用EEMD對電子舌采集信號進(jìn)行分解，EEMD添加白噪聲幅值設(shè)置α=0.2，執(zhí)行EMD次數(shù)J=100，標(biāo)準(zhǔn)差ε=0.25。經(jīng)EEMD分解后的信號如圖3所示，信號被自動分解為10個IMF分量及一個殘余分量，其中IMF信號強(qiáng)度呈逐漸減弱趨勢，為了減少信號中冗余信息，提高不同類別間的可區(qū)分度，因此采用奇異譜熵對IMF分量進(jìn)行有效性篩選。

2.2.2 IMF奇異譜熵 奇異譜熵是一種時(shí)域信號熵，熵值越大，則不確定性越高，包含的特征信息越多。對4種產(chǎn)地枸杞分別計(jì)算其10階IMF奇異譜熵，設(shè)置分段長度L=3 000，熵值分布情況如圖4所示。從圖4中可以看出，不同枸杞奇異譜熵整體分布趨勢大致相似，但熵值的大小存在差異，其值隨著IMF階數(shù)的增加呈逐漸遞減趨勢，說明隨著分解次數(shù)的增多，其所含信息量逐漸減少。由于第9、10階IMF分量熵值變化微弱，因此去除此兩階IMF分量僅保留1～8階IMF分量的奇異譜熵值作為特征向量，則一個信號樣本可得到1×8的特征向量。

圖3 原始溶液信號EEMD分解結(jié)果

圖4 不同枸杞的IMF奇異譜熵分布

2.2.3 Hilbert邊際譜 對1～8階IMF分量疊加后進(jìn)行Hilbert變換可以得到其Hilbert譜。Hilbert譜反映了電子舌信號時(shí)頻分布特性，對Hilbert譜在時(shí)間上進(jìn)行積分，可以獲得如圖5所示的Hilbert邊際譜。從圖5中可以看出，4種枸杞的Hilbert邊際譜變化趨勢大致相同，但幅值差異較大，新疆樣本最大值達(dá)到20，而甘肅樣本最大值只有7.6，可能是由于電子舌對不同產(chǎn)地枸杞的敏感性不同所致，因此產(chǎn)生的頻率累計(jì)幅值也存在較大差異。4種枸杞信號在10 Hz左右均有較為明顯的突起，其中甘肅樣本在0～10 Hz有較大的起伏，其他3種則較為平穩(wěn)。由于邊際譜的能量主要集中在0～50 Hz的區(qū)間范圍內(nèi)，且到達(dá)25 Hz左右幅值變化趨于平緩，為了兼顧信息的區(qū)分與壓縮效果，在0～25 Hz范圍內(nèi)以每次1 Hz的步長對Hilbert邊際譜提取特征點(diǎn)，可獲得25維的特征點(diǎn)。

圖5 4種枸杞信號的Hilbert邊際譜

2.3 LDA分類識別

將提取得到的奇異譜熵與Hilbert邊際譜特征點(diǎn)共33個特征向量輸入LDA模型進(jìn)行分類識別。為了證實(shí)基于HHT特征提取方法的有效性，分別采用FPE-LDA、PCA-LDA、DWT-LDA 3種模式識別方法與該方法進(jìn)行對比。其中FPE方法主要對電子舌響應(yīng)信號的極值點(diǎn)和拐點(diǎn)進(jìn)行提取，分別在每個脈沖信號提取2個拐點(diǎn)，1個極大值和1個極小值，最終得到320個特征點(diǎn)。采用PCA方法對電子舌信號進(jìn)行數(shù)據(jù)降維，當(dāng)累計(jì)貢獻(xiàn)率達(dá)到95%時(shí)，將數(shù)據(jù)作為模型輸入，則電子舌信號最終壓縮至10個數(shù)據(jù)點(diǎn)。為了找出DWT最優(yōu)的小波基和分解層數(shù)，對4個小波基(Symlet、Daubecges、Haar、Coiflet)進(jìn)行了不同階次和不同分解層數(shù)(5～9層)的測試，將分解后的數(shù)據(jù)進(jìn)行重構(gòu)，并求取與原始信號的相似系數(shù)。試驗(yàn)結(jié)果表明，Sym6母小波在分解8層相似系數(shù)達(dá)到最大值，此時(shí)壓縮效果最好，8 000個數(shù)據(jù)點(diǎn)最終壓縮至43個，波形相似系數(shù)為0.975 53。LDA分類識別結(jié)果如圖6所示。

圖6(a)為基于FPE-LDA的分類模型，其累計(jì)貢獻(xiàn)率(LD1、LD2)為90.1%。4種樣本整體距離較為接近，且各樣本類內(nèi)離散程度較大，新疆和青海樣本基本散落在相同區(qū)域，無法很好地區(qū)分，因此僅提取信號的極值點(diǎn)及拐點(diǎn)特征，雖然能夠較為直觀的體現(xiàn)不同類別間的差異，但同時(shí)也容易造成部分細(xì)節(jié)變化趨勢的缺失。圖6(b)為基于PCA-LDA的分類模型，累計(jì)貢獻(xiàn)率為89.4%，其中寧夏和新疆樣本區(qū)分度較好，甘肅和青海2種樣本有部分重合，可能是PCA算法屬于無監(jiān)督降維算法，在對信號進(jìn)行投影時(shí)僅選取貢獻(xiàn)率較大的主成分，會丟失一些重要的分類特征，但整體分類效果較FPE-LDA方法略有提高。圖6(c)為基于DWT-LDA的分類模型，其累計(jì)貢獻(xiàn)率僅為80.2%，略低于85%，新疆與甘肅樣本仍存在部分重疊現(xiàn)象，可能是DWT算法對信號進(jìn)行分解時(shí)會丟棄部分高頻分量，造成部分有用信息缺失導(dǎo)致。圖6(d)為基于HHT-LDA的分類模型，其累計(jì)貢獻(xiàn)率達(dá)到了91.6%，HHT算法主要對信號的時(shí)頻特征進(jìn)行分析，能夠很好地對非線性非平穩(wěn)信號進(jìn)行特征提取。從圖6(d)中可以看出，4種枸杞樣本被很好的區(qū)分開，且樣本類內(nèi)間距較小。試驗(yàn)結(jié)果表明，基于HHT-LDA分類模型比FPE-LDA、PCA-LDA、DWT-LDA 3種模型分類效果更好，表明在枸杞電子舌信號特征提取過程中，HHT算法能夠提取到更優(yōu)的鑒別特征。4種分類模型都能夠?qū)幭臉颖竞芎玫貐^(qū)分，說明了寧夏枸杞與其他地區(qū)枸杞物質(zhì)含量有著較大的差異，該結(jié)果驗(yàn)證了伏安電子舌結(jié)合HHT-LDA模型對不同產(chǎn)地枸杞分類的可行性和有效性，為枸杞產(chǎn)地快速分類提供了新的方法與思路。

圖6 4種模式識別模型

2.4 未知產(chǎn)地預(yù)測

為驗(yàn)證該模型對未知產(chǎn)地枸杞的分類識別能力，利用以上4種模型對未知產(chǎn)地的枸杞進(jìn)行產(chǎn)地預(yù)測。試驗(yàn)中每種枸杞分別制備20個樣本，從中隨機(jī)選取50個樣本(其中9個寧夏樣本，13個新疆樣本，16個甘肅樣本，12個青海樣本)進(jìn)行產(chǎn)地預(yù)測，其預(yù)測結(jié)果混淆矩陣如表1 所示。使用總體精度(Overall accuracy，OA)和Kappa系數(shù)作為評價(jià)指標(biāo)，經(jīng)過統(tǒng)計(jì)計(jì)算可知，基于FPE-LDA、PCA-LDA、DWT-LDA、HHT-LDA的分類模型總體分類精度分別為84%，92%，94%，98%，Kappa系數(shù)分別為0.784，0.891，0.919，0.973。基于HHT-LDA模型中只有一個甘肅樣本被錯分為了新疆樣本，總體識別效果較好，試驗(yàn)結(jié)果驗(yàn)證了該模型的有效性。

表1 預(yù)測結(jié)果混淆矩陣?

? 1、2、3、4分別為寧夏、新疆、甘肅、青海樣本。

3 結(jié)論

采用伏安型電子舌對不同產(chǎn)地枸杞進(jìn)行區(qū)分辨識，針對伏安型電子舌響應(yīng)信號數(shù)據(jù)量大，冗余信息較多的特點(diǎn)，本試驗(yàn)提出一種基于HHT-LDA的伏安電子舌模式識別方法，首先使用EEMD對信號進(jìn)行多尺度分解，通過選取有效IMF分量計(jì)算奇異譜熵，再對信號求取HHT邊際譜，組合成33維特征向量，最后將其輸入LDA模型進(jìn)行分類識別。將本特征提取方法與常用的FEP、PCA和DWT方法進(jìn)行對比，試驗(yàn)結(jié)果表明，基于HHT-LDA的分類模型可以很好地區(qū)分4種地區(qū)的枸杞樣本，表明HHT算法的特征提取性能要優(yōu)于以上3種方法。隨后建立枸杞產(chǎn)地預(yù)測模型，通過對未知產(chǎn)地樣本分類，從而驗(yàn)證該模型有效性，基于FPE-LDA、PCA-LDA、DWT-LDA和HHT-LDA 4種模型的總體分類精度分別為84%，92%，94%，98%，其中HHT-LDA模型僅有一個樣本被錯分，分類精度較高。試驗(yàn)結(jié)果驗(yàn)證了該模型的有效性，為快速鑒別枸杞產(chǎn)地提供了新的方法與思路。本方法僅對4種不同產(chǎn)地枸杞進(jìn)行了試驗(yàn)，在后續(xù)工作中將對該方法進(jìn)行更深的理論研究，并將其應(yīng)用于其他食品的檢測中。