基于Hadamard 矩阵优化排序的快速单像素成像*

李明飞 阎璐 杨然 刘院省

1) (中国航天科技集团有限公司量子工程研究中心,北京 100094)

2) (北京航天控制仪器研究所,北京 100039)

(2018 年10 月22 日收到; 2019 年1 月2 日收到修改稿)

为提升单像素成像速度,提出了基于Hadamard 矩阵优化排序的压缩采样解决方案. 利用数值仿真和室外实验对提出的5 种排序方法进行了对比分析. 研究结果表明: 按Haar 小波变换系数绝对值排序时单像素成像效果最优,排序对应到Walsh 序后可利用快速变换重建图像,速度达300 帧/秒@64 × 64 像素; 最优排序下,采样率25%仍可重建图像,采样速度可提升4 倍. 针对排序方法与成像信噪比关系,从关联成像角度给出了其物理解释: 测量基矩阵元邻域数值相等的区域面积等效于光场二阶相干面积,当光场二阶相干面积随测量基由大到小排序时成像效果最优. 本文研究成果可用于提升单像素成像速度,具有实用价值.

1 引 言

单像素成像,区别于传统阵列探测器成像,其图像信息的获取方式与采集效率不同. 在硬件的复杂度、工业成本方面,单像素成像占有优势,这对未来成像系统的简化和集成具有潜在应用前景,特别是在红外、太赫兹波段,阵列传感器件不成熟,在三维成像、光谱成像对分辨和灵敏度要求较高的情况下,单像素成像技术的出现为这些问题带来了新的解决方案. 单像素成像技术,最初采用的算法是压缩感知算法. 压缩感知理论由Candès 等[1]于2006 年正式提出,2008 年美国Rice 大学在实验室实现了单像素成像,验证了压缩感知理论在单像素成像方面的应用潜力[2]. 压缩感知理论的核心思想是: 对稀疏或可压缩的信号进行少量非适应性的线性测量,可以接近百分之百地恢复出原始信号,一定意义上突破了Shannon 采样定理的限制,还具有超灵敏感知等优势[3]. 目前用于单像素成像的压缩感知算法,测量次数仅为奈奎斯特采样极限的三分之一左右即可重建出目标图像[4],甚至某些条件下1%的奈奎斯特极限采样,就可重建包括目标绝大多数信息的图像[5,6]. 因此,压缩感知算法在单像素成像、计算层析成像等研究领域得到广泛研究和应用. 压缩感知算法特点是: 在运算中迭代求解一类凸优化问题; 重建效果与所选稀疏基密切相关.然而,算法需要迭代和对稀疏基的有限等距性质要求,使得计算时间和成像质量的稳定性受到极大挑战. 单像素成像中的另一常用算法是关联算法和矩阵求逆算法,其具有硬件要求低、成像速度快的特点,但其需要大于或等于奈奎斯特采样极限的测量次数,采样时间也较长. 如基于Walsh-Hadamard变换[7,8]、快速傅里叶变换的方法[9],提供了实现稳定、快速的成像方案. 特别是Walsh-Hadamard 算法,具有多种优势: 算法效率高,可并行和硬件加速; 调制矩阵为二值,易于在高速空间光调制器件上实现; 调制矩阵生成速度快且无需存储; Hadamard 矩阵作为测量矩阵,其均值不变且具有正交特性,便于利用差分测量降低成像误差等.

尽管存在上述优点,但基于Walsh-Hadamard变换方法要求测量次数与被测图像的像素数相等才能完美恢复图像,无法实现压缩采样. 例如,要得到 1 28×128 像素的图像,测量基和测量次数为16384 次. 随图像增大,测量次数成二次方指数增加,相应地测量时间随之增长. 虽然在结合压缩感知算法后,随机抽取Hadamard 矩阵可实现压缩测量,但算法运行时间仍在秒的量级. 实际上,高速单像素成像核心的问题仍是测量基选取问题. 学者们已就不同测量基选取及其对单像素成像质量的影响开展了大量研究[1−10],但针对同一组测量基的不同排序对成像结果影响的研究较少[11−15].

针对上述问题,本文研究了基于Hadamard 矩阵测量基的按5 种不同表象变换的稀疏度值排序并找到了解决方案,为单像素成像测量基的选取提供了一种压缩采样理论和技术路线. 本文提供的方法具有三方面特点: 第一,可解决在正交测量基选取以及按多尺度分辨率排序中遇到的一类问题,为更高效的单像素成像技术中的测量基选取提供理论和技术支撑; 第二,本文排序结果与文献[12]提供的方法类似,但本文方法的实现手段更加普适和简单,且提供了排序的理论依据; 第三,可实现文献[11—15]中的压缩测量功能,且重建图像时可采用快速Walsh-Hadamard 变换,无需存储测量矩阵,在工程应用方面更加实用化.

本文接下来的内容主要为: 第2 部分介绍Hadamard 矩阵排序理论及方法; 第3 部分介绍数值仿真结果; 第4 部分为室外实验结果及分析; 第5部分为全文总结.

2 单像素成像理论

2.1 实验系统

单像素成像的基本原理示意图如图1 所示,待成像物体通过主透镜成像到数字微镜器件(digital micro-mirror device,DMD),微镜受数字电路控制系统的控制,按加载的调制矩阵Φ,每个微镜发生+12◦或−12◦摆动,摆动方向由矩阵Φ的矩阵元决定,例如矩阵元“1”对应 + 12◦方向翻转,矩阵元“–1”对应−12◦方向翻转. 物体的像被测量矩阵调制经由中继透镜汇聚到光电探测器转为电信号,重复刷新DMD 调制矩阵将得到因光强变化引起的一系列电信号,该过程等价于调制图案与物体实像作内积的过程,模拟的电信号经A/D 进行模数转换变为数字信号,数字信号将被采集到计算设备经过算法计算重建图像.

2.2 压缩感知

压缩感知理论,是实现单像素成像的常用算法. 基于压缩感知的单像素成像理论可用简洁的矩阵求解来描述. 设待测图像的向量形式为TMN×1,其中M,N为待测图像行与列. 当观测次数为K时,设观测值为一维向量YK×1,测量矩阵为ΦK×MN,则有下列关系YK×1=ΦK×MNTMN×1+Noise,TMN×1=Ψt′MN×1,求解方法如下:

图1 单像素成像原理示意图,其中待成像物体经主透镜成像到DMD 上,经编码矩阵调制的光束被反射,反射光由中继透镜汇聚到光电探测器,探测器将光信号转为电信号,电信号经过A/D 转换由模拟信号转为数字信号,重复刷新DMD 调制矩阵将得到一系列数字信号,将数字信号采集到计算设备即可利用算法重构图像Fig. 1. Schematic diagram of the experimental setup. Lena is the object to be recovered; the main lens imaging Lena onto the DMD,which is modulated by matrices,then the light is reflected and collected by the relay lens into the photo-detector. As DMD modulation matrices refreshing continuously,the analog-to-digital converter (A/D) connected with the photo-detector receives a series of digital signals,which are finally sent to recover the image by the computer.

式中,ΦK×MN为测量矩阵;Noise为环境和探测器电信号随机噪声项;Ψ=[Ψ1,Ψ1,Ψ1,···,ΨMN]为特定稀疏基(如Noiselet,Dct 等);M,N为待恢复图像行和列的像素数;K为总的测量次数;||∗||1和||∗||2分别表示ℓ1范数和ℓ2范数;ε为噪声振幅的边界. 可见,压缩感知特点是利用图像的自然稀疏特征,稀疏基Ψ将图像变换到稀疏空间,即少数值非零,然后求解方程(1)得到最优解,所需测量次数K<

2.3 单像素Walsh-Hadamard 变换成像

单像素Walsh-Hadamard 变换成像,是一种像素复用的成像方式. Hadamard 矩阵起源较早[16],应用极其广泛. Hadamard 矩阵(简称H矩阵)具有如下特征: 首先,H矩阵是只有±1 元素构成的方阵,产生速度快,可用直积方式产生. 其次,矩阵任意两行(或列)正交; 第三,H矩阵的逆矩阵是其自身,因此可构成正交归一的完备基:

其中,K为Hadamard 矩阵元个数,I为单位矩阵.Hadamard 变换不能直接利用快算法,需按Walsh-Hadamard 序实现快速变换[17],实现手段与快速傅里叶变换类似. 在单像素成像中,也可以利用逆矩阵计算图像T,对于Hadamard 矩阵而言,二者等价.

理论上,列相机与单像素相机的探测器灵敏度相同,相比于传统成像,单像素Walsh-Hadamard 变换成像可以在信噪比(signal-to-noise ratio,SNR)方面提升单像素Walsh-Hadamard 变换成像需对目标进行全采样,按Walsh 序或Hadamard序均不能在压缩采样下快速成像,本文方法将解决这一问题.

2.4 测量基矩阵排序理论

对单像素成像的传统研究方法中,主要研究在压缩感知理论框架下,如何根据图像稀疏性来构造出一组测量基,使得图像重建效果和效率最优,即高SNR 和高压缩比. 文献[12]提出了一种Hadamard 测量基排序方法,并将该方法与两种压缩感知算法重建结果进行了对比. 与其研究方法和思路不同,本文重点研究同一组测量基的不同观测顺序对成像结果的影响,在只采用快速变换的算法条件下如何实现高SNR 和快速成像是本文的研究目标,目前还未见有文章进行过相同研究. Hadamard矩阵排序基本思路是,保留Hadamard 矩阵的正交优势和Walsh-Hadamard 可快速变换的优势,但利用不同的表象变换,将测量基作多尺度近似,然后按分辨率的尺度,即按测量基矩阵从低频到高频顺序排序. 二维测量基矩阵在变换域的空间频率,相当于其在相应表象下的稀疏性,因此按测量基疏稀度进行排序具有合理性. 定义第i个测量基矩阵的稀疏度S(i) 为

其中,Ψ为表象变换矩阵,为抽取Hadamard矩阵的第i行(或列)构成的测量基矩阵. 本文接下来主要讨论表象变换Ψ为不同变换方法时,按稀疏度S(i) 值由小到大排序,测量基排序后压缩采样次数m与单像素成像质量的关系. 在仿真实验时,由于有原始图像作为对比,评价标准选为峰值信噪比(peak signal to noise ratio,PSNR)和结构相似性(structural similarity,SSIM),从而判断最优的采样次序. 在理论研究时选取文献[18]采用的图像SNR 作为评价依据.

2.5 测量基排序对成像信噪比影响

针对测量基不同排序,为给出相应重建图像SNR,需借鉴传统鬼成像(关联成像)中赝热光成像理论. 若将Hadamard 矩阵产生的测量基矩阵元理解为光场空间强度值,则每个最小的调制单元,即测量基矩阵元邻域数值相同的区域面积等效为光场的相干面积Acoh,整个调制矩阵大小M×N相当于光场的照射面积Abeam. 在本文中M=N=64 ,最小调制单元对应一个像素,即最小相干面积Acoh=1 ,在无压缩采样时,测量次数m=4096 . 根据Ferri 的理论,文献[18]中鬼成像SNR 公式为

其中,Ncoh=Abeam/Acoh,T为待测物体.Acoh的值在光学上通常采用HBT 算法来确定,即采用自相关算法求得各测量基矩阵元任意点自相关峰值半高宽. 由于(5)式中,当待测物体确定后,Δ和为固定值,因此可将上述模型简化,

根据(6)式可得不同排序方法SNR 曲线,结果如图2 所示. 研究发现Walsh 序和Haar 小波排序SNR 曲线重合,与理论常识不符. 经仔细分析发现: Ferri 提出的SNR 理论默认散斑为圆形,符合其实验装置,但实际应用中散斑形状由激光光斑决定,例如半导体激光产生长方形光斑,对应散斑是矩形,因此只利用圆形相干面积无法准确描述最终成像的SNR 和分辨率,然而这一点在以往文献中均按圆斑处理,未进行充分研究. 而光场横向相干面积的不对称,将对成像分辨率和SNR 造成较大影响. 本实验中调制图案的平均相干面积可以是正方形或长方形,将造成实际成像过程中图像的水平和竖直两个方向上分辨率的差异,这在实际成像应用中对分辨率影响非常关键. 针对这一现象,本文处理方法是将二阶自相关峰值半高宽在水平和竖直方向分别计算,然后以二者平方的平均值作为最终相干面积. 经上述处理后,按(6)式,该理论结果显示出不同排序重建图像SNR 的差异性,注意此时的理论曲线与测量图像内容无关,直接对应测量基不同排序对重建图像SNR 影响的规律. 至此,即解释了不同排序方法导致压缩采样后重建图像SNR 变化的原因: 测量基图像排序导致测量基二阶关联时相干面积发生变化,只有随测量次数增加时,相干面积从大到小变化,才使得图像SNR 快速上升. 而按Haar 小波变换系统绝对值排序后的测量基,完全按小波尺度函数从大的相干面积到小的相干面积排序,理论上将得到在相同压缩比条件下的最佳效果. 实际上,本文的后续仿真与实验研究均证实了这一点.

图2 不同排序方法下图像SNR 与测量次数m的关系,其中排序方法包括Haar 小波变换序、Db2 小波变换序、Dct 序、Walsh 排序和随机排序Fig. 2. Image SNR versus mth measurements under different ordering methods: Haar wavelet order,Db2 wavelet order and Dct order,Walsh order and random order.

3 数值仿真

3.1 压缩采样成像结果

为验证排序效果,利用数字图像作为测试图像进行了5 种排序下压缩成像数值仿真. 测试图像采用灰度较为丰富的Lena 图像,图像大小为64×64像素,如图3(a)所示. 实验分为三个步骤: 首先,生 成 4 096×4096 大 小 的Hadamard 矩 阵; 其 次,将Hadamard 矩阵的每一行重新按顺序排列成4096 个 6 4×64 大 小 的 矩 阵; 第 三 步,对4096 个64×64大小的矩阵按2.4 节所述的5 种规则排序:Haar 小波变换、Db2 小波变换和Dct 的S(i) 值由小到大排序,以及按Walsh 序排序和按1 到4096 随机(random)序排序. 图3(b)为Db2 小波变换排序测量基25%采样下重建图像; 图3(c)为Dct 排序测量基25%采样下重建图像; 图3(d)为Walsh 排序测量基25%采样下重建图像; 图3(e)为随机排序测量基25%采样下重建图像; 图3(f)为Haar 小波变换排序测量基25%采样下重建图像. 对图3(b)—(f)从视觉上定性判断,可明显看出Haar 小波变换排序重建图像质量显著优于其他方法. 上述算法运行在Intel 酷睿i3-3240 双核3.40 GHz 中央处理器,2GB 内存的联想台式机上,软件平台为Matlab 2012a,测试图像重建算法平均计算时间 3.0 ms±0.5 ms ,该计算速度可达到300 帧/秒的单像素成像速度.

3.2 压缩采样与重建图像质量关系

为定量分析测量基排序方法下压缩采样次数m与单像素成像质量的关系,评价标准选用PSNR 和SSIM 共同进行像质评价,两种不同SNR 评价方法的引入可一定程度上确保评价的有效性. 测试图像选用图3(a),对其进行50%降采样,图像大小变为 32×32 像素,对应全测量数为1024 个数据点,测量次数从m=128到m=1024 ,间隔128 个数据点. 实验中,Hadamard 测量基的排序规则如下: 1)按3 种表象变换Ψ,分别为Haar 小波变换(或Db1)、Db2 小波变换和Dct,测量基按稀疏度S(i) 由小到大排序; 2)按2 种已知的Hadamard 测量基排序,即Walsh 顺序排序和随机排序,其中随机排序采用计算机产生的随机顺序作为排序依据. 上述 5 种排序方法下,图像的PSNR 与SSIM 值与测量次数m的关系如图4(a)和图4(b)所示. 图4 的数值仿真实验结果表明,对于灰度图像,不同排序方法得出的单像素压缩采样成像的重建图像质量不同,PSNR 变化趋势与SSIM 变化趋势具有一致性; 当测量次数增加到全测量时,不同方法的重建图像PSNR 和SSIM 均趋于相同值,与理论相符; 在压缩采样条件下,Haar 小波变换排序对应的重建图像PSNR值和SSIM 值均优于其他方法,这与上一节25%采样下图像重建的视觉判断结果一致. 综合本节的数值仿真实验结果可知,对于灰度图像,按Haar小波变换排序的Hadamard 测量基,在相同压缩采样下可获得较高的重建图像质量,且压缩采样后数据可进行重排序和快速变换,计算速度可满足 300 帧/秒的单像素成像速度.

图3 测试图像及不同排序方法在25%采样下的图像重建结果 (a) 原始测试图像64 × 64 像素; (b)—(f) 排序方法为Db2 小波序、Dct 序、Walsh 序、随机排序和Haar 小波排序下重图像重建结果Fig. 3. Original and recovered images under 25% full sampling: (a) The Lena,with 64 × 64 pixels; (b)−(f) images recovered corresponding to the ordering method of Db2 wavelet order,Dct order,Walsh order,random permutation order and Haar wavelet order,respectively.

图4 测试图3(a)在不同排序下仿真结果 (a) PSNR 值与测量次数m的关系; (b) SSIM 值与测量次数m的关系Fig. 4. Simulation results of the Fig. 3(a): (a) The PSNR versusmth measurement; (b) the SSIM versusmth measurement.

4 实验验证

为验证上述的数值模拟仿真实验及分析结论,开展了800m 距离单像素凝视成像实验. 单像素原理样机调制器DMD 为美国德州仪器公司生产的DLP4100,微镜总数及排列为1920 × 1080,微镜可±12◦翻转,每个微镜10.6 μ m . 实验时从DMD中间区域选取的调制区域为1024 × 1024,16 ×16 个微镜作为一个像素调制单元,因此共64 × 64个像素调制单元,微镜的调制速度设定为8 kHz.首先,计算机按Hadamard 编码原则生成4096 ×4096 大小的Hadamard 矩阵; 其次,将Hadamard矩阵的每一行重新按顺序排列成4096 个64 ×64 大小的矩阵; 然后将4096 个64 × 64 大小的矩阵利用直积的方式变为1024 × 1024 像素,再填0 扩充到1920 × 1080,将扩充后的序列加载到数字微镜阵列DMD的板载内存. DMD 的调制速率设置为8 kHz,采集得到一副图像需要0.512 s. 样机成像主镜为卡塞格林式天文望远镜,口径为280 mm,焦距为2800 mm. 探测器采用日本滨松H10723-20 型PMT,有效波段 350—900 nm. 数据采集卡采用16 位精度的A/D 转换模块. 为获得高SNR 成像结果,实验进行了100 次全测量,对应曝光时间51.2 s. 累加求平均后,分别按5 种规则排序: 按Haar 小波变换、Db2 小波变换和Dct 的S(i) 值由小到大排序,以及按Walsh 序排序和按1 到4096 随机排序,分别截取数据的25%用于图像重建,截取的数据如图5所示. 从截取的数据图5 可见,对于相同图像,不同测量基在数据表征上有显著区别,利用数据的特征分布规律有望进行不成像识别,此处不展开讨论. 数据虽有区别但仍不能直接给出重建效果优劣. 因此仍需利用这5 种排序的1024 个数据重建图像,其结果如图6所示. 图6(a)为成像区域,利用商用测距仪测得目标距离约800 m. 对该区域放大,得到的相机拍摄图像如图6(b)所示,该区域为选定的成像目标区域. 图6(c) 为利用单像素成像实验装置在全采样和测量次数m=4096 条件下获得的成像结果.图6(d)—(h)分别对应排序方法为Db2 小波序、Dct 序、Walsh 序、随机排序和Haar 小波排序下的图像重建结果. 通过对比发现,只有图6(h),即Haar 小波变换排序方法可在25%、测量次数m=1024条件下依然能够成像,且图像大部分信息得到还原. Haar 小波变换排序相当于提升4 倍采样速度和图像重建速度,使得单幅图像采集与重建时间为0.125 s,本实验中对应8 帧/秒@64×64像素成像. 而实际上,算法计算速度上可满足至少300 帧/秒@ 6 4×64 像素的单像素成像,此时,高速单像素成像主要受调制器调制速度的限制.

图5 不同排序方法25%采样下的实验数据,即Haar 小波变换、Db2 小波变换、Dct、Walsh 排序和随机随机排序对应的实验数据Fig. 5. Experimental data of different ordering methods:Haar wavelet,Db2 wavelet and Dct,Walsh order and random order.

图6 室外实验结果和不同排序方法在25%采样下的重建图像 (a)目标区域,对应距离800 m; (b)相机拍摄的目标图像;(c) 无压缩单像素成像,64 × 64 像素,100 幅图像累加结果; (d)—(h) 分别对应排序方法为Db2 小波序、Dct 序、Walsh 序、随机排序和Haar 小波排序下重图像重建结果Fig. 6. Outdoor experiment and recovered images under 25% full sampling with different ordering methods: (a) The target region,with the distance 800 meters; (b) target image captured by camera; (c) image recovered by single-pixel camera,with 64 × 64 pixels,with 100 recovered images averaged; (d)−(h) images recovered corresponding to the ordering method of Db2 wavelet order,Dct order,Walsh order,random permutation order and Haar wavelet order,respectively.

5 结 论

本文为提升单像素成像速度,提出了基于Hadamard 矩阵优化排序的解决方案. 首先,提出了5 种排序方法,分别利用理论分析、数值仿真和800 m 距离单像素成像室外实验进行了研究. 在压缩采样条件下,本文进行了25%全采样的不同测量基排序对单像素成像质量影响的研究. 结果表明,按Haar 小波变换系数绝对值排序时单像素成像效果最优,该方法可利用Walsh-Hadamard 快速变换重建图像,计算速度达300 帧/秒@ 6 4×64 像素,采样速度提升4 倍. 针对排序方法引起成像SNR 变化的原因,从关联成像角度给出了其物理解释: 测量基矩阵元邻域数值相等的区域面积等效于光场二阶相干面积,当光场二阶相干面积随测量基由大到小排序时成像效果最优,数值仿真和室外实验均验证了该理论的正确性. 测量基在Haar 小波变换下排序,实际上是空间频率由低频到高频的过程,本文所提出的研究方法,可用于研究单像素成像中最优测量基选取的一类科学问题,例如,利用Haar 小波变换对测量基图像稀疏度进行控制,有望寻找到除Hadamard 矩阵以外的其他最优观测矩阵,为高SNR、快速单像素成像技术提供了有效的研究思路和方法. 此外,本文提出的物理解释对单像素成像测量基设计过程中参数的选取具有指导意义,给出的Harr 小波变换排序方法可用于提升单像素成像速度,具有实际应用价值.