具体实施方式

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

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

步骤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.45 m，第一个峰与第三个峰相差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:计算 出的初始反射光子的直方图。Gaussian fit:高斯拟合得到的信号曲线拟合。Gammafit:背 景信号估计。Signal:直方图与Gamma拟合的差，即提取的目标信号光子。Target:目标 信号估计。Model fit: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不同算法重构的目标距离与真实距离的误差对比情况