CN114723902A - 基于概率膜计算的煤矿井下无人直升机slam - Google Patents

基于概率膜计算的煤矿井下无人直升机slam Download PDF

Info

Publication number
CN114723902A
CN114723902A CN202210272080.0A CN202210272080A CN114723902A CN 114723902 A CN114723902 A CN 114723902A CN 202210272080 A CN202210272080 A CN 202210272080A CN 114723902 A CN114723902 A CN 114723902A
Authority
CN
China
Prior art keywords
helicopter
matrix
slam
probability
unmanned helicopter
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.)
Pending
Application number
CN202210272080.0A
Other languages
English (en)
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.)
Anhui University of Science and Technology
Original Assignee
Anhui University of Science and Technology
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 Anhui University of Science and Technology filed Critical Anhui University of Science and Technology
Priority to CN202210272080.0A priority Critical patent/CN114723902A/zh
Publication of CN114723902A publication Critical patent/CN114723902A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05Geographic models
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/18Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T19/00Manipulating 3D models or images for computer graphics
    • G06T19/003Navigation within 3D models or images
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/20Analysis of motion
    • G06T7/246Analysis of motion using feature-based methods, e.g. the tracking of corners or segments
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/20Analysis of motion
    • G06T7/277Analysis of motion involving stochastic approaches, e.g. using Kalman filters
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/70Determining position or orientation of objects or cameras
    • G06T7/73Determining position or orientation of objects or cameras using feature-based methods

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Databases & Information Systems (AREA)
  • Multimedia (AREA)
  • Geometry (AREA)
  • Algebra (AREA)
  • Computer Graphics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Operations Research (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Computer Hardware Design (AREA)
  • Computing Systems (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明涉及无人机监测技术领域,具体涉及基于概率膜计算的煤矿井下无人直升机SLAM。本发明在分析得出LiDAR、MU和深度相机传感器数学模型的基础上,设计了概率膜系统模型和膜算法,并对建立的模型和算法进行了验证,通过理论分析并与实验相结合,将本发明提出的方法与单一、两种和三种传感器融合的建图方法进行比较,验证了概率膜计算在建图效果上具备很好的性能。

Description

基于概率膜计算的煤矿井下无人直升机SLAM
技术领域
本发明涉及无人机监测技术领域,具体涉及基于概率膜计算的煤矿井下无人直升机SLAM。
背景技术
在未知环境下,无人机可以对自身位置和姿态进行估计,是能够有效完成任务的关键技术。当无人机在低空,低速行驶时,使用全球定位系统(Global PositioningSystem,GPS)和惯性导航系统(Inertial Navigation System,INS)等方法可以获得飞行所需的位置、速度等信息,满足在城市中的飞行需求。但是GPS和INS存在使用的局限性,GPS的误差范围为3-10米,而且在室内环境和高山、峡谷等恶劣环境下由于信号弱导致无法正常工作,尤其是近年来无人机发展迅速在开采领域也取得了一定成就,尤其是露天矿场和地下掘进巷道,无人机具有速度快、精度高等优点,可以在复杂环境下快速完成数据采集,地形勘探,但是由于矿场环境复杂,传统的定位装置精度难以保证,因此设计一种多数据融合,可利用无人机自带的传感器进行定位和飞行姿态调整的方法就显得尤为重要。
发明内容
针对上述技术背景提到的不足,本发明的目的在于提供基于概率膜计算的煤矿井下无人直升机SLAM。
本发明的目的可以通过以下技术方案实现:
基于概率膜计算的煤矿井下无人直升机SLAM,其主要包括传感器模型、基于概率稀疏扩展SLAM、概率膜计算SLAM,所述传感器模型是在分析各传感器工作原理的基础上,给出相应传感器数学模型;其次,设计膜内概率膜算法,构建概率膜系统,基于稀疏滤波实现概率SLAM地图构建。
进一步的,所述传感器模型包括惯性测量单元、激光雷达、深度相机,惯性测量单元由三个线性、正交的加速度计和陀螺仪构成,可测量正交坐标轴下的角速率和加速度;激光雷达作为利用激光脉冲实现位置测量的传感器,通过控制激光束反射角度来描述方位角及俯仰角,深度相机可测量像素的深度。
进一步的,所述基于概率稀疏扩展SLAM,具体扩展方法如下:
当无人直升机不清楚自身位置,无法获取环境地图,所有的数据集中于测量数据和控制数据,从概率学角度,SLAM问题被划分为全SLAM和在线SLAM,在全SLAM中,除了直升机位姿,还需要计算路径和后验地图,时间t时刻的概率可表达为:
p(x1~t,m|y1~t,u1~t) (13)
在线SLAM,考虑即时位姿和地图后验,时间上的概率表达为:
p(xt,m|y1~t,u1~t) (14)
式(13)和式(14)中,xt表示t时刻直升机位姿,m表示地图,y1~t,u1~t为测量和控制数据,其分别表示的模型如图4所示,
在实际应用过程中,对全SLAM的过去状态积分可实现在线SLAM,积分呈现连续状态,依次进行,如式(15)所示,
p(xt,m|y1~t,u1~t)=∫∫···∫p(x1~t,m|y1~t,u1~t)dx1dx2···dxt-1 (15)
基于连续和离散的问题,连续问题包含地图中无人机位姿和定位,在特征表示过程中物体以信标等表示,而离散特性与一致性关联,检测到物体时,算法会计算前后检测到物体之间的关联,表现出离散特征,即”0”或”1”状态,因此,一致性的明确很必要,结合式(13)和式(14),加入一致性变量的在线和全SLAM如下:
Figure BDA0003553918260000031
式(16)中,ft为一致性变量下对应的矢量,
稀疏扩展滤波:与扩展卡尔曼滤波相比,出于在线运行和计算高效率目的,稀疏扩展滤波表征了信息的高效性,继承无人直升机位姿和地图后验证,即通过非零元素相连接的维持稀疏矩阵,稀疏扩展滤波的计算过程包含测量更新、运动更新、稀疏化和估计四部分,
运动更新使用矩阵的稀疏性,在时间上不依赖地图的规模,把信息矩阵Ω和矢量ζ通过更新来完成对控制的处理,根据卡尔曼滤波给出:
Figure BDA0003553918260000032
Figure BDA0003553918260000033
式(17)和式(18)中,Σ表示协方差矩阵,Fx表示直升机状态矢量矩阵,Gt表示雅可比矩阵表示对时间t的导数,
Figure BDA0003553918260000037
表示时间t时刻的估计均值,其中Gt,Fx,δ的表示如下:
Figure BDA0003553918260000034
Figure BDA0003553918260000035
Figure BDA0003553918260000036
根据式(18)和式(19)可推导出:
Figure BDA0003553918260000041
式(22)中信息矩阵Ω维数随机且在有限时间实现,假设信息矩阵Ω稀疏,则更新效率增强,定义:
Figure BDA0003553918260000042
将式(23)代入式(22)可得:
Figure BDA0003553918260000043
由矩阵逆引理可进一步得到:
Figure BDA0003553918260000044
假设有限时间内根据Ω计算得到Φt,则存在时间有限条件下的计算可行性,利用直升机位姿和地图特征的矩阵元素(非零),稀疏条件下不依赖Ω的大小,考虑Gt的逆
Figure BDA0003553918260000045
可以由以下计算得出:
Figure BDA0003553918260000046
式(26)中,对应地图特征元素非零,
测量更新考虑直升机飞行过程中的滤波更新,通过扩展卡尔曼滤波实现:
Figure BDA0003553918260000047
其中,Qt为噪声协方差矩阵。
进一步的,所述稀疏扩展滤波对信息矩阵Ω稀疏化是必要的,通过稀疏化的表征,保障其后验分布处于稀疏状态,基于此,剔除直升机位姿和地图特征之间的关联,进一步限制了特征间的数量,为实现上述思想,引入两种新的连接,首先,依靠非活动连接激活此特征,在无人直升机位姿和特征之间引入新连接;其次,直升机运动引入活动特征的两个新连接,限制活动特征的数量,避免出现两个非稀疏边界,此时的稀疏通过少的活动特征获取,
稀疏化定义过程中将特征集合分为3个子集(不相交):
m=m0+m1+m2 (28)
式(28),m1为活动继续的特征集合,m0为即将激励的活动特征,m2为非活动特征,在稀疏化步骤中继续非活动状态,同时删除直升机位姿和m0之间的连接,将稀疏化引入后验,因m1和m0包含当前所有特征,后验p(yt|z1~t,u1~t,f1~t)可表征如下:
Figure BDA0003553918260000051
式(29)中,xt不依赖非活动特征m2的前提是m0和m1已知,故m2可取随机值,利用通用项稀疏化规约,降低对m0依赖,取m2=0,则:
Figure BDA0003553918260000052
式(30)中,z1~t表示直到t时刻的测量值,u1~t在线SLAM控制量。
进一步的,所述概率膜计算SLAM具体步骤如下:
确定无人直升机的实时位置,膜控制器开始执行的每个周期都接收机载传感器数据(x,y,θ)T,位置更新输出数据(x′,y′,θ′)T,根据膜计算的分布并行特征,建立度为4的概率膜系统,
Π=(M,μ,w1,w2,w3,R,{cr}r∈R) (31)
1)M={xij,yij,θij,Error:i,j∈[1,2]};M表示以xij,yijij,Error为对象构成的集合,其中:xij表示无人直升机不同位置的横坐标,yij表示无人直升机不同位置的纵坐标,θij表示无人直升机对应xij,yij坐标下的方位角,Error为无人直升机运动产生的误差;
2)μ=[[[]2[]3]4]1,μ表示无人直升机组织膜系统膜结构,其中:[]表示单个膜,1,2,3,4分别表示膜的标记,即第i个膜(i=1,2,3,4);
3)w1=p(xt|ut,xt-1),w1表示解算密度函数,其中:xt为无人直升机t时刻姿态,ut为无人直升机控制输入,xt+1为t+1时刻无人直升机位姿;
4)w2=p(xt′|ut′,xt-1′),w2表示噪声干扰后密度计算,其中:xt′为无人直升机在噪声干扰后t时刻姿态,ut′为噪声干扰控制输入,xt+1′为无人直升机在噪声干扰后t+1时刻姿态;
5)
Figure BDA0003553918260000061
w3表示理想密度评价,其中:
Figure BDA0003553918260000062
为无人直升机t时刻理想姿态,
Figure BDA0003553918260000063
为无人直升机理想位姿对应的控制输入,
Figure BDA0003553918260000064
为t+1时刻无人直升机理想位姿;
6)R为规则集,转移过程中保持能量守恒;
7)cr表示基于R规则集下的对象产生进化的概率,
概率膜计算模型:结合概率膜计算方法,对排除m2的所有变量服从p(xt,m0,m1|m2=0)分布计算矩阵提取状态变量子矩阵:
Figure BDA0003553918260000065
式(32)中,Fx表示直升机状态矢量矩阵,m0为活动继续的特征集合,m1为即将激励的活动特征,根据矩阵引理,对p(m1|m2,z1~t,u1~t,f1~t)和p(xt,m0,m1,m2|z1~t,u1~t,f1~t)矩阵信息进行如下定义:
Figure BDA0003553918260000066
式(33)中的F为把全部状态投影成包含变量子集状态的投影阵,同理,p(m0,m1,m2|z1~t,u1~t,f1~t)可转换成矩阵:
Figure BDA0003553918260000067
将式(32)和式(34)合并可得到信息矩阵
Figure BDA0003553918260000071
Figure BDA0003553918260000072
信息矢量分别为:
Figure BDA0003553918260000073
Figure BDA0003553918260000074
进一步的,所述膜计算模型中个算法执行过程如下:
传感数据更新算法:
Figure BDA0003553918260000075
直升机运动更新算法:
Figure BDA0003553918260000076
Figure BDA0003553918260000081
直升机状态估计更新算法:
Figure BDA0003553918260000082
根据无人直升机后验输入Ωt和ξt,位移d和角向量α确定坐标系下的位置,同时将位移和向量映射到机体坐标系,基于旋转和平移,直升机位姿和地图特征分别表示为:
Figure BDA0003553918260000091
Figure BDA0003553918260000092
式(39)和式(38)包含矩阵和向量的旋转与平移,下一步给出两式的证明,定义
Figure BDA0003553918260000093
δ:
Figure BDA0003553918260000094
δ=[dx dy α]T (40)
将式(39)和式(38)进行状态向量推导得:
Figure BDA0003553918260000095
根据空间坐标变换相似性原理,进一步通过平移和旋转,信息矩阵和向量定义直升机随时间的后验:
Figure BDA0003553918260000096
由于
Figure BDA0003553918260000097
Figure BDA0003553918260000098
有:
Figure BDA0003553918260000099
融合过程中存在数据等价问题,通过增加约束进一步控制,增加地图特征惩罚矩阵C,其值愈大约束愈强,通过联合后验概率实现地图融合,融合算法如下:
Figure BDA0003553918260000101
本发明的有益效果:
本发明在分析得出LiDAR、MU和深度相机传感器数学模型的基础上,设计了概率膜系统模型和膜算法,并对建立的模型和算法进行了验证,通过理论分析并与实验相结合,将本发明提出的方法与单一、两种和三种传感器融合的建图方法进行比较,验证了概率膜计算在建图效果上具备很好的性能。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图;
图1是概率膜系统SLAM流程图;
图2是激光雷达示意图;
图3是深度相机模型;
图4是加固连接装置的结构示意图;
图5是膜计算模型P-Lingua文件框架;
图6是不同SLAM算法的角度误差比较图;
图7是不同SLAM算法绝对误差比较图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例,都属于本发明保护的范围。
基于概率膜计算的煤矿井下无人直升机SLAM(Probability Membrane ComputingSLAM,PMC-SLAM),其主要包括传感器模型、基于概率稀疏扩展SLAM、概率膜计算SLAM。
其中传感器模型包括惯性测量单元(IMU)、激光雷达(LiDAR)、深度相机(Camera)。
惯性测量单元由三个线性、正交的加速度计和陀螺仪构成,可测量正交坐标轴下的角速率和加速度。对惯性测量下的测量值基于传感器自身坐标系,非无人机机体坐标系。假设机体状态向量位姿:{cvi,rvi},角速度:wvi,加速度:
Figure BDA0003553918260000111
进一步假设传感器和机体坐标系位置关系确定,rsv表示旋转,csv表示平移。角速度w与在传感器坐标系下机体的速率关系表示如下:
w=wvicsv (1)
根据加速度测量原理,加速度测量值表示为:
Figure BDA0003553918260000121
Figure BDA0003553918260000122
表示传感器坐标原点s的加速度,g为重力加速度。由于公式2表示的值并不是机体坐标系下的值,需要将传感器和机体坐标转换,因为:
Figure BDA0003553918260000123
利用泊松公式对式(3)两次求导可得:
Figure BDA0003553918260000124
式(4)右侧由初始参数构成,
Figure BDA0003553918260000125
为传感器原点的关系加速度。将式(4)代入式(2)得:
Figure BDA0003553918260000126
式(5)中,a为加速度计模型,将式(1)和式(5)模型进一步表示成矩阵形式下的IMU传感器模型为:
Figure BDA0003553918260000127
式(6)中,g表示惯性坐标系下重力加速度。
激光雷达作为利用激光脉冲实现位置测量的传感器,通过控制激光束反射角度来描述方位角及俯仰角,距离通过飞行时间测得,其表示如图2所示。图2中,P点坐标表示为:
p=rps=[x,y,z]T (7)
以矩阵形式表示为:
p=cT(α)cT(ε)[r,0,0]T (8)
式(8)中,r表示距离,ε表示俯仰角,α表示方位角,c为沿坐标轴的旋转矩阵。为进一步表示式(7)和式(8)之间的关系,根据旋转求逆可得:
Figure BDA0003553918260000131
式(9)表明了激光雷达观测点位置(x,y,z)与距离r,俯仰角ε和方位α之间的映射关系。
深度相机可测量像素的深度,测量方法如下:P点坐标(x,y,z),左右相机的模型分别定义为:
Figure BDA0003553918260000132
式(10)中,K为相机参数,fu为水平像素焦距,fv为垂直焦距。由于两相机参数矩阵相同,当左相机获取观测值时,右侧相机沿图3中的实直线便可获取右侧相机观测值(垂直),因此,式(10)进一步组合可得到整体相机的模型为:
Figure BDA0003553918260000133
据式(11),由于参数矩阵N存在相同的两行,所以N不可逆,约束条件下左右两相机存在cl.r=1的姿态关系,可得到vl=vr,其中vl和vr为坐标点x到像素坐标下的映射。
信息融合方法,本发明以三种不同的传感器信息为基础分别构建三张地图,再对构建的地图进行融合,由于多传感器之间的数据不相互依赖,假设3种传感器生成的地图表示为map3,利用德摩根定律分解可得:
p(map)=1-Π(1-p(map3)) (12)
如果式(12)中的激光雷达、IMU或深度相机传感器构建的地图被划分单元占用,则融合地图呈现相同的占用状态。
基于概率稀疏扩展SLAM:
当无人直升机不清楚自身位置,无法获取环境地图,所有的数据集中于测量数据和控制数据。从概率学角度,SLAM问题被划分为全SLAM和在线SLAM。在全SLAM中,除了直升机位姿,还需要计算路径和后验地图,时间t时刻的概率可表达为:
p(x1~t,m|y1~t,u1~t) (13)在线SLAM,考虑即时位姿和地图后验,时间上的概率表达为:
p(xt,m|y1~t,u1~t) (14)
式(13)和式(14)中,xt表示t时刻直升机位姿,m表示地图,y1~t,u1~t为测量和控制数据,其分别表示的模型如图4所示。
在实际应用过程中,对全SLAM的过去状态积分可实现在线SLAM,积分呈现连续状态,依次进行,如式(15)所示。
p(xt,m|y1~t,u1~t)=∫∫···∫p(x1~t,m|y1~t,u1~t)dx1dx2···dxt-1 (15)
基于连续和离散的问题,连续问题包含地图中无人机位姿和定位,在特征表示过程中物体以信标等表示。而离散特性与一致性关联,检测到物体时,算法会计算前后检测到物体之间的关联,表现出离散特征,即”0”或”1”状态。因此,一致性的明确很必要,结合式(13)和式(14),加入一致性变量的在线和全SLAM如下:
Figure BDA0003553918260000141
式(16)中,ft为一致性变量下对应的矢量。
稀疏扩展滤波:与扩展卡尔曼滤波相比,出于在线运行和计算高效率目的,稀疏扩展滤波表征了信息的高效性,继承无人直升机位姿和地图后验证,即通过非零元素相连接的维持稀疏矩阵。稀疏扩展滤波的计算过程包含测量更新、运动更新、稀疏化和估计四部分。
运动更新使用矩阵的稀疏性,在时间上不依赖地图的规模,把信息矩阵Ω和矢量ζ通过更新来完成对控制的处理,根据卡尔曼滤波给出:
Figure BDA0003553918260000151
Figure BDA0003553918260000152
式(17)和式(18)中,Σ表示协方差矩阵,Fx表示直升机状态矢量矩阵,Gt表示雅可比矩阵表示对时间t的导数,
Figure BDA0003553918260000153
表示时间t时刻的估计均值。其中Gt,Fx,δ的表示如下:
Figure BDA0003553918260000154
Figure BDA0003553918260000155
Figure BDA0003553918260000156
根据式(18)和式(19)可推导出:
Figure BDA0003553918260000157
式(22)中信息矩阵Ω维数随机且在有限时间实现,假设信息矩阵Ω稀疏,则更新效率增强。定义:
Figure BDA0003553918260000161
将式(23)代入式(22)可得:
Figure BDA0003553918260000162
由矩阵逆引理可进一步得到:
Figure BDA0003553918260000163
假设有限时间内根据Ω计算得到Φt,则存在时间有限条件下的计算可行性。利用直升机位姿和地图特征的矩阵元素(非零),稀疏条件下不依赖Ω的大小,考虑Gt的逆,
Figure BDA0003553918260000164
可以由以下计算得出:
Figure BDA0003553918260000165
式(26)中,对应地图特征元素非零。
测量更新考虑直升机飞行过程中的滤波更新,通过扩展卡尔曼滤波实现:
Figure BDA0003553918260000166
其中,Qt为噪声协方差矩阵。
稀疏扩展滤波对信息矩阵Ω稀疏化是必要的,通过稀疏化的表征,保障其后验分布处于稀疏状态,基于此,剔除直升机位姿和地图特征之间的关联,进一步限制了特征间的数量。为实现上述思想,引入两种新的连接。首先,依靠非活动连接激活此特征,在无人直升机位姿和特征之间引入新连接;其次,直升机运动引入活动特征的两个新连接,限制活动特征的数量,避免出现两个非稀疏边界,此时的稀疏通过少的活动特征获取。
稀疏化定义过程中将特征集合分为3个子集(不相交):
m=m0+m1+m2 (28)
式(28),m1为活动继续的特征集合,m0为即将激励的活动特征,m2为非活动特征,在稀疏化步骤中继续非活动状态,同时删除直升机位姿和m0之间的连接。将稀疏化引入后验,因m1和m0包含当前所有特征,后验p(yt|z1~t,u1~t,f1~t)可表征如下:
Figure BDA0003553918260000171
式(29)中,xt不依赖非活动特征m2的前提是m0和m1已知,故m2可取随机值。利用通用项稀疏化规约,降低对m0依赖,取m2=0,则:
Figure BDA0003553918260000172
式(30)中,z1~t表示直到t时刻的测量值,u1~t在线SLAM控制量。
概率膜计算SLAM:在概率膜计算中,规则以概率方式执行,可以快速处理数据,结合膜计算的并行、分布式特征,构成一个独立、协作的数据处理膜计算模型。采用概率膜系统对煤矿井下无人直升机环境感知地图构建流程如图1所示。
确定无人直升机的实时位置,膜控制器开始执行的每个周期都接收机载传感器数据(x,y,θ)T,位置更新输出数据(x′,y′,θ′)T。根据膜计算的分布并行特征,建立度为4的概率膜系统。
Π=(M,μ,w1,w2,w3,R,{cr}r∈R) (31)
式(31)中:
1)M={xij,yij,θij,Error:i,j∈[1,2]};M表示以xij,yijij,Error为对象构成的集合。其中:xij表示无人直升机不同位置的横坐标,yij表示无人直升机不同位置的纵坐标,θij表示无人直升机对应xij,yij坐标下的方位角,Error为无人直升机运动产生的误差;
2)μ=[[[]2[]3]4]1,μ表示无人直升机组织膜系统膜结构,其中:[]表示单个膜,1,2,3,4分别表示膜的标记,即第i个膜(i=1,2,3,4);
3)w1=p(xt|ut,xt-1),w1表示解算密度函数。其中:xt为无人直升机t时刻姿态,ut为无人直升机控制输入,xt+1为t+1时刻无人直升机位姿;
4)w2=p(xt′|ut′,xt-1′),w2表示噪声干扰后密度计算。其中:xt′为无人直升机在噪声干扰后t时刻姿态,ut′为噪声干扰控制输入,xt+1′为无人直升机在噪声干扰后t+1时刻姿态;
5)
Figure BDA0003553918260000181
w3表示理想密度评价。其中:
Figure BDA0003553918260000182
为无人直升机t时刻理想姿态,
Figure BDA0003553918260000183
为无人直升机理想位姿对应的控制输入,
Figure BDA0003553918260000184
为t+1时刻无人直升机理想位姿;
6)R为规则集,转移过程中保持能量守恒;
7)cr表示基于R规则集下的对象产生进化的概率。
概率膜计算模型:结合概率膜计算方法,对排除m2的所有变量服从p(xt,m0,m1|m2=0)分布计算矩阵提取状态变量子矩阵:
Figure BDA0003553918260000185
式(32)中,Fx表示直升机状态矢量矩阵,m0为活动继续的特征集合,m1为即将激励的活动特征,根据矩阵引理,对p(m1|m2,z1~t,u1~t,f1~t)和p(xt,m0,m1,m2|z1~t,u1~t,f1~t)矩阵信息进行如下定义:
Figure BDA0003553918260000186
式(33)中的F为把全部状态投影成包含变量子集状态的投影阵。同理,p(m0,m1,m2|z1~t,u1~t,f1~t)可转换成矩阵:
Figure BDA0003553918260000187
将式(32)和式(34)合并可得到信息矩阵
Figure BDA0003553918260000191
Figure BDA0003553918260000192
信息矢量分别为:
Figure BDA0003553918260000193
Figure BDA0003553918260000194
由于概率膜算法循环时间上的有限性,不依赖地图本身规模,同时把控制加入计算估计,通过高效率的计算特征,在建图过程中,矩阵(矢量)中的元素更新过程中表征直升机位置和特征向量,概率膜计算模型p-Lingua文件框架如图5所示:
图5膜计算模型中各算法执行过程如下:
传感数据更新算法:
Figure BDA0003553918260000195
Figure BDA0003553918260000201
直升机运动更新算法:
Figure BDA0003553918260000202
直升机状态估计更新算法:
Figure BDA0003553918260000203
Figure BDA0003553918260000211
根据无人直升机后验输入Ωt和ξt,位移d和角向量α确定坐标系下的位置,同时将位移和向量映射到机体坐标系,基于旋转和平移,直升机位姿和地图特征分别表示为:
Figure BDA0003553918260000212
Figure BDA0003553918260000213
式(39)和式(38)包含矩阵和向量的旋转与平移,下一步给出两式的证明。定义
Figure BDA0003553918260000214
δ:
Figure BDA0003553918260000215
δ=[dx dy α]T (40)
将式(39)和式(38)进行状态向量推导得:
Figure BDA0003553918260000216
根据空间坐标变换相似性原理,进一步通过平移和旋转,信息矩阵和向量定义直升机随时间的后验:
Figure BDA0003553918260000221
由于
Figure BDA0003553918260000222
Figure BDA0003553918260000223
有:
Figure BDA0003553918260000224
融合过程中存在数据等价问题,通过增加约束进一步控制,增加地图特征惩罚矩阵C,其值愈大约束愈强。通过联合后验概率实现地图融合,融合算法如下:
Figure BDA0003553918260000225
算法4的执行过程:无人机坐标系间的相对位姿由[dx dy α]T确定,维持算法稀疏的同时,包含信息矩阵和向量的局部旋转和平移;通过构建联合后验地图实现地图融合,该步骤包含不同地图中对应的特征;执行等价约束,对于相同的两个特征,添加信息矩阵中对这两个特征的连接。
利用无人机搭载三种传感器进行信息采集,并将获得的数据进行建图分析。如图6、7所示,其中图6、图7中最下方的线为基于膜概率计算的三数据融合建图后无人机飞行的误差,可以发现利用惯性测量单元(IMU)、激光雷达(LiDAR)和深度相机(Camera)进行膜概率计算的三数据融合建图效果明显优于单一传感器或两两组合的传感器数据建图,膜概率计算的三数据融合建图得到的角度误差和绝对误差随着均明显减小,且随着无人机飞行时间和飞行里程数的增加,角度误差和绝对误差均是最小。
在本说明书的描述中,参考术语“一个实施例”、“示例”、“具体示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中,对上述术语的示意性表述不一定指的是相同的实施例或示例。而且,描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解,本发明不受上述实施例的限制,上述实施例和说明书中描述的只是说明本发明的原理,在不脱离本发明精神和范围的前提下,本发明还会有各种变化和改进,这些变化和改进都落入要求保护的本发明范围内。

Claims (6)

1.基于概率膜计算的煤矿井下无人直升机SLAM,其主要包括传感器模型、基于概率稀疏扩展SLAM、概率膜计算SLAM,其特征在于,所述传感器模型是在分析各传感器工作原理的基础上,给出相应传感器数学模型;其次,设计膜内概率膜算法,构建概率膜系统,基于稀疏滤波实现概率SLAM地图构建。
2.根据权利要求1所述的基于概率膜计算的煤矿井下无人直升机SLAM,其特征在于,所述传感器模型包括惯性测量单元、激光雷达、深度相机,惯性测量单元由三个线性、正交的加速度计和陀螺仪构成,可测量正交坐标轴下的角速率和加速度;激光雷达作为利用激光脉冲实现位置测量的传感器,通过控制激光束反射角度来描述方位角及俯仰角,深度相机可测量像素的深度。
3.根据权利要求1所述的基于概率膜计算的煤矿井下无人直升机SLAM,其特征在于,所述基于概率稀疏扩展SLAM,具体扩展方法如下:
当无人直升机不清楚自身位置,无法获取环境地图,所有的数据集中于测量数据和控制数据,从概率学角度,SLAM问题被划分为全SLAM和在线SLAM,在全SLAM中,除了直升机位姿,还需要计算路径和后验地图,时间t时刻的概率可表达为:
p(x1~t,m|y1~t,u1~t) (13)
在线SLAM,考虑即时位姿和地图后验,时间上的概率表达为:
p(xt,m|y1~t,u1~t) (14)
式(13)和式(14)中,xt表示t时刻直升机位姿,m表示地图,y1~t,u1~t为测量和控制数据,其分别表示的模型如图4所示,
在实际应用过程中,对全SLAM的过去状态积分可实现在线SLAM,积分呈现连续状态,依次进行,如式(15)所示,
p(xt,m|y1~t,u1~t)=∫∫···∫p(x1~t,m|y1~t,u1~t)dx1dx2···dxt-1 (15)
基于连续和离散的问题,连续问题包含地图中无人机位姿和定位,在特征表示过程中物体以信标等表示,而离散特性与一致性关联,检测到物体时,算法会计算前后检测到物体之间的关联,表现出离散特征,即”0”或”1”状态,因此,一致性的明确很必要,结合式(13)和式(14),加入一致性变量的在线和全SLAM如下:
Figure FDA0003553918250000021
式(16)中,ft为一致性变量下对应的矢量,
稀疏扩展滤波:与扩展卡尔曼滤波相比,出于在线运行和计算高效率目的,稀疏扩展滤波表征了信息的高效性,继承无人直升机位姿和地图后验证,即通过非零元素相连接的维持稀疏矩阵,稀疏扩展滤波的计算过程包含测量更新、运动更新、稀疏化和估计四部分,
运动更新使用矩阵的稀疏性,在时间上不依赖地图的规模,把信息矩阵Ω和矢量ζ通过更新来完成对控制的处理,根据卡尔曼滤波给出:
Figure FDA0003553918250000022
Figure FDA0003553918250000023
式(17)和式(18)中,Σ表示协方差矩阵,Fx表示直升机状态矢量矩阵,Gt表示雅可比矩阵表示对时间t的导数,
Figure FDA0003553918250000024
表示时间t时刻的估计均值,其中Gt,Fx,δ的表示如下:
Figure FDA0003553918250000025
Figure FDA0003553918250000031
Figure FDA0003553918250000032
根据式(18)和式(19)可推导出:
Figure FDA0003553918250000033
式(22)中信息矩阵Ω维数随机且在有限时间实现,假设信息矩阵Ω稀疏,则更新效率增强,定义:
Figure FDA0003553918250000034
将式(23)代入式(22)可得:
Figure FDA0003553918250000035
由矩阵逆引理可进一步得到:
Figure FDA0003553918250000036
假设有限时间内根据Ω计算得到Φt,则存在时间有限条件下的计算可行性,利用直升机位姿和地图特征的矩阵元素(非零),稀疏条件下不依赖Ω的大小,考虑Gt的逆,
Figure FDA0003553918250000037
可以由以下计算得出:
Figure FDA0003553918250000038
式(26)中,对应地图特征元素非零,
测量更新考虑直升机飞行过程中的滤波更新,通过扩展卡尔曼滤波实现:
Figure FDA0003553918250000041
其中,Qt为噪声协方差矩阵。
4.根据权利要求3所述的基于概率膜计算的煤矿井下无人直升机SLAM,其特征在于,所述稀疏扩展滤波对信息矩阵Ω稀疏化是必要的,通过稀疏化的表征,保障其后验分布处于稀疏状态,基于此,剔除直升机位姿和地图特征之间的关联,进一步限制了特征间的数量,为实现上述思想,引入两种新的连接,首先,依靠非活动连接激活此特征,在无人直升机位姿和特征之间引入新连接;其次,直升机运动引入活动特征的两个新连接,限制活动特征的数量,避免出现两个非稀疏边界,此时的稀疏通过少的活动特征获取,
稀疏化定义过程中将特征集合分为3个子集(不相交):
m=m0+m1+m2 (28)
式(28),m1为活动继续的特征集合,m0为即将激励的活动特征,m2为非活动特征,在稀疏化步骤中继续非活动状态,同时删除直升机位姿和m0之间的连接,将稀疏化引入后验,因m1和m0包含当前所有特征,后验p(yt|z1~t,u1~t,f1~t)可表征如下:
Figure FDA0003553918250000042
式(29)中,xt不依赖非活动特征m2的前提是m0和m1已知,故m2可取随机值,利用通用项稀疏化规约,降低对m0依赖,取m2=0,则:
Figure FDA0003553918250000043
式(30)中,z1~t表示直到t时刻的测量值,u1~t在线SLAM控制量。
5.根据权利要求1所述的基于概率膜计算的煤矿井下无人直升机SLAM,其特征在于,所述概率膜计算SLAM具体步骤如下:
确定无人直升机的实时位置,膜控制器开始执行的每个周期都接收机载传感器数据(x,y,θ)T,位置更新输出数据(x′,y′,θ′)T,根据膜计算的分布并行特征,建立度为4的概率膜系统,
Π=(M,μ,w1,w2,w3,R,{Cr}r∈R) (31)
1)M={xij,yij,θij,Error:i,j∈[1,2]};M表示以xij,yijij,Error为对象构成的集合,其中:xij表示无人直升机不同位置的横坐标,yij表示无人直升机不同位置的纵坐标,θij表示无人直升机对应xij,yij坐标下的方位角,Error为无人直升机运动产生的误差;
2)μ=[[[]2[]3]4]1,μ表示无人直升机组织膜系统膜结构,其中:[]表示单个膜,1,2,3,4分别表示膜的标记,即第i个膜(i=1,2,3,4);
3)w1=p(xt|ut,xt-1),w1表示解算密度函数,其中:xt为无人直升机t时刻姿态,ut为无人直升机控制输入,xt+1为t+1时刻无人直升机位姿;
4)w2=p(xt′|ut′,xt-1′),w2表示噪声干扰后密度计算,其中:xt′为无人直升机在噪声干扰后t时刻姿态,u′t为噪声干扰控制输入,xt+1′为无人直升机在噪声干扰后t+1时刻姿态;
5)
Figure FDA0003553918250000051
w3表示理想密度评价,其中:
Figure FDA0003553918250000052
为无人直升机t时刻理想姿态,
Figure FDA0003553918250000053
为无人直升机理想位姿对应的控制输入,
Figure FDA0003553918250000054
为t+1时刻无人直升机理想位姿;
6)R为规则集,转移过程中保持能量守恒;
7)cr表示基于R规则集下的对象产生进化的概率,
概率膜计算模型:结合概率膜计算方法,对排除m2的所有变量服从p(xt,m0,m1|m2=0)分布计算矩阵提取状态变量子矩阵:
Figure FDA0003553918250000061
式(32)中,Fx表示直升机状态矢量矩阵,m0为活动继续的特征集合,m1为即将激励的活动特征,根据矩阵引理,对p(m1|m2,z1~t,u1~t,f1~t)和p(xt,m0,m1,m2|z1~t,u1~t,f1~t)矩阵信息进行如下定义:
Figure FDA0003553918250000062
式(33)中的F为把全部状态投影成包含变量子集状态的投影阵,同理,p(m0,m1,m2|z1~t,u1~t,f1~t)可转换成矩阵:
Figure FDA0003553918250000063
将式(32)和式(34)合并可得到信息矩阵
Figure FDA0003553918250000064
Figure FDA0003553918250000065
信息矢量分别为:
Figure FDA0003553918250000066
Figure FDA0003553918250000067
6.根据权利要求5所述的基于概率膜计算的煤矿井下无人直升机SLAM,其特征在于,所述膜计算模型中个算法执行过程如下:
传感数据更新算法:
Figure FDA0003553918250000068
Figure FDA0003553918250000071
直升机运动更新算法:
Figure FDA0003553918250000072
直升机状态估计更新算法:
Figure FDA0003553918250000081
根据无人直升机后验输入Ωt和ξt,位移d和角向量α确定坐标系下的位置,同时将位移和向量映射到机体坐标系,基于旋转和平移,直升机位姿和地图特征分别表示为:
Figure FDA0003553918250000082
Figure FDA0003553918250000083
式(39)和式(38)包含矩阵和向量的旋转与平移,下一步给出两式的证明,定义
Figure FDA0003553918250000084
δ:
Figure FDA0003553918250000085
δ=[dx dy α]T (40)
将式(39)和式(38)进行状态向量推导得:
Figure FDA0003553918250000091
根据空间坐标变换相似性原理,进一步通过平移和旋转,信息矩阵和向量定义直升机随时间的后验:
Figure FDA0003553918250000092
由于
Figure FDA0003553918250000093
Figure FDA0003553918250000094
有:
Figure FDA0003553918250000095
融合过程中存在数据等价问题,通过增加约束进一步控制,增加地图特征惩罚矩阵C,其值愈大约束愈强,通过联合后验概率实现地图融合,融合算法如下:
Figure FDA0003553918250000096
CN202210272080.0A 2022-03-18 2022-03-18 基于概率膜计算的煤矿井下无人直升机slam Pending CN114723902A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210272080.0A CN114723902A (zh) 2022-03-18 2022-03-18 基于概率膜计算的煤矿井下无人直升机slam

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210272080.0A CN114723902A (zh) 2022-03-18 2022-03-18 基于概率膜计算的煤矿井下无人直升机slam

Publications (1)

Publication Number Publication Date
CN114723902A true CN114723902A (zh) 2022-07-08

Family

ID=82237026

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210272080.0A Pending CN114723902A (zh) 2022-03-18 2022-03-18 基于概率膜计算的煤矿井下无人直升机slam

Country Status (1)

Country Link
CN (1) CN114723902A (zh)

Similar Documents

Publication Publication Date Title
CN112347840B (zh) 视觉传感器激光雷达融合无人机定位与建图装置和方法
US10295365B2 (en) State estimation for aerial vehicles using multi-sensor fusion
CN111156998B (zh) 一种基于rgb-d相机与imu信息融合的移动机器人定位方法
Kim et al. Autonomous airborne navigation in unknown terrain environments
Doer et al. An ekf based approach to radar inertial odometry
CN109991636A (zh) 基于gps、imu以及双目视觉的地图构建方法及系统
CN108195376B (zh) 小型无人机自主导航定位方法
CN111508282B (zh) 低空无人机农田作业飞行障碍物冲突检测方法
CN110849360B (zh) 面向多机协同编队飞行的分布式相对导航方法
CN111189442A (zh) 基于cepf的无人机多源导航信息状态预测方法
Goppert et al. Invariant Kalman filter application to optical flow based visual odometry for UAVs
CN111238469A (zh) 一种基于惯性/数据链的无人机编队相对导航方法
Park et al. Design and performance validation of integrated navigation system based on geometric range measurements and GIS map for urban aerial navigation
Campbell et al. Vision-based geolocation tracking system for uninhabited aerial vehicles
Lopes et al. Attitude determination of highly dynamic fixed-wing uavs with gps/mems-ahrs integration
Gupta et al. Terrain‐based vehicle orientation estimation combining vision and inertial measurements
Chen et al. Aerial robots on the way to underground: An experimental evaluation of VINS-mono on visual-inertial odometry camera
Wang et al. Micro aerial vehicle navigation with visual-inertial integration aided by structured light
Dubey et al. Droan—disparity-space representation for obstacle avoidance
CN112923934A (zh) 一种适用于非结构化场景中结合惯导的激光slam技术
CN114723902A (zh) 基于概率膜计算的煤矿井下无人直升机slam
Xu et al. Probabilistic membrane computing-based SLAM for patrol UAVs in coal mines
Lazarus et al. Unmanned aerial vehicle navigation and mapping
Hosen et al. A vision-aided nonlinear observer for fixed-wing UAV navigation
Azizi et al. 3D inertial algorithm of SLAM for using on UAV

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