基于多目标优化约束独立成分分析方法的fMRI数据分析研究

石玉虎 曾卫明 邓 金 王倪传

1(上海海事大学信息工程学院，上海 201306) 2(淮海工学院计算机工程学院, 江苏 连云港 222023)

引言

利用功能磁共振成像(functional magnetic resonance imaging, fMRI)技术来研究人脑的功能连通性及相关的认知活动特性，可以更好地揭示人脑相关的认知活动规律，促进fMRI技术在神经信息处理以及神经认知应用等领域发挥重要作用。目前，关于fMRI脑功能连通性的研究方法主要分为基于模型驱动的方法和基于数据驱动的方法。其中，基于模型驱动的方法主要有互相关分析方法[1]、一致性分析方法[2]以及统计参数图分析方法[3]等，但这些方法均受限于预先给定的先验信息。为了克服此种局限性，便产生了诸如聚类分析[4-7]、稀疏成分分析[8-10]及独立成分分析(independent component analysis, ICA)[11-12]等数据驱动方法。与模型驱动方法相比，数据驱动方法因不需要任何的模型假设条件而被广泛应用于fMRI信号分析的研究。这些方法在fMRI信号处理上展现出较高的准确性与实用性，在具体应用中有较大的优势，但同时也存在不足之处。

ICA是一种基于高阶统计量的盲源信号分离方法，其目标是利用高阶统计量的概率统计特性，在假定源信号相互独立非高斯信号的基础上，将采集到的混合信号分解为若干相互独立的成分[13]。自从1998年McKeown等首次将该方法应用于fMRI数据分析以来[14]，ICA已经成为该领域最为流行的方法之一，被广泛应用于fMRI脑功能连通性检测研究[15-21]。作为一种数据驱动分离技术，ICA无需依赖于任何先验信息，因而具有独特的优势。但同时研究也发现，ICA在fMRI分析过程中具有其自身的局限性，如成分个数不确定、成分输出无序的问题等，由此可能导致大量计算时间和计算资源的浪费，特别是当待分析数据量比较大时问题就更加明显。

针对上述问题，研究人员提出通过在ICA模型中加入可利用先验信息来加以克服，而研究结果表明：在ICA的计算过程中引入先验知识，确实能够提高ICA的fMRI数据分析性能[22-26]。迄今为止，学者们已经提出了许多方法，以解决关于ICA中包含先验信息的问题。例如：Lu和Rajapakse在ICA的基础上，率先提出了一种约束独立成分分析(constrained ICA, CICA))来解决ICA融合先验信息的问题[27]，并进一步提出了基于增广拉格朗日方法和牛顿迭代算法的ICA with reference (ICA-R)方法[28]。随后，Lin等在CICA和ICA-R方法的基础上，选用不同的度量函数将空间先验信息引入ICA模型，并采用快速不动点算法进行求解，提出了所谓的半盲空间ICA方法，并在实际中表现出良好的分析能力[29]。与ICA方法相比，加入先验信息的CICA方法可以只提取所需要的源信号信息，这样不仅有利于后续的应用分析，而且还能减少算法计算过程中所需要的时间和存储容量[30]。另外，加入先验信息能够提高ICA算法的计算性能，使它更加适应于实时数据的分析环境[31-32]，并提高感兴趣源信号检测的质量和精度[33]。

然而，目前的CICA方法都需要一个阈值参数来约束输出信号与先验参考信号之间的相似性，使得算法的输出结果是与参考信号最为接近的唯一解。而源信号的不可知性，使得确定一个合适的阈值参数来约束输出信号与先验参考信号之间的近似度十分不容易，不合适的阈值往往会导致截然不同的结果，因此发展一种不依赖于阈值参数选择的包含先验信息的ICA方法显得十分必要。最近，Du等提出了一种多目标的优化策略来实现传统CICA方法的求解问题，通过将传统CICA方法中包含先验信息的不等式约束条件转变为一个优化目标函数，与原始ICA的独立性目标度量函数构成多目标优化策略，这样就规避了传统CICA中阈值参数选择不当导致分解结果不理想的问题[34]。

此外，目前能够嵌入的先验信息通常是一些具体且可直接利用的先验知识，如实验中激励任务的刺激模式、研究比较成熟的脑网络模板(如视觉网络)、源信号的某些特定属性等。但是，这些先验信息并不总是存在的，或在某些情况下可能存在却不易获得，比如在静息状态或复杂认知任务状态下的fMRI数据就难以获取这样的先验信息。正是由于这些先验信息难以获取，使得一些研究不得不放弃CICA方法，转而使用传统的ICA方法恢复所有的独立源信号，然后进行后期选择。由此看来，挖掘数据隐含的先验信息，不仅是为了克服CICA中先验信息难以获取的问题，而且对于CICA方法的广泛应用也十分有益。因此，如何从已有条件中获取可利用的隐含先验信息就显得尤为重要。

无论是从实验中直接获取的先验信息，还是通过数据驱动方法挖掘的隐含先验信息，它们都只是提供了一些有关源信号的粗糙信息，并不可能百分之百正确，而不同精度的先验信息对提高ICA算法的分析能力也是不同。因此，如何减少CICA方法对参考信号的精度要求，用尽量少的先验信息来最大程度提高ICA的分析能力，是一个值得研究且非常有意义的问题。此外，在有些可利用的先验信息里面，可能还会或多或少含有一些错误信息，从而会误导算法优化过程中的收敛方向，使得输出的结果不是所期望的结果。因此，如何提高包含先验信息ICA算法的信息容错能力，减少先验信息中错误信息对算法的不利影响，是一个亟待解决的关键问题。Wang等研究发现，同时利用时间和空间先验信息，要比单独使用其中某种先验信息能够使CICA算法更具有鲁棒性，减少了CICA方法对单一先验信息的准确度要求[35]。

因此，本研究在多目标优化框架的基础上，提出了一种同时融合时空先验信息的多目标优化CICA (multi-objective optimization based CICA, MOPCICA)方法，有效规避了CICA方法中阈值参数难以选择的问题。通过引入了时间和空间先验信息，不仅降低了CICA方法对先验信息的精度要求，而且提高了对先验信息中错误信息的容错能力。实验结果表明，MOPCICA方法在时间和空间源信号恢复能力方面均展现出了良好性能，取得了不错的脑功能连通性检测效果。同时，提出了一种自适应挖掘本真参考信号的方法，通过从多被试fMRI数据中挖掘可利用的先验信息，指导fMRI数据组被试水平上的ICA分析。实验结果表明，该方式在fMRI组分析中具有优越性能。总之，MOPCICA方法不仅克服和解决了现有CICA中存在的一些关键问题，同时也使得研究人员能够更加方便地使用CICA方法，使CICA方法能够更加广泛地应用于fMRI的各项研究，进而提高ICA在fMRI数据分析方面的能力。

1 材料和方法

1.1 材料

1.1.1模拟数据

本研究中的模拟数据由SimTB工具包(http://mialab.mrn.org/spftware)生成，该工具包由Vince D. Calhoun博士领导下的医学影像实验室发布，能够提供模拟信号的空间图谱与相应的时间过程或时间脑动力模式，方便不同算法进行性能对比评测。具体来说，采用SimTB生成10例模拟数据，每个被试的切片大小为100体素×100体素，共模拟了20个源信号(这些源信号直接从SimTB提供的信号模板中选取)，120个时间点，TR为2 s。该模拟数据中基准幅值设定为800，基准图谱如图1(a)所示，20个源信号的空间分布如图1(b)所示。为了模拟不同被试之间的空间差异，每个被试的源信号都通过正态偏差加入了少量的平移、旋转和伸缩变换。其中，平移和旋转变换分别服从正态分布N(0, 0.5)和N(0, 1)，而伸缩变换则服从正态分布N(1, 0.1)。

图1 模拟数据集。(a)空间基准图谱；(b)20个模拟源信号的空间分布；(c)所有模拟被试对应源信号s6的空间图谱和相应的时间过程(其中空间图谱中的横坐标和纵坐标均表示模拟的体素个数，而时间过程图中的横坐标和纵坐标则分别代表时间点数和去中心化后的幅度大小)Fig.1 Simulated dataset. (a) Spatial baseline map；(b)Spatial distribution of the 20 simulated sources；(c)Spatial map and corresponding time course of source 6 in all simulated subjects，where the abscissa and ordinate in the spatial map represent the number of simulated voxels, while the abscissa and ordinate in the time course subgraph represent the number of time points and the amplitude after decentralization, respectively

源信号s6除了混合特定的事件刺激之外，其时间过程还共享了特定区块事件刺激模式的影响，且区块事件相关的信号强度(定义为区块事件幅值与特定事件幅值之比)为2。区块事件刺激模式为(OFF, ON)×3，其中OFF或ON区块持续的时间长度均为40 s。通过将该刺激模式与典型的血氧动力学响应函数进行卷积，生成源信号s6的时间过程。剩下的19个源信号均未受到区块事件调制模式的影响，其活动模式仅由特定的血氧动力学波动信号(即通过特定事件与血氧动力学响应函数进行卷积生成)来描述。

所有模拟信号均受特定事件的影响，特定事件在每个TR内随机发生的概率为0.2。然而，对于受区块事件调制影响的源信号s6而言，特定事件的活动幅值略小，为0.4；对于不受区块事件调制影响的源信号而言，特定事件的幅值为1。对于所有的模拟信号，信号在空间域的变化强度服从正态分布N(0.03, 0.25)；与此同时，模拟信号中引入了加性噪声，其信噪比(signal-noise-ratio, SNR)服从0.2～1.2的均匀分布。最后，所有模拟被试对应源信号s6的空间图谱和时间过程如图1(c)所示。

1.1.2任务态数据

任务态数据共包括5名被试在视觉刺激任务下扫描获取的BOLD-fMRI数据。在进行数据采集之前，所有被试均被告知此次采集实验数据的目的，并获得了所有被试的数据采集许可，本次任务实验中的视觉刺激为具有放射状蓝黄团的旋转棋盘，其相应的视觉刺激为(OFF, ON)×3，每个区块(OFF或者ON状态)持续时间均为40 s；当刺激模式为ON时，刺激呈现屏中旋转棋盘；当刺激模式为OFF时，则要求被试眼睛注视刺激呈现屏的十字中心。该数据集中所有被试数据均在上海市磁共振重点实验室采集完成，采用磁场强度为3.0 T的西门子(Siemen)磁共振仪以及梯度回波平面成像技术进行覆盖全脑扫描获得，共扫描125个时间点。数据获取的其他参数为：TR=2 s，切片数=36片，扫描分辨率=64×64，片内分辨率=3.75 mm×3.75 mm，片厚=4 mm，切片间隔=1 mm。

1.1.3静息态fMRI数据

实验中使用的静息态fMRI数据集来自国际神经影像公共数据库(http://www.nitrc.org/projects/fcon_1000/)，包括23名认知能力正常、身体健康被试的静息态BOLD-fMRI扫描数据，被试年龄在20～40岁之间。所有被试数据通过磁场强度3.0T的飞利浦磁共振仪，以及单次激发敏感梯度回波平面成像技术覆盖全脑扫描获得，共扫描123个时间点。数据采集的其他参数为：TR=2.5 s，切片数=47，敏感加速因子=2.0，扫描分辨率=96×96，片内分辨率=2.67 mm×2.67 mm，片厚度=3 mm。

1.2 方法

1.2.1数据预处理

采用DPARSF(http://rfmri.org/DPARSF)软件对实验中的任务态和静息态fMRI数据进行预处理，包括时间层校正、头动校正、空间标准化和空间平滑，其中高斯平滑核设定为4 mm。特别地，在模拟数据和任务态数据实验中，仅实现了单被试水平上的ICA、包含时间先验信息的CICA (CICA with temporal reference, CICA-tR)[27]、包含空间先验信息的CICA (CICA with spatial reference, CICA-sR)[27]，以及本研究提出的MOPCICA方法；而在静息态数据中，则同时实现了组被试水平上的ICA、包含空间先验信息的基于牛顿迭代法的CICA(Newton′s method based CICA，CICA-nR)[28]和基于不动点迭代法的CICA(fixed point iteration based CICA，CICA-fR)[29]，以及MOPCICA方法和单被试水平上的ICA方法，并且在组分析中使用时空双回归方式来获取每个被试对应的时空信息。另外，采用最小描述长度方法[36]来对真实fMRI数据中的源信号个数进行估计。同时，为了在分析过程中获得更加稳定的输出结果，在ICA计算过程中使用ICASSO[37]来降低输出结果的随机性，其参数设置为运行100次，采用随机初始化和自助法相结合的采样方式。而对于其他方法则在计算过程取100次运行结果的平均值作为最终结果。最后，使用MRIcro软件，对各种方法输出结果中的空间网络和位置进行展示。

1.2.2MOPCICA方法优化

为了克服经典CICA方法中阈值参数难以选择的问题，同时考虑到时空先验信息的问题，本研究在多目标优化框架上建立CICA模型如下：

(1)

本研究采用G(v)=ln(cosh(v))，v是一个零均值单位方差的高斯随机变量，因此

(2)

在多目标优化问题式(1)中，每个解都对应一个最优分离向量wi，并且满足‖wi‖2=1。

在诸多求解多目标优化问题的方法中，线性加权求和方法是一种最为简单和有效的方法，只要权重参数严格为正且和为1，就能通过对目标函数的线性加权来进行单目标优化求解。因此，下面将采用这种方法，对所研究的多目标优化问题进行求解。只要在加权和函数中使用不同的权重值，就能找到不同解。为了避免优化过程被具有较大幅值的成本函数所控制，使用反正切函数来对式(1)中的目标函数J(wi)进行标准化，即

(3)

式中，ci是一个可以自动确定的常数，这样f1(wi)就和ε1(wi)、ε2(wi)一样，其可能的取值范围从0～1变化。

重新整理后，线性加权目标函数为

f(wi)=a1×f1(wi)+a2×ε1(wi)+a3×ε2(wi)

(4)

式中，ai(i=1,2,3)表示权重参数，且a1+a2+a3=1，通常是根据经验来决定的。

那么，优化f(wi)的迭代算法推导如下：

(5)

(6)

(7)

a3×E[Z-1×RTi]

(8)

式中，g(·)表示G(·)的导数，因此g(v)=tanh(v)，E(·)可以用所有样本的均值来进行估计，具体如下：

(9)

当目标函数式(4)的梯度被计算出来之后，就能建立最速上升的迭代公式，即

wi(k+1)=wi(k)+μ(k)×di(k)

(10)

对于P个参考信号RSp和RTp(p=1,2,…,P)，可以使用相同的方法逐个估计。最后，当获得独立成分Si后，其相应的时间过程可通过下列公式估计得到，即

Ti=Z-1×wi

(11)

1.2.3自适应先验信息挖掘

为了克服先验信息难以获取的问题，提出一种从组被试fMRI数据中获取空间先验信息的方法。假设在一组fMRI数据中有K个被试，且所有被试在经过标准化预处理后包含T个时间点和V个体素。

首先，进行单被试水平上的ICA分析，对于被试i，ICA可以表示为

Xi=Mi×Si(i=1,2,…,K)

(12)

式中：Xi表示T×V阶的观测数据；Mi表示T×Ni阶的混合矩阵；Si=(Si1,Si2,…,SiNi)'表示Ni×V阶的源信号矩阵，并且每行代表一个被试i在ICA分解过程中得到的独立成分，Sij(j=1,…,Ni)是大小为V×1的列向量。

接下来，记Sini(i=1,2,…,K)表示与感兴趣组独立成分对应被试i的第ni个独立成分，不同被试间的成分对应性可通过空间相关性进行度量[37]。

下面利用主成分分析(principal component analysis, PCA)，从由Sini(i=1,2,…,K)组成的K×v阶矩阵R来提取隐含参考信号，有

R=[S1n1,S2n2,…,SKnK]T

(13)

1.2.4算法验证

为了验证MOPCICA在fMRI数据分析中的有效性能，接下来分别通过模拟数据、任务态和静息态fMRI数据中的分析结果进行评估。其中，模拟数据和任务态fMRI数据主要用于评测MOPCICA的算法性能，而静息态fMRI数据则同时用于评估本真先验信息挖掘算法的有效性。虽然MOPCICA方法可以同时包含时间和空间先验信息，但是为了不同类型数据中前后分析的一致性，在分析中仅考虑了MOPCICA包含空间先验信息的情形。

图2 ICA、CICA-tR、CICA-sR和MOPCICA在每个模拟被试上对s6的源信号恢复性能评估。(a)空间相关系数；(b)空间AUC；(c)时间相关系数Fig.2 The performance evaluation of source signal recovery of s6 by ICA, CICA-tR, ICA-sR and MOPCICA in each simulated subject. (a)Spatial correlation coefficients；(b)Spatial AUC；(c)Temporal correlation coefficients

同时，在数据分析过程中，MOPCICA分别实现了空间权重参数从0.1～0.9、以步长为0.1取值时对应的9种情况，而在结果部分仅呈现了参数取值为0.5时的结果，主要是因为该参数在3种类型数据中均表现出良好的分析性能。另外，为了比较目的，在模拟数据和任务态fMRI数据中，同时呈现了ICA、CICA-tR和CICA-sR等方法的分析结果，并且三者均通过牛顿迭代算法实现。而在静息态fMRI数据中，则同时呈现了ICA、 CICA-nR和CICA-fR等方法的分析结果，并通过牛顿迭代算法实现ICA，还通过牛顿迭代算法实现CICA-nR和快速不动点算法实现CICA-fR。

最后，在3种类型的数据实验中，分别采用多种评价指标来评测比较上述不同方法的数据分析性能，而不同方法所对应相同指标之间的对比结果均通过置信水平为95%的双样本t检验统计获得。其中，在模拟数据实验中，采用空间AUC和时空相关系数来评测ICA、CICA-tR、CICA-sR和MOPCICA的源信号恢复性能。在任务态fMRI数据实验中，同样采用时空相关系数来评测ICA、CICA-tR、CICA-sR和MOPCICA的源信号恢复性能。同时，使用峭度和负熵来评测它们所得源信号的空间独立性。类似地，在静息态fMRI数据实验中，利用互信息、峭度和负熵来度量ICA、CICA-nR、CICA-fR和MOPCICA所得源信号的空间独立性。进一步，根据它们计算得到的组空间成分，通过时空双回归方式获得的组中每个被试对应空间成分之间的平均相关系数，评测其组分析性能。

2 结果

下面主要给出MOPCICA在模拟数据、任务态和静息态fMRI数据中的分析结果，以及与其他方法结果之间的比较情况。

2.1 模拟数据结果

为了评估上述方法的空间源信号恢复性能，分别计算了它们根据每个模拟被试数据获得的源信号成分与对应真实源信号之间的相关系数(即空间相关系数)，如图2(a)所示。可以发现，虽然与ICA相比无明显差异，但是MOPCICA在大部分被试上计算得到的相关系数都要高于CICA-tR和CICA-sR计算得到的相关系数，除了CICA-tR中被试2、3和10以及CICA-sR中被试1～3和10之外，这些有可能是由于噪声SNR不同而造成的结果。进一步，计算了它们在所有被试上的平均空间相关系数，对应ICA、CICA-tR、CICA-sR和MOPCICA分别为0.84±0.07、0.78±0.08、0.81±0.08和0.84±0.09。通过置信水平为95%的双样本t检验可知，MOPCICA显著高于CICA-tR和CICA-sR(P<0.05)，从而说明MOPCICA的空间源信号恢复能力要优于CICA-tR和CICA-sR的相关能力。

接着，ROC曲线中的能量分析也被用来评估这些方法的空间源信号恢复能力，表示为ROC曲线下的面积(AUC)。一般地，AUC越大，表示该方法对应的源信号恢复能力越强。图2(b)分别给出了ICA、CICA-tR、CICA-sR和MOPCICA在不同模拟被试上计算得到与源信号s6相应成分的AUCs。可以看出，MOPCICA在绝大部分被试上计算得到的AUCs都要高于ICA、CICA-tR和CICA-sR方法所得到的AUCs。它们在所有被试上的平均AUC分别为0.62±0.02、0.72±0.03、0.71±0.06和0.75±0.05，由双样本t检验可知，MOPCICA显著高于ICA、CICA-tR和CICA-sR(P<0.05)，表明MOPCICA具有更强的空间源信号恢复能力。

图3 ICA、CICA-tR、CICA-sR、MOPCICA和GLM计算得到的每个被试对应视觉相关网络的脑空间图谱及其时间过程(GLM除外)。其中，脑空间图谱通过z-score变换在阈值为2的条件下计算得到，对应的MNI坐标为(-10, -77, 6)；而每个时间过程小图中的横坐标和纵坐标分别代表时间点数和去中心化后的幅度大小，红线和蓝线则分别表示先验时间过程和实际计算得到的时间过程Fig.3 The spatial brain maps of the visual-related networks calculated by ICA, CICA-nR, CICA-fR, MOPCICA and GLM in each subject, as well as their time courses (except for GLM). The maps of brain networks are obtained with threshold |z| ≥ 2 after z-scored, and the corresponding MNI coordinates are (-10, -77, 6). The abscissa and ordinate of each time course subplot respectively represent the time points and the amplitude after decentralization, while the red line and the blue line respectively represent the prior time course and the time course obtained by actual calculation

此外，根据上述方法在每个模拟被试中计算得到的时间过程与相应真实时间过程之间的相关性，度量它们在时间域上的源信号恢复能力，如图2(c)所示。分别给出了ICA、CICA-tR、CICA-sR和MOPCICA在每个模拟被试上计算得到的时间过程与该被试对应真实时间过程之间的相关系数，即时间相关系数。从图中可以看出，除了少数被试之外，MOPCICA在大部分被试上的相关系数都要高于ICA、CICA-tR和CICA-sR计算得到的相关系数，这意味着MOPCICA在时间源信号恢复性能方面也要优于ICA、CICA-tR和CICA-sR的相关能力。而且，它们在所有被试上的平均时间相关系数分别为0.67±0.04、0.74±0.09、0.77±0.13和0.81±0.13，由双样本t检验可知，MOPCICA显著高于ICA、CICA-tR和CICA-sR(P<0.05)，从而进一步证明了MOPCICA在源信号检测方面的优越性能。

2.2 任务态数据结果

图3分别给出了ICA、CICA-tR、CICA-sR、MOPCICA和GLM检测得到每个被试对应视觉相关成分的脑空间图谱和时间过程(GLM除外), 以及该任务实验中的先验时间过程。考虑到大脑对外界刺激响应的延迟效应，这里的先验时间过程由实验中的Block刺激模式与延迟参数为2的伽马函数卷积得到。从图3中可以看到，MOPCICA检测得到的视觉区域要明显优于CICA-tR和ICA检测得到的视觉区域，但与CICA-sR检测得到的视觉区域相比没有明显差异，需要进一步量化分析，这意味着MOPCICA的空间源信号恢复能力要优于CICA-tR和ICA的空间源信号恢复能力。

图4 ICA、CICA-tR、CICA-sR和MOPCICA在每个被试上对视觉相关成分的源信号恢复性能评估。(a)空间相关系数；(b)峭度；(c)负熵；(d)时间相关系数Fig.4 The performance evaluation of source signal recovery of visual-related component by ICA, CICA-tR, CICA-sR and MOPCICA in each subject. (a)Spatial correlation coefficients；(b)Kurtosis；(c)Negentropy；(d)Temporal correlation coefficients

为了进一步评估上述方法的空间源信号恢复能力，分别计算了ICA、CICA-tR、CICA-sR和MOPCICA从每个被试fMRI数据中计算得到的视觉相关成分与相应由GLM方法计算得到的视觉相关成分之间的相关系数，如图4(a)所示。可以发现，除了少数被试差异不明显之外，MOPCICA在所有被试上计算得到的相关系数都要显著高于其他方法计算得到的相关系数。进一步计算了它们在所有被试上的平均相关系数，分别为0.73±0.11、0.28±0.05、0.72±0.03和0.84±0.04。由双样本t检验可知MOPCICA显著高于ICA、CICA-tR和CICA-sR(P<0.05)，意味着MOPCICA的空间源信号恢复能力要优于ICA、CICA-tR和CICA-sR的相应能力，与图3中的结论相一致。

同时，给出了上述方法计算得到每个被试对应视觉相关成分的峭度和负熵，以度量它们所得源信号的空间独立性，如图4(b)和4(c)所示。从图可以看出，除了ICA之外，MOPCICA在所有被试上计算获得的峭度和负熵都要高于CICA-tR和CICA-sR的相应情况。同时，计算了它们在所有被试上的平均峭度和负熵，分别为93.97±50.39、17.60±13.22、36.71±13.43和69.20±23.36以及0.030 2±0.004 9、0.003 7±0.002 1、0.018 4±0.004 5和0.031 2±0.007 7。根据双样本t检验可知，MOPCICA显著高于CICA-tR和CICA-sR，而平均峭度则低于ICA(P<0.05)，这说明MOPCICA估计所得源信号的独立性要优于CICA-tR和CICA-sR相应的独立性。

与模拟数据的结果分析类似，由计算得到的每个被试对应视觉相关成分的时间过程与先验时间过程之间的相关性被用来度量它们在时间域上的源信号恢复性能，图4(d)分别给出了ICA、CICA-tR、CICA-sR和MOPCICA在每个被试上计算得到的时间过程与先验时间过程之间的相关系数。可以看出，通过MOPCICA计算得到的相关系数要明显高于ICA、CICA-tR和CICA-sR计算得到的相关系数，而且它们在所有被试上的平均时间相关系数分别为0.66±0.07、0.90±0.01、0.85±0.05和0.91±0.03。

图5 ICA、CICA-nR、CICA-fR和MOPCICA计算得到的8种经典静息态脑网络对应的脑空间图谱及其相应的MNI坐标(脑网络图谱均通过z-score变换在阈值为2的条件下计算获得)Fig.5 The spatial maps of nine resting-state brain networks and their MNI coordinates including DMN, VIN, LVN, AUN, SMN, ECN, RWMN, LWMN and CCN calculated by ICA、CICA-nR、CICA-fR and MOPCICA (The maps of brain networks are obtained with threshold |z| ≥ 2 after z-scored)

由双样本t检验可知，MOPCICA显著高于ICA和CICA-sR(P<0.05)，说明MOPCICA在时间域上的分析性能要优于ICA、CICA-tR和CICA-sR等方法的相应性能。

2.3 静息态数据结果

图5分别给出了ICA、CICA-nR、CICA-fR和MOPCICA检测得到的空间独立成分对应的脑网络图谱及其相应的MNI坐标，包括默认网络(default mode network, DMN)、主视觉网络(visual network, VIN)、两侧视觉网络(bilateral visual network, BVN)、听觉网络(auditory network, AUN)、执行控制网络(executive control network, ECN)、突显网络(salience network, SAN)、右侧工作记忆网络(right working memory network, RWMN)和左侧工作记忆网络(left working memory network, LWMN)等8种经典静息态脑网络，其中所有脑网络图谱均通过z-score变换在阈值为2的条件下计算获得。由图5可知，上述方法均获得了明显的静息态脑功能网络，但无法看出明显差异，需要进一步定量分析。

图6 ICA、CICA-tR、CICA-sR和MOPCICA在组水平上对8种静息态脑网络的源信号恢复性能评估。(a)互信息；(b)峭度；(c)负熵；(d)组成分与通过双回归方式计算的组中每个被试对应独立成分之间的相关系数Fig.6 The performance evaluation of source signal recovery of eight resting-state brain networks by ICA, CICA-tR, ICA-sR and MOPCICA at group level. (a)Mutual information；(b)Kurtosis；(c)Negentropy；(d) The correlation coefficients between the group independent components and the corresponding independent components of each subject in the group obtained by dual-regression

图6给出了ICA、CICA-nR、CICA-fR和MOPCICA检测得到的8种静息态脑网络对应成分的空间独立性比较结果，其中独立性度量指标包括互信息(mutual information)、峭度(kurtosis)和负熵(negentropy)等3种情形分别如(a)～(c)所示，而(d)则分别给出了这4种方法计算得到的组空间成分与通过时空双回归方式获得的组中每个被试对应空间独立成分之间的相关系数。

由图6(a)可知，除了AUN对应的空间成分之外，MOPCICA计算得到的其他7种网络对应空间成分的互信息要低于其他方法计算得到的空间独立成分的互信息；而图6(b)和(c)的结果表明，除了少数几个脑网络之外，MOPCICA脑网络对应空间成分的峭度和负熵均要高于其他方法的值。进一步，计算了它们在所有被试上对应8种静息态脑网络的平均互信息、平均峭度和平均负熵，分别为0.44±0.04、0.44±0.04、0.44±0.04和0.43±0.04，8.68±3.28、8.85±3.32、8.78±3.25和9.15±3.64，以及5.09±5.10、4.23±4.30、4.00±4.50和4.88±4.81，根据t检验可知，除了ICA负熵之外，MOPCICA均显著优于其他3种方法(P<0.05)，从而说明MOPCICA计算得到的空间成分具有更好的独立性，更加符合ICA分析的前提假设条件，因而结果的可信性更高。

由图6(d)可知，相比于其他方法，MOPCICA计算得到的组成分与组中每个被试相应独立成分之间具有更高的相关性，而且它们在所有被试上对应8种静息态脑网络的平均相关系数为0.44±0.02、0.45±0.02、0.44±0.02和0.46±0.03。根据t检验可知，MOPCICA显著高于ICA、CICA-nR和CICA-fR(P<0.05)，说明MOPCICA计算得到的组成分更能代表组中被试的共性，具有更好的fMRI脑功能连通性检测性能。

3 讨论

在经典的CICA方法中，先验参考信号是通过不等式约束条件g(y)=ε(y,r)-ξ≤0的方式引入ICA算法的估计过程，y表示输出信号，r表示参考信号，ε(y,r)是距离测度，而ξ表示阈值参数，用于约束感兴趣信号是唯一满足不等式约束条件的输出信号[27-28, 30]。然而，在实际应用中，往往很难预先确定阈值参数ξ，因为独立成分是未知的，所以ξ的选择完全依赖于CICA的应用经验。不合适的ξ会导致两种可能的结果：当ξ大于可行域的上限时，输出信号可能是不感兴趣的信号；反之，当ξ小于可行域的下限时，可能无法输出任何的信号[34]。因此，往往需要额外的努力来决定一个合适的参数ξ。本研究采用了多目标优化策略，很好地规避了阈值参数ξ的选择问题；并且Shi等也采用了相同的策略，实验结果表明该策略的优越性能[38]。

在利用线性加权求和方法来求解式(1)中的多目标优化问题时，式(4)中的权重参数a1、a2和a3分别反映了和函数f(wi)中对应目标函数f1(wi)、ε1(wi)和ε2(wi)的重要性。加权求和方法的目标是在输出信号独立性和参考信号相似性之间寻求一个平衡，从而获得一个与参考信号最接近且独立性最大的源信号。根据Klamroth等提出的应用加权求和方法[39]，求解多目标优化问题的理论，只要权重参数满足它们严格为正且和为1的条件，那么每一组这样的权重参数都能获得Pareto最优集中的一个解。

在MOPCICA算法验证过程中，在模拟数据和任务态fMRI数据分析过程中，均只考虑一种源信号的情形。例如，在模拟数据中，仅考虑与特定区块事件刺激模式相关的源信号s6，主要是因为它在所有模拟被试上均有一个类似于区块事件刺激模式Block形状的时间过程[40]。特别地，这里的Block刺激模式为 (OFF, ON)×3，其中OFF或ON区块持续的时间长度均为40 s，这样可方便分析结果的评测；并且在计算过程中，分别利用该Block刺激模式和相应源信号s6的真实空间模板，作为时间和空间参考信号。同时，在分析过程中，模拟数据成分个数设置为21(包含20个源信号和1个背景信号)[41]。

任务态fMRI数据主要考虑与视觉任务相关的源信号成分，其对应的空间参考信号利用wfu_pickatlas软件(http://www.rad.wfubmc.edu/fmri)并根据布罗德曼模板产生，包括BA17和BA18两个脑区[42]。同时，利用先验的Block刺激模式作为时间参考信号。为了保证不同方法结果之间比较的一致性，结果分析中选取了GLM计算得到的视觉空间成分作为基准模板来进行比较验证，其中 GLM通过SPM(http://www.fil.ion.ucl.ac.uk/spm/)软件包实现，参数设定为使用FWE校正，p值取0.01，体素阈值个数设为10[3]。另外，在数据分析过程中，5名被试对应的源信号个数分别通过最小描述长度估计为50、52、54、54和44。

但是，这并不意味着MOPCICA仅适用于一种源信号的情形，因此在静息态fMRI数据分析过程中，考虑了多种源信号(对应于多种脑网络)的情形。由于静息态fMRI数据无法获得类似于模拟数据和任务态fMRI数据的时间先验信息，因此在数据分析过程中，首先通过本真先验信息挖掘算法获得几种经典脑网络(见图5)对应的空间先验信息，其中每个被试对应的源信号个数通过最小描述长度方法估计得到；然后利用ICA以及包含空间先验信息的CICA-sR和MOPCICA方法进行组水平数据分析，其中组水平源信号个数通过最小描述长度估计为26(为所有被试源信号估计个数的平均值)。上述方法在组分析过程中均采用时间级联方式，即将所有被试数据根据时间维度进行连接分析[43]，所以在实验数据预处理过程中，需要进行单被试和组水平两次主成分分析，且主成分个数分别为23和26。

在获取空间信息过程中，由于个体自身的差异性以及每个被试在采集数据过程中所受到的噪声影响不同，使得组中每个被试的成分不仅包含相同的有用信息，同时也包含不同的噪声信息[44]。因此，选择一种合适的方法，从这些成分中挖掘出最主要的信息十分关键。作为一种经典的多元数据处理方法，PCA不仅能够降低数据维度，提取出原数据中包含的主要信息，而且还能起到降噪的作用。所以，本研究采用PCA技术，从组被试的成分中挖掘出可利用的先验信息。因为第一主成分包含有原数据中最大量的信息，所以选取它为相应的参考信号。特别地，当特征值λ1=λ2时，选择它们对应主成分的平均值作为相应的参考信号。除此之外，其他方式(如平均等)也可被用来从组被试数据中求解隐含的参考信号，值得后续进一步研究。

4 结论

本研究在多目标优化的基础上，提出了一种同时包含时间和空间参考信号的MOPCICA方法，克服了经典CICA方法中阈值参数难以选择的问题，降低了CICA方法对参考信号的准确性要求，提高了它对错误信息的容错能力。同时，利用从fMRI自身挖掘的空间先验信息来指导组数据分析，使得MOPCICA获得的空间源信号成分能更好地代表组中被试的共性。