时空联合相关单光子计数三维成像方法

1.本发明属于数据处理领域，具体为一种时空联合相关单光子计数三维成像方法。


背景技术：

2.短积分时间和高噪声的探测背景下，远距离光子计数激光雷达系统通常记录大数据量但极其稀疏的数据立方体，信号光子被淹没在大量的噪声光子中，传统的图像估计方法得到的深度图与真实值偏差较大，难以准确获得探测目标信息，因此提出了以较少的内存要求、较低的计算复杂度和较快的成像速度，来估计光子缺乏情况下的准确深度图像的挑战。
3.交替方向乘子法，是一种求解优化问题的计算框架,适用于求解分布式的函数凸优化问题，特别是统计学习问题。交替方向乘子法通过分解协调过程，将大的全局问题分解为多个较小、较容易求解的局部子问题，并通过协调子问题的解而得到大的全局问题的解。
4.利用光子飞行时间原理的光子计数激光雷达相较于传统的激光雷达，能更加精确的获取探测目标的深度信息与反射率图像。在军事和民生领域，三维光子计数激光雷达的发展趋于成熟，应用范围愈加广泛。因此，在三维激光激光雷达系统研发的同时，研究人员提出了许多有效提升重建图像质量的方法，与硬件系统形成互补。
5.传统的最大似然估计(maximun likelihood estimation，mle)算法，通常用于数据量较大或者信噪比较高的情况，但在短积分时间或者大背景噪声的情况下，图像成像效果不佳。为了解决短积分时间导致的测量误差大和长积分时间导致的测量速度慢的矛盾，ahmed kirmani等人提出低光子计数情况下使用第一个探测到的光子进行成像的方法[kirmani a,venkatraman d,shin d,et al. first-photon imaging[j].science,2014,343(6166):58-61.]；dongeek shin提出适用于阵列单光子探测器的方法，利用横向像素和纵向像素的空间特性和时间特性结合的去噪方法[d shin,xu f,d venkatraman,et al.photon-efficient imagingwith a single-photon camera[j].nature communications,2016,7:12046.]。feng等人提出了时间相关光子快速去噪的方法[快速扫描光子计数激光雷达成像关键技术研究[d].南京理工大学,2018.].这些方法通过降低探测虚警率或提高图像重建正确率取得了较好的性能提升效果，但是在探测环境噪声影响较大的情况下，性能和效率会受到限制。


技术实现要素：

[0006]
本发明的目的在于提出了一种时空联合相关单光子计数三维成像方法。
[0007]
实现本发明目的的技术方案为：一种时空联合相关单光子计数三维成像方法，步骤如下：
[0008]
步骤1，利用远距离光子计数激光雷达系统，采集点云数据，并对图像数据进行预处理；
[0009]
步骤2，对预处理后的数据，进行自适应获取每个像素点的光子飞行时间，进一步
过滤噪声；
[0010]
步骤3，根据邻域信息补全空像素；
[0011]
步骤4，将步骤3得到的深度信息加入全变分正则化项构建成本函数，利用交替方向乘子法迭代更新得到最优解，作为深度信息估计值。
[0012]
优选地，对图像数据进行预处理包括：根据先验信息设置门控阈值，滤除图像中明显的噪声点。
[0013]
优选地，对预处理后的数据，进行自适应获取每个像素点的光子飞行时间，进一步过滤噪声的具体方法为：
[0014]
在像素点(i,j)，对光子的飞行时间进行排序，并将l个光子分组组成光子单元，对n-邻域中相邻光子飞行时间进行作差排序，求相邻光子单元偏差δi，找出偏差取最小值δi
min
时对应的光子单元：
[0015][0016]
式中，为在像素点(i,j)每个光子单元的飞行时间；
[0017]
由相邻光子单元偏差δi
min
对应的位置l的单元的飞行时间作为目标深度的初始估计深度t
initial
：
[0018][0019]
对光子飞行时间，判定下一组光子是否为信号光子，若下一组信号光子的飞行时间与初始估计深度之差小于阈值t
p
：
[0020][0021]
则称该类光子单元为信号光子单元，并将符合条件的光子飞行时间，记录在集合中，选定阈值k，如果像素点(i,j)处集合内的光子数达到标定值k，或者光子飞行时间与初始深度信息之差大于阈值t
p
范围，此时该像素点的判定结束，进入下一个像素点的光子判定。
[0022]
优选地，根据邻域信息补全空像素的具体方法为：对图像进行逐像素处理，在像素点(i,j)处，根据其周围8个点的邻域信息，均值填补空像素点。
[0023]
优选地，若像素(i,j)，最近周围8个像素点(i-1,j-1)、(i-1,j)、(i-1,j 1)、 (i,j-1)、(i,j 1)、(i 1,j-1)、(i 1,j)、(i 1,j 1)中标记的非空像素点不小于3，则记作小空像素点；反之，记作大空洞点，小像素点采用8-邻域非空像素均值填补，大空像素点，则视为空像素点。
[0024]
优选地，利用交替方向乘子法迭代更新得到最优解的具体方法为：
[0025]
在像素(i,j)，光子计数服从泊松分布，进行联合似然估计：
[0026][0027]
其中，s
i,j,t
＝r
i,j
g0(t-t
i,j
) b
i,j
d
i,j
，s
i,j,t
表示探测器接收到的所有光子信号， y
i,j,t
为激光雷达观测模型，t
i,j
≥0表示探测器和目标表面的距离，即深度信息， r
i,j
≥0表
示目标的强度，b
i,j
≥0表示背景和暗光子能级，g0表示系统的脉冲响应， t,r,b是n
×
1的向量，其中n＝nc×
nr，t是时间仓，通过最小化负对数似然函数得到成本函数：
[0028][0029]
其中，const表示不影响成本函数估计的常数，假设背景噪声b
i,j
＝0，仪器脉冲响应近似高斯分布其中c1是一个常数，σ是一个系统的超参数，并且是一个常数，根据最大似然估计法得到图像的深度和强度估计，和成本函数简化为：
[0030][0031]
其中，表示不影响成本函数估计的常数项，采用全变分正则化项来改进图像估计，根据邻域空间相关性可得全变分后的成本函数c
tv
(t,r)：
[0032]ctv
(t,r)＝l(t,r) τ1tv(t) τ2tv(r)
[0033]
其中τ1，τ2为正则化参数,其中和为像素n处的横坐标和纵坐标的一阶差分，为了最小化成本函数，通过交替方向乘子法最小化成本函数，交替乘子法的凸优化问题为：
[0034][0035]
其中，对深度和强度联合估计设为一个参数h
(j)
是一个选择矩阵，用来挑选联合估计中深度和强度的数据；
[0036]
在此假设一个满秩矩阵m为：
[0037][0038]
其中，j＝1,2,3,4，成本函数拆分为得：
[0039]h(1)
＝[0n,k]
[0040]h(2)
＝[k,0n]
[0041]
g3(u
(3)
＝τ1tv(z1) τ2tv(z2)，h
(3)
＝i
2n
[0042]h(4)
＝[in,0n]
[0043]
其中，g1、g2、g3、g4分别成本函数拆分后对应图像的强度、深度、tv 正则化、不必要的常数项，其中u
(3)
、u
(4)
分别为每个更新的解，r
iml
、为最似然估计方法计算得，in,kn为单位矩阵，分别对每一项进行交替方向乘子法迭代更新：
[0044][0045][0046][0047]
其中，s为求得最小值的变量值，k为迭代更新的次数，即当更新后的小于阈值threshold，跳出迭代循环，得到目标的深度估计。
[0048]
本发明与现有技术相比，其显著优点为：本发明利用探测数据间的时间相关性与被测目标的空间相关性来降低重构目标深度图像时的误差，根据信号光子和噪声光子在探测器响应在时域上呈现不同的特性，提取出所需的信号光子，分离噪声光子，再根据单光子的泊松概率分布模型结合目标空间相关性，通过提取的信号光子，重构目标深度图像；本发明成像效果优于传统最大似然估计方法和时间相关光子快速去噪方法，处理后得到的深度图像误差小，信噪比高，提高三维成像质量。
[0049]
下面结合附图对本发明作进一步详细描述。
附图说明
[0050]
图1为基于光子计数激光雷达的自适应图像去噪重建的方法总体流程示意图。
[0051]
图2为本发明中支撑实验方法的远距离光子计数激光雷达三维信息感知系统逻辑示意图。
[0052]
图3本发明中支撑实验方法的激光雷达系统示意图。
[0053]
图4为本发明利用时间相关性获取目标深度的流程示意图。
[0054]
图5为本发明深度信息加入全变分正则化项构成成本函数，利用交替方向乘子法迭代更新得到最优解的流程示意图。
[0055]
图6为用于验证本发明的一种时空联合相关单光子计数三维成像方法的实验结果，图6(a)为用于验证本发明的实验目标场景图，图6(b)为用于验证本发明的实验目标原始点云图，图6(c)为经过人工处理的深度图，作为不同方法比较的参考值，图6(d)为传统的最大似然估计算法的成像方法在采样积分时间为20ms时，处理得到的目标深度图像，图6(e)为feng提出的时间域去噪的成像方法在采样积分时间为20ms时，处理得到的目标深度图像，图6(f) 为本发明提出的成像方法在采样积分时间为20ms时，处理得到的目标深度图像。
具体实施方式
[0056]
下面将对本发明的一些技术细节进行进一步详细说明。下面结合附图对本发明作进一步详细描述。
[0057]
如图1所示，一种时空联合相关单光子计数三维成像方法，使用远距离光子计数激光雷达系统，采集一组点云数据；根据采集到点云数据，对图像数据进行预处理，应用一种开窗的方法，设置信号光子飞行时间门控，过滤大量噪声信息；根据预处理后的数据，进行自适应获取每个像素点的光子飞行时间，进一步过滤噪声；根据得到的光子飞行时间，根据邻域信息补全空像素；根据得到的深度信息加入全变分正则化项构建成本函数，利用交替方向乘子法迭代更新得到最优解，该解作为最终的深度信息估计值，具体步骤为：
[0058]
步骤1：通过如图2所示的激光雷达系统示意图，设置采集图像像素大小 nr×
nc，各像素积分时间等参数，采集到实验所需的点云数据nr×
nc×
t，对采集到的点云数据进行开窗处理。根据时间相关光子计数器探测的先验信息设置阈值门控，滤除掉图像中明显的噪点，极大的压缩了数据量，提升了数据处理效率。
[0059]
步骤2：根据步骤1初步去噪后得到的点云数据，进行下一步处理，结合图 4的流程图。在远距离探测的情况下，系统所接收的噪声光子量远远大于信号光子量，但信号光子比噪声光子更倾向于聚集在一起，由此，一种判定信号光子的方法是寻找光子最大的聚集簇。对数据进行逐像素处理，在像素点(i,j)，先对光子的飞行时间(time of flight，tof)进行排序，并将n个光子组成光子单元，假设该像素内的光子个数为l，对n-邻域中相邻光子tof进行作差排序，求相邻光子单元偏差δi，找出δi的最小值δi
min
时对应的光子单元：
[0060][0061]
式中，为在像素点(i,j)每个光子单元的飞行时间，由相邻光子单元偏差δi
min
对应的位置l的单元tof作为目标深度的初始估计深度t
initial
：
[0062][0063]
对光子飞行时间，判定下一组光子是否为信号光子，若下一组信号光子的飞行时间与初始估计深度之差小于阈值t
p
：
[0064][0065]
则称该类光子单元为信号光子单元，并将符合条件的光子飞行时间，记录在集合中，选定阈值k，如果像素点(i,j)处集合内的光子数达到标定值k，或者光子飞行时间与初始深度信息之差大于阈值t
p
范围，此时该像素点的判定结束，进入下一个像素点的光子判定。
[0066]
步骤3：在低积分时间的情况下，探测器获取的光子数减少。当使用非常低的采集时间时，传感器接收光子数急剧降低。在这些条件下，会产生很多空像素。如果没有额外的信息，就很难进行深度估计。根据步骤2所得到的数据，由于存在像素缺失的情况，根据空间相关性，像素点的领域信息填补空像素。若像素 (i,j)，最近周围8个像素点(i-1,j-1)、(i-1,j)、(i-1,j 1)、(i,j-1)、(i,j 1)、 (i 1,j-1)、(i 1,j)、(i 1,j 1)中标记的非空像素点
不小于3，则记作小空像素点；反之，记作大空洞点。对不同的空像素点匹配不同的插值模板进行填补，小像素点采用8-邻域非空像素均值填补，大空像素点，则视为空像素点。
[0067][0068]
步骤4：根据步骤3所得的tof数据，进行图像的深度信息迭代更新，结合图5的流程图，通过交替方向乘子法，得到目标较为准确的深度信息。
[0069]
激光雷达观测模型y
i,j,t
，其中(i,j)∈{1,
…
,nr}
×
{1,
…
,nc}，表示在像素(i,j)， t时间段内的光子计数，服从泊松分布：
[0070]yi,j,t
～poisson(s
i,j,t
)
[0071]
其中，s
i,j,t
＝r
i,j
g0(t-t
i,j
) b
i,j
d
i,j
，t
i,j
≥0表示探测器和目标表面的距离，表示深度信息，r
i,j
≥0表示目标的强度，b
i,j
≥0表示一个常数，表示背景和暗光子能级，g0表示系统的脉冲响应。假设观察像素之间的独立性导致联合似然：
[0072][0073]
其中t,r,b是n
×
1的向量，其中n＝nrnc，t是时间仓。
[0074]
通过最小化负对数似然函数得到：
[0075][0076]
其中，const表示不影响成本函数估计的常数，为简化成本函数，经典的估计方法设定没有背景噪声(b
i,j
＝0)，并且设定仪器脉冲响应近似高斯分布：其中c1是一个常数，σ是一个系统的超参数，并且是一个常数，根据mle得到，和在此设定的基础上，成本函数为(去除不必要的常数后)：
[0077][0078]
其中，表示不影响成本函数估计的常数项。
[0079]
在极低光子情况下，由于最大似然估计法估计得到和的质量很差，因此对估计图像的进行正则化图像恢复处理，通过考虑下面的成本函数来解释：
[0080]
c(t,r)＝l(t,r) φ(t,r)
[0081]
其中，φ(t,r)为正则项。进一步采用全变分正则化项来改进图像估计，根据邻域空间相关性，由此成本函数可得：
[0082]ctv
(t,r)＝l(t,r) τ1tv(t) τ2tv(r)
[0083]
其中τ1，τ2为正则化参数,其中和为像素n处的横坐标和纵坐标的一阶差分。
[0084]
通过交替方向乘子法去最小化成本函数，交替乘子法的凸优化问题为：
[0085][0086]
其中，对深度和强度联合估计设一个参数z，h
(j)
是一个选择矩阵，用来挑选联合估计中深度和强度的数据。假设一个满秩矩阵：
[0087][0088]
其中，j＝1,2,3,4，成本函数可拆分为得：
[0089]h(1)
＝[0n,k]
[0090]h(2)
＝[k,0n]
[0091]
g3(u
(3)
)＝τ1tv(z1) τ2tv(z2)，h
(3)
＝i
2n
[0092]h(4)
＝[in,0n]
[0093]
其中，g1、g2、g3、g4分别成本函数拆分后对应图像的强度、深度、tv 正则化、不必要的常数项，其中u
(3)
、u
(4)
分别为每个更新的解，r
iml
、为最似然估计方法计算得，in,kn为单位矩阵。分别对每一项进行交替方向乘子法迭代更新：
[0094][0095][0096][0097]
其中，s为求得最小值的变量值，k为迭代更新的次数。等于当更新后的小于阈值threshold，跳出迭代循环，得到目标的深度估计。