图1 ICA处理流程
Fig.1 Process of ICA method
摘 要：针对提取局部场电位(LFP)用于运动意图解码的特征时，存在LFP信噪比低、编码时间窗难以确定等问题，提出了一种结合独立成分分析(ICA)与小波分解的特征提取方法，用于动物转向行为的神经信息解码. 首先结合动物运动行为视频与LFP信号时频分析方法，确定编码时间窗的范围；然后用ICA对时间窗内的LFP进行去噪处理，提高LFP信噪比；接着利用小波分解进一步确定LFP编码频带，并通过滑窗方法计算频带内的时序能量，构建编码特征；最后采用k近邻方法对编码特征进行分类，验证其解码性能. 实验结果表明，利用提出的特征提取方法，经过1 000次交叉互验证，分类正确率达到(92.35±5.87)%，能够准确稳定地解码动物的转向行为.
关键词：局部场电位；独立成分分析；小波分解；时频分析
局部场电位(Local Field Potential, LFP)是由大量神经元突触活动引起的，反映了大脑局部区域神经元集群的综合电活动，具有较好的长期解码稳定性，而将其应用于运动解码的研究具有很大的应用潜力，正逐渐成为神经科学研究新的热点．Arjun等[1]用在猴子初级运动皮层和腹侧前运动皮层采集到的低频局部场电位信号，解码出了猴子抓握运动轨迹．Tomislav等[2]对猴子抓握运动参数进行解码研究发现，在运动计划和执行阶段LFP多个频带共同编码了运动信息．然而，LFP十分微弱，而且具有较强的非平稳特性，难以提取出有效的解码特征，因此确定出编码时间窗和编码频带，进而提取出可靠的编码特征是准确解码的关键．
在对LFP进行特征提取时，常用50 Hz陷波器陷波处理的方法对工频干扰进行去除，之后用经典滤波器滤波得到特定频带的信号，再提取能量特征，然而处理效果并不理想．龙飞等[3]将独立成分分析(Independent Component Analysis, ICA)方法应用到脑电中眼电伪迹的去除，消噪的同时能对有用信号的细节成分做到很好的保留，使目标信息的特征更为明显，可以进一步对相关信息进行准确提取．钟伯成等[4]利用小波变换对脑电信号中的瞬态弱电信号进行了有效提取．Camara等[5]利用小波变换对从帕金森病人身上采集到的LFP进行特征提取，实现了对不同类型病例的准确分类．LFP与皮层脑电信号较为相似，同为频率较低的非平稳信号，ICA与小波方法对LFP的特征提取具有很大的应用潜力．
笔者以鸽子为实验对象，采集了其转向运动时与运动决策相关的NCL(Nidopallium Caudolaterale)区[6]LFP信号．结合ICA和小波变换方法对特征进行提取，进而用k近邻[7](k-Nearest Neighbor, kNN)方法进行分类，并与常用的LFP特征提取方法进行了对比．
1.1 独立成分分析
ICA[8]是一种盲源分离算法，是指在不知道源信号和传输通道参数的情况下根据输入源信号的统计特性，仅由观测信号恢复出源信号各个独立成分的过程．ICA方法的处理流程如图1所示．
观测信号X=AS，在混合矩阵A和原始信号未知的情况下，只根据观测数据向量X确定分离矩阵W，使得变换后的输出Y=WX．
图1 ICA处理流程
Fig.1 Process of ICA method
1.2 小波分解
小波分解对信号具有很好的时域和频域分辨率，常用于非平稳信号的时频特性分析．用Mallat算法[9]对信号f(t)进行离散小波分解，可表示为：
(1)
式中：L为分解层数；Dj为不同尺度下的细节分量；AL为近似分量．如果信号采样率为fs，则AL,DL,DL-1,…,D1分量对应的子频段为[0,fs/2L+1],[fs/2L+1,fs/2L],[fs/2L,fs/2L-1],…,[fs/22,fs/2]．
1.3 ICA-小波特征提取
为了有效提取用于解码鸽子转向行为的LFP特征，提出一种结合ICA与小波的特征提取方法．从神经信息处理的角度出发，采用基于负熵最大的FastICA算法[10]．它以负熵最大作为一个搜寻方向，实现顺序地提取独立源，该算法采用了定点迭代的优化算法，使得收敛更加快速、稳健．
FastICA算法的基本步骤如下：
(1)对数据进行去均值、白化处理，将观测数据X变为Z．
(2)确定要分离的独立分量个数，一般与观测数据的数目相同，随机选取初始矢量Wp．
(3)由负熵最大的原则推出Wp的迭代公式
(2)
式中：Wp为第p个要估计的分离矩阵的矢量
是Wp的新值;Z为白化后的数据;g为选取的非线性函数，笔者选取为正切函数．
(4)估算出分离矩阵，FastICA算法结束．
根据局部场电位信号的特征，选取与局部场电位较为相似的Daubechines小波[11]． ICA-小波特征提取的流程如图2所示，其具体步骤如下：
(1)用ICA方法分解原始信号，对分解出的独立分量进行分析，作出分解出的各分量信号的功率谱，将明显的干扰信号分量置零，逆变换得到去除噪声的信号．
(2)根据LFP信号的常用采样频率的范围(1 000～2 000 Hz)及可能包含编码信息频带的范围(1～100 Hz)，确定出小波分解的层次；对ICA去噪后的信号用小波方法进行多层小波分解，分解出多个小波子带，对各层的细节系数进行单支重构．由时频分析作出不同方向的时频图，确定特征提取的小波频带．
(3)对选取的m个通道信号小波子带信号，在其转向发生时刻t秒内，滑动时间窗求取该时间窗内的能量值，并对其进行归一化处理， 每个通道提取n个特征，以m×n维特征作为解码特征，之后进行分类．
(4)利用kNN方法对编码特征进行分类.kNN方法不需要估计样本的概率密度函数，找到被分类对象在测试样本中k个最近的样本，然后根据样本分类属性进行投票，将预测值的分类属性赋于被分类对象，解码出鸽子运动方向．
为了获得鸽子转向运动行为的神经信号，笔者采用十字迷宫对动物进行转向训练并同步进行神经信号检测．鸽子转向训练范式如图3所示．图中有标号为1,2,3,4的4对红外装置，用于对同步记录到的神经信号数据的对应时间进行标记．
在电极植入前，首先对鸽子进行训练，训练过程如下：实验开始前将鸽子放入等待区，当实验开始提示音响起后，自动门打开，同时左转、直行、右转三个方向的食槽随机打开一个，鸽子到达十字路口通过观察食槽位置朝相应方向前进，得到食物奖励后返回到起始点，等待下次实验．
图2 ICA-小波编码特征提取及解码流程
Fig.2 Process of feature extracting by ICA-wavelet and decoding
图3 鸽子训练范式及神经信号示意图
Fig.3 Training mode and neural signals
当鸽子转向正确率达到90%以上时，将2×8通道的Microprobe微电极阵列植入到鸽子的NCL脑区，术后7 d左右利用CerebusTM的128通道数据采集系统进行信号采集．采集到的神经信号幅值较小，一般是μV级别的需要经过放大装置进行放大；数据传输采用光纤通信，采样率为2 kHz，16位模数转换精度；采用通带为0～250 Hz的二阶Butterworth滤波器对采集到的信号进行滤波得到LFP，之后降采样至1 kHz，实验全程通过视频监控，以便对鸽子运动状态进行分析．对采集到的信号用去均值的方法去除基线漂移，之后用于特征提取与解码．
为对鸽子转向的LFP进行整体认知，首先结合鸽子运动行为视频与LFP时频分析，确定编码时间窗．观察到信号时频变化如图4的原始信号所示，信号在40～60 Hz的频带内1 s的时间内能量有变化，但是工频干扰比较明显，导致信号特征不是很明显，影响解码效果，需要对信号进一步处理．用50 Hz陷波方法对原始信号进行处理，去除工频干扰的同时，将50 Hz可能包含的信号也去除了．用ICA进行处理，去除噪声的同时且不损害有用信号，使特征更为明显．对3个方向的原始信号用不同方法进行去噪后的时频图如图4所示．
由于信号的能量发生变化，对处理后的信号，提取其能量特征．用经典的Chebyshev滤波器方法提取7个通道1 s内40～60 Hz的信号，每200 ms为一个时间窗，计算该段时间内的能量值，共得到35个能量值，将其归一化(均值为0，方差为1)后为作为分类的特征向量．小波方法具有很好的局部特征提取作用，对ICA处理后的数据用db7小波进行7层小波分解，提取第4层信号(频带范围为31.25～62.5 Hz)，同样的处理方法得到分类的特征向量，各方向各通道变化如图5所示．可以看出左转，直行，右转信号的能量有所差异，小波方法提取的特征可分性更强．取左转15组，直行15组，右转15组信号，对其进行ICA-小波特征提取，左转、直行、右转各方向的信号能量变化又具有一致性，如图6所示．
图4 ICA和陷波器去噪效果对比时频图
Fig.4 Comparison of de-noising effect time-frequency diagram by ICA and notch filter method
为了评估所提出算法的解码效果，笔者比较了Chebyshev滤波器特征提取算法(原始信号陷波滤波器去除工频干扰后用Chebyshev滤波器滤波得到40～60 Hz频带)，小波滤波特征提取，ICA-Chebyshev滤波特征提取，和ICA-小波滤波特征提取四种特征提取算法的解码效果．采用kNN方法对信号特征进行分类．对4只鸽子采集的数据，用2/3的数据用作训练，1/3的数据用于测试分类，进行1 000次测试，其分类结果如表1所示．
从表1可以看出，笔者提出的方法对运动转向的分类效果最好．而经过ICA方法处理，使得原本混合在同一通道的相互独立的信号源分离开，不同通道采集到的同一信号源发出的信号加强，将明显的干扰信号去除，使解码时的左转、直行和右转的可分性更强．小波分析方法具有良好的时频特性，使得我们在时频域获得更多的信息，对于非平稳信号的特征提取具有很好的适用性．而Chebyshev等数字滤波器，容易造成频谱泄露，而且过渡带较宽，时频分辨率较差．经过小波函数提取的特征能明显反应不同转向时的特征变化，达到了很好的解码效果．
图5 Chebyshev与小波方法提取的能量特征对比
Fig.5 Comparison of energy features extracted by Chebyshev and wavelet method
图6 不同转向多个试次的特征一致性分析
Fig.6 Consistency analysis of the features in different direction
表1 ICA-小波与其它方法提取特征解码正确率
Tab.1 Decoding accuracy of extracted features by ICA-wavelet and other methods %
鸽子编号Chebyshev方法解码正确率均值解码正确率标准差小波方法解码正确率均值解码正确率标准差ICA-Chebyshev方法解码正确率均值解码正确率标准差ICA-小波方法解码正确率均值解码正确率标准差P60083.478.4486.467.2688.286.9292.355.87P60972.0714.5373.3814.1776.3213.5882.6211.60P61965.736.3866.836.2670.078.6674.018.33P60366.7510.1468.709.0677.647.3580.555.12
笔者采集到动物转向时的信号，对其进行视频分析、去噪、时频分析、提取特征，并进行分类解码．结果表明，鸽子进行运动转向时，其LFP发生了变化，说明LFP与运动决策相关，可以用于对运动方向的解码．利用ICA方法去除信号的噪声，提高了信号的信噪比，通过对数据进行时频分析，得到时频域的变化信息；并将不同实验场景发生改变的信号准确定位，确定出特征提取的时间窗与编码频带，进而结合小波分析的方法对特征进行提取．经过ICA方法处理后再用小波方法进行提取的特征，能解码鸽子运动转向行为，解码正确率较其他方法有较大提高，而且分类方差较小，证明了该方法的有效性．
参考文献：
[1] ARJUN K, CARLOS E, WILSON T. Relationships among low-frequency local field potentials, spiking activity, and three-dimensional reach and grasp kinematics in primary motor and ventral premotor cortices [J]. Journal of neurophysiology, 2011, 105(4): 1603-1619.
[2] TOMISLAV M, WILSON T, SONJA G. Local field potentials in primate motor cortex encode grasp kinetic parameters[J]. NeuroImage, 2015, 114(1): 338-355.
[3] 龙飞,吴小培,范羚. 基于独立分量分析的脑电噪声消除[J]. 生物医学工程学杂志, 2003, 20(3): 479-483.
[4] 钟伯成,吴小培. 基于小波变换的脑电信号瞬态特征提取[J]. 模式识别与人工智能, 2000, 13(2): 218-221.
[5] CAMARA C, LSASI P, WARWICK K, et al. Resting tremor classification and detection in Parkinson’s disease patients [J]. Biomedical signal processing and control, 2015, 16: 88-97.
[6] DANIEL L, ROLAND P, ONUR G. Neurons in the pigeon nidopallium caudolaterale signal the selection and execution of perceptual decisions [J]. European journal of neuroscience, 2014, 40(9): 3316-3327.
[7] 刘新玉,尚志刚,万红. 神经元锋电位检测中大幅值干扰的去除[J]. 中国科技论文, 2013, 8(1):46-50.
[8] COMON P. Independent component analysis, a new concept? [J]. Signal processing, 1994, 36(3):287-314.
[9] SHENSA M. The discrete wavelet transform: wedding the a trous and mallat algorithms[J]. IEEE transactions on signal processing, 1992, 40(10): 2464-2482.
[10]HYYARINEN A. Fast and robust fixed-point algorithms for independent component analysis [J]. IEEE Transactions on neural networks, 1999, 10(3): 626-634.
[11]尚志刚,冯平艳,刘新玉,等. 局部场电位γ频带能量对朝向调谐特性研究[J]. 郑州大学学报(工学报), 2012, 33(6): 5-9.
Abstract：In order to overcome the problems as low signal-to-noise ratio (SNR) of LFP and difficulty in identifying encoding time window when extract the features of motion intention, a method that combines independent component analysis (ICA) with Wavelet was presented to extract the features of turning. Firstly, the motion videos of animals were analyzed and the time-frequency diagrams of LFP were plotted to determine the time window of signal which encoded the motion information. Then, ICA was used to increase the SNR of LFP. Thirdly, the encode bands of LFP were extracted by wavelet method as well as the encode features were extracted by sliding time window method. Lastly, k-nearest neighbor method was used to classify the encode features. And via 1 000 times cross validation the precision was(92.35±5.87)%, the results showed that it could decode reliably the motion intention of animals.
Key words：local field potentials; independent component analysis; wavelet decomposition; time-frequency analysis
收稿日期：2016-07-01；
修订日期：2016-09-07
基金项目：国家自然科学基金资助项目(U1304602)；河南省科技攻关计划项目(122102210102)
文章编号：1671-6833(2017)03-0039-05
中图分类号：R318.04
文献标志码：A
doi：10.13705/j.issn.1671-6833.2016.06.006