CN107123097A - 一种基于优化的测量矩阵的成像方法 - Google Patents

Abstract

本发明提供一种基于优化的测量矩阵的成像方法，方法包括：步骤01、对目标图像进行傅里叶变换得到K空间数据；步骤02、采用采样矩阵对K空间数据进行采样，得到用于传输的采样信号；步骤03、对目标图像进行稀疏变换，以使重构信号时K‑稀疏的；步骤04、用采样信号获取稀疏系数；步骤05、采用快速迭代阈值算法对稀疏系数进行处理，以获得重建的图像；其中，采样矩阵为通过对随机径向采样矩阵、类圆环采样矩阵进行叠加获取的，适于与稀疏矩阵之间满足不相干特性的变密度径向类圆环矩阵。上述方法中变密度径向类圆环采样矩阵具有高度的随机性与稀疏矩阵之间容易满足不相干特性，进而恢复图像的效果较好。

Description

一种基于优化的测量矩阵的成像方法

技术领域

本发明涉及磁共振成像技术，特别是一种基于优化的测量矩阵的成像方法。

背景技术

在压缩感知理论中，测量矩阵的设计在信号的采集和重建过程中占据着十分重要的地位。测量矩阵的选取是否合理，直接影响压缩感知采样信号的稀疏性，从而影响了重建算法对原始信号的观测次数以及重构图像的质量。

在压缩感知-磁共振成像(CS-MRI)模型中，测量矩阵又称为采样矩阵即为磁共振成像系统中的k空间采样矩阵，采样矩阵能够极大的影响压缩传感的性能以及重建图像的质量。常用的采样矩阵包括螺旋采样矩阵、笛卡尔采样矩阵以及二维变密度随机采样矩阵等。

然而，采用笛卡尔采样矩阵的磁共振成像系统恢复出的图像效果不好，含有伪影。另外，采用单个的螺旋采样矩阵或者二维变密度随机采样矩阵均在磁共振成像系统恢复出的图像效果不理想，不能较好的采集图像的边缘信息且不能保证信息的完整性。为此，提供一种可以使重建算法收敛速度加快，得到的最优解更加精确的磁共振成像中的测量矩阵成为当前需要解决的问题。

发明内容

针对现有技术中的缺陷，本发明提供一种基于优化的测量矩阵的成像方法，该变密度径向类圆环采样矩阵具有高度的随机性与稀疏矩阵之间容易满足不相干特性，进而恢复图像的效果较好。

第一方面，本发明提供一种基于优化的测量矩阵的成像方法，包括：

步骤01、对目标图像进行傅里叶变换得到K空间数据；

步骤02、采用采样矩阵对K空间数据进行采样，得到用于传输的采样信号；

步骤03、对目标图像进行稀疏变换，以使重构信号时K-稀疏的；

步骤04、用采样信号获取稀疏系数；

步骤05、采用快速迭代阈值算法对稀疏系数进行处理，以获得重建的图像；

其中，采样矩阵为通过对随机径向采样矩阵、类圆环采样矩阵进行叠加获取的，适于与稀疏矩阵之间满足不相干特性的变密度径向类圆环矩阵。

可选地：稀疏变换为小波变换、离散余弦变换以及奇异值稀疏变换中的一种或多种。

可选地，所述步骤02之前，所述方法还包括：

依据预先定义的径向采样矩阵和径向采样矩阵的采样密度函数、类圆环采样矩阵及该类圆环采样矩阵的采样密度函数生成变密度径向类圆环采样矩阵。

可选地，所述径向采样矩阵和径向采样矩阵的采样密度函数包括：

给定一个256*256，且每一元素值为0的采样矩阵；

给定一个第一采样密度函数d(r)，

判断采样矩阵中的任意元素与采样矩阵中心点(x,y)之间的关系是否满足预设的第一条件；

若满足，将满足第一条件的采样矩阵中的元素置为1，否则为0，获得径向采样矩阵D；

其中，r为采样位置到采样矩阵中心点的距离，k为直线距离的斜率tan r，R为采样矩阵的点到采样矩阵中心点的最大距离；

第一条件为|y-kx|≤1/2。

可选地，类圆环采样矩阵及该类圆环采样矩阵的采样密度函数包括：

给定一个第二采样密度函数f(r)，f(r)＝1-(k′r)^p/R；k′为调整系数，0＜k′＜1,P为大于0的指数；

判断采样矩阵中的任意元素与采样矩阵中心点(x,y)之间的关系是否满足预设的第二条件；

若满足，将满足第二条件的采样矩阵中的元素置为1，否则为0，获得类圆环采样矩阵R；

其中，r为采样位置到采样矩阵中心点的距离，k′为调整稀疏，0＜k′＜1，p是指数，p>0，R为采样矩阵的电到采样矩阵中心点的最大距离；

第二条件为

可选地，所述步骤03，包括：

对目标图像依次进行小波变换，离散余弦变换以及奇异值稀疏变换。

可选地，所述步骤05包括：

通过快速迭代阈值算法解稀疏重建模型中的稀疏系数，最后重建出图像。

可选地，通过快速迭代阈值算法解稀疏重建模型中的稀疏系数，包括：

根据公式一，获得稀疏系数；

α_k＝p_L(β_k) 公式一

k＝1开始迭代，若||α_k-α_k-1||₂＞ε,,k←k+1，转入执行公式一的步骤；否则，结束；

其中，α为稀疏系数，ε为误差界

L为函数f的梯度lipschitz常数初始值β₁＝α₀∈R^N，

本发明实施例中的采样矩阵能够极大的影响压缩传感的性能以及重建图像的质量，因此对于本发明中改进的采样矩阵，必须满足与稀疏基之间的不相干特性，从而保证重构过程能够实现高精度的恢复。

另外，径向采样矩阵的采样轨迹满足高斯分布，满足高斯分布轨迹的采样矩阵与任何稀疏基之间的相关性都非常低，在磁共振采样的同时，可以选择少部分的稀疏数据来表示整个图像，将高维信号投影到低维空间中，即考虑到硬件设施的限制情况，且能够满足磁共振采样轨迹基本上在相对平滑的直线或者曲线上进行。

本实施例的优化的测量矩阵属于径向矩阵的衍生，且包含径向矩阵的各种特性，故可以使重建算法收敛速度加快，得到的最优解更加精确。

附图说明

图1为本发明的基于优化的测量矩阵的成像方法的流程图；

图2(a)为k空间数据；

图2(b)为现有技术中径向采样轨迹模式的示意图；

图2(c)为现有技术中类圆环采样矩阵的轨迹模式示意图；

图2(d)本发明中改进采样矩阵的轨迹模式示意图；

图3(a)原是真实图像；

图3(b)为径向采样方式下，用迭代阈值算法重建的图像；

图3(c)为类圆环采样方式下，用迭代阈值算法重建的图像；

图3(d)为改进采样方式下，用迭代阈值算法重建的图像。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。

实施例一

如图1所示，本实施例提供一种基于优化的测量矩阵的成像方法，该方法包括：

步骤01：对目标图像进行傅里叶变换得到K空间数据。

步骤02：采用采样矩阵对K空间数据进行采样，得到用于传输的采样信号。

在本实施例中，变密度径向类圆环采样矩阵可为通过对径向采样矩阵、类圆环采样矩阵进行叠加获取的，适于与稀疏矩阵之间满足不相干特性的采样矩阵。

也就是说，依据预先定义的径向采样矩阵和径向采样矩阵的采样密度函数、类圆环采样矩阵及该类圆环采样矩阵的采样密度函数生成变密度径向类圆环采样矩阵。

步骤03：对目标图像进行稀疏变换，以使重构信号时K-稀疏的。

举例来说，对目标图像依次进行小波变换、离散余弦变换及奇异值稀疏变换。

步骤04、用采样信号获取稀疏系数；

步骤05、采用快速迭代阈值算法对稀疏系数进行处理，以获得重建的图像。

针对上述步骤可理解的是：测量得到的欠采样的k空间数据表示为y＝Φx，其中x表示待重建的图像，Φ是采样矩阵，y是采样值。x在变换域下稀疏表示为x＝Ψα，其中Ψ为稀疏基。即y＝Φx＝ΦΨα＝Aα，故提出的稀疏重建模型为：

其中正则化参数λ用于权衡λ||α||₁和||y-Aα||²两项的重要性，α表示稀疏系数。

图像重建：解稀疏重建模型中的稀疏系数α通过快速迭代软阈值算法，最后重建出图像x。

举例来说，步骤02可为：y＝Φx

步骤03可为：x＝Ψα

则步骤04可为：y＝ΦΨα＝Aα，α为稀疏系数，A为传感矩阵，如上说明的采样矩阵和稀疏基的乘积。

具体地，针对步骤05，采用快速迭代阈值算法获取稀疏系数。

例如α_k＝p_L(β_k) 公式一

k＝1开始迭代，若||α_k-α_k-1||₂＞ε,k←k+1，转入执行公式一的步骤；否则，结束；

其中，α为稀疏系数，ε为误差界，公式一种的β只是α的线性组合，快速迭代阈值算法与迭代阈值算法的区别就在于它的更新值不只与它的前一次值有关，而是前两次或多次值的线性组合。

L为函数f的梯度lipschitz常数初始值β₁＝α₀∈R^N，

可理解的是，快速迭代阈值算法(FISTA)的迭代步骤如下：

迭代步1：给定初始值β₁＝α₀∈R^N，(的lipschitz常数)，k＝1，t₁＝1，误差界ε。

迭代步2：根据下面公式计算

α_k:＝p_L(β_k)

迭代步3：若||α_k-α_k-1||₂＞ε,k←k+1，转入迭代步骤2；否则，结束。

上述迭代步骤中：f(α):＝||y-Aα||²,g(α)＝λ||α||₁，F(α):＝f(α)+g(α)，L为函数f的梯度lipschitz常数，F(α)在点β处的二次逼近函数为：

Q_L(α,β)具有唯一的极小值点，记为：

p_L(β):＝argmin[Q_L(α,β):α∈Rⁿ]

进一步地，针对上述的步骤02，可具体说明如下。下述生成变密度径向类圆环采样矩阵中的参数与上述的快速迭代阈值算法中的参数的含义无任何关联。

第一步：给定一个256*256，且每一元素值为0的采样矩阵；

给定一个第一采样密度函数d(r)，

判断采样矩阵中的任意元素与采样矩阵中心点(x,y)之间的关系是否满足预设的第一条件；

若满足，将满足第一条件的采样矩阵中的元素置为1，否则为0，获得径向采样矩阵D；

其中，r为采样位置到采样矩阵中心点的距离，k为直线距离的斜率tan r，R为采样矩阵的电到采样矩阵中心点的最大距离；

第一条件为|y-kx|≤1/2。

第二步：给定一个第二采样密度函数f(r)，f(r)＝1-(k′r)^p/R；k′为调整系数，0＜k′＜1,P为大于0的指数；

判断采样矩阵中的任意元素与采样矩阵中心点(x,y)之间的关系是否满足预设的第二条件；

若满足，将满足第二条件的采样矩阵中的元素置为1，否则为0，获得类圆环采样矩阵R；

其中，r为采样位置到采样矩阵中心点的距离，k′为调整稀疏，0＜k′＜1，p是指数，p>0，R为采样矩阵的电到采样矩阵中心点的最大距离；

第二条件为

第三步：将径向采样矩阵D和类圆环采样矩阵R叠加，生成变密度径向类圆环采样矩阵。

在本实施例中，采样矩阵具有现有技术中随机矩阵的随机性，很容易与稀疏矩阵满足不相干特性，能高概率的恢复图像。

由于，采样矩阵能够极大的影响压缩传感的性能以及重建图像的质量，因此对于本发明实施例中改进的采样矩阵，必须满足与稀疏基之间的不相干特性，从而保证重构过程能够实现高精度的恢复。

另外，径向采样矩阵的采样轨迹满足高斯分布，满足高斯分布轨迹的采样矩阵与任何稀疏基之间的相关性都非常低，在磁共振采样的同时，可以选择少部分的稀疏数据来表示整个图像，将高维信号投影到低维空间中，即考虑到硬件设施的限制情况，且能够满足磁共振采样轨迹基本上在相对平滑的直线或者曲线上进行。

本实施例的优化的测量矩阵属于径向矩阵的衍生，且包含径向矩阵的各种特性，故可以使重建算法收敛速度加快，得到的最优解更加精确。

实施例二

本实施例中可基于CS-MRI的测量矩阵的优化，实现对图像重建的方法，具体实施步骤如下：

第一步：选择合适的稀疏变换基(即稀疏矩阵)对目标图像进行稀疏变换。

例如，可以选择小波变换。离散余弦变换以及奇异值稀疏变换作为稀疏矩阵。

第二步：通过变密度径向类圆环采样矩阵分别为各种稀疏变换基下的采样方法。

为了更直观清晰地表示改进采样矩阵Φ(即变密度径向类圆环采样矩阵)的信息。下面取采样率为30％，对改进采样矩阵进行表示，如图2所示。

从图2可以看出，改进的采样矩阵对K空间中心部分进行高度的采集，而边缘部分相比于中心位置逐渐采集密度降低，且优于传统的径向采样，边缘采集较原有采样方式多，保证信息的完整性。

第三步：最后通过快速迭代阈值算法对采样信号进行恢复重构，以获得重构的图像。

需要说明的是，本实施例中为了对比说明，对径向采样矩阵、类圆环采样矩阵、改进的采样矩阵(即变密度径向类圆环矩阵)分别进行举例说明，并对重构的图像进行说明。

特别地，本实施例中下述公式二计算改进的变密度径向类圆环矩阵与不同的稀疏变换矩阵--小波变换、DCT离散余弦变换以及奇异值分解稀疏变换矩阵之间不同的相干系数值，并与常用的径向采样矩阵和类圆环采样矩阵的相干系数进行对比。

相干度计算公式二：

其中，若Φ,Ψ之间包含相关元素，那么相关度就越大，否则就越小。当μ的值为1时，说明两者不相干，采样矩阵与稀疏矩阵之间不相干性越强，图像恢复的概率就越高，重建精度也越好；反之，相干性越强，图像恢复的概率越低。

如下述表1是采样矩阵与小波稀疏变换基之间相干性的对比，可以看出改进的变密度径向类圆环采样矩阵与小波稀疏变换基之间的相干性在五种不同的采样率下都是最小的，其次是类圆采样环矩阵，说明改进的采样矩阵用于恢复脑磁共振图像时，恢复效果也是最好的，并且采样数据减少了。

表1

表2是采样矩阵与DCT稀疏变换基之间相关性的对比，可以看出改进的变密度径向类圆环采样矩阵与DCT稀疏变换基之间的相干性在五种不同的采样率下都是最小的，其次是类圆采样环矩阵，说明改进的采样矩阵用于恢复脑磁共振图像时，恢复效果也是最好的，并且用较少的欠采样数据去完整恢复原始图像。

表2

表3是采样矩阵与奇异值分解稀疏变换基之间相关性的对比，可以看出改进的变密度径向类圆环采样矩阵与奇异值分解稀疏变换基之间的相干性在五种不同的采样率下都是最小的，其次是类圆环采样矩阵，而径向采样矩阵的相干性最大。

表3

采样矩阵设计的优劣直接关系到图像重建质量的好坏，因此用本发明采样矩阵与径向采样、类圆环采样同时恢复脑共振图像，根据恢复出的图像效果，来评判采样矩阵的优劣。从两方面进行比较，一是客观数据分析指标，峰值信噪比以及均方误差；二是主观视觉恢复效果图。其中峰值信噪比和均方误差定义如下：

表4是不同采样矩阵下图像恢复峰值信噪比，可以看出改进的采样矩阵的峰值信噪比在不同采样率下的值均高于其他两种常用采样矩阵的值。

表4

表5是不同采样矩阵下图像恢复均方误差，可以看出改进的采样矩阵在不同采样率下重建图像的均方误差均低于其他三种采样矩阵下重建的均方误差。

表5

以下采用小波变换作为稀疏变换，用径向采样、类圆环采样、以及本发明提出的改进变密度径向类圆环采样矩阵分别对大脑图片在K空间中进行采样，采样率均为30％，重构算法均选用快速迭代阈值算法(FISTA)。

实验结果如图3所示。从实验结果可以看出在改进的变密度径向类圆环采样矩阵下恢复的图像与真实的图像最为接近，不仅能完整的重构图像，而且还能很好的保留图像的细节信息，达到更好的恢复效果。

最后应说明的是：以上所述的各实施例仅用于说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或全部技术特征进行等同替换；而这些修改或替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (8)

1.一种基于优化的测量矩阵的成像方法，其特征在于，包括：

步骤01、对目标图像进行傅里叶变换得到K空间数据；

步骤02、采用采样矩阵对K空间数据进行采样，得到用于传输的采样信号；

步骤03、对目标图像进行稀疏变换，以使重构信号时K-稀疏的；

步骤04、用采样信号获取稀疏系数；

步骤05、采用快速迭代阈值算法对稀疏系数进行处理，以获得重建的图像；

其中，采样矩阵为通过对随机径向采样矩阵、类圆环采样矩阵进行叠加获取的，适于与稀疏矩阵之间满足不相干特性的变密度径向类圆环矩阵。

2.根据权利要求1所述的方法，其特征在于：稀疏变换为小波变换、离散余弦变换以及奇异值稀疏变换中的一种或多种。

3.根据权利要求1所述的方法，其特征在于，所述步骤02之前，所述方法还包括：

依据预先定义的径向采样矩阵和径向采样矩阵的采样密度函数、类圆环采样矩阵及该类圆环采样矩阵的采样密度函数生成变密度径向类圆环采样矩阵。

4.根据权利要求3所述的方法，其特征在于，所述径向采样矩阵和径向采样矩阵的采样密度函数包括：

在目标图像为256*256的图像时，给定一个256*256，且每一元素值为0的采样矩阵；

给定一个第一采样密度函数d(r)，

判断采样矩阵中的任意元素与采样矩阵中心点(x,y)之间的关系是否满足预设的第一条件；

若满足，将满足第一条件的采样矩阵中的元素置为1，否则为0，获得径向采样矩阵D；

其中，r为采样位置到采样矩阵中心点的距离，k为直线距离的斜率tan r，R为采样矩阵的点到采样矩阵中心点的最大距离；

第一条件为|y-kx|≤1/2。

5.根据权利要求3或4所述的方法，其特征在于，类圆环采样矩阵及该类圆环采样矩阵的采样密度函数包括：

给定一个第二采样密度函数f(r)，f(r)＝1-(k′r)^p/R；k′为调整系数，0＜k′＜1,P为大于0的指数；

判断采样矩阵中的任意元素与采样矩阵中心点(x,y)之间的关系是否满足预设的第二条件；

若满足，将满足第二条件的采样矩阵中的元素置为1，否则为0，获得类圆环采样矩阵R；

其中，r为采样位置到采样矩阵中心点的距离，k′为调整稀疏，0＜k′＜1，p是指数，p>0，R为采样矩阵的电到采样矩阵中心点的最大距离；

第二条件为

6.根据权利要求2所述的方法，其特征在于，所述步骤03，包括：

对目标图像依次进行小波变换，离散余弦变换以及奇异值稀疏变换。

7.根据权利要求1所述的方法，其特征在于，所述步骤05包括：

通过快速迭代阈值算法解稀疏重建模型中的稀疏系数，最后重建出图像。

8.根据权利要求7所述的方法，其特征在于，通过快速迭代阈值算法解稀疏重建模型中的稀疏系数，包括：

根据公式一，获得稀疏系数；

α_k＝p_L(β_k) 公式一；

k＝1开始迭代，若||α_k-α_k-1||₂＞ε,,k←k+1，转入执行公式一的步骤；否则，结束；

其中，α为稀疏系数，ε为误差界；

L为函数f的梯度lipschitz常数初始值β₁＝α₀∈R^N，