压缩感知磁共振成像脉冲序列的设计与实现*

黄敏，林颖，陈军波

(中南民族大学 生物医学工程学院，武汉 430074)

1 引 言

磁共振成像技术是一种无电离辐射，多层面多参数的重要医学成像方式，已广泛应用于临床和科研[1]。但MRI成像数据采集扫描过长，使该项检查的适应症大为减少。因此，提高磁共振成像速度一直是磁共振成像领域要解决的核心目标之一。

传统的快速成像方法是通过提高梯度场性能来加快 k 空间数据采集速度。然而, 受制于硬件制约和生理因素,这类方法加速效果已经接近极限。近年诞生的压缩感知(compressed sensing,CS)新理论[2-3],提供了一种不需要遵循采样定理的欠采样新方式，可明显缩短数据采集时间，并可从极少的采样数据中恢复出原始信号。

MRI图像在适当的变换域中(如小波变换域)是可压缩的，且 MRI 采集的原始数据不是直接的图像信息，而是k空间频率域信息，因此，MRI可以采用压缩感知的采集方式[4]。

要在MR扫描仪上实现压缩感知MRI新技术，需要解决两个关键问题：一是在扫描仪中实现压缩感知序列，二是对欠采样数据进行压缩感知重建。针对后一个问题，国内外已经进行了大量的图像重建方法的仿真。但在扫描仪上最终实现该技术的关键问题——压缩感知序列的设计和实现，还很少涉及。序列的设计与实现也是目前国内MR扫描仪自主研发的重点。

本研究主要就压缩感知磁共振脉冲序列的设计、编程和实现方法进行研究，设计不同的加速因子和序列，对扫描得到的数据进行离线重建，用重建结果来验证序列的正确性。

2 压缩感知MRI序列原理

压缩感知MRI与传统MRI的信号采样不同， 它只需测量少部分随机线性组合信号，通过非线性压缩感知方法从测量值中高准确度地重建出图像。将压缩感知理论应用到MRI序列中，可通过随机欠采样序列的实现，采集少量的K空间数据，提高成像速度，减少患者在扫描过程中的不适，更好地服务临床。

压缩感知MRI序列主要是要确定欠采样的方式。现有的全采样方式主要有二维笛卡尔轨迹、三维笛卡尔轨迹、二维螺旋轨迹等。在这些基础上进行欠采可分为一维随机欠采和二维随机欠采。随机欠采的相关研究有很多，如Lustig[4-6], Wang Haifeng[7]等人提出的1-D随机采样[4]，2-D随机采样[4,7]，螺旋随机采样[6]等。一维随机欠采主要是在相位编码方向随机欠采，二维随机欠采是在相位编码方向和频率编码方向同时进行随机欠采。

设计和实现压缩感知MRI序列可以有多种实现”基”。可以在自旋回波(SE)序列基础上实现，也可以在螺旋序列及其它序列基础上实现。国内MR扫描仪由于梯度系统及其它硬件配备很难实现螺旋等快速采集方式，从这种现状考虑，我们基于自旋回波序列的基础上进行设计。

在实际仪器上，对于笛卡儿采集的二维成像：由于谱仪采集数据都是等时间间隔采集数据，要在频率编码上实现欠采样必须重新调整采样数据的例程；另外，由于TR的限制，频率编码方向上即便实行了欠采样，对扫描时间也没有实质上的减少，只有相位编码方向的欠采样才能明显减少采集时间。因此，本文进行1D随机欠采的设计及实现。

对于1D加速的压缩感知MRI序列,与常规自旋回波的序列不同的是：相位编码梯度幅值不再均匀变化，而是施加‘伪随机’的变间隔梯度。压缩感知MRI脉冲序列及回波信号采集见图1。

图1 压缩感知MRI序列图

若1D欠采样加速因子为F，实际采集数据时间为全采集的1/F。一维相位编码方向上随机采样使K空间的中心附近密集采集，周围稀疏采集[5]。变密度采样的概率密度函数为：

(1)

其中r的取值范围为(0，1)，p为多项式的系数，r为采样点到k空间中心的距离；

通过以上概率密度函数设计出欠采样轨迹，根据设计的欠采样模式即可得到相应的相位编码次序，从而在序列中实现。而CS-SE序列即将SE序列的相位编码次序改变成设计的随机相位编码次序即可。

3 压缩感知MRI序列设计

本设计是在深圳贝斯达医疗器械有限公司磁共振设备上完成的，由于各方面的局限性，只能实现笛卡尔轨迹采集，并且只能在相位编码方向上欠采。

由CS 理论可知：对 K 空间数据在相位编码方向上进行随机欠采，再采用非线性重建即能重构出图像。我们先通过MATLAB设计相位编码随机欠采模式，得到不同分辨率和不同加速因子条件下，K空间欠采样的相位编码序号，并将其做成表格存储，然后将对应采样行写入序列中。

CS-SE序列是利用MR Solutions公司的磁共振谱仪对应的PPL (Pulse Program Language) 语言，来设计和编写的，主要对射频脉冲 RF，选层梯度 Gs，相位编码梯度 Gp 和频率编码梯度 Gr 进行设计，包括各脉冲的形状，脉冲宽度，时序等的控制[8]。

3.1 射频脉冲选择

磁共振2D成像要选择性激发具有一定厚度的层面，需要采用软脉冲来实现选层。选层需要在频域上为一理想窄带的矩形形状最佳，时域无限长sinc信号可满足此要求。但实际使用时通常选择一定数量的主旁瓣，以激发不同区域，常用3瓣或5瓣sinc信号。

MR Solutions公司的磁共振谱仪底层软件中，已经存在着多种固定类形的RF脉冲，如sinc型、gauss型及正切型的波形，并且各种不同波形的RF脉冲也根据带宽大小的不同而进行分类编号。设计RF脉冲需要调用这些脉冲类型。如果序列开发者觉得已有RF脉冲不满足应用，也可根据自己的需求在波形编辑器中编辑生成想要的波形并保存到RF波形文件RFstd.seq中再进行调用即可。在程序开始处需要进行RF脉冲波形调用的声明，再用对应的语句调出可能使用到的波形，如下所示：

#use RF1 "c:smisseqlibRFstd.seq" pf1

NEWSHAPE_MAC(1,pf1,"3lobe_sinc_3kHz",1332, 3000);NEWSHAPE_MAC(7, pf7,"gauss_240Hz" , 7500, 240);NEWSHAPE_MAC(8,pf8,"5sinc_7_5cos_6ms",6000, 1500)

3.2 梯度脉冲设计

MR仪器上使用的梯度脉冲是梯形状的正脉冲或者负脉冲，正负脉冲的上升和下降速度相同。为了输出不同的梯度，通常是改变梯度的强度及持续时间，各个方向上的梯度脉冲的强度大小需要定义，并存储到对应的梯度矩阵中，而持续时间也是在程序中定义的。梯度强度大小使用CREATE_MATRIX(mat, slice, phase, read)命令来定义，mat是定义的矩阵名称，slice, phase, read是固定的三个方向上的梯度值，以相位编码梯度矩阵pe_mat为例，定义如下：

CREATE_MATRIX(pe_mat, gs_on*gs_comp_s, gp_on*gp_var, gr_on*gr_comp_s)

delay( caldelay, us );

gs_on、gp_on、gr_on只能取0或者1，取0则对应方向上这段时间内的梯度值为0，取1则输出对应的梯度值，caldelay是程序中创建矩阵及计算矩阵所需要的等待时间。梯度方向上的梯度输出同样需要定义，用到MR3040_InitList()。谱仪底层软件中已有定义好的正负梯度波形，可以直接调用POSPULSE_HOLD、NEGPULSE_HOLD来组成各个方向上的梯度，也可用 MR3040_OUTPUT来分别输出各个梯度，具体代码如下：

MR3040_SetListAddress(0);

slice_list_1 = MR3040_InitList();

POSPULSE_HOLD()

if (slice_rephase = 0 & flow_comp_on == 0) {NEGPULSE(tref, clock)}

if (slice_rephase == 0 & flow_comp_on == 1) {NEGPULSE((2*tref + tramp), clock)

POSPULSE_SEC((tsel90/2-tramp/2), clock)}

slice_list_2 = MR3040_InitList();

POSPULSE_HOLD()

if ( flow_comp_on) {NEGPULSE_HOLD()}

ref_list = MR3040_InitList();

if ( flow_comp_on == 0) {POSPULSE(tref, clock)}

else {NEGPULSE_SEC_HOLD()

POSPULSE_SEC(tsel90/2, clock)}

read_list = MR3040_InitList();

POSPULSE_HOLD()

phase_list = MR3040_InitList();

POSPULSE(tref, clock)

slice_list_1和slice_list_2分别定义了90°和180°射频脉冲输出时对应的选层梯度的脉冲波形；read_list和ref_list 定义了频率编码通道输出的梯度波形；phase_list定义了相位编码梯度的输出波形。

3.3 梯度脉冲设计

为了找到适合不同场强扫描仪的压缩感知欠采样加速因子F，我们实现了F分别为1.5、2、2.5、3、4、5、8等的CS-MRI序列，实际采样矩阵大小为256×256/F，并取整数。实现过程分为以下两步：(1)用MATLAB产生1D变密度随机采样模板的概率密度函数及对应的采样模式；(2)用PPL语言实现(1)中用MATLAB设计的采样模式。

以欠采加速因子F=3为例，实际采样矩阵大小则为256×256/3(即256×84)，取r=0.15，p=6及r=0.1，p=11得到两种采样模式分别见图2(a)、(b)，图中水平方向是频率编码，垂直方向是相位编码，也就是相应的k空间填充顺序，白色表示采集，黑色不采集。

图2 压缩感知欠采样的K空间填充顺序

为了在尽可能节省时间的情况下得到较好的图像结果，就要对不同部位进行多种不同欠采样加速因子选取的实验。我们先在MATLAB中得到欠采样加速因子对应的采样模式，再将采样模式的相位编码行号写入到SMIS谱仪可读取的pph文件中，记为gp.pph。在程序开始时对该文件进行声明，并在程序中进行相应的调用。在PPL程序中具体的实现如下：

#include "C:smisincludegp.pph"

phase_encode_loop:

current_view = current_view + 1;

if (no_views=256){gp_var=gp_var-gp_inc;}

if(no_views=84){gp_var=(-gp_inc*(gp_mul2[current_view]-128));}

if (no_views=103){gp_var=(-gp_inc*(gp_mul3[current_view]-128));}

if(no_views=128){gp_var=(-gp_inc*(gp_mul5[current_view]-128));}

if(no_views=171){gp_var=(-gp_inc*(gp_mul4[current_view]-128));}

if (current_view < no_views) goto phase_encode_loop;

3.4 时序设计

MRI序列中施加的各个脉冲的顺序及各脉冲间的延迟时间的长短都会对最终的图像产生一定的影响，为了保证数据的准确性，需要对脉冲时序进行严格要求，下面以读出梯度的时序为例进行说明。

starttimer();

ret = tramp * 10 ;

MR3040_SelectMatrix( aq_mat );

MR3040_Start(CHANNEL_R);

frequency_buffer(1);

reset_frequency();

waittimer(ret);

acquire( sample_period, no_samples );

MR3040_Continue(CHANNEL_R);

delay(tramp,us);

用starttimer来提示计时开始，ret为定义的计时时间的大小，tramp是梯度脉冲上升沿及下降沿的持续时间，ret、tramp都是以us为单位的，MR3040_SelectMatrix用来选择即将输出的梯度强度矩阵，aq_mat是当前选择的读梯度强度的大小矩阵，CHANNEL_R 指的是读梯度通道，MR3040_Start(CHANNEL_R)则是执行开启读梯度通道命令，frequency_buffer频率缓存器，可存储多达到10个频率，buffer 0指的是发射频率，buffer 1指的是接收机频率，当前选择存储了接收机频率，并需要用reset_frequency()来重置一下，才能将频率切换过来，starttimer()与waittimer()是配套使用的，两者之间的持续时间即是定义的ret，包括了代码执行时间，用acquire来提示进行数据的采集，sample_period是两个数据点之间的采样间隔，no_samples 是采样点数，在参数设置时设定的，MR3040_Continue用来结束梯度的输出。

时序的精确性体现在程序开始时用MR3040_Clock(clock) 开启定时器，计时开始时，配套运用starttimer和waittimer进行时间的精确控制，及运用delay命令进行延时，及程序代码执行时间也被考虑计算进去。

程序编写完成后进行编译，生成的参数界面反映了在脉冲序列程序编写过程中设置的参数范围，一旦设置超出范围，就会进行错误提示。同时界面上的各个参数的数值与实际数据之间有各自不同的对应关系，有的是程序编写者人为定义的，有的是谱仪软件里固定的，若需要使用这些值，则需要进行一个转换。序列设计者可进一步根据自己的需求修改参数，如频率编码数、相位编码数、层数、层厚、层间距、平均次数等等，来获得自己想要的序列。

4 结果与讨论

我们在深圳贝斯达医疗器械有限公司的MR扫描仪(包括0.5T永磁型及1.5T超导型)上进行了序列的编写和调试，并进行了CS-SE序列采集试验。

对于0.5T MR扫描仪试验数据，设计了欠采样率为1/5、1/3、1/2.5、1/2、2/3、1共六种情况下的2种不同采样模式，采集得到膝盖和头部的原始数据。成像参数为：(1)头部数据：TR=520 ms，TE=20 ms，Number of Samples=256，SLICE_THICKNESS=5 mm，BW=25 KHz，SLICE_SEPARATION=1 mm，FOV=224 mm，NEX=2，TA=4 min26s； (2)膝盖数据：TR=400 ms，TE=20 ms，Number of Samples=256，SLICE_THICKNESS=5 mm，BW=25 KHz， SLICE_SEPARATION=1 mm，FOV=224 mm，NEX=3，TA=5 min7s。CS-SE序列中的相位编码步Number of Views，分别取52、84、103、128、171、256。

对全采集数据(256行)，直接采用快速傅立叶变换进行重建。对各组压缩感知欠采数据(小于256行)，按照MR扫描仪MRD文件的存储格式进行数据读取，采用本实验室在VC6.0环境下编写的压缩感知重建程序进行离线重建[9]，该程序是基于M. Lustig的非线性共轭梯度下降法，采用迭代方式求取最优值进行图像重建的。

1D随机欠采样采用了k-空间中心附近密集，周围稀疏的变密度采集方式，而k空间中心半径r及多项式次数p的选取对成像效果影响很大。对不同的欠采模式1和2下的头部数据，采用相同的Db小波结合全变差分的压缩感知方法进行重建，当迭代次数取72次，重建出一幅256×256图像需要时间为203 s，得到重建结果见图3。从左至右欠采样加速因子分别为 5、3、2.5、2、1.5、1。由于采样模式2采用的中心半径r=0.1比采样模式1的r=0.15要小的多，低频区域损失信息相对较多，加上仪器本身的噪声，在同欠采样率下，重建效果要差的多。由此可见：为了能达到最优的结果，首先要考虑中心区域应取k空间最大值的多少比例，以及各区域的采样密度的分配，才能设计出合理的欠采样模式。

图3头部CS-SE序列成像结果

(a)采样模式1；(b)采样模式2

Fig3ImagingresultofCS-SEsequenceonbrain

(a)Model 1; (b)Model 2

图像稀疏变换基有很多，如有限差分、小波基、curvelet基、Harr基、Contourlets基等等。不同部位的成像数据的稀疏度是不同的，选用不同的稀疏变换基，对图像重建速度和效果有很大的影响。我们对同一采样模式下加速因子为3的头部两层的数据以及膝盖数据，采用多种不同的稀疏变换后重建图像，结果见图4。从左至右稀疏方式分别为全变差分、5层Db小波、5层Symmlet小波、5层Battle小波、8层Db小波。迭代次数选用45，重建出256×256图像所需要时间分别为：头部13、94、92、118 s，膝盖13、92、100、118 s。注：最右边的是全采样，采用快速傅立叶变换重建。由此可见：对于稀疏程度高的膝盖数据，只需选择简单的差分运算便可得到很快且较好的结果；对于稀疏程度相对较低的头部数据，采用5层Db4的稀疏变换方式较优。

当然，迭代次数同样影响着图像重建效果。迭代次数越多，重建结果也就越接近全采样值。不同的稀疏变换下的迭代时间也各不相同，如：全变差分变换每增加一次迭代时间增加0.29 s，而5层的Db4变换一次迭代要2 s。所以，在不同稀疏变换的选择下，若是重建效果差不多，则尽可能选择重建时间短的变换；而且对于不同部位的数据，最佳的迭代次数也会不同。我们还在1.5T 超导MR扫描仪上进行了序列调试和采集数据的试验，成像参数同0.5T扫描仪。当加速因子分别为5、4、3、2.5、2，1时(即采样行数为52、64、84、103、128，256)，头部两个部位以及膝盖的重建结果见图5。根据上述扫描数据及重建结果，我们认为：在贝斯达公司的扫描仪上运用CS-SE序列后，综合考虑扫描时间及图像效果两个方面，头部合适的欠采样加速因子为2，膝盖则为4。此时将采集并重建图像的时间分别缩短到SE序列的1/2及1/4，扫描CS-SE序列获得膝盖横断位数据需要1 min 17 s，头部横断位数据需要2 min 23 s，全采膝盖需要5 min 7 s，全采头部需要4 min 26 s，CS-SE大大节省了数据采集时间。不同部位，由于图像本身的稀疏程度不同，膝盖部位可采用比头部高的欠采样加速因子。

图4不同稀疏变换的重建结果

(a) 头部层面1 ；(b) 头部层面2；(c) 膝盖
Fig4Reconstructionresultfromdifferentsparsetransforms

(a)Slice 1 on brain; (b)Slice 2 on brain;(c)On knee

图5 1.5T超导MR扫描仪CS-SE序列成像结果

由此可见：实际仪器采集CS数据可以实施的加速因子，与很多文献中采用‘仿真’的理想数据得到的可欠采样十倍甚至几十倍的结论相差甚远。因此，我们必须考虑实际采集数据中不可预见的各种噪声，不能以仿真数据作为推断。压缩感知成像技术要真正用于商业扫描仪上，还得进行大量的试验。

5 结论

目前，国内MR扫描仪的核心部分(包括谱仪，脉冲序列，线圈等)的自主研发还处于薄弱环节，压缩感知序列更没有应用于商用扫描仪上。本研究介绍的基于压缩感知的磁共振序列，对序列的设计和编程，图像的重建进行了介绍，可为国内MR开发领域在脉冲序列设计及实现方面提供参考。

CS-SE序列最终要嵌入到商用扫描仪，必须对于不同部位分析出最佳的欠采样因子，并直接固定重建程序中不同部位对应的稀疏方法以及最佳迭代次数，这样才能既实现加速扫描，获得重建质量好的MR图像。我们下一步将研究更优的压缩感知重建算法来适应更大的欠采样加速因子，并将CS重建算法集成到MR扫描仪上，进行在线重建，争取将压缩感知成像应用到临床成像序列中。

感谢深圳贝斯达医疗器械公司、湖北省医学信息分析及肿瘤诊疗重点实验室，以及国家民委脑认知实验室的支持！