梅 戬 王小宇 刘 杨 李景琦 刘克洪 杨 勇* 丛丰裕*
1(杭州电子科技大学自动化学院,杭州 310018)
2(大连理工大学生物医学工程学院,辽宁 大连 116024)
3(浙江音乐学院研究生院,杭州 310024)
4(杭州明州脑康康复医院神经康复科,杭州 310000)
5 (武警浙江省总队杭州医院康复科,杭州 310051)
大脑对音乐的加工是一个从提取声音低级特征再抽象出高级音乐结构的层级神经加工过程[1]。近年来,音乐感知与意识的关系成为认知科学领域的关注点之一,对音乐感知与意识关系的探索有助于揭示人类脑认知活动的神经机理,同时对音乐治疗的临床应用也具有重要意义。
以往对音乐感知的研究都基于意识清晰状态,如利用清醒被试的功能磁共振(functional magnetic resonance imaging, fMRI)研究音乐激活的大脑区域[2]、利用事件相关电位(event-related potential, ERP)研究音乐感知的ERP响应在波形和节律上的变化[3],以及研究意识清晰脑损伤的患者听音乐时的神经活动[4]等。随着微意识状态(minimally consciousness state, MCS)、植物状态(vegetative state, VS)概念的提出,作为天然的意识障碍模型,这两种意识状态为音乐感知与意识关系探索提供了新的模型和研究范式。以往对不同意识状态被试施加声音刺激的研究发现,植物状态被试、微意识状态被试和正常人,在受声音刺激时的脑电信号(electroencephalogram, EEG)特征[5]、失匹配负波(mismatch negative, MMN)波形[6]和功能磁共振反应[7]的脑皮层反应均有差异。
近年来涌现了许多新颖的研究方法,包括脑电信号处理方法、音乐信息学特征分析方法等,为探索音乐感知与意识关系提供了新视角。丛丰裕等利用张量分解,提取了听音乐时的脑电信号的多维特征,并验证了可行性[8]。脑电信号是大量神经信号叠加而成的混合信号,因此在研究音乐特征与脑电信号的关系之前,必须首先从脑电中提取出与音乐特征相关的神经活动成分。张量分解是一种多维的盲源分离方法,它将脑电信号分解成多个独立的具有多维特征的源成分,每个独立的源成分就是一种独立的神经活动[9]。相对其他盲源分离方法(如独立成分分析),张量分解可以在多维的特征空间利用数学模型算法,将采集到的混合信号分离为多种独立的张量成分,同时保留信号的多维特征,如时域、频域、空间分布[10-12]。最后,通过计算成分和音乐的时域特征相关性,筛选与音乐特征相关的神经活动,进而获得响应音乐特征的脑电信号成分。
此外,以往研究音乐特征多采用声学特征(音强、音高、音色)。2012年,从信号角度,基于音乐的局部旋律和音域特性,Alluri等提出了音乐时域信号特征的概念,包括fluctuation centroid、fluctuation entropy、pulse clarity、key clarity和 mode[13](详见本文附录介绍)。Alluri结合这些音乐特征,研究被试听音乐时的fMRI,发现大脑在处理这些音乐特征时,一个涉及认知、运动和脑边缘回路的大规模脑活动被激活[13]。这将音乐感知的研究从传统的声学范式引入到了信号学层面,并为音乐加工的脑机制研究范式提供了新的视角。
基于脑电信号分析的新技术及音乐信息学(信号)特征,本研究开展不同意识状态对音乐感知的差异性研究,揭示意识状态如何对大脑加工音乐产生影响。通过设计音乐刺激实验,采集相关脑电信号,利用基于层次交替最小二乘法(hierarchical alternating least squares, HALS)的非负CP张量分解算法,从采集的脑电中分离与音乐相关的张量成分。最后,对每个脑电张量成分进行比率功率谱分析和空间分布分析,从而了解不同意识水平对大脑加工音乐特征的影响。
本研究的微意识状态组和植物状态组被试的脑电信号均采自于杭州明州脑康康复医院和武警浙江省总队杭州医院。所有入选被试均无需插管由呼吸机辅助呼吸,无心肺复苏史,无精神疾病史。此外,为了在研究期间保证被试意识状态稳定及脑电特征的一致性,选用处于微意识状态或植物状态一个月以上且处于慢性期的被试[14]。同时,剔除以下病例特征:存在中度及以上听力损失患者、闭锁综合症患者,以及患有可能导致脑神经病理性缺失疾病的患者。本研究的意识水平评估,采用国际对微意识状态和植物状态区分效果较好的CRS-R量表作为意识水平量化工具,筛选出其6项得分对应小于4-5-6-3-2-3的被试。被试为3组正常被试组、微意识状态组、植物状态组。正常被试组7人,年龄分布20~30岁。微意识状态组17人,植物状态组19人,年龄分布20~55岁。微意识状态组、植物状态组被试的详细信息如表1所示。
表1 微意识状态组和植物状态组被试信息统计Tab.1 Information about MCS and VS participants
实验采用BioSemi公司生产的64导联ActiveTwoSystem脑电采集仪器,采样率为2 048 Hz。数据采集采用单极导联方式,选用A1、A2作为参考电极。根据国际10/20标准导联系统,将电极布置在被试头皮表面,在电极和头皮间注入导电凝胶,使电阻小于5 kΩ。
被试进入安静状态后,开始脑电采集实验。首先采集被试在静息态时的脑电信号60 s,然后采集被试处于音乐刺激状态下的脑电信号120 s,最后再采集被试刺激后的静息态信号60 s。本研究的刺激源选取了《茉莉花》的高潮副歌部分,时长为120 s,播放的声级为70 dB。实验中的被试是闭眼听音乐,对于不能闭眼的意识障碍被试,用毛巾遮住他的眼睛。采集过程中,无噪声和大功率用电设备干扰。被试呈静卧姿势,房间温度控制在25℃。实验的过程如图1所示。
图1 脑电信号采集流程Fig.1 Time-course of EEG acquisition.
本研究在Matlab 2016a平台上,利用EEGLAB version 13_5_4b对脑电信号进行分析。首先,对脑电进行去工频滤波,根据我国电流频率,设置去工频50 Hz;其次,根据脑电的节律特点,将带通滤波的截止频率设为低通80 Hz、高通0.5 Hz。根据EEGLAB观察滤波后的波形,剔除较大扰动和肌肉伪迹,从而得到降噪后的脑电信号。
降噪后,信号依然会存在伪迹,如眼电(electro-oculogram, EOG)等。基于Matlab 2016a,利用ICASSO工具箱做基于InfomaxICA的信号分解算法,对降噪后的脑电实施独立成分分析(independent component analysis, ICA),以去除伪迹。
利用Matlab2016a编译环境,加载MIRToolbox Version1.7.1工具箱,提取刺激源《茉莉花》的5个音乐特征,包括fluctuation centroid、fluctuation entropy、pulse clarity、key clarity和 mode[13, 15]。首先采用移窗法,利用窗宽3 s,与前后窗重叠各2 s,从零秒开始对音乐取样[16];然后利用MIRToolbox,计算每个音乐片段的音乐特征(fluctuation centroid、fluctuation entropy、pulse clarity、key clarity和 mode)的值。
选用脑电信号的空间分布(channel)、时域(time)和频域(frequency),作为三阶张量的特征。为获取每个被试脑电信号各导联的时频域特征作为张量第三阶分量,选用短时傅里叶变换对信号做时频域分析[11]。利用汉明窗对各导联脑电信号采样,窗宽3 s,与前后窗重叠各2 s。由此生成每个脑电信号的三阶张量大小为channel×frequency×time,其中channel=64,frequency=158,time=120-伪迹时间,伪迹时间为预处理时信号中因伪迹较大而被剔除部分的时长。
提取信号高维张量成分的张量分解算法一般有两种模式:canonical polyadic(CP)分解和tucker分解。CP分解是将一个给定的N阶张量X分解为一系列秩一张量的和的过程。
根据阶数不同,张量可以被划分为一阶向量、二阶矩阵和三阶及以上的高阶张量。与矩阵分解类似,一个N阶张量X可以分解为R个N阶秩一张量的和[17]。每个秩一张量就是张量的一个成分。N阶秩一张量等于N个正交单位向量的外积乘以能量系数。由此,得到一个N阶张量X的CP分解如下:
(1)
式中,张量X为一个N阶张量,λr为能量系数,ar1,ar2,…,arN为正交单位向量,符号∘表示张量外积。
本研究选用非负CP张量分解,用来分解在本文第3.1节中获得的三阶张量。在此过程中,不断趋近一系列最优非负正交向量组合,使得剩余张量范数‖E‖2趋于零。剩余张量二范数‖E‖2的计算过程如下:
(2)
为了加快剩余张量二范数的收敛,减少对算力的依赖,本研究选用了由层次交替最小二乘法(hierarchical alternating least squares, HALS)优化趋近过程的基于HALS的非负 CP张量分解[18]。
本文第3.2节介绍了选用的基于HALS的非负 CP张量分解算法。下面将利用该算法分解本文第3.1节中得到的每个被试的三阶信号张量,以获得脑电活动的成分。从信号张量提取出张量成分的个数,通过Smooth DIFFIT判断得出[19]。
以上方法提取出的张量成分都含有3个分量,分别是空间分布系数、时域包络及频谱。根据张量成分的空间分布系数画出脑地形图,然后筛选出满足dipolar形式的张量成分,即为脑电张量成分[10, 20]。然后将张量成分中的时域包络分量依次与方法2中获得的5种音乐时域特征值做Pearson相关系数分析,用Monte Carlo 法计算相关成立的相关系数阈值大小,选出与音乐特征相关的成分(P<0.05),即为响应该音乐特征变化的脑电活动[8, 11, 13]。
在本文第3.3节中得到的脑电张量成分具有空间分布系数、时域包络及频谱3个分量,其中脑电张量成分的频域分量体现了该成分的频谱特性。根据频域分量可计算其比率功率谱,即各节律的占比,用来分析不同意识水平被试的脑电信号对音乐特征的响应在节律上是否存在差异性。各个节律的比率功率谱计算如下:
(3)
式中,f1和f2分别为frequency band of interest的频域下限与上限,即delta (1~4 Hz)、theta(4~8 Hz)、alpha (8~13 Hz) 和 beta (13~30 Hz)。
脑电张量成分的空间分布分量表征了该成分在脑区的分布,系数越高说明该采样点离产生该脑电张量成分的源成分的距离越近,从而可通过脑地形图判断源成分发生的响应脑区位置。
为了验证正常被试、微意识状态组和植物状态组之间的响应差异,使用SPSS v.20 (SPSS Inc., Chicago, IL) 软件,对脑电张量成分的功率谱比率做ANOVA单因素检验和Scheffe事后检验,P<0. 05被认为统计上具有显著差异。
最后比较与音乐特征不相关和相关的脑电张量成分在比率功率谱和响应脑区分布的差异,排除3组被试在功率谱和脑地形图的差异与音乐感知无关,可能仅是本身意识水平的不同而引起的。
通过本文第3.2节中的张量分解算法和第3.3节中的张量成分筛选方法,从7例正常人组中筛选出18个成分,从17例微意识状态组筛选出41个成分,从19例植物状态组筛选出32个成分,结果如表2所示。
表2 不同意识水平组的脑电与5种音乐特征相关张量成分个数统计Tab.2 Statistics on the number of EEG tensor components correlated with five musical features among different consciousness level groups
按本文第3.4节的方法,对所有音乐特征相关的脑电张量成分做比率功率谱分析。表3~7分别为与fluctuation centroid、fluctuation entropy、pulse clarity、key clarity和 mode音乐特征相关的脑电张量成分的比率功率谱分析结果。
表3为与音乐特征fluctuation centroid相关的脑电张量成分比率功率谱的平均值。其中,正常人组脑电张量成分的alpha波和beta波占比高于theta波和delta波,而微意识状态组脑电张量成分和植物状态组脑电张量成分theta波和delta波占比高于alpha波和beta波。3组被试的脑电张量成分的alpha波占比和beta波占比存在统计差异(alpha波:F(2,21)=104.838,P<0.001;beta波:F(2,21)=10.418,P=0.001)。正常人组脑电张量成分的alpha波占比高于微意识状态组(P<0.001)和植物状态组(P<0.001),正常人组脑电张量成分的beta波占比高于微意识状态组(P<0.001)和植物状态组(P<0.001),微意识状态组脑电张量成分和植物状态组脑电张量成分在各节律波段占比没有差异。
表3 与音乐特征fluctuation centroid相关的脑电张量成分比率功率谱平均值Tab.3 The mean value of the power spectrum ratio of EEG tensor components correlated with the music feature fluctuation centroid
表4为与音乐特征fluctuation entropy相关的脑电张量成分比率功率谱的平均值。由于正常人组和植物状态组与fluctuation entropy相关的脑电张量成分各只有两个,所以不对其做独立样本统计分析。但是可以观察发现,正常人组脑电张量成分的alpha波和beta波占比高于theta波和delta波,微意识状态组和植物状态组的脑电张量成分相反。正常人组脑电张量成分的alpha波占比高于微意识状态组和植物状态组的alpha波占比,微意识状态组脑电张量成分和植物状态组脑电张量成分在各节律波段占比没有差异。
表4 与音乐特征fluctuation entropy相关的脑电张量成分比率功率谱平均值Tab.4 The mean value of the power spectrum ratio of EEG tensor components correlated with the music feature fluctuation entropy
表5为与音乐特征key clarity相关的脑电张量成分比率功率谱的平均值。其中,正常人组脑电张量成分的alpha波和beta波占比高于theta波和delta波,而微意识状态组脑电张量成分和植物状态组脑电张量成分theta波和delta波占比高于alpha波和beta波。3组被试的脑电张量成分的alpha波占比和beta波占比存在统计差异(alpha波:F(2,18)=27.349,P<0.001;beta波:F(2,18)=6.758,P=0.006)。正常人组脑电张量成分的alpha波占比高于微意识状态组(P<0.001)和植物状态组(P<0.001),正常人组脑电张量成分的beta波占比高于微意识状态组(P=0.010) 和植物状态组(P=0.022),微意识状态组和植物状态组的脑电张量成分在各节律波段占比没有差异。
表5 与音乐特征key clarity相关的脑电张量成分比率功率谱平均值Tab.5 The mean value of the power spectrum ratio of EEG tensor components correlated with the music feature key clarity
表6为与音乐特征pulse clarity相关的脑电张量成分比率功率谱的平均值。其中,常人组脑电张量成分的alpha波和beta波占比高于theta波和delta波,而微意识状态组脑电张量成分和植物状态组脑电张量成分theta波和delta波占比高于alpha波和beta波。3组被试的脑电张量成分的alpha波占比和beta波占比存在统计差异(alpha波:F(2,21)=27.349,P<0.001;beta波:F(2,18)=6.758,P=0.006)。正常人组脑电张量成分的alpha波占比高于微意识状态组(P<0.001)和植物状态组(P<0.001),正常人组脑电张量成分的beta波占比高于微意识状态组(P=0.01)和植物状态组(P=0.022),微意识状态组和植物状态组的脑电张量成分在各节律波段占比没有差异。
表6 与音乐特征pulse clarity相关的脑电张量成分比率功率谱平均值Tab.6 The mean value of the power spectrum ratio of EEG tensor components correlated with the music feature pulse clarity
表7为与音乐特征mode相关的脑电张量成分比率功率谱平均值。其中,正常人组脑电张量成分的alpha波和beta波占比高于theta波和delta波,而微意识状态组脑电张量成分和植物状态组脑电张量成分的theta波和delta波占比高于alpha波和beta波。3组被试的脑电张量成分的alpha波占比存在统计差异(alpha波:F(2,6)=5.79,P=0.040)。微意识状态组的脑电张量成分和植物状态组的脑电张量成分在各节律波段占比没有差异。
表7 与音乐特征mode相关的脑电张量成分比率功率谱平均值Tab.7 The mean value of the power spectrum ratio of EEG tensor components correlated with the music feature mode
根据本文第3.3节中得到的脑电张量的空间分布系数分量,可画出脑地形图,并判断该脑电张量成分的响应脑区。表8~10分别为正常人组、微意识状态组和植物状态组与音乐特征相关的脑电张量成分在脑区的分布统计。正常人组的脑电张量成分除2个位于颞叶,其余16个都分布于前额叶;微意识状态组的脑电张量只有6个分布于前额叶,其余的主要位于颞叶;植物状态组的脑电张量成分只有1个分布于前额叶,其余的主要分布于颞叶。如图2所示,从左到右分别为正常人组响应脑区是前额叶的脑电张量成分的脑地形图、微意识状态组响应脑区是颞叶的脑电张量成分的脑地形图和植物状态组响应脑区是颞叶的脑电张量成分的脑地形图。
多种抗心律失常药物及电复律治疗对左前分支室速效果欠佳,射频消融治疗是唯一根治性方法[7]。部分患者长时间心动过速可出现临床症状,甚至可诱发心动过速性心肌病[8]。尤其对于儿童患者,由于年龄小、体质量也小,射频消融风险相应升高[9-10],在选择射频消融根治术前,食管调搏作为一种无创的电生理检查方法,对于特发性室速的鉴别与心动过速的终止均有重要的意义,因此,食管调搏可作为紧急终止分支型室速的方式之一。
图2 三组被试的脑电信号与音乐特征相关的张量成分脑地形图。(a)正常人组;(b)微意识状态组;(c)植物状态组Fig.2 Spatial component of EEG tensor components correlated with musical feature. (a) Normal group; (b) MCS group; (c) VS group
表8 正常人组与音乐特征相关的脑电张量成分响应脑区分布统计Tab.8 The number of location of EEG tensor components correlated with musical feature in normal group
本节把3组被试脑电中与音乐特征不相关的脑电张量成分进行比率功率谱和响应脑区分布统计,以排除3组被试的功率谱和脑地形图只是由本身意识水平的不同引起而与音乐感知无关的可能性。正常人组145个脑电张量成分,微意识状态组202个脑电张量成分,植物状态组307个脑电张量成分。
表9 微意识状态组与音乐特征相关的脑电张量成分响应脑区分布统计Tab.9 The number of location of EEG tensor components correlated with musical feature in MCS group
表10 植物状态组与音乐特征相关的脑电张量成分响应脑区分布统计Tab.10 The number of location of EEG tensor components correlated with musical feature in VS group
表11为3组被试与音乐特征不相关的脑电张量成分比率功率谱平均值。delta波占比:植物状态组(0.629±0.372),高于微意识状态组(0.436±0.385),高于正常人组(0.183±0.265)(F(2,640)=77.179,P<0.001);theta波占比:微意识状态组(0.240±0.298)和植物状态组(0.213±0.293)高于正常人组(0.124±0.157)(F(2,640)=8.494,P<0.001);alpha波占比:正常人组(0.483±0.320)高于微意识状态组(0.081±0.121)和植物状态组(0.056±0.092)(F(2,640)=308.538,P<0.001);beta波占比:正常人组(0.159±0.105)高于微意识状态组(0.025±0.021)和植物状态组(0.029±0.036)(F(2,640)=299.129,P<0.001)。
表11 与音乐特征不相关的脑电张量成分比率功率谱平均值Tab.11 The mean value of the power spectrum ratio of EEG tensor components not correlated with the music feature
表12为3组被试与音乐特征不相关的脑电张量成分的比率功率谱和响应脑区分布。正常人组145个脑电张量成分中,有97个的响应脑区位于前额叶;微意识状态组202个脑电张量成分中有44个的响应脑区位于前额叶,植物状态组307个脑电张量成分中有51个的相应区域位于前额叶。
表12 与音乐特征不相关的脑电张量成分响应脑区分布统计Tab.12 The number of location of EEG tensor components not correlated with musical feature
将本文第4.3节中与音乐特征不相关的脑电张量成分的比率功率谱与响应脑区,和第4.2节中与音乐特征相关的脑电张量成分的比率功率谱与响应脑区相比较。首先,微意识状态组中与音乐特征相关的脑电张量成分的delta波占比和植物状态组没有差别,而微意识状态组中与音乐特征不相关的脑电张量成分delta波占比低于植物状态组的占比;其次,正常组与音乐特征相关的脑电张量成分的theta波占比低于微意识状态组与植物状态组的占比,而正常组与音乐特征不相关的脑电张量成分theta波占比和微意识状态组与植物状态组存在差别;最后,微意识状态组和植物状态组的与音乐特征不相关的脑电张量成分的响应脑区位于前额叶比例高于与音乐特征相关的脑电张量成分,而正常人组与音乐的时域特征相关的脑电张量成分的响应脑区位于前额叶比例高于与音乐的时域特征不相关的脑电张量成分。
因此,对比与音乐时域特征不相关和相关的脑电张量成分可知,这两组成分的响应脑区分布和功率谱在3组被试间的统计结果是存在差异的。基于此,根据实验结果得出结论,3组被试在功率谱和脑地形图的差异与音乐感知是有关的,并不仅仅是本身意识水平的不同而引起的。
首先,本研究对不同意识水平的被试在音乐刺激下产生的脑电信号做降噪和去伪迹预处理,然后,基于HALS的非负 CP张量分解,对处理后的信号做张量分解,获得脑电的张量成分。将这些张量成分的时域分量分别与5种音乐特征(fluctuation centroid、fluctuation entropy、pulse clarity、key clarity和 mode)做相关分析,提取出与音乐特征相关的脑电张量成分,分析功率谱和响应脑区分布,最后探究不同意识水平对音乐特征加工的差异。
由3组被试脑电中与音乐特征相关的脑电张量成分的功率谱和响应脑区分布的统计结果可知,被试在音乐刺激下的脑活动在节律和脑响应区域存在差异。在节律方面,通过对张量分解得到的成分做比率功率谱分析发现:除特征mode外,正常意识状态的响应脑电张量成分alpha波和beta波的占比均高于微意识状态和植物状态。微意识状态和植物状态之间的各节律占比没有明显差异。正常意识状态的响应脑电张量成分的alpha波和beta波的占比均高于delta波和theta波,而微意识状态和植物状态则是delta波和theta波占比高于alpha波和beta波。基于以往对脑电节律的研究,注意力集中(attention)和高强度认知活动(cognition)会引起beta波的出现[21],并且本研究发现,微意识状态和植物状态脑电中与音乐特征相关的脑电张量成分中beta节律明显少于正常意识状态组,可以推测这两组意识状态对音乐特征的感知和注意力更弱。
在响应脑区方面,正常意识状态时的响应脑区主要集中在前额叶,而意识降为微意识状态和植物状态时,响应主要发生于颞叶。这意味着随着意识水平降低,前额叶响应下降,而颞叶基本保持,佐证了前额叶对意识的重要性。针对音乐加工的脑机制,Koelsch等[1]提出了层级音乐加工模型。赖寒等[22]结合该模型,对大脑对音乐旋律更精细加工的第3阶段音程分析和第4阶段结构分析的神经活动做元分析,发现第3阶段音程分析激活的主要区域为双侧颞上回,而第4阶段结构分析激活的主要区域为前额叶,包括额下回、额中回和额极。根据这个理论,可推断正常意识状态对音乐进行的第4阶段结构加工并没有发生在微意识状态和植物状态中,因为处于音乐加工第四阶段的结构加工的激活区域主要在前额叶,但是微意识状态和植物状态被试的前额叶并没有被激活[22]。
其次,本研究对比正常意识状态、微意识状态和植物状态发生于颞叶的成分的频域特征,发现3组被试的响应音乐特征的脑电张量成分在节律上分布不同。正常意识状态被试对音乐特征响应的脑电张量成分的alpha波和beta波依然高于delta波和theta波,而微意识状态和植物状态对音乐特征响应的脑电张量成分以theta波和delta波为主。这说明,随着意识水平的降低,颞叶对音乐特征的响应产生的脑电节律也会降低。
研究发现,正常被试、微意识状态被试和植物状态被试对音乐特征的响应在节律和响应区域存在差异。在脑激活区域分布上,正常被试对音乐特征的响应主要集中在前额叶,而微意识状态被试和植物状态被试对音乐的响应集中在颞叶。在节律分布上,正常被试对音乐特征的脑电响应频率集中在alpha和beta波段,微意识状态被试和植物状态被试对音乐特征的脑电响应频率集中在theta和delta波段。
附录:音乐特征
5个音乐特征Fluctuation centroid、fluctuation entropy、pulse clarity、key clarity和mode是由Alluri等在2012年提出的[13],其中前3个反映音乐的节奏特征(rhythmic feature),后2个反映音乐的音域特征(tonal feature)。
这5个音乐特征的物理意义如下:fluctuation centroid是波动频谱的几何平均数,用来计算音乐0~10 Hz内的平均频率;fluctuation entropy是波动频谱的香农熵,用来表征波动频谱的混乱度,值越大包含的旋律越多;pulse clarity用来衡量旋律的清晰度,表明重音能量的大小;key clarity用来衡量音阶的清晰度;mode用来表示大调或小调的主导率。