CN102389309B - 基于压缩感知理论的磁共振图像重建的方法 - Google Patents

Abstract

本发明提出一种基于压缩感知理论的磁共振随机采样K空间数据图像重建方法，运用contourlet变换与迭代软阈值法，实现磁共振图像重建。包括首先在磁共振成像扫描仪中按照预设的观测矩阵φ来采集K空间数据，生成测量值保留y。然后从磁共振成像扫描仪的线圈中获得y并传送给计算机。最后构造同样的φ，构造任一种正交变换ψ，根据重构即从y恢复采用基于压缩传感理论的磁共振随机采样K空间数据图像重建方法，节省了扫描时间，实现了快速成像，为医学核磁共振成像检测提供了高质量的可靠图像信息，也为医学成像检测技术的进一步发展和大范围普及应用奠定了坚实的理论和实践基础。

Description

基于压缩感知理论的磁共振图像重建的方法

技术领域

本发明涉及医学成像检测技术领域，特别涉及磁共振成像技术领域，具体是指一种基于压缩感知的磁共振成像重构的方法

背景技术

核磁共振成像(MRI)作为一种无损伤的诊断手段包含了丰富的信息，具有X-CT等成像方法无法比拟的优点.然而，常规的核磁共振成像的时间太长(大约需要数分钟)，导致其成本过高，被验者身体中的生理性运动会产生影像的模糊和对比度的失真，限制了它在临床应用的范围。在保证一定MR图像质量的前提下，加快MR成像速度一直是推动MR技术发展的动力。若单纯从完成检查同一病例的成像时间作比较，MRI的成像时间比CT肯定要长。如从单一层面的成像方面考虑，常规SE序列平均每完成一层图像的扫描时间为15～30s，而常规CT每扫描一层图像所需时间为3s左右。

一般来看，MR成像速度包括2个方面：①原始数据处理速度，②扫描速度。扫描速度才是决定MR成像速度的主要方面。扫描时间由下式确定：扫描时间＝重复时间(TR)×相位编码数×平均次数，所以，影响MR成像速度的主要因素为重复时间，相位编码数及MR信号的平均次数。以上所述的内容是建立在传统自旋回波(SE)序列的基础上进行讨论的。为了加快扫描速度，可以采用缩短TR值和减少相位编码数等方法，但前提是必须保证一定的图像质量。传统的SE序列不必使所有成像区域的组织结构达到充分的弛豫而成像，就使得采用梯度回波的方法对组织成像成为可能。常采用的加快成像速度的方法有(1)梯度回波(GRE)序列(2)快速自旋回波(TSE)序列(3)快速梯度自旋回波(TGSE)序列，而MR快速成像技术的代表——平面回波成像(EPI)，其平均每一层面成像所花费的时间仅为0.1s左右。

磁共振成像技术与其它医学成像技术一个最主要的区别是成像方法的多样性和易控性。磁共振信号空间(原始数据空间)称为K空间，即为傅里叶变换空间，K空间采样得到信号经过傅里叶反变换后，即得到核磁共振图像。K空间概念被广泛地应用于分析与保证MR图像质量和提高成像技术等方面。二维傅立叶变换成像扫描时间为Ts＝Ny×TR×NEX，Ny：相位编码采样数，TR：脉冲重复时间，NEX采样平均数；若得到256×256矩阵图像，TR为2s，NEX为1，即用256条傅立叶线填充K空间，需用8.5min才得到一幅图像；为克服磁共振成像扫描时间长的缺陷，可以利用K空间性质及填充技术减少扫描时间及提高成像质量的方法。如部分K空间(Segmented K-Space)技术、快速自旋回波(Fast Spin Echo：FSE)技术、钥孔(Key Hole)成像、非线性K空间采样(圆形采样(Circular Sampling)、螺旋采样(SpiralSampling)等。

在传统采样过程中，为了避免信号失真，采样频率不得低于信号最高频率的2倍，。然而对于数字图像、视频的获取，依照香农(Shannon)定理会导致海量采样数据，大大增加了存储和传输的代价。若不满足香农定理，重建图像会产生各种伪影。

近年来，一种新兴的压缩感知理论为数据采集技术带来了新的突破，得到了研究人员的广泛关注。压缩感知理论首先由Candes、Romberg、Tao和Donoho等人在2004年提出，Candes证明了只要信号在某一个正交空间具有稀疏性，就能以较低的频率(M＜＜N)采样信号，而且可以以高概率重构该信号。压缩感知采用非自适应线性投影来保持信号的原始结构，能通过数值最优化问题准确重构原始信号，如图1所示基于压缩感知信号数据处理过程。压缩感知以远低于奈奎斯特频率进行采样，在压缩成像系统、模拟/信息转换、生物传感等领域有着广阔的应用前景。

发明内容

本发明的目的是为了克服磁共振成像上述现有技术中的缺点，提供一种能够减少磁共振成像的扫描时间，加快成像速度，从而降低成本，减少被验者身体中的生理性运动会产生影像的模糊和对比度的失真，提高磁共振图像质量，扩大了它在临床应用的范围。基于压缩感知磁共振成像过程如图2所示。

为了实现上述目的，本发明基于压缩感知理论，我们使用Contourlet变换和l₁最优化来重构图像从很少的测量值。磁共振图像重建方法如下：

该基于压缩感知理论的磁共振图像重建方法，其特征在于，所述的方法包括以下步骤：

(1)在磁共振成像扫描仪中按照预设的观测矩阵φ来采集，得到部分K空间数据。构造φ，生成测量值保留y。观测矩阵采用横向正态随机分布来选择K空间轨迹编码线；

(2)从磁共振成像扫描仪的线圈中获得y并传送给计算机。

(3)构造同样的φ，选择一种稀疏变换ψ，根据

重构

即从y恢复该基于压缩感知理论的磁共振图像重建方法，所述的系统的数学模型为：通过一个线性变换

φ是一个m×n维的矩阵，且m＜n；使得图像信号

被编码成更小的向量y。显然，y比

包含更少信息，因此它是图像信号

的压缩。由于

这样可得

也是信号

的压缩，在许多情况下，线性变换φ不是通过计算出来的，φ可以随机性(Randomness)选择，它可以是满足高斯分布的白噪声矩阵，或贝努里分布的±1矩阵(也称作Noiselet)等等。需要注意的是在编码过程

是未知的，φ可以独立于

进行选择。

该基于压缩感知理论的磁共振图像重建方法，所述的图像重建数学模型为：这时如果测量矩阵A满足约束等距条件，可以通过求解一个类似式的最优l₁范数问题(1)来重构稀疏信号

min‖x‖₁s.t‖Ax-y‖₂≤σ                    (1)

建立目标函数：

$f\left(x\right)={\left|\right|\mathrm{Ax}-y\left|\right|}_2^2+\mathrm{\λ}{\left|\right|{\mathrm{\ψ}}^{T}x\left|\right|}_1---\left(2\right)$

获得最小解

$\stackrel{\‾}{x}=\arg \min {\left|\right|\mathrm{Ax}-y\left|\right|}_2^2+\mathrm{\λ}{\left|\right|{\mathrm{\ψ}}^{T}x\left|\right|}_1---\left(3\right)$

其中参数λ是正则化系数。

该基于压缩感知理论的磁共振图像重建方法，所述的目标函数求解的数学模型为：通过迭代阈值法

$x^{i+1}=\mathrm{soft}(x^i+\frac{1}{\mathrm{\α}}{\mathrm{\ψ}}^T{\mathrm{\φ}}^T(y-\mathrm{\φ\ψ}{x}^i),\mathrm{\λ}/2\mathrm{\α})---\left(4\right)$

其中上标i是迭代指数，soft(x，T)是软阈值函数，阈值为T，

soft(x，T)：＝sign(x)max(0，|x|-T)            (5)

为了得到式(2)最小值要求式(4)收敛，我们可以设置α≥max eig(ψ^Tφ^Tφψ)，

该基于压缩感知理论的磁共振图像重建方法，所述的进行解码图像重建包括以下步骤：具体实现步骤如下：

1)输入：原始MRI图像x，Fourier域随机欠采样矩阵φ，最大迭代次数M，重构精度ε，初始权重α1、λ1及步长s。

2)随机测量：对原始图像x进行随机抽样，得到测量值y＝φx；

3)初始化：初始化重构图像x⁰＝φ^Ty，迭代次数i，且i＝0；

4)进行软阈值法迭代，利用公式

5)迭代次数i＝i+1；

6)更新参数：α1＝s×α1；

7)判断迭代停止条件；若i＞M或者

结束迭代，执行步骤8。否则，令i＝i+1返回步骤4继续迭代。

8)获得最优解：

输出重构图像

采用了该发明的基于压缩感知的磁共振数据图像重建方法，由于首先从实际磁共振设备中采集到相位的部分K数据，即找到它的M个线性测量，然后根据这M个测量和Φ，求解最优化问题，得到变换域的数据，然后反变换，便可以得到时域的图像。采用该种压缩感知模型进行磁共振图像重建方法，从而保证图像信噪比和精确度条件下，节省了扫描时间，实现了快速成像；而且相比较现有技术中的部分K空间数据图像重建方法，能够有效降低图像误差，精确显示原图像，为医学磁共振成像检测提供了高质量的可靠图像信息；同时，本发明的方法高效实用，工作性能稳定可靠，适用范围较为广泛，给人们的工作和生活带来了很大的便利，并且也为医学成像检测技术的进一步发展和大范围普及应用奠定了坚实的理论和实践基础。

该发明同时采用了Contourlet变换，由于Contourlet变换稀疏性比小波较好，在一些基于小波的压缩感知重构，重构可能出现高频振铃现象，采用Contourlet变换可以提高信噪比和视觉效果，保护图像细节。我们相信本文的压缩感知方法可以应用到别的成像领域。

附图说明

图1.压缩感知理论框架

图2 MR压缩感知成像系统

图3.Contourlet变换

图4.Contourlet变换系数图

图5.仿真K空间轨迹图

图6.各种K空间填充轨迹(直线型、S型、射线型、螺旋型)

图7.a图是满足采样定理的重构图像，b图是只采样K空间中心区域的重构图像C图是不满足采样定理重构图像，d基于压缩感知理论重建图像

图8.使用小波和Contourlet变换恢复图像结果。(a)脑部图像；(b)胰胆管图像；恢复图像分别利用5％、30％和50％最大变换系数。

具体实施方式

为了能够更清楚地理解本发明的技术内容，特举以下实施例详细说明。

本发明是从实际磁共振设备中采集部分K数据，其数据量是完整K空间数据的1/10至1/5，因而称为压缩感知图像重建方法。

在阐述本发明的整体工作过程及工作原理之前，为了更加明确其技术含义，首先需要介绍Contourlet变换。

Contourlet变换是一种图像的多尺度几何分析工具，能有效地表示轮廓和纹理丰富的图像，非线性逼近能力很强。图3所示为Contourlet变换的结构图。

Contourlet变换将多尺度分析和多方向分析分开进行，总共分两步进行，首先用拉普拉斯金字塔LP(Laplacian pyramid)变换进行多尺度分析捕获点奇异性，接着使用方向性滤波器组DFB(directional filter bank)将高频子带中分布在同一方向上的奇异点连接成线状结构，合成为一个系数。其低频子带则重复上述过程，直到达到预定的分解级数为止。LP和DFB二者的结合就构成了Contourlet变换，称为“塔式方向滤波器组”(PDFB，pyramidal directional filterbank)。图7所示为一测试图像经过Contourlet变换后的系数图，Contourlet变换中LP采用‘9-7’滤波器，DFB则采用‘pkva’方向性滤波器，LP的分解级数为5级。从图4中可以看出，由于LP的冗余性，Contourlet变换系数具有4/3的冗余度。

只有Contourlet与图像的轮廓与边界的位置与方向匹配，从而产生有意义的系数。从系数图解释了变换系数比图像本身更稀疏。当然离散余弦变换(DCT)与小波变换也可以进行稀疏变换，但当重构图像只保留5-10％的系数时候Contourlet变换相比更有些优势。

图像可以压缩，要求图像本身是稀疏的或变换后是稀疏的。首先用向量

表示一图像，寻找稀疏变换ψ，例如小波变换，那么我们可以用公式

来稀疏表示

然后保存向量

中非零元素的数值和它的对应位置。当要恢复

时，我们可以简单的使用同一个稀疏变换ψ，通过公式

来重构，这里

只有很少的K个非零元素

K＜＜N这样传统的图像压缩过程分为两步：第一步编码，即构造ψ，做正变换

保留

中最重要的K个分量，和其对应的位置。第二步解码，即把K个分量放回到对应的位置，其它位置填0，构造ψ，反变换

合理地选择稀疏基ψ，使得信号的稀疏系数个数尽可能少，不仅有利于提高采集信号的速度，而且有利于减少存储、传输信号所占用的资源。本发明用的稀疏基是Contourlet变换。

在图像压缩感知系统中，随机测量抽样并不需要任何的先验知识，但是在图像重建时需要利用图像表示的稀疏性，稀疏性越强，K越小，达到相同重构效果需要的随机投影观测值也越少。

压缩感知的处理过程可以分为两步，首先第一步编码，通过一个线性变换φ，φ一个m×n维的矩阵，m＜n；信号

被编码成更小的向量

显然，y比

包含更少信息，因此它是信号

的压缩。由于

因此

也是信号

的压缩，在许多情况下，线性变换φ不是通过计算出来的，φ可以随机性(Randomness)选择，它可以是满足高斯分布的白噪声矩阵，或贝努里分布的±1矩阵(也称作Noiselet)等等；需要注意的是在编码过程

是未知的，φ可以独立于

进行选择。

第二步解码，即从y恢复

这时如果测量矩阵A满足约束等距条件，可以通过求解一个类似式的最优l₀范数问题(1)来重构稀疏信号最小l₀范数问题是一个NP-hard问题，需要穷举x中非零值的所有CKN种排列可能，因而无法求解[27]，即使找到了它，也并不能肯定这个结果就是对的。鉴于此，求解一个更加简单的l₁优化问题会产生同等的解。

min‖x‖₁s.t‖Ax-y‖₂＜σ                       (1)

建立目标函数：

$f\left(x\right)={\left|\right|\mathrm{Ax}-y\left|\right|}_2^2+\mathrm{\λ}{\left|\right|{\mathrm{\ψ}}^{T}x\left|\right|}_1---\left(2\right)$

获得最小解

$\stackrel{\‾}{x}=\arg \min {\left|\right|\mathrm{Ax}-y\left|\right|}_2^2+\mathrm{\λ}{\left|\right|{\mathrm{\ψ}}^{T}x\left|\right|}_1---\left(3\right)$

压缩感知的编解码过程就是编码：构造φ，生成测量

保留y。解码：构造同样的φ，构造任一种正交变换ψ，根据y重构

本发明基于压缩感知理论，从很少的测量值中，使用Contourlet变换和求l₁最优化解来重构图像。

请参阅图2所示，本发明的基于压缩感知理论的磁共振图像重建方法，包括以下步骤：

(1)在磁共振成像扫描仪中按照预设的观测矩阵φ来采集，得到部分K空间数据。构造φ，生成测量值

其中φ一个m×n维的矩阵，m＜n，m满足O(k·logN)压缩测量就能以大概率精确重构原始信号。观测矩阵φ采用横向正态随机分布选择K-空间轨迹编码线，如图5。

(2)从磁共振成像扫描仪的线圈中获得y并传送给计算机。

(3)构造同样的φ，选择一种稀疏变换ψ，根据y重构

即从

恢复

包括以下步骤：

在测量矩阵A满足约束等距条件，可以通过求解一个类似式的最优l₁范数问题(1)来重构稀疏信号

min‖x‖₁s.t‖Ax-y‖₂＜σ                        (1)

建立式(1)的目标函数：

$f\left(x\right)={\left|\right|\mathrm{Ax}-y\left|\right|}_2^2+\mathrm{\λ}{\left|\right|{\mathrm{\ψ}}^{T}x\left|\right|}_1---\left(2\right)$

获得最小解

$\stackrel{\‾}{x}=\arg \min {\left|\right|\mathrm{Ax}-y\left|\right|}_2^2+\mathrm{\λ}{\left|\right|{\mathrm{\ψ}}^{T}x\left|\right|}_1---\left(3\right)$

下一步求解目标函数。所述的目标函数求解的数学模型为通过迭代阈值法求解式(3)，首先建立具体数学模型为：

$x^{i+1}=\mathrm{soft}(x^i+\frac{1}{\mathrm{\α}}{\mathrm{\ψ}}^T{\mathrm{\φ}}^T(y-\mathrm{\φ\ψ}{x}^i),\mathrm{\λ}/2\mathrm{\α})---\left(4\right)$

其中上标i是迭代指数，soft(x，T)是软阈值函数，阈值为T，

soft(x，T)：＝sign(x)max(0，|x|-T)               (5)

为了得到式(2)最小值要求式(4)收敛，我们可以设置α≥max eig(ψ^Tφ^Tφψ)。

根据建立好的求解目标函数的数学模型，进行解码图像重建，具体包括以下步骤：

1)输入：原始MRI图像x，Fourier域随机欠采样矩阵φ，最大迭代次数M，重构精度ε，初始权重α1、λ1及步长s。

2)随机测量：对原始图像x进行随机抽样，得到测量值y＝φx；

3)初始化：初始化重构图像x⁰＝φ^Ty，迭代次数i，且i＝1；

4)进行软阈值法迭代，利用公式

5)迭代次数i＝i+1；

6)更新参数：α1＝s×α1；

7)判断迭代停止条件；若i＞M或者结束迭代，执行步骤8。否则，令i＝i+1返回步骤4继续迭代。

8)获得最优解：输出重构图像

如图6所示，常见K空间填充轨迹主要分为：标准直线型、圆型、螺旋型及辐射型。在传统采样过程中，如前所述无论哪种快速成像方法，在填充K空间的时候，为了避免信号失真，采样频率都要满足奈奎斯特采样定理，即采样频率不得低于信号最高频率的2倍才能正确恢复图像，如图7(a)所示。否则图像就会产生混叠伪影，如图7(b)(c)所示。

压缩感知理论与传统奈奎斯特采样定理不同，它指出，只要信号是可压缩的或在某个变换域是稀疏的，那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上，然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原图像，如图7(d)所示，可以证明这样的投影包含了重构信号的足够信息。在该理论框架下，采样速率不决定于信号的带宽，而决定于信息在信号中的结构和内容。

采用了上述的基于压缩感知理论的磁共振压缩感知数据图像重建方法，由于首先从实际磁共振设备中采集到相位的部分K数据，然后从部分数据信息的进行模型参数估计，最后根据模型参数估计的结果由压缩感知模型进行磁共振图像重建，采用该种压缩感知模型进行磁共振图像重建方法，在保证图像信噪比分辨率和精确度条件下，节省扫描时间，实现快速成像，为医学核磁共振成像检测提供了高质量的可靠图象信息；本发明的方法高效实用，工作性能稳定可靠、适用范围较为广泛，给人们的工作和生活带来了很大的便利，并且也为医学成像检测技术的进一步的发展和大范围普及应用奠定了坚实的理论和实践基础。

本发明对重构图像非常有效，如图8所示，(a)脑部和(b)痍胆管的核磁图像，尽管输入图像的复杂程度不一样，但重构图像的质量是一样的。特别是在采样率为50％时候，对恢复一个满意可靠的图像总是有效的。30％采样率的时候也没出现明显的伪影，可惜在5％采样率的时候恢复图像有明显的伪影。这些结果能够进一步改善，如果在我们程序中允许迭代次数更多的话。

也可以看到在底低稀疏率情况，小波变换的RMS值小于或等于Contourlet变换。然而在高稀疏率情况，Contourlet变换要比小波变换更有优势。因此，通过Contourlet变换采集更少数据，一个更好的图像重建是可能的。

尽管Contourlet变换比小波变换产生更多原始图像的变换系数，、但Contourlet变换子带的熵比小波变换更小；除此之外，Contourlet是个新的多尺度、多方向的几何分析工具，它可以比小波变换更好的稀疏图像的边界和纹理，因此从这两方面因素对于作为压缩感知的稀疏变换，Contourlet变换总的来说比比小波更好一些。

不同稀疏率下的峰值信噪比，可以看出在较低的稀疏率情况下，用小波变换重构后获得的图像的峰值信噪比好于或者几乎等于Contourlet变换，然而在较高稀疏率的情况下，Contourlet变换的表现要比小波好。因此，通过Contourlet变换利用更少的采样量获得更好的图像重建是可能的。

基于压缩感知理论，我们使用Contourlet变换和l₁最优化来重构图像从很少的测量值。Contourlet变换稀疏性比小波较好，在一些基于小波的压缩感知重构，重构可能出现高频振铃现象，采用Contourlet变换可以提高信噪比和视觉效果，保护图像细节。我们相信本文的压缩感知方法可以应用到别的成像领域。

Claims (3)

1.一种基于压缩感知理论的磁共振图像重建方法，其特征在于，所述的方法包括以下步骤：

(1)在磁共振成像扫描仪中按照预设的观测矩阵φ来采集，得到随机分布的K空间数据，即构造φ，采集图像信号得到测量值向量

其中

为原始图像信号，然后保留y；观测矩阵φ采用满足高斯分布的白噪声矩阵，它是一个m×n维的矩阵，且m＜n；利用观测矩阵φ进行线性变换，使得信号

被编码成更小的测量值向量

显然，y比

包含更少信息，因此y是信号

的压缩；由于

因此

也是信号

的压缩，其中A＝φψ，即观测矩阵φ与稀疏变换ψ的乘积；

(2)从磁共振成像扫描仪的线圈中获得y并传送给计算机；

(3)构造同样的φ，选取一种稀疏变换ψ，根据

重构

是原始图像信号

在稀疏变换ψ上的等价表示，在矩阵A满足约束等距条件下，则通过求解一个类似式的最优l₁范数公式(1)来重构信号

重建误差为σ，

min||x||₁s.t||Ax-y‖₂≤σ         (1)

建立目标函数f(x)：

$f\left(x\right)={\left|\right|\mathrm{Ax}-y\left|\right|}_2^2+\mathrm{\λ}{\left|\right|{\mathrm{\ψ}}^{T}x\left|\right|}_1---\left(2\right)$

获得最优解

$\stackrel{\‾}{x}=\arg \min {\left|\right|\mathrm{Ax}-y\left|\right|}_2^2+\mathrm{\λ}{\left|\right|{\mathrm{\ψ}}^{T}x\left|\right|}_1---\left(3\right)$

公式(2)、(3)中ψ^T是ψ的转置矩阵，参数λ是正则化系数，由此即从y恢复

2.根据权利要求1所述的基于压缩感知理论的磁共振图像重建方法，其特征在于，所述的目标函数求解的数学模型为：通过迭代阈值法

$x^{i+1}=\mathrm{soft}(x^i+\frac{1}{\mathrm{\α}}{\mathrm{\ψ}}^T{\mathrm{\φ}}^T(y-\mathrm{\φ\ψ}{x}^i),\mathrm{\λ}/2\mathrm{\α})---\left(4\right)$

其中上标i是迭代次数，soft(x，T)是软阈值函数，其中T(T＝λ/2α)为迭代阈值，参数λ是正则化系数，sign(x)是取符号函数，

soft(x，T)：＝sign(x)max(0，|x|-T)      (5)

为了得到式(2)最小值要求式(4)收敛，式(4)中ψ^T是ψ的转置矩阵，φ^T是φ的转置矩阵，设置λ＞0，α≥maxeig(ψ^Tφ^Tφψ)，maxeig()是求最大特征值函数。

3.根据权利要求1或权利要求2所述的基于压缩感知理论的磁共振图像重建方法，其特征在于，所述的图像重建包括以下步骤：

1)输入：原始MRI图像

随机采样矩阵φ，最大迭代次数M，重构精度ε，初始参数α₁、λ₁及步长s；

2)随机测量：对原始MRI图像

进行随机采样，得到测量值

3)初始化：初始化重构图像x⁰＝ψ^Tφ^Ty，迭代次数i＝0；

4)进行软阈值法迭代，利用公式 $x^{i+1}=\mathrm{soft}(x^i+\frac{1}{\mathrm{\α}}{\mathrm{\ψ}}^T{\mathrm{\φ}}^T(y-\mathrm{\φ\ψ}{x}^i),\mathrm{\λ}/2\mathrm{\α}),$ 式中ψ^T是ψ的转置矩阵，φ^T是φ的转置矩阵；

5)迭代次数i＝i+1；

6)更新参数：α₁＝s×α₁；

7)判断迭代停止条件；若i＞M或者

结束迭代，执行步骤8，否则，令i＝i+1返回步骤4继续迭代；

8)获得最优解：

输出重构图像

Images