CN108562274B - 一种基于标志器的非合作目标位姿测量方法 - Google Patents
一种基于标志器的非合作目标位姿测量方法 Download PDFInfo
- Publication number
- CN108562274B CN108562274B CN201810359727.7A CN201810359727A CN108562274B CN 108562274 B CN108562274 B CN 108562274B CN 201810359727 A CN201810359727 A CN 201810359727A CN 108562274 B CN108562274 B CN 108562274B
- Authority
- CN
- China
- Prior art keywords
- arc
- marker
- coordinate system
- ellipse
- point
- 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
Links
- 239000003550 marker Substances 0.000 title claims abstract description 132
- 238000000691 measurement method Methods 0.000 title claims abstract description 21
- 238000003032 molecular docking Methods 0.000 claims abstract description 21
- 239000011159 matrix material Substances 0.000 claims description 54
- 239000013598 vector Substances 0.000 claims description 32
- 238000001514 detection method Methods 0.000 claims description 16
- 230000009466 transformation Effects 0.000 claims description 15
- 238000013519 translation Methods 0.000 claims description 15
- 238000000034 method Methods 0.000 claims description 11
- 238000012216 screening Methods 0.000 claims description 9
- 230000003247 decreasing effect Effects 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000000605 extraction Methods 0.000 claims description 4
- 238000001914 filtration Methods 0.000 claims description 4
- 210000001503 joint Anatomy 0.000 claims description 4
- 230000003287 optical effect Effects 0.000 claims description 4
- 238000007781 pre-processing Methods 0.000 claims description 4
- 238000011156 evaluation Methods 0.000 claims description 3
- 230000001186 cumulative effect Effects 0.000 claims description 2
- 238000003708 edge detection Methods 0.000 claims description 2
- 238000003384 imaging method Methods 0.000 claims description 2
- 238000012805 post-processing Methods 0.000 claims description 2
- 238000012545 processing Methods 0.000 claims description 2
- 230000002441 reversible effect Effects 0.000 claims description 2
- 238000012163 sequencing technique Methods 0.000 claims description 2
- 229920000535 Tan II Polymers 0.000 claims 3
- FPIGOBKNDYAZTP-UHFFFAOYSA-N 1,2-epoxy-3-(4-nitrophenoxy)propane Chemical compound C1=CC([N+](=O)[O-])=CC=C1OCC1OC1 FPIGOBKNDYAZTP-UHFFFAOYSA-N 0.000 claims 2
- 108091026890 Coding region Proteins 0.000 claims 1
- 238000005259 measurement Methods 0.000 abstract description 11
- 230000007547 defect Effects 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 8
- 238000005516 engineering process Methods 0.000 description 6
- 230000000007 visual effect Effects 0.000 description 4
- 238000013459 approach Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 229910002056 binary alloy Inorganic materials 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 230000000670 limiting effect Effects 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000002829 reductive effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- 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
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/80—Analysis of captured images to determine intrinsic or extrinsic camera parameters, i.e. camera calibration
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Theoretical Computer Science (AREA)
- Multimedia (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种基于标志器的非合作目标位姿测量方法,首先通过对不同的标志器进行识别,解算标志器位姿,得到各标志器坐标系与相机坐标系之间的相对位姿;进行特征点三维坐标恢复,对标志器进行定位,解算各标志器坐标系与目标器坐标系间的位姿参数;精确估算出追踪航天器与目标航天器之间的相对位姿;通过引入预先设计的标志器可以更精确地测量位姿信息,解决非合作目标位姿测量准确度低的技术问题,同时克服了交会对接近距离阶段,追踪航天器无法获取完整目标航天器特定部件图像以识别定位的缺点。
Description
技术领域
本发明属于空间交会与对接技术中的视觉测量领域,具体涉及一种基于标志器的非合作目标位姿测量方法。
背景技术
空间交会与对接技术是载人航天领域的三大基本技术之一,是实现航天器装配、回收、补给、维修、航天员交换及营救等在轨服务的首要条件。在交会对接最终逼近阶段,一般通过视觉测量提供目标航天器的相对位姿。目前,空间交会对接任务中普遍采用的都是基于合作目标的视觉测量技术,应用相对成熟。然而,空间绝大多数在轨航天器均属于非合作目标,其中包括故障或失效卫星、太空垃圾以及非己方航天器等,这就涉及到针对非合作目标的视觉测量关键技术。
非合作目标位姿测量的研究方向之一是以非合作目标的形状特征如三角架、星箭对接环、发动机喷嘴、矩形太阳能帆板、长方体目标基座等为测量对象,建立合适的参考坐标系,实现非合作目标的位姿解算。但其面临的挑战是如何从三维特征库中建立起与二维图像特征之间的特征点对应关系,并在此基础上设计目标位姿鲁棒估计方法。
发明内容
本发明提出一种基于标志器的非合作目标位姿测量方法,实现目标的精确位姿测量,解决非合作目标位姿测量准确度低的技术问题。
本发明采用如下技术方案,一种基于标志器的非合作目标位姿测量方法,采用二进制平方标记作为标志器,在交会对接逼近阶段向目标航天器抛射多个标志器,基于单目视觉测量原理对不同的标志器进行识别与定位,解算相机坐标系与标志器坐标系之间的相对位姿;同时通过识别非合作目标表面的星箭对接环,利用多幅图像的相机位姿与特征点信息确定标志器在目标器表面的相对位置,从而在对接的近距离阶段通过定位标志器即可解算目标器坐标系与相机坐标系之间的相对位姿,本发明提出了一种基于标志器的非合作目标位姿测量方法,包括以下步骤:
(1)建立相机坐标系、图像坐标系和目标器坐标系;
(2)离线标定,获取相机内参和畸变系数;
(3)对图像进行预处理,得到二值图像;
(4)标志器识别:首先在二值图像中进行轮廓检测,根据约束条件选出候选标志器;随后进行编码提取,对候选标志器的四个顶点进行逆时针排序,通过透视变换获得四边形区域的正面视图,基于最大类间方差阈值法OTSU将四边形区域划分为只包含黑白像素的均匀网格,通过识别四边形区域内部的海明编码确定标志器的序号以及初始顶点的位置;
(5)标志器位姿解算:构建标志器坐标系,利用高效N点透视相机位姿估计算法EPNP解算各标志器坐标系与相机坐标系之间的相对位姿;
(6)椭圆识别:首先进行弧段检测,提取整幅图像的边缘点信息,将边缘点分为梯度大于零和梯度小于零两组集合,即递增组和递减组,然后将边缘点合并为弧段,构造包围盒,去除不满足设定条件的弧段;随后进行弧段选择,将得到的弧段划分到四个象限,基于共圆锥曲线六点特征量CNC准则判定弧段是否属于同一椭圆,并基于象限约束和坐标约束获得有效的三弧段组合;再进行参数计算,基于椭圆上平行弦中点的连线过椭圆中心的几何定理,利用三弧段组合得到穿过椭圆中心的四条直线,取所有交点的代数均值作为椭圆中心;对椭圆参数空间进行降维处理,基于投票原则计算出椭圆的长短半轴和偏转角参数;最后进行后处理,去除三弧段中满足椭圆方程的边缘点占比小于设定值或三弧段长度和与椭圆长短半轴和之比小于设定值的候选椭圆;通过聚类算法将属于同一椭圆的多个检测结果进行合并;选取椭圆识别结果中半径最小的同心椭圆作为最终检测结果,即星箭对接环内环对应的椭圆;
(7)特征点三维坐标恢复:根据椭圆拟合的参数,在图像中构建感兴趣区域ROI,并在ROI区域内部利用累计概率霍夫变换进行直线检测,提取出星箭对接环内部相互垂直的两条直线轮廓;计算两条直线与椭圆边界的交点,将椭圆中心及四个交点作为特征点;利用标志器坐标系与相机坐标系间的位姿参数,基于最小二乘迭代的三角测量算法恢复五个特征点在各标志器坐标系下的三维坐标;
(8)标志器定位:根据三角测量恢复出的特征点三维坐标推算出各特征点在目标器坐标系下的三维坐标;基于最近点迭代ICP算法解算各标志器坐标系与目标器坐标系间的位姿参数;
(9)目标器位姿解算:将相机坐标系到标志器坐标系的变换矩阵和标志器坐标系到目标器坐标系的变换矩阵连乘,得到相机坐标系与目标器坐标系间的位姿参数。
发明达到的有益效果:本发明提出一种基于标志器的非合作目标位姿测量方法,实现目标位姿的精确测量,解决非合作目标位姿测量准确度低的技术问题;通过对不同的标志器进行识别与定位,同时基于单目视觉测量原理确定标志器在目标航天器表面的相对位置,精确估算出相机坐标系与目标器坐标系之间的相对位姿;通过引入预先设计的标志器可以更精确地测量位姿信息,同时克服了交会对接近距离阶段,追踪航天器无法获取完整目标航天器特定部件图像以识别定位的缺点。
附图说明
本发明从下面结合附图对实施例的描述中将变得明显和容易理解,其中:
图1为本发明的基于标志器的非合作目标位姿测量方法的流程图;
图2为本发明实施例中采用的基于二进制平方的标志器示意图;
图3为本发明实施例中对属于第一、三象限的弧段构建包围盒的示意图;
图4为本发明实施例中将弧段划分到四个象限的示意图;
图5为本发明实施例中构造共圆锥曲线六点特征量CNC示意图;
图6为本发明实施例中基于象限约束和坐标约束筛选有效三弧段组示意图;
图7为本发明实施例中椭圆参数方程示意图;
图8为本发明实施例中确定椭圆中心的几何特性示意图;
图9为本发明实施例中利用三弧段组确定椭圆中心示意图;
图10为本发明实施例中椭圆筛选流程图;
图11为本发明实施例中三角测量原理示意图;
图12为本发明实施例中追踪航天器与目标航天器间六自由度位姿测量值与真实值的比较结果图。
具体实施方式
下面详细描述本发明的实施例,所述实施例的示例在附图中示出,其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的,仅用于解释本发明,而不能解释为对本发明的限制。
本发明提出一种基于标志器的非合作目标位姿测量方法。基于标志器的非合作目标位姿测量方法的流程图如图1所示,该方法包括以下步骤:
步骤1:建立相机坐标系、图像坐标系和目标器坐标系,基于相机透视投影模型,以相机光心作为相机坐标系原点,X轴、Y轴分别平行于图像坐标系的u轴、v轴,光轴方向作为Z轴,构建相机坐标系;以星箭对接环中心作为目标器坐标系原点,对接环表面法向量方向作为Z轴,X轴、Y轴分别平行于太阳能帆板的长边和短边,建立目标器坐标系;在本方案中由于相机是固连在追踪航天器上的,因此以相机坐标系替代追踪器坐标系。
步骤2:利用张正友的棋盘格标定法对单目相机进行离线标定,获取相机的内参数,即CCD单目相机在相机坐标系X轴和Y轴的归一化焦距fx和fy、CCD相机的主点像素坐标(u0,v0)、径向畸变系数k1和k2以及切向畸变系数p1和p2。
步骤3:图像预处理,得到二值图像,步骤如下:
步骤3-1:高斯滤波平滑,作用为抑制高频噪声,滤波核满足二维高斯分布:
其中,(x,y)为像素点坐标,σ为高斯核的标准差;
步骤3-2:图像灰度化,求出每个像素点R,G,B三个分量的平均值,并分别赋给该像素点,得到灰度图;
步骤3-3:局部自适应阈值化,根据每个像素的邻域块内的像素值分布来确定该像素位置上的二值化阈值,将灰度图转化为二值图像。
步骤4:图2为采用的基于二进制平方的标志器,标志器识别的步骤如下:
步骤4-1:轮廓检测,基于Suzuki和Abe算法获得包含较多噪声的轮廓集合;
步骤4-2:多边形逼近,对集合中每一条轮廓运用Douglas-Peucker算法,获得近似的多边形轮廓及其顶点信息;
步骤4-3:多边形约束,通过设置约束条件筛选出候选标志器,其中约束条件包括多边形的角点数量是否为四、是否为凸多边形、四边形的边长是否满足设定值、轮廓离图像边界距离是否满足设定值以及四边形集合中四顶点之间距离和是否满足设定值。
步骤4-4:对候选标志器顶点按逆时针排序:对于四个顶点零、顶点一、顶点二和顶点三,根据顶点零、顶点一和顶点零、顶点二构成的向量计算有向面积,如果有向面积为负数,即顶点为顺时针排序,交换顶点一和顶点三的位置,使四边形的四个顶点按逆时针排序;
步骤4-5:计算变换矩阵来去除透视投影,获得四边形区域的正面视图;
步骤4-6:对正面视图进行最大类间方差OTSU阈值化:
步骤4-7:将阈值化后的区域划分为均匀网格,统计方格内非零像素值个数,若方格内非零像素个数超过方格内像素个数的一半,则该方格为白色,否则为黑色;
步骤4-8:按行遍历所有轮廓方格,若轮廓方格中存在白色方格,则舍弃该轮廓所属的候选标志器;
步骤4-9:识别内部编码区域,构造与标志器内部网格大小一致的矩阵,遍历所有网格,将黑色方格对应为数值0,白色方格对应为数值1,依次赋值给矩阵的相应元素,则n×n的网格对应于n×n的0-1矩阵;将矩阵看作由n个n维行向量构成,每一个行向量由数据位和校验位组成,以5×5规格的序号为156的标志器为例,每一个行向量由两个数据位和三个校验位组成,其中数据位00对应于编码10000,数据位01对应于编码10111,数据位10对应于编码01001,数据位11对应于编码01110,通过将该特定标志器的每一行向量与候选标志器的对应行向量做异或运算,统计计算结果中值为1的个数之和作为海明距离;利用平衡二叉树搜索,找出候选标志器与字典(特定标志器集合)中海明距离最小的标志器作为匹配结果,由此可知该候选标志器的序号;
步骤4-10:确定候选标志器的序号后,判断它的旋转状态,将标志器分为初始状态、顺时针旋转90°、顺时针旋转180°、顺时针旋转270°四种状态,分别计算各状态下的标志器与字典中该序号对应的标志器间的海明距离,将海明距离为0的状态作为正确旋转状态;以正确旋转状态下标志器左上角的顶点作为顶点零,按逆时针方向确定顶点一、顶点二和顶点三;
步骤4-11:通过亚像素提取算法对顶点位置进一步细化。
步骤5:标志器位姿解算的步骤如下:
步骤5-1:对于每一个标志器,在正确旋转状态下,以标志器中心作为标志器坐标系原点Om,以顶点零到顶点三的向量方向作为Xm轴方向,顶点一到顶点零的向量方向作为Ym轴方向,Zm轴方向由右手定则确立,构建标志器坐标系Om-XmYmZm;
步骤5-2:标志器的实际尺寸为s×s,确定正确旋转状态下的标志器顶点零到顶点三在标志器坐标系下的空间坐标:(-s/2,s/2,0)、(-s/2,-s/2,0)、(s/2,-s/2,0)、(s/2,s/2,0);
步骤5-3:结合相机内参,利用高效N点透视相机位姿估计(EPNP)算法求解相机坐标系与标志器坐标系间的位姿关系,即旋转矩阵Rcm和平移向量tcm。
步骤6:椭圆识别步骤如下:
步骤6-1:通过Canny边缘检测算子提取图像边缘点,确定每个边缘点的位置坐标(xi,yi),利用Sobel算子计算每个边缘点的梯度τi,得到边缘点信息ei=(xi,yi,τi),其中,i=1,2,...,n,τi=dyi/dxi,n为边缘点数量;
步骤6-2:根据边缘点的梯度方向不同,将边缘点分为两组,即由第二象限弧段组ArcII和第四象限弧段组ArcIV构成的递增组、由第一象限弧段组ArcI和第三象限的弧段组ArcIII构成的递减组:
其中,τi表示第i个像素的梯度,ei表示第i个像素点,ArcI、ArcII、ArcIII和ArcIV分别表示属于第一象限、第二象限、第三象限和第四象限的弧段组,∪表示并集运算。
步骤6-3:对边缘点的八连通区域进行检测,将边缘点合并为弧段;
步骤6-4:对每一条弧段构建包围盒,如图3所示,起始点和末尾点分别为e1和et的弧段,弧段长度为t,顶点(e1(x),e1(y))、(et(x),e1(y))、(et(x),et(y))和(e1(x),et(y))构成包含弧段在内的包围盒,e1(x)、e1(y)分别表示边缘点e1的横坐标和纵坐标,et(x)、et(y)分别表示边缘点et的横坐标和纵坐标,设定最短弧长Thlength,若弧段长度t<Thlength,则舍弃该弧段;
步骤6-5:计算共线三点特征量(CNL)去除直线噪声,根据弧段的起始点e1、中间点ei和末尾点et,利用下式计算CNL值:
其中,|·|表示计算行列式;
该行列式的几何解释为三角形Δe1eiet的面积,使用面积与弧段长度之比来判断e1,ei,et三点是否共线,设t表示弧段长度,Th0为给定阈值,即若CNL/t<Th0,则判定该弧段为直线段,并舍弃该弧段。
步骤6-6:将弧段划分到四个象限,如图4所示,根据弧段上下的像素数量不同,对递增组和递减组的弧段再进行划分:
对于递减组ArcI∪ArcIII的弧段,令δ表示包围盒中,弧段上方像素个数与下方像素个数之差,当弧段上部的像素数多于下部时(δ>0),将弧段划分到ArcIII,否则划分到ArcI;
对于递增组ArcII∪ArcIV的弧段,当弧段上部的像素数小于下部时(δ<0),将弧段划分到ArcII,否则划分到ArcIV;
步骤6-7:利用共圆锥曲线六点特征量(CNC)准则判定弧段是否属于同一椭圆,如图5所示,对于两段圆弧和其中分别是的中点和两个端点,是的中点和两个端点,连接得到两条线段交于点P1,连接得到两条线段交于点P2,连接得到两条线段交于点P3,由于共线三点中的任一点可由其余两点线性表示,可以得到下式:
设置CNC最小阈值为ThCNC,若CNC(P,Q)-1<ThCNC,则两条弧段属于同一椭圆;
步骤6-8:在象限约束以及坐标约束下得到四种有效弧段组合,为了减少无效的弧段组合,通过设置象限约束来选择位于相邻象限的弧段,即弧段属于一、二、四象限,弧段属于二、一、三象限,弧段属于三、二、四象限,弧段属于四、三、一象限;通过设置坐标约束,根据弧段顶点之间的相对位置关系,挑选出可能属于同一椭圆的三弧段组,结果下表1所示。
表1有效三弧段分组
结合CNC判定准则和坐标约束,从每一个有效弧段组合中筛选出属于同一椭圆的三弧段组,以属于一、二、三象限的有效弧段组合为例,三弧段组筛选算法的伪代码如下所示:
步骤6-9:椭圆中心估计
如图7所示,平面中任意的椭圆可以用椭圆中心位置(xc,yc),椭圆长半轴a,短半轴b,以及偏转角θ来表示。其数学表达式可写为:
如图8所示,对于弧段组pab,La,Lb分别是两条弧段的左顶点,Ra,Rb分别是上述两条弧段右顶点,Ma,Mb分别为两条弧段的中点,作nd条平行于LaMb的弦,斜率为r1,作nd条平行于MaRb的弦,斜率为r2,点集分别是两组弦的中点集合,其中近似位于直线l1上,斜率为t1,近似位于直线l2上,斜率为t2,表示点集的中间点,表示点集的中间点;利用改进的Theil-Sen算法得到t1和t2,算法伪代码如下:
伪代码中,GetSlope()为求直线斜率的函数,输入为一组中点集合midpoints[],输出为最佳的拟合斜率,middle表示中点个数的一半,slope为由两点计算的斜率,S[]表示斜率集合,Median()为求中值函数。
直线l1和l2的交点C可由下式计算得到:
分别为点的横坐标和纵坐标,分别为点的横坐标和纵坐标,C.x、C.y分别为交点C的横坐标和纵坐标,如图8所示,利用三弧段组筛选算法得到的有效的三弧段αa,αb,αc可以计算出四条直线,产生至多六个交点,取六个交点的代数均值作为椭圆中心位置;
步骤6-10:计算长短半轴及偏转角
如图9所示,首先将包含椭圆剩余参数即长半轴a,短半轴b以及偏转角θ的参数空间降维到半轴比R=b/a以及偏转角θ上,半轴比R、偏转角θ可通过下式计算:
q1为弧段组(αa,αb)的平行弦斜率,q3为弧段组(αd,αc)的平行弦斜率,q2为弧段组(αa,αb)的平行弦中点连线的斜率,q4为弧段组(αd,αc)的平行弦中点连线的斜率,R+为初始半轴比,K+为偏转角的初始斜率,γ和β为简化式;由上式可知,输入一组参数q1,q2,q3,q4便可以求出相应的R,θ。如图8所示,设r1 ab,为弧段组(αa,αb)的平行弦斜率,r1 dc,为弧段组(αd,αc)的平行弦斜率,直线直线为弧段组(αa,αb)的平行弦中点集合所在的直线,直线直线为弧段组(αd,αc)的平行弦中点集合所在的直线,根据上述的Theil-Sen算法可以得到直线的斜率集合直线的斜率集合直线的斜率集合直线的斜率集合q1,q2,q3,q4的取值见下表:
表2用于计算半轴比R、偏转角θ的q1,q2,q3,q4赋值表
其中,确定q1,q3的取值后,通过从斜率集合和中取不同的值赋给q2,q4,得到不同的q1,q2,q3,q4组合,对每一组合通过上式计算R,θ,得到R,θ的一维累加器,根据投票原则取累加器的峰值作为最终的R,θ。
长半轴a由下式计算:
a=ax/cos(θ)
上式中,ax为长半轴在x轴上的投影,θ为偏转角,(xc,yc)为椭圆中心坐标,(xi,yi)为三条弧段αa,αb,αc上的每一个边缘点的坐标,R为半轴比,K为偏转角θ对应的正切值,x0和y0为简化式;长半轴a在一维累加器中计算,取累加器的峰值作为a。
短半轴b由下式计算:
b=a·R
至此得到椭圆拟合的五个参数,如图9所示。
步骤6-11:椭圆评价,定义两条椭圆评价准则,第一条:三弧段中有多少边缘点满足拟合的椭圆方程,即计算满足椭圆方程的边缘点个数与边缘点总数的比值,该值越大则椭圆评分越高;第二条:三弧段的弧长总和应大于拟合出的椭圆长半轴与短半轴之和,即计算三弧段长度和与椭圆长、短半轴之和的比值,该值越大则椭圆评分越高;最终剔除评分低于设定阈值的候选椭圆;
步骤6-12:比较两个椭圆εi,εj的中心距离、半轴距离和偏转角差值来判断椭圆相似性:
δa=(|εi.a-εj.a|/max(εi.a,εj.a))<0.1
δb=(|εi.b-εj.b|/min(εi.b,εj.b))<0.1
式中,δc表示两椭圆的中心距离,δa表示两椭圆的长半轴距离,δb表示两椭圆间的短半轴距离,δθ表示两椭圆间的偏转角距离,εi.a、εi.b分别表示椭圆εi的长半轴和短半轴,εj.a、εj.b分别表示椭圆εj的长半轴和短半轴,εi.xc、εi.yc分别表示椭圆εi中心的横坐标和纵坐标,εj.xc、εj.yc分别表示椭圆εj中心的横坐标和纵坐标,εi.θ、εj.θ分别表示椭圆εi,εj的偏转角;
当以上条件成立时,椭圆εi,εj归为同一聚类,选择聚类中心作为最终检测出的椭圆,则所有的聚类中心组成椭圆集合;
步骤6-13:目标器表面包含具有同心圆结构的星箭对接环以及尺寸较小的圆形喷嘴,因此选取椭圆集合中同心椭圆中半径小的作为最终检测结果。该步骤流程图如图10所示。
步骤7:特征点三维坐标恢复步骤如下,
步骤7-1:在图像中提取出椭圆感兴趣区域(ROI),即以拟合出的椭圆中心作为矩形边界的中心点,以椭圆的长轴和短轴分别作为矩形边界的长和宽,椭圆的偏转角作为矩形边界的偏转角,生成内切于矩形边界的椭圆边界,并以椭圆中心作为种子点,基于漫水填充算法提取椭圆边界内部的图像区域。
步骤7-2:在ROI区域中,基于累计概率霍夫变换算法进行直线检测,提取出星箭对接环内部相互垂直的两条直线轮廓,并计算直线与椭圆边界的四个交点,与椭圆中心共组成五个椭圆特征点。为了方便不同视角下的特征点匹配,将特征点按固定顺序保存,即椭圆中心、上顶点、下顶点、左顶点和右顶点。
步骤7-3:计算单个标志器在目标器表面的位置,如图11所示,令三维空间点P=[x,y,z,1]T,对应的二维投影点为p=[u,v,1]T,由透视投影成像模型可得:
式中,ρ为非零常数因子,K为相机内参矩阵,R和t分别表示标志器坐标系到相机坐标系的旋转矩阵和平移向量;M=K[R t]称为相机的投影矩阵,因此空间三维点到二维投影点的变换关系可用投影矩阵M来描述。以两幅视图为例,三维点P对应的投影矩阵分别表示为:
M1=K[Rcm1 tcm1]
M2=K[Rcm2 tcm2]
式中,Rcm1,tcm1和Rcm2,tcm2分别为两幅视图下相机坐标系相对于第i个标志器坐标系Omi-XmiYmiZmi的旋转矩阵和平移向量,已通过高效N点透视相机位姿估计(EPNP)算法求解得到。
其中A为P左侧的系数矩阵,通过最小二乘法(LSM)可得到点P在该标志器坐标系下的三维坐标。同理,计算椭圆上的五个特征点在各个标志器下的三维坐标。
步骤8:标志器定位的步骤如下;
步骤8-1:根据三角测量恢复出的特征点在第i个标志器坐标系下的三维坐标,计算各特征点在目标器坐标系下的三维坐标,即:
其中,和分别表示椭圆中心、上顶点、下顶点、左顶点和右顶点在目标器坐标系下的三维坐标,Pkm(k=1,2,3,4)和分别表示上下左右四顶点和椭圆中心在第i个标志器坐标系下的三维坐标,Dis(·)表示计算两个三维点的欧式距离,s表示对接环的半径大小;
步骤8-2:迭代最近点求解位姿,经过上述步骤可以得到星箭对接环上五个特征点在第i个标志器坐标系下的三维坐标Pi m(i=1,2,3,4,5)和目标器标志器下的三维坐标Pi t(i=1,2,3,4,5),基于最近点迭代(ICP)算法计算使得下列目标函数达到最小值时的旋转矩阵R和平移向量t:
式中,J为目标函数,反映累计的重投影误差大小,||·||2表示求取二范数,Rmt和tmt分别表示标志器坐标系到目标器坐标系的旋转矩阵和平移向量,Pkm(k=1,2,3,4,5)表示特征点在标志器坐标系下的三维坐标,Pkt(k=1,2,3,4,5)表示特征点在目标器坐标系下的三维坐标。
步骤9:目标器位姿解算,相机坐标系到标志器坐标系的变换矩阵为Rcm,tcm分别为相机坐标系相对于第i个标志器坐标系Omi-XmiYmiZmi的旋转矩阵和平移向量,标志器坐标系到目标器坐标系的变换矩阵为Rmt和tmt分别表示第i个标志器坐标系到目标器坐标系的旋转矩阵和平移向量,因此相机坐标系到目标器坐标系的变换矩阵为:Rct,tct分别为相机坐标系相对于目标器坐标系的旋转矩阵和平移向量。
X偏移量、Y偏移量和Z偏移量分别为tct的三个分量。
欧拉角中比较常用的一种是用“偏航-俯仰-滚转”(yaw-pitch-roll)3个角度来描述一个旋转,它等价于ZYX轴的旋转,即
a绕物体的Z轴旋转,得到偏航角yaw;
b绕旋转后的Y轴旋转,得到俯仰角pitch;
c绕旋转后的X轴旋转,得到滚转角roll。
旋转矩阵表示为:
上式中,φ为偏航角,θ为俯仰角,ψ为滚转角,Rz(φ)表示绕z轴的旋转矩阵,Ry(θ)表示绕y轴的旋转矩阵,Rx(ψ)表示绕x轴的旋转矩阵,rij(i=1,2,3;j=1,2,3)表示旋转矩阵R中的各分量;
由上式可得目标航天器相对于追踪星的姿态参数,即相机的姿态参数:
ψ=atan2(r32,r33)
φ=atan2(r21,r11)
其中,φ为偏航角,θ为俯仰角,ψ为滚转角,atan2(y,x)为反正切函数,等价于atan(y/x),r11,r21,r31,r32,r33为旋转矩阵R对应下标的分量。
至此,便得到了目标航天器相对追踪星的姿态参数:偏航角φ、俯仰角θ、翻滚角ψ以及X偏移量、Y偏移量、Z偏移量。本发明提出的基于标志器的非合作目标位姿测量方法的测量值与真实值的比较结果见图12。
Claims (10)
1.一种基于标志器的非合作目标位姿测量方法,其特征在于,包括以下步骤:
(1)建立相机坐标系、图像坐标系和目标器坐标系;
(2)离线标定,获取相机内参和畸变系数;
(3)对图像进行预处理,得到二值图像;
(4)标志器识别:首先在二值图像中进行轮廓检测,根据约束条件选出候选标志器;随后进行编码提取,对候选标志器的四个顶点进行逆时针排序,通过透视变换获得四边形区域的正面视图,基于最大类间方差阈值法OTSU将四边形区域划分为只包含黑白像素的均匀网格,通过识别四边形区域内部的海明编码确定标志器的序号以及初始顶点的位置;
(5)标志器位姿解算:构建标志器坐标系,利用高效N点透视相机位姿估计算法EPNP解算各标志器坐标系与相机坐标系之间的相对位姿;
(6)椭圆识别:首先进行弧段检测,提取整幅图像的边缘点信息,将边缘点分为梯度大于零和梯度小于零两组集合,即递增组和递减组,然后将边缘点合并为弧段,构造包围盒,去除不满足设定条件的弧段;随后进行弧段选择,将得到的弧段划分到四个象限,基于共圆锥曲线六点特征量CNC准则判定弧段是否属于同一椭圆,并基于象限约束和坐标约束获得有效的三弧段组合;接着对三弧段组合进行参数计算,基于椭圆上平行弦中点的连线过椭圆中心的几何定理,利用三弧段组合得到穿过椭圆中心的四条直线,取所有交点的代数均值作为椭圆中心;对椭圆参数空间进行降维处理,基于投票原则计算出椭圆的长短半轴和偏转角参数;最后进行后处理,去除三弧段中满足椭圆方程的边缘点占比小于设定值或三弧段长度和与椭圆长短半轴和之比小于设定值的候选椭圆;通过聚类算法将属于同一椭圆的多个检测结果进行合并;选取椭圆识别结果中半径最小的同心椭圆作为最终检测结果,即星箭对接环内环对应的椭圆;
(7)特征点三维坐标恢复:根据椭圆拟合的参数,在图像中构建感兴趣区域ROI,并在ROI区域内部利用累计概率霍夫变换进行直线检测,提取出星箭对接环内部相互垂直的两条直线轮廓;计算两条直线与椭圆边界的交点,将椭圆中心及四个交点作为特征点;利用标志器坐标系与相机坐标系间的位姿参数,基于最小二乘迭代的三角测量算法恢复五个特征点在各标志器坐标系下的三维坐标;
(8)标志器定位:根据三角测量恢复出的特征点三维坐标推算出各特征点在目标器坐标系下的三维坐标;基于最近点迭代ICP算法解算各标志器坐标系与目标器坐标系间的位姿参数;
(9)目标器位姿解算:将相机坐标系到标志器坐标系的变换矩阵和标志器坐标系到目标器坐标系的变换矩阵连乘,得到相机坐标系与目标器坐标系间的位姿参数。
2.根据权利要求1所述的一种基于标志器的非合作目标位姿测量方法,其特征在于,在步骤(1)中基于相机透视投影模型,以相机光心作为相机坐标系原点,X轴、Y轴分别平行于图像坐标系的u轴、v轴,光轴方向作为Z轴,构建相机坐标系;以星箭对接环中心作为目标器坐标系原点,对接环表面法向量方向作为Z轴,X轴、Y轴分别平行于太阳能帆板的长边和短边,建立目标器坐标系。
3.根据权利要求1所述的一种基于标志器的非合作目标位姿测量方法,其特征在于,步骤(2)中利用张正友的棋盘格标定法对单目相机进行离线标定,获取相机的内参数,即CCD单目相机在相机坐标系X轴和Y轴的归一化焦距fx和fy、CCD相机的主点像素坐标(u0,v0)、径向畸变系数k1和k2以及切向畸变系数p1和p2。
5.根据权利要求1所述的一种基于标志器的非合作目标位姿测量方法,其特征在于,步骤(4)中标志器识别的步骤如下:
41)轮廓检测:基于Suzuki和Abe算法获得轮廓集合;
42)多边形逼近:对轮廓集合中每一条轮廓运用Douglas-Peucker算法,获得多边形轮廓及其顶点信息;
43)多边形约束:通过设置约束条件筛选出候选标志器,其中约束条件包括多边形的角点数量是否为四、是否为凸多边形、四边形的边长是否满足设定值、轮廓离图像边界距离是否满足设定值以及四边形集合中四顶点之间距离和是否满足设定值;
44)对候选标志器顶点按逆时针排序:对于四个顶点零、顶点一、顶点二和顶点三,根据顶点零、顶点一和顶点零、顶点二构成的向量计算有向面积,如果有向面积为负数,即顶点为顺时针排序,交换顶点一和顶点三的位置,使四边形的四个顶点按逆时针排序;
45)计算变换矩阵去除透视投影,获得四边形区域的正面视图;
46)对正面视图进行最大类间方差OTSU阈值化:
47)将阈值化后的区域划分为均匀网格,统计每个方格内非零像素值个数,若方格内非零像素个数超过方格内像素个数的一半,则该方格为白色,否则为黑色;
48)按行遍历所有轮廓方格,若轮廓方格中存在白色方格,则舍弃该轮廓所属的候选标志器;
49)识别内部编码区域:构造与标志器内部网格大小一致的矩阵,遍历所有网格,将黑色方格对应为数值0,白色方格对应为数值1,依次赋值给矩阵的相应元素,n×n的网格对应于n×n的0-1矩阵;将矩阵看作由n个n维行向量构成,每一个行向量由数据位和校验位组成,将特定标志器的每一行向量与候选标志器的对应行向量做异或运算,统计计算结果中值为1的个数之和作为海明距离;利用平衡二叉树搜索,找出候选标志器与字典,即特定标志器集合中海明距离最小的标志器作为匹配结果,得到候选标志器的序号;
410)判断候选标志器的旋转状态:将标志器分为初始状态、顺时针旋转90°、顺时针旋转180°和顺时针旋转270°四种状态,分别计算各状态下的标志器与字典中该序号的标志器间的海明距离,将海明距离为0的状态作为正确旋转状态;以正确旋转状态下标志器左上角的顶点作为顶点零,按逆时针方向确定顶点一、顶点二和顶点三;
411)通过亚像素提取算法对顶点位置进一步细化。
6.根据权利要求1所述的基于标志器的非合作目标位姿测量方法,其特征在于,步骤(5)中标志器位姿解算的步骤如下:
51)对于每一个标志器,在正确旋转状态下,以标志器中心作为标志器坐标系原点Om,以顶点零到顶点三的向量方向作为Xm轴方向,顶点一到顶点零的向量方向作为Ym轴方向,Zm轴方向由右手定则确立,构建标志器坐标系Om-XmYmZm;
52)标志器的实际尺寸为s×s,确定正确旋转状态下的标志器顶点零到顶点三在标志器坐标系下的空间坐标:(-s/2,s/2,0)、(-s/2,-s/2,0)、(s/2,-s/2,0)、(s/2,s/2,0);
53)利用高效N点透视相机位姿估计EPNP算法求解相机坐标系到标志器坐标系的旋转矩阵Rcm和平移向量tcm。
7.根据权利要求1所述的基于标志器的非合作目标位姿测量方法,其特征在于,步骤(6)中椭圆识别的步骤如下:
61)通过Canny边缘检测算子提取图像边缘点,确定每个边缘点的位置坐标(xi,yi),利用Sobel算子计算每个边缘点的梯度τi,得到边缘点信息ei=(xi,yi,τi),其中,i=1,2,...,n,τi=dyi/dxi,n为边缘点数量;
62)根据边缘点的梯度方向不同,将边缘点分为两组,即由第二象限弧段组ArcII和第四象限弧段组ArcIV构成的递增组、由第一象限弧段组ArcI和第三象限的弧段组ArcIII构成的递减组:
其中,τi表示第i个边缘点像素的梯度,ei表示第i个边缘点,ArcI、ArcII、ArcIII和ArcIV分别表示属于第一象限、第二象限、第三象限和第四象限的弧段组,∪表示并集运算;
63)对边缘点的八连通区域进行检测,将边缘点合并为弧段;
64)对每一条弧段构建包围盒:起始点和末尾点分别为e1和et的弧段,弧段长度为lt,顶点(e1(x),e1(y))、(et(x),e1(y))、(et(x),et(y))和(e1(x),et(y))构成包含弧段在内的包围盒,e1(x)、e1(y)分别表示边缘点e1的横坐标和纵坐标,et(x)、et(y)分别表示边缘点et的横坐标和纵坐标,设定最短弧长Thlength,若弧段长度lt<Thlength,则舍弃该弧段;
65)基于共线三点特征量CNL准则来去除直线噪声:根据弧段的起始点e1、中间点ei和末尾点et,利用下式计算CNL值:
其中,|·|表示计算行列式;
该行列式的几何解释为三角形Δe1eiet的面积,使用面积与弧段长度之比判断e1、ei、et三点是否共线,lt表示弧段长度,Th0为给定阈值,即若CNL/lt<Th0,则判定该弧段为直线段,舍弃该弧段;
66)将弧段划分到四个象限:根据弧段上下的像素数量不同,对递增组和递减组的弧段进行再划分:
对于递减组ArcI∪ArcIII的弧段,δ表示每一弧段的包围盒中弧段上方像素个数与下方像素个数之差,当弧段上部的像素数多于下部时,即δ>0,将弧段划分到ArcIII,否则划分到ArcI;
对于递增组ArcII∪ArcIV的弧段,当弧段上部的像素数小于下部时,即δ<0,将弧段划分到ArcII,否则划分到ArcIV;
67)利用共圆锥曲线六点特征量CNC准则判定弧段是否属于同一椭圆:对于两段圆弧和其中分别是的中点和两个端点,是的中点和两个端点,连接得到两条直线交于点P1,连接 得到两条直线交于点P2,连接得到两条直线交于点P3,由此得到下式:
设置CNC最小阈值为ThCNC,若CNC(P,Q)-1<ThCNC,则两条弧段属于同一椭圆;
68)在象限约束以及坐标约束下得到三弧段组:设置象限约束选择位于相邻象限的弧段,即有效弧段组合包括:弧段属于一、二、四象限,弧段属于二、一、三象限,弧段属于三、二、四象限,弧段属于四、三、一象限;结合CNC判定准则和弧段端点的相对位置约束,即坐标约束,从每一个有效弧段组合中筛选出属于同一椭圆的三弧段组;
69)确定椭圆中心:对于弧段组pab,La,Lb分别是两条弧段的左顶点,Ra,Rb分别是上述两条弧段右顶点,Ma,Mb分别为两条弧段的中点,作nd条平行于LaMb的平行弦,斜率为r1,作nd条平行于MaRb的平行弦,斜率为r2,点集分别是两组弦的中点集合,其中近似位于直线l1上,斜率为t1,近似位于直线l2上,斜率为t2,表示点集的中间点,表示点集的中间点;利用改进的Theil-Sen算法得到斜率t1和t2;
直线l1和l2的交点C可由下式计算得到:
分别为点的横坐标和纵坐标,分别为点的横坐标和纵坐标,C.x、C.y分别为交点C的横坐标和纵坐标,根据有效的三弧段组的弧段αa,αb,αc可以计算出四条直线,产生至多六个交点,取六个交点的代数均值作为椭圆中心位置;
610)计算长短半轴及偏转角:将包含长半轴a、短半轴b和偏转角θε的参数空间降维到半轴比Rh=b/a以及偏转角θε上,半轴比Rh、偏转角θε可通过下式计算:
上式中,q1为弧段组(αa,αb)的平行弦斜率,q3为弧段组(αd,αc)的平行弦斜率,q2为弧段组(αa,αb)的平行弦中点连线的斜率,q4为弧段组(αd,αc)的平行弦中点连线的斜率,R+为初始半轴比,K+为偏转角的初始斜率,γ和β为简化式;
r1 ab,为弧段组(αa,αb)的平行弦斜率,r1 dc,为弧段组(αd,αc)的平行弦斜率,直线直线为弧段组(αa,αb)的平行弦中点集合所在的直线,直线直线为弧段组(αd,αc)的平行弦中点集合所在的直线,根据Theil-Sen算法可以得到直线的斜率集合直线的斜率集合直线的斜率集合直线的斜率集合确定q1,q3的取值后,通过从斜率集合和中取不同的值赋给q2,q4,得到不同的q1,q2,q3,q4组合,对每一组合通过上式计算半轴比Rh、偏转角θε,得到半轴比Rh、偏转角θε的一维累加器,根据投票原则取累加器的峰值作为最终的半轴比Rh、偏转角θε;
长半轴a可表示为:
a=ax/cos(θε)
上式中,ax为长半轴在x轴上的投影,θε为偏转角,(xc,yc)为椭圆中心坐标,(xi,yi)为三条弧段αa,αb,αc上的每一个边缘点的坐标,tan(θε)为偏转角θε对应的正切值,x0和y0为简化式;长半轴a在一维累加器中计算,取累加器的峰值作为a;
短半轴b由下式计算:
b=a·Rh
得到椭圆拟合的五个参数;
611)椭圆评价:计算三弧段中满足拟合的椭圆方程的边缘点与边缘点总数的比值,比值越大则椭圆评分越高;计算三弧段的弧长总和与拟合的椭圆长半轴与短半轴之和的比值,比值越大则椭圆评分越高;最终剔除评分低于设定阈值的候选椭圆;
612)椭圆聚类:比较两个椭圆εi,εj的中心距离、半轴距离和偏转角差值判断椭圆相似性:
δa=(|εi.a-εj.a|/max(εi.a,εj.a))<0.1
δb=(|εi.b-εj.b|/min(εi.b,εj.b))<0.1
式中,δc表示两椭圆的中心距离,δa表示两椭圆的长半轴距离,δb表示两椭圆间的短半轴距离,表示两椭圆间的偏转角之差,εi.a、εi.b分别表示椭圆εi的长半轴和短半轴,εj.a、εj.b分别表示椭圆εj的长半轴和短半轴,εi.xc、εi.yc分别表示椭圆εi中心的横坐标和纵坐标,εj.xc、εj.yc分别表示椭圆εj中心的横坐标和纵坐标,εi.θε、εj.θε分别表示椭圆εi,εj的偏转角;∠(εi.θε-εj.θε)为两椭圆偏转角之差的弧度制表达;
当以上条件成立时,椭圆εi,εj归为同一聚类,选择聚类中心作为检测出的椭圆,则所有的聚类中心组成椭圆集合;
613)椭圆筛选:选取椭圆集合中同心椭圆中半径小的椭圆作为最终检测结果。
8.根据权利要求1所述的基于标志器的非合作目标位姿测量方法,其特征在于,步骤(7)中特征点三维坐标恢复的步骤如下:
71)在图像中提取出椭圆感兴趣区域ROI,即以拟合的椭圆中心作为矩形边界的中心点,以椭圆的长轴和短轴分别作为矩形边界的长和宽,椭圆的偏转角作为矩形边界的偏转角,生成内切于矩形边界的椭圆边界,并以椭圆中心作为种子点,基于漫水填充算法提取椭圆边界内部的图像区域;
72)在ROI区域中,基于累计概率霍夫变换算法进行直线检测,提取出星箭对接环内部相互垂直的两条直线轮廓,并计算直线与椭圆边界的四个交点,与椭圆中心共组成五个椭圆特征点,将特征点按固定顺序保存,即按椭圆中心、上顶点、下顶点、左顶点和右顶点的顺序保存;
73)计算单个标志器在目标器表面的位置,三维空间点P=[x,y,z,1]T,对应的二维投影点为p=[u,v,1]T,由透视投影成像模型得到:
式中,ρ为非零常数因子,K为相机内参矩阵,R和t分别表示标志器坐标系到相机坐标系的旋转矩阵和平移向量;M=K[R t]表示相机的投影矩阵,在两幅视图中,三维空间点P对应的投影矩阵分别表示为:
M1=K[Rcm1 tcm1]
M2=K[Rcm2 tcm2]
式中,Rcm1,tcm1和Rcm2,tcm2分别为两幅视图下相机坐标系相对于第i个标志器坐标系Omi-XmiYmiZmi的旋转矩阵和平移向量;
其中A为P左侧的系数矩阵,通过最小二乘法LSM得到三维空间点P在第i个标志器坐标系下的三维坐标,计算椭圆上的五个特征点在各标志器坐标系下的三维坐标。
9.根据权利要求1所述的基于标志器的非合作目标位姿测量方法,其特征在于,步骤(8)中标志器定位的步骤如下:
81)根据三角测量恢复出的特征点在第i个标志器坐标系下的三维坐标,计算各特征点在目标器坐标系下的三维坐标,即:
其中,和分别表示椭圆中心、上顶点、下顶点、左顶点和右顶点在目标器坐标系下的三维坐标,和分别表示上下左右四顶点和椭圆中心在第i个标志器坐标系下的三维坐标,Dis(·)表示计算两个三维点的欧式距离,sr表示对接环的半径大小;
82)迭代最近点求解位姿,由星箭对接环上五个特征点在第i个标志器坐标系下的三维坐标和目标器坐标系下的三维坐标,基于最近点迭代ICP算法计算使得下列目标函数达到最小值时的旋转矩阵R和平移向量t:
10.根据权利要求1所述的基于标志器的非合作目标位姿测量方法,其特征在于,步骤(9)中目标器位姿解算,相机坐标系到标志器坐标系的变换矩阵为Rcm,tcm分别为相机坐标系相对于第i个标志器坐标系Omi-XmiYmiZmi的旋转矩阵和平移向量,标志器坐标系到目标器坐标系的变换矩阵为Rmt和tmt分别表示第i个标志器坐标系到目标器坐标系的旋转矩阵和平移向量,因此相机坐标系到目标器坐标系的变换矩阵为:Rct,tct分别为相机坐标系相对于目标器坐标系的旋转矩阵和平移向量;
X偏移量、Y偏移量和Z偏移量分别为tct的三个分量;
旋转矩阵可表示为:
上式中,φ为偏航角,θ为俯仰角,ψ为滚转角,Rz(φ)表示绕z轴的旋转矩阵,Ry(θ)表示绕y轴的旋转矩阵,Rx(ψ)表示绕x轴的旋转矩阵,rij(i=1,2,3;j=1,2,3)表示旋转矩阵R中的各分量;
由上式得到目标航天器相对于追踪航天器,即相机的姿态参数:
ψ=a tan2(r32,r33)
φ=a tan2(r21,r11)
上式中,φ为偏航角,θ为俯仰角,ψ为滚转角,a tan2(y,x)为反正切函数,等价于a tan(y/x),r11,r21,r31,r32,r33为旋转矩阵R对应下标的分量;
追踪航天器相对目标航天器的姿态参数:偏航角φ、俯仰角θ、翻滚角ψ以及X偏移量、Y偏移量、Z偏移量。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810359727.7A CN108562274B (zh) | 2018-04-20 | 2018-04-20 | 一种基于标志器的非合作目标位姿测量方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810359727.7A CN108562274B (zh) | 2018-04-20 | 2018-04-20 | 一种基于标志器的非合作目标位姿测量方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108562274A CN108562274A (zh) | 2018-09-21 |
CN108562274B true CN108562274B (zh) | 2020-10-27 |
Family
ID=63535877
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810359727.7A Active CN108562274B (zh) | 2018-04-20 | 2018-04-20 | 一种基于标志器的非合作目标位姿测量方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108562274B (zh) |
Families Citing this family (38)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109269544B (zh) * | 2018-09-27 | 2021-01-29 | 中国人民解放军国防科技大学 | 中低速磁浮车辆悬浮传感器巡检*** |
CN109102567B (zh) * | 2018-10-11 | 2023-02-24 | 北京理工大学 | 一种基于重建误差最小化的位姿参数高精度求解方法 |
CN109521404B (zh) * | 2018-10-12 | 2020-09-04 | 上海交通大学 | 基于fmcw雷达的振动测量的准确度评估方法和*** |
CN109613923A (zh) * | 2018-11-06 | 2019-04-12 | 武汉华中天经通视科技有限公司 | 一种无人直升机的着舰控制方法 |
US11361466B2 (en) * | 2018-11-30 | 2022-06-14 | Casio Computer Co., Ltd. | Position information acquisition device, position information acquisition method, recording medium, and position information acquisition system |
CN109949361A (zh) * | 2018-12-16 | 2019-06-28 | 内蒙古工业大学 | 一种基于单目视觉定位的旋翼无人机姿态估计方法 |
CN109827578B (zh) * | 2019-02-25 | 2019-11-22 | 中国人民解放军军事科学院国防科技创新研究院 | 基于轮廓相似性的卫星相对姿态估计方法 |
CN110197509B (zh) * | 2019-04-30 | 2023-07-11 | 上海理工大学 | 一种基于彩色人工标识的相机位姿求解法 |
CN110472651B (zh) * | 2019-06-17 | 2022-11-29 | 青岛星科瑞升信息科技有限公司 | 一种基于边缘点局部特征值的目标匹配与定位方法 |
CN110189375B (zh) * | 2019-06-26 | 2022-08-23 | 中国科学院光电技术研究所 | 一种基于单目视觉测量的图像目标识别方法 |
CN110390696A (zh) * | 2019-07-03 | 2019-10-29 | 浙江大学 | 一种基于图像超分辨率重建的圆孔位姿视觉检测方法 |
CN110349207B (zh) * | 2019-07-10 | 2022-08-05 | 国网四川省电力公司电力科学研究院 | 一种复杂环境下的视觉定位方法 |
CN110608739B (zh) * | 2019-08-21 | 2021-07-27 | 深圳市人工智能与机器人研究院 | 干扰环境下运动目标的定位方法、***及电子装置 |
CN110647156B (zh) * | 2019-09-17 | 2021-05-11 | 中国科学院自动化研究所 | 基于目标物对接环的对接设备位姿调整方法、*** |
CN111091121B (zh) * | 2019-11-22 | 2022-08-26 | 重庆大学 | 一种基于图像处理的椭圆表盘检测矫正的方法 |
CN110879048A (zh) * | 2019-12-10 | 2020-03-13 | 南昌航空大学 | 一种基于标记点检测的桨叶扭转角实时监测方法 |
CN111413995B (zh) * | 2020-03-24 | 2021-03-05 | 北京科技大学 | 双刚体特征点间相对位置跟踪与姿态同步控制方法及*** |
CN111445533B (zh) * | 2020-03-27 | 2023-08-01 | 广东博智林机器人有限公司 | 一种双目相机标定方法、装置、设备及介质 |
CN111680685B (zh) * | 2020-04-14 | 2023-06-06 | 上海高仙自动化科技发展有限公司 | 一种基于图像的定位方法、装置、电子设备及存储介质 |
CN111536981B (zh) * | 2020-04-23 | 2023-09-12 | 中国科学院上海技术物理研究所 | 一种嵌入式的双目非合作目标相对位姿测量方法 |
WO2022036478A1 (zh) * | 2020-08-17 | 2022-02-24 | 江苏瑞科科技有限公司 | 一种基于机器视觉的增强现实盲区装配引导方法 |
CN112233176A (zh) * | 2020-09-27 | 2021-01-15 | 南京理工大学 | 一种基于标定物的目标位姿测量方法 |
CN112378383B (zh) * | 2020-10-22 | 2021-10-19 | 北京航空航天大学 | 基于圆和线特征非合作目标相对位姿双目视觉测量方法 |
CN113066050B (zh) * | 2021-03-10 | 2022-10-21 | 天津理工大学 | 一种基于视觉的空投货台航向姿态解算方法 |
CN113504543B (zh) * | 2021-06-16 | 2022-11-01 | 国网山西省电力公司电力科学研究院 | 无人机LiDAR***定位定姿***及方法 |
TWI787885B (zh) * | 2021-06-25 | 2022-12-21 | 凌華科技股份有限公司 | 彈出視窗按鈕的非侵入式偵測方法及裝置 |
CN113706621B (zh) * | 2021-10-29 | 2022-02-22 | 上海景吾智能科技有限公司 | 基于带标记图像的标志点定位及姿态获取方法和*** |
CN114596355B (zh) * | 2022-03-16 | 2024-03-08 | 哈尔滨工业大学 | 一种基于合作目标的高精度位姿测量方法及*** |
CN114715447A (zh) * | 2022-04-19 | 2022-07-08 | 北京航空航天大学 | 一种细胞航天器模块对接装置及视觉对准方法 |
CN114963981B (zh) * | 2022-05-16 | 2023-08-15 | 南京航空航天大学 | 一种基于单目视觉的筒状零件对接非接触式测量方法 |
CN114926526B (zh) * | 2022-05-23 | 2023-05-05 | 南京航空航天大学 | 一种基于变焦相机的位姿测量方法 |
CN115330272B (zh) * | 2022-10-13 | 2023-01-31 | 北京理工大学 | 一种复杂作战区域条件下多飞行器目标协同攻击方法 |
CN115597569B (zh) * | 2022-10-31 | 2024-05-14 | 上海勃发空间信息技术有限公司 | 利用断面扫描仪测定桩与船相对位置关系的方法 |
CN115609591B (zh) * | 2022-11-17 | 2023-04-28 | 上海仙工智能科技有限公司 | 一种基于2D Marker的视觉定位方法和***、复合机器人 |
CN116105694B (zh) * | 2022-12-09 | 2024-03-12 | 中国科学院上海技术物理研究所 | 一种多手段光学载荷复合的空间目标三维视觉测量方法 |
CN116051629B (zh) * | 2023-02-22 | 2023-11-07 | 常熟理工学院 | 面向自主导航机器人的高精度视觉定位方法 |
CN117036489B (zh) * | 2023-10-10 | 2024-02-09 | 泉州装备制造研究所 | 基于人工标识和四目全景相机的机器人定位方法及设备 |
CN117115242B (zh) * | 2023-10-17 | 2024-01-23 | 湖南视比特机器人有限公司 | 标志点的识别方法、计算机存储介质和终端设备 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104677340A (zh) * | 2013-11-30 | 2015-06-03 | 中国科学院沈阳自动化研究所 | 一种基于点特征的单目视觉位姿测量方法 |
CN105509733A (zh) * | 2015-11-30 | 2016-04-20 | 上海宇航***工程研究所 | 非合作空间圆目标的相对位姿测量方法 |
CN105806315A (zh) * | 2014-12-31 | 2016-07-27 | 上海新跃仪表厂 | 基于主动编码信息的非合作目标相对测量***及测量方法 |
CN107449402A (zh) * | 2017-07-31 | 2017-12-08 | 清华大学深圳研究生院 | 一种非合作目标的相对位姿的测量方法 |
-
2018
- 2018-04-20 CN CN201810359727.7A patent/CN108562274B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104677340A (zh) * | 2013-11-30 | 2015-06-03 | 中国科学院沈阳自动化研究所 | 一种基于点特征的单目视觉位姿测量方法 |
CN105806315A (zh) * | 2014-12-31 | 2016-07-27 | 上海新跃仪表厂 | 基于主动编码信息的非合作目标相对测量***及测量方法 |
CN105509733A (zh) * | 2015-11-30 | 2016-04-20 | 上海宇航***工程研究所 | 非合作空间圆目标的相对位姿测量方法 |
CN107449402A (zh) * | 2017-07-31 | 2017-12-08 | 清华大学深圳研究生院 | 一种非合作目标的相对位姿的测量方法 |
Non-Patent Citations (1)
Title |
---|
基于视觉技术的非合作航天器相对位姿测量方法;孙增玉 等;《宇航计测技术》;20170830;第37卷(第4期);第2-5页 * |
Also Published As
Publication number | Publication date |
---|---|
CN108562274A (zh) | 2018-09-21 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108562274B (zh) | 一种基于标志器的非合作目标位姿测量方法 | |
CN108122256B (zh) | 一种逼近状态下旋转目标位姿测量的方法 | |
CN107063228B (zh) | 基于双目视觉的目标姿态解算方法 | |
Zhang et al. | Vision-based pose estimation for textureless space objects by contour points matching | |
Kang et al. | Automatic targetless camera–lidar calibration by aligning edge with gaussian mixture model | |
JP6216508B2 (ja) | 3dシーンにおける3d物体の認識および姿勢決定のための方法 | |
JP5705147B2 (ja) | 記述子を用いて3dオブジェクトまたはオブジェクトを表す方法 | |
CN111862201B (zh) | 一种基于深度学习的空间非合作目标相对位姿估计方法 | |
JP2004516533A (ja) | 合成開口レーダーおよび前向き赤外線イメージ重ね合わせ方法 | |
CN108225319B (zh) | 基于目标特征的单目视觉快速相对位姿估计***及方法 | |
CN103136525B (zh) | 一种利用广义Hough变换的异型扩展目标高精度定位方法 | |
CN104040590A (zh) | 用于估计物体的姿态的方法 | |
CN111524168A (zh) | 点云数据的配准方法、***、装置及计算机存储介质 | |
CN106127258A (zh) | 一种目标匹配方法 | |
JP4859061B2 (ja) | 画像の補正方法、補正プログラムおよび画像歪み補正装置 | |
CN113295171B (zh) | 一种基于单目视觉的旋转刚体航天器姿态估计方法 | |
CN114358166B (zh) | 一种基于自适应k均值聚类的多目标定位方法 | |
CN104484647B (zh) | 一种高分辨率遥感图像云高度检测方法 | |
CN115112098B (zh) | 一种单目视觉的一维二维测量方法 | |
CN115760984A (zh) | 一种立方星基于单目视觉的非合作目标位姿测量方法 | |
CN115131433A (zh) | 一种非合作目标位姿的处理方法、装置及电子设备 | |
CN114972451A (zh) | 一种基于旋转不变的SuperGlue匹配的遥感影像配准方法 | |
WO2018057082A1 (en) | Curvature-based face detector | |
Martinez | USV Attitude Position Estimation by a Hovering UAV Using Monocular Images of Deck-Mounted Lights | |
Li et al. | New environmental line feature-based vision navigation: Design and analysis |
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 | ||
CB02 | Change of applicant information |
Address after: 210023, 66 new model street, Gulou District, Jiangsu, Nanjing Applicant after: NANJING University OF POSTS AND TELECOMMUNICATIONS Address before: 210023 Jiangsu city of Nanjing province Ya Dong new Yuen Road No. 9 Applicant before: NANJING University OF POSTS AND TELECOMMUNICATIONS |
|
CB02 | Change of applicant information | ||
GR01 | Patent grant | ||
GR01 | Patent grant |