技术领域
[0001] 本
发明属于数字图像领域，尤其涉及一种基于子空间特征学习的丘脑功能分区方法。
背景技术
[0002] 近年来，医学成像技术为无创地构建丘脑功能分区提供了新的途径。
弥散张量成像基于测量大脑组织
水分子弥散性质得到组织信息，并且可以重建出
纤维束，表征不同区域之间的结构连接。弥散是指分子的随机性不规则运动，基于分子在空间中沿特定方向扩散的程度情况，又分为各向同性和
各向异性。各向同性指分子的运动向各个方向的机会是均等的，例如在纯水中，水分子的弥散为各向同性弥散以及脑灰质中水分子的弥散也近似各向同性。各向异性指分子的弥散具有方向性，如脑白质中，水分子的弥散方向总与纤维束走行方向一致。
[0003] 目前，学者们则更倾向于基于结构连接特征生成丘脑核团分区。Behrens等借助弥散张量成像中水分子的弥散特性检测出具有方向性的脑组织，同时使用概率纤维束成像的方法获得丘脑的结构连接模式，采用少数皮层功能分区定义丘脑
体素的结构连接特征，生成丘脑核团分区。这种方法获得的信息损失严重、特征区分度低，只能获得粗糙的丘脑核团分区结果。与rsfMRI的聚类方法类似，最近几年则通过全脑其他区域体素描述丘脑体素的结构连接特征。伦敦国王学院O'Muircheartaigh与Toulmin等学者，在获得丘脑体素到全脑其他区域的概率
纤维跟踪结果之后，使用独立成分分析聚类，生成个体的丘脑核团分区。伦敦大学圣乔治医学院Lambert等学者使用与相同的特征，计算丘脑体素之间的距离特征图，生成个体的丘脑核团分区。然而利用全脑其他体素计算丘脑体素的连接特征数据维度过高，而且容易受到噪声的影响。
发明内容
[0004] 发明目的：由于使用全脑体素作为丘脑特征会得到极高维度的特征，容易受到噪声影响，而使用解剖结构定义的粗糙皮层分区则使得特征区分度较低。同时丘脑皮层间的结构连接特征是复杂非线性的，常用的
降维学习手段不能有效学习出体素的特征表达。针对以上问题，本发明提出一种基于子空间特征学习的丘脑功能分区方法，可以高效地获得丘脑功能分区。
[0005] 技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种基于子空间特征学习的丘脑功能分区方法，包括以下步骤：
[0006] 步骤1，对弥散张量图像原始数据预处理，提取图像特征，获得丘脑-皮层结构连接特征；
[0007] 步骤2，将丘脑-皮层结构连接特征作为输入数据，搭建深度子空间网络模型，训练网络模型得到输入数据的低维子空间特征以及自表达特征；所述低维子空间特征即为输入数据的潜在特征；
[0008] 步骤3，使用无监督聚类方法，对步骤2中得到的输入数据的潜在特征以及自表达特征进行聚类，并对聚类结果进行分割，形成个体丘脑功能分区。
[0009] 通过对比实验，使用不同方法对个体丘脑进行功能分区，并组合起来形成群组图谱，检验效果。
[0010] 进一步，步骤1所述对弥散张量图像原始数据预处理，提取图像特征，获得丘脑-皮层结构连接特征；方法如下：
[0011] 1-1，采用纤维束示踪成像技术获得原始图像数据的纤维概率结果；
[0012] 1-2，采用精细皮层分区方法提取出体素的成像特征，减小弥散张量图像数据和纤维追踪的噪声；
[0013] 1-3，根据步骤1-1和步骤1-2，获得丘脑-皮层结构连接特征。
[0014] 进一步，步骤1-2所述精细皮层分区方法，即考虑精细度与噪声，使用Shen Atlas脑图谱作为先验脑图谱，提取出体素的成像特征。所述Shen Atlas脑图谱将脑区分为不同数目的分区，特征维度为268。
[0015] 进一步，步骤2所述搭建深度子空间网络模型，训练网络模型得到输入数据的低维子空间特征以及自表达特征；步骤如下：
[0016] 2-1，将丘脑-皮层结构连接特征作为
特征向量输入，搭建子空间自
编码器网络模型用于提取低维子空间特征。子空间自编码器网络模型包含三个部分，即编码结构、自表达层、解码结构。
[0017] 在编码结构中学习特征的潜在表示，编码结构包含两个全连接层，每个全连接层后使用了参数化的线性修正单元，对输入数据进行非线性激活。
[0018] 为了避免用于编码的神经网络参数过于稀疏化，导致学习到的信息少而不能有效表达原始特征的子空间结构。本发明使用改进的带参数的线性修正单元，避免了大范围的参数抑制，利于丘脑体素结构连接特征学习。
[0019] 带参数的线性修正单元定义为：
[0020]
[0021] 其中yi代表非线性激活函数f在第i个通道的数据输入，ai是第i个通道的控制输入数据非负部分的斜率；当ai为0时，上式为线性修正单元。
[0022] 在解码结构将特征的潜在表示进行回映射，得到解码输出数据；在编码结构和解码结构之间加入自表达层，自动学习输入数据之间的相关性。
[0023] 对于给定的图像中的体素数据点集合{xi}i＝1,2,...,M，M代表体素总数，这些体素数据点分布在K个线性子空间{Si}i＝1,2,...,K中，如果属于某个子空间中的点能够被同一子空间中的其它点通过线性组合的方法表达，这种性质称为自表达属性。
[0024] 将这一系列体素数据点组成数据矩阵X，有如下公式：
[0025] X＝XC (1)
[0026] 其中C为系数矩阵，即自表达属性的相关性矩阵。
[0027] 有研究表明，在数据子空间互相独立的假设下，通过最小化矩阵C的特定范数就能使矩阵C保持分
块对
角结构，此时若数据点xi和xj属于同一子空间，那么相关性矩阵C的元素cij≠0，因此能使用矩阵C构建体素数据点的亲和度矩阵用于聚类处理。
[0028] 加入特定范式约束C后的自表达公式表示为：
[0029]
[0030] 上式‖·‖p表示任一矩阵范数，根据对矩阵C的不同约束可以选择不同的范数，如稀疏约束的L1范数、低秩约束的核范数。
[0031] 在实际应用中，将数据的噪声考虑进去，公式改写为：
[0032]
[0033] 其中λ表示自表达损失的权重参数，控制自表达损失在公式(3)中的权重占比。
[0034] 2-2，定义子空间自编码器网络模型的损失函数，在深度子空间聚类(DSC)中为了通过子空间自编码器网络学习数据的自表达，引入全新的损失函数：
[0035]
[0036] 式中，X表示原始特征输入， 表示输入数据经过
解码器重构后的输出，Θ表示权值参数，C为系数矩阵，表示自表达层的权值，ZΘ表示编码阶段的特征输出；λ1、λ2为权重权衡系数；
[0037] 公式(4)的第一项表示数据的重建损失；第二项表示自表达层的正则项，通过最小化系数矩阵C的特定范数确保系数矩阵C维持分块对角结构，从而使得ZΘ具有可分割特性；第三项表示数据的自表达项，数据空间ZΘ中的每个体素点ZΘ,i通过系数矩阵C和其他数据点{ZΘ,j},j＝1,...,N,j≠i表示，N为数据空间ZΘ的点个数；使用全连接层表示自表达项，系数矩阵C即为自表达层的权值，通过系数矩阵C构建数据空间ZΘ的亲和度矩阵。
[0038] 通过优化上述损失函数，使得到的数据空间ZΘ低维表示输入数据，同时使得系数矩阵C达到自表达效果；使用反向传播联合求解网络模型中的权值参数Θ和系数矩阵C；所述系数矩阵C作为输入数据的自表达特征。
[0039] 2-3，采用两步训练策略，训练子空间自编码器网络模型：
[0040] 第一步，使用不添加自表达层的传统自编码器进行子空间自编码器网络模型预训练，同时模型的损失函数采用输入与输出均方误差函数，优化所述损失函数，使得网络初步学习丘脑特征的潜在表示，所述丘脑特征的潜在表示即为编码器的输出；
[0041] 第二步，使用第一步得到的预训练模型里的神经网络权重参数，初始化深度子空间自编码器网络模型结构中对应的编码与解码部分，将输入数据按批次投入训练，将同一丘脑的丘脑-皮层结构连接特征加入一个批次进行网络模型训练；
[0042] 同时，因左右丘脑具有差异性，将左右丘脑分开训练，对于单独的左丘脑或右丘脑，全部体素的特征向量在一个批次中进行训练，最小化深度子空间自编码器网络模型的损失函数，每个特征向量由编码器映射得到潜在表示，即提取低维子空间特征。
[0043] 本发明采用基于梯度下降的自适应矩估计方法优化子空间自编码器网络模型中的权重参数以及系数矩阵，所述权重参数指不属于自表达层的神经元参数。
[0044] 进一步，步骤3所述使用无监督聚类方法，对步骤2中得到的输入数据的潜在特征以及自表达特征进行聚类，并对聚类结果进行分割，形成个体丘脑功能分区；步骤如下：
[0045] 3-1，不同体素间的亲和度代表彼此间的特征相关性，依据此相关性可以找出体素间的差异性，构建体素的亲和度矩阵。
[0046] 为了解决基于自表达系数构建亲和度矩阵在个体水平上聚类结果差异大的问题，本发明使用基于特征潜在表示的亲和度矩阵。
[0047] 根据步骤2中得到的特征潜在表示作为构建体素亲和度矩阵的衡量数据；
[0048] 体素点i,j对应的潜在特征向量为ZΘ,i与ZΘ,j，体素点i,j之间的相似性由它们的潜在特征向量之间的欧式距离表示；公式如下：
[0049] Sij＝||ZΘ,i-ZΘ,j||2
[0050] 式中，Sij表示体素点i,j的相似性，Sij取值范围在0到1之间，越接近0表示体素点i,j相似度越高，越接近1表示体素点i,j差异性越大。
[0051] 3-2，将所有体素构成一个邻接矩阵，每个体素只与邻域内的体素相连，对该邻接矩阵，采用融合测地距离(即加入丘脑体素空间约束)以及体素间欧式距离的方法来衡量体素间的相似性，能够反映体素空间的拓扑结构，丰富对空间信息的提取；
[0052] 所述邻接矩阵为一个无向相似度图，邻接矩阵元素的定义为：
[0053]
[0054] 其中aij表示邻接矩阵的元素值，Sij表示对应体素i,j两点间的亲和度值，E为体素构成的邻接矩阵(无向相似度图)的边集，vi,vj表示体素i,j，属于邻接矩阵(无向相似度图)的
顶点集，ε为定义的邻域半径；
[0055] 所述测地距离通过下式计算：
[0056]
[0057] i,j表示邻接矩阵(无向相似度图)中的两个顶点，P(*)为权重函数，即特征的欧式距离，P(x)>0，γ(t)表示从顶点i到顶点j的所有路径的集合。
[0058] 3-3，使用归一化割方法对弥散张量图像进行分割，根据分割结果形成个体丘脑的分区。
[0059] 有益效果：与
现有技术相比，本发明的技术方案具有以下有益的技术效果：
[0060] (1)本发明采用精细皮层分区能获得连接，减小了弥散张量成像数据和纤维追踪的噪声，保证特征的高分区度，并且能够降低噪声的影响。
[0061] (2)本发明选择Shen Atlas作为先验脑图谱，考虑到精细度与噪声的权衡，在特征维度为268的时候能提取出体素的成像特征，并且含有较少的噪声，有利于后续处理。
[0062] (3)本发明提出了融合测地距离与欧式距离的距离度量方式度量体素相似度，更好地反映体素空间的拓扑结构，丰富对空间信息的提取，利于丘脑体素的有效聚类。
附图说明
[0064] 图2是深度子空间网络图。
具体实施方式
[0065] 下面结合附图和
实施例对本发明的技术方案作进一步的说明。
[0066] 本发明所述的一种基于子空间特征学习的丘脑功能分区方法，如图1所示，首先使用弥散张量成像进行纤维追踪以获得活体大脑内部结构连接信息；其次针对复杂非线性的丘脑皮层结构连接特征，本发明使用深度子空间网络学习特征的隐藏子空间映射(潜在表示)；最后，对体素特征加以空间约束，降低噪声的影响，更好地反映空间拓扑结构，丰富对空间信息的提取，利于丘脑体素的有效聚类。具体如下：
[0067] 步骤1，对弥散张量图像原始数据预处理，提取图像特征，获得丘脑-皮层结构连接特征；
[0068] 1-1，对原始数据进行预处理，采用纤维束示踪成像技术获得图像的纤维概率结果。
[0069] 1-2，采用精细皮层分区方法提取出体素的成像特征，从而减小弥散张量成像数据和纤维追踪的噪声，获得丰富的连接保证特征的高分区度。
[0070] 1-3，根据步骤1-1和步骤1-2，获得丘脑-皮层结构连接特征。
[0071] 步骤2，将丘脑-皮层结构连接特征作为输入数据，搭建深度子空间网络模型，训练网络模型得到输入数据的低维子空间特征以及自表达特征，所述低维子空间特征即为输入数据的潜在特征；具体包括：
[0072] 2-1，将丘脑-皮层结构连接特征作为特征向量输入，搭建子空间自编码器网络模型用于提取低维子空间特征；子空间自编码器网络模型包含三个部分，即编码结构、自表达层、解码结构。
[0073] 在编码部分学习特征的潜在表示，编码部分包含两个全连接层，每个全连接层后使用了参数化的线性修正单元，对输入数据进行非线性激活。
[0074] 带参数的线性修正单元定义为：
[0075]
[0076] 其中yi代表非线性激活函数f在第i个通道的数据输入，ai是第i个通道的控制输入数据非负部分的斜率。
[0077] 在解码结构将特征的潜在表示进行回映射，得到解码输出数据；在编码结构和解码结构之间加入自表达层，自动学习输入数据之间的相关性。
[0078] 对于给定的图像中的体素数据点集合{xi}i＝1,2,...,M，M代表体素总数，这些体素数据点分布在K个线性子空间{Si}i＝1,2,...,K中，如果属于某个子空间中的点能够被同一子空间中的其它点通过线性组合的方法表达，这种性质称为自表达属性。
[0079] 将这一系列体素数据点组成数据矩阵X，有如下公式：
[0080] X＝XC (1)
[0081] 其中C为系数矩阵，即自表达属性的相关性矩阵。
[0082] 有研究表明，在数据子空间互相独立的假设下，通过最小化矩阵C的特定范数就能使矩阵C保持分块对角结构，此时若数据点xi和xj属于同一子空间，那么相关性矩阵C的元素cij≠0，因此能使用矩阵C构建体素数据点的亲和度矩阵用于聚类处理。
[0083] 加入特定范式约束C后的自表达公式表示为：
[0084]
[0085] 上式‖·‖p表示任一矩阵范数，根据对矩阵C的不同约束可以选择不同的范数，如稀疏约束的L1范数、低秩约束的核范数。
[0086] 在实际应用中，将数据的噪声考虑进去，公式改写为：
[0087]
[0088] 其中λ表示自表达损失的权重参数，控制自表达损失在公式(3)中的权重占比。
[0089] 2-2，定义子空间自编码器网络模型的损失函数，在深度子空间聚类(DSC)中为了通过网络学习数据的自表达，引入全新的损失函数：
[0090]
[0091] 式中，X表示原始特征输入， 表示输入数据经过解码器重构后的输出，Θ表示权值参数，C为系数矩阵，表示自表达层的权值，ZΘ表示编码阶段的特征输出；λ1、λ2为权重权衡系数；
[0092] 公式(4)的第一项表示数据的重建损失；第二项表示自表达层的正则项，通过最小化系数矩阵C的特定范数确保系数矩阵C维持分块对角结构，从而使得ZΘ具有可分割特性；第三项表示数据的自表达项，数据空间ZΘ中的每个体素点ZΘ,i通过系数矩阵C和其他数据点{ZΘ,j},j＝1,...,N,j≠i表示，N为数据空间ZΘ的点个数；使用全连接层表示自表达项，系数矩阵C即为自表达层的权值，通过系数矩阵C构建数据空间ZΘ的亲和度矩阵。
[0093] 通过优化上述损失函数，使得到的数据空间ZΘ低维表示输入数据，同时使得系数矩阵C达到自表达效果；使用反向传播联合求解网络模型中的权值参数Θ和系数矩阵C；所述系数矩阵C作为输入数据的自表达特征。
[0094] 2-3，采用两步训练策略，训练子空间自编码器网络模型：
[0095] 第一步，使用不添加自表达层的传统自编码器进行子空间自编码器网络模型预训练，同时模型的损失函数采用输入与输出均方误差函数，优化所述损失函数，使得网络初步学习丘脑特征的潜在表示，所述丘脑特征的潜在表示即为编码器的输出；
[0096] 第二步，使用第一步得到的预训练模型里的神经网络权重参数，初始化深度子空间自编码器网络模型结构中对应的编码与解码部分，即预训练模型直接作为接下来子空间自编码器网络模型的初始状态；将输入数据按批次投入训练，将同一丘脑的丘脑-皮层结构连接特征加入一个批次进行网络模型训练；
[0097] 同时，因左右丘脑具有差异性，将左右丘脑分开训练，对于单独的左丘脑或右丘脑，全部体素的特征向量在一个批次中进行训练，最小化深度子空间自编码器网络模型的损失函数，每个特征向量由编码器映射得到潜在表示，即提取低维子空间特征。
[0098] 采用基于梯度下降的自适应矩估计方法优化子空间自编码器网络模型中的权重参数以及系数矩阵。
[0099] 步骤3，使用无监督聚类方法，对步骤2中得到的输入数据的潜在特征以及自表达特征进行聚类，并对聚类结果进行分割，形成个体丘脑功能分区；具体包括：
[0100] 3-1，不同体素间的亲和度代表彼此间的特征相关性，依据此相关性可以找出体素间的差异性，构建体素的亲和度矩阵。
[0101] 为了解决基于自表达系数构建亲和度矩阵在个体水平上聚类结果差异大的问题，本发明使用基于特征潜在表示的亲和度矩阵。
[0102] 在特征学习的过程中，利用深度子空间聚类模型进行丘脑体素特征非线性映射的学习，在网络收敛以后，提取编码器输出即原始特征在进行子空间约束后的潜在表示。
[0103] 根据特征的潜在表示可以构建体素亲和度矩阵。获取编码器输出ZΘ，体素点i,j对应的潜在特征向量为ZΘ,i与ZΘ,j，体素点i,j之间的相似性由它们的潜在特征向量之间的欧式距离表示；公式如下：
[0104] Sij＝||ZΘ,i-ZΘ,j||2
[0105] 式中，Sij表示体素点i,j的相似性，Sij取值范围在0到1之间，越接近0表示体素点i,j相似度越高，越接近1表示体素点i,j差异性越大。
[0106] 3-2，将所有体素构成一个邻接矩阵，每个体素只与邻域内的体素相连，对该邻接矩阵，采用融合测地距离(即加入丘脑体素空间约束)以及体素间欧式距离的方法来衡量体素间的相似性，能够反映体素空间的拓扑结构，丰富对空间信息的提取；
[0107] 所述邻接矩阵为一个无向相似度图，邻接矩阵元素的定义为：
[0108]
[0109] 其中aij表示邻接矩阵的元素值，Sij表示对应体素i,j两点间的亲和度值，E为体素构成的邻接矩阵(无向相似度图)的边集，vi,vj表示体素i,j，属于邻接矩阵(无向相似度图)的顶点集，ε为定义的邻域半径；
[0110] 所述测地距离通过下式计算：
[0111]
[0112] i,j表示邻接矩阵(无向相似度图)中的两个顶点，P(*)为权重函数，即特征的欧式距离，P(x)>0，γ(t)表示从顶点i到顶点j的所有路径的集合。
[0113] 3-3，使用归一化割方法对弥散张量图像进行分割，根据分割结果形成个体丘脑的分区。
[0114] 下面通过对比实验，不同方法得到的个体分区组合起来构建丘脑功能群组分区图谱。
[0115] 本实施例的数据集采用了人类连接组计划项目数据，总共100个正常人的弥散张量成像模态数据，参见https://www.humanconnectome.org/，包括T1结构
磁共振成像与弥散张量成像数据。T1成像的参数为：repetition time(TR)为1000ms，echo time(TE)为2.14ms，采集矩阵为320×320，层数为256，层厚为0.7mm，
视野为224×224mm^2。弥散张量成像的采集参数如下：TR为5520ms，TE为89.5ms，视野为210×210mm^2，层数为145，层厚为
1.25mm，b值1000s/m^2、2000s/m^2、3000s/m^2分别获得288个梯度方向的图像，b等于0获得
18个基准图像。
[0116] 本发明
算法是基于标准空间2mm下，其成像大小应为90*109*90。在样本的原始空间下，其丘脑体素数量为6000个左右，而在标准空间下只有约2000个，因此能降低受弥散张量成像噪声的影响。概率性纤维追踪算法得到的结果是基于每个体素的三维图像，其图像里每个体素的值表示了此体素与当前体素的纤维连接强度，因此在原始空间下，每个样本都有约6000个三维的纤维追踪图像。
[0117] 本发明提取的丘脑体素精细皮层分区特征的维度为268，而左右丘脑体素的数目均为1164，因此子空间自编码器模型的输入数据特征构成的矩阵大小为1164*268。对于矩阵的每一行代表当前体素与268个脑区的纤维连接强度，使用min-max归一化方法作为特征预处理手段，将特征矩阵的每一行元素值变换到-1到1的区间内。
[0118] 可复现性广泛应用于验证脑功能分区方法的鲁棒性。鲁棒的算法在充足的样本上会得到高度相似的分区结果。通过基于深度子空间聚类的丘脑功能分区算法得到群组丘脑功能分区图谱之后，依次在个体水平上把个体丘脑功能分区结果与群组丘脑功能分区图谱进行比对，以Dice相似度为指标，个体丘脑功能分区结果与群组丘脑功能分区图谱越相似则表示算法鲁棒性强，可复现性高。可复现性的验证分为两个方案：一是基于全部数据集；二是基于奇偶分半数据集。
[0119] 基于全部数据集的可复现性实验过程为：在标准空间对个体丘脑进行功能分区，然后通过多数表决法构建群组图谱，依次将个体功能分区结果与群组图谱进行比对得到Dice相似度。
[0120] 表1是基于脑连接的丘脑功能分区方法结果对比，表中(a)为本发明方法；(b)为直接使用深度子空间聚类网络自表达层系数构建亲和度矩阵的方法；(c)为使用AAL图谱进行特征提取的方法；(d)为使用粗糙皮层图谱进行特征提取的方法。可以发现在全部样本集、奇数样本集、偶数样本集以及奇数样本集群组图谱与偶数样本集群组图谱对比上本发明方法均得到最高的Dice系数值，表示本发明方法可复现性显著强于对比实验方法。
[0121] 表1
[0122]
[0123] 以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以
权利要求书的保护范围为准。
