CN110411459A - 一种使用图像序列的空间非合作目标自旋参数的辨识方法 - Google Patents

Abstract

本发明公开的一种使用图像序列的空间非合作目标自旋参数的辨识方法，首先拍摄目标航天器的图像，其次将得到的每副图像进行特征峰检测处理，得到帆板的图像坐标，然后，根据得到的帆板的顶点图像坐标反解出目标航天器的可能姿态，最后将多副图像处理得到的姿态序列拟合得到自旋轴和转速。本发明公开的方法适用于有较少先验几何信息(仅需要特征的长宽比)的非合作目标，限于其帆板和本体不存在相对转动的情形。

Description

一种使用图像序列的空间非合作目标自旋参数的辨识方法

技术领域

本发明属于航天导航技术领域，具体涉及一种使用图像序列的空间非合作目标自旋参数的辨识方法。

背景技术

根据运动稳定性理论，空间失效目标大多数趋于绕最大惯量轴旋转状态。对于开展太空垃圾清理，故障卫星维修等空间操作/在轨维护任务来说，要先确定目标的自旋轴和转速，任务航天器才能够通过姿态控制实现与目标本体的相对静止，从而实现进一步的抓捕、对接、消旋及拖曳等动作。目前通过光学图像辨识目标自旋参数的技术，一般只能通过图像的变化规律获得非合作目标的转动周期，而周期信息不足以支持对目标的姿态捕获等控制。也有文献使用图像进行目标的姿态估计，但计算中需要较多的先验尺寸信息，不适用于信息不完整的非合作目标。

发明内容

本发明的目的是提供一种使用图像序列的空间非合作目标自旋参数的辨识方法，通过应用光学图像的成像投影关系，可在仅有目标尺寸长宽比的情况下计算目标的几种可能姿态，解决了先验信息不足时的姿态估计问题；通过引入多幅连续图像和差分进化算法，解决了从可能姿态中选取真实姿态的问题，同时获得目标的自旋参数。

本发明所采用的技术方案是，一种使用图像序列的空间非合作目标自旋参数的辨识方法，具体实施过程包括如下步骤：

步骤1，拍摄目标航天器的图像；

控制携带相机的航天器接近目标至相机拍摄范围内，拍摄过程中航天器利用自身姿态反馈控制使得相机光轴对准目标，设置相机按照小于目标转动周期的五分之一，拍摄得到清晰的目标航天器的图像序列；

步骤2，将步骤1得到的每副图像进行特征峰检测处理，得到帆板的图像坐标；

对于每一副图像，首先使用Sobel算子检测边缘，然后利用随机Hough变换检测所有的直线，最后将检测出的直线进行平行、闭合等匹配，选择出最大的平行四边形即认为是帆板，计算获得帆板4个顶点的图像坐标

步骤3，根据步骤2得到的帆板的顶点图像坐标反解出目标航天器的可能姿态；

步骤4，将多副图像处理得到的姿态序列拟合得到自旋轴和转速。

本发明的其他特点还在于，

步骤3的具体操作过程如下：

定义成像坐标系O-XYZ，其中，相机中心为原点，光轴为Z轴，相机水平面为X轴，Y轴根据X轴和Z轴叉乘获得，则成像坐标系到J2000惯性系O-X_gY_gZ_g的坐标转换矩阵M根据相机所在航天器的轨道和姿态计算得到：

其中，r和v分别为航天器惯性系位置和速度，M_cam为已知的相机安装矩阵，M_body为航天器本体姿态矩阵，q＝[q₁,q₂,q₃,q₄]为四元数，即相机所在航天器的轨道系姿态，e是单位矢量符号；

定义目标本体系O-X_bY_bZ_b，其中Y_b轴为帆板长度方向，X_b轴为帆板宽度方向，Z_b为帆板法向，设帆板实际长为h，宽为w，则目标航天器帆板顶点在目标本体系的坐标p_i,i＝1,2,3,4如式(2)所示：

由于成像距离较远且不需要绝对尺度信息，将相机成像投影关系简化为平行投影，则有

其中，p_i'为计算中间变量，p_i ^p表示某一姿态下帆板顶点通过投影关系计算出的在图像上坐标，与步骤2中的通过图像处理获得的坐标在没有误差的理想情况下应当相同，k为单位长度对应的像素行数，A为成像坐标系到目标本体系的转换矩阵，使用四元数方式表述为q＝[q₁,q₂,q₃,q₄]，则转换矩阵A的表达式如式(4)所示：

从式(3)中得到，经过投影后丢失了Z方向的信息，帆板相邻边的夹角会发生变化，设分别为θ_4-1-2,θ_1-2-3,θ_2-3-4,θ_3-4-1，cosθ_1-2-3计算公式如式(5)所示，其余三个可按照下标类推：

其中，p_i ^p表示某一姿态下帆板顶点通过投影关系计算出的在图像上坐标，与步骤2中的通过图像处理获得的坐标在没有误差的理想情况下应当相同；

投影后帆板与成像坐标系X轴和Y轴不再重合，定义每条边与X轴夹角分别为θ_{1_2}，θ_{2_3}，θ_{3_4}，θ_{4_1}，tanθ_{1_2}计算公式如式(6)所示，其余三个可按照下标类推：

其中，p_i ^p表示某一姿态下帆板顶点通过投影关系计算出的在图像上的坐标；

将步骤2中得到的与式(2)计算得到的p_i ^p带入式(4)和式(5)中，得到8个方程构成的方程组，将方程组中的k消除，仅需给定h/w，然后求解方程组可得对应的目标航天器的姿态四元数q；

求解方程组的具体过程如下：

首先定义3维求解变量trial[3]，则目标航天器的姿态四元数q表示为：

代入上述关系式(5)-式(6)，使得按照投影及姿态关系计算出的角度和图像检测计算的值尽可能接近，优化的目标函数如式(8)所示：

应用差分进化算法，多次设置不同的初值，即可获得目标航天器在成像坐标系下的所有四种可能姿态，将获得的姿态四元数按照式(3)转为矩阵A，则惯性系的姿态矩阵如式9所示：

A_e＝M·A (9)

将式(1)中的M带入式(14)，再按照式(2)逆向计算得到惯性系下目标航天器的姿态四元数q。

步骤4的具体过程如下：

针对绕固定轴匀速转动的目标，按照一定时间间隔拍摄得到的图像序列，对应的姿态理论上应当满足式(10)所示的关系：

其中，q(t_m)为某一时刻姿态，q(t_m+1)为下一时刻姿态，[n_x,n_y,n_z]为转轴，ω为转速，Δt为时间间隔，由于图像拍摄和处理的误差，实际计算得到的姿态并不能严格满足该关系式，只能使拟合偏差最小；

待优化的变量如下：设图像序列个数为N，则需要辨识N个图像对应的四元数在可能取值中的序号，即[1，2，3，4]之间的整数，同时还有待拟合的转轴和转速，转轴为单位矢量2维，转速1维，以及序列的初始姿态q(t₀)，序列的初始姿态为3维，总的优化变量X维数为N+6，优化目标为使得实际图像解算姿态与根据初值和匀速转动计算的理论姿态的夹角偏差α_m的平方和最小，目标函数如式(11)所示：

α_m按照四元数的运算规则计算：

其中q(t_m)为图像解算的姿态取值，q'(t_m)为根据初始姿态和转速计算的值：

为实现统一实数编码，以实数变量取整的方式确定真实姿态四元数在可能姿态中的取值，设可能姿态为则用来拟合计算的姿态如式(14)所示：

其中，floor为向下取整函数，由于X_i定义域为[0,1)，floor(2X_i+2)取值是0，1，2，3；

将其代入式(13)和式(12)，最后再带入式(11)，使用混合差分进化算法进行求解，优化使得实际图像解算姿态与根据初值和匀速转动计算的理论姿态的夹角偏差α_m的平方和最小，此时[n_x,n_y,n_z]和w即分别为优化出的惯性坐标系的转轴和转速。

本发明的有益效果是，一种使用图像序列的空间非合作目标自旋参数的辨识方法，适用于有较少先验几何信息(仅需要特征的长宽比)的非合作目标，限于其帆板和本体不存在相对转动的情形。根据对目标拍摄的图像，利用图像处理获得帆板的特征像素，根据姿态转换和投影关系反解出可能姿态，最后将多个图像解算的姿态拟合为目标的自旋参数，获得转轴和转速。

附图说明

图1是本发明的一种使用图像序列的空间非合作目标自旋参数的辨识方法的流程图；

图2是图像处理过程中的目标航天器的示意图；

图3是目标航天器计算的多个姿态的示意图；

图4是实施例中对目标拍摄的仿真图像。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明的一种使用图像序列的空间非合作目标自旋参数的辨识方法，如图1所示，具体实施过程包括如下步骤：

步骤1，拍摄目标航天器的图像；

控制携带相机的航天器接近目标至相机拍摄范围内，拍摄过程中航天器利用自身姿态反馈控制使得相机光轴对准目标，设置相机按照小于目标转动周期的五分之一，拍摄得到清晰的目标航天器的图像序列；

步骤2，将步骤1得到的每副图像进行特征峰检测处理，得到帆板的图像坐标；

对于每一副图像，首先使用Sobel算子检测边缘，然后利用随机Hough变换检测所有的直线，最后将检测出的直线进行平行、闭合等匹配，选择出最大的平行四边形即认为是帆板，计算获得帆板4个顶点的图像坐标如图2所示；

Sobel算子如下，也可使用其它边缘检测算子；

随机Hough变换的步骤为：

1)从边缘图像中随机选择两个点，计算对应的直线参数，即y＝ax+b；

2)计算剩下的像素点与该直线的距离，以判断匹配程度(距离小于δ)，如果匹配的像素点较多(大于给定阈值)则说明图中存在一条直线，否则返回1)重新选择随机点；

3)从图像中剔除已检测为直线的像素点，返回1)重新选择随机点继续检测下一条直线，直到剩下的所有像素点都不能匹配为直线；

获得图像中的直线集合以后，再按照排列组合将所有4条直线的组合进行判断，判据为：

1)直线两两平行(夹角小于给定阈值)

2)不平行的直线端点相接(端点距离小于给定阈值)

步骤3，根据步骤2得到的帆板的顶点图像坐标反解出目标航天器的可能姿态；

步骤3的具体操作过程如下：

定义成像坐标系O-XYZ，其中，相机中心为原点，光轴为Z轴，相机水平面为X轴，Y轴根据X轴和Z轴叉乘获得，则成像坐标系到J2000惯性系O-X_gY_gZ_g的坐标转换矩阵M根据相机所在航天器的轨道和姿态计算得到：

其中，r和v分别为航天器惯性系位置和速度，M_cam为已知的相机安装矩阵，M_body为航天器本体姿态矩阵，q＝[q₁,q₂,q₃,q₄]为四元数，即相机所在航天器的轨道系姿态，e是单位矢量符号；

定义目标本体系O-X_bY_bZ_b，其中Y_b轴为帆板长度方向，X_b轴为帆板宽度方向，Z_b为帆板法向，设帆板实际长为h，宽为w，则目标航天器帆板顶点在目标本体系的坐标p_i,i＝1,2,3,4，如式(2)所示：

由于成像距离较远且不需要绝对尺度信息，将相机成像投影关系简化为平行投影，则有

其中，p_i'为计算中间变量，p_i ^p表示某一姿态下帆板顶点通过投影关系计算出的在图像上坐标，与步骤2中的通过图像处理获得的坐标在没有误差的理想情况下应当相同，k为单位长度对应的像素行数，A为成像坐标系到目标本体系的转换矩阵，使用四元数方式表述为q＝[q₁,q₂,q₃,q₄]，则转换矩阵A的表达式如式(4)所示：

从式(3)中得到，经过投影后丢失了Z方向的信息，帆板相邻边的夹角会发生变化，设分别为θ_4-1-2,θ_1-2-3,θ_2-3-4,θ_3-4-1，cosθ_1-2-3计算公式如式(5)所示，其余三个可按照下标类推：

其中，p_i ^p表示某一姿态下帆板顶点通过投影关系计算出的在图像上坐标，与步骤2中的通过图像处理获得的坐标在没有误差的理想情况下应当相同；

投影后帆板与成像坐标系X轴和Y轴不再重合，定义每条边与X轴夹角分别为θ_{1_2}，θ_{2_3}，θ_{3_4}，θ_{4_1}，tanθ_{1_2}计算公式如式(6)所示，其余三个可按照下标类推：

其中，p_i ^p表示某一姿态下帆板顶点通过投影关系计算出的在图像上坐标；

将步骤2中得到的与式(3)计算得到的p_i ^p带入式(4)和式(5)中，得到8个方程构成的方程组，将方程组中的k消除，仅需给定h/w，然后求解方程组可得对应的目标航天器的姿态四元数q；

如图3所示，由于三维坐标到二维坐标的映射损失了目标相对于相机的距离信息，因此导致有四种可能姿态都满足式(4)和式(5)构成的方程组。

使用差分进化算法求解方程组的具体过程如下：

首先，定义3维求解变量trial[3]，则目标航天器的姿态四元数q表示为：

代入上述关系式(5)和式(6)，使得按照投影及姿态关系计算出的角度和图像检测计算的值尽可能接近，优化的目标函数如式(8)所示：

应用差分进化算法，多次设置不同的初值，即可获得目标航天器在成像坐标系下的所有四种可能姿态，将获得的姿态四元数按照式(3)转为矩阵A，则惯性系的姿态矩阵如式9所示：

A_e＝M·A (9)

将式(1)中的M带入式14，再按照式(2)逆向计算得到惯性系下目标航天器的姿态四元数q。

步骤4，将多副图像处理得到的惯性系姿态序列拟合得到自旋轴和转速；

步骤4的具体过程如下：

针对绕固定轴匀速转动的目标，按照一定时间间隔拍摄得到的图像序列，对应的姿态理论上应当满足式(10)所示的关系：

其中，q(t_m)为某一时刻姿态，q(t_m+1)为下一时刻姿态，[n_x,n_y,n_z]为转轴，ω为转速，Δt为时间间隔，由于图像拍摄和处理的误差，实际计算得到的姿态并不能严格满足该关系式，只能使拟合偏差最小；

待优化的变量如下：设图像序列个数为N，则需要辨识N个图像对应的四元数在可能取值中的序号，即[1，2，3，4]之间的整数，同时还有待拟合的转轴和转速，转轴为单位矢量2维，转速1维，以及序列的初始姿态q(t₀)，序列的初始姿态为3维，总的优化变量X维数为N+6，目标函数如式(11)所示：

α_m按照四元数的运算规则计算，如式12所示：

其中q(t_m)为图像解算的姿态取值，q'(t_m)为根据初始姿态和转速计算的值：

为实现统一实数编码，以实数变量取整的方式确定真实姿态四元数在可能姿态中的取值，设可能姿态为则用来拟合计算的姿态如式(14)所示：

其中，floor为向下取整函数，由于X_i定义域为[0,1)，floor(2X_i+2)取值就是0，1，2，3；

将其代入式(13)和式(12)，最后再带入式(11)得到目标函数，使用混合差分进化算法进行求解，优化使得实际图像解算姿态与根据初值和匀速转动计算的理论姿态的夹角偏差α_m的平方和最小，此时[n_x,n_y,n_z]和w即分别为优化出的惯性坐标系的转轴和转速。

本发明通过应用光学图像的成像投影关系，可在仅有目标尺寸长宽比的情况下计算目标的几种可能姿态，解决了现有方法先验信息不足时的姿态估计难题；根据目标为自旋状态这一约束，通过引入多幅连续图像处理数据和差分进化算法，解决了单幅图像计算姿态存在多解的问题。

实施例

首先获得对目标拍摄的图像，如图4所示，本发明中使用仿真生成航天器连续拍摄的自旋目标图像，时间间隔1s，已知帆板长宽比为310/56；

根据步骤2，进行图像处理，识别的帆板顶点坐标见表1；

根据步骤4，应用混合编码差分进化算法，总的优化变量X维数为12，种群数设为100，进化代数设为800，解得转轴方向为：n_x＝-0.825179；n_y＝0.067050；n_z＝-0.560878；转速ω＝1.0472rad/s，初始时刻姿态四元数为q＝[0.109931,0.455528,0.276139,0.83914]。

表1图像处理解得的帆板顶点坐标

根据步骤3，计算得到的可能姿态四元数见表2；

表2姿态四元数计算结果

Claims (3)

1.一种使用图像序列的空间非合作目标自旋参数的辨识方法，其特征在于，具体实施过程包括如下步骤：

步骤1，拍摄目标航天器的图像；

控制携带相机的航天器接近目标至相机拍摄范围内，拍摄过程中航天器利用自身姿态反馈控制使得相机光轴对准目标，设置相机按照小于目标转动周期的五分之一，拍摄得到清晰的目标航天器的图像序列；

步骤2，将步骤1得到的每副图像进行特征峰检测处理，得到帆板的图像坐标；

对于每一副图像，首先使用Sobel算子检测边缘，然后利用随机Hough变换检测所有的直线，最后将检测出的直线进行平行、闭合等匹配，选择出最大的平行四边形即认为是帆板，计算获得帆板4个顶点的图像坐标

步骤3，根据步骤2得到的帆板的顶点图像坐标反解出目标航天器的可能姿态；

步骤4，将多副图像处理得到的姿态序列拟合得到自旋轴和转速。

2.如权利要求1所述的一种使用图像序列的空间非合作目标自旋参数的辨识方法，其特征在于，所述步骤3的具体操作过程如下：

定义成像坐标系O-XYZ，其中，相机中心为原点，光轴为Z轴，相机水平面为X轴，Y轴根据X轴和Z轴叉乘获得，则成像坐标系到J2000惯性系O-X_gY_gZ_g的坐标转换矩阵M根据相机所在航天器的轨道和姿态计算得到：

其中，r和v分别为航天器惯性系位置和速度，M_cam为已知的相机安装矩阵，M_body为航天器本体姿态矩阵，q＝[q₁,q₂,q₃,q₄]为四元数，即相机所在航天器的轨道系姿态，e是单位矢量符号；

定义目标本体系O-X_bY_bZ_b，其中Y_b轴为帆板长度方向，X_b轴为帆板宽度方向，Z_b为帆板法向，设帆板实际长为h，宽为w，则目标航天器帆板顶点在目标本体系的坐标p_i,i＝1,2,3,4如式(2)所示：

由于成像距离较远且不需要绝对尺度信息，将相机成像投影关系简化为平行投影，则有

其中，p_i'为计算中间变量，p_i ^p表示某一姿态下帆板顶点通过投影关系计算出的在图像上坐标，与步骤2中的通过图像处理获得的坐标在没有误差的理想情况下应当相同，k为单位长度对应的像素行数，A为成像坐标系到目标本体系的转换矩阵，使用四元数方式表述为q＝[q₁,q₂,q₃,q₄]，则转换矩阵A的表达式如式(4)所示：

从式(3)中得到，经过投影后丢失了Z方向的信息，帆板相邻边的夹角会发生变化，设分别为θ_4-1-2,θ_1-2-3,θ_2-3-4,θ_3-4-1，cosθ_1-2-3计算公式如式(5)所示，其余三个可按照下标类推：

其中，p_i ^p表示某一姿态下帆板顶点通过投影关系计算出的在图像上坐标，与步骤2中的通过图像处理获得的坐标在没有误差的理想情况下应当相同；

投影后帆板与成像坐标系X轴和Y轴不再重合，定义每条边与X轴夹角分别为θ_{1_2}，θ_{2_3}，θ_{3_4}，θ_{4_1}，tanθ_{1_2}计算公式如式(6)所示，其余三个可按照下标类推：

其中，p_i ^p表示某一姿态下帆板顶点通过投影关系计算出的在图像上坐标；

将步骤2中得到的与式(2)计算得到的p_i ^p带入式(4)和式(5)中，得到8个方程构成的方程组，将方程组中的k消除，仅需给定h/w，然后求解方程组可得对应的目标航天器的姿态四元数q；

求解方程组的具体过程如下：

首先定义3维求解变量trial[3]，则目标航天器的姿态四元数q表示为：

代入上述关系式(5)-式(6)，使得按照投影及姿态关系计算出的角度和图像检测计算的值尽可能接近，优化的目标函数如式(8)所示：

应用差分进化算法，多次设置不同的初值，即可获得目标航天器在成像坐标系下的所有四种可能姿态，将获得的姿态四元数按照式(3)转为矩阵A，则惯性系的姿态矩阵如式9所示：

A_e＝M·A (9)

将式(1)中的M带入式(14)，再按照式(2)逆向计算得到惯性系下目标航天器的姿态四元数q。

3.如权利要求1所述的一种使用图像序列的空间非合作目标自旋参数的辨识方法，其特征在于，所述步骤4的具体过程如下：

针对绕固定轴匀速转动的目标，按照一定时间间隔拍摄得到的图像序列，对应的姿态理论上应当满足式(10)所示的关系：

其中，q(t_m)为某一时刻姿态，q(t_m+1)为下一时刻姿态，[n_x,n_y,n_z]为转轴，ω为转速，Δt为时间间隔，由于图像拍摄和处理的误差，实际计算得到的姿态并不能严格满足该关系式，只能使拟合偏差最小；

待优化的变量如下：设图像序列个数为N，则需要辨识N个图像对应的四元数在可能取值中的序号，即[1，2，3，4]之间的整数，同时还有待拟合的转轴和转速，转轴为单位矢量2维，转速1维，以及序列的初始姿态q(t₀)，序列的初始姿态为3维，总的优化变量X维数为N+6，优化目标为使得实际图像解算姿态与根据初值和匀速转动计算的理论姿态的夹角偏差α_m的平方和最小，目标函数如式(11)所示：

α_m按照四元数的运算规则计算：

其中q(t_m)为图像解算的姿态取值，q'(t_m)为根据初始姿态和转速计算的值：

为实现统一实数编码，以实数变量取整的方式确定真实姿态四元数在可能姿态中的取值，设可能姿态为则用来拟合计算的姿态如式(14)所示：

其中，floor为向下取整函数，由于X_i定义域为[0,1),floor(2X_i+2)取值是0，1，2，3；

将其代入式(13)和式(12)，最后再带入式(11)，使用混合差分进化算法进行求解，优化使得实际图像解算姿态与根据初值和匀速转动计算的理论姿态的夹角偏差α_m的平方和最小，此时[n_x,n_y,n_z]和w即分别为优化出的惯性坐标系的转轴和转速。