一种改进的线性预测滤波算法

张 靖,张晓君,江万寿,王建超,郭大海

(1.武汉大学测绘遥感信息工程国家重点实验室,武汉 430079;2.中国国土资源航空物探遥感中心,北京 100083)

一种改进的线性预测滤波算法

张 靖1,张晓君1,江万寿1,王建超2,郭大海2

(1.武汉大学测绘遥感信息工程国家重点实验室,武汉 430079;2.中国国土资源航空物探遥感中心,北京 100083)

对传统的线性预测算法进行改进,通过粗差剔除、初始点选取、特殊地形分析及光滑度检验等步骤提高了滤波精度。利用 ISPRS发布的 6组标准测试数据与原方法进行实验比较,结果证明,改进后的算法提高了滤波总精度以及地面点分类精度。

L iDAR;滤波;线性预测;粗差;分块

0 引言

机载激光雷达 (LightDetection And Ranging,Li-DAR)是一种快速获取高精度三维信息的新型遥感技术。LiDAR点云的滤波就是要从原始 LiDAR数据中滤掉非地面点,保留有效的地形信息。目前国际上已提出了多种滤波算法[1],如 Linderberger[2]、Killian等[3]和 Zhang等[4]提出的数学形态学方法,Axelsson[5]提出的渐进式 TIN滤波算法以及 Kraus等[6,7]提出的线性预测滤波算法等。

现有滤波方法的基本原理都是基于局部地形的连续性变化这一假设,算法设计的关键在于如何快速准确地检测出局部的高程突变,这需要解决两个关键问题:①如何确定激光点之间的邻域关系;②选择什么标准作为高程变化的判别条件。由于线性预测算法不需要确定激光点之间的严格位置关系,也不需要人为设定精确阈值,算法的适应性好,处理精度较高,因此被广泛地采用[8-10],并成为一种标准的L iDAR滤波方法。

但是,传统的线性预测算法也存在着一些缺陷:①分块大小固定,未考虑分块方式的科学性、地形形态的适应性及相邻分块的光滑程度;②算法易受低点误差影响。假设高程较低的点为地面点,当数据中存在低于地面的奇异点时,会导致计算结果出错;③算法处理过程中需进行多次迭代。迭代过程中用整个数据集拟合趋势面,计算效率较低;④在森林和城镇地区地面点较少时,很难估计出准确的地面[8]。

针对上述问题,本文提出了一种改进的线性预测滤波方法。具体改进如下:①对数据进行分块时,依照点密度的大小设置不同的分块大小,并通过高程信息的丰富程度划分不同的层次;②在滤波前通过对分块内直方图进行分析,去除奇异点;③初次拟合时,选取高程分布在 [15%,65%]范围内的点作初始值,再逐渐加点,以减少拟合的计算量;④通过分析特殊地形的高程直方图,对仅含有房屋点和少量地面点的分块进行拟合初值设定,提高了算法的适用性。实验证明,改进的线性滤波算法与原算法比较,具有更强的适应性和更高的滤波精度。

1 算法描述

1.1 传统线性预测算法

线性预测滤波是在滤波的同时内插DEM。对于给定的一系列离散点集,为了更好地进行线性估计,首先将原始数据划分为小块,用一个移动曲面来逼近局部区域的趋势面,一般设为二次曲面,其解析式为

式中,i为分块号;a1i等为第 i块曲面的参数;X,Y为平面坐标。

然后,用原始数据中每一点的高程值减去这点趋势面的高程值,即为拟合的残差。利用拟合残差确定该点在下一次曲面拟合中的权重,权函数定义为

式中,vi为第 i点的拟合残差;ρi为第 i点的权值;g为偏移值;w为地形最大起伏;a,b决定了权函数的陡峭程度,其计算方法为

式中,h为权值函数半宽值;s为半宽处的斜率;h和 s的含义如图 1所示。在一般情况下,可以取a=1,b=4。

图 1 权函数示意图Fig.1 Diagram of the weight function

参数 g可以根据残差直方图统计得到,主要有3种方法:①用预期的地面点精度计算,将残差直方图不断向左移,直至计算出的地面点精度与预期精度相同时,向左移动的距离即为 g(但当存在大的粗差时,该方法就会失效);②计算出向左移动每一个位置后点的标准差,找出最小的标准差,对应的向左移动的距离为 g(图 2);③通过估计的地面点的比例来计算。如果初始估计地面点比例为 40%,残差直方图要移动到有 20%地面点的地方[6]。

图 2 移动距离 g的计算方法示意图Fig.2 Calculate moving distanceg

1.2 改进后的算法流程

图 3为改进的线性预测滤波的算法流程图,其中实线框标明的部分表示传统的线性预测滤波流程,虚线框部分为改进的部分。

图 3 改进的线性滤波流程Fig.3 Flow chart of progressive linear prediction filtering

1.3 改进的关键技术

针对上述线性预测滤波算法的几个缺点,本文做出如下改进:

(1)粗差剔除。一般的滤波算法都认为在局部范围内最低点为地面点,而机载 L iDAR在飞行过程中都会出现粗差点 (即极低点或极高点),许多滤波算法都没有去除极低点的能力。本文从计算高程直方图的角度简单去除粗差。无论是极低点还是极高点,都表现为其高程与其他点的高程差异较大,且个数较少。在计算分块内高程直方图时,若把份数划分得较多,粗差点与正常点之间就会存在大量的空白频数,因此可通过去除空白频数之前的点,达到去除低点误差的目的。

(2)合理分块。线性预测滤波是用局部拟合的方法逐渐逼近地面的,因此分块不可太大;但为了更好地滤除建筑物,又要求分块不可太小。本文依据经验估算出最大建筑物的尺寸,一般设置为40m×40m,并根据点密度的不同自动调整分块大小。若窗口内点云的密度超过测区的平均点云密度,则将窗口分为 4个小窗口。这样,在保证分块大小相似的同时,确保了每个分块内有近似数量的激光点,拟合时可以得到更稳定的结果。

(3)初始点选取。Kraus等[7]提出的迭代线性最小二乘内插法,其计算始终是面向整个数据集,运算速度较慢。后来,有人提出首先计算出整个区域内所有点的高程平均值 Zi,设一个水平面 Z=Zi和一个高程阈值 a,将位于 Z=Zi-a与 Z=Zi+a之间的点作为初始点进行拟合,再加入符合条件的点,不断逼近 (其中 a的域值选取多为分块内最大高程值与最小高程值差值的 1/3～1/4)。这种方法减少了初始拟合点数,提高了运算速度,但在有些情况下这一范围的点太少,会影响平面拟合的精度。基于此,本文利用高程直方图统计的方法来选择初始点。一般选取高程数目分布在[15%,65%]之间的点数进行初步拟合,这样既可以提高运算速度,也可以选到足够数量的可靠点 (如图 4所示,可以选取高程范围在[334,356]范围内的点作为初始点)。

图 4 分块内高程直方图Fig.4 Histogram of elevation

(4)特殊地形分析。线性迭代最小二乘内插法是通过逐步拟合平面来逼近地面的,但只有分块内地面点相对较多时,拟合出来的平面才能代表正确的地面。有时房屋面积远大于分块面积,由于分块内地面点数太少,拟合出的最终平面是介于地面点和房屋顶面中间的平面,使得该方法失效。但从高程直方图的角度来看,分块的选择一般需满足以下条件:①高程有两个峰值;②峰值间高程相差较大;③峰值相邻的高程信息较为简单;④第一个峰值的点数较少。因此可以先判断分块是否满足特殊地形条件,若满足则选取第一个峰值附近的点为地面初始点。这样拟合出来的平面为真正地面,可以得到较为准确的滤波结果。

(5)光滑度检验。最小二乘内插法利用局部拟合平面的方法来逐步滤出地面点。由于地表被划分为一个个的小块,因此不仅小块内部要满足光滑、连续的条件,分块间也要满足光滑的条件。本文先对每一小块进行拟合滤波,得到该块的地面表达式后(如图 4中的阴影块),检验它与 4个邻域小块 (图 4中的 1,2,3,4)的光滑程度,可以按照不同的要求选择 0阶光滑 (即 f(x,y)=g(x,y))或一阶光滑 (即f′(x,y)=g′(x,y))等判断标准。若不满足光滑条件,可联合阴影部分与该块共同拟合平面,以提高滤波准确度。

图 5 相邻分块示意图Fig.5 Ad jacen t b lock

2 实验与讨论

2.1 试验数据及滤波结果

本文选取国际摄影测量与遥感学会 (ISPRS)在线发布的参考数据作为试验数据。对每个参考数据都事先做好了人工分类,标上地面点或者非地面点,以利于检验滤波的正确性。共选取了 6套测试数据,这些数据中既有城区地形也有农村地形,其中测区 1(Sample 24)的典型地物为陡坎、人字型房屋和植被,测区 2(Sample 31)的典型地物为一个方形复杂的建筑物,测区 3(Sample 12)的典型地物为城区密集建筑,测区 4(Sample 41)的典型地物为形状不规则层次多的建筑,测区 5(Sample 51)的典型地物为斜坡及斜坡上的植物,测区 6(Sample 54)为密集低矮植物。图 6显示了测区 3和测区 6的原始数据及其滤波结果。

图 6 测试数据与滤波结果Fig.6 Raw data and filtering results

2.2 滤波结果分析

从图 6可以看出,该算法能够较好地滤除密集的规则房屋 (测区 3)等人工建筑物,对于平地上低矮的植物和斜坡 (测区 6)也有较好的效果。但是,该算法对于地表断裂处 (测区 1)存在一定的误差。

表 1列出 6个测区的原方法与改进方法的滤波结果,包括每个测区中标准数据和两种滤波方法滤出的正确地面点及地物点点数。

表 1 各测区滤波结果比较Tab.1 Comparasion of filtering results (个)

从表 1可以看出,线性迭代最小二乘滤波可以较好地滤出地物点,这样可以保证得到的DEM有较高的精度;但有时原算法的地面点分类正确率相对较低 (如测区 4),而改进的方法在尽量保证地物点滤波正确性的同时,提高了地面点的正确率。

图 7示出 6个测区的误差统计图。假定 a表示数据中地面点的点数,b表示数据中地物点点数,c表示分类正确的地面点点数,d表示分类正确的地物点点数,则正确分类的地面点比例为 100×%,正确分类的地物点的比例为 100×%,总分类正确点的比例为 100×%。

图 7 滤波结果正确率比较Fig.7 Comparison of filter ing accuracy

从图 7(a)可以看出,原滤波算法的总点数精度在 90%左右,但波动较大,说明该方法对于不同的地形有着一定的依赖性。例如对于含有较密集建筑的第 3测区,滤波结果就不如滤除植被的第 6测区好。而改进的方法,由于提前对特殊地形进行了分析,在处理房屋时就提高了滤波精度。图 7(b)改进的方法比原方法的地面正确率更高。从单个测区看 ,第 2、4、6测区的正确率高于第 1、3、5测区 ,说明该方法在陡坎、斜坡密集的地方,滤出的地面点精度会降低。图 7(c)除了第 1测区外,其他测区精度都较高。由于第 1测区中有一个陡坎,陡坎有高程突变,滤波时会误检测为地物点,造成一定误差。

总的来看,无论是城区的复杂建筑、密集房屋,还是农村的斜坡、低矮植被,改进的滤波算法都有较好的滤波结果。从地面点和地物点分类正确率来看,改进的算法在尽量保证地物分类正确的同时,提高了地面点分类的正确率。

3 结论

(1)在总结各种滤波算法基础上,通过对高程直方图的分析,在粗差剔除、初始点选取、特殊地形分析方面对线性迭代最小二乘滤波方法进行了改进,并在滤波结束后检验了相邻分块的光滑程度。

(2)采用总点数分类正确率、地面点正确率和地物点正确率 3个评价指标对改进的方法和原线性迭代最小二乘滤波结果进行了比较。结果表明,改进的方法在处理密集房屋等方面提高了滤波精度。

[1] Vosselman S G.Experimental Comparison of Filter Algorithms for Bare-Earth Extraction from Airborne Laser Scanning Point Clouds[J].ISPRS Journal of Photogrammetry and Remote Sensing,2004,59(1/2):85-101.

[2] Lindenberger J.Laser-profilmenssungen zur Topographischen Gelandeaufnahm e[D].Stuttgart:Universitat Stuttgart,Verleg der Bayerischen Akademie der Wissenschaften,1993.

[3] Kilian J,Haala N,EnglichM.Cap ture and Evaluation of Airborne Laser ScannerData[J].IAPRS,1996(31-B3):383.

[4] Zhang K Q,Chen S,Whitman D,et al.Progressive Morphological Filter for Removing Nonground Measurements from Airborne L i-DAR Data[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(4):872-882.

[5] Axelsson P.DEM Generation from Laser Scanner Data Using Adaptive TIN Models[J].International Archives of the Photogrammetry,Remote Sensing and Spatial Information Sciences,2000,33(B4/1):110-117.

[6] Kraus K,Pfeifer N.Determination of Terrain Models in Wooded Areas with Airborne Laser Scanner Data[J].ISPRS Journal of Photogrammetry and Remote Sensing,1998,54(4):193-203.

[7] Kraus K,Pfeifer N.Advance DTM Generation from LiDAR Data[J].International Archives of Photogrammetry and Remote Sensing,2001,34(3/W 4):23-35.

[8] 刘经南,许晓东,张小红,等.机载激光扫描测高数据分层迭代选权滤波方法及其质量评价[J].武汉大学学报 (信息科学版),2008,33(6):551-555.

[9] 王 刃,徐 青,朱新慧,等.用分层稳健线性估计法从机载 L i-DAR数据中获取DEM[J].测绘科学技术学报,2008,25(3):188-191.

[10]曾齐红,毛建华,李先华,等.机载激光雷达点云的阶层式分类[J].测绘科学,2008,33(1):103-105.

(责任编辑:刘心季)

Progressive Linear Prediction Fitting for Extracting DTM from Airborne LiDAR Data

ZHANG Jing1,ZHANG Xiao-jun1,JIANGW an-shou1,WANG Jian-chao2,GUO Da-hai2
(1.LIESMARS,Wuhan University,Wuhan 430079,China;2.China Aero Geophysical Survey&Remote Sensing Center for Land and Resousces,Beijing 100083,China)

A progressive linear prediction filtering algorithm is proposed for extracting DEM from LiDAR data.Some processes are inserted in the ordinary linear prediction algorithm,such as gross error detection,initial value selection,land form analysis,and smoothness detection.The authors used this algorithm to process 6 datasets published by ISPRS as standard filtering test data.The results show that the improvement in the traditional methods can increase the precision of DEM.

LiDAR;Filtering;Linear prediction;Outlier;Partitioning

江万寿 (1967-),男,研究员,博士生导师。主要研究方向为数字摄影测量及遥感数据处理。联系电话:(+86)27-68778092-8321(O),E-mail:jw s@lm ars.whu.edu.cn。

TP 75

A

1001-070X(2011)01-0052-05

2010-04-21;

2010-06-10

国家“863”计划项目 (编号:2006AA 06A 208)和国家自然科学基金项目 (编号:40671159)共同资助。

张 靖 (1982-),男,博士研究生,主要从事 LiDAR数据处理方面的研究。