【概述】

分辨率是成像系统的一个重要参数, 获得高分辨率图像一直是鬼成像系统的一个目标. 本文提出了以成 像系统点扩散函数作为先验知识, 基于稀疏测量的超分辨压缩感知鬼成像重建模型. 搭建了一套计算鬼成像 实验装置, 用于验证该模型对于提高鬼成像系统分辨率的有效性, 并与传统的鬼成像计算模型进行了对比.实验表明, 利用该模型可突破成像系统衍射极限分辨率的限制, 得到超分辨鬼成像.

【关键词:】鬼成像, 压缩感知, 超分辨, 稀疏测量

1、引言

鬼成像(ghost imaging, GI), 是一种通过光场 强度关联测量恢复物体信息的技术. 由于其具有 区别于传统成像的特殊性质, 受到很多研究者的关 注 [1−7] . 在热光鬼成像中, 光源发出的光被分成两 路: 一路作为参考光, 一路作为信号光. 参考光不 经过物体, 在经过一定距离的自由传播后, 其强度 分布信息被一具有空间分辨能力的探测器所探测.信号光照射到待成像物体, 透射光或反射光被一不 具有空间分辨能力的桶探测器收集并记录总光强 值. 通过对两路光场的强度值进行关联计算, 即可 恢复出物体的信息. 由于鬼成像具有抗湍流扰动 能力以及可实现无透镜成像等优点 [8,9] , 因而其在 对地观测 [10,11]、雷达成像 [12]、生命科学 [13]、保密通 信 [14] 等领域具有广泛的应用价值.

分辨率是评价鬼成像性能的一个关键参数, 实 现高分辨率的成像, 对鬼成像投入实际应用具有重 要的促进作用, 如何提高鬼成像的分辨率也成为研 究的热点 [15−17] . 热光鬼成像的分辨率由光场的相 干长度决定, 光场的相干长度越小, 成像分辨率越高[18]. 由于真热光的相干时间很短, 探测难度较 大[19−21], 热光鬼成像较多地使用激光经过旋转的 毛玻璃等散射体产生的赝热光作为光源. 有研究表 明, 对于赝热光源, 在散射体表面光场的相干长度 最小, 随着光场的传播, 相干长度逐渐增大[22,23].因此要实现高分辨率的鬼成像, 可使用计算鬼成像 的光路系统 [24−26] , 将散射体表面的光场成像到物 体表面. 然而, 此时成像系统的衍射极限决定了物 体表面光场的最小相干长度, 限制了鬼成像分辨率 的进一步提高.

近年来, Donoho等[27]和Candès等[28]提出的 压缩感知理论提供了一种新型的信号重建方法. 该 理论指出, 只要信号是可压缩的或在某个变换域下 是稀疏的, 可通过随机的观测矩阵将高维信号投影 到低维空间上, 并通过求解一个优化问题从少量的 投影值重建原始信号. 韩申生等 [29−33] 提出通过引 入稀疏约束, 利用压缩感知理论可实现超分辨鬼 成像.

本文提出使用稀疏测量的探测方法, 并引入成 像系统的先验知识, 利用压缩感知理论实现超越衍 射极限的鬼成像, 得到高于传统稀疏约束鬼成像的 分辨率. 我们利用数字微镜阵列 (digital mircomir-ror device, DMD) 搭建了计算鬼成像的光路系统,并通过实验验证了这一方法提高鬼成像分辨率的 作用.

2、理论2.1 计算鬼成像

图1 为计算鬼成像原理图. 激光照射在空间 光调制器(spatial light modulator, SLM) 上, 空间 光调制器对照在其上的光进行强度调制, 成像透镜L1 将调制后的具有光强分布 A(x, y) 的光场照射到 物体表面, 透镜L2对透过物体的所有光进行收集,由探测器Dt 探测其总光强. 与传统鬼成像相比,计算鬼成像光路的特殊之处在于, 参考光路光强分 布通过计算预先得知, 因而可省略参考光路. 计算 鬼成像只需一无空间分辨能力的桶探测器收集总 光强, 通过计算光场的光强分布和桶探测器得到的 总光强值的二阶关联, 即可得到对物体的成像.

图 1：计算鬼成像原理图

SLM：空间光调制器

L1、L2：透镜

Dt：桶探测器

成像时, SLM 进行 N 次随机调制, 每次调制产 生强度分布为Ai(x, y)的散斑场, Ai(x, y) 表示第i 次探测时坐标为 (x, y) 处的光强值. 桶探测器 Dt记录与每次随机调制对应的N个总光强值Yi. 将 散斑场的强度分布A 与光强测量值Y 进行关联运 算, 得到二阶关联函数:

G(2)(x, y) = ⟨Ai(x, y), Yi⟩− ⟨Ai(x, y)⟩ · ⟨Yi⟩, (1)

其中⟨⟩表示取平均值运算; i = 1,2,··· ,N,G(2) (x,y)即可给出物体的像.

由于光源的涨落性质, GI的信噪比较低. Ferri等 [34] 提出差分鬼成像 (differential GI, DGI) 模型 用于提高成像信噪比, 在桶探测器后用差分值代替 光强值, 去除光场的平均强度涨落, 可以明显地抑制背景噪声, 提取有效的物体信息.G(2) (x, y)= ⟨Ai(x, y), Yi⟩ − ⟨Yi⟩ · ⟨Ai(x, y), Ri⟩, (2)

其中

表示每次探测时，探测矩阵值的总和, i = 1,2,··· ,N, ⟨⟩表示取平均值运算.

2.2 基于压缩感知的鬼成像

压缩感知理论是由Donoho, Candès和Tao等 提出的全新的数学理论. 压缩感知理论指出, 长度 为 n 的信号 X 在某个已知的变换域 Ψ 下稀疏, 利用 一个维度为m × n(m ≪ n)的与变换基底Ψ 不相关 的观测矩阵Φ, 对信号在变换域下的系数进行线性 观测, 得到维度为m × 1的观测值Y . 由于观测值 的维度小于未知信号的维度, 在传统信号恢复理论 中, 信号无法精确求解. 而根据压缩感知理论, 利 用观测值 Y , 变换基底 Ψ 和观测矩阵 Φ, 通过求解 优化问题可以精确重构出原始信号X[35]. 数学模型如下:

如果原始信号在原域上是稀疏的, 则不需要通过变换, 模型简化为

将压缩感知重建模型引入到鬼成像系统中, 得 到基于压缩感知的鬼成像模型:

其中

i为探测次数, T为待成像物体的透射函数, A为探 测矩阵.

2.3 基于稀疏测量的超分辨鬼成像

有研究表明, 测量矩阵的稀疏性对压缩感知图像恢复的效果有影响 [36] . 测量矩阵越稀疏, 图像恢 复的精确度越高, 这是由于使用高稀疏度矩阵进行 测量时, 各次测量获得的信息重叠度低, 图像恢复时各像素点之间的干扰变小. 同时稀疏测量能够更清晰地体现出图像本身的稀疏性, 更有利于压缩感 知算法进行稀疏信号的求解.

实际上, 高分辨率图像的获取也来源于图像的 精确恢复. 要清晰地分辨图像的细节, 就要求图像 中相邻像素点的干扰减小, 并且对图像中细节的分 辨往往意味着恢复出的图像具有更高的稀疏性. 因 此, 可以通过提高测量矩阵的稀疏度来提高鬼成像 的分辨率. 传统压缩感知鬼成像使用 0-1 等概率分 布的伯努利矩阵, 并不具有很高的稀疏性. 为了获 得尽可能高的成像分辨率, 我们提出使用扫描的方 式对物体进行测量. 这是因为扫描时每一次的测量 矩阵中理论上只有一个像素的值为1, 其他像素均 为0, 达到了最高的测量矩阵稀疏度.

另一方面, 在计算鬼成像中, 成像透镜分辨率 受瑞利衍射极限限制, 照射在物体表面的真正实现 测量矩阵功能的光场强度分布与 SLM 上预设的光 场强度调制并不完全相同. 因此, 在重建模型中引 入成像系统的点扩散函数h(x, y; x0, y0), 可以在一 定程度上消除衍射受限成像系统对测量矩阵的影 响, 实现超分辨鬼成像. 基于稀疏测量和点扩散函 数的鬼成像模型:

i 为探测次数, T 为待成像物体的透射函数, SLM 预置的稀疏测量矩阵, SH 为经过光学系统点 扩散函数作用的真实测量矩阵, ⊗表示卷积运算.

3、实验与讨论

为验证我们提出的超分辨鬼成像方法, 利用DMD 的可控性, 搭建了基于稀疏测量的计算超分 辨鬼成像实验装置, 如图 2 所示.

卤素灯与 DMD 组合构成具有时间和空间涨落 特性的热光源. 卤素灯经过中心波长为 550 nm 的 滤光片后, 经过透镜L1汇聚到DMD上32×32的像素区域内, 每个像素尺寸为 27.36 μm × 27.36 μm. DMD 对每个像素上的光进行独立调制, 按照预 设矩阵将部分像素的光反射到光路系统中. 经过DMD 调制后的反射光经过 2f -2f 成像系统成像 到物体表面, 成像透镜L2的直径d = 2 mm, 焦距f = 400 mm. 透过物体的光被汇聚透镜 L3 收集,利用 CCD (charge-coupled devices) 作为桶探测器 探测总光强. 需要强调的是, 透镜 L3 只是用来汇聚 所有光, 没有成像功能, 而 CCD 只是用来测量总光 强, 并没有利用其空间分辨能力, 对物体的成像只 需要桶探测器即可实现.

实验中的待成像物体为双缝, 如图 3 (a1) 所示.双缝缝宽为100 μm, 缝中心间距为200 μm缝高度 为780 μm. 将图中各列求和, 得到物体的一维图 像, 归一化处理之后, 如图 3 (a2) 所示. 根据成像透 镜L2的孔径和焦距, 通过计算可以得到成像系统 的瑞利衍射极限分辨率为268.4 μm, 成像物体双缝 在此成像系统中无法被区分. DMD上一个像素经 过成像系统后到达物体表面的衍射光斑如图 3 (b1)所示, 截面图归一化之后, 如图 3 (b2) 所示, 其半高 宽为259.1 μm, 与计算得到的极限分辨率基本一 致. 为获得准确的点扩散函数形式, 我们根据点扩 散函数的半高宽, 利用 sinc 函数对其进行拟合, 得 到成像系统点扩散函数h(x, y; x0, y0)的具体表达 形式.

传统计算鬼成像实验结果如图4 所示, 测量 矩阵为0-1等概率分布的二值随机矩阵, 采样次数 为1024次, 得到1024个随机矩阵A 以及对应的强 度值 Y . 分别利用 DGI 和压缩感知的算法对数据 进行恢复, 得到结果如图 4 (a1), (b1), 图 4 (a2) 和(b2) 分别为将各列求和得到的一维图像, 图 4 (a2)和 (b2) 纵坐标值已归一化处理. 图中可见, 两种算 法恢复的图像都无法分辨出双缝. 这是由于受成 像系统衍射受限的影响, 物体表面测量矩阵的分辨率超出了双缝的中心间距, 虽然在计算中已经引入 了系统的点扩散函数, 但在普通的非稀疏测量下,各次测量值包含的信息重叠度过高, 并不能准确获 取物体每个部分的信息, DGI 和压缩感知都无法对 物体图像进行准确的恢复. 有研究表明, 在单像素 相机实验中, 通过普通的非稀疏测量, 利用压缩感 知和点扩散函数依然可以实现超分辨成像[31]. 在 单像素相机光路中, 物体经过衍射受限系统成像到DMD上, DMD上的测量矩阵是准确且不受系统干 扰的. 而本文的鬼成像光路中, DMD上的测量矩 阵经过衍射受限系统到达物平面时, 已严重模糊 平滑, 无法有效提取物体信息. 本文受衍射受限系 统影响的是测量矩阵, 而不是待成像的物体, 因而 无法实现超分辨. 由此说明, 与模糊化物体相比,压缩感知成像的分辨率对测量矩阵的模糊程度更 敏感.

为实现稀疏测量, 将DMD设置为每次调制时 只将DMD上1个像素的光向光路中反射, 依次反 射32 × 32像素区域内每个像素的光, 采样次数为1024次. 由于此时进行了扫描探测, 将测量得到的 强度值按像素位置排列, 即可直接得到物体的图 像, 如图 5 (a1), (a2), 由于成像系统分辨率的限制,显然扫描结果是无法分辨双缝的. 引入点扩散函 数, 使用DGI算法对扫描数据进行恢复, 得到结果 如图 5 (b1), (b2). 由于传统鬼成像的二阶关联计 算是一种基于统计平均的算法, 而不是对原始图像 的精确求解, 在引入稀疏测量和点扩散函数的情况 下, DGI 仍然无法获得较高分辨率的成像结果. 使用压缩感知算法对扫描数据进行恢复, 得到结果如 图 5 (c1), (c2). 可明显看出, 使用压缩感知算法进 行恢复的双缝清晰分辨, 实现了超越衍射极限分辨 率的成像. 与图 4 (b1), (b2) 相比可以看出, 稀疏测 量在实现超分辨鬼成像上具有重要的作用, 通过稀 疏测量可以得到比传统压缩感知测量成像结果更 高的分辨率.

根据本文提出的基于压缩感知的超分辨方法,在稀疏测量方式下可实现超分辨鬼成像. 通过最稀 疏的稀疏测量方式—–逐点扫描, 实验验证了该方 法. 随机稀疏矩阵同样满足稀疏测量的方式, 根据 本文方法依然可以实现超分辨鬼成像. 我们通过 仿真实验来验证该问题, 并探讨随机稀疏矩阵的稀 疏性与超分辨成像质量之间的关系. 有研究表明,基于压缩感知的超分辨重建算法对于灰度物体依 然有效 [37] . 本文提出的基于稀疏测量的压缩感知 超分辨方法并没有对物体的透过函数进行限制, 针 对灰度物体同样适用. 因此, 在仿真实验中, 我们 设计有灰度的双缝, 透过率分别为1和0.5, 缝宽为28 μm, 中心间距为 140 μm, 如图 6 (a). 系统点扩 散函数半高宽为182.5 μm, 利用sinc 函数对其进 行模拟. 设计随机矩阵, 比较五种不同稀疏度下压 缩感知超分辨重建的结果, 稀疏度(每次测量中非 零点的个数)依次为1, 5, 50, 250, 500. 仿真结果如 图6所示.

由图6可知, 随机稀疏测量方式下, 利用本文 的方法可以在一定程度上实现超分辨成像, 对于灰 度物体同样适用. 并且, 随机测量矩阵的稀疏性对超分辨成像质量有一定的影响. 测量矩阵越稀疏,采样点之间的距离大于系统点扩散函数宽度的概 率越大, 每次测量采集到的信息重叠度较低, 保证 了测量矩阵对物体信息的准确采集, 因而能够较好 地实现超分辨鬼成像. 当测量矩阵稀疏度较低, 如图6 (b), (c)中为1或5时, 能够较好地满足这一要 求, 因而超分辨成像质量较好. 当稀疏度增大时,重建质量明显变差, 如图 6 (d), (e), (f) 所示. 由此 说明, 测量矩阵越稀疏, 基于压缩感知的超分辨鬼 成像重建质量越高.

4、结论

本文在以成像系统点扩散函数作为先验知识 的条件下, 提出了基于稀疏测量的压缩感知鬼成像 重建模型. 基于 DMD 搭建了计算超分辨鬼成像光 路系统, 实验验证了在稀疏测量方式下, 该模型可 突破成像系统瑞利衍射极限, 实现超分辨鬼成像,促进了鬼成像的实用化进展. 考虑到采样次数对于 实现实时应用的限制, 今后的工作将围绕着压缩采 样超分辨成像进行, 研究亚采样情况下基于压缩感 知的超分辨成像.

参考文献

1. [1]  Pittman T B, Shih Y H, Strekalov D V, Sergienko A V1995 Phy. Rev. A 52 R3429

2. [2]  Gatti A, Brambilla E, Bache M, Lugiato L A 2004 Phy. Rev. Lett. 93 093602

3. [3]  CaiYJ,ZhuSY2005Phy.Rev.E71056607

4. [4]  Valencia A, Scarcelli G, Angelo M D, Shih Y H 2005

   Phy. Rev. Lett. 94 063601

5. [5]  CaoDZ,XiongJ,WangKG2005Phy.Rev.A71

   013801

6. [6]  LiuXF,YaoXR,LiMF,YuWK,ChenXH,Sun ZB,WuLA,ZhaiGJ2013ActaOpt.Sin.62184205(in Chinese) [刘雪峰, 姚旭日, 李明飞, 俞文凯, 陈希浩, 孙 志斌, 吴令安, 翟光杰 2013 物理学报 62 184205]

7. [7]  Zhao S M, Zhuang P 2014 Chin. Phys. B 23 054203

8. [8]  Cheng J 2009 Opt. Express 17 7916

9. [9]  Cheng J, Han S S 2004 Phy. Rev. Lett. 92 093903

10. [10]  ShiDF,FanCY,ZhangPF,ZhangJH,ShenH,Qiao C H, Wang Y J 2012 Opt. Express 20 27992

11. [11]  YuWK,LiS,YaoXR,LiuXF,WuLA,ZhaiGJ2013 Appl. Opt. 52 7882

12. [12]  ZhaoCQ,GongWL,ChenML,LiWR,WangH,Xu W D, Han S S 2012 Appl. Phys. Lett. 101 141123

13. [13]  LiuXF,YaoXR,ChenXH,WuLA,ZhaiGJ2012

    J. Opt. Soc. Am. A 29 1922

[14] YuWK,LiMF,YaoXR,LiuXF,WuLA,ZhaiGJ2014 Opt. Express 22 7133

[15] Ferri F, Magatti D, Gatti A, Bache M, Brambilla E, Lugiato L A 2005 Phys. Rev. Lett. 94 183602

[16] ZhangPL,GongWL,ShenX,HuangDJ,HanSS2009 Opt. Lett. 34 1222

[17] ZhouY,LiuJB,SimonJ,ShihYH2012J.Opt.Soc. Am. B 29 377

[18] Gatti A, Bache M, Magatti D, Brambilla E, Ferri F, Lugiato L A 2006 J. Mod. Opt. 53 739

[19] ZhangD,ZhaiYH,WuLA,ChenXH2005Opt.Lett.30 2354

[20] ChenXH,LiuQ,LuoKH,WuLA2009Opt.Lett.34695

[21] LiuXF,ChenXH,YaoXR,YuWK,ZhaiGJ,Wu L A 2014 Opt. Lett. 39 2314

[22] Gatti A, Magatti D, Ferri F 2008 Phys. Rev. A 78063806

[23] Magatti D, Gatti A, Ferri F 2009 Phys. Rev. A 79053831

[24] Shapiro J H 2008 Phys. Rev. A 78 061802            [25] Bromberg Y, Katz O, Silberberg Y 2009 Phys. Rev. A

79 053840            [26] LiuXF,LiMF,YaoXR,YuWK,ZhaiGJ2013AIP

Adv. 3 052121
            [27] Donoho D L 2004 IEEE Trans. Inf. Theory 52 1289 [28] Candès E J, Romberg, Tao T 2006 IEEE Trans. Inf.

Theory 52 489            [29] DuJ,GongWL,HanSS2012Opt.Lett.371067            [30] Gong W L, Han S S 2012 Phys. Lett. A 376 1519            [31] ChenJ,GongWL,HanSS2013Phys.Lett.A377

1844             [32] Bai X, Li Y Q, Zhao S M 2013 Acta Phys. Sin. 62 044209

(in Chinese) [白旭, 李永强, 赵生妹 2013 物理学报 62

044209]
            [33] ChengJ,HanSS,YanYJ2006Chin.Phys.B152002[34] Ferri F, Magatti D, Lugiato L A, Gatti A 2010 Phys.

Rev. Lett. 104 253603            [35] Candès E J 2006 Proceedings of the International

Congress of Mathematicians 3 1433
            [36] ChenML,LiER,WangH,ShenX,GongWL,HanS

S 2012 Acta Opt. Sin. 32 0503001            [37] Gong W L, Han S S 2014 Appl. Phys. Express 7 062503

原创信息，侵权删！

1)(中国科学院空间科学与应用研究中心,中国科学院复杂航天系统电子信息技术重点实验室,北京 100190)

2)(中国科学院大学,北京 100190)

李龙珍1)2) 姚旭日1)2) 刘雪峰1) 俞文凯1)2) 翟光杰1)†