-
地面激光扫描技术可以获得高精度、高密度的三维TLS点云，进而构建出精细的三维模型。相对于机载激光扫描技术，TLS具有显著的精度优势和时效优势, 因此在精密地形测绘与建筑景观建模中逐渐得到广泛的应用。基于TLS点云的三维模型构建往往需要各点的法向量信息以辅助点云的分割、滤波与三维表面重建等运算[1-2]。现有的点云处理软件(如CloudCompare、PCL、TerraSolid等)主要采用基于局部平面拟合的方法计算大规模点云的法向量。该方法由Hoppe等提出[3]，假设待建模表面处处光滑，则可通过主成分分析法(principal component analysis，PCA)计算邻居点集的最小二乘平面作为局部拟合平面，该平面的法向量即为该点处的点云法向量，所以也称基于PCA的点云法向量计算方法。
由于TLS点云中缺乏点与点之间的拓扑信息且采样点分布不均，因此对其法向量的准确估计具有一定的复杂性；并且在实际的TLS点云中，受地面覆盖物、多站匹配误差、采样密度差异等因素的影响，不可避免地存在着影响地形建模的粗差。基于PCA的点云法向量求解方法只能消除一般的随机误差，所以在实际应用中需要克服粗差对法向量计算的影响[4]。在基于PCA的法向量研究方法基础上，有学者又提出了基于加权最小二乘、高阶多项式拟合以及移动最小二乘的法向量求解方法[5-8]。文献[9]指出，由于少量的粗差点就可以对法向量的计算造成极大的干扰，所以拟合函数的改进对提高法向量计算结果的可靠性较为有限。对此，最理想的方法是通过点云滤波与分类算法将与地形重建无关的粗差剔除；但由于实际自然地形的复杂性以及TLS数据本身的局限，目前还没有广泛适用的方法能够完全有效地剔除影响地形建模的粗差。
为了克服粗差的影响，文献[10]将计算机视觉中的多结构检测思想引入到点云的处理中，采用随机抽样一致(random sample consensus，RANSAC)方法进行点云的平面检测、拟合与法向量的计算。该方法在含有噪声的点云中实现了大片连续平面的高效检测，但对点云中局部微小平面拟合的效率不高，且参数难以自适应设定。文献[11]采用统计学中的集成方法，通过在多个邻域尺度下进行点云法向量的估计，剔除差异较大的值后再由剩余值计算平均法向量, 以避免粗差的影响，但其计算量也随着邻域尺度的增大而成倍增加。文献[12]结合抗差噪声强度估计与核密度估计方法，提出了一种能够克服噪声与尖锐特征影响的抗差法向量计算方法，但主要针对的还是采样密度均匀的点云数据。现有的抗差法向量求解方法针对的多是人造模型与建筑物点云中的法向量计算，所依赖的先验条件与假设较多，不适用于场景复杂、密度变化较大、含有大量粗差的TLS点云。因此，有必要研究可靠高效的TLS点云抗差法向量求解方法。
本文通过对局部平面拟合的法向量计算方法的分析，基于邻居点集的协方差矩阵，提出了一种TLS点云的抗差法向量求解方法。该方法采用确定型最小广义方差估计方法(deterministic minimum covariance determinant，DetMCD)[13]与多元马氏距离对影响协方差计算的粗差进行识别与剔除，最终得到法向量的抗差估计。
HTML
-
对于点云P={pi∈R3, i=1, 2…n}中的某点pi=[xi yi zi]T，可以采用局部平面拟合的方法得到其法向量[3]。首先，寻找点pi的邻居点集NB(pi)。根据各点间的欧氏距离，采用最近k邻居(k-nearest neighbors，kNN)点集作为其定义[14]，即NB(pi)=kNN(pi, k)。根据各邻居点pj∈NB(pi)，计算该邻居点集的协方差矩阵Σi=(1/k)∑(pj－μi)(pj－μi)T，其中μi=(1/k)∑pj为该点集的中心向量。对协方差矩阵Σi进行特征值分解, 可以得到Σivj=λjvj(j=1, 2, 3), 其中，λj是协方差矩阵的特征值，vj为相应的特征向量。假设λ1 < λ2 < λ3，则v1即为邻居点集拟合平面的法向量，也就是点pi处的法向量。受到采样点密度分布、采样表面曲率、采样噪声等的影响，v1与真实表面的法向量并不完全一致。
基于局部平面拟合的法向量计算的关键是对邻居点集协方差矩阵Σi的求解。文献[15]研究了随机噪声对Σi的影响，给出了法向量计算误差的估计。但在真实的TLS点云中，除了一般的设备测量误差等随机误差外, 还存在着对协方差矩阵的计算有较大影响的粗差。对此，可以采用最小广义方差估计(minimum covariance determi nant，MCD)对其进行抗差估计[16]。对于一个大小为k的样本，MCD方法通过在所有的Ckh种大小为h的子样本组合中, 选取协方差矩阵的行列式(广义方差)最小的一个子样本组合, 以确保其中不含粗差。但是基于直接枚举的MCD方法计算复杂度很大，直到Fast-MCD的出现该方法才得以实际应用, 它通过随机抽取初始样本，用优化迭代获得总体样本的MCD估计。文献[17]采用基于Fast-MCD的抗差PCA方法实现了点云法向量的抗差估计, 然而Fast-MCD采用的随机化算法时间效率依然较低，因此难以满足大规模TLS点云数据的计算。
-
根据基于局部平面拟合的法向量求解方法与DetMCD方法，通过计算邻居点集协方差矩阵的抗差估计，得到TLS点云的抗差法向量。
1) 采用DetMCD方法计算最近k邻居点集合NB(pi, k)的中心向量和协方差矩阵的估计$ \left( {{{\mathit{\boldsymbol{\hat \mu }}}_{ir}}, {{\mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }}}_{ir}}} \right)$。作为Fast-MCD方法的改进，DetMCD采用了与其相同的优化迭代方法，但通过用确定型算法代替随机算法进行迭代初值的选取，提高了算法的收敛速度与稳定性。DetMCD方法的主要参数是优化迭代中保留的抗差样本比例H=h/N∈[0, 1]，它需要对总体样本N中的正常样本数量h进行估计，但对于实际的点云数据, 该参数难以预先确定。因此本文首先设定DetMCD的参数H=0.5，从而得到中心向量和协方差矩阵的保守估计，再结合多元马氏距离进行粗差剔除，确保算法对不同TLS点云的适应性。
2) 据多元统计学，假设各邻居点符合多元正态分布，则邻居点马氏距离的平方服从自由度为p的卡方(Chi-square)分布[18]。但传统的马氏距离直接根据邻居点集的原始中心向量与协方差矩阵进行计算，本身也受到粗差的干扰；而根据DetMCD得到的中心向量与协方差矩阵的保守估计$ \left( {{{\mathit{\boldsymbol{\hat \mu }}}_{ir}}, {{\mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }}}_{ir}}} \right) $可以排除大部分粗差的影响，由其与各邻居点pj∈NB(pi, k)的坐标可以计算得到各点的抗差马氏距离(robust Mahalanobis distance, RMD)：
根据RMD2(pj)~χ2(p)取显著性水平为α，可得到置信度为1-α的粗差判别限差：$ L = \sqrt {\chi _{3, 1 - \alpha }^2} $。由各邻居点的抗差马氏距离RMD(pj), 从统计意义出发，剔除影响协方差计算的粗差。
图 1为本文方法计算得到的一组邻居点的抗差马氏距离与原始马氏距离的比较。图 1中各点的横坐标为各邻居点的原始马氏距离，纵坐标为结合DetMCD得到的抗差马氏距离。两条红色直线为粗差判别限差：$ L = \sqrt {\chi _{3,0.975}^2} \approx 3.057\;5 $(置信度为0.975)。由图 1可知，原始马氏距离受粗差点的影响，只能识别出小部分粗差(红色虚线右侧的点)；而本文方法能够克服粗差影响，有效识别出大部分粗差(红色实线上方的点)。
Figure 1. Gross Error Identification Based on Original Mahalanobis Distance and Robust Mahalanobis Distance of Each Neighbor Point
剔除粗差后，将保留下的邻居点构成集合：
其协方差矩阵$ {\mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }}_{rs}} $就是邻居点集NB(pi, k)的抗差协方差估计。根据PCA方法，可由$ {\mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }}_{rs}} $得到点pi的抗差法向量Vnormal(pi)。
3) 对各点的法向量方向进行全局统一。根据TLS点云的特点，采用基于视点的法向量方向调整方法统一点云法向量方向。设TLS点云P的采样视点为vp=[xp yp zp]T，则可对各点抗差法向量ni∈Vnormal (P)的方向进行统一，即：
调整后各点的法向量方向相对于视点保持一致，从而实现了法向量方向的全局统一。对于从多个测站得到的点云数据，需要根据其各自的视点分别调整统一法向量方向。若没有对应的视点坐标，则需要采用法向量方向传播等相关算法对全局法向量方向进行统一[19]。
-
1) 读取TLS点云P：P=[p1 p2 … pn]，pi=[xi yi zi]T。
2) 遍历点pi∈P，搜索其最近k邻居点集NB(pi, k), 主要参数为邻居点集大小k。本文采用快速近似最近邻居搜索算法(fast library for approximate nearest neighbors，FLANN)进行最近k邻居点集的计算[14]。它通过随机k维森林算法与k均值优先搜索树算法对点云数据建立结构化索引，实现最近k邻居点集的快速查找。
3) 对邻居点集NB(pi, k)，采用DetMCD方法计算其中心向量与协方差矩阵的保守估计$ \left( {{{\mathit{\boldsymbol{\hat \mu }}}_{ir}}, {{\mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }}}_{ir}}} \right) $。DetMCD优化迭代过程中保留的抗差子集比例H=0.5。
4) 根据各邻居点的坐标pj∈NB(pi, k)与DetMCD计算得到的$ \left( {{{\mathit{\boldsymbol{\hat \mu }}}_{ir}}, {{\mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }}}_{ir}}} \right) $，由式(1)计算各邻居点的抗差马氏距离RMD(pj)。
5) 参照式(2)设定粗差判别限差L，识别并剔除粗差点。输入为各邻居点的抗差马氏距离RMD(NB(pi, k))，输出为剔除粗差点后的邻居点集NBR(pi)，参数为粗差判定的显著性水平α。
6) 根据剔除粗差点后的邻居点集NBR(pi)计算得到抗差协方差矩阵$ {{{\mathit{\boldsymbol{ \boldsymbol{\widehat \varSigma} }}}_{rs}}} $; 根据PCA方法得到点pi处的抗差法向量Vnormal(pi)。
7) 点云法向量的全局统一。对于计算得到的法向量Vnormal(pi)，根据地面激光扫描仪的视点vp和各点空间位置pi，使用式(3)统一各点的法向量方向。
该方法的参数主要有邻居点集的大小k与粗差判别的显著性水平α。邻居点集的大小需根据实际TLS点云数据的特点和应用需求确定，显著性水平α一般可取0.025。
2.1. 基于DetMCD的TLS点云抗差法向量计算原理
2.2. 基于DetMCD的TLS点云抗差法向量求解流程
-
将本文方法分别与基于PCA、鲁棒PCA(robust PCA，ROBPCA)[20]和RANSAC的法向量计算方法进行比较。通过构建模拟TLS点云分析各方法的抗差性与时间效率，并通过对真实TLS点云数据进行泊松表面重建，以评估该方法的实际应用效果。
各算法均使用Matlab语言实现，在Xeon E5工作站上进行实验。ROBPCA算法与RANSAC算法详细实现过程分别参见文献[20]和[21]。算法的精度由计算得到的各点法向量与真实法向量之间的夹角来度量。设点pi的真实法向量为n0，计算得到的法向量为ni，则点pi的法向量角度计算误差为：
由于实际TLS点云中各点的真实法向量难以获得，因此实验中真实法向量由构建的模拟点云数据预先确定。
-
为了检验本文方法的抗差性，构建包含粗差的平面TLS点云模拟模型进行实验。首先，根据TLS点云的采样密度和测量精度，构建一块2 m×2 m的正方形平面点云P={pi| pi=[pxi pyi pzi]T, i=1, 2…NP}。其具体定义为：
其中，NP=round(N·(1-G))为模拟TLS点云中平面点的规模，round(·)为取整函数；N为TLS点云总体规模(设为12 000)；G为粗差点比例，等于粗差点与点云总体规模的比值，取值范围为[0, 1]；S为模拟平面TLS点云边长(设为2 m)；rand (0, 1)是随机数，取值范围为[0, 1]；C为测量精度(设为0.01 m)。由于各模拟平面TLS点位于理想平面上，所以真实法向量均为n0=[0 0 1]T。
然后，在模拟平面TLS点云中添加随机分布的模拟粗差点P′={p′i| p′i=[p′xi p′yi p′zi]T, i=1, 2…NG)}。其定义为：
其中，NG=N-NP为粗差点规模；H为粗差分布高度(设为0.2 m)。最终由模拟平面点云P和模拟粗差点P′得到模拟TLS点云模型PT(N，G，S，C，H)=P+P′。根据该模型参数构建的模拟TLS点云PT(12 000，G，2，0.01，0.2)的密度和精度与一般TLS点云相符。但由于构建的模拟TLS点云内部的采样密度基本一致，为了更好地分析各方法在密度差异与粗差干扰下的法向量计算效果，本文随机选取模拟TLS点云模型边缘的点进行法向量的计算与比较。图 2为模拟TLS点云(粗差比例为30%)边缘选取的一组邻居点(k=70)。由图 2可见，受粗差和采样密度差异的影响，各方法计算得到的法向量相对真实法向量都有误差，但本文方法的法向量计算误差相对较小。
Figure 2. Normals of a Neighbor Point Set in Simulated TLS Point Cloud by Different Methods (with 30% Gross Error)
进一步以法向量角度θ的计算误差均值mean(θ)为指标，对比各方法对不同粗差比例的点云数据法向量计算误差。首先由模拟TLS点云模型PT的定义分别构建粗差点比例G={0%，10%…70%}共8组不同模拟TLS点云{PT0，PT1…PT7}=PT(12 000，G，2，0.01，0.2)。然后在生成的8组点云的边缘随机抽取各1 000个点构成测试点集{S0，S1…S7}。分别采用PCA方法、ROBPCA方法、RANSAC方法(参数σ分别设为0.01、0.02)、原始DetMCD方法(H=0.5)以及本文提出的方法计算{S0，S1…S7}中各点的法向量，各方法的邻居点集计算尺度k设为70。实验得到各方法在不同粗差比例下的法向量角度计算误差均值，如图 3所示。
Figure 3. Impact of Different Gross Error Rates on Mean of Normal Angle Errors Calculated by Different Methods
由图 3可见：①PCA方法在粗差比例大于10%的情况下，法向量角度误差均值已大于2°，而在粗差比例达到50%时，角度误差均值可达16°；②本文方法在粗差比例小于60%情况下的角度误差均值仍小于1°，在粗差比例大于65%后，角度误差均值才略有增加；③原始DetMCD方法由于保留的抗差子集大小固定，在粗差比例较低时误差相对偏大，不能适用于不同粗差情况的实际TLS点云；④RANSAC方法在参数σ=0.01时，抗差性能与ROBPCA方法接近，但当参数σ=0.02时误差急剧增大。在粗差比例较大时，甚至比PCA方法的误差还大，说明RANSAC方法的抗差性能对参数设定较为敏感，在实际应用中，参数σ的设定存在较大的困难；⑤ROBPCA方法与本文方法的抗差性能接近，但在粗差比例较大时，本文方法更优。
-
为了检验和比较各方法的时间效率，对不同规模的点云集分别测试各方法的运行时间。模拟TLS点云模型定义为：ST=PT(300 000，0.3，10，0.01，0.1)，在其边缘分别随机提取规模为{1 000，2 000…10 000}的10组测试点集{T1，T2…T10}。分别采用PCA方法、RANSAC方法、本文方法以及本文方法的并行优化改进方法进行比较实验(最近邻居点集尺度k均设为70)，得到各方法对不同规模点云的法向量计算时间，如图 4所示。由于ROBPCA方法对规模为1 000的点云数据的法向量计算时间已达到921 s，远远大于其他方法，难以适用于大规模点云数据，所以不在图 4中标示。由图 4可见：①本文方法的运行时间相对于RANSAC方法具有显著优势，但高于PCA方法，约为其19倍；②本文方法通过并行优化改进后，运行时间可降为PCA方法的2.6倍左右。因此，对于数据大小为100万的点云，采用并行优化改进的方法可以在约77 min内完成计算，能够满足大部分应用的时效需要。
Figure 4. Comparison of Running Time with Different Point Cloud Scales by Different Methods
-
为了验证本文方法的实际效果，采用野外采集的真实TLS点云数据进行实验验证。数据来自甘肃天水桥子沟的一条典型切沟。该切沟断面呈典型的Ⅴ字型，90%的地表被低矮豆科类植物和荒草覆盖。利用Leica HDS6100型三维激光扫描仪对切沟进行扫描并存储点云数据。该TLS点云包含激光点2 573 105个，覆盖面积为720.476 m2，平均密度为3 571.4个/m2。
对该TLS点云采用本文方法进行法向量计算，取邻居点集合大小k为20，置信度为0.975，得到的点云抗差法向量如图 5所示(为了便于显示，从点云中随机抽取了1 000个点)。分别根据本文方法得到的抗差法向量与PCA方法得到的法向量对该点云进行泊松表面重建，结果如图 6所示。由图 6(a)可见，由PCA方法求得的法向量重建得到的泊松表面在地形变化剧烈的区域存在较大的变形(如图 6(a)中4个标黑框的区域)；而由本文方法求得的法向量重建的泊松表面基本反映了实际地面的形态(图 6(b))，说明本文方法通过获取TLS点云的抗差法向量能够提高后续处理的可靠性。
Figure 5. Robust Normal Estimation of Real TLS Point Cloud by the Proposed Method
Figure 6. Poisson Surface Reconstruction Results Based on Normal Estimation by PCA and the Proposed Method
3.1. 模拟TLS点云实验1：抗差性能分析
3.2. 模拟TLS点云实验2：时间效率比较
3.3. 真实TLS点云实验：泊松表面重建
-
本文针对TLS点云中存在的粗差干扰与采样密度不均的问题, 根据局部平面拟合的原理，从邻居点集的协方差矩阵出发，提出了一种基于最小广义方差估计的TLS点云抗差法向量求解方法。该方法依赖的假设条件与先验知识较少，适用于场景复杂、采样密度差异较大的TLS点云抗差法向量计算。实验结果表明：①本文方法相比于PCA方法具有显著的抗差性优势，也优于ROBPCA方法和RANSAC方法；②本文方法的并行优化改进方法具有较高的时间效率，可以应用于大规模TLS点云的计算；③通过对真实TLS点云的法向量计算与泊松表面重建，进一步证明了其对实际TLS点云处理的有效性。为了更好地对复杂自然场景TLS点云的法向量进行求解，下一步将研究最优法向量计算尺度的自适应选择。
DownLoad:
