基于在體植入電極的生物嗅覺傳感系統設計及其氣味識別*

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

(浙江大學生物醫學工程與儀器科學學院，生物傳感器國家專業實驗室，生物醫學工程教育部重點實驗室，杭州310027)

生物傳感技術被廣泛應用于感覺、運動系統生理機制研究等領域，具有實時、快速、靈敏、簡便等優點。動物的嗅覺十分靈敏，利用動物嗅覺構建的傳感系統在辨識氣味方面具有獨特的優勢。嗅球是哺乳動物氣味信息編碼的重要器官，它接收來自嗅上皮嗅覺感受細胞的電興奮，而后者將氣味刺激這一化學信號轉換成電信號用于編碼[1]。研究者發現氣味刺激前后嗅球僧帽細胞層神經電活動十分強烈，除了神經元個體的脈沖式放電行為之外，大量神經元動作電位疊加表現出較大的電位波動，稱之為場電位。目前的研究表明，氣味刺激后嗅球內的局部場電位振蕩會明顯增加，其幅值和波形的變化受到氣味種類的調控[2－3]。

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

本文介紹了我們采用微電極陣列傳感器植入大鼠嗅球構建新型嗅覺傳感系統的探索性研究，分析了該傳感系統在不同氣味刺激下嗅球場電位信號的特征，并對采集得到的信號通過多窗譜估計與K最鄰近分類，實現了氣味的初步識別。

1 大鼠嗅覺傳感系統的設計及信號采集

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

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

圖1 嗅覺傳感系統及采集到的嗅球場電位信號

2 基于頻譜的氣味識別算法

2.1 多窗譜估計算法

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

多窗譜定義如下:

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

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

2.2 功率譜時頻圖

通過選擇移動窗寬度與每次譜估計的時間步長，對每個窗的數據進行多窗法(MTM)譜估計，構建以時間為橫坐標，頻率為縱坐標的功率譜時頻圖(Spectrogram)。

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

2.3 K最鄰近分類算法

K最鄰近分類算法(KNN)是一種基于向量空間模型和實例學習的對象有監督的分類方法，其基本思想是:給定一個經過分類的訓練樣本集合，在對新對象進行分類時，首先從訓練樣本中找出與新對象最相似的K個訓練樣本，并根據這K個訓練樣本所屬的類別，判定待定對象所屬的類別。

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

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

3 實驗結果及討論

實現氣味識別的主要步驟為:(1)對不同氣味刺激中的嗅球僧帽細胞層信號進行低通濾波，保留200 Hz以下的低頻場電位信號。(2)對低頻場電位信號使用多窗譜估計的方法計算功率譜。(3)選擇合適的移動窗和移動步長建立功率譜時頻圖。(4)對不同氣味刺激產生的功率譜時頻圖分成訓練集與測試樣本集，使用K最鄰近分類算法識別測試樣本集對應的刺激氣味。

實驗中，每種氣味刺激記錄的數據時長為30 s，其中包括氣味刺激前10 s，連續氣味刺激過程4 s和氣味刺激結束后16 s的僧帽細胞層場電位信號。圖2分別是異丁醇、苯甲醚、香芹酮和檸檬醛4種氣味刺激相關的30 s嗅球僧帽細胞層低頻信號的功率譜估計。可以看到，多窗法相比于傳統的周期圖法在功率譜估計結果上更顯平滑。分析功率譜的能量分布，信號能量主要集中在低頻段，且4種氣味刺激的能量分布差異度不大。高能量的低頻信號掩蓋了氣味刺激過程中在時域中某些頻率段能量出現的變化情況。因此，使用30 s信號的譜估計僅能給出信號全局功率譜特征，沒有凸顯出不同氣味刺激下短時間內的功率譜能量變化的特征。

圖2 4種氣味刺激下嗅球信號的多窗法譜估計

為了放大氣味刺激后的局部功率譜特征，我們采用多窗法功率譜估計時頻圖分析時長為4 s的氣味刺激過程數據。首先對分析數據段選取移動時間窗寬為0.5 s，步長為 0.05 s，正交數據窗數量選擇為5個，然后使用基于移動窗的多窗法譜估計，得到與時間、頻率信息相關的信號功率譜矩陣。值得注意的是，正交數據窗數量選擇越多，功率譜估計越平滑，但得到的頻率分辨率越低。以時間為橫坐標，頻率為縱坐標，可以得到功率譜能量隨時間變化的時頻圖。功率譜時頻圖可理解為每個時間點都進行了數據長度為0.5 s的多窗法譜估計來計算功率譜，并且選取了正交數據窗數量的平均值作為該時間點上功率隨頻率的分布。由于嗅球僧帽細胞層的信號在低頻段(0～20 Hz)能量較大，為顯示氣味刺激相關的功率譜變化信息，每一幅功率譜圖都選取了氣味刺激前的3個時間窗(1.5 s)數據平均值作為基線信號用于該段功率譜的規范化處理。圖3所示為在4種氣味刺激下典型的功率譜時頻圖。圖中的時間信息提高了氣味刺激后4 s內嗅球僧帽細胞層響應信息的分辨率。從圖中觀察發現，氣味刺激過程中功率譜能量改變主要集中在gamma頻段(40 Hz～120 Hz)，生理研究也表明氣味刺激主要引起gamma頻段振蕩[11]。圖3中顯示的4幅時頻圖表明對于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。在不同次氣味刺激中，功率譜能量變化出現的時間點略有差異，但能量主要分布頻率范圍基本一致。

圖3 4種氣味刺激下的功率譜時頻圖

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

圖4 功率譜各子頻段用于氣味識別的準確率統計

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

4 結論與展望

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

[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]平建峰，吳堅，應義斌.基于短時傅立葉變換的雞蛋破損檢測技術的研究[J].傳感技術學報，2009，22(7):1055－1060.

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

[6]楊廣映，羅志增.肌電信號的功率譜分析方法[J].傳感技術學報，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].科學計算及信息處理，2008，12:150－152.

[10]張寧，賈自艷，史忠植.使用KNN算法的文本分類[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.