周 俊，董 琪，莊柳靜，吳春生，劉清君，王 平

(浙江大學(xué)生物醫(yī)學(xué)工程與儀器科學(xué)學(xué)院，生物傳感器國(guó)家專(zhuān)業(yè)實(shí)驗(yàn)室，生物醫(yī)學(xué)工程教育部重點(diǎn)實(shí)驗(yàn)室，杭州310027)

生物傳感技術(shù)被廣泛應(yīng)用于感覺(jué)、運(yùn)動(dòng)系統(tǒng)生理機(jī)制研究等領(lǐng)域，具有實(shí)時(shí)、快速、靈敏、簡(jiǎn)便等優(yōu)點(diǎn)。動(dòng)物的嗅覺(jué)十分靈敏，利用動(dòng)物嗅覺(jué)構(gòu)建的傳感系統(tǒng)在辨識(shí)氣味方面具有獨(dú)特的優(yōu)勢(shì)。嗅球是哺乳動(dòng)物氣味信息編碼的重要器官，它接收來(lái)自嗅上皮嗅覺(jué)感受細(xì)胞的電興奮，而后者將氣味刺激這一化學(xué)信號(hào)轉(zhuǎn)換成電信號(hào)用于編碼[1]。研究者發(fā)現(xiàn)氣味刺激前后嗅球僧帽細(xì)胞層神經(jīng)電活動(dòng)十分強(qiáng)烈，除了神經(jīng)元個(gè)體的脈沖式放電行為之外，大量神經(jīng)元?jiǎng)幼麟娢化B加表現(xiàn)出較大的電位波動(dòng)，稱(chēng)之為場(chǎng)電位。目前的研究表明，氣味刺激后嗅球內(nèi)的局部場(chǎng)電位振蕩會(huì)明顯增加，其幅值和波形的變化受到氣味種類(lèi)的調(diào)控[2－3]。

多窗譜估計(jì)算法是基于FFT的一種譜估計(jì)算法，與傳統(tǒng)的周期圖法只用一個(gè)數(shù)據(jù)窗[4]不同，該方法采用了多個(gè)數(shù)據(jù)窗進(jìn)行平滑，因此在譜估計(jì)結(jié)果上具有良好的方差性能[5]。同時(shí)由于數(shù)據(jù)窗相互正交，可有效阻止頻譜泄露，保證該方法的譜估計(jì)分辨率。在多窗譜估計(jì)中加入合適的移動(dòng)窗可以得到頻率隨時(shí)間分布的功率譜時(shí)頻圖，克服了頻譜分析中無(wú)法描述信號(hào)局部特征[6]的不足，從而突出在特定事件前后不同頻率段的信號(hào)能量變化情況[7]。利用這一特征，功率譜時(shí)頻圖可以對(duì)不同氣味刺激后的嗅覺(jué)信號(hào)特征進(jìn)行表征。K最鄰近分類(lèi)算法K-NNA(K-Nearest Neighbor Algorithm)是理論較為成熟的方法，將該分類(lèi)算法用于分類(lèi)不同氣味刺激下信號(hào)的功率譜時(shí)頻圖可以區(qū)分引起不同嗅球振蕩頻率特征的氣味。

本文介紹了我們采用微電極陣列傳感器植入大鼠嗅球構(gòu)建新型嗅覺(jué)傳感系統(tǒng)的探索性研究，分析了該傳感系統(tǒng)在不同氣味刺激下嗅球場(chǎng)電位信號(hào)的特征，并對(duì)采集得到的信號(hào)通過(guò)多窗譜估計(jì)與K最鄰近分類(lèi)，實(shí)現(xiàn)了氣味的初步識(shí)別。

1 大鼠嗅覺(jué)傳感系統(tǒng)的設(shè)計(jì)及信號(hào)采集

大鼠嗅覺(jué)傳感系統(tǒng)以大鼠嗅覺(jué)感受細(xì)胞作為氣味敏感傳感單元，通過(guò)記錄和分析具有氣味刺激特征的嗅球電位信號(hào)實(shí)現(xiàn)氣味識(shí)別。傳感系統(tǒng)主要由氣味刺激裝置、多通道電極陣列傳感器、多通道信號(hào)采集裝置及計(jì)算機(jī)數(shù)據(jù)分析模塊4部分組成。氣味刺激裝置由濃度為10 mmol/L液體形式存儲(chǔ)的氣味瓶、微型真空泵和電磁閥組成，能輸出0.2 L/min的頂空氣體。實(shí)驗(yàn)中記錄信號(hào)的多通道電極陣列傳感器采用自制8通道，直徑為65 μm鎳鎘合金微絲(AM System，USA)的電極陣列，電極阻抗為300 kΩ～500 kΩ(1 kHz)。多通道信號(hào)采集裝置為OmniPlex數(shù)據(jù)采集系統(tǒng)(Plexon，USA)。該信號(hào)采集器具有20倍的前置放大器，16 bit的A/D轉(zhuǎn)換器，能以最高40 kHz頻率5 000倍信號(hào)增益實(shí)現(xiàn)128個(gè)通道實(shí)時(shí)同步數(shù)據(jù)采集，并帶有用于低頻場(chǎng)電位和高頻鋒電位的頻帶隔離可編程硬件濾波器。整個(gè)實(shí)驗(yàn)臺(tái)上架設(shè)了由120目銅網(wǎng)搭建的金屬屏蔽罩用于屏蔽環(huán)境周?chē)赡軒?lái)的噪聲，如工頻噪音(50 Hz)、射頻噪音等。信號(hào)采集裝置通過(guò)網(wǎng)線將采集到的數(shù)據(jù)發(fā)送到計(jì)算機(jī)，利用Matlab編譯環(huán)境實(shí)現(xiàn)數(shù)據(jù)的處理和分析。

實(shí)驗(yàn)對(duì)象SD大鼠購(gòu)于浙江省醫(yī)學(xué)科學(xué)院，質(zhì)量是200 g～250 g。對(duì)大鼠實(shí)施電極植入手術(shù)，流程是在腹腔注射10%的水合氯醛(4 mL/kg)，完全麻醉后移除單側(cè)嗅球組織上覆蓋的頭骨，利用油壓微推進(jìn)器(Narishige，Japan)將電極植入到嗅球僧帽細(xì)胞層，平均深度約為 200 μm ～ 300 μm[8]。實(shí)驗(yàn)期間大鼠體下放置加熱毯，以保證其體溫在37℃左右。使用濃度為10 mmol/L的異丁醇、苯甲醚、香芹酮和檸檬醛4種氣體作為刺激物，單次刺激給予體積為13 mL、時(shí)間為4 s的氣味，刺激間隔為60 s。整個(gè)實(shí)驗(yàn)都是在通風(fēng)的環(huán)境中進(jìn)行，并間隔一定時(shí)間用風(fēng)扇加速空氣流動(dòng)，以避免氣味刺激后仍殘余在實(shí)驗(yàn)臺(tái)附周?chē)?。信?hào)采集裝置同時(shí)記錄低頻場(chǎng)電位信號(hào)和刺激開(kāi)始結(jié)束兩個(gè)標(biāo)記。由于嗅球僧帽細(xì)胞放電幅值較大較明顯，能說(shuō)明電極植入到目標(biāo)細(xì)胞層，故我們選取氣味刺激時(shí)能記錄到細(xì)胞鋒電位信號(hào)的通道，分析該通道的低頻場(chǎng)電位信號(hào)。圖1所示為大鼠嗅覺(jué)傳感系統(tǒng)和采集到的嗅球低頻場(chǎng)電位原始信號(hào)。信號(hào)采樣頻率為1 000 Hz，信號(hào)幅值通常為200 mV～300 mV，信噪比大于5∶1，能保證長(zhǎng)時(shí)間穩(wěn)定的記錄。

圖1 嗅覺(jué)傳感系統(tǒng)及采集到的嗅球場(chǎng)電位信號(hào)

2 基于頻譜的氣味識(shí)別算法

2.1 多窗譜估計(jì)算法

多窗譜估計(jì)MSE(Multitaper Spectrum Estimation)算法是由Thomson在1982年提出的。該算法的基本思路是對(duì)同一段隨機(jī)信號(hào)使用多組互相正交的數(shù)據(jù)窗分別直接求功率譜，然后進(jìn)行平均得到功率譜估計(jì)[9]。

多窗譜定義如下:

其中K是使用的窗個(gè)數(shù)，Sk(ω)是第k個(gè)功率譜，計(jì)算方法如下[7]:

式(2)中N是數(shù)據(jù)序列的長(zhǎng)度，hk(n)是第k個(gè)數(shù)據(jù)窗序列。一般hk(n)序列是選取一組相互正交的離散橢球序列(DPSS序列)，且該序列滿(mǎn)足:

2.2 功率譜時(shí)頻圖

通過(guò)選擇移動(dòng)窗寬度與每次譜估計(jì)的時(shí)間步長(zhǎng)，對(duì)每個(gè)窗的數(shù)據(jù)進(jìn)行多窗法(MTM)譜估計(jì)，構(gòu)建以時(shí)間為橫坐標(biāo)，頻率為縱坐標(biāo)的功率譜時(shí)頻圖(Spectrogram)。

假設(shè)大鼠在不同氣味刺激下選取連續(xù)氣味刺激時(shí)長(zhǎng)為4 s的低頻場(chǎng)電位信號(hào)(角標(biāo)i表示刺激氣味種類(lèi)，分別對(duì)應(yīng)異丁醇、苯甲醚、香芹酮和檸檬醛4種氣味，上標(biāo)k表示第i種氣味的第k次刺激)。對(duì)信號(hào)選取時(shí)間窗為win，步長(zhǎng)為step，進(jìn)行多窗法譜估計(jì)，并選取氣味刺激前的3個(gè)時(shí)間窗數(shù)據(jù)平均值作為基線信號(hào)用于該段功率譜的規(guī)范化處理，得到隨時(shí)間變化的功率譜向量。

2.3 K最鄰近分類(lèi)算法

K最鄰近分類(lèi)算法(KNN)是一種基于向量空間模型和實(shí)例學(xué)習(xí)的對(duì)象有監(jiān)督的分類(lèi)方法，其基本思想是:給定一個(gè)經(jīng)過(guò)分類(lèi)的訓(xùn)練樣本集合，在對(duì)新對(duì)象進(jìn)行分類(lèi)時(shí)，首先從訓(xùn)練樣本中找出與新對(duì)象最相似的K個(gè)訓(xùn)練樣本，并根據(jù)這K個(gè)訓(xùn)練樣本所屬的類(lèi)別，判定待定對(duì)象所屬的類(lèi)別。

本文采用KNN算法分類(lèi)氣味刺激信息相關(guān)的功率譜時(shí)頻圖。4種氣味刺激共有120組時(shí)頻圖向量，完成二重交叉檢驗(yàn)。具體算法步驟如下:(1)從每種氣味相關(guān)的30組功率譜向量中隨機(jī)選取其中的一半數(shù)據(jù)，共60組向量的集合作為訓(xùn)練樣本 x={x1，x2，…，x30}，并建立對(duì)應(yīng)的氣味標(biāo)識(shí)矩陣G。(2)用剩余的另一半向量作為待分類(lèi)的新樣本y={y1，y2，…，y30};(3)計(jì)算新樣本和訓(xùn)練樣本之間的相似度，采用兩者之間的歐式距離:

在訓(xùn)練集中選出與待分類(lèi)對(duì)象最接近的K個(gè)樣本;(4)在K個(gè)鄰居中通過(guò)投票的方式選取投票最多的類(lèi)別，并將對(duì)象歸于此類(lèi)別。將訓(xùn)練集與樣本集對(duì)換，完成二重交叉檢驗(yàn)，取兩者識(shí)別準(zhǔn)確率的平均值作為最終識(shí)別結(jié)果。KNN算法的優(yōu)點(diǎn)是這類(lèi)算法對(duì)于訓(xùn)練集中的噪聲不敏感，有很好的魯棒性。此外，它的學(xué)習(xí)過(guò)程能最大限度地利用樣本和樣本之間的關(guān)系，減少了類(lèi)別特征選擇不當(dāng)對(duì)分類(lèi)造成的不利影響[10]。

3 實(shí)驗(yàn)結(jié)果及討論

實(shí)現(xiàn)氣味識(shí)別的主要步驟為:(1)對(duì)不同氣味刺激中的嗅球僧帽細(xì)胞層信號(hào)進(jìn)行低通濾波，保留200 Hz以下的低頻場(chǎng)電位信號(hào)。(2)對(duì)低頻場(chǎng)電位信號(hào)使用多窗譜估計(jì)的方法計(jì)算功率譜。(3)選擇合適的移動(dòng)窗和移動(dòng)步長(zhǎng)建立功率譜時(shí)頻圖。(4)對(duì)不同氣味刺激產(chǎn)生的功率譜時(shí)頻圖分成訓(xùn)練集與測(cè)試樣本集，使用K最鄰近分類(lèi)算法識(shí)別測(cè)試樣本集對(duì)應(yīng)的刺激氣味。

實(shí)驗(yàn)中，每種氣味刺激記錄的數(shù)據(jù)時(shí)長(zhǎng)為30 s，其中包括氣味刺激前10 s，連續(xù)氣味刺激過(guò)程4 s和氣味刺激結(jié)束后16 s的僧帽細(xì)胞層場(chǎng)電位信號(hào)。圖2分別是異丁醇、苯甲醚、香芹酮和檸檬醛4種氣味刺激相關(guān)的30 s嗅球僧帽細(xì)胞層低頻信號(hào)的功率譜估計(jì)?？梢钥吹剑啻胺ㄏ啾扔趥鹘y(tǒng)的周期圖法在功率譜估計(jì)結(jié)果上更顯平滑。分析功率譜的能量分布，信號(hào)能量主要集中在低頻段，且4種氣味刺激的能量分布差異度不大。高能量的低頻信號(hào)掩蓋了氣味刺激過(guò)程中在時(shí)域中某些頻率段能量出現(xiàn)的變化情況。因此，使用30 s信號(hào)的譜估計(jì)僅能給出信號(hào)全局功率譜特征，沒(méi)有凸顯出不同氣味刺激下短時(shí)間內(nèi)的功率譜能量變化的特征。

圖2 4種氣味刺激下嗅球信號(hào)的多窗法譜估計(jì)

為了放大氣味刺激后的局部功率譜特征，我們采用多窗法功率譜估計(jì)時(shí)頻圖分析時(shí)長(zhǎng)為4 s的氣味刺激過(guò)程數(shù)據(jù)。首先對(duì)分析數(shù)據(jù)段選取移動(dòng)時(shí)間窗寬為0.5 s，步長(zhǎng)為 0.05 s，正交數(shù)據(jù)窗數(shù)量選擇為5個(gè)，然后使用基于移動(dòng)窗的多窗法譜估計(jì)，得到與時(shí)間、頻率信息相關(guān)的信號(hào)功率譜矩陣。值得注意的是，正交數(shù)據(jù)窗數(shù)量選擇越多，功率譜估計(jì)越平滑，但得到的頻率分辨率越低。以時(shí)間為橫坐標(biāo)，頻率為縱坐標(biāo)，可以得到功率譜能量隨時(shí)間變化的時(shí)頻圖。功率譜時(shí)頻圖可理解為每個(gè)時(shí)間點(diǎn)都進(jìn)行了數(shù)據(jù)長(zhǎng)度為0.5 s的多窗法譜估計(jì)來(lái)計(jì)算功率譜，并且選取了正交數(shù)據(jù)窗數(shù)量的平均值作為該時(shí)間點(diǎn)上功率隨頻率的分布。由于嗅球僧帽細(xì)胞層的信號(hào)在低頻段(0～20 Hz)能量較大，為顯示氣味刺激相關(guān)的功率譜變化信息，每一幅功率譜圖都選取了氣味刺激前的3個(gè)時(shí)間窗(1.5 s)數(shù)據(jù)平均值作為基線信號(hào)用于該段功率譜的規(guī)范化處理。圖3所示為在4種氣味刺激下典型的功率譜時(shí)頻圖。圖中的時(shí)間信息提高了氣味刺激后4 s內(nèi)嗅球僧帽細(xì)胞層響應(yīng)信息的分辨率。從圖中觀察發(fā)現(xiàn)，氣味刺激過(guò)程中功率譜能量改變主要集中在gamma頻段(40 Hz～120 Hz)，生理研究也表明氣味刺激主要引起gamma頻段振蕩[11]。圖3中顯示的4幅時(shí)頻圖表明對(duì)于4種氣味刺激在gamma頻段能量分布略有不同，異丁醇能量主要分布在60 Hz～70 Hz，苯甲醚能量主要分布在100 Hz～120 Hz，香芹酮能量主要分布在40 Hz～50 Hz和80 Hz～90 Hz，檸檬醛能量主要分布在0.3 Hz～5 Hz和40 Hz～50 Hz。在不同次氣味刺激中，功率譜能量變化出現(xiàn)的時(shí)間點(diǎn)略有差異，但能量主要分布頻率范圍基本一致。

圖3 4種氣味刺激下的功率譜時(shí)頻圖

根據(jù)以上分析得到的功率譜能量分布特征，我們對(duì)功率譜時(shí)頻圖按頻率分為δ(0.3 Hz～5 Hz)，θ－α(5 Hz～15 Hz)，β(15 Hz～40 Hz)，γ(40 Hz～120 Hz)和γ2(120 Hz～200 Hz)共5個(gè)子頻段分別進(jìn)行氣味識(shí)別。在120組氣味刺激功率譜時(shí)頻圖中對(duì)每種氣味刺激選取15組時(shí)頻圖，使訓(xùn)練集達(dá)到60組向量，剩余另60組向量作為測(cè)試集。采用二重交叉檢驗(yàn)的檢驗(yàn)的方法，將測(cè)試集與訓(xùn)練集交換計(jì)算平均識(shí)別準(zhǔn)確率。如圖4所示，采用K最鄰近分類(lèi)法能夠在γ頻段實(shí)現(xiàn)較好地氣味識(shí)別，得到的總體氣味分辨準(zhǔn)確率為77.4%，高于其他子頻段的識(shí)別準(zhǔn)確率，這與氣味刺激后信號(hào)在gamma頻段的功率譜能量分布較大且差異性較明顯結(jié)果一致。對(duì)于每種單一氣體在gamma頻段的識(shí)別率，結(jié)果顯示異丁醇90%、苯甲醚83.3%的識(shí)別準(zhǔn)確率相對(duì)高于檸檬醛66.7%、香芹酮70%。單一氣味的識(shí)別準(zhǔn)確率可能與嗅球中氣味感受域有一定的關(guān)系[12]。當(dāng)氣味的感受域較大時(shí)，在嗅球中響應(yīng)細(xì)胞分布更為廣泛，功率譜能量變化較為明顯，更容易通過(guò)場(chǎng)電位信號(hào)進(jìn)行氣味識(shí)別。

圖4 功率譜各子頻段用于氣味識(shí)別的準(zhǔn)確率統(tǒng)計(jì)

由于嗅球低頻場(chǎng)電位信號(hào)穩(wěn)定且幅值相對(duì)較大，便于獲得較高信噪比的信號(hào)，故被記錄用于分析。實(shí)驗(yàn)設(shè)計(jì)中刺激氣味選擇了4種具有不同官能團(tuán)的小分子氣體，希望能得到差異度較大的信號(hào)。在信號(hào)解讀方面，時(shí)域信號(hào)較難得到氣味刺激相關(guān)的特征性信息。多窗譜估計(jì)法用于生物信號(hào)的功率譜分析具有一定優(yōu)勢(shì)[13]，特別是加入了移動(dòng)窗分析后能體現(xiàn)出功率譜隨時(shí)間變化的具體信息。移動(dòng)窗寬0.5 s的選取與大鼠呼吸節(jié)律一致，一般認(rèn)為這個(gè)時(shí)間窗精度可以包含嗅覺(jué)信號(hào)中的各頻段振蕩信息[14]。功率譜時(shí)間信息加入氣味識(shí)別中能引入更多特征信息，但由于生物信號(hào)本身存在著一定的時(shí)變特性，是否能顯著地提高氣味識(shí)別率仍有待進(jìn)一步分析。利用功率譜時(shí)頻圖γ頻段的識(shí)別結(jié)果達(dá)到基本識(shí)別要求。為提高識(shí)別精度，可引入單細(xì)胞鋒電位變化信號(hào)協(xié)同解碼[15]。同時(shí)，也可以嘗試其他分類(lèi)識(shí)別算法如支持向量機(jī)、偏最小二乘回歸等，實(shí)現(xiàn)嗅覺(jué)傳感系統(tǒng)氣味的準(zhǔn)確識(shí)別。

4 結(jié)論與展望

本文利用哺乳動(dòng)物嗅覺(jué)系統(tǒng)中的嗅球作為氣味信息的生物傳感單元，采用微電極陣列傳感器植入的手段記錄嗅球中的生物電信號(hào)并解讀，初步實(shí)現(xiàn)氣味識(shí)別，是生物學(xué)與工程學(xué)結(jié)合的一種氣體傳感系統(tǒng)設(shè)計(jì)。目前氣味識(shí)別較多采用化學(xué)傳感器陣列組成電子鼻的方法，構(gòu)建以生物嗅覺(jué)系統(tǒng)為初級(jí)傳感單元的氣味識(shí)別系統(tǒng)具有新穎性，同時(shí)也具有一定的挑戰(zhàn)性。利用生物嗅覺(jué)系統(tǒng)對(duì)于氣味識(shí)別所特有的高特異性和高靈敏性可以構(gòu)造出更具有吸引力的新一代生物電子鼻。另一方面，由于生物系統(tǒng)的不穩(wěn)定性以及氣味相關(guān)信號(hào)的解碼仍面臨著一些困難，因此該嗅覺(jué)傳感系統(tǒng)仍需要進(jìn)行進(jìn)一步的探索和改進(jìn)，以實(shí)現(xiàn)氣味更準(zhǔn)確、更有效地識(shí)別。

[1]Laurent G.A Systems Perspective on Early Olfactory Coding[J].Science，1999，286:723－728.

[2]Rosero M A，Aylwin M L.Sniffing Shapes the Dynamics of Olfactory Bulb Gamma Oscillations in Awake Behaving Rats[J].Eur.J.Neurosci.，2011，34:787－799.

[3]Kay L M，Beshel J.A Beta Oscillation Network in the Rat Olfactory System During a 2-Alternative Choice Odor Discrimination Task[J].J.Neurophysiol.，2010，104:829－839.

[4]平建峰，吳堅(jiān)，應(yīng)義斌.基于短時(shí)傅立葉變換的雞蛋破損檢測(cè)技術(shù)的研究[J].傳感技術(shù)學(xué)報(bào)，2009，22(7):1055－1060.

[5]Thomson D J.Spectrum Estimation and Harmonic Analysis[J].Proc.IEEE，1982，70:1055－1096.

[6]楊廣映，羅志增.肌電信號(hào)的功率譜分析方法[J].傳感技術(shù)學(xué)報(bào)，2004，3:355－358.

[7]Frederik J S，Maria T Z，Jun K.Isostatic Response of the Australian Lithosphere:Estimation of Effective Elastic Thickness and Anisotropy Using Multitaper Spectral Analysis[J].J.Geophys.Res.，2000，105:19163－19184.

[8]Rinberg D，Koulakov A，Gelperin A.Sparse Odor Coding in Awake Behaving Mice[J].J.Neurosci.，2006，26:8857－8865.

[9]武鵬鵬，趙剛，鄒明.基于多窗譜估計(jì)的改進(jìn)譜減法[J].科學(xué)計(jì)算及信息處理，2008，12:150－152.

[10]張寧，賈自艷，史忠植.使用KNN算法的文本分類(lèi)[J].計(jì)算機(jī)工程，2005，31(8):171－173.

[11]Cenier T，Amat C，Litaudon P，et al.Odor Vapor Pressure and Quality Modulate Local Field Potential Oscillatory Patterns in the Olfactory Bulb of the Anesthetized Rat[J].Eur.J.Neurosci.，2008，27:1432－1440.

[12]Luo M，Katz L C.Response Correlation Maps of Neurons in the Mammalian Olfactory Bulb[J].Neuron，2001，32(6):1165－1179.

[13]Mollazadeh M，Aggarwal V，Singhal G，et al.Spectral Modulation of LFP Activity in m1 During Dexterous Finger Movements[J].Conf.Proc.IEEE Eng.Med.Biol.Soc.，2008，5314－5317.

[14]Cindy P，Jeffry S I.Odor Representations in Olfactory Cortex:“Sparse”Coding，Global Inhibition，and Oscillations[J].Neuron，2009，62(6):850－861.

[15]Zhou J，Dong Q，Zhuang L J.Odor Discrimination by Mitral Cells in Rat Olfactory Bulb Using Microwire Array Recording[J].J.Innovative Opt.Health Sci.，2012;5(1):97－103.