一种激光雷达全波形信号分析方法

摘要

本发明提出一种激光雷达全波形信号分析方法。本发明方法在每次更新操作时在回波信号的数量、入射到激光雷达探测器的噪声信号、每个回波信号的峰值位置和每个回波信号的峰值幅中随机选择任意一个参量继续进行更新操作；参量采用Metropolis算法思想进行更新，参量采用Metropolis算法或者模拟回火算法思想进行更新，参量和分别采用模拟回火算法思想进行更新。本发明能够快速处理特性未知目标的全波形信号，且不受全波形信号信噪比的影响。

说明书

技术领域

本发明属于激光雷达技术领域，具体涉及一种激光雷达全波形信号分析方法。 

背景技术

传统的激光雷达通常只能提取激光光斑内的一个距离特征，无法处理光斑内多个距离特征目标的回波信号。采用具有完整地记录接收信号功能的全波形激光雷达可以得到这种包含多个回波信号的全波形信号，为全波形信号分析提供了基础，所谓全波形信号由激光光斑印记内所有特征距离的回波信号共同产生。为了有效解析单个激光光斑印记内包含的距离特征，需要采用全波形信号分析技术提取出其中包含的回波信号的峰值位置、数量、峰值幅度、噪声共4个关键的全波形参量，分别标记为t₀、k、β、B。其中，回波数量和峰值位置是全波形参量是最为重要的两个参数，用于标记距离特征的数量和相互位置关系并完成全波形信号分解，峰值位置即准确距离上的光子飞行时间。 

全波形信号y可以看成是接收信号F的后验抽样结果，以能量在时间轴上的能量分布形式来表示。全波形信号分析的目的就是利用全波形信号y求解全波形参量t₀、k、β、B，其过程可以看作是利用后验信息对先验信息进行估计。其中，后验信息是全波形信号y，而相对应的先验信息是全波形参量t₀、k、β、B。通常可以采用后验分布概率密度函数来评价全波形参量估计值的准确程度。 

典型的全波形分析策略主要包括采用LM算法的非线性最小二乘波形拟合方法和采用EM算法及其改进算法的全波形参量最大似然值估计方法。前者需要事先预知目标的有效数据模型，只适用于单一特性的场景，例如：植被、建筑等；后者的收敛性能强烈依赖于初始值的预设，无法实现完全的非人工干预全波形参量提取，有时存在失效的风险。 

近来，一种基于贝叶斯估计框架的可逆跳跃马尔可夫蒙特卡罗算法（RJMCMC）（详见IEEE Transactions on Pattern Analysis and Machine Intelligence2007,29(12),2170-2180）被提出作为一种全波形分析策略。在此策略中，其关键之处是采用RJMCMC算法按照固定的顺序以循环形式逐个进行待解全波形参量t₀、k、β和B，并由后验分布概率密度函数评估所有由参量近似解构成的全局近似解（t₀、k、β、B），挑选最为合适的全 局近似解（t₀、k、β、B）作为全局最优解。该方法允许在更新过程中待解全波形参量可以在不同尺度的预设分布模型之间进行跳跃完成对分布空间的探测而不依赖先验信息，从而实现对目标特性未知的全波形信号分析。但该方法对分布空间的快速探测主要依赖于有效的预设分布，而且各参量的转移概率通常较低，从而导致较慢的运算速度。 

发明内容

本发明提出一种激光雷达全波形信号分析方法，能够快速处理特性未知目标的全波形信号，且不受全波形信号信噪比的影响。 

为了解决上述技术问题，本发明提供一种激光雷达全波形信号分析方法，包括以下步骤： 

步骤一、输入一个激光雷达全波形信号y，明确待求参量为：t₀、k、β、B，其中： 

k为回波信号的数量，其取值范围为：{0,1,2Kk_max}，k_max是预设的k的最大值； 

B为入射到激光雷达探测器的噪声信号，其取值范围为：(0,max(y))，max(y)是全波形信号y的最大值； 

t₀为每个回波信号中的光子准确飞行时间构成的k维向量，其元素对应k个回波信号的峰值位置，t₀的所有元素的取值范围为：(0,i_max)，i_max是全波形信号y的长度； 

β为每个回波信号的峰值幅度构成的k维向量，其元素对应k个回波信号的峰值幅度，β的所有元素的取值范围为：(0,max(y))； 

从各待求参量的取值范围内选取任意值分别赋值给参量t₀、k、β和B，完成各待求参量的初始化；在完成初始化的待求参量中选取任意一个参量，以被选中的参量开始进行一次更新操作； 

步骤二、将更新操作后的所有待求参量组成一个全局近似解（t₀、k、β、B），在每次更新操作后组成的全局近似解（t₀、k、β、B）中随机选择任意一个参量继续进行更新操作；重复本步骤，直至更新操作次数的总和达到预先设定的更新操作次数的最大值后停止更新操作； 

步骤三、将每一次更新操作后获得的所有全局近似解分别代入后验分布概率密度函数，求出每一个全局近似解的后验分布概率值，将值最大的后验分布概率值对应的全局近似解作为全局最优解；所述后验分布概率密度函数如公式（1）所示： 

$\pi (t_0,k,\beta ,B|y)=\left(\underset{i=1}{\overset{i_{\max }}{\Pi }}\exp (-(\underset{j=1}{\overset{k}{\Sigma }}f(i,{\beta }_j,t_{0j})+B))\frac{{(\underset{j=1}{\overset{k}{\Sigma }}f(i,{\beta }_j,t_{0j})+B)}^{y_i}}{y_i!}\right)$

$\times \frac{1}{k_{\max }}{\left(\frac{1}{i_{\max }}\right)}^k{f}_G\left(B\right|c,d)\underset{n=1}{\overset{k}{\Pi }}{f}_G\left({\beta }_n\right|a,b)---\left(1\right)$

式（1）中，π(t₀,k,β,B|y)为后验分布概率值；f_G为伽玛分布概率密度函数，a和c为f_G的尺度参量，b和d为f_G的形状参量，且a、b、c、d均为预设的正数；y_i是y在时间i上的离散表示；f(i,β_j,t_0j)为第j个回波信号在时间i上的能量分布，j是不大于k_max的非负整数，i通常取0至i_max间的整数，t_0j和β_j分别为第j个回波信号的光子准确飞行时间和峰值幅度、也是t₀和β的第j个元素，f(i,β_j,t_0j)的函数表达式如公式（2）所示： 

$f(i,{\beta }_j,t_{0j})={\beta }_j\left(\begin{array}{cc}{e}^{-{(t_{1j}-t_{0j})}^2/{2\sigma }^2}{e}^{(i-t_{1j})/{\tau }_1},& i<t_{1j}\\ e^{-{(i-t_{0j})}^2/{2\sigma }^2},& t_{1j}\le i<t_{2j}\\ e^{-{(t_{2j}-t_{0j})}^2/{2\sigma }^2}{e}^{-(i-t_{2j})/{\tau }_2},& t_{2j}\le i<t_{3j}\\ e^{-{(t_{2j}-t_{0j})}^2/{2\sigma }^2}{e}^{-(t_{3j}-t_{2j})/{\tau }_2}{e}^{-(i-t_{3j})/{\tau }_3},& i\ge t_{3j}\end{array}\right)---\left(2\right)$

式（2）中，t_0j、t_1j、t_2j、t_3j间存在稳定的差值大小关系，即t_0j-t_1j＝c₁、t_0j-t_2j＝c₂、t_0j-t_3j＝c₃，其中c₁、c₂、c₃为预设的常数；τ₁、τ₂、τ₃和σ均为预设的常数，且0＜τ₁＜τ₂＜τ₃、σ＝0。 

本发明与现有技术相比，其显著优点在于：（1）可以随机选择参量进行参量更新，一定程度上减小了运算量，提升了求解速度；（2）峰值位置、峰值幅度、噪声三个参量更新中增加了主动回火步骤，实现了近似解空间的临时扩大，实现了参量的快速更新，提高了近似解的探测效率，达到了待解全波形参量近似解快速更新的目的并保持了鲁棒性；（3）不同性质的参量采取不同的更新算法，有效减小了运算量，提高了运行速度。 

附图说明

图1是本发明一种激光雷达全波形信号分析方法流程图。 

图2是本发明实验一中探测目标为毛玻璃与纸箱，毛玻璃与纸箱不同间距平行放置， 获取其七个不同间距的全波形信号的信号分解图，其中，图2(a)是间距为5cm的全波形信号的信号分解图，图2(b)是间距为8cm的全波形信号的信号分解图，图2(c)是间距为10cm的全波形信号的信号分解图，图2(d)是间距为20cm的全波形信号的信号分解图，图2(e)是间距为30cm的全波形信号的信号分解图，图2(f)是间距为50cm的全波形信号的信号分解图，图2(g)是间距为100cm的全波形信号的信号分解图。 

图3是本发明在实验二中获得的低信噪比全波形信号的信号分解图，探测目标为毛玻璃与纸箱，毛玻璃与纸箱不同间距平行放置，其中，图3(a)是间距为10cm的低信噪比全波形信号的信号分解图、图3(b)是间距为20cm的低信噪比全波形信号的信号分解图。 

图4是本发明在实验三中重建的目标在距离上的分布，探测目标为等距平行放置的毛玻璃、N形字母和纸箱，达到了厘米级的水平。 

具体实施方式

如图1所示，本发明激光雷达全波形信号分析方法，包括以下步骤： 

步骤一、输入一个激光雷达全波形信号y，明确待求参量为：t₀、k、β、B，其中： 

k为回波信号的数量，其取值范围为：{0,1,2Kk_max}，k_max是预设的k的最大值，根据预先了解的目标全波形信号中的回波数量可以设定合适的k_max值，以减小计算量； 

B为入射到激光雷达探测器的噪声信号，其取值范围为：(0,max(y))，max(y)是全波形信号y的最大值； 

t₀为每个回波信号中的光子准确飞行时间构成的k维向量，其元素对应k个回波信号的峰值位置，t₀的所有元素的取值范围为：(0,i_max)，i_max是全波形信号y的长度； 

β为每个回波信号的峰值幅度构成的k维向量，其元素对应k个回波信号的峰值幅度，β的所有元素的取值范围为：(0,max(y))； 

从各待求参量的取值范围内选取任意值分别赋值给参量t₀、k、β和B，完成各待求参量的初始化；在完成初始化的待求参量中选取任意一个参量，以被选中的参量开始 进行一次更新操作，即根据选中的不同参量对应的不同更新操作方法进行本次更新操作，其它未被选中的参量参与本次更新操作，所有参量在一次更新操作后有可能被更新也有可能保持不变； 

步骤二、将更新操作后的所有待求参量组成一个全局近似解（t₀、k、β、B），在每次更新操作后组成的全局近似解（t₀、k、β、B）中随机选择任意一个参量继续进行更新操作；重复本步骤，直至更新操作次数的总和达到预先设定的更新操作次数的最大值后停止更新操作； 

步骤三、将每一次更新操作后获得的所有全局近似解分别代入后验分布概率密度函数，求出每一个全局近似解的后验分布概率值，将值最大的后验分布概率值对应的全局近似解作为全局最优解；所述后验分布概率密度函数如公式（1）所示： 

$\pi (t_0,k,\beta ,B|y)=\left(\underset{i=1}{\overset{i_{\max }}{\Pi }}\exp (-(\underset{j=1}{\overset{k}{\Sigma }}f(i,{\beta }_j,t_{0j})+B))\frac{{(\underset{j=1}{\overset{k}{\Sigma }}f(i,{\beta }_j,t_{0j})+B)}^{y_i}}{y_i!}\right)$

$\times \frac{1}{k_{\max }}{\left(\frac{1}{i_{\max }}\right)}^k{f}_G\left(B\right|c,d)\underset{n=1}{\overset{k}{\Pi }}{f}_G\left({\beta }_n\right|a,b)---\left(1\right)$

式（1）中，π(t₀,k,β,B|y)为后验分布概率值；f_G为伽玛分布概率密度函数，a和c为f_G的尺度参量，b和d为f_G的形状参量，且a、b、c、d均为预设的正数；y_i是y在时间i上的离散表示；f(i,β_j,t_0j)为第j个回波信号在时间i上的能量分布，j是不大于k_max的非负整数，i通常取0至i_max间的整数，t_0j和β_j分别为第j个回波信号的光子准确飞行时间和峰值幅度、也是t₀和β的第j个元素，f(i,β_j,t_0j)的函数表示可以参见文献（S.Pellegrini,G.Buller,J.Smith,A.Wallace,S.Cova.Laser-Based Distance Measurement Using Picosecond Resolution TCSPC,Measurement Science and Technology,vol.11,pp.712-716,2000），f(i,β_j,t_0j)的具体表示如公式（2）所示： 

$f(i,{\beta }_j,t_{0j})={\beta }_j\left(\begin{array}{cc}{e}^{-{(t_{1j}-t_{0j})}^2/{2\sigma }^2}{e}^{(i-t_{1j})/{\tau }_1},& i<t_{1j}\\ e^{-{(i-t_{0j})}^2/{2\sigma }^2},& t_{1j}\le i<t_{2j}\\ e^{-{(t_{2j}-t_{0j})}^2/{2\sigma }^2}{e}^{-(i-t_{2j})/{\tau }_2},& t_{2j}\le i<t_{3j}\\ e^{-{(t_{2j}-t_{0j})}^2/{2\sigma }^2}{e}^{-(t_{3j}-t_{2j})/{\tau }_2}{e}^{-(i-t_{3j})/{\tau }_3},& i\ge t_{3j}\end{array}\right)---\left(2\right)$

式（2）中，t_0j、t_1j、t_2j、t_3j间存在稳定的差值大小关系，即t_0j-t_1j＝c₁、t_0j-t_2j＝c₂、t_0j-t_3j＝c₃，其中c₁、c₂、c₃在运算过程中均是常数，需要在处理全波形信号y前设置；τ₁、τ₂、τ₃和σ均为预设的常数，且0＜τ₁＜τ₂＜τ₃、σ＝0，需要在处理全波形信号y前设置。在实际应用中，t_0j、t_1j、t_2j、t_3j、τ₁、τ₂、τ₃和σ会在一定程度上影响本发明方法的运算速度和运算准确性，因此需要根据由单个回波构成的全波形信号的波形调整t_0j、t_1j、t_2j、t_3j、τ₁、τ₂、τ₃和σ以适应各种形状的回波波形，这些参量均可以在处理全波形信号y前设置。 

如前所述，在每次更新操作前需要在参量t₀、k、β和B之间任意随机选择一个参量开始本次更新操作，不同参量对应的不同更新操作方法。 

随机选中参量k后的更新方法： 

参量k在四个待求参量中是一个特例，k经过一次更新操作后，有可能达到状态k′，也有可能达到状态k′′，其中k′＝k+1，k′′＝k-1（k存在于{0，1，2Kk_max}，同时k是参量t₀和β的维数，若k发生更新，参量t₀和β必定发生更新，所以k采用较慢渐变的方式以减小计算量。）。当一次更新操作选定的参量为k时，可以进一步在k′和k′′这两种状态中随机选择k经过本次更新操作后可能达到的状态，然后根据随机选择的结果（选中k′或者k′′的概率相等），按照k′和k′′分别对应的不同的更新方法进行本次更新操作。另外也可以将k′和k′′与参量t₀、β和B放在一起选择本次更新操作的方法，即在k′、k′′、t₀、β和B中选择任意一个参量，将该参量对应的更新方法作为本次更新操作的方法，其中，选中k′或者k′′的概率相等。 

无论选择状态k′还是状态k′′，k均采用Metropolis算法思想进行更新，但具体步骤不同。 

如果随机选择的k经过本次更新操作后可能达到的状态为k′时，本次更新操的计算 过程为： 

（1）、若k＝k_max，则各参量t₀、k、β和B保持不变，跳出本次更新操作，重新选择参量进行下一次更新操作；否则继续步骤（2）； 

（2）、设t₀′表示t₀经过本次更新操作后可能到达的状态，且t₀′为k′维向量，t₀′的前k个元素依次为t₀的前k个元素； 

设β′表示β经过本次更新操作后可能到达的状态，且β′为k′维向量，β′的前k个元素依次为t₀的前k个元素； 

在区间(0,i_max)中选择任意一个服从均匀分布的随机数赋值给t₀′的第k′个元素，在区间(0,max(y))中选择任意一个服从均匀分布的随机数赋值给β′的第k′个元素； 

（3）、将y、t₀和β代入公式（2）求解依据y、t₀和β预估的回波信号的能量分布{f(i,β_j,t_0j)}，j是不大于k的任意非负整数，β_j和t_0j分别是向量β和t₀的第j个元素；将y、t₀′和β′代入公式（2）求解依据y、t₀′和β′预估的回波信号的能量分布{f(i,β_m′,t_0m′)}，m是不大于k′的任意非负整数，β_m′和t_0m′分别是向量β′和t₀′的第m个元素； 

然后将{f(i,β_j,t_0j)}、k和B代入公式（3）求解入射到激光雷达探测器的信号F(i,β,t₀,k,B)，将{f(i,β_m′,t_0m′)}、k′和B代入公式（3）求解预估的入射到激光雷达探测器的信号F(i,β′,t₀′,k′,B)； 

再将F(i,β,t₀,k,B)和y代入公式（4）求解全波形信号y在信号F(i,β,t₀,k,B)的联合分布概率L(y|t₀,k,β,B)，将F(i,β′,t₀′,k′,B)和y代入公式（4）求解全波形信号y在信号F(i,β′,t₀′,k′,B)的联合分布概率L(y|t₀′,k′,β′,B)； 

最后将L(y|t₀,k,β,B)、L(y|t₀′,k′,β′,B)代入公式（5）求解k′的转移概率α(k′,k)，公式（3）、（4）、（5）具体如下所示： 

$F(i,\beta ,t_0,k,B)=\underset{j=1}{\overset{k}{\Sigma }}f(i,{\beta }_j,t_{0j})+B---\left(3\right)$

$L\left(y\right|t_0,k,\beta ,B)=\underset{i=1}{\overset{i_{\max }}{\Pi }}\exp (-F\left(i\right))\frac{F{\left(i\right)}^{y_i}}{y_i!}---\left(4\right)$

式（4）中，y₁,y₂K为y在时间上的离散表示， 

$\alpha (k^{\prime },k)=L\left(y\right|t_0^{\prime },k^{\prime },{\beta }^{\prime },B)-L\left(y\right|t_0,k,\beta ,B)---\left(5\right)$

（4）、若α(k′,k)＞0，则将t₀′、β′和k′分别赋值给t₀、β和k，参量B保持不变，完成本次更新操作，在由t₀、k、β和B组成的新的全局近似解（t₀、k、β、B）中随机选择任意参量进行下一次更新操作； 

（5）若α(k′,k)≤0，则进一步判断此时是否到达最大尝试更新次数N_k（预先设置），其中步骤（2）至步骤（4）的过程是一次尝试更新；若到达最大尝试更新次数N_k，则各参量t₀、k、β和B保持不变，跳出本次更新操作，重新选择参量进行下一次更新操作；若未到达最大尝试更新次数N_k，则跳到步骤（2）继续本次更新操作。 

如果随机选择的k经过本次更新操作后可能达到的状态为k′′时，本次更新操的计算过程为： 

（1）、若k＝0，则各参量t₀、k、β和B保持不变，跳出本次更新操作，重新选择参量进行下一次更新操作；否则，则继续步骤（2）； 

（2）、假设此次是第l次尝试更新（步骤（2）至步骤（4）是一次尝试更新），设t₀′、β′分别是本次更新操作后t₀、β可能更新到的状态，且t₀′、β′为k′′维向量，将剔除第l个元素后向量t₀和β的剩余元素按顺序依次分别赋值给向量t₀′和β′的k′′个元素； 

（3）、将y、t₀和β代入公式（2）求解依据y、t₀和β预估的回波信号的能量分布{f(i,β_j,t_0j)}，j是不大于k的任意非负整数，β_j和t_0j分别是向量β和t₀的第j个元素；将y、t₀′和β′代入公式（2）求解依据y、t₀′和β′预估的回波信号的能量分布{f(i,β_m′,t_0m′)}，m是不大于k′′的任意非负整数，β_m′和t_0m′分别是向量β′和t₀′的第m个元素； 

然后将{f(i,β_j,t_0j)}、k和B代入公式（3）求解依据{f(i,β_j,t_0j)}、k和B预估的入射到激光雷达探测器的信号F(i,β,t₀,k,B)，将{f(i,β_j′,t_0j′)}、k′′和B代入公式（3）求解依据{f(i,β_m′,t_0m′)}、k′′和B预估的入射到激光雷达探测器的信号F(i,β′,t₀′,k′′,B)； 

再将F(i,β,t₀,k,B)和y代入公式（4）求解y在信号F(i,β,t₀,k,B)的联合分布概率L(y|t₀,k,β,B)，将F(i,β′,t₀′,k′′,B)和y代入公式（4）求解y在信号F(i,β′,t₀′,k′′,B)的联合分布概率L(y|t₀′,k′′,β′,B)； 

最后将L(y|t₀,k,β,B)、L(y|t₀′,k′′,β′,B)代入公式（6）求解k′′的转移概率α(k′′,k)； 

$\alpha (k^{\prime \prime },k)=L\left(y\right|t_0^{\prime },k^{\prime \prime },{\beta }^{\prime },B)-L\left(y\right|t_0,k,\beta ,B)---\left(6\right)$

（4）、若α(k′′,k)＞0，则将t₀′、β′和k′′的值分别赋值给t₀、β和k，参量B保持不变，完成本次更新操作，在由t₀、k、β和B组成的新的全局近似解（t₀、k、β、B）中随机选择任意参量进行下一次更新操作； 

（5）若α(k′′,k)≤0，则进一步判断是否到达最大尝试更新次数k，若到达最大尝试更新次数k，则各参量t₀、k、β和B保持不变，跳出本次更新操作，重新选择参量进行下一次更新操作；若未到达最大尝试更新次数k，则跳到步骤（2）继续本次更新操作。 

随机选中参量B后的更新方法： 

如果在参量t₀、k、β和B之间任意随机选择的是参量B，则由参量B开始本次更新操作，参量B的更新操作可以采取两种不同的方法进行更新，该两种不同的方法分别基于Metropolis算法和模拟回火算法， 

基于Metropolis算法的更新操作过程为： 

（1）、设B′是B可能更新到的状态，在区间(0，max(y))选择一个服从伽玛分布的随机数赋值给B′，伽玛分布的尺度参量不能过大或者过小，过大或者过小会导致过慢的运算速度，根据经验伽玛分布的尺度参量选择1.5，伽玛分布的形状参量选择为B值； 

（2）、首先将y、t₀和β代入公式（2）求解依据y、t₀和β预估的回波信号的能量分布{f(i,β_j,t_0j)}，j是不大于k的任意非负整数，β_j和t_0j是β和t₀的第j个元素； 

然后将{f(i,β_j,t_0j)}、k和B代入公式（3）求解依据{f(i,β_j,t_0j)}、k和B预估的入射到激光雷达探测器的信号F(i,β,t₀,k,B)，将{f(i,β_j,t_0j)}、k和B′代入公式（3）求解依据{f(i,β_j,t_0j)}、k和B′预估的入射到激光雷达探测器的信号F(i,β,t₀,k,B′)； 

再将F(i,β,t₀,k,B)和y代入公式（4）求解y在信号F(i,β,t₀,k,B)的联合分布概率L(y|t₀,k,β,B)，将F(i,β,t₀,k,B′)和y代入公式（4）求解y在信号F(i,β,t₀,k,B′)的联合分布概率L(y|t₀,k,β,B′)； 

最后将L(y|t₀,k,β,B)、L(y|t₀,k,β,B′)代入公式（7）求解B′的转移概率α(B′,B)； 

α(B′,B)＝L(y|t₀,k,β,B′)-L(y|t₀,k,β,B)    （7） 

（3）、若α(B′,B)＞0，则将B′赋值给B，其它参量t₀、k、β保持不变，完成本次更新操作，在由t₀、k、β和B组成的新的全局近似解（t₀、k、β、B）中随机选择 任意参量进行下一次更新操作； 

（4）若α(B′,B)≤0，则进一步判断是否到达最大尝试更新次数N_B1（预先设置），其中步骤（1）至步骤（3）是一次尝试更新，若达到最大尝试更新次数N_B1，则所有参量t₀、k、β和B保持不变，跳出本次更新操作，重新选择参量进行下一次更新操作；若未到达最大更新次数N_B1，则跳到步骤（1）继续本次更新操作。 

基于模拟回火算法的更新操作过程为： 

（1）、根据模拟回火算法的要求，设定函数T＝g(T₀,u)，1＜T＜T_max，且T同时满足公式（8）和（9）； 

T′＞0    （8） 

(T′)²-2T(T′)²+T′′(T)²≤0    （9） 

公式（8）和（9）中，T′和T′′分别为T关于u的一次导函数和二次导函数，初始化此函数为T＝T₀(T₀＞1)； 

（2）、设T_v＝g(T₀,v)，v是u的范围内的任意值，设置在T_v＝g(T₀,v)下的最大尝试更新次数N_B2（预先设置）； 

（3）、设B′是B可能更新到的状态，在区间(0，max(y))选择一个服从伽玛分布的随机数赋值给B′，伽玛分布的尺度参量不能过大或者过小，过大或者过小会导致较慢的运算速度，根据经验伽玛分布的尺度参量选择1.5，伽玛分布的形状参量选择为B值； 

（4）、首先将y、t₀和β代入公式（2）求解依据y、t₀和β预估的回波信号的能量分布{f(i,β_j,t_0j)}，j是不大于k的任意非负整数，β_j和t_0j是β和t₀的第j个元素； 

然后将{f(i,β_j,t_0j)}、k和B代入公式（3）求解预估的入射到激光雷达探测器的信号F(i,β,t₀,k,B)，将{f(i,β_j,t_0j)}、k和B′代入公式（3）求解预估的入射到激光雷达探测器的信号F(i,β,t₀,k,B′)； 

再将F(i,β,t₀,k,B)和y代入公式（4）求解y在信号F(i,β,t₀,k,B)的联合分布概率L(y|t₀,k,β,B)，将F(i,β,t₀,k,B′)和y代入公式（4）求解y在信号F(i,β,t₀,k,B′)的联合分布概率L(y|t₀,k,β,B′)； 

最后将L(y|t₀,k,β,B)、L(y|t₀,k,β,B′)代入公式（7）求解B′的转移概率α(B′,B)； 

（5）、若α(B′,B)＞0，则将B′赋值给B，其它参量t₀、k和β保持不变，完成本次 更新操作，在由t₀、k、β和B组成的新的全局近似解（t₀、k、β、B）中随机选择任意参量进行下一次更新操作； 

（6）、若α(B′,B)≤0，则计算提高转移概率η(B′,B)，计算方式如公式（10）所示， 

η(B′,B)＝exp(α(B′,B)/T_v)    （10） 

（7）、若η(B′,B)＞ε，ε为在区间(0，1)的一个满足均匀分布的随机数，则将B′赋值给B，其它参量t₀、k、β保持不变，完成本次更新操作，在由t₀、k、β和B组成的新的全局近似解（t₀、k、β、B）中随机选择任意参量进行下一次更新操作； 

（8）、若η(B′,B)≤ε，则进一步判断是否到达在T_v＝g(T₀,v)条件下最大尝试更新次数N_B2，其中步骤（3）至步骤（7）是一次尝试更新，若未到达在T_v＝g(T₀,v)条件下最大尝试更新次数N_B2，则跳到步骤（3）继续本次更新操作；若到达在T_v＝g(T₀,v)条件下最大尝试更新次数N_B2，则按照公式（11）所示方式调节T： 

T_v+1＝g(T₀,v+1)    （11） 

式（11）中，T_v+1为调节后的T； 

（9）、若T_v+1＜T_max，转到步骤（2）继续本次更新操作；若T_v+1≥T_max，则参量t₀、k、β和B均保持不变，跳出本次更新操作，重新选择参量进行下一次更新操作。 

随机选中参量t₀后的更新方法： 

如果在参量t₀、k、β和B之间任意随机选择的是参量t₀，则由参量t₀开始本次更新操作，参量t₀的更新操作基于模拟回火算法进行更新，其步骤为： 

（1）、根据模拟回火算法的要求，设定函数T＝g(T₀,u)，且满足公式（8）（9）；其中，限制u的取值范围使得1＜T＜T_max，T′和T′′分别为T关于u的一次导函数和二次导函数，初始化此函数T＝T₀(T₀＞1)； 

（2）、设T_v＝g(T₀,v)，v是u范围内的任意值，设置在T_v＝g(T₀,v)下的最大尝试更新次数（预先设置）； 

（3）、设k维向量t₀′是t₀可能更新到的状态，在区间(0,i_max)选择k个服从均匀分布的随机数赋值给k维向量t₀′的所有元素； 

（4）、将y、t₀和β代入公式（2）求解依据y、t₀和β预估的回波信号的能量分布 {f(i,β_j,t_0j)}，将y、t₀′和β代入公式（2）求解依据y、t₀′和β预估的回波信号的能量分布{f(i,β_j,t_0j′)}，j是不大于k的任意非负整数，t_0j、β_j和t_0j′分别是向量t₀、β和_t0′的第j个元素； 

然后将{f(i,β_j,t_0j)}、k和B代入公式（3）求解预估的入射到激光雷达探测器的信号F(i,β,t₀,k,B)，将{f(i,β_j,t_0j′)}、k和B代入公式（3）求解预估的入射到激光雷达探测器的信号F(i,β,t₀′,k,B)； 

再将F(i,β,t₀,k,B)和y代入公式（4）求解y在信号F(i,β,t₀,k,B)的联合分布概率L(y|t₀,k,β,B)，将F(i,β,t₀′,k,B)和y代入公式（4）求解y在信号F(i,β,t₀′,k,B)的联合分布概率L(y|t₀′,k,β,B)； 

最后将L(y|t₀,k,β,B)、L(y|t₀′,k,β,B)代入公式（12）求解t₀′的转移概率α(t₀′,t₀)； 

$\alpha (t_0^{\prime },t_0)=L\left(y\right|t_0^{\prime },k,\beta ,B)-L\left(y\right|t_0,k,\beta ,B)---\left(12\right)$

（5）、若α(t₀′,t₀)＞0，则将t₀′赋值给t₀，其他参量k、β和B保持不变，完成本次更新操作，在由t₀、k、β和B组成的新的全局近似解（t₀、k、β、B）中随机选择任意参量进行下一次更新操作； 

（6）、若α(t₀′,t₀)≤0，则计算提高转移概率η(t₀′,t₀)，计算方式如公式（13）所示， 

$\eta (t_0^{\prime },t)=\exp (\alpha (t_0^{\prime },t)/T_v)---\left(13\right)$

（7）、若η(t₀′,t)＞ε，ε为区间(0,1)的一个满足均匀分布随机数，则将t₀′赋值给t₀，其它参量k、β和B保持不变，完成本次更新操作，在由t₀、k、β和B组成的新的全局近似解（t₀、k、β、B）中随机选择任意参量进行下一次更新操作； 

（8）若η(t₀′,t)≤ε，则进一步判断是否到达在T_v＝g(T₀,v)条件下最大尝试更新次数 其中步骤（3）至步骤（7）是一次尝试更新，若未到达在T_v＝g(T₀,v)条件下最大尝试更新次数则转到步骤（3）继续本次更新操作；若到达在T_v＝g(T₀,v)条件下最大尝试更新次数则按照公式（11）所示方式调节T； 

（9）、若T_v+1＜T_max，其中T_v+1是调节后的T，转到步骤（2）继续本次更新操作；若T_v+1≥T_max，参量t₀、k、β和B均保持不变，且跳出本次更新操作，重新选择参量进 行下一次更新操作。 

随机选中参量β后的更新方法： 

如果在参量t₀、k、β和B之间任意随机选择的是参量β，则由参量β开始本次更新操作，参量β的更新操作基于模拟回火算法进行更新，其步骤为： 

（1）、根据模拟回火算法的要求，设定函数T＝g(T₀,u)，且满足公式（8）（9）；其中，限制u的取值范围使得1＜T＜T_max，T′和T′′分别为T关于u的一次导函数和二次导函数，初始化此函数T＝T₀(T₀＞1)； 

（2）、设T_v＝g(T₀,v)，v是u范围内任意值，设置在T_v＝g(T₀,v)下的最大尝试更新次数N_β（预先设置）； 

（3）、设k维向量t₀′是t₀可能更新到的状态，在区间(0,max(y))选择k个服从均匀分布的随机数赋值给β′的所有个元素； 

（4）、首先将y、t₀和β代入公式（2）求解依据y、t₀和β预估的回波信号的能量分布{f(i,β_j,t_0j)}，将y、t₀和β′代入公式（2）求解依据y、t₀和β′预估的回波信号的能量分布{f(i,β_j′,t_0j)}，j是不大于k的任意非负整数，t_0j、β_j和β_j′分别是向量t₀、β和β′的第j个元素； 

然后将{f(i,β_j,t_0j)}、k和B代入公式（3）求解入射到激光雷达探测器的信号F(i,β,t₀,k,B)，将{f(i,β_j′,t_0j)}、k和B代入公式（3）求解入射到激光雷达探测器的信号F(i,β′,t₀,k,B)； 

再将F(i,β,t₀,k,B)和y分别代入公式（4）求解y在信号F(i,β,t₀,k,B)的联合分布概率L(y|t₀,k,β,B)，将F(i,β′,t₀,k,B)和y分别代入公式（4）求解y在信号F(i,β′,t₀,k,B)的联合分布概率L(y|t₀,k,β′,B)； 

最后将L(y|t₀,k,β,B)、L(y|t₀,k,β′,B)代入公式（15）求解β′的转移概率α(β′,β)； 

α(β′,β)＝L(y|t₀,k,β′,B)-L(y|t₀,k,β,B)    （15） 

（5）、若α(β′,β)＞0，则将β′赋值给β，其它参量t₀、k和B保持不变，完成本次更新操作，在由t₀、k、β和B组成的新的全局近似解（t₀、k、β、B）中随机选择任意参量进行下一次更新操作； 

（6）、若α(β′,β)≤0，则计算提高的转移概率η(β′,β)，计算方式如公式（16）所 示， 

η(β′,β)＝exp(α(β′,β)/T_v)    （16） 

（7）、若η(β′,β)＞ε，ε为区间[(0,1)]的一个满足均匀分布随机数，则将β′赋值给β，其它参量t₀、k和B保持不变，完成本次更新操作，在由t₀、k、β和B组成的新的全局近似解（t₀、k、β、B）中随机选择任意参量进行下一次更新操作； 

（8）、若η(β′,β)≤ε，则进一步判断是否到达在T_v＝g(T₀,v)条件下最大尝试更新次数N_β，其中步骤（3）至步骤（7）是一次尝试更新，若未到达在T_v＝g(T₀,v)条件下最大尝试更新次数N_β，则转到步骤（3）继续本次更新操作；若到达在T_v＝g(T₀,v)条件下最大尝试更新次数N_β，则按照公式（11）所示方式调节T； 

（9）、若T_v+1＜T_max，其中T_v+1是调节后的T，转到步骤（2）继续本次更新操作；若T_v+1≥T_max，参量t₀、k、β和B均保持不变，跳出本次更新操作，重新选择参量进行下一次更新操作。 

本发明所述Metropolis算法的原理可参见文献（Journal of Computer and System Science 1998 57(1):20-36），模拟回火算法的原理可参见文献（Europhysics Letters 2007 29(12):2170-2180）。 

本发明的有益效果可以通过以下实验进一步说明： 

本发明在同一实验环境下作了四个实验，实验环境为：使用自行研制的光子计数雷达获取全波形信号y，系统最小时间分辨为24ps，使用本发明方法处理此全波形信号，提取全波形参量中的回波信号数量k和峰值位置t₀，将运算得到的峰值位置t₀和回波信号峰值幅度β代入公式（2）所示的回波信号在的能量分布函数，获得依据峰值位置t₀和回波信号峰值幅度β的回波，即完成全波形信号的信号分解，从而实现激光光斑内的距离的提取。 

根据本发明方法要求，设置的初始条件为：c₁＝10、c₂＝-5、c₁＝-20，τ₁＝5，τ₂＝10，τ₃＝20，伽玛分布概率密度函数f_G的尺度参量a、c和形状参量b、d分别取值为：a=10、b=1、c=1.1、d=100；β＝max(y)；B＝0；k=1；选择（0,i_max）中max(y)所在位置赋值给t₀；t₀和β采用本发明所述的基于模拟回火算法的更新操作过 程，其中，T＝T₀×1.2^v，T₀＝2，通过控制v改变T，并且在任意T下的尝试更新次数均是20，即和N_β均为20；k和B采用本发明所述的基于Metropolis算法的更新操作过程，并且尝试更新次数均为20，即N_k和N_B1均为20；设置所有参量更新次数总和为600；用matlab工具编写程序，处理相关数据。 

实验一： 

本实验探测目标是毛玻璃及其后与其平行放置的纸箱。探测目标距离激光雷达约为40米。调整纸箱与毛玻璃的间距分别为：5厘米、8厘米、10厘米、20厘米、30厘米、50厘米以及100厘米，获得7组全波形实验数据。将使用本发明方法得到的t₀和β的对应元素代入公式（2），还原出回波波形，全波形和还原出的回波波形如图2所示。图2中，带毛刺的曲线为全波形号，其它两个曲线是通过本发明方法还原出的目标回波信号。图2（a）是间距为5cm的全波形信号的信号分解图、图2（b）是间距为8cm的全波形信号的信号分解图、图2（c）是间距为10cm的全波形信号的信号分解图、图2（d）是间距为20cm的全波形信号的信号分解图、图2（e）是间距为30cm的全波形信号的信号分解图、图2（f）是间距为50cm的全波形信号的信号分解图、图2（g）是间距为100cm的全波形信号的信号分解图。 

在图2可以看出，采用相同的光子计数激光雷达，针对不同距离的探测目标，相同的算法初始设置，均能实现特征距离的提取，实现信号的自动化处理，这是其它一些算法所不具有的优势。LM算法的非线性最小二乘波形拟合方法只适用于单一特性的场景；EM算法及其改进算法的全波形参量最大似然值估计方法的运算结果收敛性强烈依赖于初始值的预设，无法实现完全的非人工干预全波形参量提取，有时存在失效的风险。 

从图2（a）、图2（b）、图2（c）、图2（d）中可以看出，毛玻璃与纸箱间距越小，波形重叠越严重，其分解越困难，这也证明本算法在精度上的优势。从图2（e）、图2（f）、图2（g）可以看出，毛玻璃与纸箱间距越大，回波之间的重叠越小，算法运算结果计算获得的的峰值位置t₀与粗略观察到的探测到的目标信号峰值位置非常接近，这与实际的距离相当接近，可以证明本算法的准确性。 

实验二： 

为了进一步检验本发明处理全波形信号的能力，本实验降低激光发射能量以获取低信噪比的全波形信号。本实验的探测目标是毛玻璃及其后的与其平行放置的纸箱。探测目标距激光雷达约40米。调整纸箱与毛玻璃的间距粗略测量分别为10cm和20cm。经 本发明对全波形信号的处理，将通过运算得到t₀和β的对应元素代入到公式（2）得到分解出的回波波形，将同一组数据分解出的波形与原始波形显示在一张图，其中带毛刺的曲线为全波形号，其它两个曲线是通过本发明方法还原出的目标回波信号。图3（a）是间距为10cm的低信噪比全波形信号的分解图，通过运算结果得到其间距为9.76cm；图3（b）是间距为20cm的低信噪比全波形信号的分解图，通过运算结果得到其间距为21.6cm。 

实验三： 

将本发明处理激光雷达通过点扫描得到的面阵目标的全波形信号，还原出目标在距离上的分布。如图4所示，探测目标为依次平行放置的毛玻璃、N形字母及纸箱，纸箱及N形字母与毛玻璃间的距离分别为20cm和10cm（图4（a）为目标的俯视图，图4（b）为目标的主视图）。每一束激光均穿过毛玻璃照射到N形字母或者背后的纸箱上，处理这样的每一束激光回波信号，并重建三维形态如图4（c）。图4(c)给出了还原的目标形态，图中黑色区域均是激光可以照射到的位置，黑色块状区域对应实际的毛玻璃，中间的黑色N形字母即N形字母硬纸板，剩下的黑色区域即纸箱上激光可以照射到的区域，后方的白色N形字母即纸箱上激光照射不到的区域，从左侧的距离坐标可以看出还原出的目标在距离上的信息，与真实值相近。 

实验四： 

本实验将本发明方法与RJMCMC方法做了比较实验，进一步对本发明进行评价。两种方法在相同条件下对同样的试验数据进行处理。两种方法对实验一中获取的7组试验数据分别处理100次（主要是提取两个距离特征），RJMCMC算法按照固定顺序循环更新所有参量，每个参量更新600次，本发明的所有参量的更新次数总和为600。分别统计平均运行时间以及100次运行结果，如表1所示统计了本发明方法与RJMCMC方法平均运行时间以及100次运行结果（各间距均值，及标准差）。在间距均值和间距标准差这两个指标上，RJMCMC方法运算结果的稳定性要比本发明好；而在程序运行时间均值上的比较，本发明在速度上具有巨大的优势。 

表1 本发明与RJMCMC方法在距离提取及运算时间的统计表