[0001] 本发明涉及激光雷达成像技术领域，具体涉及一种基于分数阶全变分正则化与非局部正则化的Gm‑APD激光雷达距离像压缩感知重构方法。

[0002] 盖革模式雪崩光电二极管(Geiger‑Mode ofAvalanche Photodiodes,Gm‑APD)激光雷达基于Gm‑APD的雪崩特性，可探测单光子级的回波信号，具备作用距离远和距离分辨率高等优势，在目标检测、遥感、军事制导、安防监控等领域具有重大应用价值。但随着Gm‑APD探测灵敏度的提高，源于日光、大气后向散射以及探测器暗计数等噪声的影响增强，目标信号极易淹没在噪声中，极大提升了回波信号距离信息提取难度；此外，由于Gm‑APD为概率统计型器件，且Gm‑APD激光雷达系统存在死时间和暗计数，导致输出距离像包含大批失落信息和距离反常噪声，使距离像可靠性大大降低。为了实现高质量距离像的重构，现有研究通常使用局部滤波或全局滤波方法进行距离像复原，这类方法不仅需要对维度极高的图像数据进行采样，并且算法复杂度高，算法运行时间长，不利于工程应用。因此亟需开展一种适用于工程应用的低信背比条件下Gm‑APD激光雷达距离像重构方法研究。

[0003] 目前针对Gm‑APD激光雷达距离像重构方法的研究较多，2015年，在欧洲信号处理会议《利用稀疏单光子数据构建的强度和深度图像的恢复》一文中，赫瑞瓦特大学团队发表了首达光子距离像重构算法，在少数光子计数的限制下，通过不同像素所消耗的时间作为空间反射率特性的表征，再采用全变分(Total Variation,TV)正则化联合空间邻域重建目标信息，但重建时对采集数据量要求很大。同年，在IEEE计算成像汇刊《单光子探测器的光子高效计算三维反射成像》一文中，麻省理工学院将赫瑞瓦特大学的算法扩展为稀疏正则化重建框架，并于2016年将该算法应用于基于盖革雪崩焦平面阵列探测系统中。2016年，在IEEE信号处理通讯《基于子空间联合模型的单光子深度成像》一文中，该团队进一步提出了基于联合子空间模型的重建算法，该算法将距离信息求解转化为约束优化问题，提出了一种改进的压缩采样匹配追踪算法以解决该欠定线性逆问题，但是由于实际应用环境下目标回波光子数量很少，使得算法对噪声和干扰的敏感度较高，对于有强噪声干扰的情况下，该方法重构精度会大幅下降。2017年，在IEEE计算成像汇刊《众多光子中的少数光子：光子高效有源成像的信号和噪声解混》一文中，波士顿大学团队提出了UA(Unmixing Algorithm)算法，该算法首先基于噪声分布的随机性与信号在回波时刻分布上的一致性，构建超像素数据直方图消除噪声，再利用TV正则化对图像进行平滑，但是该算法计算复杂度高，所需的成像时间相对较长。2019年，在IEEE图像处理汇刊《以非局部空间相关性恢复稀疏的3D单光子数据》一文中，陈松懋等提出了一种OPN3DR(Optimized  based  nonlocal 3D restoration)算法，用于在光子匮乏状态(每个像素少于一个光子)情况下获取的单光子三维激光雷达图像的空间相关性和非局部恢复。该算法通过Laplacian矩阵引入像素间的非局部空间关系以提升重建质量，同时，还提出了一种压缩计算量的非均匀采样方法，不仅能应对ppp小于1等极端情况，并且提升了算法的运行速度。但是该算法在强噪声环境下重建效果较差，并且算法复杂度高，难以进一步工程化实现。2020年，在光子学研究《通过混合惩罚进行首光子成像》一文中，北京航空航天大学团队提出了一种基于小波域的1‑范数正则化与TV正则化的混合正则化的图像重建算法。1‑范数正则化可以在一定程度上保持图像的稀疏性，TV正则化可以保持图像的纹理特征。这两种方法的结合可以相互补充，平衡了去噪效果与图像细节的保持，重建效果优于单正则项方法。但是在信噪比较低的条件下，重建图像的质量会受到较大影响。同年，在光子学《基于首信号光子单元法的高效抗噪光子计数成像》一文中，中国科学院空间光电精密测量科技重点实验室设计了首序列成像方法，该方法对噪声较高的像素累积时间较长，直至可以在回波中分辨出信号簇再扫描下一个像素，并通过中值滤波与TV正则化对逐点重建的结果进行处理，以进一步提升成像质量。虽然该算法可以实现强噪声环境下的精确重构，但对于大规模复杂的成像任务，算法运行可能会进一步加大。2022年，在艾斯维尔《基于光子数的Gm‑APD激光雷达少帧检测距离像重构系统算法》一文中，哈尔滨工业大学团队提出了一种二维双阈值近似去噪算法，通过设定子像素邻域阈值实现噪声的去除，并采用邻域平滑算法恢复噪声像素的真实值，该方法仅构建子像素与周围像素的联系进行去噪和复原，复原的目标深度信息偏差较大，且相比于全局滤波算法，在目标整体平滑效果上较差。

[0004] 为实现低信背比条件下Gm‑APD激光雷达距离像的精确重构，本发明提出了一种基于分数阶全变分正则化与非局部正则化的Gm‑APD激光雷达距离像压缩感知重构方法。

[0036] 结合图1至图14说明本实施方式，基于分数阶全变分正则化与非局部正则化的Gm‑APD激光雷达距离像压缩感知重构方法，该方法首先基于峭度因子峰值法在低信背比条件下提取目标距离像；其次构建分数阶全变分正则化和非局部正则化的Gm‑APD距离像压缩感知重构模型；然后利用增广拉格朗日乘子法求解该模型，实现目标距离像精准重构，最后通过仿真和实验对本文所提算法的有效性进行验证。其具体步骤如下：

[0037] 步骤一、引入峭度因子，实现低信背比目标的距离像提取，具体实现过程为：

[0038] 步骤一一、计算滑窗峭度因子Qd；

[0039] 确定信号的总体均值μ和总体标准差σ，计算滑窗峭度因子Qd，用下式表示为：

[0040]

[0041] 式中，xi为接收到的含噪回波信号，n为正整数；

[0042] 步骤一二、如图2所示，生成含噪回波信号的触发概率直方统计结果；

[0043] 对接收到的含噪回波信号的时序特征以及触发概率，构建出含噪回波信号的触发概率直方统计结果Y：

[0044] Y＝Histogram(xi)

[0045] 式中，Histogram(·)表示对输入量进行直方统计。

[0046] 步骤一三、如图3所示，令含噪回波信号触发概率直方统计结果Y与滑窗峭度因子Qd进行卷积操作，获得卷积操作后的直方统计结果Y'：

[0047] Y'＝Y*Qd

[0048] 根据所述直方统计结果Y'即可以得到凸特征增强、平滑特征被抑制的曲线。

[0049] 步骤一四、完成对目标距离像的提取，通过对各个阵元处回波信号进行寻址，能够提取对应的距离信息，最终完成对目标距离像的提取：用下式表示为：

[0050] Upos＝argmax[Y′(ipos)]

[0051] 式中，Upos是图像pos位置处的距离信息，ipos表示图像pos位置处卷积操作后的直方统计结果Y'的第i个位置，arg max[·]表示寻找输入量最大值，并返回其对应的位置信息。

[0052] 步骤二、对信号进行压缩采样，构建基于分数阶全变分正则化和非局部正则化的距离像压缩感知重构模型，具体步骤如下：

[0053] 步骤二一、目标距离像分割；

[0054] 将目标距离像分割等分成ns×ns个相等且不重叠的图像块，ns为给定正整数；由图像所有位置处距离信息构成的N维距离矩阵u，设定矩阵u中给定图像块的中心像素为ui，ui周围各图像块的中心像素为uj；

[0055] 步骤二二、计算权重；计算周围各图像块的中心像素uj相对给定图像块的中心像素ui的权重wij；

[0056]

[0057] 式中，h为高斯核的控制因子，ci为归一化常数。

[0058] 步骤二三、基于非局部均值滤波思想，构造非局部正则化惩罚项NR(u)；

[0059]

[0060] 式中， 为给定图像块周围各图像块的中心像素uj各自的权重系数组成的列向量的转置形式，ki为给定图像块周围各图像块的中心像素距离值组成的列向量，非局部正则化惩罚项的矩阵形式为：

[0061]

[0062] 式中，W为权重系数矩阵，wij为各中心像素的权重系数， 表示L2正则化项。

[0063] 步骤二四、将步骤二三构造的非局部正则化惩罚项引入基于全变分正则化和非局部正则化的距离像压缩感知重构模型：

[0064] s.t.Au＝b

[0065] 式中，D是全变分算子。

[0066] 步骤二五、基于分数阶微积分形式，将步骤二四构建的基于全变分正则化和非局部正则化的距离像压缩感知重构模型拓展至分数域。

[0067] s.t.Au＝b

[0068] 式中，u是由图像所有位置处距离信息构成的N维距离矩阵，Dv是分数阶全变分算子，v是分数阶阶次，表示L1正则化项，指输入量中所有元素绝对值之和，NR(u)为关于u非局部正则化惩罚项，α为权重系数，A为观测矩阵，b为M维观测向量。

[0069] 步骤三、利用增广拉格朗日乘子法对距离像重构模型求解，具体步骤为：

[0070] 步骤三一、将所述距离像重构模型等价转化为最优化问题；引入辅助变量和，将距离像重构模型进行等价转化：

[0071] s.t.Dvu＝w,u＝x,Au＝b

[0072] 步骤三二、将该问题表达为增广拉格朗日函数LA(w,u,x)；

[0073]

[0074] 式中，β，θ，μ是与惩罚项相关的正则化参数，T表示对向量进行转置。

[0075] 步骤三三、利用迭代求解该模型。

[0076]

[0077]

[0078] 式中，k表示迭代次数，ω、γ、λ分别是约束条件Dvu＝w，u＝x，Au＝b的拉格朗日乘子。

[0079] 给定未知量、u、x中的任意两个，可以求解余下的量，结果如下：

[0080] 给定u、x，与w相关的子问题为：

[0081]

[0082] 给定w、x，与u相关的子问题为：

[0083]

[0084] 给定w、u，与x相关的子问题为：

[0085]

[0086] 本实施方式中，采用蒙特卡罗方法对回波数据进行模拟，同时利用客观评价指标，即：利用目标还原度及峰值信噪比对本发明所提算法进行评价，并与其他算法进行比较。模型的模拟深度图像如图4所示，图中横纵坐标表示像素点位置，颜色值代表目标距离相机的相对距离。

[0087] 为验证低信背比情况下基于峭度因子法距离参数提取精度，将分别采用峰值法(Peak)和本发明采用的峭度因子法(KFM)对SBR为0.15，0.18，0.2的回波数据进行处理，当统计帧数为40，50，60，70，80时，两种方法的成像结果如图5所示。

[0088] 从图5可知，当SBR为0.20时，峰值法仅能准确提取出瓶子的个别距离值，但本发明采用的峭度因子法，能获取瓶子的距离像，且随着统计帧数的增加，杯子的距离像越清晰。当SBR为0.15和0.18时，峰值法不能获取距离像，而本发明采用的峭度因子法仍能较清晰的获得杯子的距离像。不过，利用峭度因子法获取的距离像仍然需要进一步去噪以提高成像质量。

[0089] 在不同SBR与不同统计帧数下分别进行1000组蒙特卡罗重复实验，采用目标还原度K、峰值信噪比PSNR两个指标对距离像质量进行评价，如图6和图7所示。

[0090] 从图6、图7可知，本发明方法目标还原度优于峰值法，且当SBR等于0.2，成像帧数等于40，本发明方法的目标还原度达到0.8245，峰值信噪比(PSNR)达到24.8462。

[0091] 为探究采样率对本发明方法重构性能的影响，设置统计帧数为40帧，采用本发明重构方法分别对在SBR为0.15、0.18和0.2时利用本发明采用的峭度因子法获取的距离像进行重构，分别进行1000组蒙特卡罗重复实验，取目标还原度均值对距离像重构质量进行对比评价，结果如图8所示。

[0092] 由图8可知，随着采样率增加，本发明方法的目标还原度K递增，当采样率为0.6时，不同信背比下的目标还原度均高于0.75，此时已经能够基本实现距离像的重构，故本发明选择采样率为0.6。

[0093] 为验证不同信背比下本发明方法对低信背比目标的距离像重构性能，将经过峭度因子法(KFM)提取之后不加去噪算法的距离像作为参照组，分别采用全变分(TV)去噪算法，分数阶全变分(FOTV)去噪算法和本发明的方法(Proposed)，对在统计帧数为50，SBR为0.15、0.18、0.2时利用本发明采用的峭度因子法获取的距离像进行重构，结果如图9所示。

[0094] 从图9可知，经本发明方法重构后的目标区域更为平滑，并且去除目标位置大量距离反常值噪声，采用量化指标对图像质量进行评价，在不同SBR下分别进行1000次蒙特卡罗重复实验，得到各个指标均值，结果如图10和图11所示。表1为各算法的量化评价指标。

[0095] 表1

[0096]

[0097] 结合表1可知，本发明方法与对比算法均可以提升目标还原度K、峰值信噪比PSNR，且相较于对比算法，在帧数为50时，SBR为0.2时，本发明方法的K至少提升2.66％，PSNR至少提升6.69％。

[0098] 本发明搭建了64×64阵列Gm‑APD激光雷达对所提算法性能进行验证，系统原理如图12所示。该系统使用Gm‑APD作为探测器，首先通过发射光学系统将激光照射到目标区域，同时脉冲激光器产生触发信号来启动Gm‑APD的计时功能。当激光在目标表面产生漫反射后，接收光路将光信号聚焦到Gm‑APD的焦平面，此时Gm‑APD停止计时，之后读出电路将激光光子的飞行时间进行时间‑数字转换，并将数据传递给上位机进行信号处理，然后经过算法处理后显示出深度图像。

[0099] 成像实验验证:为了验证本发明方法的图像重构性能，在白天对463.2‑446.1m的居民楼进行成像实验。目标场景图如图13所示，图中方框内为探测区域。

[0100] 当统计帧数为200帧时，分别采用TV和FOTV算法对目标深度图像进行重构，结果如图14中的(b)和(c)所示。

[0101] 从图14中的(d)可知，本发明的去噪算法比对比方法更能恢复目标距离像中的噪声，产生更平滑的目标区域，表2为成像实验量化评价指标。

[0102] 表2

[0103]

[0104] 结合表2可知，本发明的去噪方法，目标还原度至少提升了4.29％，峰值信噪比为4.6969，两项指标均领先于对比算法，可见本发明方法提升了Gm‑APD距离像重构精度。

[0105] 上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

[0106] 以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。