CN105912851A - 一种利用PCA和非环形特性估计复数fMRI数据模型阶数的方法 - Google Patents

Abstract

一种利用PCA和非环形特性估计复数fMRI数据模型阶数的方法，属于生物医学信号处理领域。先对复数fMRI数据进行连续主成分数的PCA消噪，再计算消噪数据的非环形度量DOI，得到DOI曲线并进行必要的调整，最后采用SORTE准则检测DOI曲线的拐点，该拐点对应的PCA成分数即为所估计模型阶数。本发明利用了完整的复数fMRI数据，能估计出更高更准确的模型阶数，进而分离更多更好的空时成分。在敲击手指任务下采集的复数fMRI数据的ICA分析中，本发明估计阶数下获取的单被试SM和TC与参考信号的相关系数最大可提高202.42％和51.89％以及123.15％和431.30％(DMN)。

Description

一种利用PCA和非环形特性估计复数fMRI数据模型阶数的方法

技术领域

本发明涉及生物医学信号处理领域，特别是涉及一种复数功能磁共振成像(functional magnetic resonance imaging,fMRI)数据的模型阶数估计方法。

背景技术

由于具有空间分辨率高、非侵入等优点，fMRI技术已成为脑功能分析和脑疾病诊断的有力工具。原始采集到的fMRI数据是复数，包括幅值数据和相位数据。幅值数据和相位数据都含有独特的脑功能信息，但相位数据的噪声水平较之幅值数据高得多，这导致复数fMRI数据的噪声水平较之幅值fMRI数据高得多。尽管如此，辅以适当的消噪环节，人们已经从复数fMRI数据中提取了更多更完整的脑功能信息，表明了复数fMRI数据分析在脑科学研究中的明显优势。

复数fMRI数据属于混合数据，其中包含着数十个脑空间激活区(spatial map，SM)成分及其相应的时间过程(time course，TC)成分，比幅值fMRI数据所包含的成分数多。在已知成分数的情况下，可以采用独立成分分析(independentcomponent analysis,ICA)、张量分解(tensor decomposition)等数据驱动方法提取出SM成分和TC成分。然而，当成分数过大时，有些成分会分裂成多个子成分；而当成分数过小时，会出现多个成分的混合。因此，成分数估计，也称模型阶数估计，是基于数据驱动方法分析fMRI数据的重要前提。

目前，用于估计复数fMRI数据模型阶数的方法主要分为两类，一类为基于信息论的Akaike’s information criterion(AIC)、Kullback-Leibler information criterion(KIC)和minimum description length(MDL)等方法，另一类是传统的寻找协方差矩阵特征值曲线拐点的方法。其中，第一类信息论方法要求输入数据满足独立同分布，而fMRI数据间存在相关性，所以在代入信息论准则计算模型阶数之前，需要对fMRI数据进行下采样以获取独立同分布子集，这便导致了该类方法的信息损失缺限。换句话说，应用信息论方法对复数fMRI数据进行模型阶数估计时，会因为数据下采样致使所估计成分数偏少，模型阶数偏低。此外，由于实际复数fMRI信号的高噪声水平，第二类特征值曲线拐点方法在估计模型阶数时难于得到正确的结果，而且由于设定阈值的需要，该方法具有很大的不确定性。总之，与成熟的幅值fMRI数据模型阶数估计方法相比，复数fMRI数据的模型阶数估计方法面临着较大的挑战。因此，探索适于高噪声完整复数fMRI数据的模型阶数估计方法，是今后深入开展复数fMRI数据分析的重要前提。

发明内容

本发明的目的在于，提供一种适于高噪声完整复数fMRI数据的模型阶数估计方法，通过引入PCA(principle component analysis)进行噪声抑制，通过直接利用复数fMRI成分的非环形特性提高对复数fMRI数据分析的针对性，通过利用碎石图(scree plot)元素二阶统计量的方法提高非环形度量拐点(即模型阶数)估计的准确性，进而显著提高ICA分析的性能。

本发明的技术方案是，首先对复数fMRI数据进行连续主成分数的PCA消噪，然后计算消噪数据的非环形度量DOI(degree of impropriety)，得到DOI曲线并进行必要的调整，最后采用SORTE(second order statistic of eigenvalue)准则检测DOI曲线的拐点，该拐点对应的PCA成分数即为所估计模型阶数。具体实现步骤如下：

第一步：输入复数fMRI数据X∈C^J×V。J表示时间维的全脑扫描次数，V表示空间维的脑内体素数目。

第二步：设置PCA主成分个数M＝1,…,N，N<J。由于复数fMRI数据较之幅值fMRI数据增加了相位信息，最大主成分数N可取为幅值fMRI数据所估计模型阶数的三倍。

第三步：进行PCA消噪，得到N组主成分数连续的消噪数据。设U＝[u₁,...,u_M]∈C^J×M包含X协方差矩阵特征分解前M个最大特征值对应的特征矢量，则Y＝U^HX∈C^M×V为PCA消噪数据，含有M个主成分。

第四步：计算N组PCA消噪数据的DOI。设置循环M＝1,…,N，分别计算各组消噪矩阵Y的DOI如下：

$\mathrm{\&rho;}=\frac{\det \left({\tilde{R}}_{YY}{R}_{YY}^{-*}{\tilde{R}}_{YY}^{*}\right)}{{\mathrm{detR}}_{YY}}---\left(1\right)$

式中，det表示求行列式值，R_YY＝E(YY^H)，E表示数学期望，上标H表示共轭转置，上标T表示转置，上标*表示共轭，上标-表示求逆；得到N组PCA消噪数据的DOI序列ρ₁,...,ρ_N。

第五步：检测和调整DOI曲线。理想情况下ρ₁≥ρ₂≥...≥ρ_N，即DOI曲线为递降曲线。但是，由于实际复数fMRI数据含有高水平噪声，DOI曲线可能出现个别波动点，不满足递减规律。为此，从前至后逐点检测DOI值，当DOI序列中某个DOI值高于前一个DOI值时，令其等于前一个DOI值，即当ρ_i>ρ_i-1,i＝2,...,N时，令ρ_i＝ρ_i-1。

第六步：计算SORTE序列：

$SORTE\left(p\right)=\left\{\begin{array}{cc}\frac{\mathrm{var}\&lsqb;{\{\&dtri;{\mathrm{\&rho;}}_i\}}_{i=p+1}^{N-1}\&rsqb;}{\mathrm{var}\&lsqb;{\{\&dtri;{\mathrm{\&rho;}}_i\}}_{i=p}^{N-1}\&rsqb;},& \mathrm{var}\&lsqb;{\{\&dtri;{\mathrm{\&rho;}}_i\}}_{i=p}^{N-1}\&rsqb;>0\\ +\mathrm{\&infin;},& \mathrm{var}\&lsqb;{\{\&dtri;{\mathrm{\&rho;}}_i\}}_{i=p}^{N-1}\&rsqb;=0\end{array}\right.---\left(2\right)$

其中，为调整后DOI曲线的差值序列，p＝1,...,(N-2)，

$\mathrm{var}\&lsqb;{\{\&dtri;{\mathrm{\&rho;}}_i\}}_{i=p}^{N-1}\&rsqb;=\frac{1}{N-p}{\mathrm{\&Sigma;}}_{i=p}^{N-1}{(\&dtri;{\mathrm{\&rho;}}_i-\frac{1}{N-p}{\mathrm{\&Sigma;}}_{i=p}^{N-1}\&dtri;{\mathrm{\&rho;}}_i)}^2$

第七步：计算模型阶数。根据SORTE准则，计算模型阶数P如下：

$P=\arg \underset{p=1,...,(N-3)}{min}SORTE\left(p\right)---\left(3\right)$

第八步：输出模型阶数估计结果P。

本发明所达到的效果和益处是，较之含有独立同分布下采样的信息论方法，本发明利用了完整的复数fMRI数据，能估计出更高的模型阶数，进而分离更多的空时成分。而且，通过利用PCA降噪和复数fMRI数据的非环形特性，本发明估计的模型阶数在ICA分析中体现出明显的性能改善，表明了估计的准确性。例如，在16被试敲击手指任务下采集的复数fMRI数据分析中，信息论MDL方法估计的模型阶数平均值为27，本发明估计的模型阶数平均值为35。应用复数ICA的EBM算法在两种模型阶数下进行ICA分解，以任务相关成分和默认网络(defaultmode network,DMN)为感兴趣信号，较之信息论结果，本发明获取的单被试SM与参考SM的相关系数最大可提高202.42％(任务相关成分)和123.15％(DMN)，单被试TC与Model TC的相关系数最大可提高51.89％(任务相关成分)和431.30％(DMN)。因此，本发明为复数fMRI数据分析提供了有力的保障。

附图说明

图1是本发明估计单被试复数fMRI数据模型阶数的工作流程图。

具体实施方式

下面结合技术方案和图1，详细叙述本发明的一个具体实施例。

现有单被试在执行敲打手指任务下采集的复数fMRI数据。时间维度上进行了J＝165次扫描，每次扫描都获得了53×63×46的全脑数据，脑内体素数V＝59610。应用信息论MDL准则估计该单被试幅值fMRI数据的模型阶数为20。采用本发明进行模型阶数估计的步骤如附图所示。

第一步：输入单被试复数fMRI数据X∈C^165×59610。

第二步：设置PCA主成分个数M＝1,…,N。已知幅值fMRI数据的模型阶数为20，故设置最大主成分数N＝60。

第三步：进行PCA消噪，得到60组主成分数连续的消噪数据Y∈C^M×59610，M＝1,…,60。

第四步：代入公式(1)，循环计算60组PCA消噪数据的DOI，得到DOI序列ρ₁,...,ρ₆₀。

第五步：检测和调整DOI曲线。当ρ_i>ρ_i-1,i＝2,...,60时，令ρ_i＝ρ_i-1。

第六步：代入公式(2)，计算SORTE序列SORTE(p)，p＝1,...,58。

第七步：代入公式(3)，计算模型阶数P。SORTE(1),…,SORTE(57)最小值对应的p＝38，故模型阶数P＝38。

第八步：输出模型阶数估计结果P＝38。

Claims (2)

1.一种利用PCA和非环形特性估计复数fMRI数据模型阶数的方法，首先对复数fMRI数据进行连续主成分数的PCA消噪，然后计算消噪数据的非环形度量DOI，得到DOI曲线并进行必要的调整，最后采用SORTE准则检测DOI曲线的拐点，该拐点对应的PCA成分数即为模型阶数；具体实现步骤如下：

第一步：输入复数fMRI数据X∈C^J×V；J表示时间维的全脑扫描次数，V表示空间维的脑内体素数目；

第二步：设置PCA主成分个数M＝1,…,N，N<J；

第三步：进行PCA消噪，得到N组主成分数连续的消噪数据；设U＝[u₁,...,u_M]∈C^J×M包含X协方差矩阵特征分解前M个最大特征值对应的特征矢量，Y＝U^HX∈C^M×V为PCA消噪数据，含有M个主成分；

第四步：计算N组PCA消噪数据的DOI；设置循环M＝1,…,N，分别计算各组消噪矩阵Y的DOI如下：

$\mathrm{\&rho;}=\frac{\det \left({\tilde{R}}_{YY}{R}_{YY}^{-*}{\tilde{R}}_{YY}^{*}\right)}{{\mathrm{detR}}_{YY}}---\left(1\right)$

式中，det表示求行列式值，R_YY＝E(YY^H)，E表示数学期望，上标H表示共轭转置，上标T表示转置，上标*表示共轭，上标-表示求逆；得到N组PCA消噪数据的DOI序列ρ₁,...,ρ_N；

第五步：检测和调整DOI曲线；从前至后逐点检测DOI值，当ρ_i>ρ_i-1,i＝2,...,N时，令ρ_i＝ρ_i-1；

第六步：计算SORTE序列：

$SORTE\left(p\right)=\left\{\begin{array}{cc}\frac{\mathrm{var}\&lsqb;{\left\{\&dtri;{\mathrm{\&rho;}}_i\right\}}_{i=p+1}^{N-1}\&rsqb;}{\mathrm{var}\&lsqb;{\left\{\&dtri;{\mathrm{\&rho;}}_i\right\}}_{i=p}^{N-1}\&rsqb;},& \mathrm{var}\&lsqb;{\left\{\&dtri;{\mathrm{\&rho;}}_i\right\}}_{i=p}^{N-1}\&rsqb;>0\\ +\mathrm{\&infin;},& \mathrm{var}\&lsqb;{\left\{\&dtri;{\mathrm{\&rho;}}_i\right\}}_{i=p}^{N-1}\&rsqb;=0\end{array}\right.---\left(2\right)$

其中，i＝1,...,(N-1)为调整后DOI曲线的差值序列，p＝1,...,(N-2)，

$\mathrm{var}\&lsqb;{\{\&dtri;{\mathrm{\&rho;}}_i\}}_{i=p}^{N-1}\&rsqb;=\frac{1}{N-p}{\mathrm{\&Sigma;}}_{i=p}^{N-1}{(\&dtri;{\mathrm{\&rho;}}_i-\frac{1}{N-p}{\mathrm{\&Sigma;}}_{i=p}^{N-1}\&dtri;{\mathrm{\&rho;}}_i)}^2$

第七步：计算模型阶数；根据SORTE准则，计算模型阶数P如下：

$P=\arg \underset{p=1,...,(N-3)}{min}SORTE\left(p\right)---\left(3\right)$

第八步：输出模型阶数估计结果P。

2.根据权利要求1所述的一种利用PCA和非环形特性估计复数fMRI数据模型阶数的方法，最大主成分数N取为幅值fMRI数据所估计模型阶数的三倍。