一种适于多被试复数fMRI数据分析的自适应定点IVA算法

摘要

一种适于多被试复数fMRI数据分析的自适应定点IVA算法，属于生物医学信号处理领域。采用基于多维广义高斯分布(MGGD)的非线性函数估计复数fMRI数据的源向量成分(SCV)分布；采用最大似然法自适应地估计MGGD的形状参数，与变化SCV分布自动匹配；在SCV主导子空间更新基于MGGD的非线性函数，实现对复数fMRI数据的消噪；在算法更新过程中加入输入数据的伪协方差阵，直接利用复数fMRI数据的非环形特性，进一步提高IVA分析复数fMRI数据的针对性。本发明能够有效分析高噪声水平但脑功能信息最为全面的多被试复数fMRI数据，在被试间差异性大且信噪比低的不利情况下，能为脑功能研究和脑疾病诊断提供更好的依据。

主权项

一种适于多被试复数fMRI数据分析的自适应定点IVA算法，采用基于多维广义高斯分布MGGD的非线性函数估计复数fMRI数据的SCV分布；采用最大似然估计MLE方法自适应地估计MGGD的形状参数，与不断变化的SCV分布自动匹配；在SCV主导子空间更新基于MGGD的非线性函数，实现对复数fMRI数据的消噪；在算法更新过程中加入输入数据的伪协方差阵，直接利用复数fMRI数据的非环形特性，进一步提高IVA对复数fMRI数据分析的针对性；包括以下步骤：第一步：输入多被试复数fMRI数据K表示被试数目；J表示时间维的全脑扫描次数；M表示空间维的脑内体素数目；第二步：对各被试复数fMRI数据X^(k)分别进行PCA压缩和白化；设各被试的SM和TC的成分数为N，将被试k的复数fMRI数据X^(k)压缩并白化为为压缩阵，N＜J，将时间维由J降为N，为白化阵，使得的方差为1；第三步：初始化；随机初始化解混矩阵作用于计算得到第n个SCV成分的初始值是W^(k)的第n列，n＝1,…,N，是的第m列，m＝1,…,M，上标T表示转置，上标H表示共轭转置；忽略m，将y_n(m)、和x^(k)(m)简写为y_n、和x^(k)；利用式(1)计算自适应定点IVA算法的代价函数初始值其中，E表示数学期望，p(y_n)是y_n的概率密度函数；G(·)为实值非线性函数；|·|表示取复数的模值；设λ_n是协方差矩阵的最大特征值，为λ_n对应的特征向量，令式(1)中G(·)采用式(2)所示基于MGGD的非线性函数：$G\left(q_n\right)={\left(q_n\right)}^{{\beta }_n}={\left[{\lambda }_n{\left({\Sigma }_{k=1}^K{v}_{kn}\right|y_n^{\left(k\right)}\left|\right)}^2\right]}^{{\beta }_n}---\left(2\right)$式中，β_n为MGGD的形状参数，也是各SCV分布的形状参数，初始化β_n＝β₀；第四步：更新解混矩阵W^(k)；对W^(k)的每一列分别按照式(3)进行更新：$\begin{array}{c}{w}_n^{\left(k\right)}\leftarrow -E\left\{y_n^{\left(k\right)*}{G}^{\prime }\left(q_n\right)x^{\left(k\right)}\right\}\\ +E\{G^{\prime }\left(q_n\right)+{\left|y_n^{\left(k\right)}\right|}^2{G}^{\prime \prime }\left(q_n\right)\}{w}_n^{\left(k\right)}\\ +E\left\{x^{\left(k\right)}{\left(x^{\left(k\right)}\right)}^T\right\}E\left\{{\left(y_n^{\left(k\right)*}\right)}^2{G}^{\prime \prime }\left(q_n\right)\right\}{w}_n^{\left(k\right)*}\end{array}---\left(3\right)$上标*表示共轭，G′(·)和G″(·)分别为G(·)的一阶导数和二阶导数：$G^{\prime }\left(q_n\right)={\beta }_n{\left(q_n\right)}^{{\beta }_n-1}---\left(4\right)$$G^{\prime \prime }\left(q_n\right)={\beta }_n({\beta }_n-1){\left(q_n\right)}^{{\beta }_n-2}---\left(5\right)$W^(k)各列全部更新之后，对W^(k)再进行式(6)所示去相关操作：W^(k)←(W^(k)(W^(k))^H)^‑1/2W^(k)    (6)第五步：采用基于Newton‑Raphson优化的最大似然估计方法，估计并更新各SCV分布的形状参数β_n：${\beta }_n\leftarrow {\beta }_n-\frac{\log L({\beta }_n;y_n(1),...,y_n(M\left)\right)}{\partial \log L({\beta }_n;y_n(1),...,y_n(M\left)\right)/\partial {\beta }_n}---\left(7\right)$其中，β_n的最大似然估计为式中，Γ(·)为伽马函数；Σ_n的逆阵为第六步：计算本次迭代的代价函数；记本次迭代次数为iter，采用式(1)计算本次迭代的代价函数第七步：迭代终止条件判断；计算本次迭代代价函数与上次迭代代价函数之差当的绝对值小于预设阈值ε或者达到最大迭代次数iter_max时，自适应定点IVA算法结束迭代，否则返回第四步；第八步：计算各被试的N个SM成分和N个TC成分如下：$Y^{\left(k\right)}=W{\tilde{X}}^{\left(k\right)}---\left(10\right)$$R^{\left(k\right)}={\left(U_2^{\left(k\right)}{U}_1^{\left(k\right)}\right)}^{-1}{\left(W^{\left(k\right)}\right)}^H---\left(11\right)$其中，第九步：进行相位校正与SM相位消噪；首先应用各TC成分进行相位校正，然后去除对应SM成分中相位范围在‑π/4～π/4之外的体素，得到相位消噪的SM成分；第十步：输出各被试相位消噪后的N个SM成分以及相位校正后的N个TC成分。