CN106772430A - 基于多分辨率小波逼近的单像素光子计数三维成像系统及方法 - Google Patents

Abstract

本发明公开了一种基于多分辨率小波逼近的单像素光子计数三维成像系统及方法，所述成像系统包括结构光投影系统、光子接收系统、同步控制及信号处理系统，采用了单像素光子探测器和数字微镜器件组合的单光子相机结构对目标场景进行成像。根据多分辨率小波逼近原理，从初始分辨率图像开始，逐渐获取目标场景由低分辨率到高分辨率的细节信息，用于重构最终分辨率三维图像。本发明有效减少了采样次数，缩短了成像时间，适用于高分辨率三维成像应用；同时，避免了CS算法所需的计算开销，减少了重构所需的时间；利用单像素光子探测器和DMD组合的单光子相机结构，减小了系统尺寸，简化了系统结构，具有结构简单、可靠性高、成本低的特点。

Description

基于多分辨率小波逼近的单像素光子计数三维成像系统及 方法

技术领域

本发明属于光子计数三维成像技术领域，特别是一种使用单像素探测器的光子计数三维成像系统，以及应用于该系统的基于多分辨率小波逼近的光子计数三维成像方法。

背景技术

光子计数三维成像具有探测灵敏度高、距离分辨率高的特点，可在极暗光照条件下获得目标场景反射率和三维结构信息，受到广泛关注。

在光子计数三维成像中，系统通过测量光子从发射到经目标场景反射被单光子探测器接收的飞行时间获得距离信息。目前，传统的光子计数三维成像系统可以达到亚毫米级的距离了分辨率，但是由于受限于其单点扫描成像机制，获得高空间分辨率需要很长的成像时间([1]卡塞格林式激光雷达结构，CN201320519712.[2]一种时间分辨光子计数成像系统及方法，CN201110152839.3)。

目前，为解决上述问题，已知有两种技术途径。一种方法仍然采用单点扫描成像机制，在成像所需测量次数不变的情况下，通过优化测量过程，减少单点的测量时间，从而减少成像时间。例如，Ahmed Kirmani等人提出的首光子成像方法([3]Ahmed kirmani,et al,First-Photon Imaging,Science,2014,vol.343,pp:58-61.)，通过建立单光子探测过程的概率统计模型，结合目标相邻像素的空间相关性，利用第一个探测到的回波光子信息，获取目标的三维信息。另一种方法使用单像素相机架构，其中点探测器选用具有单光子灵敏度的盖革模式雪崩二极管(Geiger-mode avalanche photoelectric diodes,GM-APD)或光电倍增管(photomultiplier tube,PMT)。这种方法利用目标场景强度图像的稀疏性，通过压缩感知(Compressed Sensing，CS)算法获取目标空间信息，极大地减少了成像所需测量次数。例如，Howland等人([4]G.A.Howland,P.B.Dixon,and J.C.Howell,Photon-countingcompressive sensing laser radar for 3D imaging,Appl.Opt.50,5917-5920,2011.)通过将CS算法和单像素探测器相结合获取目标场景空间信息，并使用距离门控获取距离信息。此外，中国科学院上海技术物理研究所的马彦鹏等人研制了基于压缩感知的激光雷达成像系统([5]一种基于压缩感知的激光雷达成像系统CN201410403458.1)。

第一种方法在目标场景的低反射区域，特别是边缘细节部分成像不准确，容易出现平滑现象。这两种技术途径均需要较复杂的后续算法处理，尤其是第二种基于CS的方法需要迭代计算解最优化问题重构图像，计算开销巨大，重构时间随着成像分辨率而成指数增长，成像速度缓慢，不适用于实时应用场景。

发明内容

本发明的目的在于提供一种快速、准确、高分辨率的基于多分辨率小波逼近的单像素光子计数三维成像系统及方法，解决传统光子计数三维成像方法受限于空间分辨率、成像时间长的问题。

实现本发明目的的技术解决方案为：一种基于多分辨率小波逼近的单像素光子计数三维成像系统及其方法，由结构光投影系统、光子接收系统、同步控制及信号处理系统构成三维成像系统，所述结构光投影系统包括激光光源、扩束镜、数字微镜器件DMD以及投影透镜；光子接收系统包括窄带滤光片、成像透镜以及光电倍增管PMT；同步控制及信号处理系统包括信号发生器、时间相关单光子模块TCSPC以及计算机PC；

所述信号发生器生成同步脉冲，驱动激光光源发出脉冲激光，经扩束镜扩束，照亮DMD微镜平面；根据PC生成的调制图案，DMD对脉冲激光进行空间调制，经投影透镜照射目标场景；目标场景反射光子经窄带滤光片滤除杂散光后，经成像透镜收集，由PMT测量得到光子探测信号；TCSPC将信号发生器输出的激光器发射信号与PMT光子探测信号相关联得到该脉冲的光子飞行时间TOF。

本发明与现有技术相比，其显著优点：(1)与传统光子计数三维成像方法相比，本发明采用的多分辨率小波逼近方法自适应采集目标场景重要信息，有效减少了测量次数，缩短了成像时间，适用于高分辨率三维成像应用。(2)与现有的基于CS的光子计数三维成像方法相比，本发明避免了CS算法所需的计算开销，减少了重构所需的时间。(3)本发明减小了系统尺寸，简化了系统结构，具有结构简单、可靠性高、成本低的特点。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是本发明基于多分辨率小波逼近的单像素光子计数三维成像系统示意图，包括超短脉冲激光光源(图1(a))、扩束镜(图1(b))、数字微镜器件(图1(c))、投影透镜的投影镜头(图1(d))、窄带滤光片(图1(e))、成像透镜(图1(f))、光电倍增管(图1(g))、信号发生器(图1(h))、时间相关单光子模块(图1(i))以及计算机(图1(j))。

图2是扩展小波树的结构示意图，数字64、128、256/512代表四个分辨率，箭头表示小波系数相关性。扩展小波树结构是一种小波分解提供的一种多分辨率的数据组织结构，不同分辨率对应相同空间区域的小波系数之间具有相关性。本发明所述基于多分辨率小波逼近的单像素光子计数三维成像方法利用了扩展小波树结构，由以获取的上一分辨率重要系数位置估计当前分辨率所需采样的重要系数位置，从而生成所需的DMD调制图案。

图3是本发明基于多分辨率小波逼近的单像素光子计数三维成像方法中某一分辨率的流程图，其中，图3(a)是结构光投影脉冲，图3(b)是目标场景反射光子，图3(c)是光子计数直方图，图3(d)是小波系数立方，图3(e)是图像立方，图3(f)是像素飞行时间直方图，图3(g)和图3(h)分别是当前分辨率深度图像和强度图像。图3描述了由已获取的上一分辨率深度图像小波系数估计当前重要小波系数位置，生成DMD调制图案，得到当前分辨率深度图像与强度图像，并估计下一分辨重要小波系数位置的多分辨率小波逼近过程。

图4是本发明基于多分辨率小波逼近的单像素光子计数三维成像系统及方法的成像结果，以及与文献[4]方法得到的结果的比较图。图4(a)是目标场景照片，图4(b)和图4(c)分别是采样率为5％时使用本发明所述的系统和方法得到的强度图像和深度图像，图4(d)和图4(e)是文献[4]方法得当的强度图像和深度图像。

具体实施方式

如图1所示，本发明基于多分辨率小波逼近的单像素光子计数三维成像系统由结构光投影系统、光子接收系统、同步控制及信号处理系统构成，所述结构光投影系统包括激光光源(图1(a))、扩束镜(图1(b))、数字微镜器件(digital micromirror device,DMD)(图1(c))以及投影透镜(图1(d))，其中脉冲激光经扩束后由DMD进行空间调制，实现结构光投影。光子接收系统包括窄带滤光片(图1(e))、成像透镜(图1(f))以及单像素光子计数型的光电倍增管(photomultiplier tube，PMT)(图1(g))，其中回波光子由PMT测量，得到光子探测信号，窄带滤光片的作用是阻挡杂散光，提高接收信噪比。同步控制及信号处理系统包括信号发生器(图1(h))、时间相关单光子模块(time-corrected single-photon counting，TCSPC)(图1(i))以及计算机(personal computer,PC)(图1(j))，其中信号发生器生成激光器驱动脉冲，TCSPC将激光器发射信号与PMT光子探测信号相关联得到光子飞行时间(time-of-flight，TOF)，PC将一段时间内的TOF综合得到光子计数直方图，用于三维图像重建。

本发明基于多分辨率小波逼近的单像素光子计数三维成像系统中，对于一次测量，信号发生器生成同步脉冲，驱动激光光源发出脉冲激光，经扩束镜扩束，照亮DMD微镜平面；根据PC生成的调制图案，DMD对脉冲激光进行空间调制，经投影透镜照射目标场景；目标场景反射光子经窄带滤光片滤除杂散光后，经成像透镜收集，由PMT测量得到光子探测信号；TCSPC将信号发生器输出的激光器发射信号与PMT光子探测信号相关联得到该脉冲的光子飞行时间(time-of-flight,TOF)。上述测量过程持续一段时间，称为积分时间，得到多个脉冲的TOF，PC将积分时间内的TOF累计得到本次测量的光子计数直方图。对目标场景的测量由多次测量组成，根据这些测量结果，重建目标场景三维信息。

本发明基于多分辨率小波逼近的单像素光子计数三维成像方法，从初始分辨率图像开始，利用多分辨率小波逼近逐渐获取目标场景由低分辨率到高分辨率的细节信息，用于重构最终分辨率三维图像，具体步骤如下：

第一步，DMD使用扫描调制图案，对目标场景进行逐点扫描。根据逐点扫描机制原理，扫描所需DMD调制图案个数与初始分辨率图像像素数相同，并且，每个DMD调制图案积分时间光子计数直方图记录了到达该像素的光子的飞行时间，该像素的深度值即为光子计数直方图最大值位置对应深度，强度值即为光子计数直方图的平均值，从而得到初始分辨率深度图像和强度图像。

第二步，由于深度图像比强度图像在小波域中更加稀疏，将上一分辨率的深度图像作为当前分辨率重要小波系数位置估计的依据，以减少所需采样次数。具体地，首先对上一分辨率深度图像进行单层小波分解，得到其小波系数矩阵。然后，根据扩展小波树结构(图2)中相邻分辨率对应相同空间位置的小波系数间的相关性，由上一分辨率深度小波系数矩阵中重要小波系数位置估计当前分辨率需要采集的重要小波系数的位置。

第二步中所述重要小波系数，指的是对应图像边缘、纹理等重要信息的绝对值较大的小波系数，只占自然图像小波系数矩阵中很少一部分。使用少量重要小波系数可以在不改变视觉质量的前提下重建图像。这个原理已广泛应用于图像和视频压缩算法，如JPEG2000。小波分解提供了一种分层的多分辨率图像数据表示结构，称为小波树结构。小波树相邻分辨率对应相同空间位置的小波系数具有相关性，如果某一低分辨率的小波系数是重要小波系数，其相邻高分辨率对应空间位置的小波系数也有很大可能是重要小波系数。戴慧东等人将同一分辨率小波系数相关性加入小波树，构成扩展小波树结构(具体可参加文献[6]，H.Dai,et al.Adaptive compressed sampling based on extended wavelettrees.2014,Applied Optics 53(29):6619-6628.)，如图2所示。

第三步，根据第二步估计结果，对当前分辨率每个重要小波系数生成两个对应小波基调制图案，并由DMD使用这些调制图案对脉冲激光进行调制，构成结构光(图3(a))，对目标场景进行探测。目标场景反射的光子由成像透镜收集，被PMT记录(图3(b))。TCSPC将信号发生器输出的激光器发射信号与PMT光子探测信号相关联得到该脉冲的光子飞行时间。对每个DMD调制图案积分时间内光子飞行时间进行累计，得到一组光子计数直方图对(图3(c))。将直方图时间轴均匀离散为多个时间点，并计算各个时间点的深度，对于每个深度，将计算得到的当前分辨率重要小波系数与上一分辨率小波系数相结合，得到该位置的小波系数逼近结果，再综合各个深度，得到小波系数立方(图3(d))，再通过分别对每个深度的小波系数矩阵做小波逆变换，得到图像立方(图3(e))，提取每个像素的光子飞行时间直方图(图3(f))，直方图最大值位置对应深度即为该像素的深度值，光子计数直方图的平均值即为该像素的强度值，从而得到当前分辨率深度图像(图3(g))和强度图像(图3(h))。

第三步中所述由重要小波系数位置生成DMD调制图案，构成结构光，以及由结构光测量结果重建深度图像和强度图像的方法如下：

以常见的Haar小波为例，任一Haar小波水平子带系数可由表示为：

式中，f表示最终分辨率图像，x₁,x₂是图像像素索引，j是小波尺度(决定了该系数的分辨率)，k＝(k₁,k₂)是该系数在系数矩阵中的位置，e表示三种小波子带系数类型：水平子带系数(e＝1)、垂直子带系数(e＝2)以及对角子带系数(e＝3)。公式(1)中，两个空间域积分对应两个小波基图案。因此，在多分辨率小波逼近中，任一小波系数可由两次结构光测量结果计算得到，其中DMD调制图案为构成该系数的小波基图案。公式(1)的结论可以推广到垂直子带系数和对角子带系数。

第四步，重复第二步和第三步，获取更高分辨率三维图像，直到获得最终分辨率三维图像。

实施例

本发明基于多分辨率小波逼近的单像素光子计数三维成像系统选用PicoQuant公司LDH系列830nm波长脉冲激光器作为结构光光源，具有4mW平均发射功率、10MHz脉冲重复频率以及300ps脉冲宽度。所用DMD是Vialux ALP 4.2系列，具有1024×768分辨率和4G片上内存，微镜最高翻转频率约为22.7kHz。本实施例中，积分时间设置为1ms。投影透镜的投影镜头和成像透镜分别选用50mm和35mm尼康标准镜头。单光子探测器选用滨松H7422P-50光子计数型PMT。窄带滤光片透过波长宽度为10nm。TCSPC为PicoQuant公司的PicoHarp300。

如图4(a)所示，目标场景是由字母“NJ”、“U”、“ST”硬纸板，以及黑色橡胶垫背景组成的，各个字母硬纸板尺寸如图4(a)中的标注所示。目标物放在距离本发明所述系统1.5～1.8m范围内。本次测量中，采样率设为全部像素(512×512)的5％。

首先，通过DMD使用扫描调制图案，获得64×64分辨率的初始分辨率深度图像和强度图像。

然后，根据扩展小波树结构，由该深度图像小波系数估计128×128分辨率重要小波系数位置，生成一组DMD调制图案，对目标场景进行测量。经目标场景反射的光子由成像透镜收集，被PMT记录，平均接收功率为0.125pW。TCSPC将发生脉冲与接收脉冲相关，记录每个DMD调制图案积分时间内的光子飞行时间，并综合成光子计数直方图。经深度离散，小波系数矩阵计数，小波逆变换得到空间分辨率为128×128的图像立方。对于每个空间像素，提取其距离轴信息，找出该像素的深度值和强度值，得到128×128分辨率的强度图像和深度图像。

接下来，利用多分辨率小波逼近，依次获取128×128、256×256分辨率三维图像，直至最终获取如图4(b)和图4(c)所示的512×512分辨率强度图像和深度图像。图4(d)和图4(e)是文献[4]所提出的方法在5％采样率下的重建结果。可以看出，本发明得到的强度图像和深度图像的质量均优于文献[4]等现有技术所提出的方法。此外，本发明方法三维图像重建用时17.5秒，而文献[4]的方法重建用时超过了20分钟。

Claims (4)

1.一种基于多分辨率小波逼近的单像素光子计数三维成像系统，其特征在于由结构光投影系统、光子接收系统、同步控制及信号处理系统构成三维成像系统，所述结构光投影系统包括激光光源、扩束镜、数字微镜器件DMD以及投影透镜；光子接收系统包括窄带滤光片、成像透镜以及光电倍增管PMT；同步控制及信号处理系统包括信号发生器、时间相关单光子模块TCSPC以及计算机PC；

所述信号发生器生成同步脉冲，驱动激光光源发出脉冲激光，经扩束镜扩束，照亮DMD微镜平面；根据PC生成的调制图案，DMD对脉冲激光进行空间调制，经投影透镜照射目标场景；目标场景反射光子经窄带滤光片滤除杂散光后，经成像透镜收集，由PMT测量得到光子探测信号；TCSPC将信号发生器输出的激光器发射信号与PMT光子探测信号相关联得到该脉冲的光子飞行时间TOF。

2.一种利用权利要求1所述的基于多分辨率小波逼近的单像素光子计数三维成像系统进行成像方法，其特征在于过程如下：

对于一次测量，信号发生器生成同步脉冲，驱动激光光源发出脉冲激光，经扩束镜扩束，照亮DMD微镜平面；

根据PC生成的调制图案，DMD对脉冲激光进行空间调制，经投影透镜照射目标场景；

目标场景反射光子经窄带滤光片滤除杂散光后，经成像透镜收集，由PMT测量得到光子探测信号；

TCSPC将信号发生器输出的激光器发射信号与PMT光子探测信号相关联得到该脉冲的光子飞行时间TOF；

上述测量过程持续一段时间构成积分时间，得到多个脉冲的TOF，PC将积分时间内的TOF累计得到本次测量的光子计数直方图，对目标场景的测量由多次测量组成，根据这些测量结果，重建目标场景三维信息。

3.根据权利要求2所述的成像方法，其特征在于具体步骤如下：

第一步，DMD使用扫描调制图案，对目标场景进行逐点扫描，根据逐点扫描机制原理，扫描所需DMD调制图案个数与初始分辨率图像像素数相同，并且每个DMD调制图案积分时间光子计数直方图记录了到达该像素的光子的飞行时间，该像素的深度值即为光子计数直方图最大值位置对应深度，强度值即为光子计数直方图的平均值，从而得到初始分辨率深度图像和强度图像；

第二步，将上一分辨率的深度图像作为当前分辨率重要小波系数位置估计的依据，即首先对上一分辨率深度图像进行单层小波分解，得到其小波系数矩阵；然后，根据扩展小波树结构中相邻分辨率对应相同空间位置的小波系数间的相关性，由上一分辨率深度小波系数矩阵中重要小波系数位置估计当前分辨率需要采集的重要小波系数的位置；

第三步，根据第二步估计结果，对当前分辨率每个重要小波系数生成两个对应小波基调制图案，并由DMD使用这些调制图案对脉冲激光进行调制，构成结构光，对目标场景进行探测，目标场景反射的光子由成像透镜收集，被PMT记录，TCSPC将信号发生器输出的激光器发射信号与PMT光子探测信号相关联得到该脉冲的光子飞行时间，并对每个DMD调制图案积分时间内光子飞行时间进行累计，得到一组光子计数直方图对，将直方图时间轴均匀离散为多个时间点，并计算各个时间点的深度，对于每个深度，将计算得到的当前分辨率重要小波系数与上一分辨率小波系数相结合，得到该位置的小波系数逼近结果，再综合各个深度，得到小波系数立方，再通过分别对每个深度的小波系数矩阵做小波逆变换，得到图像立方，提取每个像素的光子飞行时间直方图，直方图最大值位置对应深度即为该像素的深度值，光子计数直方图的平均值即为该像素的强度值，从而得到当前分辨率深度图像和强度图像；

第四步，重复第二步和第三步，获取更高分辨率三维图像，直到获得最终分辨率三维图像。

4.根据权利要求3所述的成像方法，其特征在于第三步中所述由重要小波系数位置生成DMD调制图案，构成结构光，以及由结构光测量结果重建深度图像和强度图像的方法如下：

以常见的Haar小波为例，任一Haar小波水平子带系数可由表示为：

$w_{j,k}^e=<f,{\mathrm{\&psi;}}_{j,k}^1>=2^{-j}\&lsqb;{\&Integral;}_{2^j{k}_1}^{2^j(k_1+1)}{\&Integral;}_{2^j{k}_2}^{2^j(k_2+1/2)}f(x_1,x_2){\mathrm{dx}}_1{\mathrm{dx}}_2-{\&Integral;}_{2^j{k}_1}^{2^j(k_1+1)}{\&Integral;}_{2^j(k_2+1/2)}^{2^j(k_2+1)}f(x_1,x_2){\mathrm{dx}}_1{\mathrm{dx}}_2\&rsqb;$

式中，f表示最终分辨率图像，x₁,x₂是图像像素索引，j是小波尺度，k＝(k₁,k₂)是该系数在系数矩阵中的位置，e表示三种小波子带系数类型：e＝1是指水平子带系数、e＝2是指垂直子带系数以及e＝3是指对角子带系数，公式中，两个空间域积分对应两个小波基图案；因此，在多分辨率小波逼近中，任一小波系数可由两次结构光测量结果计算得到，其中DMD调制图案为构成该系数的小波基图案。