基于在体植入电极的生物嗅觉传感系统设计及其气味识别*

周 俊，董 琪，庄柳静，吴春生，刘清君，王 平

(浙江大学生物医学工程与仪器科学学院，生物传感器国家专业实验室，生物医学工程教育部重点实验室，杭州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.