-
机载激光雷达(light detection and ranging，LiDAR)作为一种主动遥感技术，可以直接量测激光扫描仪与地形之间的距离，获取地面反射点的三维坐标，已成为对地空间信息获取的一种重要技术手段[1]。2004年，RIEGL公司首先推出了商用的全波形LiDAR系统，可以以纳秒为单位记录后向散射信号，得到一个连续的高分辨率回波波形，即全波形数据，该波形反映了地物的垂直分布、几何和物理特征[2]。相比传统离散多回波记录方式，全波形记录能够提供更加丰富的地形细节，提高测距精度和地表分辨率，而且用户可以通过分析后向散射波形得到地物更多的物理特征(振幅、宽度等信息)，这在植被提取和地物分类方面有很大的潜力，例如，低矮植被的分离以及冠层高度的精确提取[3]。波形数据分解处理已经成为LiDAR数据处理研究中的热点问题之一。目前对于波形数据的处理方法包括脉冲阈值检测法、相关性技术分析法和信号分解方法[4-8]等。由于激光脉冲回波信号的分布形态近似高斯分布，基于高斯函数的波形信号分解方法逐渐成为LiDAR波形分解的主流方法。Hofton等[9]用两个相邻拐点确定一个高斯分量，提出了基于高斯模型的波形分解方法。Wagner等[6]论述了将波形建模为一系列高斯脉冲的理论基础。Chauve等[10]考虑到波形的细节特征，采用广义高斯函数和对数高斯函数进行波形分解。波形的高斯分解通常包括两个步骤：①确定高斯分量的个数以及其初始参数，②特征参数的优化及后处理。初始参数的估计方法主要包括局部极值、重心法和拐点法等，这些方法通过计算波形数据的一阶导数零点、重力中心、二阶导数零点等[11-13]信息确定初始参数，高斯分量个数可以通过迭代方式获得。波形参数的优化方法通常有两种：①基于概率模型的期望最大化(expectation maximization，EM)及其改进算法[14-15]；②基于非线性最小二乘的相关算法[16-17]，如列文伯格-马夸尔特(Levenber Marquard，LM)优化算法。由于高斯分量的初始值估计容易受噪声影响，参数优化结果依赖比较精确的初始值，且并不能优化高斯分量的个数[12]，因此，针对叠加波等复杂波形，如何更好地估计高斯分量个数及其初始参数仍是目前基于高斯函数分解法需要解决的关键问题。
本文提出一种横向高斯分解方法，该方法在波形数据去噪的基础上，检测波峰类型并估计其初始参数，采用横向迭代分解估计初始高斯分量，利用非线性最小二乘算法优化参数，并解算得到分解点云，利用Leica ALS60记录的波形数据进行波形分解实验和分析，比较了分解点云和系统点云在三维信息表达上的能力。
HTML
-
激光雷达回波波形可以看作激光脉冲和目标表面散射函数的卷积结果，两者都可以用高斯函数来描述，假设接收到的回波信号满足混合高斯分布，则可以用混合高斯模型来描述:
式中，Ai是第i个高斯分量的振幅；μi是第i个高斯分量的波峰位置；σi为第i个高斯分量的标准差；n为高斯函数个数；b为背景噪声。
波形分解流程如图 1所示，由于波形数据以扩展变长记录形式存储在LAS文件末尾，为了便于利用分解参数解算点云，在对波形数据处理前，需要建立波形数据与点云之间的索引。
Figure 1. Waveform Decomposition Technical Procedure
-
波形预处理可以消除原始波形信息中的背景噪声和随机噪声。本文采用塞威兹克-戈利(Savitzky-Golay，S-G)算法去除随机噪声。S-G滤波是一种在时域内基于局域多项式最小二乘法拟合的滤波方法，该滤波器最大的特点在于在滤除噪声的同时可以确保信号的形状、宽度不变[18]。图 2显示了使用S-G滤波和高斯滤波对实测波形数据的滤波效果对比。通过对比发现，S-G滤波在对波形平滑的同时很好地保留了信号的形状、宽度信息，而高斯滤波后波形在波峰处与原始波形相比会出现明显的变形。
Figure 2. Return Waveform Before and After Smoothing Using Gaussian Filter and S-G Filter
-
估计高斯分量的个数及其参数是波形高斯分解的关键步骤，优化过程需要尽可能精确地估计初始参数。对于多次明显叠加波，波峰位置由于叠加的原因会偏离原始位置，因此可利用左右拐点位置确定波峰位置和宽度；对于弱叠加波，左右拐点位置难以判断，选择半波宽公式求高斯分量宽度。对预处理后的波形，由局部最大值确定明显波峰位置。对于检测到的明显波峰，若相邻波峰间距小于4σ0(σ0为系统发射脉冲宽度)，则标记为明显叠加波；对于明显叠加波，在波峰左右一定范围内寻找左右拐点，拐点通过离散数据的二阶五点数值微分求得，并通过式(2)确定该高斯分量的初始参数。
式中，pl、pr分别为波峰μi左右拐点位置；σ为该高斯分量的宽度；μ为波峰位置。对于弱叠加波和独立波，在波峰左右寻找为波峰值一半处的波形数据的位置hl、hr，如果第i-1个波形数据值大于波峰值的一半且第i+1个波形数据值小于波峰值的一半，则第i个波形数据值的位置就作为hl、hr的值。波峰位置为峰值Xi在波形中的位置，宽度由半波宽公式求得。算法通过迭代方式进行波形分解，每次从当前波形中减去波形宽度最大的波形以减少波形叠加对初始参数估计的影响。
-
由于噪声影响，实际的波形分布形态并不完全符合高斯函数的形状，完成初始参数的估计后还需要去掉无意义的高斯分量。实验表明，如果两个高斯分量波峰间距小于一个脉冲宽度时，叠加之后只会有一个主要波峰存在，且脉冲形状与标准的高斯分布相似。因此，对于相邻波形分量波峰位置间距如果小于一个脉冲长度，则把波峰值较小的高斯分量视为无效分量，其判别条件为：
式中，Gi为初始估计的任一高斯分量；DGi, j代表任意两个高斯脉冲之间的波峰距离；Ai为高斯分量的振幅；At为设定的弱波检测阈值；δi为当前高斯分量的宽度；δ0为发射脉冲的宽度。对剩余的回波按弱波检测阈值分为明显波和弱波，对明显波进行拟合，如果拟合后的均方根误差大于预定的精度，则添加一个弱波到拟合过程中，直到满足预定的精度为止。LM算法同时具有梯度法和牛顿法的优点[19]，本文采用LM优化算法对确定初值的高斯分量进行波形拟合，得到的波形分量可用于后续的点云坐标解算以及植被提取等方面的应用。
-
系统点云是LiDAR系统产生的有限次回波点，点位获取方法对用户是不透明的，无法评价点位所在位置的可靠性，点云密度也有待提高。因此，对全波形回波信号进行分析，并解算出更具物理意义的点位信息对于提高点位的可靠性以及辅助后续点云应用具有重要意义。本文基于波形分解得到的高斯分量的波峰位置作为分解点云的响应位置，由波形数据对应的系统点云及相应的转换参数可以得到分解点云，如式(4)：
式中，X、Y、Z为根据波形参数推算点的空间位置；X0、Y0、Z0为起始点的位置；X(t)、Y(t)、Z(t)为激光信号的速度矢量在各坐标方向上的分量；t为相对起始点的时间，单位为皮克秒。
1.1. 波形预处理
1.2. 横向高斯分解
1.3. 参数优化及后处理
1.4. 点云提取
-
本文的实验数据来自Leica ALS 60系统在2012年采集的某地数据，点云和波形信息以LAS1.3标准格式存储在LAS文件中，波形数据采样间隔为1 ns，采样数为256，发射脉冲宽度为5 ns。测区内地物丰富，主要包含建筑物、树木和道路等。
-
图 3显示了原始波形信号进行横向高斯分解及优化的过程，这里背景噪声为回波信号最后5%的数据的平均值。图 4显示了算法处理各种复杂回波的效果，处理的波形包含了弱叠加波，明显叠加波和左右叠加波等类型的回波，图 4(a)、4(b)给出了初始估计的高斯分量，可以看出，对于复杂回波，初始估计的高斯分量中包含有无效的高斯分量。图 4(c)、4(d)显示了优化后的结果，优化过程很好地确定了最佳的高斯分量个数以及参数。相比系统方法和一阶导数零点法，本文方法能够有效地检测出叠加波和弱波等复杂波形，可以分解得到更多的高斯分量以及对应的波形特征信息。
Figure 3. Lateral Gaussian Decomposition Procedure of Waveform and Optimization
Figure 4. Gaussian Decomposition of Complicated Waveform
-
为了对比分解点云与系统点云的差异，本文将波形分解的点云响应位置与系统记录点云在波形中的位置相近的点作为同一地物点，图 5为分解点与系统点在波形上的响应位置对比，可以看出，分解点响应位置更加稳定和可靠。结合分解点响应位置以及公式(4)，可以解算得到分解点坐标，表 1显示了分解点云与系统点云在点位坐标上的差异，由于与分解点与系统点在波形中的响应位置不同，点云坐标在X、Y、Z轴向均存在一定程度的偏差，X、Y偏差程度为厘米级。Z值偏差较大，其原因主要在于，激光脉冲在传播时Z方向被赋予较大的速度矢量。
Figure 5. Comparison of Response Location Between Decomposition Point and System Point
点号 系统点坐标/m 分解点坐标/m 坐标差/m 距离差/m X=483 233.673 X=483 233.642 0.031 1 Y=240 462.434 Y=240 462.451 -0.017 0.225 Z=86.026 Z=85.772 0.254 X=483 234.766 X=483 234.744 0.022 2 Y=240 461.359 Y=240 461.371 -0.012 0.179 Z=85.902 Z=85.723 0.179 X=483 245.717 X=483 245.697 0.020 3 Y=240 450.802 Y=240 450.811 -0.009 0.169 Z=88.426 Z=88.257 0.169 X=483 248.772 X=483 248.734 0.038 4 Y=240 447.607 Y=240 447.626 -0.019 0.334 Z=84.380 Z=84.042 0.338 Table 1. Coordinate Error Between Decomposition Point and System Record Point
表 2给出了分解点云与系统记录点云在回波次数分布情况上的统计对比分析。由表看出，系统最多只能检测4次回波，本文的分解方法可以检测出多达8次回波，分解点数量比系统点增加了61.8%，多次回波分解得到的点云个数明显增加。图 6显示了分解点云与系统点云的差异，图 6(a)为系统点云与分解点云局部细节对比图，可以看出低矮植被可以从地表分离出来，分解点云显示了更丰富的地形细节，图 6(b)为分解点与系统点的差值图，显示出分解的额外点主要集中在植被(包括树木、草地、绿化带)地区，分解的多次回波点明显多于系统点，这主要是由于激光具有很强的植被穿透能力，同时说明了分解点云在林业领域的应用前景。
回波次数 分解点数 系统点数 差异 1 152 299 356 507 -204 208 2 381 256 34 170 347 086 3 74 472 1 275 73 197 4 20 924 20 20 904 5 4 485 0 4 485 6 606 0 606 7 77 0 77 8 8 0 8 总和 634 127 391 927 242 200 Table 2. Comparison of Points Derived from Decomposition and the Hardware System
Figure 6. Comparison of Decomposition Point and System Point
2.1. 波形分解实验数据
2.2. 波形分解及优化
2.3. 点云生成及分析
-
针对传统LiDAR波形数据分解方法易受噪声影响，高斯组分个数及叠加波和弱波参数估计不精确等问题，本文提出并实现了一种横向高斯分解方法，可用于小光斑LiDAR波形数据精细分解。实验表明，该方法能够处理各种类型的回波信号，具有良好的鲁棒性和适应性，分解的点云不仅在数量上具有优势，在点位可靠性方面也有所提高。本文主要采用S-G滤波算法既对波形进行了平滑，保持了波形的形状、宽度不变；对于不同类型的回波分别采用合适的方法进行参数估计；通过剔除无效高斯分量，确定了最佳高斯分量个数；从点云数量、信号位置以及坐标精度等多方面定量比较了波形分解点云与系统点云的差异。由于缺少地面实测数据，本文只对比分析了分解点云与系统点云的相对位置精度。进一步的研究包括点位绝对精度验证以及波形特征参数在地物分割、分类方面的应用。
DownLoad:
