有监督的多模态脑影像融合方法

摘要

本发明公开了一种有监督的多模态脑影像融合方法。其中，该方法包括步骤1：计算所述各模态的特征；步骤2：对所述各模态的特征进行矩阵化和归一化；步骤3：利用奇异值分解算法，对所述各模态特征进行降维；步骤4：基于所述步骤S3得到的降维后的各模态特征，同时最大化所述各模态典型变量之间以及所述典型变量与先验信息之间的相关性平方和，并进行迭代循环，直至收敛；步骤5：将所述步骤S4得到的各模态成分串联，使用联合独立成分分析算法，计算得到所述各模态与所述先验信息显著相关的独立成分和混合矩阵，从而实现有监督的多模态脑影像融合。本发明实施例具有较好的鲁棒性，可以揭示复杂脑疾病认知功能损伤的生理病理机制。

说明书

技术领域

本发明实施例涉及医学图像处理技术领域，具体涉及一种有监督的多模态脑影像融合方法。

背景技术

近年来，使用多种非侵入式的图像技术(功能磁共振影像fMRI、结构磁共振影像sMRI、弥散张量成像DTI等)对同一个被试收集不同模态的数据已经被研究人员广泛采用。模态就是一台磁共振机器可以实现的多种成像。每一种模态从不同角度反映了大脑的功能或者结构。如fMRI基于血氧水平依赖性(blood oxygen level dependent，BOLD)信号，来反映大脑在做某种任务或者静息态时，对应脑区的神经元活动；sMRI则提供大脑的组织结构信息：如灰质(gray matter，GM)、白质(white matter，WM)、脑脊液(cerebrospinal fluid，CSF)等；弥散张量成像DTI可以研究脑白质的完整性和连接性。单独使用某一种模态的图像都不能同时获得这几种模态的图像的完整信息，而在实际应用中它们又都无法代替对方。越来越多的研究表明，多模态融合(multimodal fusion)能够结合多种模态的脑影像从多个角度理解人脑，并且具有一定的互补性，可以发现脑疾病的多模态影像共变模式[1](Sui J.,Huster R.,and Yu Q.B.et al.Function–structure associations of the brain:Evidence from multimodal connectivity and covariance studies.NeuroImage,102:11–23,2014)。

现有的多模态融合方法分为：基于数据(first-level)的多模态融合方法，如：独立成分分析(ICA，independent component analysis)[2-3](Hyvarinen A.,and Oja E.Independent component analysis:algorithms and applications.Neural Networks,13:411–430,2000；Du Y.,Fan Y.Group information guided ICA for fMRI data analysis.NeuroImage,69:157-197,2013)，该方法无法获得多模态共变成分，只能单独分析各个模态。基于特征(second-level)的多模态融合方法，如：联合独立成分分析(joint ICA)[4](Calhoun V.D.,Adali T.,and Kiehl K.A.et al.A Method for Multitask fMRI Data Fusion Applied to Schizophrenia，Human Brain Mapping,27:598–610,2006)，该方法得到的各模态的混淆矩阵是相同的。MCCA+jICA[4](Sui J.,He H.,and Pearlson G.D.et al.Three-way(N-way)fusion of brain imaging data based on mCCA+jICA and its application to discriminating schizophrenia.NeuroImage,2:119–132,2013),Parallel ICA[5](pICA)(Liu J.,Pearlson G.,and Windemuth A.et al.Combining fMRI and SNP data to investigate connections between brain function and genetics using parallel ICA.Human Brain Mapping,30:241-55,2009),该方法同时保证成分的独立性以及不同模态之间的对应性，然而以上方法都是基于无监督的多模态融合。pICA-R(pICA with references)[6](Chen J.,Calhoun V.D.,and Pearlson G.D.et al.Guided exploration of genomic risk for gray matter abnormalities in schizophrenia using parallel independent component analysis with reference.Neuroimage,83:384–396,2013)，利用基因作为先验信息指导多模态的融合，然而该方法只能融合两种模态，其中必须包含基因，并且融合的效果严重依赖于基因的精确度。

另一方面，利用脑影像数据来识别病人的认知损伤或症状改变的神经结构基础是一项极具临床价值的研究课题，在理解精神疾病的发病机制上具有重要意义。在临床医学中，精神分裂症(schizophrenia)是一种常见、慢性、高致残性的重型精神障碍。多起病于青壮年，常有感知、思维、情感、行为等多方面的障碍和精神活动的不协调，也有包括注意、记忆和执行功能障碍等在内的认知功能缺陷，患病率约1％。目前，我国重性精神疾病患者约为1600万，而罹患精神分裂症者多达780万。精神分裂症患者的社会功能降低，给患者、家庭及国家带来沉重的经济负担。因此如何利用现代医学成像技术，研究精神分裂症的病理机制，并对精神分裂症做出早期预警和干预具有重要意义。值得注意的是，认知功能缺失被认为是精神分裂症患者的核心特征，精神分裂症患者在认知方面的功能状态，相较于其他精神疾病的诊断，有明显的落差。虽然许多精神疾病都会出现认知功能缺失，但是精神分裂症患者的认知功能缺失，较其他类别诊断更为严重。随着药物科技的不断精进，患者的症状可以快速受到控制，让患者出院的机率大幅提高。然而，精神分裂症患者往往不能重新回归社会和适应日常生活，这主要归因于其认知功能的缺失是持久的，而认知功能缺失的严重程度，恰恰影响了精神分裂症患者心理复健的成效和适应社会生活的程度。认知能力包含多个认知子领域，比如：处理速度、注意力、工作记忆、语言学习、视觉学习、推理能力和社会认知等等。目前针对患者的症状评分或者特定认知领域评分的研究很大一部分是临床评分和药物测评，与人脑磁共振影像尤其是多模态磁共振影像相关的研究非常有限，大都是基于独立的分析然后求解相关关系，而基于有监督学习的、探索多模态神经影像(包括磁共振影像、脑电图、脑磁图)与特定临床指标相互关系的研究尚未展开。

有鉴于此，特提出本发明。

发明内容

鉴于上述问题，提出了本发明以便提供一种克服上述问题或者至少部分地解决上述问题的一种有监督的多模态脑影像融合方法。

为了实现上述目的，根据本发明的一个方面，提供了以下技术方案：

一种有监督的多模态脑影像融合方法，所述方法至少包括：

步骤1：计算所述各模态的特征；

步骤2：对所述各模态的特征进行矩阵化和归一化；

步骤3：利用奇异值分解算法，对所述各模态特征进行降维；

步骤4：基于所述步骤S3得到的降维后的各模态特征，同时最大化所述各模态典型变量之间以及所述典型变量与先验信息之间的相关性平方和，并进行迭代循环，直至收敛；

步骤5：将所述步骤S4得到的各模态成分串联，使用联合独立成分分析算法，计算得到所述各模态与所述先验信息显著相关的独立成分和混合矩阵，从而实现有监督的多模态脑影像融合。

与现有技术相比，上述技术方案至少具有以下有益效果：

本发明实施例通过加入先验信息(如认知评分)来指导多模态融合，同时最大化各模态典型变量之间以及典型变量与先验信息reference之间的相关性平方和，并利用jICA算法，得到各模态与reference显著相关的成分。不同于目前大多数无监督多模态融合，本发明实施例采用认知评分作为先验信息，能够有目的地提取出与特定兴趣指标显著相关的独立成分。一旦确定某一认知评分作为先验引导，便能够有监督的、稳定地提取出多中心数据的多模态共变的目标成分，从而挖掘出认知与多模态神经影像在精神疾病中的交互关系，具有较好的鲁棒性。另外，本发明实施例能够将认知映射到多模态神经影像中，从而揭示复杂脑疾病认知功能损伤的生理病理机制。

当然，实施本发明的任一产品不一定需要同时实现以上所述的所有优点。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其它优点可通过在所写的说明书、权利要求书以及附图中所特别指出的方法来实现和获得。

附图说明

附图作为本发明的一部分，用来提供对本发明的进一步的理解，本发明的示意性实施例及其说明用于解释本发明，但不构成对本发明的不当限定。显然，下面描述中的附图仅仅是一些实施例，对于本领域普通技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。在附图中：

图1为根据一示例性实施例示出的有监督的多模态脑影像融合方法的流程示意图；

图2a为根据另一示例性实施例示出的maxvar、ssqcor、MCCAR、MCCA+jICA、MCCAR+jICA相对于成分个数求平均的混合矩阵的精确度的示意图；

图2b为根据另一示例性实施例示出的maxvar、ssqcor、MCCAR、MCCA+jICA、MCCAR+jICA相对于成分个数求平均的独立成分的精确度的示意图；

图2c为根据另一示例性实施例示出的maxvar、ssqcor、MCCAR、MCCA+jICA、MCCAR+jICA相对于16个噪声水平求平均的混合矩阵的精确度的示意图；

图2d为根据另一示例性实施例示出的maxvar、ssqcor、MCCAR、MCCA+jICA、MCCAR+jICA相对于16个噪声水平求平均的独立成分的精确度的示意图；

图3a为根据一示例性实施例示出的maxvar、ssqcor、MCCAR、MCCA+jICA、和MCCAR+jICA方法在峰值信噪比为7的情况下提取出与先验信息显著相关的联合独立成分能力的对比结果示意图；

图3b为根据一示例性实施例示出的maxvar、ssqcor、MCCAR、MCCA+jICA、MCCAR+jICA方法对各模态的目标成分的混合矩阵的精确度的对比结果示意图；

图3c为根据一示例性实施例示出的maxvar、ssqcor、MCCAR、MCCA+jICA、MCCAR+jICA方法对各模态的目标成分的精确度的对比结果示意图；

图3d为根据一示例性实施例示出的maxvar、ssqcor、MCCAR、MCCA+jICA、MCCAR+jICA和3ICA方法估计得到的各模态目标成分之间的相关性对比示意图；

图3e为根据一示例性实施例示出的maxvar、ssqcor、MCCAR、MCCA+jICA、MCCAR+jICA和3ICA方法估计得到的各模态目标成分之间的相关性与各模态目标成分之间真实的相关性差的绝对值的对比示意图；

图4a为根据一示例性实施例示出的本发明实施例提出的方法检测到的与CMINDS认知评分相关的多模态共变联合成分的空间视图；

图4b为根据一示例性实施例示出的各模态的组间差异图；

图4c为根据一示例性实施例示出的各模态与认知评分的相关性拟合图。

这些附图和文字描述并不旨在以任何方式限制本发明的构思范围，而是通过参考特定实施例为本领域技术人员说明本发明的概念。

具体实施方式

下面结合附图以及具体实施例对本发明实施例解决的技术问题、所采用的技术方案以及实现的技术效果进行清楚、完整的描述。显然，所描述的实施例仅仅是本申请的一部分实施例，并不是全部实施例。基于本申请中的实施例，本领域普通技术人员在不付出创造性劳动的前提下，所获的所有其它等同或明显变型的实施例均落在本发明的保护范围内。本发明实施例可以按照权利要求中限定和涵盖的多种不同方式来具体化。

需要说明的是，在下面的描述中，为了方便理解，给出了许多具体细节。但是很明显，本发明的实现可以没有这些具体细节。

需要说明的是，在没有明确限定或不冲突的情况下，本发明中的各个实施例及其中的技术特征可以相互组合而形成技术方案。

图1为本发明有监督的多模态脑影像融合方法的流程图，如图1所示，该方法可以包括步骤S1至步骤S5。

S1：计算各模态的特征。

在本步骤中，模态可以是静息态功能脑影像fMRI、弥散张量图像dMRI、结构磁共振影像sMRI。在实际实施过程中，可以选取其中的任意两个或者其他更多的模态进行实施。

本实施例中仅以精神分裂症的多模态数据融合为示例进行说明，但是本发明的保护范围并不仅限于此。各模态的数据可以是所有精神障碍疾病的多模态数据，包括精神分裂症，抑郁症等精神疾病，以及双性情感障碍等非精神疾病的多模态数据。

对于fMRI，可以进行如下处理：(1)头动校正；(2)层间时间校正；(3)标准化到蒙特利尔标准空间，优选地，重采样成3×3×3mm；(4)回归掉6个头动参数和白质(WM)、灰质(GM)、脑脊液(CSF)；(5)用8mm full width half max(FWHM)Gaussian filter(8毫米半最大值全宽高斯滤波核)做平滑；(6)计算分数低频波动振幅(fALFF)特征。

对于dMRI，(1)质量检测，其包括去除梯度方向大于头动的方向以及振动的伪影；(2)运动和涡流校正；(3)校正梯度方向；(4)计算分数各项异性(FA)；(5)用8毫米半最大值全宽高斯滤波核做平滑。

对于sMRI，(1)标准化到蒙特利尔标准空间，优选地，重采样成3×3×3mm；(2)分割成白质、灰质、脑脊液；(3)用8毫米半最大值全宽高斯滤波核对灰质做平滑，(4)检测被试野点，以确保所有被试被正确分割。

经过上述处理后，各模态的特征为3D图像(也即3D矩阵)。对于静息态功能脑影像，模态特征可以为分数低频波动振幅特征。对于弥散张量图像，模态特征可以为分数各项异性特征。对于结构磁共振影像，模态特征可以为灰质密度特征。

S2：对各模态的特征进行矩阵化和归一化。

具体地，矩阵化的处理过程可以为：将经过步骤S1处理后得到的3D矩阵转换成行向量(1×L)，L为体素个数。将所有被试的3D矩阵都进行转换，则得到N×L的矩阵X_k(k＝1,2,3)，N表示被试个数，k表示模态数目。

对3个模态dMRI、dMRI和sMRI都进行上述转换，则得到3个模态的特征矩阵。

接着，对上述特征矩阵进行归一化，使其含有相同的平均平方和，以确保3个模态的特征数值在相同的范围内。

S3：利用奇异值分解算法，对各模态特征进行降维。

具体地，根据以下公式对各模态特征进行降维：

Y_k＝X_kE_k

其中，Y_k表示降维后的特征矩阵；X_k表示被试个数乘以体素个数(N×L)的特征矩阵，N表示被试个数，L表示体素个数；E_k表示包含了前M个最大的奇异值对应的特征向量，M取正整数。

S4：基于步骤S3得到的降维后的各模态特征，同时最大化各模态典型变量之间以及典型变量与先验信息之间的相关性平方和，并进行迭代循环，直至收敛。

其中，先验信息用于指导多模态融合，其包括但不限于认知评分、症状评分、某一位点的基因。先验信息可以根据实际研究目的进行选取。

具体地，根据以下公式进行同时最大化所述各模态典型变量之间以及所述典型变量与先验信息之间的相关性平方和：

$\{w_1^{\left(1\right)},w_2^{\left(1\right)},\mathrm{...},w_n^{\left(1\right)}\}=\arg \underset{w}{\max }\{{\Sigma }_{k,j=1}^n{\left|r_{k,j}^{\left(1\right)}\right|}^2+\lambda ·{\Sigma }_{k=1}^n{\left|r_{k,ref}^{\left(1\right)}\right|}^2\}---\left(1\right)$

$\{w_1^{\left(i\right)},w_2^{\left(i\right)},\mathrm{...},w_n^{\left(i\right)}\}=\arg \underset{w}{\max }\{{\Sigma }_{k,j=1}^n{\left|r_{k,j}^{\left(i\right)}\right|}^2+\lambda ·{\Sigma }_{k=1}^n{\left|r_{k,ref}^{\left(i\right)}\right|}^2\}---\left(2\right)$

$\left(\begin{array}{ccc}s.t.& w_k^{\left(i\right)}\perp \{w_k^{\left(1\right)},w_k^{\left(2\right)},\mathrm{...},w_k^{(i-1)}\};& k,j=1,2,\mathrm{...},n\end{array}\right)$

其中，在公式(1)中，i取1，在公式(2)中，i＝2......M；M表示融合的成分个数，优选地，M可以参考最小描述长度准则(minimum description length criterion)进行选取，例如，本实施例中M取20；n表示模态个数；是矩阵w的第i列；表示A_k与A_j的第i列的相关性；表示A_k与先验信息的第i列的相关性；A_k和A_j表示不同模态的典型变量；λ表示控制和之间的权重的正则化参数以及算法优化的收敛性。

求解能量函数公式(1)和公式(2)的最大化问题具体可以包括：公式(1)为第一阶段的求解，将相关性平方和(SSQCOR)作为目标函数，通过将目标函数相对于每个求偏导数，并令之等于0，得到含有的偏微分方程。由于目标函数相关性平方和SSQCOR是的二次方程，因此相关性平方和SSQCOR对的偏微分是的线性函数，所以可以求得近似解。给定一个初始点，每一个向量将通过公式(1)逐步迭代被更新，并且保证目标函数SSQCOR的值不断增加，当满足目标函数的收敛准则时，迭代停止，此时的作为公式(1)的最优解。公式(2)为第二阶段的求解。公式(2)以及更高阶段的求解，通过将典型相关系数的正交约束作为拉格朗日乘子加入到目标函数中，使用类似于公式(1)的求解方式求得稳定解，即公式(2)的最优解。

基于上述最大化算法，同时可得到各模态的典型变量A_k，其满足：

$E\left\{A_k^T{A}_k\right\}=I,E\left\{A_k^T{A}_k\right\}\approx diag(r_{k,j}^{\left(1\right)},r_{k,j}^{\left(2\right)},...,r_{k,j}^{\left(M\right)}),k\ne j,k,j\in \{1,2,3\}$

其中，表示典型相关系数，此时就可以得到各模态整个第一阶段的典型变量A_k和成分C_k。

S5：将步骤S4得到的各模态成分串联，使用联合独立成分分析算法，计算得到各模态与先验信息显著相关的独立成分和混合矩阵，从而实现有监督的多模态脑影像融合。

具体地，根据以下公式得到独立成分和混合矩阵：

W[C₁,C₂,…,C_n]＝[S₁,S₂,…,S_n]

X_k＝(A_k·W^-1)·S_k,k＝1,2,…,n

其中，[C₁,C₂,…,C_n]为步骤S4得到各模态成分的串联，S_k为最终求解得的各模态独立成分，A_k·W^-1为混合矩阵。对于jICA的求解，这里采用目前广泛使用的Infomax(信息极大)算法。

本实施例中虽然将各个步骤按照上述先后次序的方式进行了描述，但是本领域技术人员可以理解，为了实现本实施例的效果，不同的步骤之间不必按照这样的次序执行，其可以同时(并行)执行或以颠倒的次序执行，这些简单的变化都在本发明的保护范围之内。

图2a、图2b、图2c、图2d对比了各方法分离得到的混合矩阵A_k·W^-1和独立成分S_k的精确度。

在本文中，各方法为：maxvar是指多组别典型相关分析-最大特征值；ssqcor是指多组别典型相关分析-相关性平方和；MCCAR是指有监督的多组别典型相关分析；MCCA+jICA是指多组别典型相关分析+联合独立成分分析；MCCAR+jICA是指有监督的多组别典型相关分析+联合独立成分分析；3ICA是指单个模态的独立成分分析。另外，real corr是指真实的相关性。

精确度是指各方法得到的混合矩阵、独立成分与真实的混合矩阵、成分的相关性，两者的相关性越高则表明精确度越高。

图2a示例性地示出了各方法相对于成分个数求平均的混合矩阵A_k·W^-1的精确度。图2b示例性地示出了各方法相对于成分个数求平均的独立成分S_k的精确度。图2c示例性地示出了各方法相对于16个噪声水平(峰值信噪比PSNR＝[1>k·W^-1的精确度。图2d示例性地示出了各方法中相对于16个噪声水平(峰值信噪比PSNR＝[1>k的精确度。通过图2c和图2d可以看出，本发明实施例提出的方法在混合矩阵A_k·W^-1和独立成分S_k的精确度上都是最优的。同时，也可以看出本发明实施例中添加jICA步骤后的精确度要明显高于MCCAR，这说明了使用联合独立成分分析jICA的重要性和必要性。

图3a示例性地示出了各方法在峰值信噪比PSNR＝7的情况下提取出与reference显著相关的联合独立成分能力的对比结果。其中，real corr为真实的相关性；如果相同索引的成分在不止一个模态中有组间差异，则这个成分称为多模态共变的联合独立成分。从图3a中可以看出只有MCCAR和MCCAR+jICA可以得到联合独立成分。图3b示例性地示出了各方法对各模态的目标成分的混合矩阵A_k·W^-1的精确度的对比结果。图3c示例性地示出了各方法对各模态的目标成分S_k的精确度的对比结果。这里，盒型图是对于峰值信噪比PSNR＝[134]统计的，可以看出本发明实施例提出的方法所分离出的目标成分的精确度是最高的。

图3d和图3e为各方法估计得到的各模态目标成分之间的相关性、各模态目标成分之间的相关性与各模态目标成分之间真实的相关性差的绝对值的对比图；其中，3ICA指separate ICA，即对3个模态分别做ICA；real表示各模态目标成分之间真实的相关性。从中可以看出本发明实施例提出的方法估计得到的目标成分之间的相关性与真实的相关性最接近。

下面对本发明示例里提出的方法进行的测试进行详细说明。

本发明实施例使用Function Biomedical Informatics Research Network(FBIRN)Phase III study数据集(包含147个精神分裂症患者和147个正常对照被试的3个模态数据)。利用这些数据对本发明实施例提出的方法进行了测试。这里，我们将总的认知评分作为先验信息。

图4a、图4b和图4c示例性地示出了本发明实施例在数据集FBIRN上的结果。其中，图4a是本发明实施例提出的方法检测到的与CMINDS认知评分相关的多模态共变成分的空间视图，可以看出其包含的脑区位置以及激活程度与现有的研究基本一致；图4b是各模态的组间差异图；图4c是各模态与认知评分CMINDS的相关性拟合图。综上所述，可以看出本发明提出的方法可以分离出既具有组检差异又与CMINDS认知评分显著相关的联合独立成分。

以上对本发明实施例所提供的技术方案进行了详细的介绍。虽然本文应用了具体的个例对本发明的原理和实施方式进行了阐述，但是，上述实施例的说明仅适用于帮助理解本发明实施例的原理；同时，对于本领域技术人员来说，依据本发明实施例，在具体实施方式以及应用范围之内均会做出改变。

需要说明的是，本文中涉及到的流程图不仅仅局限于本文所示的形式，其还可以进行划分和/或组合。

需要说明的是：附图中的标记和文字只是为了更清楚地说明本发明，不视为对本发明保护范围的不当限定。

术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法、物品或者设备/装置所固有的要素。

本发明的各个步骤可以用通用的计算装置来实现，例如，它们可以集中在单个的计算装置上，例如：个人计算机、服务器计算机、手持设备或便携式设备、平板型设备或者多处理器装置，也可以分布在多个计算装置所组成的网络上，它们可以以不同于此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。因此，本发明不限于任何特定的硬件和软件或者其结合。

本发明提供的方法可以使用可编程逻辑器件来实现，也可以实施为计算机程序软件或程序模块(其包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件或数据结构等等)，例如根据本发明的实施例可以是一种计算机程序产品，运行该计算机程序产品使计算机执行用于所示范的方法。所述计算机程序产品包括计算机可读存储介质，该介质上包含计算机程序逻辑或代码部分，用于实现所述方法。所述计算机可读存储介质可以是被安装在计算机中的内置介质或者可以从计算机主体上拆卸下来的可移动介质(例如：采用热插拔技术的存储设备)。所述内置介质包括但不限于可重写的非易失性存储器，例如：RAM、ROM、快闪存储器和硬盘。所述可移动介质包括但不限于：光存储介质(例如：CD－ROM和DVD)、磁光存储介质(例如：MO)、磁存储介质(例如：磁带或移动硬盘)、具有内置的可重写非易失性存储器的媒体(例如：存储卡)和具有内置ROM的媒体(例如：ROM盒)。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。