CN106683168A - 一种极低照度条件下光子计数激光三维计算成像方法 - Google Patents

Abstract

本发明公开了一种极低照度条件下光子计数激光三维计算成像方法，包括：同时测量得到所有像素的光子计数数据；构建光子数量与反射率的测量模型，并构建第一负对数似然函数；根据总变分半范数和第一负对数似然函数，构建第一目标函数；逐点迭代求解第一目标函数的凸优化问题，得到反射率图像；根据空间相关性构造二元假设检验统计量，筛选背景光子噪声；建立光子到达时间与测量距离的测量模型，并构建第二负对数似然函数；根据总变分半范数和第二负对数似然函数，构建第二目标函数；逐点迭代求解第二目标函数的凸优化问题，得到深度图像。通过使用本发明所提供的方法，可以在光子计数模式下同时获取高质量的激光三维距离像和强度像。

Description

一种极低照度条件下光子计数激光三维计算成像方法

技术领域

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

背景技术

现有技术中的激光三维成像探测技术既可实时获取目标的深度图像和激光反射率图像，从而同时生成目标场景的三维图像，又可实时自主测量相对姿态与运动参数，为视觉操控提供实时的三维测量与感知信息，是一种新型的自主三维视觉测量技术。现有技术中的激光三维成像遥感技术，可以直接快速生成目标场景的高分辨率三维影像或地形图，不需要地面控制点的配合，是一种新型定量遥感技术。

然而，现有技术中的激光三维成像探测和激光三维成像遥感技术等存在的共性难点问题在于：

1)远距离探测导致激光回波能量很低，激光探测器阵列接收到的光辐射只有光子量级照度；

2)需要在极低照度条件下同时生成深度图像和激光反射率图像，或三维立体测绘和光谱图像；

3)非扫描激光凝视成像，需要满足动态条件下可靠性、适装性要求。

现有技术中的激光显微断层成像仪存在的主要难点问题是：需要在极低照度条件下同时获取高质量的立体图像和激光反射率图像(反射率信息)，满足对活体组织结构动态、定量、三维的显微观测需求，也要求非扫描激光凝视成像。

目前，现有技术中已有的激光成像探测体制主要可以分为以下几类：

1)基于盖革模式雪崩光电二极管(Avalanche Photodiode in Geiger Mode,GM-APD)激光探测器阵列的激光凝视成像探测体制。该体制是解决此类激光成像探测的有效方法，但是，通过时间相关方法检测激光回波的飞行时间，只能产生三维距离像，而不能产生强度像；而且，当背景光不连续时，即使有光子计数过程中内在的泊松噪声，每个像素检测到几十个光子，也足以精确的距离成像，当由背景光产生的光子检测导致伪距离时，距离成像需要与反射成像相当的光子数量，需要几百个光子才能产生精确的强度像。

2)基于像增强CCD的激光强度成像体制仅接收激光回波的强度，不测量激光脉冲的飞行时间，只能产生强度像；在没有背景光照射的条件下，精确的反射成像需要几百个光子才能产生精确的强度像，当存在背景光照射时精确的反射成像需要的激光回波强度更大。

3)基于线性雪崩光电二极管(APD)激光探测器阵列激光凝视成像探测体制，既可以测量激光回波的飞行时间，还可测量激光回波的强度。但是，该体制强度成像需要连续的光子流，不能工作在光子计数模式，灵敏度很低，激光回波利用效率不高，不能满足极低照度条件下成像。

4)激光扫描成像体制，通过光机扫描或微扫描机构实现激光束逐点成像，从而获得距离像、强度像。但是，该技术要求激光重复频率很高，存在运动光机结构，平台适应性差，应用范围受限，因此不适用激光三维成像遥感，激光三维成像探测。

由此可知，现有技术中所要解决的技术问题是：在极低照度条件下、在很短驻留时间内，每个GM-APD像素只能检测到少量、不连续的光子，这些光子远不足以清晰成像；在没有背景光照射的条件下，清晰反射率成像需要几百个连续光子流，清晰距离成像则需要几十个光子流，因此难以同时生成激光三维距离像和反射率图像。

发明内容

有鉴于此，本发明提供一种极低照度条件下光子计数激光三维计算成像方法，从而可以在光子计数模式下同时获取高质量的激光三维距离像和强度像。

本发明的技术方案具体是这样实现的：

一种极低照度条件下光子计数激光三维计算成像方法，该方法包括如下步骤：

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

构建光子数量与反射率的测量模型，并构建给定光子数量条件下反射率的第一负对数似然函数；

根据总变分半范数和第一负对数似然函数，构建第一目标函数G_α(x,y)；

逐点迭代求解第一目标函数G_α(x,y)的凸优化问题，得到并输出反射率图像；

根据空间相关性构造二元假设检验统计量，筛选背景光子噪声；

建立光子到达时间与测量距离的测量模型，并构建光子达到时间条件下测量距离的第二负对数似然函数；

根据总变分半范数和第二负对数似然函数，构建第二目标函数J_R(x,y)；

逐点迭代求解第二目标函数J_R(x,y)的凸优化问题，得到并输出深度图像。

较佳的，所述光子计数数据包括：

每个像素接收到的光子数量和每个光子的到达时间。

较佳的，所述同时测量得到所有像素的光子计数数据包括：

设定激光探测器阵列的驻留时间T_w；在所述驻留时间内，使用泛光激光脉冲持续照射某一预设目标场景；

激光探测器阵列中的每个像素同时检测接收到的光子数量，并测量每个光子到达时间，从而得到所有像素的光子计数数据。

较佳的，所述第一负对数似然函数为：

其中，log(·)为自然对数，P₀(x,y)为没有检测到光子的概率，N为激光脉冲照射的次数，n为检测接收到的光子数量，α(x,y)为反射率，A表示平均光子数量，B表示背景光子平均到达速率。

较佳的，所述逐点迭代求解第一目标函数G_α(x,y)的凸优化问题，得到并输出反射率图像包括：

采用序列二阶Taylor级数逼近可分解的对数似然函数H(α)，迭代求解第一目标函数G_α(x,y)的最优解的最小值，即迭代求解如下公式的最小值：

其中，α_ML(x,y)是第一目标函数G_α(x,y)的最优解，β_α是惩罚系数，惩罚函数pl(α)是凸函数；H、W分别表示图像的行和列；

惩罚函数pl(α)定义为反射率{α(x,y)}的总变分半范数：

在k次迭代，H(α)的二阶Taylor级数逼近H^k(z):

其中，为负对数似然函数H(·)在z^k处的一阶导数，参数σ_k＞0，为向量z与迭代值z^k的二范数的平方；

k+1次迭代公式如下：

其中，z≥0，惩罚函数pl(α)为总变分半范数，β是惩罚系数；

凸目标函数G_α(x,y)的最优解简化为：

然后，通过如下的步骤得到反射率图像：

步骤A1，进行初始化，包括：

步骤A1a，设每个像素反射率初始值为

步骤A1b，对第一负对数似然函数求导，计算梯度

步骤A1c，选择参数σ_k∈[σ_min,σ_max]、梯度下降步长η＞0、最大迭代次数M∈Z₊；

步骤A2，开始进行迭代，k＝0,······,M；

步骤A3，迭代求解简化后的凸目标函数G_α(x,y)的最优解，得到α_k+1；

步骤A4，判断α_k+1是否满足迭代下降准则：

如果否，则α_k+1←ηα_k+1，然后再返回执行步骤A3；

如果是，则k←k+1，然后进行下一次算法迭代，直到算法迭代满足收敛准则：

其中，η＞0是很小的常数。

较佳的，所述根据空间相关性构造二元假设检验统计量包括如下步骤：

对于每一个像素位置(x,y)，根据八邻域像素检测到的光子到达时间的中值的平均值构造二元假设检验统计量H(x,y)；

根据二元假设检验统计量H(x,y)设置筛选光子准则。

较佳的，所述二元假设检验统计量H(x,y)通过如下的公式计算得到：

其中，八邻域像素检测到的光子到达时间分别为：

{t_j(x-1,y-1),j＝1,…,k₁}、{t_j(x-1,y),j＝1,…,k₂}、{t_j(x-1,y+1),j＝1,…,k₃}、

{t_j(x,y-1),j＝1,…,k₄}、{t_j(x,y+1),j＝1,…,k₅}、{t_j(x+1,y-1),j＝1,…,k₆}、

{t_j(x+1,y),j＝1,…,k₇}、{t_j(x+1,y+1),j＝1,…,k₈}；

对应到达时间中值分别为tm₁、tm₂、tm₃、tm₄、tm₅、tm₆、tm₇、tm₈；

k₁，k₂，k₃，k₄，k₅，k₆，k₇，k₈分别表示八邻域像素接收到的光子数量。

较佳的，所述筛选光子准则为：

其中，t_j(x,y)为光子的到达时间，H(x,y)为二元假设检验统计量，α_ML(x,y)是反射率的最大似然估计，T_p是激光脉冲宽度的均方根。

较佳的，所述第二负对数似然函数为：

其中，L[R(x,y)|t_j(x,y),j＝1,…,n(x,y)]为第二负对数似然函数，当n(x,y)＝0时，L[R(x,y)|t_j(x,y),j＝1,…,n(x,y)]＝0；g(·)为激光脉冲的波形。

较佳的，所述逐点迭代求解第二目标函数J_R(x,y)的凸优化问题，得到并输出深度图像包括如下步骤：

采用序列二阶Taylor级数逼近可分解的对数似然函数H(R)，迭代求解第二目标函数J_R(x,y)的最优解的最小值，即迭代求解如下公式的最小值：

其中，R_ML(x,y)是第二目标函数J_R(x,y)的最优解，β_R是惩罚系数，惩罚函数pl(R)是凸函数，H、W分别表示图像的行和列；

惩罚函数pl(R)定义为深度图像{R(x,y)}的总变分半范数：

在k次迭代，H(R)的二阶Taylor级数逼近H^k(z):

其中，为负对数似然函数H(·)在z^k处的一阶导数，参数σ_k＞0，为向量z与迭代值z^k的二范数的平方；

k+1次迭代公式如下：

其中，z≥0，惩罚函数pl(α)为总变分半范数，β是惩罚系数；

凸目标函数J_R(x,y)的最优解简化为：

然后，通过如下所述的步骤得到深度图像：

步骤B1，进行初始化，包括：

步骤B1a，经过筛选后的像素，初始深度值为其中，d为距离分辨率，δ为测量偏差；未经过筛选后的像素

步骤B1b，对第二负对数似然函数求导，计算梯度

步骤B1c，选择参数σ_k∈[σ_min,σ_max]、梯度下降步长η＞0、最大迭代次数M∈Z₊；

步骤B2，开始进行迭代，k＝0,······,M；

步骤B3，迭代求解简化后的凸目标函数J_R(x,y)的最优解，得到R_k+1；

步骤B4，判断R_k+1是否满足迭代下降准则：

如果否，则R_k+1←ηR_k+1，然后再返回执行步骤B3；

如果是，则k←k+1，然后进行下一次算法迭代，直到算法迭代满足收敛准则：

其中，η＞0是很小的常数。

如上可见，本发明所提供的极低照度条件下光子计数激光三维计算成像方法，可以精确测量在很短的驻留时间内每个GM-APD像素检测到的光子数量以及每个光子到达时间，对测量光子计数数据进行筛选，剔除噪声光子数据，根据经过筛选的光子数量和到达时间数据分别重构目标场景的反射率图像和三维距离像，实现了GM-APD探测器阵列同时生成三维距离像和反射率图像，从而可以在光子计数模式下同时获取高质量激光三维距离像和强度像。

附图说明

图1为本发明实施例中的极低照度条件下光子计数激光三维计算成像方法的流程示意图。

具体实施方式

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

本实施例提供了一种极低照度条件下光子计数激光三维计算成像方法。

图1为本发明实施例中的极低照度条件下光子计数激光三维计算成像方法的流程示意图。如图1所示，本发明实施例中的极低照度条件下光子计数激光三维计算成像方法主要包括如下所述的步骤：

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

在本步骤中，将同时测量所有像素的光子计数数据，从而获取驻留时间内的光子数据。

具体来说，较佳的，在本发明的具体实施例中，所述步骤11可以包括如下所述的步骤：

步骤111，设定激光探测器阵列的驻留时间T_w；在所述驻留时间内，使用泛光激光脉冲持续照射某一预设目标场景；

另外，较佳的，在本发明的具体实施例中，所述激光探测器阵列可以是：GM-APD激光探测器阵列。

步骤112，激光探测器阵列中的每个像素同时检测接收到的光子数量，并测量每个光子到达时间，从而得到所有像素的光子计数数据。

例如，较佳的，在本发明的具体实施例中，所述像素可以是：GM-APD激光探测器阵列中的GM-APD像素。

另外，在本发明的较佳实施例中，所述光子计数数据包括：每个像素接收到的光子数量和每个光子的到达时间，因此，在本发明的具体实施例中，可以设每个像素同时检测接收到的光子数量为：{n(x,y)}，每个光子的到达时间为：{t_j(x,y),j＝1,…,n(x,y)}；其中，t_j(x,y)是相对于最近一个激光脉冲的时间。

在上述的步骤111、112中，由于为激光探测器阵列设定了驻留时间T_w，相当于为激光探测器阵列中的每个像素设定了相同的驻留时间，从而可以保证面阵探测器中的所有像素可以同步工作，适用于面阵探测器的并行成像。

当然，并不是每个像素都能产生有效的光子计数数据，有的像素可能没有检测到任何光子，对于这类数据缺失，可以根据场景目标的空间相关性和反射率相关性来推断。

步骤12，构建光子数量与反射率的测量模型，并构建给定光子数量条件下反射率的第一负对数似然函数。

在本步骤中，将先构建光子数量n(x,y)与反射率α(x,y)的测量模型，从而建立了光子数量n(x,y)与反射率α(x,y)的数量关系；然后，即可根据测量模型构建给定光子数量{n(x,y)}条件下反射率α(x,y)的负对数似然函数L[α(x,y)|n(x,y)]，可称之为第一负对数似然函数。

在相同的驻留时间内，相同数量的激光脉冲照射目标时，反射率越高的目标返回的激光回波能量越大，单光子探测器检测到的光子数量越多；反射率越低的目标返回的激光回波能量越小，单光子探测器检测到的光子数量越小。因此，目标像素点的激光反射率α(x,y)可由检测到的光子数量n(x,y)进行编码，即每个像素的激光反射率α(x,y)是该像素检测到的光子数量n(x,y)的函数。

根据单光子激光探测器的统计特性可知，一次激光脉冲照射，光子计数探测器最多产生一次计数脉冲，其发生概率与背景光、信号光和反射率大小有关，且分布服从泊松分布模型：

没有检测到光子的概率：P₀(x,y)＝exp{-[a(x,y)·A+B]} (1)

检测到光子的概率：P₁(x,y)＝1-P₀(x,y) (2)

其中，A表示平均光子数量，B表示背景光子平均到达速率，α(x,y)为反射率。

对于N次激光脉冲照射，每个GM-APD像素检测到的光子数量n(x,y)服从二项分布，该分布是反射率α(x,y)的函数：

因此，在本发明的技术方案中，可以构建光子数量n(x,y)与反射率α(x,y)的测量模型。然后，构建给定光子数量{n(x,y)}条件下反射率α(x,y)的第一负对数似然函数L[α(x,y)|n(x,y)]，将反射率图像{α(x,y)}估计转换为逐点求解负对数似然函数与惩罚函数构成的目标函数G_α(x,y)的最优解。

较佳的，在本发明的具体实施例中，所述第一负对数似然函数可以表示为：

其中，log(·)为自然对数。去掉与α(x,y)无关的常数，L[α(x,y)|n(x,y)]是反射率α(x,y)的凸函数。

步骤13，根据总变分半范数和第一负对数似然函数，构建第一目标函数G_α(x,y)。

在本步骤中，可以由总变分半范数和所有像素位置的激光反射率{α(x,y)}的第一负对数似然函数L[α(x,y)|n(x,y)]构建第一目标函数G_α(x,y)。

单光子探测器阵列获取的光子测量数据均为非负数，第一负对数似然函数L[α(x,y)|n(x,y)]是光子数量n(x,y)的凸函数。因此，可以选择总变分半范数作为惩罚函数。第一目标函数G_α(x,y)也是光子数量n(x,y)的凸函数，则估计反射率图像即是逐点求解凸函数G_α(x,y)的最优解。

但是，测量数据的非负性在目标函数最优化中引入了一组不等式约束，导致最优化求解变得困难。因此，较佳的，在本发明的具体实施例中，针对非负约束的目标函数的凸优化问题，可以使用一种二阶Taylor级数逼近二次序列目标函数的迭代数值求解算法。

较佳的，在本发明的具体实施例中，可以考虑反射率估计平滑约束，在第一目标函数中选择总变分半范数作为惩罚函数，从所有像素检测到的光子数量{n(x,y)}估计反射率图像{α(x,y)}，即是逐点迭代求解负对数似然函数与惩罚函数构成的目标函数G_α(x,y)的凸优化问题的最优解。

例如，较佳的，在本发明的具体实施例中，可以通过如下所述的公式得到第一目标函数G_α(x,y)的最优解：

其中，α_ML(x,y)是最优解，β_α是惩罚系数，惩罚函数pl(α)是凸函数。约束反射率估计的非平滑性，H、W分别表示图像的行和列。

另外，惩罚函数pl(α)可以定义为反射率{α(x,y)}的总变分半范数：

步骤14，逐点迭代求解第一目标函数G_α(x,y)的凸优化问题，得到并输出反射率图像。

测量数据的非负性在目标函数最优化中引入了一组不等式约束，导致最优化求解变得困难。因此，较佳的，在本发明的具体实施例中，针对非负约束的目标函数的凸优化问题，可以使用一种二阶Taylor级数逼近二次序列目标函数的迭代数值求解算法。

例如，较佳的，在本发明的具体实施例中，所述步骤14可以通过如下所述的方法实现：

采用序列二阶Taylor级数逼近可分解的对数似然函数H(α)，迭代求解第一目标函数G_α(x,y)的最优解的最小值，即迭代求解公式(5)的最小值；

在k次迭代，H(α)的二阶Taylor级数逼近H^k(z):

其中，为负对数似然函数H(·)在z^k处的一阶导数，参数σ_k＞0，为向量z与迭代值z^k的二范数的平方；

k+1次迭代公式如下：

其中，z≥0，惩罚函数pl(α)为总变分半范数，β是惩罚系数；

凸目标函数G_α(x,y)的最优解简化为：

然后，通过如下所述的步骤得到反射率图像：

步骤141，进行初始化，包括：

步骤141a，设每个像素反射率初始值为

步骤141b，对第一负对数似然函数，即公式(4)求导，计算梯度

步骤141c，选择参数σ_k∈[σ_min,σ_max]、梯度下降步长η＞0、最大迭代次数M∈Z₊；

步骤142，开始进行迭代，k＝0,······,M；

步骤143，迭代求解简化后的凸目标函数G_α(x,y)的最优解，即迭代求解公式(17)、(18)，得到α_k+1；

步骤144，判断α_k+1是否满足迭代下降准则：

如果否(即α_k+1不满足迭代下降准则)，则α_k+1←ηα_k+1，然后再返回执行步骤143；

如果是，则k←k+1，然后进行下一次算法迭代，直到算法迭代满足收敛准则：

其中，η＞0是很小的常数。

通过上述的方法，即可逐点迭代求解第一目标函数G_α(x,y)的凸优化问题，得到反射率图像。

另外，较佳的，在本发明的具体实施例中，在得到反射率图像之后，即可输出该反射率图像{α(x,y)}，x＝1,…,W，j＝1,…,H。其中，H、W分别表示图像的行和列。

步骤15，根据空间相关性构造二元假设检验统计量，筛选背景光子噪声。

由于单光子探测器无法区分检测到的光子是来自目标返回的激光回波、背景光还是暗计数，而背景光产生的光子以及探测器本身的暗计数导致光子到达时刻的负对数似然函数呈现多模态，表现为非凸函数，是导致逐点最大似然估计失败的原因，因此检测到的光子到达时间t_j(x,y)必须在逐点似然估计距离之前进行背景噪声筛选，剔除噪声光子，即剔除来自背景光的光子和探测器本身的暗计数，才能用于后续的三维结构估计。

在物理场景中三维空间结构的高度空间相关性隐含信号光产生的光子到达时刻有非常小的条件方差，它们在时间上集中，在空间上相关；噪声光子的到达时刻是在时间区间[0,T_r]上均匀分布，在空间位置互相独立分布，有较高的方差。

因此，在本发明的技术方案中，可以根据空间相关性原理，构造二元假设检验统计量H(x,y)来筛选背景光子噪声。

例如，较佳的，在本发明的具体实施例中，所述根据空间相关性构造二元假设检验统计量可以包括如下所述的步骤：

步骤151，对于每一个像素位置(x,y)，根据八邻域像素检测到的光子到达时间的中值的平均值构造二元假设检验统计量H(x,y)；

例如，较佳的，在本发明的具体实施例中，所述二元假设检验统计量H(x,y)可通过如下所述的公式计算得到：

其中，八邻域像素检测到的光子到达时间分别为：

{t_j(x-1,y-1),j＝1,…,k₁}、{t_j(x-1,y),j＝1,…,k₂}、{t_j(x-1,y+1),j＝1,…,k₃}、

{t_j(x,y-1),j＝1,…,k₄}、{t_j(x,y+1),j＝1,…,k₅}、{t_j(x+1,y-1),j＝1,…,k₆}、

{t_j(x+1,y),j＝1,…,k₇}、{t_j(x+1,y+1),j＝1,…,k₈}。

对应到达时间中值分别为tm₁、tm₂、tm₃、tm₄、tm₅、tm₆、tm₇、tm₈；

k₁，k₂，k₃，k₄，k₅，k₆，k₇，k₈分别表示八邻域像素接收到的光子数量。

因此可知，上述的二元假设检验统计量H(x,y)为八个邻域像素光子到达时间的中值的平均值。

对于缺失光子数据的像素，H(x,y)＝0。

步骤152，根据二元假设检验统计量H(x,y)设置筛选光子准则；

例如，较佳的，在本发明的具体实施例中，所述筛选光子准则可以是：

其中，α_ML(x,y)是反射率的最大似然估计，T_p是激光脉冲宽度的均方根。

因此，满足上述公式(7)的光子属于信号光，否则，该光子属于背景光。

在低照度条件下，单光子激光探测器测量大光子计数数据包含大量的背景光子或暗计数，二元统计量H(x,y)是真实光子达到时间的较好近似值。单光子激光成像来说，α(x,y)A+B<<1，激光脉冲宽度在纳秒级，公式(8)右边项是很小的稳定值，不会导致二元统计量H(x,y)筛选光子准则性能降低。

步骤16，建立光子到达时间与测量距离的测量模型，并构建光子达到时间条件下测量距离的第二负对数似然函数。

在本步骤中，将建立光子到达时间t_j(x,y)与测量距离R(x,y)的测量模型，即建立光子到达时间t_j(x,y)与深度图像R(x,y)的数量关系；然后，即可根据测量模型构建光子到达时间t_j(x,y)条件下测量距离R(x,y)的负对数似然函数L[R(x,y)|t_j(x,y)]，可称之为第二负对数似然函数。

在真实世界中三维空间结构的高度空间相关性隐含信号光产生的光子到达时刻具有非常小的条件方差，因此基于空间相关性构造二元假设检验统计量，审查光子到达时刻tj(x,y)是否满足时间条件约束。

单光子探测器在一次激光脉冲持续期间只发生一次光子计数事件，只提供第一个检测到的光子的计时信息，此光子计数过程具有独立增量特性。因此，光子到达时间t_j(x,y)可用返回光信号的时间平移波形来表征，而平移时间由测量距离R(x,y)确定，从而可以建立光子到达时间t_j(x,y)与测量距离R(x,y)的测量模型。

例如，较佳的，在本发明的具体实施例中，对于经过筛选后的光子数据，在激光脉冲持续期间光子到达时间t(x,y)∈[0,Tr)的概率为：

其中，第一项中g(·)表示光信号的时间平移波形，系数为常数，第二项为常数。因此，光子到达时间t_j(x,y)的测量模型是测量距离R(x,y)的函数。

然后，可根据光子到达时间t_j(x,y)与测量距离R(x,y)的测量模型，构造光子达到时间t_j(x,y)与距离R(x,y)的第二负对数似然函数L[R(x,y)|t_j(x,y)]，将三维距离像{R(x,y))}估计转换为逐点求解由负对数似然函数与惩罚函数构成的目标函数J_R(x,y)的最优解。

例如，较佳的，在本发明的具体实施例中，可以根据驻留时间内光子到达时间的测量模型，给定光子到达时间数据{t_j(x,y),j＝1,…,n(x,y)}条件下测量距离R(x,y)的第二负对数似然函数L[R(x,y)|t_j(x,y),j＝1,…,n(x,y)]：

当n(x,y)＝0时，L[R|t_j,j＝1,…,n](x,y)＝0，对距离估计没有贡献。

激光脉冲的波形g(·)通常近似为：

g(t)∝exp[-v(t)] (11)

其中，t是激光脉冲持续时间，v(t)是激光脉冲波形的近似函数。

第二负对数似然函数L[R(x,y)|t_j(x,y)]可以记为H(R)。

另外，在本发明的技术方案中，上述的步骤16可以与步骤12同时执行，也可以按照预设执行顺序执行(例如，先执行步骤12再执行步骤16，或者相反)，本发明对此并不做限定。

步骤17，根据总变分半范数和第二负对数似然函数，构建第二目标函数J_R(x,y)。

单光子探测器阵列获取的光子测量数据均为非负数，第二负对数似然函数L[R(x,y)|t_j(x,y)]是光子到达时间t_j(x,y)的凸函数。因此，可以选择总变分半范数作为惩罚函数。第二目标函数J_R(x,y)也是光子到达时间t_j(x,y)的凸函数，则估计三维距离像即是逐点求解凸函数J_R(x,y)的最优解。

经过二元假设检验筛选后的光子数据，才可用于深度图像{R(x,y)}估计。由筛选后的像素位置的负对数似然函数，构建第二目标函数J_R(x,y)，考虑空间相关性约束，在目标函数中选择总变分半范数作为惩罚函数。从所有像素检测到的光子到达时间{t_j(x,y)}估计深度图像{R(x,y)}，即是逐点迭代求解负对数似然函数与惩罚函数构成的第二目标函数J_R(x,y)的凸优化问题的最优解。

例如，较佳的，在本发明的具体实施例中，可以通过如下所述的公式得到第二目标函数J_R(x,y)的最优解：

其中，R_ML(x,y)是最优解，β_R是惩罚系数，惩罚函数pl(R)是凸函数。约束空间距离估计的非平滑性，H、W分别表示图像的行和列。

另外，惩罚函数pl(R)可以定义为深度图像{R(x,y)}的总变分半范数：

步骤18，逐点迭代求解第二目标函数J_R(x,y)的凸优化问题，得到并输出三维距离像(即深度图像)。

测量数据的非负性在目标函数最优化中引入了一组不等式约束，导致最优化求解变得困难。因此，较佳的，在本发明的具体实施例中，针对非负约束的目标函数的凸优化问题，可以使用一种二阶Taylor级数逼近二次序列目标函数的迭代数值求解算法。具体过程与步骤14中的相关处理过程相类似。

例如，较佳的，在本发明的具体实施例中，所述步骤18可以通过如下所述的方法实现：

采用序列二阶Taylor级数逼近可分解的对数似然函数H(R)，迭代求解第二目标函数J_R(x,y)的最优解的最小值，即迭代求解公式(13)的最小值；

在k次迭代，H(R)的二阶Taylor级数逼近H^k(z):

其中，为负对数似然函数H(·)在z^k处的一阶导数，参数σ_k＞0，为向量z与迭代值z^k的二范数的平方；

k+1次迭代公式如下：

其中，z≥0，惩罚函数pl(α)为总变分半范数，β是惩罚系数；

凸目标函数J_R(x,y)的最优解简化为：

然后，通过如下所述的步骤得到深度图像：

步骤181，进行初始化，包括：

步骤181a，经过筛选后的像素，初始深度值为其中，d为距离分辨率，δ为测量偏差；未经过筛选后的像素

步骤181b，对第二负对数似然函数，即公式(12)求导，计算梯度

步骤181c，选择参数σ_k∈[σ_min,σ_max]、梯度下降步长η＞0、最大迭代次数M∈Z₊；

步骤182，开始进行迭代，k＝0,······,M；

步骤183，迭代求解简化后的凸目标函数J_R(x,y)的最优解，即迭代求解公式(17)、(18)，得到R_k+1；

步骤184，判断R_k+1是否满足迭代下降准则：

如果否(即R_k+1不满足迭代下降准则)，则R_k+1←ηR_k+1，然后再返回执行步骤183；

如果是，则k←k+1，然后进行下一次算法迭代，直到算法迭代满足收敛准则：

其中，η＞0是很小的常数。

综上可知，本发明所提供的极低照度条件下光子计数激光三维计算成像方法，可以精确测量在很短的驻留时间内每个GM-APD像素检测到的光子数量以及每个光子到达时间，对测量光子计数数据进行筛选，剔除噪声光子数据，根据经过筛选的光子数量和到达时间数据分别重构目标场景的反射率图像和三维距离像，实现了GM-APD探测器阵列同时生成三维距离像和反射率图像，从而可以在光子计数模式下(极低照度条件、很短驻留时间)同时获取高质量激光三维距离像和强度像。

该方法在激光成像探测、激光遥感和激光显微断层成像仪等领域有相当大的应用价值。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明保护的范围之内。

Claims (10)

1.一种极低照度条件下光子计数激光三维计算成像方法，其特征在于，该方法包括如下步骤：

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

构建光子数量与反射率的测量模型，并构建给定光子数量条件下反射率的第一负对数似然函数；

根据总变分半范数和第一负对数似然函数，构建第一目标函数G_α(x,y)；

逐点迭代求解第一目标函数G_α(x,y)的凸优化问题，得到并输出反射率图像；

根据空间相关性构造二元假设检验统计量，筛选背景光子噪声；

建立光子到达时间与测量距离的测量模型，并构建光子达到时间条件下测量距离的第二负对数似然函数；

根据总变分半范数和第二负对数似然函数，构建第二目标函数J_R(x,y)；

逐点迭代求解第二目标函数J_R(x,y)的凸优化问题，得到并输出深度图像。

2.根据权利要求1所述的方法，其特征在于，所述光子计数数据包括：

每个像素接收到的光子数量和每个光子的到达时间。

3.根据权利要求2所述的方法，其特征在于，所述同时测量得到所有像素的光子计数数据包括：

设定激光探测器阵列的驻留时间T_w；在所述驻留时间内，使用泛光激光脉冲持续照射某一预设目标场景；

激光探测器阵列中的每个像素同时检测接收到的光子数量，并测量每个光子到达时间，从而得到所有像素的光子计数数据。

4.根据权利要求1所述的方法，其特征在于，所述第一负对数似然函数为：

$\begin{array}{c}L\&lsqb;\mathrm{\&alpha;}|n\&rsqb;(x,y)=-\log C_n^N{P}_0{(x,y)}^{N-n}{\&lsqb;1-P_0(x,y)\&rsqb;}^n\\ =\left(N-n\right)a\left(x,y\right)\&CenterDot;A-n\log \left\{1-e^{-\&lsqb;a\left(x,y\right)\&CenterDot;A+B\&rsqb;}\right\}\end{array};$

其中，log(·)为自然对数，P₀(x,y)为没有检测到光子的概率，N为激光脉冲照射的次数，n为检测接收到的光子数量，α(x,y)为反射率，A表示平均光子数量，B表示背景光子平均到达速率。

5.根据权利要求4所述的方法，其特征在于，所述逐点迭代求解第一目标函数G_α(x,y)的凸优化问题，得到并输出反射率图像包括：

采用序列二阶Taylor级数逼近可分解的对数似然函数H(α)，迭代求解第一目标函数G_α(x,y)的最优解的最小值，即迭代求解如下公式的最小值：

${\mathrm{\&alpha;}}_{ML}(x,y)=\underset{a(x,y)\&GreaterEqual;0}{\mathrm{argmin}}(1-{\mathrm{\&beta;}}_{\mathrm{\&alpha;}})\underset{x=1}{\overset{W}{\&Sigma;}}\underset{y=1}{\overset{H}{\&Sigma;}}L\&lsqb;\mathrm{\&alpha;}|n\&rsqb;(x,y)+{\mathrm{\&beta;}}_{\mathrm{\&alpha;}}pl\left(\mathrm{\&alpha;}\right);$

其中，α_ML(x,y)是第一目标函数G_α(x,y)的最优解，β_a是惩罚系数，惩罚函数pl(α)是凸函数；H、W分别表示图像的行和列；

惩罚函数pl(α)定义为反射率{α(x,y)}的总变分半范数：

$pl\left(\mathrm{\&alpha;}\right)=\underset{x=1}{\overset{W-1}{\&Sigma;}}\underset{y=1}{\overset{H}{\&Sigma;}}|{\mathrm{\&alpha;}}_{x,y}-{\mathrm{\&alpha;}}_{x+1,y}|+\underset{x=1}{\overset{W}{\&Sigma;}}\underset{j=1}{\overset{H-1}{\&Sigma;}}|{\mathrm{\&alpha;}}_{x,y}-{\mathrm{\&alpha;}}_{x,y+1}|;$

在k次迭代，H(α)的二阶Taylor级数逼近H^k(z):

$H^k\left(z\right)=H\left(z^k\right)+{(z-z^k)}^T\&dtri;H\left(z^k\right)+\frac{{\mathrm{\&sigma;}}_k}{2}\left|\right|z-z^k|{|}_2^2;$

其中，为负对数似然函数H(·)在z^k处的一阶导数，参数σ_k＞0，为向量z与迭代值z^k的二范数的平方；

k+1次迭代公式如下：

$z^{k+1}=\underset{z\&Element;R^{W\&times;H}}{\arg \min }\{H^k\left(z\right)+\mathrm{\&beta;}\&CenterDot;pl\left(z\right)\};$

其中，z≥0，惩罚函数pl(α)为总变分半范数，β是惩罚系数；

凸目标函数G_α(x,y)的最优解简化为：

$z^{k+1}=\underset{z\&Element;R^{WH}}{\mathrm{argmin}}\{\mathrm{\&Phi;}\left(z\right)=||z-s^k|{|}_2^2/2+\mathrm{\&beta;}\&CenterDot;pl\left(z\right)/{\mathrm{\&sigma;}}_k\};$

$s^k=z^k-\&dtri;H\left(z^k\right)/{\mathrm{\&sigma;}}_k;$

然后，通过如下的步骤得到反射率图像：

步骤A1，进行初始化，包括：

步骤A1a，设每个像素反射率初始值为

步骤A1b，对第一负对数似然函数求导，计算梯度

$\&dtri;H\left({\mathrm{\&alpha;}}_0\right)(x,y)=(N-n)A+nA\frac{\exp \&lsqb;-a_0(x,y)\&CenterDot;A-B\&rsqb;}{1-\exp \&lsqb;-a_0(x,y)\&CenterDot;A-B\&rsqb;};$

步骤A1c，选择参数σ_k∈[σ_min,σ_max]、梯度下降步长η＞0、最大迭代次数M∈Z₊；

步骤A2，开始进行迭代，k＝0,……,M；

步骤A3，迭代求解简化后的凸目标函数G_α(x,y)的最优解，得到α_k+1；

步骤A4，判断α_k+1是否满足迭代下降准则：

$\mathrm{\&Phi;}\left({\mathrm{\&alpha;}}_{k+1}\right)\&le;\underset{1\&le;i\&le;k}{max}\{\mathrm{\&Phi;}\left({\mathrm{\&alpha;}}_{k+1}\right)-\frac{{\mathrm{\&epsiv;\&sigma;}}_k}{2}||{\mathrm{\&alpha;}}_{k+1}-{\mathrm{\&alpha;}}_k|{|}_2^2\};$

如果否，则α_k+1←ηα_k+1，然后再返回执行步骤A3；

如果是，则k←k+1，然后进行下一次算法迭代，直到算法迭代满足收敛准则：

$\left|\right|{\mathrm{\&alpha;}}_{k+1}-{\mathrm{\&alpha;}}_k|{|}_2^2/|\left|{\mathrm{\&alpha;}}_{k+1}\right|{|}_2^2\&le;\mathrm{\&eta;};$

其中，η＞0是很小的常数。

6.根据权利要求5所述的方法，其特征在于，所述根据空间相关性构造二元假设检验统计量包括如下步骤：

对于每一个像素位置(x,y)，根据八邻域像素检测到的光子到达时间的中值的平均值构造二元假设检验统计量H(x,y)；

根据二元假设检验统计量H(x,y)设置筛选光子准则。

7.根据权利要求6所述的方法，其特征在于，所述二元假设检验统计量H(x,y)通过如下的公式计算得到：

$H(x,y)=\frac{{\mathrm{tm}}_1+{\mathrm{tm}}_2+{\mathrm{tm}}_3+{\mathrm{tm}}_4++{\mathrm{tm}}_5++{\mathrm{tm}}_6++{\mathrm{tm}}_7++{\mathrm{tm}}_8}{8};$

其中，八邻域像素检测到的光子到达时间分别为：

{t_j(x-1,y-1),j＝1,…,k₁}、{t_j(x-1,y),j＝1,…,k₂}、{t_j(x-1,y+1),j＝1,…,k₃}、{t_j(x,y-1),j＝1,…,k₄}、{t_j(x,y+1),j＝1,…,k₅}、{t_j(x+1,y-1),j＝1,…,k₆}、{t_j(x+1,y),j＝1,…,k₇}、{t_j(x+1,y+1),j＝1,…,k₈}；

对应到达时间中值分别为tm₁、tm₂、tm₃、tm₄、tm₅、tm₆、tm₇、tm₈；

k₁，k₂，k₃，k₄，k₅，k₆，k₇，k₈分别表示八邻域像素接收到的光子数量。

8.根据权利要求7所述的方法，其特征在于，所述筛选光子准则为：

$|t_j(x,y)-H(x,y)|<2T_p\frac{B}{{\mathrm{\&alpha;}}_{ML}(x,y)S+B},j=1,...,n(x,y);$

其中，t_j(x,y)为光子的到达时间，H(x,y)为二元假设检验统计量，α_ML(x,y)是反射率的最大似然估计，T_p是激光脉冲宽度的均方根。

9.根据权利要求8所述的方法，其特征在于，所述第二负对数似然函数为：

$L\&lsqb;R|t_j,j=1,...,n\&rsqb;(x,y)=-\underset{j=1}{\overset{n(x,y)}{\mathrm{\&Sigma;}}}log\&lsqb;g\&lsqb;t_j(x,y)-2R(x,y)/c\&rsqb;\&rsqb;;$

其中，L[R(x,y)|t_j(x,y),j＝1,…,n(x,y)]为第二负对数似然函数，当n(x,y)＝0时，L[R(x,y)|t_j(x,y),j＝1,…,n(x,y)]＝0；g(·)为激光脉冲的波形。

10.根据权利要求9所述的方法，其特征在于，所述逐点迭代求解第二目标函数J_R(x,y)的凸优化问题，得到并输出深度图像包括如下步骤：

采用序列二阶Taylor级数逼近可分解的对数似然函数H(R)，迭代求解第二目标函数J_R(x,y)的最优解的最小值，即迭代求解如下公式的最小值：

$R_{ML}(x,y)=\underset{R(x,y)\&Element;\&lsqb;0,cTr/2)}{\mathrm{argmin}}(1-{\mathrm{\&beta;}}_R)\underset{x=1}{\overset{W}{\&Sigma;}}\underset{y=1}{\overset{H}{\&Sigma;}}L\&lsqb;R|t_j,j=1,...,n\&rsqb;(x,y)+{\mathrm{\&beta;}}_{R}pl\left(R\right);$

其中，R_ML(x,y)是第二目标函数J_R(x,y)的最优解，β_R是惩罚系数，惩罚函数pl(R)是凸函数，H、W分别表示图像的行和列；

惩罚函数pl(R)定义为深度图像{R(x,y)}的总变分半范数：

$pl\left(R\right)=\underset{x=1}{\overset{W-1}{\&Sigma;}}\underset{y=1}{\overset{H}{\&Sigma;}}|R_{x,y}-R_{x+1,y}|+\underset{x=1}{\overset{W}{\&Sigma;}}\underset{j=1}{\overset{H-1}{\&Sigma;}}|R_{x,y}-R_{x,y+1}|;$

在k次迭代，H(R)的二阶Taylor级数逼近H^k(z):

$H^k\left(z\right)=H\left(z^k\right)+{(z-z^k)}^T\&dtri;H\left(z^k\right)+\frac{{\mathrm{\&sigma;}}_k}{2}\left|\right|z-z^k|{|}_2^2;$

其中，为负对数似然函数H(·)在z^k处的一阶导数，参数σ_k＞0，为向量z与迭代值z^k的二范数的平方；

k+1次迭代公式如下：

$z^{k+1}=\underset{z\&Element;R^{W\&times;H}}{\arg \min }\{H^k\left(z\right)+\mathrm{\&beta;}\&CenterDot;pl\left(z\right)\};$

其中，z≥0，惩罚函数pl(α)为总变分半范数，β是惩罚系数；

凸目标函数J_R(x,y)的最优解简化为：

$z^{k+1}=\underset{z\&Element;R^{WH}}{\mathrm{argmin}}\{\mathrm{\&Phi;}\left(z\right)=||z-s^k|{|}_2^2/2+\mathrm{\&beta;}\&CenterDot;pl\left(z\right)/{\mathrm{\&sigma;}}_k\};$

$s^k=z^k-\&dtri;H\left(z^k\right)/{\mathrm{\&sigma;}}_k;$

然后，通过如下所述的步骤得到深度图像：

步骤B1，进行初始化，包括：

步骤B1a，经过筛选后的像素，初始深度值为其中，d为距离分辨率，δ为测量偏差；未经过筛选后的像素

步骤B1b，对第二负对数似然函数求导，计算梯度

$\&dtri;H\left(R_0\right)(x,y)=-\underset{j=1}{\overset{n(x,y)}{\mathrm{\&Sigma;}}}dv\&lsqb;t_j(x,y)-2R(x,y)/c\&rsqb;/{\mathrm{dt}}_j;$

步骤B1c，选择参数σ_k∈[σ_min,σ_max]、梯度下降步长η＞0、最大迭代次数M∈Z₊；

步骤B2，开始进行迭代，k＝0,……,M；

步骤B3，迭代求解简化后的凸目标函数J_R(x,y)的最优解，得到R_k+1；

步骤B4，判断R_k+1是否满足迭代下降准则：

$\mathrm{\&Phi;}\left(R_{k+1}\right)\&le;\underset{1\&le;i\&le;k}{max}\{\mathrm{\&Phi;}\left(R_{k+1}\right)-\frac{{\mathrm{\&epsiv;\&sigma;}}_k}{2}||R_{k+1}-R_k|{|}_2^2\};$

如果否，则R_k+1←ηR_k+1，然后再返回执行步骤B3；

如果是，则k←k+1，然后进行下一次算法迭代，直到算法迭代满足收敛准则：

$\left|\right|R_{k+1}-R_k|{|}_2^2/|\left|R_{k+1}\right|{|}_2^2\&le;\mathrm{\&eta;};$

其中，η＞0是很小的常数。