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

杨晓兰，朱永贵，丛佳，刘平

(中国传媒大学理学院，北京 100024)

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

杨晓兰，朱永贵，丛佳，刘平

(中国传媒大学理学院，北京 100024)

由于一些器官的边界信息在大多数核磁共振图像中都是稀疏的，所以利用压缩感知从数量非常有限的观测数据集合中重构同样的核磁共振图像并且大大减少核磁共振图像的扫描磨损成为可能。然而，为了能够做到这一点，我们必须要解决定义在大量数据集合上的非光滑函数的最小化这一困难问题。为了解决这一问题，我们给出了一个有效算法，它克服以往求解l1问题的计算复杂性，提出β范数近似逼近l1范数的思想，由于β范数具有光滑性，可采用Bregman迭代正则化方法进行求解。数值实验证明，核磁共振图像可以从全部数据的40%抽样中几乎精确重构。

压缩感知；核磁共振成像；重构算法；Bregman方法

1 引言

对于压缩感知成像，其目的是尽可能减少探测单元数量，更精确的重构出图像。因此在感知信息获取投影方面，该方面的研究需要沿着以下几个方面进行：(a)投影方法的选择，投影方法需要满足以下三个条件：1)在信息采集方面，能够使需要进行的测量尽量少；2)便于投影方法的硬件实现和优化算法的迭代运算；3)能够适用于大多数的稀疏情况。(b)寻找稳定性更高，重建精度更好的重构算法，使得获得的图像的精度满足实际应用需要。

如今，成像在医疗诊断方面起到了着非常重要的作用。事实表明减少扫描时间非常重要，即减少核磁共振成像所需要的时间，那就意味着在保证图像质量的情况下，尽可能少地从频域中收集观测数据。然而，这看上去直接违背了长期以来建立的奈奎斯特标准：即所获得的观测数据的数量必须至少与恢复图像所需要的信息数量相匹配。这就意味着完美的图像重构将是不可能的，但压缩感知允许我们在没有相关噪声的情况下重构图像。

我们给出了一个有效算法，它克服以往求解l1问题的计算复杂性，提出β范数近似逼近范数的思想，由于β范数具有光滑性，可采用Bregman迭代正则化方法进行求解。数值实验证明，核磁共振图像可以从全部数据的40%抽样中几乎精确重构。

2 基于压缩感知的核磁共振图像重建的Bregman方法

通过压缩感知我们得到b=Ru，其中传感矩阵R是大小为m×n的矩阵。现在若想通过已经得到的b恢复出原信号u，则通过直接求解上述方程组是无法得到，因为上述方程组的未知数个数n大于方程个数m。如果一个方程组中未知数个数多于方程个数，则一般有无穷解，通常的求解方法为[1]：

(1)寻求更多的方程；

(2)选自由变量，给一组线性无关的值后，进行求解。

如果未知u的其它信息，则上述逆问题很难求得唯一解。但若信号u是k-稀疏的，设m≥k·log(n)，则上述方程组可以求解，因为此时未知数个数事实上只有k个。尽管如此，由于我们并不知道u中哪些分量为零，所以不能直接求解，而是通过下式求最稀疏的向量，从而获得u：

(2.1)

其中‖·‖0表示u中非零元素的个数。然而，尽管从理论上来讲，这种方法可以实现，但是在实际中不可行，因为这是一个组合问题(NP难问题)，计算量非常大。因此我们需要寻求合适的范数来近似求解上述问题。采用l1范数最小化求得的解，为稀疏向量，非常接近最小化l0范数所得的真实解。因此，我们可以选择l1范数最小化来近似求解上述优化问题，即

(2.2)

刚才讨论的信号u本身是稀疏的，接下来讨论信号u本身并不稀疏，但是可以压缩的情况。我们需要将信号变换到变换域考虑。设Ф是压缩基(如小波基)或紧框架，满足规范正交，即Ф*Ф=I，作变换Фu=x，显然x是稀疏的。于是在这种情况下的求解方法为：

(2.3)

其中A=RФ-1，并且要求m≥μ2(R，Ф-1)·k·logn。

以上考虑的都是等式约束，然而实际中，测量过程可能会引入噪声，这时约束条件式中的Ax=b必须被放松，从而引出问题

(2.4)

或者问题

(2.5)

其中σ与μ是参数。从最优化理论上来说，问题(2.4)和(2.5)在某种意义上是等价的，解决了其中一个问题就能决定另一个问题的参数，使得它俩能够给出相同的解[2]。

我们将提出β范数近似逼近l1范数的思想，他克服了以往求解l1问题的计算复杂性的问题，由于β范数具有光滑性，使得模型可采用Bregman迭代正则化方法进行求解。

(2.6)

令J(u)=μ|u|β，基于凸函数J(u)的u，v之间的Bregman距离定义为：

p∈∂J(v)是J在v点的次微分中的某一个次梯度。

可通过求解：

(2.8)

k=0，1，2，3…

u0=0，p0=0(k=0为原始问题(2.6))来替代对(2.6)问题的求解，其证明参见文献[3]。

因为J(u)不是处处可微，故J(u)的次微分可能包含不止一个元素，从(2.8)中uk+1的最优性，对(2.8)关于u求偏导再以uk+1替换可以得到：

0∈∂J(uk+1)-pk+A*(Auk+1-b)

因此令：

pk+1=pk+A*b-A*Auk+1。

(2.8)与(2.6)的区别在于正则化的应用，(2.6)通过直接最小化u的β-范数将其正则化，而(2.8)是通过最小化u的基于β-范数的Bregman距离将其正则化。

文献[3]中的数值结果表明，当u足够大，这一简单迭代过程明显改进了原始模型(2.6)的去噪能力。

除了对于迭代过程k=0，对于所有的k，(2.8)可以演绎为问题(2.6)只要满足：

bk+1：=b+(bk-Auk)
b0=u0=0

其初始条件为：b0=u0=0，因此迭代过程(2.8)可以等价于

(2.9)

(2.9)可以通过任意一种已存在的求解(2.6)的方法求解。

方法一：

u0←0，p0←0

(2.81)

Fork=0，1，…do

(2.82)

pk+1←pk-A*(Auk+1-b)

(2.83)

方法二：

b0←0，u0←0

(2.91)

Fork=0，1，…，N，do

bk+1=b+(bk-Auk)

(2.92)

(2.93)

在方法一中，给出uk，pk，uk+1满足一阶最优条件：

0∈∂J(uk+1)-pk+▽H(uk+1)
=∂J(uk+1)-pk+A*(Auk+1-b)

因此：

pk+1=pk-A*(Auk+1-b)∈∂J(uk+1)

定理二表示Bregman迭代过程的方法一(2.81)-(2.83)，方法二(2.91)-(2.93)在(2.82)及(2.93)有相同的目标函数(达到同一个常量)的情况下是等价的。

由(2.91)可以得出：b1=b因此当k=0时，(2.82)(2.93)求解相同的最优化问题：

A*(b-Au)
=wF*P*(b-PFw-1u)

故：

根据(2.83)p0=0，b=b1可以得到：

试证明：

(i) (2.82)(2.93)最优化问题的第k次迭代是等价的

关于(i)的证明：由假设可知：

其中C1，C2，C3为常数，因此(2.82)(2.93)有相同的目标函数；

由(i)及[4]的结论可以得到；

关于(iii)的证明：根据假设以及(2.83)(2.92)，及(ii)可以得到：

pk+1=pk-A*(Auk+1-b)

=pk-wF*P*(PFw-1uk+1-b)

注：如果J不是严格凸函数，方法一、二也许会求的不止一个解，上面的证明过程告诉我们，即使方法一及方法二的某一步迭代产生了不同的中间值，他们之后还是会产生相同的结果。

(2.93)的每次迭代都是(2.6)的一个实例，都可以通过FPC(Fixed-point continuation Methed)[5]来求解，尽管对任何严格正的μ都能得到收敛的结果，我们还是要选取能使得(2.6)可以用FPC有效求解的μ，并且使得Bregman迭代的总时间是最优的。

下面我们将求解模型：

(2.9)

bk+1=b+(bk-Auk)

此时u是近似稀疏信号，且A=PF指的是部分傅利叶变换，如果令x代表原始信号，w代表小波变换，那么可以得到u=wx，同样的x=w-1u，因为在实际计算时，能直接得到并且可以做初值的是原始信号，所以需要对模型中的u做一个替换，用wx替换u，那么我们可以得到：

由于A=PF而w代表小波变换，对原始图进行完小波变换后，不能直接对稀疏信号进行福利叶变换，故要对稀疏信号再进行一次逆小波变换再进行傅利叶变换，此时为使模型简单，继续用u=wx替换wx但同时要加入逆小波变换那么得到：

但是此时使A=PFw-1，下面我们就对上述模型进行计算推导。

令

令导数为0，那么式子变为：

需要求解上式中的u，我们使用逐步递归的方法，用u(k+1)代替u得到：

在分母部分用u(k)代替u(k+1)，并且去分母，得到：

利用

b(k+1)←b+(b(k)-Au(k))

得到：

μu(k+1)+A*(Au(k+1)-b-(b(k)-Au(k)))

合并同类项得到：

由

A=PFw-1

A*A=(PFw-1)*(PFw-1)=wF*P*PFw-1

A*Au(k)=wF*P*PFw-1u(k)

A*=wF*P*

A*b=wF*P*b

那么就得到：

对两边同时进行逆小波变换、傅里叶变换。得到：

具体算法如下：

b0←0，u0←0
For：k=0，1，…，N，do
bk+1=b+(bk-Auk)

(3.0)

bk=bk+1
kk=uk+1

3 数值实验

在压缩的核磁共振成像中，采样矩阵A是由A=RФ-1给出的，其中Ф是一个小波变换，R是一个部分二维离散傅里叶变换。假设一个核磁共振图像有n个像素。在我们的算法中，R是从相应于一个完整的二维离散傅里叶变换的n×n阶矩阵中随机抽取m行组成的，即R=PF，其中P是从n×n阶单位矩阵中随机抽取m行组成的矩阵，F是二维离散傅里叶变换矩阵，m▯n。所选的m行指出所选择的频率，在这些频率中，b中的观测数据被收集。m值越小，通过一个核磁共振扫描器来获取b所需要的时间就越少。在核磁共振成像中，我们有一定的自由来选择行(然而，实际的限制可能会影响我们的选择，但是这已经超出了这篇论文的讨论范围)。在我们的实验中，由R估量得出的傅里叶系数并不是在随机的频率中选取的。我们是通过以下方式来选择的：在k-space中，我们分别采取了两种不同采样方法进行试验：一种采样是沿着一定数量的从中心散开的呈辐射状的直线来选取，即半径抽样。例如图4.1显示了在一个k-space中的22*5条辐射状直线，即这是22*5 views抽样频域图；另一种则是矩形采样，例如图4.2显示了在一个中采样的区域(白色部分为采样区域)。我们发现这种选择允许我们从数量较少的采样数据中，而不是整体随机采样选择来恢复方形的核磁共振图像。实际上，在一个核磁共振成像扫描中，频率的集合以及采样的速度都是受物理和心理极限限制的，所以我们的采样方法是理想化的。

白色显示的是采样位置，这时的采样率是图1：42.38%，图2：42.95%。

图1 图2

(3.1)

我们在256×256的心脏原始图像核磁共振图像上检测我们的编码。在所有的检测问题中，采用的噪声是均值为0方差为0.01的高斯白噪声。对256×256的心脏原始图像进行110 views、66views、22 views频域抽样，重构的图像如图3每组图中图(a)为原始图像图(b)、(c)、(d)为恢复后图像。然后，我们再对256×256的心脏原始图像进行矩形区域抽样，重构的图像如图4示，每组图中图(a)为原始图像，图(b)、(c)、(d)为恢复后图像。

在数值试验中，相对误差(Relative Error)和信噪比(Signal to Noise Rations，简称SNR)用来评估重构图像的质量。相对误差和信噪比的定义如下：

(3.2)

(3.3)

(a) (b)

(c) (d)

图3(a)是原始心脏图像；(b)、(c)和(d)分别是恢复后的图像，它们的采样率分别是42.18%、26.85%和9.36%。

(a) (b)

(c) (d)

图4 (a)是原始心脏动脉图像；(b)、(c)和(d)是恢复后的图像，它们的采样率分别是42.95%、24.29%和14.44%．

图5 半径采样下相对误差、信噪比CPU时间与采样率之间的关系

图6 矩形采样下相对误差、信噪比、CPU时间与采样率之间的关系

表1半径采样下核磁共振图像重建实验结果：

表1

表2矩形采样下核磁共振图像重建实验结果

表2

4 结论

在这篇论文中，基于压缩感知理论，我们通过较少数量的观测数据来恢复核磁共振图像，对此，我们研发了一个有效算法，它克服以往求解问题的计算复杂性，提出β范数近似逼近l1范数的思想，由于β范数具有光滑性，可采用Bregman迭代正则化方法进行求解，这在l1问题的求解方面是一种新颖的方法。在真实的核磁共振图像上进行的数值实验表明，这种算法可以在采样率相对较小的情况下，使用不到一分钟的时间来恢复忠实于原图的正方形图像。通过对相对误差、信噪比和CPU时间的对比，我们知道通过使用最优化技巧，例如光滑以及更有效的线性搜索等，我们的算法速度仍然可以加快。然而，本论文未能就Bregman法的收敛速度进行系统分析，对这种方法的理论分析也缺乏详细的研究，这些问题要留待以后再解决。另外，我们认为，这篇论文中所阐述的压缩感知的算法可以被延拓到其他相关的图像以及可视化的应用中。

[1]Osher S，Burger M，Goldfarb D，Xu J，Yin W.An iterated regularization method for total variation-based image restoration，Multiscale Model.Simul[J].2005，4：460-489.

[2]赵瑞珍.压缩传感与稀疏重构的理论及应用[EB/OL].中国科技论文在线http：//www.paper.edu.cn

[3]Nocedal J，Wright S.Numerical Optimization.Springer[M].New York，2nd edition，2006.

[4]Rudin L I ，Stanley Osher.Total Variation Based Image Reatoration with free local constratints[C]. Image Processing，ICIP-94，1994，1：31-35.

[5]Gilbert A C，Muthukrishnan S，Strauss M J. Improved time bounds for nearoptimal sparse Fourier representation [ A] .Proceedings of SPIE，Wavelets XI [ C] .Bellingham WA：International Society for Optical Engineering，2005，5914：1215.

(责任编辑：王 谦)

TheBregmanMethedResearchonImageReconstructionBasedonCompressiveSensing

YANG Xiao-lan，ZHU Yong-gui，CONG Jia，LIU Ping

(School of Sciences，Communication University of China，Beijing 100024，China)

Because information such as boundaries of organs is very sparse in most MR images，compressed sensing makes it possible to reconstruct the same MR images from a very limited set of measurements significantly reducing the MRI scan duration.In order to do that，however，one has to solve the difficult problem of minimizing nonsmooth functions on large data sets.To handle this，we propose an efficient algorithm that It has overcome the computational complexity of solving thel1problem，this paper puts forwardβnorm approximationl1norm thought，becauseβnorm with slickness，one can use Bregman iterative regularization method to solve this problem.The numerical experiments demonstrate that original MR images can be reconstructed exactly from the mere 40 percent of the complete set of measurements by our approach.

compressed sensing；magnetic resonance imaging；reconstruction algorithm；The Bregman methed

2012-06-06

中国传媒大学理科规划项目(XNL1105)

杨晓兰(1985-)，女(满族)，新疆昌吉人，中国传媒大学理学院硕士研究生.E-mail：yxlan_2010@163.com

TP391

A

1673-4793(2012)04-0018-09