第7卷第1期 智能系统学报 Vol.7 No.1 2012年2月 CAAI Transactions on Intelligent Systems Feh.2012 D0I:10.3969/i.issn.16734785.201110010 网络出版t地址:htp://www.cnki.net/kcma/detail/23.1538.TP.20120216.2036.002.html 压缩感知理论及其在成像技术中的应用 赵春晖,刘巍 (哈尔滨工程大学信息与通信工程学院,黑龙江哈尔滨150001) 摘要:在传统的Shannon/Nyquist采样定理指导下,信号处理往往面临两大难题:A/D转换器技术的限制和海量采 样数据的处理压力.压缩感知(CS)理论表明当信号具有稀疏性或可压缩性时,可以通过全局非自适应线性投影的方 式,用远低于Shannon/Nyquist采样定理要求的颍率获取信号的全部信息.以CS理论为基础的压缩成像(CI)技术在 近年来得到了快速的发展.在概述CS理论的基础上,着重介绍了CI技术的原理及其发展状况,并从稀疏分解、观测 矩阵的构造和重建算法3个方面对其关键问题进行了分析.CI系统可以显著节省感光器件的数量,避免了传统成像 技术“先采样再压缩”方式带来的资源浪费. 关键词:压缩感知:压缩成像:稀疏分解:观测矩阵:重建算法 中图分类号:TP391.41;TN911.73文献标识码:A文章编号:1673-4785(2012)01-0021-12 Compressive sensing theory and its application in imaging technology ZHAO Chunhui,LIU Wei College of Information and Communication Engineering,Harbin Engineering University,Harbin 150001,China) Abstract:Under the guidance of the traditional Shannon/Nyquist sampling theorem,signal processing often faces two problems:technology limitation of the A/D converter and processing pressure caused by a mass of sampled da- ta.Compressive sensing (CS)theory suggests that when the signal is sparse or compressible,by means of global non-adaptive linear projection,all the signal information can be obtained with the samples much less than the sam- pling theorem required.CS theory based compressive imaging(CI)technology has been developed significantly in recent years.This paper first reviewed the principles of CS,and on this basis,discussed the theory and develop- ment of CI technology.The key issues of CI were also analyzed from three aspects of sparse decomposition,con- struction of measurement matrix,and the reconstruction algorithm.The CI system can significantly cut down on the number of photosensitive devices to avoid resource waste caused by a traditional"sample-then-compress"frame- work. Keywords:compressive sensing (CS);compressive imaging (CI);sparse decomposition;measurement matrix; reconstruction algorithm 作为信号处理领域的一块基石,传统的Shan 宽带通信等领域中,Shannon/Nyquist采样定理所要 non/Nyquist采样定理给出了一种将模拟信号离散 求的采样频率超出了一般A/D转换器的处理能力, 化的局部采样方案并确定了采样间距的下限,它作 导致信号包含的信息难以在采样过程中得到完整的 为信号处理领域中的一个核心原理指导着现今几乎 保存并且极易引入混叠和各种噪声,其导致的另一 所有的信息采集工作;但是在诸如雷达、医学成像、 个结果就是接收端被淹没在海量数据之中.而为了 存储和传输这些海量的采样数据,又必须额外引入 收稿日期:2011-10-28.网络出版时间:201202-16. 基金项目:国家自然科学基金资助项目(61077079);高等学校博士学 复杂的压缩算法.随着社会的进步,人们对于信息的 科点专项基金资助项目(20102304110013);哈尔滨市优秀 需求正在成倍的增长,传统的采样方式已经成为了 学科带头人基金资助项目(2009 RFXXG034). 通信作者:赵春晖.E-mail:zhaochunhui(@hrbeu.edu.cn. 当前信息技术发展的主要瓶颈之一,迫使人们开始
·22 智能系统学报 第7卷 在信息获取和处理等方面寻找新的出路 式中:更是一个以(*表示共轭转置)为行的 近些年来,一种名为“压缩感知”(compressive M×W矩阵,表示一个普遍意义下的采样行为,不妨 sensing,CS)[1s]的理论开始频繁出现在人们的视线之 称其为观测矩阵,这样可以简单地将信号与波形 中,该理论证明了当信号具有稀疏性或可压缩性时, P()联系起来:如果P(t)是冲击函数,那么y就 可以通过全局非自适应线性投影的方式,用远低于 是理想的信号时域抽样序列;如果9:()是正弦波, Shannon/Nyquist采样定理要求的频率获取信号的全 那么y就是信号的傅里叶系数,这也正是核磁共振 部信息,并且采用非线性优化方法重建该信号.2004 成像(magnetic resonance imaging,MRI)等医学类断 年加州理工大学的E.Candes等在研究医学图像重建 层成像技术的信号获取模型。利用这种采样机制可 的过程中发现了这一现象4]:2006年D.Donoho所撰 以通过M个向量P.(t)∈R"来获取信号f∈R“的 写的一篇以“压缩感知”为名的论文发表在EEE 信息.令人真正感兴趣的是M≤N的情况,也就是说 Transactions on Information Theory上,该理论因此得 当观测值∫的长度远远小于信号值∫的长度时是否 名;随后E.Candes等又进一步证明了采用非自适 有可能从少量的测量值中精确恢复信号呢?是否可 应的随机线性投影方法,可以在观测值高度不完备和 以人为地设计M(M≤N)个P:(t)使其获得信号的 不精确的情况下以极大的概率实现信号重建6).在 全部信息呢? 此基础上D.Donoho等又在多篇文章中对这一技术进 从数学的角度来看上述想法几乎是不可能的, 行了深人的分析与进一步的探索[8],使得该理论的 要解决这样的问题就涉及到求解一个条件远远少于 框架日趋完善.作为一种全局采样方式,CS采样过程 未知数的“欠定”方程组,而这种方程组的解应当有 需要同时操作信号的全部采样值.由于潜在像素时间 无穷多个,但是如果信号具有某种特殊的形式,比如 上的同时性,以及自然图像的成像恰好符合这一要 一个信号的值在大部分时间里都是0,只有S个位 求,因此以CS理论为基础的压缩成像(compressive 置有值,那么只要有M≥S个观测样本,方程组中的 imaging,CI)技术241可以显著节省传感器的数量, 方程数量相对于未知数来讲就显得“充足”了,这就 避免了“先采样再压缩”带来的资源浪费,该项技术的 涉及到了信号中普遍存在的稀疏性和可压缩性, 发展将会对诸多领域产生革命性的影响.例如在非可 1,2信号的稀疏性和可压缩性 见光谱成像中,感光器件的价格往往是极其昂贵的, 如果将信号∫的支撑集的势定义为sup(f):= 若要使用传统方法在某一波段对特定的场景成像,那 {i:f≠0},那么包含在集合s:={fup()≤S}中的 么一次就要根据成像精度的要求使用几百万或上干 信号就称为S稀疏信号,换句话来说S稀疏信号就是 万个针对该频谱的感光器件,但是如果调整成像波 指信号中非零值元素个数最多不超过S个的信号. 段,那么就要将这几百万个感光器件全部更换,假如 信号的稀疏性并不仅仅体现在时域上,考虑一 能够通过CS理论将感光器件的数量降低到一个,那 个长度为N的实值一维离散时间信号f∈R“(高维 么将大大节省成本 信号可以向量化为一维信号).对于一组N×1维正 CS理论的魅力在于其暗示了一种全新的信号 交基向量{中:1,可以将信号表示为这组正交基的 获取准则,引领人们从另外一个角度审视信息与载 线性组合: 体信号之间的深层联系,在此基础上发展起来的CI 技术在降低成像成本、优化成像方式等方面拥有巨 (1) 大的发展潜力,在医学成像、空间探测、非可见光成 式中:x为正交变换的系数.如果以少作为列,构成 像等领域具有广阔的应用前景.本文在概述CS理 一个N×N维的正交基矩阵平的话,式(1)可以方 论的基础上,着重介绍了CI技术的原理及其发展状 便地表示为 况,并对其关键问题进行了分析。 f=Vx. (2) 1 式中:x是以x:为元素的N×1系数向量 压缩感知的基本原理 显然时域形式∫或频域形式x在表示同一个信 1.1常规采样机制 号时是等价的.从这个角度上来说,如果系数矩阵x 为了便于后续的讨论,首先要介绍一下信号处 的支撑集的势为sup(x):={i:x:≠0},那么当 理中常规的采样机制.一个信号f(t)的获取可以表 sup(x)≤S时,依然可以将信号f称为S-稀疏信号, 示为 只是该信号的稀疏性表现在了频域上. y=(f,p〉,k=1,2,…,M 在实际的处理过程中很难遇到严格的稀疏信 或者等价的 号,待处理信号的能量不会仅仅集中在时域或频域 y =f. 中的若干个点上,大多数情况下它们是散布于整个
第1期 赵春晖,等:压缩感知理论及其在成像技术中的应用 ·23· 信号持续的范围之内的.一般来讲信号的频域系数 于信号∫,即整个观测过程不必是自适应的,它是独立 会遵循一种指数规律516,例如对于一个信号∫ 于信号的形式而存在的;另外业可以认为是由任意一 R“,用一组正交基{少:将其分解为f=∑1x:, 组将信号∫分解为稀疏形式的正交基(或紧框架)组成 如果将系数:按照从大到小的递减顺序重新排列 的矩阵,它只与信号的重建有关 的话,即x)≥xl2≥…≥xlwM,那么对于每个1≤ 在CS理论中,观测矩阵的构造问题实际上就 n≤V来说,各系数应该满足6] 等同于编码模式的选择问题,规测系统的设计与实 xl(am≤R·np 现都是围绕着观测矩阵展开的.观测矩阵在对信号 式中:0≤p≤∞.在这个前提下可以知道,一般信号 进行采样观测的同时需要保持原信号的结构,因为 的系数x中只有少数几个较大的系数,其他系数相 一旦观测过程破坏了信号中的信息,那么精确重建 对较小,称这样的信号为“可压缩信号”,信号的“可 就是不可能的了].观测矩阵的构造问题同时涉及 压缩性”是普遍存在于自然界中的,在此基础上可 到观测矩阵和压缩感知矩阵,因为对于一个严格稀 以在不造成太大感官损失的情况下抛弃小系数,仅 疏的信号来说,观测矩阵Φ就是压缩感知矩阵Q. 仅保留较大的系数。 在CS过程中,压缩感知矩阵⊙必须保证长度为N 1.3压缩感知的编解码原理 的x可以通过M(MN)个观测值y恢复出来.由 1.3.1编码观测原理 于M<N,这个问题似乎是病态的.然而如果x是S 对于一个长度为N的信号f∈R",假设其在时 稀疏的,并且x中的S个非零系数的位置已知,那么 域不具有稀疏性,考虑一个线性观测过程y=,其 只要令M≥S就可以解决这个问题.将这个问题简 中更为一个M×N矩阵.假定MgN,得到一个M× 化成为良态问题的一个充要条件是,对任意x与具 1向量y.显然观测向量y就是信号∫在中上的一个 有相同非零元素的向量y来说,压缩感知矩阵⊙需 线性投影,如果∫对于一组正交基是稀疏的(或可压 要满足一定的原则. 缩的),那么根据式(2),可以用变换域系数x的S 在文献[8]中E.Candes等介绍了一种“一致不 稀疏形式(或近似稀疏形式)代替其时域形式∫,则y 确定原则”(uniform uncertainty principle,UUP),即 可以表示成 压缩感知矩阵⊙必须满足 y=时=Φ%=⊙x, (3) 式中:压缩感知矩阵⊙=ΦD亚是一个M×N矩阵,观 %11≤161≤1I 测向量y由⑧中对应于信号稀疏值位置的S列线 UUP在本质上阐明了®具有一种“受限等距特 性叠加而成.以上构成了CS的基本编码步骤, 性”(restricted isometry property,RP).之后在文献 如图1所示. [11]中这一原则又被进一步修改为如式(4)形式, 若从⑨中提取T列(TC{1,2,…,N},IT1≤S)构成 子矩阵⊙,则⑨,满足 1-d≤I8,≤1+d 2 (4) 式中:δ为⊙的S-受限等距常数.这一原则本质上 要求压缩感知矩阵⑨表现为近似正交的(由于M (a)时域形式 N,所以不可能完全正交).式(4)所表示的就是IP 的具体形式,IP从本质上说明了一个问题:要想完 全重建原始信号,必须确保观测矩阵不会把2个不 NXI 同的S-稀疏信号映射到同一个观测集合当中,这就 要求从观测矩阵中抽取的每T个列向量构成的矩 阵都是非奇异的. 尽管RP理论特性完美,但是它只是观测矩阵 (b)变换域形式 的一个充分条件,而不是必要条件.换句话说,满足 图1压缩感知编码原理 P条件的矩阵可以作为观测矩阵,但是可以作为 Fig.1 Encoding schematic diagram of compressive 观测矩阵的矩阵未必都满足RP条件.另外RP条 sensing 件不具有可计算矩阵结构的简单易行的公式,这样 这里称M×N矩阵Φ为观测矩阵,N×N矩阵亚 一来很难用它来判断某一观测矩阵是否拥有这种特 为稀疏分解矩阵此处需要强调的是Φ的形式不依赖 性,并且也很难用它来指导观测矩阵的设计.R.G
24 智能系统学报 第7卷 Baraniuk给出的RP的等价条件是观测矩阵Φ和 的一个高度非线性空间.解集空间H=N(⊙)+x由 稀疏矩阵业不相关,其具体内容就是Φ的行不能 于矩阵®的随机性而具有一个随机的角度,如 由业的列来表示,反之亦然。 图2(a)所示(在实际使用中M,S≥3,此处为了直观 1.3.2解码重建原理 选择三维),整个优化过程可以想象成是一个“弱球 考虑S-稀疏系数向量x,如果⑨x=y,那么对于 体”从原点开始慢慢膨胀,这个“球”首次与解集空 任意位于⊙的空集空间N(⊙)中的向量r(⑨=0) 间H接触的点即为最终的解,由于所期望得到的解 来说,都有O(x+r)=y成立.这样在式(3)中,当 是稀疏形式的,所以“球”和H首次接触的点应该位 M<N时应该有无穷多个解x'满足⊙x'=y.因此, 于坐标轴(或若干坐标轴组成的超平面)上.2-最小 重建算法的目标就是在解集空间H=N(⊙)+x中 化得出的解x是在H上最接近原点的点,这个点可 找到信号的稀疏系数向量x,欲在众多候选向量之 以通过一个与H相切的超球面(1,2-球)来找到.考虑 中找到最符合稀疏条件的解,自然要求助于最优化 到H方向的随机性,最近点x有很高的可能性出现 算法.若将向量x的,范数定义为 在远离坐标轴的位置,因此可能既不稀疏又不接近 正确答案x.与此相对地,图2(b)中的L-球具有坐 (IxI)P=∑I:IP (5) 标轴上的点:因此,当构建球从原点开始膨胀的 这里有3种范数优化方式可供选择: 时候,它会首先与解集空间丑相交于一个靠近坐标 1)最小化2范数重建.一般来讲,解决这类优 轴的点,这个点正是稀疏向量x的位置. 化问题的经典方法是通过求解优化方程: x=min‖x'‖2满足r'=y, 找到变换空集空间中,2范数(能量)最小的向量x. 这种最优化方法有方便的近似解x=Q(⑨0)-y, 遗憾的是,范数最小化方法几乎难以找到一个真 正S稀疏的解,取而代之的是含有大量非0值的x. 这是因为从能量的角度上讲,2范数最小化是在众 多信号可能的形式中寻找使其能量最小的解,而将 能量散布于整个信号范围内往往是将能量降低的有 效形式,显然这不是希望见到的形式 2)最小化范数重建.由于2范数最小化表示 (a)l,2范数最小化视觉效果 的信号能量不具有稀疏性,所以考虑针对x中非零 值计算的,范数(S-稀疏向量的,范数等于S).将 优化方程修改为 x=min‖x'‖o满足r'=y, (6) 求解这一优化方程可以精确恢复一个S稀疏信号. 但一般来讲,求解式(6)既是数值不稳定的同时也 是NP-hard的,即无法设计有效的算法进行计算,只 能通过穷举x中非零值位置的全部可能性来一一验 证,直到找到满足条件的解为止 3)最小化1范数重建.相比以上2种方法,求 解基于最小化11范数的优化方程: (b),范数最小化视觉效果 x=min‖x'‖1满足⊙r'=y, (7) 图2重建算法的几何解释 可以精确恢复S稀疏信号和近似可压缩信号.这是 Fig.2 Visualization of reconstruction algorithm 一个凸优化问题1,可以简化为计算复杂度约为 1.4模拟信号的压缩感知 O(N3)的基追踪(basis pursuit,BP)问题9T, CS是从离散信号处理中发展起来的,当然关心它 文献[2-3]分别在二维和三维空间中对最小化 是否可以推广到模拟信号处理当中.在实时处理摸拟 1:范数重建信号的有效性做出了几何解释.通过图 信号时,通常采用的模数(A/D)转换技术往往受到 2可以直观地了解为什么2-最小化所无法求得的稀 Shannon/Nyquist采样率的限制,为了摆脱这种束缚人 疏解,1-最小化却可以找到.在R中所有的S-稀疏 们开始转向更进一步的深层次转换—模拟/信息 向量x的集合是由坐标系中全部S-维超平面组成 (analog-to-information,A/I)转换n2].A/I转换针对模
第1期 赵春晖,等:压缩感知理论及其在成像技术中的应用 ·25· 拟信号本身处理,不必将其转换为数字形式就可以通 下限首先要引入矩阵相关性的概念.观测矩阵Φ和 过一系列的步骤得到离散形式的观测值, 稀疏分解矩阵业之间的相关性可以采用矩阵间的 对于模拟信号f(t),利用一组连续的正交基 相关系数(Φ,)表示, {中.(t)}=可以将其分解为 u(D,)=Nmw(,,〉l. ft)=∑,.(t). 式中:(中,业)∈[1,√n].简单来说,相关系数表征 这里假设系数向量x是严格S稀疏的.尽管每一个 了2个矩阵所包含的元素之间最大的相关关系。 业(t)都可能具有极宽的带宽,但是信号f八t)却可以 (中,)越大表示2个矩阵之间的相关性越大,反 只用很少的自由度确定, 之则越小 A/I转换的处理过程大致分为3步:调制、滤波 定理1】假设信号fR"在稀疏分解矩阵 和采样21,如图3所示.首先利用元素为±1的伪 亚下其系数x是S稀疏的,在Φ域中随机选择M 随机序列P。(t)对模拟信号f(t)进行调制(即相 个规测值满足 乘),这样做的目的是使信号的频域成分被分散开, M≥C×u2(Φ,)×S×logN, 不受后续滤波处理的破坏,此处的调制速率必须大 式中:C为一个正值常数那么式(7)表示的1-最小 于或者等于输入信号f(t)的Shannon/Nyquist采样 化将以极大的概率取得精确解 率;之后用冲击响应为(t)的低通滤波器对其进行 对定理1有如下4点说明: 滤波;最后使用传统的模数(A/D)转换器以频率 1)相关性的作用是十分明显的,相关性越小则 采样.尽管输入信号是模拟量,但该系统得到的观测 需要的采样数就越少, 值y依然是离散形式的 2)每一个观测值中都包含着原始信号的一部 分信息,可以通过任意M个观测值的组合获得信号 y[m]= v.(r)p.(r)h(mg-)dr. 的全部信息 (8) 3)可以通过解一个凸优化问题,从任意M个观 测值中重建原始信号∫,这个重建过程是独立于信 榄拟 号的,不需要知道任何与原始信号有关的先验知识, 滤波器) 只要原始信号足够稀疏就可以得到精确的结果。 量化器 p0{-1,lH 4)“以极大的概率取得精确解”暗示了例外情 况存在的可能性,某些经过人为精巧安排的不遵守 图3A/I转换 此定理的信号是可能存在的,但是它们几乎不可能 Fig.3 Schematic diagram of A/I conversion 在自然情况下出现 为了与离散形式的CS统一起来,依然希望y是 2.2对可压缩信号的重建 由f(t)通过一个M×N压缩感知矩阵⊙来得到.根 在前面介绍的CS技术都是针对严格稀疏的信 据式(8)可以如此构造压缩感知矩阵,令其第m行、 号的,但是实际应用中遇到的信号往往不是严格稀 第n列的元素0m,n为 疏的,而是可压缩的.那么自然要问这样的问题:CS 0-v.(r)p.(r)h(mg-r)dr. 对于可压缩信号的重建性能如何?显然精确重建是 显然压缩感知矩阵⊙仍然由两部分组成:稀疏分解 不可能的;但是考虑到可压缩信号能够由严格稀疏 矩阵亚将模拟信号分解为频域系数;观测矩阵Φ 信号近似,那么重建信号也应该是精确解的近似.假 从模拟信号中提取离散的观测值.如此一来,模拟信 设对于任意向量x∈R“,保留其中最大的S个元素 号的CS就被纳入到离散CS的轨道中来了,离散形 并将其余置0得到xs,下面将以xs为标准来衡量可 式中的理论可以简单地推广到模拟信号的处理当 压缩信号的重建效果, 中,以实现对更加复杂的模拟信号进行灵活的处 定理20假设6s<2-1,那式(7)的解x满足: 理23241 Ix-xI2≤Cx-‖山 2压缩感知的性能 √3 lx-x‖1≤C川x-x3‖… 2.1压缩能力 观测值向量y作为最终的编码结果,其长度与 式中:C为常数.如果x是可压缩的并满足指数规律 观测矩阵Φ的某一维长度M相关,M的大小显然 xl(m≤R×n,即x属于一个半径为R的弱l。球 不可能是任意的,它必然有一个下限为了找到这个 时,根据经典计算可知: