近几年，随着商业机载小光斑地形测量设备的不断涌现，一种新兴的机载全波形激光探测与测量(light detection and ranging, LiDAR)系统受到广泛关注^[1-2]，它具备激光脉冲信号的全数字化记录能力。相比传统机载LiDAR系统，它主要提供给用户的不再是离散三维点云数据，而是由许多一维脉冲采样集合组成的波形数据，通过对波形数据进行处理不仅得到离散点云，还能够分析出目标地物的返回特性信息(也称为波形信息)^[3]。全波形LiDAR技术不仅对植被及生物量信息有着显著作用，对激光地形测量也有着非常大的帮助，其中，高质量数字高程模型(digital elevation model, DEM)生产是比较典型的应用，利用全波形LiDAR可以得到精度更高的DEM产品。

利用全波形LiDAR数据生产高质量DEM包括波形数据分解和点云滤波两个关键步骤。现有波形分解算法主要包括非线性最小二乘法^[4-5]、期望最大化法^[6-7]与蒙特卡洛马尔科夫链算法^[8]，其中基于Levenberg Marquardt(LM)的非线性最小二乘法应用广泛，它灵活吸收了梯度法与牛顿法的优点，然而该算法对初始值的要求较高，容易陷入局部最优解。基于离散点云几何特性的滤波算法可以分为形态学方法^[9-10]、基于内插的方法^[11-12]与基于曲面约束的方法^[13]这3类^[14]。曲面约束法中，由于几何信息的局限，地面种子点的筛选会出现一定偏差，由此会带来对曲面拟合滤波效果的不利影响；利用波形信息的点云滤波研究还很少，国外学者通过简单设定脉冲半宽阈值进行了有关滤波尝试^[15-17]，验证了利用波形参数区分地面与非地面的可行性，但是将波形信息融入到滤波过程的算法仍有待进一步研究。

针对以上问题，本文提出了一种利用波形信息的加权曲面拟合LiDAR点云滤波方法。该方法利用全局收敛LM算法解决波形分解过程中参数解陷入局部最优的问题，利用波形参数与抗差估计理论筛选可靠的地面种子点，在曲面拟合过程中根据波形参数值设置合理权重，综合考虑窗口尺寸与曲面拟合中误差的影响，设置适用于各种地形的自适应高差阈值，从而能较好地克服上述不足, 如图 1所示。

图  1  算法基本流程图

Figure  1.  BasicFlowchart of Algorithm

下载: 全尺寸图片 幻灯片

1 波形分解

波形数据一般由两部分组成:一部分记录波形属性信息，包括发射脉冲数、激光波长、发射脉冲宽度以及首末采样点的三维坐标等；另一部分记录所有脉冲的采样数据。波形分解将波形分成多个回波，从而对脉冲下的不同目标进行特征描述，其中波形模拟与波形拟合是两个关键技术。通常，波形数据被看作为若干高斯分量的叠加^[18]，因此采用高斯函数模型对其进行模拟；波形拟合用以求解最优化参数，如前所述，在非线性最小二乘法的研究中均采用LM算法，该算法往往得到局部最优解，造成波形分解结果不准确。本文引入全局收敛LM算法来解决这一问题，基于全局收敛LM的波形分解步骤如下。

1) 波形预处理。激光脉冲扫描过程中由于天气等方面影响会产生背景噪声，这些噪声在波形中呈现小振幅抖动形态，如果不将其去除，波形误分解的情况很容易发生。

2) 波峰检测。计算每个采样点的二阶导数，将零点位置视为波峰初始位置，检测出的波峰数量决定波形参数个数。

3) 参数初始化。对于高斯函数模型，需要初始化的参数包括脉冲振幅A_k、脉冲距离μ_k和脉冲半宽ω_k。

4) 参数最优化。采用全局收敛LM算法实现波形的拟合，即得到波形参数的全局最优解，算法原理详见本文1.1节。

5) 计算波形残差。计算拟合波形后的残差ξ，若残差值小于前一次的计算结果，则对残差波形进行步骤2)~5)的处理，如此迭代计算直至残差值停止减小。残差ξ的计算公式为：

其中，n是采样点的个数；m为参数个数；y是目标波形中各个采样点的振幅值；f(x)是拟合波形中采样点的振幅值。

6) 生成三维点云。迭代完成后得到的是组合高斯函数的模型参数，它们也就是描述波形特征的波形参数。每个波形分量对应一个地物反射回波，因此由该波形分量参数能够计算对应地物的空间三维坐标，计算方法见1.2节。

1.1   全局收敛LM算法

全局收敛LM是传统LM的改进算法，该算法能够弥补传统LM局部收敛的问题，其在燃气轮机部件特性曲线的拟合优化中有显著效果^[19]。本文将全局收敛LM算法引入到激光雷达全波形数据的波形拟合中，不仅拓展了全局收敛LM算法的应用领域，也为激光雷达全波形数据分解提供了更加可靠的解决方案。

传统LM算法的计算式为：

其中，f(x_i, p)为待定系数(p₁, p₂…p_m)的函数，其观测数据为(x_i, y_i), i=1, 2…n；p⁽⁰⁾=(p₁⁽⁰⁾, p₂⁽⁰⁾…p_m⁽⁰⁾)为函数参数的初始值；λ为阻尼系数，当λ=0时，即为高斯牛顿法的形式；J(x, p)为函数f(x, p)的一阶偏导数组成的雅克比矩阵H(x, p)为函数f(x, p)的海塞矩阵。在解算过程中，若${\hat{p}}$与p⁽⁰⁾中各元素之差的绝对值很小，则迭代结束；若‖${\hat{p}}$-p⁽⁰⁾‖较大，则将解算得到的${\hat{p}}$作为新的p⁽⁰⁾重新进行计算，如此反复迭代直到‖${\hat{p}}$-p⁽⁰⁾‖达到收敛条件或迭代次数达到最大为止。

假设解集非空且记为${\hat{p}}$，若J(x, ${\hat{p}}$)非奇异且初始点p⁽⁰⁾与${\hat{p}}$充分靠近时，LM算法产生的迭代点列二阶收敛于${\hat{p}}$。然而，在实际运用中，雅克比矩阵非奇异这一条件往往过于严格，而且非奇异性条件表明该方法得到的解是局部唯一的。因此，采用一种基于信赖域技巧的全局收敛LM算法使得解具有全局收敛性。算法选取阻尼系数λ_k=α_k[θ‖F_k‖+(1-θ)‖J_k^TF_k‖], θ∈[0, 1]其中目标函数F_k=y-f(x, p^(k))；k表示当前迭代次数；α_k利用信赖域技巧来修正。第k步迭代的实际下降量A_k和预估下降量P_k分别为：

其中，Δp^(k)=p^(k+1)-p^(k)，用实际下降量A_k和预估下降量P_k的比值r_k决定是否接受试探步Δp^(k)以及调整迭代参数中α_k因子的大小。一般来说，r_k越大，说明‖F‖²下降得越多，因此接受Δp^(k)，期望下一试探步Δp^(k+1)更长，故减小α_k；反之，r_k越小，则考虑拒绝接受Δp^(k)，从而增大α_k。

1.2   生成三维点云

利用波形分解出的脉冲距离μ_k与波形数据属性信息可以计算对应回波的地物三维坐标^{[7, 20]}，公式为：

其中(E₀，N₀，H₀)是起始脉冲波形首个采样的三维坐标；(dE，dN，dH)为波形采样单元的坐标微分；W是从起始脉冲波形首个采样到返回脉冲波形首个采样的偏移量。

2 利用波形信息的滤波

点云滤波是指将分布在地形表面不同目标上、不规则分布的离散点进行处理和识别，即区分地面点云(裸地与道路)和非地面点云(建筑物、车辆与植被等)^[10]。传统滤波算法往往基于离散点云的空间几何特征，它们对于低矮植被的滤除收效甚微，利用波形分解得到的波形参数将有利于提升滤波效果。国外学者利用波形信息进行了简单滤波尝试，在山地区域验证了利用脉冲半宽区分地面与植被部分的可行性，然而，更充分地利用波形信息的滤波方法还有待进一步挖掘，对城市等其他区域的实验也有待进一步验证。传统曲面拟合滤波算法^[13]通过选取窗口内的高程最低点作为地面种子点，并由多个窗口种子点拟合区块地形，当窗口内包含大型建筑物时，建筑物脚点会被误选为地面种子点，区块地形的拟合会出现偏差。

针对以上问题，本文采用波形分解得到的波形参数与抗差估计理论筛选可靠的地面种子点，设计了基于波形参数的加权曲面拟合与高差阈值自适应确定方法来弥补传统方法的不足。该方法的具体流程如图 2所示，根据传统曲面约束法的多级窗口滤波思想，它也是一个多级窗口的迭代处理过程，随着窗口的不断增大，逐步实现滤波处理，直至达到最大窗口尺寸为止。关键技术包括异常种子点检测、加权曲面拟合以及自适应高差阈值确定。

图  2  滤波流程图

Figure  2.  Flowchart of Filtering

下载: 全尺寸图片 幻灯片

2.1   异常种子点检测

文献[16]指出，脉冲半宽对植被与建筑边缘具备较强的识别能力。例如，对于Riegl LMS-Q560全波形数据，脉冲半宽大于1.70 ns的脚点属于植被或建筑边缘的可能性大；反之，属于平坦表面的可能性大。脉冲半宽阈值ω由系统规格参数中给定的脉冲宽度Ω决定，计算公式为$\omega =\mathit{\Omega /}\rm{2}\sqrt{2\ln 2}$。

本文首先利用脉冲半宽对种子点进行初步筛选。接着，引入抗差估计理论对异常种子点作进一步检测。抗差迭代的初值要保证高崩溃污染率^[21]，因此，利用高崩溃污染率的抗差估计方法进行趋势面模型的解算，并根据抗差迭代后每个种子点对应的等价权来检测异常。具体步骤^[22]如下。

1) 首先计算强淘汰权函数，计算式为：

其中，p_i⁰为第i个种子点的强淘汰权函数值；残差v_i表示第i个种子点高程与所有种子点高程中位数之差，即v_i=z_i-med{z_i}；${{{\hat{\sigma }}}_{i}}$是v_i的均方差因子抗差解，${{{\hat{\sigma }}}_{i}}$=med{|v_i|}/0.674 5；c₀为权函数临界值，通常设置为1.0或1.5。

2) 利用式(6)获得的p_i⁰计算趋势面模型的初值。设有误差方程V = AX -L，根据抗差估计理论，A^TPV =0，系数向量的第(k+1)步迭代估值为：

式中，迭代等价权矩阵为：

其中，p_i^(k)为各种子点对应的等价权，利用IGGΙΙΙ权函数确定：

其中，标准化残差v′_i^(k)=v_i^(k)/$\hat{\sigma }_{i}^{\left( k \right)}$，$\hat{\sigma }_{i}^{\left( k \right)}$=med{|v_i^(k)|}/0.674 5；临界值k₀与k₁的取值范围分别为1.0~2.0和2.5~4.0。

3) 迭代计算的停止条件为max|x^(k+1)-x^(k)|＜ε。迭代完成后，将等价权值为0的种子点标记为异常。

抗差估计在迭代过程中，通过判断标准化残差的大小，不断减小不可信种子点的权函数值，从而提高了检测异常种子点的准确性。

2.2   加权曲面拟合

异常种子点检测后，利用区块中保留的地面种子点进一步拟合地形。为获取可靠的地形模拟效果，在曲面拟合中依据脉冲半宽设置相应种子点权值。波形分解得到的脉冲半宽ω反映照射目标表面的平整程度，ω越小，该点位于平坦表面的可能性越大，因此，设计权值计算式为：

其中，p_j为区域内第j个种子点的权值，其脉冲半宽为ω_j；n是区域内种子点总数。按照这种方法设置权重的意义在于增大脉冲半宽小的点对曲面拟合的贡献，降低脉冲半宽大的点对曲面拟合的影响，从而得到尽可能真实的地形拟合结果。

2.3   自适应高差阈值

区块内离散点与拟合曲面的高程差异将最终决定其属性，当高差大于阈值T时，该点标记为非地面点；反之，记录为地面点。然而，不同区块地形各异，高差阈值设置为常数显然是不合适的。因此，本文综合考虑窗口尺寸与曲面拟合中误差的影响，提出了自适应确定高差阈值的方法，如式(11)所示。基本思想是窗口尺寸与中误差越小，认为拟合结果越能够表达真实地形，拟合精度越高，此时高差阈值可以设置相对严格，以保证低矮地物脚点被有效滤除；反之，随着窗口尺寸与中误差增大，其表现真实地形的能力逐步减弱，说明拟合地形偏差增大，因此为避免地面脚点的错误滤波，高差阈值应设置相对宽松。

式中，d_normalized为归一化的窗口尺寸；m_normalized是归一化的中误差；窗口权重P_d与中误差权重P_m均设置为0.5。如果高差阈值设置得过小，将会牺牲一部分地面脚点，使得地形不连续；如果设置过大，则会保留部分非地面点，从而达不到滤波效果，因此本文将高差阈值限定在0.3~2 m之间，即T_min为0.3，T_max为2。

3 点云滤波实验与分析

采用中国黑河综合遥感联合实验中的小光斑激光雷达波形数据进行实验，该数据集由RIEGL公司的LMS-Q560机载全波形LiDAR测量系统获取，激光波束发散角为0.5 mrad，激光波长1 064 nm，发射脉冲宽度为10 ns。在中国张掖-盈科飞行区选取城市区域与耕地区域的实验数据；在中国大野口森林区选取山地区域实验数据，数据如表 1所示。

表  1  实验数据基本情况

Table  1.  BasicInformation of Test Data

实验数据	类型	脉冲数	描述
测区1	城市区域	188 305	主要包含建筑物、车辆与少量植被
测区2	耕地区域	117 065	耕地与低矮植被
测区3	山地区域	70 201	较大地形起伏，高植被

下载: 导出CSV 

| 显示表格

3.1   波形分解实验与分析

图 3为单条脉冲的分解对比结果，其中图 3(a)为传统LM算法的分解效果，图 3(b)是全局收敛LM算法得到的结果。图 3中黑色实心圆点是原始波形的采样数据，红色虚线表示分解出的波形分量，蓝色实线为模拟波形结果，横轴与纵轴分别表示波形采样时间与振幅值。表 2给出了对应波形分量参数，并分别计算了波形分解后的残差ξ。由图 3可知，全局收敛LM算法较传统LM算法更加可靠，模拟波形更准确。表 2中，传统LM对该脉冲仅分解出2个波形分量，这是局部收敛造成的结果，相比之下，全局收敛LM算法得到4个波形分量，并且残差由3.612减小至0.158。

图  3  脉冲2波形分解结果

Figure  3.  Decomposing Results of Pulse 2

下载: 全尺寸图片 幻灯片

表  2  脉冲2波形分量参数与残差

Table  2.  WaveformParameters and Residual of Pulse 2

算法	波形分量数	脉冲振幅Ak	脉冲距离μk	脉冲半宽ωk	残差ξ
传统LM	2	29.368	14.156	1.746	3.612
		8.142	19.143	1.562
全局收敛LM	4	29.412	14.167	1.753	0.158
		5.935	18.325	1.075
		6.852	22.343	2.435
		6.088	24.214	1.245

下载: 导出CSV 

| 显示表格

3.2   点云生成实验与分析

点云生成效果是评价波形分解的另一重要指标。分解出的波形分量越多，点云的层次感越强。图 4与图 5分别为两个不同剖面的对比结果。图 4剖面1是城市数据中一颗树木的剖面，图 5剖面2是山区中一个地形条带的剖面。从图 4剖面1的对比结果可以看出，标记A的是树冠部分，此部分在返回脉冲中位置靠前；标记B的为下部枝条部分，该部分在返回脉冲中位置靠后。可以看出，全局收敛LM方法在A部分有非常明显的改进，波形分解生成的脚点更加密集，在B区域也分解出了更多的脚点。对于图 5剖面2，标记A与标记B分别为山坡的下部与上部区域，全局收敛方法在这两个标记区域同样分解出更多激光点云。这进一步反映出了全局收敛LM方法在点云层次表达上的优势，同时也验证了该方法在波形分解结果上的先进性与可靠性。

图  4  剖面1点云对比

Figure  4.  Comparison of Cross 1

下载: 全尺寸图片 幻灯片

图  5  剖面2点云对比

Figure  5.  Comparison of Cross 2

下载: 全尺寸图片 幻灯片

3.3   滤波实验与分析

对这3个测区(表 1)的滤波结果进行精度统计，分别统计传统方法、渐进加密三角网滤波方法以及本文方法的Ⅰ类错误率、Ⅱ类错误率以及总错误率。Ⅰ类错误率是指地面点误分为非地面点的比率；Ⅱ类错误率为非地面点误判断为地面点的比率；总错误率是Ⅰ、Ⅱ类错误率的加权求和。Ⅰ、Ⅱ类错误率表现算法的适应性，总错误率反映算法的可行性^[14]。参考滤波结果由人工编辑获取。

无论是曲面拟合法还是渐进加密三角网滤波方法，最大窗口尺寸都应根据地物具体情况进行设置。城市数据中将最大窗口尺寸设置为90 m；耕地数据的最大窗口尺寸设为50 m；山地数据设置为30 m。为保证对比实验的公正性，在同一测区均使用相同的最大窗口尺寸。具体实验结果如表 3~5所示。

表  3  城市数据点云滤波结果

Table  3.  Filtering Results of Urban Area

方法	参考结果	滤波结果		参考总计	错误率(%)
		地面点	非地面点
传统方法	地面点	76 941	30 192	107 133	Ⅰ类	28.18
	非地面点	2 378	110 853	113 231	Ⅱ类	2.10
	滤波总计	79 319	141 045	220 364	总错误率	14.78
渐进加密三角网滤波法	地面点	95 956	11 177	107 133	Ⅰ类	10.43
	非地面点	340	112 891	113 231	Ⅱ类	0.30
	滤波总计	96 296	124 068	220 364	总错误率	5.23
本文方法	地面点	100 377	6 756	107 133	Ⅰ类	6.31
	非地面点	283	112 948	113 231	Ⅱ类	0.25
	滤波总计	100 660	119 704	220 364	总错误率	3.19

下载: 导出CSV 

| 显示表格

表  4  耕地数据点云滤波结果

Table  4.  Filtering Results of Farmland Area

方法	参考结果	滤波结果		参考总计	错误率(%)
		地面点	非地面点
传统方法	地面点	79 078	24 873	103 951	Ⅰ类	23.93
	非地面点	3 670	41 064	44 734	Ⅱ类	8.20
	滤波总计	82 748	65 937	148 685	总错误率	19.20
渐进加密三角网滤波法	地面点	94 752	9 199	103 951	Ⅰ类	8.85
	非地面点	1 588	43 146	44 734	Ⅱ类	3.55
	滤波总计	96 340	52 345	148 685	总错误率	7.25
本文方法	地面点	94 821	9 130	103 951	Ⅰ类	8.78
	非地面点	1 538	43 196	44 734	Ⅱ类	3.44
	滤波总计	96 359	52 326	148 685	总错误率	7.17

下载: 导出CSV 

| 显示表格

表  5  山地数据点云滤波结果

Table  5.  Filtering Results of Mountain Area

方法	参考结果	滤波结果		参考总计	错误率(%)
		地面点	非地面点
传统方法	地面点	52 401	2 077	54 478	Ⅰ类	3.81
	非地面点	21 205	68 853	90 058	Ⅱ类	23.55
	滤波总计	73 606	70 930	144 536	总错误率	16.11
渐进加密三角网滤波法	地面点	45 544	8 934	54 478	Ⅰ类	16.40
	非地面点	13 232	76 826	90 058	Ⅱ类	14.69
	滤波总计	58 776	85 760	144 536	总错误率	15.34
本文方法	地面点	45 667	8 811	54 478	Ⅰ类	16.17
	非地面点	12 752	77 306	90 058	Ⅱ类	14.16
	滤波总计	58 419	86 117	144 536	总错误率	14.92

下载: 导出CSV 

| 显示表格

实验结果表明，本文方法在城市数据以及耕地数据的Ⅰ类、Ⅱ类和总错误率相比传统方法均有大幅度提高。这2个测区均属于平坦地形，传统方法中，建筑物与密集植被上的错误种子点使得拟合曲面与实际地形不符，从而错误滤除了大量地面脚点，而部分非地面点却无法滤除，因此其Ⅰ类与Ⅱ类错误率均较高，本文方法通过异常种子点检测很好地弥补了以上不足。在山地数据中，由于地形起伏较大，传统方法设置的高差阈值也会相应较大，因此，其滤波结果会保留很多非地面点，而Ⅰ类错误率小只是牺牲Ⅱ类错误率所造成的“假象”，本文方法通过波形信息的融入以及高差阈值的自适应设置使得该问题得以改善。对比渐进加密三角网滤波方法的结果可知，本文方法在城市数据方面存在一定优势，耕地数据与山地数据的滤波精度与渐进加密三角网滤波方法相当。

4 结语

全波形LiDAR是当前摄影测量与遥感领域的新技术，相比传统LiDAR，它的最大特点在于对回波脉冲的完整记录，波形处理得到的信息有利于高质量DEM的生成。目前，利用波形信息的滤波方法研究还不深入，传统LM波形分解算法容易陷入局部最优，曲面约束法在选取地面种子点拟合曲面时存在缺陷。针对以上问题，本文提出了一种利用波形信息的加权曲面拟合滤波方法。相比传统方法，它具有波形分解更可靠、滤波精度更高的特点，能够较好解决波形分解局部收敛、地面种子点错误选取、地形曲面拟合偏差以及高差阈值难以设定的问题，具有较强实用价值。波形信息不仅可用于DEM生产，也可用于点云数据精细分类，深入挖掘波形参数特征将进一步提高分类精度，这也是未来研究的重点。