基于地形断裂线约束的机载激光点云高精度滤波方法

杨宇妍,臧玉府,肖雄武,管海燕,彭代峰

1. 南京信息工程大学遥感与测绘工程学院,江苏 南京 210044; 2. 武汉大学测绘遥感信息工程国家重点实验室,湖北 武汉 430079

机载激光扫描系统可高效、准确地获取大范围区域的空间点云数据,其中地面点是建立高精度数字高程模型的重要数据源,已在地形地貌、水土保持、土壤侵蚀分析等领域中得到广泛应用[1]。机载激光点云滤波作为识别地面点的主要手段,对地面信息提取、道路规划及目标精准识别等应用具有重要作用。但由于点云数量巨大,且地物易被遮挡,使得处理效率低下,目标分类困难,对现有滤波方法提出了挑战[2]。

目前,机载激光点云的滤波方法主要包括4类[3]:①基于坡度变化的方法。逐点计算点与相邻点之间的坡度,通过比较两点间实际坡度与坡度阈值滤除地物点[4-6]。但这类方法在处理陡坡、高程突变区域时会错误滤除大量地面点,且坡度越陡滤波精度越低。②数学形态学方法。对地形进行腐蚀膨胀,构造开运算,利用开运算后非地面点高程变化大这一特点,将高差大于阈值的点归为非地面点并进行滤除[7-9]。这类方法得到的滤波结果依赖于最佳窗口尺寸的选择,较小的窗户尺寸会保留较大的建筑物,较大的尺寸会忽视地形细节。③基于曲面的方法。通过迭代从原始数据集中选择地面测量值,逐渐逼近地面,创建一个近似裸露的地球表面[1,9-13]。文献[12]利用种子点创建稀疏三角网,并在迭代过程中加密,用来处理不连续的表面。文献[13]提出基于渐进三角网滤波方法,它在经典三角网滤波的基础上,获取更具可靠性的种子点,以增强对不同地形的适应性。不同于传统基于曲面的方法,布料模拟滤波(cloth simulation filter,CSF)[14]将点云翻转,假设有一块布料受到重力落下,最终落下的布代表当前地形。此类方法在地形变化平缓的地区效果较好、所需设置的参数少。④基于采样线的滤波方法。此类方法因能顾及整个区域的地形特征,高效地处理大范围场景区域,为机载激光点云滤波提供了一种思路。文献[15]中用球面坐标代替笛卡儿坐标,利用方位角和径向距离球面坐标作为判断标准,通过球坐标中径向距离曲线的断点和转折点分割地面点和非地面点,可有效避免过度分割和欠分割问题。文献[16]针对每条采样线,利用半径可变的圆环滚过,保留采样线上被圆环滚过的点,采用B样条拟合曲面,通过比较拟合高程与实际高程的差值得到地面点。文献[17]提出一种基于采样线处理的地面滤波算法,利用一维迭代样条插值算法改进地面点的提取。

现有滤波模型虽能基于一维地形特征实现快速滤波,但不能避免断裂线周围的大量地面点被错误滤除,难以适用于由断裂线引起的高程突变区。针对以上问题,本文基于采样线滤波优势,构建融合断裂线约束的采样线滤波模型,使其在滤波过程中较大程度保留地形断裂处的地面点,适用于多种复杂地面场景的地物点滤除,提高滤波的完整度和准确度。

1 基于断裂线的激光点云滤波

本文通过建立断裂线约束模型提出一种高效的采样线滤波方法。首先,基于点曲率及邻域范围内距离之和创建约束能量项,改进Snake能量模型[18],并通过能量函数最小化提取断裂线特征;然后,通过构建正三角形模型,计算点云平坦度及融合断裂线约束进行采样线滤波;最后,根据采样线上的地面点利用改进最小二乘法拟合全局曲面,精准滤除地物点。该方法主要流程如图1所示。

1.1 地形断裂处点云滤波问题

采集区域内通常场景复杂,存在大量的高程突变区域,严重影响了现有滤波方法的性能,一直是点云滤波中的重难点。研究发现点云中高程突变区域附近往往存在显著的断裂线特征,为充分利用该信息,本文按地形表面是否连续的特点,将断裂线特征分为以下两类:山谷、山脊等处地形表面连续变化且存在明显线特征的Ⅰ类断裂线(图2(a)),以及在人工建筑、陡坎、冲沟等地表不连续且存在明显高程突变的Ⅱ类断裂线(图2(b))。图2中黄色曲线为地形断裂线特征。现有滤波方法大多没有综合考虑上述地形断裂线的影响,使得地形突变区域的地面点或地物点被错误地滤除或保留。针对上述问题,本文提出一种融合断裂线约束模型的采样线高效滤波方法,提高滤波精度、避免高程突变区域的地形失真。

图2 断裂线Fig.2 Breaklines

1.2 地形断裂线提取

曲率是反映地形表面起伏程度的重要度量,断裂线上的点大多具有局部曲率极值。本文基于点曲率及邻域范围内距离之和创建约束能量项,改进Snake能量模型,并通过能量函数最小化提取断裂线特征。对点云中任意一点及其邻域点,采用最小二乘拟合法实现其局部曲面参数化,设该点局部曲面函数为f(xsur,ysur),其一阶、二阶偏导数分别表示为rx、ry、rxx、rxy、ryy,则该点处的平均曲率H(s)[19]为

(1)

(2)

式中,Eint表示模型内部能量,控制曲线的平滑性和连续性;Eext表示模型外部能量,吸引初始线向着曲线延伸或伸缩;Econ表示模型约束能量,使初始断裂线向外部能量较弱的地方继续演化。内部能量Eint、外部能量Eext和约束能量Econ构建如下

(3)

Eext[l(s)]=ω(s)·|u(s)|

(4)

(5)

通过循环迭代求解能量函数最小值[22],排除噪声点的影响、优化Snake曲线,得到最终的断裂线,如图3所示。

图3 部分断裂线提取结果Fig.3 Results of partial breaklines extraction

1.3 融合断裂线约束的滤波模型

为满足大范围点云滤波需求,本文由测区的机载点云生成多条采样线,将空间滤波问题转换为一维滤波处理。在点云场景范围内沿纵坐标方向等间距对点云进行采样,获取点序列构成的采样线(图4),再将每根采样线沿纵轴方向上给定距离d内的点归入对应的采样线,将采样线上的点按横坐标值升序排列。距离d的大小会影响采样线准确性:若给定距离过大,则采样线上的点数量增多,造成数据冗余,继而影响采样线滤波的效率;反之,采样线上的点数量过少,采样线丧失连续性,不但降低采样线滤波精度,同时也使得曲面拟合丧失真实性,故本文通过式(6)计算距离d[15]

(6)

图4 采样线Fig.4 Sampling line

式中,num为点云数量;S和W分别表示点云所表示地形范围的长和宽。

由图4可知,本文求得的采样线连贯、准确,能保证后续操作顺利进行。

1.3.1 预处理

通过构建等边三角形模型、计算点平坦度及邻域范围内高差进行采样线粗滤波(图5),滤除采样线上的车辆、植被、建筑物等非地面点,再选取粗滤波后采样线上的最左端点为采样中心点,取距采样中心点最近的点为基准点。

图5 采样线粗过滤Fig.5 Coarse filtering of a sampling line

(7)

(8)

在△O1P2O2中,运用三角余弦定理计算O1P2和P2O2的夹角θ,判断θ与夹角阈值θthreshold的关系:若θ大于θthreshold,则P3为非地面点,继续取采样线上与P2相邻下一点作为目标点,重复上述步骤;反之,P3为地面点,用P2更新P1、P3更新P2继续判断后续点,直至采样线终点。

滤除采样线上大部分非地面点后,需要利用点云平坦度λ[23]进一步过滤植被、车辆和人员等地物。本文根据上述处理后采样线点局部邻域R内原始点Xi(0

(9)

(10)

计算点云平坦度滤除采样线上的植被、车辆和人员等地物后,需要考虑建筑物等局部邻域平坦的非地面点对采样中心选取的影响[24]。为剔除采样线上的建筑物点,对点云中任一点,计算统计高程大于该点的邻域点数与全部邻域点数的比值,比值越小,该点为地面点的可能性越大。对每条采样线进行粗过滤后,采样线上最左端点为采样中心,距采样中心最近的点为基准点。

1.3.2 基于断裂线的采样线点云遍历

对任一条采样线精细滤波,基于断裂线约束,结合采样中心点、基准点和各采样点构建有向向量计算夹角[25],滤除地物点。确定采样中心点和基准点后,从原始采样线最左端开始遍历,由近及远按点云顺序计算坡度ρ,得到采样线上的地面点。坡度ρ的计算公式为

(11)

图6为基于断裂线的采样线点云遍历示意图,具体的遍历步骤如下。

图6 采样线点云遍历Fig.6 A sampling line point cloud traversal

∥输入:采样中心点O、基准点B、采样点数n

fori←0 ton

ifPi==断裂点

保留Pi

ifPi!=断裂点

连接BPi、OB,根据式(11)计算OB和BPi的夹角ρ

ifρ≤ρthreshold(阈值)

保留Pi,Pi=B,i++

else

滤除Pi,i++

∥输出:一维地面点(GRPts)

基于断裂线约束的采样线滤波结果如图7所示,由图7(b)可知,本文方法能准确区分采样线上蓝色地面点与红色非地面点,确保后续处理顺利执行。

图7 采样线滤波前后对比结果Fig.7 Comparison before and after sampling lines filtering

1.4 改进的二次曲面拟合

由于地表复杂多样,简单的二次曲面拟合方法不能真实表示滤波后的地表形态,在连续性上存在较大缺陷。为较大程度还原地面,保证地表的真实性和连续性,本文格网化整个测区(格网单元尺寸原则上应大于最大地物的边长),且各格网的边界区域相交重叠(图8)。图8中黄色区域表示任一格网的不重叠区域,绿色区域表示两个格网的重叠部分,蓝色区域表示4个格网的重叠部分。先根据落在单元格(图8中虚线矩形)内的处理后的m个采样线上的地面点,运用最小二乘原理计算误差平方和σ2达到最小时的近似地面参数A、B、C、D、E、F,依据地面参数计算点云的拟合高程z(x,y)。σ2和z(x,y)的计算公式为

图8 格网划分Fig.8 Grids division

(12)

z(x,y)=Ax2+By2+Cxy+Dx+Ey+F

(13)

式中,σ2为计算值与真实值之间的误差平方和;z(x,y)为点云在二次曲面上的高程。判断原始点云中每个点所在的格网位置:若原始点落入非重叠区域,则依据式(13)计算点云的拟合高程z(x,y);若点落入格网重叠区域,则计算点在不同格网曲面拟合高程z(x,y)的平均值。比较拟合值与实际高程的差值,将差值大于阈值的点从候选地形点中去除。

2 试验结果与分析

为验证本文方法的有效性,采用国际摄影测量与遥感学会(ISPRS)工作组提供的数据集和国内成都山区点云两套数据进行测试。参考数据集是通过手动过滤点云数据集生成的,地形样本中的每个点被分类为地面点或地物点,后续将手动分类结果作为真值检验方法精度。

2.1 ISPRS试验数据

对于ISPRS提供的数据,本文选择了4组具有不同地形特征的数据集对方法性能进行测试,详细信息见表1。

表1 地形特征

2.1.1 试验结果

为了评估方法的性能,本文计算了4个样本的Ⅰ型误差、Ⅱ型误差和总误差[1]。此外,还通过Kappa系数[26]衡量分类结果的总体一致性。各型误差、Kappa系数和处理时间见表2。结果表明,对于不同的地形特征,本文方法均能有效提取地面点,并将总滤波误差控制在5%左右。其中地形2的总误差最小,这表明本文方法在存在大型建筑和桥梁、且地形平坦的情形下,性能依然优异。本文方法在处理Ⅰ型误差时,表现尤为突出,均控制在2%左右,这表明该方法能较大程度地保留地面点,还原地面形态。4个地形的滤波时间分别为22.6、18.3、25.4、8.3 s,可见本文方法兼顾了精度和效率。总体而言,本文方法对地形的起伏较为敏感,能够处理各种地貌类型,并且在斜坡、建筑边缘起伏处表现良好。

表2 场景误差分析

试验中4个典型地形滤波结果如图9所示。图9(a)、(c)、(e)、(g)为滤波前地形,图9(b)、(d)、(f)、(h)展示了各类地形中附着在DEM表面上的Ⅱ型误差(红色点),DEM由地面点生成。从图9中对不同地形的滤波结果看,疏密程度不同的地表植被和各种形状的建筑物基本被滤除,同时保留了地形特征和地面细节。

图9 滤波结果Fig.9 Filtering results

为了全面分析本文方法的准确性,将本文方法的总误差和Kappa系数与现有方法[9,14,27-29]进行了比较。这些方法和本文方法的总误差与Kappa系数如图10—图11所示。总体而言,对于不同特征的地形,本文方法通用性好并且滤波误差更小。除地形3外,本文方法与传统方法相比滤波精度提高了3%～5%。Kappa系数是统计学中度量一致性的指标。Kappa系数越大,表示模型预测结果与实际滤波结果一致性越强;若Kappa系数在[0.81,1]区间,则说明模拟结果与实际结果几乎完全一致,本文方法计算得到的Kappa系数均高于88%。总误差和Kappa系数的比较结果,说明本文方法准确性、精度与现有方法比较均表现优异。

图10 本文方法与其他方法的总误差比较Fig.10 Total error compared between the proposed method and other reported methods

图11 本文方法与其他方法的Kappa系数比较Fig.11 Kappa coefficient compared between the proposed method and other reported methods

2.1.2 参数设置

本文将通过大量测试确定的普适参数设置为固定值,在此基础上,针对不同样本确定采样线间距、坡度、格网边长及曲率计算时的邻域范围等参数的最优值。为了定量评估采样线间距的影响,本文测试了具有不同间距的样本(从6～20 m,步长为2 m)。每组总误差(取决采样线间距)如图12所示。数据表明,所有组的总误差基本随着采样线间距增大而增大,只有地形1在间距最大(20 m)的时候,取得最好的滤波效果。若采样线间距太小,会得到大量采样线点云,将失去简化数据、节省时间的意义。若采样线间距太大,后续生成的格网中可能不存在处理后的地面点,导致曲面拟合失败,无法对全域进行滤波。因此,大多样本均选择6 m作为采样线间距值,因为它对所有样本都产生了相对较好的结果,并且可以应用于许多情况而无须调整。

滤波过程中采样线上地物点滤除效果与坡度值设置有较强的关系。图13显示了不同坡度值下的总误差变化,可以看出,随着坡度的增加,滤波总误差也随之增大。所以所有样本的最佳坡度值均为20°。

图13 坡度阈值对应的总误差Fig.13 Total errors for each slope threshold

曲面拟合过程中格网边长将影响拟合效果(图14)。如图14所示,若格网边长设置过小,格网内采样线上的地面点数可能为0,将无法求解多项式系数,导致拟合曲面失败;若格网边长设置过大,拟合的曲面将不能反映真实的地表情况。因此,大多地形选择边长为40 m的格网,此时的拟合效果最佳。

图14 格网边长对应的总误差Fig.14 Total errors for each grid size

除上述参数外,断裂线的完整性与滤波结果紧密相关。本文采用改进Snake模型提取断裂线,曲率作用于约束能量,故计算曲率时邻域的选取制约着断裂线的提取(图15)。由图15可知,若邻域范围选取过大,造成断裂线冗余,将会影响采样线滤波结果,从而导致拟合曲面过高、错误滤除地面点;相反,断裂线提取不完整将不能滤除高程突变区的地物点。

图15 曲率计算时邻域范围对应的总误差Fig.15 Total error corresponding to different neighborhood in curvature calculation

2.1.3 消融研究

对于场景中包含许多陡坡或建筑的地形,断裂线是影响本文方法准确性的重要因素。假设不考虑断裂线对滤波结果的影响,那么在陡坡、建筑物边缘处的地面点将被错误地移除,影响最终的滤波精度。解决此问题的一种方法是提取整块区域的断裂线,并在随后的采样线滤波中保留断裂点,正确处理陡坡区域。加入断裂线约束前后滤波结果见表3,可以看出4类地形的总误差均减小,其中地形2的总误差减少近10%;Ⅰ型误差也大幅减小,这表明在加入断裂线约束后,在陡峭地表边缘的地面点会被正确保留;与此同时紧贴地表的地物点会被视为地面点而被错误保留,因此Ⅱ型误差略有增大。

表3 添加断裂线约束前后的误差比较

传统Snake模型利用内部能量控制曲线的平滑性和连续性,利用外部能量限制Snake的演变,而初始轮廓线距断裂线较远时,外部能量的作用不明显,故本文加入曲率约束和距离约束构成能量约束项,用于辅助外部能量对Snake演变的限制。Snake模型改进前后的滤波结果见表4,可以看出,采用传统的Snake模型提取断裂线能保证滤波精度在90%以上,这表明本文基于断裂线约束的滤波方法颇具成效。在Snake模型中加入约束能量项后,地形2的滤波精度能达到97.18%,其余地形的总误差也控制在5%左右,表明改进Snake模型的有效性。

采用传统二次曲面拟合方法对测区进行整体拟合,不能较大程度逼近真实地面,故本文按测区点云的x、y坐标划分格网,进行格网局部拟合,且各格网的边界区域相交重叠,提高了滤波精度。二次曲面拟合方法改进前后的滤波结果见表5。由表5可知,除地形3的Ⅱ型误差在改进二次曲面拟合后略微增加2%外,其余地形的各类误差均不同程度降低,其中地形4的滤波效果改善最为明显,滤波精度提高了8.64%,表明改进二次曲面拟合方法能真实展现地面形态,保证了地表的真实性和连续性。

表5 二次曲面拟合改进前后的误差比较

2.2 成都山区试验数据

数据来源于成都市山区,区域面积约28.9万m2,共约157.5万点,海拔为408.55～525.17 m,包括陡坡、茂密植被、梯田、河流和低矮建筑,如图16所示。图17显示了本文方法的滤波结果,其中扫描线间距为12 m、坡度阈值为20°、格网边长为40 m、计算曲率时的邻域范围为6 m。手动检查表明,本文方法在复杂区域表现良好,只有少数错误点。由于分割不足和树木遮挡,地面点被错分为植被点从而造成少量Ⅰ型误差,如图17的矩形A、B、C所示。此外,一些附着地面的低矮植被点被错误提取,如图17的矩形D、E所示;加入断裂线约束后,断裂线处紧贴地表的地物点被视为地面点而错误保留,因此Ⅱ型误差略有增大,如图17的矩形F所示。

图16 成都山区机载点云数据Fig.16 Point clouds from the city of Chengdu

图17 成都山区点云滤波结果Fig.17 The filtering results of the Chengdu data

比较本文方法与商业软件Terrasolid TerraScan中布料滤波方法,CC(CloudCompare)中布料滤波、基于坡度滤波方法的各类误差和Kappa系数(表6),可知本文方法的滤波性能更优、表现更好。

表6 本文方法、商业软件TerraScan和CloudCompare的误差比较

图18显示了研究区域的断裂线及基于断裂线约束生成的DEM。由图18中成都山区的DEM可知,低矮建筑、茂密植被等被滤除,山区地面的细节特征也得以保留。为探究断裂线对DEM精度的影响,选取图18中的3个区域(A1、A2、A3)分析有、无断裂线约束的DEM与参考地面之间的差异。

图18 断裂线提取与DEM生成Fig.18 Extracting brealines and generating DEMs with breaklines

图19比较了图18中A1、A2、A3区域有无断裂线约束DEM与参考地面点的二维横截面,可以看出,具有断裂线的DEM横截面更接近参考地面点。此外,由图19(a)、图19(b)可知,在无断裂线约束时,地形起伏大的地面会被错误滤除, 说明本文方法在Ⅰ类断裂线和Ⅱ类断裂线处都能准确地还原地表。在图19(c)这类地形起伏较小的Ⅱ类断裂线处,有无断裂线约束都能生成高质量DEM,而基于断裂线约束生成的DEM横截面更贴合参考地面。

图19 3个区域的有、无断裂线约束的DEM二维横截面比较Fig.19 2D cross-sections comparisons between the DEMs with/without breaklines and the ground truths in three areas

为进一步探讨断裂线约束的重要性,本文分别计算图18中A1、A2、A3区域参考DEM与有无断裂线约束DEM之间的高差,定量对比有无断裂线约束下的滤波后DEM精度,结果如图20所示。误差主要分布在断裂线周围,且对于Ⅰ类断裂线和Ⅱ类断裂线所在区域,无断裂线约束的DEM与参考DEM的高差更大,表明基于断裂线约束生成的DEM质量更高,在地形陡峭处表现更好。

3 结论与展望

基于地形表面崎岖不平、存在高程突变的特点,本文提出了一种基于断裂线约束和改进二次曲面拟合的点云滤波方法。不同于将单一滤波问题转换为一维采样线问题,本文改进Snake能量模型,并通过能量函数最小化提取断裂线特征,在采样线滤波基础上融合断裂线约束,阻止了高程突变处的大量地面点被错误滤除;根据采样线上的地面点利用改进二次曲面拟合过滤非地面点,较大程度还原地貌,提高了点云滤波精度,增加了点云滤波模型适用的场景类型。在4份ISPRS标准数据中,滤波的准确率分别达到94.74%、97.18%、94.96%和95.26%;在成都山区数据中,滤波精度达到97.61%。试验结果表明,本文方法能从机载激光点云中精确提取地面点,不仅适用于平坦的地形表面,而且对复杂崎岖地面也能取得较好的结果。目前本文滤波模型的参数设置还需要一定的先验知识,因此未来的一个工作方向是改进模型自适应参数选择,进一步提高模型自动化程度。