融合点、对象、关键点等3种基元的点云滤波方法

㊀２０１６年１１月A c t aG e o d a e t i c ae tC a r t o g r a p h i c aS i n i c a N o v e m b e r,２０１６㊀㊀第４５卷㊀第１１期测㊀绘㊀学㊀报V o l．４５,N o．１１引文格式:林祥国,张继贤,宁晓刚,等．融合点㊁对象㊁关键点等３种基元的点云滤波方法[J]．测绘学报,２０１６,４５(１１):１３０８Ｇ１３１７．D O I:１０．１１９４７/j．A G C S．２０１６．２０１６０３７２．L I N X i a n g g u o,Z HA N GJ i x i a n,N I N G X i a o g a n g,e ta l．F i l t e r i n g o fP o i n tC l o u d sU s i n g F u s i o no fT h r e eT y p e so fP r i m i t i v e sI n c l u d i n g P o i n t s,O b j e c t sa n d K e y P o i n t s[J]．A c t aG e o d a e t i c ae tC a r t o g r a p h i c aS i n i c a,２０１６,４５(１１):１３０８Ｇ１３１７．D O I:１０．１１９４７/j．A G C S．２０１６．２０１６０３７２．融合点㊁对象㊁关键点等３种基元的点云滤波方法林祥国,张继贤,宁晓刚,段敏燕,臧㊀艺中国测绘科学研究院,北京１００８３０F i l t e r i n g o fP o i n tC l o u d sU s i n g F u s i o no fT h r e eT y p e so fP r i m i t i v e s I n c l u d i n g P o i n t s,O b j e c t sa n dK e y P o i n t sL I NX i a n g g u o,Z H A N GJ i x i a n,N I N GX i a o g a n g,D U A N M i n y a n,Z A N GY iC h i n e s eA c a d e m y o f S u r v e y i n g a n dM a p p i n g,B e i j i n g１００８３０,C h i n aA b s t r a c t:P r i m i t i v e,b e i n g t h eb a s i c p r o c e s s i n g u n i t,i s o n eo f t h e k e y f a c t o r s t od e t e r m i n e t h ea c c u r a c y a n de f f i c i e n c y o f p o i n t c l o u d f i l t e r i n g．T r i a n g u l a r i r r e g u l a r n e t w o r k(T I N)p r o g r e s s i v ed e n s i f i c a t i o n(T P D) a n do b j e c tＧb a s e dT I N p r o g r e s s i v ed e n s i f i c a t i o n(O T P D)a r et w oe x i s t i n g f i l t e r i n g m e t h o d s,b u ts i n g l e p r i m i t i v e i s e m p l o y e db y t h e m．Am u l t i p l eＧp r i m i t i v e sＧb a s e dT I N p r o g r e s s i v ed e n s i f i c a t i o n(M P T P D)f i l t e r i n g m e t h o d i s p r o p o s e d．I t i s c o m p o s e do f t h r e e k e y s t a g e s,i n c l u d i n gp o i n t c l o u d s e g m e n t a t i o n,e x t r a c t i o no f k e y p o i n t s o f o b j e c t s,t h e k e yＧp o i n t sＧb a s e d j u d g i n g o f t h eo b j e c t s．S p e c i f i c a l l y,p o i n t,o b j e c t a n d t h ek e y p o i n t sa r et h e p r i m i t i v eo ft h e a b o v et h r e es t a g e sr e s p e c t i v e l y．F o u rt e s t i n g d a t a s e t s,i n c l u d i n g t w o a i r b o r n e L i D A Ra n d t w o p h o t o g r a m m e t r i c p o i n t c l o u d s,a r eu s e d t ov e r i f y t h eo v e r a l l p e r f o r m a n c e so f t h e a b o v e t h r e e f i l t e r i n g m e t h o d s．E x p e r i m e n t a l r e s u l t s s u g g e s t t h a t t h e p r o p o s e dM P T P Dh a s t h eb e s t o v e r a l l p e r f o r m a n c e．I nt h ev i e w p o i n to fa c c u r a c y,M P T P D a n d O T P D h a v et h es i m i l a ra c c u r a c y．M o r e o v e r, c o m p a r e dw i t h t h eT P D,M P T P D i sa b l e t o r e d u c eo m i s s i o ne r r o r sa n d t o t a l e r r o r sb y２２．０７％a n d８．４４％r e s p e c t i v e l y．I n t h ev i e w p o i n to fe f f i c i e n c y,u n d e rm o s to f t h ec a s e s,T P Di s t h eh i g h e s t,M P T P Di s t h e s e c o n d,a n dO T P D i s t h es l o w e s t．M o r e o v e r,t h e t o t a l t i m ec o s t o fM P T P D i so n l y５７．９３％o f t h eo n eo f O T P D．K e y w o r d s:f i l t e r i n g;L i D A R p o i n t c l o u d;p h o t o g r a m m e t r i c p o i n t c l o u d;o b j e c t s;t r i a n g u l a r i r r e g u l a r n e t w o r k F o u n d a t i o n s u p p o r t:T h eN a t i o n a l N a t u r a l S c i e n c e F o u n d a t i o n s o f C h i n a(N o．４１３７１４０５);T h e F o u n d a t i o n f o r R e m o t eS e n s i n g Y o u n g T a l e n t sb y t h e N a t i o n a l R e m o t e S e n s i n g C e n t e r o f C h i n a;T h e B a s i c R e s e a r c h F u n d o f t h eC h i n e s eA c a d e m y o f S u r v e y i n g a n dM a p p i n g(N o．７７７１６１１０３)摘㊀要:基元是影响点云滤波精度和效率的关键因素之一.本文提出了一种基于多基元的三角网渐进加密(M P T P D)滤波方法.它包括点云分割㊁对象关键点提取㊁基于关键点的对象类别判别３个主要阶段,且３个阶段的基元分别为点㊁对象㊁关键点.使用了４景机载激光雷达和摄影测量点云数据对M P T P D㊁三角网渐进加密(T P D)㊁基于对象的三角网渐进加密(O T P D)３种滤波方法进行了性能测试.试验表明,M P T P D方法具有整体上最优的性能:在精度方面,M P T P D与O T P D两种方法的精度相当, M P T P D方法的一类误差I㊁总误差T比T P D的相应误差分别低约２２．０７％和８．４４％;在效率方面,多数情况下T P D㊁M P T P D㊁O T P D方法的效率依次降低,且M P T P D的平均耗时是O T P D平均耗时的５７．９３％.关键词:滤波;激光雷达点云;摄影测量点云;对象;三角网中图分类号:P２３７㊀㊀㊀㊀文献标识码:A㊀㊀㊀㊀文章编号:１００１Ｇ１５９５(２０１６)１１Ｇ１３０８Ｇ１０基金项目:国家自然科学基金(４１３７１４０５);遥感青年科技人才创新资助计划;中国测绘科学研究院基本第１１期林祥国,等:融合点㊁对象㊁关键点等３种基元的点云滤波方法科研业务费(７７７１６１１０３)㊀㊀随着激光雷达(l i g h td e t e c t i o na n dr a n g i n g, L i D A R)[１]测量㊁多视影像密集匹配[２]技术的完善和行业应用的深入,点云滤波的重要性日益突出.本文的点云涉及机载L i D A R点云和航空㊁航天多视立体影像密集匹配的点云等３种类型.在点云处理和信息提取领域,滤波是指区分点云中的地面点和非地面点的过程[１,３],它是生成数字高程模型(D E M)㊁分类㊁目标识别和三维重建的基础和必经的步骤[４].文献[１,４Ｇ５]对目前众多点云滤波方法进行了系统的介绍.其中,有代表性的方法有三角网(t r i a n g u l a r i r r e g u l a r n e t w o r k,T I N)渐进加密(T I N p r o g r e s s i v e d e n s i f i c a t i o n,T P D)[６Ｇ７]㊁分层稳健线性内插[８]㊁坡度滤波[９]㊁数学形态学滤波[１０]㊁基于聚类/对象的滤波[１１Ｇ１２]等.已有方法中涉及的基元(基本处理单元)有点[６Ｇ７]㊁对象[１１Ｇ１２]㊁体素[１３]或剖面[１４]等多种类型;且后３种基元具有一定的共性,本质上是点基元的一种集合和再组织方式,本文仅关注其中的对象.由于点易受粗差㊁地形断裂的负面影响,而对象比点更能增强点云处理效果[１５],因此基于对象的点云滤波方法[１６Ｇ２１]是研究的一个热点.然而,与基于点的滤波方法相比,尽管基于对象的滤波方法可以在一定程度上提高滤波精度,但是也存在效率低下的问题[２１].文献[２２]提出基于多实体的点云分类方法,在分类的不同阶段使用不同的实体以实现更优的分类效果.借鉴上述策略,本文设计一种既能继承基于对象方法的优势㊁又不显著降低基于点方法的效率的滤波技术,即同时提高基于点的T P D 方法[６]精度和基于对象的三角网渐进加密(o b j e c tＧb a s e dT P D,O T P D)方法[２１]效率,本文称之为基于多基元的三角网渐进加密(m u l t i p l eＧp r i m i t i v e sＧb a s e dT P D,M P T P D)方法.它有３个创新点:①使用多基元㊁而非单一的基元参与运算,其中多基元包括点㊁对象㊁关键点等３种类型,且在不同阶段使用不同类型的基元;②使用关键点代替对象参与判别,即在核心判别步骤中,使用对象的关键点替代对象进行运算以提高效率;③提出一种简单㊁快捷的关键点检测算法.特别指出,本文的一个 对象 指 点云分割后具有同一标号的点集 , 关键点 又是对象点集的一个子集,即关键点本质上仍然是原始点云中的点,而非额外创造的,但是关键点具有特殊性.另外,处理一个 对象 ,可以通过处理该对象包含的点集来实现,也可以通过处理 关键点 来实现.１㊀本文的滤波方法M P T P D方法包括基于表面生长的点云分割㊁对象关键点提取㊁基于关键点的对象类别判别等３个主要步骤.整体技术框架如图１所示;图２展示了某一点云各个处理步骤的效果,文中数字 １ 代表 非地面点类 ,数字 ２ 代表 地面点类.图１㊀本文方法的整体技术流程图F i g．１㊀T h ew h o l ew o r k f l o wo f t h e p r o p o s e dm e t h o d １．１㊀基于表面生长的点云分割点云分割是对点云数据中每个点按照一定的判别规则进行标号的过程.分割后,满足同一规则的点集被赋予同一标号,且每一点集称为一个对象.本文判别的是３D空间中邻近且共平面的点.另外,不满足上述判别规则的孤立点㊁邻近点数目不足的点㊁共平面性差的点亦会被标号.鉴于表面生长[２３]算法具有所需参数少㊁分割效果好㊁普适性好的特点,本文使用它对点云进行分割,其主要步骤如下.第１步,估计法向量和残差,处理过程如下.(１)加载点云数据,并将所有点的类别标记为 １ ㊁标号状态记为 未分割 ,设共计有n 个点.(２)建立点云的三维k dＧt r e e[２４]空间索引.(３)逐一处理每一个点,即对第i(i＝０,２, ,n－１)个点,首先利用k dＧt r e e求取其k个最临近点,然后利用特征值法[２５]求当前点及k个邻近点构成点集的拟合平面方程,即可确定第i个９０３１N o v e m b e r ２０１６V o l ．４５N o ．１１A G C Sh t t p :ʊx b ．s i n o m a ps ．c o m 点的法向量φi 及其残差λi .第２步,进行区域生长.图２展示了表面生长过程.该过程需要设定两个参数,包括法向量间角度差异阈值α㊁距离阈值r .其中,α的使用规则为:对于当前要扩张的种子点,利用k d Ｇt r e e 查询其某一邻接点作为待处理的点,设当前种子点的法向量为αi ,待处理的邻接点法向量为αj ,且αi 和αj 均为单位法向量;如果待处理的邻接点与种子点之间法向量的点积αiαj 小于阈值c o s α,则将该待处理的邻接点加入到种子点队列中,并作为新的种子点,继续判别其他的未处理点;反之,则跳过该待处理的邻接点,处理其他的未处理点.r 的使用规则为:如果待处理的邻接点到种子点确定的种子平面的距离小于给定距离阈值r ,则将该待处理的邻接点加入到种子点队列中,并作为新的种子点,继续判别其他的未处理点;反之,则跳过该待处理的邻接点,处理其他的未处理点.具体处理如下.图２㊀表面生长过程的示意图F i g ．２㊀I l l u s t r a t i o no f t h e s u r f a c e g r o w i n gpr o c e s s (１)记 区域标记号 从０开始.(２)检查 未分割 点集中点的数量,如果数量为０,则转到步骤(６);否则,接着从 未被分割 的点集中,寻找出残差λ最小的点,以该点为种子点并将该点压入一个种子点的队列,且将该点的处理状态标记为 未处理 ,开始进行区域生长.(３)取种子点队列中第一个 未处理 的种子点,利用k d Ｇt r e e 求取该种子点的k 个最临近点.(４)逐一对于每一个临近点进行判别.如果临近点已经被赋予分割号,则不予以处理;另外,若临近点已经在种子点的队列中,则不予以处理;接着,分别按照法向量间角度差异和距离差异的规则进行当前种子点和该邻近点的相似性判别.如果该邻近点同时满足两个相似性的条件,则将该邻近点加入到种子点队列中;反之,如果该邻近点没有同时满足两个相似性的条件,则不予以处理.邻近点判别完毕后,将该种子点的处理状态标记为 已处理.(５)顺序检查种子点队列中是否有 未处理 的点.如果有,返回步骤(３);否则,将种子点队列中的点集的标号记为 区域标记号 ,状态记为 已分割 ,同时 区域标记号 自增１,清空种子点队列后返回步骤(２).(６)结束.经过上述分割后,任意一点被划到一个对象,但是部分对象的点的数量较少.个别情况下,一个对象仅仅包含一个点.图３(a)展示了某机载L i D A R 点云数据,图３(b)展示了其分割效果,其中地面点被聚为若干个对象,多数地面对象包含点的数量较多;一个建筑物可以被聚为一个或者若干个对象,这与其类型㊁点云密度㊁精度等多个因素有关;孤立的植被点㊁粗差点也往往被分割为一个对象.图３㊀本文提出滤波方法的关键步骤处理效果示意图F i g ．３㊀I l l u s t r a t i o no ft h e p r o c e s so ft h e p r o po s e d f i l t e r i n g me t h o d １．２㊀对象关键点的提取本文的关键点包括外轮廓点㊁内特征点㊁最高点和最低点.第１步,计算每个对象包含点的数量.如果数量不大于经验阈值４,则将该对象的点集记为对象的关键点.反之,进入第２步.第２步,识别每个对象的关键点.图４展示了某对象关键点检测的主要过程.提取的基本原理是仅利用某一对象点集的水平坐标信息生成T I N .该T I N 中,处于边缘的三角形仅有两个三角形通过边相邻;而处于非边缘的三角形有３个０１３１第１１期林祥国,等:融合点㊁对象㊁关键点等３种基元的点云滤波方法三角形通过边相邻.因此,可以通过一个三角形的一边为邻边的邻接三角形的数量来判断该三角形是否处于边缘.处于边缘的三角形涉及的３个顶点记为 外轮廓点 ,如图４(b )所示.接着,删除该T I N 中的短边,并以删除短边后的T I N 为索引进行连通区域分析(c o n n e c t e dc o m p o n e n t a n a l ys i s ).则该对象被分割为若干子对象.若子对象包含的点的数量大于经验阈值４,则该子对象的点集被认为是 内特征点 ,如图４(c)所示.注意上述关键点无需有序排列,这与A l ph a s h a p e [２６]算法有着显著的差别.第２步具体内容如下.(１)建立某一对象的二维T I N .(２)检测 外轮廓点 .通过上述邻接三角形数量的规则识别 外轮廓点 ,如图４(b )所示.(３)删除短边.假设输入原始点云的平均点间距为g (单位:m ),且g 已知.删除T I N 中二维边缘长度小于经验阈值３g 的边.(４)检测 内特征点 .通过上述连通区域分析获取的子对象包含点的数量的规则识别 内特征点 ,如图４(c)所示.(５)检测最高点㊁最低点.另外,外轮廓点㊁内特征点㊁最高点㊁最低点不可重复.如有重复,则只保留其中一个.图４(a )展示了某一地面对象包含３２１９９８个点,其外轮廓点㊁内特征点和关键点分别如图４(b )㊁(c )和(d )所示,图４(d )只包含１９８７５个点.从图４(a )和(d )中的D E M 看,尽管构建D E M 的点的数量差别悬殊,但是D E M 的表达效果却趋于一致.对两个D E M ,不仅最高㊁最低值一致,且相应像素值之差的绝对值的平均值和标准差分别为０．０５和０．０１,这反映了提取的关键点既能显著地减少点的数量㊁又能逼近真实的对象原始形态.另外,图２(c )展示了图２(a )中点云的关键点检测结果,其中原始点云包含８２６４１６个点㊁而关键点只包含６４２５７个点,关键点数量只占原始点数量的７．７８％.图４㊀对象关键点检测的示意图F i g ．４㊀I l l u s t r a t i o no f d e t e c t i o no f t h ek e yp o i n t s o f a no b je c t １．３㊀基于关键点的对象类别判别经典T P D 方法的运算过程中,T I N 构建和点类别判别占整个滤波时间的比重很大[２７].本文采用关键点替代对象的目的是能同时显著地减少参与T I N 构建的㊁参与判别的点的数量以提高效率,又可使构建的T I N 尽可能地逼近区域真实的D E M 以确保精度.本节是一个迭代过程.共涉及４个参数:最大建筑物长度b (单位:m )㊁最大角度阈值θ(单位:ʎ)㊁最大距离阈值d (单位:m )和最大地形角度阈值t (单位:ʎ).具体过程如下.第１步:格网划分.求点云在X O Y 平面上的最小外包矩形,并在X O Y 平面上对该最小外包矩形进行格网划分㊁且格网的尺寸为b ˑb .１１３１N o v e m b e r２０１６V o l．４５N o．１１A G C S h t t p:ʊx b．s i n o m a p s．c o m第２步:将全部对象处理状态均标记为 未处理 .第３步:选择初始地面种子点.逐一选择每一个格网的地面种子点.即对每一个格网,找到格网中高程值最低的点所在的对象.如果对象的面积小于４．００m２,则继续找到高程值次低的点所在的对象直至找到面积大于４．００m２的对象.将该对象的关键点作为该格网的地面种子点,且该对象的类别被标记为 ２ ㊁处理状态被标记为 已处理 .第４步:构建初始地面种子点的T I N.该T I N代表该区域初始的D E M.第５步:迭代的判别对象类别.子步骤包括: (１)迭代次数记为０.(２)以对象为基本处理单元,逐一通过每一个 未处理 对象的关键点的判别,实现该对象类别的判别.对每一个 未处理 对象逐一判别其 未处理 关键点,找到该关键点落入的三角形,计算该关键点到三角形构成的平面的距离及该关键点到三角形３点的夹角,并找出３个夹角中的最大角.进行下述判别:如果同时满足距离小于d㊁最大夹角小于θ,则认为该关键点是地面点,将该关键点的类别号标记为 ２ ,处理状态标记为 已处理 ,继续处理下一个 未处理 关键点;否则,检查三角形构成的平面的倾角,进行下述判别:如果倾角小于,继续处理下一个 未处理 关键点;反之,将当前关键点以所在三角形的最高点为中心做一个镜像点(参考文献[６,１２]),且该镜像点的高等于该关键点.对该镜像点进行类似的判别.如果该镜像点被判别为 ２ ,则将该关键点的类别号标记为 ２ ,处理状态标记为 已处理 ;继续处理下一个 未处理 关键点.判别完毕,统计该对象的关键点的数量㊁关键点属于地面点的比例.如果该比例大于５０％,将该对象的处理状态标记为 已处理 ,类别为 ２ ;否则,将该对象的处理状态重新标记为 未处理 ,类别为 １ .(３)利用新识别的属于地面点的关键点更新T I N,同时迭代次数自增１.(４)重复上述步骤(１)至(３),继续执行直至迭代次数达到经验阈值５,或者没有新识别的地面点则停止迭代.２㊀试验与分析基于V i s u a l S t u d i o２０１０C＋＋集成开发环境实现了本文提出的M P T P D方法,同时对T P D[６]㊁O T P D[２１]两种方法进行性能比较,上述３种方法均采用串行计算,未采用并行计算技术.其中,T P D包括低位粗差点剔除和１．３节描述的５个主要步骤,其基元为点;O T P D包括１．１节和１．３节两个相似部分,其基元为对象.为了增加效率的可比性,尽管３种滤波方法的基元的不同,但相似步骤涉及的算法一致.试验平台的配置: T h i n k P a d W５２０笔记本,C P U为I n t e l酷睿i７Ｇ２７６０QM２．４G H z,内存２．９８G B,装配W i n d o w s X P系统.２．１㊀试验数据及结果本文共使用了４个场景的点云开展试验(图５),它们的基本信息见表１.前两个场景的点云为开放的机载L i D A R数据;后两个为摄影测量点云,其中,第３个场景的原始影像由T r i m b l e G e r m a n y G m b H公司免费提供㊁点云由中国测绘科学研究院的P i x e l G r i d软件生成,第４个场景的点云由德国宇航局免费提供.另外,试验数据１为国际摄影测量与遥感协会第三委员会提供的测试数据C S i t e１,该数据位于德国,如图５(a)所示;试验数据２由I E E E G R S LF u s i o nC o n t e s t２０１３提供,该数据位于美国休斯敦大学附近,本文截取了原始点云的一部分,如图５(b)所示;试验数据３对应的斜影像由天宝公司A O S系统获取,该数据位于德国柏林市和波茨坦市附近,本文截取了由影像生成的摄影测量点云的一部分,如图５(c)所示;试验数据４由G e o E y eＧ１的立体影像对生成,该数据位于德国的慕尼黑市,本文截取了该摄影测量点云的一部分,如图５(d)所示.对４个试验数据进行滤波时,使用的相关参数的值见表２.其中,T P D㊁O T P D和M P T P D ３种方法对１．３节的参数b㊁θ㊁d㊁t采用了相同的参数值,O T P D和M P T P D两种方法对１．１节的参数k㊁d㊁r采用了相同的参数值.鉴于篇幅的原因未展示３种方法的滤波结果.后续试验分析表明３种滤波方法均能正确地区分多数的地面点和非地面点,且O T P D和M P T P D两种方法的滤波效果相当㊁并优于T P D方法的滤波效果.２．２㊀精度评价本文采用文献[４]中的一类误差I(将地面点２１３１第１１期林祥国,等:融合点㊁对象㊁关键点等３种基元的点云滤波方法错分为非地面点的数量占地面点数量的比例)㊁二类误差I I(将非地面点错分为地面点的数量占非地面点数量的比例)和总误差T(错分点数量占全部点数量的比例)３个指标定量衡量滤波精度.同时,使用了人工半自动解译的方式识别了４个试验数据的地面点和非地面点,并将每个试验数据的人工识别结果作为真值计算滤波方法的误差.４个试验数据的３类误差值见表３.图５㊀４个试验数据F i g．５㊀T h e f o u r t e s t i n g d a t a s e t s㊀㊀表３中,T P D方法的Ⅰ㊁Ⅱ㊁T３类误差的平均值分别为２７．０９％㊁２．２７％和１２．５２％,O T P D方法相关误差的平均值分别为４．９６％㊁２．７５％和４．０４％,M P T P D方法相关误差的平均值分别为５．０３％㊁２．７８％和４．０８％.数字说明O T P D方法与M P T P D方法的各类误差均十分接近.且,一类误差I和总误差T均呈现T P D＞O T P DʈM P T P D的趋势,而二类误差I I呈现T P D＜O T P DʈM P T P D的趋势.另外,M P T P D和O T P D两种方法的一类误差I㊁总误差T的平均值比T P D的分别低约２２．１０％㊁８．４６％,而两种方法的二类误差I I的平均值比T P D的高约０．４９％.可见,与T P D相比,M P T P D和O T P D两者会轻微地提高二类误差I I,但会显著地降低一类误差I㊁总误差T.可知,O T P D和M P T P D两种方法的滤波精度相当㊁且高于T P D方法的滤波精度.表１㊀４个试验数据的基本信息T a b．１㊀T h e f u n d a m e n t a l i n f o r m a t i o no f t h e f o u r t e s t i n g d a t a s e t s试验数据点云类型点数/个长/mˑ宽/m平均点间距/m是否有粗差地形类型１机载L i D A R１３６６４０８１１００．７１ˑ７０１．４６０．６９有城区与森林２机载L i D A R２３３５０４８１０８７．４７ˑ８８９．５１０．６４有城区３航空影像匹配２０４４４４７６５５．１２ˑ６２４．３９０．４４有城区４卫星影像匹配２４９７６３８１０１１．８５ˑ１０３５．８５０．５０有城区表２㊀４个试验数据中３种滤波方法的相关参数取值T a b．２㊀T h e v a l u e s o f t h e i n p u t p a r a m e t e r s o f t h e t h r e e f i l t e r i n g m e t h o d s f o r t h e f o u r t e s t i n g d a t a s e t s方法T P D O T P D和M P T P D参数b/m t/(ʎ)θ/(ʎ)d/m k/点α/(ʎ)r/m 试验数据１２０８０．００６．００１．４０２０３０．０００．５０试验数据２㊁３㊁４６０８８．００６．００１．４０２０１０．０００．５０３１３１N o v e m b e r ２０１６V o l ．４５N o ．１１A G C Sh t t p :ʊx b ．s i n o m a ps ．c o m 表３㊀４个试验数据中３种滤波方法３类误差的统计值T a b ．３㊀T h ev a l e so ft h r e et y pe so fe r r o r so ft h et h r e ef i l t e r i ng m e th o d s f o r t h e f o u r t e s ti n g da t a s e t s 试验数据误差类型T P D /(％)O T P D /(％)M P T P D /(％)１Ⅰ３４．７４１１．３５１１．３６Ⅱ２．３９４．２１４．１６T１７．９２７．６４７．６２２Ⅰ２２．４６６．８１６．９８Ⅱ３．０２０．７９０．７７T１３．３５４．０４４．０７３Ⅰ２３．７１１．２９１．３３Ⅱ１．８０５．７４５．８１T９．４０４．２０４．２６４Ⅰ２７．４５０．３９０．４３Ⅱ１．８７０．２４０．３６T９．４００．２８０．３８２．３㊀效率评价本文采用时间花费来衡量滤波效率.为此,分别记录了３种方法对４个试验数据进行滤波的各个阶段耗时及总耗时,具体的统计数据见表４.其中,将滤波过程划分为５个阶段:低位粗差点剔除㊁基于表面生长的点云分割(１．１节)㊁对象关键点提取(１．２节)㊁基于关键点的对象类别判别(１．３节)前３步和后两步,在表４中分别称为阶段１㊁阶段２㊁阶段３㊁阶段４㊁阶段５.表４的统计数据表明,在４个试验数据中,有３个出现了３种滤波方法的总效率均呈现T P D＞M P T P D＞O T P D 的规律.以试验数据３为例,T P D ㊁M P T P D ㊁O T P D３种方法的总耗时由少到多依次为约２９９s ㊁４１５s ㊁８６６s .但在第１个试验数据中,M P T P D 的效率高于T P D ㊁O T P D .整体上,T P D 效率最高,M P T P D 次之,O T P D 最低.表４还表明T P D ㊁M P T P D ㊁O T P D３种方法在每个试验数据上的总耗时平均分别为２４８．７６s㊁３５５．２６s ㊁６１３．３０s .如果以效率最慢的O T P D 的基准,T P D 和M P T P D 的效率分别是O T P D 的２．４７倍㊁１．７３倍.而且,表４还表明每种滤波方法的各个阶段的耗时的比例也有着显著的差别,表现为:(１)T P D 方法中,阶段１㊁阶段４和阶段５这３个阶段的耗时占总耗时比例的平均值分别为８４．６１％㊁６．４４％㊁８．９５％,可见,阶段１和阶段５这两个阶段占了T P D 方法总耗时的绝大部分比例,其中粗差剔除的相关比例最大㊁且点云中的粗差越复杂相应的比例越大.例如,试验数据１的粗差多且多样,粗差剔除的时间花费占总时间花费的９５．１９％,是４个试验数据中比例最大的,这是在第１个试验中M P T P D 的效率高于T P D 的原因.(２)O T P D 方法中,阶段２㊁阶段４㊁阶段５等３个阶段的耗时占总耗时比例的平均值分别为４５．６２％㊁１０．１９％㊁４４．８７％.可见,O T P D 的３个阶段中,除了第２个阶段外,其他两个阶段均比较耗时㊁且占总耗时的比例相当.(３)M P T P D 方法中,阶段２㊁阶段３㊁阶段４㊁阶段５等４个阶段的耗时占总耗时比例的平均值分别为７５．４７％㊁２０．８６％㊁２．８７％㊁０．８０％.可见,M P T P D 的４个阶段中,前两个阶段的累计耗时占了总耗时的绝大部分,而后两个阶段的累计耗时占总耗时的比例极小.这是M P T P D 与O T P D尽管在阶段２的耗时相同㊁但是M P T P D 的总耗时显著的低于O T P D 的原因.另外,表４展示了一个很有兴趣的现象.３种滤波方法均有阶段４㊁阶段５两个阶段,但由于基元的不同导致两个阶段(尤其是阶段５)的效率有着显著差别.例如,第４个试验数据中,T P D 和O T P D 两种方法的阶段４的耗时分别为９．４５s㊁１０９．７０s ,但M P T P D 的相应耗时仅为７．７８s ;T P D ㊁O T P D ㊁M P T P D３种方法在阶段５的耗时分别为１３．７９s ㊁３４２．６４s ㊁２．０９s ,即M P T P D 方法在此阶段的耗时仅为T P D 的１５．１６％㊁O T P D 的０．６１％.其他３个试验数据亦表现出类似的规律.这证明了基元对滤波效率有着显著的影响,M P T P D 方法的关键点显著提高了其后续阶段的效率.但是,与T P D 方法相比,M P T P D 方法的点云分割和关键点提取耗费了较多的时间,因此该方法的整体效率低于T P D 方法.２．４㊀分析与讨论文献[４]指出T P D 方法对具有不同场景复杂度的点云数据均具有较高的滤波精度.本文的统计数据表明,T P D 方法的总误差T 平均值约１２．５２％,精度较高,符合既有结论.另外,表２表明T P D 方法所需的４个参数的取值对场景的变化不是很敏感,且具有显著的物理意义,根据实际情况微调参数取值即可.但T P D 方法存在对低位粗差㊁地形断裂敏感的问题,因此该滤波方法仍然存在一定的误差.表３表明T P D 方法的一类误差I 平均值约２７．０９％,显著地高于其他两种滤波方法的相关误差值.O T P D 方法是对T P D 方法的改进,具有对４１３１第１１期林祥国,等:融合点㊁对象㊁关键点等３种基元的点云滤波方法低位粗差㊁地形断裂不敏感的优势,但是比较耗时.统计数据表明,与T P D 相比,O T P D 方法的一类误差I ㊁总误差T 比T P D 的分别低约２２．１３％㊁８．４８％,但耗时是T P D 的２．４７倍.本文提出的M P T P D 方法,既有与O T P D 方法相当的滤波精度,又有更高的效率.统计数据表明,M P T P D 方法与O T P D 方法的一类误差Ⅰ㊁二类误差Ⅱ㊁总误差T 的差值绝对值分别为０．０７％㊁０．０３％㊁０．０４％,但是M P T P D 方法的平均耗时却仅有O T P D 方法的约５８％.与O T P D 方法的滤波精度相当㊁但效率显著提升的原因在于,M P T P D 方法中的关键点既能逼近原始点云的三维形态㊁又能显著地减少参与后续判别的计算量.另外,与T P D 方法相比,M P T P D 和O T P D 两种方法需要额外的３个参数.但是,表２表明４个试验数据中k 和r 两个参数的取值可相同,而α的取值在３个试验数据亦相同.这表明额外的３个参数亦 具有显著的物理意义,根据实际情况微调参数取值即可.表４㊀４个试验数据中３种滤波方法的时间花费T a b ．４㊀T h e t i m e c o s t s o f t h e t h r e e f i l t e r i n g m e t h o d s f o r t h e f o u r t e s t i n g da t a s e t s 试验数据方法阶段１耗时/s 阶段２耗时/s 阶段３耗时/s 阶段４耗时/s 阶段５耗时/s 总耗时/s １T P D２６１．１３ＧＧ５．３４７．８５２７４．３２O T P D Ｇ１７１．２６Ｇ５９．２７１３５．２９３５６．８２M P T P D Ｇ１７１．２６４７．７８４．２１１．３７２２４．６２２T P D １７９．４０ＧＧ３４．７７３３．８４２４８．０４O T P D Ｇ２３７．４５Ｇ３５．８８１６１．４８４３３．８１M P T P D Ｇ２３７．４５５５．３４１７．１３２．９４３１２．８６３T P D２５１．７０ＧＧ１３．０３３３．９７２９８．７０O T P D Ｇ３１６．２０Ｇ１８．３９５３１．０６８６５．６５M P T P D Ｇ３１６．２０８４．０５１０．２３４．９８４１５．４６４T P D１５０．７２ＧＧ９．４５１３．７９１７３．９６O T P DＧ３４４．５７Ｇ１０９．７０３４２．６４７９６．９１M P T P D Ｇ３４４．５７１１３．６３７．７８２．０９４６８．０７３㊀结㊀论点云是一种新型的数据源,其数据处理方法亟待研究.滤波是点云数据处理的一个必要的关键环节.目前,多数滤波方法采用单一的基元,但采用单一基元的滤波方法很难平衡滤波精度和滤波效率.为此,本文提出了M P T P D 方法,该方法在滤波的不同阶段采用了不同的基元.其涉及的基元包括点㊁对象㊁关键点等３种,且滤波的基本原理与T P D ㊁O T P D 方法相似.采用４个有代表性的点云数据进行了试验.试验表明,M P T P D方法具有整体上最优的性能.其中,精度方面,M P T P D 与O T P D 两种方法的精度相当,M P T P D 方法的一类误差I ㊁总误差T 比T P D 的相应误差低分别约２２．０７％㊁８．４４％;效率方面,多数情况下T P D ㊁M P T P D ㊁O T P D 方法的效率依次降低,但少数情况下M P T P D 的效率最高,且M P T P D 的平均耗时是O T P D 平均耗时的５７．９３％.笔者下一步的研究围绕两个方面开展:①采用并行计算技术提高滤波效率;②探索多基元的点云分类.参考文献:[１]㊀黄先锋,李卉,王潇,等．机载L i D A R 数据滤波方法评述[J ]．测绘学报,２００９,３８(５):４６６Ｇ４６９．HU A N G X i a n f e n g ,L I H u i ,WA N G X i a o ,e ta l ．F i l t e r A l go r i t h m s o fA i r b o r n eL i D A R D a t a :R e v i e w a n dP r o s Ｇp e c t s [J ]．A c t aG e o d a e t i c ae tC a r t o g r a ph i c aS i n i c a ,２００９,３８(５):４６６Ｇ４６９．[２]㊀王竞雪,朱庆,王伟玺．多匹配基元集成的多视影像密集匹配方法[J ]．测绘学报,２０１３,４２(５):６９１Ｇ６９８．WA N G J i n g x u e ,Z HU Q i n g,WA N G W e i x i ．A D e n s e M a t c h i n g A l g o r i t h mo fM u l t i Ｇv i e wI m a geB a s e d o n t h e I n Ｇt e g r a t e d M u l t i p l e M a t c h i n g P r i m i t i v e s [J ]．A c t a G e o d a e t i c a e t C a r t o g r a ph i c a S i n i c a ,２０１３,４２(５):６９１Ｇ６９８．[３]㊀张小红．机载激光雷达测量技术理论与方法[M ]．武汉:武汉大学出版社,２００７．Z HA N G X i a o h o n g ．T h e o r i e s ,M e t h o d s o f A i r b o r n e L i D A R T e c h n i q u e [M ]．W u h a n :W u h a n U n i v e r s i t y P r e s s ,２００７．５１３１。