CN115406447A - 拒止环境下基于视觉惯性的四旋翼无人机自主定位方法 - Google Patents
拒止环境下基于视觉惯性的四旋翼无人机自主定位方法 Download PDFInfo
- Publication number
- CN115406447A CN115406447A CN202211348263.2A CN202211348263A CN115406447A CN 115406447 A CN115406447 A CN 115406447A CN 202211348263 A CN202211348263 A CN 202211348263A CN 115406447 A CN115406447 A CN 115406447A
- Authority
- CN
- China
- Prior art keywords
- frame
- image
- imu
- unmanned aerial
- aerial vehicle
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 60
- 230000000007 visual effect Effects 0.000 title claims description 39
- 238000005259 measurement Methods 0.000 claims abstract description 55
- 230000010354 integration Effects 0.000 claims abstract description 37
- 238000000605 extraction Methods 0.000 claims abstract description 21
- 238000012545 processing Methods 0.000 claims abstract description 14
- 238000001514 detection method Methods 0.000 claims abstract description 13
- 230000003287 optical effect Effects 0.000 claims description 47
- 239000013598 vector Substances 0.000 claims description 37
- 238000005457 optimization Methods 0.000 claims description 36
- 239000011159 matrix material Substances 0.000 claims description 13
- 238000010168 coupling process Methods 0.000 claims description 10
- 230000004044 response Effects 0.000 claims description 9
- 238000012795 verification Methods 0.000 claims description 8
- 230000002159 abnormal effect Effects 0.000 claims description 6
- 238000011156 evaluation Methods 0.000 claims description 6
- 238000003672 processing method Methods 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 5
- 230000001360 synchronised effect Effects 0.000 claims description 5
- 230000008878 coupling Effects 0.000 claims description 4
- 238000005859 coupling reaction Methods 0.000 claims description 4
- 238000005516 engineering process Methods 0.000 claims description 4
- 230000000750 progressive effect Effects 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 abstract 1
- 230000006870 function Effects 0.000 description 10
- 238000010586 diagram Methods 0.000 description 6
- 230000001133 acceleration Effects 0.000 description 3
- 239000000654 additive Substances 0.000 description 2
- 230000000996 additive effect Effects 0.000 description 2
- 239000000284 extract Substances 0.000 description 2
- 230000005484 gravity Effects 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000009827 uniform distribution Methods 0.000 description 2
- 230000000903 blocking effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000005764 inhibitory process Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000008447 perception Effects 0.000 description 1
- 230000008092 positive effect Effects 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/20—Instruments for performing navigational calculations
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C11/00—Photogrammetry or videogrammetry, e.g. stereogrammetry; Photographic surveying
- G01C11/04—Interpretation of pictures
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
- G01C21/12—Navigation; 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/16—Navigation; 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/165—Navigation; 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/1656—Navigation; 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
- G06T7/246—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
- G06T7/269—Analysis of motion using gradient-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/70—Determining position or orientation of objects or cameras
- G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Multimedia (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Theoretical Computer Science (AREA)
- Automation & Control Theory (AREA)
- Image Analysis (AREA)
Abstract
本发明涉及无人机定位技术领域,尤其涉及一种拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,首先设计前端数据处理模块,预处理双目相机图像信息和IMU测量信息,对图像进行改进的快速特征提取与跟踪,并将图像帧与IMU帧对齐,给出IMU预积分形式;其次基于滑动窗口模型构建位姿估计最小化误差函数,减少位姿估计的累积误差;最后利用词袋模型进行回环检测,消除漂移,构建全局一致的运动轨迹,实现完整的无人机自主定位。本发明方法不依赖于传统的GPS系统与繁琐的地面站控制,仅靠重量较轻的双目相机与IMU,并通过机载计算机进行实时处理,能够有效提高无人机的高机动自主定位性能,具备高精度、高可靠性的优点。
Description
技术领域
本发明涉及无人机定位技术领域,尤其涉及一种拒止环境下基于视觉惯性的四旋翼无人机自主定位方法。
背景技术
四旋翼无人机由于其低成本、小尺寸和优越的机动性,正逐渐成为室内外环境中广泛应用的理想机器人平台,在四旋翼无人机的应用中,其自主定位能力在导航运动规划中至关重要。如GPS的大量商业化应用为室外开阔区域的设备提供了较为简单的定位方案,凭借全方位、全天候、全时段和高精度的优势,以GPS为代表的卫星定位系统被广泛应用于各类机器人设备。而卫星信号被干扰、遮挡以及欺骗攻击等导致GPS定位设备无法正常输出参数的情形被称为拒止环境,如深山峡谷、密集森林、高层楼宇间、室内环境以及战场区域。在GPS拒止环境中,四旋翼无人机由于无法获得精确的先验定位信息,而无法在上述区域执行自主的高机动复杂飞行任务。
为了实现自主导航和高效安全地执行任务,无人机需要能够根据环境对自身进行定位。在定位领域,同步定位与建图是一种在未知环境下估计传感器运动和重建环境信息的技术,而仅利用相机视觉和IMU(Inertial Measurement Unit,惯性测量单元)信息进行的同步定位与建图被称为视觉惯性里程计(VIO)。相比较激光雷达,相机和IMU具有重量轻、尺寸小、成本低的优点,是无人机低负载情况下良好的传感器设备。VIO融合了相机图像信息和IMU数据信息,用于定位和环境感知,其经典VIO算法框架主要包括前端数据处理、后端优化、回环检测和建图。
发明内容
本发明公开一种拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,实现无人机在无先验信息的未知环境中完成对自身位姿的精确估计。
为了实现本发明的目的,所采用的技术方案是:拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,四旋翼无人机在GPS拒止环境下的复杂场景中进行快速的高机动飞行,定位方法基于VIO算法框架,融合双目相机视觉图像信息和IMU测量信息进行同步定位,并使用基于优化的紧耦合方法,实现四旋翼无人机位姿估计;
定位方法具体包括如下步骤:
S1、在前端数据处理部分,对于安装在四旋翼无人机前部的双目相机所实时采集到的图像帧序列,依靠图像处理方法对视觉图像进行特征提取与跟踪,并剔除其中特征跟踪匹配失败的异常点;
S2、在前端数据处理部分,针对IMU测量信息,将图像观测频率(双目相机的观测频率)与IMU测量频率进行对齐,在两个连续图像帧之间利用IMU预积分技术进行处理,减少IMU积分次数,将图像帧之间的IMU测量数据的积分形式约束在相对运动状态之外;避免因迭代优化造成重复积分,简化后端位姿优化模型;
S3、在后端优化估计部分,基于滑动窗口的后端优化方法来减少VIO算法框架中位姿估计的累积误差,通过非线性优化,用紧耦合的方法构建滑动窗口位姿估计最小化误差函数;其中融合了视觉测量和IMU测量;
S4、四旋翼无人机由于长时间的运动,对于返回已访问的区域难免会有累积漂移,运用基于描述符的词袋模型进行回环检测,将视觉图像的特征描述符转换成词袋向量,计算相似度,通过验证,完成回环检测。
作为本发明的优化方案,基于优化的方法通过使用双目相机视觉图像信息和IMU测量信息联合最小化残差来获得最优估计;而紧耦合则将双目相机的全局状态和IMU的全局状态融合到一个观测方程中进行状态估计。
作为本发明的优化方案,在步骤S1中,对于双目图像信息,基于图像特征提取与跟踪以获得图像与三维空间的2D-3D对应关系,采用改进的基于FAST角点提取与改进的KLT光流跟踪的图像处理方法,并基于PROSAC方法剔除异常点。
作为本发明的优化方案,改进的FAST角点提取方法,针对候选角点,定义为中心像素p,像素灰度值I p ,按顺时针方向计算像素圆中像素灰度值I i ,若由16个像素组成的像素圆中存在12个连续的像素i满足公式1,则中心像素p为角点,否则筛除像素p;
其中,T表示为灰度阈值,circle(p)表示为围绕中心像素p的像素圆;
将角点与其对应的连续12个像素组成的像素弧之间灰度差绝对值的总和作为特征响应值函数S,如公式2所示,保留特征响应最大的角点作为唯一角点;
作为本发明的优化方案,改进后的KLT光流跟踪以连续图像帧中的前一帧I (x, y)和后一帧J (x, y)的原始图像作为L0层,建立图像金字塔模型,依次生成L1-L4层,前一层图像的大小是后一层图像的4倍,图像I (x, y)和J (x, y)的金字塔模型记为:
在金字塔模型中自顶而下迭代计算第L层的光流场矢量,建立以残差光流矢量 L d为变量的最小化匹配跟踪误差函数,
其中,x表示为图像帧中角点的x像素坐标值,y表示为图像帧中角点的y像素坐标值, L e表示为第L层的最小化匹配跟踪误差函数, L d表示为第L层迭代的角点残差光流矢量, L d x 表示为第L层中x方向上迭代的角点残差光流矢量, L d y 表示为第L层中y方向上迭代的角点残差光流矢量, L u x 表示为第L层迭代的角点像素二维坐标x点, L u y 表示为第L层迭代的角点像素二维坐标y点,w x = w y = 2,{ L u x - w x , L u x + w x , L u y – w y , L u y + w y }组成邻域为(2w x + 1) × (2w y + 1)的像素窗口, L g表示为第L层迭代的角点光流估计值, L g x 表示为第L层中x方向上迭代的角点光流估计值, L g y 表示为第L层中y方向上迭代的角点光流估计值;
基于PROSAC方法为基于PROSAC渐进一致采样算法,降序排列在前后图像帧中特征跟踪点对(因为后一帧图像上的特征点要跟踪前一帧图像上的特征点,所以这两个点形成了一对特征跟踪点对)的最小欧氏距离,将其作为评价判断,选取一定量评价值最高的数据,从中随机采样4组数据通过迭代,寻找到最优的参数模型使得能匹配上的特征点对最多,记为模型M,剔除与模型M的误差大于阈值的点对。
作为本发明的优化方案,在步骤S2中,以连续图像帧i的时间刻度为基准,针对IMU测量数据冗余的问题,将双目相机图像帧与IMU测量帧进行帧对齐(将图像观测频率与IMU测量频率进行对齐),并在对齐帧之间进行IMU预积分,计算IMU数据在图像帧之间的积分,运用预积分方法,将图像帧之间的IMU测量数据的积分形式约束在相对运动状态之外,避免因迭代优化造成重复积分
有如下IMU预积分形式:
其中,表示为机体坐标系b中第i帧和第i + 1帧之间位置的预积分量,表
示为机体坐标系b中第i帧和第i + 1帧之间速度的预积分量,表示为机体坐标系b中第i帧和第i + 1帧之间旋转的预积分量,t表示为第i帧和第i + 1帧之间的时间间隔,表
示为机体坐标系b中t时刻下相对于第i帧的旋转矩阵,表示为机体坐标系b中第i帧和t
时刻下帧之间旋转的预积分量,表示为t时刻下IMU的加速度计原始测量值,表示为t
时刻下IMU的陀螺仪原始测量值,表示为t时刻下IMU的加速度计偏置,表示为t时刻
下IMU的陀螺仪偏置。
作为本发明的优化方案,在步骤S3中,通过VINS-Mono算法中基于滑动窗口的后端优化方法来减少VIO算法框架中位姿估计的累积误差,通过非线性优化,用紧耦合的方法构建滑动窗口误差函数E,最小化边缘化先验信息、IMU测量残差和视觉测量残差,三种残差均用马氏距离(与量纲无关)来表示,即:
其中,
其中,表示为视觉和IMU的全状态向量,r p -H p X 表示为边缘化先验信息,表示为第i帧和第i + 1帧之间的IMU测量残差,表示为第i帧图像中
首次观测到的第l个特征点在第j帧图像中特征点观测的视觉测量残差,B是所有IMU测量值
的集合,C是在当前滑动窗口中观察到至少两次的特征集,表示为IMU预积分噪声项的
协方差矩阵, 表示为视觉观测噪声的协方差矩阵,i表示为第i帧图像,l表示为第l个
特征点,j表示为第j帧图像,ρ(s)表示为Huber范数,s为中间变量。
作为本发明的优化方案,在步骤S4中,基于普通BRIEF描述符的DBoW词袋模型进行回环检测,根据BRIEF描述符创建二进制特征向量,用以描述特征点周围的像素信息,并将其转换成词袋向量,计算双目相机当前帧与词袋模型的相似度,通过外观特征验证和几何一致性验证判断两个不同图像位置的角点是否相似,完成回环检测,重新估计四旋翼无人机位姿。
本发明具有积极的效果:1)本发明根据VIO算法框架设计了一种拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,在GPS信号被干扰、遮挡或是欺骗攻击时,能够在无先验信息的未知环境中完成对自身位姿的精确估计;
2)本发明所提出的方法在保证四旋翼无人机高自主性的飞行性能下,不需要依靠传统的GPS系统与冗余的地面站控制,最大程度减少了对传统的、低精度的、不可靠设备的依赖;
3)本发明所使用的外部传感器仅靠重量较轻、精度较高的双目相机与IMU,运用较为先进的算法进行处理,并通过搭载在无人机上的高性能机载计算机进行实时运算,能够有效提高无人机的高机动自主定位性能,具备高精度、高可靠性的优点。
附图说明
下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1为本发明实施例的拒止环境下基于视觉惯性的四旋翼无人机自主定位方法总体系统框图;
图2为本发明实施例中改进的FAST角点提取示意图;
图3为本发明实施例中改进的FAST角点提取步骤结果对比图;
图4为本发明实施例中四叉树金字塔结构KLT稀疏光流法示意图;
图5为本发明实施例中特征跟踪算法对比图;
图6为本发明实施例中图像帧与IMU帧对齐示意图。
具体实施方式
下面结合附图和具体实施例对本发明的方法进行进一步说明。
如图1所示,本发明设计了拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,四旋翼无人机在GPS拒止环境下的复杂场景中进行快速的高机动飞行,定位方法基于VIO算法框架,融合双目相机视觉图像信息和IMU测量信息进行同步定位,并使用基于优化的紧耦合方法,实现四旋翼无人机位姿估计;基于优化的方法通过使用双目相机视觉图像信息和IMU测量信息联合最小化残差来获得最优估计;而紧耦合则将双目相机的全局状态和IMU的全局状态融合到一个观测方程中进行状态估计。
定位方法具体包括如下步骤:
S1、在前端数据处理部分,设计一种轻量级的视觉处理系统,对于安装在四旋翼无人机前部的双目相机所实时采集到的图像帧序列,依靠图像处理方法对视觉图像进行特征提取与跟踪,并剔除其中特征跟踪匹配失败的异常点,图像特征提取与跟踪用以获得图像与三维空间的2D-3D对应关系;
在步骤S1中,采用改进的基于FAST角点提取与改进的KLT光流跟踪的图像处理方法,并基于PROSAC方法剔除异常点,为后续融合IMU测量信息估算四旋翼无人机位姿建立基础。
FAST是一种用于图像检测中较典型的角点提取方法,主要用于检测局部像素灰度水平的明显变化。
本发明通过改进的FAST角点提取法,自适应高低阈值对均匀划分的图像帧提取角点,对不同亮度和纹理的图像区域均能提取到一定量的角点,再通过四叉树分裂原理对特征分布不均匀的图像进行分块切割,采取类似非极大值抑制的方法制定特征响应值函数,保留高响应特征,剔除低响应特征,实现图像中特征的均匀分布,进一步改善特征点的鲁棒性。FAST角点提取示意图如图2所示,其中,像素p表示为待提取的中心像素,编号1-16的像素表示为围绕中心像素p的像素圆,由方框标记编号1、 5、9、 13的像素表示为算法中首先与中心像素p作比较的像素。
以公共数据集The EuRoC MAV Dataset中的双目视觉图像为例,具体改进的FAST角点提取方法步骤如下:
(1)首先确定最终图像帧中均匀分布的特征点数量需保持在90-100之间;
(2)将图像帧按照30 × 30的固定尺寸划分网格,对每个网格提取FAST角点,设待提取像素为中心像素p,其灰度值设为I p ,设高阈值T h = 0.2I p 和低阈值T l = 0.1I p ,围绕中心像素p的16个像素组成的像素圆的灰度值设为I i (i∈1, 2, …, 16);
(3)在高阈值T h 条件下,对由16个像素组成的像素圆进行采样,首先计算像素1、5、9、13的灰度值I i 与中心像素p的灰度值I p 之间灰度差的绝对值,若至少存在3个像素满足公式(1),则中心像素p为候选角点,继续下一步考量,否则筛除像素p;
(4)针对候选角点,按顺时针方向继续计算像素圆中其余的像素灰度值I i ,若由16个像素组成的像素圆中存在12个连续的像素i满足公式(1),则中心像素p为角点,否则筛除像素p;
其中,T表示为灰度阈值,circle(p)表示为围绕中心像素p的像素圆。
(5)若使用高阈值T h 未在该网格区域内提取到角点,则自动降低阈值条件,使用低阈值T l 参照步骤(3)-(4)重新提取,使得在弱纹理区域也能提取到一定量的角点,而若使用低阈值也未提取到该网格区域内的角点时,则放弃提取,避免提取到质量差的角点,非均匀化特征点如图3(a)所示;
(6)在由上述步骤所得到的非均匀分布特征点图像帧中构建四叉树网络节点。以初始图像帧为四叉树的单位根节点,将1个根节点分裂为4个子节点,如图3(b)所示,分别统计每个子节点内所包含的FAST角点数量。若节点中角点数量为0,则删除该节点;若节点中角点数量为1,则停止分裂该节点;判断此时子节点总数是否满足特征点需求数量90-100,若不满足,则继续分裂角点数量大于1的节点,且优先分裂角点数量较多的节点,直至子节点总数满足期望数量90-100,如图3(c)-(d)所示;
(7)在上一步骤所得到的全部子节点中,若存在角点数量大于1的节点,采取类似非极大值抑制的方法,将角点与其对应的连续12个像素组成的像素弧之间灰度差绝对值的总和作为特征响应值函数S (Score Function),如公式(2)所示,保留特征响应最大的角点作为该节点的唯一角点。由此便得到了数量适中且分布均匀的特征点,同时也保证了所得特征点的鲁棒性,如图3(e)所示。
其中,arc(p)表示为围绕中心像素p的像素弧。
图3(f)所示即为最终FAST角点的提取效果,与图3(a)进行对比,可以明显地看出改进后的FAST角点具有显著的均匀分布状态,避免了特征点的严重扎堆,能够为后续的特征匹配跟踪提供良好的特征性能。
光流表示一个空间物体在成像平面上像素运动的瞬时速度矢量,其运动信息是由相邻帧之间密集像素的对应关系来计算的,而传统的KLT稀疏光流法是一种通过稀疏像素点处理微小运动的方法。本发明与上文四叉树均匀化特征点相结合,运用包含四叉树和金字塔结构的KLT稀疏光流方法,建立分层金字塔模型,利用其分层缩小化像素运动,计算光流场矢量,层次化迭代执行四叉树均匀特征点跟踪,如图4所示。传统的光流跟踪不包括异常值过滤过程,本发明在所提出的特征跟踪方案中添加PROSAC方法优化匹配跟踪结果,该方法将所有匹配跟踪结果按照质量进行采样排序,更加快捷有效地剔除误匹配。在保持精度的同时,提高视觉数据处理系统的运行速度,节省计算资源。
总体而言,改进的四叉树金字塔结构KLT稀疏光流特征跟踪(简称改进后的KLT光流跟踪)遵从光流法的三个假设条件:(A)亮度恒定性,图像中物体的像素灰度在连续帧之间不会发生变化;(B)空间一致性,在一定邻域内相邻像素具有相似的运动;(C)微小运动,相邻帧之间的时间足够短,且物体运动较小,由金字塔模型的缩放保证该假设成立。具体改进后的算法步骤流程如下所示:
(1)以连续图像帧中的前一帧I (x, y)和后一帧J (x, y)的原始图像作为L0层,建立图像金字塔模型,依次生成L1-L4层,前一层图像的大小是后一层图像的4倍,图像I(x, y)和J (x, y)的金字塔模型记为:
(2)使用改进后的FAST角点提取方法,利用四叉树结构提取每层图像的FAST角点。设金字塔图像中第L层迭代的角点光流估计值 L g = [ L g x L g y ]T,设第L层迭代的角点残差光流矢量 L d = [ L d x L d y ]T,金字塔模型中迭代的光流场矢量由光流估计值与残差光流矢量构成。在前一帧图像I (x, y)中设某个待匹配跟踪的角点像素二维坐标u = [u x u y ]T,其第L层迭代的角点像素二维坐标 L u = [ L u x L u y ]T;
(3)根据假设条件(B),在一定邻域内相邻像素具有相似的运动,设该邻域为(2w x + 1) × (2w y + 1)的像素窗口,w x 和w y 的取值一般为2, 3, 4, 5, 6, 7像素,此处取w x =w y = 2,在金字塔模型中自顶而下迭代计算第L层的光流场矢量,建立以残差光流矢量 L d为变量的最小化匹配跟踪误差函数:
(4)由经典Lucas-Kanade光流公式对邻域内像素定义可知,对于金字塔模型内的像素窗口进行简化,存在定义:
由Lucas-Kanade光流公式可得第L层迭代的角点残差光流矢量为:
其中:
(5)初始化最顶层L m 层光流估计值,令:
由于金字塔每层级是上一层级的4倍,即长宽尺寸是上一层级的两倍,故第L - 1层迭代的角点光流估计值为:
(6)故最终迭代得到的光流场矢量为:
其中, 0 g表示为第0层光流估计值, 0 d表示为第0层迭代的角点残差光流矢量;
则在后一帧图像J (x, y)中匹配跟踪的角点像素二维坐标为:
v = u + d (15)
其中,v表示为在后一帧图像J (x, y)中匹配跟踪的角点像素二维坐标,u表示为前一帧图像I (x, y)中待匹配跟踪的角点像素二维坐标,d表示为最终迭代得到的光流场矢量。
(7)基于PROSAC渐进一致采样算法,降序排列在前后图像帧中特征跟踪点对的最小欧氏距离,将其作为评价判断,选取一定量评价值最高的数据,从中随机采样4组数据通过迭代,寻找到最优的参数模型使得能匹配上的特征点对最多,记为模型M,剔除与模型M的误差大于阈值的点对。
以EuRoC数据集为例,使用本节特征匹配跟踪方法前后的具体对比结果分别如图5(a)、图5(b)所示,可以看出本发明的算法能够更加有效地跟踪稀疏特征点,并合理剔除误匹配点。
S2、在前端数据处理部分,相机的观测频率一般为几十HZ,而IMU的测量频率能达到几百HZ,在双目相机拍摄的两帧图像之间的时间间隔内会有大量的IMU测量数据,且对每个IMU数据进行积分用以更新位姿状态会占用大量的运算资源。针对IMU测量数据冗余的问题,将图像观测频率与IMU测量频率进行对齐,在两个连续图像帧之间利用IMU预积分技术进行处理,减少IMU积分次数,将图像帧之间的IMU测量数据的积分形式约束在相对运动状态之外,避免因迭代优化造成重复积分,简化后端位姿优化模型。
以连续图像帧i的时间刻度为基准,将相机图像帧与IMU测量帧进行帧对齐,并在对齐帧之间进行IMU预积分,如图6所示。
定义机体坐标系b与世界坐标系w。在机体坐标系b中,IMU的加速度计和陀螺仪原始测量值分别如下式所示:
其中,表示为t时刻下IMU的加速度计原始测量值,表示为t时刻下IMU的陀螺
仪原始测量值,a t 表示为t时刻下IMU的加速度真实值,ω t 表示为t时刻下IMU的角速度真实
值,表示为t时刻下IMU的加速度计偏置,表示为t时刻下IMU的陀螺仪偏置,表示
为在t时刻下世界坐标系w相对于IMU的旋转矩阵,g w 表示为在世界坐标系w中的重力加速
度,n a 表示为IMU加速度计的附加噪声,n ω 表示为IMU陀螺仪的附加噪声。
在机体坐标系b中给定对齐帧中的两帧,即第i帧和第i + 1帧,记为b i 和b i + 1。在世界坐标系w中对时间间隔[t i , t i+ 1]内的IMU测量数据积分可得第i + 1帧的位置、速度和旋转,此处旋转采用四元数的表示形式,分别如下所示:
其中:
其中,表示为机体坐标系b中第i帧相对于世界坐标系w的位置,表示为机
体坐标系b中第i+1帧相对于世界坐标系w的位置,表示为机体坐标系b中第i帧相对于世
界坐标系w的速度,表示为机体坐标系b中第i+1帧相对于世界坐标系w的速度,Δt i 表示
为时间间隔[t i , t i+ 1]之内的连续时间;表示为t时刻下相对于世界坐标系w的旋转矩
阵,表示为机体坐标系b中第i帧相对于世界坐标系w的旋转,表示为机体坐标系b中
第i+1帧相对于世界坐标系w的旋转,表示为t时刻下相对于机体坐标系b中第i帧的旋
转,t表示为第i帧和第i + 1帧之间的时刻,表示为角速度,表示为x轴上的角速度,表示为y轴上的角速度,表示为z轴上的角速度。
由式(18)可知,第i + 1帧的IMU积分依赖于初始帧第i帧的初始位置、速度和旋转,在后端优化时,需要不断迭代求解初始帧第i帧机体坐标系到世界坐标系的初始旋转变量,即t i 时刻在世界坐标系下的初始姿态,以此来更新第i帧的位置、速度和旋转,故每次迭代后都需根据优化变量重新进行积分,这将非常损耗运算资源。
其中:
其中,表示为在第i帧时世界坐标系w到机体坐标系b的旋转矩阵,表示为
机体坐标系b中第i+1帧相对于世界坐标系w的位置,表示为机体坐标系b中第i帧相对于
世界坐标系w的位置,表示为机体坐标系b中第i帧相对于世界坐标系w的速度, 表示
为机体坐标系b中第i帧相对于世界坐标系w的旋转,表示为机体坐标系b中第i+1帧相对
于世界坐标系w的旋转,表示为世界坐标系w相对于机体坐标系b中第i帧的旋转,t表示
为第i帧和第i + 1帧之间的时刻,g w 表示为在世界坐标系w中的重力加速度,表示为机
体坐标系b中第i帧和第i + 1帧之间位置的预积分量,表示为机体坐标系b中第i帧和
第i + 1帧之间速度的预积分量,表示为机体坐标系b中第i帧和第i + 1帧之间旋转的
预积分量,表示为机体坐标系b中t时刻下相对于第i帧的旋转矩阵,表示为机体坐
标系b中第i帧和t时刻下帧之间旋转的预积分量,表示为t时刻下IMU的加速度计原始测
量值,表示为t时刻下IMU的陀螺仪原始测量值,表示为t时刻下IMU的加速度计偏置,表示为t时刻下IMU的陀螺仪偏置。
由式(5)可知,积分项中的旋转矩阵变为了,而在初始帧t i 时刻的初值,这样便消除了对初值旋转优化估计的依赖。由此便得到预积分项式(5),此积分可以通过在机体坐标系b中IMU的原始测量值直接单独获得,且只与第i帧和第i + 1帧之间的IMU偏置有关,与其他状态量无关。
S3、在后端优化估计部分,基于滑动窗口的后端优化方法来减少VIO算法框架中位姿估计的累积误差,通过非线性优化,用紧耦合的方法构建滑动窗口位姿估计最小化误差函数,其中融合了视觉测量和IMU测量。
在滑动窗口内定义n个对齐帧内IMU的全状态向量与m个特征点的全状态向量,即:
其中,表示为全状态向量,x i 表示为第i个对齐帧下IMU的状态向量,i=0,1,……
n;表示为双目相机的状态向量,λ m 表示为三维特征点的逆深度信息,m=0,1,……m;ba表
示为在机体坐标系b中IMU的加速度计偏置,bg表示为在机体坐标系b中IMU的陀螺仪偏置,表示为相机坐标系c到机体坐标系b的位移,表示为相机坐标系c到机体坐标系b的旋
转。
在步骤S3中,通过VINS-Mono算法中基于滑动窗口的后端优化方法来减少VIO算法框架中位姿估计的累积误差,通过非线性优化,用紧耦合的方法构建滑动窗口误差函数E,最小化边缘化先验信息、IMU测量残差和视觉测量残差,三种残差均用马氏距离(与量纲无关)来表示,即:
其中,
其中,表示为视觉和IMU的全状态向量,r p -H p X 表示为边缘化先验信息,表示为第i帧和第i + 1帧之间的IMU测量残差,表示为第i帧图像中
首次观测到的第l个特征点在第j帧图像中特征点观测的视觉测量残差,B是所有IMU测量值
的集合,C是在当前滑动窗口中观察到至少两次的特征集,表示为IMU预积分噪声项的
协方差矩阵, 表示为视觉观测噪声的协方差矩阵,i表示为第i帧图像,l表示为第l个
特征点,j表示为第j帧图像。
紧耦合则将双目相机的全局状态和IMU的全局状态融合到一个观测方程中进行状态估计,该观测方程即为公式(6)。
S4、无人机由于长时间的运动,对于返回已访问的区域难免会有累积漂移,本发明基于普通BRIEF描述符的DBoW词袋模型进行回环检测。根据BRIEF描述符创建二进制特征向量,用以描述特征点周围的像素信息,并将其转换成词袋向量,计算双目相机当前帧与词袋模型的相似度,通过外观特征验证和几何一致性验证判断两个不同图像位置的角点是否相似,完成回环检测,重新估计位姿信息。
本发明基于机载计算机进行开发和设计,所用的机载计算机为Intel NUC7系列套件,搭载Ubuntu 18.04系统,运行内存8GB,CPU处理器为Intel Core i5-7260U @ 2.20GHz× 4,运行ROS Melodic版本。
以上所述的具体实施例,对本发明的目的、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限制本发明,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (8)
1.拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,其特征在于:四旋翼无人机在GPS拒止环境下的复杂场景中进行快速的高机动飞行,定位方法基于VIO算法框架,融合双目相机视觉图像信息和IMU测量信息进行同步定位,并使用基于优化的紧耦合方法,实现四旋翼无人机位姿估计;
定位方法具体包括如下步骤:
S1、在前端数据处理部分,对于安装在四旋翼无人机前部的双目相机所实时采集到的图像帧序列,依靠图像处理方法对视觉图像进行特征提取与跟踪,并剔除其中特征跟踪匹配失败的异常点;
S2、在前端数据处理部分,针对IMU测量信息,将图像观测频率与IMU测量频率进行对齐,在两个连续图像帧之间利用IMU预积分技术进行处理,将图像帧之间的IMU测量数据的积分形式约束在相对运动状态之外;
S3、在后端优化估计部分,基于滑动窗口的后端优化方法来减少VIO算法框架中位姿估计的累积误差,通过非线性优化,用紧耦合的方法构建滑动窗口位姿估计最小化误差函数;
S4、四旋翼无人机由于长时间的运动,对于返回已访问的区域难免会有累积漂移,运用基于描述符的词袋模型进行回环检测,将视觉图像的特征描述符转换成词袋向量,计算相似度,通过验证,完成回环检测。
2.根据权利要求1所述的拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,其特征在于:基于优化的方法通过使用双目相机视觉图像信息和IMU测量信息联合最小化残差来获得最优估计;而紧耦合则将双目相机的全局状态和IMU的全局状态融合到一个观测方程中进行状态估计。
3.根据权利要求1所述的拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,其特征在于:在步骤S1中,采用改进的基于FAST角点提取与改进的KLT光流跟踪的图像处理方法,并基于PROSAC方法剔除异常点。
5.根据权利要求3所述的拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,其特征在于:改进的KLT光流跟踪以连续图像帧中的前一帧I (x, y)和后一帧J (x, y)的原始图像作为L0层,建立图像金字塔模型,依次生成L1-L4层,前一层图像的大小是后一层图像的4倍,图像I (x, y)和J (x, y)的金字塔模型记为:
在金字塔模型中自顶而下迭代计算第L层的光流场矢量,建立以残差光流矢量 L d为变量的最小化匹配跟踪误差函数,
其中,x表示为图像帧中角点的x像素坐标值,y表示为图像帧中角点的y像素坐标值, L e表示为第L层的最小化匹配跟踪误差函数, L d表示为第L层迭代的角点残差光流矢量, L d x 表示为第L层中x方向上迭代的角点残差光流矢量, L d y 表示为第L层中y方向上迭代的角点残差光流矢量, L u x 表示为第L层迭代的角点像素二维坐标x点, L u y 表示为第L层迭代的角点像素二维坐标y点,w x = w y = 2,{ L u x - w x , L u x + w x , L u y – w y , L u y + w y }组成邻域为(2w x +1) × (2w y + 1)的像素窗口, L g表示为第L层迭代的角点光流估计值, L g x 表示为第L层中x方向上迭代的角点光流估计值, L g y 表示为第L层中y方向上迭代的角点光流估计值;
基于PROSAC方法为基于PROSAC渐进一致采样算法,降序排列在前后图像帧中特征跟踪点对的最小欧氏距离,将其作为评价判断,选取一定量评价值最高的数据,从中随机采样4组数据通过迭代,寻找到最优的参数模型使得能匹配上的特征点对最多,记为模型M,剔除与模型M的误差大于阈值的点对。
6.根据权利要求1所述的拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,其特征在于:在步骤S2中,以连续图像帧i的时间刻度为基准,将双目相机图像帧与IMU测量帧进行帧对齐,并在对齐帧之间进行IMU预积分,有如下IMU预积分形式:
7.根据权利要求1所述的拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,其特征在于:在步骤S3中,通过VINS-Mono算法中基于滑动窗口的后端优化方法来减少VIO算法框架中位姿估计的累积误差,通过非线性优化,用紧耦合的方法构建滑动窗口误差函数E,最小化边缘化先验信息、IMU测量残差和视觉测量残差,三种残差均用马氏距离来表示,即:
其中,
8.根据权利要求1所述的拒止环境下基于视觉惯性的四旋翼无人机自主定位方法,其特征在于:在步骤S4中,基于普通BRIEF描述符的DBoW词袋模型进行回环检测,根据BRIEF描述符创建二进制特征向量,用以描述特征点周围的像素信息,并将其转换成词袋向量,计算双目相机当前帧与词袋模型的相似度,通过外观特征验证和几何一致性验证判断两个不同图像位置的角点是否相似,完成回环检测,重新估计四旋翼无人机位姿。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211348263.2A CN115406447B (zh) | 2022-10-31 | 2022-10-31 | 拒止环境下基于视觉惯性的四旋翼无人机自主定位方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211348263.2A CN115406447B (zh) | 2022-10-31 | 2022-10-31 | 拒止环境下基于视觉惯性的四旋翼无人机自主定位方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN115406447A true CN115406447A (zh) | 2022-11-29 |
CN115406447B CN115406447B (zh) | 2023-03-24 |
Family
ID=84167777
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202211348263.2A Active CN115406447B (zh) | 2022-10-31 | 2022-10-31 | 拒止环境下基于视觉惯性的四旋翼无人机自主定位方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN115406447B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115523920A (zh) * | 2022-11-30 | 2022-12-27 | 西北工业大学 | 一种基于视觉惯性gnss紧耦合的无缝定位方法 |
CN116442248A (zh) * | 2023-06-19 | 2023-07-18 | 山东工程职业技术大学 | 一种基于目标检测的机器人视觉定位模块及方法 |
CN117647263A (zh) * | 2023-12-06 | 2024-03-05 | 中山大学 | 基于非线性优化的单光子相机视觉惯性里程计方法及系统 |
CN117726678A (zh) * | 2023-12-12 | 2024-03-19 | 中山大学·深圳 | 无人系统位姿估计方法、装置、计算机设备和存储介质 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108051002A (zh) * | 2017-12-04 | 2018-05-18 | 上海文什数据科技有限公司 | 基于惯性测量辅助视觉的运输车空间定位方法及系统 |
CN109520497A (zh) * | 2018-10-19 | 2019-03-26 | 天津大学 | 基于视觉和imu的无人机自主定位方法 |
CN109900265A (zh) * | 2019-03-15 | 2019-06-18 | 武汉大学 | 一种camera/mems辅助北斗的机器人定位算法 |
CN112233177A (zh) * | 2020-10-10 | 2021-01-15 | 中国安全生产科学研究院 | 一种无人机位姿估计方法及系统 |
CN112484725A (zh) * | 2020-11-23 | 2021-03-12 | 吉林大学 | 一种基于多传感器融合的智能汽车高精度定位与时空态势安全方法 |
CN112509044A (zh) * | 2020-12-02 | 2021-03-16 | 重庆邮电大学 | 一种基于点线特征融合的双目视觉slam方法 |
CN115112123A (zh) * | 2022-06-27 | 2022-09-27 | 华东理工大学 | 基于视觉-imu融合的多移动机器人协同定位方法及系统 |
-
2022
- 2022-10-31 CN CN202211348263.2A patent/CN115406447B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108051002A (zh) * | 2017-12-04 | 2018-05-18 | 上海文什数据科技有限公司 | 基于惯性测量辅助视觉的运输车空间定位方法及系统 |
CN109520497A (zh) * | 2018-10-19 | 2019-03-26 | 天津大学 | 基于视觉和imu的无人机自主定位方法 |
CN109900265A (zh) * | 2019-03-15 | 2019-06-18 | 武汉大学 | 一种camera/mems辅助北斗的机器人定位算法 |
CN112233177A (zh) * | 2020-10-10 | 2021-01-15 | 中国安全生产科学研究院 | 一种无人机位姿估计方法及系统 |
CN112484725A (zh) * | 2020-11-23 | 2021-03-12 | 吉林大学 | 一种基于多传感器融合的智能汽车高精度定位与时空态势安全方法 |
CN112509044A (zh) * | 2020-12-02 | 2021-03-16 | 重庆邮电大学 | 一种基于点线特征融合的双目视觉slam方法 |
CN115112123A (zh) * | 2022-06-27 | 2022-09-27 | 华东理工大学 | 基于视觉-imu融合的多移动机器人协同定位方法及系统 |
Non-Patent Citations (3)
Title |
---|
刘全攀 等: "《基于双目视觉-惯性导航的轻型无人机导航算法》", 《兵工学报》 * |
杨震洲: "《无人机自主定位与建图方法研究及应用》", 《万方数据知识服务平台》 * |
陈文杰: "《无人机室内自主定位与穿越多重门框关键技术研究》", 《中国优秀硕士学位论文全文数据库 工程科技II辑》 * |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115523920A (zh) * | 2022-11-30 | 2022-12-27 | 西北工业大学 | 一种基于视觉惯性gnss紧耦合的无缝定位方法 |
CN115523920B (zh) * | 2022-11-30 | 2023-03-10 | 西北工业大学 | 一种基于视觉惯性gnss紧耦合的无缝定位方法 |
CN116442248A (zh) * | 2023-06-19 | 2023-07-18 | 山东工程职业技术大学 | 一种基于目标检测的机器人视觉定位模块及方法 |
CN117647263A (zh) * | 2023-12-06 | 2024-03-05 | 中山大学 | 基于非线性优化的单光子相机视觉惯性里程计方法及系统 |
CN117647263B (zh) * | 2023-12-06 | 2024-07-09 | 中山大学 | 基于非线性优化的单光子相机视觉惯性里程计方法及系统 |
CN117726678A (zh) * | 2023-12-12 | 2024-03-19 | 中山大学·深圳 | 无人系统位姿估计方法、装置、计算机设备和存储介质 |
Also Published As
Publication number | Publication date |
---|---|
CN115406447B (zh) | 2023-03-24 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN115406447B (zh) | 拒止环境下基于视觉惯性的四旋翼无人机自主定位方法 | |
CN109520497B (zh) | 基于视觉和imu的无人机自主定位方法 | |
CN109211241B (zh) | 基于视觉slam的无人机自主定位方法 | |
US10665115B2 (en) | Controlling unmanned aerial vehicles to avoid obstacle collision | |
KR101725060B1 (ko) | 그래디언트 기반 특징점을 이용한 이동 로봇의 위치를 인식하기 위한 장치 및 그 방법 | |
CN108230361B (zh) | 用无人机探测器和追踪器融合来增强目标追踪方法及系统 | |
KR101708659B1 (ko) | 이동 로봇의 맵을 업데이트하기 위한 장치 및 그 방법 | |
Panahandeh et al. | Vision-aided inertial navigation based on ground plane feature detection | |
KR101784183B1 (ko) | ADoG 기반 특징점을 이용한 이동 로봇의 위치를 인식하기 위한 장치 및 그 방법 | |
CN111288989B (zh) | 一种小型无人机视觉定位方法 | |
CN112219087A (zh) | 位姿预测方法、地图构建方法、可移动平台及存储介质 | |
CN112240768A (zh) | 基于Runge-Kutta4改进预积分的视觉惯导融合SLAM方法 | |
CN110726406A (zh) | 一种改进的非线性优化单目惯导slam的方法 | |
Zhang et al. | Vision-aided localization for ground robots | |
KR20200075727A (ko) | 깊이 맵 산출 방법 및 장치 | |
CN114001733B (zh) | 一种基于地图的一致性高效视觉惯性定位算法 | |
US20160055646A1 (en) | Method for estimating the angular deviation of a mobile element relative to a reference direction | |
CN112802096A (zh) | 实时定位和建图的实现装置和方法 | |
Mostafa et al. | A smart hybrid vision aided inertial navigation system approach for UAVs in a GNSS denied environment | |
CN108827287B (zh) | 一种复杂环境下的鲁棒视觉slam系统 | |
CN112945233B (zh) | 一种全局无漂移的自主机器人同时定位与地图构建方法 | |
CN112731503A (zh) | 一种基于前端紧耦合的位姿估计方法及系统 | |
WO2023030062A1 (zh) | 一种无人机的飞行控制方法、装置、设备、介质及程序 | |
Han | Optical flow/ins navigation system in four-rotor | |
Yuan et al. | LIWO: LiDAR-Inertial-Wheel Odometry |
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 |