CN113866736B - 激光雷达回波信号高斯拐点选择分解方法 - Google Patents

Abstract

本发明公开了激光雷达回波信号高斯拐点选择分解方法，包括：(1)收集相应的激光雷达回波信号数据；(2)对激光雷达回波信号数据进行预处理；(3)构建高斯模型并求解；(4)高斯模型精度验证。本发明操作简便且效率高，能够有效的分解复杂激光雷达回波信号，特别是只存在一个峰值的复杂激光雷达回波信号，为激光雷达回波信号分解提供了便捷的分解算法，在激光雷达遥感探测中具有重要意义。

Description

激光雷达回波信号高斯拐点选择分解方法

技术领域

本发明涉及激光雷达领域，特别涉及到激光雷达数据处理，提出了激光雷达回波信号高斯选择拐点分解方法。

背景技术

激光雷达已经广泛应用于测绘、地质、城市规划、农业开发、水利水电、土木工程、环境检测等领域。为了获得高精度的激光雷达点云数据，有效的对全波形激光雷达回波信号进行波形分解，成为全波形激光雷达数据处理的一个热点。

全波形激光雷达回波信号过程是：激光从激光器中发射而出，在激光照射到地物后，接收器接收到地物的反射光并将数据进行存储，再将存储的数据转换为数字的形式，所得到的数据就是全波形激光雷达回波信号数据。波形分解需要计算全波形激光雷达回波信号在时间轴上的相应幅值位置，从而得到精确的时间间隔，通过光速与时间间隔相乘就可以得到激光器到目标物的距离。因此，计算准确的时间间隔是全波形激光雷达回波信号数据处理工作的重中之重。

近二十年来，我国许多科研人员从事对激光雷达回波信号分解的研究，主要方法有反卷积法和数学模拟法，由于反卷积法容易受到噪声影响，并且其分解之后的信号变窄，无法通过波形的形状、宽度、强度和偏度提供了更多的物体信息，所以常使用数学模拟法应用于全波形激光雷达回波信号的数据处理中。在数学模拟法中最常用的方法就是高斯分解法。高斯分解法的原理是通过多个高斯函数进行建模，模拟全波形激光雷达回波信号。高斯分解法需要估计回波信号中高斯分量的个数和高斯分量的参数，再对拟合回波信号的高斯估计参数进行优化。武汉大学马洪超等使用拐点的方法计算高斯参数，拐点法计算全波形激光雷达回波信号的拐点，再通过拐点和峰值点计算出高斯函数需要的三个参数，从而完成参数估计。拐点法虽然能够寻找到高斯参数，但是并不适用于分解复杂的回波信号，在参数拟合中容易使得高斯波形产生负向震荡的波形，导致距离提取时数值失真。南京大学徐帆等提出峰值拐点的方法计算高斯参数。峰值拐点法是寻找所有波形的峰值，并计算所有拐点位置，通过拐点和峰值判断高斯波形的个数。但是如果有多个回波信号叠加的复杂回波信号，这些回波信号非常接近并且只产生一个峰值点时，峰值拐点法就会分解出一个高斯波形。

针对上述全波形激光雷达波形分解方法遇到的问题，本发明公开了高斯拐点选择分解方法，该方法能够解决分解复杂回波信号的问题，并且能够处理只有一个峰值点并存在多个拐点的复杂叠加波形。

发明内容

为了达到上述的目的，本发明提供了如下技术方案：

激光雷达高斯拐点选择分解方法，包括如下步骤：

(1)收集相应的激光雷达回波信号数据；

(2)对激光雷达回波信号数据进行预处理；

(3)构建高斯模型并求解；

(4)高斯模型精度验证。

优选地，步骤(1)中提供全波形激光雷达回波数据。

优选地，步骤(2)所述回波信号预处理包括Savitzky-Golay算法进行信号去噪与五点三次平滑法处理回波信号。

优选地，步骤(3)所述的高斯模型为：

A_i，第i个高斯函数幅值；

μ_i，第i个高斯函数的中心位置；

σ_i，第i个高斯函数的标准偏差；

Noise，背景噪声；

优选地，对多高斯函数模型进行求解时采用高斯拐点选择分解方法，根据峰值和拐点数目的计算出高斯参数，将各参数带入高斯函数模型，得到多组高斯分量。接着将不符合要求的高斯分量进行剔除。然后利用传统的Levenberg-Marquarelt算法对高斯分量的标准偏差σ_i进行优化。再利用全局Levenberg-Marquarelt算法，优化所有的高斯分量的参数。最后将优化后高斯分量区分为明显的高斯分量，隐藏的高斯分量和无效的高斯分量。

优选地，步骤(4)中采用原始波形数据与计算得到的最优参数高斯函数模型的均方根误差，平均拟合优度R²，波形分解数，波形分解类型来实现高斯模型精度的验证。

本发明采用上述的技术方案，具体有以下的增益效果：

本发明利用获取到的全波形激光雷达回波信号数据，能够根据不同拐点数判断出全波形激光雷达回波信号的幅值，中心位置和标准偏差，能够为激光雷达回波信号数据处理提供相应的较为准确的分解方法，在全波形激光雷达回波信号分解应用中具有重要的意义。

附图说明

图1是高斯拐点选择分解算法的流程

图2高斯拐点选择分解算法波形半宽选择原理

图3是高斯拐点选择分解算法的分解原理

图4a是部分波形分解结果

图4b是部分波形分解统计的结果

图4c是回波信号高斯拐点选择分解算法分解的RMS统计图

图4d是回波信号高斯拐点选择分解算法分解的波形类型数统计图

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下举出优选实施例，对本发明进一步详细说明。然而，需要说明的是，说明书中列出的许多细节仅仅是为了使读者对本发明的一个或多个方面有一个透彻的理解，即便没有这些特定的细节也可以实现本发明的这些方面。这些实例及其说明并不限制本发明的权利要求，任何在本发明的启示下得出的与本发明相同或想近似的方法，均在保护范围之内。

实施例1

以美国密西西比州伊萨奎纳县(34°07′,N,88°28′,E)的16384条回波信号数据为例，本实例提供了高斯拐点选择分解算法，包括如下步骤：

(1)收集全波形激光雷达回波信号，例如，从LVIS网站获取https:// lvis.gsfc.nasa.gov/Data/Maps/GEDI2019Map.html；

(2)对步骤1的数据进行预处理；

(2-1)采用Savizkg-Golag算法对原始数据进行信号去噪；

(2-2)采用五点三次平滑算法再对原始数据进行平滑；

(2-3)计算原始数据首尾各10％的采样点的平均值作为背景噪声Noise；

(2-4)使用以下公式计算迭代停止的阈值σ_th:

σ_th＝3σ_n+Noise (3)

Y(t)为原始数据，f(t)为滤波后的采样点数据，n为采样点数。

(3)使用高斯拐点选择分解算法，找到高斯参数的估计值

高斯分解实际上是使用由多个高斯分量组成的高斯模型来尽可能地拟合回波信号。高斯拐点选择分解算法利用回波信号峰值附近拐点的个数来判断波形半宽的位置，选择合适的波形半宽对回波信号进行迭代分解，从而得到高斯分量。具体步骤如下：

(3-1)计算采样点的幅值，并将局部最大的幅值作为公式(1)中的A，并将局部最大的幅值对应的采样点的时间位置t作为中心位置当满足下了公式时就是局部最大的幅值。

f(t)＞f(t+1)＞f(t+2)&&f(t)＞f(t-1)＞f(t-2) (4)

其中t表示时间。

(3-2)接着需要对振幅A左右两侧的拐点数目进行判断，使用以下式子寻找波形拐点。

(f(t-2)+y(t)-2f(t-1))(f(t-1)+f(t+1)-2f(t))＜0 (5)

(3-3)当寻找到拐点之后会出现下面的4种情况，根据不同的情况选择不同的半波宽。

当幅值A左侧的拐点数等于1，右边拐点数大于1的情况。计算幅值A的一半为A/2，在幅值A左侧的采样点中找到最接近A/2的采样点，并将这个采样点对应点时间位置记为tgl，波形半宽σ则由公式(6)得到。如图2(a)所示。

当振幅A左侧的拐点数大于1，右边拐点数等于1的情况。计算幅值A的一半为A/2，在幅值A右侧的采样点中找到最接近A/2的采样点，并将这个采样点对应点时间位置记记为tgr，波形半宽σ则由公式(6)得到。如图2(b)所示。

当幅值A左右拐点数都大于等于1的情况，在这种情况下存在一个回波或者多个叠加的回波，波形振幅最大值在叠加回波的中间，只存在一个波峰。首先找到接近振幅A左侧的两个拐点对应的时间位置，再找到接近振幅A右侧的两个拐点对应的时间位置。计算左侧两个拐点的时间位置的平均值a₁，找到最接近a₁的采样点的时间位置记为S₁。计算振幅A右侧两个拐点的时间位置的平均值a₂，找到最接近a₂的采样点的时间位置记为S₂。寻找幅值A两侧最接近幅值A一半的采样点，记录其对应的时间位置记为S₃。然后将S₁,S₂和S₃与中心位置μ进行相减并取绝对值，去掉最大值和最小值后所得到的数值记为|μ-mid_dis|。最后将|μ-mid_dis|带入公式(6)中进行计算得到σ。如图2(c)所示。

当上述三种情况没有出现。寻找幅值A两侧最接近幅值A一半位置的采样点，记录其对应的时间位置为dis，将dis带入公式(6)中计算得到σ。

当找到初始参数A、μ、σ时就可以带入公式(1)中，然后使用滤波后的回波信号和高斯函数进行相减，为了得相减后得到振幅值非负，使原始波形信号振幅低于0mv的部分等于0mv。完成迭代(3-1)-(3-3)寻找剩下的高斯参数，随着迭代振幅越来越小，当振幅小于阈值时停止迭代。

公式(6)中LINF表示幅值A左侧的拐点数，RINF表示幅值A右侧的拐点数。

(3-4)波形分解结果如附图3所示，分解顺序为(a)—(j)，(k)为分解结果。

(4)由于迭代分量的分解会产生一些无意义的高斯分量，因此在迭代时和迭代结束之后要对这些高斯分量进行剔除。任意两个高斯分量如果小于一个发射脉冲宽度，则可以将波峰值较小的分量进行剔除。

(5)对参数的波形半宽进行优化，传统Levenberg-Marquardt算法(LM算法)将幅值A_i和中心位置μ_i当成常数，迭代优化σ_i。

p＝p⁽⁰⁾+[H(t,p⁽⁰⁾)+λE]^-1J^T(t,p⁽⁰⁾)[f(t)-f(t,p⁽⁰⁾)] (7)

其中f(t,p)为待定系数p₁,p₂,p₃,...,p_n的函数，p⁽⁰⁾表示初始参数p⁽⁰⁾＝(p₁,p₂,p₃,...,p_n)，H(t,p)表示f(t,p)的海塞矩阵，J^T(t,p)表示f(t,p)的转置雅克比矩阵，f(t)表示为预处理后的数据，λ为阻尼系数。若||p-p⁽⁰⁾||相差很大，则将p作为新的p⁽⁰⁾，直到||p-p⁽⁰⁾||相差很小或者循环达到所设定的值时，迭代结束。

(6)采用全局LM算法对所有的参数进行优化，主要通过定义增益比ρ，用于调节公式(7)中λ的大小，达到全局优化参数的目的。

g是LM算法梯度，表示为：

g＝J^T(t,p⁽⁰⁾)[f(t)-f(t,p⁽⁰⁾)] (9)

dp是LM算法步长，表示为：

dp＝[H(t,p⁽⁰⁾)+λE]^-1*g (10)

v是常量。

分解的结果如图4a，图4b和图4c所示。

(7)精度评估：采用原始波形数据与计算得到的最优参数多高斯函数模型的均方根误差，平均拟合优度R²，波形分解类型来实现高斯模型精度的验证。

波形分解类型分为明显波形分量，隐藏波形分量和无效波形分量。首先判断出拟合回波信号的所有峰值点和拐点，每两个拐点的位置作为一组，分为N组区间。将距离峰值点位置最近的波形分量判断为明显波形分量。明显波形分量区分出来后，判断剩下的波形分量的中心位置是否在N组拐点位置的区间内，如果在这N组区间内，则该波形分量就是隐藏的波形分量。上述条件都不满足的波形分量分为无效的波形分量。分解结果如图4d所示。

本实例各公式的参数说明见表1。

表1

本发明未详尽描述的技术内容均为公知技术。

Claims (4)

1.激光雷达回波信号高斯拐点选择分解方法，其特征在于包括以下步骤：

步骤S1，采用Savizkg-Golag算法与五点三次平滑算法对回波信号进行预处理，得到预处理后的回波信号数据；

步骤S2，采用高斯拐点选择分解方法对预处理后的回波数据进行波形分解，得到高斯分量的特征参数；

步骤S3，采用传统Levenberg-Marquarelt算法与全局Levenberg-Marquarelt对所述的高斯分量的特征参数进行优化，以获得特征参数的最优解；

步骤S4，将优化后的高斯分量分类，判断为明显波形分量，隐藏波形分量和无效波形分量；包括以下子步骤：判断出拟合回波信号的所有峰值点和拐点，每两个拐点的位置作为一组，分为N组区间，将距离峰值点位置最近的波形分量判断为明显波形分量；判断剩下的波形分量的中心位置是否在N组拐点位置的区间内，如果在这N组区间内，则该波形分量就是隐藏波形分量；上述条件都不满足的波形分量分为无效波形分量。

2.根据权利要求1所述方法，其特征在于，步骤S1具体包括以下子步骤：

步骤S1.1采用Savizkg-Golag算法对回波信号进行初次处理；

步骤S1.2采用五点三次平滑算法对初次处理的到的回波信号数据进行平滑。

3.根据权利要求1所述方法，其特征在于，步骤S2所述包括以下子步骤：

步骤S2.1对回波信号的采样点f(t)采用以下公式寻找回波信号的最大值记为幅值A，并将时间t记为中心位置μ：

f(t)＞f(t+1)＞f(t+2)并且f(t)＞f(t-1)＞f(t-2)

步骤S2.2对回波信号的采样点f(t)采用以下公式寻找回波信号拐点，在幅值A左侧的拐点记为LINF，在幅值A右侧的拐点记为RINF：

(f(t-2)+y(t)-2f(t-1))(f(t-1)+f(t+1)-2f(t))＜0

步骤S2.3根据幅值A左右侧的拐点数采用以下公式判断波形半宽σ：

其中tgl为振幅左侧一半位置的时间点，tgr为振幅右侧一半位置时间点，mid_dis为与中心位置中间的振幅点的距离；dis为与中心位置距离最短的振幅点的距离；

步骤S2.4当找到初始参数A、μ、σ时就可以带入以下公式中，

其中A_i第i个高斯函数幅值，μ_i第i个高斯函数的中心位置，σ_i为第i个高斯函数的标准偏差，Noise为背景噪声，在带入参数计算时可以忽略其影响；

然后使用滤波后的回波信号和高斯函数进行相减，为了得相减后得到振幅值非负，使原始波形信号振幅低于0mv的部分等于0mv；重复S2.1-S2.4,寻找剩下的高斯分量当振幅小于阈值时停止迭代。

4.根据权利要求1所述方法，其特征在于，步骤S3中包括以下子步骤：

步骤S3.1，采用传统Levenberg-Marquardt算法(LM算法)将幅值A_i和中心位置μ_i当成常数，对波形半宽σ_i迭代优化；

p＝p⁽⁰⁾+[H(t,p⁽⁰⁾)+λE]^-1J^T(t,p⁽⁰⁾)[f(t)-f(t,p⁽⁰⁾)]

其中为待定系数p₁,p₂,p₃,...,p_n的函数，/>表示初始参数p⁽⁰⁾＝(p₁,p₂,p₃,...,p_n)，H(t,p)表示f(t,p)的海塞矩阵，J^T(t,p)表示f(t,p)的转置雅克比矩阵，f(t)表示为预处理后的数据，λ为阻尼系数，E为单位矩阵；若||p-p⁽⁰⁾||相差很大，则将p作为新的p⁽⁰⁾，直到||p-p⁽⁰⁾||相差很小或者循环达到所设定的值时，迭代结束；

步骤S3.2，采用全局Levenberg-Marquardt算法对上述的参数进行再次优化，主要通过定义增益比ρ，用于调节公式上述式子中λ的大小，达到全局优化参数的目的；定义增益比ρ通过下面式子定义：

其中f(t,p⁽⁰⁾)为待定系数p₁,p₂,p₃,...,p_n的函数，g是LM算法的梯度，表示为：

g＝J^T(t,p⁽⁰⁾)[f(t)-f(t,p⁽⁰⁾)]，

dp是LM算法的步长，表示为：

dp＝[H(t,p⁽⁰⁾)+λE]^-1*g，

v是常量。