CN113866736A - 激光雷达回波信号高斯拐点选择分解方法 - Google Patents
激光雷达回波信号高斯拐点选择分解方法 Download PDFInfo
- Publication number
- CN113866736A CN113866736A CN202110956522.9A CN202110956522A CN113866736A CN 113866736 A CN113866736 A CN 113866736A CN 202110956522 A CN202110956522 A CN 202110956522A CN 113866736 A CN113866736 A CN 113866736A
- Authority
- CN
- China
- Prior art keywords
- gaussian
- echo signal
- laser radar
- waveform
- radar echo
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000354 decomposition reaction Methods 0.000 title claims abstract description 49
- 238000000034 method Methods 0.000 title claims abstract description 37
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 21
- 238000007781 pre-processing Methods 0.000 claims abstract description 6
- 210000003127 knee Anatomy 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000009499 grossing Methods 0.000 claims description 4
- 238000012795 verification Methods 0.000 claims description 3
- 238000001514 detection method Methods 0.000 abstract 1
- 238000005070 sampling Methods 0.000 description 13
- 238000012545 processing Methods 0.000 description 5
- 238000002592 echocardiography Methods 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 238000003672 processing method Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 102100036167 CXXC-type zinc finger protein 5 Human genes 0.000 description 1
- 101000947154 Homo sapiens CXXC-type zinc finger protein 5 Proteins 0.000 description 1
- 238000012897 Levenberg–Marquardt algorithm Methods 0.000 description 1
- 238000013016 damping Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
- G01S7/418—Theoretical aspects
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Radar Systems Or Details Thereof (AREA)
- Optical Radar Systems And Details Thereof (AREA)
Abstract
本发明公开了激光雷达回波信号高斯拐点选择分解方法，包括：(1)收集相应的激光雷达回波信号数据；(2)对激光雷达回波信号数据进行预处理；(3)构建高斯模型并求解；(4)高斯模型精度验证。本发明操作简便且效率高，能够有效的分解复杂激光雷达回波信号，特别是只存在一个峰值的复杂激光雷达回波信号，为激光雷达回波信号分解提供了便捷的分解算法，在激光雷达遥感探测中具有重要意义。
Description
技术领域
本发明涉及激光雷达领域，特别涉及到激光雷达数据处理，提出了激光雷达回波信号高斯选择拐点分解方法。
背景技术
激光雷达已经广泛应用于测绘、地质、城市规划、农业开发、水利水电、土木工程、环境检测等领域。为了获得高精度的激光雷达点云数据，有效的对全波形激光雷达回波信号进行波形分解，成为全波形激光雷达数据处理的一个热点。
全波形激光雷达回波信号过程是：激光从激光器中发射而出，在激光照射到地物后，接收器接收到地物的反射光并将数据进行存储，再将存储的数据转换为数字的形式，所得到的数据就是全波形激光雷达回波信号数据。波形分解需要计算全波形激光雷达回波信号在时间轴上的相应幅值位置，从而得到精确的时间间隔，通过光速与时间间隔相乘就可以得到激光器到目标物的距离。因此，计算准确的时间间隔是全波形激光雷达回波信号数据处理工作的重中之重。
近二十年来，我国许多科研人员从事对激光雷达回波信号分解的研究，主要方法有反卷积法和数学模拟法，由于反卷积法容易受到噪声影响，并且其分解之后的信号变窄，无法通过波形的形状、宽度、强度和偏度提供了更多的物体信息，所以常使用数学模拟法应用于全波形激光雷达回波信号的数据处理中。在数学模拟法中最常用的方法就是高斯分解法。高斯分解法的原理是通过多个高斯函数进行建模，模拟全波形激光雷达回波信号。高斯分解法需要估计回波信号中高斯分量的个数和高斯分量的参数，再对拟合回波信号的高斯估计参数进行优化。武汉大学马洪超等使用拐点的方法计算高斯参数，拐点法计算全波形激光雷达回波信号的拐点，再通过拐点和峰值点计算出高斯函数需要的三个参数，从而完成参数估计。拐点法虽然能够寻找到高斯参数，但是并不适用于分解复杂的回波信号，在参数拟合中容易使得高斯波形产生负向震荡的波形，导致距离提取时数值失真。南京大学徐帆等提出峰值拐点的方法计算高斯参数。峰值拐点法是寻找所有波形的峰值，并计算所有拐点位置，通过拐点和峰值判断高斯波形的个数。但是如果有多个回波信号叠加的复杂回波信号，这些回波信号非常接近并且只产生一个峰值点时，峰值拐点法就会分解出一个高斯波形。
针对上述全波形激光雷达波形分解方法遇到的问题，本发明公开了高斯拐点选择分解方法，该方法能够解决分解复杂回波信号的问题，并且能够处理只有一个峰值点并存在多个拐点的复杂叠加波形。
发明内容
为了达到上述的目的，本发明提供了如下技术方案：
激光雷达高斯拐点选择分解方法，包括如下步骤：
(1)收集相应的激光雷达回波信号数据；
(2)对激光雷达回波信号数据进行预处理；
(3)构建高斯模型并求解；
(4)高斯模型精度验证。
优选地，步骤(1)中提供全波形激光雷达回波数据。
优选地，步骤(2)所述回波信号预处理包括Savitzky-Golay算法进行信号去噪与五点三次平滑法处理回波信号。
优选地，步骤(3)所述的高斯模型为：
Ai，第i个高斯函数幅值；
μi，第i个高斯函数的中心位置；
σi，第i个高斯函数的标准偏差；
Noise，背景噪声；
优选地，对多高斯函数模型进行求解时采用高斯拐点选择分解方法，根据峰值和拐点数目的计算出高斯参数，将各参数带入高斯函数模型，得到多组高斯分量。接着将不符合要求的高斯分量进行剔除。然后利用传统的Levenberg-Marquarelt算法对高斯分量的标准偏差σi进行优化。再利用全局Levenberg-Marquarelt算法，优化所有的高斯分量的参数。最后将优化后高斯分量区分为明显的高斯分量，隐藏的高斯分量和无效的高斯分量。
优选地，步骤(4)中采用原始波形数据与计算得到的最优参数高斯函数模型的均方根误差，平均拟合优度R2，波形分解数，波形分解类型来实现高斯模型精度的验证。
本发明采用上述的技术方案，具体有以下的增益效果：
本发明利用获取到的全波形激光雷达回波信号数据，能够根据不同拐点数判断出全波形激光雷达回波信号的幅值，中心位置和标准偏差，能够为激光雷达回波信号数据处理提供相应的较为准确的分解方法，在全波形激光雷达回波信号分解应用中具有重要的意义。
附图说明
图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右侧的两个拐点对应的时间位置。计算左侧两个拐点的时间位置的平均值a1，找到最接近a1的采样点的时间位置记为S1。计算振幅A右侧两个拐点的时间位置的平均值a2，找到最接近a2的采样点的时间位置记为S2。寻找幅值A两侧最接近幅值A一半的采样点，记录其对应的时间位置记为S3。然后将S1,S2和S3与中心位置μ进行相减并取绝对值，去掉最大值和最小值后所得到的数值记为|μ-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算法)将幅值Ai和中心位置μi当成常数，迭代优化σi。
p＝p(0)+[H(t,p(0))+λE]-1JT(t,p(0))[f(t)-f(t,p(0))] (7)
其中f(t,p)为待定系数p1,p2,p3,...,pn的函数，p(0)表示初始参数p(0)＝(p1,p2,p3,...,pn)，H(t,p)表示f(t,p)的海塞矩阵，JT(t,p)表示f(t,p)的转置雅克比矩阵，f(t)表示为预处理后的数据，λ为阻尼系数。若||p-p(0)||相差很大，则将p作为新的p(0)，直到||p-p(0)||相差很小或者循环达到所设定的值时，迭代结束。
(6)采用全局LM算法对所有的参数进行优化，主要通过定义增益比ρ，用于调节公式(7)中λ的大小，达到全局优化参数的目的。
g是LM算法梯度，表示为：
g＝JT(t,p(0))[f(t)-f(t,p(0))] (9)
dp是LM算法步长，表示为：
dp＝[H(t,p(0))+λE]-1*g (10)
v是常量。
分解的结果如图4a，图4b和图4c所示。
(7)精度评估：采用原始波形数据与计算得到的最优参数多高斯函数模型的均方根误差，平均拟合优度R2，波形分解类型来实现高斯模型精度的验证。
波形分解类型分为明显波形分量，隐藏波形分量和无效波形分量。首先判断出拟合回波信号的所有峰值点和拐点，每两个拐点的位置作为一组，分为N组区间。将距离峰值点位置最近的波形分量判断为明显波形分量。明显波形分量区分出来后，判断剩下的波形分量的中心位置是否在N组拐点位置的区间内，如果在这N组区间内，则该波形分量就是隐藏的波形分量。上述条件都不满足的波形分量分为无效的波形分量。分解结果如图4d所示。
本实例各公式的参数说明见表1。
表1
本发明未详尽描述的技术内容均为公知技术。
Claims (5)
1.高斯拐点选择分解方法，其特征在于，包括如下步骤：
(1)收集相应的激光雷达回波信号数据；
(2)对激光雷达回波信号数据进行预处理；
(3)构建高斯模型并求解；
(4)高斯模型精度验证。
2.根据权利要求1所述高斯拐点选择分解方法，其特征在于，步骤(1)中提供全波形激光雷达回波数据。
3.根据权利要求1所述高斯拐点选择分解方法，其特征在于，步骤(2)所述回波信号预处理包括Savitzky-Golay算法进行信号去噪与五点三次平滑法处理回波信号。
4.根据权利要求1所述高斯拐点选择分解方法，其特征在于，步骤(3)所述的高斯模型为：
Ai，第i个高斯函数幅值；
μi，第i个高斯函数的中心位置；
σi，第i个高斯函数的标准偏差；
Noise，背景噪声；
采用高斯拐点选择分解方法，根据峰值和拐点数目计算出高斯参数，将各参数带入高斯函数模型，得到多组高斯分量；接着将不符合要求的高斯分量进行剔除；然后利用传统的Levenberg-Marquarelt算法对高斯分量的标准偏差σi进行优化；再利用全局Levenberg-Marquarelt算法，优化所有的高斯分量的参数；最后将优化后高斯分量区分为明显的高斯分量，隐藏的高斯分量和无效的高斯分量。
5.根据权利要求1所述高斯拐点选择分解方法，其特征在于步骤(4)中采用原始波形数据与计算得到的最优参数高斯函数模型的均方根误差，平均拟合优度R2，波形分解数，波形分解类型来实现高斯模型精度的验证。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110956522.9A CN113866736B (zh) | 2021-08-19 | 2021-08-19 | 激光雷达回波信号高斯拐点选择分解方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110956522.9A CN113866736B (zh) | 2021-08-19 | 2021-08-19 | 激光雷达回波信号高斯拐点选择分解方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113866736A true CN113866736A (zh) | 2021-12-31 |
| CN113866736B CN113866736B (zh) | 2023-12-19 |
Family
ID=78990734
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110956522.9A Active CN113866736B (zh) | 2021-08-19 | 2021-08-19 | 激光雷达回波信号高斯拐点选择分解方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113866736B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2023134396A1 (zh) * | 2022-01-14 | 2023-07-20 | 华为技术有限公司 | 设备感知方法及装置 |
| CN117665744A (zh) * | 2023-09-25 | 2024-03-08 | 中国人民解放军国防科技大学 | 对高重叠高光谱激光雷达回波波形进行处理的方法和系统 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110673109A (zh) * | 2019-11-01 | 2020-01-10 | 自然资源部国土卫星遥感应用中心 | 一种星载大光斑激光雷达全波形数据分解方法 |
-
2021
- 2021-08-19 CN CN202110956522.9A patent/CN113866736B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110673109A (zh) * | 2019-11-01 | 2020-01-10 | 自然资源部国土卫星遥感应用中心 | 一种星载大光斑激光雷达全波形数据分解方法 |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2023134396A1 (zh) * | 2022-01-14 | 2023-07-20 | 华为技术有限公司 | 设备感知方法及装置 |
| CN117665744A (zh) * | 2023-09-25 | 2024-03-08 | 中国人民解放军国防科技大学 | 对高重叠高光谱激光雷达回波波形进行处理的方法和系统 |
| CN117665744B (zh) * | 2023-09-25 | 2024-06-28 | 中国人民解放军国防科技大学 | 对高重叠高光谱激光雷达回波波形进行处理的方法和系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113866736B (zh) | 2023-12-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104614718B (zh) | 基于粒子群算法的激光雷达波形数据分解的方法 | |
| CN113866736A (zh) | 激光雷达回波信号高斯拐点选择分解方法 | |
| CN106500671B (zh) | 一种基于lm算法分解激光雷达波形确定海水深度的方法 | |
| CN110134976B (zh) | 一种机载激光测深信号提取方法及系统 | |
| CN105676205A (zh) | 一种机载LiDAR波形数据高斯分解方法 | |
| CN103217679A (zh) | 一种基于遗传算法的全波形激光雷达回波数据的高斯分解法 | |
| CN105137498A (zh) | 一种基于特征融合的地下目标探测识别系统及方法 | |
| Schwarz et al. | Exponential decomposition with implicit deconvolution of lidar backscatter from the water column | |
| CN102521835A (zh) | 一种基于空间三维模板的点云数据树高提取方法 | |
| CN106199537B (zh) | 基于逆高斯纹理海杂波幅度分布参数的分位点估计方法 | |
| CN102722640B (zh) | 一种顾及邻近波形信息的机载激光波形数据分解算法 | |
| CN113608195B (zh) | 一种激光雷达全波形数据分解测深方法、装置及电子设备 | |
| Zhang et al. | Denoising for satellite laser altimetry full-waveform data based on EMD-Hurst analysis | |
| Wang et al. | Online multi-target laser ranging using waveform decomposition on FPGA | |
| CN114779355A (zh) | 基于发射电流全波形的地面瞬变电磁法反演方法及装置 | |
| CN110133680B (zh) | 一种机载激光测深接收波形有效信号初值确定方法及系统 | |
| Berntsen et al. | Internal pressure gradient errors in σ-coordinate ocean models in high resolution fjord studies | |
| CN114488168B (zh) | 基于最大正向偏差的卫星激光测距全波形高斯拟合方法 | |
| Mader et al. | Volumetric nonlinear ortho full-waveform stacking in airborne LiDAR bathymetry for reliable water bottom point detection in shallow waters | |
| CN109655910A (zh) | 基于相位校正的探地雷达双参数全波形反演方法 | |
| CN115128680A (zh) | 一种磁性源多波形组合的瞬变电磁靶向测量方法 | |
| Richter et al. | Determination of 3D Water turbidity parameter fields from LiDAR bathymetry data by volumetric data analysis | |
| Chen et al. | The improved kriging interpolation algorithm for local underwater terrain based on fractal compensation | |
| CN113589364A (zh) | 基于佐布里兹方程约束的地震数据规则化处理方法 | |
| CN117310696B (zh) | 一种探地雷达自聚焦后向投影成像方法及装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |