搜索

x

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

时空域联合编码扩频单光子计数成像方法

沈姗姗 顾国华 陈钱 何睿清 曹青青

引用本文:
Citation:

时空域联合编码扩频单光子计数成像方法

沈姗姗, 顾国华, 陈钱, 何睿清, 曹青青

Time-space united coding spread spectrum single photon counting imaging method

Shen Shan-Shan, Gu Guo-Hua, Chen Qian, He Rui-Qing, Cao Qing-Qing
PDF
HTML
导出引用
  • 本文结合空间编码单像素成像技术和扩频时间编码扫描成像技术, 提出一种时空域联合编码扩频单光子计数成像方法. 该方法具备可避免距离模糊、大时宽带宽积的优势, 并且在噪声干扰下, 能够准确恢复距离像. 本文推导了基于单光子探测的时空域联合相关非线性探测模型、成像正向模型和信噪比模型, 并通过凸优化反演算法恢复深度图像, 理论模型和仿真实验均证明, 与传统的基于空间编码的单像素成像方法相比, 本方法提高了重建的质量. 其中, 采用m序列作为时间编码, 成像的噪声鲁棒性更高. 和传统的空间编码单像素成像技术相比, 本文提出方法的成像均方误差降低了4/5, 引入二次相关方法后, 成像均方误差降低了9/10. 本文所提出的成像架构为非扫描激光雷达成像方法提供了新思路.
    In this paper, we demonstrate a new imaging architecture called time-space united coding spread spectrum single photon counting imaging technique by combining the space coding based single-pixel imaging technology and spread spectrum time coding based scanning imaging technology. This method has the advantages of range ambiguity-free and large time-bandwidth product. Under the interference of noise, this method can accurately restore depth images. In this work, the time-space united correlation nonlinear detection model based on single photon detection, forward imaging model and signal-to-noise ratio model is derived, and the depth image is restored by convex optimization inversion algorithm. The theoretical model and simulation experiments show that compared with the traditional single pixel imaging method based on spatial coding, this method improves the quality of scene reconstruction. Using m-sequence as time coding, imaging has higher noise robustness. In addition, compared with the traditional space coding single pixel imaging technology, the imaging mean square error of the proposed method is reduced by 4/5 and the imaging mean squared error is reduced by 9/10 after introducing the second correlated method. The proposed imaging architecture in this paper may provide a new path for non-scanning lidar imaging methods.
      通信作者: 沈姗姗, ssssoner@niit.edu.cn
    • 基金项目: 国家自然科学基金(批准号: 62205144, 61905108, 62005128)、江苏省高等学校自然科学研究项目(批准号: 19KJB140010, 21KJB510049)和南京工业职业技术大学人才引进基金项目(批准号: YK-20-03-02, YK-20-03-03)资助的课题.
      Corresponding author: Shen Shan-Shan, ssssoner@niit.edu.cn
    • Funds: Project supported by the National Natural Science Foundation of China (Grant Nos. 62205144, 61905108, 62005128), the Natural Science Research Project of Higher Education Institutions of Jiangsu Province, China (Grant Nos. 19KJB140010, 21KJB510049), and the Talent Introduction Foundation of Nanjing Vocational University of Industry Technology, China (Grant Nos. YK-20-03-02, YK-20-03-03).
    [1]

    Rehain P, Sua Y M, Zhu S Y, Dickson I, Muthuswamy B, Ramanathan J, Shahverdi A, Huang Y P 2020 Nat. Commun. 11 921Google Scholar

    [2]

    Tachella J, Altmann Y, Mellado N, McCarthy A, Tobin R, Buller G S, Tourneret J Y, McLaughlin S 2019 Nat. Commun. 10 4984Google Scholar

    [3]

    Brock J C, Purkis S J 2009 J. Coastal Res. 25 1Google Scholar

    [4]

    Howland G A, Lum D J, Ware M R, Howell J C 2013 Opt. Express 21 23822Google Scholar

    [5]

    Zhao Y N, Hou H Y, Han J C, Liu H C, Zhang S H, Cao D Z, Liang B L 2021 Opt. Lett. 46 4900Google Scholar

    [6]

    Radwell N, Johnson S D, Edgar M P, Higham C F, Murray-Smith R, Padget M J 2019 Appl. Phys. Lett. 115 231101Google Scholar

    [7]

    Li F Q, Chen H J, Pediredla A, Yeh C, He K, Veeraraghavan A, Cossairt O 2017 Opt. Express 11 31096

    [8]

    Liu X L, Shi J H, Sun L, Li Y H, Fan J P, Zeng G H 2020 Opt. Express 28 8132Google Scholar

    [9]

    Asmann A, Mota J F C, Stewart B D, Wallace A M 2019 IEEE Trans. Comput. Imaging 8 385

    [10]

    Misra P, Hu W, Yang M R, Duarte M, Jha S 2017 IEEE Trans. Mob. Comput. 16 2037Google Scholar

    [11]

    Krichel N J, McCarthy A, Buller G S 2010 Opt. Express 18 9192Google Scholar

    [12]

    Zhang Y F, He Y, Yang F, Luo Y, Chen W B 2016 Chin. Opt. Lett. 14 111101Google Scholar

    [13]

    Yu Y, Liu B, Chen Z 2018 Appl. Opt. 57 7733Google Scholar

    [14]

    Norman D M, Gardner C S 1988 Appl. Opt. 27 3650Google Scholar

    [15]

    Ding Y C, Wu H X, Gao X L, Wu B, Shen Y H 2022 J. Opt. Soc. Am. A 39 206

    [16]

    Hiskett P A, Parry C S, McCarthy A, Buller G S 2008 Opt. Express 16 13685Google Scholar

    [17]

    Yu Y, Liu B, Chen Z, Kang Li Z K 2020 Sensors 20 2204Google Scholar

    [18]

    Wu F, Yang L, Chen X L, Li Z H, Wu G 2022 Chin. Opt. Lett. 20 021202Google Scholar

    [19]

    Bioucas-Dias J M, Figueiredo M A T 2007 IEEE Trans. Image Process. 16 2992Google Scholar

    [20]

    Mandel L, Wolf E 1995 Optical Coherence and Quantum Optics (Cambridge: Cambridge University Press) p123

    [21]

    冒添逸, 陈钱, 何伟基, 庄佳衍, 邹云浩, 戴慧东, 顾国华 2016 物理学报 65 08427Google Scholar

    Mao T Y, Chen Q, He W J, Zhuang J Y, Zou Y H, Dai H D, Gu G H 2016 Acta Phys. Sin. 65 08427Google Scholar

    [22]

    Gatt P, Johnson S, Nichols T 2009 Appl. Opt. 48 3261Google Scholar

    [23]

    沈姗姗, 陈钱, 何伟基, 陈云飞, 尹文也, 戴慧东 2014 光学学报 34 1012001Google Scholar

    Shen S S, Chen Q, He W J, Chen Y F, Yin W Y, Dai H D 2014 Acta Opt. Sin. 34 1012001Google Scholar

    [24]

    Dai H D, Gu G H, He W J, Liao F J, Zhuang J Y, Liu Xiao J, Chen Q 2014 Appl. Opt. 53 6619Google Scholar

    [25]

    Zierler N, Siam J 1996 Electron. Commun. Eng. J. 1996 79

    [26]

    Molisch A F 2011 Wireless Communications (2nd Ed.) (New York: John Wiley & Sons Ltd.) p334

    [27]

    Shen S S, Chen Q, He W J, Gu G H 2020 Opt. Quantum Electron. 52 2020Google Scholar

    [28]

    Shen S S, Chen Q, He W J, Wang Y Q 2017 Chin. Opt. Lett. 15 090101Google Scholar

  • 图 1  时空域联合编码扩频单光子计数成像架构

    Fig. 1.  Time-space united coding spread spectrum single photon counting imaging method architecture.

    图 2  不同压缩比的时空域联合编码扩频单光子计数成像 (a) “飞机”深度真值; (b)−(g)不同压缩比的深度成像 (b) 5%, (c)15%, (d) 25%, (e) 45%, (f) 65%, (g) 85%. 颜色条表示深度数值, 单位为m

    Fig. 2.  Time-space united coding spread spectrum single photon counting image with different compression proportions: (a) Depth ground truth of ‘airplane’; (b)−(g) depth maps with different compression proportions of (b) 5%, (c)15%, (d) 25%, (e) 45%, (f) 65% and (g) 85%. Colorbars show the depth with unit of meter.

    图 3  时空域联合编码扩频单光子计数成像性能模拟 (a) 1码概率为0.1的伪随机序列深度图像; (b) 1码概率为0.3的伪随机序列深度图像; (c) Gold序列的深度图像; (d) m序列的深度图像; (e) 死时间对成像性能的影响理论模拟; (f) 死时间为200 ns; (g) 死时间为20 ns; (h) 死时间为1ns. 颜色条表示深度数值, 单位为m

    Fig. 3.  Time-space united coding image spread spectrum single photon counting imaging performance simulation: (a) Depth maps by pseudo-random sequences with ‘1’ bit of 0.1; (b) depth maps by pseudo-random sequences with ‘1’ bit of 0.3; (c) depth maps by Gold sequences; (d) depth maps by m-sequences; (e) theoretical simulation of dead time influence on imaging performance; (f) dead time is 200 ns; (g) dead time is 20 ns; (d) dead time is 1 ns. Colorbars show the depth information with unit of meter.

    图 4  不引入二次相关法的时空域联合编码扩频单光子计数深度成像, 码长为 (a) 2048, (b) 4096, (c) 8192, (d) 16384. 引入二次相关法的时空域联合编码扩频单光子计数深度图像, 码长为(e) 2048, (f) 4096, (g) 8192, (h)16384; (i) 传统的基于空间编码的单像素深度图像; (j) 引入二次相关法(点虚线)和不引入二次相关法(虚线)的深度MSE

    Fig. 4.  Time-space united coding image spread spectrum single photon counting imaging without second correlated method by code length of (a) 2048, (b) 4096, (c) 8192, (d)16384. Depth maps simulations with second correlated method by code length of (e) 2048, (f) 4096, (g) 8192 and (h) 16384. (i) Traditional single pixel imaging method based on space coding. (j) The depth MSE with second correlated method (dot dashed line) and without second correlated method (dashed line).

    图 5  成像边缘检测 (a) 不引入二次相关法的时空域联合编码扩频单光子成像边缘检测; (b) 引入二次相关法的时空域联合编码扩频单光子计数成像边缘检测; (c) 传统的基于空间编码的单像素成像边缘检测

    Fig. 5.  Image edge detection: (a) Time-space united coding spread spectrum single photon counting imaging without second correlated method image edge detection; (b) time-space united coding spread spectrum single photon counting imaging with second correlated method image edge detection; (c) traditional single pixel imaging method based on space coding image edge detection.

    图 6  单光子计数扫描成像测试图像深度重建仿真 (a) 玩偶深度真值; (b) $ {m_n} = 10 $, $ {\text{MSE = 0}}{\text{.176 m}} $, 传统的基于空间编码的单像素成像深度图; 在(c) ${m_{\rm{n}}} = 5$, $ {\text{MSE = 0}}{\text{.019 m}} $, (d) ${m_{\rm{n}}} = 8$, ${\text{MSE = 0}}{\text{.03\; m}}$, (e) $ {m_n} = 10 $, ${\text{MSE = 0}}{\text{.035 \;m}}$条件下, 不引入二次相关法的时空域联合编码扩频单光子计数深度图; 在(f) ${m_{\rm{n}}} = 5$, ${\text{MSE = 0}}{\text{.005\; m}}$, (g) ${m_{\rm{n}}} = 8$, ${\text{MSE = 0}}{\text{.011\; m}}$, (h) ${m_{\rm{n}}} = 10$, ${\text{MSE = 0}}{\text{.017 \;m}}$条件下, 引入二次相关法的时空域编码扩频单光子计数深度图, 深度单位为m

    Fig. 6.  Depth reconstruction simulation by using single photon counting scanning imaging test image: (a) Ground truth of ‘doll’; (b) the depth map by traditional single pixel imaging method based on space coding image when ${m_{\rm{n}}} = 10$, ${\text{MSE = 0}}{\text{.176\; m}}$. Depth maps by time-space united coding spread spectrum single photon counting imaging without second correlated method when (c) ${m_{\rm{n}}} = 5$, ${\text{MSE = 0}}{\text{.019\; m}}$, (d) ${m_{\rm{n}}} = 8$, ${\text{MSE = 0}}{\text{.03\; m}}$ and (e) ${m_{\rm{n}}} = 10$, ${\text{MSE = 0}}{\text{.035\; m}}$. Depth maps by time-space united coding spread spectrum single photon counting imaging with second correlated method when (f) ${m_{\rm{n}}} = 5$, ${\text{MSE = 0}}{\text{.005\; m}}$, (g) ${m_{\rm{n}}} = 8$, ${\text{MSE = 0}}{\text{.011\; m}}$, and (h) ${m_{\rm{n}}} = 10$, ${\text{MSE = 0}}{\text{.017\; m}}$. The depth unit is m.

    图 7  成像边缘检测 (a) 不引入二次相关法的时空域联合编码扩频单光子成像边缘检测; (b) 引入二次相关法的时空域联合编码扩频单光子成像边缘检测; (c) 传统的基于空间编码的单像素成像边缘检测

    Fig. 7.  Image edge detection: (a) Time-space united coding spread spectrum single photon counting imaging without second correlated method image edge detection; (b) time-space united coding spread spectrum single photon counting imaging with second correlated method image edge detection; (c) traditional single pixel imaging method based on space coding image edge detection.

    图 8  噪声$ {m_n} $=10不同的凸优化反演方法的图像重建 (a) TVAL3; (b) GPSR; (c) FISTA

    Fig. 8.  Different convex optimization inversion methods imaging reconstruction when$ {m_n} $=10: (a) TVAL3; (b) GPSR; (c) FISTA.

    表 1  不同重建方法的性能统计

    Table 1.  Performance statistical record of different reconstruction methods.

    MSE/m Time/s
    Method${m_{\rm{n} } } = 5$ ${m_{\rm{n}}} = 10$ ${m_{\rm{n}}} = 5$ ${m_{\rm{n}}} = 10$
    TVAL30.033 0.011 3557
    GPSR0.060 0.018 1324
    FISTA0.042 0.019 1624
    下载: 导出CSV
  • [1]

    Rehain P, Sua Y M, Zhu S Y, Dickson I, Muthuswamy B, Ramanathan J, Shahverdi A, Huang Y P 2020 Nat. Commun. 11 921Google Scholar

    [2]

    Tachella J, Altmann Y, Mellado N, McCarthy A, Tobin R, Buller G S, Tourneret J Y, McLaughlin S 2019 Nat. Commun. 10 4984Google Scholar

    [3]

    Brock J C, Purkis S J 2009 J. Coastal Res. 25 1Google Scholar

    [4]

    Howland G A, Lum D J, Ware M R, Howell J C 2013 Opt. Express 21 23822Google Scholar

    [5]

    Zhao Y N, Hou H Y, Han J C, Liu H C, Zhang S H, Cao D Z, Liang B L 2021 Opt. Lett. 46 4900Google Scholar

    [6]

    Radwell N, Johnson S D, Edgar M P, Higham C F, Murray-Smith R, Padget M J 2019 Appl. Phys. Lett. 115 231101Google Scholar

    [7]

    Li F Q, Chen H J, Pediredla A, Yeh C, He K, Veeraraghavan A, Cossairt O 2017 Opt. Express 11 31096

    [8]

    Liu X L, Shi J H, Sun L, Li Y H, Fan J P, Zeng G H 2020 Opt. Express 28 8132Google Scholar

    [9]

    Asmann A, Mota J F C, Stewart B D, Wallace A M 2019 IEEE Trans. Comput. Imaging 8 385

    [10]

    Misra P, Hu W, Yang M R, Duarte M, Jha S 2017 IEEE Trans. Mob. Comput. 16 2037Google Scholar

    [11]

    Krichel N J, McCarthy A, Buller G S 2010 Opt. Express 18 9192Google Scholar

    [12]

    Zhang Y F, He Y, Yang F, Luo Y, Chen W B 2016 Chin. Opt. Lett. 14 111101Google Scholar

    [13]

    Yu Y, Liu B, Chen Z 2018 Appl. Opt. 57 7733Google Scholar

    [14]

    Norman D M, Gardner C S 1988 Appl. Opt. 27 3650Google Scholar

    [15]

    Ding Y C, Wu H X, Gao X L, Wu B, Shen Y H 2022 J. Opt. Soc. Am. A 39 206

    [16]

    Hiskett P A, Parry C S, McCarthy A, Buller G S 2008 Opt. Express 16 13685Google Scholar

    [17]

    Yu Y, Liu B, Chen Z, Kang Li Z K 2020 Sensors 20 2204Google Scholar

    [18]

    Wu F, Yang L, Chen X L, Li Z H, Wu G 2022 Chin. Opt. Lett. 20 021202Google Scholar

    [19]

    Bioucas-Dias J M, Figueiredo M A T 2007 IEEE Trans. Image Process. 16 2992Google Scholar

    [20]

    Mandel L, Wolf E 1995 Optical Coherence and Quantum Optics (Cambridge: Cambridge University Press) p123

    [21]

    冒添逸, 陈钱, 何伟基, 庄佳衍, 邹云浩, 戴慧东, 顾国华 2016 物理学报 65 08427Google Scholar

    Mao T Y, Chen Q, He W J, Zhuang J Y, Zou Y H, Dai H D, Gu G H 2016 Acta Phys. Sin. 65 08427Google Scholar

    [22]

    Gatt P, Johnson S, Nichols T 2009 Appl. Opt. 48 3261Google Scholar

    [23]

    沈姗姗, 陈钱, 何伟基, 陈云飞, 尹文也, 戴慧东 2014 光学学报 34 1012001Google Scholar

    Shen S S, Chen Q, He W J, Chen Y F, Yin W Y, Dai H D 2014 Acta Opt. Sin. 34 1012001Google Scholar

    [24]

    Dai H D, Gu G H, He W J, Liao F J, Zhuang J Y, Liu Xiao J, Chen Q 2014 Appl. Opt. 53 6619Google Scholar

    [25]

    Zierler N, Siam J 1996 Electron. Commun. Eng. J. 1996 79

    [26]

    Molisch A F 2011 Wireless Communications (2nd Ed.) (New York: John Wiley & Sons Ltd.) p334

    [27]

    Shen S S, Chen Q, He W J, Gu G H 2020 Opt. Quantum Electron. 52 2020Google Scholar

    [28]

    Shen S S, Chen Q, He W J, Wang Y Q 2017 Chin. Opt. Lett. 15 090101Google Scholar

计量
  • 文章访问数:  238
  • PDF下载量:  10
  • 被引次数: 0
出版历程
  • 收稿日期:  2022-07-19
  • 修回日期:  2022-09-29
  • 上网日期:  2022-10-27
  • 刊出日期:  2023-01-20

时空域联合编码扩频单光子计数成像方法

  • 1. 南京工业职业技术大学航空工程学院, 南京 210023
  • 2. 南京理工大学电子工程与光电技术学院, 智能感知和微光成像实验室, 南京 210094
  • 3. 南京工程学院信息与通信工程学院, 南京 211167
  • 通信作者: 沈姗姗, ssssoner@niit.edu.cn
    基金项目: 国家自然科学基金(批准号: 62205144, 61905108, 62005128)、江苏省高等学校自然科学研究项目(批准号: 19KJB140010, 21KJB510049)和南京工业职业技术大学人才引进基金项目(批准号: YK-20-03-02, YK-20-03-03)资助的课题.

摘要: 本文结合空间编码单像素成像技术和扩频时间编码扫描成像技术, 提出一种时空域联合编码扩频单光子计数成像方法. 该方法具备可避免距离模糊、大时宽带宽积的优势, 并且在噪声干扰下, 能够准确恢复距离像. 本文推导了基于单光子探测的时空域联合相关非线性探测模型、成像正向模型和信噪比模型, 并通过凸优化反演算法恢复深度图像, 理论模型和仿真实验均证明, 与传统的基于空间编码的单像素成像方法相比, 本方法提高了重建的质量. 其中, 采用m序列作为时间编码, 成像的噪声鲁棒性更高. 和传统的空间编码单像素成像技术相比, 本文提出方法的成像均方误差降低了4/5, 引入二次相关方法后, 成像均方误差降低了9/10. 本文所提出的成像架构为非扫描激光雷达成像方法提供了新思路.

English Abstract

    • 单光子计数三维成像技术[1-3]广泛应用于机器视觉、遥感和无人驾驶等领域. 其中, 利用目标场景的可压缩性, 从欠采样线性投影序列中重建三维图像是近年来广泛研究的单光子计数三维成像方法. 该方法在空间光调制器(digital micro mirror device, DMD)上加载一系列的随机编码. 进而在空间域调制光信号, 单光子探测器(single photon avalanche detector, SPAD)用于测量场景和随机编码的内积, 从远小于像素大小的投影中获取目标信息, 并在稀疏约束下根据欠采样测量值重建图像. Howland等[4]从当前测量的信号中减去前一次测量的参考信号以削弱噪声能量, 实现32×32分辨率高帧速动态目标跟踪. Zhao等[5]采用Hadamard基提出一种高效的单像素全彩成像方法. 在微光环境下, Radwell等[6]提出了一种基于深度学习的单像素成像优化算法. 李凤强等[7]利用L2范数进行约束, 并采用TwIST(two-step iterative shrinkage/thresholding)反演方法重建图像, 所提出的100万像素的压缩感知(compressed sensing, CS)激光雷达的成像空间分辨率提高了4倍. Liu等[8]提出一种SP3I技术以提高系统的接收效率和输出信噪比(signal-to-noise ratio, SNR). Asmann等[9]提出一种阵列的压缩感知激光雷达系统, 成像信噪比逼近全帧互相关方法的成像信噪比, 并分析了全变分增强拉格朗日交替方向法(total variation augmented Lagrangian alternation direction algorithm, TVAL3)和梯度投影稀疏重建法(gradient projection for sparse reconstruction, GPSR)的重建结果. 传统的基于空间编码的单像素激光雷达通常以固定频率(10 kHz至10 MHz)发射激光脉冲, 限制了最大不模糊距离和信号探测速率. 此外, 分辨率和带宽之间的矛盾难以调和, 成像结果易受噪声影响, 使得此类方法成像信噪比低, 仅限于室内和短距离成像的应用场合.

      扩频通信技术广泛应用于无线传感器网络[10]、全球定位系统和扫描式激光雷达[11-18]等领域, 其中, 基于扩频通信时间编码的单光子计数扫描成像方法[11-18], 通过接收端的匹配滤波器实现脉冲压缩, 高效提高输出信噪比并减少信号衰减的概率, 具有可避免距离模糊、大时宽带宽积、低功耗和易于单片集成的优势. 该方法通常采用伪随机序列进行时间编码, Ding等[15]采用加权的多距离时间相关函数作为目标函数, 通过模拟退火方法优化伪随机序列结构, 提高时间相关函数的距离分辨能力. Hiskett等[16]设计“1”码概率为0.1的序列, 并指出若“1”码概率增加, 噪声相关水平也会增加.

      本文结合空间编码单像素成像技术和扩频时间编码扫描成像技术, 提出了一种基于时空域联合编码的扩频单光子计数三维成像架构. 本方法采用伪随机序列在时间域调制激光脉冲, 替代激光脉冲的周期性调制方法, 并将调制后的光脉冲投影到DMD, 在空间域调制激光脉冲. SPAD用于探测返回的光脉冲信号, 最后将接收信号和参考信号互相关, 从而建立时空域联合相关非线性探测的数学模型, 该模型将场景点多路复用到单光子探测器上, 进而将压缩感知凸优化反演技术(如TwIST[19])引入到新的成像架构中, 由此, 建立成像正向模型和输出信噪比模型, 并研究死时间对成像质量的影响规律. 仿真实验和理论模型共同研究表明, 通过优化时域编码的平衡度, 可降低环境噪声和暗计数的干扰, 采用更长的序列可以获得更高的信噪比和深度确度. 此外, 针对一般的时域编码序列, 理论模型和仿真实验表明, 引入的二次相关重建方法 [20,21]可降低噪声影响, 提高成像质量.

    • 扩频时间编码单光子计数成像系统中包括伪随机序列发生器, 其用于生成0, 1码流, 码为1时触发激光器发射光脉冲, 激光脉冲经目标后漫反射得到回波信号, 经由光学系统接收, 到达SPAD感光面. SPAD启动一次探测并驱动光子到达时间记录仪(time-to-digital converter, TDC)产生时间值. 经过上述发射接收过程, 完成单个像素的有效探测. 定义成像场景的反射率为$ r(x, y) $, 激光脉冲由SPAD和TDC接收后, 通过匹配滤波器进行脉冲压缩, 也即将接收序列$ s(t) $与发送副本序列$ f(t) $互相关, 则像素$ {(x, y)^{{\rm{th}}}} $的互相关信号为

      $ \begin{split} & C(x,y,\tau ) = \displaystyle\sum\limits_{t = 0}^T {{P_R}r(x,y)s(t)f(t + \tau )} \\ =\;&\left\{ \begin{aligned} &r(x,y){P_{\rm{R}}}\left[1 - \dfrac{{T + 1}}{T}\dfrac{{\left| \tau \right|}}{{\Delta t}}\right], \quad |\tau| \leqslant \Delta t , \\ &- \dfrac{1}{T}, \quad\qquad \qquad \qquad \dfrac{{T\Delta t}}{{\text{2}}} \geqslant \left| \tau \right| > \Delta t, \end{aligned} \right. \end{split}$

      其中$ T $表示一个周期的序列长度; $ {P_{\rm{R}}} $为信号光功率; $ \Delta t $为码元宽度. 选取适当的距离门信号可以消除后向散射峰值, (1)式中第一个最高的互相关峰值对应的横坐标$ \tau (x, y) $为从场景目标反射的光子到达时间的估计值. 因此, 可采用最大似然估计法(maximum likelihood estimation, ML)[14]估计$\tau {(x, y)_{{\rm{ML}}}}$:

      $ \tau {(x,y)_{\rm ML}} = \arg \max \left[C\left(x,y,\tau \right)\right] . $

      目标的距离值为

      $ d(x,y) = ({{\tau (x,y) \times c}})/{2} \text{, } $

      其中$ c $为光速.

    • 本文提出的时空域联合编码扩频单光子技术成像架构如图1所示. 采用基于可编程逻辑器件的序列发生器发送扩频伪随机序列. 并驱动激光器发射调制后的激光脉冲, 在时间域完成激光脉冲的编码; 采用空间光调制器调制返回的激光脉冲, 在空间域完成激光脉冲的编码. 采用偏振分束镜将从空间光调制器反射的激光光束与从目标反射的激光光束分开, 透镜收集时空域联合编码的激光光束并导入SPAD, 信号处理器记录光子到达时间点, 并且和模板序列做互相关运算, 通过后续信号处理重建目标的三维图像. 设先后被第$ i $次时域扩频伪随机序列和空域随机矩阵调制后的光束, 依次通过偏振分束镜和透镜, 到达单光子探测器的总光子计数率记为$ {R_i}(t) $:

      图  1  时空域联合编码扩频单光子计数成像架构

      Figure 1.  Time-space united coding spread spectrum single photon counting imaging method architecture.

      $ \begin{split} \;\\ {R_i}(t) =\;& \bigg[{P_{R}}\displaystyle\sum\limits_{y = 1}^N \displaystyle\sum\limits_{x = 1}^N {{\boldsymbol{\varPhi}} _i}(x,y)r(x,y){s_i}(t) \\ &+{B_{i1}}(x,y,t) \bigg] + {B_{i2}}(t) , \end{split} $

      式中$ N $表示$ x $轴或$ y $轴方向的像素个数; $ {{\boldsymbol{\varPhi}} _i}(x, y) $表示第$ i $次调制的随机矩阵; $ {B_{i1}}(x, y, t) $表示空域外部环境噪声光子计数率; $ {B_{i2}}(t) $表示时域外部环境噪声光子计数率. 第$ i $次时空域联合调制后的光束被SPAD探测后产生响应, 该过程为非齐次泊松过程. 假设探测器的量子效率为$ \eta $, 则探测到的光子计数率为$ \eta {R_i}(t) $. 探测器暗电流产生响应的过程为齐次泊松过程, $ d $表示暗电流响应速率, 可得探测到的光子计数率为

      $\begin{split} \lambda (t) =\; &\eta {R_i}(t) + d = \eta \bigg\{ \bigg[\displaystyle\sum\limits_{y = 1}^N \displaystyle\sum\limits_{x = 1}^N {P_R}{\varPhi _i}(x,y)r(x,y){s_i}(t) \\ &+{B_{i1}}(x,y,t) \bigg] + {B_{i2}}(t)\bigg\} + d .\\[-16pt] \end{split}$

      采用相关接收技术, 将(5)式和模板序列$ f(t) $互相关, 时空域联合相关的非线性探测模型为

      $ \begin{split} {C_i}(\tau ) =\;& \displaystyle\sum\limits_{t = 0}^T {\lambda (t)} f(t + \tau ) = \sum\limits_{t = 0}^{T} {[\eta {R_i}(t) + d} ]f(t + \tau ) = \eta \Bigg[\displaystyle\sum\limits_{t = 0}^T {\sum\limits_{x = 1}^N {\sum\limits_{y = 1}^N {{P_{\rm{R}}}{{\boldsymbol{\varPhi}} _i}(x,y)r(x,y){s_i}(t)f(t + \tau )} } } \\ &+\underbrace {\displaystyle\sum\limits_{t = 0}^T {\displaystyle\sum\limits_{x = 1}^N {\displaystyle\sum\limits_{y = 1}^N {{B_{i1}}(x,y,t)f(t + \tau )} } } }_{{B_1'}(\tau )} + \underbrace {\displaystyle\sum\limits_{t = 0}^T {{B_{i2}}(t)f(t + \tau )} }_{{B_2'}(\tau )}\Bigg] + \underbrace {\displaystyle\sum\limits_{t = 0}^T {df(t + \tau )} }_{{B_d}(\tau )}\\ =\;&\eta \displaystyle\sum\limits_{t = 0}^T {\displaystyle\sum\limits_{x = 1}^N {\displaystyle\sum\limits_{y = 1}^N {{\varPhi _i}(x,y)r(x,y){s_i}(t)f(t + \tau )} } } + b(\tau ), \end{split} $

      (6)式中, 第一项为时空相关信号, 外部环境噪声$ {B_{i1}}(x, y, t) $$ {B_{i2}}(t) $与模板序列互相关后分别得到噪声项$B'_1(\tau )$$B'_2(\tau )$, 暗计数噪声$ d $和模板序列互相关得到噪声项$ {B_d}(\tau ) $.

      $b(\tau ) = B'_1(\tau ) + B'_2(\tau ) + {B_d}(\tau )$, 图1中返回的目标信号向量${\boldsymbol{p}} \in {\mathbb{R}^{Q \times 1}}$先被投影至空间光调制器, 并被随机矩阵调制, 空间调制后的光信号被SPAD接收. 本项目考虑多次调制的成像过程, 将(6)式写为矩阵形式, 基于时空域联合编码的扩频单光子计数成像正向模型为

      $ {\boldsymbol{C}}(\tau ) = {\boldsymbol{\varPhi}} {\boldsymbol{p}}(\tau ) + {\boldsymbol{b}}(\tau ) \text{, } $

      其中

      $ {\boldsymbol{p}}(\tau ) = \left\{ \begin{aligned} &r\left[1 - \dfrac{{T + 1}}{T}\dfrac{{\left| \tau \right|}}{{\Delta t}}\right],\quad \left| \tau \right| \leqslant {\Delta t} , \\ &- \dfrac{1}{T},\qquad \qquad\qquad \;\; \dfrac{{T\Delta t}}{{\text{2}}} \geqslant \left| \tau \right| > {\Delta t} {\text{. }} \end{aligned} \right.{\kern 1pt} $

      测量向量为 $ {\boldsymbol{C}} \in {\mathbb{R}^{M \times 1}} $, 空间光调制器加载的测量矩阵为$ {\boldsymbol{\varPhi}} \in {\mathbb{R}^{M \times Q}} $, 投影次数为M, 噪声随机向量为 ${\boldsymbol{ b}} \in {\mathbb{R}^{Q \times 1}} $, 像素总数为$Q = {N^2}$. (6)式表明$ {C_i}(\tau ) $包括信号光子$ {C_{i{\text{s}}}} $和噪声光子$ {C_{i{\text{n}}}} $, 本文通过给出$ {C_i}(\tau ) $中变量的统计特征来建立信噪比模型. 在微光环境下, 单光子探测服从泊松过程[22], 假设噪声和信号是独立的, 且期望和方差相等. 则第$ i $次调制得到的$ {C_i}(\tau ) $中包含的信号光子$ {C_{i{\text{s}}}} $的数学期望为

      $ E\left[{P_R}\displaystyle\sum\limits_{t = 0}^T {\displaystyle\sum\limits_{x = 1}^N {\displaystyle\sum\limits_{y = 1}^N {{{\boldsymbol{\varPhi}} _i}(x,y)r(x,y){s_i}(t)f(t + \tau )} } } \right] =T{R_a}{M_2}\overline r \int\nolimits_{{t_1}}^{{t_2}} {\frac{{{P_R}}}{{hv}}} {\rm{d}}t{\text{ = }}T{R_a}{M_2}\overline r {m_{\text{s}}} \text{, } $

      其中$ {R_a} $为伪随机序列“1”码的概率. (6)式中第$ i $次调制得到的$ {C_i}(\tau ) $中包含噪声光子${C_{i{\rm{n}}}}$的方差为

      $\begin{split} &{\rm{Var}}\left[ \displaystyle\sum\limits_{t = 0}^T {\displaystyle\sum\limits_{x = 1}^N {\displaystyle\sum\limits_{y = 1}^N {B_{i1}}(x,y,t)f(t + \tau )} } \right] + {\rm{Var}}\left[ \displaystyle\sum\limits_{t = 0}^T {B_{i2}}(t)f(t + \tau )\right] + {\rm{Var}}\left[\displaystyle\sum\limits_{t = 0}^T d f(t + \tau ) \right] \\ =\;& {N^2}\displaystyle\sum\limits_{t = 0}^T {E[{B_{i1}}(t)]} f(t + \tau ) + \displaystyle\sum\limits_{t = 0}^T {E[{B_{i2}}(t)]f(t + \tau )} +\displaystyle\sum\limits_{t = 0}^T {E[d]f(t + \tau )} \\ =\;& {N^2}\displaystyle\sum\limits_{t = 0}^T {\int\nolimits_{{t_1}}^{{t_2}} {\dfrac{{({P_{B1}})}}{{hv}}} {\rm{d}}tf(t + \tau )} +\displaystyle\sum\limits_{t = 0}^T {\int\nolimits_{{t_1}}^{{t_2}} {\dfrac{{({P_{B2}})}}{{hv}}} {\rm{d}}tf(t + \tau )} +\displaystyle\sum\limits_{t = 0}^T {\int\nolimits_{{t_1}}^{{t_2}} {\dfrac{{({P_{B3}})}}{{hv}}} {\rm{d}}t f(t + \tau )} \\ = &\;{N^2}{m^{(1)}_n}\displaystyle\sum\limits_{t = 0}^T {f(t + \tau )} + {m^{(2)}_n}\displaystyle\sum\limits_{t = 0}^T {f(t + \tau )} +{m^{(3)}_n}\displaystyle\sum\limits_{t = 0}^T {f(t + \tau )} . \end{split}$

      设在空间编码矩阵中的采样像素用“1”表示, DMD的总像素为M2. 将平均反射率$ \overline r $代替$ r(x, y) $. $ E $表示变量的数学期望, $ {\text{Var}} $表示变量的方差, $ {m_{\text{s}}} $表示平均信号光子计数值, $ {P_{B1}} $, $ {P_{B2}} $$ {P_{B3}} $分别表示$ {B_{i1}}(x, y, t) $, $ {B_{i2}}(t) $, $ d $三种噪声的功率. $m_{\rm{n}}^{(1)}$, $m_{\rm{n}}^{(2)}$$m_{\rm{n}}^{(3)}$分别表示三种噪声的平均光子计数值. 文献[22]给出平均光子计数值与激光功率之间的关系. $ h $是普朗克常数, $ v $是激光的光学频率, $ hv $是光子能量. 一个周期T内第$ i $个随机矩阵调制下的SNR模型表示为[22]

      $ {\rm{SNR}} = \dfrac{{E{{[{C_{i{\rm{s}}}}]}^2}}}{{{\text{Var}}[{C_{i{\rm{s}}}} + {C_{i{\text{n}}}}]}} =\dfrac{{{{(T{R_{\rm{a}}}{M_2}\overline r {m_{\rm{s}}})}^2}}}{{T{R_{\rm{a}}}{M_2}\overline r {m_{\rm{s}}} + {N^2}m_{\rm{n}}^{(1)}\displaystyle\sum\limits_{t = 0}^T {f(t + \tau )} + m_{\rm{n}}^{(2)}\displaystyle\sum\limits_{t = 0}^T {f(t + \tau )} + m_{\rm{n}}^{(3)}\displaystyle\sum\limits_{t = 0}^T {f(t + \tau )} }} . $

      假设单光子探测器工作在平稳状态下, 探测到光子后, 在一段时间内无法再次启动探测光子, 这段时间称为“死时间”. 死时间导致雪崩触发概率[22,23]发生改变, 影响光子探测总概率, 进而影响光子计数率. 根据团队前期工作, 将探测到的信号光子或噪声光子的总概率[23] ${\eta _{{\rm{s}},{\rm{ n}}}} = \dfrac{{\dfrac{1}{{1 + {d_{\text{t}}} \times {\varphi _{{\text{s}}, {\text{n}}}}(t)}}(1 - {{\rm{e}}^{ - {m_{{\text{s}}, {\text{n}}}}}})}}{{\eta {m_{\rm s, n}}}}$ 引入(10)式, 得到(11)式. 其中$ \dfrac{1}{{1 + {d_{\text{t}}} \times {\varphi _{{\text{s}}, {\text{n}}}}(t)}} $为雪崩触发概率, $ \eta $为量子效率, dt为探测器死时间, 信号光子流速率为$ {\varphi _{\text{s}}}(t) = R_{\text{a}} \times Sp \times {m_{\text{s}}} $, 噪声光子流速率为$ {\varphi _{\text{n}}}(t) = Sp \times {m_{\text{n}}} $, $ Sp $表示伪随机序列速率.

      $ {\rm{SNR}} = \dfrac{{E{{[{C_{i{\rm{s}}}}]}^2}}}{{{\rm{Var}}[{C_{i{\text{s}}}} + {C_{i{\text{n}}}}]}} = \dfrac{{{{({\eta _{\text{s}}}T{R_{\text{a}}}{M_2}\overline r {m_{\text{s}}})}^2}}}{{{\eta _{\text{s}}}T{R_{\text{a}}}{M_2}\overline r {m_{\text{s}}} + {\eta _{\text{n}}}{N^2} m^{(1)}_{\text{n}}\displaystyle\sum\limits_{t = 0}^T {f(t + \tau )} + {\eta _{\text{n}}} m^{(2)} _{\text{n}}\displaystyle\sum\limits_{t = 0}^T {f(t + \tau )} + {\eta _{\text{n}}} m^{(3)}_{\text{n}}\displaystyle\sum\limits_{t = 0}^T {f(t + \tau )} }} . $

      (11)式表明输出信噪比与序列长度成正比. 分母的后三项均包括$ f(t + \tau ) $, 如果序列中“1”码的概率几乎等于“–1”码的概率, 则噪声能量可忽略不计. 因此, 通过合理配置时域编码序列正负比特的平衡性, 可降低噪声能量. 其次, 考虑探测器的死时间, 引入探测到的信号光子或噪声光子的总概率, 死时间增大, 雪崩触发概率降低, 探测到信号光子的总概率降低, 输出信噪比降低.

    • 根据(6)式, 互相关函数$ {C_i}(\tau ) $中包含信号和噪声项, 其中, 第一项是信号项. 只有当目标平面和调制平面之间在空间上满足严格的一对一相关性, 以及接收到的光子到达时间点和发送的副本序列之间在时间上满足严格的一对一的相关性, 才能生成无噪声干扰的信号项. 噪声主要来自系统和环境, 这两类噪声不相关, 将两者相乘不会生成直流分量. 因此考虑通过二次相关来降低噪声. 假设矩阵${\boldsymbol{ \varPhi}} $的每一列表示为$ {\varphi _\textit{z}} \in {\mathbb{R}^{M \times 1}} $, $ M = 1, 2, \cdots , Q $, 则加性噪声 $ b $ 可以等效为随机矩阵$ {\boldsymbol{V}} $与信号$ p $的乘积, 表达为$ b = {\boldsymbol{V}}p $, 其中${\boldsymbol{ V }}$$ {\boldsymbol{\varPhi}} $非相关. 对(7)式做二次相关:

      $ \begin{split} &\varDelta (\tau) =\langle {C( \tau ),{\boldsymbol{\varPhi}} ( \tau )} \rangle - \langle {C( \tau )} \rangle \langle {{\boldsymbol{\varPhi}} ( \tau )} \rangle \\ =\;&\langle {{\boldsymbol{\varPhi}} p( \tau ) + {\boldsymbol{V}}p( \tau ),\Phi ( \tau )} \rangle - \langle {{\boldsymbol{\varPhi}} p( \tau ) + {\boldsymbol{V}}p( \tau )} \rangle \langle {{\boldsymbol{\varPhi}} ( \tau )} \rangle \\ =\;& \langle {{\boldsymbol{\varPhi}} p( \tau ),{\boldsymbol{\varPhi}}( \tau )} \rangle - \langle {{\boldsymbol{\varPhi}} p( \tau )} \rangle \langle {{\boldsymbol{\varPhi}} ( \tau )} \rangle \\ &+\langle {{\boldsymbol{V}}p( \tau ),{\boldsymbol{\varPhi}} ( \tau )} \rangle - \langle {{\boldsymbol{V}}p( \tau )} \rangle \langle {{\boldsymbol{\varPhi}} ( \tau )} \rangle \\ =\;&\dfrac{{M\displaystyle\sum\limits_{j = 1}^Q {( {{\varphi _\tau } \times {\varphi _j}} )} {p_j} - \displaystyle\sum\limits_{i = 1}^M {{\varphi _{i\tau }}\displaystyle\sum\limits_{i = 1}^M {\displaystyle\sum\limits_{j = 1}^Q {{\varphi _{ij}}{p_j}} } } }}{{{M^2}}} \\ &+\dfrac{{M\displaystyle\sum\limits_{j = 1}^Q {( {{\varphi _\tau } \times {v_j}} )} {p_j} - \displaystyle\sum\limits_{i = 1}^M {{\varphi _{i\tau }}\displaystyle\sum\limits_{i = 1}^M {\displaystyle\sum\limits_{j = 1}^Q {{v_{ij}}{p_j}} } } }}{{{M^2}}}. \end{split} $

      根据统计理论:

      $ {\varphi _\tau } \times {\varphi _j} = M\left[{\mu _\tau }{\mu _j} + {\rm cov} \left( {{\varphi _\tau },{\varphi _j}} \right)\right] \text{, } $

      $ \displaystyle\sum\limits_{i = 1}^M {{\varphi _{i\tau }}} \displaystyle\sum\limits_{i = 1}^M {{\varphi _{ij}} = {M^2}} {\mu _\tau }{\mu _j} \text{, } $

      $ {\varphi _\tau } \times {v_j} = M\left[{\mu _\tau }{\mu '_j} + {\rm cov} \left( {{\varphi _\tau },{v_j}} \right)\right] \text{, } $

      $ \displaystyle\sum\limits_{i = 1}^M {{\varphi _{i\tau }}} \displaystyle\sum\limits_{i = 1}^M {{v_{ij}} = {M^2}} {\mu _\tau }{\mu '_j} \text{, } $

      其中$ {\mu _\tau } $$ {\mu _j} $$ {\mu '_j} $分别是 $ {\varphi _\tau } $$ {\varphi _j} $${\varphi '}_j$的均值; ${\rm cov} (\cdot)$表示协方差运算.

      $ \varDelta \left( \tau \right) = \displaystyle\sum\limits_{j = 1}^Q {{\rm cov} \left( {{\varphi _\tau },{\varphi _j}} \right)} {p_j} + \displaystyle\sum\limits_{j = 1}^Q {{\rm cov} \left( {{\varphi _\tau },{v_j}} \right)} {p_j} . $

      $ {\boldsymbol{V}} $$ {\boldsymbol{\varPhi}} $非相关, 则协方差矩阵 ${\rm cov} \left( {{\varphi _\tau }, {v_j}} \right)$趋于0. (17)式表示为

      $ \varDelta = {\boldsymbol {Hp}} \text{, } $

      其中 ${\varDelta} \in {\mathbb{R}^{Q \times 1}}$, ${\boldsymbol{H}} = {\rm cov} ({\boldsymbol{\varPhi}} , {\boldsymbol{\varPhi}} )$. 压缩感知充分利用信号的可压缩性, 通过优化的方法恢复稀疏表示的测量值. 如果${\boldsymbol{ p}} $在稀疏基$ {\boldsymbol{\varPsi}} $的表示下是足够稀疏的, 则有${\boldsymbol{p}} = {\boldsymbol{\varPsi }}{ Z}$, 其中${ Z}$是稀疏系数. (18)式等效为$\varDelta = {{\boldsymbol{H}}}{\boldsymbol{\varPsi}} {\boldsymbol{Z}}{\text{ = }}{\boldsymbol{AZ}}$, 其中${\boldsymbol{A }}$为感知矩阵且${\boldsymbol{A}} = {{\boldsymbol{H}}}{\boldsymbol{\varPsi}}$, 采用以下凸优化反演方法重建图像$\hat{Z}$:

      $ \hat{ Z } = \arg \mathop {\min }\limits_Z \dfrac{{\text{1}}}{{\text{2}}}{\left\| {\varDelta - {\boldsymbol{A}}Z} \right\|^{\text{2}}} + \lambda \kappa (Z) \text{, } $

      其中$ \lambda $是正则化参数; $ \kappa (Z) $是正则项. 实验中, 采用TwIST[19]方法求解(19)式. 为了平衡计算的鲁棒性和复杂性, 引入最大似然估计(ML)来优化重建.

    • 为了验证本文提出成像方法的有效性, 采用MATLAB进行成像仿真实验. 由MATLAB生成10 GHz伪随机序列驱动激光器发射激光脉冲, 伪随机序列中1码的概率由随机数阈值确定, 设置为0.1, 长度为2048. 灰度图像“飞机”的深度范围为20—20.8 m, 如图2(a)所示. 将生成的M个64×64大小的高斯矩阵依次加载至DMD上, 接收端数字化时间编码和空间编码的信号后, 结合TwIST和ML法重建图像. 图2(b)(g)给出采用不同压缩比(5%—85%)进行重建的结果, 压缩比表示压缩的程度, 定义为$ \dfrac{测量次数}{总像素数} $, 压缩程度越高, 测量次数越少, 重建速度越快, 但重建质量越差. 采用不同压缩比探测和重建图像的总时间在8—78 s之间. 图2(e)(g)中45%或更大的压缩比可实现高分辨率和高对比度的重建. 由于ML的稳定性, 本方法的重建结果与参考文献[7,24]中强度图像重建工作所得出的结论一致, 即25%的压缩比也可清晰呈现目标场景, 如图2(d)所示. M序列和Gold序列是CDMA系统上行链路中最常用的伪随机序列[25,26]. 为了比较不同序列的抗噪声能力, 在仿真中引入泊松噪声. 平均噪声光子计数之和定义为三种噪声光子计数值$ m_{\rm{n}}^{(1)}$, $ m_{\rm{n}}^{(2)} $$ m_{\rm{n}}^{(3)} $之和, 记为${m_{\rm{n}}}$. 假设平均信号光子计数$ {m_{\text{s}}} $不变且$ {m_{\text{s}}} $=5. 噪声干扰由泊松随机数发生器生成, 将${m_{\rm{n}}} = 1$的噪声引入成像仿真, 由图3(a)图3(b), 1码概率为0.3的伪随机序列成像效果优于1码概率为0.1的序列成像效果. 由于m序列中1码的概率和–1码的概率基本相同, 图3(d)中所示的m序列的成像效果优于其他序列的成像效果, 与(11)式序列平衡度和成像信噪比的关系相互验证. 考虑单光子探测器死时间, 根据(11)式, 数值模拟死时间和成像信噪比的制约关系, 如图3(a)所示. 随着死时间的增大, 信号光子探测概率降低, 信噪比降低; 相同死时间下, 码流速度越大, 在死时间内到达探测器的信号光脉冲数量越多, 探测到的信号光子数目减少, 信噪比降低. 16384码长下1码概率为0.3, 引入不同死时间后的成像仿真结果如图3(e)(h)所示, 实验结果表明死时间减小, 信噪比增大, 成像质量提高, 变化规律和理论模型数值模拟结果一致. 定义深度图像重建质量的均方误差(mean squared error, MSE)为

      图  2  不同压缩比的时空域联合编码扩频单光子计数成像 (a) “飞机”深度真值; (b)−(g)不同压缩比的深度成像 (b) 5%, (c)15%, (d) 25%, (e) 45%, (f) 65%, (g) 85%. 颜色条表示深度数值, 单位为m

      Figure 2.  Time-space united coding spread spectrum single photon counting image with different compression proportions: (a) Depth ground truth of ‘airplane’; (b)−(g) depth maps with different compression proportions of (b) 5%, (c)15%, (d) 25%, (e) 45%, (f) 65% and (g) 85%. Colorbars show the depth with unit of meter.

      图  3  时空域联合编码扩频单光子计数成像性能模拟 (a) 1码概率为0.1的伪随机序列深度图像; (b) 1码概率为0.3的伪随机序列深度图像; (c) Gold序列的深度图像; (d) m序列的深度图像; (e) 死时间对成像性能的影响理论模拟; (f) 死时间为200 ns; (g) 死时间为20 ns; (h) 死时间为1ns. 颜色条表示深度数值, 单位为m

      Figure 3.  Time-space united coding image spread spectrum single photon counting imaging performance simulation: (a) Depth maps by pseudo-random sequences with ‘1’ bit of 0.1; (b) depth maps by pseudo-random sequences with ‘1’ bit of 0.3; (c) depth maps by Gold sequences; (d) depth maps by m-sequences; (e) theoretical simulation of dead time influence on imaging performance; (f) dead time is 200 ns; (g) dead time is 20 ns; (d) dead time is 1 ns. Colorbars show the depth information with unit of meter.

      $ {\text{MSE = }}\sqrt {\dfrac{1}{{{N^2}}}\displaystyle\sum\limits_{n = 1}^{{N^2}} {{{{\text{(}}{{\text{Z}}_n} - {G_n})}^2}} } . $

      设第n个像素重建的深度为$ {{\text{Z}}_n} $; $ {{\text{G}}_n} $为第n个像素的深度真值; N 2为像素数. 该项指标用于衡量成像准确程度, 误差越小, 确度越高.

      当噪声${m_{\rm{n}}} = 10$时, 图4(a)—(d)给出m序列长度分别为2048, 4096, 8192和16384重建的深度图. 随着序列长度的增加, 图像的信噪比增强, 与(11)式一致. 图4(j)给出码长1024到32768重建图像的MSE, 32768码长的MSE为0.033 m, 1024码长的MSE为0.1 m. 由于m序列良好的自相关性和ML的稳定性[18], 在环境噪声和暗计数噪声的干扰下, 仿真设定基于空间编码的单像素成像方法的积分时间与本文所提出的16384码长的积分时间相同, 图4(d)所示本文提出的成像方法优于图4(i)所示传统的成像方法, 也表明了大多数单像素激光雷达在远距离和户外成像质量不高的原因.

      图  4  不引入二次相关法的时空域联合编码扩频单光子计数深度成像, 码长为 (a) 2048, (b) 4096, (c) 8192, (d) 16384. 引入二次相关法的时空域联合编码扩频单光子计数深度图像, 码长为(e) 2048, (f) 4096, (g) 8192, (h)16384; (i) 传统的基于空间编码的单像素深度图像; (j) 引入二次相关法(点虚线)和不引入二次相关法(虚线)的深度MSE

      Figure 4.  Time-space united coding image spread spectrum single photon counting imaging without second correlated method by code length of (a) 2048, (b) 4096, (c) 8192, (d)16384. Depth maps simulations with second correlated method by code length of (e) 2048, (f) 4096, (g) 8192 and (h) 16384. (i) Traditional single pixel imaging method based on space coding. (j) The depth MSE with second correlated method (dot dashed line) and without second correlated method (dashed line).

      其次, 在本文提出的成像方法基础上, 引入二次相关重建方法, 重建的深度图4(e)—(h) MSE比不采用二次相关重建的深度图4(a)—(d) MSE小, 准确度高, 和图4(j)一致. 传统的基于空间编码的单像素成像的MSE为0.201 m, 引入二次相关法后成像的MSE为0.02 m. 本文提出的成像方法的MSE是传统成像方法的1/10,成像误差降低了9/10. 最后, 采用Sobel算子对图4(d)图4(h)图4(i)进行边缘检测, 结果分别如图5(a)图5(b)图5(c)所示. 引入二次相关方法后, 边缘检测的噪声鲁棒性更高, 边界更加平滑和清晰.

      图  5  成像边缘检测 (a) 不引入二次相关法的时空域联合编码扩频单光子成像边缘检测; (b) 引入二次相关法的时空域联合编码扩频单光子计数成像边缘检测; (c) 传统的基于空间编码的单像素成像边缘检测

      Figure 5.  Image edge detection: (a) Time-space united coding spread spectrum single photon counting imaging without second correlated method image edge detection; (b) time-space united coding spread spectrum single photon counting imaging with second correlated method image edge detection; (c) traditional single pixel imaging method based on space coding image edge detection.

      最后, 采用本课题组单光子计数扫描成像系统进行实际场景的模拟测试[27,28]. 噪声光子计数值为100 (counts/s), 单个像素积分时间为1 s, 成像系统的时间分辨率为8 ps, 图6(a)为采用该系统测量的高精度深度图, 将其作为真实深度图的模拟. 采用本文提出的方法进行目标探测和重建, 结果如图6(b)(h)所示. 由于所采用的系统深度分辨率比图像“飞机”高, 因此, 在相同条件下, 图6(b)图6(e)图6(h)中的MSE分别比图4(i)图4(d)图4(h)中的MSE小. m序列长度为16384, 当$ {m_{\rm{n}}} = 10 $时, 传统的基于空间编码的单像素深度成像如图6(b)所示. 在不引入二次相关法的情况下, 当$ {m_{\rm{n}}} = 5 $到10时, 时空域联合编码扩频单光子计数成像分别如图6(c)图6(d)图6(e)所示. 时空域联合编码扩频单光子计数成像性能优于基于空间编码的单像素成像性能. 基于空间编码的单像素成像的MSE为0.176 m, 而提出的时空域联合编码扩频单光子计数成像的MSE为0.035 m. 因此, 本文提出方法的MSE降低了4/5. 当$ {m_{\rm{n}}} = 5 $—10时, 采用二次相关法重建的深度图分别如图6(f)图6(g)图6(h)所示, 与不采用二次相关法的成像性能相比, 二次相关法重建的深度图MSE平均减少量约为0.018 m. 与图6(b)传统空间编码单像素成像方法相比, 引入二次相关法后MSE降低了9/10, 图7(b)的成像边缘更清晰, 对噪声的鲁棒性更强, 相比其他边缘检测结果更加平滑. 针对(19)式采用TVAL3, GPSR, FISTA重建的深度图如图8所示, 重建结果表明TVAL3的性能最好.

      图  6  单光子计数扫描成像测试图像深度重建仿真 (a) 玩偶深度真值; (b) $ {m_n} = 10 $, $ {\text{MSE = 0}}{\text{.176 m}} $, 传统的基于空间编码的单像素成像深度图; 在(c) ${m_{\rm{n}}} = 5$, $ {\text{MSE = 0}}{\text{.019 m}} $, (d) ${m_{\rm{n}}} = 8$, ${\text{MSE = 0}}{\text{.03\; m}}$, (e) $ {m_n} = 10 $, ${\text{MSE = 0}}{\text{.035 \;m}}$条件下, 不引入二次相关法的时空域联合编码扩频单光子计数深度图; 在(f) ${m_{\rm{n}}} = 5$, ${\text{MSE = 0}}{\text{.005\; m}}$, (g) ${m_{\rm{n}}} = 8$, ${\text{MSE = 0}}{\text{.011\; m}}$, (h) ${m_{\rm{n}}} = 10$, ${\text{MSE = 0}}{\text{.017 \;m}}$条件下, 引入二次相关法的时空域编码扩频单光子计数深度图, 深度单位为m

      Figure 6.  Depth reconstruction simulation by using single photon counting scanning imaging test image: (a) Ground truth of ‘doll’; (b) the depth map by traditional single pixel imaging method based on space coding image when ${m_{\rm{n}}} = 10$, ${\text{MSE = 0}}{\text{.176\; m}}$. Depth maps by time-space united coding spread spectrum single photon counting imaging without second correlated method when (c) ${m_{\rm{n}}} = 5$, ${\text{MSE = 0}}{\text{.019\; m}}$, (d) ${m_{\rm{n}}} = 8$, ${\text{MSE = 0}}{\text{.03\; m}}$ and (e) ${m_{\rm{n}}} = 10$, ${\text{MSE = 0}}{\text{.035\; m}}$. Depth maps by time-space united coding spread spectrum single photon counting imaging with second correlated method when (f) ${m_{\rm{n}}} = 5$, ${\text{MSE = 0}}{\text{.005\; m}}$, (g) ${m_{\rm{n}}} = 8$, ${\text{MSE = 0}}{\text{.011\; m}}$, and (h) ${m_{\rm{n}}} = 10$, ${\text{MSE = 0}}{\text{.017\; m}}$. The depth unit is m.

      图  7  成像边缘检测 (a) 不引入二次相关法的时空域联合编码扩频单光子成像边缘检测; (b) 引入二次相关法的时空域联合编码扩频单光子成像边缘检测; (c) 传统的基于空间编码的单像素成像边缘检测

      Figure 7.  Image edge detection: (a) Time-space united coding spread spectrum single photon counting imaging without second correlated method image edge detection; (b) time-space united coding spread spectrum single photon counting imaging with second correlated method image edge detection; (c) traditional single pixel imaging method based on space coding image edge detection.

      图  8  噪声$ {m_n} $=10不同的凸优化反演方法的图像重建 (a) TVAL3; (b) GPSR; (c) FISTA

      Figure 8.  Different convex optimization inversion methods imaging reconstruction when$ {m_n} $=10: (a) TVAL3; (b) GPSR; (c) FISTA.

      表1中计算了不同噪声水平下TVAL3, GPSR和FISTA的均方误差和重建时间, 以充分研究凸优化反演方法的性能. TwIST方法的MSE如图6所示. 噪声分别为${m_{\rm{n}}} = 5$${m_{\rm{n}}} = 10$, TwIST方法的运行时间分别为20 s和28 s. FISTA和GPSR反演的图像质量几乎与TwIST相同, FISTA和GPSR运行时间更短. TVAL3需要对参数进行更多调整, 虽然重建质量高, 但耗费时间长.

      MSE/m Time/s
      Method${m_{\rm{n} } } = 5$ ${m_{\rm{n}}} = 10$ ${m_{\rm{n}}} = 5$ ${m_{\rm{n}}} = 10$
      TVAL30.033 0.011 3557
      GPSR0.060 0.018 1324
      FISTA0.042 0.019 1624

      表 1  不同重建方法的性能统计

      Table 1.  Performance statistical record of different reconstruction methods.

    • 本文提出了一种结构紧凑、非扫描的时空域联合编码成像架构. 基于单光子探测理论, 建立了时空域联合相关的探测模型、成像正向模型和输出信噪比模型. 理论模型和实验结果证明, 通过配置时域编码的平衡度可优化成像质量. 将二次相关法引入到成像模拟中, 实验结果与理论模型相符, 均证明该技术具有探测距离远、抗干扰能力强、准确度高的优势. 今后的工作中, 课题组将构建时空域联合编码扩频单光子计数无扫描成像系统. 希望该项工作能够提高机器视觉、远程探测和无人驾驶等应用的成像性能, 为传统的单光子计数单像素激光雷达的研究打开新的局面.

参考文献 (28)

目录

    /

    返回文章
    返回