CN114623817B - 基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法 - Google Patents

基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法 Download PDF

Info

Publication number
CN114623817B
CN114623817B CN202210158652.2A CN202210158652A CN114623817B CN 114623817 B CN114623817 B CN 114623817B CN 202210158652 A CN202210158652 A CN 202210158652A CN 114623817 B CN114623817 B CN 114623817B
Authority
CN
China
Prior art keywords
frame
feature
bundles
characteristic
track
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
CN202210158652.2A
Other languages
English (en)
Other versions
CN114623817A (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.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
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 Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202210158652.2A priority Critical patent/CN114623817B/zh
Priority to PCT/CN2022/080346 priority patent/WO2023155258A1/zh
Publication of CN114623817A publication Critical patent/CN114623817A/zh
Application granted granted Critical
Publication of CN114623817B publication Critical patent/CN114623817B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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
    • 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/10Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/165Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments
    • G01C21/1656Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments with passive imaging devices, e.g. cameras
    • 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/10Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/18Stabilised platforms, e.g. by gyroscope
    • 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/10Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/183Compensation of inertial measurements, e.g. for temperature effects

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Automation & Control Theory (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Image Analysis (AREA)
  • Navigation (AREA)

Abstract

本发明公开了一种基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法,属于多传感器融合导航定位领域。传统的滤波方法会按时间推移不断的去除老旧的状态量,在退化运动的情况下,保留的状态量对应的帧束之间没有足够的视差,所以难以约束运动并导致漂移。本发明方法包括图像特征提取、基于关键帧的特征关联、滤波器初始化、基于IMU的状态推算、利用特征观测更新滤波器和基于关键帧的状态量管理这几个步骤,可以实时的对传感器的几何参数进行估计。本发明方法基于关键帧组织状态量,在退化运动时,这些关键帧束对应的状态量不会被去除,所以仍然可以保证良好的观测,避免了漂移,本发明是第一个支持自校准的基于关键帧的滑窗滤波方法。

Description

基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法
技术领域
本发明属于多传感器融合导航定位领域,具体涉及基于单目或双目相机与惯性测量单元(IMU)组合系统的位置、速度与姿态的推算和系统中传感器参数的实时自标定。
背景技术
如何实时估计设备在场景中的位置和姿态(即位姿),是增强现实(AR)、移动机器人、无人机等应用领域的核心问题之一。视觉惯性里程计(visual inertial odometry)利用设备上固定安装的相机和IMU采集的数据,能够实时追踪设备的位姿。该技术融合了计算机三维视觉和惯性航位推算,利用低成本的传感器,能够达到厘米级的实时定位精度。
很多AR相关的企业对该技术进行了研发。谷歌公司在2014年Tango项目中展示了含有鱼眼相机和IMU的手机可以通过融合这些传感器的信息来估计手机的位姿并进行场景的三维重建,证明了移动设备可以在不依赖卫星导航系统的情况下拥有位姿和空间感知能力。随后,谷歌公司协同手机厂商推出了基于Tango的智能手机。该手机通过配备的鱼眼相机和IMU,可以跟踪自身的三维运动轨迹。2017年,苹果公司推出了针对含相机和IMU的iOS手持设备的ARKit框架,用于开发增强现实的应用,其中的运动跟踪技术是一种典型的视觉惯性里程计[1],[2]。同年,谷歌公司推出了适用于Android手机的含视觉惯性里程计[3],[4]的ARCore框架。
视觉惯性里程计也广泛应用于无人机和自动驾驶等行业。大疆新型的无人机通过视觉惯性里程计技术[5]可以在复杂场景中精确的估计位姿,从而实现避障和穿越等一系列飞行动作。多种自动驾驶套件中使用了视觉惯性里程计方法辅助自动泊车。
目前的视觉惯性里程计方法仍然存在一些问题,比如如何应对瑕疵的参数标定和如何应对退化运动(如静止和近似纯旋转)。对于大量的消费级产品,低成本的传感器通常没有精确的标定,这些不精确的标定会明显的降低视觉惯性里程计的位姿估计精度。另一方面,这些设备(如手机)在实际使用中会经常经历退化运动,现有的实时里程计方法在这些情况下经常会产生明显的位姿漂移。
发明内容
为了解决上述两个问题,本发明提出了一种含自标定的基于关键帧滑窗滤波(Keyframe-based Sliding Window Filter,KSWF)的适用于相机和IMU组合系统的视觉惯性里程计方法。该方法基于扩展卡尔曼滤波(EKF)进行状态估计,所以具有很好的实时性。该方法能够实时利用场景特征标定(即自标定)相机和IMU的几何参数,并且基于关键帧的概念进行针对滑动窗口中变量的滤波,来解决退化运动。在下文描述中,我们将该滤波器估计的位置、速度、姿态以及传感器的参数等状态量总称为状态向量。该滤波器也估计了这些状态量的不确定度,即状态向量的协方差。为简化起见,位置和姿态简称为位姿,在同一时刻t的位姿和速度简称为导航状态量。本文中传感器泛指相机和IMU,设备指承载这些传感器的装置,如手机,滤波器指融合数据进行状态估计的程序。KSWF适用于由一个或两个相机和一个IMU固定组成的系统,基于滑窗滤波器融合他们采集的数据以估计该系统(或承载设备)的状态向量和协方差。KSWF可以实时的标定IMU的零偏、尺度因子、错位和加速度敏感性,以及相机的投影几何参数、外部参数、相机的时间延迟和卷帘快门效应。其特征关联前端从单目相机获取一系列帧或从双目相机获取一系列大致同步的帧对作为输入,为了简化,后文将单目的帧或双目的帧对统一称为帧束。该前端在帧束中提取特征点,并基于关键帧的概念在帧束间匹配特征点,从而得到特征点观测的轨迹(即特征轨迹),并根据当前帧束与历史帧束视野重叠的比例选择关键帧束。后端滤波器利用IMU数据进行位姿推算,然后利用特征轨迹中的观测更新滤波器的状态向量和协方差,最后基于关键帧的概念管理滤波器中的状态量以保证实时运算。
因为KSWF在状态向量中包含了相机和IMU的几何标定参数,所以可以实时估计传感器的标定参数,以应对标定不准确的情况。在有准确标定参数的情况下,状态向量中的标定参数可以被固定,以减少计算量。
本方法为了避免计算量随时间不断增加,基于关键帧的概念选取并去除多余的状态量。去除多余状态量的环节也叫边缘化。KSWF根据导航状态量是否与关键帧束相关联来确定其边缘化的顺序,以延长关键帧束存在的时长。这样在设备进行退化运动时,由于关键帧束此时基本不发生变化,所以该方法可以避免传统方法的位姿漂移。
本发明针对现有相机—IMU融合状态估计技术难以处理标定误差和退化运动的不足,提供一种含自标定的基于关键帧滑窗滤波的视觉惯性里程计方法。本发明适用的设备上固定安装有至少一个IMU(含三轴陀螺仪和三轴加速度计)和至少一个相机(例如宽视角相机或鱼眼相机),相机的快门方式可以为卷帘快门或全局快门。本发明假设该设备上的相机和IMU有过粗略的标定。这种设备的实例有智能手机、AR眼镜、快递机器人等。该设备通过N个相机采集图像帧束的序列(每个帧束含由N个相机拍摄的N张时间相近的图像,N≥1)和一个IMU采集3轴加速度计数据和3轴陀螺仪数据,然后利用本方法融合相机和IMU的数据来估计设备的状态和参数,以及不确定度,将位置、速度、姿态以及相机和IMU的参数总称为状态向量,这些状态向量的不确定度由状态向量的协方差来描述;为简化起见,在同一个时刻t的位姿和速度简称为导航状态量;
本发明的实现方法包含以下步骤:
步骤1,当第j帧束到来时,在其中的N张图像上提取特征点和对应的描述子,对于图像k,得到其特征点的集合
进一步的,步骤1中特征点的检测方法包括ORB和BRISK,生成描述子的方法包括BRISK、FREAK和基于深度学习的特征描述子生成方法。
步骤2,将帧束j与先前的K个关键帧束进行匹配,以扩充特征点的特征轨迹。两个帧束的匹配包含N对图像的匹配,对每对图像,根据特征描述子匹配两帧图像中的特征点。两个匹配的特征点对应一个场景中的地标点Lm,m=0,1,…,该地标点Lm在多帧图像(j=0,1,…)中的观测(/>其中u,v是特征点在第j帧束的第k张图像中的像素坐标)形成一个特征轨迹/>当前帧束j与K个关键帧束匹配完成后,根据帧束j中的特征匹配个数与特征点个数的比例决定是否把帧束j选为关键帧束。
接着,如果上一帧束j-1不是关键帧束,还将匹配帧束j和帧束j-1,以扩充特征点的特征轨迹。
对于N>1的情况,还会对帧束j中的N帧图像进行两两匹配,以扩充特征点的特征轨迹。
两帧图像间的匹配通过寻找Hamming距离相近的特征描述子来实现。搜索相近的特征描述子的实现方式包括暴力搜索、k最近邻或基于窗口的搜索。
步骤3,如果滑窗滤波器尚未初始化,将通过几个帧束(如两个关键帧束)中的特征轨迹和IMU数据初始化位姿、速度和传感器参数等状态向量以及他们的协方差。
步骤3中状态向量X(t)包含当前导航状态量π(t)、IMU参数ximu,N个相机的参数 过去的Nkf个关键帧束对应的导航状态量/>和最近Ntf+1个帧束对应的导航状态量/>以及m个地标点L1,L2,…,Lm;其中,ximu包括IMU的零偏、尺度因子、错位和加速度敏感性,xC包含每个相机的投影几何参数、外部参数、相机的时间延迟和卷帘快门效应;过去的Nkf个关键帧束和最近的Ntf个帧束形成一个导航状态量的滑动窗口。
步骤4,在滤波器完成初始化的情况下,对于完成特征匹配的当前帧束j,根据IMU数据,从上一帧束j-1对应的位姿和速度递推估算帧束j对应时刻tj的位姿和速度,并克隆推算的位姿和速度以扩充滤波器的状态向量和协方差。
步骤5,用帧束j中的多个特征轨迹作为观测信息,更新滤波器的状态向量和协方差。更新的状态向量和协方差可以发布给下游的应用,例如AR中虚拟场景的投影和机器人路径规划等。根据特征轨迹是否在当前帧束j中消失或者对应的地标点是否在状态向量中,确定特征轨迹的处理方式。进而,利用这些特征轨迹中的观测来更新状态向量和协方差。因为传感器的几何参数可以置于状态向量中以得到更新,所以本方法支持实时自标定。
具体而言,针对每条特征轨迹可能发生的三种情况,采取不同的方式处理。这三种情况包括该特征轨迹在当前帧束j中消失,该特征轨迹对应的地标点在状态向量中,以及该特征轨迹对应不在状态向量中但观测良好的地标点。对于第一种情况,首先用该特征轨迹三角化一个地标点,并计算特征轨迹中所有观测的重投影误差和雅可比矩阵,然后利用矩阵零空间投影以消除地标点参数的雅可比矩阵,再通过马氏检验去除异常值。对于第二种情况,首先计算特征轨迹在当前帧束j中观测的重投影误差和雅可比矩阵,然后通过马氏检验去除异常值。对于第三种情况,首先用该特征轨迹三角化一个地标点,然后计算特征轨迹中所有观测的重投影误差和雅可比矩阵,通过马氏检验去除异常值,接着添加该地标点到状态向量和协方差矩阵。最后,按以上三种观测情况把所有的特征轨迹分为三类,对每类特征轨迹,叠加上述的重投影误差和雅可比矩阵,对状态向量和协方差矩阵进行更新,更新的方式与经典EKF相同。。
步骤6,为了控制计算量,在更新完毕后,会检测当前状态向量中冗余的导航状态量。这些冗余的状态量根据其是否与关键帧束对应来选取,接着会被移出滤波器。
具体而言,当滤波器状态向量中导航状态量的数量超过允许的关键帧束数目Nkf和最近帧束数目Ntf之和时,则从中选取冗余帧束并将其边缘化。每次边缘化操作,至少选择Nr个冗余帧束(对于单目相机,Nr为3,对于双目相机,Nr为2)。为了达到这个要求,首先在最近的非关键帧束中选择冗余帧束,同时排除最近的Ntf个帧束,然后在最早的关键帧束中选择冗余帧束。在找到Nr个冗余帧束后,可以利用其中长度超过两次观测的特征轨迹来更新滤波器——即如果这样的一条特征轨迹可以利用其所有的观测成功三角化一个地标点,则使用该特征轨迹上冗余帧束中的特征点观测进行EKF更新,类似步骤5第一种情况;冗余帧束中属于其他特征轨迹的特征点将被丢弃。EKF更新后,删除这些冗余帧束对应的导航状态量和协方差矩阵的行和列。
对下一帧束j+1,将重复执行步骤1-6。每次循环都将发布滤波器估计的状态向量和协方差,以服务下游的应用。这些发布的信息描述了设备的位姿、速度和传感器参数等状态量。
与现有技术相比,本发明的优点和有益效果如下:
1.由于该里程计方法在滤波过程中可以估计IMU内参、相机内参、相机外参、时间延迟和卷帘快门效应,所以能够适用于初始标定不准确或使用卷帘快门的相机-IMU组合,例如手机上的相机-IMU组合。
2.该里程计方法采用了基于关键帧的特征关联和状态量管理,所以能够稳健的处理退化运动的情形,例如静止、悬停、近似纯旋转等。
附图说明
图1为基于关键帧滑窗滤波的视觉惯性里程计方法示意图,灰底模块运用了关键帧的概念。
图2为基于关键帧滑窗滤波的视觉惯性里程计方法的多线程实现示意图,灰底模块运用了关键帧的概念。
图3为基于关键帧的图像特征点匹配的步骤①-⑦的流程图,以双目相机为例。
图4为状态向量中所含状态量的两个实例示意图。
图5为针对图3的情况,对当前帧束j进行导航状态量克隆之前和之后的状态向量示意图,πj是新增的状态量。
图6为针对图3的情况,在对当前帧束j导航状态量克隆之前的协方差矩阵示意图。
图7为针对图3的情况,在对当前帧束j导航状态量克隆之后的协方差矩阵示意图,内填斜线的方格对应新增的行和列,注意内填斜线的行列与内填棋盘格的行列具有相等的元素。
图8为地标点特征轨迹用于状态向量和协方差更新的三种方式示意图,如果地标点在当前帧束消失且不在状态向量中,则投影以去除其雅可比系数矩阵,然后再用于滤波更新;如果地标点在状态向量中,则直接用于滤波更新;如果地标点在当前帧束被观测到并且可以被成功三角化,则添加到状态向量中,然后用于滤波更新。
图9为冗余帧束选取的方法示意图。图中显示了关键帧束(灰色底)和普通帧束(白色底)对应的导航状态量。注意在边缘化冗余帧束后需要保留至少一个关键帧束。
图10为冗余帧束选取的两个例子示意图,假设Ntf=3,Nr=2,灰色底标记了关键帧束,白色底标记了普通帧束,虚线标记了将被边缘化的冗余帧束。左图:针对图3的情况,对当前帧束完成滤波更新后,需要选取的冗余帧束。右图:针对另一种假设情况,需要选取的冗余帧束。
具体实施方式
本发明利用相机和IMU采集的数据,考虑低端传感器标定的残余误差和静止等退化运动导致的问题,提出了一种含自标定的基于关键帧的滑窗滤波器进行实时位姿估计。该滤波器适用于单目或双目相机和IMU的组合系统,相机的快门方式可以为卷帘快门或全局快门。通过实时估算传感器的几何标定参数,可以提高位姿估计的精度;通过运用关键帧的概念,可以避免退化运动时的位姿漂移现象。
本发明提供的方法能够用计算机软件技术实现,实施主要步骤分为:图像特征提取、基于关键帧的特征关联、滤波器初始化、基于IMU的状态推算、利用特征观测更新滤波器和基于关键帧的状态量管理,参见图1。实施例以双目相机和IMU组合系统为例对本发明的流程进行一个具体阐述,如下:
步骤1,对于由N个相机捕获的含N帧时间相近图像的帧束j,检测每张图像中的特征点,并对这些特征点生成描述子(即特征向量),特征点的检测方法包括ORB[6]或BRISK[7],生成描述子的方法有BRISK[7]或FREAK[8]。
实施例具体的实施过程说明如下(例如N=2,特征检测和描述方法都为BRISK):
对由2个相机捕获的帧束j,检测和提取图像特征。为了对每帧图像检测特征点,首先构建含4个图层和4个内图层的尺度空间,得到该帧图像8个不同采样率的图像。然后利用FAST[9]检测算子定位特征点,并在多尺度空间中对定位的特征点进行非极大值抑制,即每个特征点与该特征点的图像空间(8邻域)和尺度空间(上下两层9×2=18邻域)的共26个邻域点比较,保留在邻域中FAST响应值取得最大值的特征点,并去除邻域中的非极大值。最后,基于二次拟合得到特征点的亚像素级精确位置和尺度空间的精确尺度。为了描述每个特征点,首先以该特征点为中心,构建不同半径的同心圆,在每个圆上获取一定数目的等间隔采样点,这些采样点包括特征点本身共有60个,对所有采样点进行高斯滤波以消除混叠效应,并计算这些采样点形成的灰度主方向;然后旋转该特征点周围的采样区域至主方向,得到新的采样区域,并构建512个采样点对,根据点对的灰度差异形成512维的二进制编码,即每个特征向量为64Bytes。
步骤2,利用步骤1提取的图像特征描述子进行基于关键帧的特征关联(也称为特征匹配)。如图3所示,匹配过程将当前帧束j的特征与几个先前帧束的特征相关联,包括K个最近的关键帧束和上一帧束j-1。在与K关键帧束的匹配完成后,如果当前帧束中特征匹配个数与特征点个数的比例<Tr,则将帧束j选为关键帧束。如果有多个相机,也可以在当前帧束j的N张图像之间进行特征关联。两帧图像之间特征描述子的匹配通过暴力搜索、k最近邻或窗口搜索的方法完成。
实施例具体的实施过程说明如下(例如N=2,K=2,暴力搜索,Tr=0.2):
图3展示了对于一个含双目相机的设备采集的数据,进行特征点匹配的过程,其中灰底的帧束表示关键帧束,白底帧束表示普通帧束。特征匹配包含①-⑦等步骤,①是当前帧束j与上一帧束j-1左相机图像的匹配,②是当前帧束j与关键帧束j-2左相机图像的匹配,③是当前帧束j与关键帧束j2左相机图像的匹配,④是当前帧束j与上一帧束j-1右相机图像的匹配,⑤是当前帧束j与关键帧束j-2右相机图像的匹配,⑥是当前帧束j与关键帧束j2右相机图像的匹配,⑦是当前帧束j中左右两个相机的两帧图像之间的匹配。特征匹配首先执行当前帧束j与前两个关键帧束j-2,j2的匹配—②③⑤⑥,然后判断是否选定当前帧束j为关键帧束。如果当前帧束中特征匹配的个数小于其特征点个数的20%,则选定j为关键帧束。接着,执行当前帧束j与上一帧束j-1的匹配—①④。最后,执行当前帧束j中2帧图像之间的匹配—⑦。每两帧图像的匹配含两个步骤,首先根据特征描述子间的Hamming距离通过暴力搜索寻找最近的特征匹配,然后通过内含5点法的RANSAC方法去除异常匹配。
步骤3,如果滑窗滤波器尚未初始化,利用几个帧束(如2个关键帧)中的特征轨迹和IMU数据初始化状态向量和协方差。在时刻t所有状态量组成的状态向量X(t)一般情况下如图4所示,其中包含当前时刻t的导航状态量π(t)、IMU参数ximu,N个相机的参数 过去的Nkf个关键帧束对应的导航状态量/>和最近Ntf+1个帧束对应的导航状态量/>以及m个地标点L1,L2,…,Lm(在图4中,Nkf=2,Ntf=3)。其中,ximu可以包括IMU的零偏、尺度因子、错位和加速度敏感性(图4实例2),xC可以包含每个相机的投影几何参数、外部参数、相机的时间延迟和卷帘快门效应(图4实例2)。这些过去的Nkf个关键帧束和最近的Ntf个帧束形成一个导航状态量的滑动窗口。
实施例具体实施方案为:
等到步骤2给出两个关键帧束j1,j2和相应的IMU数据后,对滤波器进行初始化。滤波器初始的位置和速度设为零,初始的姿态通过起初的加速度计数据来确定,以使得滤波器的参考世界坐标系{W}的z轴沿着负重力方向。对于IMU的零偏,如果帧束j2相对于j1的像素移动(即光流)较小,则通过平均IMU数据来设置陀螺仪和加速度计的零偏。否则,陀螺仪和加速度计的零偏将被设置为零。相机相对IMU的时间延迟初始化为0。其他的传感器参数和状态量的协方差根据数据表或经验设定。如果传感器事先已经做了很好的标定,可以设定这些传感器参数对应的协方差矩阵的条目为零,来固定这些参数,即这些参数不被更新。
步骤4,对于当前帧束j,通过IMU数据将导航状态量π(t)和协方差矩阵传播到该帧束对应的时刻tj。然后,克隆导航状态量π(tj)=π(t)并添加至状态向量,并且在协方差矩阵中添加对应π(tj)的行和列,这些行和列的值从对应π(t)的行和列复制。
实施例具体实施方案为:
对于帧束j,将导航状态量π(t)和协方差矩阵通过IMU数据传播到该帧束对应的时刻tj。然后,将导航状态量π(t)克隆为π(tj)并添加至状态向量(如图5所示,注意πj∶=π(tj)),并且对协方差矩阵(如图6所示)的相应行和列也进行扩增(如图7所示)。
步骤5,利用当前帧束j中的特征轨迹更新状态向量和协方差。针对特征轨迹可能发生的三种情况,采取不同的方式执行EKF更新。这三种情况包括该特征轨迹在当前帧束中消失、该特征轨迹对应的地标点在状态向量中,以及该特征轨迹对应不在状态向量中但观测良好的地标点。观测良好的地标点需要满足两个条件,一是有足够长的观测,比如在7帧图像中被观测到,二是其对应的地表点可以成功的三角化,即这些观测有足够的视差。对于第一种情况,首先利用该特征轨迹三角化一个地标点,然后计算该特征轨迹上所有观测的重投影误差和雅可比矩阵,之后利用矩阵零空间投影以消除地标点参数的雅可比矩阵,再通过马氏检验去除异常值。对于第二种情况,首先计算特征轨迹在当前帧束j中的观测的重投影误差和雅可比矩阵,然后通过马氏检验去除异常值。对于第三种情况,首先利用该特征轨迹三角化一个地标点,然后计算该特征轨迹上所有观测的重投影误差和雅可比矩阵,通过马氏检验去除异常值,接着添加该地标点到状态向量和协方差矩阵。最后,按以上三种观测情况把所有的特征轨迹分为三类,对每类特征轨迹,叠加上述的重投影误差和雅可比矩阵,对状态向量和协方差矩阵进行更新,更新的方式与经典EKF相同。
实施例具体实施方案为:
如图8所示,对于在当前帧束j中消失的特征轨迹,若其包含至少3个观测,则可作为一个完整的特征轨迹来三角化一个地标点。如果三角化成功,将计算该特征轨迹上所有观测的重投影误差和雅可比矩阵,之后利用矩阵零空间投影以消除地标点参数的雅可比矩阵,进而通过马氏检验去除异常值。对于在状态向量中的地标点的特征轨迹,首先计算当前帧束j中该地标点观测的重投影误差和雅可比矩阵,然后通过马氏检验去除异常值。对于不在状态向量中的地标点的足够长的特征轨迹(如观测个数>7),首先三角化地标点,如果三角化成功,将计算该特征轨迹上所有观测的重投影误差和雅可比矩阵,再通过马氏检验去除异常值,进而将该地标点添加到状态向量中并相应的扩展协方差矩阵。最后,把所有的特征轨迹按以上情况,分为三类,将每一类中所有观测轨迹对应的重投影误差和雅可比矩阵叠在一起,执行一次EKF更新步骤(共三次更新)。
步骤6,为了限制计算量,当滑动窗口中导航状态量的数量超过允许的关键帧束数目Nkf和最近帧束数目Ntf之和时,则从滑动窗口中选取冗余帧束并将其边缘化。因为对一个未知地标点的两个观测不能对位姿提供有效约束,所以在一次边缘化操作中至少选择Nr个冗余帧束(对于单目相机,Nr为3,对于双目相机,Nr为2)。为了满足这个要求,首先在最近的非关键帧束中选择冗余帧束,同时排除最近的Ntf个帧束,然后在最早的关键帧束中选择冗余帧束。在找到Nr个冗余帧束后,可以利用其中长度超过两次观测的特征轨迹来更新滑窗滤波器,具体而言如果这样一条特征轨迹可以利用其所有的观测成功三角化一个地标点,则使用该特征轨迹上冗余帧束内的特征点观测进行EKF更新,类似步骤5第一种情况;冗余帧束中属于其他特征轨迹的特征点将被丢弃。EKF更新后,删除这些冗余帧束对应的导航状态量和协方差矩阵的行和列。
实施例具体实施方案为:
为了选取Nr个冗余帧束(对于单目相机,Nr为3,对于双目相机,Nr为2),如图9所示,先排除最近的Ntf个帧束,然后在最近的非关键帧束中选择冗余帧束(假设选取nr个),然后在最早的关键帧束中选择Nr-nr个冗余帧束。注意去除这些冗余帧束后,至少要保留一个关键帧束。对于图3中Nkf=2和Ntf=3的情况,冗余帧束为关键帧束j1和普通帧束j-3(如图10左图虚线框所示)。对于另外一种假设的情况,如图10右图所示,冗余帧束为最旧的两个帧束。选定冗余帧束后,将利用在冗余帧束中观测次数超过两次的地标点的观测来更新滤波器。如果该地标点可以利用其所有的观测成功三角化,则使用其特征轨迹上冗余帧束内的观测进行EKF更新,类似步骤5第一种情况;冗余帧束中其他的特征点观测将被丢弃。EKF更新后,删除这些冗余帧束对应的导航状态量和协方差矩阵的行和列。
本发明提供的方法,也可以利用多线程设计实现为相应程序,如图2所示。基于关键帧滑窗滤波的视觉惯性里程计的多线程程序包含特征提取线程、基于关键帧的特征关联线程、相机-IMU同步线程和滑窗滤波线程。该多线程实现可以显著提高程序的吞吐率。
特征提取线程,用于检测相机拍摄的每一帧图像的特征点,并用一种特征描述方法对特征点生成描述子。对于由N个相机捕获的含N个时间戳接近图像的一个帧束j,使用如BRISK[7]方法检测特征点并提取描述子。BRISK利用FAST[9]算法进行特征点检测,通过构造尺度空间中的图像金字塔来检测特征。为了生成特征点的描述子,通过比较特征点邻域内512个像素对的灰度值,得到描述子的二进制编码。
基于关键帧的特征关联线程,用于关联当前帧束与先前帧束的图像特征,根据特征提取线程提取的图像特征进行基于关键帧的特征点匹配。具体的匹配步骤包括匹配当前帧束与K个关键帧束,匹配当前帧束与上一帧束,以及两两匹配当前帧束的N张图像。在与K个关键帧束匹配完成后,将根据特征匹配的比例判断是否选定当前帧束为关键帧束。两帧图像特征匹配的一种实现方法是通过暴力搜索的方式找到距每个特征描述子最近的描述子,最后这些匹配会经过一个包含5点法的RANSAC方法来去除外点。
相机-IMU同步线程,用于同步相机和IMU数据。该线程可以通过条件变量来锁住其他线程,以等待对应当前帧束的IMU数据到来。
滑窗滤波线程,利用特征关联线程得到的特征轨迹和IMU数据更新滤波器。如果滤波器尚未初始化,将利用帧束间的特征轨迹、IMU数据、数据表和经验值来初始化滤波器的状态向量和协方差。对于已经初始化的情形,当帧束j到达时,将导航状态量π(t)和协方差矩阵通过IMU数据传播到j帧束的时刻tj,然后克隆导航状态量π(tj)=π(t)并添加到状态向量,并且对协方差矩阵也进行相应扩展。接着,根据特征轨迹的属性,分三种情况执行EKF更新:在当前帧束j中消失的特征轨迹,在状态向量中的地标点对应的特征轨迹,以及不在状态向量中但观测足够长的地标点的特征轨迹。三种情况在如何准备更新方面有所不同。对于第一种情况,利用矩阵零空间投影以消除地标点参数的雅可比矩阵。对于第三种情况,三角化的新地标点将被添加到状态向量和协方差矩阵。这些特征轨迹按这三种情况分为三类。每类执行一次经典的EKF更新。当状态向量中导航状态量的数目超过允许的关键帧束数目Nkf和最近帧束数目Ntf之和时,从中选择冗余帧束并将其边缘化。每次至少选择Nr个冗余帧束(对于单目设备,Nr为3,对于双目设备,Nr为2)。这些帧束选取时,将排除最近的Ntf个帧束,接着在非关键帧束中选择冗余帧束,最后在最早的关键帧束中选择冗余帧束。选取冗余帧束后,可以利用在冗余帧束中观测超过两次的地标点的观测来更新滤波器。找出这些地标点后,选取其中可以成功三角化的地标点,利用他们在冗余帧束中的观测进行一次EKF更新,类似于上述更新的第一种情况。冗余帧束中的其他特征点观测将被丢弃。更新完成后,删除这些冗余帧束对应的导航状态量和协方差矩阵的条目。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各式修改、补充或采用类似的方式替代,但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。
[1]A.Flint,O.Naroditsky,C.P.Broaddus,A.Grygorenko,S.Roumeliotis,andO.Bergig,“Visual-based inertial navigation,”US9424647B2,Aug.23,2016Accessed:Apr.26,2021.[Online].Available:https://patents.google.com/patent/US9424647/nl
[2]S.I.Roumeliotis and K.J.Wu,“Inverse sliding-window filters forvision-aided inertial navigation systems,”US9658070B2,May 23,2017Accessed:Apr.26,2021.[Online].Available:https://patents.google.com/patent/US9658070B2/en?oq=US9658070B2-patent-Inverse+sliding-window+filters+for+vision-aided+inertial+navigation
[3]S.I.Roumeliotis andA.I.Mourikis,“Vision-aided inertialnavigation,”US9766074B2,Sep.19,2017Accessed:Apr.26,2021.[Online].Available:https://patents.google.com/patent/US9766074B2/en
[4]M.Li and A.Mourikis,“Real-time pose estimation system usinginertial and feature measurements,”US9798929B2,Oct.24,2017Accessed:Feb.18,2022.[Online].Available:https://patents.google.com/patent/US9798929B2/en
[5]T.Qin,P.Li,and S.Shen,“VINS-Mono:A robust and versatile monocularvisual-inertial state estimator,”IEEE Trans.Robot.,vol.34,no.4,pp.1004–1020,Aug.2018,doi:10.1109/TRO.2018.2853729.
[6]E.Rublee,V.Rabaud,K.Konolige,and G.Bradski,“ORB:An efficientalternative to SIFT or SURF,”in 2011International Conference on ComputerVision,Nov.2011,pp.2564–2571.doi:10.1109/ICCV.2011.6126544.
[7]S.Leutenegger,M.Chli,and R.Y.Siegwart,“BRISK:Binary robustinvariant scalable keypoints,”in Intl.Conf.on Computer Vision(ICCV),Barcelona,Spain,Nov.2011,pp.2548–2555.doi:10.1109/ICCV.2011.6126542.
[8]A.Alahi,R.Ortiz,and P.Vandergheynst,“FREAK:Fast retina keypoint,”in 2012IEEE Conference on Computer Vision and Pattern Recognition,Jun.2012,pp.510–517.doi:10.1109/CVPR.2012.6247715.
[9]E.Rosten and T.Drummond,“Machine learning for high-speed cornerdetection,”in European conference on computer vision,2006,pp.430–443.

Claims (6)

1.一种基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法,其特征在于,包含以下步骤:
首先在设备上固定安装一个IMU和N个相机,相机用于采集图像帧束的序列,IMU用于采集3轴加速度计数据和3轴陀螺仪数据,然后通过滤波器执行以下步骤估计设备的状态和参数,以及不确定度;以下将位置、速度、姿态以及相机和IMU的参数总称为状态向量,这些状态向量的不确定度由状态向量的协方差来描述;为简化起见,在同一个时刻t的位姿和速度简称为导航状态量;
步骤1,对于由N个相机捕获的含N帧时间相近的图像的帧束j,检测每张图像中的特征点,并对这些特征点生成描述子,即特征向量;
步骤2,利用步骤1提取的特征描述子进行基于关键帧的特征匹配;
步骤3,滤波器尚未初始化时,将通过帧束中的特征轨迹和IMU数据初始化状态向量和协方差;
步骤4,在滤波器完成初始化的情况下,对于完成特征匹配的帧束j,根据IMU的数据,从帧束j-1对应的位姿和速度递推估算帧束j对应时刻tj的位姿和速度,并克隆推算的位姿和速度以扩充滤波器的状态向量和协方差;
步骤5,利用特征点的特征轨迹更新滤波器的状态向量和协方差;
步骤6,更新完毕后,检测是否存在冗余帧,并删除冗余帧对应的导航状态量和协方差矩阵相应的行和列;
对下一帧束j+1,将重复执行步骤1-6;
每次循环都发布滤波器估计的状态向量和协方差,以服务下游的应用。
2.如权利要求1所述的一种基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法,其特征在于:步骤1中特征点的检测方法包括ORB和BRISK,生成描述子的方法包括BRISK、FREAK和基于深度学习的特征描述子生成方法。
3.如权利要求1所述的一种基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法,其特征在于:步骤2中特征匹配过程是将当前帧束j的特征与几个先前帧束的特征相关联,包括K个最近的关键帧和上一帧束j-1,在与K个关键帧的匹配完成后,如果当前帧中特征匹配个数与特征点个数的比例小于Tr,则将帧束j选为关键帧;如果每个帧束含多个相机的N张图像,则在当前帧束j的N张图像之间进行特征关联,两帧图像之间特征描述子的匹配通过暴力搜索、k最近邻或窗口搜索的方法来完成。
4.如权利要求1所述的一种基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法,其特征在于:步骤3中状态向量X(t)包含当前导航状态量π(t)、IMU参数ximu,N个相机的参数过去的Nkf个关键帧束对应的导航状态量/>和最近Ntf+1个帧束对应的导航状态量/>以及m个地标点L1,L2,…,Lm;其中,ximu包括IMU的零偏、尺度因子、错位和加速度敏感性,xC包含每个相机的投影几何参数、外部参数、相机的时间延迟和卷帘快门效应;过去的Nkf个关键帧束和最近的Ntf个帧束形成一个导航状态量的滑动窗口。
5.如权利要求1所述的一种基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法,其特征在于:步骤5中针对每条特征轨迹可能发生的三种情况,采取不同的方式执行EKF更新,这三种情况包括特征轨迹在当前帧中消失、特征轨迹对应的地标点在状态向量中,以及特征轨迹对应不在状态向量中但观测良好的地标点;对于第一种情况,首先利用特征轨迹三角化一个地标点,然后计算特征轨迹上所有观测的重投影误差和雅可比矩阵,之后利用矩阵零空间投影以消除地标点参数的雅可比矩阵,再通过马氏检验去除异常值;对于第二种情况,首先计算特征轨迹在当前帧束j中的观测的重投影误差和雅可比矩阵,然后通过马氏检验去除异常值;对于第三种情况,首先利用特征轨迹三角化一个地标点,然后计算该特征轨迹上所有观测的重投影误差和雅可比矩阵,通过马氏检验去除异常值,接着添加该地标点到状态向量和协方差矩阵;最后,按以上三种观测情况把所有的特征轨迹分为三类,对每类特征轨迹,利用上述的重投影误差和雅可比矩阵对状态向量和协方差矩阵进行更新,更新的方式与经典EKF相同。
6.如权利要求5所述的一种基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法,其特征在于:步骤6的具体实现方式如下;
当滤波器状态向量中导航状态量的数量超过允许的关键帧数目Nkf和最近帧数目Ntf之和时,则从中选取冗余帧束并将其边缘化;每次边缘化操作,至少选择Nr个冗余帧束,为了达到这个要求,首先在最近的非关键帧束中选择冗余帧,同时排除最近的Ntf个帧束,然后在最早的关键帧束中选择冗余帧,在找到Nr个冗余帧后,利用其中长度超过两次观测的特征轨迹来更新滑窗滤波器,即如果这种特征轨迹可以利用其所有的观测成功三角化一个地标点,则使用该轨迹上冗余帧中的特征轨迹进行EKF更新,更新方法和步骤5中的第一种情况相同;冗余帧中属于其他特征轨迹的特征点将被丢弃,EKF更新后,删除这些冗余帧对应的导航状态量和协方差矩阵的行和列。
CN202210158652.2A 2022-02-21 2022-02-21 基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法 Active CN114623817B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202210158652.2A CN114623817B (zh) 2022-02-21 2022-02-21 基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法
PCT/CN2022/080346 WO2023155258A1 (zh) 2022-02-21 2022-03-11 基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210158652.2A CN114623817B (zh) 2022-02-21 2022-02-21 基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法

Publications (2)

Publication Number Publication Date
CN114623817A CN114623817A (zh) 2022-06-14
CN114623817B true CN114623817B (zh) 2024-04-26

Family

ID=81899165

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210158652.2A Active CN114623817B (zh) 2022-02-21 2022-02-21 基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法

Country Status (2)

Country Link
CN (1) CN114623817B (zh)
WO (1) WO2023155258A1 (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116824463B (zh) * 2023-08-31 2023-12-19 江西啄木蜂科技有限公司 视频的关键帧提取方法、计算机可读存储介质及电子设备
CN117826141A (zh) * 2023-12-29 2024-04-05 广东工业大学 一种复杂环境下分布式无人机群的协同定位方法
CN117705107B (zh) * 2024-02-06 2024-04-16 电子科技大学 基于两阶段稀疏舒尔补的面向视觉惯性定位方法
CN117760428B (zh) * 2024-02-22 2024-04-30 西北工业大学 一种基于多立体视觉惯性紧耦合的自主定位方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109166149A (zh) * 2018-08-13 2019-01-08 武汉大学 一种融合双目相机与imu的定位与三维线框结构重建方法与系统
CN111739063A (zh) * 2020-06-23 2020-10-02 郑州大学 一种基于多传感器融合的电力巡检机器人定位方法
CN111750864A (zh) * 2020-06-30 2020-10-09 杭州海康机器人技术有限公司 一种基于视觉地图的重定位方法和装置
CN114001733A (zh) * 2021-10-28 2022-02-01 浙江大学 一种基于地图的一致性高效视觉惯性定位算法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9709404B2 (en) * 2015-04-17 2017-07-18 Regents Of The University Of Minnesota Iterative Kalman Smoother for robust 3D localization for vision-aided inertial navigation
CN111024066B (zh) * 2019-12-10 2023-08-01 中国航空无线电电子研究所 一种无人机视觉-惯性融合室内定位方法
CN111707261A (zh) * 2020-04-10 2020-09-25 南京非空航空科技有限公司 一种微型无人机高速感知和定位方法
CN111780754B (zh) * 2020-06-23 2022-05-20 南京航空航天大学 基于稀疏直接法的视觉惯性里程计位姿估计方法
CN113066127B (zh) * 2021-04-02 2024-04-19 视辰信息科技(上海)有限公司 一种在线标定设备参数的视觉惯性里程计方法和系统
CN113865584B (zh) * 2021-08-24 2024-05-03 知微空间智能科技(苏州)有限公司 一种基于视觉惯性里程计的uwb三维寻物方法和装置

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109166149A (zh) * 2018-08-13 2019-01-08 武汉大学 一种融合双目相机与imu的定位与三维线框结构重建方法与系统
CN111739063A (zh) * 2020-06-23 2020-10-02 郑州大学 一种基于多传感器融合的电力巡检机器人定位方法
CN111750864A (zh) * 2020-06-30 2020-10-09 杭州海康机器人技术有限公司 一种基于视觉地图的重定位方法和装置
CN114001733A (zh) * 2021-10-28 2022-02-01 浙江大学 一种基于地图的一致性高效视觉惯性定位算法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
动态场景下视觉里程计算法研究与实现;王子国;cnki硕士论文电子期刊;20220515;全文 *
基于Realsense传感器的机器人视觉里程计研究;廖萱;陈锐志;李明;;地理空间信息;20200228(第02期);全文 *

Also Published As

Publication number Publication date
WO2023155258A1 (zh) 2023-08-24
CN114623817A (zh) 2022-06-14

Similar Documents

Publication Publication Date Title
CN114623817B (zh) 基于关键帧滑窗滤波的含自标定的视觉惯性里程计方法
CN112197770B (zh) 一种机器人的定位方法及其定位装置
CN111024066B (zh) 一种无人机视觉-惯性融合室内定位方法
JP7326720B2 (ja) 移動体位置推定システムおよび移動体位置推定方法
CN112734852B (zh) 一种机器人建图方法、装置及计算设备
US9071829B2 (en) Method and system for fusing data arising from image sensors and from motion or position sensors
CN110084832B (zh) 相机位姿的纠正方法、装置、系统、设备和存储介质
CN111795686B (zh) 一种移动机器人定位与建图的方法
CN110411476B (zh) 视觉惯性里程计标定适配及评价方法和系统
CN110726406A (zh) 一种改进的非线性优化单目惯导slam的方法
CN107167826B (zh) 一种自动驾驶中基于可变网格的图像特征检测的车辆纵向定位系统及方法
CN108519102B (zh) 一种基于二次投影的双目视觉里程计算方法
CN110553648A (zh) 一种用于室内导航的方法和系统
CN112747750B (zh) 一种基于单目视觉里程计和imu融合的定位方法
CN110207693B (zh) 一种鲁棒立体视觉惯性预积分slam方法
CN111986261B (zh) 一种车辆定位方法、装置、电子设备及存储介质
Michot et al. Bi-objective bundle adjustment with application to multi-sensor slam
CN111882602B (zh) 基于orb特征点和gms匹配过滤器的视觉里程计实现方法
Zhang et al. Vision-aided localization for ground robots
Cai et al. Mobile robot localization using gps, imu and visual odometry
CN110751123B (zh) 一种单目视觉惯性里程计系统及方法
CN111623773B (zh) 一种基于鱼眼视觉和惯性测量的目标定位方法及装置
US20160055646A1 (en) Method for estimating the angular deviation of a mobile element relative to a reference direction
Seiskari et al. Hybvio: Pushing the limits of real-time visual-inertial odometry
Chen et al. Stereo visual inertial pose estimation based on feedforward-feedback loops

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