留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

高斯混合聚类对移动曲面拟合滤波分类的应用

邢承滨 龚声胜 于晓亮 李易馨

邢承滨, 龚声胜, 于晓亮, 李易馨. 高斯混合聚类对移动曲面拟合滤波分类的应用[J]. 红外与激光工程, 2021, 50(10): 20200501. doi: 10.3788/IRLA20200501
引用本文: 邢承滨, 龚声胜, 于晓亮, 李易馨. 高斯混合聚类对移动曲面拟合滤波分类的应用[J]. 红外与激光工程, 2021, 50(10): 20200501. doi: 10.3788/IRLA20200501
Xing Chengbin, Gong Shengsheng, Yu Xiaoliang, Li Yixin. Application of Gaussian Mixture Clustering to moving surface fitting filter classification[J]. Infrared and Laser Engineering, 2021, 50(10): 20200501. doi: 10.3788/IRLA20200501
Citation: Xing Chengbin, Gong Shengsheng, Yu Xiaoliang, Li Yixin. Application of Gaussian Mixture Clustering to moving surface fitting filter classification[J]. Infrared and Laser Engineering, 2021, 50(10): 20200501. doi: 10.3788/IRLA20200501

高斯混合聚类对移动曲面拟合滤波分类的应用

doi: 10.3788/IRLA20200501
基金项目: 中央级公益性科研院所基本科研业务费专项资金(TKS20200318,SJY20200101)
详细信息
    作者简介:

    邢承滨,男,硕士生,主要从事机载雷达点云算法方面的研究

    通讯作者: 龚声胜,男,高级工程师,硕士,主要从事海洋测绘技术方面的研究。
  • 中图分类号: P237

Application of Gaussian Mixture Clustering to moving surface fitting filter classification

  • 摘要: 为了提高激光雷达点云滤波算法的精度和自适应性,对移动曲面滤波算法进行改进。采用格网边界点构建曲面约束条件,检验格网内是否全部为建筑物点。利用区域拟合求解地形的起伏,引入机器学习中高斯混合模型(GMM)对地形起伏进行滤波分类,将移动曲面中的种子点作为聚类算法中的靶向点参与分类学习。实验数据为雷达飞行的自测区,对于自测区采用随机抽样的方式,检验判断滤波效果。同时为检验GMM算法的准确性,在三类误差检验方式的基础上,增加了Kappa系数作为检验方式。通过与谱系聚类分类算法对比,证明所提算法能取得较好的滤波效果。
  • 图  1  测区格网划分和种子点选取

    Figure  1.  Grid division of survey area and selection of seed points

    图  2  点云真实高程与拟合高程分布

    Figure  2.  Point cloud real elevation and fitting elevation distribution

    图  3  角度阈值法检验建筑点

    Figure  3.  Angle threshold method to inspect building points

    图  4  三维数据的降维处理

    Figure  4.  Dimensionality reduction processing of 3D data

    图  5  EM-GMM模型的移动曲面算法

    Figure  5.  Moving surface algorithm of EM-GMM model

    图  6  工业区第28格网数据展示

    Figure  6.  Data display of the 28th grid in the industrial zone

    图  7  工业区第56格网数据展示

    Figure  7.  Data display of the 56th grid in the industrial zone

    图  8  两类聚类算法三类误差和 Kappa 系数分布

    Figure  8.  Three types of errors and Kappa coefficient distributions of the two types of clustering algorithms

    图  9  两类聚类算法对于样本Samp53数据处理展示

    Figure  9.  Two types of clustering algorithms show the data processing of Samp53

    图  10  两类聚类算法对于样本Samp53中指定格网数据处理展示

    Figure  10.  Two types of clustering algorithms show the processing of the specified grid data in Samp53

    表  1  3类误差与Kappa系数判定方法

    Table  1.   Three types of errors and Kappa coefficient judgment method

    Test sampleSample nameTopographic features
    Site1 Samp11 Vegetation and buildings on steep slopes
    Samp12 Small features
    Site2 Samp21 Narrow bridge
    Samp22 Bridge/passage
    Samp23 Complex and huge buildings, intermittent terrain
    Samp24 Slope
    Site3 Samp31 Cluster low-value points (multi-path effect)
    Site4 Samp41 Discontinuous terrain
    Samp42 Buildings, high-frequency undulating features
    Site5 Samp51 Vegetation on the slope
    Samp52 Low vegetation, fractured steep ridges
    Samp53 Intermittent terrain
    Samp54 Low-resolution building points
    Site6 Samp61 Low-resolution building points
    Site7 Samp71 Bridge, discontinuous terrain
    下载: 导出CSV

    表  2  测区地形特征与误差分布

    Table  2.   Topographic features and error distribution of the survey area

    Total numberType I errorType II errorTotal errorPAPCKappa
    pe=pa+pb+pc+pdpb/(pa+pb)pc/(pc+pd)(pb+pc)/pe(pa+pd)/pe[(pa+pb)(pa+pc)+(pc+pd)(pb+pd)]/pe2(PA-PC)/(1-PC)
    下载: 导出CSV

    表  3  两类聚类算法对于数据样本误差统计

    Table  3.   Two types of clustering algorithms for data sample error statistics

    Algorithm typeGaussian clustering model classification algorithmPedigree clustering model classification algorithm
    Sample nameType I error Type II error Total error Kappa coefficient Type I errorType II error Total errorKappa coefficient
    Samp113.43%14.38%5.40%81.79%22.25%6.38%13.82%72.00%
    Samp123.68%6.25%5.15%89.56%10.06%6.01%8.36%83.01%
    Samp211.59%4.52%2.54%94.19%8.89%8.64%8.81%80.37%
    Samp227.31%6.15%6.48%84.16%10.81%22.41%13.09%61.71%
    Samp239.16%9.81%9.56%80.13%15.38%16.82%16.26%66.47%
    Samp243.70%14.2%7.32%83.45%28.82%10.15%20.12%60.15%
    Samp314.04%6.24%5.22%89.51%9.18%12.50%10.00%74.48%
    Samp411.95%2.17%2.04%95.81%5.61%3.53%4.77%90.18%
    Samp425.56%3.87%4.21%87.53%7.42%7.46%7.45%84.46%
    Samp514.91%4.84%4.90%86.23%16.39%8.51%12.96%74.01%
    Samp529.36%10.53%9.72%77.84%14.67%7.27%11.54%76.76%
    Samp535.75%15.79%6.48%60.69%10.27%34.04%14.39%52.59%
    Samp545.34%1.30%3.16%93.62%5.13%12.23%8.98%82.03%
    Samp612.67%5.88%2.85%76.53%10.94%0.00%9.15%72.69%
    Samp713.23%7.25%4.26%88.60%5.29%12.90%7.80%82.26%
    下载: 导出CSV
  • [1] Yang Bisheng, Liang Fuxun, Huang Ronggang. Progress, challenges and perspectives of 3D LiDAR point cloud processing [J]. Acta Geodatica et Cartographica Sinica, 2017, 46(10): 1509-1516. (in Chinese) doi:  10.11947/j.AGCS.2017.20170351
    [2] Yan W Y, Shaker A, El-Ashmawy N. Urban land cover classification using airbone LiDAR data: a review [J]. Remote Sensing of Environment, 2015, 158: 295-310. doi:  http://dx.doi.org/10.1016/j.rse.2014.11.001
    [3] Li Deren, Ma Jun, Shao Zhenfeng. Innovation in the census and monitoring of geographical national conditions [J]. Geomatics and Information Science of Wuhan University, 2018, 43(1): 1-9. (in Chinese) doi:  10.13203/j.whugis20170356
    [4] Yang Bisheng, Dong Zhen. Progress and perspective of point cloud intelligence [J]. Acta Geodaetica et Cartographica Sinica, 2019, 48(12): 1575-1585. (in Chinese) doi:  10.11947/j.AGCS.2019.20190465
    [5] Li Deren. From geomatics to geospatial intelligent service science [J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(10): 1207-1212. (in Chinese) doi:  10.11947/j.AGCS.2017.20170263
    [6] Vosselman G, Mass H G. Airborne and terrestrial laser scanning[M]. Boca Ration, FL: CRC Press, 2010.
    [7] Hu Han, Ding Yulin, Zhu Qing, et al. An adaptive surface filter for airborne laser scanning point clouds by means of regularization and bending energy [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 92: 98-111.
    [8] Sui Lichun, Zhang Yibin, Liu Yan, et al. Filtering of airborn LiDAR point cloud data based on the adaptive mathematical morphology [J]. Acta Geodaeticaet Cartographica Sinica, 2010, 39(4): 390-396. (in Chinese)
    [9] Sui Lichun, Yang Yun. Filtering of airborne LiDAR point cloud data based on car(p, q) model and mathematical morphology [J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(2): 219-224. (in Chinese)
    [10] Zhu Xiaoxiao, Wang Cheng, Xi Xiaohuan, et al. Research progress of ICESat-2/ATLAS data processing and applications [J]. Infrared and Laser Engineering, 2020, 49(11): 20200259. (in Chinese) doi:  10.3788/IRLA20200259
    [11] Qin Wei, Chen Xi, Ma Yuanyuan, et al. The analysis of edge detection algorithms based on mathematical morphology [J]. Information Technology, 2019, 43(11): 33-36. (in Chinese)
    [12] Zhao Chuan, Zhang Baoming, Guo Haitao, et al. Airborne LiDAR point cloud classification using transfer learning [J]. Optics and Precision Engineering, 2019, 27(7): 1601-1612. (in Chinese) doi:  10.3788/OPE.20192707.1601
    [13] Liu Zhiqing, Li Pengcheng, Chen Xiaowei, et al. Classification of airborne LiDAR point cloud data based on information vector machine [J]. Optics and Precision Engineering, 2016, 24(1): 210-219. (in Chinese) doi:  10.3788/OPE.20162401.0210
    [14] Du Ruijian, Ge Baozhen, Chen Lei. Texture mapping of multi-view high-resolution images and binocular 3D point clouds [J]. Chinese Optics, 2020, 13(5): 1055-1064. (in Chinese) doi:  10.37188/CO.2020-0034
    [15] Zong Wenpeng, Li Guangyun, Li Minglei, et al. A survey of laser scan matching methods [J]. Chinese Optics, 2018, 11(6): 914-930. (in Chinese) doi:  10.3788/CO.20181106.0914
    [16] Chen Ying, Zhu Ming, Li Zhaoze. Remote sensing digital image enhancement based on Gaussian mixture modeling [J]. Chinese Journal of Lasers, 2014, 41(12): 1209002. (in Chinese)
    [17] Tao Zhiyong, Liu Xiaofang, Wang Hezhang. Clustering algorithm of Gaussian mixture model based on density peaks [J]. Journal of Computer Applications, 2018, 38(12): 3433-3437, 3443. (in Chinese)
    [18] Feng Fajie, Ding Yazhou, Li Junping. Airborne LiDAR point cloud filtering using saliency division [J]. Infrared and Laser Engineering, 2020, 49(8): 20190439. (in Chinese) doi:  10.3788/IRLA20190439
    [19] Sun Hao, Guo Yingqing, Zhao Wanli. Improved model for on-board real-time by constructing empirical model via GMM clustering method [J]. Journal of Northwestern Polytechnical University, 2020, 38(3): 507-514. (in Chinese) doi:  10.1051/jnwpu/20203830507
    [20] Sithole G, Vosselman G. A comparison of existing automatic filters[R].Techniacal Report, ISPRS test on extracting DEMs from point clouds, 2003.
  • [1] 范伟, 刘博, 蒋贇.  基于体光栅窄带光学滤波的激光雷达收发波长动态匹配技术研究 . 红外与激光工程, 2022, 51(7): 20210639-1-20210639-7. doi: 10.3788/IRLA20210639
    [2] 靳辰飞, 田小芮, 唐勐, 王峰, 杨杰, 乔凯, 史晓洁, 张思琦.  非视域三维成像激光雷达的研究进展 . 红外与激光工程, 2022, 51(3): 20210471-1-20210471-16. doi: 10.3788/IRLA20210471
    [3] 郭静菁, 费晓燕, 葛鹏, 周安然, 王磊, 李正琦, 盛磊.  基于全光纤光子计数激光雷达的高分辨率三维成像 . 红外与激光工程, 2021, 50(7): 20210162-1-20210162-10. doi: 10.3788/IRLA20210162
    [4] 张欣, 李思宁, 孙剑峰, 姜鹏, 刘迪, 王鹏辉.  基于Mean Shift点法向量分类的目标三维姿态估计 . 红外与激光工程, 2020, 49(S2): 20200109-20200109. doi: 10.3788/IRLA20200109
    [5] 张楠, 孙剑峰, 姜鹏, 刘迪, 王鹏辉.  激光雷达场景三维姿态点法向量估计方法 . 红外与激光工程, 2020, 49(1): 0105004-0105004(8). doi: 10.3788/IRLA202049.0105004
    [6] 龚道然, 李思宁, 姜鹏, 刘迪, 孙剑峰.  激光雷达三维距离像超分辨重构方法研究 . 红外与激光工程, 2020, 49(8): 20190511-1-20190511-7. doi: 10.3788/IRLA20190511
    [7] 曹杰, 郝群, 张芳华, 徐辰宇, 程阳, 张佳利, 陶禹, 周栋, 张开宇.  APD三维成像激光雷达研究进展 . 红外与激光工程, 2020, 49(9): 20190549-1-20190549-10. doi: 10.3788/IRLA20190549
    [8] 胡善江, 贺岩, 陶邦一, 俞家勇, 陈卫标.  基于深度学习的机载激光海洋测深海陆波形分类 . 红外与激光工程, 2019, 48(11): 1113004-1113004(8). doi: 10.3788/IRLA201948.1113004
    [9] 李小路, 曾晶晶, 王皓, 徐立军.  三维扫描激光雷达系统设计及实时成像技术 . 红外与激光工程, 2019, 48(5): 503004-0503004(8). doi: 10.3788/IRLA201948.0503004
    [10] 李晶, 车英, 宋暖, 翟艳男, 陈大川, 李君.  三维激光雷达共光路液体透镜变焦光学系统设计 . 红外与激光工程, 2019, 48(4): 418002-0418002(9). doi: 10.3788/IRLA201948.0418002
    [11] 邓迁, 吴德成, 况志强, 刘东, 谢晨波, 王英俭.  用于水汽混合比自标定的532 nm/660 nm双波长激光雷达 . 红外与激光工程, 2018, 47(12): 1230004-1230004(5). doi: 10.3788/IRLA201847.1230004
    [12] 马跃, 张文豪, 张智宇, 马昕, 李松.  基于半解析模型的激光测高回波海水海冰波形分类方法 . 红外与激光工程, 2018, 47(5): 506005-0506005(7). doi: 10.3788/IRLA201847.0506005
    [13] 张欣婷, 安志勇, 亢磊.  三维激光雷达发射/接收共光路光学系统设计 . 红外与激光工程, 2016, 45(6): 618004-0618004(5). doi: 10.3788/IRLA201645.0618004
    [14] 史风栋, 刘文皓, 汪鑫, 丁娟, 史屹君, 修春波.  室内激光雷达导航系统设计 . 红外与激光工程, 2015, 44(12): 3570-3575.
    [15] 曹开法, 黄见, 胡顺星.  边界层臭氧差分吸收激光雷达 . 红外与激光工程, 2015, 44(10): 2912-2917.
    [16] 姜成昊, 杨进华, 张丽娟, 李祥.  新型多普勒成像激光雷达原理设计与仿真 . 红外与激光工程, 2014, 43(2): 411-416.
    [17] 李小珍, 吴玉峰, 郭亮, 曾晓东.  合成孔径激光雷达下视三维成像构型及算法 . 红外与激光工程, 2014, 43(10): 3276-3281.
    [18] 杨文秀, 付文兴, 周志伟, 余巍, 马杰.  基于投影降维的激光雷达快速目标识别 . 红外与激光工程, 2014, 43(S1): 1-7.
    [19] 曹文勇, 马明, 赵彬, 何幸锴, 侯天晋, 周鼎富, 邓华荣.  滤波器组在相干激光测风雷达信号处理中的应用 . 红外与激光工程, 2013, 42(5): 1179-1183.
    [20] 孙崇利, 苏伟, 武红敢, 刘睿, 刘婷, 黄健熙, 朱德海, 张晓东, 刘峻明.  改进的多级移动曲面拟合激光雷达数据滤波方法 . 红外与激光工程, 2013, 42(2): 349-354.
  • 加载中
图(10) / 表(3)
计量
  • 文章访问数:  231
  • HTML全文浏览量:  96
  • PDF下载量:  32
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-12-18
  • 修回日期:  2021-02-24
  • 刊出日期:  2021-10-20

高斯混合聚类对移动曲面拟合滤波分类的应用

doi: 10.3788/IRLA20200501
    作者简介:

    邢承滨,男,硕士生,主要从事机载雷达点云算法方面的研究

    通讯作者: 龚声胜,男,高级工程师,硕士,主要从事海洋测绘技术方面的研究。
基金项目:  中央级公益性科研院所基本科研业务费专项资金(TKS20200318,SJY20200101)
  • 中图分类号: P237

摘要: 为了提高激光雷达点云滤波算法的精度和自适应性,对移动曲面滤波算法进行改进。采用格网边界点构建曲面约束条件,检验格网内是否全部为建筑物点。利用区域拟合求解地形的起伏,引入机器学习中高斯混合模型(GMM)对地形起伏进行滤波分类,将移动曲面中的种子点作为聚类算法中的靶向点参与分类学习。实验数据为雷达飞行的自测区,对于自测区采用随机抽样的方式,检验判断滤波效果。同时为检验GMM算法的准确性,在三类误差检验方式的基础上,增加了Kappa系数作为检验方式。通过与谱系聚类分类算法对比,证明所提算法能取得较好的滤波效果。

English Abstract

    • 激光雷达(Light Detection and Ranging, LiDAR)是当前流行的一种主动获取点云数据的手段。经过20多年的发展,激光雷达的硬件技术已经取得长足的进步,当前流行的机载平台有Leica3(ALS80-UP),Optech4(Titan),北京天绘5(AP-3500)等,都采用1 064 nm的波长。激光雷达主要面向中小型区域地形植被、建筑物、水下地形等三维数据获取[1-3]

      激光雷达系统测量原理是利用发射激光束测量与目标之间距离,结合发射器位置和姿态联合解算,得到目标区域的三维点云(X,Y,Z)。相较于传统航空摄影,激光雷达不仅能够全天候作业,并且抗干扰能力较强,点云数据量丰富。当前点云数据是继三大传统地理空间产业(GNSS,GIS,RS)后的第四大产业。但是点云数据也存在缺陷,例如数据量庞大,高冗余,采样粒度分布不均等[4-7]

      点云数据的后处理一直是难题,其中雷达点云滤波(雷达点云地面点和非地面点的分离称为滤波)更是数据处理的关键。点云滤波技术发展已有几十年的历史,众多学者提出了多种滤波算法,并且不断完善。如今滤波算法主要分为四类:(1)形态学滤波算法[8-11];(2)基于坡度的滤波算法[12];(3)基于曲面拟合的滤波算法[13-14];(4)基于聚类分割的滤波算法[15]

      移动曲面滤波算法最早由武汉大学的张小红提出,文中基于传统的移动曲面算法,引入机器学习中高斯混合模型(Gaussian Mixture Model, GMM)聚类算法对地形进行分类,为检验文中算法的分类效果,与谱系聚类算法(另一种常见的机器学习聚类方法)进行实验对比,判断两种机器学习算法的滤波效果。

    • 移动曲面算法的原理将测区划分为规则格网,采用移动窗口的方式对所有数据点进行遍历。对规则格网进行种子点选择,利用种子点拟合地面模型,求解真实数据点的高程与拟合高程之间的高差。

    • 文中采用直方图方法对数据进行粗差剔除,初步剔除粗差有利于种子点的准确选择。对于测区进行格网划分,如图1(a)所示,ij表示测区格网建立的行列号。横纵步长相等,即dx=dy。建立格网的目的是保证每一个数据点都有唯一索引号,同一格网内数据点索引号一致的。

      图  1  测区格网划分和种子点选取

      Figure 1.  Grid division of survey area and selection of seed points

      文中采用大小$3 \times 3$的移动格网,单个格网内最低点作为种子点(再经过剔除粗差的数据区域,格网内最低点通常确定为地面点)。移动窗口能较好覆盖建筑区域,并且格网尺寸不宜太大,尽量简化格网内地物种类,为后续建立分类做准备。如图1(b)所示,星号点为移动格网内种子的地面点,利用种子点建立初始地面。种子点的选择对于文中算法至关重要,不仅初期进行曲面拟合,后期更会作为靶向点参与地面点的分类。

    • 文中采用二次曲面拟合方式建立初始地面模型。

      $$ {\textit{z}} = a{x^2} + b{y^2} + cxy + dx + ey + f $$ (1)

      式中:z为种子点的高程;xy为种子点的平面坐标;af为拟合方程的参数。利用观测值建立误差方程:

      $$ \left\{ {\begin{array}{*{20}{c}} {{{\textit{z}}_1} = {a_0} + {a_1}{x_1} + {a_2}{y_1} + {a_3}x_1^2 + {a_4}{x_1}{y_1} + {a_5}y_1^2} \\ {{{\textit{z}}_2} = {a_0} + {a_1}{x_2} + {a_2}{y_2} + {a_3}x_2^2 + {a_4}{x_2}{y_2} + {a_5}y_2^2} \\ \vdots \\ {{{\textit{z}}_n} = {a_0} + {a_n}{x_n} + {a_n}{y_n} + {a_n}x_n^2 + {a_n}{x_n}{y_n} + {a_5}y_n^2} \end{array}} \right. $$ (2)

      利用最小二乘原理,求解上述二次曲面的6个参数,再利用求解的拟合参数建立地面模型。如图2(a)所示,为某格网内数据真实高程与拟合高程的分布。图2(b)表示真实高差与拟合高差的分布。通常,地面点拟合高程与真实高程近似,地物点拟合高程与真实高程差异较大,高差为地物实际高度,这是为后续分类提供重要理论。

      图  2  点云真实高程与拟合高程分布

      Figure 2.  Point cloud real elevation and fitting elevation distribution

    • 多种滤波算法将格网内最低点定义为种子点,但是存在部分特殊情况,例如格网内地形都为建筑的特殊地形,后续的分类将建筑点全部误分类为地面点,造成对数据的错误判断。对于这种特殊地形,最低点判断为种子点就会存在误差。建筑地区有着明显的特征,边界存在明显高差,建筑区域内高程起伏较小,这也是与地面点容易混淆的特征。文中基于建筑地区的特征,提出角度阈值检验算法判断格网内是否全为建筑点。

      图3所示,abcd为同一格网中的边界点,判断最低点与模拟平面的角度,经过众多实验判断出$\max $$ \left\{ {{\theta _1},{\theta _2},{\theta _3},{\theta _4}} \right\} \leqslant {5^ \circ }$,表示平面地形较为平滑。对于较为平坦的地形,进行二次判断。二次判断采用梯度差方式进行。建筑顶端表面平整,将格网内平均高程与周围临近4个格网内平均高程逐次比较。经过大量实验计算,平均高程与周围邻近格网平均高程阈值为8 cm,对建筑地区滤波效果最佳,阈值范围内该格网定义为建筑格网,格网内全部定义为地面点。当前还没获取一种高效的科学计算手段判断角度阈值与梯度阈值,判断方式还停留在统计计算。阈值设定主要判断测区内是否完全为建筑点,这种验证方式只对存在大量建筑的区域效果明显,是一种先验条件下的判别方式;对于非密集建筑区意义不大。通常所选测区都比较复杂,先验一下测区数据,避免格网全部为建筑点造成的误差。

      图  3  角度阈值法检验建筑点

      Figure 3.  Angle threshold method to inspect building points

    • 上文提出地面点拟合高程与真实高程近似,地物点拟合高程与真实高程之间高差近似为地物高度。如图4所示,A1为地面点,B2为相邻地物点(非地面点),A1B2水平距离为$\Delta d$$ \Delta {h}_{2} $A点与拟合值高差,$ \Delta {h}_{1} $B点与拟合值高差。AB两点的地形起伏为$ \Delta {h}_{1}+\Delta {h}_{2} $。相邻点的水平距离和高差距离作为GMM聚类算法的数据基础。将三维数据降维为二维起伏数据,利用GMM聚类算法进行学习判断。

      图  4  三维数据的降维处理

      Figure 4.  Dimensionality reduction processing of 3D data

    • GMM算法是一种在机器学习中广泛使用的聚类方式。GMM算法原理是将事物分解为若干的基于高斯概率密度的函数,每一个混合成分对应一个高斯分布,任意形状的概率分布都可以用多个高斯分布函数近似表达[16-17]。GMM算法原理概率模型为:

      $$ P\left( {{{y}}|\theta } \right) = \sum\limits_{k = 1}^K {{a_k}} \phi \left( {y|{\theta _k}} \right) $$ (3)

      式中:K为最大高斯概率密度函数的个数;${a_k}$为每一个高斯分布的权重,${a_k} \geqslant 0$$\phi \left( {y|{\theta _k}} \right)$为第k个高斯分布的概率密度函数;y为观测值。其中$ \theta = ( {{a_1},{a_2}, \cdots {a_k};{\theta _1}}, $$ {{\theta _2}, \cdots {\theta _k}}),{\theta _k} = \left( {{\mu _k},a_k^2} \right),k = 1,2, \cdots k $$\theta $为未知参数,包括权重${a_k}$${\theta _{{k}}}$${\theta _{{k}}}$包括期望${\;\mu _{\rm{k}}}$和方差$a_k^2$,高斯混合k个模型需要估计3 k个参数。

      $$ \phi \left( {y|{\theta _k}} \right) = \frac{1}{{\sqrt {2\pi } {\sigma _k}}}\exp \left( { - \frac{{{{\left( {y - {\mu _k}} \right)}^2}}}{{2\sigma _k^2}}} \right) $$ (4)

      公式(4)中,每个观测样本出现的概率表示为k个高斯分布的加权。聚类的核心思想是对于某个样本${y_j}$,把该样本代入到k个高斯分布中求出属于每个类别的概率,公式(5)表示第k个类别概率:

      $$ \frac{{{a_k}\phi \left( {{y_j}|{\theta _k}} \right)}}{{\displaystyle\sum\limits_{k = 1}^{{a_k}\phi \left( {{y_j}|{\theta _k}} \right)} {{a_k}\phi \left( {{y_j}|{\theta _k}} \right)} }} $$ (5)
    • EM(Expectation-Maximization)算法又称最大期望算法,要估计的参数$\theta =( {{a_1},{a_2}, \cdots {a_k};{\theta _1}}, {\theta _2}, \cdots {{\theta _k}}), $$ {\theta _k} = \left( {{\mu _k},a_k^2} \right),k = 1,2, \cdots k$。采用极大似然估计参数$\theta $,使得观测值数据y的对数似然函数$L\left( \theta \right) = \log p\left( {y|\theta } \right)$极大化[18-19]。公式(6)、(7)为建立的似然方程:

      $$ {\theta ^*} = \arg \mathop {\max }\limits_\theta L\left( \theta \right) $$ (6)
      $$ L\left( \theta \right) = \log P\left( {y|\theta } \right) = \sum\limits_{j = 1}^N {\left[ {\log \left( {\sum\limits_{k = 1}^k {{a_k}\phi \left( {{y_j}|{\theta _k}} \right)} } \right)} \right]} $$ (7)

      根据EM算法求解参数$\theta $。可以获取观测值y,但是确定观测数据来自第k个模型的数据却是未知,用隐函数变量z表示。EM算法的迭代过程分为E步和M步:

      (1) 选择参数的初始值${\theta _{\rm{0}}}$开始进行迭代;

      (2) E(Expection)步:记第n次迭代参数为${\theta _{{n}}}$,因此n+1的E步为:

      $$ Q\left( {\theta |{\theta _n}} \right) = \sum\limits_{i = 1}^n {\sum\limits_{\textit{z}} {p{\theta _n}\left( {{\textit{z}}|{y_i}} \right)} } {\log ^{p\theta \left( {{y_i},{\textit{z}}} \right)}} $$ (8)

      $Q$函数是指在给定观测数据y和第n次迭代的参数${\theta _n}$时,对数似然函数的$\log p\theta \left( {y,{\textit{z}}} \right)$的期望;

      (3) M(Maximization)步:求得$Q\left( {\theta |{\theta _n}} \right)$最大化的$\theta \left( {n + 1} \right)$,即确定第n+1次的模型参数;

      $$ {\theta _{n + 1}} = \arg \mathop {\max }\limits_\theta Q\left( {\theta |{\theta _n}} \right) $$ (9)

      (4) 重复步骤(2)、(3)直到收敛;

      利用EM-高斯混合模型聚类求解出参数$\theta $$\theta $参数中包含均值向量$\mu _i'$和协方差矩阵。

      计算期望值$\mu _i'$

      $$ u_i' = \frac{{\displaystyle\sum\limits_{k = 1}^m {p\left( {{{\textit{z}}_k} = i|{y_k}} \right) \cdot {{{y}}_k}} }}{{\displaystyle\sum\limits_{k = 1}^m {p\left( {{{\textit{z}}_k} = i|{y_k}} \right)} }} $$ (10)

      计算协方差矩阵$\sigma _{{k}}^{2'}$为:

      $$ \sigma _k^{2'} = \frac{{\displaystyle\sum\limits_{k = 1}^m {p\left( {{{\textit{z}}_k} = i|{y_k}} \right){{\left( {{y_k} - u_i'} \right)}^{\rm T}}} }}{{\displaystyle\sum\limits_{{\rm{k}} = 1}^m {p\left( {{{\textit{z}}_k} = i|{y_k}} \right)} }} $$ (11)

      计算混合系数$a_k'$为:

      $$ a_k' = \frac{{\displaystyle\sum\limits_{k = 1}^m {p\left( {{{\textit{z}}_k} = i|{y_k}} \right)} }}{m} $$ (12)

      式中:Zk为第k个高斯模型。利用最大期望迭代求解出高斯混合模型的权重和期望方差,再利用高斯混合模型判断点的聚类结果。

    • 高斯聚类结果不能获取数据的属性,即聚类只是表示同一类型点的集合,并不能显示地面点为哪一类。因此对于地面点类型需要新的判断方式,上文提及格网内种子点为地面点,种子点所在类型就是地面点。种子点作为靶向点参与滤波分类,种子点所在聚类结果中的类型判定为地面点类型。EM-高斯混合模型对于移动曲面算法改进算法流程如图5所示。

      图  5  EM-GMM模型的移动曲面算法

      Figure 5.  Moving surface algorithm of EM-GMM model

    • 文中采用定性和定量方式检验滤波分类效果,定性分析中采取随机抽样的方式检验滤波效果。

    • 定性分析主要通过目视方式对滤波效果进行判断,对比滤波前后的地形检验滤波效果。文中实验采用定性分析的数据由华测AS300便携机载激光雷达扫描系统获取。测域为天津某一工业区域,内部设有建筑厂房、主办公楼、湖泊、树木、汽车等小型地物。采用随机抽样的方式判断滤波效果。测区内数据点约为两千万,随机抽取某一格网数据进行展示,实验抽取第28格网为例。因文中测区数据平面坐标经过重心化处理,XYZ会有部分呈现负值。

      图6(a)、(b)表示测区滤波分类前后点云分布,图6(c)表示聚类结果。经过多次实验,通常将高斯混合模型定义为4个。树木、建筑地形、地面点、独立小地物,比如汽车等,此格网为一个包含部分建筑区域的格网。树木呈现为渐变式高差,高差逐渐起伏;建筑物呈现两极分化,分为距离陡变和逐渐变化;地面点由于点云数据量十分庞大,地形起伏较为平缓,绝大多数数据量集中在0附近,由于每个格网内存在靶向点,靶向点的数据类型可以直接确定地面点的类型;独立地物变化区域较小。如图6(a)、(b)所示,平均高程下降七八米左右,正是测区内建筑物的高度,区域内较高建筑地区点云被剔除。图6(d)、(e)表示点云滤波前后的地面模型,对于图6(d)中标注的建筑地物很好过滤,地面模型保留较好。图6(f)、 (g)为等高线表示,经过滤波处理,图6(g)中不存在图6(f)中明显的建筑物闭合,边界密集的等高线分布。GMM算法对于建筑物有着较好的过滤效果。

      图  6  工业区第28格网数据展示

      Figure 6.  Data display of the 28th grid in the industrial zone

      抽取第56个格网进行滤波效果分析。由等高线结合点云分布,测区为一个存在独立地物的斜坡地形。图7(a)为聚类后4类地物,图7(b)为滤波前的地面模型,图7(c)为滤波后的地面模型。图7(b)中地形可以判断存在独立地物,比如旗杆、汽车等。这类地物变化区域较小,独立地物大小相对有限,文中算法可以很好地过滤其中的独立地物,将图7(b)中两处独立地物很好地剔除。根据图7(d)中等值线的分布可以判断独立地物的存在。两处相对独立的等高线,而图7(e)中不存在这两处区域等高线,只显示逐渐起伏的地形。

      图  7  工业区第56格网数据展示

      Figure 7.  Data display of the 56th grid in the industrial zone

      虽然随机抽样对于滤波效果判断存在片面性,但是随机抽样方式是对点云类型检验滤波效果的有效方式,抽样数据能一定程度反映全部滤波效果。

    • 此次定量实验数据采用2003年摄影测量与遥感协会(ISPRS)公布的雷达点云数据。该数据为Optech ALTM机载激光扫描仪获取的Vaihingen/Enz和Stuttgart市的点云数据。其中城市数据:Site1~Site4;村庄数据:Site5~Site8,表1对每个样本的属性进行详细解释。

      表 1  3类误差与Kappa系数判定方法

      Table 1.  Three types of errors and Kappa coefficient judgment method

      Test sampleSample nameTopographic features
      Site1 Samp11 Vegetation and buildings on steep slopes
      Samp12 Small features
      Site2 Samp21 Narrow bridge
      Samp22 Bridge/passage
      Samp23 Complex and huge buildings, intermittent terrain
      Samp24 Slope
      Site3 Samp31 Cluster low-value points (multi-path effect)
      Site4 Samp41 Discontinuous terrain
      Samp42 Buildings, high-frequency undulating features
      Site5 Samp51 Vegetation on the slope
      Samp52 Low vegetation, fractured steep ridges
      Samp53 Intermittent terrain
      Samp54 Low-resolution building points
      Site6 Samp61 Low-resolution building points
      Site7 Samp71 Bridge, discontinuous terrain

      滤波误差类型分为3种:I类误差(Type I error),将地面点误分为非地面点;II类误差(Type II error),将非地面点误分为地面点;总误差(Total error),所有误分点数/数据点总数[20]。为增加数据的稳健性测度,文中同时引用了Kappa系数。3类误差与Kappa系数的计算公式如表2所示。

      此次定量实验采用对比分析的方式,引入另一种机器学习算法:谱系聚类算法。两种聚类算法对同一样本进行分类分析,判断哪种算法滤波分类效果更好。如表3所示。

      表 2  测区地形特征与误差分布

      Table 2.  Topographic features and error distribution of the survey area

      Total numberType I errorType II errorTotal errorPAPCKappa
      pe=pa+pb+pc+pdpb/(pa+pb)pc/(pc+pd)(pb+pc)/pe(pa+pd)/pe[(pa+pb)(pa+pc)+(pc+pd)(pb+pd)]/pe2(PA-PC)/(1-PC)

      表 3  两类聚类算法对于数据样本误差统计

      Table 3.  Two types of clustering algorithms for data sample error statistics

      Algorithm typeGaussian clustering model classification algorithmPedigree clustering model classification algorithm
      Sample nameType I error Type II error Total error Kappa coefficient Type I errorType II error Total errorKappa coefficient
      Samp113.43%14.38%5.40%81.79%22.25%6.38%13.82%72.00%
      Samp123.68%6.25%5.15%89.56%10.06%6.01%8.36%83.01%
      Samp211.59%4.52%2.54%94.19%8.89%8.64%8.81%80.37%
      Samp227.31%6.15%6.48%84.16%10.81%22.41%13.09%61.71%
      Samp239.16%9.81%9.56%80.13%15.38%16.82%16.26%66.47%
      Samp243.70%14.2%7.32%83.45%28.82%10.15%20.12%60.15%
      Samp314.04%6.24%5.22%89.51%9.18%12.50%10.00%74.48%
      Samp411.95%2.17%2.04%95.81%5.61%3.53%4.77%90.18%
      Samp425.56%3.87%4.21%87.53%7.42%7.46%7.45%84.46%
      Samp514.91%4.84%4.90%86.23%16.39%8.51%12.96%74.01%
      Samp529.36%10.53%9.72%77.84%14.67%7.27%11.54%76.76%
      Samp535.75%15.79%6.48%60.69%10.27%34.04%14.39%52.59%
      Samp545.34%1.30%3.16%93.62%5.13%12.23%8.98%82.03%
      Samp612.67%5.88%2.85%76.53%10.94%0.00%9.15%72.69%
      Samp713.23%7.25%4.26%88.60%5.29%12.90%7.80%82.26%

      定量实验可知,在15组样本中,高斯聚类I,II类误差均在10%,表明高斯聚类分类准确度较高,相较于谱系聚类中,部分样本中I,II类误差超过10%,在Samp24中,谱系聚类I类误差28.82%,相较于高斯聚类的3.7%,存在较大误差。由图8直方图可以更加直观判断,两类聚类算法的误差分布,另外通过Kappa系数对比,高斯聚类稳健性较高。

      图  8  两类聚类算法三类误差和 Kappa 系数分布

      Figure 8.  Three types of errors and Kappa coefficient distributions of the two types of clustering algorithms

      从15组样本中不同样本数据可知,文中聚类算法对于平坦地区分类效果较好,稳健性较高。对于个别间断样本,例如样本Samp53,误差分布较大并且稳健性较低,高斯聚类在60.69%,谱系聚类低至52.59%。

      此次算法重点分析两算法统一存在误差最大和稳健性最差的样本Samp53。对于Samp53样本,采用高斯聚类与谱系聚类进行对比分析。图9(a)(d)分别为高斯聚类和谱系聚类结果显示;图9(b)(e)为两类算法误差分布与地面模型;图9(c)(f)为两类算法分类后的地面模型。由上述可知,Samp53为间断地形,如图9(b)所示,I,II类误差主要集中在地物点与地面点的边界处。I类误差过大的主要原因是陡变的地面点,拟合高程与真实高程高差过大,相邻点起伏过大被分为地物点,而边界点作为地面点,造成错误的分类。真实地形与拟合地形在间断地区存在较大差异。

      图  9  两类聚类算法对于样本Samp53数据处理展示

      Figure 9.  Two types of clustering algorithms show the data processing of Samp53

      尽管两类滤波对于间断地形都略有不足,但是高斯混合模型与谱系聚类算法针间断地形分类效果存在较大差异,以图9(b)中II类格网数据为例。图10(a)~(c)为高斯聚类分类:图10(a)为根据格网内靶向点判断地面类型,其他属性合并为地物点;图10(b)为分类后数据属性,图10(c)为高斯聚类滤波效果。同理,图10(d)~(f)为谱系聚类分类。

      图  10  两类聚类算法对于样本Samp53中指定格网数据处理展示

      Figure 10.  Two types of clustering algorithms show the processing of the specified grid data in Samp53

      通过对图10(c)图10(f)进行对比发现,高斯聚类中误差数目,尤其是II类误差明显小于谱系聚类。II类误差主要分布在山脊线,I类误差主要分布在山体边缘线。

      高斯聚类和谱系聚类同样是采用相邻点地形起伏作为聚类依据,造成分类差异的主要原因为高斯聚类里利用高斯曲线模拟地物类型,利用不同类别高斯曲线差异完成分类。高斯聚类的优势在于大量数据分类噪声不敏感,对于距离和密度分类,适合高维特征,同时高斯聚类可以完成较多形状的聚类。

      而谱系聚类对于噪声非常敏感,并且需要在测试前准备判断类别的个数。陡峭地区存在个别噪声,对于谱系聚类算法分类效果影响较大,同时在山体地区,主要地物不再是建筑,而是树木、斜坡等,地物种类不固定于4类,谱系聚类按照同一4类特征处理,也会造成部分误差。

    • 文中通过移动曲面方式,利用相邻数据点高程起伏作为高斯混合模型分类的依据,对于测区内数据进行分类,利用种子点作为靶向点进行分类属性判断。通过定性和定量实验分析,该算法可以整体取得较好的滤波效果。传统滤波算法与新型机器学习技术的结合为点云滤波技术的发展提供了新型思路。引入高斯混合模型可以充分利用数据点之间的相关性,对于点云数据拓扑关系不明确的问题有较好改善。对于其中仍存在的问题,例如定义格网数据类型固定为4种,部分特殊地形不能用高斯曲线很好拟合等,后续工作将会逐步完善。

参考文献 (20)

目录

    /

    返回文章
    返回