快速有效地提取和比较不同神经信号中所含信息的相同和区别之处是研究者关注的问题。本文介绍了一种基于小波变换的时频建模方法。该方法用少量的半椭球模型来表征神经信号在时域和频域上的变化,克服了传统时频分析中背景干扰大、参数多的缺陷。将该方法应用于嗅球场电位的研究,与支持向量机(SVM)算法相结合,可初步实现气味的分类。
引用本文： 董琪, 胡靓, 庄柳静, 周俊, 王平. 基于小波变换的时频图建模及其在嗅觉场电位分析中的应用. 生物医学工程学杂志, 2014, 31(3): 481-486. doi: 10.7507/1001-5515.20140089 复制
引言
近年来，随着脑-机接口技术的发展，低频神经信号(＜200 Hz)，包括非侵入式检测(如脑电图，electroencephalogram，EEG)[1] 和侵入式检测(如局部场电位，local field potential，LFP)[2]，得到越来越多研究者的重视。和单神经元放电(spike)信号相比，这些低频神经信号反映了一定区域内所有神经元电生理活动的总体特征，因此具有更高的稳定性，并被广泛应用于神经解码、信号分析、疾病检测等领域 [3-8]。
时频联合分析是上世纪90年代发展起来的神经信号处理方法。该方法能够表征神经信号的频率成分组成以及相关频率成分的能量分布随时间的变化，从而克服了传统时域和频域分析不能完整准确地表征信号的缺陷[1]，在各种神经信号研究中得到了大量的应用[5-7, 9]。在EEG和LFP的分析中，主要的难点在于：(1)与事件相关的信息往往被淹没在强大的背景噪声中，很难从原始信号中有效提取有用特征[3]；(2)基于短时傅里叶变换等手段的传统时频分析方法会产生大量的参数，不利于统计归纳。
本文介绍了一种新的低频神经信号的处理方法。该方法将小波时频图和椭球方程建模方法有机地结合起来，用时频图表征神经信号的时频特征，并通过建模的方式从杂乱的背景噪声中有效提取其中的有用信息，大大减少了时频图原有的参数量，在一定程度上克服了传统时频分析中存在的缺陷。将该方法应用于大鼠嗅球LFP的研究，可提取与气味刺激相关的事件特征，与支持向量机(support vector machine,SVM)算法相结合，初步实现了三种气味刺激的分类。
1 分析方法
该方法的核心思想是用尽可能少的方程参数来表征原始时频图，以达到去除噪声干扰，提取有用信息供进一步分析的目的。主要包括两个步骤：① 构建神经信号基于小波变换的功率谱时频图；② 使用半椭球方程对时频图进行建模。
1.1 基于小波变换的功率谱时频图
小波变换是以某些特殊函数为基础，将数据过程或数据系列变换为级数系列，以发现它的类似频谱的特征，从而实现数据处理。这种方法是由法国工程师Morlet 在1974 年首先提出，经过之后数十年的研究，该方法建立了完整的数学化体系和扎实的理论基础，并被广泛应用于各种神经生理信号的分析[1-2, 10-13]。函数f(x)的连续小波变换可定义为
| ${{W}_{f}}\left( a,b \right)=\frac{1}{\sqrt{\left| a \right|}}\int_{-\infty }^{+\infty }{f\left( x \right)\bar{\Psi }}\left( \frac{x-b}{a} \right)dx,$ |
式中Ψ(x)是母小波，(x)是Ψ(x)的共轭形式，参数a 和b 分别被称为尺度因子和平移因子。由式(1)可知，连续小波变换是空间(时间)和尺度的局部变换，因而能有效地从信号中提取信息。通过伸缩和平移等运算功能可对函数或信号进行多尺度的细化分析，解决了傅里叶变换不能解决的许多困难。
由连续小波变换得到的小波功率谱可定义为
| $TF\left( a,b \right)={{\left| {{W}_{f}}\left( a,b \right) \right|}^{2}}$ |
式中TF(a,b)是函数f(x)在尺度a和位移b时所具有的能量。严格地说，连续小波变换提供了时间尺度分析，而不是时间频率分析。然而，使用适当的尺度频率变换，依然可以得到一个近似于时间频率分析的工具。与短时傅里叶变换相比，小波时频图的窗函数具有可变的特性。对于较高的频率，时间窗会变小，使时间分析精度提高，而对于较低的频率，时间窗会变大，保留了完整的低频信息。
小波变换的关键是要选择合适的母小波[1]。复Morlet小波是一种使用广泛的基本小波,已被众多研究证明，适用于电生理信号的时频分析[1-3, 10-11, 13]。复Morlet小波w(t,f)对于某频率的表达式为
| $w\left( t,f \right)=Aexp(2i\pi ft)exp({{t}^{2}}/2{{\sigma }_{i}}^{2})$ |
式中A为归一化因子，其目的是使小波能量为1；f与σt(时域分辨率)是需要适当选择的参数，两者的乘积决定了小波的周期数，同一复Morlet小波家族具有恒定的f与σt乘积，本研究所使用的小波家族定义为2πfσt=7(见图 1a)。复Morlet小波w(t,f)是一种复值调制的高斯(Gaussian)函数，在时域和频域上都具有高斯分布，其傅里叶频谱的标准差满足
| ${{\sigma }_{f}}=\frac{1}{2\pi {{\sigma }_{t}}},$ |
特别的，在本研究中，
| ${{\sigma }_{f}}=\frac{1}{2\pi {{\sigma }_{t}}}=\frac{f}{7},$ |
因此，不同频率f所对应的σt和σf不相同，即在整个时频平面上具有可变的时频分辨率：在高频区能提供高的时间分辨率，在低频区能提供高的频率分辨率(见图 1b)。
图1
复Morlet小波的特性
(a)复Morlet小波的实部(公式3)，左：f/
(a) real part of equation 3,with
由于神经信号在低频段具有的能量比高频段大很多，为了在同一比例下观察场电位信号在低频段和高频段能量随时间的变化，由连续小波变换获得的信号功率时频图需要标准化
| ${{z}_{ft}}=\frac{{{c}_{ft}}-{{m}_{f}}}{{{\delta }_{f}}},$ |
式中cft为f(t)在频率f及时间t时的小波功率谱，mf和δf分别为刺激前1 s内信号在频率f的小波功率谱的平均值和标准差。
图 2为一段人造信号的原始波形及其基于复Morlet小波变换的时频图。该人造信号(见图 2a)模拟了神经信号在各频率段的成分组成，具有1个5 Hz的基础振荡，并在0.5~0.7 s和1.5 ~1.7 s范围内分别具有一个β段(15~40 Hz)和σ段(40~90 Hz)的突变。图 2(b)为该信号的原始时频图，由于低频段基础振荡的能量大，使得高频段的突变不明显，经公式(6)归一化之后，能量的突变点被凸显出来，而恒定的能量波段被弱化。
图2
人造信号及其小波时频图
(a)人造信号原始波形；(b)人造信号原始小波时频图；(c)经过归一化之后的小波时频图
Figure2. An artificial signal and its wavelet time-frequency maps(a) the waveform of an artificial signal； (b) the raw wavelet time-frequency map； (c) the normalized wavelet time-frequency map
1.2 基于半椭球方程的时频图建模
原始时频图包含的数据量非常大，很难直接提取信号的时频特征。在本研究中，采用了建模的方法进行信息提取和降纬处理。该方法的核心思想是利用事先预定的参数方程φ对原始时频图进行拟合，主要包含时频图窗口化、求和与拟合三个步骤。
首先将原始时频图分解成一系列重叠的窗口，以确定需要建模的区域。在经行窗口化时，窗口的尺寸必须满足三点要求：① 由于原始时频图是通过连续小波变换获得，在进行窗口化时必须使窗口的尺寸能够自动适应小波变换的时频分辨率，即L/σt=H/σf。其中L和H分别为窗口的长度和高度，σt和σf分别为小波变换的时域标准差和频域标准差。 ② 要选择合适的窗口宽度，在充分体现频率特征的基础上凸显能量的突变，一般选择的窗口需要包含当前频率3~5 个振荡周期。计算中，我们选取的值为4。 ③ 为了使位于不同位置的窗口具有相同的评判标准，窗口的面积必须相同。
因此，对于特定的频率f，窗口的长度L=4/f，即一个窗口包含当前频率的4个振荡周期，则L/σt=8π/7，由此换算出窗口的高度H=8πf/49。用这个值去计算原始时频图上每个点，就能获得一系列面积为3249π的重叠窗口。然后逐一计算每个窗口中的能量总和，窗口总能量值体现了该窗口内场电位的振荡强度。对窗口wi，总能量值Si计算公式如式7所示，然后找出所有大于预设阈值的窗口。
| ${{S}_{i}}=\sum\limits_{t,f\in {{w}_{i}}}{{{z}_{ft}}}$ |
最后，在选中的窗口中使用参数方程进行拟合。本研究选择半椭球方程φ作为拟合方程并定义为
| $\varphi \left( f,t \right)=\left\{ \begin{align} & a\sqrt{1-v},0 ＜v＜1 \\ & 0,v>1 \\ \end{align} \right.$ |
其中v=(f/μf)2/lf2+(t/μt)/lt2，μf和μt为该椭球在时频图上的位置，a为椭球的极半径(幅度)，lf和lt为半椭球的赤道半径。图 3为一个半椭球方程的示意图。该方程共有5个参数，并且需要满足以下条件：
(1)μf ＞ 0，μt ＞ 0，且使半椭球位于所选窗口的中心；
(2)0 ＜ lt ＜ L,0 ＜ lf ＜ H。
图3
由5个参数进行表征的半椭球方程示意图
Figure3.
Half-ellipsoid model with its five parameters
在拟合的过程中，要使建模误差函数C达到最小，函数C定义为
| $C=\frac{1}{2}\sum\limits_{t,f\in {{w}_{i}}}{({{z}_{ft}}-\varphi \left( f,t \right)),}$ |
图 4为采自麻醉大鼠在香芹酮刺激下嗅球的局部场电位β段信号的小波时频图以及建模后的结果。从图中我们可以发现，用5个半椭球模型就能将原始时频图的主要能量变化从大量参数和背景噪声中提取出来。
图4
嗅球场电位信号β段时频图和半椭球方程建模结果
左：归一化后的小波时频图；中：半椭球建模的2D表示；右：半椭球建模的3D表示
Figure4. Time-frequency map (left) of β range of signal from rat’s olfactory bulb and its modeling results in two (middle) and three (right) dimensions2 大鼠嗅球场电位记录和气味识别分析
气味分子的识别是由嗅黏膜上的一千多种气味受体完成的，一种气味受体只能与特定分子结构的气味相结合[15-17]。表达同种气味受体的感觉神经元投射到嗅球表面的一个或几个嗅小球中，因此特定的气味能够使嗅球产生特定时空特性的活动[18-22]。研究者们发现，气味刺激前后嗅球僧帽细胞层神经电活动十分强烈，LFP振荡会明显增加，其幅值和波形的变化受到气味种类的调控 [9, 23-25]。
2.1 实验设计和信号记录
实验对象SD大鼠购于浙江省医学科学院，质量是200~250 g。大鼠进行电极植入手术流程是在腹腔注射10% 的水合氯醛(4 mL/kg)麻醉后，移除单侧嗅球上覆盖的头骨，利用油压微推进器(Narishige,Japan)将自制的16通道直径为65 μm镍镉合金微丝(AM system,USA)电极阵列，植入到嗅球僧帽细胞层，平均深度约为200~300 μm。由于嗅球僧帽细胞放电幅值较大，能提示电极已植入到目标细胞层，故我们选取气味刺激时能记录到细胞锋电位信号的通道，分析该通道的低频场电位信号。
使用浓度为10 mmol/L的异丁醇、苯甲醚和香芹酮三种气味作为刺激物，单次刺激给予体积为13 mL、时间为5 s的气味，刺激间隔为60 s。多通道信号采集装置为商用OmniPlex数据采集系统(Plexon,USA)，具有20倍的前置放大器，16位的A/D转换器，能以最高频率为40 kHz，5 000倍模拟放大，实现256个通道实时数据采集，并带有用于低频场电位和高频锋电位的频带隔离可编程硬件滤波器。信号采集装置通过网线将采集到的数据发送到计算机，利用Matlab编译环境实现数据的处理和分析。实验的整体示意图如图 5所示。
图5
实验整体示意图
Figure5.
Schematic of the experiment
2.2 特征提取和气味分类
考虑到了小波变换和半椭球方程建模时可能造成的边界效应，我们将记录的嗅球场电位信号分成三个独立的频段分别进行变换。频率的步进分别为0.25 Hz(θ段、β段)和1 Hz(γ段)。1.5~18.5 Hz：用于研究θ段的信号(5~12 Hz)；9~37.5 Hz：用于研究β段的信号(12~30 Hz)；25~85 Hz：用于研究γ段的信号(30~75 Hz)。经过半椭球方程建模之后，原始信号可以用一系列的半椭球的参数来表征。我们定义了以下参数作为分类的特征值：F1，椭球的个数；F2，椭球的最大幅值; F3，椭球的幅值总和。
因此，每个频段能获得3个特征值，每段信号能获得9个特征值。这些特征值经过归一化之后，用于后续的基于SVM算法的气味分类。 本文采用了SVM算法[26]对气味信号行进分类识别。三种气味共218组时频图向量(苯甲醚73组，异丁醇72组，香芹酮73组)，随机选取一半作为样本A，另一半作为样本B，进行分类测试，最终识别结果如表 1所示。二重交叉检验得到的总体准确率为79.4%，其中苯甲醚、异丁醇和香芹酮的检测准确率分别为82.2%、80.6%和75.3%，其准确度和特性性均高于在以往研究中基于短时傅里叶变换的分析和分类方法[9]。
表1
气味分类结果
Table1.
Results of odor discrimination
3 结论
本文介绍了一种新型的低频神经信号的处理方法。该方法将基于小波变换的时频图分析方法与半椭球建模方法有机结合起来，能够快速有效地提取神经信号能量变化的有用信息，也可用于比较大量数据在相同与不同刺激下的相同与不同之处，克服了传统时频分析方法背景噪声大、参数多的缺陷。该方法目前存在的主要问题是建模算法比较复杂，当待处理数据段增长时，对计算机硬件配置要求较高，数据处理时间大幅增加，有待进一步的优化。
将该方法用于麻醉大鼠嗅球场电位的分析中，可用于提取不同气味刺激下，嗅球的电生理响应特征，通过模式识别的方法，初步实现不同气味刺激的分类，其准确率和特异性均高于基于短时傅里叶变换的分析和分类方法。
引言
近年来，随着脑-机接口技术的发展，低频神经信号(＜200 Hz)，包括非侵入式检测(如脑电图，electroencephalogram，EEG)[1] 和侵入式检测(如局部场电位，local field potential，LFP)[2]，得到越来越多研究者的重视。和单神经元放电(spike)信号相比，这些低频神经信号反映了一定区域内所有神经元电生理活动的总体特征，因此具有更高的稳定性，并被广泛应用于神经解码、信号分析、疾病检测等领域 [3-8]。
时频联合分析是上世纪90年代发展起来的神经信号处理方法。该方法能够表征神经信号的频率成分组成以及相关频率成分的能量分布随时间的变化，从而克服了传统时域和频域分析不能完整准确地表征信号的缺陷[1]，在各种神经信号研究中得到了大量的应用[5-7, 9]。在EEG和LFP的分析中，主要的难点在于：(1)与事件相关的信息往往被淹没在强大的背景噪声中，很难从原始信号中有效提取有用特征[3]；(2)基于短时傅里叶变换等手段的传统时频分析方法会产生大量的参数，不利于统计归纳。
本文介绍了一种新的低频神经信号的处理方法。该方法将小波时频图和椭球方程建模方法有机地结合起来，用时频图表征神经信号的时频特征，并通过建模的方式从杂乱的背景噪声中有效提取其中的有用信息，大大减少了时频图原有的参数量，在一定程度上克服了传统时频分析中存在的缺陷。将该方法应用于大鼠嗅球LFP的研究，可提取与气味刺激相关的事件特征，与支持向量机(support vector machine,SVM)算法相结合，初步实现了三种气味刺激的分类。
1 分析方法
该方法的核心思想是用尽可能少的方程参数来表征原始时频图，以达到去除噪声干扰，提取有用信息供进一步分析的目的。主要包括两个步骤：① 构建神经信号基于小波变换的功率谱时频图；② 使用半椭球方程对时频图进行建模。
1.1 基于小波变换的功率谱时频图
小波变换是以某些特殊函数为基础，将数据过程或数据系列变换为级数系列，以发现它的类似频谱的特征，从而实现数据处理。这种方法是由法国工程师Morlet 在1974 年首先提出，经过之后数十年的研究，该方法建立了完整的数学化体系和扎实的理论基础，并被广泛应用于各种神经生理信号的分析[1-2, 10-13]。函数f(x)的连续小波变换可定义为
| ${{W}_{f}}\left( a,b \right)=\frac{1}{\sqrt{\left| a \right|}}\int_{-\infty }^{+\infty }{f\left( x \right)\bar{\Psi }}\left( \frac{x-b}{a} \right)dx,$ |
式中Ψ(x)是母小波，(x)是Ψ(x)的共轭形式，参数a 和b 分别被称为尺度因子和平移因子。由式(1)可知，连续小波变换是空间(时间)和尺度的局部变换，因而能有效地从信号中提取信息。通过伸缩和平移等运算功能可对函数或信号进行多尺度的细化分析，解决了傅里叶变换不能解决的许多困难。
由连续小波变换得到的小波功率谱可定义为
| $TF\left( a,b \right)={{\left| {{W}_{f}}\left( a,b \right) \right|}^{2}}$ |
式中TF(a,b)是函数f(x)在尺度a和位移b时所具有的能量。严格地说，连续小波变换提供了时间尺度分析，而不是时间频率分析。然而，使用适当的尺度频率变换，依然可以得到一个近似于时间频率分析的工具。与短时傅里叶变换相比，小波时频图的窗函数具有可变的特性。对于较高的频率，时间窗会变小，使时间分析精度提高，而对于较低的频率，时间窗会变大，保留了完整的低频信息。
小波变换的关键是要选择合适的母小波[1]。复Morlet小波是一种使用广泛的基本小波,已被众多研究证明，适用于电生理信号的时频分析[1-3, 10-11, 13]。复Morlet小波w(t,f)对于某频率的表达式为
| $w\left( t,f \right)=Aexp(2i\pi ft)exp({{t}^{2}}/2{{\sigma }_{i}}^{2})$ |
式中A为归一化因子，其目的是使小波能量为1；f与σt(时域分辨率)是需要适当选择的参数，两者的乘积决定了小波的周期数，同一复Morlet小波家族具有恒定的f与σt乘积，本研究所使用的小波家族定义为2πfσt=7(见图 1a)。复Morlet小波w(t,f)是一种复值调制的高斯(Gaussian)函数，在时域和频域上都具有高斯分布，其傅里叶频谱的标准差满足
| ${{\sigma }_{f}}=\frac{1}{2\pi {{\sigma }_{t}}},$ |
特别的，在本研究中，
| ${{\sigma }_{f}}=\frac{1}{2\pi {{\sigma }_{t}}}=\frac{f}{7},$ |
因此，不同频率f所对应的σt和σf不相同，即在整个时频平面上具有可变的时频分辨率：在高频区能提供高的时间分辨率，在低频区能提供高的频率分辨率(见图 1b)。
图1
复Morlet小波的特性
(a)复Morlet小波的实部(公式3)，左：f/
(a) real part of equation 3,with
由于神经信号在低频段具有的能量比高频段大很多，为了在同一比例下观察场电位信号在低频段和高频段能量随时间的变化，由连续小波变换获得的信号功率时频图需要标准化
| ${{z}_{ft}}=\frac{{{c}_{ft}}-{{m}_{f}}}{{{\delta }_{f}}},$ |
式中cft为f(t)在频率f及时间t时的小波功率谱，mf和δf分别为刺激前1 s内信号在频率f的小波功率谱的平均值和标准差。
图 2为一段人造信号的原始波形及其基于复Morlet小波变换的时频图。该人造信号(见图 2a)模拟了神经信号在各频率段的成分组成，具有1个5 Hz的基础振荡，并在0.5~0.7 s和1.5 ~1.7 s范围内分别具有一个β段(15~40 Hz)和σ段(40~90 Hz)的突变。图 2(b)为该信号的原始时频图，由于低频段基础振荡的能量大，使得高频段的突变不明显，经公式(6)归一化之后，能量的突变点被凸显出来，而恒定的能量波段被弱化。
图2
人造信号及其小波时频图
(a)人造信号原始波形；(b)人造信号原始小波时频图；(c)经过归一化之后的小波时频图
Figure2. An artificial signal and its wavelet time-frequency maps(a) the waveform of an artificial signal； (b) the raw wavelet time-frequency map； (c) the normalized wavelet time-frequency map
1.2 基于半椭球方程的时频图建模
原始时频图包含的数据量非常大，很难直接提取信号的时频特征。在本研究中，采用了建模的方法进行信息提取和降纬处理。该方法的核心思想是利用事先预定的参数方程φ对原始时频图进行拟合，主要包含时频图窗口化、求和与拟合三个步骤。
首先将原始时频图分解成一系列重叠的窗口，以确定需要建模的区域。在经行窗口化时，窗口的尺寸必须满足三点要求：① 由于原始时频图是通过连续小波变换获得，在进行窗口化时必须使窗口的尺寸能够自动适应小波变换的时频分辨率，即L/σt=H/σf。其中L和H分别为窗口的长度和高度，σt和σf分别为小波变换的时域标准差和频域标准差。 ② 要选择合适的窗口宽度，在充分体现频率特征的基础上凸显能量的突变，一般选择的窗口需要包含当前频率3~5 个振荡周期。计算中，我们选取的值为4。 ③ 为了使位于不同位置的窗口具有相同的评判标准，窗口的面积必须相同。
因此，对于特定的频率f，窗口的长度L=4/f，即一个窗口包含当前频率的4个振荡周期，则L/σt=8π/7，由此换算出窗口的高度H=8πf/49。用这个值去计算原始时频图上每个点，就能获得一系列面积为3249π的重叠窗口。然后逐一计算每个窗口中的能量总和，窗口总能量值体现了该窗口内场电位的振荡强度。对窗口wi，总能量值Si计算公式如式7所示，然后找出所有大于预设阈值的窗口。
| ${{S}_{i}}=\sum\limits_{t,f\in {{w}_{i}}}{{{z}_{ft}}}$ |
最后，在选中的窗口中使用参数方程进行拟合。本研究选择半椭球方程φ作为拟合方程并定义为
| $\varphi \left( f,t \right)=\left\{ \begin{align} & a\sqrt{1-v},0 ＜v＜1 \\ & 0,v>1 \\ \end{align} \right.$ |
其中v=(f/μf)2/lf2+(t/μt)/lt2，μf和μt为该椭球在时频图上的位置，a为椭球的极半径(幅度)，lf和lt为半椭球的赤道半径。图 3为一个半椭球方程的示意图。该方程共有5个参数，并且需要满足以下条件：
(1)μf ＞ 0，μt ＞ 0，且使半椭球位于所选窗口的中心；
(2)0 ＜ lt ＜ L,0 ＜ lf ＜ H。
图3
由5个参数进行表征的半椭球方程示意图
Figure3.
Half-ellipsoid model with its five parameters
在拟合的过程中，要使建模误差函数C达到最小，函数C定义为
| $C=\frac{1}{2}\sum\limits_{t,f\in {{w}_{i}}}{({{z}_{ft}}-\varphi \left( f,t \right)),}$ |
图 4为采自麻醉大鼠在香芹酮刺激下嗅球的局部场电位β段信号的小波时频图以及建模后的结果。从图中我们可以发现，用5个半椭球模型就能将原始时频图的主要能量变化从大量参数和背景噪声中提取出来。
图4
嗅球场电位信号β段时频图和半椭球方程建模结果
左：归一化后的小波时频图；中：半椭球建模的2D表示；右：半椭球建模的3D表示
Figure4. Time-frequency map (left) of β range of signal from rat’s olfactory bulb and its modeling results in two (middle) and three (right) dimensions2 大鼠嗅球场电位记录和气味识别分析
气味分子的识别是由嗅黏膜上的一千多种气味受体完成的，一种气味受体只能与特定分子结构的气味相结合[15-17]。表达同种气味受体的感觉神经元投射到嗅球表面的一个或几个嗅小球中，因此特定的气味能够使嗅球产生特定时空特性的活动[18-22]。研究者们发现，气味刺激前后嗅球僧帽细胞层神经电活动十分强烈，LFP振荡会明显增加，其幅值和波形的变化受到气味种类的调控 [9, 23-25]。
2.1 实验设计和信号记录
实验对象SD大鼠购于浙江省医学科学院，质量是200~250 g。大鼠进行电极植入手术流程是在腹腔注射10% 的水合氯醛(4 mL/kg)麻醉后，移除单侧嗅球上覆盖的头骨，利用油压微推进器(Narishige,Japan)将自制的16通道直径为65 μm镍镉合金微丝(AM system,USA)电极阵列，植入到嗅球僧帽细胞层，平均深度约为200~300 μm。由于嗅球僧帽细胞放电幅值较大，能提示电极已植入到目标细胞层，故我们选取气味刺激时能记录到细胞锋电位信号的通道，分析该通道的低频场电位信号。
使用浓度为10 mmol/L的异丁醇、苯甲醚和香芹酮三种气味作为刺激物，单次刺激给予体积为13 mL、时间为5 s的气味，刺激间隔为60 s。多通道信号采集装置为商用OmniPlex数据采集系统(Plexon,USA)，具有20倍的前置放大器，16位的A/D转换器，能以最高频率为40 kHz，5 000倍模拟放大，实现256个通道实时数据采集，并带有用于低频场电位和高频锋电位的频带隔离可编程硬件滤波器。信号采集装置通过网线将采集到的数据发送到计算机，利用Matlab编译环境实现数据的处理和分析。实验的整体示意图如图 5所示。
图5
实验整体示意图
Figure5.
Schematic of the experiment
2.2 特征提取和气味分类
考虑到了小波变换和半椭球方程建模时可能造成的边界效应，我们将记录的嗅球场电位信号分成三个独立的频段分别进行变换。频率的步进分别为0.25 Hz(θ段、β段)和1 Hz(γ段)。1.5~18.5 Hz：用于研究θ段的信号(5~12 Hz)；9~37.5 Hz：用于研究β段的信号(12~30 Hz)；25~85 Hz：用于研究γ段的信号(30~75 Hz)。经过半椭球方程建模之后，原始信号可以用一系列的半椭球的参数来表征。我们定义了以下参数作为分类的特征值：F1，椭球的个数；F2，椭球的最大幅值; F3，椭球的幅值总和。
因此，每个频段能获得3个特征值，每段信号能获得9个特征值。这些特征值经过归一化之后，用于后续的基于SVM算法的气味分类。 本文采用了SVM算法[26]对气味信号行进分类识别。三种气味共218组时频图向量(苯甲醚73组，异丁醇72组，香芹酮73组)，随机选取一半作为样本A，另一半作为样本B，进行分类测试，最终识别结果如表 1所示。二重交叉检验得到的总体准确率为79.4%，其中苯甲醚、异丁醇和香芹酮的检测准确率分别为82.2%、80.6%和75.3%，其准确度和特性性均高于在以往研究中基于短时傅里叶变换的分析和分类方法[9]。
表1
气味分类结果
Table1.
Results of odor discrimination
3 结论
本文介绍了一种新型的低频神经信号的处理方法。该方法将基于小波变换的时频图分析方法与半椭球建模方法有机结合起来，能够快速有效地提取神经信号能量变化的有用信息，也可用于比较大量数据在相同与不同刺激下的相同与不同之处，克服了传统时频分析方法背景噪声大、参数多的缺陷。该方法目前存在的主要问题是建模算法比较复杂，当待处理数据段增长时，对计算机硬件配置要求较高，数据处理时间大幅增加，有待进一步的优化。
将该方法用于麻醉大鼠嗅球场电位的分析中，可用于提取不同气味刺激下，嗅球的电生理响应特征，通过模式识别的方法，初步实现不同气味刺激的分类，其准确率和特异性均高于基于短时傅里叶变换的分析和分类方法。
