一、什么是压缩感知(CS)?
compressed sensing又称compressed sampling,似乎后者看上去更加直观一些。没错,CS是一个针对信号采样的技术,它通过一些手段,实现了“压缩的采样”,准确说是在采样过程中完成了数据压缩的过程。
因此我们首先要从信号采样讲起:
1. 我们知道,将模拟信号转换为计算机能够处理的数字信号,必然要经过采样的过程。问题在于,应该用多大的采样频率,即采样点应该多密多疏,才能完整保留原始信号中的信息呢?
---------------------------------------
2. 奈奎斯特给出了答案——信号最高频率的两倍。一直以来,奈奎斯特采样定律被视为数字信号处理领域的金科玉律。
---------------------------------------
3. 至于为什么是两倍,学过信号处理的同学应该都知道,时域以τ为间隔进行采样,频域会以1/τ为周期发生周期延拓。那么如果采样频率低于两倍的信号最高频率,信号在频域频谱搬移后就会发生混叠。
4. 然而这看似不容置疑的定律却受到了几位大神的挑战。Candes最早意识到了突破的可能,并在不世出的数学天才陶哲轩以及Candes的老师Donoho的协助下,提出了压缩感知理论,该理论认为:如果信号是稀疏的,那么它可以由远低于采样定理要求的采样点重建恢复。
5. 而突破的关键就在于采样的方式。当我们说“采样频率”的时候,意味着做的是等间距采样,数字信号领域通常都是做等间距采样,也服从奈奎斯特采样定律。
但是如果是不等间距采样呢?依然必须要服从采样定理吗?
6. 答案是,随机的亚采样给了我们恢复原信号的可能。
上图非常关键,它可以简单直观地表述压缩感知的思路。 如图b、d为三个余弦函数信号叠加构成的信号,在频域的分布只有三条线(图a)。 如果对其进行8倍于全采样的等间距亚采样(图b下方的红点),则频域信号周期延拓后,就会发生混叠(图c),无法从结果中复原出原信号。
7. 而如果采用随机亚采样(图b上方的红点),那么这时候频域就不再是以固定周期进行延拓了,而是会产生大量不相关(incoherent)的干扰值。如图c,最大的几个峰值还依稀可见,只是一定程度上被干扰值覆盖。这些干扰值看上去非常像随机噪声,但实际上是由于三个原始信号的非零值发生能量泄露导致的(不同颜色的干扰值表示它们分别是由于对应颜色的原始信号的非零值泄露导致的)
P.S:为什么随机亚采样会有这样的效果?
这可以理解成随机采样使得频谱不再是整齐地搬移,而是一小部分一小部分胡乱地搬移,频率泄露均匀地分布在整个频域,因而泄漏值都比较小,从而有了恢复的可能。
8. 接下来的关键在于,信号该如何恢复? 下面讲一种典型的算法(匹配追踪):
(1) 由于原信号的频率非零值在亚采样后的频域中依然保留较大的值,其中较大的两个可以通过设置阈值,检测出来(图a)。
(2) 然后,假设信号只存在这两个非零值(图b),则可以计算出由这两个非零值引起的干扰(图c)。
(3) 用a减去c,即可得到仅由蓝色非零值和由它导致的干扰值(图d),再设置阈值即可检测出它,得到最终复原频域(图e)
(4) 如果原信号频域中有更多的非零值,则可通过迭代将其一一解出。
以上就是压缩感知理论的核心思想——以比奈奎斯特采样频率要求的采样密度更稀疏的密度对信号进行随机亚采样,由于频谱是均匀泄露的,而不是整体延拓的,因此可以通过特别的追踪方法将原信号恢复。
二、压缩感知的前提条件
接下来我们总结一下,能实现压缩感知的关键在于什么,即需要哪些前提条件。
9. 在刚才的讲述中大家可以感受到,这个例子之所以能够实现最终信号的恢复,是因为它满足了两个前提条件:
1. 这个信号在频域只有3个非零值,所以可以较轻松地恢复出它们。
2. 采用了随机亚采样机制,因而使频率泄露均匀地分布在整个频域。
这两点对应了CS的两个前提条件——稀疏性(sparsity)、不相关性(incoherence)。
10. 关于稀疏性可以这样简单直观地理解:若信号在某个域中只有少量非零值,那么它在该域稀疏,该域也被称为信号的稀疏域。
因此,第一个前提条件要求信号必须在某一个变换域具有稀疏性。比如例子中,信号在频域是稀疏的,因而可以通过所述的重建方法轻松地在稀疏域(频域)复原出原信号。
然而通常信号在变换域中不会呈现完全的稀疏性。其实只要它近似满足稀疏性,即大部分值趋于零,只有少量大的非零值,就可以认为它是可压缩信号,可以对它进行CS亚采样。
对于之前讲的例子,如果它在频域中不稀疏,我们可以做DWT、DCT等,找到它的稀疏变换。
---------------------------------------
11. 这里针对信号的稀疏性和信号压缩额外补充一下:其实,信号的稀疏性已经在图像压缩领域有了很广泛的应用。利用信号的稀疏性,可以对信号进行压缩。如图像压缩领域的JPEG格式,就是将图像变换到离散余弦域,得到近似稀疏矩阵,只保留较大的值,从而实现压缩。
---------------------------------------
12. 比如这个例子中,仅用原图像6.9%的点就复原了和原图像基本相同的图像。我们还可以采用小波变换,即为JPEG-2000,压缩效果更好。
---------------------------------------
13. 这里需要指出,图像压缩和压缩感知这两个概念很容易弄混,大家一定要分清。
它们其实有着本质上的区别。图像压缩是先进行了全采样,然后再变换域丢弃小系数,完成压缩;
而压缩感知不同,它的思想其实从图像压缩中借鉴了很多:既然全采样了还要再丢弃,我们为什么不能直接少采样一些点?因此,压缩感知直接进行了亚采样,然后再用算法消除亚采样导致的伪影。可以说,压缩感知直接在采样时就完成了压缩。
---------------------------------------
14. 接下来,在将第二个前提条件之前,还是需要引入必要的数学表达的。上图是一个大家在压缩感知相关的书籍文献中会经常看到的一张示意图。很多文章试图用这张图给大家讲清楚什么是压缩感知,结果导致大家看得一头雾水,混淆在各种“矩阵”当中。。不过相信有了我之前的讲解,现在这张图会好理解很多。这张图也就是把亚采样的过程用矩阵的方式表达出来而已:
如图,x是为长度N的一维信号,也就是原信号,稀疏度为k。此刻它是未知的。
Φ为观测矩阵,对应着亚采样这一过程。它将高维信号x投影到低维空间,是已知的。
y=Φx为长度M的一维测量值,也就是亚采样后的结果。显然它也是已知的。
因此,压缩感知问题就是在已知测量值y和测量矩阵Φ的基础上,求解欠定方程组y=Φx得到原信号x。
然而,一般的自然信号x本身并不是稀疏的,需要在某种稀疏基上进行稀疏表示。令x=Ψs,Ψ为稀疏基矩阵,s为稀疏系数。
于是最终方程就变成了:y=ΦΨs。已知y、Φ、Ψ,求解s。
15. 对应一开始的例子大家就能明白:x就是三个正弦信号叠加在一起的原信号;稀疏矩阵Ψ就是傅里叶变换,将信号变换到频域S;而观测矩阵Φ就对应了我们采用的随机亚采样方式;y就是最终的采样结果。
16. y=ΦΨs有点长,我们把ΦΨ合并成一个矩阵,称之为传感矩阵。即令Θ=ΦΨ ,则y=ΘS。
问题即为,已知y和Θ,求解S。
求解出S后,由x=Ψs即可得到恢复出的原信号x。
然而在正常情况下,方程的个数远小于未知数的个数,方程是没有确定解的,无法重构信号。但是,由于信号是K稀疏,如果上式中的Φ满足有限等距性质(RIP),则K个系数就能够从M个测量值准确重构(得到一个最优解)。
17.接下来的数学内容可以简短略过:陶大神和Candès大神证明了RIP才是观测矩阵要满足的准确要求。但是,要确认一个矩阵是否满足RIP非常复杂。于是Baraniuk证明:RIP的等价条件是观测矩阵和稀疏表示基不相关(incoherent)。
这就是压缩感知的第二个前提条件。
18. 那怎样找到不相关的观测矩阵呢?陶哲轩和Candès又证明: 独立同分布的高斯随机测量矩阵可以成为普适的压缩感知测量矩阵。
于是满足高斯分布的随机测量矩阵就成了CS最常用的观测矩阵。
对于二维信号,往往就采用如右上图所示的采样矩阵对图像进行亚采样。
对于一维信号,采用前文提到的随机不等间距的亚采样即可。
------------------------------------------------------------------------------
到这里,我们可以这样用一句话概括地描述什么是压缩感知:
如果一个信号在某个变换域是稀疏的,那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上,然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号。
以上可以算作是压缩感知的定义吧。但是如果要再简洁一点呢?
在我看来,压缩感知可以用这样一句话来表述:
直接采集出一个JPEG
——之前图像压缩的方法是全采样之后再压缩,抛弃稀疏变换域中的一些小系数;而CS直接减少了采样点,采集完后、经过重建的图像,就是一副在某变换域稀疏的压缩图像,比如JPEG。
那这么做有什么优势呢?
对于很多情形,比如照相机拍摄照片,这样减少采样点并没有优势。因为所有像素的采集在一瞬间就都完成了。
但是对于一些采集比较慢的情形,比如核磁共振成像,CS就可以发挥巨大优势。原本一副MRI图像常常需要几十秒,速度慢也是MRI的一大缺陷。而应用CS技术后,只需要采集全采样几分之一的数据,就可以重建出原图。这样就可以把成像速度提高好几倍,同时对图像质量影响不大。
另一个应用是Rice大学开发的单像素相机,也就是说这种相机只需要一个像素,非常有趣。感兴趣的朋友可以自己去调查。
三、压缩感知的重建方法
如前文所述,CS的重建也就是求解欠定方程组y=ΘS的方法。这是一个零范数(l0)最小化问题,是一个NP完全问题(没有快速解法的问题),因此往往转换成一范数(l1)最小化的求解,或者用一些近似估计的算法。这部分的具体内容在这里就不再详述了。