2.793

2018影响因子

(CJCR)

  • 中文核心
  • EI
  • 中国科技核心
  • Scopus
  • CSCD
  • 英国科学文摘

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

分层向量自回归的多通道脑电信号的特征提取研究

王金甲 陈春

王金甲, 陈春. 分层向量自回归的多通道脑电信号的特征提取研究. 自动化学报, 2016, 42(8): 1215-1226. doi: 10.16383/j.aas.2016.c150461
引用本文: 王金甲, 陈春. 分层向量自回归的多通道脑电信号的特征提取研究. 自动化学报, 2016, 42(8): 1215-1226. doi: 10.16383/j.aas.2016.c150461
WANG Jin-Jia, CHEN Chun. Multi-channel EEG Feature Extraction Using Hierarchical Vector Autoregression. ACTA AUTOMATICA SINICA, 2016, 42(8): 1215-1226. doi: 10.16383/j.aas.2016.c150461
Citation: WANG Jin-Jia, CHEN Chun. Multi-channel EEG Feature Extraction Using Hierarchical Vector Autoregression. ACTA AUTOMATICA SINICA, 2016, 42(8): 1215-1226. doi: 10.16383/j.aas.2016.c150461

分层向量自回归的多通道脑电信号的特征提取研究


DOI: 10.16383/j.aas.2016.c150461
详细信息
    作者简介:

    陈春 燕山大学信息科学与工程学院硕士研究生.主要研究方向为信号处理.E-mail:xkz1124357698@sina.cn

    通讯作者: 王金甲 燕山大学信息科学与工程学院教授.主要研究方向为信号处理和模式识别.E-mail:wjj@ysu.edu.cn
  • 基金项目:

    河北省博士后专项基金 B2014010005

    中国博士后科学基金 2014M561202

    国家自然科学基金 61473339

    首批“河北省青年拔尖人才”项目 [2013]17

Multi-channel EEG Feature Extraction Using Hierarchical Vector Autoregression

More Information
    Author Bio:

    Master student at the School of Information Science and Engineering, Yanshan University. Her main research interest is signal processing.E-mail:

    Corresponding author: WANG Jin-Jia Professor at the School of Information Science and Engineering, Yanshan University. His research interest covers signal processing and pattern recognition. Corresponding author of this paper.E-mail:wjj@ysu.edu.cn
  • Fund Project:

    Hebei Province Postdoctoral Special Foundation B2014010005

    China Postdoctoral Science Foundation 2014M561202

    National Natural Science Foundation of China 61473339

    Hebei Province Top Young Talents [2013]17

  • 摘要: 有效的特征提取方法能提高脑机接口(Brain-computer interface,BCI)系统对脑电(Electroencephalogram,EEG)信号的识别率.因脑电信号都是多通道的,本文将分层向量自回归(Hierarchical vector autoregression,HVAR)模型用于脑电信号的特征提取,并结合传统的线性支持向量机(Support vector machine,SVM)用于脑电信号识别.该模型不仅克服了自回归(Autoregression,AR)模型只能用来提取单通道特征的局限性,而且不再采用传统VAR(Vector autoregression)模型所有通道共用一个时滞的处理方法.创新之处在于在传统的VAR模型基础上添加正则化思想,有效地压缩参数空间,实现合理的分层结构.本文首次将HVAR模型用于由Keirn等采集并整理的脑电数据中.实验结果证明HVAR模型在阶数较小的情况下(2阶)与阶数较大(6阶)的AR模型效果相当,可见低阶的HVAR能很好地刻画脑电信号的时空关联关系,这说明HVAR可能是刻画EEG信号的一种新颖的方法,这对其他多通道时间序列分析都有借鉴意义.
  • 图  1  ${\rm{LASSO}}-{\rm{VA}}{{\rm{R}}_2}(4)$ 得到的稀疏模式示例图

    Fig.  1  The sparse pattern for ${\rm{LASSO}}-{\rm{VA}}{{\rm{R}}_2}(4)$

    图  2  ${\rm{HVAR}}{{\rm{C}}_{\rm{3}}}(4)$ 的分层时滞结构示例图

    Fig.  2  A componentwise hierarchical lag structure: ${\rm{HVAR}}{{\rm{C}}_{\rm{3}}}(4)$

    图  3  ${\rm{HVAR}}{{\rm{O}}_{\rm{3}}}(4)$ 的分层时滞结构示例图

    Fig.  3  An own-other hierarchical lag structure: ${\rm{HVAR}}{{\rm{O}}_{\rm{3}}}(4)$

    图  4  ${\rm{HVAR}}{{\rm{E}}_{\rm{3}}}(4)$ 的分层时滞结构示例图

    Fig.  4  An elementwise hierarchical lag structure: ${\rm{HVAR}}{{\rm{E}}_{\rm{3}}}(4)$

    图  5  采样电极放置图

    Fig.  5  Electrode placement

    图  6  受试者1的分类任务6个通道的多元时序图

    Fig.  6  The sequence diagram of each channel for Subject 1

    图  7  第5 种任务组合方式对应的不同方法的正确率箱线图

    Fig.  7  The boxplot of the fifth task combination fordifferent methods

    图  8  第3 种任务组合方式对应的不同方法的正确率箱线图

    Fig.  8  The boxplot of the third task combination fordifferent methods

    图  9  不同受试者的分类结果

    Fig.  9  The results for all subjects

    图  10  阶数为 6 时不同方法的分类结果箱线图

    Fig.  10  The boxplots of all methods using six order for allsubjects

    表  1  受试者的样本数

    Table  1  Trails of subjects

    受试者样本数
    110
    25
    310
    410
    515
    610
    75
    下载: 导出CSV

    表  2  任务组合方式

    Table  2  Patterns of task combinations

    任务编号组合方式任务编号组合方式
    1基准、乘法计算6乘法计算、几何图旋转
    2基准、字母组合7乘法计算、视觉计算
    3基准、几何图旋转8字母组合、几何图旋转
    4基准、视觉计算9字母组合、视觉计算
    5乘法计算、字母组合10几何图旋转、视觉计算
    下载: 导出CSV

    表  3  VAR模型不同时滞的平均正确率

    Table  3  Average accuracy rate using VAR model with differentorder

    12345678910
    20.830.730.850.850.810.910.840.770.740.78
    30.790.740.810.820.800.880.800.730.730.75
    40.750.750.800.820.800.860.800.710.710.70
    50.730.710.780.780.810.860.770.700.650.70
    60.730.720.750.760.780.840.760.680.680.69
    70.710.670.710.730.770.810.720.650.650.65
    下载: 导出CSV

    表  4  LASSO-VAR模型不同时滞的平均正确率

    Table  4  Average accuracy rate using LASSO-VAR model with different order

    12345678910
    20.790.700.830.790.810.880.830.700.720.67
    30.800.730.790.800.800.860.780.650.650.65
    40.790.730.820.790.780.850.760.650.630.64
    50.770.690.790.760.760.830.710.670.610.64
    60.760.710.780.730.750.840.750.660.620.63
    70.760.660.760.730.730.840.750.630.580.61
    下载: 导出CSV

    表  5  HVARC模型不同时滞的平均正确率

    Table  5  Average accuracy rate using HVARC model with different order

    12345678910
    20.820.730.850.800.860.900.840.730.700.70
    30.810.750.820.840.840.900.840.700.680.69
    40.800.760.840.820.830.890.820.690.660.68
    50.780.730.810.790.850.870.780.690.640.64
    60.800.730.810.790.820.880.810.700.650.66
    70.780.730.810.790.800.870.810.670.630.66
    下载: 导出CSV

    表  6  HVARO模型不同时滞的平均正确率

    Table  6  Average accuracy rate using HVARO model with different order

    12345678910
    20.810.720.840.810.840.900.860.710.720.72
    30.820.760.830.820.850.910.830.690.690.69
    40.820.740.840.800.830.880.820.670.680.67
    50.810.710.820.810.830.870.820.670.640.65
    60.810.720.810.790.810.870.820.670.670.66
    70.790.710.810.800.810.870.820.650.620.64
    下载: 导出CSV

    表  7  HVARE模型不同时滞的平均正确率

    Table  7  Average accuracy rate using HVARE model with different order

    12345678910
    20.790.700.830.800.820.870.830.690.710.68
    30.770.720.800.820.810.870.790.640.680.69
    40.780.730.800.800.790.860.770.650.650.64
    50.780.710.780.790.800.850.750.660.640.64
    60.780.710.760.790.800.850.770.680.640.66
    70.780.680.760.800.760.860.770.640.640.65
    下载: 导出CSV

    表  8  不同特征提取方法的结果总结

    Table  8  Summary of classification results for all subjects

    AR-BGVARLASSO-VARHVARCHVAROHVARE
    平均值0.780.810.770.790.790.77
    最大值0.830.910.880.910.900.87
    最佳任务乘法计算、乘法计算、乘法计算、乘法计算、乘法计算、乘法计算、
    1组合方式几何图旋转几何图旋转几何图旋转几何图旋转几何图旋转几何图旋转
    平均值0.740.730.670.680.690.67
    最大值0.890.820.760.740.770.76
    最佳任务字母组合、字母组合、字母组合、字母组合、字母组合、字母组合、
    2组合方式视觉计算视觉计算视觉计算视觉计算视觉计算视觉计算
    平均值0.670.730.680.710.700.68
    最大值0.770.840.780.820.770.81
    最佳任务字母组合、几何图旋转、几何图旋转、几何图旋转、字母组合、几何图旋转、
    3组合方式几何图旋转视觉计算视觉计算视觉计算几何图旋转视觉计算
    平均值0.770.820.750.770.780.75
    最大值0.930.910.860.890.910.86
    最佳任务乘法计算、乘法计算、乘法计算、乘法计算、乘法计算、乘法计算、
    4组合方式视觉计算视觉计算视觉计算视觉计算视觉计算视觉计算
    下载: 导出CSV
  • [1] 王行愚, 金晶, 张宇, 王蓓. 脑控:基于脑——机接口的人机融合控制. 自动化学报, 2013, 39(3):208-221 doi:  10.1016/S1874-1029(13)60023-3

    Wang Xing-Yu, Jin Jing, Zhang Yu, Wang Bei. Brain control:human-computer integration control based on brain-computer interface. Acta Automatica Sinica, 2013, 39(3):208-221 doi:  10.1016/S1874-1029(13)60023-3
    [2] 伏云发, 王越超, 李洪谊, 徐保磊, 李永程. 直接脑控机器人接口技术. 自动化学报, 2012, 38(8):1229-1246 doi:  10.3724/SP.J.1004.2012.01229

    Fu Yun-Fa, Wang Yue-Chao, Li Hong-Yi, Xu Bao-Lei, Li Yong-Cheng. Direct brain-controlled robot interface technology. Acta Automatica Sinica, 2012, 38(8):1229-1246 doi:  10.3724/SP.J.1004.2012.01229
    [3] McFarland D J, Wolpaw J R. Brain-computer interfaces for communication and control. Communications of the ACM, 2011, 54(5):60-66 doi:  10.1145/1941487
    [4] Yang B H, Yan G Z, Yan R G, Wu T. Adaptive subject-based feature extraction in brain-computer interfaces using wavelet packet best basis decomposition. Medical Engineering & Physics, 2007, 29(1):48-53 http://cn.bing.com/academic/profile?id=2153925452&encoded=0&v=paper_preview&mkt=zh-cn
    [5] Shannon M, Zen H, Byrne W. Autoregressive models for statistical parametric speech synthesis. IEEE Transactions on Audio, Speech, and Language Processing, 2013, 21(3):587-597 doi:  10.1109/TASL.2012.2227740
    [6] 孙会文, 伏云发, 熊馨, 杨俊, 刘传伟, 余正涛. 基于HHT运动想象脑电模式识别研究. 自动化学报, 2015, 41(9):1686-1692 http://www.aas.net.cn/CN/abstract/abstract18742.shtml

    Sun Hui-Wen, Fu Yun-Fa, Xiong Xin, Yang Jun, Liu Chuan-Wei, Yu Zheng-Tao. Identification of EEG induced by motor imagery based on Hilbert-Huang transform. Acta Automatica Sinica, 2015, 41(9):1686-1692 http://www.aas.net.cn/CN/abstract/abstract18742.shtml
    [7] 伏云发, 徐保磊, 李永程, 李洪谊, 王越超, 余正涛. 基于运动相关皮层电位握力运动模式识别研究. 自动化学报, 2014, 40(6):1045-1057 http://www.aas.net.cn/CN/abstract/abstract18374.shtml

    Fu Yun-Fa, Xu Bao-Lei, Li Yong-Cheng, Li Hong-Yi, Wang Yue-Chao, Yu Zheng-Tao. Recognition of actual grip force movement modes based on movement-related cortical potentials. Acta Automatica Sinica, 2014, 40(6):1045-1057 http://www.aas.net.cn/CN/abstract/abstract18374.shtml
    [8] Mousavi E A, Maller J J, Fitzferald P B, Lithgow B J. Wavelet common spatial pattern in asynchronous offline brain computer interfaces. Biomedical Signal Processing and Control, 2011, 6(2):121-128 doi:  10.1016/j.bspc.2010.08.003
    [9] Li P Y, Wang X R, Li F L, Zhang R, Ma T, Peng Y H, Lei X, Tian Y, Guo D Q, Liu T J, Yao D Z, Xu P. Autoregressive model in the Lp norm space for EEG analysis. Journal of Neuroscience Methods, 2015, 240:170-174 doi:  10.1016/j.jneumeth.2014.11.007
    [10] Jain N, Dandapat S. Constrained autoregressive (CAR) model. In:Proceedings of 2005 Annual IEEE India International Conference. Chennai, India:IEEE, 2005.255-257
    [11] Huan N J, Palaniappan R. Neural network classification of autoregressive features from electroencephalogram signals for brain-computer interface design. Journal of Neural Engineering, 2004, 1(3):142-150 doi:  10.1088/1741-2560/1/3/003
    [12] Lawhern V, Hairston W D, McDowell K, Westerfield M, Robbins K. Detection and classification of subject-generated artifacts in EEG signals using autoregressive models. Journal of Neuroscience Methods, 2012, 208(2):181-189 doi:  10.1016/j.jneumeth.2012.05.017
    [13] Chen L L, Madhavan R, Rapoport B I, Anderson W S. Real-time brain oscillation detection and phase-locked stimulation using autoregressive spectral estimation and time-series forward prediction. IEEE Transactions on Biomedical Engineering, 2013, 60(3):753-762 doi:  10.1109/TBME.2011.2109715
    [14] Anderson C W, Stolz E A, Shamsunder S. Multivariate autoregressive models for classification of spontaneous electroencephalographic signals during mental tasks. IEEE Transactions on Biomedical Engineering, 1998, 45(3):277-286 doi:  10.1109/10.661153
    [15] Pei X M, Zheng C X. Feature extraction and classification of brain motor imagery task based on MVAR model. In:Proceedings of 2004 International Conference on Machine Learning and Cybernetics. Shanghai, China:IEEE, 2004.3726-3730
    [16] Hu X, Nenov V. Multivariate AR modeling of electromyography for the classification of upper arm movements. Clinical Neurophysiology, 2004, 115(6):1267-1287 http://cn.bing.com/academic/profile?id=1970693885&encoded=0&v=paper_preview&mkt=zh-cn
    [17] Wang J, Xu G Z, Wang L, Zhang H Y. Feature extraction of brain-computer interface based on improved multivariate adaptive autoregressive models. In:Proceedings of the 3rd International Conference on Biomedical Engineering and Informatics. Yantai, China:IEEE, 2010.895-898
    [18] Zhao C L, Zheng C X, Zhao M, Tu Y L, Liu J P. Multivariate autoregressive models and kernel learning algorithms for classifying driving mental fatigue based on electroencephalographic. Expert Systems with Applications, 2011, 38(3):1859-1865 doi:  10.1016/j.eswa.2010.07.115
    [19] Heger D, Terziyska T, Schultz T. Connectivity based feature-level filtering for single-trial EEG BCIS. In:Proceedings of the 2014 IEEE International Conference on Acoustics, Speech and Signal Processing. Florence, Italy:IEEE, 2014.2064-2068
    [20] Faes L, Erla S, Porta A, Nollo G. A framework for assessing frequency domain causality in physiological time series with instantaneous effects. Philosophical Transactions of the Royal Society A-Mathematical, Physical and Engineering Sciences, 2013, 371(1997):20110618 doi:  10.1098/rsta.2011.0618
    [21] Varotto G, Fazio P, Rossi Sebastiano D, Duran D, D'Incerti L, Parati E, Sattin D, Leonardi M, Franceschetti S, Panzica F. Altered resting state effective connectivity in long-standing vegetative state patients:an EEG study. Clinical Neurophysiology, 2014, 152(1):63-68 http://cn.bing.com/academic/profile?id=2032694363&encoded=0&v=paper_preview&mkt=zh-cn
    [22] Panzica F, Camafoglia L, Framceschetti S. EEG-EMG information flow in movement-activated myoclonus in patients with Unverricht-Lundborg disease. Clinical Neurophysiology, 2014, 125(9):1803-1808 doi:  10.1016/j.clinph.2014.01.005
    [23] Wang J J, Zhang Y N. A novel method of multi-channel feature extraction combining multivariate autoregression and multiple-linear principal component analysis. Journal of Biomedical Engineering, 2015, 32(1):19-24 http://cn.bing.com/academic/profile?id=2403326835&encoded=0&v=paper_preview&mkt=zh-cn
    [24] Tibshirani R. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, Series B, 1996, 58(1):267-288 http://cn.bing.com/academic/profile?id=2135046866&encoded=0&v=paper_preview&mkt=zh-cn
    [25] Yuan M, Lin Y. Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society, Series B, 2006, 68(1):49-67 doi:  10.1111/rssb.2006.68.issue-1
    [26] Zhao P, Rocha G, Yu B. The composite absolute penalties family for grouped and hierarchical variable selection. The Annals of Statistics, 2009, 37(6A):3468-3497 doi:  10.1214/07-AOS584
    [27] Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems.SIAM Journal on Image Sciences, 2009, 2(1):183——202 doi:  10.1137/080716542
    [28] Jenatton R, Mairal J, Obozinski G, Obozinski G, Bach F. Proximal methods for hierarchical sparse coding. The Journal of Machine Learning Research, 2011, 12:2297-2334 http://cn.bing.com/academic/profile?id=1539012881&encoded=0&v=paper_preview&mkt=zh-cn
    [29] Keirn Z A, Aunon J I. A new mode of communication between man and his surroundings. IEEE Transactions on Biomedical Engineering, 1990, 37(12):1209-1214 doi:  10.1109/10.64464
  • [1] 张建朋, 裴雨龙, 刘聪, 李邵梅, 陈鸿昶. 基于因子图模型的动态图半监督聚类算法[J]. 自动化学报, 2020, 46(4): 670-680. doi: 10.16383/j.aas.c170363
    [2] 孙小棋, 李昕, 蔡二娟, 康健楠. 改进模糊熵算法及其在孤独症儿童脑电分析中的应用[J]. 自动化学报, 2018, 44(9): 1672-1678. doi: 10.16383/j.aas.2018.c170334
    [3] 杨默涵, 陈万忠, 李明阳. 基于总体经验模态分解的多类特征的运动想象脑电识别方法研究[J]. 自动化学报, 2017, 43(5): 743-752. doi: 10.16383/j.aas.2017.c160175
    [4] 刘志勇, 孙金玮, 卜宪庚. 单通道脑电信号眼电伪迹去除算法研究[J]. 自动化学报, 2017, 43(10): 1726-1735. doi: 10.16383/j.aas.2017.c160191
    [5] 张毅, 尹春林, 蔡军, 罗久飞. Bagging RCSP脑电特征提取算法[J]. 自动化学报, 2017, 43(11): 2044-2050. doi: 10.16383/j.aas.2017.c160094
    [6] 耿志强, 张怡康. 一种基于胶质细胞链的改进深度信念网络模型[J]. 自动化学报, 2016, 42(6): 943-952. doi: 10.16383/j.aas.2016.c150727
    [7] 俞斌峰, 季海波. 稀疏贝叶斯混合专家模型及其在光谱数据标定中的应用[J]. 自动化学报, 2016, 42(4): 566-579. doi: 10.16383/j.aas.2016.c150255
    [8] 张靖, 周明全, 张雨禾, 耿国华. 基于马尔科夫随机场的散乱点云全局特征提取[J]. 自动化学报, 2016, 42(7): 1090-1099. doi: 10.16383/j.aas.2016.c150627
    [9] 唐朝辉, 朱清新, 洪朝群, 祝峰. 基于自编码器及超图学习的多标签特征提取[J]. 自动化学报, 2016, 42(7): 1014-1021. doi: 10.16383/j.aas.2016.c150736
    [10] 孟明, 朱俊青, 佘青山, 马玉良, 罗志增. 多类运动想象脑电信号的两级特征提取方法[J]. 自动化学报, 2016, 42(12): 1915-1922. doi: 10.16383/j.aas.2016.c160122
    [11] 倪鼎, 马洪兵. 基于近邻协同的高光谱图像谱-空联合分类[J]. 自动化学报, 2015, 41(2): 273-284. doi: 10.16383/j.aas.2015.c140043
    [12] 孙会文, 伏云发, 熊馨, 杨俊, 刘传伟, 余正涛. 基于HHT运动想象脑电模式识别研究[J]. 自动化学报, 2015, 41(9): 1686-1692. doi: 10.16383/j.aas.2015.c150007
    [13] 余航, 焦李成, 刘芳. 基于上下文分析的无监督分层迭代算法用于SAR图像分割[J]. 自动化学报, 2014, 40(1): 100-116. doi: 10.3724/SP.J.1004.2014.00100
    [14] 王行愚, 金晶, 张宇, 王蓓. 脑控: 基于脑——机接口的人机融合控制[J]. 自动化学报, 2013, 39(3): 208-221. doi: 10.3724/SP.J.1004.2013.00208
    [15] 林玉娥, 顾国昌, 刘海波, 沈晶, 赵靖. 适用于小样本问题的具有类内保持的正交特征提取算法[J]. 自动化学报, 2010, 36(5): 644-649. doi: 10.3724/SP.J.1004.2010.00644
    [16] 詹宇斌, 殷建平, 刘新旺. 基于大间距准则和图像矩阵双向投影的人脸特征提取方法[J]. 自动化学报, 2010, 36(12): 1645-1654. doi: 10.3724/SP.J.1004.2010.01645
    [17] 高全学, 谢德燕, 徐辉, 李远征, 高西全. 融合局部结构和差异信息的监督特征提取算法[J]. 自动化学报, 2010, 36(8): 1107-1114. doi: 10.3724/SP.J.1004.2010.01107
    [18] 徐科, 李文峰, 杨朝霖. 基于幅值谱与不变矩的特征提取方法及应用[J]. 自动化学报, 2006, 32(3): 470-474.
    [19] 杜恩祥, 李科杰. 基于多重分形和小波变换的声目标信号特征提取[J]. 自动化学报, 2004, 30(5): 742-746.
    [20] 谭枫, 曾小明. 基于类别可分离性的遥感图象特征提取方法[J]. 自动化学报, 1990, 16(2): 174-178.
  • 加载中
图(10) / 表(8)
计量
  • 文章访问数:  1289
  • HTML全文浏览量:  247
  • PDF下载量:  819
  • 被引次数: 0
出版历程
  • 收稿日期:  2015-07-20
  • 录用日期:  2016-02-18
  • 刊出日期:  2016-08-01

分层向量自回归的多通道脑电信号的特征提取研究

doi: 10.16383/j.aas.2016.c150461
    基金项目:

    河北省博士后专项基金 B2014010005

    中国博士后科学基金 2014M561202

    国家自然科学基金 61473339

    首批“河北省青年拔尖人才”项目 [2013]17

    作者简介:

    陈春 燕山大学信息科学与工程学院硕士研究生.主要研究方向为信号处理.E-mail:xkz1124357698@sina.cn

    通讯作者: 王金甲 燕山大学信息科学与工程学院教授.主要研究方向为信号处理和模式识别.E-mail:wjj@ysu.edu.cn

摘要: 有效的特征提取方法能提高脑机接口(Brain-computer interface,BCI)系统对脑电(Electroencephalogram,EEG)信号的识别率.因脑电信号都是多通道的,本文将分层向量自回归(Hierarchical vector autoregression,HVAR)模型用于脑电信号的特征提取,并结合传统的线性支持向量机(Support vector machine,SVM)用于脑电信号识别.该模型不仅克服了自回归(Autoregression,AR)模型只能用来提取单通道特征的局限性,而且不再采用传统VAR(Vector autoregression)模型所有通道共用一个时滞的处理方法.创新之处在于在传统的VAR模型基础上添加正则化思想,有效地压缩参数空间,实现合理的分层结构.本文首次将HVAR模型用于由Keirn等采集并整理的脑电数据中.实验结果证明HVAR模型在阶数较小的情况下(2阶)与阶数较大(6阶)的AR模型效果相当,可见低阶的HVAR能很好地刻画脑电信号的时空关联关系,这说明HVAR可能是刻画EEG信号的一种新颖的方法,这对其他多通道时间序列分析都有借鉴意义.

English Abstract

王金甲, 陈春. 分层向量自回归的多通道脑电信号的特征提取研究. 自动化学报, 2016, 42(8): 1215-1226. doi: 10.16383/j.aas.2016.c150461
引用本文: 王金甲, 陈春. 分层向量自回归的多通道脑电信号的特征提取研究. 自动化学报, 2016, 42(8): 1215-1226. doi: 10.16383/j.aas.2016.c150461
WANG Jin-Jia, CHEN Chun. Multi-channel EEG Feature Extraction Using Hierarchical Vector Autoregression. ACTA AUTOMATICA SINICA, 2016, 42(8): 1215-1226. doi: 10.16383/j.aas.2016.c150461
Citation: WANG Jin-Jia, CHEN Chun. Multi-channel EEG Feature Extraction Using Hierarchical Vector Autoregression. ACTA AUTOMATICA SINICA, 2016, 42(8): 1215-1226. doi: 10.16383/j.aas.2016.c150461
  • 通过建立计算机与脑电波之间的联系从而实现大脑与外界信息的交换已经成为多年的一个研究热点[1-2]. 脑机接口(Brain-computerinterface,BCI) 是一种新的人机交互方式,只需要通过大脑的活动信号将命令传给计算机或其他外部设备[3].这种不依赖外围神经肌肉的通信系统为完全瘫痪的人提供全新的生活方式.

    为了能够正确有效地控制BCI系统,必须通过一定的技术将不同的大脑活动转换为系统可以识别的模式从而发送不同的命令到外部设备中.BCI系统的核心之一是脑电图(Electroencephalogram,EEG)的信号分类,包括预处理、特征提取和选择、分类三个过程.由此可见,在EEG识别过程中,特征提取起着非常重要的作用,提取出来的特征会直接影响系统的分类正确率. 目前,脑信号的特征提取方法主要包括:小波变换和小波包变换[4]、自回归(Autoregression,AR)模型参数估计[5]、时频分析(Time-frequencyanalysis)[6-7]、共空间模式 (Common spatial patterns,CSP)[8]等. 大致可以将 上述方法分为参数法和非参数法两类. 非参数法需要大量的样本,脑电信号潜在信息的泄露不可避免; 而参数法是从有限的记录数据中得到估计参数来描述随机信号,基于该方法不需要对原始数据的产生做任何假设使得它在脑信号研究中得到广泛的应用[9]. 其中,AR模型因其表达式简单、 适用范围广等优点[10]受到众多学者青睐. Huan等使用3种不同的算法求解AR模型参数,将得到的系数作为神经网络分类器的输入实现EEG信号分类[11];Lawhern等利用AR模型对连续的脑电信号进行特征提取,与支持向量机(Support vector machine,SVM)结合进行分类,在几个不同的人工制造的脑电信号上取得了接近94%的正确率[12]; Chen等将AR模型应用到提取EEG的相位和频率信息中[13];Li等在AR模型的基础上添加了Lp范数的约束条件来压缩噪声等外界干扰对脑电信号的分析[9]. 但是,AR模型是针对单个通道进行分析的,所以在对多通道脑电信号的处理上会忽略各个通道之间的相互作用,不能有效地挖掘出通道之间的互信息. 相比之下,针对多通道信号建立的多变量自回归模型(Multivariateautoregressive,MVAR) 更适合用于脑信号的研究.Anderson等对4个受试者6个通道的两类任务进行四种方法的特征提取,结果表明使用MVAR系数表示不同的行为模式能取得更好的分类性能,且结果更为稳定[14]; Pei等在对左右手运动两类任务进行分类时首先用MVAR进行特征提取,然后使用基于Mahalanobis距离的判别分类器进行分类,实验结果证明,在分类正确率、 信噪比和互信息等方面,MVAR模型都优于AR模型[15]. Hu等将AR、MVAR模型运用于基于胳膊运动的分类任务上,得到了同样的结论,即MVAR模型因能够挖掘出肌电图通道之间的互信息可以提高分类正确率[16];Wang等提出了新的多元自适应自回归(Multivariate adaptiveautoregressive,MVAAR) 模型对左右手和脚三类任务进行特征提取,结果表明MVAAR模型对于在线BCI系统是一种更为有效的特征提取方法[17];Zhao等结合MVAR模型和核主成分分析(Kernel principal componentanalysis,KPCA) 降维方法用于EEG信号检测驾驶精神疲劳状态,不仅提高了正确率,而且加快了收敛速度[18]. Heger等用MVAR模型提取特征,分析脑信号信息流之间的联系[19]; Faes等提出了扩充MVAR (Extended MVAR,eMVAR) 框架,并将其成功地应用在心血管变异时间序列和多通道的脑电图中[20];Varotto等运用MVAR模型分析 使用频域线性索引连接的脑信号[21];Panzica等使用MVAR模型对皮质性肌阵挛患者的脑电图和肌电图进行分析[22];Wang等结合MVAR与MPCA (Mutilinear principal componentanalysis)降维方法 对UCI数据库中的脑电数据进行分类,分类器选用贝叶斯分类器,取得了很好的分类效果[23].但MVAR模型系数是矩阵形式,在进行分析时很容易出现高维小样本问题.

    本文将正则化思想添加在MVAR模型中,并提出分层结构,不仅能克服传统MVAR模型 共用一个时滞的缺点,而且由于正则化罚的使用有效地压缩了模型的参数空间,将MVAR模型成功推广到高维小样本数据的应用上,此外,采用加速近邻梯度法求解模型加快了算法的收敛速度,降低了时间复杂度.

    • 1980年,Sims首先提出向量自回归(Vector autoregression,VAR) 模型,即上面提到的MVAR模型.这种模型通过采用多方程联立的形式挖掘各个变量之间的线性关系.定义 $\left\{ {{{ y}_{t}} \in {{\bf R}^k}} \right\}_{{t =1}}^{T}$ 是长度为Tk维时间序列,p阶的VAR模型具有如下形式:

      $$\begin{align} & {{y}_{t}}=v+{{B}^{(1)}}{{y}_{t-1}}+\cdots +{{B}^{(p)}}{{y}_{t-p}}+{{u}_{t}}, \\ & t=1,\cdots ,T \\ \end{align}$$ (1)

      其中,p是预先给定的最大时滞, ${ v} \in {{\bf R}^{k}}$ 是k维截距向量, $\left\{ {{B^{(l)}} \in {{\bf R}^{k \times k}}} \right\}_{l = 1}^p$ 是时滞为l 的系数矩阵, $B_{ij}^l$ 表示的是第j 个时间序列前l 个时刻的值对第i个时间序列当前时刻值的影响. $\left\{ {{{ u}_{t}} \in {{\bf R}^k}} \right\}_{t = 1}^T$ 定义的是均值为0、同期协方差矩阵为 ${\Sigma _u}$ 的白噪声,即 ${{ u}_t}$ 满足以下两个条件: $1){\rm E}({{ u}_t}) = 0;\;2){\rm E}({{ u}_t}{{ u}_\tau }) =\left\{\begin{array}{l}0,\;\;\;t e \tau \\{\Sigma_u},\;t = \tau\end{array} \right.$ .

      最小二乘估计(Ordinary least square estimation,OLSE)是求解线性回归模型 参数估计最常用的方法.使用最小二乘估计VAR模型的参数,则参数估计可以描述为:

      $$({ v},B){\rm{ = }}\mathop {\arg} \min \limits_{({ v},B)}\frac{1}{2}\sum\limits_{t = 1}^{T} {\left\| {{{ y}_t} - { v}- \sum\limits_{l = 1}^p {{B^{(l)}}{{ y}_{t - l}}} } \right\|}_2^2$$ (2)

      对于 ${\rm{VA}}{{\rm{R}}_k}(p)$ 模型,需要估计的参数个数为 ${k^2}p$ ,也就是说,参数空间随着时间序列的增加呈平方增长,这对低频高维的时序数据很容易出现参数任意化问题.所以要想将VAR模型应用在高维数据上,必须采取一定的正则化技术压缩参数空间.

      Tibshirani 提出的LASSO模型[24]是在最小二乘的基础上对系数的1-范数施加约束,不仅能够达到压缩系数的目的,而且能够产生稀疏解,从而降低参数空间.

      为表示方便,定义: $ k \times T$ 维的响应矩阵 $ Y = {y_1},\cdots,{y_T} $ ,表示所有时间序列不同时刻的值; $k \times kp $ 维的参数矩阵 $B =[{B^{(1)}} \cdots {B^{(p)}}]$ ,表示不同序列间的相互影响; $kp\times 1$ 维的向量 ${{ z}_t} = {\left[{{ y}_{t-1}^{\rm{T}}\cdots { y}_{t - p}^{\rm{T}}} \right]^{\rm{T}}}$ ,表示所有时间序列滞后当前时刻1,2, $\cdots$ ,p 个时刻的值; $kp \times T$ 维的协变量矩阵 $Z = \left[{{z_1} \cdots {z_{T}}}\right] $ ,此处,我们需要用到 $\left\{ {{{ y}_{ - (p - 1)}},\cdots ,{{ y}_0}} \right\} $ 来初始化 Z; $k \times T$ 维的扰动矩阵 $U = \left[{{{ u}_1}\cdots{{ u}_{{T}}}}\right] $ 以及 $T \times 1$ 维的全 1 向量 ${\bf 1} = {\left[{1\cdots1} \right]^{\rm{T}}}$ . 在最小二乘的基础上添加 ${L_1}$ 罚,也就是系数的绝对值之和. 这就是我们所说的 LASSO-VAR 模型:

      $$\min \frac{1}{2}\left\| {Y - { v}{{\bf 1}^{\rm T}} - BZ}\right\|_2^2 + \lambda {\left\| B \right\|_1}$$ (3)

      其中, $\left\| A \right\|_2^{} = \sqrt {\sum_{i = 1}^m {\sum_{j =1}^n {{a_{ij}}^2} } }$ 是矩阵 ${A_{m \times n}}$ 的F-范数, ${\left\| A \right\|_1} = \sum_i {\sum_j {\left| {{a_{ij}}} \right|}} $ 是矩阵 ${A_{m \times n}}$ 的1-范数, $\lambda $ 是正则化参数.LASSO-VAR与Bayesian相比,不仅可以使最小二乘解的某些项收缩至0,实现变量选择,而且使得在回归参数个数 $k \times kp$ 接近甚至超过样本大小 $k \times T$ 的情况下也能得到可靠估计.下一节我们将采用坐标下降法来求解该模型.

    • 记 $f(B,{ v}) = \frac{1}{2}\left\| {Y - { v}{{\bf 1}^{\rm T}} - BZ} \right\|_2^2 + \lambda {\left\| B \right\|_1}$ ,将每一个时间序列的回归模型看作一个子问题,那么可将目标函数分解为k个子问题. 第j个子问题的形式为:

      $${f_j}({B_{j.}},{v_j}) = \frac{1}{2}\left\| {{Y_{jt}} - {v_j}{{\bf 1}^{\rm T}} - {B_{j.}}Z} \right\|_2^2 + \lambda {\left\|{{B_{j.}}} \right\|_1}$$ (4)

      式中, ${v_j}$ 和 ${B_{j.}}$ 为未知变量.首先通过求 ${f_j}({B_{j.}},{v_j})$ 关于 ${v_j}$ 的偏导数并令其为0得到截距值 ${v_j}$ . 即: $0 = { abla_v}{f_j}({B_{j.}},{v_j}) = ({Y_{jt}} - {v_j}{{\bf 1}^{\rm T}} -{B_{j.}}Z){\bf 1}$ ,从而求得: ${v_j} = {\overline Y _j} -{\widehat B_{j.}}\overline { Z} $ . 所以截距向量 ${ v} =\overline { Y} - \widehat B\overline { Z}$ . 其中, $\overline { Y} $ 是 $k \times 1$ 维的响应矩阵的行均值向量, $\overline { Z} $ 是 $kp \times 1$ 维的协变量矩阵的行均值向量.接下来用坐标下降法求 ${B_{j.}}$ .坐标下降法的基本思想是每次只优化一维变量,其余变量看作常量,这样优化问题可以很快完成,且相关方程可以在变量循环中更新;而通常情况下,极小化的变量在众多参数中不随循环改变,因此整个迭代过程将很快完成.

      假设YZ都是中心化的, 则目标函数中不再包含截距向量.使用坐标下降法逐个更新系数 矩阵B的元素.假设当前估计是Bjr列的元素Bjr, 那么将式(4)转换为单变量形式:

      \begin{align} f_{jr}({B_{jr}}) = & \dfrac{1}{2}\sum\limits_t {{{\left( {{Y_{jt}} - \sum\limits_{l e r} {{B_{jl}}{Z_{lt}} - {B_{jr}}{Z_{rt}}} } \right)}^2}} + \\ & \lambda \sum\limits_{l e r} {\left| {{B_{jl}}} \right|} + \lambda \left| {{B_{jr}}} \right| \end{align} (5)

      分 ${B_{jr}} > 0, {B_{jr}} < 0, {B_{jr}} = 0$ 三种情况对Bjr 求导,令 ${R_t} = {Y_{jt}} - \sum_{l e r}{{B_{jl}}{Z_{lt}}_{}} $ , 结果如下:

      \begin{align} {B_{jr}} =~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\\ \left\{ \begin{array}{l} \dfrac{{\sum\limits_t {{R_t}{Z_{rt}} - \lambda } }}{{\sum\limits_t {{Z_{rt}}^2} }},\quad \sum\limits_t {{R_t}{Z_{rt}} > 0,\;\left| {\sum\limits_t {{R_t}{Z_{rt}}} } \right| > \lambda } \\[8mm] \dfrac{{\sum\limits_t {{R_t}{Z_{rt}} + \lambda } }}{{\sum\limits_t {{Z_{rt}}^2} }},\quad \sum\limits_t {{R_t}{Z_{rt}} < 0,\;\left| {\sum\limits_t {{R_t}{Z_{rt}}} } \right| > \lambda } \\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left| {\sum\limits_t {{R_t}{Z_{rt}}} } \right|< \lambda \end{array} \right. \end{align} (6)

      可见Bjr的坐标更新为软阈值形式: ${B_{jr}} = \frac{{S(\sum_t {{R_t}{Z_{rt}},\lambda )} }}{{\sum_t {{Z_{rt}}^2} }}$ . 其中, $S(x,\phi ) = {\mathop{\rm sgn}} (x){(\left| x \right| - \phi )_ +}, {(\left| x \right| - \phi )_ + }=\max (\left| x \right| - \phi ,0)$ .

      从系数矩阵的求解过程中可以看出LASSO-VAR模型是将所有变量看作无关变量,这也使得虽然通过添加L1 罚能得到稀疏解,但是这种稀疏模式没有固定结构. 举例来说: 最大时滞为4的二元时间序列,通过LASSO-VAR2(4) 得到的系数的可能稀疏模式如图 1所示.

      图  1  ${\rm{LASSO}}-{\rm{VA}}{{\rm{R}}_2}(4)$ 得到的稀疏模式示例图

      Figure 1.  The sparse pattern for ${\rm{LASSO}}-{\rm{VA}}{{\rm{R}}_2}(4)$

    • HVAR (Hierarchical vector autoregression) 模型在LASSO-VAR的基础上,通过对系数矩阵B 的不同约束得到不同的分层结构,可以生成一个具有较高预测性能的可解释模型. 在讨论之前做以下定义:

      $$\begin{array}{l}{B^{(l:p)}} = [{B^{(l)}} \cdots {B^{(p)}}] \in {{\bf R}^{k \times k(p - l + 1)}}\\{B_i}^{(l:p)} = [{B_i}^{(l)} \cdots {B_i}^{(p)}] \in {{\bf R}^{1 \times k(p - l + 1)}}\\{B_{ij}}^{(l:p)} = [{B_{ij}}^{(l)} \cdots {B_{ij}}^{(p)}] \in{{\bf R}^{1 \times (p - l + 1)}}\end{array}$$

      此外,我们还定义了一个 $k \times k$ 维的最大时滞矩阵L,L中的每个元素定义如下: ${L_{ij}} = \left\{\begin{array}{l}\max \left\{ l \right\},\;\;\;B_{ij}^l e 0\\{\rm{0}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;B_{ij}^l{\rm{ = }}0\end{array} \right.,\;l = 1,\cdots,p$ . 因此, ${L_{ij}}$ 表示的是时间序列j对时间序列i的最长作用时间. 若L中的所有元素都相等即为普通的 ${\rm{VA}}{{\rm{R}}_k}(p)$ 模型,本文详细介绍三种分层结构.

      1) 分量方式HVARC: 将每一个时间序列 ${y_{jt}}$ 看成一个边际方程,每个方程有自己的最大时滞,即: ${L_{ij}} = {L_i},\;\forall j,i = 1,\cdots,k$ . 若定义时滞矩阵: $L = \left[\begin{array}{l}4\;4\;4\\{\rm{2}}\;{\rm{2}}\;{\rm{2}}\\{\rm{3}}\;{\rm{3}}\;{\rm{3}}\end{array} \right]$ ,则该分层结构的时滞矩阵如图 2所示. 图 2表示所有时间序列对第一个时间序列的作用时间都为4个单位,对第二个时间序列的作用时间为2个单位.

      图  2  ${\rm{HVAR}}{{\rm{C}}_{\rm{3}}}(4)$ 的分层时滞结构示例图

      Figure 2.  A componentwise hierarchical lag structure: ${\rm{HVAR}}{{\rm{C}}_{\rm{3}}}(4)$

      2) 自身--其他分层结构HVARO: 在HVARC的基础上增加假设条件:自身作用时间不短于其他时间序列的作用时间,且最多长于其他序列1个时间单位. 即: ${L_{ij}} = \left\{ \begin{array}{l}{L_i}^{other},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i e j\\{L_i}^{other}\;\mbox{或者}\;{L_i}^{other} + 1,\; i = j\end{array} \right.,i,j = 1,\cdots,k$. 表现在系数矩阵上是指: ${B}_{ij}^{\left( {\left[{{L_i}^{other} + 1} \right]:p} \right)}=0,{B}_{ij}^{\left( {{L_i}^{other}} \right)} e 0,\;i e j,\;{B}_{ii}^{\left( {{L_i}^{other}} \right)} e 0$ ,且对 ${{B}_{ii}}^{(other + 1)}$ 没有约束, ${{B}_{ii}}^{([other +2]:p)} = 0$ . 若定义时滞矩阵: ${L} = \left[\begin{array}{l}4\;{\rm{3}}\;{\rm{3}}\\{\rm{2}}\;{\rm{2}}\;{\rm{2}}\\{\rm{2}}\;{\rm{2}}\;{\rm{3}}\end{array} \right]$ ,则该分层结构的时滞矩阵如图 3所示.即: 第一个时间序列对其自身的作用时间为4个单位,而其他时间序列对第一个时间序列的作用时间为3个单位.

      图  3  ${\rm{HVAR}}{{\rm{O}}_{\rm{3}}}(4)$ 的分层时滞结构示例图

      Figure 3.  An own-other hierarchical lag structure: ${\rm{HVAR}}{{\rm{O}}_{\rm{3}}}(4)$

      图  4  ${\rm{HVAR}}{{\rm{E}}_{\rm{3}}}(4)$ 的分层时滞结构示例图

      Figure 4.  An elementwise hierarchical lag structure: ${\rm{HVAR}}{{\rm{E}}_{\rm{3}}}(4)$

      3) 单个元素的分层结构HVARE: 该分层结构比较灵活,不同的时间序列彼此之间的作用时间不同. 即: ${B}_{ij}^{\left({{L_{ij}}} \right)} e 0,\;{B}_{ij}^{\left( {[{L_{ij}} + 1]:p}\right)} = 0,\;\;\forall i,j$ . 若定义时滞矩阵: ${L} = \left[\begin{array}{l}4\;{\rm{1}}\;{\rm{2}}\\{\rm{2}}\;{\rm{0}}\;{\rm{1}}\\{\rm{4}}\;{\rm{4}}\;{\rm{2}}\end{array} \right]$ ,则时滞矩阵如图4所示. 彼此之间没有特定的约束关系.

      为了实现上述三种分层结构,我们采用分层组LASSO罚对系数加以约束.分层组LASSO罚是指具有嵌套结构的组LASSO罚.理论已经证明组LASSO罚[25]将每一组的变量看成一块,对块变量施加约束,能够实现变量的组选择;而嵌套组[26]能够实现分层选择,即某一组变量没有被选中,则嵌套于该组的另一组变量也不被选中.可以通过在VAR模型上添加不同的罚函数实现以上三种分层结构.

      1) HVARC: 前面已经提到该分层结构是指每一个边际方程的阶数相同,也就是要求满足这样的条件: 如果 $ {\widehat {B}_i}^{(l)} = {0}$ ,那么 ${\widehat B_i}^{(l')} = {0},l' > l,i = 1,\cdots,k$ .即嵌套组LASSO罚可以实现这样的分层结构,故只需在VAR模型的基础上添加嵌套组LASSO罚函数,目标函数定义如下:

      $${\min _{{ v},B}}\left\{ {\frac{1}{2}\left\| {Y - { v}{{\bf 1}^{\rm T}} - BZ} \right\|_2^2 + \lambda \sum\limits_{i = 1}^k{\sum\limits_{l = 1}^p {{{\left\| {B_i^{(l:p)}} \right\|}_2}} } }\right\}$$ (7)

      随着 $\lambda $ 的增加,i 越大,l越小,就会有 ${\widehat {B}_i}^{(l:p)} = {0}$ .这样的罚更趋向于选择一个较小的最大时滞,而不是简单地得到一个没有任何特殊结构的系数矩阵 B.

      2) HVARO: 该分层结构要求实现的是对角元素优于非对角元素,所以只需在HVARC 的基础上对非对角元素进行再一次的惩罚,故目标函数定义如下:

      $$\begin{align} & {{\min }_{v,B}}\{\frac{1}{2}\left\| Y-v{{\mathbf{1}}^{\text{T}}}-BZ \right\|_{2}^{2}+ \\ & \lambda \sum\limits_{i=1}^{k}{\sum\limits_{l=1}^{p}{\left[ {{\left\| B_{i}^{(l:p)} \right\|}_{2}}+{{\left\| B_{i,-i}^{(l)},B_{i}^{(\left[ l+1 \right]:p)} \right\|}_{2}} \right]}}\} \\ \end{align}$$ (8)

      其中, ${B}_{i,- i}^{\left( l \right)} = \left\{ {B_{ij}^{\left( l\right)}:j e i} \right\}$ ,并且约定 ${B}_i^{\left( {[p + 1]:p}\right)} = 0$ .式(8)中罚函数与HVARC的罚函数的不同之处在于增加了第二项,放宽了对对角元素 ${B}_{ii}^{\left( l \right)}$ 的约束,这样就能保证在时滞为l 时,即使序列之间的交互影响已经为0,而自身影响可能并不为0.该模型能够确保自身的影响可能长于其他时间序列的影响1个时间单位.

      3) HVARE: 该分层结构对时滞矩阵的结构没有硬性的约束,任意时间序列对的时滞都可以不同,所以只需对任意时间序列对的时滞矩阵进行约束,故目标函数定义如下:

      $${\min _{{ v},B}}\Bigg\{ \frac{1}{2}\left\| {Y - { v}{{\bf 1}^{\rm T}} - BZ} \right\|_2^2 + \lambda \times \\\sum\limits_{i = 1}^k {\sum\limits_{j = 1}^k {\sum\limits_{l =1}^p {\left[{{{\left\| {B_{ij}^{(l:p)}} \right\|}_2}} \right]} }} \Bigg\}$$ (9)

      在这样的分层结构中,因为每对时间序列都有自己的最大时滞,因此 $\widehat {B}_{ij}^{(l:p)} = {0}$ 可能出现在时序对 $\left({i,j} \right)$ 的任意 l值的情况下. 这样的分层结构是最灵活的.

      前面已经讨论过截距向量与响应矩阵、协变量矩阵的关系.此外仍假定响应矩阵与协变量矩阵都是中心化的,所以在下面的具体求解过程中将不再讨论截距向量. 式 (7) $\sim$ (9)都是基于分层组LASSO 罚的,在这儿引入函数 ${\Omega_i}(B_i^{(l:p)})$ 表示系数矩阵不同形式的范数,即将式(7) $\sim$ (9) 的罚函数统一用该函数表示,那么可将三个凸优化问题统一为以下形式:

      $${\min _B}\left\{ {\frac{1}{2}\left\| {Y - BZ} \right\|_2^2 +\lambda \sum\limits_{i = 1}^k {\sum\limits_{l = 1}^p {{\Omega_i}(B_i^{(l:p)})} } } \right\}$$ (10)

      因为每一个边际方程都是独立的,所以可以将式(10)分解成k个如下的子问题求解:

      $${\min _{{B_i}}}\left\{ {\frac{1}{2}\left\| {{Y_i} - {B_i}Z}\right\|_2^2 + \lambda \sum\limits_{l = 1}^p {{\Omega_i}(B_i^{(l:p)})} } \right\}$$ (11)
    • 下面详细讲述加速近邻梯度法求解式(11).加速近邻梯度法是由YurriNesterov 于1983年首先提出的[27].该方法以负梯度的方向作为下降方向,将上两代的迭代值做线性组合得到一个中间变量,然后由该中间变量导出当代的最优解.加速近邻梯度法通过对近邻梯度法的一个简单修正,不仅保持了梯度法的简洁性,更提高了算法的效率.理论上已经证明了加速近邻梯度法的收敛率是 $1/{k^2}$ .本文选取的相邻两次迭代值线性组合的权重为 ${\theta _m} =2/{(m +2)}$ ,步长 $t \leftarrow {\sigma _1}{\left( Z \right)^{ -2}}$ ,其中 ${\sigma _1}\left( Z \right)$ 是协变量矩阵Z的最大奇异值. 记 ${L_i}({B_i}) =\frac{1}{2}\left\| {{Y_i} - {B_i}Z} \right\|_2^2$ ,则有 $ abla{L_i}({B_i}) = \left( {{B_i}Z - {Y_i}} \right){Z^{\rm T}}$ ,式(11) 的系数矩阵B的详细更新步骤见算法1.

      算法 1. 加速近邻梯度法求解式(11)

      输入 . Y,Z, $\widehat B\left[0 \right],\lambda,\varepsilon$ ;

      步骤 1. 初始化: $\widehat B\left[1 \right]\leftarrow \widehat B\left[0 \right];\widehat B\left[2 \right]\leftarrow \widehat B\left[0 \right];t \leftarrow {\sigma_1}{\left( Z \right)^{ - 2}},i=1$ ;

      步骤 2. 迭代直至收敛:

      1) 更新中间变量 $\widehat b,\widehat b \leftarrow {\widehat B_i}\left[{m - 1} \right] + \frac{{m - 2}}{{m + 1}}\left({{{\widehat B}_i}\left[{m - 1} \right] - {{\widehat B}_i}\left[{m - 2} \right]} \right)$ ;

      2) 更新当代的系数值 $\;{\widehat B_i}\left[m \right],\;{\widehat B_i}\left[m \right] \leftarrow pro{x_{t\lambda\Omega _i^{}}}\left( {\widehat b - t\left( {\widehat bZ - {Y_i}}\right){Z^{\rm T}}} \right)$ ;

      3) 判断 $\left\| {\widehat b - {{\widehat B}_i}\left[m \right]}\right\| \le \varepsilon $ 是否成立,若成立,内层循环结束,转至步骤3; 否则, $m\leftarrow m+ 1$ ,转至1);

      步骤 3.i=k,算法结束;否则 $i \leftarrow i + 1$ ,转至步骤2;

      输出. 系数矩阵B.

      三种不同分层结构的罚函数都是基于分层组LASSO 罚,由近邻算子的定义可知:

      $$pro{x_{t\lambda {\Omega _i}}}\left( u \right) = \arg \min \left\{{\frac{1}{2}\left\| {u{\rm{ - }}v} \right\|_2^2 + t\lambda {\Omega_i}(v)} \right\}$$ (12)

      将三种范数的近邻算子统一记为:

      $$\widehat x = \arg \min \left\{ {\frac{1}{2}\left\| {x{\rm{ -}}\tilde x} \right\|_2^2 + \lambda \sum\limits_{h = 1}^H{{\omega _h}{{\left\| {{x_{{g_h}}}} \right\|}_2}} } \right\}$$ (13)

      其中, ${g_h}$ 是第h组包含的变量的索引集合,且 ${g_1} \subset\cdots \subset {g_H}$ . 我们已经知道对于组LASSO罚,组间是不重叠的,它的近邻问题按组可分离,近邻算子是广义软阈值算子. 而在分层结构中,嵌套组的出现导致不能直接得到分层组 LASSO罚的近邻算子. Jenatton等已经证明对于分层组LASSO 的近邻问题,可以通过求其对偶范数的近邻问题得到[28].本文直接使用Jenatton 等的结果. 见算法2.

      算法 2. 求解分层组 LASSO 罚的近邻算子

      输入. $\tilde x,\lambda ,{\omega _1},\cdots,{\omega _H}$ ;

      步骤 1. 初始化: $r \leftarrow \tilde x$ ;

      步骤 2. 循环 $h = 1,\cdots,H$:\\ $\quad\quad {r_{{g_h}}} \leftarrow {\left( {1 - \lambda {\omega_h}/\left\| {{r_{{g_h}}}} \right\|} \right)_ + }{r_{{g_h}}}$ ;

      输出. $\widehat x=r$ .

    • 将HVAR 模型应用到脑机接口的背景中. 实验采用的是由Keirn等采集的脑电数据,详见网址:http://www.cs.colostate.edu/eeg/main/data/1989\_Keirn\_and\_Aunon.该数据集记录的是7位受试者在执行想象任务时的EEG数据.每次实验要求受试者完成一种任务的想象,时间持续10s,采样频率是250Hz,所以每一个样本的每个通道会产生2500个采样点.在一个周期中每个动作重复执行5次,共5类任务.采样期间受试者的EEG数据通过7个位于头部的电极记录,采样电极的分布如图 5所示. 每个受试者的样本数列于表 1.

      图  5  采样电极放置图

      Figure 5.  Electrode placement

      表 1  受试者的样本数

      Table 1.  Trails of subjects

      受试者样本数
      110
      25
      310
      410
      515
      610
      75

      图 5中,电极C3,C4,P3,P4,O1,O2,EOG记录7个通道的EEG信号,电极A1,A2之间通过一个放大器(Grass7P511) 连接,该放大器的模拟滤波器的带宽为0.1Hz $\sim$ 100Hz,这两个电极的记录是电乳突的参考电位.

      表 1可知,受试者1,3,4,6 分别执行了两个周期的任务,受试者2,7 是执行了 一个周期的任务,受试者5则是执行了3个周期的任务.前面已提到每个受试者在每个周期中共执行5类任务,详细介绍如下:

      1) 基准任务: 在该任务下,要求受试者尽可能地放松,尽量不要去想任何事情.所以可将该任务看作脑电图的基准.

      2) 乘法计算任务: 要求受试者计算一个复杂的乘法问题,在计算过程中,受试者不能发声或有明显的动作,每次的乘法问题都是不重复的,因此受试者不可能直接给出答案,而且也已经证实受试者在每次的执行过程中都没有在10s之内给出答案,所以在采样期间受试者一直处于计算状态.

      3) 字母组合任务: 要求受试者不发声地给朋友或亲人写一封信.因为任务重复5次,所以在后边几次的任务中要求受试者回想在先前的信中是否有遗漏的字母.

      4) 几何图旋转任务: 在该任务下,首先给受试者30s的时间学习画一幅复杂的3维图,30s之后,画板撤掉,然后要求受试者想象该物体沿着某个轴旋转.

      5) 视觉计算任务: 要求受试者想象将数字按顺序写在一块黑板上.

    • 电极EOG是用来检测受试者是否有眨眼行为,已经有相关实验证明该通道对任务分类正确率没有影响[29],因此在实验中我们只用了6 个通道的值. 图 6 是受试者1 的5 类任务的6个通道的时序图.

      图  6  受试者1的分类任务6个通道的多元时序图

      Figure 6.  The sequence diagram of each channel for Subject 1

      图 6 的每一幅子图都是长度为2500 的6 元时间序列,由图可见,这6元时间序列基本上都是平稳的. 实验中用到的是各个任务的10 个样本,从表 1 可知,受试者2 和7 每一类任务只有5 个样本,而受试者4的数据中包含了未知的无效值,故可用的是受试者1,3,5,6的数据.对于受试者5,我们用的是他的前10 个样本,重新对1,3,5,6 编号为1,2,3,4. 此外,我们用到了0.5s的时间窗,所以将原始长度为 2500 的时间序列分割成20 段,每个时间段中包含了125 个采样点. 考虑到每个人的思维方式不同,所以对每个受试者单独进行实验.

      首先对这20段分别用VAR、LASSO-VAR、HVARC、HVARO、HVARE5种方法进行特征提取,那么每类任务的每一个样本就会得到20种特征模式,因此每类想象任务共得到200种特征

      模式,也就是说每类任务都会有200 个样本.然后将不同的任务进行两两组合,不同的任务组合方式见表 2;接着将组合后的 400个样本随机分成大小相等的训练集和测试集,共进行10次试验; 最后将10 个测试集上得到的正确率求平均.分类器我们选择线性支持向量机.实验中每种方法的阶数都选择2 $\sim$ 7.

      表 2  任务组合方式

      Table 2.  Patterns of task combinations

      任务编号组合方式任务编号组合方式
      1基准、乘法计算6乘法计算、几何图旋转
      2基准、字母组合7乘法计算、视觉计算
      3基准、几何图旋转8字母组合、几何图旋转
      4基准、视觉计算9字母组合、视觉计算
      5乘法计算、字母组合10几何图旋转、视觉计算
    • 实验中使用网格法选择正则化参数 $\lambda $ . 网格的深度与长度分别设置为30 和10. 首先我们详细分析受试者1 的实验结果. 受试者1 的不同任务组合方式,在5种不同特征提取方法的不同时滞下得到的10次试验的平均分类正确率如表 3~7所示.

      表 3~7可见,不同的任务组合方式所对应的最大时滞是不同的,但总的看来,当阶数等于2 时,不同的特征提取方法至少在5种任务组合方式下能获得最高正确率.由此我们可认为受试者1的最大时滞是2. 为了对10次实验的分类结果有更直观的认识,本文以第5种和第3种任务组合方式为例进行详细分析. 图 7图 8表示在阶数为2 时,五种不同方法分别在这两种情况下得到的分类结果的箱线图.

      表 3  VAR模型不同时滞的平均正确率

      Table 3.  Average accuracy rate using VAR model with differentorder

      12345678910
      20.830.730.850.850.810.910.840.770.740.78
      30.790.740.810.820.800.880.800.730.730.75
      40.750.750.800.820.800.860.800.710.710.70
      50.730.710.780.780.810.860.770.700.650.70
      60.730.720.750.760.780.840.760.680.680.69
      70.710.670.710.730.770.810.720.650.650.65

      表 4  LASSO-VAR模型不同时滞的平均正确率

      Table 4.  Average accuracy rate using LASSO-VAR model with different order

      12345678910
      20.790.700.830.790.810.880.830.700.720.67
      30.800.730.790.800.800.860.780.650.650.65
      40.790.730.820.790.780.850.760.650.630.64
      50.770.690.790.760.760.830.710.670.610.64
      60.760.710.780.730.750.840.750.660.620.63
      70.760.660.760.730.730.840.750.630.580.61

      表 5  HVARC模型不同时滞的平均正确率

      Table 5.  Average accuracy rate using HVARC model with different order

      12345678910
      20.820.730.850.800.860.900.840.730.700.70
      30.810.750.820.840.840.900.840.700.680.69
      40.800.760.840.820.830.890.820.690.660.68
      50.780.730.810.790.850.870.780.690.640.64
      60.800.730.810.790.820.880.810.700.650.66
      70.780.730.810.790.800.870.810.670.630.66

      表 6  HVARO模型不同时滞的平均正确率

      Table 6.  Average accuracy rate using HVARO model with different order

      12345678910
      20.810.720.840.810.840.900.860.710.720.72
      30.820.760.830.820.850.910.830.690.690.69
      40.820.740.840.800.830.880.820.670.680.67
      50.810.710.820.810.830.870.820.670.640.65
      60.810.720.810.790.810.870.820.670.670.66
      70.790.710.810.800.810.870.820.650.620.64

      表 7  HVARE模型不同时滞的平均正确率

      Table 7.  Average accuracy rate using HVARE model with different order

      12345678910
      20.790.700.830.800.820.870.830.690.710.68
      30.770.720.800.820.810.870.790.640.680.69
      40.780.730.800.800.790.860.770.650.650.64
      50.780.710.780.790.800.850.750.660.640.64
      60.780.710.760.790.800.850.770.680.640.66
      70.780.680.760.800.760.860.770.640.640.65

      图  7  第5 种任务组合方式对应的不同方法的正确率箱线图

      Figure 7.  The boxplot of the fifth task combination fordifferent methods

      图  8  第3 种任务组合方式对应的不同方法的正确率箱线图

      Figure 8.  The boxplot of the third task combination fordifferent methods

      从上边两个箱线图来看,不同任务组合对应的最佳特征提取方法是不一样的. 对于第5种任务组 合,较好的特征提取方法是HVAR方法,以HVARC为最佳,在10次实验中,无论是最小值、最大值、上下四分位点还是中位数都远远高于其他四种方法,但是HVARC方法得到的正确率的分布 不如VAR均匀. 而对于第3种任务组合,则是VAR的效果更好些,HVARC次之.两种方法50%的值的分布范围大致相同,但总的来看仍属VAR的分类正确率较高. 在两种情况下,没有固定稀疏结构LASSO-VAR的分类性能都不好,而HVARO和HVARE都在不同程度上出现了离散点.

      对其他受试者进行同样分析,得到受试者2、3、4的最大时滞也都为2,将各个受试者在不同方法下得到的分类结果绘于图 9,每个受试者都选用各自对应的阶数.图中横坐标代表 10种不同任务组合方式,纵坐标代表不同方法10次实验结果均值.

      图  9  不同受试者的分类结果

      Figure 9.  The results for all subjects

      图 9我们可得到以下结论: 1) 因每个人的思维方式不同,所以不同的受试者对应的最佳的特征提取方法不同; 2) 对于不同的受试者,最优的任务组合方式不同; 3)在传统的VAR模型的基础上添加合适的分层结构在某些识别任务中能提高分类性能,像受试者1的第5种任务组合方式以及受试者3的第4种任务组合方式.结合4位受试者来看,在阶数为2时,具有分层结构的HVAR模型并没有表现出明显的优势,甚至在绝大多数情况下得到的分类结果 都不如普通的VAR模型高. 此外,将EEG信号的不同通道不同时刻的值看作单一变量的LASSO-VAR方法,虽然能有效地压缩参数空间,但是作为EEG信号的特征提取方法不能得到较为满意的识别效果.

      图 9不能很好地看出各种方法的差异性,我们进一步采用重复测量方差分析的方法来探究不同特征提取方法的显著性,并重点分析受试者1的实验结果.首先对不同任务组合方式的分类结果进行方差齐次性检验,得到的p值分别为:0.1696、0.002571、0.08715、0.3105、0.2482、0.5554、0.6208、0.2263、0.5239、0.913.由此可知第2种任务组合方式的结果不满足方差齐性的要求,所以不再对其进行重复性方差分析,对其余9种任务组合方式进行分析得到的p值为:1.5E-08、6.59E-08、8.87E-05、2.18E-09、1.52E-05、0.00013、0.0307、0.000102、0.00258.可见差异都较为显著.

      表 8是本文用到的五种特征提取方法与文献[11]中方法的分类正确率的比较.文献[11]的分类器选择神经网络. 表中AR-BG表示用Burg算法求解AR系数.

      表 8  不同特征提取方法的结果总结

      Table 8.  Summary of classification results for all subjects

      AR-BGVARLASSO-VARHVARCHVAROHVARE
      平均值0.780.810.770.790.790.77
      最大值0.830.910.880.910.900.87
      最佳任务乘法计算、乘法计算、乘法计算、乘法计算、乘法计算、乘法计算、
      1组合方式几何图旋转几何图旋转几何图旋转几何图旋转几何图旋转几何图旋转
      平均值0.740.730.670.680.690.67
      最大值0.890.820.760.740.770.76
      最佳任务字母组合、字母组合、字母组合、字母组合、字母组合、字母组合、
      2组合方式视觉计算视觉计算视觉计算视觉计算视觉计算视觉计算
      平均值0.670.730.680.710.700.68
      最大值0.770.840.780.820.770.81
      最佳任务字母组合、几何图旋转、几何图旋转、几何图旋转、字母组合、几何图旋转、
      3组合方式几何图旋转视觉计算视觉计算视觉计算几何图旋转视觉计算
      平均值0.770.820.750.770.780.75
      最大值0.930.910.860.890.910.86
      最佳任务乘法计算、乘法计算、乘法计算、乘法计算、乘法计算、乘法计算、
      4组合方式视觉计算视觉计算视觉计算视觉计算视觉计算视觉计算

      表 8可知,采用多通道的VAR相关方法在多数情况下的分类性能都优于单通道的AR模型.此外,文献[11]中AR相关方法使用的阶数是6,而本文中所有方法选择的阶数都是2,由此说明,多通道VAR的相关方法因其能挖掘出通道之间的互信息,可以在较短的时间范围内就得到最优的分类结果.对于受试者1、2、4,本文的方法得到的最佳任务组合方式与AR相同.

      在阶数为2的时候,VAR模型的参数空间还不太大,在这种情况下,我们并没有看出HVAR相关方法的优势.下面我们分析阶数为6时各种方法的实验结果. 图 10是4个受试者在阶数为6时五种方法对不同任务的分类性能的箱线图.

      图  10  阶数为 6 时不同方法的分类结果箱线图

      Figure 10.  The boxplots of all methods using six order for allsubjects

      图 10,我们可以发现: 当参数空间增大时,对于受试者1、3、4,在传统的VAR模型的基础上添加合适的分层结构有助于分类性能的提高,但是所有特征提取方法在受试者1 的实验中都出现了正偏的情况,也就是说对于不同的识别任务各种方法都没有太好的稳定性;从受试者3的实验结果来看,HVAR相关方法与VAR的性能基本持平,就最大值看HVAR的三种方法都要远高于VAR;对于受试者4,HVAR的性能远优于VAR,尤其以HVARC效果最好,值分布比较均匀,也就是说该方法对于不同的任务组合都有较好的识别正确率,稳定性很好. 总的看来,HVARO和HVARC两种方法的适用范围较广,在所有受试者上都能取得较好的实验结果,这说明当参数空间增大的时候,具有分层结构的VAR模型更有优势,这也反映出在多通道的信号处理中,不同通道的信号间的最长相互作用时间是不一样的,且自影响大于交互影响.但正如前面提到,每个人的思维方式不同,对于受试者2,HVAR模型的分类正确率反而不如VAR高. 此外,将EEG信号的不同通道不同时刻的值看作单一变量的LASSO-VAR方法,虽然能有效地压缩参数空间,但是在一定程度上会降低分类性能.

      图 10表明在参数空间大的时候,在传统VAR模型上添加罚函数的HVAR模型与VAR模型的性能相当,这是因为在阶数大的情况下,两种模型的参数求解误差都比较大,造成分类效果相当. 但是,我们知道VAR模型只有在 $T >kp$ 的情况下才能求解参数,故当时间序列之间的相互作用时间比较长的时候必然会造成VAR模型不再适用,而HVAR模型并不受该条件的约束,因此我们猜想虽然高阶HVAR模型对本文的脑电数据并没有实用性,但是基于其不受 $T > kp$ 条件的约束,该方法对处理多通道时间序列的分析会有一定的借鉴意义.

    • 在EEG识别过程中,特征提取起着非常重要的作用,提取出来的特征会直接影响系统的分类正确率.本文对多通道的EEG信号建立HVAR模型,将不同的分层结构模型得到的系数作为特征,并用Keirn等采集并整理的脑电数据进行检验,结果表明,在阶数远低于文献[11]给出的单通道AR的最优阶数时,多通道VAR相关方法就可以得到与AR模型相近的分类结果.在高阶的时候,参数求解误差较大,分类效果一般,但HVAR模型与VAR模型相比优势在于其不受 $T > kp$ 条件的约束,由于正则化思想的加入,该模型能够实现变量选择,剔除无关变量,压缩参数空间,在一定程度上降低参数求解误差,故对于相互作用时间较长的多通道时间序列,HVAR模型不失为一种新颖的方法.

参考文献 (29)

目录

    /

    返回文章
    返回