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

Abstract

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

Description

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

技术领域

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

背景技术

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

发明内容

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

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

1.1数据的采集

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

1.2数据的预处理

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

采用非局部的均值滤波（NonLocalMeansFilter）方法，对每个像素点x_i处理如下式（5-1）：

$\mathrm{NL}\left(v\right)\left(x_i\right)={\mathrm{\Σ}}_{x_i\∈I}w(x_i,y_i)v\left(x_j\right)(5-1)$

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

$w(x_i,y_j)=\frac{1}{z\left(i\right)}{e}^{\frac{d(v\left(N_i\right),v\left(N_j\right))}{(\mathrm{h\θ}{)}^2}}(5-2)$

这里，Z(i)是标准化的常数Z(i)＝∑_jw(i，j)，θ是使用残差方法估计的噪声标准偏差，h是一个滤波参数，距离d表示为如下： $d(v\left(N_i\right),v\left(N_j\right))=\sqrt{\frac{1}{N}{\mathrm{\Σ}}_k^N\mathrm{\Δ}(v\left(y_k\right),v\left(z_k\right))},$ 这里N＝card(Ni)＝card(N_j)，y_k和z_k分别代表x_i和y_j领域N_i和N_j的第k个点，对于灰度级图像，Δ(v(y_k)，v(z_k))＝|(v(y_k)，v(z_k))||。

张量的拟合方案：

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

$S\left(b\right)=S\left(0\right)\exp (-{\mathrm{bD}}_{\mathrm{app}}+\frac{1}{6}{b}^2{D_{\mathrm{app}}}^2{K}_{\mathrm{app}}).(5-3)$

这里，

$D_{\mathrm{app}}=g^T\mathrm{Dg}={\mathrm{Dg}}^2={\mathrm{\Σ}}_{i=1}^3{\mathrm{\Σ}}_{j=1}^3{D}_{\mathrm{ij}}{g}_i{g}_j\·(5-4-1)$

$K_{\mathrm{app}}=\frac{{\mathrm{MD}}^2}{{D_{\mathrm{app}}}^2}W{g}^4={\mathrm{\Σ}}_{i=1}^3{\mathrm{\Σ}}_{j=1}^3{\mathrm{\Σ}}_{k=1}^3{\mathrm{\Σ}}_{l=1}^3{W}_{\mathrm{ijkl}}{g}_i{g}_j{g}_k{g}_l.(5-4-2)$

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

上述式（5-3）通过两个不同b值的信号S(b1)和S(b2)，b1＝1000s/mm²、b2＝2000s/mm²以及一个b＝0的信号；联立得到对应b1和b2的关于未知变量D_app和K_app的二元二次方程组如式（5-5），可解得对应30个方向上的D_app和K_app；

$\left\{\begin{array}{c}\ln [S\left(b_1\right)/S\left(0\right)]=-b_1{D}_{\mathrm{app}}+\frac{1}{6}{b}_1^2{D_{\mathrm{app}}}^2{K}_{\mathrm{app}}\\ \ln [S\left(b_2\right)/S\left(0\right)]=-b_2{D}_{\mathrm{app}}+\frac{1}{6}{b}_2^2{D_{\mathrm{app}}}^2{K}_{\mathrm{app}}\end{array}\right.\·(5-5)$

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

$D_{\mathrm{app}}=D_{11}{g}_1^2+D_{22}{g}_2^2+D_{33}{g}_3^2+2(D_{12}{g}_1{g}_2+D_{23}{g}_2{g}_3+D_{13}{g}_1{g}_3).(5-6)$

${\mathrm{Wg}}^4=W_{1111}{g}_1^4+W_{2222}{g}_2^4+4W_{1112}{g}_1^3{g}_2+4W_{1113}{g}_1^3{g}_3+4W_{2221}{g}_2^3{g}_1+{4W}_{2223}{g}_2^3{g}_3$

$+{4W}_{3331}{g}_3^3{g}_1+{4W}_{3332}{g}_3^3{g}_2+12W_{1123}{g}_1^2{g}_2{g}_3+12W_{2213}{g}_2^2{g}_1{g}_3+12W_{3312}{g}_3^2{g}_1{g}_2$

$+6W_{1122}{g}_1^2{g}_2^2+{6W}_{1133}{g}_1^2{g}_3^2+{6W}_{2233}{g}_2^2{g}_3^2.(5-7)$

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

特征提取方案：

张量的分解算法

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

$D={\left[v_1{v}_2{v}_3\right]}^T\left[\begin{array}{ccc}{\mathrm{\λ}}_1& 0& 0\\ 0& {\mathrm{\λ}}_2& 0\\ 0& 0& {\mathrm{\λ}}_3\end{array}\right]\left[v_1{v}_2{v}_3\right]\·$

这里，λ_i和v_i是对应的特征值和特征向量（λ₁≥λ₂≥λ₃），且3个v_i两两正交，分别对应三个特征值的方向；其中作为特征值的λ_i代表组织内三分方向上的扩散系数大小，且最大特征值对应方向极为组织内神经纤维的主要走向，另两个特征值反映神经纤维髓鞘的特性；

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

$\left\{\begin{array}{c}\max ,W{x}^4\\ s.t.x^{T}x=1\end{array}\right.$

即得到如下方程组，

$\left\{\begin{array}{c}W{x}^3=\mathrm{\λx}\\ x^{T}x=1\end{array}\right.$

带入x＝(x₁，x₂，x₃)展开如下，

$\left\{\begin{array}{c}{W}_{1111}{x}_1^3+W_{1222}{x}_2^3+W_{1333}{x}_3^3+{3W}_{1112}{x}_1^2{x}_2+{3W}_{1113}{x}_1^2{x}_3+{3W}_{1223}{x}_2^2{x}_3+{3W}_{1122}{x}_1{x}_2^2\\ +{3W}_{1133}{x}_1{x}_3^2+{3W}_{1233}{x}_2{x}_3^2={\mathrm{\λ}}_{x_1}\\ W_{2111}{x}_1^3+W_{2222}{x}_2^3+W_{2333}{x}_3^3+{3W}_{2112}{x}_1^2{x}_2+{3W}_{2113}{x}_1^2{x}_3+{3W}_{2223}{x}_2^2{x}_3+{3W}_{2122}{x}_1{x}_2^2\\ +3W_{2133}{x}_1{x}_3^2+{3W}_{2233}{x}_2{x}_3^2={\mathrm{\λ}}_{x_2}\\ W_{3111}{x}_1^3+W_{3222}{x}_2^3+W_{3333}{x}_3^3+{3W}_{3112}{x}_1^2{x}_2+{3W}_{3113}{x}_1^2{x}_3+{3W}_{3223}{x}_2^2{x}_3+{3W}_{3122}{x}_1{x}_2^2\\ +{3W}_{3133}{x}_1{x}_3^2+{3W}_{3233}{x}_2{x}_3^2={\mathrm{\λ}}_{x_3}\\ x_1^2+x_2^2+x_3^2=1\end{array}\right.$

其解的情况如下，

如果W₂₁₁₁＝W₃₁₁₁＝0，即x₂＝x₃＝0，则λ＝W₁₁₁₁，且对应特征向量为x＝(±1，0，0)；若x₂≠0，x₃＝0，则另t＝x₁/x₂，消元得到如下方程，

$\left\{\begin{array}{c}-W_{2111}{t}^4+(W_{1111}-{3W}_{2112})t^3+3(W_{1112}-W_{2122})t^2+({3W}_{1122}-W_{2222})t+W_{1222}=0\\ W_{3111}{t}^3+{3W}_{3112}{t}^2+{3W}_{3122}t+W_{3222}=0\end{array}\right.$

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

$\begin{array}{c}x=\±\frac{1}{\sqrt{t^2+1}}{(t,\mathrm{1,0})}^T\\ \mathrm{\λ}=W{x}^4\end{array}$

若x₃≠0，则另u＝x₁/x₃，v＝x₂/x₃，消元得到如下方程，

$\left\{\begin{array}{c}-W_{3111}{u}^4-{3W}_{3112}{u}^{3}v+(W_{1111}-{3W}_{3113})u^3-{3W}_{3122}{u}^2{v}^2+3(W_{1112}-{6W}_{3123})u^{2}v\\ +({3W}_{1113}-{3W}_{3133})u^2-{3W}_{3223}u{v}^2-W_{3222}{\mathrm{uv}}^3+{3W}_{1122}{\mathrm{uv}}^2+({6W}_{1123}-{3W}_{3233})\mathrm{uv}+\\ +({3W}_{1133}-W_{3333})u+W_{1222}{v}^3+{3W}_{1223}{v}^2+{3W}_{1233}v+W_{1333}=0\\ -W_{3111}{u}^{3}v+W_{2111}{u}^3-{3W}_{3112}{u}^2{v}^2+({3W}_{2112}-{3W}_{3113})u^{2}v+{3W}_{2113}{u}^2-{3W}_{3122}{\mathrm{uv}}^3\\ +({3W}_{2122}-{6W}_{3123}){\mathrm{uv}}^2+({6W}_{2123}-{3W}_{3133})\mathrm{uv}+{3W}_{2133}u+{3W}_{2223}{v}^2-W_{3222}{v}^4\\ +\left(W_{2222}-{3W}_{3223}\right)v^3-{3W}_{3233}{v}^2+({3W}_{2233}-W_{3333})v+W_{2333}=0\end{array}\right.$

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

$\begin{array}{c}x=\frac{\±{(u,v,1)}^T}{\sqrt{u^2+v^2+1}}\\ \mathrm{\λ}=W{x}^4\end{array}$

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

特征参数计算

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

平均峭度系数MK，meankurtosiscoefficient，

MK＝mean(λ_i)

峭度各向异性KA，kurtosisanisotropy，

$\mathrm{KA}=\sqrt{\frac{n}{n-1}}.\sqrt{\frac{{\mathrm{\Σ}}_{i=1}^n{({\mathrm{\λ}}_i-\mathrm{MK})}^2}{{\mathrm{\Σ}}_{i=1}^n{\mathrm{\λ}}_i^2}}$

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

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

K_i′＝mean(λ_Xi′)(i′＝1，2，3)

径向的峭度系数RK，radialkurtosiscoefficient，

RK＝K₁

轴向的峭度系数AK，axialkurtosiscoefficient，

AK＝(K₂+K₃)/2。

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

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

附图说明

图1图像降噪示意图。

图2方法步骤流程图。

图3提取参数图像。

具体实施方式

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

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

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

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

1.1数据的采集方案

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

1.2数据的预处理方案

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

$\mathrm{NL}\left(v\right)\left(x_i\right)={\mathrm{\Σ}}_{x_i\∈I}w(x_i,y_i)v\left(x_j\right)(5-1)$

这里，I代表待处理图像；x_i，x_j均代表I中的像素点，x_i为待处理像素点，0＜i≤card(I)，0＜j≤card(I)，且i≠j，i、j∈Z（Z代表整数，card(N)代表域N的像素点个数）；v(x_j)代表像素x_j的灰度值；w(x_i，x_j)代表恢复v(x_i)对应于v(x_j)的权值。

$w(x_i,y_j)=\frac{1}{z\left(i\right)}{e}^{\frac{d(v\left(N_i\right),v\left(N_j\right))}{(\mathrm{h\θ}{)}^2}}(5-2)$

这里，Z(i)是标准化的常数Z(i)＝∑_jw(i，j)，是使用为残差方法估计的噪声标准偏差，h是一个滤波参数。距离d表示为如下： $d(v\left(N_i\right),v\left(N_j\right))=\sqrt{\frac{1}{N}{\mathrm{\Σ}}_k^N\mathrm{\Δ}(v\left(y_k\right),v\left(z_k\right))},$ 这里N_i和N_j分别代表像素点x_i和x_j对应的领域；v(N_i)，v(N_j)分别代表领域N_i和N_j的所有像素点的灰度值矩阵，N＝card(N_i)＝card(N_j)；y_k和z_k分别代表x_i和x_j领域N_i和N_j的第k个点，则v(y_k)，v(z_k)分别代表其对应灰度值；对于灰度级图像，Δ(v(y_k)，v(z_k))＝|(v(y_k)，v(z_k))|。

1.3张量的拟合方案

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

$S\left(b\right)=S\left(0\right)\exp (-{\mathrm{bD}}_{\mathrm{app}}+\frac{1}{6}{b}^2{D_{\mathrm{app}}}^2{K}_{\mathrm{app}}).(5-3)$

这里，

$D_{\mathrm{app}}=g^T\mathrm{Dg}={\mathrm{Dg}}^2={\mathrm{\Σ}}_{i=1}^3{\mathrm{\Σ}}_{j=1}^3{D}_{\mathrm{ij}}{g}_i{g}_j\·(5-4-1)$

$K_{\mathrm{app}}=\frac{{\mathrm{MD}}^2}{{D_{\mathrm{app}}}^2}W{g}^4={\mathrm{\Σ}}_{i=1}^3{\mathrm{\Σ}}_{j=1}^3{\mathrm{\Σ}}_{k=1}^3{\mathrm{\Σ}}_{l=1}^3{W}_{\mathrm{ijkl}}{g}_i{g}_j{g}_k{g}_l.(5-4-2)$

其中，D_app和K_app分别为组织在某方向(这里指数据采集的梯度敏感方向g)上的表观扩散系数和表观峭度系数；D_ij代表张量矩阵D的第i行第j列的元素；W_ijkl代表张量矩阵W的第一维第i，第二维第j，第三维k和第四维第l的元素，对应的张量即为D和W，分别对应一个二维三阶的实对称矩阵D和一个四维三阶的实对称矩阵W；MD是对应组织部位的平均扩散系数，由扩散张量矩阵D计算得来的；g依然是扩散敏感梯度方向，g_i、g_j、g_k、g_l是_g的第i、j、k、l个分量元素，同上标T表示转置；

上述式（5-3）通过两个不同b值的信号S(b₁)和S(b₂)(b₁＝1000s/mm²、b₂＝2000s/mm²)以及一个b₀＝0的信号，这里b₀，b₁，b₂为b值的三个不同的取值；联立得到对应b₁和b₂的关于未知变量D_app和K_app的二元二次方程组如式（5-5），可解得对应30个方向上的D_app和K_app。

$\left\{\begin{array}{c}\ln [S\left(b_1\right)/S\left(0\right)]=-b_1{D}_{\mathrm{app}}+\frac{1}{6}{b}_1^2{D_{\mathrm{app}}}^2{K}_{\mathrm{app}}\\ \ln [S\left(b_2\right)/S\left(0\right)]=-b_2{D}_{\mathrm{app}}+\frac{1}{6}{b}_2^2{D_{\mathrm{app}}}^2{K}_{\mathrm{app}}\end{array}\right.\·(5-5)$

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

$D_{\mathrm{app}}=D_{11}{g}_1^2+D_{22}{g}_2^2+D_{33}{g}_3^2+2(D_{12}{g}_1{g}_2+D_{23}{g}_2{g}_3+D_{13}{g}_1{g}_3).(5-6)$

${\mathrm{Wg}}^4=W_{1111}{g}_1^4+W_{2222}{g}_2^4+4W_{1112}{g}_1^3{g}_2+4W_{1113}{g}_1^3{g}_3+4W_{2221}{g}_2^3{g}_1+{4W}_{2223}{g}_2^3{g}_3$

$+{4W}_{3331}{g}_3^3{g}_1+{4W}_{3332}{g}_3^3{g}_2+12W_{1123}{g}_1^2{g}_2{g}_3+12W_{2213}{g}_2^2{g}_1{g}_3+12W_{3312}{g}_3^2{g}_1{g}_2$

$+6W_{1122}{g}_1^2{g}_2^2+{6W}_{1133}{g}_1^2{g}_3^2+{6W}_{2233}{g}_2^2{g}_3^2.(5-7)$

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

1.4特征提取方案

1.4.1张量的分解算法

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

$D={\left[v_1{v}_2{v}_3\right]}^T\left[\begin{array}{ccc}{\mathrm{\λ}}_1& 0& 0\\ 0& {\mathrm{\λ}}_2& 0\\ 0& 0& {\mathrm{\λ}}_3\end{array}\right]\left[v_1{v}_2{v}_3\right]\·$

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

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

$\left\{\begin{array}{c}\max ,W{x}^4\\ s.t.x^{T}x=1\end{array}\right.$

即得到如下方程组，

$\left\{\begin{array}{c}W{x}^3=\mathrm{\λx}\\ x^{T}x=1\end{array}\right.$

带入x＝(x₁，x₂，x₃)展开如下，

$\left\{\begin{array}{c}{W}_{1111}{x}_1^3+W_{1222}{x}_2^3+W_{1333}{x}_3^3+{3W}_{1112}{x}_1^2{x}_2+{3W}_{1113}{x}_1^2{x}_3+{3W}_{1223}{x}_2^2{x}_3+{3W}_{1122}{x}_1{x}_2^2\\ +{3W}_{1133}{x}_1{x}_3^2+{3W}_{1233}{x}_2{x}_3^2={\mathrm{\λ}}_{x_1}\\ W_{2111}{x}_1^3+W_{2222}{x}_2^3+W_{2333}{x}_3^3+{3W}_{2112}{x}_1^2{x}_2+{3W}_{2113}{x}_1^2{x}_3+{3W}_{2223}{x}_2^2{x}_3+{3W}_{2122}{x}_1{x}_2^2\\ +3W_{2133}{x}_1{x}_3^2+{3W}_{2233}{x}_2{x}_3^2={\mathrm{\λ}}_{x_2}\\ W_{3111}{x}_1^3+W_{3222}{x}_2^3+W_{3333}{x}_3^3+{3W}_{3112}{x}_1^2{x}_2+{3W}_{3113}{x}_1^2{x}_3+{3W}_{3223}{x}_2^2{x}_3+{3W}_{3122}{x}_1{x}_2^2\\ +{3W}_{3133}{x}_1{x}_3^2+{3W}_{3233}{x}_2{x}_3^2={\mathrm{\λ}}_{x_3}\\ x_1^2+x_2^2+x_3^2=1\end{array}\right.$

这里诸如W₁₂₂₂均为张量矩阵的某一个独立元素值，即W_ijkl各小标在（1,2,3）取值得到；上式其解的情况如下，

如果W₂₁₁₁＝W₃₁₁₁＝0，即x₂＝x₃＝0，则λ＝W₁₁₁₁，且对应特征向量为x＝(±1，0，0)。

若x₂≠0，x₃＝0，则另t＝x₁/x₂，消元得到如下方程，

$\left\{\begin{array}{c}-W_{2111}{t}^4+(W_{1111}-{3W}_{2112})t^3+3(W_{1112}-W_{2122})t^2+({3W}_{1122}-W_{2222})t+W_{1222}=0\\ W_{3111}{t}^3+{3W}_{3112}{t}^2+{3W}_{3122}t+W_{3222}=0\end{array}\right.$

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

$\begin{array}{c}x=\±\frac{1}{\sqrt{t^2+1}}{(t,\mathrm{1,0})}^T\\ \mathrm{\λ}=W{x}^4\end{array}$

λ＝Wx⁴

若x₃≠0，则另u＝x₁/x₃，v＝x₂/x₃，消元得到如下方程，

$\left\{\begin{array}{c}-W_{3111}{u}^4-{3W}_{3112}{u}^{3}v+(W_{1111}-{3W}_{3113})u^3-{3W}_{3122}{u}^2{v}^2+3(W_{1112}-{6W}_{3123})u^{2}v\\ +({3W}_{1113}-{3W}_{3133})u^2-{3W}_{3223}u{v}^2-W_{3222}{\mathrm{uv}}^3+{3W}_{1122}{\mathrm{uv}}^2+({6W}_{1123}-{3W}_{3233})\mathrm{uv}+\\ +({3W}_{1133}-W_{3333})u+W_{1222}{v}^3+{3W}_{1223}{v}^2+{3W}_{1233}v+W_{1333}=0\\ -W_{3111}{u}^{3}v+W_{2111}{u}^3-{3W}_{3112}{u}^2{v}^2+({3W}_{2112}-{3W}_{3113})u^{2}v+{3W}_{2113}{u}^2-{3W}_{3122}{\mathrm{uv}}^3\\ +({3W}_{2122}-{6W}_{3123}){\mathrm{uv}}^2+({6W}_{2123}-{3W}_{3133})\mathrm{uv}+{3W}_{2133}u+{3W}_{2223}{v}^2-W_{3222}{v}^4\\ +\left(W_{2222}-{3W}_{3223}\right)v^3-{3W}_{3233}{v}^2+({3W}_{2233}-W_{3333})v+W_{2333}=0\end{array}\right.$

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

$\begin{array}{c}x=\frac{\±{(u,v,1)}^T}{\sqrt{u^2+v^2+1}}\\ \mathrm{\λ}=W{x}^4\end{array}$

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

1.4.2特征参数计算

经典的扩散张量成像（DTI）算法的常用参数有各向异性FA（fractionalanisotropy），平均扩散系数MD（meandiffusivity），轴向扩散系数AD（axialdiffusivity）和径向扩散系数RD（radialdiffusivity）。本发明主要根据5.4.1的峭度张量特征分解结果给出峭度张量的提取特征，

平均峭度系数MK，meankurtosiscoefficient，

MK＝mean(λ_i)

峭度各向异性KA，kurtosisanisotropy，

$\mathrm{KA}=\sqrt{\frac{n}{n-1}}.\sqrt{\frac{{\mathrm{\Σ}}_{i=1}^n{({\mathrm{\λ}}_i-\mathrm{MK})}^2}{{\mathrm{\Σ}}_{i=1}^n{\mathrm{\λ}}_i^2}}$

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

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

K_i′＝mean(λ_Xi′)(i’＝1，2，3)

径向的峭度系数RK，radialkurtosiscoeffocient，

RK＝K₁

轴向的峭度系数AK，axialkurtosiscoefficient，

AK＝(K₂+K₃)/2

1.5特征的验证方案

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

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

Claims (4)

1.一种扩散峭度张量成像的高阶张量特征参数提取方法，其特征是，包括如下步骤：通过受试者在磁共振扫描仪上采集组织沿多个方向的扩散加权信号，拟合得到反映组织内水分子扩散分布概率密度函数特征的二阶及四阶张量，通过一定的张量分析得到扩散系数和峭度系数两个不同角度的参数；对应张量矩阵的特征值和特征向量，采集组织沿多个方向的扩散加权信号进一步细化为：

1.1数据的采集

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

1.2数据的预处理

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

采用非局部的均值滤波(NonLocalMeansFilter)方法，对每个像素点X_p处理如下式(5-1)：

$NL\left(v\right)\left(X_p\right)={\Σ}_{x_p\∈I}w(X_p,Y_q)v\left(X_p\right)---(5-1)$

这里，I代表待处理图像；v(X_p)代表像素X_p的灰度值；NL(v)(X_p)代表像素X_p的灰度值v(X_P)经过非局部均值滤波的结果；w(X_p，Y_q)代表恢复v(X_p)对应于v(X_q)的权值；其中p和q均代表像素的索引编号：

这里，Z(p)是标准化的常数Z(p)＝∑_qw(p，q)，是使用残差方法估计的噪声标准偏差，h是一个滤波参数，距离d表示为如下：这里N＝card(N_p)＝card(N_q)，和分别代表像素X_p和像素Y_q领域N_p和N_q中的第个像素点，对于灰度级图像，

2.如权利要求1所述的扩散峭度张量成像的高阶张量特征参数提取方法，其特征是对应张量矩阵的特征值和特征向量，张量的拟合方案：

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

$S\left(b\right)=S\left(0\right)\exp (-{\mathrm{bD}}_{app}+\frac{1}{6}{b}^2{D_{app}}^2{K}_{app})---(5-3)$

这里，

$D_{app}=g^{T}Dg={\mathrm{Dg}}^2={\Σ}_{i=1}^3{\Σ}_{j=1}^3{D}_{ij}{g}_i{g}_j---(5-4-1)$

$K_{app}=\frac{{\mathrm{MD}}^2}{{D_{app}}^2}{\mathrm{Wg}}^4={\Σ}_{i=1}^3{\Σ}_{j=1}^3{\Σ}_{k=1}^3{\Σ}_{l=1}^3{W}_{ijkl}{g}_i{g}_j{g}_k{g}_l---\left(5-4-2\right)$

其中，D_app和K_app分别为组织在梯度敏感方向g上的表观扩散系数和表观峭度系数；对应的张量即为D和W，分别对应一个二维三阶的实对称矩阵D和一个四维三阶的实对称矩阵W；MD是对应组织部位的平均扩散系数，由扩散张量矩阵D计算得来的；g依然是扩散敏感梯度方向，g_i、g_j、g_k、g_l是向量g的第i、j、k、1个分量元素，i/j/k/l＝1，2，3，上标T表示转置；

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

$\left\{\begin{array}{c}\ln \[S\left(b_1\right)/S\left(0\right)\]=-b_1{D}_{app}+\frac{1}{6}{b}_1^2{D_{app}}^2{K}_{app}\\ \ln \[S\left(b_2\right)/S\left(0\right)\]=-b_2{D}_{app}+\frac{1}{6}{b}_2^2{D_{app}}^2{K}_{app}\end{array}\right.---\left(5-5\right)$

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

$D_{app}=D_{11}{g}_1^2+D_{22}{g}_2^2+D_{33}{g}_3^2+2(D_{12}{g}_1{g}_2+D_{23}{g}_2{g}_3+D_{13}{g}_1{g}_3)---(5-6)$

$\begin{array}{c}{\mathrm{Wg}}^4=W_{1111}{g}_1^4+W_{2222}{g}_2^4+W_{3333}{g}_3^4+4W_{1112}{g}_1^3{g}_2+4W_{1113}{g}_1^3{g}_3+4W_{2221}{g}_2^3{g}_1+4W_{2223}{g}_2^3{g}_3\\ +4W_{3331}{g}_3^3{g}_1+4W_{3332}{g}_3^3{g}_2+12W_{1123}{g}_1^2{g}_2{g}_3+12W_{2213}{g}_2^2{g}_1{g}_3+12W_{3312}{g}_3^2{g}_1{g}_2\\ +6W_{1122}{g}_1^2{g}_2^2+6W_{1133}{g}_1^2{g}_3^2+6W_{2233}{g}_2^2{g}_3^2\end{array}---\left(5-7\right)$

D₁₁、D₁₂、D₁₃、D₂₂、D₂₃、D₃₃表示矩阵D的元素，其中D₁₂＝D₂₁、D₂₃＝D₃₂、D₁₃＝D₃₁，W₃₃₂₁与之类似；通过(5-4-1)、(5-4-2)、(5-6)、(5-7)式采用非线性最小二乘算法即可得到最优的张量矩阵D和W，即对应的张量D和W。

3.如权利要求2所述的扩散峭度张量成像的高阶张量特征参数提取方法，其特征是对应张量矩阵的特征值和特征向量，特征提取方案：

张量的分解算法

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

$D={\[v_1{v}_2{v}_3\]}^T\left[\begin{array}{ccc}{L}_1& 0& 0\\ 0& L_2& 0\\ 0& 0& L_3\end{array}\right]\[v_1{v}_2{v}_3\]$

这里，和是对应D的特征值和特征向量，L₁≥L₂≥L₃，且3个两两正交，分别对应三个特征值的方向；其中作为特征值的代表组织内三分方向上的扩散系数大小，且最大特征值对应方向极为组织内神经纤维的主要走向，另两个特征值反映神经纤维髓鞘的特性；

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

$\left\{\begin{array}{c}\begin{array}{cc}\max & {\mathrm{Wx}}^4\end{array}\\ s.t.x^{T}x=1\end{array}\right.$

即得到如下方程组，

$\left\{\begin{array}{c}{\mathrm{Wx}}^3=\mathrm{\λ}x\\ x^{T}x=1\end{array}\right.$

带入x＝(x₁，x₂，x₃)展开如下，

$\left\{\begin{array}{c}{W}_{1111}{x}_1^3+W_{1222}{x}_2^3+W_{1333}{x}_3^3+3W_{1112}{x}_1^2{x}_2+3W_{1113}{x}_1^2{x}_3+3W_{1223}{x}_2^2{x}_3+3W_{1122}{x}_1{x}_2^2\\ +3W_{1133}{x}_1{x}_3^2+3W_{1233}{x}_2{x}_3^2={\mathrm{\λx}}_1\\ W_{2111}{x}_1^3+W_{2222}{x}_2^3+W_{2333}{x}_3^3+3W_{2112}{x}_1^2{x}_2+3W_{2113}{x}_1^2{x}_3+3W_{2223}{x}_2^2{x}_3+3W_{2122}{x}_1{x}_2^2\\ +3W_{2133}{x}_1{x}_3^2+3W_{2233}{x}_2{x}_3^2={\mathrm{\λx}}_2\\ W_{3111}{x}_1^3+W_{3222}{x}_2^3+W_{3333}{x}_3^3+3W_{3112}{x}_1^2{x}_2+3W_{3113}{x}_1^2{x}_3+3W_{3223}{x}_2^2{x}_3+3W_{3122}{x}_1{x}_2^2\\ +3W_{3133}{x}_1{x}_3^2+3W_{3233}{x}_2{x}_3^2={\mathrm{\λx}}_3\\ x_1^2+x_2^2+x_3^2=1\end{array}\right.$

其解的情况如下，

如果W₂₁₁₁＝W₃₁₁₁＝0，即x₂＝x₃＝0，则λ＝W₁₁₁₁，且对应特征向量为x＝(±1，0，0)；若x₂≠0，x₃＝0，则另t＝x₁/x₂，消元得到如下方程，

$\left\{\begin{array}{c}-W_{2111}{t}^4+(W_{1111}-3W_{2112})t^3+3(W_{1112}-W_{2122})t^2+(3W_{1122}-W_{2222})t+W_{1222}=0\\ W_{3111}{t}^3+3W_{3112}{t}^2+3W_{3122}t+W_{3222}=0\end{array}\right.$

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

$x=\±\frac{1}{\sqrt{t^2+1}}{(t,1,0)}^T$

λ＝Wx⁴

若x₃≠0，则另u＝x₁/x₃，v＝x₂/x₃，消元得到如下方程，

$\left\{\begin{array}{c}-W_{3111}{u}^4-3W_{3112}{u}^{3}v+\left(W_{1111}-3W_{3113}\right)u^3-3W_{3122}{u}^2{v}^2+3\left(W_{1112}-6W_{3123}\right)u^{2}v\\ +\left(3W_{1113}-3W_{3133}\right)u^2-3W_{3223}{\mathrm{uv}}^2-W_{3222}{\mathrm{uv}}^3+3W_{1122}{\mathrm{uv}}^2+\left(6W_{1123}-3W_{3233}\right)uv+\\ +\left(3W_{1133}-W_{3333}\right)u+W_{1222}{v}^3+3W_{1223}{v}^2+3W_{1233}v+W_{1333}=0\\ -W_{3111}{u}^{3}v+W_{2111}{u}^3-3W_{3112}{u}^2{v}^2+\left(3W_{2112}-3W_{3113}\right)u^{2}v+3W_{2113}{u}^2-3W_{3122}{\mathrm{uv}}^3\\ +\left(3W_{2122}-6W_{3123}\right){\mathrm{uv}}^2+\left(6W_{2123}-3W_{3133}\right)uv+3W_{2133}u+3W_{2223}{v}^2-W_{3222}{v}^4\\ +\left(W_{2222}-3W_{3223}\right)v^3-3W_{3233}{v}^2+\left(3W_{2233}-W_{3333}\right)v+W_{2333}=0\end{array}\right.$

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

$x=\frac{\±{(u,v,1)}^T}{\sqrt{u^2+v^2+1}}$

λ＝Wx⁴

由此得到W的特征向量x_w和特征值λ_w，1≤w≤n，且3≤n≤13代表特征值的个数。

4.如权利要求3所述的扩散峭度张量成像的高阶张量特征参数提取方法，其特征是对应张量矩阵的特征值和特征向量，特征参数计算：

根据5-4-2的峭度张量特征分解结果给出峭度张量的提取特征，

平均峭度系数MK(meankurtosiscoefficient)：

MK＝mean(λ_w)，

峭度各向异性KA(kurtosisanisotropy)：

$KA=\sqrt{\frac{n}{n-1}}\·\sqrt{\frac{{\Σ}_{i=1}^n{\left({\mathrm{\λ}}_w-MK\right)}^2}{{\Σ}_{i=1}^n{\mathrm{\λ}}_w^2}}$

对每个w，将W的特征向量x_w按着三个正交方向v_w’进行分组，即寻找w其对应的最优w’，

得到三组向量V₁、V₂、V₃，分别对每组V_w’向量对应的λ_Vw’求平均，则

K_w‘＝mean(λ_Vw’)

径向的峭度系数RK(radialkurtosiscoefficient)：RK＝K₁

轴向的峭度系数AK(axialkurtosiscoefficient)：AK＝(K₂+K₃)/2。