CN102005031B - 一种mri系统中k空间采样数据的运动伪影消除方法及装置 - Google Patents

Abstract

本发明涉及一种MRI系统中K空间采样数据的运动伪影消除方法及装置，该方法包括数据采集步骤；数据预处理步骤，保证回波信号强度的一致性，并将回波最大值移到中心；旋转运动估计步骤，用频域相似性寻找旋转数据带和基础数据带之间的回波参考对，计算参考对之间的角度相对偏移量，建立旋转运动方程组并引入约束条件求解旋转运动参数；平移运动估计步骤，根据数据一致性原理等量方程和上步骤中所得旋转角度，建立平移运动方程组求解平移运动参数；运动参数补偿步骤，根据步骤三、四中求出的运动参数进行运动补偿；滤波反投影加权重建步骤，采用滤波反投影方法加权重建本发明提供的方法或装置能有效抑制运动伪影。

Description

一种MRI系统中K空间采样数据的运动伪影消除方法及装置

技术领域

本发明涉及一种MRI系统中K空间采样数据的运动伪影消除方法及装置。

背景技术

磁共振成像(Magnetic Resonance Imaging，MRI)是当今最先进的医学成像方法之一，在临床上和科学研究中得到了越来越广泛的应用。磁共振成像过程中，病人的自主性和生理性运动常常难以避免，它将破坏数据的采集过程，并在所成的图像上形成伪影，使医生难以做出正确的诊断。由于运动的不确定性以及难以获得运动的先验知识，对运动伪影的校正也就十分困难，严重影响和阻碍了磁共振成像技术的发展和应用。因此，对磁共振运动伪影消除方法的研究，引起了国内外学者的广泛关注，是当前医学成像领域中的研究热点和重要难题之一。

由于病人运动对采样数据的破坏，通过后处理的方式对传统采样方式获取的MR数据进行运动伪影校正是非常困难的，因此研究人员把刚性运动分解成独立的平移和旋转运动进行校正。平移运动伪影的校正主要有基于凸集投影的校正算法、基于运动熵的校正算法和基于傅立叶投影的校正算法。旋转伪影校正主要有基于图像能量和基于磁共振点扩散函数的校正算法。但是病人运动往往同时包含平移和旋转运动，这些算法都不能有效的解决问题。直到1999年，James G.Pipe提出按PROPELLER(Periodically RotatedOverlapping Parallel Lines with Enhanced Reconstruction)方式采样磁共振数据和成像算法，才在磁共振成像过程中可以有效的消除病人刚体运动所形成的伪影。但PROPELLER采样方式需要过采样数据，延长了采集时间。另外，PROPELLER采样数据的重建过程中，需要将大量的位于不规则坐标系上的数据点网格化到笛卡尔坐标系下，以便于快速傅里叶重建，网格化过程需要消耗大量的计算时间。

发明内容

本发明所要解决的第一个技术问题是针对上述现有技术提供一种MRI系统中K空间采样数据的运动伪影消除方法，该方法能有效抑制运动伪影，避免了数据过采样导致采集时间长、数据量大以及重建过程计算代价大的问题。

本发明所要解决的第二个技术问题是针对上述现有技术提供一种MRI系统中K空间采样数据的运动伪影消除装置，运用该装置能有效抑制运动伪影，避免了数据过采样导致采集时间长、数据量大以及重建过程计算代价大的问题。

本发明解决上述第一个技术问题所采用的技术方案为：该MRI系统中K空间采样数据的运动伪影消除方法，其特征在于：包括以下步骤：

步骤一、数据采集步骤

旋转采集K空间数据，获得具有放射状的K空间数据，具体方法为：

(1-1)、首先从零度开始以Δθ＝2π/N的间隔均匀获取共

个回波，称为基础数据，记为M_b；

(1-2)、然后从开始，仍以Δθ＝2π/N的间隔均匀获取另外

个回波，称为旋转数据，记为M_r；

基础数据和旋转数据共同组成本步骤采集的K空间原始采样数据；

步骤二、数据预处理步骤

(2-1)、将步骤一采集到的K空间原始采样数据的每个回波做一维傅里叶变换，然后将变换后的数据取模，即得到投影数据；

(2-2)、接着根据投影数据中任意一行回波数据中所有采样点对应的投影数据的总和相同的特征归一化经过步骤(2-1)处理后的投影数据，去除不同回波信号强度不均的影响，校正公式为：其中校正后的投影数据为P，P₀为经过步骤(2-1)处理后得到的投影数据中任意一行回波中任意一个采样点对应的投影数据，l为经过步骤(2-1)处理后得到的投影数据中任意一行回波的采样点数，将基础数据校正后的投影数据记为P_b，将旋转数据校正后的投影数据记为P_r；

(2-3)、对旋转数据校正后的投影数据P_r在角度方向进行T倍上采样得到旋转数据校正后的T倍采样投影数据P_r0；

(2-4)、最后，将基础数据校正后的投影数据为P_b和旋转数据校正后的T倍采样投影数据P_r0做一维傅里叶逆变换并取模，从而得到本步骤预处理后的基础数据F_b和旋转数据F_r；

步骤三、旋转运动估计步骤

(3-1)、利用数据相关性在经过步骤二预处理后的基础数据F_b及经过步骤二预处理后的旋转数据F_r中寻找回波参考对，一对回波参考对是指在基础数据和旋转数据中相同旋转角度位置所采集到的一对回波，当扫描过程中被检测物体发生旋转运动时，回波参考对在K空间位置会发生偏移：

对于经过步骤二预处理后的基础数据F_b中的第n行数据F_b(k，n)，k是每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，n为自然数，其取值从1到经过步骤二预处理后的基础数据F_b数据的总行数；通过如下公式计算与经过步骤二预处理后的旋转数据F_r中第m行数据F_r(k，m)的相关性C(n，m)，k为每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，m为自然数，其取值从1到经过步骤二预处理后的旋转数据F_r数据的总行数：

$C(n,m)=\frac{\underset{k=1}{\overset{S}{\mathrm{\Σ}}}(F_b(k,n)-{\stackrel{\‾}{F}}_b\left(n\right))(F_r(k,m)-{\stackrel{\‾}{F}}_r\left(m\right))}{{\mathrm{\σ}}_n{\mathrm{\σ}}_m}$

其中C(n，m)表示两个回波之间的相关性，k为每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，分别为F_b(k，n)和F_r(k，m)的均值，σ_n、σ_m分别为F_b(k，n)和F_r(k，m)的方差，S为单个回波最大的采样点数；

经过步骤二预处理后的旋转数据F_r中第m行数据F_r(k，m)如果与经过步骤二预处理后的基础数据F_b中第n行数据F_b(k，n)具有最大相关性，即C(n，m)最大，则它们构成一对回波参考对；

(3-2)、计算经过步骤二预处理后的基础数据F_b及经过步骤二预处理后的旋转数据F_r中回波参考对之间的旋转角度偏移Δθ，然后根据这个旋转角度偏移Δθ，列出被检测物体的旋转运动方程组：

设回波参考对F_b(k，n)和F_r(k，m)的采样时刻分别为t_n和t_m，用φ(t)描述被检测物体的旋转运动方程，可得到如下等式：

φ(t_m)-φ(t_n)＝Δθ

被检测物体的旋转运动方程φ(t)用采样时刻t的Q阶多项式表示为：

φ(t)＝r₁t+r₂t²+…+r_Qt^Q其中，r₁到r_Q为关于t的Q阶多项式的系数；

这样，就可得到包含Q个旋转运动参数的方程：

r₁(t_m-t_n)+r₂(t_m ²-t_n ²)+…+r_Q(t_m ^Q-t_n ^Q)＝Δθ

设可以在经过步骤二预处理后的基础数据F_b及经过步骤二预处理后的旋转数据F_r中总共能找到M对回波参考对，且因此共有M个方程，将上述方程组写成矩阵形式为：

AR＝G

其中

R＝(r₁，r₂，…，r_Q)^T，T表示矩阵转置；

G＝(Δθ₁，Δθ₂，…，Δθ_M)^T，T表示矩阵转置，Δθ_M为第M对回波参考对之间的旋转角度偏移；

(3-3)、然后根据被检测物体在检测过程中可能的最大旋转角度作为约束条件，利用二次规划的方法求解被检测物体的旋转向量；

根据受检物体在检测过程中可能的最大旋转角度作为约束条件，将被检测物体的旋转运动方程组转化为二次规划问题：

$\left\{\begin{array}{cc}\min & \frac{1}{2}{R}^T\mathrm{HR}+\mathrm{fR}\\ \mathrm{subject\; to}& \left|\mathrm{\φ}\left(t\right)\right|\≤{\mathrm{\φ}}_{\max }\end{array}\right.$

上式中，H＝2A^TA，f＝-2G^TA，R＝(r₁，r₂，…，r_Q)^T，φ_max为被检测物体在检测过程中可能的最大旋转角度；利用二次规划方法求解上式可得到旋转参数R，从而得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]，N为步骤一中提到的采集到的回波总个数；

步骤四、平移运动估计步骤

(4-1)、根据数据一致性原理(Helgason-Ludwig Consistency Condition，HLCC)，被检物体发生运动后存在如下等量关系：

$\underset{-\∞}{\overset{\∞}{\∫}}{(l-d_1\cos \left({\mathrm{\θ}}_c\left(n\right)\right)-d_2\sin \left({\mathrm{\θ}}_c\left(n\right)\right))}^{k}P(n,l)\mathrm{dl}=\underset{r=0}{\overset{k}{\mathrm{\Σ}}}\left(\begin{array}{c}k\\ r\end{array}\right)m_{r,k-r}{\cos }^r\left({\mathrm{\θ}}_c\left(n\right)\right){\sin }^{k-r}\left({\mathrm{\θ}}_c\left(n\right)\right)$

上式中，l是经过步骤(2-1)得到的某一投影数据的位置坐标，d₁是被检测物体在x方向平移运动分量，d₂是受检物体在y方向平移运动分量，θ_c为利用步骤三中得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]进行修正后的旋转角度，θ_c等于K空间旋转向量减去步骤三中得到的被检测物体的旋转向量，k为HLCC等式的阶数，可取k＝0，1，2......，P(n，l)为步骤二中(2-2)校正后的第n行投影数据，m_r，k-r为HLCC定义的被检测物体被检断层图象的几何动量，对于给定的一个受检物体，m_r，k-r为常量；

(4-2)、将物体的平移运动在x方向和y方向的分量用如下关于采样时刻t的Q阶多项式描述：

d_i(t)＝a_i1t+a_i2t²+…+a_iQt^Q，其中参数a_i1到a_iQ为多项式系数；

将上式代入(4-1)中的公式，其中θ_c为利用步骤三中得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]进行修正后的旋转角度，并取k＝1，得到

$\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}{C}_{1q}\left(t\right)a_{1q}+\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}{C}_{2q}\left(t\right)a_{2q}+m_{\mathrm{0,1}}\sin \left({\mathrm{\θ}}_c\left(n\right)\right)+m_{\mathrm{1,0}}\cos \left({\mathrm{\θ}}_c\left(n\right)\right)=C_0\left(t\right)$

其中 $\left\{\begin{array}{c}{C}_0\left(t\right)=\∫l\×P(n,l)\mathrm{dl}\\ C_{\mathrm{lq}}\left(t\right)=\cos \left({\mathrm{\θ}}_c\left(n\right)\right)\×t^q\×P(n,l)\mathrm{dl}\\ C_{2q}\left(t\right)=\sin \left({\mathrm{\θ}}_c\left(n\right)\right)\×t^q\×\∫P(n,l)\mathrm{dl}\end{array}\right.$

对于步骤一中采集的每个回波都可以得到一个方程，共N个方程，利用最小二乘可以求出平移运动分量的多项式参数a_1q，a_2q，从而得到被检测物体的平移运动向量；

步骤五、运动参数补偿步骤

(5-1)、将K空间旋转向量减去步骤三中得到的被检测物体的旋转向量，得到修正后的K空间旋转向量θ_c；

(5-2)、被检测物体的旋转运动导致了投影数据在角度方向上分布不均，在角度方向进行加权补偿，首先对θ_c进行升序排序得到θ_co，然后计算各个旋转角度对应的加权因子

其中W₁(n)为每个角度对应的加权因子；

(5-3)、对步骤一中采集的基础数据和旋转数据补偿一个与步骤四得到的被检测物体的平移动运动向量相关的相位偏移来消除平移的影响，补偿公式为M′(n，k)＝e^{i×2π×k×d(n)}×M(n，k)，M(n，k)为步骤一中采集到的原始回波数据，M′(n，k)即为相位补偿后的回波数据，k为每个回波上的采样点的编号，d为被检测物体平移运动在K空间径向方向的运动分量，其数学表达式为：

d(n)＝d₁cos(θ_c(n))+d₂sin(θ_c(n))；

(5-4)、将相位补偿后的回波数据M′(n，k)中每个回波做一维傅里叶变换，然后将变换后的数据取模，得到投影数据P₁；

(5-5)、根据投影数据中任意一行回波数据中所有采样点对应的投影数据的总和相同的特征归一化经过步骤(5-4)处理后的投影数据，校正公式为

其中校正后的投影数据为P_c，P₁为经过步骤(5-4)处理后得到的投影数据中任意一个回波中任意一个采样点对应的投影数据，l为经过步骤(5-4)处理后得到的投影数据中任意一行回波的采样点数；

步骤六、滤波反投影加权重建步骤

(6-1)、计算投影向量加权因子，进一步减少未能校正到的运动对图像的影响，C为权重调节因子，C∈[0，∞]；

$W_2\left(n\right)=\left\{\begin{array}{cc}\frac{1}{\mathrm{\π}}\mathrm{\θ}\left(n\right)+C& n\≤\frac{1}{2}N\\ 1-\frac{1}{\mathrm{\π}}\mathrm{\θ}\left(n\right)+C& n>\frac{1}{2}N\end{array}\right.$

上式中，W₂(n)投影向量加权因子，θ(n)即为第n行回波数据采集时的旋转角度；

(6-2)、对步骤五中经过步骤(5-5)校正后的投影数据P_c进行加权

P_w(n，l)＝W₁(n)×W₂(n)×P_c(n，l)

其中P_W(n，l)为经过加权后的投影数据，W₁(n)为第n行回波数据对应的加权因子，

W₂(n)为投影向量加权因子，P_c(n，l)为经过步骤(5-5)处理后的第n行回波的投影数据；

(6-3)、以加权后的投影数据P_w作为重建数据，采用滤波反投影方法重建即可得到校正后的图像。

本发明解决上述第二个技术问题所采用的方法为：该MRI系统中K空间采样数据的运动伪影消除装置，其特征在于：包括

数据采集装置，用于旋转采集K空间数据，获得具有放射状的K空间数据，具体方法为：

(1-1)、首先从零度开始以Δθ＝2π/N的间隔均匀获取共

个回波，称为基础数据，记为M_b；

(1-2)、然后从

Δθ开始，仍以Δθ＝2π/N的间隔均匀获取另外

个回波，称为旋转数据，记为M_r；

基础数据带和旋转数据带共同组成本数据采集装置采集的K空间原始采样数据；

数据预处理装置，用于以下功能的处理：

(2-1)、将数据采集装置采集到的K空间原始采样数据的每个回波做一维傅里叶变换，然后将变换后的数据取模，即得到投影数据；

(2-2)、接着根据投影数据中任意一行回波数据中所有采样点对应的投影数据的总和相同的特征归一化经过步骤(2-1)处理后的投影数据，去除不同回波信号强度不均的影响，校正公式为：

其中校正后的投影数据为P，P₀为经过步骤(2-1)处理后得到的投影数据中任意一行回波中任意一个采样点对应的投影数据，l为经过步骤(2-1)处理后得到的投影数据中任意一行回波的采样点数，将基础数据校正后的投影数据记为P_b，将旋转数据校正后的投影数据记为P_r；

(2-3)、对旋转数据校正后的投影数据P_r在角度方向进行T倍上采样得到旋转数据校正后的T倍采样投影数据P_r0；

(2-4)、最后，将基础数据校正后的投影数据为P_b和旋转数据校正后的T倍采样投影数据P_r0做一维傅里叶逆变换并取模，从而得到本步骤预处理后的基础数据F_b和旋转数据F_r；

旋转运动估计装置，用于以下功能的处理：

(3-1)、利用数据相关性在经过数据预处理装置预处理后的基础数据F_b及经过数据预处理装置预处理后的旋转数据F_r中寻找回波参考对，一对回波参考对是指在基础数据和旋转数据中相同旋转角度位置所采集到的一对回波，当扫描过程中被检测物体发生旋转运动时，回波参考对在K空间位置会发生偏移：

对于经过数据预处理装置预处理后的基础数据F_b中的第n行数据F_b(k，n)，k是每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，θ(n)为是第n个回波被采集时的旋转角度，n为自然数，其取值从1到经过数据预处理装置预处理后的基础数据F_b数据的总行数；通过如下公式计算与经过数据预处理装置预处理后的旋转数据F_r中第m行数据F_r(k，m)的相关性C(n，m)，k为每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，θ(m)为第m个回波被采集时的旋转角度，m为自然数，其取值从1到经过数据预处理装置预处理后的旋转数据F_r数据的总行数；：

$C(n,m)=\frac{\underset{k=1}{\overset{S}{\mathrm{\Σ}}}(F_b(k,n)-{\stackrel{\‾}{F}}_b\left(n\right))(F_r(k,m)-{\stackrel{\‾}{F}}_r\left(m\right))}{{\mathrm{\σ}}_n{\mathrm{\σ}}_m}$

其中C(n，m)表示两个回波之间的相关性，k为每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，分别为F_b(k，n)和F_r(k，m)的均值，σ_n、σ_m分别为F_b(k，n)和F_r(k，m)的方差，S为单个回波最大的采样点数；

经过数据预处理装置预处理后的旋转数据F_r中第m行数据F_r(k，m)如果与经过数据预处理装置预处理后的基础数据F_b中第n行数据F_b(k，n)具有最大相关性，即C(n，m)最大，则它们构成一对回波参考对；

(3-2)、计算经过数据预处理装置预处理后的基础数据F_b及经过数据预处理装置预处理后的旋转数据F_r中回波参考对之间的旋转角度偏移Δθ，然后根据这个旋转角度偏移Δθ，列出被检测物体的旋转运动方程组：

设回波参考对F_b(k，n)和F_r(k，m)的采样时刻分别为t_n和t_m，用φ(t)描述被检测物体的旋转运动方程，可得到如下等式：

φ(t_m)-φ(t_n)＝Δθ

被检测物体的旋转运动方程φ(t)用采样时刻t的Q阶多项式表示为：

φ(t)＝r₁t+r₂t²+…+r_Qt^Q其中，r₁到r_Q为关于t的Q阶多项式的系数；

这样，就可得到包含Q个旋转运动参数的方程：

r₁(t_m-t_n)+r₂(t_m ²-t_n ²)+…+r_Q(t_m ^Q-t_n ^Q)＝Δθ

设可以在经过数据预处理装置预处理后的基础数据F_b及经过数据预处理装置预处理后的旋转数据F_r中总共能找到M对回波参考对，且因此共有M个方程，将上述方程组写成矩阵形式为：

AR＝G

其中

R＝(r₁，r₂，…，r_Q)^T，T表示矩阵转置；

G＝(Δθ₁，Δθ₂，…，Δθ_M)^T，T表示矩阵转置，Δθ_M为第M对回波参考对之间的旋转角度偏移；

(3-3)、然后根据被检测物体在检测过程中可能的最大旋转角度作为约束条件，利用二次规划的方法求解被检测物体的旋转向量；

根据被检测物体在检测过程中可能的最大旋转角度作为约束条件，将被检测物体的旋转运动方程组转化为二次规划问题

$\left\{\begin{array}{cc}\min & \frac{1}{2}{R}^T\mathrm{HR}+\mathrm{fR}\\ \mathrm{subject\; to}& \left|\mathrm{\φ}\left(t\right)\right|\≤{\mathrm{\φ}}_{\max }\end{array}\right.$

上式中，H＝2A^TA，f＝-2G^TA，R＝(r₁，r₂，…，r_Q)^T，φ_max为被检测物体在检测过程中可能的最大旋转角度；利用二次规划方法求解上式可得到旋转参数R，从而得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]，N为数据采集装置中提到的采集到的回波总个数；

平移运动估计装置，用于以下功能的处理：

(4-1)、根据HLCC数据一致性原理，被检物体发生运动后存在如下等量关系：

$\underset{-\∞}{\overset{\∞}{\∫}}{(l-d_1\cos \left({\mathrm{\θ}}_c\left(n\right)\right)-d_2\sin \left({\mathrm{\θ}}_c\left(n\right)\right))}^{k}P(n,l)\mathrm{dl}=\underset{r=0}{\overset{k}{\mathrm{\Σ}}}\left(\begin{array}{c}k\\ r\end{array}\right)m_{r,k-r}{\cos }^r\left({\mathrm{\θ}}_c\left(n\right)\right){\sin }^{k-r}\left({\mathrm{\θ}}_c\left(n\right)\right)$

上式中，l是经过步骤(2-1)得到的某一投影数据的位置坐标；d₁是被检测物体在x方向平移运动分量；d₂是受检物体在y方向平移运动分量；θ_c为利用旋转运动估计装置中得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]进行修正后的旋转角度，θ_c等于K空间旋转向量减去旋转运动估计装置中得到的被检测物体的旋转向量；k为HLCC等式的阶数，可取k＝0，1，2......，P(n，l)为旋转运动估计装置中(2-2)校正后的第n行投影数据；m_r，k-r为HLCC定义的被检测物体被检断层图象的几何动量，对于给定的一个受检物体，m_r，k-r为常量；

(4-2)、将物体的平移运动在x方向和y方向的分量用如下关于采样时刻t的Q阶多项式描述：

d_i(t)＝a_i1t+a_i2t²+…+a_iQt^Q  ，其中参数a_i1到a_iQ为多项式系数；

将上式代入(4-1)中的公式，其中θ_c为利用旋转运动估计装置中得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]进行修正后的旋转角度，θ_c等于K空间旋转向量减去旋转运动估计装置中得到的被检测物体的旋转向量，并取k＝1，得到

$\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}{C}_{1q}\left(t\right)a_{1q}+\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}{C}_{2q}\left(t\right)a_{2q}+m_{\mathrm{0,1}}\sin \left({\mathrm{\θ}}_c\left(n\right)\right)+m_{\mathrm{1,0}}\cos \left({\mathrm{\θ}}_c\left(n\right)\right)=C_0\left(t\right)$

其中 $\left\{\begin{array}{c}{C}_0\left(t\right)=\∫l\×P(n,l)\mathrm{dl}\\ C_{\mathrm{lq}}\left(t\right)=\cos \left({\mathrm{\θ}}_c\left(n\right)\right)\×t^q\×P(n,l)\mathrm{dl}\\ C_{2q}\left(t\right)=\sin \left({\mathrm{\θ}}_c\left(n\right)\right)\×t^q\×\∫P(n,l)\mathrm{dl}\end{array}\right.$

对于数据采集装置中采集的每个回波都可以得到一个方程，共N个方程，利用最小二乘可以求出平移运动分量的多项式参数a_1q，a_2q，从而得到被检测物体的平移运动向量；

运动参数补偿装置，用于以下功能的处理：

(5-1)、将K空间旋转向量减去旋转运动估计装置中得到的被检测物体的旋转向量，得到修正后的K空间旋转向量θ_c；

(5-2)、被检测物体的旋转运动导致了投影数据在角度方向上分布不均，在角度方向进行加权补偿，首先对θ_c进行升序排序得到θ_co，然后计算各个旋转角度对应的加权因子

其中W₁(n)为每个角度对应的加权因子；

(5-3)、对数据采集装置中采集的基础数据和旋转数据补偿一个与平移运动估计装置得到的被检测物体的平移动运动向量相关的相位偏移来消除平移的影响，补偿公式为M′(n，k)＝e^{i×2π×k×d(n)}×M(n，k)，M(n，k)为数据采集装置中采集到的原始回波数据，M′(n，k)即为相位补偿后的回波数据，k为每个回波上的采样点的编号，d为被检测物体平移运动在K空间径向方向的运动分量，其数学表达式为：

d(n)＝d₁cos(θ_c(n))+d₂sin(θ_c(n))；

(5-4)、将相位补偿后的回波数据M′(n，k)中每个回波做一维傅里叶变换，然后将变换后的数据取模，得到投影数据P₁；

(5-5)、根据投影数据中任意一行回波数据中所有采样点对应的投影数据的总和相同的特征归一化经过步骤(5-4)处理后的投影数据，校正公式为

其中校正后的投影数据为P_c，P₁为经过步骤(5-4)处理后得到的投影数据中任意一个回波中任意一个采样点对应的投影数据，l为经过步骤(5-4)处理后得到的投影数据中任意一行回波的采样点数；

滤波反投影加权重建装置，用于以下功能的处理：

(6-1)、计算投影向量加权因子，进一步减少未能校正到的运动对图像的影响，C为权重调节因子，C∈[0，∞]；投影向量加权因子W₂(n)表示为

$W_2\left(n\right)=\left\{\begin{array}{cc}\frac{1}{\mathrm{\π}}\mathrm{\θ}\left(n\right)+C& n\≤\frac{1}{2}N\\ 1-\frac{1}{\mathrm{\π}}\mathrm{\θ}\left(n\right)+C& n>\frac{1}{2}N\end{array}\right.$

上式中，θ(n)即为第n行回波数据采集时的旋转角度；

(6-2)、对运动参数补偿装置中经过步骤(5-5)校正后的投影数据P_c进行加权

P_w(n，l)＝W₁(n)×W₂(n)×P_c(n，l)

其中P_W(n，l)为经过加权后的投影数据，W₁(n)为第n行回波数据对应的加权因子，

W₂(n)为投影向量加权因子，P_c(n，l)为经过步骤(5-5)处理后的第n行回波的投影数据；

(6-3)、以加权后的投影数据P_w作为重建数据，采用滤波反投影方法重建即可得到校正后的图像。

与现有技术相比，本发明的优点在于：运用本发明提供的方法或装置能有效抑制运动伪影，避免了数据过采样导致采集时间长、数据量大以及重建过程计算代价大的问题。

附图说明

图1为本发明实施例方法的流程图；

图2为本发明实施例方法步骤一采用的放射状K空间采样模式图；

图3为本发明实施例方法步骤二数据预处理流程图；

图4为基础数据和旋转数据带中的回波参考对的示意图；

图5为本发明实施例实验中旋转角度估计结果示意图；

图6为本发明实施例实验中平移运动估计结果示意图；

图7为本发明实施例试验中直接重建成像示意图；

图8为本发明实施例中采用本发明提供方法进行运动伪影消除后的图像示意图。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

本发明首先提供了一种MRI系统中K空间采样数据的运动伪影消除方法，其包括以下步骤，请参见图1所示：

步骤一、数据采集步骤，基于特殊的放射状K空间数据采集方式采集数据；

步骤二、数据预处理步骤，保证回波信号强度的一致性，并将回波最大值移到中心；

步骤三、旋转运动估计步骤，用频域相似性寻找旋转数据带和基础数据带之间的回波参考对，计算参考对之间的角度相对偏移量，建立旋转运动方程组并引入约束条件求解旋转运动参数；

步骤四、平移运动估计步骤，根据HLCC等量方程和步骤3中所得旋转角度，建立平移运动方程组求解平移运动参数；

步骤五、运动参数补偿步骤，根据步骤三、四中求出的运动参数进行运动补偿；

步骤六、滤波反投影加权重建步骤，采用滤波反投影方法加权重建。

上述步骤一数据采集步骤具体如下：旋转采集K空间数据，获得具有放射状的K空间数据，具体方法为，请参见图2所示。

(1-1)、首先从零度开始以Δθ＝2π/N的间隔均匀获取共

个回波，称为基础数据，记为M_b；

(1-2)、然后从

开始，仍以Δθ＝2π/N的间隔均匀获取另外

个回波，称为旋转数据，记为M_r；

基础数据和旋转数据共同组成本步骤采集的K空间原始采样数据。

本发明实施例中共旋转采集了180条K空间数据，首先采集基础数据带，旋转角度为[0°，2°，4°......176°，178°]，共90条K空间线；接着采集旋转数据带，旋转角度为[1°，3°，5°......177°，179°]。

上述步骤二的数据预处理步骤具体如下，请参看图3所示：

由于梯度系统的非线性，主磁场的不理想以及涡流的影响等因素，采集到的K空间数据存在信号强度不均和最大值不在中心的问题，因此需要

(2-1)、将步骤一采集到的K空间原始采样数据的每个回波做一维傅里叶变换，然后将变换后的数据取模，即得到投影数据；

(2-2)、接着根据投影数据中任意一行回波数据中所有采样点对应的投影数据的总和相同的特征归一化经过步骤(2-1)处理后的投影数据，去除不同回波信号强度不均的影响，校正公式为：

其中校正后的投影数据为P，P₀为经过步骤(2-1)处理后得到的投影数据中任意一行回波中任意一个采样点对应的投影数据，l为经过步骤(2-1)处理后得到的投影数据中任意一行回波的采样点数，将基础数据校正后的投影数据记为P_b，将旋转数据校正后的投影数据记为P_r；

(2-3)、对旋转数据校正后的投影数据P_r在角度方向进行T倍上采样得到旋转数据校正后的T倍采样投影数据P_r0，本实施例中，采样倍数T为4；

(2-4)、最后，将基础数据校正后的投影数据为P_b和旋转数据校正后的T倍采样投影数据P_r0做一维傅里叶逆变换并取模，从而得到本步骤预处理后的基础数据F_b和旋转数据F_r。

上述步骤三旋转运动估计步骤具体如下：

(3-1)、利用数据相关性在经过步骤二预处理后的基础数据F_b及经过步骤二预处理后的旋转数据F_r中寻找回波参考对，一对回波参考对是指在基础数据和旋转数据中相同旋转角度位置所采集到的一对回波，当扫描过程中被检测物体发生旋转运动时，回波参考对在K空间位置会发生偏移：

对于经过步骤二预处理后的基础数据F_b中的第n行数据F_b(k，n)，k是每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，n为自然数，其取值从1到经过步骤二预处理后的基础数据F_b数据的总行数；通过如下公式计算与经过步骤二预处理后的旋转数据F_r中第m行数据F_r(k，m)的相关性C(n，m)，k为每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，m为自然数，其取值从1到经过步骤二预处理后的旋转数据F_r数据的总行数：

$C(n,m)=\frac{\underset{k=1}{\overset{S}{\mathrm{\Σ}}}(F_b(k,n)-{\stackrel{\‾}{F}}_b\left(n\right))(F_r(k,m)-{\stackrel{\‾}{F}}_r\left(m\right))}{{\mathrm{\σ}}_n{\mathrm{\σ}}_m}$

其中C(n，m)表示两个回波之间的相关性，k为每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，

分别为F_b(k，n)和F_r(k，m)的均值，σ_n、σ_m分别为F_b(k，n)和F_r(k，m)的方差，S为单个回波最大的采样点数；

经过步骤二预处理后的旋转数据F_r中第m行数据F_r(k，m)如果与4经过步骤二预处理后的基础数据F_b中第n行数据F_b(k，n)具有最大相关性，即C(n，m)最大，则它们构成一对回波参考对，请参见图4所示；

(3-2)、计算经过步骤二预处理后的基础数据F_b及经过步骤二预处理后的旋转数据F_r中回波参考对之间的旋转角度偏移Δθ，然后根据这个旋转角度偏移Δθ，列出被检测物体的旋转运动方程组：

设回波参考对F_b(k，n)和F_r(k，m)的采样时刻分别为t_n和t_m，用φ(t)描述被检测物体的旋转运动方程，可得到如下等式：

φ(t_m)-φ(t_n)＝Δθ

被检测物体的旋转运动方程φ(t)用采样时刻t的Q阶多项式表示为：

φ(t)＝r₁t+r₂t²+…+r_Qt^Q其中，r₁到r_Q为关于t的Q阶多项式的系数；

这样，就可得到包含Q个旋转运动参数的方程：

r₁(t_m-t_n)+r₂(t_m ²-t_n ²)+…+r_Q(t_m ^Q-t_n ^Q)＝Δθ

设可以在经过步骤二预处理后的基础数据F_b及经过步骤二预处理后的旋转数据F_r中总共能找到M对回波参考对，且

因此共有M个方程，将上述方程组写成矩阵形式为：

AR＝G

其中

R＝(r₁，r₂，…，r_Q)^T，T表示矩阵转置；

G＝(Δθ₁，Δθ₂，…，Δθ_M)^T，T表示矩阵转置，Δθ_M为第M对回波参考对之间的旋转角度偏移；

(3-3)、然后根据被检测物体在检测过程中可能的最大旋转角度作为约束条件，利用二次规划的方法求解被检测物体的旋转向量；

根据被检测物体在检测过程中可能的最大旋转角度作为约束条件，将被检测物体的旋转运动方程组转化为二次规划问题

$\left\{\begin{array}{cc}\min & \frac{1}{2}{R}^T\mathrm{HR}+\mathrm{fR}\\ \mathrm{subject\; to}& \left|\mathrm{\φ}\left(t\right)\right|\≤{\mathrm{\φ}}_{\max }\end{array}\right.$

上式中，H＝2A^TA，f＝-2G^TA，R＝(r₁，r₂，…，r_Q)^T，φ_max为被检测物体在检测过程中可能的最大旋转角度；利用二次规划方法求解上式可得到旋转参数R，从而得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]，N为步骤一中提到的采集到的回波总个数。

本实施例中，设病人的旋转运动幅度最大为15度。因此取基础数据带中属于[15°，165°]之间的K空间线在旋转数据带中寻找参考对。具体方法为，假设基础数据带中的K空间线旋转角度为a，则在旋转数据带中属于[a-15°，a+15°]的所有K空间线计算相关度，相关度最大的即为参考对。计算出参考对的相对偏移角度并列出运动参数方程组，并根据病人的最大旋转角度作为约束条件，利用二次规划的方法求解旋转参数。实验中，旋转曲线φ(t)的阶数Q取3。

上述步骤四的平移运动估计步骤，具体为：

(4-1)、根据HLCC数据一致性原理，被检物体发生运动后存在如下等量关系：

$\underset{-\∞}{\overset{\∞}{\∫}}{(l-d_1\cos \left({\mathrm{\θ}}_c\left(n\right)\right)-d_2\sin \left({\mathrm{\θ}}_c\left(n\right)\right))}^{k}P(n,l)\mathrm{dl}=\underset{r=0}{\overset{k}{\mathrm{\Σ}}}\left(\begin{array}{c}k\\ r\end{array}\right)m_{r,k-r}{\cos }^r\left({\mathrm{\θ}}_c\left(n\right)\right){\sin }^{k-r}\left({\mathrm{\θ}}_c\left(n\right)\right)$

上式中，l是经过步骤(2-1)得到的某一投影数据的位置坐标，d₁是被检测物体在x方向平移运动分量，d₂是受检物体在y方向平移运动分量，θ_c为利用步骤三中得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]进行修正后的旋转角度，θ_c等于K空间旋转向量减去步骤三中得到的被检测物体的旋转向量，k为HLCC等式的阶数，可取k＝0，1，2......，P(n，l)为步骤二中(2-2)校正后的第n行投影数据，m_r，k-r为HLCC定义的被检测物体被检断层图象的几何动量，对于给定的一个受检物体，m_r，k-r为常量；

(4-2)、将物体的平移运动在x方向和y方向的分量用如下关于采样时刻t的Q阶多项式描述：

d_i(t)＝a_i1t+a_i2t²+…+a_iQt^Q，其中参数a_i1到a_iQ为多项式系数；

将上式代入步骤(4-1)中的公式，其中θ_c为利用步骤三中得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]进行修正后的旋转角度，并取k＝1，得到

$\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}{C}_{1q}\left(t\right)a_{1q}+\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}{C}_{2q}\left(t\right)a_{2q}+m_{\mathrm{0,1}}\sin \left({\mathrm{\θ}}_c\left(n\right)\right)+m_{\mathrm{1,0}}\cos \left({\mathrm{\θ}}_c\left(n\right)\right)=C_0\left(t\right)$

其中 $\left\{\begin{array}{c}{C}_0\left(t\right)=\∫l\×P(n,l)\mathrm{dl}\\ C_{\mathrm{lq}}\left(t\right)=\cos \left({\mathrm{\θ}}_c\left(n\right)\right)\×t^q\×P(n,l)\mathrm{dl}\\ C_{2q}\left(t\right)=\sin \left({\mathrm{\θ}}_c\left(n\right)\right)\×t^q\×\∫P(n,l)\mathrm{dl}\end{array}\right.$

对于步骤一中采集的每个回波都可以得到一个方程，共N个方程，利用最小二乘可以求出平移运动分量的多项式参数a_1q，a_2q，从而得到被检测物体的平移运动向量。

根据HLCC数据一致性原理建立平移运动参数方程，共180个方程，其中病人平移运动曲线的阶数Q取3；求取方程组的最小二乘解即可得到病人运动曲线。

上述步骤五运动参数补偿步骤，具体如下：

(5-1)、将K空间旋转向量减去步骤三中得到的被检测物体的旋转向量，得到修正后的K空间旋转向量θ_c；

(5-2)、被检测物体的旋转运动导致了投影数据在角度方向上分布不均，在角度方向进行加权补偿，首先对θ_c进行升序排序得到θ_co，然后计算各个旋转角度对应的加权因子

其中W₁(n)为每个角度对应的加权因子；

(5-3)、对步骤一中采集的基础数据和旋转数据补偿一个与步骤四得到的被检测物体的平移动运动向量相关的相位偏移来消除平移的影响，补偿公式为M′(n，k)＝e^{i×2π×k×d(n)}×M(n，k)，M(n，k)为步骤一中采集到的原始回波数据，M′(n，k)即为相位补偿后的回波数据，k为每个回波上的采样点的编号，d为被检测物体平移运动在K空间径向方向的运动分量，其数学表达式为：

d(n)＝d₁cos(θ_c(n))+d₂sin(θ_c(n))；

(5-4)、将相位补偿后的回波数据M′(n，k)中每个回波做一维傅里叶变换，然后将变换后的数据取模，得到投影数据P₁；

(5-5)、根据投影数据中任意一行回波数据中所有采样点对应的投影数据的总和相同的特征归一化经过步骤(5-4)处理后的投影数据，校正公式为，其中校正后的投影数据为P_c，P₁为经过步骤(5-4)处理后得到的投影数据中任意一个回波中任意一个采样点对应的投影数据，l为经过步骤(5-4)处理后得到的投影数据中任意一行回波的采样点数。

上由于磁场、梯度等的不理想导致采集到的数据与理想数据有偏差，前面的校正步骤很难完全消除运动对数据的影响，因此通过减小存在运动偏差可能性较大的那部分K空间数据对重建图像的贡献来进一步抑制运动伪影；述步骤六中滤波反投影加权重建步骤具体如下：

(6-1)、计算投影向量加权因子，进一步减少未能校正到的运动对图像的影响，C为权重调节因子，C∈[0，∞]；投影向量加权因子W₂(n)表示为

$W_2\left(n\right)=\left\{\begin{array}{cc}\frac{1}{\mathrm{\π}}\mathrm{\θ}\left(n\right)+C& n\≤\frac{1}{2}N\\ 1-\frac{1}{\mathrm{\π}}\mathrm{\θ}\left(n\right)+C& n>\frac{1}{2}N\end{array}\right.$

上式中，θ(n)即为第n行回波数据采集时的旋转角度；

(6-2)、对步骤五中经过步骤(5-5)校正后的投影数据P_c进行加权

P_w(n，l)＝W₁(n)×W₂(n)×P_c(n，l)

其中P_W(n，l)为经过加权后的投影数据，W₁(n)为第n行回波数据对应的加权因子，

W₂(n)为投影向量加权因子，P_c(n，l)为经过步骤(5-5)处理后的第n行回波的投影数据；

(6-3)、以加权后的投影数据P_w作为重建数据，采用滤波反投影方法重建即可得到校正后的图像。

为了检验本发明实施例中提供的运动伪影的去除效果，本实施例中的在0.4T永磁型磁共振设备中采集水模数据进行试验，扫描过程中不断移动水模，直接重建结果参见图7所示，可见图像产生了严重伪影，分辨率小孔完全模糊；用本算法校正后，图像轮廓和细节都得到很大程度恢复，小孔清晰可见，参见图8所示。图5～图8中，Fov为220mm，放射状k空间线数N为180条，每个回波采样点数为256点。

本发明进一步提供了一种MRI系统中K空间采样数据的运动伪影消除装置，该装置包括数据采集装置、数据预处理装置、旋转运动估计装置、平移运动估计装置、运动参数补偿装置、滤波反投影加权重建装置；其中数据采集装置、数据预处理装置、旋转运动估计装置、平移运动估计装置、运动参数补偿装置、滤波反投影加权重建装置的作用和处理步骤与上述运动伪影方法中各步骤相对应一致，这里不再累述。

Claims (2)

1.一种MRI系统中K空间采样数据的运动伪影消除方法，其特征在于：包括以下步骤：

步骤一、数据采集步骤：旋转采集K空间数据，获得具有放射状的K空间数据，具体方法为：

(1-1)、首先从零度开始以Δθ＝2π/N的间隔均匀获取共 

个回波，称为基础数据，记为M_b；

(1-2)、然后从 

开始，仍以Δθ＝2π/N的间隔均匀获取另外 

个回波，称为旋转数据，记为M_r；

基础数据和旋转数据共同组成本步骤采集的K空间原始采样数据；

步骤二、数据预处理步骤：

(2-1)、将步骤一采集到的K空间原始采样数据的每个回波做一维傅里叶变换，然后将变换后的数据取模，即得到投影数据；

(2-2)、接着根据投影数据中任意一行回波数据中所有采样点对应的投影数据的总和相同的特征归一化经过步骤(2-1)处理后的投影数据，去除不同回波信号强度不均的影响，校正公式为： 

其中校正后的投影数据为P，P₀为经过步骤(2-1)处理后得到的投影数据中任意一行回波中任意一个采样点对应的投影数据，l为经过步骤(2-1)处理后得到的投影数据中任意一行回波的采样点数，将基础数据校正后的投影数据记为P_b，将旋转数据校正后的投影数据记为P_r；

(2-3)、对旋转数据校正后的投影数据P_r在角度方向进行T倍上采样得到旋转数据校正后的T倍采样投影数据P_r0；

(2-4)、最后，将基础数据校正后的投影数据为P_b和旋转数据校正后的T倍采样投影数据P_r0做一维傅里叶逆变换并取模，从而得到本步骤预处理后的基础数据F_b和旋转数据F_r；

步骤三、旋转运动估计步骤， 

(3-1)、利用数据相关性在经过步骤二预处理后的基础数据F_b及经过步骤二预处理后的旋转数据F_r中寻找回波参考对，一对回波参考对是指在基础数据和旋转数据中相同旋转角度位置所采集到的一对回波，当扫描过程中被检测物体发生旋转运动时，回波参考对在K空间位置会发生偏移：

对于经过步骤二预处理后的基础数据F_b中的第n行数据F_b(k，n)，k是每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，n为自然数，其取值从1到经过步骤二预处理后的基础数据F_b数据的总行数；通过如下公式计算与经过步骤二预处理后的旋转数据F_r中第m行数据F_r(k，m)的相关性C(n，m)，k为每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，m为自然数，其取值从1到经过步骤二预处理后的旋转数据F_r数据的总行数：

其中C(n，m)表示两个回波之间的相关性，k为每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S， 

分别为F_b(k，n)和F_r(k，m)的均值，σ_n、σ_m分别为F_b(k，n)和F_r(k，m)的方差，S为单个回波最大的采样点数；

经过步骤二预处理后的旋转数据F_r中第m行数据F_r(k，m)如果与经过步骤二预处理后的基础数据F_b中第n行数据F_b(k，n)具有最大相关性，即C(n，m)最大，则它们构成一对回波参考对；

(3-2)、计算经过步骤二预处理后的基础数据F_b及经过步骤二预处理后的旋转数据F_r中回波参考对之间的旋转角度偏移Δθ，然后根据这个旋转角度偏移Δθ，列出被检测物体的旋转运动方程组：

设回波参考对F_b(k，n)和F_r(k，m)的采样时刻分别为t_n和t_m，用φ(t)描述被检测物体的旋转运动方程，可得到如下等式：

φ(t_m)-φ(t_n)＝Δθ

被检测物体的旋转运动方程φ(t)用采样时刻t的Q阶多项式表示为： 

φ(t)＝r₁t+r₂t²+…+r_Qt^Q其中，r₁到r_Q为关于t的Q阶多项式的系数；

这样，就可得到包含Q个旋转运动参数的方程：

r₁(t_m-t_n)+r₂(t_m ²-t_n ²)+…+r_Q(t_m ^Q-t_n ^Q)＝Δθ

设可以在经过步骤二预处理后的基础数据F_b及经过步骤二预处理后的旋转数据F_r中总共能找到M对回波参考对，且 

因此共有M个方程，将上述方程组写成矩阵形式为：

AR＝G

其中 

R＝(r₁，r₂，…，r_Q)^T，T表示矩阵转置；

G＝(Δθ₁，Δθ₂，…，Δθ_M)^T，T表示矩阵转置，Δθ_M为第M对回波参考对之间的旋转角度偏移；

(3-3)、然后根据被检测物体在检测过程中可能的最大旋转角度作为约束条件，利用二次规划的方法求解被检测物体的旋转向量；

根据被检测物体在检测过程中可能的最大旋转角度作为约束条件，将被检测物体的旋转运动方程组转化为二次规划问题

上式中，H＝2A^TA，f＝-2G^TA，R＝(r₁，r₂，…，r_Q)^T，φ_max为被检测物体在检测过程中可能的最大旋转角度；利用二次规划方法求解上式可得到旋转参数R，从而得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]，N为步骤一中提到的采集到的回波总个数；

步骤四、平移运动估计步骤：

(4-1)、根据HLCC数据一致性原理，被检物体发生运动后存在如下等量关系：

上式中，l是经过步骤(2-1)得到的某一投影数据的位置坐标，d₁是被检测物体在x方向平移运动分量，d₂是受检物体在y方向平移运动分量，θ_c为利用步骤三中得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]进行修正后的旋转角度，θ_c等于K空间旋转向量减去步骤三中得到的被检测物体的旋转向量，k为HLCC等式的阶数，可取k＝0，1，2......，P(n，l)为步骤二中(2-2)校正后的第n行投影数据，m_r，k-r为HLCC定义的被检测物体被检断层图象的几何动量，对于给定的一个受检物体，m_r，k-r为常量；

(4-2)、将物体的平移运动在x方向和y方向的分量用如下关于采样时刻t的Q阶多项式描述：

d_i(t)＝a_i1t+a_i2t²+…+a_iQt^Q，其中参数a_i1到a_iQ为多项式系数；

将上式代入(4-1)中的公式，其中θ_c为利用步骤三中得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]进行修正后的旋转角度，并取k＝1，得到

其中

对于步骤一中采集的每个回波都可以得到一个方程，共N个方程，利用最小二乘可以求出平移运动分量的多项式参数a_1q，a_2q，从而得到被检测物体的平移运动向量；

步骤五、运动参数补偿步骤：

(5-1)、将K空间旋转向量减去步骤三中得到的被检测物体的旋转向量，得到修正后的K空间旋转向量θ_c；

(5-2)、被检测物体的旋转运动导致了投影数据在角度方向上分布不均，在角度方向进行加权补偿，首先对θ_c进行升序排序得到θ_co，然后计算各个旋转角度对应的加权因 子 

其中W₁(n)为每个角度对应的加权因子；

(5-3)、对步骤一中采集的基础数据和旋转数据补偿一个与步骤四得到的被检测物体的平移动运动向量相关的相位偏移来消除平移的影响，补偿公式为M′(n，k)＝e^{i×2π×k×d(n)}×M(n，k)，M(n，k)为步骤一中采集到的原始回波数据，M′(n，k)即为相位补偿后的回波数据，k为每个回波上的采样点的编号，d为被检测物体平移运动在K空间径向方向的运动分量，其数学表达式为：

d(n)＝d₁cos(θ_c(n))+d₂sin(θ_c(n))；

(5-4)、将相位补偿后的回波数据M′(n，k)中每个回波做一维傅里叶变换，然后将变换后的数据取模，得到投影数据P₁；

(5-5)、根据投影数据中任意一行回波数据中所有采样点对应的投影数据的总和相同的特征归一化经过步骤(5-4)处理后的投影数据，校正公式为 ，其中校正后的投影数据为P_c，P₁为经过步骤(5-4)处理后得到的投影数据中任意一个回波中任意一个采样点对应的投影数据，l为经过步骤(5-4)处理后得到的投影数据中任意一行回波的采样点数；

步骤六、滤波反投影加权重建步骤：

(6-1)、计算投影向量加权因子，进一步减少未能校正到的运动对图像的影响，C为权重调节因子，C∈[0，∞]；

上式中，W₂(n)投影向量加权因子，θ(n)即为第n行回波数据采集时的旋转角度；

(6-2)、对步骤五中经过步骤(5-5)校正后的投影数据P_c进行加权

P_w(n，l)＝W₁(n)×W₂(n)×P_c(n，l)

其中P_W(n，l)为经过加权后的投影数据，W₁(n)为第n行回波数据对应的加权因 子， 

W₂(n)为投影向量加权因子，P_c(n，l)为经过步骤(5-5)处理后的第n行回波的投影数据；

(6-3)、以加权后的投影数据P_w作为重建数据，采用滤波反投影方法重建即可得到校正后的图像。

2.一种MRI系统中K空间采样数据的运动伪影消除装置，其特征在于：包括

数据采集装置，用于旋转采集K空间数据，获得具有放射状的K空间数据，具体方法为：

(1-1)、首先从零度开始以Δθ＝2π/N的间隔均匀获取共 

个回波，称为基础数据，记为M_b；

(1-2)、然后从 

开始，仍以Δθ＝2π/N的间隔均匀获取另外 

个回波，称为旋转数据，记为M_r；

基础数据带和旋转数据带共同组成本数据采集装置采集的K空间原始采样数据；

数据预处理装置，用于以下功能的处理：

(2-1)、将数据采集装置采集到的K空间原始采样数据的每个回波做一维傅里叶变换，然后将变换后的数据取模，即得到投影数据；

(2-2)、接着根据投影数据中任意一行回波数据中所有采样点对应的投影数据的总和相同的特征归一化经过步骤(2-1)处理后的投影数据，去除不同回波信号强度不均的影响，校正公式为： 

其中校正后的投影数据为P，P₀为经过步骤(2-1)处理后得到的投影数据中任意一行回波中任意一个采样点对应的投影数据，l为经过步骤(2-1)处理后得到的投影数据中任意一行回波的采样点数，将基础数据校正后的投影数据记为P_b，将旋转数据校正后的投影数据记为P_r；

(2-3)、对旋转数据校正后的投影数据P_r在角度方向进行T倍上采样得到旋转数据校正后的T倍采样投影数据P_r0；

(2-4)、最后，将基础数据校正后的投影数据为P_b和旋转数据校正后的T倍采样投影数据P_r0做一维傅里叶逆变换并取模，从而得到本步骤预处理后的基础数据F_b和旋转 数据F_r；

旋转运动估计装置，用于以下功能的处理：

(3-1)、利用数据相关性在经过数据预处理装置预处理后的基础数据F_b及经过数据预处理装置预处理后的旋转数据F_r中寻找回波参考对，一对回波参考对是指在基础数据和旋转数据中相同旋转角度位置所采集到的一对回波，当扫描过程中被检测物体发生旋转运动时，回波参考对在K空间位置会发生偏移：

对于经过数据预处理装置预处理后的基础数据F_b中的第n行数据F_b(k，n)，k是每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，θ(n)为是第n个回波被采集时的旋转角度，n为自然数，其取值从1到经过数据预处理装置预处理后的基础数据F_b数据的总行数；通过如下公式计算与经过数据预处理装置预处理后的旋转数据F_r中第m行数据F_r(k，m)的相关性C(n，m)，k为每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S，θ(m)为第m个回波被采集时的旋转角度，m为自然数，其取值从1到经过数据预处理装置预处理后的旋转数据F_r数据的总行数：

其中C(n，m)表示两个回波之间的相关性，k为每个回波上的采样点的编号，取值从1到单个回波最大的采样点数S， 

分别为F_b(k，n)和F_r(k，m)的均值，σ_n、σ_m分别为F_b(k，n)和F_r(k，m)的方差，S为单个回波最大的采样点数；

经过数据预处理装置预处理后的旋转数据F_r中第m行数据F_r(k，m)如果与经过数据预处理装置预处理后的基础数据F_b中第n行数据F_b(k，n)具有最大相关性，即C(n，m)最大，则它们构成一对回波参考对；

(3-2)、计算经过数据预处理装置预处理后的基础数据F_b及经过数据预处理装置预处理后的旋转数据F_r中回波参考对之间的旋转角度偏移Δθ，然后根据这个旋转角度偏 移Δθ，列出被检测物体的旋转运动方程组：

设回波参考对F_b(k，n)和F_r(k，m)的采样时刻分别为t_n和t_m，用φ(t)描述被检测物体的旋转运动方程，可得到如下等式：

φ(t_m)-φ(t_n)＝Δθ

被检测物体的旋转运动方程φ(t)用采样时刻t的Q阶多项式表示为：

φ(t)＝r₁t+r₂t²+…+r_Qt^Q其中，r₁到r_Q为关于t的Q阶多项式的系数；

这样，就可得到包含Q个旋转运动参数的方程：

r₁(t_m-t_n)+r₂(t_m ²-t_n ²)+…+r_Q(t_m ^Q-t_n ^Q)＝Δθ

设可以在经过数据预处理装置预处理后的基础数据F_b及经过数据预处理装置预处理后的旋转数据F_r中总共能找到M对回波参考对，且 

因此共有M个方程，将上述方程组写成矩阵形式为：

AR＝G

其中 

R＝(r₁，r₂，…，r_Q)^T，T表示矩阵转置；

G＝(Δθ₁，Δθ₂，…，Δθ_M)^T，T表示矩阵转置，Δθ_M为第M对回波参考对之间的旋转角度偏移；

(3-3)、然后根据被检测物体在检测过程中可能的最大旋转角度作为约束条件，利用二次规划的方法求解被检测物体的旋转向量；

根据被检测物体在检测过程中可能的最大旋转角度作为约束条件，将被检测物体的旋转运动方程组转化为二次规划问题

上式中，H＝2A^TA，f＝-2G^TA，R＝(r₁，r₂，…，r_Q)^T，φ_max为被检测物体在检测过程中可能的最大旋转角度；利用二次规划方法求解上式可得到旋转参数R，从而得到被 检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]，N为数据采集装置中提到的采集到的回波总个数；

平移运动估计装置，用于以下功能的处理：

(4-1)、根据HLCC数据一致性原理，被检物体发生运动后存在如下等量关系：

上式中，l是经过步骤(2-1)得到的某一投影数据的位置坐标；d₁是被检测物体在x方向平移运动分量；d₂是受检物体在y方向平移运动分量；θ_c为利用旋转运动估计装置中得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1，φ_N]进行修正后的旋转角度，θ_c等于K空间旋转向量减去旋转运动估计装置中得到的被检测物体的旋转向量；k为HLCC等式的阶数，可取k＝0，1，2......，P(n，l)为旋转运动估计装置中(2-2)校正后的第n行投影数据；m_r，k-r为HLCC定义的被检测物体被检断层图象的几何动量，对于给定的一个受检物体，m_r，k-r为常量；

(4-2)、将物体的平移运动在x方向和y方向的分量用如下关于采样时刻t的Q阶多项式描述：

d_i(t)＝a_i1t+a_i2t²+…+a_iQt^Q，其中参数a_i1到a_iQ为多项式系数；

将上式代入(4-1)中的公式，其中θ_c为利用旋转运动估计装置中得到被检测物体的旋转向量[φ₁，φ₂，φ₃......φ_N-1 ，φ_N]进行修正后的旋转角度，θ_c等于K空间旋转向量减去旋转运动估计装置中得到的被检测物体的旋转向量，并取k＝1，得到

其中

对于数据采集装置中采集的每个回波都可以得到一个方程，共N个方程，利用最小二乘可以求出平移运动分量的多项式参数a_1q，a_2q，从而得到被检测物体的平移运动向量；

运动参数补偿装置，用于以下功能的处理：

(5-1)、将K空间旋转向量减去旋转运动估计装置中得到的被检测物体的旋转向量，得到修正后的K空间旋转向量θ_c；

(5-2)、被检测物体的旋转运动导致了投影数据在角度方向上分布不均，在角度方向进行加权补偿，首先对θ_c进行升序排序得到θ_co，然后计算各个旋转角度对应的加权因子 

其中W₁(n)为每个角度对应的加权因子；

(5-3)、对数据采集装置中采集的基础数据和旋转数据补偿一个与平移运动估计装置得到的被检测物体的平移动运动向量相关的相位偏移来消除平移的影响，补偿公式为M′(n，k)＝e^{i×2π×k×d(n)}×M(n，k)，M(n，k)为数据采集装置中采集到的原始回波数据，M′(n，k)即为相位补偿后的回波数据，k为每个回波上的采样点的编号，d为被检测物体平移运动在K空间径向方向的运动分量，其数学表达式为：

d(n)＝d₁cos(θ_c(n))+d₂sin(θ_c(n))；

(5-4)、将相位补偿后的回波数据M′(n，k)中每个回波做一维傅里叶变换，然后将变换后的数据取模，得到投影数据P₁；

(5-5)、根据投影数据中任意一行回波数据中所有采样点对应的投影数据的总和相同的特征归一化经过步骤(5-4)处理后的投影数据，校正公式为 

，其中校正后的投影数据为P_c，P₁为经过步骤(5-4)处理后得到的投影数据中任意一个回波中任意一个采样点对应的投影数据，l为经过步骤(5-4)处理后得到的投影数据中任意一行回波的采样点数；

滤波反投影加权重建装置，用于以下功能的处理：

(6-1)、计算投影向量加权因子，进一步减少未能校正到的运动对图像的影响，C为权重调节因子，C∈[0，∞]； 

上式中，W₂(n)投影向量加权因子，θ(n)即为第n行回波数据采集时的旋转角度；

(6-2)、对运动参数补偿装置中经过步骤(5-5)校正后的投影数据P_c进行加权

P_w(n，l)＝W₁(n)×W₂(n)×P_c(n，l)

其中P_W(n，l)为经过加权后的投影数据，W₁(n)为第n行回波数据对应的加权因子， 

W₂(n)为投影向量加权因子，P_c(n，l)为经过步骤(5-5)处理后的第n行回波的投影数据；

(6-3)、以加权后的投影数据P_w作为重建数据，采用滤波反投影方法重建即可得到校正后的图像。 

Images