一种顾及邻近波形信息的机载激光波形数据分解算法

摘要

本发明公开了一种顾及邻近波形信息的机载激光波形数据分解方法，包括A、波形数据的预处理；B、波形数据的脉冲波峰探测；C、基于广义高斯模型的显著脉冲波峰提取；D、顾及邻近波形的候选脉冲波峰提取；E、基于区域生长策略的待分析波形列表生成；F、波形分解结果输出。本发明采用高、低双阈值，可以克服单阈值难以合适选取的缺陷，有效考虑波形数据中的微弱脉冲信号；波形分解过程中，通过区域生长的方式，不断利用邻近波形显著脉冲波峰所提供的可靠信息，对待分析波形的候选脉冲波峰进行分析，可以实现这类微弱脉冲波峰信息的稳健提取，获取更加完整的点云信息。

说明书

技术领域

本发明涉及测绘科学与技术领域，涉及一种顾及邻近波形信息的机载激光波形 数据分解算法，该方法能对波形数据中微弱的脉冲信号进行有效检测，以获取更加 完整的点云，非常适用于植被覆盖严重的地区。

背景技术

机载LiDAR(Light Detection And Ranging)，又称机载激光扫描技术，作为一种 三维空间信息的实时获取手段，可以快速、精确的获取大范围区域的地表信息。根 据回波记录方式的不同，机载激光雷达系统可以分为离散激光雷达系统和全波形激 光雷达系统。前者记录有限个离散的回波信号，而后者以很小的采样间隔对来自目 标的激光反射信号进行采样记录，形成一个随时间变化的回波信号。相对于离散激 光雷达而言，全波形雷达可以提供更多的细节信息，具有更强的地物区分能力和更 优的植被穿透特性等。同时，也给数据处理和信息提取提出了更高的要求。

如何从波形采样数据中提取出高质量的点云是数据处理中的一个关键问题。 Hofton等提出机载激光雷达波形数据高斯分解算法，Wagner等从理论角度阐述了机 载全波形激光雷达的成像机理，并使用高斯模型提取波形数据峰值点。这些方法大 多数情况下可以获取较好结果，但是由于各种干扰因素，回波信息并不是高斯函数 的精确表达。Chauve等提出利用广义高斯模型来拟合波形数据，并通过迭代的方式 实现波形数据的分解。马洪超等采用Expectation-Maximization算法(EM)进行机载雷 达波形数据高斯混合模型参数的最大似然估计，取得较好结果但计算量较大。通过 分析发现，目前算法中主要存在如下几点缺陷：

第一，噪声阈值难以确定。目前，为排除信号噪声的影响，大多算法采用经验 取值的方法，并指定为较高的数值。

第二，微弱信号难以可靠检测。当利用较高阈值对波形数据进行预处理后，一 些微弱脉冲信号由于信号强度不够而无法得到考虑，从而不能进行检测。这种现象 在植被茂密地区尤为严重，由于植被遮挡造成地面点反射脉冲信号比较微弱，给地 面信息的提取造成较大挑战。

发明内容

本发明的目的是在于针对现有技术存在的上述问题，提供一种顾及邻近波形信 息的机载激光波形数据分解方法，利用高、低双阈值数据预处理，同时通过顾及邻 近波形显著脉冲波峰信号所提供的信息，对比较微弱的脉冲信息进行分析，融入区 域生长策略，最终实现波形数据中脉冲波峰信号的完整提取，可以得到更加精确、 完整的点云数据。

本发明的上述目的通过以下技术方案实现：

一种顾及邻近波形信息的机载激光波形数据分解方法，包括以下步骤：

步骤1、波形数据的预处理：利用一维高斯模板对原始波形数据进行平滑处理； 利用一阶微分算子计算经过平滑处理后的波形数据的一阶导数；计算每一个波形数 据对应的近似物方平面坐标，并根据近似物方平面坐标建立k-d树索引；

步骤2、波形数据的脉冲波峰探测：通过计算步骤1中得到的一阶导数零交叉 点，确定波形的脉冲波峰中心的初始位置，采用高、低双阈值对波峰强度进行分析， 确定出显著脉冲波峰和候选脉冲波峰，并剔除伪脉冲波峰；

步骤3、基于广义高斯模型的显著脉冲波峰提取：对当前所确定的显著脉冲波 峰，采用广义高斯模型和非线性最小二乘Levenburg-Marquardt算法进行波形数据的 拟合，提取显著脉冲波峰的参数，并统计波形数据拟合的单位权中误差；

步骤4、顾及邻近波形的候选脉冲波峰提取：以含有候选脉冲波峰的波形为中 心，通过k-d树索引搜索设定半径范围内的邻近波形，并计算邻近波形的显著脉冲 波峰提供的高程值与设置为中心的波形的候选脉冲波峰的高程值的高程差，若高程 差小于给定阈值，则认为该候选脉冲波峰与邻域环境相容，将该候选脉冲波峰重新 认定为显著脉冲波峰，采用广义高斯模型和非线性最小二乘Levenburg-Marquardt 算法进行重新认定的显著脉冲波峰的参数的提取，并将重新认定为显著脉冲波峰的 对应的波形加入到已更新波形列表中；

步骤5、基于区域生长策略的待分析波形列表生成：以已更新波形列表中的波 形为中心，通过k-d树搜索设定半径范围内的邻近波形，如果邻近波形中含有候选 脉冲波峰，将该邻近波形加入到待分析波形列表中，当待分析波形列表不为空时， 返回步骤4；当待分析波形列表为空时，则进入步骤6；

步骤6、波形分解结果输出：利用所获取的显著脉冲波峰的参数，计算其对应 激光点的三维空间坐标，输出点云。

如上所述的步骤1中一维高斯模板，其连续情况下的数字形式如下：

$f\left(x\right)=\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{x^2}{2{\sigma }^2})$

其中，σ为高斯模板的标准偏，π为圆周率。

如上所述的σ是基于以下公式：

$\sigma =\frac{1}{2\sqrt{2\ln 2}}\mathrm{FWHM}$

其中，FWHM为激光器波形的半宽参数。

如上所述的步骤1中建立k-d树索引包括以下步骤：

步骤4.1、计算波形对应的近似物方平面坐标，第k个波形近似物方平面坐标 的计算，采用如下公式：

$E_k^{\mathrm{appr}}=E_0+\mathrm{dE}\times \mathrm{WFOFFSET}$

$E_k^{\mathrm{appr}}=N_0+\mathrm{dN}\times \mathrm{WFOFFSET}$

其中，E₀,N₀为第k个波形第一个采样点所对应的平面坐标，dE,dN为第k个波形的 平面坐标微分，WFOFFSET为第k个波形的回波数据第一个采样点到发射波形第一 次采样点的偏移量；

步骤4.2、利用波形对应的近似物方平面坐标建立k-d树索引。

如上所述的步骤2中采用高、低双阈值对波峰强度进行分析，确定出显著脉冲 波峰和候选脉冲波峰，并剔除伪脉冲波峰包括以下步骤：

设定高阈值为HighThred，低阈值为LowThred，脉冲波峰位置i处信号强度为 WVal[i]；

当WVal[i]＞HighThred时，该脉冲波峰为显著脉冲波峰；

当LowThred≤WVal[i]≤HighThred时，该脉冲波峰为候选脉冲波峰，将候选脉冲波 峰所在的波形作为待分析波形加入待分析波形列表中；

当WVal[i]＜LowThred时，该脉冲波峰为伪脉冲波峰，剔除伪脉冲波峰；

若波形中候选脉冲波峰个数不为零，则将该波形加入待分析波形列表中。

如上所述的步骤4中高程值是基于以下公式：

H＝H₀+dH×(WFOFFSET+t)

其中，H为高程值，H₀为显著脉冲波峰/候选脉冲波峰所在波形第一个采样点 所对应的高程坐标，dH为波形数据单元的高程坐标微分，WFOFFSET为回波数据 第一个采样点到发射波形第一次采样点的偏移量，t为显著脉冲波峰/候选脉冲波峰 对应的位置。

如上所述的步骤6中点云的计算基于以下公式：

E＝E₀+dE×(WFOFFSET+t)

N＝N₀+dN×(WFOFFSET+t)

H＝H₀+dH×(WFOFFSET+t)

其中，E、N、H为波形点云三维坐标，E₀、N₀、H₀为波形第一次采样点所 对应的平面坐标，dE、dN、dH为波形的平面坐标微分，WFOFFSET为回波数据 第一次采样点到发射波形第一次采样点的偏差量，t为显著脉冲波峰对应的位置。

本发明与现有技术相比，具有以下优点：

1、采用高、低双阈值，可以克服单阈值难以合适选取的缺陷，有效考虑波形数 据中的微弱脉冲信号；

2、波形分解过程中，通过区域生长的方式，不断利用邻近波形显著脉冲波峰所 提供的可靠信息，对待分析波形的候选脉冲波峰进行分析，可以实现这类微弱脉冲 波峰信息的稳健提取，获取更加完整的点云信息。

3、利用该算法提取的点云相比于相关商业软件生成的点云，提取的有效地面点 云信息可达到30%以上，成效显著。

附图说明

图1是本发明的流程示意图。

具体实施方式

以下结合附图对本发明的技术方案作进一步详细描述。

实施例1

如图1，一种顾及邻近波形信息的机载激光波形数据分解方法，包括以下步骤：

步骤1、波形数据的预处理：利用一维高斯模板对原始波形数据进行平滑处理； 利用一阶微分算子计算经过平滑处理后的波形数据的一阶导数；计算每一个波形数 据对应的近似物方平面坐标，并根据近似物方平面坐标建立k-d树索引；

基于高斯模板的波形数据平滑处理：利用数字图像处理中常用的一维高斯模板， 其连续情况下的数字形式如下：

$f\left(x\right)=\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{x^2}{2{\sigma }^2})$

其中，σ为高斯模板的标准偏，π为圆周率。σ与信号滤波平滑的程度密切相关，σ 越大，高斯滤波器的频带就越宽，平滑的程度就越强；反之，平滑的程度就越弱。 为了在平滑噪声的同时尽量保留信号特征，本实施例中σ的选取根据激光器波形的 半宽参数(full width at half maximum,FWHM)确定，如果激光器波形的半宽参数厂商 没有提供，可以利用平地处所反射激光波形的半宽参数来代替。当得到半宽参数 FWHM后，依据公式(2)可以计算得到高斯模板的标准偏差σ：

$\sigma =\frac{1}{2\sqrt{2\ln 2}}\mathrm{FWHM}\approx 0.42466\mathrm{FWHM}---\left(2\right)$

高斯滤波模板窗口的半径取为3σ，按照一定间隔，对连续的一维高斯函数进行离散 化和归一化处理，即可获取所需的高斯模板。之后，利用高斯模板对原始波形数据 进行卷积运算，即可实现波形数据中噪声的有效抑制。

计算波形数据的一阶导数值：假设高斯滤波后的波形数据为WVal[i](i＝0,…,m)， 其中m为波形数据的采样数，对于位置i处波形数据对应的一阶导数Grad[i]，采用 如下公式计算：

Grad[i]＝WVal[i+1]-WVal[i]    (3)

建立波形数据的k-d树索引：首先计算波形对应的近似物方平面坐标，第k个 波形近似物方平面坐标的计算，采用如下公式：

$E_k^{\mathrm{appr}}=E_0+\mathrm{dE}\times \mathrm{WFOFFSET}$(4)

$E_k^{\mathrm{appr}}=N_0+\mathrm{dN}\times \mathrm{WFOFFSET}$

式中，(E₀,N₀)为第k个波形第一个采样点所对应的平面坐标，(dE,dN)为第k个波 形的平面坐标微分，WFOFFSET为第k个波形的回波数据第一个采样点到发射波形 第一次采样点的偏移量，这些参数均可以从后缀名为lgc的文件中读取。

然后利用波形对应的平面坐标建立k-d树索引，k-d树是由Bentley于1975年提 出并逐渐发展成为一种比较好的多维空间索引，特别适合于空间点状目标。它通过 超平面将一个空间递归划分为两个子空间来实现搜索，其实现算法比较稳定和成熟， 本专利中采用了马里兰大学David M.Mount和香港理工大学Sunil Arya提供的最近 邻算法软件库(Approximate Nearest Neighbor，ANN)实现，该软件包提供了最临近查 询和范围查询两种方式，其中最临近查询方式为查询给定点周围若干最临近点，范 围查询方式为查询给定点给定半径大小范围内的所有点。至此，可以利用k-d树索 引和波形对应的近似平面坐标完成邻域波形的查询。

步骤2、波形数据的脉冲波峰探测：通过计算步骤1中得到的一阶导数零交叉 点，确定波形的脉冲波峰中心的初始位置，采用高、低双阈值对波峰强度进行分析， 确定出显著脉冲波峰和候选脉冲波峰，并剔除伪脉冲波峰；主要包括以下步骤：

设定高阈值为HighThred，低阈值为LowThred，脉冲波峰位置i处信号强度为 WVal[i]；

当WVal[i]＞HighThred时，该脉冲波峰为显著脉冲波峰；

当LowThred≤WVal[i]≤HighThred时，该脉冲波峰为候选脉冲波峰，将候选脉冲波 峰所在的波形作为待分析波形加入待分析波形列表中；

当WVal[i]＜LowThred时，该脉冲波峰为伪脉冲波峰，剔除伪脉冲波峰；

若波形中候选脉冲波峰个数不为零，则将该波形加入待分析波形列表中。遍历 每一个波形数据，并进行上述处理。

步骤3、基于广义高斯模型的显著脉冲波峰提取：对当前所确定的显著脉冲波 峰，采用广义高斯模型和非线性最小二乘Levenburg-Marquardt算法进行波形数据的 拟合，提取显著脉冲波峰的参数，并统计波形数据拟合的单位权中误差；

当波形的显著脉冲波峰个数不为零时，利用当前确定的显著脉冲波峰的近似参 数，通过广义高斯模型和非线性最小二乘算法Levenberg-Marquardt对波形进行分解。

波形分解的目的就是要将波形分解成多个脉冲信号的集合，通过记录各个脉冲 信号的特征来描述不同地物的位置和反射强度等信息。本实施例中将每个脉冲信号 视为广义高斯函数，则波形的形式如下：

$y=f\left(x\right)=\underset{i=1}{\overset{n}{\Sigma }}{f}_i\left(x\right)---\left(5\right)$

式(5)中，y为波形值，n代表脉冲信号个数，取值为当前所确定的显著脉冲波峰个 数，f_i(x)为广义高斯函数，其具体表达形式如下：

$f_i\left(x\right)=P_i\exp (-\frac{{(x-t_i)}^{r_i^2}}{2s_i^2})---\left(6\right)$

式(6)中，P_i为脉冲波形的振幅，t_i为脉冲波形对应波峰的位置，s_i为脉冲波形的标 准方差，r_i为形状参数。容易发现，当时，相当于标准的高斯函数；广义高 斯模型通过在标准高斯函数的基础上增加了形状参数r_i，达到改变波形形状的目的， 从而可以更好的拟合波形的形状。

在建立波形分解的数学模型后，波形分解的过程其实就是各个脉冲波形参数的 优化计算过程，具体采用Levenberg-Marquardt算法实现。该算法是数学中最常见的 优化算法，是一种介于牛顿法与梯度下降法之间的一种非线性优化算法，其优点在 于对过参数化问题不敏感，能够有效处理冗余参数带来的强相关问题，是迭代优化 陷于局部极小值的机会大大减少。

利用Levenberg-Marquardt算法实现波形数据分解的步骤如下：

1、参数初始值的确定：对于广义高斯函数参数，在首次迭代中，t_i为利用波形 一阶导数的零交叉点所确定的波峰位置，P_i为波峰位置处信号强度值WVal[t_i]，缺省 认为脉冲为标准高斯函数形状，故r_i取为此外，脉冲信号半宽参数FWHM一 般为3，故s_i取为1.5；在随后的迭代中，参数初始值取为前一次迭代平差所计算的 新参数值。对于Levenberg-Marquardt算法中所使用的变量，k为迭代次数，λ_k为阻 尼系数，v为阻尼系数调整常数，本实施例中k初始化为0，λ_k初始值λ₀设为0.01， v取值为10。

2、计算拟合函数的雅可比矩阵J_k：

式（7）中，n为当前探测出的显著波峰个数，m为所使用的波形数据点个数， p_k＝[P₁ t₁ s₁ r₁…P_n t_n s_n r_n]^T为待解算的未知数，为各脉冲信号的波形 参数。

然后，计算海塞矩阵H:

$H=J_k^T{J}_k+{\lambda }_{k}I---\left(8\right)$

并构造增量正规方程

$H·{\delta }_k=J_k^T{ϵ}_k---\left(9\right)$

式中，δ_k表示待解算的未知数p_k的增量向量，ε_k表示利用未知数当前值所计 算的波形值与其对应观测值的误差向量。利用公式(9)进行求解，可以得到δ_k。

3、阻尼系数调整：

3.1、利用δ_k更新参数向量并计算新的误差向量，如果‖y-f(p_k+δ_k)‖＜||ε_k||(||·||表 示向量的L2-范数)，即误差向量变小，则令p_k+1＝p_k+δ_k；如果||ε_k‖小于为给定收敛 阈值，满足收敛条件，则停止迭代，输出结果，否则令λ_k+1＝λ_k/v，转到步骤2中 进行迭代。

3.2、如果‖y-f(p_k+δ_k)‖＞||ε_k‖，则令λ_k+1＝λ_k·v，返回到步骤1重新解算。

计算拟合中误差σ：当迭代满足收敛条件后，利用最终得到的误差向量ε_k计算 拟合中误差

步骤4、顾及邻近波形的候选脉冲波峰提取：以含有候选脉冲波峰的波形为中 心，通过k-d树索引搜索设定半径范围内的邻近波形，并计算邻近波形的显著脉冲 波峰提供的高程值与设置为中心的波形的候选脉冲波峰的高程值的高程差，若高程 差小于给定阈值，则认为该候选脉冲波峰与邻域环境相容，将该候选脉冲波峰重新 认定为显著脉冲波峰，采用广义高斯模型和非线性最小二乘Levenburg-Marquardt 算法进行重新认定的显著脉冲波峰的参数的提取，并将重新认定为显著脉冲波峰的 对应的波形加入到已更新波形列表中；

首先，利用步骤1中波形数据预处理中所建立的k-d树索引，和待分析波形的 近似平面位置，查询待分析波形一定半径大小范围内的所有波形数据，半径取值为 点云平均距离的3倍。然后，利用如下公式计算邻近波形中显著脉冲波峰所对应的 地面高程H：

H＝H₀+dH×(WFOFFSET+t)    (10)

式中，H₀为显著脉冲波峰所在波形第一个采样点所对应的高程坐标，dH为波形数 据单元的高程坐标微分，WFOFFSET为回波数据第一个采样点到发射波形第一次采 样点的偏移量，t为显著脉冲波峰对应的位置，这些参数均可以从后缀名为lgc文件 获取得到。

同样，采用公式(10)，计算待分解波形中候选脉冲波峰对应的高程值，并将其 与邻近波形中显著脉冲波峰对应的高程值求差，如果差值小于给定高差阈值，则认 为该候选脉冲与周围环境是相容的，并将其确定为显著脉冲波峰。

如果有候选脉冲波峰被认定为显著脉冲波峰，则重新对该波形进行基于广义高 斯模型的波形分解，精确提取当前显著脉冲波峰的参数信息。波形分解成功后，将 该波形加入到已更新波形列表中。

步骤5、基于区域生长策略的待分析波形列表生成：以已更新波形列表中的波 形为中心，通过k-d树搜索设定半径范围内的邻近波形，如果邻近波形中含有候选 脉冲波峰，将该邻近波形加入到待分析波形列表中，当待分析波形列表不为空时， 返回步骤4；当待分析波形列表为空时，则进入步骤6；

对于已更新波形列表，由于增加了新的显著脉冲波峰，将对其所在邻域中含有 候选脉冲波峰的波形数据造成影响。为此，对已更新波形列表中的波形，逐一通过 前面建立的k-d树索引，搜索给定半径大小范围内的邻近波形，半径取值为点云平 均距离的3倍。如果邻近波形中存在候选脉冲波峰，将该波形加入到待分析波形列 表中。

当得到待分析波形列表后，采用步骤4中的顾及邻近波形的候选脉冲波峰提取 方法，逐一对待分析波形的候选脉冲波峰进行分析和探测。波形分解完成后，重新 生成已更新波形列表。同样，通过建立的k-d树索引，根据新产生的已更新波形列 表，生成新的待分析波形列表。可以发现，这种方式非常类似于区域生长策略，仅 不断考虑有改变的波形的邻域，这样不仅减少了计算量，而且实现候选脉冲波峰的 递进提取。

步骤6、波形分解结果输出：利用所获取的显著脉冲波峰的参数，计算其对应 激光点的三维空间坐标，输出点云。

点云的计算基于以下公式：

E＝E₀+dE×(WFOFFSET+t)

N＝N₀+dN×(WFOFFSET+t)

H＝H₀+dH×(WFOFFSET+t)

其中，E、N、H为波形点云三维坐标，E₀、N₀、H₀为波形第一次采样点所 对应的平面坐标，dE、dN、dH为波形的平面坐标微分，WFOFFSET为回波数据 第一次采样点到发射波形第一次采样点的偏差量，t为显著脉冲波峰对应的位置。 此外，根据实际需要，可以输出激光点所对应脉冲波峰的强度、波形宽度等信息。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术 领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的 方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。