首页> 中国专利> 扩散峭度张量成像的高阶张量特征参数提取方法

扩散峭度张量成像的高阶张量特征参数提取方法

摘要

本发明涉及医疗器械技术领域,为提出一种新的实用的人脑白质结构异常的磁共振成像参数提取方法,通过一定的特征提取方法对组织的各向异性进行分析,达到客观评价病变程度的目的,本发明采取的技术方案是,扩散峭度张量成像的高阶张量特征参数提取方法,包括如下步骤:通过受试者在磁共振扫描仪上采集组织沿多个方向的扩散加权信号,拟合得到反映组织内水分子扩散分布概率密度函数特征的二阶及四阶张量,通过一定的张量分析得到由于组织结构异常引起的水分子扩散分布异常的病变特征,主要有扩散系数和峭度系数两个不同角度的参数。本发明主要应用于医疗器械设计制造。

著录项

  • 公开/公告号CN103142229A

    专利类型发明专利

  • 公开/公告日2013-06-12

    原文格式PDF

  • 申请/专利权人 天津大学;

    申请/专利号CN201310056881.4

  • 申请日2013-02-22

  • 分类号A61B5/055(20060101);

  • 代理机构12201 天津市北洋有限责任专利代理事务所;

  • 代理人刘国威

  • 地址 300072 天津市南开区卫津路92号

  • 入库时间 2024-02-19 17:57:55

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2015-12-02

    授权

    授权

  • 2013-07-17

    实质审查的生效 IPC(主分类):A61B5/055 申请日:20130222

    实质审查的生效

  • 2013-06-12

    公开

    公开

说明书

技术领域

本发明涉及医疗器械技术领域,具体讲,涉及扩散峭度张量成像的高阶张量特征参数提 取方法。

背景技术

阿尔兹海默症(Alzheimer’s disease,AD)是老年痴呆症中最普遍最重要的一种,也是严重 影响老年人健康的三大疾病之一。随着目前我国乃至直接全球社会老龄化程度的加剧,此类 疾病给社会带来的经济负担和社会负担也随之大大的加重。此类疾病依然还没有能治愈的药 物或者其它方法,而且该疾病也是有着发病周期长具有隐匿性,早期的诊断对于该疾病的早 处理有着十分重要的意义。目前针对阿尔兹海默症的诊断研究开展较多,也有一定的研究进 展,但是依然缺乏统一标准,国内外学者也积极采用了多种方法进行研究。

发明内容

本发明旨在克服现有技术的不足,提出一种新的实用的人脑白质结构异常的磁共振成像 参数提取方法,通过一定的特征提取方法对组织的各向异性进行分析,达到客观评价病变程 度的目的,为此,本发明采取的技术方案是,扩散峭度张量成像的高阶张量特征参数提取方 法,包括如下步骤:通过受试者在磁共振扫描仪上采集组织沿多个方向的扩散加权信号,拟 合得到反映组织内水分子扩散分布概率密度函数特征的二阶及四阶张量,通过一定的张量分 析得到由于组织结构异常引起的水分子扩散分布异常的病变特征,主要有扩散系数和峭度系 数两个不同角度的参数。

采集组织沿多个方向的扩散加权信号进一步细化为:

1.1数据的采集

在3.0T的磁共振扫面仪器上,采用单次激发平面回波(SE-EPI)序列进行30个敏感梯度方 向3个扩散敏感因子b值(b=0,b1000和b2000s/mm2)的扩散加权DWI信号采集。

1.2数据的预处理

首先采用FSL(the FMRIB Software Library)软件的FDT对数据进行涡流矫正和头动矫 正操作,以及SUSAN对数据进行初步的去噪处理,其次在张量拟合再对数据进行非局部均 值滤波操作。

采用非局部的均值滤波(Non Local Means Filter)方法,对每个像素点xi处理如下式(5-1):

NL(v)(xi)=ΣxiIw(xi,yi)v(xj)(5-1)

这里,I代表待处理图像;v(xj)代表像素xj的灰度值;NL(v)(xi)代表像素v(xj)经过非局部 均值滤波的结果;w(xi,xj)代表恢复v(xi)对应于v(xj)的权值;其中i和j均代表像素的索引编号。

w(xi,yj)=1z(i)ed(v(Ni),v(Nj))()2(5-2)

这里,Z(i)是标准化的常数Z(i)=∑jw(i,j),θ是使用残差方法估计的噪声标准偏差, h是一个滤波参数,距离d表示为如下:d(v(Ni),v(Nj))=1NΣkNΔ(v(yk),v(zk)),这里 N=card(Ni)=card(Nj),yk和zk分别代表xi和yj领域Ni和Nj的第k个点,对于灰度级图像, Δ(v(yk),v(zk))=|(v(yk),v(zk))||。

张量的拟合方案:

S(b)是添加了梯度磁场的信号强度,即水分子在某个方向上的扩散情况,b值是扩散磁敏 感因子;S(0)是未加梯度磁场的信号强度,及水分子的背景扩散情况;g是梯度敏感方向的单 位向量,D是对应组织的扩散张量,是一个3×3的二维矩阵;基于量化水分子真实扩散分布 偏离高斯扩散分布程度的峭度的扩散峭度成像(DKI,diffusion kurtosis imaging)有如下模型, 如式(5-3)

S(b)=S(0)exp(-bDapp+16b2Dapp2Kapp).(5-3)

这里,

Dapp=gTDg=Dg2=Σi=13Σj=13Dijgigj·(5-4-1)

Kapp=MD2Dapp2Wg4=Σi=13Σj=13Σk=13Σl=13Wijklgigjgkgl.(5-4-2)

其中,Dapp和Kapp分别为组织在梯度敏感方向g上的表观扩散系数和表观峭度系数;对应 的张量即为D和W,分别对应一个二维三阶的实对称矩阵D和一个四维三阶的实对称矩阵W; MD是对应组织部位的平均扩散系数,由扩散张量矩阵D计算得来的;g依然是扩散敏感梯度 方向,gi、gj、gk、gl是g的第i、j、k、l个分量元素,同上标T表示转置;

上述式(5-3)通过两个不同b值的信号S(b1)和S(b2),b1=1000s/mm2、b2= 2000s/mm2以及一个b=0的信号;联立得到对应b1和b2的关于未知变量Dapp和Kapp的二元二 次方程组如式(5-5),可解得对应30个方向上的Dapp和Kapp

ln[S(b1)/S(0)]=-b1Dapp+16b12Dapp2Kappln[S(b2)/S(0)]=-b2Dapp+16b22Dapp2Kapp·(5-5)

D的9个元素中仅有6个未知量,而W的81个元素中只有15个位置变量,分别对(5-4-1) 和(5-4-2)展开,可得到如下式,

Dapp=D11g12+D22g22+D33g32+2(D12g1g2+D23g2g3+D13g1g3).(5-6)

Wg4=W1111g14+W2222g24+4W1112g13g2+4W1113g13g3+4W2221g23g1+4W2223g23g3

+4W3331g33g1+4W3332g33g2+12W1123g12g2g3+12W2213g22g1g3+12W3312g32g1g2

+6W1122g12g22+6W1133g12g32+6W2233g22g32.(5-7)

D11表示矩阵D的元素,其下标表示该元素在矩阵中的位置,与之类似W3321。通过(5-4)、 (5-6)、(5-7)式采用非线性最小二乘算法即可得到最优的张量矩阵D和W,即对应的张量D和 W。

特征提取方案:

张量的分解算法

首先是扩散张量矩阵D的分解,依然采用对角化分解方法,如下式,

D=[v1v2v3]Tλ1000λ2000λ3[v1v2v3]·

这里,λi和vi是对应的特征值和特征向量(λ1≥λ2≥λ3),且3个vi两两正交,分别对应三 个特征值的方向;其中作为特征值的λi代表组织内三分方向上的扩散系数大小,且最大特征 值对应方向极为组织内神经纤维的主要走向,另两个特征值反映神经纤维髓鞘的特性;

其次对于高阶的峭度张量W,转化为解决此优化问题,

max,Wx4s.t.xTx=1

即得到如下方程组,

Wx3=λxxTx=1

带入x=(x1,x2,x3)展开如下,

W1111x13+W1222x23+W1333x33+3W1112x12x2+3W1113x12x3+3W1223x22x3+3W1122x1x22+3W1133x1x32+3W1233x2x32=λx1W2111x13+W2222x23+W2333x33+3W2112x12x2+3W2113x12x3+3W2223x22x3+3W2122x1x22+3W2133x1x32+3W2233x2x32=λx2W3111x13+W3222x23+W3333x33+3W3112x12x2+3W3113x12x3+3W3223x22x3+3W3122x1x22+3W3133x1x32+3W3233x2x32=λx3x12+x22+x32=1

其解的情况如下,

如果W2111=W3111=0,即x2=x3=0,则λ=W1111,且对应特征向量为x=(±1,0,0);若 x2≠0,x3=0,则另t=x1/x2,消元得到如下方程,

-W2111t4+(W1111-3W2112)t3+3(W1112-W2122)t2+(3W1122-W2222)t+W1222=0W3111t3+3W3112t2+3W3122t+W3222=0

取其实根t,则对应的特征向量和特征值为,

x=±1t2+1(t,1,0)Tλ=Wx4

若x3≠0,则另u=x1/x3,v=x2/x3,消元得到如下方程,

-W3111u4-3W3112u3v+(W1111-3W3113)u3-3W3122u2v2+3(W1112-6W3123)u2v+(3W1113-3W3133)u2-3W3223uv2-W3222uv3+3W1122uv2+(6W1123-3W3233)uv++(3W1133-W3333)u+W1222v3+3W1223v2+3W1233v+W1333=0-W3111u3v+W2111u3-3W3112u2v2+(3W2112-3W3113)u2v+3W2113u2-3W3122uv3+(3W2122-6W3123)uv2+(6W2123-3W3133)uv+3W2133u+3W2223v2-W3222v4+(W2222-3W3223)v3-3W3233v2+(3W2233-W3333)v+W2333=0

取其实根对(u,v),则特征向量和特征值为,

x=±(u,v,1)Tu2+v2+1λ=Wx4

由此得到W的特征向量xi和特征值λi,1≤i≤n,且3≤n≤13代表特征值的个数。

特征参数计算

根据5.4.1的峭度张量特征分解结果给出峭度张量的提取特征,

平均峭度系数MK,mean kurtosis coefficient,

MK=mean(λi)

峭度各向异性KA,kurtosis anisotropy,

KA=nn-1.Σi=1n(λi-MK)2Σi=1nλi2

对每个i,将xi按着三个正交的vi′方向进行分组,即寻找其对应的最优i‘,

得到三组X,分别对每组Xi′向量对应的λXi′求平均,则

Ki′=mean(λXi′)(i′=1,2,3)

径向的峭度系数RK,radial kurtosis coefficient,

RK=K1

轴向的峭度系数AK,axial kurtosis coefficient,

AK=(K2+K3)/2。

本发明的技术特点及效果:

本发明提出了基于峭度信息的表征脑白质异变特征的参数提取方法。在采集MR信号并 经过与处理之后,采用本发明所述的方法进行计算,可以得出包括扩散张量参数在内的以及 峭度张量峭度特征参数的多种参数图像。利用这些参数,通过ROI统计分析、分类算法等可 以获取反映白质退化程度的明显表征,对阿尔兹海默症早期的病变特征具有明显的更敏感的 反映。如附图,及为各种参数反映阿尔兹海默症病人病变部位的图像。

附图说明

图1图像降噪示意图。

图2方法步骤流程图。

图3提取参数图像。

具体实施方式

本发明旨在提出一种新的实用的人脑白质结构异常的磁共振成像参数提取方法,通过多 方向检测的DWI信号,拟合得到扩散位移分布的方差张量和峭度张量,通过一定的特征提取 方法对组织的各向异性进行分析,达到客观评价病变程度的目的。

本发明提出了基于磁共振扩散峭度张量成像(MR-DKI,MR-diffusion kurtosis tensor  imaging)理论的反映组织结构特征的参数图像提取方法。其技术流程是:通过受试者在磁共振 扫描仪上采集组织沿多个方向的扩散加权信号,拟合得到反映组织内水分子扩散分布概率密 度函数特征的二阶及四阶张量,通过一定的张量分析得到由于组织结构异常引起的水分子扩 散分布异常的病变特征,主要有扩散系数和峭度系数两个不同角度的参数。通过上述特征参 数的分析,可以用于临床AD病人的诊断依据。

磁共振扫描仪可以对生物组织呈现多种模态的检测和信息获取,水分子的扩散特性是生 物组织内的十分普遍的现象也能充分反映组织的围观结构特性。磁共振扩散加权成像 (MR-DWI,diffusion weighted imaging)是检测活体组织内水分子扩散情况的技术,DKI技术 就是通过多个方向上的DWI信号采集,获取组织内水分子的扩散分布情况,进而研究组织的 微观结构特性。同时,本专利利用了DKI理论模型中的高阶峭度张量(KT,kurtosis tensor), 通过敏感的四阶张量,能更加敏感和精细的得到组织微观结构的异常,为阿尔兹海默症的进 一步诊断提供有效的依据。

本发明方法步骤流程图如附图2。

1.1数据的采集方案

在3.0T的磁共振扫面仪器上,采用单次激发平面回波(SE-EPI)序列进行30个敏感梯度方 向3个扩散敏感因子b值(b=0,b=1000和b=2000s/mm2)的扩散加权DWI信号采集。另外其 它参数设置为:重复时间TR(repetition time)=10500ms,回波时间TE(ehco time)=103ms, 重复次数为Average(or NEXT)=1,数据获取时间TA(Acquisition Time)=11’14’’,噪声水平 为30,获取图像矩阵大小为128×128,视野FOV=230×230mm2,层厚为1.8mm,全脑采集 73层,层间无间隔。因为对每个被试,将采集得到一个大小为128×128×73×61的矩阵数 据集,采集时间为约11分钟。

1.2数据的预处理方案

首先采用英国牛津大学的FSL(the FMRIB Software Library)软件的FDT对数据进行涡 流矫正和头动矫正操作,以及SUSAN对数据进行初步的去噪处理。其次在张量拟合再对数 据进行非局部均值滤波操作,本发明采用非局部的均值滤波方法,对每个像素点xi处理如下 式(5-1):(见图1)

NL(v)(xi)=ΣxiIw(xi,yi)v(xj)(5-1)

这里,I代表待处理图像;xi,xj均代表I中的像素点,xi为待处理像素点,0<i≤ card(I),0<j≤card(I),且i≠j,i、j∈Z(Z代表整数,card(N)代表域N的像素点个数); v(xj)代表像素xj的灰度值;w(xi,xj)代表恢复v(xi)对应于v(xj)的权值。

w(xi,yj)=1z(i)ed(v(Ni),v(Nj))()2(5-2)

这里,Z(i)是标准化的常数Z(i)=∑jw(i,j),是使用为残差方法估计的噪声标准偏差, h是一个滤波参数。距离d表示为如下:d(v(Ni),v(Nj))=1NΣkNΔ(v(yk),v(zk)),这里Ni和 Nj分别代表像素点xi和xj对应的领域;v(Ni),v(Nj)分别代表领域Ni和Nj的所有像素点的灰度值 矩阵,N=card(Ni)=card(Nj);yk和zk分别代表xi和xj领域Ni和Nj的第k个点,则 v(yk),v(zk)分别代表其对应灰度值;对于灰度级图像,Δ(v(yk),v(zk))=|(v(yk),v(zk))|。

1.3张量的拟合方案

扩散峭度成像(DKI,diffusion kurtosis imaging)技术是基于传统的扩散张量成像(D TI, diffusion tensor imaging)推广发展而来的,经典的扩散张量成像模型如式S(b)=S(0)e-b(xTDx), 其中S(b)是添加了梯度磁场的信号强度(b=1000,2000s/mm2),即水分子在某个方向上的扩散 情况,b值是扩散磁敏感因子;S(0)是未加梯度磁场的信号强度,及水分子的背景扩散情况; g是梯度敏感方向的单位向量,D是对应组织的扩散张量,是一个3×3的二维矩阵。但是一 个张量D获取的组织结构信息对于微观病变特性很难识别,所以新提出的基于量化水分子真 实扩散分布偏离高斯扩散分布程度的峭度的扩散峭度成像(DKI)技术有如下模型,如式(5-3)

S(b)=S(0)exp(-bDapp+16b2Dapp2Kapp).(5-3)

这里,

Dapp=gTDg=Dg2=Σi=13Σj=13Dijgigj·(5-4-1)

Kapp=MD2Dapp2Wg4=Σi=13Σj=13Σk=13Σl=13Wijklgigjgkgl.(5-4-2)

其中,Dapp和Kapp分别为组织在某方向(这里指数据采集的梯度敏感方向g)上的表观扩散 系数和表观峭度系数;Dij代表张量矩阵D的第i行第j列的元素;Wijkl代表张量矩阵W的第 一维第i,第二维第j,第三维k和第四维第l的元素,对应的张量即为D和W,分别对应一个 二维三阶的实对称矩阵D和一个四维三阶的实对称矩阵W;MD是对应组织部位的平均扩散系 数,由扩散张量矩阵D计算得来的;g依然是扩散敏感梯度方向,gi、gj、gk、glg的第 i、j、k、l个分量元素,同上标T表示转置;

上述式(5-3)通过两个不同b值的信号S(b1)和S(b2)(b1=1000s/mm2、b2= 2000s/mm2)以及一个b0=0的信号,这里b0,b1,b2为b值的三个不同的取值;联立得到对 应b1和b2的关于未知变量Dapp和Kapp的二元二次方程组如式(5-5),可解得对应30个方向上 的Dapp和Kapp

ln[S(b1)/S(0)]=-b1Dapp+16b12Dapp2Kappln[S(b2)/S(0)]=-b2Dapp+16b22Dapp2Kapp·(5-5)

因为张量D和W,对应的矩阵是实对称的,D的9个元素中仅有6个未知量,而W的81 个元素中只有15个位置变量。因此后分别对(5-4-1)和(5-4-2)展开,可得到如下式,

Dapp=D11g12+D22g22+D33g32+2(D12g1g2+D23g2g3+D13g1g3).(5-6)

Wg4=W1111g14+W2222g24+4W1112g13g2+4W1113g13g3+4W2221g23g1+4W2223g23g3

+4W3331g33g1+4W3332g33g2+12W1123g12g2g3+12W2213g22g1g3+12W3312g32g1g2

+6W1122g12g22+6W1133g12g32+6W2233g22g32.(5-7)

通过(5-4)、(5-6)、(5-7)式采用非线性最小二乘算法即可得到最优的张量矩阵D和W, 即对应的张量D和W。

1.4特征提取方案

1.4.1张量的分解算法

首先是扩散张量矩阵D的分解,依然采用对角化分解方法,如下式,

D=[v1v2v3]Tλ1000λ2000λ3[v1v2v3]·

这里,λi和vi是对应D的特征值和特征向量(λ1≥λ2≥λ3),vi是特征分解λi所对应的1×3 的特征向量,这里i=1,2,3,且3个vi两两正交,分别对应三个特征值的方向;其中 其中作为特征值的λi代表组织内三分方向上的扩散系数大小,且最大特征值对应方向极为组 织内神经纤维的主要走向,另两个特征值反映神经纤维髓鞘的特性。

其次对于高阶的峭度张量W,我们转化为解决此优化问题,

max,Wx4s.t.xTx=1

即得到如下方程组,

Wx3=λxxTx=1

带入x=(x1,x2,x3)展开如下,

W1111x13+W1222x23+W1333x33+3W1112x12x2+3W1113x12x3+3W1223x22x3+3W1122x1x22+3W1133x1x32+3W1233x2x32=λx1W2111x13+W2222x23+W2333x33+3W2112x12x2+3W2113x12x3+3W2223x22x3+3W2122x1x22+3W2133x1x32+3W2233x2x32=λx2W3111x13+W3222x23+W3333x33+3W3112x12x2+3W3113x12x3+3W3223x22x3+3W3122x1x22+3W3133x1x32+3W3233x2x32=λx3x12+x22+x32=1

这里诸如W1222均为张量矩阵的某一个独立元素值,即Wijkl各小标在(1,2,3)取值得到; 上式其解的情况如下,

如果W2111=W3111=0,即x2=x3=0,则λ=W1111,且对应特征向量为x=(±1,0,0)。

若x2≠0,x3=0,则另t=x1/x2,消元得到如下方程,

-W2111t4+(W1111-3W2112)t3+3(W1112-W2122)t2+(3W1122-W2222)t+W1222=0W3111t3+3W3112t2+3W3122t+W3222=0

取其实根t,则对应的特征向量和特征值为,

x=±1t2+1(t,1,0)Tλ=Wx4

λ=Wx4

若x3≠0,则另u=x1/x3,v=x2/x3,消元得到如下方程,

-W3111u4-3W3112u3v+(W1111-3W3113)u3-3W3122u2v2+3(W1112-6W3123)u2v+(3W1113-3W3133)u2-3W3223uv2-W3222uv3+3W1122uv2+(6W1123-3W3233)uv++(3W1133-W3333)u+W1222v3+3W1223v2+3W1233v+W1333=0-W3111u3v+W2111u3-3W3112u2v2+(3W2112-3W3113)u2v+3W2113u2-3W3122uv3+(3W2122-6W3123)uv2+(6W2123-3W3133)uv+3W2133u+3W2223v2-W3222v4+(W2222-3W3223)v3-3W3233v2+(3W2233-W3333)v+W2333=0

取其实根对(u,v),则特征向量和特征值为,

x=±(u,v,1)Tu2+v2+1λ=Wx4

由此得到W的特征向量xi和特征值λi,1≤i≤n,且3≤n≤13代表特征值的个数。

1.4.2特征参数计算

经典的扩散张量成像(DTI)算法的常用参数有各向异性FA(fractional anisotropy),平 均扩散系数MD(mean diffusivity),轴向扩散系数AD(axial diffusivity)和径向扩散系数RD (radial diffusivity)。本发明主要根据5.4.1的峭度张量特征分解结果给出峭度张量的提取特 征,

平均峭度系数MK,mean kurtosis coefficient,

MK=mean(λi)

峭度各向异性KA,kurtosis anisotropy,

KA=nn-1.Σi=1n(λi-MK)2Σi=1nλi2

对每个i,将xi对应三个正交的vi′方向进行分组,即求夹角寻找其对应的最优i‘,这里 i‘=1,2,3,代表扩散张量D的三个特征值或特征向量,参见1.4.1,

得到特征向量的三个分组X,分别对每组Xi’向量对应的λXi′求平均,则

Ki′=mean(λXi′)(i’=1,2,3)

径向的峭度系数RK,radial kurtosis coeffocient,

RK=K1

轴向的峭度系数AK,axial kurtosis coefficient,

AK=(K2+K3)/2

1.5特征的验证方案

分别对更具现有阿尔兹海默症最新诊断标准,获取阿尔兹海默症(AD)、轻度认知障碍 (MCI)和控制组(ND)三组病人MR数据,采用本发明所述的方法进行特征提取,得到扩 散信息参数和峭度信息参数若干。然后对阿尔兹海默症发病的几个显著部位(如颞叶、额叶 等)的ROI区域或者体积进行统计病理分析等。

本发明针对阿尔兹海默症发病周期长,具有隐匿性,且诊断标准主观化的特点,提出了 利用磁共振扩散结构成像技术进行客观诊断的新方法。基于峭度张量得到峭度信息相关参数, 一次作为脑白质病变程度的输入参数,实现对患阿尔兹海默症高风险的老年人进行客观有效 的诊断识别,具有重大的社会意义。

去获取专利,查看全文>

相似文献

  • 专利
  • 中文文献
  • 外文文献
获取专利

客服邮箱:kefu@zhangqiaokeyan.com

京公网安备:11010802029741号 ICP备案号:京ICP备15016152号-6 六维联合信息科技 (北京) 有限公司©版权所有
  • 客服微信

  • 服务号