本发明公开了一种基于序列图像的跟踪精度反演方法，其步骤如下：一、基于现已知的观测卫星的轨/姿数据、目标卫星的轨/姿数据、连续的图像序列，采用坐标变换方法求得选取目标点的理想位置；二、基于同面特征点搜索获取基础变换矩阵方法获取目标点的实际位置；三、求得目标点的理想位置和实际位置之差，结合给定相机参数，可得到以角度形式表示的反演误差。本发明使用所拍摄连续的图像序列，反演追踪观测卫星对目标卫星的指向误差，该方法能够涵盖悬停、接近和绕飞三种情况。
1.一种基于序列图像的跟踪精度反演方法，其特征在于所述方法步骤如下：一、基于现已知的观测卫星的轨/姿数据、目标卫星的轨/姿数据、连续的图像序列，采用坐标变换方法求得选取目标点的理想位置；
二、基于同面特征点搜索获取基础变换矩阵方法获取目标点的实际位置；
三、求得目标点的理想位置和实际位置之差，结合给定相机参数，可得到以角度形式表示的反演误差；
所述光学中心衍生中心包括目标轮廓中心、目标边框中心、有效配准中心；
所述采用坐标变换方法求得选取目标点的理想位置的步骤如下：(1)建立坐标系：为了分析地球、观测卫星、目标卫星三者之间的位置随时间变化的关系，建立J2000.0惯性坐标系、轨道坐标系与卫星本体坐标系；
(2)坐标系变换：若知道观测卫星和目标卫星的轨道根数，由开普勒定律可得观测卫星和目标卫星在J2000.0惯性坐标系中的矢量坐标：p＝Rz(-Ω)Rx(-i)Rz(-ω)p0，式中：i—轨道倾角；Ω—升交点赤经；ω—近地点辐角；Rx—绕x轴旋转变换矩阵；Rz—绕z轴旋转变换矩阵；p0—近焦坐标系下的矢量坐标；
将各量在J2000.0惯性坐标系中的矢量坐标转换到轨道坐标系，则为：qo＝L001Rz(ω+v)Rx(i)Rz(Ω)p，式中：qo—各量在轨道坐标系中的矢量坐标；L001——坐标轴的反向变换矩阵；
为了便于目标特性分析，将各量在轨道坐标系中的矢量坐标转换到本体坐标系，则为：q＝Aqo，
对于三轴稳定卫星，常采用zxy转序，则姿态矩阵由旋转变换矩阵表示为：式中：Ψ—卫星绕偏航轴转动的角度，即偏航角；θ—卫星绕俯仰轴转动的角度，即俯仰角；—卫星绕滚动轴转动的角度，即滚动角；Ψ，θ与 称为卫星姿态描述的欧拉角；
焦面上任意一点p坐标为(yc,zc)，则相机焦平面转换到图像坐标系下的坐标为：[(yc-y0)/D，(zc-z0)/D]；
在本体坐标系下最终结果可表示成q＝[Xc,YC,ZC]的形式，将其中的YC,ZC,带入[(yc-y0)/D，(zc-z0)/D]中，其中D为像元尺寸，即可得到在图像中的坐标形式；
因此，给定观测目标、目标卫星的轨道参数、连续图像的时间间隔、像元尺寸及目标卫星姿态特性的欧拉角后，即可算得目标点的理想位置；
所述J2000.0惯性坐标系FI：OXYZ，原点O位于地球的质心，OX轴和OZ轴分别指向春分点和北极，并与OY轴构成右手直角坐标系；
所述轨道坐标系FO：SXoYoZo，原点S位于卫星的质心，SXo轴与轨道速度方向一致，SZo轴指向地心，SYo轴垂直轨道平面并构成右手直角坐标系；
所述本体坐标系FB：SXbYbZb，原点S位于卫星的质心，三轴为卫星的三个惯性主轴，其中SXb为滚动轴，SYb为俯仰轴，SZb为偏航轴；
所述基于同面特征点搜索获取基础变换矩阵方法获取目标点的实际位置的方法如下：第一步，采用SIFT算法来获取图像的特征点；
第二步，对经过SIFT算法筛选的特征点进行初步筛选，筛选过程如下：对于各特征点Pi，获取其在前帧图上的坐标位置Pi(x,y)before以及本帧图像上的坐标位置Pi(x,y)now，计算其位置移动的变化量：剔除移动距离超过设定阀值Dk的特征点，阀值Dk设定可取：式中，n为配准特征点的总数，τ为阀值系数；
剔除移动距离较大点后，剩余特征点进行如下变换：剩余特征点与目标点的距离 PT(x,y)为目标点在本帧图像中的位置坐标，则有：取 满足 的特征点，其中m为剩余特征点的总数，η为阀值系数，经过与标记点的距离筛选后，获得较近的配准特征点；
对筛选后的配准特征点，取离标记点最近的k个点，以周边三个像素为半径，取各配准特征点周围的小图像，并计算特征点周围对应小图的直方图，对各特征点的直方图，与其他的配准点所对应的直方图进行相关系数计算：其中：H′k(i)＝Hk(i)-(1/N)(∑jHk(j))，N等于直方图中数据的个数；
随后进行求和，并取最大值对应的配准特征点为基准点Po(x,y)：然后以基准点Po(x,y)对其他各配准特征点对应的直方图求相关系数，取相关系数大于阀值ξ＞0.9的认为是同面配准特征点；经过直方图系数计算获取的基准点，认定为同面配准特征点；
当确定了共面配准特征点后，通过RANSAC方法实现两帧图像间的基础变换矩阵，再将前帧的目标位置通过基础变换矩阵变换到本帧上，获取目标点在本帧图像中的坐标；
所述误差用角度表示的形式为： 其中f为相机焦距，a为理想点距离光轴距离，b为实际点距离光轴距离。
一种基于序列图像的跟踪精度反演方法
技术领域
[0001] 本发明属于航天卫星在轨实时跟踪技术领域，涉及一种基于序列图像的跟踪精度反演方法。
背景技术
[0002] 航天领域中，对于卫星在轨实时跟踪领域，国内研究集中在跟踪系统的基本原理、结构、关键技术及跟踪过程，而对于跟踪精度方向也主要集中在误差产生原因、误差在成像链路各个环节的定量分析与计算及误差补偿系统的研究，而对于基于图像进行反演得到跟踪精度的定量计算研究较少，且这是近年来刚刚提出的方向，可以查阅的资料较少。
发明内容
[0003] 本发明的目的在于提供一种基于序列图像的跟踪精度反演方法，使用所拍摄连续的图像序列，反演追踪观测卫星对目标卫星的指向误差，该方法能够涵盖悬停、接近和绕飞三种情况。
[0004] 本发明的目的是通过以下技术方案实现的：
[0005] 一、基于现已知的观测卫星的轨/姿数据、目标卫星的轨/姿数据、连续的图像序列，采用坐标变换方法求得选取目标点的理想位置；
[0006] 二、基于同面特征点搜索获取基础变换矩阵方法获取目标点的实际位置；
[0007] 三、求得目标点的理想位置和实际位置之差，结合给定相机参数，可得到以角度形式表示的反演误差。
[0008] 本发明可以较快速的反演出目标点与实际点的距离之差，并且误差精度在可接受的范围之内，为后续的整个卫星链路的反演提供数据支撑，并为将来发射卫星做一些先验经验。
附图说明
[0009] 图1为坐标系；
[0010] 图2为焦平面像元在相机本体系下的坐标；
[0011] 图3为单幅图像的处理过程；
[0012] 图4为基于同面配准特征点获取基础变换矩阵流程；
[0013] 图5为SIFT算法获取特征点过程；
[0014] 图6为像平面上的角度表示方法。
具体实施方式
[0015] 下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。
[0016] 本发明提供了一种基于序列图像的跟踪精度反演方法，具体包括以下内容：
[0017] 1、问题的提出和假设条件
[0018] 本问题的提出在于：使用所拍摄连续的图像序列，反演追踪观测卫星对目标卫星的指向误差，能够涵盖悬停、接近和绕飞三种情况。
[0019] 本问题的前提包括：
[0020] 1)目标卫星为非合作目标，其形貌和内在参数均未知。
[0021] 2)目标卫星在观测期间，其星载部件可能发生改变。
[0022] 3)可获取的参数主要为拍摄的连续图像序列，并有观测卫星的轨道和姿态数据及目标卫星的轨道和姿态数据。
[0023] 2、指向误差反演研究思路
[0024] 2.1 指向误差反演分析
[0025] 悬停工况：悬停状态包括无需主动控制的同轨追踪悬停以及需要主动控制的正向悬停和侧向悬停。卫星悬停状态下，当相对指向不发生改变时，如果目标卫星维持姿态不变，则在观测卫星的视场中，目标卫星的图像基本维持不变，仅受到定轨误差的影响，发生小幅的漂移。当目标卫星的姿态改变或所载部件运动时，在观测卫星的视场中，目标卫星的成像形状随姿态的变化或部件运动发生改变。
[0026] 接近工况：观测卫星向目标卫星接近的情况下，即使观测卫星视线始终对准目标卫星，由于观测卫星轨道与目标卫星的轨道角均变化，在观测卫星的视场中，目标由小及大，在图中的位置和角度也随之变化。当目标卫星的姿态改变或所载部件运动时，在观测卫星的视场中，目标卫星的成像形状随姿态的变化或部件运动发生改变。
[0027] 绕飞工况：观测卫星与目标卫星处于相近轨道，典型的周期型绕飞有椭圆绕飞、摆型绕飞、空间圆绕飞等。在绕飞工况下，目标卫星在视场内成像形状随绕飞的状态发生改变，当目标卫星的姿态改变或所载部件运动时，将进一步增加成像形状变化的复杂性。
[0028] 对象无论是悬停、接近或者绕飞工况下，目标的成像形状均可能发生变化，还存在成像大小变化的耦合情况。由上述分析可知，在进行指向误差反演时，需考虑的是通用性算法，同时涵盖目标在成像面上的缩放、三维旋转效应。
[0029] 2.2 指向误差反演思路
[0030] 基于现已知的观测卫星的轨/姿数据、目标卫星的轨/姿数据、连续的图像序列，拟采用坐标变换方法求得选取目标点的理想位置，再通过图像匹配的方法找到目标点的实际位置，求得两者之差(ΔX，ΔY)，结合给定相机参数，转化成为角度误差，当图像序列的每一幅图均获取了该角度误差后，则会形成指向误差的变化序列。
[0031] 该问题进一步分解为三个问题：
[0032] 1)目标点如何选取？
[0033] 2)如何确定图像中目标点的理想位置？
[0034] 3)如何确定图像中目标点的实际位置？
[0035] 对于目标点的确定问题，与观测卫星的指向对准策略相关。理论上来说，有几种预期目标点：
[0036] 1)质心点：目标卫星的质心位置作为预期对准点具有最佳的唯一性和不变性，但是由于目标形貌、内部参数等信息未知，因此质心位置仅理论存在，但无法获知。
[0037] 2)几何中心：几何中心对于无大部件移动的对象而言，也具有较好的唯一性和不变性，但建立和获取几何中心至少需要在前期对对象的外形形状和尺寸开展深入的分析，以确定不同角度下的几何中心位置，显然本问题中无达成该条件。
[0038] 3)某个典型形貌：基于典型形貌的对准广泛用于光学指向任务中，但其有两个前提，第一，典型形貌需要提前预知；第二典型形貌需要在始终存在于视场中；显然对于本问题，该方法所需的前提均不成立，既没有对象的外形，更不用提典型形貌，而且由于姿态改变或绕飞使得观测角度的变化，也不能保证典型形貌始终存在于视场中。
[0039] 4)光学中心：从工程实现角度来看，光学中心的可实现性最高，然而光学中心随成像的形状改变而变化，不存在唯一性，不变性也较差，这对反演指向误差是不利的。对于红外成像，由于目标的细节难以准确反映，除了光学中心外，也没有更好的方法。
[0040] 对于目标点的确定，本发明给出的方法是在光学中心及其衍生中心为基础下，发展由人介入进行标记的方法。所述的光学衍生中心，包括如下几种：
[0041] 1)目标轮廓中心：即以目标边缘轮廓为像集的光学中心。
[0042] 2)目标边框中心：即以包括目标边缘的矩形框的几何中心。
[0043] 3)有效配准中心：即前帧图像与当帧图像的有效配准点为像集的光学中心。
[0044] 上述光学中心及其衍生中心均可由计算机自动处理获得，但由于这类中心不存在物理意义，均无不变性，由此引入人介入进行标记。人的介入是必要也是有可实现性的，首先，由于误差反演是获取图像序列数据后的后处理，因此具备人处理的条件，人可以通过小部件的认知，主体的认知对质心点、几何中心点或者是表面的特征部件进行先期学习和基本标记。
[0045] 对于目标点在各幅图像中的位置确定问题，人难以做到的是对每帧图像均进行标记和分析，因此本发明的方法是人首先对序列图像进行按时序分组，对于每一组的首帧图像，选取目标点，该目标点可以任意选取，而该组的其他帧图像，以第一帧为基准，通过坐标变换确定第一帧的目标点在后续各帧的位置，再找到目标点的实际位置，以实现序列图的误差反演。
[0046] 3、指向误差反演算法
[0047] 3.1 基于坐标变换方法确定目标点的理论位置
[0048] 3.1.1 轨道要素：卫星在空间轨道上的位置，可用6个参数描述。这些参数分别为：
[0049] (1)升交点赤经Ω：轨道升交点与地球惯性坐标系指向春分点轴的夹角。
[0050] (2)轨道倾角i：轨道面与地球赤道面的夹角。
[0051] (3)近地点幅角ω：轨道面内由升交点到近地点拱线的夹角。
[0052] (4)轨道长半轴a：椭圆长轴的一半。
[0053] (5)轨道偏心率e：椭圆两焦点之间的距离与长轴的比值，该值决定了椭圆偏心大小。
[0054] (6)卫星过近地点时刻τ：卫星飞过近地点的时刻，一般用年、月、日、时、分、秒表示，它是描述卫星运动时间的起始点。
[0055] 在六个要素中，升交点赤经Ω和轨道倾角i决定了轨道面在空间的位置；近地点幅角ω决定了轨道面的指向；轨道长半轴a和轨道偏心率e决定了轨道大小和形状；卫星飞过近地点的时刻τ决定了卫星在轨道上的位置。
[0056] 3.1.2 空间位置确定：为了确定目标表面(网格单元)光线的入射方向矢量及探测器的观测方向矢量，本发明建立参考坐标系，基于目标卫星、观测卫星轨道参数及目标卫星的姿态特性进行坐标变换确定它们之间的相互位置关系。
[0057] (1)建立坐标系：为了分析地球、观测卫星、目标卫星三者之间的位置随时间变化的关系，首先建立三个坐标系：J2000.0惯性坐标系、轨道坐标系与卫星本体坐标系，如图1所示。
[0058] J2000.0惯性坐标系FI：OXYZ。原点O位于地球的质心，OX轴和OZ轴分别指向春分点和北极，并与OY轴构成右手直角坐标系。
[0059] 轨道坐标系FO：SXoYoZo。原点S位于卫星的质心，SXo轴与轨道速度方向一致，SZo轴指向地心，SYo轴垂直轨道平面并构成右手直角坐标系。
[0060] 本体坐标系FB：SXbYbZb。原点S位于卫星的质心，三轴为卫星的三个惯性主轴。其中SXb为滚动轴，SYb为俯仰轴，SZb为偏航轴。
[0061] (2)坐标系变换：观测卫星和目标卫星相对于地球的运动遵循开普勒定律，常用椭圆运动方程表示为：
[0062]
[0063] (1)
[0064] v＝M+e(2-e2/4+5e4/96)sinM+e2(5/4-11e2/24)sin2M+ (2)
[0065] e3(13/12-43e2/64)sin3M+103e4sin4M/96+1097e5sin5M/960
[0066]
[0067] (3)
[0068] 式中：a——轨道半长轴(km)；
[0069] e——轨道偏心率；
[0070] t——观测时刻；
[0071] τ——过近地点时刻；
[0072] μ——开普勒常数，μ＝3.986032×1014m3·s-2。
[0073] t与τ均采用J2000.0标准历元起算的儒略日计时。
[0074] 若知道观测卫星和目标卫星的轨道根数，由开普勒定律可得观测卫星和目标卫星在J2000.0惯性坐标系中的矢量坐标：
[0075] p＝Rz(-Ω)Rx(-i)Rz(-ω)p0 (4)式中：i——轨道倾角(deg)；
[0076] Ω——升交点赤经(deg)；
[0077] ω——近地点辐角(deg)。
[0078] Rx——绕x轴旋转变换矩阵，
[0079] Rz——绕z轴旋转变换矩阵，
[0080] p0——近焦坐标系下的矢量坐标，p0＝(rcosv,rsinv,0)。
[0081] 将各量在J2000.0惯性坐标系中的矢量坐标转换到轨道坐标系，则为：
[0082] qo＝L001Rz(ω+v)Rx(i)Rz(Ω)p (5)
[0083] 式中：qo——各量在轨道坐标系中的矢量坐标；
[0084] L001——坐标轴的反向变换矩阵，
[0085] 为了便于目标特性分析，将各量在轨道坐标系中的矢量坐标转换到本体坐标系，则为：
[0086] q＝Aqo (6)
[0087] 式中：A——卫星姿态矩阵。
[0088] 对于三轴稳定卫星，常采用zxy转序，则姿态矩阵由旋转变换矩阵表示为：
[0089] (7)
[0090] 式中：Ψ——卫星绕偏航轴转动的角度，即偏航角。向右偏转为正，向左偏转为负；
[0091] θ——卫星绕俯仰轴转动的角度，即俯仰角。向上为正又称抬头，向下为负又称低头；
[0092] ——卫星绕滚动轴转动的角度，即滚动角。逆时针为正，顺时针为负。
[0093] Ψ，θ与 称为卫星姿态描述的欧拉角。
[0094] 为了使本体坐标系与图像坐标有效转换，可进行如下变换(图2)：
[0095] q＝[Xc,YC,ZC]。
[0096] 定义OYZ为图像坐标系，O'YCZC为相机焦平面坐标系，O’在OTZ下的坐标为(y0,zo)。
[0097] 焦面上任意一点p坐标为(yc,zc)，则相机焦平面转换到图像坐标系下的坐标为：
[0098] [(yc－yo)/D，(zc－zo)/D] (8)。
[0099] 前一步最终的结果可整理的q＝[Xc,YC,ZC]，其中D为像元尺寸，带入公式(8)，即可得到在图像中的坐标形式。
[0100] 因此，给定观测目标、目标卫星的轨道参数、连续图像的时间间隔、像元尺寸及目标卫星姿态特性的欧拉角后，即可算的目标点的理想位置。
[0101] 3.2 基于同面特征点搜索获取基础变换矩阵方法获取目标点的实际位置[0102] 3.2.1 配准区的获取
[0103] 为了提高后续配准算法的准确度及提高效率，需要首先获得算法的实施区域。为此，可对单幅图像进行如下过程的处理：二值化、开运算、闭运算、提取轮廓，确定配准区等步骤，如图4所示。
[0104] 通过二值化和开运算，可以剔除星光背景。
[0105] 通过闭运算形成连续闭合的图像，通过轮廓提取和最长的轮廓搜索选取外轮廓。
[0106] 3.2.2 获取基础变换矩阵
[0107] 为了找到同面配准特征点，需要将本帧图像与下一帧图像进行配准。对同面配准特征点进行多次筛选获得同面特征点，最后基于同面特征点获取基础变换矩阵，如图3所示。
[0108] 第一步是采用较成熟的SIFT算法来获取图像的特征点，它具有良好的尺度、光照和空间旋转的不变性。是通过对原始图像进行高斯滤波，建立高斯金字塔，获得高斯差分来找到那感兴趣的关键点，也就是在以后的图像匹配中可能发挥作用的特征点。而这些特征点需要经过稳定性的筛选才能成为真正的关键点。再根据每一个关键点邻域内点的梯度方向，为关键点分配梯度方向和梯度的模，最终特征点以多维向量的形式描述出来。图4显示了一个特征点的选取过程。
[0109] 第二步，对经过SIFT算法筛选的特征点进行初步筛选，本筛选主要分为两个过程。
[0110] 一、对于各特征点Pi，获取其在前帧图上的坐标位置Pi(x,y)before，以及本帧图像上的坐标位置Pi(x,y)now，计算其位置移动的变化量：
[0111]
[0112] (9)
[0113] 剔除移动距离超过设定阀值Dk的特征点，阀值Dk设定可取：
[0114]
[0115] (10)
[0116] 式中，n为配准特征点的总数，t为阀值系数，一般t>1。
[0117] 剔除移动距离较大点后，剩余特征点进行如下变换：
[0118] 剩余特征点与目标点的距离 ，PT(x,y)为目标点在本帧图像中的位置坐标。
[0119]
[0120] (11)
[0121] 取 满足 的特征点，其中m为剩余特征点的总数，h为阀值系数，一般h<1。经过与标记点的距离筛选后，获得较近的配准特征点。
[0122] 对筛选后的配准特征点，取离标记点最近的k个点，以周边三个像素为半径，取各配准特征点周围的小图像，并计算特征点周围对应小图的直方图，对各特征点的直方图，与其他的配准点所对应的直方图进行相关系数计算：
[0123]
[0124] 其中： N等于直方图中数据的个数。
[0125] 随后进行求和，并取最大值对应的配准特征点为基准点Po(x,y)。
[0126]
[0127] 然后以基准点Po(x,y)对其他各配准特征点对应的直方图求相关系数，取相关系数大于阀值x>0.9的认为是同面配准特征点。经过直方图系数计算获取的基准点，认定为同面配准特征点。
[0128] 通过直方图的相关系数搜索同面配准特征点的物理基础在于：对于同面特征点的纹理相近。方法先设法找到与所有配准特征点纹理最相近的点作为基准点，然后寻找与该基准点纹理接近的点作为同面点。结合之前的筛选，该方法对于同面点搜索具有一定的适应性。
[0129] 当确定了共面配准特征点后，通过成熟的RANSAC方法实现两帧图像间的基础变换矩阵，再将前帧的目标位置通过基础变换矩阵变换到本帧上，获取目标点在本帧图像中的坐标。
[0130] 4、角度误差的求取
[0131] 经过前面的过程，即可求得图像中目标点的实际位置与理想位置，在已知相机的焦距的情况下，可得到以角度形式表示的反演误差，如图5所示。
[0132] 可见，误差用角度表示的形式为：
[0133] 其中f为相机焦距，a为理想点距离光轴距离，b为实际点距离光轴距离。
