CN112184736B - 一种基于欧式聚类的多平面提取方法 - Google Patents

一种基于欧式聚类的多平面提取方法 Download PDF

Info

Publication number
CN112184736B
CN112184736B CN202011078122.4A CN202011078122A CN112184736B CN 112184736 B CN112184736 B CN 112184736B CN 202011078122 A CN202011078122 A CN 202011078122A CN 112184736 B CN112184736 B CN 112184736B
Authority
CN
China
Prior art keywords
grid
plane
point
ground
map
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011078122.4A
Other languages
English (en)
Other versions
CN112184736A (zh
Inventor
张雪波
刘皓月
宋锐
方勇纯
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nankai University
Original Assignee
Nankai University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nankai University filed Critical Nankai University
Priority to CN202011078122.4A priority Critical patent/CN112184736B/zh
Publication of CN112184736A publication Critical patent/CN112184736A/zh
Application granted granted Critical
Publication of CN112184736B publication Critical patent/CN112184736B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C15/00Surveying instruments or accessories not provided for in groups G01C1/00 - G01C13/00
    • G01C15/002Active optical surveying means
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/005Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 with correlation of navigation data from several sources, e.g. map or contour matching
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/23Clustering techniques
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/27Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T15/003D [Three Dimensional] image rendering
    • G06T15/005General purpose rendering architectures
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/136Segmentation; Edge detection involving thresholding
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10028Range image; Depth image; 3D point clouds
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20024Filtering details

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Automation & Control Theory (AREA)
  • Medical Informatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)
  • Image Analysis (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)

Abstract

本发明公开一种基于欧式聚类的多平面提取方法,可以适用于多种野外复杂环境。首先,针对树冠、隧道、矿洞等一类具有悬挂物的环境,提出基于欧式聚类的混合高度图方法。通过将映射到二维栅格中三维点云进行聚类,并利用混合高度图法区分地面点和悬挂物点,实现了对地面上方障碍物的准确检测和下方对应可通行区域的正确提取。除此之外,针对野外崎岖环境,本发明还提出了基于随机采样一致性的多平面提取方法。通过在gazebo虚拟仿真环境和真实校园环境中进行了实验验证。实验结果表明,本发明方法与已有算法相比在点云分割的精度和效率方面均有明显提高。

Description

一种基于欧式聚类的多平面提取方法
技术领域
本发明属于点云分割和移动机器人可通行域提取的技术领域,特别是涉及一种基于欧式聚类的多平面提取方法。
背景技术
目前自主移动机器人[1]被广泛应用到各个领域中,其相关技术也成为研究热点。其中,环境感知技术是移动机器人实现自主导航的关键。户外机器人的应用领域包括无人驾驶汽车、勘测机器人、救援机器人等[2-4],而正确的检测出周围环境中的障碍物,识别出机器人能够安全穿越的区域是此类机器人的一项最基本的要求。由于户外非结构化环境下的地形变化较大、障碍物形状不规则等原因,实现机器人的自主导航[5]具有复杂性。除了检测树木、石块、陡坡等常规障碍物以避免与它们碰撞外,陷坑、洞穴等更是增加了检测难度。这些都使得可通行区域的检测成为机器人导航的首要任务[2,6]。
户外环境的非结构化特性要求机器人在三维空间中感知周围环境。近年来,多线三维激光扫描仪(如Velodyne)因其在获取距离数据方面的高速性和准确性而得到了广泛的应用。由于三维(three-dimensional,3D)距离扫描每帧返回的数据点数量庞大,而占用网格法可以实现数据的快速压缩和高效处理,因此该方法成为最常用的一种环境表示的方法。占用网格[7-9]是一种2.5D表示方法,它通过将三维点云映射到固定在地面上的二维网格来实现。Thrun等人提出了一种基于网格的方法[10],它根据网格单元内点云高度之间的最大绝对差将网格单元划分为地面单元和非地面单元。在2007年美国国防部高级研究计划局城市挑战赛中,为了分割点云和探测赛道上的其他车辆,许多团队都采用了这种方法,取得了巨大的成功。虽然这种表示方法的计算效率比较高,但仍存在如下一些问题。首先,受栅格大小影响,该方法容易出现欠分割,即障碍物与附近的地面一起被视为不可通行区域。其次,由于将三维点云映射到2.5D网格结构中,引入了降维技术,而并没有建模任意两点之间的自由垂直空间,造成该方法并不能识别悬挂障碍物。例如,树顶和树下的地面由于在同一栅格中而将被视为不可通行的物体,这给点云分割带来了新的挑战。针对悬挂物的处理,文[11]提出了一种扩展高度图的方法。除此之外,由P.Pfaff等人提出的多尺度高度图[12]是用于处理悬挂和垂直结构的高度图的一种变种,它由两维栅格组成,每一个栅格包含一系列平面块。但是该方法由于保存的数据较多,处理的效率并不高。
除了上述基于栅格的方法,对地面和障碍物点的分割还包括基于模型的方法[13-16]。Himmelsbach等人[17,18]将三维点云映射到极坐标栅格地图中,并从中选取高度最低的点作为候选地面点。利用候选点的分布在柱坐标系中进行直线拟合,其中满足一定坡度阈值的直线段上的点被认为是地面点。Douillard等人[19]将64线激光扫描仪数据投影到栅格地图中,采用二维高斯过程回归算法直接对整个地图进行地面拟合和障碍物的提取,该算法可以用于描述起伏路面,但是由于该算法计算比较复杂,只能获得近似实时的效果。Chen等人[20]使用非平稳协方差函数确定适应地形起伏的地面点,使其对粗糙环境具有鲁棒性。然而,在检测较小的障碍物时会出现错误的分类结果,例如,将台阶视为可通行地面而将上坡路视为障碍物。除此之外,多个文献[21,22]中都用到了RANSAC的方法来提取平面,但该方法仍然欠缺处理悬挂物体的能力。文[23]由多个不相交平面区域拟合不完全平面的地面。该方法在实时计算复杂度较高的RANSAC方案下可以近似地找到最优的区域划分和平面假设。
虽然上述方法取得了很大的进展,但针对野外环境下的点云分割的研究还远未成熟。基于网格的方法虽然可以实现点云的快速分割,但在具有悬挂物的环境中难以准确区分地面点和悬挂障碍物点,造成分割的精度下降。而基于模型的方法难以适应特征多变的野外环境。总结文献可以得知,现有方法在点云分割的精度上存在欠缺,影响移动机器人可通行域的提取和导航任务的实施。
发明内容
本发明的目的是为了克服现有技术中的不足,提供一种基于欧式聚类的多平面提取方法。
本发明的目的是通过以下技术方案实现的:
一种基于欧式聚类的多分辨率混合高度图方法,包括以下步骤:
(1)采用极坐标栅格表示地图;
通过二维极坐标栅格将地图按如下方式进行划分:以机器人所在位置为中心、前进方向为起始边界线,按逆时针方向将地图均匀划分为m个扇形块,对于每一个扇形块si,沿半径方向等间距划分为n个子块;将机器人上三维激光扫描仪得到的三维点云数据投影到二维极坐标栅格;
(2)遍历地图中的所有栅格,运用基于欧式聚类分割的混合高度图方法对栅格属性进行初步判断;具体如下:
当某一栅格的高度差小于给定阈值T1时,将该栅格中所有的点云都加入到准地面点集合中;如果栅格的高度差大于给定阈值T1,则判定对应栅格为不确定栅格;
通过欧式聚类分割方法对不确定栅格中的点云进行分组并计算每组中点云的高度平均值,并从各组点云中选取高度平均值最小的一组点云作为基准组,其他组作为比较组;然后,将各个比较组中点云高度最小的值与基准组中点云高度最大的值进行比较;如果所有比较组和基准组之间的高度差值均大于可通行间隙阈值T2,则将基准组中的点云加入到准地面点集合中,而将其余所有的比较组中的点云加入到障碍物点集合中;
(3)通过随机采样一致性的多平面提取方法对基于欧式聚类分割方法得到的准地面点进行再次滤波处理,具体如下:
(301)将基于欧式聚类分割方法得到的准地面点作为随机采样一致性方法的输入测量点序列进行平面模型提取,获取平面参数并建立平面模型如式(3)所示:
ax+by+cz+d=0 (3)
得到拟合平面的法向量为:
Figure BDA0002717226900000031
以机器人车体坐标系z轴方向的单位向量
Figure BDA0002717226900000032
作为参考法向量;计算拟合平面的法向量
Figure BDA0002717226900000037
与参考法向量
Figure BDA0002717226900000033
之间的夹角θ,夹角θ的计算方式如式(4)所示:
Figure BDA0002717226900000034
若夹角θ小于最大坡度阈值θth,则将拟合平面的内点加入最终地面点集合中;
(302)而拟合平面的外点则继续作为输入测量点序列进行平面模型提取,如果满足约束条件则按(301)中步骤处理,否则将本次拟合平面的内点加入障碍物点集中;
(303)剩余的外点则按照步骤(302)重复进行操作,直到剩余的外点数量小于设定的阈值Tnum
进一步的,步骤(1)中,假设当前时刻t下三维激光扫描仪得到的无序的三维点云数据为Pt={p1,…,pm,…,pN},即共包含N个激光数据点;将上述三维点云数据投影到二维极坐标栅格中,对每一个激光数据点pm=(xm,ym,zm)T,其映射方式如下:
Figure BDA0002717226900000035
Figure BDA0002717226900000036
其中,θ是三维数据点pm在投影平面内与X轴正方向之间的夹角,大小范围为[0,2π];Δth和Δr分别表示角度分辨率和半径分辨率,二者共同决定栅格的大小;而floor是一个向下取整函数,用于计算投影栅格在地图中的索引位置(si,bj);
进一步的,步骤(2)中,首先在低分辨率栅格下运用最大最小高度差法进行判断,如果满足对应栅格的高度差小于给定阈值T1则直接将该栅格内的点云加入到准地面点集合中,否则转入高分辨率的栅格下重新进行判断;如果直至最高分辨率层,栅格内点云的最大最小高度差值仍然大于给定阈值T1,则判定对应栅格为不确定栅格;对于不确定栅格则按照步骤(2)所提的基于欧式聚类的混合高度图方法进一步确定栅格属性。
与现有技术相比,本发明的技术方案所带来的有益效果是:
本发明方法克服了之前方法无法处理悬挂物的困难,运用所提出的基于欧式聚类的混合高度图方法将点云进行了分组,实现了地面点云和悬挂物点云的分离,提高了点云的分割精度,达到了很好的点云分割效果。而且该方法在点云稀疏的情况下也可以采用,这不同于文献[25]只能适用于密集点云的处理。除此之外,上述方法并没有以通常的单分辨率地图作为地图表达,而是采用了新的多分辨率地图的方式,这在一定程度上解决了由于栅格大小造成的欠分割与时间复杂度之间的矛盾的问题。最后,本发明还采用了一种新的基于随机采样一致性的多平面提取方法,通过迭代的方式自适应地提取符合地形约束条件的平面点,可以用于复杂环境(尤其是具有起伏的不平坦地形)中的点云分割。该方法是根据地势自适应的,能容忍地面的不平整。该方法解决了目前大多数主要针对城市以及乡村平坦道路的点云分割而在富有挑战的野外崎岖地形下仍然存在技术空缺的问题。通过在搭建的虚拟环境和实验环境中与多种已有算法进行比较,该算法在点云分割的精度和效率方面均有明显提高。
附图说明
图1为本发明整体算法流程示意图;
图2为极坐标栅格表示地图。通过将三维点云投影到二维极坐标栅格进行环境表示;
图3a至图3d为基于欧式聚类分割的混合高度图算法示意图。其中图3a表示运用最大最小高度图法对栅格属性进行初步判断,图3b用来计算悬挂障碍物与位于其下方的地面之间的距离,图3c表示经过欧式聚类分割后可能只存在一组点云的特殊情况,图3d用来判断地面处是否存在阻碍机器人通过的障碍物;
图4为基于多分辨率的极坐标栅格地图。
图5为采用随机采样一致性方法提取平面特征示意图。
图6a至图6c为基于随机采样一致性的多平面提取算法过程示意图。其中图6a表示基于欧式聚类的混合高度图方法存在的一些误检情况,图6b表示通过RANSAC方法进行滤波处理,图6c表示通过基于RANSAC的多平面提取方法进行滤波处理;
图7a至图7d为husky移动平台在gazebo仿真环境中的仿真结果。其中图7a表示第一类仿真环境示意图;图7b至图7d分别表示最大最小高度差方法、均值高度图方法以及本发明所提的基于欧式聚类的混合高度图方法;
图8a至图8d为summit-XL机器人在校园环境中的实验结果。其中图8a表示第一类实验环境示意图;图8b至图8d分别表示最大最小高度差方法、均值高度图方法以及本发明所提的基于欧式聚类的混合高度图方法;
图9a至图9c为husky移动平台在gazebo仿真环境中的仿真结果。其中图9a表示第二类仿真环境示意图;图9b和图9c分别表示文献[17]所用的三段拟合方法以及本发明所提的基于RANSAC的自适应多平面提取方法;
图10a至图10c为summit-XL机器人在校园环境中的实验结果。其中图10a表示第二类实验环境示意图;图10b和图10c分别表示文献[17]所用的三段拟合方法以及本发明所提的基于RANSAC的自适应多平面提取方法;
图11a和图11b关于单分辨率栅格地图和多分辨率栅格地图在时间性能上的对比示意图。
注:说明书附图中所涉及的点云图中,浅色点均代表地面点,深色点均代表障碍物点。
具体实施方式
以下结合附图和具体实施例对本发明作进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
本发明提供一种基于欧式聚类的多平面点云分割方法。通过将三维点云投影到二维极坐标栅格中,利用栅格内点云之间的高度差初步识别出可通行区域。对于不满足高度差约束条件的栅格,采用基于欧式聚类的方法进一步检查是否存在自由垂直空间。即通过将栅格内的点云进行分组,从中选出属于地面的部分,而位于上方的悬挂物部分与地面之间的高度决定了该栅格的可通行性。为了提高点云分割的效率,本发明采用了多分辨率栅格地图的方法,可以有效降低计算存储。除此之外,通过迭代地提取具有平面特征的点,基于随机采样一致性的多平面提取方法可以很好的适应不平坦的地面。整体算法流程如附图1所示。通过虚拟环境仿真和真实校园环境的实验结果表明,相比于现有方法,本发明具有更高的点云分割精度和准确度。
具体的,本发明提供的基于欧式聚类的多平面点云分割方法包括:
步骤一:采用极坐标栅格表示地图;
三维激光扫描仪的单帧数据通常包含几万甚至几十万的测量点,如果直接对这些测量点进行操作会产生巨大的时间消耗。所以在对点云进行分割之前,通常需要对数据进行合理的表示以提高分割效率。基于栅格地图的表示是最常用的一种三维数据表示方法。多线激光扫描仪的数据特性表现为近处稠密、远处较为稀疏。为了适应上述特点,本实施例采用极坐标栅格的形式来表示地图。该方法在一定程度上可以克服随着距离增加,激光扫描仪的测量点变得越来越稀疏所带来的数据分布不均匀的问题。极坐标栅格将地图按如下方式进行划分:以机器人所在位置为中心、前进方向为起始边界线,按逆时针方向将地图均匀划分为m个扇形块,对于每一个扇形块si,沿半径方向等间距划分为n个子块,具体划分结果如附图2所示。
假设当前时刻t下三维激光扫描仪得到的无序点云序列为Pt={p1,…,pm,…,pN},即共包含N个激光数据点。将这些三维点云投影到二维极坐标栅格中,对每一个激光数据点pm=(xm,ym,zm)T,其映射方式如下:
Figure BDA0002717226900000061
Figure BDA0002717226900000062
其中,θ是三维数据点pm在投影平面内与X轴正方向之间的夹角,其大小范围为[0,2π]。Δth和Δr分别表示角度分辨率和半径分辨率,二者共同决定栅格的大小。而floor是一个向下取整函数,用于计算投影栅格在地图中的索引位置(si,bj)。
步骤二,基于欧式聚类分割的混合高度图方法;
基于栅格的方法只保留了每个栅格中的部分数据信息,这使得需要处理的数据量变少,处理效率变高。目前,该方法成为最常用的一种三维点云分割方法。基于栅格的方法根据每个栅格中所包含的数据信息的不同,又可以分为均值高度图、最大最小值高度图等。其中,均值高度图的每个栅格中只存储一个值,即该地形位置的高度信息平均值。通过将该平均值与阈值进行比较,来判断每个栅格的状态是障碍物还是非障碍物。而在最大最小值高度图中,每个栅格仅包含两个值——所有投影到同一栅格中的测量点高度的最大值和最小值。栅格高度差法是通过计算最大高度值与最小高度值之间的差值,并将其与设定的阈值进行比较,如果大于阈值,则认为该栅格内的点均属于障碍物点,反之则判断为地面点。尽管上述方法在点云分割技术中被广泛应用,并在一定程度上表现出了良好的分割效果。但是他们都不具备处理悬挂结构(树冠、隧道、矿洞等)的能力。
算法1基于欧式聚类分割的混合高度图方法
Figure BDA0002717226900000071
因此该步骤基于欧式聚类分割的混合高度图方法区分地面点和悬挂物点,实现对地面上方障碍物的准确检测和下方对应可通行区域的正确提取。欧式聚类分割的混合高度图方法的流程图如算法1所示。
首先遍历极坐标地图中的所有栅格,运用最大最小高度图法对栅格属性进行初步判断(如附图3a所示)。当某一栅格的高度差小于给定阈值T1时,将该栅格中所有的点云都加入到准地面点集合中;如果栅格的高度差大于给定阈值T1,则判定为不确定栅格,需要对该不确定栅格的属性进行进一步的判断。由于悬挂物体与位于其下方的地面之间存在明显的间距,所以本发明根据这一特性来判断是否存在使机器人从中安全通行的垂直自由空间。在算法1的第6行,通过欧式聚类分割方法[24]对“不确定栅格”中的点云进行分组。函数Group()用于计算每组中点云的高度平均值,并从中选取平均值最小的一组点云作为基准组,其他组作为比较组。然后,将各个比较组中点云高度最小的值与基准组中点云高度最大的值进行比较。如果所有比较组和基准组之间的高度差值均大于可通行间隙阈值T2,则将基准组中的点云加入到准地面点集合中,而将其余所有的比较组中的点云加入到障碍物点集中。在算法1的第11-16行中,本实施例采用了最大最小高度图方法,用来计算悬挂障碍物与位于其下方的地面之间的距离,如附图3b所示。
但是,在使用上述基于欧式聚类的混合高度图方法对栅格属性进行进一步判断时,也会出现一些特殊情况。例如,经过欧式聚类分割后可能只存在一组点云,这说明该栅格中点云分布较为密集。出现这种情况的原因可能是因为多线激光扫描仪并未同时从同一地点的地面位置和其上方的悬挂物位置获得激光测量点。此时,将这一组点云全部划分到障碍物点集中,原因如下:(1)在采用最大最小高度差法进行初步判断时,这组点云的高度差大于阈值T1。如果其由地面处返回,则说明该位置存在地面障碍物。(2)如果该组点云由悬挂物位置返回,则其同样属于障碍物点集(如附图3c所示)。除此之外,还会出现如下一类情况,即栅格内地面与悬挂物之间的距离大于阈值T2,但是地面处存在阻碍机器人通过的障碍物。对此,除了需要计算对比组与基准组之间的距离之外,还需要对基准组的点云进行检测。即在基准组内部再次采用最大最小高度差法判断是否存在地面障碍物(如附图3d所示)。只要上述测得的两个距离中有一个不满足阈值条件,则该栅格内的点云就全部划分到障碍物点集中,见算法1中第14-15行。本实施例中,阈值T1设置为0.2m,T2设置为1m。
本发明在欧式聚类分割方法的基础上结合了多种栅格高度图法,实现了在具有悬挂结构的环境中点云的正确分割。该方法不仅可以有效区分地面部分和悬挂物部分,而且对地面上的障碍物也具有准确的检测识别能力。
步骤三:多分辨率栅格地图法;
基于栅格地图的方法虽然在一定程度上可以实现对点云的初步分割,但是栅格大小的选择会直接影响点云分割的精度。显然,栅格越大,就会出现欠分割的情况,影响点云分割的准确度;但过小的栅格又会增加时间耗费,影响判断的实时性。由此可以发现,点云分割的准确度和其所花费的时间代价在单一分辨率下是矛盾的。
因此,本实施例还提供了一种基于多分辨率的极坐标栅格地图方法(如附图4所示)。即首先在低分辨率栅格下运用最大最小高度差法进行判断,如果满足阈值条件T1则直接将该栅格内的点云加入到准地面点集合中,否则继续划分该栅格转入较高分辨率的栅格下重新进行判断。如果直至最高分辨率层,栅格内点云的最大最小高度差值仍然大于阈值,则转入步骤二中所提出的基于欧式聚类分割的混合高度图方法进行判断。
当栅格较大时,一些平坦的区域可以直接被判断为可通行地面,这不仅不会影响点云分割的准确度而且减少了不必要的时间耗费;而当栅格较小时,又能够对一些障碍物区域实现精细划分,大大降低出现欠分割的概率。然而,不同层级的分辨率的选择也很重要。本文选定角度分辨率为2.5°/半径分辨率为0.25m作为栅格的最高分辨率,而角度分辨率为20°/半径分辨率为2m作为栅格的最低分辨率。表1所示为各层栅格的分辨率大小。
表1栅格各层分辨率大小
Figure BDA0002717226900000091
步骤四,基于随机采样一致性方法(Random Sample Consensus,RANSAC)的自适应多平面提取方法;
步骤一中的点云分割的方法,不仅可以处理具有悬挂结构的物体,而且避免了欠分割的问题。但是上述基于欧式聚类分割的多分辨率混合高度图方法仍存在如下问题:(1)这种基于栅格的方法只反映局部特性,对于环境中出现的高空水平物体,由于局部的高度未发生明显变化,所以很容易将其划分为地面点,出现错分的情况;(2)栅格高度图法只体现垂直高度上的变化,并不能反映水平特性,这会导致一些具有较大倾角的坡面由于也具有平面特征而被错误地划分为可行驶路面。如附图6a所示,实线框所标记的点云均出现错误分割。针对上述情况,本实施例采用随机采样一致性方法对基于欧式聚类分割方法得到的准地面点进行再次滤波处理。
随机采样一致性方法(Random Sample Consensus,RANSAC)是一种模型提取方法,它采用迭代的方式从一组包含外点的观测数据集中估计出数学模型的参数。RANSAC的基本假设是:
a.观测数据集由符合某种数学模型分布的内点测量序列和不符合该分布的外点测量序列组成;
b.给定一组观测数据集,存在一种模型用来描述其中的内点,而外点作为噪声不能适应该模型。
给定一组包括内点和外点的观测序列,RANSAC方法利用一种投票程序逐步搜索到最适合给定数学模型的参数。其迭代过程包含如下两个步骤:首先假定计算给出数学模型参数的所需最小样本数量为n,则从观测数据集中随机挑选出样本数量为n的一组子序列,基于该子序列计算所给数学模型的参数;其次遍历整个观测序列,判断样本是否符合上步中得到的模型,记录符合该模型的内点数量。
重复执行以上两个步骤直到迭代次数达到所给阈值,最后从中找到符合模型的内点数量最多的一个模型并记录该模型参数。
本实施例采用RANSAC方法提取平面特征,筛选得到符合移动机器人可通行要求的地面点云,如附图5所示。首先,将第1节中得到的“准地面点”作为一组观测序列,基于该观测数据集进行平面模型提取。由于RANSAC方法得到的是内点数量最多的一个平面模型,所以该平面一般为以机器人所处位置为中心的主平面。经过滤波处理,图附6b中实线框所标注的点云被正确地划分为障碍物点。但同时,一些可通行的倾角较小的坡面点由于和主平面点距离较远而被作为噪声滤除,如附图6b虚线框标注所示。
针对上述问题,本实施例基于RANSAC提出了一种自适应的多平面提取方法,通过迭代的方式提取满足一定坡度约束的平面。首先,基于步骤二中得到的准地面点,获取平面参数并建立平面模型如式(3)所示:
ax+by+cz+d=0 (3)
平面的法向量为:
Figure BDA0002717226900000101
然后以车身z轴方向的单位向量
Figure BDA0002717226900000102
作为参考法向量。则拟合的地平面应该满足一定的坡度约束,即拟合平面的法向量
Figure BDA0002717226900000103
与参考法向量
Figure BDA0002717226900000104
之间的夹角θ应该小于最大坡度阈值θth。夹角θ的计算方式如式(4)所示:
Figure BDA0002717226900000105
基于随机采样一致性的自适应多平面提取方法的实现步骤如下:
(401)将准地面点作为随机采样一致性方法的观测序列进行平面模型提取,得到拟合平面的法向量
Figure BDA0002717226900000111
利用式(4)计算地平面法向量与参考平面法向量之间的夹角,若小于最大坡度阈值θth(本实施例中,阈值θth设置为15°),则将拟合平面的内点加入最终地面点集合中。
(402)而拟合平面的外点则继续作为观测序列进行平面模型提取,如果满足约束条件则按(1)中步骤处理,否则将本次拟合平面的内点加入障碍物点集中。
(403)剩余的外点则按照步骤(402)重复进行操作,直到剩余的外点数量小于设定的阈值Tnum(本实施例中,阈值Tnum设置为0.1N)。
基于RANSAC的自适应多平面提取方法不仅可以滤除局部噪点,还对坡度有一定的约束作用。如附图6c所示,实线框中标注的点云均被正确划分为障碍物点,而虚线框中的坡面点也由于多次迭代提取平面而被正确划分为可通行地面点,证明了该算法在起伏地面的点云分割中具有良好的性能。
具体的,为了验证本文所提方法的性能,本发明分别在虚拟仿真环境和校园实际环境中做了实验。一方面,通过机器人仿真平台Gazebo搭建不同类型的野外环境,包括树冠,隧道,山坡,岩石等。采用室外移动平台husky来模拟野外机器人的作业,可以实现真实的动力学仿真。另一方面,在校园环境中也选择了合适的场景来说明算法性能,分别为树下的人行道和湖边的草地。移动机器人平台则选用能够应用于室内外各种复杂地形的summit-XL机器人。除此之外,两种环境下所选择的三维激光扫描仪类型相同,均为Velodyne VLP-16激光扫描仪。
本发明中所提出的基于欧式聚类的混合高度图方法与最大最小高度差方法和均值高度图方法进行了比较,上述三种算法都结合了RANSAC方法进行滤波。该对比实验可以充分说明基于欧式聚类的混合高度图方法在处理悬挂障碍物方面具有良好表现。附图7a至附图8d分别是husky移动平台在gazebo仿真环境中和summit-XL机器人在校园环境中的实验结果。由图可知,最大最小高度差方法和均值高度图方法均能正确的检测出地面障碍物,但却无法处理悬挂物。上述两种方法都不能有效识别出位于树冠下的可通行区域,这将严重影响移动机器人的行进过程。图中实线框所标记的点云均属于被误检为障碍物的部分。而本发明方法无论在树冠还是隧道中,在处理悬挂物时均表现出了明显的优势。显然,本发明的方法有效的解决了上述难题,增强了路面的可通行性。
上述的实验结果表明,基于欧式聚类的混合高度图方法在处理悬挂物方面具有良好的性能,进一步通过实验验证本发明方法在不平坦路面上的性能。将本发明中所提的基于随机采样一致性的多平面提取方法与文献[17]中的提取路面的方法进行比较。文献[17]通过将初始点云沿着机器人行进方向分割成多个分段,并对每个分段拟合一个平面,从而纠正被错误标注的地面点。如附图9b和9c所示,不同分段所拟合的平面用不同灰度标注。而基于随机采样一致性的多平面提取方法是按照地形自适应的分成多个分段,并对每个分段拟合一个平面。两种方法的对比结果如附图9a至图9c和附图10a至图10c所示,显然,基于随机采样一致性的多平面提取方法分割点云的精度更高,具有更好的效果。
关于时间性能;为了说明多分辨率栅格地图方法在减少时间消耗方面的优势,本实验方案在上述四种环境下进行时间性能测试。点云分别在单分辨率栅格地图和多分辨率栅格地图下进行分割测试。附图11a和图11b显示了两种栅格分辨率之间的比较结果。由于第一类环境的复杂度较高,两种分辨率下的运行时间几乎接近。但是在第二类环境中,多分辨率下的点云分割只需要20ms,本发明方法在时间性能上提高了接近20%。
参考文献
[1].Quan Q,Jianda H.2.5-dimensional angle potential field algorithmfor the real-time autonomous navigation of outdoor mobile robots.ScienceChina Information Sciences,2011,54(10):2100-2112.
[2].Satish K R,Prabir K P.Computing an unevenness field from 3D laserrange data to obtain traversable region around a mobile robot.Robotics andAutonomous Systems,2016,84:48-63.
[3].Suger B,Steder B,Burgard W.Traversability analysis for mobilerobots in outdoor environments:A semi-supervised learning approach based on3D-lidar data.Proceedings IEEE International Conference on Robotics andAutomation,2015:3941-3946.
[4].Zhu B,Xiong G,Di H,et al.A Novel Method of Traversable AreaExtraction Fused With LiDAR Odometry in Off-road Environment.2019IEEEInternational Conference on Vehicular Electronics and Safety(ICVES),2019.
[5].Papadakis P.Terrain traversability analysis methods for unmannedground vehicles:A survey.Engineering Applications of Artificial Intelligence,2013,26(4):1373-1385.
[6].Byun J,Na K I,Seo B S,et al.Drivable Road Detection with 3D PointClouds Based on the MRF for Intelligent Vehicle.Springer Tracts in AdvancedRobotics,2015,105:49-60.
[7].Himmelsbach M,Luettel T,Wuensche H J.Real-time objectclassification in 3D point clouds using point feature histograms.IEEE/RSJInternational Conference on Intelligent Robots and Systems,2009:994-1000.
[8].Igor B,Cyrill S.Efficient Online Segmentation for Sparse 3D LaserScans.Photogrammetrie Fernerkundung Geoinformation,2017,85(1):41-52.
[9].Reina A J,Martinez J L,Mandow A,et al.Collapsible cubes:Removingoverhangs from 3D point clouds to build local navigable elevation maps.IEEE/ASME International Conference on Advanced Intelligent Mechatronics,2014.
[10].Thrun S,Montemerlo M,Dahlkamp H,et al.Stanley:The Robot that Wonthe DARPA Grand Challenge.Journal of Field Robotics,2006,23(9):661-692.
[11].Pfaff P,Triebel R,Burgard W.An Efficient Extension to ElevationMaps for Outdoor Terrain Mapping and Loop Closing.International Journal ofRobotics Research,2007,26(2):217-230.
[12].Triebel R,Pfaff P,Burgard W.Multi-Level Surface Maps for OutdoorTerrain Mapping and Loop Closing.IEEE/RSJ International Conference onIntelligent Robots and Systems,2006:2276-2282.
[13].Liu B,Chen X,Han Y,et al.Accelerating DNN-based 3D point cloudprocessing for mobile computing.Science China Information Sciences,2019,62(11):36-46.
[14].Chu P M,Cho S,Park Y W,et al.Fast point cloud segmentation basedon flood-fill algorithm.2017 IEEE International Conference on MultisensorFusion and Integration for Intelligent Systems(MFI),2017.
[15].Wen C,Sun X,Li J,et al.A deep learning framework for roadmarking extraction,classification and completion from mobile laser scanningpoint clouds.ISPRS Journal of Photogrammetry and Remote Sensing,2019,147(1):178-192.
[16].Gao B,Xu A,Pan Y,et al.Off-Road Drivable Area Extraction Using3D LiDAR Data.2019 IEEE Intelligent Vehicles Symposium(IV),2019.
[17].Zermas D,Izzat I,Papanikolopoulos N.Fast segmentation of 3Dpoint clouds:A paradigm on LiDAR data for autonomous vehicleapplications.IEEE International Conference on Robotics and Automation,2017.
[18].Himmelsbach M,Hundelshausen F V,Wuensche H J.Fast segmentationof 3D point clouds for ground vehicles.IEEE Intelligent Vehicles Symposium,2010:560-565.
[19].Douillard B,Underwood J,Kuntz N,et al.On the segmentation of 3DLIDAR point clouds.IEEE International Conference on Robotics and Automation,2011:2798-2805.
[20].Chen T,Dai B,Wang R,et al.Gaussian-Process-Based Real-TimeGround Segmentation for Autonomous Land Vehicles.Journal of Intelligent andRobotic Systems,2014,76(3-4):563-582.
[21].
Figure BDA0002717226900000141
A,Hermans A,Engelmann F,et al.Multi-Scale Object Candidatesfor Generic Object Tracking in Street Scenes.IEEE International Conference onRobotics and Automation,2016:3180-3187.
[22].Meng X,Cao Z,Liang S,et al.A terrain description method fortraversability analysis based on elevation grid map.International Journal ofAdvanced Robotic Systems,2018.
[23].Li B.On Enhancing Ground Surface Detection from Sparse LidarPoint Cloud.IEEE/RSJ International Conference on Intelligent Robots andSystems,2019.
[24].Spark D N.Euclidean cluster analysis algorithm.AppliedStatistics,1973,22(1):126-130.
[25].Douillard B,Underwood J,Melkumyan N,et al.Hybrid elevation maps:3D surface models for segmentation.IEEE/RSJ International Conference onIntelligent Robots and Systems,2010.
本发明并不限于上文描述的实施方式。以上对具体实施方式的描述旨在描述和说明本发明的技术方案,上述的具体实施方式仅仅是示意性的,并不是限制性的。在不脱离本发明宗旨和权利要求所保护的范围情况下,本领域的普通技术人员在本发明的启示下还可做出很多形式的具体变换,这些均属于本发明的保护范围之内。

Claims (3)

1.一种基于欧式聚类的多平面提取方法,其特征在于,包括以下步骤:
(1)采用极坐标栅格表示地图;
通过二维极坐标栅格将地图按如下方式进行划分:以机器人所在位置为中心、前进方向为起始边界线,按逆时针方向将地图均匀划分为m个扇形块,对于每一个扇形块si,沿半径方向等间距划分为n个子块;将机器人上三维激光扫描仪得到的三维点云数据投影到二维极坐标栅格;
(2)遍历地图中的所有栅格,运用基于欧式聚类分割的混合高度图方法对栅格属性进行初步判断;具体如下:
当某一栅格的高度差小于给定阈值T1时,将该栅格中所有的点云都加入到准地面点集合中;如果栅格的高度差大于给定阈值T1,则判定对应栅格为不确定栅格;
通过欧式聚类分割方法对不确定栅格中的点云进行分组并计算每组中点云的高度平均值,并从各组点云中选取高度平均值最小的一组点云作为基准组,其他组作为比较组;然后,将各个比较组中点云高度最小的值与基准组中点云高度最大的值进行比较;如果所有比较组和基准组之间的高度差值均大于可通行间隙阈值T2,则将基准组中的点云加入到准地面点集合中,而将其余所有的比较组中的点云加入到障碍物点集合中;
(3)通过随机采样一致性的多平面提取方法对基于欧式聚类分割的混合高度图方法得到的准地面点进行再次滤波处理,具体如下:
(301)针对基于欧式聚类分割方法得到的准地面点,获取平面参数并建立平面模型如式(3)所示:
ax+by+cz+d=0 (3)
平面的法向量为:
Figure FDA0002717226890000011
以机器人车体坐标系z轴方向的单位向量
Figure FDA0002717226890000012
作为参考法向量;拟合平面的法向量
Figure FDA0002717226890000013
与参考法向量
Figure FDA0002717226890000014
之间的夹角θ小于最大坡度阈值θth;夹角θ的计算方式如式(4)所示:
Figure FDA0002717226890000015
(302)将基于欧式聚类分割方法得到的准地面点作为随机采样一致性方法的输入测量点序列进行平面模型提取,得到拟合平面的法向量
Figure FDA0002717226890000016
利用式(4)计算地平面法向量与参考平面法向量之间的夹角,若小于最大坡度阈值θth,则将拟合平面的内点加入最终地面点集合中;
(303)而拟合平面的外点则继续作为输入测量点序列进行平面模型提取,如果满足约束条件则按(301)中步骤处理,否则将本次拟合平面的内点加入障碍物点集中;
(304)剩余的外点则按照步骤(302)重复进行操作,直到剩余的外点数量小于设定的阈值Tnum
2.根据权利要求1所述一种基于欧式聚类的多平面提取方法,其特征在于,步骤(1)中,假设当前时刻t下三维激光扫描仪得到的无序的三维点云数据为Pt={p1,…,pm,…,pN},即共包含N个激光数据点;将上述三维点云数据投影到二维极坐标栅格中,对每一个激光数据点pm=(xm,ym,zm)T,其映射方式如下:
Figure FDA0002717226890000021
Figure FDA0002717226890000022
其中,θ是三维数据点pm在投影平面内与X轴正方向之间的夹角,大小范围为[0,2π];Δth和Δr分别表示角度分辨率和半径分辨率,二者共同决定栅格的大小;而floor是一个向下取整函数,用于计算投影栅格在地图中的索引位置(si,bj)。
3.根据权利要求1所述一种基于欧式聚类的多平面提取方法,其特征在于,步骤(2)中,首先在低分辨率栅格下运用最大最小高度差法进行判断,如果满足对应栅格的高度差小于给定阈值T1则直接将该栅格内的点云加入到准地面点集合中,否则转入高分辨率的栅格下重新进行判断;如果直至最高分辨率层,栅格内点云的最大最小高度差值仍然大于给定阈值T1,则判定对应栅格为不确定栅格。
CN202011078122.4A 2020-10-10 2020-10-10 一种基于欧式聚类的多平面提取方法 Active CN112184736B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011078122.4A CN112184736B (zh) 2020-10-10 2020-10-10 一种基于欧式聚类的多平面提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011078122.4A CN112184736B (zh) 2020-10-10 2020-10-10 一种基于欧式聚类的多平面提取方法

Publications (2)

Publication Number Publication Date
CN112184736A CN112184736A (zh) 2021-01-05
CN112184736B true CN112184736B (zh) 2022-11-11

Family

ID=73947346

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011078122.4A Active CN112184736B (zh) 2020-10-10 2020-10-10 一种基于欧式聚类的多平面提取方法

Country Status (1)

Country Link
CN (1) CN112184736B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP4258078A4 (en) * 2021-01-13 2023-12-27 Huawei Technologies Co., Ltd. POSITIONING METHOD AND DEVICE AND VEHICLE
WO2022213376A1 (zh) * 2021-04-09 2022-10-13 深圳市大疆创新科技有限公司 障碍物检测方法、装置、可移动平台及存储介质
US11741621B2 (en) 2021-05-10 2023-08-29 Qingdao Pico Technology Co., Ltd. Method and system for detecting plane information
CN113240678B (zh) * 2021-05-10 2023-05-30 青岛小鸟看看科技有限公司 平面信息检测方法及***
CN113276130B (zh) * 2021-05-28 2022-10-04 山东大学 一种基于点云切片的自由曲面喷涂路径规划方法及***
CN113920134B (zh) * 2021-09-27 2022-06-07 山东大学 一种基于多线激光雷达的斜坡地面点云分割方法及***
CN114119998B (zh) * 2021-12-01 2023-04-18 成都理工大学 一种车载点云地面点提取方法及存储介质

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101877128A (zh) * 2009-12-23 2010-11-03 中国科学院自动化研究所 一种三维场景中不同物体的分割方法
CN103226833A (zh) * 2013-05-08 2013-07-31 清华大学 一种基于三维激光雷达的点云数据分割方法
CN103268729A (zh) * 2013-05-22 2013-08-28 北京工业大学 基于混合特征的移动机器人级联式地图创建方法
CN110211169A (zh) * 2019-06-06 2019-09-06 上海黑塞智能科技有限公司 基于多尺度超像素和相位相关的窄基线视差的重构方法
CN110967024A (zh) * 2019-12-23 2020-04-07 苏州智加科技有限公司 可行驶区域的检测方法、装置、设备及存储介质
CN111292275A (zh) * 2019-12-26 2020-06-16 深圳一清创新科技有限公司 基于复杂地面的点云数据滤除方法、装置、计算机设备
CN111624622A (zh) * 2020-04-24 2020-09-04 库卡机器人(广东)有限公司 障碍物检测方法、装置

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105184852B (zh) * 2015-08-04 2018-01-30 百度在线网络技术(北京)有限公司 一种基于激光点云的城市道路识别方法及装置
CN107818288B (zh) * 2016-09-13 2019-04-09 腾讯科技(深圳)有限公司 标志牌信息获取方法及装置
GB2559157A (en) * 2017-01-27 2018-08-01 Ucl Business Plc Apparatus, method and system for alignment of 3D datasets

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101877128A (zh) * 2009-12-23 2010-11-03 中国科学院自动化研究所 一种三维场景中不同物体的分割方法
CN103226833A (zh) * 2013-05-08 2013-07-31 清华大学 一种基于三维激光雷达的点云数据分割方法
CN103268729A (zh) * 2013-05-22 2013-08-28 北京工业大学 基于混合特征的移动机器人级联式地图创建方法
CN110211169A (zh) * 2019-06-06 2019-09-06 上海黑塞智能科技有限公司 基于多尺度超像素和相位相关的窄基线视差的重构方法
CN110967024A (zh) * 2019-12-23 2020-04-07 苏州智加科技有限公司 可行驶区域的检测方法、装置、设备及存储介质
CN111292275A (zh) * 2019-12-26 2020-06-16 深圳一清创新科技有限公司 基于复杂地面的点云数据滤除方法、装置、计算机设备
CN111624622A (zh) * 2020-04-24 2020-09-04 库卡机器人(广东)有限公司 障碍物检测方法、装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"Fast Planar Detection System Using a GPU-Based 3D Hough Transform for LiDAR Point Clouds";Yifei Tian et al.;《Applied Science》;20200304;全文 *
"基于LiDAR/INS 的野外移动机器人组合导航方法";宋锐等;《智能***学报》;20200730;第15卷(第4期);全文 *

Also Published As

Publication number Publication date
CN112184736A (zh) 2021-01-05

Similar Documents

Publication Publication Date Title
CN112184736B (zh) 一种基于欧式聚类的多平面提取方法
Azim et al. Detection, classification and tracking of moving objects in a 3D environment
WO2018068653A1 (zh) 点云数据处理方法、装置及存储介质
CN103268729B (zh) 基于混合特征的移动机器人级联式地图创建方法
Weng et al. Pole-based real-time localization for autonomous driving in congested urban scenarios
Liu et al. Point cloud segmentation based on Euclidean clustering and multi-plane extraction in rugged field
CN109256028B (zh) 一种用于无人驾驶的高精度路网自动生成的方法
CN108845569A (zh) 生成三维高清道路图水平弯道行车线的半自动点云方法
CN109101743B (zh) 一种高精度路网模型的构建方法
CN113345008A (zh) 一种考虑轮式机器人位姿估计的激光雷达动态障碍物检测方法
Li et al. Robust localization for intelligent vehicles based on pole-like features using the point cloud
CN115620261A (zh) 基于多传感器的车辆环境感知方法、***、设备及介质
CN113936215A (zh) 一种矿区路面凹坑的识别方法、***及无人驾驶货车
Arora et al. Static map generation from 3D LiDAR point clouds exploiting ground segmentation
Xiong et al. Road-Model-Based road boundary extraction for high definition map via LIDAR
Sun et al. Oriented point sampling for plane detection in unorganized point clouds
Cai et al. A lightweight feature map creation method for intelligent vehicle localization in urban road environments
Ai et al. A real-time road boundary detection approach in surface mine based on meta random forest
Yan et al. RH-Map: Online Map Construction Framework of Dynamic Object Removal Based on 3D Region-wise Hash Map Structure
CN113077473B (zh) 三维激光点云路面分割方法、***、计算机设备及介质
Li et al. Intersection detection and recognition for autonomous urban driving using a virtual cylindrical scanner
Hu et al. LiDAR-based road extraction for UGV in high definition map
CN115792958A (zh) 一种基于3d激光雷达的无人矿车障碍物检测方法
Rangan et al. Improved localization using visual features and maps for Autonomous Cars
Mi et al. Automatic road structure detection and vectorization Using Mls point clouds

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant