毕业论文-基于压缩感知的磁共振成像欠采样恢复方法对比与分析

本科生毕业论文（设计）题目:基于压缩感知的磁共振成像欠采样恢复方法对比与分析姓名:张作恒学院:基础医学院专业:生物医学工程班级学号:生物医学工程指导教师:王伟职称:实验师2014 年 6月 10日南京医科大学教务处制摘要摘要磁共振成像（MRI）技术属无创伤、无射线检查，具有对比分辨率高，多方位、多参数采集及功能成像等优势，然而，MRI技术是一种速度相对较慢的成像模式，数据采集的时间过长是其最大缺点。压缩感知（Compressed Sensing，CS）理论是近年来新兴的研究方向，它指出，只要信号是可压缩的或在某个变换域是稀疏的，那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上，然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号，可以证明这样的投影包含了重构信号的足够信息。在该理论框架下，采样速率不决定于信号的带宽，而决定于信息在信号中的结构和内容。目前，基于CS理论的MRI技术尚处于理论研究阶段，很多关键问题需要解决，主要包括：k空间部分采样方式的确定，不同类型成像的稀疏性表达，相关的快速优化算法，以及从医学临床的角度对重建图像质量进行评价。本文首先简要介绍了MRI技术的基本原理和CS理论的基本概念。从k空间重建原理、k空间的欠采样填充轨迹、MRI图像的稀疏化方法、相关的恢复重建算法的角度对基于CS理论的MRI技术的发展过程进行了详细说明；之后选取脑部断层MR图像进行仿真，对k空间分别进行矩形FOV、放射状、螺旋状和Propeller轨迹欠采样，并通过5种恢复算法进行恢复重建，并从均方误差、峰值信噪比、结构化相似度和运算时间等客观指标方面进行比较分析，最终得出k空间欠采样选择放射状欠采样轨迹、稀疏化方法选择6层haar小波变换，恢复算法选择线性Bregman法可以在较短恢复时间内取得较优恢复效果的结论；最后对基于CS理论的MRI技术下一步的研究和应用进行展望。关键词：压缩感知，磁共振成像，k空间欠采样，恢复算法72AbstractIIAbstractMagnetic resonance imaging (MRI)inspection has no trauma and no X-ray.With high contrast and resolution, it possesses the advantages of multiparameter, comprehensive and functional imaging.However, magnetic resonance imaging (MRI) technology is a kind of relatively slow imaging mode.Its biggest drawback is the long scan time.Theory of Compressed Sensing (CS) is proposed in recent years. It pointed out that as long as the signal is compressible or in a transform domain is sparse, a observation matrix unrelated to transform bases can project high dimensional signal to a low dimensional space, and then solving an optimization problem can reconstruct the original signal from these few projection with high probability. It can be proved that the projection contains enough information to reconstruct signal.Under the theoretical framework, the sampling rate is determined by the structure and content of information in the signal, instead of the bandwidth of the signal.So far, the MRI technology based on the theory of the CS was still in the stage of theoretical research. A lot of key problems need to be solved,including: the k-space undersampling way, sparse representation of different image, the relevant fast optimization algorithm and the clinical evaluation of reconstruction image quality.This paper first briefly introduces magnetic resonance imaging (MRI) technology and the CS theory.From the aspects of k-space reconstruction, k-space undersampling, sparse methods and reconstruction algorithms of MRI, this paper reviews the development of MRI techniques based on the theory of the CS. After performing brain MR image simulation, using the FOV,radial,spiral,Propeller undersampling trajectory of k-space and 5 recovery algorithms, the reconstructed images are compared each other and analysed from mean square error, peak signal-to-noise ratio and structural similarity and computing time to concluded that choosing radial undersampling trajectory of k space, haar wavelet transform, linear Bregman recovery algorithm can recover images well in a shorter recovery time. Finally, the research and application of MRI technology based on the theory of CS is also prospected.Key words: Compressed Sensing, magnetic resonance imaging, k-space undersa-mpling, reconstruction algorithm目录目 录摘要IAbstractII第 1 章绪论11.1磁共振成像技术与压缩感知理论11.1.1磁共振成像技术11.1.2压缩感知理论11.1.3基于压缩感知理论的磁共振成像21.2主要研究内容和章节安排3第 2 章磁共振成像原理52.1核磁共振理论52.1.1原子核的磁性52.1.2磁共振的基本原理62.1.3弛豫和弛豫时间72.2磁共振成像82.2.1磁共振成像的空间定位82.2.2k空间与磁共振图像重建92.3本章小结10第 3 章压缩感知基本理论113.1信号的稀疏表示113.2压缩观测143.2.1随机矩阵153.2.2确定性矩阵163.2.3结构随机矩阵163.3重建算法173.3.1基追踪法173.3.2匹配追踪法183.3.3正交匹配追踪法193.4本章小结22第 4 章基于CS理论的MRI技术234.1K空间的欠采样填充轨迹234.2MRI图像常用稀疏化方法254.3MRI欠采样恢复算法254.3.1直接填零法264.3.2二维OMP法264.3.3快速收缩阈值法264.3.4线性Bregman法294.3.5共轭梯度法314.4本章小结32第 5 章压缩感知效果对比335.1图像重建质量的客观评价方法335.1.1均方误差（mean square error，MSE）335.1.2峰值信噪比(peak signal-to-noise ratio，PSNR)335.1.3结构化相似度（structural similarity，SSIM）335.1.4运算时间345.2实验与分析345.3本章小结42第 6 章结论和展望43附录 A 图表索引45附录 B 评价数据47参考文献53综 述57致 谢63绪论第 1 章 绪论1.1 磁共振成像技术与压缩感知理论1.1.1 磁共振成像技术磁共振成像（MRI，magnetic resonance imaging）是根据生物体磁性核（氢核）在磁场中表现出的共振特性成像的高新技术1。MRI是以核磁共振（Nuclear Magnetic Resonance，NMR）这一物理现象为理论基础的，之所以省去“核”字，是为了突出这一检查技术不存在对人体有害的电离辐射的优点，使之区别于要使用X射线的放射科检查以及要使用放射性核素的核医学检查。从NMR的发现到MRI设备的诞生，这中间经历了几代物理学家及医学家长达数十年的辛勤努力。1946年美国物理学家布洛赫（Felix Bloch）和珀赛尔（Edward Purcell）分别独立地发现核磁共振现象，并荣获1952年诺贝尔物理学奖。1972年美国内科医生达玛迪安（Raymond Damadian）提出利用磁共振原理测定活体组织的纵向弛豫时间T1和横向弛豫时间T2的方法，以及用于医学诊断的设想。1973年，美国州立大学教授劳特伯（Paul Lauterbur）首次使用梯度磁场进行空间编码，产生了第一幅磁共振图像。从此，成像概念发生了质的飞跃。尤其近年来，随着超导技术、低温技术、电子技术、成像技术和计算机等相关技术的进步，MRI技术得到了飞速的发展。MRI系统已经成为当今医学影像领域最先进、最昂贵的诊断设备之一。目前以分子探针和谱成像为标志的MRI分子影响学使成像进入一个新的层次。然而，数据采集时间较长仍然是磁共振成像技术的最大缺点，由于成像速度慢，使该项检查的适用范围大为减少，例如不适合运动器官和危重病人的检查；对于躁动或者丧失自制能力的患者，如不使用镇静剂，也是难以成像；儿科的某些应用同样受到限制。缩短成像时间不仅可以提高效率和病人的舒适度、减少时间依赖性伪影，还是实现心血管检查、功能信息获取、实时温度检测与介入手术成像等动态成像的关键。因此缩短成像时间一直以来都是磁共振成像技术发展的重要目标之一。现在的研究者致力于通过分析k空间图像数据的冗余性，从而减少成像所需的采集数据，达到缩短成像时间的目的。本文论述的部分k空间重建就是针对部分k空间采集技术的重建方法。部分k空间采集技术就是在k空间数据时，只采集k空间数据，以达到缩短扫描时间的一种应用广泛的扫描采集技术。1.1.2 压缩感知理论近几年来出现的一种新颖的理论Compressed sensing（也称为Compressive sampling)2-7使得采样速率低于信号带宽两倍时精确恢复信号成为可能。它指出，只要信号是可压缩的或在某个变换域是稀疏的，那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上，然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号，可以证明这样的投影包含了重构信号的足够信息。在该理论框架下，采样速率不决定于信号的带宽，而决定于信息在信号中的结构和内容。事实上，压缩感知理论的某些抽象结论源于Kashin创立的范函分析和逼近论8，最近由Candes，Romberg4，Tao9和Donoho5等人构造了具体的算法并且通过研究表明了这一理论的巨大应用前景。压缩感知理论是应用数学与信号处理领域中一个非常新的研究方向，自从2006年有了正式的论文发表之后，迅速引起了国内外相关领域研究者的高度重视，在医疗成像10、光学/遥感成像11,12、无线通信13,14、目标识别15、雷达探测1619、地震数据重建20、图像融合21、冷链物流检测22、声音信号的采集与压缩23,24等领域均受到高度关注，并被美国科技评论评为“2007年度十大科技进展”，D Donoho因此还获得了“2008年IEEE IT学会最佳论文奖”。该领域的先驱者有加州大学洛杉矶分校的Terence Tao、加州理工大学的EmmanuelCandes、斯坦福大学的DavidDonoho、以及莱斯大学的RichardBaranink等。国内关于这方面的研究则刚刚起步，发表的论文甚少，但中科院电子所、西安电子科技大学、燕山大学、西南交通大学、华南理工大学、北京交通大学等单位的一些研究所已经着手开始研究。在硬件实现上，Rice大学的Baranink教授等研制出单像素相机25和A/I转换器26，吸引了国内外众多媒体的眼球。随后也有一些硬件方面的相继报道，例如，麻省理工学院Wald教授等人研制出MRIRF脉冲设备27，麻省理工学院Freeman教授等人研制出编码孔径相机28，伊利诺伊州立大学Milenkovic等人研制出DNA微阵列传感器29等等，然而，由于缺乏有效的压缩感知矩阵判别理论，除Rice大学的单像素相机（硬件成本昂贵，重建算法效率低下）外，其他硬件均缺乏严格的理论分析。经过近两年的发展，压缩感知在理论方面已经取得了许多重要的成果，许多研究者已经将之投入到实际应用当中，如信息、医学等学科。1.1.3 基于压缩感知理论的磁共振成像如今，成像在医疗诊断方面起到了着非常重要的作用。现在临床MRI的几个研究热点是加速MR功能成像、卡尔曼滤波重建磁共振实时动态成像、插值法加速二维多层面成像、分割算法快速成像、加速超极化3D MR波谱成像、加速MR电影成像、高空间分辨率并行采集成像等30。在MRI领域，由于扫描仪器所采集的不是直接的图像像素，而是由图像经过全局傅里叶变换将原始采集的时域图像转化得到的频域图像。运用压缩感知理论可以大大减少采样数据量，从而为后续数据传输、处理和存储减小压力，降低对硬件的需求；此外，该技术创新地改变医学信息的获取方式，可以将速度提高到原来的几千倍，在能保证令人满意的空间分辨率的情况下，缩短扫描时间，减少病人在扫描过程中的不适；在保持原有采样时间的基础上，可以进一步提高现有设备的空间分辨率，使医生的诊断更为可靠，故其在临床MRI的应用备受关注。磁共振成像应用范围广，对于不同成像部位以及成像的不同类型，可能会对k空间数据欠采样方法、图像稀疏化方法以及恢复重建算法的选择产生影响，因此，针对成像部位和类型的差异选择合适的压缩感知方案十分重要。1.2 主要研究内容和章节安排本文主要针对压缩感知理论在磁共振成像技术中的应用进行了研究，并且通过多种磁共振图像进行仿真。研究内容主要包括以下几个方面：（1） 从数据稀疏化、数据采样、数据恢复三个方面学习压缩感知理论的基本概念；（2） 磁共振成像技术的基本原理以及数据采样方式；（3） 分析磁共振成像数据的稀疏性并对比使用多种稀疏变换的成像效果；（4） 磁共振成像数据k空间欠采样轨迹的选择；（5） 基于压缩感知的磁共振成像的信号重构恢复算法；（6） 从均方误差、峰值信噪比、结构化相似度和运算时间四个方面分析压缩感知应用于MRI的优越性，并得出最优的稀疏化方法、k空间欠采样轨迹以及恢复算法。本文共分为六章，主要内容安排如下：第1章是绪论，对磁共振成像技术与压缩感知理论的基本原理、发展历史与现状以及本课题的研究意义做了简要介绍；第2章介绍了磁共振成像技术的基本原理，重点对空间编码技术、k空间的特性做了较为详细的介绍，为之后的k空间欠采样方法的研究做好必要准备；第3章阐述了CS理论的基本要素，CS理论中常用的欠采样方法、稀疏化方法以及恢复重建算法，结合lena图像和批量离子浓度检测分别说明稀疏化方法和恢复算法的优劣；第4章介绍基于CS理论的MRI技术的几种常用k空间欠采样方法、稀疏化方法以及恢复重建算法；第5章介绍均方误差、峰值信噪比、结构化相似度和运算时间四个客观评价指标，将多种欠采样轨迹、稀疏表示和恢复算法相互结合，在四种欠采样率情况下恢复图像，并从上述评价指标方面进行对比，得出最优压缩感知方案；第6章对本文所做的工作进行总结和展望。磁共振成像原理第 2 章 磁共振成像原理本章主要介绍了磁共振成像的基本原理31，主要包括核磁共振、磁共振信号的产生与检测、磁共振成像的空间编码与图像重建等内容。对k空间采样过程以及k空间的特性做了较为详细的介绍，为后面基于CS的MRI的K空间欠采样做好理论上的准备。2.1 核磁共振理论世界是物质的，物质是由分子构成的，分子又是由原子构成的。原子由原子核和核外电子组成，原子核由带正电荷的质子和不带电荷的中子组成。核磁共振所要研究的对象就是原子核，而且是具有磁性的原子核。磁共振成像是以核磁共振为理论基础的，为了突出这一检查技术不存在对人体有害的电离辐射的优点，使之区别于要使用X射线的放射科检查以及要使用放射性核素的核医学检查，因此将“核”字省去。2.1.1 原子核的磁性在微观世界中，电子、中子、质子和原子核等微观粒子除了具有一定的大小、电荷、质量等属性外，还有一种固有属性-自旋，用自旋角动量（spin angular momentum）描述。为便于理解，微观粒子的自旋运动可以简单地看成微观粒子的自旋，虽然实际情况并非如此。微观粒子除了具有自旋角动量外，还具有轨道角动量（orbital angular momentum）。原子核的自旋的取值是离散的、不连续的、也就是说是量子化的：(21)其中，为原子核的自旋量子数（spin quantum number），取整数和半整数：，为普朗克常数。大小取决于值，不同的原子核值不同。处于静磁场中的原子核，它的自旋在空间所取的方向也是离散的、不连续的，具有空间量子化的性质，即在静磁场方向只能有若干个特定的取向，取向的数量取决于值的大小，为种，这也使得在静磁场方向（方向）的投影值有个不同的分量，分别为：(22)其中，为核自旋磁量子数（spin magnetic quantum number）。原子核可看成是具有一定质量与体积的均匀带点球体，因此原子核的自旋运动会产生绕核心旋转的环形电流，继而会在其周围空间产生磁场，所以自旋不为零的原子核（简称自旋核）就会具有一定的磁性。为描述自旋核在其周围空间所产生的磁场特性，引入一物理量磁矩（magnetic moment）。自旋核的磁矩和自旋都是由原子核的自旋运动引起来的，它们之间存在着一定的比例关系：(23)为比例系数，称为磁旋比（guromagnetic ratio）；为一个取决于原子核种类的无量纲的数，称为朗德因子；为电子的电荷数；为光速；为质子的质量。与原子核的自旋一样，原子核的磁矩在静磁场中也存在种可能的取向，因此核磁矩在静磁场方向（方向）的投影值有个不同的分量：(24)自旋不为零的原子核都是磁性核，也只有磁性核才能和静磁场相互作用产生磁共振。在生物组织中，存在很多的磁性核，如H、C、N、F、Na、P但是由于氢原子占到生物组织原子数的2/3，氢核的磁化强度也是人体常见磁性核中最高的，所以目前能用于临床MRI的却只有氢核。当人体被置于磁体内，人体内部的磁性核就会受到静磁场的作用，使得其运动状态发生改变。在人体进入磁场之前，磁性核的磁矩所选择的方向（即所谓的取向）处于一种杂乱无章的状态，磁矩沿空间各方向等概率分布。当磁性核处于静磁场中时，就会在静磁场的作用下，只沿空间某几个特定方向分布，也就是所谓的空间量子化。例如，氢核的自旋，它在磁场中的取向就只有两种，一是顺着磁场方向，能量状态较低，为；另一是反着磁场方向，能量状态较高，为，它们之间的能量差为：(25)其中，为静磁场强度，为时氢核所具有的能量。在静磁场的作用下，核磁矩会有特定的空间取向，使得核磁矩和静磁场存在特定的夹角，因此使得原子核的磁矩以夹角在以静磁场为轴做拉摩尔旋进（Larmor precession），旋进的角速度为：(26)其中，为拉摩尔旋进频率，为磁旋比，为静磁场强度。2.1.2 磁共振的基本原理共振是自然界普遍存在的一种物理现象，如音叉的共振，磁共振等。音叉的共振是在外来声波的激励下产生的，而磁共振则是在外来电磁波的激励下产生的。处于静磁场中的氢核会有两种取向，取向不同，它们之前的能量差为。如果外界施加的电测波的能量（量子）正好等于不同取向的氢核之间的能量差，则处于低能态的氢核就会吸收电磁波能量跃迁到高能态（受激吸收），这就是所谓的核磁共振，即处于静磁场中的磁性核受电磁波的作用而产生的不同能级之间的共振跃迁现象。磁共振中，所施加的电磁波又叫射频波（radio frequence wave），简称RF波。如果射频脉冲的频率与质子绕轴进动的频率相等，那么就会发生共振。共振会导致两个同时发生的作用：一是低能级的质子吸收了射频脉冲的能量后跃迁到了高能级，纵向磁化矢量变小；二是受射频脉冲磁场的磁化作用，进动的质子趋向于射频磁场方向而变为同步、同速运动，这样在平面上叠加起来，就形成了一个新的宏观磁化矢量，即横向磁化矢量，用表示。质子同时绕磁场以频率，又绕射频波的磁矢量进动，这导致净磁化矢量由轴到平面的螺旋形运动，即章动（nutation）。在各种射频脉冲作用的情况下，的大小都不会变化，只是随着其偏离轴的程度加大，会在平面内产生越来越大的横向磁化分量。随着翻转角度的变化，逐渐增加，当到达平面后，达到在平面上的最大分量。当射频脉冲继续持续作用时，就会继续偏离平面向轴负方向翻转，其在平面上的分量又逐渐减少。如果在平面设置一接收线圈，根据法拉第电磁感应原理，通过闭合回路的磁通量发生变化时，闭合回路内产生感应电压，感应电压的大小与磁通量的变化率成正比。于是在线圈两端就会感应出交流电动势，这个电动势就是线圈接收到的磁共振信号。2.1.3 弛豫和弛豫时间弛豫（relaxation）是指当射频脉冲消失后，质子由激发状态返回到平衡状态的过程。一旦关闭射频脉冲，质子将会重新沿静磁场方向排列，并释放出多余能量。因此，在关闭射频脉冲以后，将会发生两种情况：（1）自旋的质子将回复到最低能态；（2）自旋的质子彼此之间将出现相位差。这是由两种同时发生但互相独立的过程引起的：（1）纵向弛豫（longitudinal ralaxtion），是指纵向磁化逐渐恢复为的过程；（2）横向弛豫（transverse ralaxtion），是指横向磁化逐渐衰减恢复为零的过程。Bloch从实验发现，弛豫过程中磁化强度偏离平衡状态的程度越大，则其恢复的速度就越快。这一规律在旋转坐标系中表述成如下形式：(27)(28)考虑样品受到的是90RF脉冲作用，且把90RF脉冲过后的时间点作为弛豫过程开始的起点，因此时，据此可推出和随时间的变化规律：(29)(210)其中， 、为引入的两个系数，分别表示纵向弛豫时间常数和横向弛豫时间常数，简称为纵向弛豫时间和横向弛豫时间。2.2 磁共振成像2.2.1 磁共振成像的空间定位在纵向磁场的全身MRI系统中，磁体的轴一般定义为磁体的轴向，之后通过右手法则定义轴、轴。MRI系统通过三个梯度子系统来实现成像层面的选择和空间编码。（1） 层面选择若只考虑在方向选层，沿主磁场方向叠加一个线性梯度磁场，沿轴各层面上质子的进动频率可表示为：(211)若所施加的激励RF脉冲的角频率为：，则只有这一层的核受到激发产生自由感应衰减（FID）信号，这样就选出了一层。由于激励脉冲RF的频率有一定的范围，则选中的一层将有一定的厚度。若梯度一定，则RF脉冲频宽越大层越厚。（2） 相位编码所谓相位编码（phase encoding），就是先利用相位编码梯度造成质子有规律的进动相位差，然后利用此相位差来标定体素空间位置的方法。加入相位梯度磁场后，在的作用下，相位编码方向上各体素处于不同强度的磁场中，因此该方向上各磁化强度矢量将以不同频率进动，其进动频率为：(212)进动频率的不同必然导致进动相位的不同。设相位编码梯度的持续时间为，则时间后相位编码方向上个体素的进动相位为：(213)不同的坐标处，磁化矢量旋过角度不同，其相位差为：(214)与成正比，也就是说可以通过相位差来确定自旋核所处的空间位置，这就是空间坐标用相位进行了编码。（3） 频率编码通过选层及相位编码确定了体素的坐标和坐标，但相同的体素还不能区分。经时间撤销后，加入沿方向的梯度场，磁化矢量旋进角频率为：(215)不同处旋进角频率不同。实际上频率编码与相位编码没有本质的不同，不过各体素的坐标通过频率差别来确定。采集到的MR信号是断层内所有体素自旋核产生信号的总和，要实现重建必须把信号重新按不同的频率和相位分解开，才能得到断层的每一位置自旋核所产生信号的强度。在MRI中普遍使用的方法是傅立叶变换成像图像重建。2.2.2 k空间与磁共振图像重建k空间填充技术是MRI仪生产过程中的关键技术之一，也是做成像分析的重要手段。在MRI中，相位和频率标记了位置，隐含了时间。每次信号采集的都是总信号，为了更快更方便对采集的的时域信号进行变换，人为地建立一个抽象的二维k空间。以自旋回波为例讨论，序列开始与90脉冲同时，启动选层梯度场后启动相位编码梯度场，持续时间为，90脉冲后间隔时间加180脉冲，重聚后产生回波，此时加梯度脉冲，对坐标不同的体素频率编码，持续时间为，打开的同时进行数据采集，若图像矩阵为，在时间内以一定的时间间隔如，时间段进行采集，每次采集个信号，对应一次相位编码，这个数据形成以时间为变量的数据矩阵的一行。一次采集结束后再加RF脉冲重新开始，每次大小不变采集时间不变，所以相同坐标不同坐标的频率是相同的。而相同坐标的不同体素共振频率不同，经一维傅立叶变换后用频谱区分。每次采集到的一行数据经一维傅立叶变换后，某点的频域信号等于对应该坐标点全部坐标所发出信号的总和即所有相位的叠加。把采集到的每个信号变换为以为自变量的数据填充到以一定顺序储存数据的k空间中去。可见k空间是以空间频率为坐标轴的空间。横轴代表频率编码，纵轴代表相位编码，列数等于取样点数，行数等于相位编码步数。k空间内的空间频率分布是中心频率为零，距中心越远则频率越高。k空间内信号的幅度分布是k空间中心部分所对应的MR信号幅度大，主要形成图像的对比度。k空间外围部分所对应的MR信号幅度低，主要形成图像的分辨力。由于MRI成像设备价格昂贵，没有亲自进行磁共振成像的实验条件，所以本课题通过对磁共振图像进行二维傅立叶变换，反推出相应k空间。如图21所示，其中(a)、(b)、(c)分别为脑部断层MR图像、脑部血管造影MR图像、腿部断层MR图像，(d)、(e)、(f)分别为对应的k空间图像。从图中可以看出，大部分数据集中在k空间的中心处，而k空间的外围部分的幅值远低于k空间的中心处。 (a)脑部断层图像 (b)脑部血管造影图像 (c)腿部断层图像 (d)脑部断层k空间 (e)脑部血管造影k空间 (f)腿部断层k空间图21 MRI图像及相应k空间图像2.3 本章小结本章主要介绍了磁共振成像的基本原理，磁共振成像主要包括三个过程：激发过程、弛豫过程和成像过程，对于磁共振成像的空间编码、k空间采样过程以及k空间的特性做了更为详细的介绍，并提出通过对磁共振图像进行二维傅立叶变换，反推出相应k空间，实现基于CS的MRI的k空间欠采样软件仿真模拟的方案，为后续研究做好准备工作。压缩感知基本理论第 3 章 压缩感知基本理论压缩感知（CS）理论主要包括三个方面信号的稀疏表示、观测方法以及重建算法32。所谓信号的稀疏表示，指的是当我们把原始信号根据特定变换基投影变换之后,变换后的系数中大多数的绝对值很小, 也就是说变换向量是稀疏的, 则可以认为变换向量是原始信号的另一种更简洁的表现形式33，这是压缩感知的先验条件。观测矩阵在选取时，要求所选矩阵必须满足约束等距性（RIP）条件34，再通过原始信号与该矩阵的乘积得到线性投影测量。最后，选择合适的重建算法恢复原始信号。本章首先介绍压缩感知的理论框架，对稀疏表示、观测矩阵以及几类重建算法作简要介绍。3.1 信号的稀疏表示稀疏表示作为一种有效的信号表示模型，能够用尽可能简洁稀疏的方式表示信号。表示系数中较少的非零分量揭示了信号的主要结构和本质属性，从而为信号的后续处理带来了很大的便利。假设信号可由某一组原子线性表示为：(31)将式(31)改写成矩阵形式，可以得到：(32)其中，为稀疏矩阵，是在一个稀疏变换域中的投影系数，展开稀疏系数向量。若系数向量中只有个非零分量，其余分量均为零，则称信号在变换下是-稀疏（-sparse）的。其中，称为信号的稀疏度（sparsity）。然而，实际的自然信号或图像往往都不是严格稀疏的，在信号的表示系数中，大部分的分量都接近于零，只有少部分的分量幅值较大，这时信号被称为是可压缩的（compressible）。由上述内容可知，信号的稀疏性或可压缩性依赖于所选择的变换，不同的变换对某一特定信号的稀疏表示能力通常是不同的。到目前为止，信号的稀疏变换发展大体经历了傅立叶变换、小波变换35,36、多尺度几何分析37和过完备冗余字典38-40等几个主要阶段。傅立叶变换是信号和数字图像处理的理论基础，揭示了时间函数与频谱函数之间的内在联系，反映了信号和图像在“整个”时间范围内的“全部”频谱成分。然而虽然有很强的频域局部化能力，但是却不具有时域局部化的能力，只适合于平稳信号或图像的表示和分析。小波分析将信号和数字图像处理带入到一个崭新的领域，是处理非平稳信号和图像的有力工具。它通过牺牲一定的频域定位性以获得时频局部性的折中，对一维有界变差函数表现出良好的逼近性能。但是小波变换在一维情形下所具有的优异性能并不能简单地推广到二维或高维信号，不能最优表示含线或者面奇异的高维函数。多尺度几何分析是继小波变换之后的新一代信号分析工具，具有多分辨、局部化和多方向性等优良特性，较好地克服了小波变换不能有效表示含线或面奇异的高维函数的不足，能更加有效地表示和处理图像等高维空间数据。近年来，学者们又提出了以过完备冗余字典为基础的稀疏表示理论，旨在实现图像的自适应稀疏表示38-40。在这种理论框架下，传统的基函数由“原子”所取代，表示基则更名为冗余字典。由于字典的构造具有一定的灵活性，正交基可以视为特殊的字典，多个正交基的联合也可以构成字典，甚至不同类型的函数族也可以存在于同一字典中，此外还可由某一生成函数通过参数离散化的方式构造字典，因此图像的自适应稀疏表示成为可能。信号的稀疏表示是压缩感知的理论基础，在一定程度上影响着测量值个数的选取和重建结果的精度。具体地，当重建精度一定的情况下，信号越稀疏，需要的测量值个数就越少；而当测量个数一定时，信号越稀疏，重建的精度就越高。因此在压缩感知重建过程中，信号稀疏变换的选择对重建结果的精度至关重要。图像稀疏性是变换系数满足诊断要求重建的百分比，然而诊断质量是主观的10，为了能够较客观地评价图像稀疏性，本文选取具有代表性的512512像素点的lena图像，如图31所示，分别通过离散傅立叶变换、离散余弦变换、3层离散haar小波变换和x方向有限差分变换进行稀疏表示，结果如图32所示。之后将稀疏域系数由大到小排序，分别选择前2%、5%、10%和20%的系数进行重建，重建效果对比图如图33所示，由观察可得，离散傅立叶变换、离散余弦变换和3层离散haar小波变换在利用前10%的系数恢复时已经可以获得较好的图像，而x方向有限差分变换在利用前20%的系数恢复时效果仍不理想，会产生明显的条纹状干扰，因此，x方向有限差分变换与另外三种变换相比稀疏效果较差。此外，稀疏变换原理不同对图像恢复造成的影响也不相同，离散傅立叶变换与离散余弦变换将图像域变换到由正弦函数为基底的稀疏域中，在使用较少系数恢复后，会产生较为明显的模糊，这是由于数值较小的高频被过滤的缘故；而3层离散haar小波变换是将图像利用相邻像素反复求和求差，对图像进行压缩，如图32中(c)所示，图像信息被集中在左上角的分辨率较小的低频区域内，因此在使用较少系数恢复时，会产生类似马赛克的干扰，即图像分辨率的降低；x方向有限差分可视为小波变换的高频分量，由于缺少低频分量，使用较少系数恢复时会产生明显的条纹状干扰。图31 512512像素的lena图像 (a)离散傅立叶稀疏变换 (b)离散余弦稀疏变换 (c)3层haar小波稀疏变换 (d)x方向有限差分稀疏变换图32 lena图的四种稀疏变换系数图像2%5%10%20%离散傅立叶变换离散余弦变换离散小波变换有限差分变换图33 lena图由四种稀疏变换的最大2%、5%、10%和20%系数恢复效果对比3.2 压缩观测设向量空间中的一维离散实值信号在变换下是稀疏的或是可压缩的，则可通过压缩感知测量过程表示为：(33)其中，为个非自适应的测量向量，为相应的测量值，且满足。零测量矩阵，测量值向量，则压缩感知过程可表示为如下矩阵形式：(34)由于，则式(34)可等价表示为：(35)其中，CS信息算子为，压缩感知观测向量的矩阵表示如。图34 压缩感知观测向量的矩阵表示示意图测量值向量包含了稀疏信号在变换下的个线性测量值，能否通过无失真地重建原始信号取决于测量矩阵是否能够保持的原始结构信息。显然，如果测量过程破坏或丢失了中的必要信息，精确重建是不可行的。要想使原始信号完全重建，必须保证测量矩阵不会把两个不同的K-稀疏信号映射到同一个测量值向量上，这就要求从测量矩阵中抽取的任意K个列向量构成的子矩阵是近似正交的。Candes和Tao2,41,42等指出如果测量矩阵满足限制等距性条件（Restricted Isometry Property，RIP），则原始信号可由相应的测量值向量精确重建。然而，直接构造一个满足RIP条件的测量矩阵是一个组合复杂度问题，当前的主要构造方法有：随机矩阵、确定性矩阵和结构随机矩阵43。3.2.1 随机矩阵测量矩阵构造常用的方法是通过随机的方式按照某种特定的概率分布选取测量矩阵的元素。当测量值个数与信号的稀疏度K之间满足一定的不等式关系时，相应的随机测量矩阵可高概率地满足RIP条件。此类测量矩阵主要包括随机高斯矩阵和随机伯努力矩阵。随机高斯矩阵，的每个元素均独立地采样于正态分布，如果满足，则以较高概率满足RIP条件；随机伯努力矩阵，的每个元素均独立地采样于对称伯努力（Bernoulli）分布，即，其中三者应满足的关系与高斯测量矩阵相同。这类随机测量矩阵的优点是完全重建所需的测量值个数达到了理论最优值，且几乎与任意的稀疏变换矩阵都不相关，具有一定的普适性。但缺点是矩阵是稠密的，所需的存储空间很大，并且由于其非结构化的特点导致实际的计算复杂度很高，在对图像等含有大量数据的高维信号进行测量和重建时，实际的存储和运算几乎不可能实现。3.2.2 确定性矩阵为了克服随机高斯/伯努利矩阵的相关缺点，一些学者提出了确定性矩阵。其中一类测量矩阵是以随机托普利兹和随机循环结构矩阵（Random Toeplitz and Circulant structured Matrices）作为测量矩阵。该类矩阵可表示为：(36)或(37)其中的元素为独立地选自与某同一概率分布，为随机托普利兹矩阵，为随机循环矩阵。这类矩阵的特点是：（1）只需产生个独立的随机变量，相对于随机高斯/伯努力矩阵，所需的存储空间大大减少；（2）矩阵向量乘法可通过快速傅里叶变换高效实现，从而降低了相应测量和重建过程的时间复杂度。3.2.3 结构随机矩阵由于Gauss矩阵与Bernoulli矩阵随机性较强，确定性矩阵难以证明具有阶数较好的RIP性质。我们将介绍介于确定与随机矩阵之间的一种矩阵：结构随机矩阵。与确定性矩阵相比，结构随机矩阵多了些随机性，因而可以证明其具有较好的RIP性质，同时，结构随机矩阵的随机性较弱，一般仅具有行随机。更为重要的是，在很多实际应用中，观测矩阵为一结构随机矩阵。三类压缩感知观测矩阵效果对比如表31。表31 压缩感知观测矩阵效果对比压缩感知观测矩阵随机矩阵确定性矩阵结构随机矩阵RIP特性强弱中普适性强弱中存储空间很大大大减少中构造复杂度高低中3.3 重建算法要精确重建原始信号，首先要保证观测矩阵满足RIP条件，在此前提下，我们可以先通过求解范数下的最优化问题，得到原始信号在变换域的稀疏表示系数，得到稀疏系数之后，将其变换为原始信号，完成数据重建。对原始信号的重建过程可以视作是范数下的最优化问题：(38)由于观测值的维数远小于稀疏系数的维数，所以这是一个NP-Hard问题，想要求解该问题，则必须将稀疏系数中非零值的种可能一一列举，即使用穷举法求解，然而这种方式会耗费大量的时间，换言之，我们在实际应用中不能直接根据上式来求解44。但是可以证明，最小范数和最小范数能够在一定条件下互相转化，具有等价性45。那么式(38)可转化为最小范数下的最优化问题：(39)以下将介绍几种典型的重建算法：3.3.1 基追踪法由于最小范数问题为NP-hard问题，可将问题转换为范数问题46：(310)这个问题是一个凸优化问题，我们可以把这个凸优化问题转化为一个线性规划（LP）问题加以求解，称为基追踪（Basis Pursuit，BP）47方法。线性规划问题求解思路如下48：将式(38)进行如下变换：(311)其中为线性规划的初始参数，为线性规划的系统函数，为测量值。该算法将待恢复数据分为正负两部分进行重建，重建完成后将结果相加得到恢复的稀疏系数。BP算法是一种全局最优化算法，它能够很稳定的分解信号，需要的测度最少，但缺点是复杂性高，复杂性高意味着在进行重建时需要大的计算量从而会影响重建速度，这个缺点使得BP算法在很多场合并不适用。如果考虑误差，我们可以将上述问题进一步转化如下：(312)转化后的问题可以用二阶圆锥规划解决41。3.3.2 匹配追踪法匹配追踪（MP）算法属于贪婪算法，其基本思想是在每一次的迭代过程中，先从观测矩阵中选取与信号最匹配的列向量进行稀疏逼近，同时计算出信号余量，接着继续选取与余量最为匹配的列向量进行稀疏，亦计算余量。反复迭代之后，信号便可以由这些选中的列向量进行线性表示，如此便可得到重建结果。虽然MP算法可以重建出原始信号，但是由于这种挑选原子（观测矩阵列向量）的方式有其不足之处，信号在已选定原子集合上的投影并不具有正交性，这种非正交性带来的影响是无法保证每次迭代结果都最优，因而有时要经历多次迭代才能获得好的重建质量。该算法解决的是范数下的最优化问题：(313)若考虑重建误差，并将观测矩阵与稀疏化矩阵的乘积看作压缩感知的恢复算子（即）可将式(313)转化为如下的近似形式：(314)其基本思想是：采用贪婪迭代的方法选择恢复算子的列向量（原子），对所选原子的要求是：每次选中的原子都要与当前的冗余向量具有最大相关性，按此方式进行反复迭代处理，当迭代次数达到稀疏度K时，终止迭代。算法具体步骤如下：已知条件：采样向量，观测矩阵，稀疏化矩阵，恢复算子，稀疏度；初始化：残差，索引集，重建原子集合，；（1）计算残差和恢复算子列向量的内积，找出内积最大时对应的索引，即；（2）更新索引集，将找到的原子（恢复算子列向量）添加至重建原子集合； （3）更新残差，迭代次数累加；（4）判断迭代终止条件，如果，结束迭代；否则，继续回到执行步骤（1）；输出：信号的K-稀疏系数。3.3.3 正交匹配追踪法正交匹配追踪（Orthogonal Matching Pursuit，OMP）49算法是在MP算法的基础上发展起来的，该算法有效解决了MP算法中信号在选定原子集合上的投影非正交性的问题。OMP算法依然使用MP算法挑选原子的方式，即在每次迭代时得到一个原子，不过该算法还会对选定原子集合通过递归保证正交化，通过正交化可以保证单次迭代的最优性，这样可以减少迭代次数，且性能优于MP算法。由于具有正交性，可根据选出的原子从观测矩阵的向量中去掉相关部分，按此方式