CN113406594B - 一种基于双量估计法的单光子激光透雾方法 - Google Patents

Abstract

本发明公开一种基于双量估计法的单光子激光透雾方法。步骤1：基于Gamma分布建立雾的模型；步骤2：建立目标回波信号模型；步骤3：探测总回波光子；步骤4：对步骤3中探测的总回波光子内的目标信号通过步骤1雾的模型和步骤2目标回波信号模型进行提取。本发明用以解决全参数估计法、单量参数估计法与传统峰值法目标信号提取能力弱的问题。

Description

一种基于双量估计法的单光子激光透雾方法

技术领域

本发明属于Gm-APD激光雷达探测领域；具体涉及一种基于双量估计法的单光子激光透雾方法。

背景技术

随着汽车激光雷达和无人机等新应用的发展，人们对在视觉退化环境下，特别是在浓雾环境下的目标的高分辨率成像越来越感兴趣。当激光传播环境是高散射介质时，主动光学成像的主要限制是激光的吸收和散射，这会在相对较短的传播距离内导致明显的信号衰减。在浓雾中进行深度成像时，目标回波的光子数太小导致无法从浓雾信号中提取出目标信号。因此，开展一种将目标信号从烟雾背景中分离的算法，满足室外浓密烟雾背后目标深度图像探测和短时间成像的需求，提升Gm-APD激光雷达的天候适应性能力，是有必要的。

传统峰值法(Peakselectionalgorithm，PSA)：选取每个像素点中信号的峰值点作为目标的距离位置。

单量参数估计法(Singleparameterestimationalgorithm，SPEA)：根据实测的雾衰减系数和烟雾的Gamma分布模型，采用单参数估计构造算法(SPEA)重建深度图像。

全参数估计法(Allparameterestimationalgorithm，APEA)：将目标信号近似认为是噪声，采用极大似然估计直接对光子分布直方图进行估计，进而重构深度图像。

传统峰值法(Peakselectionalgorithm，PSA)：峰值法只能对低浓度的烟雾背后目标进行信号提取。当烟雾浓度增加时，烟雾回波强度较大，目标回波强度较小，目标峰几乎淹没于烟雾回波信号的下降沿中，基本找不到目标信号。

单量参数估计法(Singleparameterestimationalgorithm，SPEA)：由于雾的稠密性、动态性和非均匀性，在测量其衰减系数时很难准确得到其真实值，因此利用其进行参数估计，误差较大，并且该方案不适合室外实验。

全参数估计法(Allparameterestimationalgorithm，APEA)：由于回波信号中包含烟雾信号和目标信号，忽略目标信号，直接对回波光子直方图进行参数估计的方法，会导致参数估计误差较大。

发明内容

本发明提供一种基于双量估计法的单光子激光透雾方法，可解决上述三个方法中目标信号提取能力弱的问题，提升了后两个算法的烟雾信号估计精度，在极端环境下能更好的提取微弱目标信号。

本发明通过以下技术方案实现：

一种基于双量估计法的单光子激光透雾方法，所述单光子激光透雾方法包括以下步骤：

步骤1：基于Gamma分布建立雾的模型；

步骤2：建立目标回波信号模型；

步骤3：探测总回波光子；

步骤4：对步骤3中探测的总回波光子内的目标信号通过步骤1雾的模型和步骤2目标回波信号模型进行提取。

进一步的，所述步骤1具体为，Gm-APD探测器的触发中，当激光脉冲能量在极其微弱时，背景噪声、暗计数噪声和回波光子产生的初始电子数均服从Poisson分布；所述回波光子包括目标回波、后向散射以及环境光0；根据Poisson分布统计，在t₁至t₂探测间隔内产生q个光电事件的概率为：

其中Nr(t₁,t₂)为在t₁至t₂探测间隔内的平均光子数，表示为：

其中λ(t)为GM-APD探测器中的初始电子的率函数；在探测时间t₁至t₂内不产生光电子的概率为P(0)＝exp(-Nr(t₁,t₂))，所以产生光电子的概率为1-exp(-Nr(t₁,t₂))；由于Gm-APD探测器在距离门内只产生一次探测，在第j个时隙产生探测的概率为：

根据Gm-APD探测器输出的光子数在时间轴上的探测概率分布曲线，反推得到第j个时隙到达探测器焦平面的平均光子数Nr(j)，表示为：

由于烟雾颗粒对激光具有强散射和吸收作用，因此在激光前向传输的过程中光子存在连续散射事件；若每次散射事件的初始光子数为N₀，发生一次散射后的光子数N(t)随时间的变化满足关系；

N(t)＝N₀exp(-αct) (5)

其中，α为激光在雾中的衰减系数，c为光速；

在一次散射事件中的光子概率密度函数为：

由于各散射事件相互独立，每次散射光子的探测时间τ_i均满足指数分布；经过多次散射后的光子探测时间为

其中k为散射次数；Gamma分布是用来解决从开始到k个随机事件都发生所需要的时间问题，k个独立指数随机变量的总和服从Gamma分布，T～GAMMA(k，β)；因此得到激光在烟雾中传输时回波光子数随时间t变化的密度函数：

其中Г(k)是Gamma函数，k和β是形状和逆尺度参数；当k＝1时，Gamma分布就是指数分布：

G_T(t|k＝1,β)＝βexp(-βt) (8)

令β＝αc，利用公式(4)(6)计算得到，经烟雾连续散射后到达探测器焦平面处的光子数为：

其中r与烟雾的后向散射系数相关。

进一步的，所述步骤2具体为，将高斯函数作为目标回波信号的拟合模型，拟合的模型表示为：

其中，中心位置t_target与目标距离有关，t_target＝2d/c，其中d为目标距离；在模型中引入FWHM来表征回波信号的特性；激光发射脉冲与目标表面作用后的FWHM为τ_p，τ_p与σ成正比，通过以下公式得到：

由公式(11)-(13)得到目标信号光子在时间t时的光子分布为：

式中s与目标的后向散射系数相关。

进一步的，所述步骤3具体为，由于烟雾距离雷达***较近，且浓度较大，导致其回波信号强度大于目标信号强度，探测的回波信号呈双峰分布，并且目标信号回波峰存在于烟雾信号的下降沿；在时间t时探测到的回波光子数如下：

进一步的，所述步骤4具体包括以下步骤：

步骤4.1：回波信号的时间剖面估计；

步骤4.2：基于步骤4.1的时间剖面估计，确定雾信号位置；

步骤4.3：基于步骤4.2的雾信号位置，进行烟雾信号估计；

步骤4.4：基于步骤4.3的烟雾信号估计，进行目标信号估计。

进一步的，所述步骤4.1具体为，光子计数雷达的每个像素点记录了接收到的光子的飞行时间，因此得到回波光子的飞行时间直方图，根据成像帧数得到在每个时间bin内的探测概率，进而由公式(4)计算得到到达探测器焦平面处的光子分布Nf(t)，然后利用高斯函数进行滤波，得到时间轮廓估计。

进一步的，所述步骤4.2具体为，对烟雾和目标回波信号的峰值位置进行估计；利用连续小波变化对估计后的时间轮廓信号进行目标和烟雾的峰值检测；

所述连续小波变换公式如下：

其中CoC_s是CWT系数，s(t)是通过高斯拟合获得的信号，

是母小波，p和q是缩放因子和移位因子，/>

是缩放和平移小波。

进一步的，所述步骤4.3具体为，利用计算得到的N_f(t)进行G_T(t|k,β)的参数估计；烟雾信号估计主要有以下步骤：

步骤4.3.1：根据CWT变换得到烟雾回波信号的起始位置τ_a和结束位置τ_b，将该范围τ_a至τ_b内的回波光子数提取出来得到N_f(τ_ab)；

步骤4.3.2：将步骤4.3.1的N_f(τ_ab)进行归一化处理，计算得到其均值E(τ_ab)和方差D(τ_ab)。由GAMMA分布函数可知，如果t～GAMMA(k，β)，均值和方差分别为E(t)＝k/β，D(t)＝k/β²。因此，在时间τ_a至τ_b内的β＝E(τ_ab)/D(τ_ab)；

步骤4.3.3：根据步骤4.3.2求得的β，利用极大似然估计得到分布参数k的估计值，最终得到估计的烟雾回波信号分布G_T(t|k,β)；

步骤4.3.4：由CWT得到的烟雾回波信号N_f(τ_ab)的峰值强度及位置，将步骤4.3.3得到是G_T(t|k,β)根据N_f(τ_ab)的最大值进行尺度统一，变化尺度即为公式(9)中的r，经峰值平移后得到烟雾回波信号G_Tr(t|k,β)。

进一步的，所述步骤4.4具体为，将计算得到的总体回波光子数N_f(t)与估计的烟雾信号G_Tr(t|k,β)两个信号相减，得到目标信号的初始光子数分布S_T(t)；利用激光发射脉冲与目标表面相互作用后的反射回波信号拟合得到FWHMτ_p，进而根据公式(14)得到高斯形状的***响应函数R；采用互相关法从每个像素的直方图中提取目标的深度信息；对于每个像素，互相关系数C_r(t)是通过目标回波的直方图S_T(t)与每个像素的***响应函数R计算得到，表达式如下：

计算结果Cr(t)中最大的点为目标的距离位置；

通过对每个像素点中目标的距离估计，最终得到整个目标的深度图像。

本发明的有益效果是：(优点，多多益善)

附图说明

附图1为无烟雾存在的目标的标准图像。

附图2为无烟雾存在的目标的标准图像的直方图分布图。

附图3为本发明实验场景及目标的位置分布图。

附图4为不同衰减系数下的目标单个像素回波信号分布及信号提取情况图。

附图5为不同衰减长度下的目标深度图像重构结果对比图，其中，(a)为利用峰值法对时间轮廓估计进行目标信号提取图，(b)为根据测量得到的衰减长度值估计烟雾回波信号图，(c)为根据测量得到的光子飞行时间的直方图原始数据估计得到的烟雾回波信号图，(d)为本发明提取的目标信号图。

附图6为不同采集时间下的目标深度图像重构结果对比，其中，F表示采集帧数(a)为利用峰值法对时间轮廓估计时采集帧数为3000、7500、12500和17500时的对比图，(b)，(c)，(d)

附图7为不同采集帧数下的深度图像重构结果评价指标对比，其中，(a)为AL＝1.2时重构图像的TR随采集帧数的变化曲线图，(b)为AL＝1.2时重构图像的RARE随采集帧数的变化曲线图，(c)为AL＝3.6时重构图像的RARE随采集帧数的变化曲线图，(d)为AL＝3.6时重构图像的TR随采集帧数的变化曲线图。

附图8为少量采集帧数下不同算法重构图像的直方图分布情况图，其中，(a)衰减长度AL＝1.2，F＝3000，A、B和Background分别对应标准图像中目标A、B和背景的位置图，(b)衰减长度AL＝2.7，F＝3000，A、B和Background分别对应标准图像中目标A、B和背景的位置图。

附图9为本发明方法流程图。

具体实施方式

下面将结合本发明实施例中的附图对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

一种基于双量估计法的单光子激光透雾方法，所述单光子激光透雾方法包括以下步骤：

步骤1：基于Gamma分布建立雾的模型；

步骤2：建立目标回波信号模型；

步骤3：探测总回波光子；

步骤4：对步骤3中探测的总回波光子内的目标信号通过步骤1雾的模型和步骤2目标回波信号模型进行提取。

进一步的，所述步骤1具体为，基于Gm-APD探测器的脉冲激光雷达通过光子飞行时间法进行目标的距离测量。作为光子计数雷达的重要组成部分，激光器向目标发射周期性的光束，同时启动时间数字转换器(TDC)开始计时，在接收到回波脉冲时触发TDC停止计时，由TDC记录的飞行时间获得目标距离；Gm-APD探测器的触发中，当激光脉冲能量在极其微弱时，背景噪声、暗计数噪声和回波光子产生的初始电子数均服从Poisson分布；所述回波光子包括目标回波、后向散射以及环境光0；根据Poisson分布统计，在t₁至t₂探测间隔内产生q个光电事件的概率为：

其中Nr(t₁,t₂)为在t₁至t₂探测间隔内的平均光子数，表示为：

其中λ(t)为GM-APD探测器中的初始电子的率函数；在探测时间t₁至t₂内不产生光电子的概率为P(0)＝exp(-Nr(t₁,t₂))，所以产生光电子的概率为1-exp(-Nr(t₁,t₂))；由于Gm-APD探测器在距离门内只产生一次探测，在第j个时隙产生探测的概率为：

根据Gm-APD探测器输出的光子数在时间轴上的探测概率分布曲线，反推得到第j个时隙到达探测器焦平面的平均光子数Nr(j)，表示为：

由于烟雾颗粒对激光具有强散射和吸收作用，因此在激光前向传输的过程中光子存在连续散射事件；若每次散射事件的初始光子数为N₀，发生一次散射后的光子数N(t)随时间的变化满足关系；

N(t)＝N₀exp(-αct) (5)

其中，α为激光在雾中的衰减系数，c为光速；

在一次散射事件中的光子概率密度函数为：

由于各散射事件相互独立，每次散射光子的探测时间τ_i均满足指数分布；经过多次散射后的光子探测时间为

其中k为散射次数；Gamma分布是用来解决从开始到k个随机事件都发生所需要的时间问题，k个独立指数随机变量的总和服从Gamma分布，T～GAMMA(k，β)；因此得到激光在烟雾中传输时回波光子数随时间t变化的密度函数：

其中Г(k)是Gamma函数，k和β是形状和逆尺度参数；当k＝1时，Gamma分布就是指数分布：

G_T(t|k＝1,β)＝βexp(-βt) (8)

令β＝αc，利用公式(4)(6)计算得到，经烟雾连续散射后到达探测器焦平面处的光子数为：

其中r与烟雾的后向散射系数相关。

进一步的，所述步骤2具体为，由于接收的回波信号是激光发射脉冲与不同距离处目标表面反射信号强度的卷积结果，并且多数雷达***的发射脉冲和目标表面反射信号近似满足高斯分布；将高斯函数作为目标回波信号的拟合模型，拟合的模型表示为：

其中，中心位置t_target与目标距离有关，t_target＝2d/c，其中d为目标距离；由于回波信号的FWHM参数反映了表面粗糙度，其他一些表面几何特性以及发射脉冲的宽度；在模型中引入FWHM来表征回波信号的特性；激光发射脉冲与目标表面作用后的FWHM为τ_p，τ_p与σ成正比，通过以下公式得到：

由公式(11)-(13)得到目标信号光子在时间t时的光子分布为：

式中s与目标的后向散射系数相关。

进一步的，所述步骤3具体为，考虑激光雷达***应用场景为烟雾存在于目标和***中间，并且烟雾浓度较大，处于一种非线性、非稳态状态。回波光子中主要包括烟雾回波光子，目标回波光子以及噪声光子；相对于烟雾强度来说噪声强度极其微弱，因此，本发明模型不考虑噪声光子；由于烟雾距离雷达***较近，且浓度较大，导致其回波信号强度大于目标信号强度，探测的回波信号呈双峰分布，并且目标信号回波峰存在于烟雾信号的下降沿；在时间t时探测到的回波光子数如下：

进一步的，由于面阵光子计数雷达的每个像素具有时间分辨能力，每个像素中的回波信号包括烟雾后向散射光和信号光，我们的目的是将目标信号从烟雾背景信号中分离出来。所述步骤4具体包括以下步骤：

步骤4.1：回波信号的时间剖面估计；

步骤4.2：基于步骤4.1的时间剖面估计，确定雾信号位置；

步骤4.3：基于步骤4.2的雾信号位置，进行烟雾信号估计；

步骤4.4：基于步骤4.3的烟雾信号估计，进行目标信号估计。

进一步的，所述步骤4.1具体为，光子计数雷达的每个像素点记录了接收到的光子的飞行时间，因此得到回波光子的飞行时间直方图，根据成像帧数得到在每个时间bin内的探测概率，进而由公式(4)计算得到到达探测器焦平面处的光子分布Nf(t)，然后利用高斯函数进行滤波，得到时间轮廓估计。

进一步的，所述步骤4.2具体为，由于目标与烟雾距离较小，导致回波信号存在波峰重叠的现象。为了将目标信号从背景烟雾信号中分离，对烟雾和目标回波信号的峰值位置进行估计；小波变换作为一种信号分析工具，在峰值检测方面取得了很大的成功；连续小波变换(CWT)能实现对激光雷达全波形回波信号分量数量和位置的检测。因此，利用连续小波变化对估计后的时间轮廓信号进行目标和烟雾的峰值检测；

所述连续小波变换公式如下：

其中CoC_s是CWT系数，s(t)是通过高斯拟合获得的信号，

是母小波，p和q是缩放因子和移位因子，/>

是缩放和平移小波。

进一步的，所述步骤4.3具体为，激光经过烟雾传输后的回波光子数随时间t满足Gamma分布，因此，利用计算得到的N_f(t)进行G_T(t|k,β)的参数估计；由于N_f(t)中包含烟雾回波和目标回波信息，忽略目标信号，直接对N_f(t)进行参数估计的方法，会导致参数估计误差较大。为提高G_T(t|k,β)参数估计的准确性，在参数估计时只针对烟雾回波信号进行计算。烟雾信号估计主要有以下步骤：

步骤4.3.1：根据CWT变换得到烟雾回波信号的起始位置τ_a和结束位置τ_b，将该范围τ_a至τ_b内的回波光子数提取出来得到N_f(τ_ab)；

步骤4.3.2：将步骤4.3.1的N_f(τ_ab)进行归一化处理，计算得到其均值E(τ_ab)和方差D(τ_ab)。由GAMMA分布函数可知，如果t～GAMMA(k，β)，均值和方差分别为E(t)＝k/β，D(t)＝k/β²。因此，在时间τ_a至τ_b内的β＝E(τ_ab)/D(τ_ab)；

步骤4.3.3：根据步骤4.3.2求得的β，利用极大似然估计得到分布参数k的估计值，最终得到估计的烟雾回波信号分布G_T(t|k,β)；

步骤4.3.4：由CWT得到的烟雾回波信号N_f(τ_ab)的峰值强度及位置，将步骤4.3.3得到是G_T(t|k,β)根据N_f(τ_ab)的最大值进行尺度统一，变化尺度即为公式(9)中的r，经峰值平移后得到烟雾回波信号G_Tr(t|k,β)。

进一步的，所述步骤4.4具体为，将计算得到的总体回波光子数N_f(t)与估计的烟雾信号G_Tr(t|k,β)两个信号相减，得到目标信号的初始光子数分布S_T(t)；由于烟雾分布的非稳态性和***的不稳定性等因素，导致各像素点的目标回波光子S_T(t)存在非均匀分布。利用激光发射脉冲与目标表面相互作用后的反射回波信号拟合得到FWHMτ_p，进而根据公式(14)得到高斯形状的***响应函数R；采用互相关法从每个像素的直方图中提取目标的深度信息；对于每个像素，互相关系数C_r(t)是通过目标回波的直方图S_T(t)与每个像素的***响应函数R计算得到，表达式如下：

计算结果Cr(t)中最大的点为目标的距离位置；

通过对每个像素点中目标的距离估计，最终得到整个目标的深度图像。

为了检验雷达***的稳定性，在实验开始前首先需要对无烟雾存在的目标进行成像，分析回波信号的直方图是否满足实验需求。成像结果如图1所示，目标A和B的轮廓信息与图3中的实际目标情况很符合，并且与背景完全分离。图1的直方图结果如图2所示，横坐标为光子飞行时间，一个时间Bin代表1ns，纵坐标为该时间Bin的光子数。曲线主要包括三个峰，其中第一个峰与第二个峰相差3个时间Bin，即为3ns。根据激光在空气中的速度约为3×10⁸m/s，可以计算得到两个峰之间的距离ΔL₁＝3×10^-9×3×10⁸/2＝0.45m，第一个峰与第三个峰相差38个时间Bin，相差距离ΔL₂约为5.7m。计算结果ΔL₁、ΔL₂与实际目标A、B之间的距离和目标A与背景板之间的距离相等。由于激光器和探测器等仪器设备的时间抖动，测量得到的目标位置与实际位置相差3个时间Bin，但是目标的相对位置是与实际符合的，因此，***的测量结果满足实验要求。本发明将图1做为后续烟雾图像处理结果评判的标准图像。图2标准图像的直方图分布，三个峰依次表示为目标A、目标B和背景板，图3为实验场景及目标的位置分布。

为了检验我们算法对透过烟雾的微弱目标信号的提取能力，***对目标B进行20000帧成像。选择目标范围内的一个像素点，分别得到衰减系数AL从1.2增加至3.6时的目标回波信号，并且根据第3节中的目标信号提取算法，提取出烟雾和目标信号，结果如图4所示，红色虚线为拟合的烟雾信号，绿色实线为目标信号。随着衰减系数的增加，由于烟雾对激光的后向散射能力的增加，透过烟雾的光子数减少，烟雾回波信号的强度逐渐增大，目标回波信号逐渐减小，叠加的回波信号F_T(t)由双峰变为单峰分布，回波脉宽逐渐减小。当衰减系数AL＝3.6时，目标信号已经完全被烟雾信号掩盖，此时的SBR值为-23.2dB，相较于AL＝1.2时减少了16.9dB，但是我们的算法仍然能够实现对微弱目标信号的提取。

图4不同衰减系数下的目标单个像素回波信号分布及信号提取情况。Histogram:计算出的初始反射光子的直方图。Gaussianfit:高斯拟合得到的信号曲线拟合。Gammafit:背景信号估计。Signal:直方图与Gamma拟合的差，即提取的目标信号光子。Target:目标信号估计。Modelfit:Gamma拟合与目标信号之和。

***对不同衰减长度下的目标A、B进行连续成像20000帧，与之对应的采集时间为1s，利用不同的算法进行目标信号提取，重构的目标深度图像如图5所示。随着衰减长度AL从1.2增加至3.6时，烟雾对目标信号的干扰能力逐渐增强，导致四种算法重构得到的目标深度图像的TR值越来越小、RARE值越来越大，其中PSA、SPEA、APEA和本发明方法重构图像的TR值分别减少了99.3％、86.9％、48.1％和65.4％。虽然本发明方法的TR值减少的百分比相对于APEA算法大了17.3％，但是相对于PSA和SPEA分别减小33.9％和21.5％。当AL＝3.6时，相对于SPEA和APEA，本发明方法得到的TR分别提升了0.2300和0.0408，RARE分别减小了0.3793和0.0231，恢复的目标轮廓更加完整，提取的目标位置更加准确。因此，本发明方法在衰减长度变化时具有较好的稳定性和目标信号提取能力。

图5不同衰减长度下的目标深度图像重构结果对比，不同的列显示的是在同一衰减长度下不同算法的重构结果，不同的行显示的是同一算法在不同衰减长度下的重构结果。TR为目标恢复度，RARE为相对平均测距误差。(a)为利用峰值法对时间轮廓估计进行目标信号提取。(b)为根据测量得到的衰减长度值估计烟雾回波信号。(c)为根据测量得到的光子飞行时间的直方图原始数据估计得到的烟雾回波信号。(d)为本发明方法提取的目标信号。

当衰减长度AL＝1.2时，对目标进行连续多帧成像，采集帧数F分别为3000、7500、12500和17500，与之对应的采集时间分别为0.150s、0.375s、0.625s和0.875s，利用不同的算法进行目标信号提取，重构的目标深度图像如图6所示。虽然本发明方法在F＝3000时重构的深度图像相对于SPEA和APEA具有更多的噪声，但是得到的TR分别提升了0.2361和0.2583，RARE分别减小了0.1299和0.1137。随着采集时间的增大，本发明方法重构的目标轮廓逐渐趋于完整，而SPEA和APEA算法与之相反。因此，本发明方法在较短采集时间下具有最好的目标恢复能力。

为了研究采集时间对目标深度图像重构结果的影响，本发明对衰减长度AL分别为1.2和3.6的图像进行目标信号提取，成像帧数从3000增加至20000，重构图像的TR和RARE随采集帧数的变化曲线如图5所示。当烟雾浓度较低时(AL＝1.2)，随着采集帧数的增加，四种算法重构图像的TR逐渐增大，RARE逐渐减小。相对于SPEA和APES算法，本发明方法重构图像的TR分别平均提升了0.3271和0.3327，RARE分别平均减小了0.1356和0.1146。当烟雾浓度较高时(AL＝3.6)，随着采集帧数的增加，APEA和本发明方法重构图像的TR逐渐增大，而SPEA的TR先增大后减小，并且其值接近于零。虽然本发明方法得到的RARE值与APEA算法的结果几乎相等，但是TR平均提升了0.0613。综合来看，本发明方法具有最大的TR值和最小的RARE值，也就是说本发明方法对高散射介质或采集时间非常短的目标深度图像进行重构时，具有较好的稳定性和目标信号提取能力。

图7不同采集帧数下的深度图像重构结果评价指标对比。(a)、(b)分别为AL＝1.2时重构图像的TR和RARE随采集帧数的变化曲线，(c)、(d)分别为AL＝3.6时重构图像的TR和RARE随采集帧数的变化曲线。

重构目标深度图像的作用就是能够将图像中的目标按距离位置进行区别。因此，本发明将少量采集帧数下重构后的目标深度图像进行直方图统计，结果如图8所示。当烟雾衰减长度AL＝1.2或AL＝2.7时，APEA和SPEA算法得到的目标距离值与标准图像的目标距离值相差较大，但是本发明方法重构的目标深度图像直方图与标准图像一致，包含三个回波峰。利用峰值法对三种算法的直方图进行目标搜索，得到重构图像中目标A、B和背景与真实目标的距离误差，具体如表1所示。当衰减长度AL＝1.2或2.7时，本发明方法的平均误差为0.33Bin，比传统算法提升至少4Bins。结果表明，本发明方法在烟雾浓密或采集时间较小的情况下，能够将多目标按距离值区分，并且在具有最高灵敏度的同时，测距精度最高。

图8少量采集帧数下不同算法重构图像的直方图分布情况，其中A、B和Background分别对应标准图像中目标A、B和背景的位置。(a)衰减长度AL＝1.2，F＝3000，(b)衰减长度AL＝2.7，F＝3000。

表1不同算法重构的目标距离与真实距离的误差对比情况

/>

Claims (5)

1.一种基于双量估计法的单光子激光透雾方法，其特征在于，所述单光子激光透雾方法包括以下步骤：

步骤1：基于Gamma分布建立雾的模型；

步骤2：建立目标回波信号模型；

步骤3：探测总回波光子；

步骤4：对步骤3中探测的总回波光子内的目标信号通过步骤1雾的模型和步骤2目标回波信号模型进行提取；

所述步骤1具体为，Gm-APD探测器的触发中，当激光脉冲能量在极其微弱时，背景噪声、暗计数噪声和回波光子产生的初始电子数均服从Poisson分布；所述回波光子包括目标回波、后向散射以及环境光0；根据Poisson分布统计，在t₁至t₂探测间隔内产生q个光电事件的概率为：

其中Nr(t₁,t₂)为在t₁至t₂探测间隔内的平均光子数，表示为：

其中λ(t)为GM-APD探测器中的初始电子的率函数；在探测时间t₁至t₂内不产生光电子的概率为P(0)＝exp(-Nr(t₁,t₂))，所以产生光电子的概率为1-exp(-Nr(t₁,t₂))；由于Gm-APD探测器在距离门内只产生一次探测，在第j个时隙产生探测的概率为：

根据Gm-APD探测器输出的光子数在时间轴上的探测概率分布曲线，反推得到第j个时隙到达探测器焦平面的平均光子数Nr(j)，表示为：

由于烟雾颗粒对激光具有强散射和吸收作用，因此在激光前向传输的过程中光子存在连续散射事件；若每次散射事件的初始光子数为N₀，发生一次散射后的光子数N(t)随时间的变化满足关系；

N(t)＝N₀ exp(-αct) (5)

其中，α为激光在雾中的衰减系数，c为光速；

在一次散射事件中的光子概率密度函数为：

由于各散射事件相互独立，每次散射光子的探测时间τ_i均满足指数分布；经过多次散射后的光子探测时间为

其中k为散射次数；Gamma分布是用来解决从开始到k个随机事件都发生所需要的时间问题，k个独立指数随机变量的总和服从Gamma分布，T～GAMMA(k，β)；因此得到激光在烟雾中传输时回波光子数随时间t变化的密度函数：

其中Г(k)是Gamma函数，k和β是形状和逆尺度参数；当k＝1时，Gamma分布就是指数分布：

G_T(t|k＝1,β)＝βexp(-βt) (8)

令β＝αc，利用公式(4)(6)计算得到，经烟雾连续散射后到达探测器焦平面处的光子数为：

其中r与烟雾的后向散射系数相关；

所述步骤4具体包括以下步骤：

步骤4.1：回波信号的时间剖面估计；

步骤4.2：基于步骤4.1的时间剖面估计，确定雾信号位置；

步骤4.3：基于步骤4.2的雾信号位置，进行烟雾信号估计；

步骤4.4：基于步骤4.3的烟雾信号估计，进行目标信号估计；

所述步骤4.2具体为，对烟雾和目标回波信号的峰值位置进行估计；利用连续小波变化对估计后的时间轮廓信号进行目标和烟雾的峰值检测；

所述连续小波变换公式如下：

其中CoC_s是CWT系数，s(t)是通过高斯拟合获得的信号，

是母小波，p和q是缩放因子和移位因子，/>

是缩放和平移小波；

所述步骤4.3具体为，利用计算得到的N _f(t)进行G_T(t|k,β)的参数估计；烟雾信号估计主要有以下步骤：

步骤4.3.1：根据CWT变换得到烟雾回波信号的起始位置

和结束位置/>

，将该范围

至/>

内的回波光子数提取出来得到N _f(/>

)；

步骤4.3.2：将步骤4.3.1的N _f(

)进行归一化处理，计算得到其均值E(/>

)和方差D(/>

)；由GAMMA分布函数可知，如果t～GAMMA(k，β)，均值和方差分别为E(t)＝k/β，D(t)＝k/β²；因此，在时间/>

至/>

内的β＝E(/>

)/D(/>

)；

步骤4.3.3：根据步骤4.3.2求得的β，利用极大似然估计得到分布参数k的估计值，最终得到估计的烟雾回波信号分布G_T(t|k,β)；

步骤4.3.4：由CWT得到的烟雾回波信号N _f(

)的峰值强度及位置，将步骤4.3.3得到的G_T(t|k,β)根据N _f(/>

)的最大值进行尺度统一，变化尺度即为公式(9)中的r，经峰值平移后得到烟雾回波信号G_Tr(t|k,β)。

2.根据权利要求1所述一种基于双量估计法的单光子激光透雾方法，其特征在于，所述步骤2具体为，将高斯函数作为目标回波信号的拟合模型，拟合的模型表示为：

其中，中心位置t_target与目标距离有关，t_target＝2d/c，其中d为目标距离；在模型中引入FWHM来表征回波信号的特性；激光发射脉冲与目标表面作用后的FWHM为τ_p，τ_p与σ成正比，通过以下公式得到：

由公式(11)-(13)得到目标信号光子在时间t时的光子分布为：

式中s与目标的后向散射系数相关。

3.根据权利要求2所述一种基于双量估计法的单光子激光透雾方法，其特征在于，所述步骤3具体为，由于烟雾距离雷达***较近，且浓度较大，导致其回波信号强度大于目标信号强度，探测的回波信号呈双峰分布，并且目标信号回波峰存在于烟雾信号的下降沿；在时间t时探测到的回波光子数如下：

4.根据权利要求1所述一种基于双量估计法的单光子激光透雾方法，其特征在于，所述步骤4.1具体为，光子计数雷达的每个像素点记录了接收到的光子的飞行时间，因此得到回波光子的飞行时间直方图，根据成像帧数得到在每个时间bin内的探测概率，进而由公式(4)计算得到到达探测器焦平面处的光子分布N_f(t)，然后利用高斯函数进行滤波，得到时间轮廓估计。

5.根据权利要求1所述一种基于双量估计法的单光子激光透雾方法，其特征在于，所述步骤4.4具体为，将计算得到的总体回波光子数N_f(t)与估计的烟雾信号G_Tr(t|k,β)两个信号相减，得到目标信号的初始光子数分布S_T(t)；利用激光发射脉冲与目标表面相互作用后的反射回波信号拟合得到FWHMτ_p，进而根据公式(14)得到高斯形状的***响应函数R；采用互相关法从每个像素的直方图中提取目标的深度信息；对于每个像素，互相关系数C_r(t)是通过目标回波的直方图S_T(t)与每个像素的***响应函数R计算得到，表达式如下：

计算结果Cr(t)中最大的点为目标的距离位置；

通过对每个像素点中目标的距离估计，最终得到整个目标的深度图像。

Images