独立成分分析的变系数统计变量模型

独立成分分析的变系数统计变量模型

1sica和fmri信号分析方法独立成分分析（ica）是一种使用高统计量的盲源分析。研究表明,ICA能够在功能磁共振(functionalMRI,fMRI)数据中分析出各种生理的或者非生理的信号。目前,ICA对fMRI数据分析方法主要分为空间ICA(spatialICA,sICA)和时间ICA(temporalICA,tICA)两种,分别检测的是脑部空间激活和时域上动态响应。由于tICA数据维数高,计算量大,对数据信噪比敏感,即使仅处理单层数据也很难保证结果收敛等因素,所以目前较多使用sICA方法处理fMRI数据。在对fMRI信号ICA分析过程中,一般假设观测数据的信源信号是线性混叠,即由ICA分解得到的激活体元稀疏且不重叠,每个体元的fMRI信号是各个成分激活信号的线性混合。然而随着研究的深入,一些研究者提出构成fMRI信号的血氧水平依赖(bloodoxygenleveldependent,BOLD)等信源信号是非线性混叠,这样就不满足ICA统计隐变量模型的初始假设条件。BOLD信号是指血氧水平依赖性的增强效应所引起的信号,反映了血氧代谢过程,该信号的变化可以看为一个能量变化过程。如果用ICA分解脑部能量变化信号即瞬时功率信号,就可以得到反映了各脑部活跃区域能量的独立成分。同时,由于瞬时功率信号是标量,从而有效解决了ICA分解fMRI信号时所面对的需转化非线性混叠信源为线性混叠信源问题。本文的研究目的就是如何将fMRI数据转化为与BOLD信号的能量变化相关的信号,并采用sICA检测脑部能量异常活跃区域。2bold信号总能量的变化在表示BOLD信号的能量形式之前,首先来看一下电路中的能量表示形式。假设R表示电路中的电阻,在电路中保持恒定。于是,电路的总能量可以表示为E=∫U(t)I(t)dt=1R∫U2(t)dt=R∫I2(t)dt(1)E=∫U(t)Ι(t)dt=1R∫U2(t)dt=R∫Ι2(t)dt(1)式中,U(t)和I(t)分别表示瞬时电压和瞬时电流,t表示时间。如果将电阻看作一个固定系数,那么总能量表示形式就与∫U2(t)dt或∫I2(t)dt相关。电路的能量变化情况可以通过瞬时功率表示。P(t)=dEdt=U(t)I(t)=U2(t)R=I2(t)RΡ(t)=dEdt=U(t)Ι(t)=U2(t)R=Ι2(t)R(2)根据式(1),可类似地将fMRI的T*2信号的能量定义为ET*2=∫(x(t))2dt(3)式中,x(t)表示一个体素的fMRI时间序列。由于BOLD信号与T*2信号的变化相关,所以可以用fMRI数据的方差来表示BOLD信号的总能量:EBOLD=∫(x(t)−x¯)2dtEBΟLD=∫(x(t)-x¯)2dt(4)式中,x¯x¯表示x(t)的均值。根据式(2)的定义,进一步得到与BOLD信号总能量相关的瞬时功率:PBOLD(t)=∂EBOLD∂t=(x(t)−x¯)2ΡBΟLD(t)=∂EBΟLD∂t=(x(t)-x¯)2(5)因为BOLD信号Δ(t)在T*2信号中可以表示为(详细推导见附录)x(t)=x¯⋅eΔ(t)x(t)=x¯⋅eΔ(t)(6)所以BOLD信号Δ(t)可以表示为Δ(t)=ln(x(t))−lnx¯Δ(t)=ln(x(t))-lnx¯(7)由于BOLD信号在1.5T磁场下的变化量大约占总的T*2信号的2%,所以Δ(t)很小。根据泰勒级数,式(7)近似地表示为Δ(t)=lnx(t)x¯=ln(1+x(t)−x¯x¯)=x(t)−x¯x¯(8)Δ(t)=lnx(t)x¯=ln(1+x(t)-x¯x¯)=x(t)-x¯x¯(8)由式(7)和式(8)得到:x(t)−x¯=x¯⋅(ln(x(t))−lnx¯)x(t)-x¯=x¯⋅(ln(x(t))-lnx¯)(9)更进一步地,将式(5)改写为P′BOLD(t)=x¯2⋅(ln(x(t))−lnx¯)2Ρ′BΟLD(t)=x¯2⋅(ln(x(t))-lnx¯)2(10)式(5)和式(10)都能够表示BOLD信号能量的瞬时变化状况,从生理角度来讲,式(10)能够更好地反映BOLD信号能量变化。当x(t)→0时(如脑外部的噪声信号),ln(x(t))趋于无穷大,这不符合实际信号的变化。可以通过加入一个正常数ε>0,使得x(t)+ε>0来解决这个问题,即P∗BOLD(t)=x¯2⋅(lnx(t)+εx¯+ε)2ΡBΟLD*(t)=x¯2⋅(lnx(t)+εx¯+ε)2(11)3瞬时功率信号sica的信号特性为了定义ICA,文献采用了一种统计“隐变量(latentvariables)”模型。该模型假设n个观测数据x1,x2,…,xn是由n个独立成分c1,c2,…,cn线性混合得到,即xi=ai1c1+ai2c2+…+aincn(12)然而,观测数据中一般都会存在复杂的非线性混合情况(如乘性混合),ICA将无法准确地分析这类复杂情况。设一个观测数据xi是由m个信源s1,s2,…,sm构成,那么xi的能量Exi可以由各信源能量的线性叠加得到:Exi=Es1+Es2+…+Esm(13)根据式(2)和式(5)的定义,xi的瞬时功率表示为Pxi(t)=∂Exi∂t=∂Es1∂t+∂Es2∂t+⋯+∂Esm∂t=Ps1(t)+Ps2(t)+⋯+Psm(t)(14)Ρxi(t)=∂Exi∂t=∂Es1∂t+∂Es2∂t+⋯+∂Esm∂t=Ρs1(t)+Ρs2(t)+⋯+Ρsm(t)(14)因为观测数据的瞬时功率可以由信源的瞬时功率线性叠加得到,所以这种瞬时功率表示形式符合式(12)定义的统计隐变量模型。因此,瞬时功率信号适合ICA进行盲源分析,并且检测得到的成分与BOLD反映的能量相关。具体的,瞬时功率信号sICA步骤如下:(1)采用统计参数图法(statisticalparametricmapping,SPM)对fMRI数据进行空间预处理,包括:①头动校正;②消除间隔扫描的层间差异;③配准到标准空间坐标(可选);④空间平滑(可选)。(2)将4维fMRI数据(3维脑功能数据×采集时间)转化为由空间和时间构成的2维矩阵Fspace×time,矩阵中每个元素表示大脑空间某个体素在时间t时的采样值。(3)根据式(5)或式(11)计算每个空间位置的瞬时功率Pspace×time,并计算瞬时功率的转置Xtime×space=PTtime×spacetime×spaceΤ。(4)采用sICA,将信号Xtime×space分解为相互独立的空间成分Ccomponent×space和与时间信息相关的混合矩阵Atime×component两部分:Xtime×space=Atime×componentCcomponent×space(15)(5)对每个独立成分Ccomponent×space的Z值的绝对值(虽然ICA得到的Z值存在正负,但是由于瞬时功率是标量,所以本文对Z值取绝对值。)取阈值,并分别映射到对应的空间位置。称此处理方法为基于fMRI瞬时功率的空间独立成分分析方法(instantaneouspowerbased-sICA,IP-sICA)。4仿真数据处理为了验证方法的有效性及可靠性,选用一名健康被试在无刺激情况下采集的单一64×64的EPIfMRI图像作为模板,重复100次,构成79×95×100的仿真数据体,然后在其上两个大小不等的ROI(regionofinterest)区域(M:13200个体素,N:30600个体素)叠加了与标准HRF(hemodynamicimpulsefunction)卷积后的事件相关序列信号(图1(a)),其幅度值为背景图像像素平均强度的2%。同时在背景数据体上叠加79×95×100个高斯白噪声,信噪比(SNR)由0递增至1(不含0),步长0.01。为了提高仿真数据的稳定性,重复以上过程5次。对仿真的fMRI数据分别用IP-sICA方法(由于式(5)和式(11)结果一致,以下结果均采用式(5))和传统sICA方法(ICA核心算法选取fastICA参数选择Approach为symmetric,Nonlinearity为gauss,stabilization为on)处理,计算出精确度、真阳性率(truepositivefraction,TPF)和假阳性率(falsepositivefraction,FPF)值,同时计算出受试者工作特征(receiveroperatingcharacteristic,ROC)power的值,绘制图形得到如下比较结果。图2给出了在4个SNR(0.22,0.46,0.68,0.90)和不同置信度p下两种方法的精确度对比图,新方法较传统sICA方法在不同的置信度下具有更高的精确性,尤其在p>0.05时。图3给出了两种方法的ROCpower值。在同一阈值下,ROCpower的值越高,检测激活的空间精确性越高,检测的效果越好。从仿真数据处理的结果来看,对不同的SNR,提出的算法准确有效地找出了模拟的空间激活区域,同时较传统sICA方法具有更高的精确性。计算出各SNR下ROCpower曲线下面积值(areaunderROCcurve,AUC),当SNR分别为0.22,0.46,0.68,0.90时,IP-sICA和传统sICA得到的AUC值分别为0.9950,0.9365,0.9805,0.9523,0.9877,0.9074和0.9697,0.9075。尤其当SNR为0.22时,IP-sICA方法得到的AUC值几乎为1,近似理想情况。同时发现在重复的5次仿真分解中,IP-sICA方法分解的所有成分中有90%以上的成分空间分布于2个ROI中,而传统sICA分解的所有成分中只有1个或2个成分空间分布于2个ROI中,从此可以看出新方法更加容易提取fMRI的sICA成分。5功能像及数据采集和分析癫痫发作是由于异常脑电生理改变所致,同时引起脑血氧能量代谢异常,所以本研究纳入南京军区南京总医院采集的3例局灶性癫痫患者。这些患者发作类型均为部分性发作,并且颅内脑实质有明确结构异常,临床已经确诊癫痫灶。结构像与功能像数据在GE-Signa1.5T磁共振仪下进行采集。结构像采用T1Flair序列,TR=2019.3ms,TE=25.3ms,层厚4mm,层间距0.5mm,平行前后联合水平进行全脑扫描,共21层。功能像采用GRE-EPI序列采集,TR=2000ms,TE=40ms,FA=80°,共200个时间点,采集方位与结构像相同。被试者平静放松并用海绵垫固定其头部,清醒状态闭眼平卧于MR仪内。为了观察结果的稳定性,每人扫描2次功能像,其间间隔3～5min。在对3例癫痫被试者fMRI数据分析前,通过主成分分析(principlecomponentanalysis,PCA)进行降维,对方差和占99%以上的主成分进行sICA分解,成分数目为50。提出方法有效提取3例患者颅内病灶周围激活区域,而传统sICA方法无法正确提取独立成分,从而限制其对能量异常的癫痫fMRI数据的应用。图4显示了一名被试者的一个断面(+32/MNI305)fMRI空间激活对比图。其中图4(a),图4(c)分别为IP-SICA分解和传统SICA分解的50个空间独立成分叠加到结构像上的激活图,可以看出IP-sICA方法分解的90%以上fMRI信号成分定位于癫痫病灶周围(右枕顶叶,白实心箭头),而传统sICA算法分解的各个成分结果分散,在没有预先假设信息的条件下,无法正确提取sICA成分。为了便于对比,图4(b),图4(d)分别为图4(a)、图4(c)中50个独立成分的Z值平均后的结果,新方法(图4(b))可以直观地得出癫痫灶。6基于信号正则化的fmri信号模型ICA盲分析方法在fMRI数据分析领域的应用日趋成熟,可在不需要任何有关时间序列的先验假设下,将多道观察信号按照统计独立的原则,通过优化算法分解为若干相互独立成分,有效地探测出一些其他传统假设驱动方法所得不到的激活,例如静息状态下fMRI激活提取。它的建立基础在于假设fMRI信源信号的统计独立性。随着研究的深入,发现fMRI实际观测信号是由各种信源信号包括BOLD、线性漂移、生理噪声等信号非线性混叠而成,从而不满足ICA“隐变量”模型(式(12))。同时由于无法对ICA分离的各分量排序及正确选取等客观因素,限制其在fMRI应用领域里的进一步深入。针对ICA对fMRI信号处理中的不足,提出了一种基于fMRI信号瞬时功率的ICA模型方法:由两种不同的正则化方法对非线性fMRI信号进行正则化,这不仅使得信号自动满足I