激光雷达是一种以激光作为光源的主动式遥感技术。激光光源发射一束短脉冲激光，经过目标的后向散射后，使用探测器接收这一回波脉冲，通过高精度的时间测量系统，获取出射脉冲和回波脉冲之间的飞行时间，即可借此计算出目标与测量系统之间的距离(Wagner 等，2006)。而相比普通的激光雷达，全波形激光雷达通过数字化采样获取了回波脉冲的完整波形(Leeuwen 等，2011)，从中不仅可以获取目标的距离信息，还可以提取出激光足印内目标的几何与物理特征，如粗糙度(Shi 等，2013)、倾斜角(Mahoney 等，2014；Li 等，2012)、森林冠盖结构(Kotchenova 等，2003；Hovi和Korpela，2014；Bye 等，2017)、水深(Pan 等，2015)等，此外，全波形回波数据也可以用来辅助地物分类(Koenig和Höfle，2016)。

全波形回波脉冲可以认为是在一个激光脚印中，不同距离的多个目标面散射后产生的单回波脉冲的叠加。因此，可以通过分解全波形回波的方式获取对应不同目标面的单脉冲，进而由分析这些单脉冲的参数得到相应目标面的距离与其他特征信息(Mallet和Bretar，2009)。分解全波形回波的方法通常包含对原始波形的滤波与拟合，滤波的同时需要进行噪声水平的估计，拟合过程则包括初始参数估计和拟合优化两部分(Hofton 等，2000)。

当前常用的全波形回波滤波法主要有低通滤波、维纳滤波(Jutzi和Stilla，2006)、高斯平滑滤波(Hofton 等，2000；Xu 等，2016)和小波分解滤波(Zhou 等，2013；Li 等，2017)等。低通滤波对于频率在截止频率附近的噪声的滤波效果较差。维纳滤波需要对输入信号与预期输出信号之间的互相关向量的预估计，实现难度较大。而高斯平滑滤波存在一个核心问题，即如何选取合适的高斯核函数的脉宽参数，脉宽的变化会对滤波的结果产生极大的影响。同样，小波分解滤波也需要人为选择小波基函数，而小波基函数的选择目前缺少一般性的选择原则。因而，本文使用了可变阈值的经验模态分解滤波法 (EMD-soft)对全波形回波进行滤波，并同时估计其噪声水平。

目前广泛使用在全波形回波拟合优化上的算法主要是LM (Levenberg-Marquard)算法(Qin 等，2012；Ma 等，2017；Marquardt，1963)与期望最大化EM(Expectation Maximization)算法(Soederman 等，2005；Dempster 等，1977)，两种方法都具有较高的拟合速度与拟合精度，但也面临着相同的问题，即两种方法都对于初始参数的准确性具有非常高的要求，如果输入的初始参数与期望值差距太大，会导致拟合收敛到局部最小值，而不是所需的全局最小值，甚至导致拟合发散。因此，本文着重完善了LM算法初始参数估计过程，得到更加准确的初始参数，进而提高拟合精度。

全波形激光雷达系统中使用的激光脉冲近似为高斯函数，因此本文选择高斯函数作为分解模型，此时获取的全波形回波可以表示为式(1)所示的高斯模型

式中，N为回波中的高斯脉冲数，a_i、t_i、f_i分别为第i个高斯脉冲的幅值、位置与半高宽， $n(t)$ 为噪声带来的残余误差，上述表达式中 ${f_i}$ 可以直接对应高斯函数的半高宽，因此使用该表达式便于参数的对应。式中，每一个高斯脉冲都由相应的目标面的散射特性决定。

由于全波形回波在获取过程中，会加入由背景光、光电探测器和模数转换器等带来的噪声，在拟合与分解之前需要对采集到的原始波形进行滤波，这里我们采用的是EMD-soft滤波法。

经验模态分解EMD(Empirical Mode Decomposition)理论认为，任何信号都是由一些不同的本征模态函数IMF(Intrinsic Mode Function)组成的

式中，X(t)表示原始信号，IMF_j(t)表示分解出来的n个本征模态函数，r_n(t)表示分解后的余项，也即信号的单调趋势项。每个IMF必须符合以下两个条件：两个零点之间只能存在一个极值点，且由其局部最大值和局部最小值定义的上下包络必须均值为零。通过反复从信号中减去上下包络的均值，直到残量满足IMF的条件的方法获取IMF，而后对被减去的部分重复上述步骤，获取更多的IMF直到最终剩余的单调趋势项，即完成了整个EMD过程(Huang 等，1998；Li 等，2014)。

从滤波的角度来看，EMD过程中，逐个提取出的IMF分别包含了不同频率的分量，先提取出的IMF包含了较高频率，后提取出的IMF包含较低频率，最终剩余一个单调趋势项，因而通过舍去含不同频率分量的IMF，即可实现不同的滤波效果，由于全波形回波中噪声的频率远高于高斯脉冲的频率，本文采用舍弃较高频率的IMF来实现低通滤波的效果，即

式中， ${F_0}(t)$ 为滤波后的波形，舍去前p个IMF。

然而直接舍去前p个IMF有可能会丢失部分有用的信息，因此，需要对EMD过程进行一些改进，从而将可能的有用信息从被舍去的IMF中提取回来，由于该改进方法采用了一个可变的阈值，或者说柔性的阈值，因而称之为EMD-soft滤波法。在EMD过程结束后，通过下述方法从被舍去的IMF中提取有用信息 ${s_j}(t)$ (Boudraa 等，2013)

式中， $j \leqslant p$ ， ${\sigma _j}$ 为对 ${{\rm{IMF}}_j}$ 的噪声估计，L为 ${{\rm{IMF}}_j}$ 的数组长度。将提取出的 ${s_j}(t)$ 加到 ${F_0}(t)$ 中即得到了经过EMD-soft滤波法处理后的回波波形 $F(t)$

滤波完成后，可以对噪声水平进行估计

将 ${\rm{Noise}}(t)$ 的均值与标准差记为 ${\mu _n}$ 与 ${\sigma _n}$ 。

由于后续的拟合优化需要提供较为准确的初始参数，即分解模型中a_i、t_i、f_i的估计，因而在优化之前需要对滤波后的回波进行初始参数的估计。通过数值微分获取 $F(t)$ 的一阶与二阶导数，并分别由一阶、二阶导数的零点得到 $F(t)$ 的极值点与拐点。

由于滤波并不能完全去除原始波形中的噪声，残余的一些噪声可能会带来错误的极值点与拐点，所以为了避免这一问题，需要设定一个噪声阈值 $t{h_n}$

${\mu _n}$ 与 ${\sigma _n}$ 为上一节中得到的噪声的均值与标准差。而后，在获取极值点与拐点的过程中，如果该点的幅值小于 $t{h_n}$ ，则认为该点是一个错误的极值点或拐点，将其忽略。此外，在计算一阶与二阶导数的过程中，还可以对导数进行EMD-soft滤波，以进一步减少残余噪声带来的误判。

全波形回波中的高斯脉冲，存在相互独立与相互交叠两种情况，相互独立的高斯脉冲可以通过极值点来确定其数量、幅值与位置，而相互交叠的高斯脉冲在一些情况下由于两个脉冲距离过于接近，其中一个脉冲的极值点被另一个脉冲掩盖，则需要通过其他的方法估计其初始参数。图1是一个含有4个高斯脉冲的全波形回波，4个高斯脉冲分别记为 ${G_1}$ 、 ${G_2}$ 、 ${G_3}$ 、 ${G_4}$ ，可以看到， ${G_1}$ 与 ${G_2}$ 为独立高斯脉冲，而 ${G_3}$ 与 ${G_4}$ 则相互交叠， ${G_3}$ 的极值点被 ${G_4}$ 掩盖，因而整个回波中只能得到3个极值点，将高斯脉冲数N暂记为3。

图 1 一个含有4个高斯脉冲的全波形回波 Figure 1 A full-waveform echo with four Gaussian pulses

对于独立脉冲例如 ${G_1}$ ，其对应的极值点为 ${P_1}({t_{p1}}, {a_{p1}})$ ，则 ${a_{p1}}$ 、 ${t_{p1}}$ 可以作为 ${G_1}$ 幅值与位置的初始参数，而对于独立脉冲，极值点 ${P_1}$ 左右两侧的单调区间均只有一个拐点，即 ${I_1}({t_{I_1}}, {a_{I_1}})$ 与 ${I_2}({t_{I_2}}, {a_{I_2}})$ ，则其在时间轴上的差 ${f_1} = {t_{I_1}} + {t_{I_2}}$ 可以直接作为 ${G_1}$ 半高宽的初始参数。

对于相互交叠的脉冲例如 ${G_3}$ 与 ${G_4}$ ，只能得到一个极值点 ${P_3}$ ，但在 ${P_3}$ 左侧的单调区间中有3个拐点 ${I_5}$ 、 ${I_6}$ 、 ${I_7}$ ，我们把这种在某一极值点一侧的单调区间中有多于1个拐点的情况作为有脉冲被交叠掩盖的标志，此时将高斯脉冲数N的值改为N+1，并由 ${P_3}({t_{p3}}, {a_{p3}})$ 得到 ${G_3}$ 的幅值与位置的初始参数。 ${G_3}$ 的半高宽与 ${G_4}$ 的初始参数则需要进一步的分析：

对于一个高斯函数

任意横坐标 $t$ 与中心位置 ${t_0}$ 的距离为 $\Delta t = \left| {t - {t_0}} \right|$ ，已知一高斯函数上任意两点的横坐标与中心位置的距离 $\Delta {t_1}$ 、 $\Delta {t_2}$ ，及其幅值比 ${a_1}/{a_2}$ ，可以求得该函数的半高宽 ${f_0}$

即

因而对于 ${G_3}$ ，其左侧拐点受到 ${G_4}$ 的影响，位置并不在半高宽处，但可以通过上述公式，由 ${P_3}$ 两侧与其最接近的两个拐点，即 ${I_7}({t_{I_7}}, {a_{I_7}})$ 与 ${I_8}({t_{I_8}}, {a_{I_8}})$ 来确定，独立脉冲的半高宽估计，是此方法在所取两个拐点恰好为半高宽对应的两点时的特殊情况。

而对于 ${G_4}$ ，由于 ${I_6}$ 并不一定恰好处在其极值点所在的位置，而且在有些情况下，由于 ${G_4}$ 与 ${G_3}$ 距离过近，会出现只能获取 ${I_5}$ 一个拐点的情况，因而对于这种极值点被掩盖的脉冲，需要采用另一种方法进行初始参数估计。认为 ${I_5}({t_{I_5}}, {a_{I_5}})$ 为 ${G_4}$ 的一个半高宽对应点，则幅值的初始参数可以由式(14)进行估计

而后在全波形回波中找到幅值最接近 ${a_4}$ 的点，其横坐标即可作为 ${G_4}$ 的位置初始参数，进而 ${I_5}$ 与该点在横轴上距离的两倍可作为 ${G_4}$ 的半高宽初始参数。

通过上述的所有方法，即可实现对回波中所有高斯脉冲的初始参数估计。

在经过前述的滤波与初始参数估计后，获取了回波中高斯脉冲的数量N以及每一个高斯脉冲对应函数表达式的幅值、位置和半高宽，即a_i、t_i、f_i的初始估计，利用这些参数对滤波后的回波 $F(t)$ 进行非线性最小二乘拟合，可以获取更加准确的高斯脉冲参数。

Levenberg-Marquardt(LM)法是一种被广泛用于激光雷达回波数据拟合优化的算法，它是对Gauss-Newton(GN)法的改进，通过加入参数 ${\mu _k}$ ，将GN法与最速下降法结合起来，在所拟合的数据接近收敛点时，算法趋于GN法，拥有GN法的快速收敛性，而当所拟合的数据远离收敛点时，算法趋于最速下降法，拥有最速下降法的强适应性，使LM法可以更加快速有效的完成拟合优化(Marquardt，1963)。

LM法的迭代公式如下

式中，F为N个高斯函数相加组成的拟合目标函数，N为之前获取的高斯脉冲数，J为F的Jacobian矩阵。当目标函数与回波 $F(t)$ 之差达到设定的阈值或迭代次数达到设定的上限，迭代停止，即可得到经过拟合优化后精度更高的高斯脉冲函数表达式，完成回波的分解。

为了检验本分解方法的有效性与分解精度，使用高斯函数之和并加入高斯白噪声的方法对全波形回波进行仿真，采样率与采样时间分别设置为5 GHz和199 ns，与后面实验中采集实测回波时的参数相同。

仿真的函数表达式如下

式中，t为采样时间，N为回波中高斯脉冲的数量，设为1—4的随机整数，a_i、t_i、f_i为第i个高斯脉冲的幅值、位置与半高宽，分别设为3—30的随机整数、40—160 ns的随机数和10—20 ns的随机整数， $n(t)$ 为高斯白噪声，仿真中按照信噪比分别为15 dB、25 dB、35 dB、40 dB设置4组，每组的仿真总次数为500次。图2是一包含4个高斯脉冲，信噪比25 dB的仿真全波形回波。

图 2 仿真全波形回波 Figure 2 A simulated full-waveform echo

分解实验中，将本方法与另外两种高斯分解方法进行对比，两种方法分别记为Hofton法(Hofton 等，2000)、Duong法(Duong，2010)，通过相同的评价参数比较3种方法的分解效果。

定义一些评价参数对仿真试验的结果进行评价。分解成功率S

式中， ${N_0}$ 为该组仿真的总回波数，此处为500， ${N_S}$ 为分解成功的回波数，这里定义满足分解后高斯脉冲数量与设定值相同，且每个高斯脉冲的位置与设定值误差小于1 ns这两个条件的回波分解为分解成功的回波。

其他评价参数有幅值、位置、半高宽与设定值之差的均值μ_a、μ_t、μ_f及其标准差σ_a、σ_t、σ_f，4种方法的分解结果如表1所示，分别列出了4种方法对于上述4组信噪比不同的仿真回波的分解结果评价参数值。

表 1(Table 1)

表 1 不同信噪比下的仿真评价参数 Table 1 Evaluation parameters of simulation under different SNRs 分解方法 信噪比/dB S/% μa μt /ns μf /ns σa σt /ns σf /ns Hofton法 15 46.5 –2.04 0.806 0.424 3.85 0.967 0.988 25 62.3 –1.32 0.592 0.309 2.27 0.411 0.610 35 77.6 –0.905 0.498 0.228 1.68 0.265 0.291 40 87.8 –0.356 0.086 0.122 0.674 0.120 0.144 Duong法 15 48.6 –1.75 0.824 0.324 3.70 0.704 1.374 25 62.9 –1.10 0.601 0.231 2.46 0.323 0.747 35 79.8 –0.851 0.502 0.170 1.69 0.102 0.409 40 88.2 –0.322 0.090 0.112 0.623 0.064 0.187 本文方法 15 78.3 –0.347 0.016 0.242 0.439 0.138 0.343 25 79.4 –0.175 –0.004 0.170 0.241 0.119 0.257 35 87.5 –0.107 0.001 0.091 0.112 0.033 0.085 40 91.5 –0.096 0.001 0.078 0.066 0.010 0.032

表 1 不同信噪比下的仿真评价参数 Table 1 Evaluation parameters of simulation under different SNRs

由表1中的仿真结果可以看出，本文方法分解成功率在78.3%以上，位置与半高宽分解误差在0.1 ns量级，对应测距误差为15 mm，即厘米量级。随着信噪比的提高，分解的成功率与精度等都会有相应的提升，当原始波形信噪比达到40 dB时，分解成功率可以达到91.5%，位置与半高宽的分解误差也会达到0.01 ns量级，对应测距误差达到毫米量级，这也说明实际使用中分解方法的有效性与精度同时也受到硬件条件的较大影响。而对比其他两种方法，本方法的分解成功率有了极大的提高，同时分解误差也有了明显的下降。从表格中的数据也可以看出，本方法在信噪比较低的情况下具有更明显的优势，这一方面是由于EMD-soft滤波算法较其他两种算法使用的高斯平滑滤波的滤波效果更好，同时也是得益于本方法中初始参数估计算法相较于其他两种算法，对于交叠脉冲的参数估计以及参数估计的准确性都有较大的提高。

使用实验室自制的全波形激光雷达实验系统获取实测的全波形回波，并使用本方法进行分解。所用的激光雷达系统结构如图3所示，脉冲激光器波长为1064 nm，发射脉冲半高宽为15 ns，脉冲出射后分为两束，一束由光电探测器(PIN)T1接收作为出射波，另一束作为探测光经目标面散射后由雪崩二极管(APD)探测器T2接收作为回波，并使用示波器或采集卡采集全波形数据。

图 3 激光雷达实验系统 Figure 3 LiDAR experimental system

实测回波的分解实验通过分解含有两个高斯脉冲的回波，计算出两个高斯脉冲对应的目标面间的距离，并对比实际目标面间距，来验证分解方法的精度。

在距实验系统一定距离处设置两块目标板，使其各有一部分处于同一激光脚印中，调整两块目标板的距离，分别在2.15 m、2.30 m、2.45 m、3.20 m、4.80 m、6.20 m距离下各采集30组回波数据，图4是实验中目标板的分布图，图5是上述6组距离下的实测波形，图5(a)—(f)分别对应这6组距离，随着两目标板距离的增加，回波中的两个高斯脉冲在时间轴上的距离也逐渐增加，由交叠状态逐渐分离。将回波数据使用本分解方法以及前述的另外两种方法进行分解，得到两个高斯脉冲的时间差并换算为距离值，30组回波的距离平均值 $\mu $ 与标准差 $\sigma $ 如图6所示。

图 4 实验中目标板的分布图 Figure 4 Distribution of the target boards setting in experiment

图 5 实测波形 Figure 5 Measured waveforms

图 6 不同距离下分解回波得到的距离平均值与标准差(红色条形和黑色线分别代表回波分解得到的距离值的平均值 $\mu $ 和标准差 $\sigma $ ，蓝色线代表实际距离) Figure 6 Means and standard deviations of the distances obtained from decomposing echoes (Red bars and black lines represent the means $\mu $ and standard deviations $\sigma $ of distances derived from decomposition of the echoes respectively, blue lines represent the real distances)

由图6的实验结果可以看出，本方法的分解结果较其他两种方法更接近实际值，且多组回波分解的标准差也明显小于其他两种方法，证明了本方法分解误差更小，稳定性更好，实验中最大分解误差小于0.1 m，也在厘米量级，经过对噪声水平的估计，实测波形的信噪比约为27 dB，而误差与25 dB信噪比下的仿真结果处于同一数量级，考虑到实际实验中存在硬件系统或环境影响导致的波形变型等问题会影响到拟合效果，可以认为实验结果与仿真结果相符。

此外，3种算法在分解实际波形时的分解误差相比于分解仿真波形时，均有了一定的增大，其产生的原因主要在于回波实验中的APD探测器，APD管与后续跨阻放大电路中用到的运放等元器件对于脉冲信号的响应均会有一定的延迟，而该延迟时间相对于标称值存在一定的误差。查阅所用APD管与运放等的数据手册可知，其延迟时间的误差之和约为0.3 ns，即该误差会为单次测距结果加入约4.5 cm的误差，这一结果也符合实测数据与仿真数据分解后误差的增大量。该误差仅为单次测量中可能产生的偏差，是一个在真值附近浮动的随机误差，在实际测量中，在测量速度允许的范围内，通过多次测量取平均值的方法即可极大的消除其影响，上述的实测实验中，30次测量平均后均值存在的误差已经与仿真实验结果相近，证明了多次测量取均值可以消除该误差。

本文实现了一种基于EMD-soft滤波与LM优化算法的全波形激光雷达回波分解方法，通过对仿真波形和实测回波波形的分解，验证了本方法的有效性，同时得到了本方法将回波分解后的精度在0.1 ns量级，对应的距离精度在厘米量级，随着波形信噪比的提高，分解精度也会相应提高，达到毫米量级，实现了对全波形激光雷达回波的高精度分解，为后续进一步的回波分析提供了可靠的原始数据，将在遥感、测绘等生产与科研领域中发挥非常重要的作用。

相较于已有的回波分解方法，本方法采用的滤波算法不需要针对不同的回波数据设定不同的滤波参数，具有更好的普适性与可靠性；提出的初始参数估计方法提高了拟合初始参数的精度，进而提高了LM算法的拟合效果。与其他两种高斯分解方法的对比也直接证明了本方法在分解成功率上具有显著的优势，在分解精度上也有了明显的提升。

本方法可以从全波形回波数据中获取所有独立高斯脉冲的幅值、位置、半高宽等参数，其中除文中已提取到的距离信息外还有很多的其他的信息如反射率、粗糙度、面型等可以进行进一步的提取与分析。回波分解方法是回波分析的基础，接下来的工作中，我们将基于本分解方法与现有的实验系统，尝试提取出更多的特征信息，建立目标特征与回波参数的对应函数关系，实现回波与复杂面型之间的反演，并通过实验系统进行验证。