[发明专利]一种利用时空信息的磁共振动脉自旋标记序列部分容积校正方法

权利要求书

1.一种利用时空信息的磁共振动脉自旋标记序列的部分容积校正方法，其特征在于，包括以下步骤：

(1)采集被试者的磁共振数据，包括结构像和动脉自旋标记序列；结构像为T1序列；

(2)将结构像和动脉自旋标记序列进行配准；

所述将结构像和动脉自旋标记序列进行配准，是利用SPM提供的MNI坐标系作为中间值来进行配准；

使用SPM读入结构像和动脉自旋标记数据时，会产生图像坐标系与MNI坐标系之间的转换矩阵，MNI坐标对应的是标准模板，实现结构像和动脉自旋标记数据的配准；

(3)利用SPM软件对结构像进行分割，获取灰质、白质和脑脊液的概率分布图像；

(4)利用线性回归方法对动脉自旋标记序列进行空间部分容积校正，并利用校正结果作为最大期望算法的初始值，在时间方向上进行部分容积校正，充分利用动脉自旋标记序列的空间和时间信息，获得准确的校正结果；

步骤(4)所述的利用校正结果作为最大期望值算法的初始值，在时间方向上进行部分容积校正是利用动脉自旋标记的四维数据，沿着时间轴，将具有相同三维空间坐标i的体素构成一个时间向量；假定该向量中所有元素相互独立，且其中的灰质和白质符合高斯分布，利用最大期望算法估计混合组织中灰质和白质灌注信号，具体求解方法如下：

对于某一空间体素i，磁矩表示为：

ΔM_i＝P_iGMΔM_iGM+P_iWMΔM_iWM(1)

其中，P_iGM和P_iWM是混合体素中灰质和白质的概率；ΔM_iGM和ΔM_iWM代表灰质和白质的灌注信号，灌注信号是通过label/control图像的差值ΔM来描述，且某一空间体素的脑血流与ΔM和M0图像的比值ΔM/M₀等价；根据动脉自旋标记成像的两腔室模型得到：f_tissue＝(ΔM/M₀)·F_tissue ；

其中，F_tissue是与血脑相关的血流灌注系数；

M0图像是与灌注图像序列分别扫描的，采用T1序列的形式，并通过与灌注图像配准，转换成与灌注图像相同大小；故假设M0图像未受到部分容积效应的影响，得到如下关系式：

$CBF=f_{GM}^P+f_{WM}^P=P_{GM}\left(\frac{{\mathrm{\ΔM}}_{iGM}}{M_0}\right)F_{GM}+P_{WM}\left(\frac{{\mathrm{\ΔM}}_{iWM}}{M_0}\right)F_{WM}---\left(2\right)$

其中，P_GM和P_WM通过动脉自旋标记序列与同一被试的结构图像进行配准后获得，F_GM和F_WM与图像的成像参数有关；ΔM_iGM和ΔM_iWM为灰质和白质的灌注信号；沿时间轴，将具有相同三维空间坐标i的体素构成一个时间向量{Y_it,t＝1,…,T}，T代表了时间向量的维数；在观察值Y_it中，灰质和白质的分量表示为：

Y_it＝X_itGM+X_itWM(3)

其中，X_itGM和X_itWM分别是以均数为和方差为和的随机变量，假设所有的T个体素是相互独立的，则有：

$p\left(X\right|{\stackrel{\‾}{X}}_{iGM},{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iGM}^2,{\mathrm{\σ}}_{iWM}^2)=\underset{t=1}{\overset{T}{\Π}}\left\{p\left(X_{itGM}\right|{\stackrel{\‾}{X}}_{iGM},{\mathrm{\σ}}_{iGM}^2)p\left(X_{itWM}\right|{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iWM}^2)\right\}---\left(4\right)$

将灌注模型和式(4)结合，得到：

$\begin{array}{cccc}{\stackrel{\‾}{X}}_{iGM}=R_{iGM}{\mathrm{\ΔM}}_{iGM};& {\stackrel{\‾}{X}}_{iWM}=P_{iWM}{\mathrm{\ΔM}}_{iWM};& {\mathrm{\σ}}_{iGM}^2=P_{iGM}{S}_{iGM};& {\mathrm{\σ}}_{iWM}^2=P_{iWM}{S}_{iWM};\end{array}$

其中，S_iGM和S_iWM分别代表了灰质和白质灌注信号的方差；

假设X_itGM和X_itWM服从高斯分布，则公式(4)变换为:

$\begin{array}{c}p\left(X\right|{\mathrm{\ΔM}}_{iGM},{\mathrm{\ΔM}}_{iWM},S_{iGM},S_{iWM})=\\ \underset{t=1}{\overset{T}{\Π}}\left\{\left(\frac{1}{\sqrt{2{\mathrm{\πP}}_{iGM}{S}_{iGM}}}{e}^{-\frac{{(X_{itGM}-P_{iGM}{\mathrm{\ΔM}}_{iGM})}^2}{2P_{iGM}{S}_{iGM}}}\right)\left(\frac{1}{\sqrt{2{\mathrm{\πP}}_{iWM}{S}_{iWM}}}{e}^{-\frac{{(X_{itWM}-P_{iWM}{\mathrm{\ΔM}}_{iWM})}^2}{2P_{iWM}{S}_{iWM}}}\right)\right\}\end{array}---\left(5\right)$

在最大期望算法中，第t个体素的观察值Y_it是一个不完整的随机变量；X_itGM和X_itWM代表的是第t个体素中完整的混合组织信息，是一个完整变量；

在以公式(1)为条件进行积分，建立不完整变量{Y_it}与完整变量{X_itGM}和{X_itWM}之间的概率分布的关系，如下式：

$\begin{array}{c}p\left(Y_{it}\right|{\mathrm{\ΔM}}_{iGM},{\mathrm{\ΔM}}_{iWM},S_{iGM},S_{iWM})\\ ={\∫}_{\{T_{it}=X_{itGM}+X_{itWM}\}}p\left(X_{itGM}\right|{\stackrel{\‾}{X}}_{iGM},{\mathrm{\σ}}_{iGM}^2)p\left(X_{itWM}\right|{\stackrel{\‾}{X}}_{iWM},{\mathrm{\σ}}_{iWM}^2)dX\end{array}---\left(6\right)$

其中，i代表某一体素三维空间位置，{Y_it}代表具有相同空间位置i的体素所构成一个时间向量；{X_itGM}和{X_itWM}分别为观察值{Y_it}中灰质和白质的分量。

2.根据权利要求1所述的利用时空信息的磁共振动脉自旋标记序列的部分容积校正方法，其特征在于，式(6)的最优解获取步骤如下：

采用最大期望算法求取完全变量条件概率分布的最大期望，E-step是用于对变量的对数似然估计，M-step用于求取期望最大值：

E-step：计算条件概率期望值p(X|Θ)，其中条件期望表示为：

$\begin{array}{c}Q\left(\mathrm{\Θ}\right|{\mathrm{\Θ}}^{\left(n\right)})=E_{Y_{it}=X_{itGM}+X_{itWM}}\[\ln \left(p\right(X|\mathrm{\Θ}\left)\right)|Y,{\mathrm{\Θ}}^{\left(n\right)}\]\\ =E_{Y_{it}=X_{itGM}+X_{itWM}}\[-\frac{1}{2}\underset{t}{\Σ}\{\ln \left(2{\mathrm{\πP}}_{iGM}{S}_{iGM}\right)+\frac{1}{P_{iGM}{S}_{iGM}}\[X_{itGM}^2-2P_{iGM}{\mathrm{\ΔM}}_{iGM}{X}_{itGM}+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}\right)}^2\]\}|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]+\\ E_{Y_{it}=X_{itGM}+X_{itWM}}\[-\frac{1}{2}\underset{t}{\Σ}\{\ln \left(2{\mathrm{\πP}}_{iWM}{S}_{iWM}\right)+\frac{1}{P_{iWM}{S}_{iWM}}\[X_{itWM}^2-2P_{iWM}{\mathrm{\ΔM}}_{iWM}{X}_{itWM}+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}\right)}^2\]\}|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]\\ =-\frac{1}{2}\underset{t}{\Σ}\left\{\begin{array}{c}\ln \left(2{\mathrm{\πP}}_{iGM}{S}_{iGM}\right)+\frac{1}{P_{iGM}{S}_{iGM}}\[E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGM}^2\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})-2P_{iGM}{\mathrm{\ΔM}}_{iGM}{E}_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}\right)}^2\]+\\ \ln \left(2{\mathrm{\πP}}_{iWM}{S}_{iWM}\right)+\frac{1}{P_{iWM}{S}_{iWM}}\[E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}^2\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})-2P_{iWM}{\mathrm{\ΔM}}_{iWM}{E}_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}\right)}^2\]\end{array}\right\}\end{array}---\left(7\right)$

根据条件期望的推导，得出：

$\begin{array}{c}{X}_{itGM}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itGM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})\\ =P_{iGM}{\mathrm{\ΔX}}_{iGM}^{\left(n\right)}+\frac{P_{iGM}{S}_{iGM}^{\left(n\right)}}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}\[Y_{it}-(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)})\]\end{array}---\left(8\right)$

$\begin{array}{c}{X}_{itWM}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\left(X_{itWM}\right|Y_{it},{\mathrm{\Θ}}^{\left(n\right)})\\ =P_{iWM}{\mathrm{\ΔX}}_{iWM}^{\left(n\right)}+\frac{P_{iWM}{S}_{iWM}^{\left(n\right)}}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}\[Y_{it}-(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)})\]\end{array}---\left(9\right)$

${\left(X_{itGM}^2\right)}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\[X_{itGMt}^2|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]={\left(X_{itGM}^{\left(n\right)}\right)}^2+\frac{\left(P_{iGM}{S}_{iGM}^{\left(n\right)}\right)\left(P_{iWM}{S}_{iWM}^{\left(n\right)}\right)}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}---\left(10\right)$

${\left(X_{itWM}^2\right)}^{\left(n\right)}=E_{Y_{it}=X_{itGM}+X_{itWM}}\[X_{itWM}^2|Y_{it},{\mathrm{\Θ}}^{\left(n\right)}\]={\left(X_{itWM}^{\left(n\right)}\right)}^2+\frac{\left(P_{iGM}{S}_{iGM}^{\left(n\right)}\right)\left(P_{iWM}{S}_{iWM}^{\left(n\right)}\right)}{P_{iGM}{S}_{iGM}^{\left(n\right)}+P_{iWM}{S}_{iWM}^{\left(n\right)}}---\left(11\right)$

M-step：通过n+1次迭代使条件概率期望值最大化，混合组织模型中的均值通过最大化条件概率的n+1次迭代求出：

$\frac{\∂Q}{\∂{\mathrm{\ΔM}}_{iGM}}{|}_{{\mathrm{\ΔM}}_{iGM}={\mathrm{\ΔM}}_{iGM}^{(n+1)}}=0\⇒{\mathrm{\ΔM}}_{iGM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T{X}_{itGM}^{\left(n\right)}}{T\·P_{iGM}}---\left(12\right)$

$\frac{\∂Q}{\∂{\mathrm{\ΔM}}_{iWM}}{|}_{{\mathrm{\ΔM}}_{iWM}={\mathrm{\ΔM}}_{iWM}^{(n+1)}}=0\⇒{\mathrm{\ΔM}}_{iWM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T{X}_{itWM}^{\left(n\right)}}{T\·P_{iWM}}---\left(13\right)$

$S_{iGM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T\[{\left(X_{itGM}^2\right)}^{\left(n\right)}-2X_{itGM}^{\left(n\right)}{P}_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}+{\left(P_{iGM}{\mathrm{\ΔM}}_{iGM}^{\left(n\right)}\right)}^2\]}{T\·P_{iGM}}---\left(14\right)$

$S_{iWM}^{(n+1)}=\frac{{\mathrm{\Σ}}_{t=1}^T\[{\left(X_{itWM}^2\right)}^{\left(n\right)}-2X_{itWM}^{\left(n\right)}{P}_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)}+{\left(P_{iWM}{\mathrm{\ΔM}}_{iWM}^{\left(n\right)}\right)}^2\]}{T\·P_{iWM}}---\left(15\right).$