CN109613556A - 基于稀疏表征的光子计数激光三维探测成像方法 - Google Patents

Abstract

本发明提供一种基于稀疏表征的光子计数激光三维探测成像方法，包括：同时测量得到所有像素的光子计数数据；根据光子计数数据，分别构建光子数量与反射率的模型以及光子到达时间与深度的测量模型；同时利用时空信息，添加全变分正则化约束，构建TV/L2最小化稀疏正则化信号重建模型；通过拉格朗日方法处理优化TV/L2最小化稀疏正则化信号重建模型，使用交替方向方法来迭代地找到子问题的解；得到并输出反射率和深度图像。本发明方法可以在极低照度条件下，对光子数据进行去噪处理，基于稀疏表征的思想实现高光子效率，在光子计数模式下获取高质量的反射率图像和深度图像，重构目标三维图像。

Description

基于稀疏表征的光子计数激光三维探测成像方法

技术领域

本发明涉及激光三维计算成像技术，特别涉及一种低照度条件下光子计数激光三维探测成像方法。

背景技术

激光三维成像技术是一项利用照射激光在目标探测所在场景中反射回来的光脉冲，从中提取强度和范围信息的技术，在目标识别与辨认等方向都有着巨大的前景。与传统成像技术相比，具有高角度分辨一种。它所具有的高探测灵敏度和时间分辨率特性，使得在诸多领域如荧光成像、荧光相关光率、高距离分辨率、小发散角的特点，可应用于地形测绘、构建虚拟环境、数字城市、目标识别与辨认等多方面领域。而时间相关单光子计数技术是基于统计采样的原理，在记录低强度、高重复频率的脉冲信号时通过统计大量的高重频信号的处理方法来减小随机误差，属于单光子探测技术的谱等都有着重要意义。将时间相关单光子计数技术应用于激光三维成像系统中，能大大的提高激光探测的时间分辨精度和探测灵敏度，减少激光脉冲能量，快速实现对目标物体的三维成像。基于时间相关单光子技术的激光三维探测技术是近年来激光三维探测的主要发展方向之一。稀疏表征是通过将图像基于过完备字典进行稀疏分解，得到相应的稀疏表示，而同时噪声在字典上没有对应的稀疏表示，从而将噪声与图像中的有用信息分离，在图像恢复领域具有巨大的潜力。在成像过程中，往往会受到各式噪声的影响，使得最终图像不能达到预期的要求，不能很好地反映目标的真实形象轮廓，如何改善图像质量，抑制噪声对其的不利影响，后续处理工作显得十分重要，基于稀疏表征的方法，对获得的反映目标信息建立模型，重建更为清晰标准的三维图像。将稀疏表征应用于光子计数激光三维探测技术中，通过对多个光子构建稀疏正则化的表示方法来实现图像去噪，抑制背景噪声给成像过程带来的影响，已成为近年的研究热门。

在传统的激光探测与测量技术中，形成精准的三维图像，每个像素则需要采集成百上千个光子。然而，很多科学研究例如远距离遥感成像、生物显微镜成像等领域内，都是在弱光环境下进行的，光通量都非常小，亦或是探测远程目标时，只能探测到有限的光子，加上背景噪声的影响，成像质量受到严重限制。如何在限制光子数的情况下，能快速捕捉到高质量的反映目标特征的图像，成为该领域研究的重点问题。

发明内容

本发明针对现有技术的不足，提供一种基于稀疏表征的光子计数激光三维探测方法，从而可以在低照度条件下以光子计数模式同时获取高质量的激光三维距离像和强度像。

本发明的技术方案为基于稀疏表征的光子计数激光三维探测方法，该方法包括以下步骤：

步骤S1，同时测量得到所有像素的光子计数数据；

步骤S2，根据光子计数数据，分别构建光子数量与反射率的模型以及光子到达时间与深度的测量模型；

步骤S3，同时利用时空信息，添加全变分正则化约束，构建TV/L2最小化稀疏正则化信号重建模型；

步骤S4，通过拉格朗日方法处理优化TV/L2最小化稀疏正则化信号重建模型，使用交替方向方法来迭代地找到子问题的解；

步骤S5，得到并输出反射率和深度图像。

进一步的，步骤S1中光子计数数据包括：每个像素接收到的光子数量k_i,j以及每个光子的到达时间其中l指光子探测数，取值从1到k_i,j。

进一步的，步骤S2中光子数量与反射率的模型为：

其中，η为探测器量子效率，N为脉冲个数，k_i,j表示每个像素(i,j)接收到的光子数量，为每个脉冲重复周期的总信号，T_r为设定的重复周期，s(t)表示单个脉冲的波形。

进一步的，步骤S2中光子到达时间与深度z的测量模型为：

其中，c为光速，t_i,j为每个光子的到达时间。

进一步的，步骤S3中构建的TV/L2最小化信号重建模型为：

其中，g为通过成像系统测量得到的光子计数数据通过测量模型估计后得到的数据，即被噪声影响待优化的信号，大小为ROWS×COLS×N，f即为所求经算法处理后的结果，大小同为ROWS×COLS×N，μ为正则化参数，可调节TV正则化对整个优化过程的影响；

所添加的TV正则化为各向异性，即表示为

其中，算子D_x，D_y，D_t分别是沿水平，垂直方向以及时间方向的正向有限差分算子，(β_x,β_y,β_t)是常数，[f]_i表示向量f的第i个分量；将运算符D定义为三个子运算符的集合其中每个子运算符的定义为：

D_xf＝vec(f(x+1,y，t)-f(x,y,t))

D_yf＝vec(f(x,y+1，t)-f(x,y,t))

D_tf＝vec(f(x,y，t+1)-f(x,y,t))

为了在控制沿每个方向的前向差异方面具有更大的灵活性，引入三个比例因子，定义标量β_x，β_y，β_t，并分别将它们与D_x，D_y，D_t相乘，使得：

通过调整β_x，β_y，β_t控制对个别项D_xf，D_yf，D_tf的相对强调。

进一步的，步骤S4的具体实现方式如下，

将核心优化TV/L2最小问题写出其拉格朗日形式为：

其中，y是与约束相关的拉格朗日乘数，其初始值为预先设定，ρ_r是指与后面相乘的二次惩罚项相关联的正则化参数；

通过交替方向方法迭代求解以下子问题：

(1)f-子问题

对上述拉格朗日形式进行对f求导，得出方程：

(μI+ρ_rD^TD)f＝μg+ρ_rD^Tu-D^Ty

I是单位矩阵，由此可得解为：

其中，F表示3D傅里叶变换算子；

(2)u-子问题

该子问题可通过收缩公式解决，定义：

v_x＝β_xD_xf+(1/ρ_r)y_x

v_y＝β_yD_yf+(1/ρ_r)y_y

v_t＝β_tD_tf+(1/ρ_r)y_t

由此得出u：

(3)y-子问题

y_k+1＝y_k-ρ_r(u_k+1-Df_k+1)

其中k为迭代次数。

进一步的，步骤S5的具体实现方式为，

步骤S51，计算并输出反射率图像，根据光子数量与反射率α的模型，选择合适的参数μ，β_x，β_y，β_t，设置初始值f₀＝g_α，u₀＝Df₀，y＝0，k＝0，然后开始进行迭代求解TV/L2信号重建模型最小化的最优解，得到f_k+1，使k＝k+1，直到满足收敛准则：

||f_k+1-f_k||₂/||f_k||₂≤tol

其中tol为指定误差；

对算法输出优化后的N幅图像组成的视频f_α(ROWS×COLS×N)取中值，最终输出结果反射率图像f_α(ROWS×COLS)；

步骤S52，背景噪声审查，背景计数不包含任何场景深度信息，它们的检测时间在空间位置上相互独立，方差为T_r ²/12；相反，由于光脉冲的持续时间为T_p<<T_r，深度在空间位置上相关，给定来自相邻位置的数据，信号计数的检测时间具有条件方差，其远低于T_r ²/12；基于该关键观察，在每个像素(i，j)处对8个相邻像素的检测时间进行排序，计算其中值为(若相邻像素都没有探测值，设)，对光子进行筛选：

其中，T_p为脉冲宽度的均方根；

步骤S53，计算并输出深度图像，根据审查后的结果以及光子到达时间与深度z的测量模型，选择合适的参数μ，β_x，β_y，β_t，设置初始值f₀＝g_z，u₀＝Df₀，y＝0，k＝0，然后开始进行迭代求解TV/L2信号重建模型最小化的最优解，得到f_k+1，使k＝k+1，直到满足收敛准则：

||f_k+1-f_k||₂/||f_k||₂≤tol

其中tol为指定误差，

对算法输出优化后的N幅图像组成的视频f_z(ROWS×COLS×N)取中值，最终输出结果深度图像f_z(ROWS×COLS)。

由上可见，本发明所提供的一种基于稀疏表征的光子计数激光三维探测方法，可以在极低照度条件下，对光子数据进行去噪处理，在过去只考虑空间相关性的基础上进行改进，同时利用时空信息，基于稀疏表征的思想实现高光子效率，在光子计数模式下获取高质量的反射率图像和深度图像，重构目标三维图像。

附图说明

图1是本发明实施例的光子计数激光三维探测方法的流程图。

图2是本发明实施例的TV/L2算法流程图。

图3是本发明实施例中对scene_man_flower_N100使用TV/L2算法重建的反射率图像。

图4是本发明实施例中对scene_man_flower_N100使用TV/L2算法重建的深度图像。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下参照附图并举实施例，对本发明进一步详细说明。

本实施例提供了一种在极低照度条件下基于稀疏表征的光子计数激光三维探测方法。

如图1所示，本发明实施例中的基于稀疏表征的光子计数激光三维探测方法主要包括如下几个步骤：

步骤1，同时测量得到所有像素的光子计数数据；

一个完整的激光三维成像系统由激光光源、光学收发扫描系统、单光子探测器、时间相关单光子技术模块以及数据处理与系统控制模块这几个部分组成，通过光栅扫描的方式使用周期性脉冲激光照亮场景，由单光子探测器探测反射信号，通过计数器进入数据处理模块，重建目标3D场景。我们假设用一系列N个脉冲照射每个像素(i,j)，为模拟来自环境光和暗计数的噪声，设在工作光学波长λ处的光通量为b_λ的背景光照射在探测器上。

具体来说，在本发明的具体实施例中，所采用的信号采集模型为以光栅扫描的方式使用周期性脉冲激光照亮场景，由单光子探测器探测反射信号，通过计数器进入数据处理模块，重建目标三维场景。

另外，在本发明实施例中，所述光子计数数据包括：每个像素接收到的光子数量k_i,j以及每个光子的到达时间其中l指光子探测数，取值从1到k_i,j。

步骤2，分别构建光子数量与反射率的模型以及光子到达时间与深度的测量模型；S21：构建光子数量与反射率α的模型。我们用周期性脉冲激光以光栅扫描的方式照亮场景，设定重复周期T_r，单个脉冲的波形由s(t)表示，c为光速。光源照射目标反射光经探测器检测的光通量为：

r_i,j(t)＝α_i,js(t-2z_i,j/c)+b_λ

由于低通量条件，忽略SPAD复位时间，由SPAD响应s(t)传输的反射光而产生的光子检测构成具有随时间变化率函数ηr_i,j(t)的非均匀泊松过程，其中η为探测器量子效率。对于这些光子检测，我们添加探测器器暗计数，它来自于速率为d的独立齐次泊松过程。将暗计数与背景生成的光子探测集合在一起产生SPAD输出，即当仅传输单个脉冲时，生成具有速率函数的非均匀泊松过程：

λ_i,j(t)＝ηr_i,j(t)+d

＝ηα_i,js(t-2z_i,j/c)+(ηb_λ+d)

这意味着在时间间隔中到达数量的概率质量函数(Probability Mass Function,PMF)为

k是指时间间隔中到达的光子数量，其中

定义B＝(ηb_λ+d)T_r作为每个脉冲重复周期的总信号和背景计数，并将其用于所有后续计算中的背景计数包括暗计数以及环境光产生的计数。我们假定B是已知的，在低光通量成像系统中，每个脉冲重复周期每个像素的光通量远小于1，即ηα_i,jS+B≤1。所以这种每脉冲重复周期的低光通量需要在许多脉冲重复周期内成像。忽略探测器复位时间，通过之前给出的速率函数，响应于单个照射脉冲检测到的光子数量的PMF由上式给出，其中一个重复周期内平均探测器次数：

而零探测的概率为：

由于这些概率对于每个脉冲是独立的，所以在N个脉冲测到的光子的数量k_i,j是以概率质量函数二项分布的：

其中，为组合数：

给定在像素(i，j)处的光子探测总数k_i,j和受约束的反射率估计是：

传统的，归一化的光子计数值用作反射率的估计值：

其中，S是定义的作为每个脉冲重复周期的总信号；N指脉冲的个数。

S22：构建光子到达时间与深度z的测量模型。该步骤与S21相类似，深度z的估计值为：

步骤3，同时利用时空信息，添加全变分正则化约束，构建TV/L2最小化稀疏正则化信号重建模型：

由于背景噪声的影响，我们同时利用单光子输出数据的空时信息，对其添加稀疏正则化约束，构建TV/L2最小化信号重建模型：

其中，g为通过成像系统测量得到的光子计数数据通过测量模型估计后得到的数据(g 可以是反射率，也可以是深度信息)，大小为ROWS×COLS×N，f即为待求解变量(g是反射率数据时，f为反射率对应的重建值；g是深度信息数据时，f就是对应的深度信息重建值)，大小同为ROWS×COLS×N。μ为正则化参数，可调节TV正则化对整个优化过程的影响。

所添加的TV正则化为各向异性，即表示为，

其中算子D_x,D_y，D_t分别是沿水平，垂直方向以及时间方向的正向有限差分算子。这里， (β_x,β_y,β_t)是常数，[f]_i表示向量f的第i个分量。我们将运算符D定义为三个子运算符的集合其中每个子运算符的定义为：

D_xf＝vec(f(x+1,y，t)-f(x,y,t))

D_yf＝vec(f(x,y+1，t)-f(x,y,t))

D_tf＝vec(f(x,y，t+1)-f(x,y,t))

为了在控制沿每个方向的前向差异方面具有更大的灵活性，我们引入三个比例因子。我们定义标量β_x，β_y，β_t，并分别将它们与D_x，D_y，D_t相乘，使得：

我们通过调整β_x，β_y，β_t，我们可以控制对个别项D_xf，D_yf，D_tf的相对强调。可知||f||_TV1等于Df上的矢量1范数，即||f||_TV1＝||Df||₁。因此，为了简单起见，我们使用||Df||₁来代替。

步骤4，通过拉格朗日方法处理优化，使用交替方向方法来迭代地找到子问题的解，得到深度图像；

分析信号重建模型TV/L2最小化，我们通过上述步骤3将其可以从无约束最小化问题转化为等价约束最小化问题。

当利用多光子数据时，我们将其视作由N帧图像组成的一段视频，这里N指每个像素取 N个光子的数值，即从时间上对N次测量光子数据进行优化，则核心优化TV/L2最小问题为：

subject to u＝Df

写出其拉格朗日形式为：

其中，y是与约束相关的拉格朗日乘数，其初始值为预先设定，ρ_r是指与后面相乘的二次惩罚项相关联的正则化参数。

S42：通过交替方向方法迭代求解以下子问题：

(1)f-子问题

对上述拉格朗日形式进行对f求导，得出方程：

(μI+ρ_rD^TD)f＝μg+ρ_rD^Tu-D^Ty

I是单位矩阵，由此可得解为：

其中，F表示3D傅里叶变换算子。

(2)u-子问题

该子问题可通过收缩公式解决，我们定义：

v_x＝β_xD_xf+(1/ρ_r)y_x

v_y＝β_yD_yf+(1/ρ_r)y_y

v_t＝β_tD_tf+(1/ρ_r)y_t

由此得出u：

(3)y-子问题

y_k+1＝y_k-ρ_r(u_k+1-Df_k+1)。

其中k为迭代次数。

步骤5：得到并输出反射率和深度图像。

S51：计算并输出反射率图像。根据光子数量与反射率α的模型，选择合适的参数μ，β_x，β_y，β_t，设置初始值f₀＝g_α(由探测得到的光子数量数据，经过光子数量与反射率α的模型得到的反射率估计数据)，u₀＝Df₀，y＝0，k＝0，然后开始进行迭代求解TV/L2信号重建模型最小化的最优解，得到f_k+1，使k＝k+1，直到满足收敛准则：

||f_k+1-f_k||₂/||f_k||₂≤tol

其中tol为指定误差。

对算法输出优化后的N幅图像组成的视频f_α(ROWS×COLS×N)取中值，最终输出结果反射率图像f_α(ROWS×COLS)。

S52：背景噪声审查。背景计数不包含任何场景深度信息，它们的检测时间在空间位置上相互独立，方差为T_r ²/12。相反，由于光脉冲的持续时间为T_p<<T_r，深度在空间位置上相关，给定来自相邻位置的数据，信号计数的检测时间具有条件方差，其远低于T_r ²/12。基于该关键观察，我们在每个像素(i，j)处对8个相邻像素的检测时间进行排序，计算其中值为(若相邻像素都没有探测值，设)，对光子进行筛选：

其中，T_p为脉冲宽度的均方根。

S53：计算并输出深度图像。根据审查后的结果以及光子到达时间与深度z的测量模型，选择合适的参数μ，β_x，β_y，β_t，设置初始值f₀＝g_z(由探测得到的光子到达时间数据，经过光子到达时间与深度z的测量模型得到的深度估计数据)，u₀＝Df₀，y＝0，k＝0，然后开始进行迭代求解TV/L2信号重建模型最小化的最优解，得到f_k+1，使k＝k+1，直到满足收敛准则：

||f_k+1-f_k||₂/||f_k||₂≤tol

其中tol为指定误差。

对算法输出优化后的N幅图像组成的视频f_z(ROWS×COLS×N)取中值，最终输出结果深度图像f_z(ROWS×COLS)。

具体的，在该实施例中，两次设置参数μ为0.2，β_x＝1，β_y＝1，β_t＝15，以及指定误差tol＝1e-3。

另外，这里使用的实测数据是由Dongeek Shin和Jeffrey H.Shapiro等人在2015年发表的Single-Photon Depth Imaging Using a Union-of-Subspaces Model中在低光照条件下进行光子深度成像中的实验数据。

对该实验数据进行TV/L2重建算法处理，得到结果如图3，图4。从图3和图4中可以看出，本发明提供的基于稀疏表征的光子计数激光三维探测方法，能够有效抑制在弱光环境下的背景噪声，在光子计数模式下短时间内实现对目标的高质量三维成像。

本文中所描述的具体实施例仅仅是对本发明作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (7)

1.一种基于稀疏表征的光子计数激光三维探测成像方法，其特征在于，包含以下步骤：

步骤S1，同时测量得到所有像素的光子计数数据；

步骤S2，根据光子计数数据，分别构建光子数量与反射率的模型以及光子到达时间与深度的测量模型；

步骤S3，同时利用时空信息，添加全变分正则化约束，构建TV/L2最小化稀疏正则化信号重建模型；

步骤S4，通过拉格朗日方法处理优化TV/L2最小化稀疏正则化信号重建模型，使用交替方向方法来迭代地找到子问题的解；

步骤S5，得到并输出反射率和深度图像。

2.如权利要求1所述的一种基于稀疏表征的光子计数激光三维探测成像方法，其特征在于：步骤S1中光子计数数据包括：每个像素接收到的光子数量k_i,j以及每个光子的到达时间其中l指光子探测数，取值从1到k_i,j。

3.根据权利要求2所述的一种基于稀疏表征的光子计数激光三维探测成像方法，其特征在于：步骤S2中光子数量与反射率的模型为：

其中，η为探测器量子效率，N为脉冲个数，k_i,j表示每个像素(i,j)接收到的光子数量，为每个脉冲重复周期的总信号，T_r为设定的重复周期，s(t)表示单个脉冲的波形。

4.根据权利要求2所述的一种基于稀疏表征的光子计数激光三维探测成像方法，其特征在于：步骤S2中光子到达时间与深度z的测量模型为：

其中，c为光速，t_i,j为每个光子的到达时间。

5.根据权利要求1所述的一种基于稀疏表征的光子计数激光三维探测成像方法，其特征在于：步骤S3中构建的TV/L2最小化信号重建模型为：

其中，g为通过成像系统测量得到的光子计数数据通过测量模型估计后得到的数据，大小为ROWS×COLS×N，f即为所求经算法处理后的结果，大小同为ROWS×COLS×N，μ为正则化参数，可调节TV正则化对整个优化过程的影响；

所添加的TV正则化为各向异性，即表示为

其中，算子D_x，D_y，D_t分别是沿水平，垂直方向以及时间方向的正向有限差分算子，(β_x,β_y,β_t)是常数，[f]_i表示向量f的第i个分量；将运算符D定义为三个子运算符的集合其中每个子运算符的定义为：

D_xf＝vec(f(x+1,y，t)-f(x,y,t))

D_yf＝vec(f(x,y+1，t)-f(x,y,t))

D_tf＝vec(f(x,y，t+1)-f(x,y,t))

为了在控制沿每个方向的前向差异方面具有更大的灵活性，引入三个比例因子，定义标量β_x，β_y，β_t，并分别将它们与D_x，D_y，D_t相乘，使得：

通过调整β_x，β_y，β_t控制对个别项D_xf，D_yf，D_tf的相对强调。

6.根据权利要求5所述一种基于稀疏表征的光子计数激光三维探测成像方法，其特征在于：步骤S4的具体实现方式如下，

将核心优化TV/L2最小问题写出其拉格朗日形式为：

其中，y是与约束相关的拉格朗日乘数，其初始值为预先设定，ρ_r是指与后面相乘的二次惩罚项相关联的正则化参数；

通过交替方向方法迭代求解以下子问题：

(1)f-子问题

对上述拉格朗日形式进行对f求导，得出方程：

(μI+ρ_rD^TD)f＝μg+ρ_rD^Tu-D^Ty

I是单位矩阵，由此可得解为：

其中，F表示3D傅里叶变换算子；

(2)u-子问题

该子问题可通过收缩公式解决，定义：

v_x＝β_xD_xf+(1/ρ_r)y_x

v_y＝β_yD_yf+(1/ρ_r)y_y

v_t＝β_tD_tf+(1/ρ_r)y_t

由此得出u：

(3)y-子问题

y_k+1＝y_k-ρ_r(u_k+1-Df_k+1)

其中k为迭代次数。

7.根据权利要求3和4和6所述的一种基于稀疏表征的光子计数激光三维探测成像方法，其特征在于：步骤S5的具体实现方式为，

步骤S51，计算并输出反射率图像，根据光子数量与反射率α的模型，选择合适的参数μ，β_x，β_y，β_t，设置初始值f₀＝g_α，u₀＝Df₀，y＝0，k＝0，然后开始进行迭代求解TV/L2信号重建模型最小化的最优解，得到f_k+1，使k＝k+1，直到满足收敛准则：

||f_k+1-f_k||₂/||f_k||₂≤tol

其中tol为指定误差；

对算法输出优化后的N幅图像组成的视频f_α(ROWS×COLS×N)取中值，最终输出结果反射率图像f_α(ROWS×COLS)；

步骤S52，背景噪声审查，背景计数不包含任何场景深度信息，它们的检测时间在空间位置上相互独立，方差为T_r ²/12；相反，由于光脉冲的持续时间为T_p<<T_r，深度在空间位置上相关，给定来自相邻位置的数据，信号计数的检测时间具有条件方差，其远低于T_r ²/12；基于该关键观察，在每个像素(i，j)处对8个相邻像素的检测时间进行排序，计算其中值为(若相邻像素都没有探测值，设)，对光子进行筛选：

其中，T_p为脉冲宽度的均方根；

步骤S53，计算并输出深度图像，根据审查后的结果以及光子到达时间与深度z的测量模型，选择合适的参数μ，β_x，β_y，β_t，设置初始值f₀＝g_z，u₀＝Df₀，y＝0，k＝0，然后开始进行迭代求解TV/L2信号重建模型最小化的最优解，得到f_k+1，使k＝k+1，直到满足收敛准则：

||f_k+1-f_k||₂/||f_k||₂≤tol

其中tol为指定误差，

对算法输出优化后的N幅图像组成的视频f_z(ROWS×COLS×N)取中值，最终输出结果深度图像f_z(ROWS×COLS)。