CN103440513B - 基于稀疏非负张量分解的大脑特定视觉认知状态判定方法 - Google Patents

Abstract

本发明公开了基于稀疏非负张量分解的大脑特定视觉认知状态判定方法,其创新点在于，采用SNTF算法，把fMRI认知数据当成一个大张量，从张量的层面上构建高阶非负张量模型，然后对fMRI认知数据每个维度上进行特征降维，从而得到维度较小的稀疏非负特征张量，最后结合支持向量机的特点，有效地实现对大脑特定视觉的认知状态的判定。本发明的方法利用稀疏非负张量分解进行降维和特征提取多方向、多角度的提取出原有数据中的潜在结构信息,由于l₁范数正则化和非负性的约束使得提取出的相关成分更稀疏，合乎大脑感知的直观体验，再结合支持向量机的特点，提高了分类判别的准确性。

Description

基于稀疏非负张量分解的大脑特定视觉认知状态判定方法

技术领域

本发明属于生物特征提取以及大脑认知状态判别领域，涉及与任务相关的功能核磁成像（functional magnetic resonance images，即fMRI）的预处理、稀疏非负张量分解（sparse nonnegative tensor factorization，即SNTF）的特征提取和支持向量机（support vector machine，即SVM）的判别分类，具体涉及一种基于稀疏非负张量分解的大脑特定视觉认知状态判定方法，可用于生物特征提取、维度约减、模式识别等领域。

背景技术

近年来，随着脑功能磁共振成像(functional Magnetic Resonance Imaging，fMRI)技术的发展，人们进行脑研究的能力大大增强，与此同时也产生了大量的数据，而仅仅利用计算机对实验数据进行简单的存储、查询和计算，根本无法挖掘出数据背后隐藏的大量信息，更不能满足研究的需求。因此，必须利用数据挖掘的方法，合理组织脑功能成像数据，把目前散乱的研究数据和结果加以有序化，整合心理学、语言学、脑科学等不同层次的研究数据，运用先进的信息学工具和方法帮助研究人员对实验数据进行分析、处理、整合和建模，从而发现新的规律，揭示脑高级认知活动的深层机制。

目前fMRI已被广泛的用来研究脑在执行某项认知功能时的激活，而已有的研究工作大多是比较不同任务的脑图或任务脑图与静息脑图之间的差异。但fMRI所采集的数据是一类维度高、数据量大，结构特别复杂的数据。因此，脑fMRI数据分析是一件困难而意义重大的工作，它直接决定了采用fMRI方法进行脑功能研究的成败。目前对fMRI数据分析的技术有多种，其中比较流行的有相关分析、方差分析、独立成分分析等。

相关分析是利用实验任务的先验知识来确定激活区和它们的强度，该方法的不足之处在于如果实验的fMRI信号变化的模型是未知的或不是常量时，例如，受试者是脑疾病患者或者正在进行复杂的学习任务，这种技术就不能探测出激活区。

方差分析需要数据满足4个假设：(1)随机变量的分布是已知的；(2)每次实验都有相同的方差和协方差；(3)不同因子的时间序列可预先被可靠地估计；(4)不同像素的信号是独立的，然而fMRI数据并不能保证满足所有的假设。

独立成分分析(Independent Component Analysis，即ICA)是一种采用高阶统计量的盲源分析方法(Blind Source Separation，即BSS)，目的在于从未知源信号的观测混合信号中分离(或抽取)相互统计独立的源信号。但该方法有一定的局限性：

(1)在单任务刺激中，ICA方法可以成功地分离出与任务相关的分量，主要是因为该任务与其它的心理活动所激活的脑区之间是空间独立的。但是，在多任务的情况下，我们不能假定不同任务所激活的脑区是空间独立的，因此，利用ICA方法进行功能定位必然受到局限。

(2)ICA方法分解的独立成分往往会有负成分的存在，即它不一定保证分解的成分是非负的。从计算的角度来看，分解的结果中可以存在负值，但负值元素在实际问题中往往缺失物理意义。

目前常用的分析方法的共同点是允许分解后的结果中出现负值，但就应用角度看负值是没有实际意义的，采用非负矩阵分解通过添加“矩阵中所有元素均为非负数”的限制条件,可以保证分解结果的可解释性,同时,它还具有实现简便和占用存储空间小的优点,从而更加贴近应用领域的特点，此外上述数据分析方法很少从机器学习的层面上进行大脑认知判定。文献“Padilla P，López M，Górriz J M，et al.NMF-SVM Based CAD ToolApplied to Functional Brain Images for the Diagnosis of Alzheimer's Disease[J].Medical Imaging，IEEE Transactions on，2012，31(2):207-216.”利用NMF（非负矩阵）对fMRI数据进行特征提取，同时结合支持向量机的特点获取了对痴呆病人的诊断。在进行特征提取时，它先将fMRI数据矢量化处理，然后利用非负矩阵进行特征分解。非负矩阵是盲源分析方法的一种，即它不需要人为设定参数或是根据任何经验模型对计算进行干预，最后的结果只与数据相关，同时它也克服了矩阵分解中负成分存在的局限性。正是由于具有非负性和稀疏性等特点，使得它对于特征提取和数据降维都具有很好的效果。在模式识别领域，它已成功运用到了人脸图像的识别当中。

但是，在现代科学技术和工程问题中常常会碰到大量的大规模数据处理问题，而非负矩阵分解在对大规模高维数据处理时，仅仅是将图像数据矢量化再根据特征值和特征向量进行特征提取，而忽略了处理对象往往是多维数据这个因素，从而造成了只在一个方向上投射而导致的降维效果不理想以及破坏了原始图像的结构和相关性，不能完全保持原始图像中的冗余和结构等问题。

作为向量和矩阵的推广，张量相当于一个高维的矩阵。由于张量模型可以处理指标集为多维的数据，所以相比矩阵模型其更能接近于实际问题的属性，因此能更好地描述实际问题。同矩阵的低秩逼近的思想一样，对于高阶张量模型人们也想借助于张量的分解模型来提取张量模型数据中潜在的内部结构信息。

发明内容

本发明针对传统矩阵分解模型和方法的不足，结合机器学习中支持向量机（SVM）的特点，提出了一种基于稀疏非负张量分解（SNTF）的大脑特定视觉认知状态判定方法。

为解决上述技术问题，本发明采用的技术方案：结合SNTF和SVM的特点，利用SNTF-SVM算法有效地进行大脑所处认知任务状态的鉴别。实现本发明目的的技术思路是，对从磁共振机器采集到的静息态数据先进行预处理分析，然后利用稀疏限制的非负张量(SNTF)方法分离出一个或多个与实验任务相关的独立成分，得到与脑部任务相关的激活区域，最后利用支持向量机(SVM)，获得对大脑认知状态的判定。

本发明的方法包括如下步骤：

步骤一：输入样本集，对fMRI认知数据进行预处理,并把样本集按照扫描时间序列以及行为实验类别重新组合，并分为训练集和测试集两部分，每部分都包含一定比例的认知数据。

对输入数据进行预处理，具体步骤如下：

a、对磁共振采集到的数据进行时间差异矫正，时间差异矫正就是矫正1个volume中层与层之间的采集时间的差异，进而保证各层之间都是从相同的时间得到的；

b、由于脑功能成像实验持续时间比较长，被试的呼吸、血流脉动等生理因素造成的头部运动在所难免，所以对时间矫正过的数据进行头动对齐，也就是将一个序列的每一帧图像都和这个序列的第一帧图像进行配准，配准到同一坐标系下，以矫正头动，对齐后再对头动矫正后的数据进行人工检查，如果平移和旋转分别超过1毫米和1度就排除掉、不予分析；

c、头动对齐后的数据空间标准化，试验中存在多个被试，被试和被试之间的脑部大小形状存在一定的差异，为了后续的统计分析，须进行大脑形状的归一化，将被试的大脑配准归一化到标准的大脑模板上；

d、采用6毫米半高全宽(FWHM，Full Width at Half Maximum)的高斯核函数平滑处理标准化后的数据，以达到降低噪声、提高信噪比及消除不同被试脑结构之间的细微差别的目的；

e、将预处理后的图像数据按照扫描时间序列以及行为实验类别进行重新组合，并将其与包括灰质、白质、脑脊液的统计模版进行匹配，去除白质区域，构成样本数据集；

f、将样本集按照扫描时间序列以及行为实验类别重新组合，并分为训练集和测试集两部分，每部分都包含一定比例的认知数据。

步骤二：计算训练样本集的特征分解，求出特征向量以及特征变换矩阵，将样本投射到特征子空间，得到特征张量，之后直接矢量化，得到训练集的特征向量集。

计算训练样本的特征分解，求出特征张量以及最优稀疏特征变换矩阵的具体步骤如下：

a、选取训练样本，构建张量空间：

稀疏编码的基本思想可以概括为将一个特征向量表示为一组从预定义的字典中选择的少量基的线性组合，因此导致了稀疏性的概念。换句话说，在给定一个基向量集合的基础上，稀疏编码提供了给定信号数据的低维近似表达。从模型拟合的角度来看，稀疏编码可以通过在评价拟合函数时去除不重要的变量的方法来避免过拟合的问题。等价地说，稀疏编码在最小化数据拟合的损失函数与模型系数的稀疏正则化之间选择一个平衡的解决。具体可参见下式:

公式(2.1)

其中，X∈R^D×N是输入信号的列矩阵，对应的稀疏系数矩阵为S∈R^P×N,字典基为B∈R^D×P。

在公式(2.1)中，最流行的损失函数选择为平方损失函数,正则化选项为l₁范数，因此公式(2.1)成为如下的l₁正则化的最小平方优化问题： $\arg \min \frac{1}{2}{\left|\right|X-\mathrm{BS}\left|\right|}_F^2+\mathrm{\β}{\left|\right|S\left|\right|}_1$

$\mathrm{subject\; to}{\mathrm{\Σ}}_i{B}_{i,j}^2\≤c,\∀j=1,...,N---\left(2.2\right)$

由于公式(2.2)中，当字典基B一定时，第一项是凸的，具有全局最小值。虽然第二项的l₁范数也是凸的，但是在原点不可导，因此需要一些特殊的策略来处理。如果不保持基B一定，那么学习到的B对不同的数据域有更好的适应性.尽管公式(2.2)中的问题并不是对B和S同时凸的，但是，我们仍可以按下面的策略获得求解。即交替地保持其中一个量一定(如B)，求解公式(2.2)所示的目标函数来更新另一个量(如S)，如此不断迭代直到收敛。经过实验证明，数据经过l1范数正则化处理主要得到了两方面的效果:一方面获得了稀疏性；另一方面更加突出的特征会被保留下来。

对于选取的训练样本，这里定义M个样本组成的集合为{χ₁,χ₂,…,χ_M}，每个样本对应的张量空间为其中I_n(1≤n≤N)是张量的n-mode（模式）维度。稀疏非负张量分解目标是寻找多维稀疏变换矩阵集U⁽ⁿ⁾代表n-mode（模式）方向上的非负特征投射矩阵，将原始张量空间同时从N个维度上投影到稀疏非负张量子空间

$R^{J_1}\⊗R^{J_2}\⊗...\⊗R^{J_N}(J_n\≤I_n,n=1,...,N):$

其中

其中，χ_m×_nU⁽ⁿ⁾表示张量χ_m与矩阵在n-mode（模式）方向上的乘，其中X_m第n模方向的维度必须与矩阵的行相等；

稀疏非负张量分解（SNTF）的目标函数是:

$\underset{U^1,...,U^N}{\min }\frac{1}{2}{\left|\right|\mathrm{\χ}-G{\×}_1{U}^{{\left(1\right)}^T}{\×}_2{U}^{{\left(2\right)}^T}...{\×}_N{U}^{{\left(N\right)}^T}\left|\right|}^2+\underset{1\≤n\≤N}{\mathrm{\Σ}}{\mathrm{\λ}}_n\left|U^{\left(n\right)}\right|$

s.t.U⁽ⁿ⁾≥0,1≤n≤N(保证分解的投影变换矩阵非负)

其中 $\mathrm{\χ}\∈R^{I_1\×I_2\×...\×I_N},U^{\left(n\right)}\∈R^{I_n\×J_n},n=1,...,N,$ 特征张量G是同一张量，即只含有0、1两个元素，主对角元素上都是1，其余的为0，λ_n表示稀疏系数，主要是用于限制投影变换矩阵U⁽ⁿ⁾稀疏。

这样该问题就转换成了一个最优化问题，为了找到全局上的最小值，可以将上面的最优化问题重新改写为：

$\underset{u_i}{\min }:\frac{1}{2}{\left|\right|B_{\left(n\right)}^T{u}_i-x_i\left|\right|}^2+{\mathrm{\λ}}_n\left|u_i\right|$

s.t.u_i≥0(代表U_n上的每一列)

其中 $B_{\left(n\right)}={(G{\×}_1{U}^{{\left(1\right)}^T}{\×}_2{U}^{{\left(2\right)}^T}...{\×}_N{U}^{{\left(N\right)}^T})}_{\left(n\right)},{\mathrm{\χ}}_{\left(n\right)}={[x_1,x_2,...,x_n]}^T,$ $U_n={[u_1,u_2,...,u_n]}^T.$

稀疏系数λ_n是用于控制变换矩阵的稀疏度，即l₁范数。通常所说的稀疏度是与非零矩阵中零元的个数所占的比重有关，而与元素的符号及顺序无关。

b、初始化：

对非负特征投射矩阵集进行随机初始化，满足非负要求。

c、求解特征变换矩阵：

在迭代优化的过程中，需要对作进一步优化。为了保证优化的始终是稀疏的、非负的，每次迭代优化时都用替代u_i：

$u_i^{*}=\frac{t-{\mathrm{\λ}}_n}{b_j{b}_j^T}(t>{\mathrm{\λ}}_n),$ 其中 $t=b_j(B_{\left(n\right)}^T{u}_i-b_j^T{a}_i)-b_j{u}_{\mathrm{ij}}$

如果t<λ_n时，则

这里即u_i代表的是U_n的每一列；当完成对的一次更新后，通过训练样本与更新的投射矩阵对应模方向上相乘，得到

这样就将原始张量空间投射到更新之后的张量子空间，通过判断是否满足这个终止条件来决定是否还需要迭代以优化 $\{U^{\left(n\right)}\&Element;R^{I_n\&times;J_n},n=1,...,N\},$ 其中η是用户自定义的比较小的阈值；如果不满足则重复迭代优化步骤继续寻优，直到满足条件退出优化迭代过程。

这里的迭代更新算法叫做列元坐标递减算法，即每次迭代时都是按U_n的列进行更新，相比传统迭代算法（每次更新一个元素），列元迭代优化算法效率会更高，具体步骤如下：

（ⅰ）随机初始化非负投影矩阵并设定相应的稀疏度λ_n。

（ⅱ）对于k=1:K(迭代次数)

对于j=1:J_n(张量阶数)

设定 $M=B_{\left(n\right)}{B}_{\left(n\right)}^T,N={\mathrm{\&chi;}}_{\left(n\right)}{B}_{\left(n\right)}^T-{\mathrm{\&lambda;}}_n,D=[M\left(\mathrm{1,1}\right)M\left(\mathrm{2,2}\right)...M(J_n,J_n)];$

令M(j,j)=0

if not convergent(不收敛)

$U_n(:,j)=\max (0,\frac{N(:,j)-U_{n}M}{D\left(j\right)})$

end(收敛)

这样不断通过迭代优化U_n(:,j)，进而求取特征投影矩阵U_n。这里的收敛条件是（η是用户自定义的比较小的阈值），如果满足收敛条件或者满足k=K，则跳出循环，得到最优非负特征投射矩阵集进入到下一步。

d、计算训练样本的特征张量集：

训练样本χ_m(m=1,…,M)向稀疏非负张量投影算子分解（SNTF）之后的非负特征张量子空间上投影，得到训练样本的特征张量集为 $\{G_m={\mathrm{\&chi;}}_m{\&times;}_1{U}^{{\left(1\right)}^T}{\&times;}_2{U}^{{\left(2\right)}^T}...{\&times;}_N{U}^{{\left(N\right)}^T},m=1,...,M\},G_m\&Element;R^{J_1\&times;J_2\&times;...\&times;J_N}.$

e、将特征张量G进行矢量化处理，得到训练集的特征向量集。

步骤三：将上述SNTF降维后的特征向量集作为SVM的输入，提取支持向量，构建最优分类器。

利用SVM构建最优超平面（分类器）的具体步骤如下：

a、对训练集中的各个样本的特征向量分类标记，根据所属的类别分别标记为C₁,C₂,C₃,…；

b、对已经标记好的特征向量进行SVM训练，得到多个最优超平面分类器。即以标记好的特征向量集作为支持向量机的输入，求取最优SVM投射方向w及偏移常量b。

最优SVM投射方向W通过计算如下优化问题求得：

$\underset{w,\mathrm{\&zeta;},b\&GreaterEqual;0}{\min }<w,w>+C\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&zeta;}}_i$

其中C为支持向量机中的惩罚因子，ζ=[ζ₁,ζ₂,…ζ_M]为支持向量机中的松弛因子，y_i∈[1,-1]为训练样本类别标记，通常我们定义1标记一类认知任务，-1标记另一类认知任务。这里由于有多个认知任务，因此得构建多个最优分类器。

步骤四：特征分类，利用构建好的分类器对测试样本进行分类判别。

利用训练得到的分类器对测试样本分类的具体步骤如下：

a、将测试样本分别输入到多个最优分类器中，得到相应的分类结果，累计计算各类的得分（分类准确的个数）；

b、统计步骤A中所得的分类结果，选择得分最高的类作为测试数据的类别，从而获得各个样本的分类。

本发明的有益效果：

（1）利用稀疏非负张量分解进行降维和特征提取，克服了传统矩阵分解因矢量化处理而破坏了原始图像的结构和相关性，保留了fMRI成像数据的空间结构信息；

（2）提取实验数据中与认知活动相关成分更稀疏。稀疏限制的非负张量分解算法由于其特有的稀疏性使得提取出的相关成分更稀疏，即用更少量的数据去描述大量数据，从降维的角度来看，它有效的实现了维度约减；

（3）对提取的相关成分的解释显得更为合理。由于稀疏限制的非负张量分解算法的非负性，使得提取出的相关成分合乎大脑感知的直观体验，并具有明确的物理含义；

（4）结合支持向量机（SVM）分类判别和有限样本相适应的特点，克服了如高维度、小样本、局部极小值等问题，提高了分类判别的准确性。

附图说明

图1为本发明对fMRI数据预处理分析示意图；

图2为本发明详细流程图，主要包括稀疏非负张量分解和支持向量机的实施示意图。

具体实施方式

以下结合附图，详细描述一个完整的基于稀疏张量分解的大脑特定视觉认知判定方法的具体流程。

一、参照图1，本发明预处理部分包括：对磁共振扫描后采集到的数据进行数据预处理（基于SPM8），然后将它们按照时间序列与实验类别进行重新组合，最后与标准模板进行匹配，得到试验数据的灰质、白质、脑脊液的部分。具体步骤如下：

步骤1：对磁共振采集到的数据进行时间差异矫正，时间差异矫正就是矫正1个volume中层与层之间的采集时间的差异，进而保证各层之间都是从相同的时间得到的。

步骤2：由于脑功能成像实验持续时间比较长，被试的呼吸、血流脉动等生理因素造成的头部运动在所难免，所以对时间矫正过的数据进行头动对齐，也就是将一个序列的每一帧图像都和这个序列的第一帧图像进行配准，配准到同一坐标系下，以矫正头动；对齐后再对头动矫正后的数据进行人工检查，如果平移和旋转分别超过1毫米和1度就排除掉、不予分析。

步骤3：头动对齐后的数据空间标准化，试验中存在多个被试，被试和被试之间的脑部大小形状存在一定的差异，为了后续的统计分析，须进行大脑形状的归一化，将被试的大脑配准归一化到标准的大脑模板上。

a）对头动对齐后的数据利用12个参数的仿射变换平均图像配准到MNI标准模板；

b）重切为3x3x3立方毫米的体素，并将MNI(Montreal Neurological Institute)坐标变换到Talairach坐标系；

步骤4：采用6毫米半高全宽(FWHM，Full Width at Half Maximum)的高斯核函数平滑处理标准化后的数据，以达到降低噪声、提高信噪比及消除不同被试脑结构之间的细微差别的目的。

步骤5：将预处理后的图像数据按照扫描时间序列以及行为实验类别进行重新组合，并将其与包括灰质、白质、脑脊液的统计模版进行匹配，去除白质区域，构成样本数据集；

步骤6：将样本集按照扫描时间序列以及行为实验类别重新组合，并分为训练集和测试集两部分，每部分都包含一定比例的认知数据。

二、参照图2，本发明稀疏非负张量分解和支持向量机部分包括：先将实验数据构建到张量空间，然后利用稀疏限制的非负张量(SNTF)方法分离出一个或多个与实验任务相关的独立成分，得到与脑部任务相关的激活区域，最后利用支持向量机(SVM)，获得对大脑认知状态的判定。具体步骤如下：

步骤1：选取训练样本，构建张量空间。

对于选取的训练样本，定义M个样本组成的集合为{χ₁,χ₂,…,χ_M}，每个样本对应的张量空间为其中I_n(1≤n≤N)是张量的n-mode（模式）维度；稀疏非负张量分解目标是寻找多维稀疏变换矩阵集U⁽ⁿ⁾代表n-mode（模式）方向上的非负特征投射矩阵，将原始张量空间同时从N个维度上投影到稀疏非负张量子空间

$R^{J_1}\&CircleTimes;R^{J_2}\&CircleTimes;...\&CircleTimes;R^{J_N}(J_n\&le;I_n,n=1,...,N):$

其中

其中，χ_m×_nU⁽ⁿ⁾表示张量X_m与矩阵在n-mode（模式）方向上的乘，其中χ_m第n模方向的维度必须与矩阵的行相等。

稀疏非负张量分解（SNTF）的目标函数是:

$\underset{U^1,...,U^N}{\min }\frac{1}{2}{\left|\right|\mathrm{\&chi;}-G{\&times;}_1{U}^{{\left(1\right)}^T}{\&times;}_2{U}^{{\left(2\right)}^T}...{\&times;}_N{U}^{{\left(N\right)}^T}\left|\right|}^2+\underset{1\&le;n\&le;N}{\mathrm{\&Sigma;}}{\mathrm{\&lambda;}}_n\left|U^{\left(n\right)}\right|$

s.t.U⁽ⁿ⁾≥0,1≤n≤N(保证分解的投影变换矩阵非负)

其中 $\mathrm{\&chi;}\&Element;R^{I_1\&times;I_2\&times;...\&times;I_N},U^{\left(n\right)}\&Element;R^{I_n\&times;J_n},n=1,...,N,$ 特征张量G是同一张量，即只含有0、1两个元素，主对角元素上都是1，其余的为0，λ_n表示稀疏系数，主要是用于限制投影变换矩阵U⁽ⁿ⁾稀疏。

这样该问题就转换成了一个最优化问题，为了找到全局上的最小值，可以将上面的最优化问题重新改写为：

$\underset{u_i}{\min }:\frac{1}{2}{\left|\right|B_{\left(n\right)}^T{u}_i-x_i\left|\right|}^2+{\mathrm{\&lambda;}}_n\left|u_i\right|$

s.t.u_i≥0(代表U_n上的每一列)

其中 $B_{\left(n\right)}={(G{\&times;}_1{U}^{{\left(1\right)}^T}{\&times;}_2{U}^{{\left(2\right)}^T}...{\&times;}_N{U}^{{\left(N\right)}^T})}_{\left(n\right)},{\mathrm{\&chi;}}_{\left(n\right)}={[x_1,x_2,...,x_n]}^T,$ $U_n={[u_1,u_2,...,u_n]}^T.$

稀疏系数λ_n是用于控制变换矩阵的稀疏度，通常所说的稀疏度是与非零矩阵中零元的个数所占的比重有关，而与元素的符号及顺序无关。对于非零向量x，记||x||₀为向量的非零个数，称为0范数，则稀疏值可简单定义如下：

${\mathrm{sp}}_0\left(x\right)=\frac{n-{\left|\right|x\left|\right|}_0}{n-1}\&Element;\left[\mathrm{0,1}\right].$

我们称为向量的0-稀疏值。这种定义比较直观，对于n维向量来说，它的值域是区间[0,1]的n个等分点（包含0和1）。注意到，如果向量的零元或者小元素受到微小的扰动(如零元变成的ε，向量的0-稀疏值会有较大的偏离(至少为)。为解决这种不确定性，通常考虑寻找一个逼近sp₀(x)的连续函数作为稀疏值的定义。

现有稀疏值的定义是

$\mathrm{\&lambda;}\left(x\right)={\mathrm{sp}}_1\left(x\right)=\frac{\sqrt{n}-1\left(x\right)}{\sqrt{n}-1}\&Element;\left[\mathrm{0,1}\right].$

这个通常称为向量的1-稀疏值，其中由于它有效的克服了0-稀疏值出现的问题，而且现在用的比较广泛，因此本文中我们选用向量的1-稀疏值。

步骤2：初始化。

对非负特征投射矩阵集进行随机初始化，满足非负要求。

步骤3：求解特征变换矩阵。

在迭代优化的过程中，需要对作进一步优化。为了保证优化的始终是稀疏的、非负的，每次迭代优化时都用替代u_i：

$u_i^{*}=\frac{t-{\mathrm{\&lambda;}}_n}{b_j{b}_j^T}(t>{\mathrm{\&lambda;}}_n),$ 其中 $t=b_j(B_{\left(n\right)}^T{u}_i-b_j^T{a}_i)-b_j{u}_{\mathrm{ij}}$

如果t<λ_n时，则

这里即u_i代表的是U_n的每一列；当完成对的一次更新后，通过训练样本与更新的投射矩阵对应模方向上相乘，得到

这样就将原始张量空间投射到更新之后的张量子空间，通过判断是否满足这个终止条件来决定是否还需要迭代以优化 $\{U^{\left(n\right)}\&Element;R^{I_n\&times;J_n},n=1,...,N\},$ 其中η是用户自定义的比较小的阈值；如果不满足则重复迭代优化步骤继续寻优，直到满足条件退出优化迭代过程。

这里的迭代更新算法叫做列元坐标递减算法，即每次迭代时都是按U_n的列进行更新，相比传统迭代算法（每次更新一个元素），列元迭代优化算法效率会更高。

其具体步骤如下：

（a）随机初始化非负投影矩阵并设定相应的稀疏度λ_n。

（b）对于k=1:K(迭代次数)

对于j=1:J_n(张量阶数)

设定 $M=B_{\left(n\right)}{B}_{\left(n\right)}^T,N={\mathrm{\&chi;}}_{\left(n\right)}{B}_{\left(n\right)}^T-{\mathrm{\&lambda;}}_n,D=[M\left(\mathrm{1,1}\right)M\left(\mathrm{2,2}\right)...M(J_n,J_n)];$

令M(j,j)=0

if not convergent(不收敛)

$U_n(:,j)=\max (0,\frac{N(:,j)-U_{n}M}{D\left(j\right)})$

end(收敛)

这样不断通过迭代优化U_n(:,j)，进而求取特征投影矩阵U_n。这里的收敛条件是（η是用户自定义的比较小的阈值），如果满足收敛条件或者满足k=K，则跳出循环，得到最优非负特征投射矩阵集 $\{U^{\left(n\right)}\&Element;R^{I_n\&times;P_n},n=1,...,N\},$ 进入到下一步。

步骤4：计算训练样本的特征张量集。

训练样本χ_m(m=1,…,M)向稀疏非负张量投影算子分解（SNTF）之后的非负特征张量子空间上投影，得到训练样本的特征张量集为 $\{G_m={\mathrm{\&chi;}}_m{\&times;}_1{U}^{{\left(1\right)}^T}{\&times;}_2{U}^{{\left(2\right)}^T}...{\&times;}_N{U}^{{\left(N\right)}^T},m=1,...,M\},G_m\&Element;R^{J_1\&times;J_2\&times;...\&times;J_N}.$

步骤5：将特征张量G进行矢量化处理，得到训练集的特征向量集。

步骤6：对训练集中的各个样本的特征向量分类标记，根据所属的类别分别标记为C₁,C₂,C₃,…。

步骤7：对已经标记好的特征向量进行SVM训练，得到多个最优超平面分类器，即以标记好的特征向量集作为支持向量机的输入，求取最优SVM投射方向W及偏移常量b。

最优SVM投射方向w通过计算如下优化问题求得：

$\underset{w,\mathrm{\&zeta;},b\&GreaterEqual;0}{\min }<w,w>+C\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&zeta;}}_i$

其中C为类似与支持向量机中的惩罚因子，ζ=[ζ₁,ζ₂,，ζ_M]等同于支持向量机中的松弛因子，y_i∈[1,-1]为训练样本类别标记，通常我们定义1标记一类认知任务，-1标记另一类认知任务。这里由于有多个认知任务，因此我们得构建多个最优分类器。

步骤8：根据步骤7，利用训练得到的分类器对测试样本分类的具体步骤如下：

a)将测试样本分别输入到多个最优分类器中，得到相应的分类结果，累计计算各类的得分（分类准确的个数）。

b)统计步骤a)中所得的分类结果，选择得分最高的类作为测试数据的类别，从而获得各个样本的分类。

本发明所述的基于稀疏张量分解的大脑特定视觉认知状态判定方法并不仅限于说明书和实施方式中的描述。凡在本发明的精神和原则之内，所做的任何修改、同等替换、改进等，均包含在本发明的权利要求范围之内。

Claims (5)

1.基于稀疏非负张量分解的大脑特定视觉认知状态判定方法，其特征在于：采用基于张量模型的稀疏非负张量分解(SNTF)的特征降维和支持向量机(SVM)的分类判定，包括如下步骤：

1)输入样本集，对fMRI认知数据进行预处理,并把样本集按照扫描时间序列以及行为实验类别重新组合，并分为训练集和测试集两部分，每部分都包含一定比例的认知数据；

2)计算训练样本集的特征分解，求出特征向量以及特征变换矩阵，将样本投射到特征子空间，得到特征张量之后直接矢量化，得到训练集的特征向量集；

3)将稀疏非负张量分解特征降维后的特征向量集作为SVM的输入，提取支持向量，构建最优分类器；

4)特征分类，利用构建好的分类器对测试样本进行分类判别。

2.根据权利要求1所述的基于稀疏非负张量分解的大脑特定视觉认知状态判定方法，其特征在于：所述的步骤1)中，对fMRI认知数据的预处理步骤如下：

a)利用SPM软件对fMRI认知数据进行预处理，包括时间校正、空间校正、配准图像、分割、图像标准化、平滑处理步骤；

b)将预处理后的图像数据按照扫描时间序列以及行为实验类别进行重新组合，并将其与包括灰质、白质、脑脊液的统计模版进行匹配，构成样本数据集。

3.根据权利要求1所述的基于稀疏非负张量分解的大脑特定视觉认知状态判定方法，其特征在于：所述的步骤2)中，计算特征张量及特征变换矩阵的步骤如下：

a)选取训练样本，构建张量空间：对于选取的训练样本，定义M个样本组成的集合为{χ₁,χ₂,…,χ_M}，每个样本对应的张量空间为其中I_n(1≤n≤N)是张量的n-mode(模式)维度；寻找多维稀疏变换矩阵集U⁽ⁿ⁾代表n-mode(模式)方向上的非负特征投射矩阵；

b)初始化：对非负特征投射矩阵集进行随机初始化，满足非负要求；

c)求解特征变换矩阵：在迭代优化的过程中，对作进一步优化，保证优化的始终是稀疏的、非负的；

d)计算训练样本的特征张量集：训练样本χ_m(m＝1,…,M)向稀疏非负张量投影算子分解(SNTF)之后的非负特征张量子空间上投影，得到训练样本的特征张量集为

e)将特征张量G进行矢量化处理，得到训练集的特征向量集。

4.根据权利要求1所述的基于稀疏非负张量分解的大脑特定视觉认知状态判定方法，其特征在于：所述的步骤3)中，利用SVM构建最优超平面分类器的步骤如下：

a)对训练集中的各个样本的特征向量分类标记，根据所属的类别分别标记为C₁,C₂,C₃,…；

b)对已经标记好的特征向量进行SVM训练，得到多个最优超平面分类器，即以标记好的特征向量集作为支持向量机的输入，求取最优SVM投射方向及偏移常量b；

5.根据权利要求1所述的基于稀疏非负张量分解的大脑特定视觉认知状态判定方法，其特征在于：所述的步骤4)中，利用训练得到的最优分类器对测试样本分类的步骤如下：

a)将测试样本分别输入到多个最优分类器中，得到相应的分类结果，累计计算各类的得分；

b)统计步骤a)中所得的分类结果，选择得分最高的类作为测试数据的类别，从而获得各个样本的分类。