CN103116764A

CN103116764A - 一种基于多线性主元分析的大脑认知状态判定方法

Info

Publication number: CN103116764A
Application number: CN2013100665640A
Authority: CN
Inventors: 李军; 甘云; 徐鑫秀; 王洪勇; 李明欣; 袁森; 曹凯; 梁继民; 秦伟
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2013-03-02
Filing date: 2013-03-02
Publication date: 2013-05-22
Anticipated expiration: 2033-03-02
Also published as: CN103116764B

Abstract

本发明公开了一种基于多线性主元分析的大脑认知状态判定方法，包括如下步骤：1)输入样本集，对输入数据进行预处理；2)计算训练样本集的特征分解，求出最优特征变换矩阵，将训练样本投射到张量特征子空间，得到训练集的特征张量集；3)将降维后的低维特征张量数据矢量化作为LDA的输入，求出LDA最优投射矩阵，并将矢量化的低维特征张量数据投射到LDA特征子空间进一步提取训练集的判别特征向量；4)特征分类，将训练图像和测试图像的投射所得判别特征向量进行特征匹配，进而分类。本发明利用多线性主元分析直接对多阶张量数据进行降维和特征提取，克服了传统PCA单纯的进行降维而破坏了原始图像数据的结构和相关性，不能完全保持原始图像中的冗余和结构的不足，保留了fMRI成像数据的空间结构信息。

Description

一种基于多线性主元分析的大脑认知状态判定方法

技术领域

本发明属于生物特征提取以及大脑认知状态判定分类领域，涉及脑血氧水平的功能核磁成像的预处理、多线性主成分的提取和线性判别分类，是基于多线性主元分析的特征提取算法和线性判别分类算法。

背景技术

认知神经科学是在认知科学和神经科学的基础上发展起来的一门学科，其核心分支学科是认知心理学和人工智能。其中，认知心理学结合信息加工理论解释人的认知过程，并采用科学的方法对人的感知觉、注意、记忆、语言等认知过程展开实验研究。与此同时，人工智能的发展引发了学界对人脑认知过程更多的关注与研究，一些新的理论成果相继出现，大大丰富了认知科学的内涵。认知神经科学旨在从分子水平、细胞水平对大脑内与学习和记忆相关的一系列神经活动进行研究。

当前，基于血氧水平依赖（Blood-oxygen-level dependent，简称BOLD）的fMRI（functional magnetic resonance images，即功能磁共振成像技术）已经成为认知神经科学中研究大脑功能使用最为广泛的研究手段。fMRI所采集的实验数据是一类维数高、数据量巨大、噪声非常强、结构特别复杂的数据。因此，脑fMRI数据分析是一件困难而意义重大的工作，它直接决定了采用fMRI方法进行脑功能研究的成败。

现有的fMRI数据分析方法主要有两类：模型驱动的参数方法和数据驱动的非参数方法。模型驱动的参数方法中的模型假设需要对脑的结构、功能和处理机制有比较深刻的认识，但这些认识本身就是脑功能研究的目标。所以，模型驱动的参数方法只是一种对假设模型的验证方法，而模型中关于生理、神经、功能等机理的假设还有待进一步研究、验证和解释。数据驱动的非参数方法主要是直接选用统计学、机器学习等学科中已有的多元统计方法对fMRI数据进行处理，只是非常有限地利用了fMRI和脑认知领域的知识。因此，非参数方法的处理结果往往不是很理想，并且合理、准确地解释这些结果也比较困难。鉴于上述情况，脑fMRI研究迫切需要能更有效、更准确地提取实验数据中与认知活动相关的信息数据分析的新方法。在模式识别里的张量对象通常都是指高维的张量空间，直接应用于这种空间的识别算法往往要受到维度的限制，即所谓的维数灾难。为了解决维数灾难，出现了很多降维算法。

主元分析法(Principalcomponentanalysis，简称PCA)是一种经典的用于降维的无监督的线性技术。PCA的核心思想是减少包含大量无关参数的数据集的维度，同时尽可能多的保留出现在原始数据集中的参数，这能够通过转化为一个新的参数集（即主成分）来实现，这些参数是不相关的、有序的，以至于第一个很少能保留原始数据参数的大部分。PCA理想化的应用于张量对象要求它们重新变形为高维的向量（向量化），这明显地造成了计算机内存的消耗。除了这些实施的问题外，很容易理解的就是向量化会重新破坏它的本来结构和原始数据间的相关性，消除了冗余，或者是对原始数据集有了更高的顺序依赖性，并且可能失去了更简洁、更有用的能从原始形式中获得的表现形式。向量化作为PCA的预处理忽略了张量对象本来就是多维的对象。因此，直接作用在张量对象上而不是在它的向量化上的降低维度算法更加可取。

线性判别分析(Linear Discriminant Analysis,LDA)，也叫做Fisher线性判别(Fisher Linear Discriminant,FLD)，是模式识别的又一经典算法，它是在1996年由Belhumeur引入模式识别和人工智能领域的。线性鉴别分析的基本思想是将高维的模式样本投影到最佳鉴别矢量空间，以达到抽取分类信息和压缩特征空间维数的效果，投影后保证模式样本在新的子空间有最大的类间距离和最小的类内距离，即模式在该空间中有最佳的可分离性。因此，它是一种有效的特征抽取方法。使用这种方法能够使投影后模式样本的类间散布矩阵最大，并且同时类内散布矩阵最小。就是说，它能够保证投影后模式样本在新的空间中有最小的类内距离和最大的类间距离，即模式在该空间中有最佳的可分离性。

主成分分析（PCA）与线性判别（LDA）有着非常近似的意思，LDA的输入数据是带标签的，而PCA的输入数据是不带标签的，所以PCA是一种无监督的学习方式（unsupervised learning）。LDA通常来说是作为一个独立的算法存在，给定了训练数据后，将会得到一系列的判别函数，之后对于新的输入，就可以进行预测了。而PCA更像是一个预处理的方法，它可以将原本的数据降低维度，而使得降低了维度的数据之间的方差最大。

PCA仅仅是将图像数据矢量化再根据特征值和特征向量进行特征提取，而忽略了张量对象往往是多阶数据这个因素，从而造成了只在一个方向上投射而导致的降维效果的不理想以及破坏了原始图像的结构和相关性，不能完全保持原始图像中的冗余和结构等问题。

发明内容

鉴于PCA的局限性，我们提出将一种改进的PCA方法应用于脑fMRI数据分析：多线性主元分析法（简称MPCA)。在此基础上，我们又结合了LDA的特点，最终形成了MPCA+LDA的方法应用于脑fMRI数据分析中，可以通过该方法来判定大脑所处的认知任务状态。

本发明针对目前生物特征提取和分类判别领域中存在的问题，提出了一种应用于大脑认知状态判定领域的多线性主元分析方法。

改进的PCA，我们称之为多线性主元分析法（Multilinear PrincipalComponent Analysis，简称MPCA)。这种方法延续了经典PCA的范式，直接对原始张量数据进行处理，利用多线性主元分析（MPCA）考虑了在一个局部领域内的图像像素的空间相关性，从多个方向上对脑fMRI数据张量对象进行降维和特征提取，克服了传统PCA单纯的进行降维而破坏了原始图像数据的结构和相关性，不能完全保持原始图像数据中的冗余和结构的不足。这种MPCA方法是一种用来在所有张量模式方向上降维并在每一个模式方向上寻求它们的基础成分（主要成分）的多线性算法，它允许仿射张量捕获出现在原始张量上的大多数变化。

结合多线性主元分析（MPCA）与线性判别分析（LDA）的特点，利用MPCA-LDA算法进行大脑所处认知任务状态的鉴别。通过MPCA算法求得训练样本的特征子空间，并在此基础上计算出LDA算法的特征子空间。训练样本与测试样本分别向上述特征子空间投影，从而得到判别特征。最后利用最近领域法分类器完成识别分类。

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

1)输入样本集，对输入数据集进行预处理；

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

3)将上述MPCA降维后的特征向量集作为LDA的输入，求出LDA最优投射矩阵，并将MPCA降维后的特征向量集投射到LDA特征子空间进一步提取训练集的判别特征向量；

4)特征分类，将训练图像和测试图像的投射所得判别特征向量进行特征匹配，进而分类。

所述步骤2)中，计算训练样本集的特征分解，求出特征张量及最优投射矩阵

的具体步骤是：

对于选取的训练样本，这里定义M个样本

组成的集合为{χ₁,χ₂,…,χ_M}，每个样本对应的张量空间为

其中I_n(1≤n≤N)是张量的n-mode（模式）维度，N表示张量的阶数。多线性主元分析的目标是寻找多线性变换

它将原始张量空间

投射到张量子空间

R^{P_{1}} &CircleTimes; R^{P_{2}} &CircleTimes; \cdot \cdot \cdot &CircleTimes; R^{P_{N}} (P_{n} < I_{n}, n = 1, \cdot \cdot \cdot, N) :

Y_{m} = {\overset{&OverBar;}{χ}}_{m} \times_{1} {U^{(1)}}^{T} \times_{2} {U^{(2)}}^{T} \cdot \cdot \cdot \times_{N} {U^{(N)}}^{T}, m = 1, \cdot \cdot \cdot, M

其中，是原始训练样本中心化后的张量数据。

MPCA的目标函数是确定使投影后的总张量离散度

最大化的N个投射矩阵

满足:

{U^{(n)}, n = 1, \cdot \cdot \cdot, N} = \arg ma x_{U^{(1)}, U^{(2)}, \cdot \cdot \cdot, U^{(N)}} Ψ_{Y}

这样

就能捕获原始样本中的大多数变化，更好的表示出各个样本的最主要特征。

2.1)对样本数据进行中心化处理

其中

μ = \frac{1}{M} Σ_{m = 1}^{M} χ_{m}

是样本均值。

2.2)初始化。构造每个mode上的协方差矩阵

Φ^{{(n)}^{*}} = Σ_{m = 1}^{M} {\overset{&OverBar;}{X}}_{m (n)} \cdot {\overset{&OverBar;}{X}}_{m (n)}^{T}

其中，

是张量

的n-mode展开矩阵，

是将中心化后的张量从第n-mode方向展开而求得的协方差矩阵。

再分别对各mode方向上的协方差矩阵

分别进行奇异值分解，求出其特征值

{λ⁽ⁿ⁾,1≤n≤N}（这里将特征值按从大到小的顺序排列）和对应的特征向量

{{\tilde{U}}^{(n)}, 1 \leq n \leq N} .

一般情况下，前10%甚至1%的奇异值的和就占了全部的奇异值总和的99%以上了，

因此，我们可以用前k（远小于矩阵的秩）大的奇异值来近似的表示所有奇异值，也就可以用来近似的描述矩阵了。我们按照保留总信息量的97%来选取特征值个数，也就是选取的特征值的和占全部特征值总和的比率为97%，即

Σ_{i_{n} = 1}^{k} λ_{i_{n}}^{(n)} / Σ_{i_{n} = 1}^{r} λ_{i_{n}}^{(n)} \leq β

其中，r表示协方差矩阵的秩，β表示在n-mode方向前k（小于r）个最大特征值的和占全部特征值总和的比值，这里设为定值0.97，其值可以根据实际需要自行设定。表示第i_n个完全投射的n-mode方向上的特征值。

这样就得到了选取的各个mode方向上特征值个数，同时从各个方向上对原始张量数据进行降维，达到了在保留原结构信息的基础上还能实现维数的下降，进而减少计算量。最后将得到的各个mode方向上各自前k（远小于矩阵的秩）个的最大特征值对应的特征向量组成的特征矩阵记作

将其赋值给

作为它的初始值。

2.3)局部优化。就是需要对

作进一步的优化。

当选择更新U⁽ⁿ⁾时，就保持原有的{U⁽¹⁾,U⁽²⁾,…U^(n-1),U⁽ⁿ⁺¹⁾…U^(N-1),U^(N)}不变，则

Φ^{(n)} = Σ_{m = 1}^{M} {\overset{&OverBar;}{X}}_{m (n)} \cdot U_{Φ_{n}} \cdot U_{Φ_{n}}^{T} {\overset{&OverBar;}{X}}_{m (n)}^{T}, 1 \leq n \leq N

其中

U_{Φ_{n}} = U^{(N)} &CircleTimes; U^{(N - 1)} &CircleTimes; \cdot \cdot \cdot {&CircleTimes; U}^{(n + 1)} &CircleTimes; U^{(n - 1)} &CircleTimes; \cdot \cdot \cdot U^{(2)} &CircleTimes; U^{(1)},

表示克罗内克乘积，求取Φ⁽ⁿ⁾的前P_n个最大特征值所对应的特征向量赋值给U⁽ⁿ⁾完成对它的更新。依次取n=1，2，…，N。则完成对

的一次寻优。之后通过中心化的训练样本与更新的投射矩阵相乘，得到

Y_{m} = {\overset{&OverBar;}{χ}}_{m} \times_{1} {U^{(1)}}^{T} \times_{2} {U^{(2)}}^{T} \cdot \cdot \cdot \times_{N} {U^{(N)}}^{T}, m = 1, \cdot \cdot \cdot, M

这样就将原始张量空间

投射到更新之后的张量子空间

R^{P_{1}} &CircleTimes; R^{P_{2}} &CircleTimes; \cdot \cdot \cdot &CircleTimes; R^{P_{N}} (P_{n} < I_{n}, n = 1, \cdot \cdot \cdot, N) .

通过判断是否满足

代表更新前的总张量离散度，

代表更新后的总张量离散度）这个终止条件来决定是否还需要迭代以更新

如果不满足则重复局部优化步骤继续寻优，直到满足条件退出局部优化迭代过程。

具体步骤如下：

(a)对m=1,…,M，计算

{Y_{m} = {\overset{&OverBar;}{χ}}_{m} \times_{1} {\tilde{U}}^{(1)}^{T} \times_{2} {\tilde{U}}^{(2)}^{T} \cdot \cdot \cdot \times_{N} {\tilde{U}}^{(N)}^{T}} \cdot

(b)计算初始张量离散度

由(a)步骤求得。

(c)对于k=1:K（K表示迭代次数）

对于n=1:N（N表示张量阶数）

令

U_{Φ_{n}} = U^{(N)} &CircleTimes; U^{(N - 1)} &CircleTimes; \cdot \cdot \cdot {&CircleTimes; U}^{(n + 1)} &CircleTimes; U^{(n - 1)} &CircleTimes; \cdot \cdot \cdot U^{(2)} &CircleTimes; U^{(1)},

按照

Φ^{(n)} = Σ_{m = 1}^{M} {\overset{&OverBar;}{X}}_{m (n)} \cdot U_{Φ_{n}} \cdot U_{Φ_{n}}^{T} {\overset{&OverBar;}{X}}_{m (n)}^{T}, 1 \leq n \leq N,

求取Φ⁽ⁿ⁾的前P_n个最大特征值对应的特征向量赋值给U⁽ⁿ⁾完成对它的更新。

计算{Y_m,m=1…,M}和

（k表示迭代次数）。如果

（η是用户定义的比较小的阈值）或者满足k=K，跳出循环，得到局部最优的N个投射矩阵

进入到下一步。

(d)投射。中心化处理后的训练图像向MPCA特征张量子空间上投影，得到训练图像的特征张量为

{{Y}_{m} = χ_{m} \times_{1} {U^{(1)}}^{T} \times_{2} {U^{(2)}}^{T} \cdot \cdot \cdot \times_{N} {U^{(N)}}^{T}, m = 1, \cdot \cdot \cdot, M},

之后直接矢量化得到训练集的特征向量集

{{\tilde{y}}_{m}, m = 1, \cdot \cdot \cdot, M} .

所述步骤3)中，将多线性主元分析(MPCA)降维后训练集的特征向量集作为线性判别分析(LDA)的输入，求取最优LDA投射子空间V_lda，具体步骤是：

运用LDA寻找使得类间散布矩阵S_B和类内散布矩阵S_W比值最大的投射子空间V_lda，最优的投射子空间

V_{lda} = \arg \max \frac{| V^{T} S_{B} V |}{| V^{T} S_{W} V |} = [\begin{matrix} v_{1} & v_{2} & \cdot \cdot \cdot & v_{H_{z}} \end{matrix}] .

其中类内散布矩阵

S_{W} = Σ_{m = 1}^{M} ({\tilde{y}}_{m} - {\overset{&OverBar;}{y}}_{c_{m}}) {({\tilde{y}}_{m} - {\overset{&OverBar;}{y}}_{c_{m}})}^{T},

{\overset{&OverBar;}{y}}_{c_{m}} = \frac{1}{N_{c}} Σ_{m}^{c_{m} = C} {\tilde{y}}_{m},

类间散布矩阵

S_{B} = Σ_{c = 1}^{C} N_{c} ({\overset{&OverBar;}{y}}_{c} - \overset{&OverBar;}{y}) {({\overset{&OverBar;}{y}}_{c} - \overset{&OverBar;}{y})}^{T},

\overset{&OverBar;}{y} = \frac{1}{M} Σ_{m} {\tilde{y}}_{m} .

这里M表示训练样本总个数，c_m表示类别标签，C表示训练样本总类别个数，N_c表示训练样本类别标签是c_m的样本总个数，

表示训练样本向最优MPCA投射矩阵投影并矢量化之后的的特征向量。

是满足

即

的前H_z(=C-1)个最大的特征值，而

是特征值

对应的广义特征向量。即得到LDA算法的最优投射矩阵

及其对应的特征值

{λ_{h_{z}}, h_{z} = 1, \cdot \cdot \cdot, H_{z}} .

所述步骤4)中，特征分类，将训练样本和测试样本的判别特征向量进行特征匹配，进而分类，具体步骤是：

首先选中心化的训练样本集和测试样本集它们分别向MPCA和LDA的最优投射子空间投影，从而得到它们各自的判别特征向量集{y_m,m=1,…,M}和{y_l,l=1,…,L}，其中M,L分别表示训练样本集和测试样本集的样本数。

下面引用的测试判别特征集表示包含所有测试样本判别特征向量的集合，训练判别特征集表示包含所有训练样本判别特征向量的集合。或者分别简称测试集和训练集。

之后利用最近领域分类算法计算出测试判别特征集中各个子集与训练判别特征集中每个子集的相似度，用匹配分数S(p,g)来表示

S (p, g) = (\frac{1}{N_{p}}) Σ_{n_{p} = 1}^{N_{p}} S (z_{n_{p}}, g) + (\frac{1}{N_{g}}) Σ_{n_{g} = 1}^{N_{g}} S (z_{n_{g}}, p)

其中，

表示测试判别特征集中类别为p的每一个测试样本特征向量

和训练序列g（即训练集中类别为g的样本组成的序列）的相似度，

表示训练判别特征集中类别为的g每一个训练样本特征向量

和测试序列p（即测试集中类别为p样本组成的序列）的相似度，S(p,g)表示测试序列p和训练序列g的匹配分数以及训练序列g和测试序列p的匹配分数的平均值的和，其最高的匹配分数就表示测试序列p匹配训练序列g，即将测试序列p中的所有测试样本判断归类为训练序列g所属的类别。

S (y, c) = - \min {d (y, z_{n_{c}}), n_{c} = 1,2, \cdot \cdot \cdot, N_{c}}

表示一个测试样本判别特征向量y和训练判别特征集中某一个类c中的N_c个训练样本判别特征向量的相似度。其修正马氏距离计算公式是其中选取的特征值个数H=H_z（H_z=C-1），

本发明的特点和优势在于：

1.本发明利用多线性主元分析（MPCA）直接对多阶张量数据进行降维和特征提取，克服了传统PCA单纯的进行降维而破坏了原始图像数据的结构和相关性，不能完全保持原始图像中的冗余和结构的不足，保留了fMRI成像数据的空间结构信息。

2.结合带类别标签的线性判别（LDA）算法，获取MPCA-LDA的特征空间，对不同类数据之间的差别进行建模，克服了单纯依赖PCA降维而可能导致的不相关特征引入、过度拟合等问题，便于后序的分类器获得更好的识别分类效果。

附图说明

图1为本发明算法的流程示意图；

具体实施方式

以下结合具体实施例，对本发明进行详细说明。

步骤1，结合附图1描述的模块1：预处理部分。

首先，采集行为学实验功能数据，对功能数据进行预处理，其目的是去除数据在采集过程中混入的与任务无关的一些影响因素，以提高图像的信噪比，更好地提高数学建模和分析的效力。这里采用SPM软件对脑fMRI数据进行预处理，其步骤包括：

1）时间片校正，目的是为了校正组成体素的层与层之间获取时间点上的差异。因为时间片的扫描方式有多种，每个切片是在不同的时间点上获取的，这种差异会对统计分析带来一定的影响。

2）空间校正，也被称为头部校正，由于脑功能磁共振成像实验持续时间较长，测量次数也较多，被试的呼吸、心跳、血流脉动等生理上的因素或者外面环境上的因素造成头部运动的情况不可避免，这样就会导致在不同的时刻采集的数据出现偏移，破坏了图像之间的对应关系。考虑到是同一被试者同一种成像方式，目前一般采用刚性变换（rigid transformation）方法求解校正参数，来进行头动校正。

3）图像融合，上述空间校正的刚性变换方法仅针对同一被试者的同一中成像方式，当对于同一被试者的不同成像方式所获取的图像时，由于它们之间没有足够的可比性，这时就需要使用图像融合的方法来进行空间校正。首先，可以根据图像信号强度，人为的选择不同的阈值将脑组织的灰质（GM）、白质（WM）和脑脊液（CSF）三个部分分割开；分割完后，同一种组织的不同成像方式所得到的图像之间就有了足够的可比性，这样就可以使用空间校正的方法来进行空间对齐了。

4）图像标准化，鉴于不同被试间脑结构的差异性，就需要将不同被试者的脑图像映射到一个统一的标准模板空间中，这时刚性变换就不再适用，而要利用整体形变的仿射变换和局部非线性变换来进行标准化操作，这就是图像标准化。目前比较常用的标准模板主要是Talairach模板和MNI（Montreal NeurologicalInstitute）模板，这两者整体上基本相同，除了一些细节上的区别外。

5）空间平滑，就是将图像数据在空间上用一个光滑函数（通常是高斯核函数）去进行卷积运算。平滑有以下几个优点：可以提高信噪比，脑功能磁共振成像检测到的是神经元激活所产生的血液动力学信号，在图像重建中主要对应着低频域部分，而噪声则对应着信号的高频部分，经过空间平滑后，噪声将会得到很大的抑制，这样就可以提高图像的信噪比了；使数据符合高斯随机场的要求，这对利用高斯随机场理论来做统计推断是很重要的，因为这样就可以提高统计推断的准确度；可以消除不同被试者间脑结构的差别，对于需要使用不同被试者平均结果的实验来说，这些差异也会带来很大的影响，而经过空间平滑后这种差异性就会被模糊化，同时不会滤除掉有意义的高频部分的信息。

然后，将SPM软件处理后的原始图像数据按照扫面时间序列以及行为实验类别进行重新组合，构成试验数据。

最后，再将数据进行归一化处理。

步骤2，结合附图1描述的模块2：特征提取部分。

本模块包含两个部分：MPCA投射空间和LDA投射空间的产生。

第一部分，MPCA投射空间的产生：

对于选取的训练样本，这里定义M个样本组成的集合为{χ₁,χ₂,…,χ_M}，样本对应的张量空间为其中I_n(1≤n≤N)是张量的n-mode（模式）维度，N表示张量的阶数。多线性主元分析的目标是寻找多线性变换

它将原始张量空间

投射到张量子空间

R^{P_{1}} &CircleTimes; R^{P_{2}} &CircleTimes; \cdot \cdot \cdot &CircleTimes; R^{P_{N}} (P_{n} < I_{n}, n = 1, \cdot \cdot \cdot, N) :

Y_{m} = {\overset{&OverBar;}{χ}}_{m} \times_{1} {U^{(1)}}^{T} \times_{2} {U^{(2)}}^{T} \cdot \cdot \cdot \times_{N} {U^{(N)}}^{T}, m = 1, \cdot \cdot \cdot, M

其中，

是原始训练样本中心化后的张量数据。

MPCA的目标函数是确定使投影后的总张量离散度

最大化的N个最优投射矩阵即满足:

{U^{(n)}, n = 1, \cdot \cdot \cdot, N} = \arg ma x_{U^{(1)}, U^{(2)}, \cdot \cdot \cdot, U^{(N)}} Ψ_{Y}

其中投影所得张量集{Y₁,Y₂,…,Y_M}的总离散度Ψ_Y定义为：

Ψ_{Y} = Σ_{m = 1}^{M} {| | Y | |}^{2}

这样

1）对样本数据进行中心化处理

其中

μ = \frac{1}{M} Σ_{m = 1}^{M} χ_{m}

是样本均值。

2）初始化。构造原始张量样本数据中心化后数据集的每个mode方向上的协方差矩阵即

Φ^{{(n)}^{*}} = Σ_{m = 1}^{M} {\overset{&OverBar;}{X}}_{m (n)} \cdot {\overset{&OverBar;}{X}}_{m (n)}^{T} (n = 1,2, \cdot \cdot \cdot N)

其中，

是张量

的n-mode展开矩阵，是

的转置矩阵，

再分别对各mode方向上展开的协方差矩阵

分别进行奇异值分解，求出其特征值{λ⁽ⁿ⁾,1≤n≤N}（这里将特征值按从大到小的顺序排列）和对应的特征向量

{{\tilde{U}}^{(n)}, 1 \leq n \leq N} .

因此，我们可以用前k（远小于矩阵的秩）个的最大奇异值来近似的表示所有奇异值，也就可以用来近似的描述矩阵了。我们按照保留总信息量的97%来选取特征值个数，也就是选取的特征值的和占全部特征值总和的比率为97%，即

Σ_{i_{n} = 1}^{P_{n}} λ_{i_{n}}^{(n)} / Σ_{i_{n} = 1}^{r} λ_{i_{n}}^{(n)} \leq β

其中，r表示协方差矩阵的秩，β表示在n-mode方向前P_n(n=1,…,N)（小于r）个最大特征值的和占全部特征值总和的比值，这里设β为定值0.97，其值可以根据实际需要自行设定。

表示第i_n个完全投射的n-mode方向上的特征值。

这样就得到了选取的各个mode方向上特征值个数，同时从各个方向上对原始张量数据进行降维，达到了在保留多数原结构信息的基础上还能实现维度的下降，进而减少计算量。最后将得到的各个mode方向上各自前k（远小于矩阵的秩）个的最大特征值对应的特征向量组成的特征矩阵记作

将其赋值给

作为它的初始值。

3)局部优化。就是需要对

作进一步的优化。

Φ^{(n)} = Σ_{m = 1}^{M} {\overset{&OverBar;}{X}}_{m (n)} \cdot U_{Φ_{n}} \cdot U_{Φ_{n}}^{T} {\overset{&OverBar;}{X}}_{m (n)}^{T}, 1 \leq n \leq N

其中

U_{Φ_{n}} = U^{(N)} &CircleTimes; U^{(N - 1)} &CircleTimes; \cdot \cdot \cdot {&CircleTimes; U}^{(n + 1)} &CircleTimes; U^{(n - 1)} &CircleTimes; \cdot \cdot \cdot U^{(2)} &CircleTimes; U^{(1)},

表示矩阵克罗内克乘积运算，求取Φ⁽ⁿ⁾的前P_n个最大特征值所对应的特征向量赋值给U⁽ⁿ⁾完成对它的更新。依次取n=1,2,…,N。则完成对

的一次寻优。之后将通过中心化的训练样本与更新的投射矩阵相乘，得到

Y_{m} = {\overset{&OverBar;}{χ}}_{m} \times_{1} {U^{(1)}}^{T} \times_{2} {U^{(2)}}^{T} \cdot \cdot \cdot \times_{N} {U^{(N)}}^{T} (m = 1, \cdot \cdot \cdot, M)

这样就将原始张量空间

投射到更新之后的张量子空间

R^{P_{1}} &CircleTimes; R^{P_{2}} &CircleTimes; \cdot \cdot \cdot &CircleTimes; R^{P_{N}} (P_{n} < I_{n}, n = 1, \cdot \cdot \cdot, N) .

通过判断是否满足

代表更新前的总张量离散度，

具体步骤如下：

(a)对m=1,…,M，计算

{Y_{m} = {\overset{&OverBar;}{χ}}_{m} \times_{1} {U^{(1)}}^{T} \times_{2} {U^{(2)}}^{T} \cdot \cdot \cdot \times_{N} {U^{(N)}}^{T}} \cdot

(b)计算初始张量离散度由(a)步骤求得。

(c)对于k=1:K（K表示迭代次数）

对于n=1:N（N表示张量阶数）

令

U_{Φ_{n}} = U^{(N)} &CircleTimes; U^{(N - 1)} &CircleTimes; \cdot \cdot \cdot {&CircleTimes; U}^{(n + 1)} &CircleTimes; U^{(n - 1)} &CircleTimes; \cdot \cdot \cdot U^{(2)} &CircleTimes; U^{(1)},

按照

Φ^{(n)} = Σ_{m = 1}^{M} {\overset{&OverBar;}{X}}_{m (n)} \cdot U_{Φ_{n}} \cdot U_{Φ_{n}}^{T} {\overset{&OverBar;}{X}}_{m (n)}^{T}, 1 \leq n \leq N,

计算{Y_m,m=1…,M}和

（k表示迭代次数）。如果

（η是用户定义的比较小的阈值）或者满足k=K，跳出循环，得到局部最优投射矩阵

{U^{(n)} &Element; R^{I_{n} \times P_{n}}, n = 1, \cdot \cdot \cdot, N},

进入到下一步。

(d)投射。训练图像数据χ_m(m=1,…,M)经过MPCA后得到投射矩阵

将其中心化后的训练样本

与最优投射矩阵

相乘:

Y_{m} = {\overset{&OverBar;}{χ}}_{m} \times_{1} {U^{(1)}}^{T} \times_{2} {U^{(2)}}^{T} \cdot \cdot \cdot \times_{N} {U^{(N)}}^{T}, (m = 1, \cdot \cdot \cdot, M) .

实现向MPCA特征子空间上投影，得到训练集的特征张量集{Y_m,m=1,…,M}，之后直接矢量化Y_m得其矢量化形式

即得到变化后的特征向量集

第二部分，LDA投射空间的产生：

将多线性主元分析(MPCA)降维后的特征向量集

作为线性判别分析(LDA)的输入，求取最优LDA投射子空间V_lda，以便进一步提取训练样本的判别特征向量集{y_m,m=1,…,M}。

运用LDA寻找使得类间散布矩阵S_B和类内散布矩阵S_W比值最大的投影子空间V_lda，最优的投射子空间

V_{lda} = \arg \max \frac{| V^{T} S_{B} V |}{| V^{T} S_{W} V |} = [\begin{matrix} v_{1} & v_{2} & \cdot \cdot \cdot & v_{H_{z}} \end{matrix}] .

其中类内散布矩阵

S_{W} = Σ_{m = 1}^{M} ({\tilde{y}}_{m} - {\overset{&OverBar;}{y}}_{c_{m}}) {({\tilde{y}}_{m} - {\overset{&OverBar;}{y}}_{c_{m}})}^{T},

{\overset{&OverBar;}{y}}_{c_{m}} = \frac{1}{N_{c}} Σ_{m}^{c_{m} = C} {\tilde{y}}_{m},

类间散布矩阵

S_{B} = Σ_{c = 1}^{C} N_{c} ({\overset{&OverBar;}{y}}_{c} - \overset{&OverBar;}{y}) {({\overset{&OverBar;}{y}}_{c} - \overset{&OverBar;}{y})}^{T},

\overset{&OverBar;}{y} = \frac{1}{M} Σ_{m} {\tilde{y}}_{m} .

表示训练样本向最优MPCA投射矩阵投影并矢量化之后的的特征向量。是满足

即

的前H_z(≤C-1)个最大的特征值，而

是特征值

对应的广义特征向量。即得到LDA算法的特征值及其对应的最优投射矩阵

V_{lda} = {v_{h_{z}}, h_{z} = 1, \cdot \cdot \cdot, H_{z}} .

步骤3，结合附图1描述的模块3：识别分类部分。

首先，将中心化的训练样本集和测试样本集分别与步骤2得到的最优投射矩阵

相乘，即投射到MPCA特征张量子空间，得到各自的特征张量

{Y_{s} = \overset{&OverBar;}{χ} \times_{1} {U^{(1)}}^{T} \times_{2} {U^{(2)}}^{T} \cdot \cdot \cdot \times_{N} {U^{(N)}}^{T}, s = 1, 2, \cdot \cdot \cdot, M, M + 1, M + 2, \cdot \cdot \cdot M + L},

之后矢量化Y_s得其特征向量集在这一模块里的M,L分别代表训练样本集和测试样本集元素的个数。

将多线性主元分析(MPCA)降维后的特征向量集作为线性判别分析(LDA)的输入与步骤2得到的V_lda相乘，即得到样本在最优LDA投射子空间V_lda的投影{y_s,s=1,2,…,M,M+1,M+2,…M+L}。实现提取训练样本判别特征向量和测试样本判别特征向量。将包含所有训练样本判别特征向量的集合{y_m,m=1,…,M}称作训练判别特征集，包含所有测试样本判别特征向量的集合{y_l,l=1,…,L}称作测试判别特征集。

然后，利用最近邻域法分类器完成识别分类。将包含所有训练样本判别特征向量的训练判别特征集和包含所有测试样本判别特征向量的测试判别特征集进行特征匹配以完成分类，匹配就是用最近邻分类算法计算出每个测试(训练)样本判别特征向量与训练(测试)判别特征集中每个类的相似度，再利用匹配分数S(p,g)来表示测试判别特征集中类别为p子集与训练判别特征集中类别为g子集的相似度。具体步骤是：

利用最近邻分类算法来计算每个测试样本判别特征向量与训练判别特征集中某一个类c的相似度，在这里用

表示一个测试样本判别特征向量y和训练判别集中某一个类c中的N_c个训练样本判别特征向量

的相似度。N_c表示训练判别特征集中属于类c的训练样本判别特征向量的个数，

表示一个测试样本判别特征向量y和训练判别特征集中某一个类c中的N_c个训练样本判别特征向量的修正马氏距离的最小值，即最近邻算法。其修正马氏距离计算公式是

其中选取的特征值个数H=H_z（H_z=C-1），

步骤2中已算得。

同理利用最近邻分类算法来计算每个训练样本判别特征向量与测试判别特征集中每一个类的相似度。

最后利用匹配分数S(p,g)来表示测试集中类别为p子集与训练集中类别为g子集的相似度：

S (p, g) = (\frac{1}{N_{p}}) Σ_{n_{p} = 1}^{N_{p}} S (z_{n_{p}}, g) + (\frac{1}{N_{g}}) Σ_{n_{g} = 1}^{N_{g}} S (z_{n_{g}}, p)

其中，

表示测试集中类别为p的每一个测试样本判别特征向量

和训练序列g（即训练集中类别为g的N_g个样本组成的序列）的相似度，

表示训练集中类别为g每一个训练样本判别特征向量

和测试序列p（即测试集中类别为p的N_g样本组成的序列）的相似度，N_p表示测试集中类别为p样本个数，N_g表示训练集中类别为g的样本个数。

假定训练集中包括了N个类别{g_i,i=1,2,…N}，分别计算类别为p测试子集与每个训练类别g_i子集的匹配分数。其最高的匹配分数就表示测试序列p匹配训练序列g_i，即将测试序列p中的所有样本判断归类为训练序列g_i所属的类别。依次确定测试集中每个序列对应的所属的类别，最终完成大脑认知状态判别。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims

1.一种基于多线性主元分析的大脑认知状态判定方法，其特征在于，包括如下步骤：

1)输入样本集，对输入数据进行预处理；

2)计算训练样本集的特征分解，求出特征变换矩阵，将训练样本到张量特征子空间，得到训练集的特征张量集；

3)将降维后的低维特征张量数据矢量化作为LDA的输入，求出LDA最优投射矩阵，并将矢量化的低维特征张量数据投射到LDA特征子空间进一步提取训练集的判别特征向量；

2.根据权利要求1所述的一种基于多线性主元分析的大脑认知状态判定方法，其特征在于：所述步骤2)中，计算训练样本集的特征分解，求出特征张量以及最优特征变换矩阵

的具体步骤是：

对于选取的训练样本，这里定义M个样本

组成的集合为{χ₁,χ₂,…,χ_M}，每个样本对应的张量空间为

其中I_n(1≤n≤N)是张量的n-mode（模式）维度，N表示张量的阶数；多线性主元分析的目标是寻找多线性变换矩阵

它将原始张量空间

投射到张量子空间

R^{P_{1}} &CircleTimes; R^{P_{2}} &CircleTimes; \cdot \cdot \cdot &CircleTimes; R^{P_{N}} (P_{n} < I_{n}, n = 1, \cdot \cdot \cdot, N) :

Y_{m} = {\overset{&OverBar;}{χ}}_{m} \times_{1} {U^{(1)}}^{T} \times_{2} {U^{(2)}}^{T} \cdot \cdot \cdot \times_{N} {U^{(N)}}^{T}, m = 1, \cdot \cdot \cdot, M

其中，

是原始样本中心化后的张量数据；

MPCA的目标函数是确定使投射后的总张量离散度

最大化的N个投射矩阵

满足:

{U^{(n)}, n = 1, \cdot \cdot \cdot, N} = \arg ma x_{U^{(1)}, U^{(2)}, \cdot \cdot \cdot, U^{(N)}} Ψ_{Y}

这样

就能捕获原始样本中的大多数变化，更好的表示出各个样本的最主要特征；

3.1)对训练样本数据进行中心化处理

其中

μ = \frac{1}{M} Σ_{m = 1}^{M} χ_{m}

是样本均值；

3.2)初始化；对训练样本每个mode上的协方差矩阵

进行奇异值特征分解，其中是张量的n-mode展开矩阵；

分别对各mode方向上的协方差矩阵

{{\tilde{U}}^{(n)}, 1 \leq n \leq N};

然后再按照保留总信息量的97%来选取特征值个数，也就是选取的特征值的和占全部特征值总和的比率为97%，即

Σ_{i_{n} = 1}^{P_{n}} λ_{i_{n}}^{(n)} / Σ_{i_{n} = 1}^{r} λ_{i_{n}}^{(n)} \leq β

其中，r表示协方差矩阵的秩，β表示在n-mode方向前P_n(n＝1,…,N)（小于r）个最大特征值的和占全部特征值总和的比值，这里设为定值0.97，其值可以根据实际需要自行设定；

表示第i_n个完全投射的n-mode方向上的特征值；

这样就得到了选取的各个mode方向上特征值个数，同时从各个方向上对原始张量数据进行降维，达到了在保留原结构信息的基础上还能实现维数的下降，进而减少计算量；最后将得到的各个mode方向上各自前k（远小于矩阵的秩）个的最大特征值对应的特征向量组成的特征矩阵记作将其赋值给

作为它的初始值；

3.3)局部优化；就是需要对

作进一步的局部优化；

Φ^{(n)} = Σ_{m = 1}^{M} {\overset{&OverBar;}{X}}_{m (n)} \cdot U_{Φ_{n}} \cdot U_{Φ_{n}}^{T} {\overset{&OverBar;}{X}}_{m (n)}^{T}, 1 \leq n \leq N

其中

U_{Φ_{n}} = U^{(N)} &CircleTimes; U^{(N - 1)} &CircleTimes; \cdot \cdot \cdot {&CircleTimes; U}^{(n + 1)} &CircleTimes; U^{(n - 1)} &CircleTimes; \cdot \cdot \cdot U^{(2)} &CircleTimes; U^{(1)},

表示矩阵克罗内克乘积运算，求取Φ⁽ⁿ⁾的前P_n个最大的特征值所对应的特征向量赋值给U⁽ⁿ⁾完成对它的更新；依次取n＝1,2,…,N；则完成对

的一次寻优；之后通过中心化的训练样本与更新的投射矩阵相乘，得到：

Y_{m} = {\overset{&OverBar;}{χ}}_{m} \times_{1} {U^{(1)}}^{T} \times_{2} {U^{(2)}}^{T} \cdot \cdot \cdot \times_{N} {U^{(N)}}^{T}, m = 1, \cdot \cdot \cdot, M

这样就将原始张量空间

投射到了更新之后的张量子空间

R^{P_{1}} &CircleTimes; R^{P_{2}} &CircleTimes; \cdot \cdot \cdot &CircleTimes; R^{P_{N}} (P_{n} < I_{n}, n = 1, \cdot \cdot \cdot, N);

通过判断是否满足（代表更新前的总张量离散度，代表更新后的总张量离散度）这个终止条件来决定是否还需要迭代以更新

如果不满足则重复局部优化步骤继续寻优，直到满足条件退出局部优化迭代过程；

具体步骤如下：

(a)对m＝1,…,M，计算

{Y_{m} = {\overset{&OverBar;}{χ}}_{m} \times_{1} {\tilde{U}}^{(1)}^{T} \times_{2} {\tilde{U}}^{(2)}^{T} \cdot \cdot \cdot \times_{N} {\tilde{U}}^{(N)}^{T}};

(b)计算初始张量离散度由(a)步骤求得；

(c)对于k＝1:K（K表示迭代次数）

对于n＝1:N（N表示张量阶数）

令

U_{Φ_{n}} = U^{(N)} &CircleTimes; U^{(N - 1)} &CircleTimes; \cdot \cdot \cdot {&CircleTimes; U}^{(n + 1)} &CircleTimes; U^{(n - 1)} &CircleTimes; \cdot \cdot \cdot U^{(2)} &CircleTimes; U^{(1)},

按照

Φ^{(n)} = Σ_{m = 1}^{M} {\overset{&OverBar;}{X}}_{m (n)} \cdot U_{Φ_{n}} \cdot U_{Φ_{n}}^{T} {\overset{&OverBar;}{X}}_{m (n)}^{T}, 1 \leq n \leq N,

求取Φ⁽ⁿ⁾的前P_n个最大特征值对应的特征向量赋值给U⁽ⁿ⁾完成对它的更新；

计算{Y_m,m＝1…,M}和

（k表示迭代次数）；如果

（η是用户定义的比较小的阈值）或者满足k＝K，跳出循环，得到局部最优投射矩阵

{U^{(n)} &Element; R^{I_{n} \times P_{n}}, n = 1, \cdot \cdot \cdot, N},

进入到下一步；

(d)投射；训练图像数据χ_m(m＝1,…,M)经过MPCA后得到投射矩阵

将中心化后的训练样本向MPCA特征子空间上投影，即得到训练图像的特征张量为

{{Y}_{m} = {\overset{&OverBar;}{χ}}_{m} \times_{1} {U^{(1)}}^{T} \times_{2} {U^{(2)}}^{T} \cdot \cdot \cdot \times_{N} {U^{(N)}}^{T}, m = 1, \cdot \cdot \cdot, M};

之后直接矢量化Y_m得到MPCA变化后的特征向量集

3.根据权利要求1所述的一种基于多线性主元分析的大脑认知状态判定方法，其特征在于：所述步骤3)中，将多线性主元分析(MPCA)降维后的特征向量集作为线性判别分析(LDA)的输入，求取最优LDA投射子空间V_lda，具体步骤是：

V_{lda} = \arg \max \frac{| V^{T} S_{B} V |}{| V^{T} S_{W} V |} = [\begin{matrix} v_{1} & v_{2} & \cdot \cdot \cdot & v_{H_{z}} \end{matrix}];

其中类内散布矩阵

S_{W} = Σ_{m = 1}^{M} ({\tilde{y}}_{m} - {\overset{&OverBar;}{y}}_{c_{m}}) {({\tilde{y}}_{m} - {\overset{&OverBar;}{y}}_{c_{m}})}^{T},

{\overset{&OverBar;}{y}}_{c_{m}} = \frac{1}{N_{c}} Σ_{m}^{c_{m} = C} {\tilde{y}}_{m},

类间散布矩阵

S_{B} = Σ_{c = 1}^{C} N_{c} ({\overset{&OverBar;}{y}}_{c} - \overset{&OverBar;}{y}) {({\overset{&OverBar;}{y}}_{c} - \overset{&OverBar;}{y})}^{T},

\overset{&OverBar;}{y} = \frac{1}{M} Σ_{m} {\tilde{y}}_{m};

表示训练样本向最优MPCA投射矩阵投影并矢量化之后的的特征向量；

是满足

即

的前H_z(＝C-1)个最大的特征值，则

是特征值

对应的广义特征向量；即得到LDA算法的最优特征矩阵

及其对应的特征值

{λ_{h_{z}}, h_{z} = 1, \cdot \cdot \cdot, H_{z}} .

4.根据权利要求1所述的一种基于多线性主元分析的大脑认知状态判定方法，其特征在于：所述步骤4)中，特征分类，将训练图像和测试图像的判别特征向量进行特征匹配，进而分类，具体步骤是：

首先将中心化后的训练样本集和测试样本集分别向MPCA和LDA的最优投射子空间投影，从而得到它们各自的判别特征向量集{y_m,m＝1,…,M}和{y_l,l＝1,…,L}，其中M,L分别表示训练样本集和测试样本集的样本数；

下面引用的测试判别特征集表示包含所有测试样本判别特征向量的集合，训练判别特征集表示包含所有训练样本判别特征向量的集合，分别简称测试集和训练集；

S (p, g) = (\frac{1}{N_{p}}) Σ_{n_{p} = 1}^{N_{p}} S (z_{n_{p}}, g) + (\frac{1}{N_{g}}) Σ_{n_{g} = 1}^{N_{g}} S (z_{n_{g}}, p)

其中，

表示测试判别特征集中类别为p的每一个测试样本特征向量

表示训练判别特征集中类别为的g每一个训练样本特征向量

和测试序列p（即测试集中类别为p样本组成的序列）的相似度，S(p,g)表示测试序列p和训练序列g的匹配分数以及训练序列g和测试序列p的匹配分数的平均值的和，其最高的匹配分数就表示测试序列p匹配训练序列g，即将测试序列p中的所有测试样本判断归类为训练序列g所属的类别；

S (y, c) = - \min {d (y, z_{n_{c}}), n_{c} = 1,2, \cdot \cdot \cdot, N_{c}}

表示一个测试样本判别特征向量y和训练判别特征集中某一个类c中的N_c个训练样本判别特征向量

的相似度；其修正马氏距离计算公式是

其中选取的特征值个数H＝H_z（H_z＝C-1），