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。

4.根据权利要求1所述的一种基于标志器的非合作目标位姿测量方法，其特征在于，步骤(3)中图像预处理的步骤如下：

31)高斯滤波平滑，滤波核满足二维高斯分布：

其中，(x,y)为像素坐标，σ为高斯核的标准差；

32)图像灰度化，求出每个像素点的R分量、G分量和B分量的平均值赋给该像素点，得到灰度图；

33)局部自适应阈值化，根据每个像素的邻域块内的像素值分布确定该像素位置上的二值化阈值，将灰度图转化为二值图像。

5.根据权利要求1所述的一种基于标志器的非合作目标位姿测量方法，其特征在于，步骤(4)中标志器识别的步骤如下：

41)轮廓检测：基于Suzuki和Abe算法获得轮廓集合；

42)多边形逼近：对轮廓集合中每一条轮廓运用Douglas-Peucker算法，获得多边形轮廓及其顶点信息；

43)多边形约束：通过设置约束条件筛选出候选标志器，其中约束条件包括多边形的角点数量是否为四、是否为凸多边形、四边形的边长是否满足设定值、轮廓离图像边界距离是否满足设定值以及四边形集合中四顶点之间距离和是否满足设定值；

44)对候选标志器顶点按逆时针排序：对于四个顶点零、顶点一、顶点二和顶点三，根据顶点零、顶点一和顶点零、顶点二构成的向量计算有向面积，如果有向面积为负数，即顶点为顺时针排序，交换顶点一和顶点三的位置，使四边形的四个顶点按逆时针排序；

45)计算变换矩阵去除透视投影，获得四边形区域的正面视图；

46)对正面视图进行最大类间方差OTSU阈值化：

*

其中，[0,L-1]为图像的灰度级范围，th为灰度阈值，th为最佳灰度阈值， 为不同灰度级的类间方差，Arg max(·)表示使目标函数取最大值时的变量值；

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，由此得到下式：其中,Pi为直线交点的像素坐标, 为弧段上点的像素坐标,为对应系数；

通过上式求出系数 并代入下式计算共圆锥曲线六点特征量CNC值：其中,CNC(P,Q)表示两段圆弧的CNC值， 为对应系数，i表示直线交点P的索引，j表示弧段上构成直线的像素点的索引；Π(·)表示累乘运算；

设置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+为偏转角的初始斜率，γ和β为简化式；

r1ab, 为弧段组(αa,αb)的平行弦斜率，r1dc, 为弧段组(α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的旋转矩阵和平移向量；

三维空间点P＝[x,y,z,1]T在两幅图像上的投影点分别为p1＝[u1,v1,1]T和p2＝[u2,v2,

1]T，由p1＝M1P,p2＝M2P， 得到：

其中A为P左侧的系数矩阵，通过最小二乘法LSM得到三维空间点P在第i个标志器坐标系下的三维坐标，计算椭圆上的五个特征点在各标志器坐标系下的三维坐标。

9.根据权利要求1所述的基于标志器的非合作目标位姿测量方法，其特征在于，步骤(8)中标志器定位的步骤如下：

81)根据三角测量恢复出的特征点在第i个标志器坐标系下的三维坐标，计算各特征点在目标器坐标系下的三维坐标，即：

其中， 和 分别表示椭圆中心、上顶点、下顶点、左顶点和

右顶点在目标器坐标系下的三维坐标， 和 分别表示上下左右四顶点和椭圆中心在第i个标志器坐标系下的三维坐标，Dis(·)表示计算两个三维点的欧式距离，sr表示对接环的半径大小；

82)迭代最近点求解位姿，由星箭对接环上五个特征点在第i个标志器坐标系下的三维坐标和目标器坐标系下的三维坐标，基于最近点迭代ICP算法计算使得下列目标函数达到最小值时的旋转矩阵R和平移向量t：式中，J为目标函数，反映累计的重投影误差大小，||·||2表示求取二范数，Rmt和tmt分别表示标志器坐标系到目标器坐标系的旋转矩阵和平移向量， 表示特征点在标志器坐标系下的三维坐标， 表示特征点在目标器坐标系下的三维坐标。

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偏移量。