CN104793160B - 一种减少欠采样磁共振成像的频率混迭效应的方法 - Google Patents

Abstract

本发明公开了一种减少欠采样磁共振成像的频率混迭效应的方法，包括以下步骤：步骤一，通过核磁共振仪得到K空间，采用随机振荡法对原始的K空间数据进行欠采样，得到欠采样数据；随机振荡法包括变密度振荡的螺旋轨迹和变密度振荡放射轨迹；步骤二，选择稀疏方法构成超完备字典；步骤三，根据超完备字典，利用恢复算法对步骤一的欠采样数据进行恢复重建，求得原始数据，通过傅里叶变换，得到成像图像，并实时成像。

Description

一种减少欠采样磁共振成像的频率混迭效应的方法

技术领域

本发明涉及一种核磁共振K空间欠采样的随机振荡方法，特别是涉及一种基于压缩感知原理的减少欠采样磁共振成像的频率混迭效应的方法。

背景技术

数据采集时间较长是磁共振成像技术的最大缺点，由于成像速度慢，使该项检查的适用范围大为减少，例如不适合运动器官和危重病人的检查；对于躁动或者丧失自制能力的患者，如不使用镇静剂，也是难以成像；儿科的某些应用同样受到限制。缩短成像时间不仅可以提高效率和病人的舒适度、减少时间依赖性伪影，还是实现心血管检查、功能信息获取、实时温度检测与介入手术成像等动态成像的关键。因此缩短成像时间一直以来都是磁共振成像技术发展的重要目标之一。

中国专利申请2014103131974提出了“一种基于预扫描和非均匀采样的薄层快速磁共振成像方法”，虽然该方案能够通过K空间数据进行分析，提取各层K空间数据中大信号的位置信息，对于传统的方法，可以更加快捷的获得更好地效果。但是这就是简单地变密度算法，很多地方还是会存在一定的混叠伪影影响，随机性效果不明显。

中国专利申请201210390285X提出了一种“磁共振成像采样轨迹优化方法”，虽然该方法在低频段和高频段采用不同的采集策略，对采样轨迹进行优化，优化效率高保证了较高的下采样因子和重建图像质量，但是实质上也是一种简单的变密度平行K空间欠采样，还是有一定的改进空间，同时改方法只用于平行法中，不能有更多的扩展。

在欠采样率较低的情况下，使用放射状与螺旋状欠采样后利用重建算法恢复的图像会产生较为明显的干涉条纹，该现象是由于欠采样导致的频率混叠所造成的。干涉条纹会严重影响MRI成像质量，产生伪影，为临床诊断工作带来干扰。干涉条纹的形状与欠采样轨迹相关，放射状轨迹产生放射状干涉条纹，螺旋状轨迹产生螺旋状干涉条纹。干涉条纹随着欠采样率的提升逐渐淡化，直至消失。

综上所述，如何减少数据采集，消除混叠伪影是核磁共振成像需要解决问题之一。

发明内容

发明目的：本发明所要解决的技术问题是针对现有技术的不足，提供一种减少欠采样磁共振成像的频率混迭效应的方法。

为了解决上述技术问题，本发明公开了一种基于压缩感知原理的K空间欠采样的随机振荡方法，包括如下具体步骤：

步骤一，通过核磁共振仪得到K空间，采用随机振荡法对原始的K空间数据进行欠采样，得到欠采样数据；随机振荡法包括变密度振荡的螺旋轨迹和变密度振荡的放射轨迹；

步骤二，选择合适的稀疏方法，6层小波变换稀疏法，构成超完备字典。

步骤三，利用bregman恢复算法对所得数据进行恢复重建，求得原始数据，通过傅里叶变换，得到成像图像，实时成像。

步骤一中，K空间是寻常空间在傅利叶转换下的对偶空间，主要应用在磁振造影的成像分析，其他如磁振造影中的射频波形设计，以及量子计算中的初始态准备亦用到K空间的概念。K和出现在波动数学中的波数相应，可说都是“频率空间频率”的概念。K-空间是一个抽象空间(三维空间)或平面(二维空间)，MR成像数据根据不同的空间频率编排在特定的K-空间位置，最后被变换成图像。由于K-空间以空间频率为单位(Hz/cm)，空间频率K又是由空间互垂的3个分量Kx、Ky、Kz来描述，Kx、Ky、Kz正好对应一个三维频率空间，所以将该抽象空间称为K-空间。

本发明步骤一可以减少数据的采集，通过构建不同的欠采样轨迹。对于螺旋状欠采样轨迹r_i(θ_i)₁，在极坐标系中的构造方法如式(1-1)，之后可通过式(1-2)转到笛卡尔直角坐标中：

其中n为螺旋状欠采样轨迹的旋转圈数，为螺旋状欠采样轨迹上对应角度的欠采样点到圆心的距离，r_c为半径-角度系数，取值为K空间边长a的一半即r_c小于等于a/2，分别为极坐标生成的螺旋状欠采样轨迹对应的笛卡尔直角坐标。为了满足螺旋形状，与应成正比，但是考虑到图像稀疏变换后大部分信息集中于二维变换域的中心，如果采用关系构造均匀分布的曲线，必然造成中心区域大量数据丢失，因此本研究采用简单变密度螺旋状欠采样轨迹，通过构造欠采样轨迹，使得集中了大量信息的中心区域的欠采样轨迹排布更为紧密，而外围区域的排布较为疏散，简单变密度螺旋状欠采样轨迹的构造方法如式(1-3)：

虽然上述方法生成的螺旋状欠采样轨迹能够欠采样较为丰富的数据，但是在数据恢复后仍然无法避免欠采样所造成的干涉条纹，因此本发明在原有基础的螺旋状欠采样轨迹延径向进行微小随机振荡，优化构造方法如下(1-4)：

其中为优化后的轨迹，dr_i为对应的调整距离，调整公式如下(1-5)：

其中i＝1,2,3,...为对应的序号，，初始i＝1时之后均与的正负有关，k_i为满足高斯分布的一个随机数，为第i₁-1次的调整距离，σ₁为高斯分布的标准差。

对于放射状欠采样轨迹，在极坐标系中的构造方法如式(1-6)，之后可通过式(1-7)转到笛卡尔直角坐标中：

其中r∈[-R,R]，为放射欠状欠采样轨迹上对应角度的欠采样点到圆心的距离，即一个角度对应的到圆心的距离为[-R,R]，分别为极坐标生成的放射状欠采样轨迹对应的笛卡尔直角坐标。

上述方法在数据恢复后无法避免欠采样所造成的干涉条纹，根据螺旋状欠采样的优化方法，本发明在原有基础的放射状欠采样轨迹延角度方向进行微小随机振荡，优化构造方法如下(1-8)：

其中为优化后的角度，为对应的调整角度，调整公式如下(1-9)：

其中j＝1,2,3,...为对应的序号，之后dθ_ij均与的正负有关，k_j为满足高斯分布的一个随机数，为上一次的调整量，σ₂为高斯分布的标准差，为每次调整的角度值，两者同时决定振荡幅度。

考虑到图像稀疏变换后大部分信息集中于二维变换域的中心，如果为恒定常数，则使得欠采样轨迹在中心处的振荡过小，而中心处的欠采样对最终成像的干涉条纹影响最大，因此导致优化效果不明显，本发明对进行缩放处理，提升在中心处的振荡幅度，缩放公式如下(1-10)：

其中Δθ为角度常数，取决于欠采样率，欠采样率越低角度越大，欠采样率跟角度的平方成反比。R为K空间的半径，为当前点到圆心的距离。

本发明与现有技术相比其显著优点在于：一是本发明是基于压缩感知原理对K空间进行欠采样的方法，大大减少了数据的采集，数据小于或等于现有数据次数的35％，即可进行重构，减少了测量过程中引入的系统误差，提高了测量精度，也减少了计算机运行，成像时间；二是本发明随机振荡法能有效的消除欠采样所带来的混叠，伪影影响，在高质量成像的基础下，减少了计算机运算次数。本发明身体各部分核磁共振成像中等。

表1：本发明与现有技术测量方法和装置的对比结果

本发明具有成像速度快、数据采集量小、成像效果好等优点，适用于核磁共振对人体进行成像处理。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1是本发明提出的螺旋轨迹结构对比图。

图2是本发明提出的放射轨迹对比图。

图3是本发明提出的螺旋轨迹结构示意图。

图4是本发明提出放射轨迹结构示意图。

图5为实施例1核磁共振图像在不同σ₁下恢复的图像。

图6为实施例2核磁共振图像在不同σ₂下恢复的图像。

具体实施方式

下面结合附图和实施例对本发明的具体实施方式做进一步的详细说明。

根据本发明提出的一种基于压缩感知原理的K空间欠采样的随机振荡方法，其特征在于包括如下具体步骤：

步骤一，通过核磁共振仪得到K空间，并设计对K空间欠采样的轨迹，采用随机振荡法，变密度振荡的螺旋轨迹和变密度振荡放射轨迹对原始的K空间进行欠采样。

步骤二，选择合适的稀疏方法，6层小波变换稀疏法，构成超完备字典。

步骤三，利用bregman恢复算法对所得数据进行恢复重建，求得原始数据，通过傅里叶变换，得到成像图像，实时成像。

结合图1螺旋轨迹结构示意图和图2放射轨迹结构示意图，本发明提出的一种基于压缩感知原理的K空间欠采样的随机振荡方法，设计的一系列欠采样轨迹函数。它包括核磁共振仪和成像计算机，所述计算机为含有matlab软件。图1中(a)是欠采样率为5％，σ₁＝1情况下的螺旋轨迹。图1中(b)是欠采样率为5％，σ₁＝2情况下的螺旋轨迹。图1中(c)是欠采样率为5％，σ₁＝4情况下的螺旋轨迹。图1中(d)是欠采样率为5％，σ₁＝8情况下的螺旋轨迹。σ₁越大轨迹振荡幅度越大。图2中(a)是欠采样率为5％，σ₂＝0.01情况下的螺旋轨迹。图2中(b)是欠采样率为5％，σ₂＝0.02情况下的螺旋轨迹。图2中(c)是欠采样率为5％，σ₂＝0.05情况下的螺旋轨迹。图2中(d)是欠采样率为5％，σ₂＝0.1情况下的螺旋轨迹。σ₂越大轨迹振荡幅度越大。

以下进一步说明本发明的具体实施例。

实施例1

以本发明应用以螺旋随机振荡为例应用于头部MR成像：

测量目的：对头部进行成像，效果是否优于其他方法，是否存在混叠伪影。

测量装置：核磁共振仪器，计算机

测量样品：江苏省人民医院核磁共振室

测量方法：本发明应用于头部MR成像的具体步骤包括如下：

步骤一，通过核磁共振仪得到K空间，并设计对K空间欠采样的轨迹，采用随机振荡法，变密度振荡的螺旋轨迹对原始的K空间进行欠采样。

对于螺旋状欠采样轨迹r_i(θ_i)₁，在极坐标系中的构造方法如式，之后转到笛卡尔直角坐标中：

其中n为螺旋状欠采样轨迹的旋转圈数，为螺旋状欠采样轨迹上对应角度的欠采样点到圆心的距离，r_c为半径-角度系数，取值为K空间边长的一半即r_c小于等于N/2，分别为极坐标生成的螺旋状欠采样轨迹对应的笛卡尔直角坐标。为了满足螺旋形状，与应成正比，但是考虑到图像稀疏变换后大部分信息集中于二维变换域的中心，如果采用关系构造均匀分布的曲线，必然造成中心区域大量数据丢失，因此本研究采用简单变密度螺旋状欠采样轨迹，通过构造欠采样轨迹，使得集中了大量信息的中心区域的欠采样轨迹排布更为紧密，而外围区域的排布较为疏散，简单变密度螺旋状欠采样轨迹的构造方法如式：

虽然上述方法生成的螺旋状欠采样轨迹能够欠采样较为丰富的数据，但是在数据恢复后仍然无法避免欠采样所造成的干涉条纹，因此本发明在原有基础的螺旋状欠采样轨迹延径向进行微小随机振荡，优化构造方法如下：

其中为优化后的轨迹，dr_i为对应的调整距离，调整公式如下：

其中i＝1,2,3,...为对应的序号，初始i＝1时之后均与的正负有关，k_i为满足高斯分布的一个随机数，为第i₁-1次的调整距离，σ₁为高斯分布的标准差。

其中欠采样率为15％。具体情况结合图1和图3。图3中(a)为欠采样率为5％的基础螺旋轨迹方法，图3中(b)为欠采样率为15％的基础螺旋轨迹方法，图3中(c)为欠采样率为30％的基础螺旋轨迹方法。图3中(d)为欠采样率5％，σ₁＝4的随机振荡螺旋欠采样轨迹。图3中(e)为欠采样率15％，σ₁＝4的随机振荡螺旋欠采样轨迹。图3中(f)为欠采样率30％，σ₁＝4的随机振荡螺旋欠采样轨迹。

步骤二，选择6层小波变换作为稀疏变换，构成超完备字典；

X＝ΨΘ

其中：Ψ＝[ψ₁,ψ₂,l,ψ_N]∈R^N×N为正交基字典矩阵(满足ΨΨ^T＝Ψ^TΨ＝I)，Θ是X在一个稀疏变换域中的投影系数，展开稀疏系数向量Θ＝[θ₁,θ₂,l,θ_N]^T；

结合对信号X的压缩观测，记CS信息算子为A^CS＝ΦΨ，可以得到：

Y＝ΦX＝ΦΨΘ＝A^CSΘ

虽然从Y中恢复Θ也是一个病态问题，但是因为系数Θ是稀疏的，这样未知数个数大大减少，使得信号重构成为可能；

常用的稀疏化方法还包括离散余弦变换、离散傅立叶变换等，由于Θ系数的稀疏度很大程度地影响最终的恢复效果。因此，为了获得最好的恢复效果，选择6层小波变换。

步骤三，在含matlab软件的计算机上，利用恢复算法(具体是bregman迭代算法，可以参见Osher S.,Burger M.,et al.An iterated regularization method for totalvariation-based image restoration[J].Multiscale Model.Simul.,2005(4):460-489.或者Darbon F.,Osher S.Fast discrete optimizations for sparse approximationsand deconvolutions.to appear,2007.)对核磁共振成像的数据进行恢复。

常用的恢复算法除正交匹配追踪法外，还可采用基追踪法、匹配追踪法或共轭梯度法等。

最终计算得出核磁共振图像结果详见图5。图5中(a)是在欠采样率为15％下基础方法恢复效果。图5中(b)是在欠采样率为15％，σ₁＝1情况下的随机振荡螺旋法恢复效果。图5中(c)是欠采样率为15％，σ₁＝2情况下的随机振荡螺旋法恢复效果。图5中(d)是欠采样率为15％，σ₁＝4情况下的随机振荡螺旋法恢复效果。图5中(e)是欠采样率为15％，σ₁＝8情况下的随机振荡螺旋法恢复效果。我们可以看到采用本发明的恢复成像比基础方法更清晰，成功消除了混叠伪影现象。

实施例2，以本发明应用以放射随机振荡为例应用于头部MR成像：

测量目的：对头部进行成像，效果是否优于其他方法，是否存在混叠伪影。

测量装置：核磁共振仪器，计算机

测量样品：江苏省人民医院核磁共振室

测量方法：本发明应用于测量头部MR成像的具体步骤包括如下：

步骤一，通过核磁共振仪得到K空间，并设计对K空间欠采样的轨迹，采用随机振荡法，变密度振荡的放射轨迹对原始的K空间进行欠采样。

生成变密度振荡的放射轨迹

其中r∈[-R,R]，为放射状欠采样轨迹上对应角度的欠采样点到圆心的距离，即一个角度对应的到圆心的距离为[-R,R]，分别为极坐标生成的放射状欠采样轨迹对应的笛卡尔直角坐标。

上述方法在数据恢复后无法避免欠采样所造成的干涉条纹，根据螺旋状欠采样的优化方法，本发明在原有基础的放射状欠采样轨迹延角度方向进行微小随机振荡，优化构造方法如下：

其中为优化后的角度，为对应的调整角度，调整公式如下：

其中j＝1,2,3,...为对应的序号，之后dθ_ij均与的正负有关，k_j为满足高斯分布的一个随机数，为上一次的调整量，σ₂为高斯分布的标准差，为每次调整的角度值，两者同时决定振荡幅度。

考虑到图像稀疏变换后大部分信息集中于二维变换域的中心，如果为恒定常数，则使得欠采样轨迹在中心处的振荡过小，而中心处的欠采样对最终成像的干涉条纹影响最大，因此导致优化效果不明显，本发明对进行缩放处理，提升在中心处的振荡幅度，缩放公式如下：

其中Δθ为角度常数，取决于欠采样率。R为K空间的半径，为当前点到圆心的距离。

其中欠采样率为15％。具体情况结合图2和图4。图4中(a)为欠采样率为5％的基础放射轨迹方法，图4中(b)为欠采样率为15％的基础放射轨迹方法，图4中(c)为欠采样率为30％的基础放射轨迹方法。图4中(d)为欠采样率5％，σ₂＝0.05的随机振荡螺旋欠采样轨迹。图4中(e)为欠采样率15％，σ₂＝0.05的随机振荡螺旋欠采样轨迹。图4中(f)为欠采样率30％，σ₂＝0.05的随机振荡螺旋欠采样轨迹。

步骤二，选择6层小波变换作为稀疏变换，构成超完备字典；

X＝ΨΘ

其中：Ψ＝[ψ₁,ψ₂,l,ψ_N]∈R^N×N为正交基字典矩阵(满足ΨΨ^T＝Ψ^TΨ＝I)，Θ是X在一个稀疏变换域中的投影系数，展开稀疏系数向量Θ＝[θ₁,θ₂,l,θ_N]^T；

结合对信号X的压缩观测，记CS信息算子为A^CS＝ΦΨ，可以得到：

Y＝ΦX＝ΦΨΘ＝A^CSΘ

虽然从Y中恢复Θ也是一个病态问题，但是因为系数Θ是稀疏的，这样未知数个数大大减少，使得信号重构成为可能；

常用的稀疏化方法还包括离散余弦变换、离散傅立叶变换等，由于Θ系数的稀疏度很大程度地影响最终的恢复效果。因此，为了获得最好的恢复效果，选择6层小波变换。

步骤三，在含matlab软件的计算机上，利用恢复算法(具体是bregman迭代算法)对核磁共振成像的数据进行恢复。

常用的恢复算法除正交匹配追踪法外，还可采用基追踪法、匹配追踪法或共轭梯度法等。

最终计算得核磁共振图像结果详见图6。图6中(a)是在欠采样率为15％下基础方法恢复效果。图6中(b)是在欠采样率为15％，σ₂＝0.01情况下的随机振荡放射法恢复效果。图6中(c)是欠采样率为15％，σ₂＝0.02情况下的随机振荡放射法恢复效果。图6中(d)是欠采样率为15％，σ₂＝0.05情况下的随机振荡放射法恢复效果。图6中(e)是欠采样率为15％，σ₂＝0.1情况下的随机振荡放射法恢复效果。我们可以看到采用本发明的恢复成像比基础方法更清晰，成功消除了混叠伪影现象。

本发明提供了一种减少欠采样磁共振成像的频率混迭效应的方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (5)

1.一种减少欠采样磁共振成像的频率混迭效应的方法，其特征在于，包括以下步骤：

步骤一，通过核磁共振仪得到K空间，采用随机振荡法对原始的K空间数据进行欠采样，得到欠采样数据；随机振荡法包括变密度振荡的螺旋轨迹和变密度振荡的放射轨迹；

步骤二，选择稀疏方法构成超完备字典；

步骤三，根据超完备字典，利用恢复算法对步骤一的欠采样数据进行恢复重建，求得原始数据，通过傅里叶变换，得到成像图像，并实时成像；

步骤一中，对于变密度振荡的螺旋轨迹，采用如下构造公式计算轨迹

$adjust\_r_{i_1}\left({\mathrm{\θ}}_{i_1}\right)=r_{i_1}\left({\mathrm{\θ}}_{i_1}\right)+{\mathrm{dr}}_{i_1},$

其中，为简单变密度螺旋状欠采样轨迹，为对应的调整距离；

简单变密度螺旋状欠采样轨迹的计算公式如下：

$r_{i_1}\left({\mathrm{\θ}}_{i_1}\right)=r_c\left|\frac{{\mathrm{\θ}}_{i_1}^3}{2n\mathrm{\π}}\right|,$

其中变密度振荡的螺旋轨迹角度n为螺旋状欠采样轨迹的旋转圈数，r_c为半径-角度系数；

调整距离dr_i的计算方法如下：

${\mathrm{dr}}_{i_1}=\left\{\begin{array}{cc}{\mathrm{dr}}_{i_1-1}+1& i_1\≠1,k_i~N_1\left(-{\mathrm{dr}}_{i_1-1},{\mathrm{\σ}}_1\right)>0\\ {\mathrm{dr}}_{i_1-1}-1& i_1\≠1,k_i~N_1\left(-{\mathrm{dr}}_{i_1-1},{\mathrm{\σ}}_1\right)\≤0\\ 0& i_1=1\end{array}\right.,$

其中i₁＝1,2,3,...为对应的序号，k_i为满足高斯分布的一个随机数，dr_i1-1为第i₁-1次的调整距离，σ₁为高斯分布的标准差。

2.根据权利要求1所述的方法，其特征在于，步骤一中，对于变密度振荡放射轨迹，采用如下构造公式计算轨迹

$r_{i_2}\left({\mathrm{\θ}}_{i_2}\right)=r,$

其中变密度振荡放射轨迹角度r∈[-R,R]，R为K空间的半径；

对进行优化，计算

$adjust\_{\mathrm{\θ}}_{i_{2}j}={\mathrm{\θ}}_{i_2}+{\mathrm{d\θ}}_{i_{2}j},$

其中，为对应的调整角度；

调整角度的计算公式如下：

${\mathrm{d\θ}}_{i_{2}j}=\left\{\begin{array}{cc}{\mathrm{d\θ}}_{i_{2}j-1}+{\mathrm{\Δ\θ}}_{i_{2}j}& j\≠1,k_j~N_2\left(-{\mathrm{d\θ}}_{i_{2}j-1},{\mathrm{\σ}}_2\right)>0\\ {\mathrm{d\θ}}_{i_{2}j-1}-{\mathrm{\Δ\θ}}_{i_{2}j}& j\≠1,k_j~N_2\left(-{\mathrm{d\θ}}_{i_{2}j-1},{\mathrm{\σ}}_2\right)\≤0\\ 0& j=1\end{array}\right.$

其中j＝1,2,3,...为对应的序号，k_j为满足高斯分布的一个随机数，为第j-1次的调整量，σ₂为高斯分布的标准差，为每次调整的角度差值；每次调整的角度差值如下：

${\mathrm{\Δ\θ}}_{i_{2}j}=\mathrm{\Δ}\mathrm{\θ}\·\frac{R}{2\left|r_{i_{2}j}\right|+R},$

其中Δθ为角度常数，R为K空间的半径，为当前点到圆心的距离。

3.根据权利要求2所述的方法，其特征在于，步骤二中，采用6层小波变换稀疏构成超完备字典。

4.根据权利要求1所述的方法，其特征在于，步骤三中，采用bregman算法对核磁共振成像的数据进行恢复。

5.根据权利要求1所述的方法，其特征在于，步骤三中，采用离散傅立叶变换对数据进行变换，得出磁共振图像。