CN112559959B - 基于特征向量的天基成像非合作目标旋转状态解算方法 - Google Patents

基于特征向量的天基成像非合作目标旋转状态解算方法 Download PDF

Info

Publication number
CN112559959B
CN112559959B CN202011440247.7A CN202011440247A CN112559959B CN 112559959 B CN112559959 B CN 112559959B CN 202011440247 A CN202011440247 A CN 202011440247A CN 112559959 B CN112559959 B CN 112559959B
Authority
CN
China
Prior art keywords
rotation
coordinate system
target
axis
coordinates
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
Application number
CN202011440247.7A
Other languages
English (en)
Other versions
CN112559959A (zh
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.)
China Xian Satellite Control Center
Original Assignee
China Xian Satellite Control Center
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 China Xian Satellite Control Center filed Critical China Xian Satellite Control Center
Priority to CN202011440247.7A priority Critical patent/CN112559959B/zh
Publication of CN112559959A publication Critical patent/CN112559959A/zh
Application granted granted Critical
Publication of CN112559959B publication Critical patent/CN112559959B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Image Analysis (AREA)

Abstract

本发明提供了一种基于特征向量的天基成像非合作目标旋转状态解算方法，以特征向量旋转状态代替目标旋转状态，通过提取目标连续图像中的特征向量投影，建立特征向量投影变换模型，解算成像目标相对相机坐标系的旋转轴指向与转速。本发明具有无需目标特征尺寸、普适性好的优点。

Description

基于特征向量的天基成像非合作目标旋转状态解算方法
技术领域
本发明属于飞行器导航与制导领域，适用于对天基成像非合作目标旋转轴指向和旋转角速度的解算。
背景技术
搭载于卫星的高分辨率光学探测设备，具有不受空间天气影响、灵活机动的先天优势，能够完成对太空目标的观察探测。随着天基成像系统应用范围的拓展，从航天器碰撞规避预警、航天器在轨故障维修需求考虑，对天基成像目标的形态识别、运动特征分析、空间状态确定的要求越来越高。综上所述，需要一种天基成像非合作目标的旋转状态解算方法。
传统的目标旋转状态解算方法需已知目标物上特征直线的尺寸及指向，通过图像识别获取特征直线在图像中的坐标，两者匹配解算旋转矩阵，由不同时刻的旋转矩阵解算目标的旋转轴和旋转速度，不适用于非合作目标旋转状态解算。
发明内容
为了克服现有技术的不足，本发明提供一种基于特征向量的天基成像非合作目标旋转状态解算方法，以特征向量旋转状态代替目标旋转状态，通过提取目标连续图像中的特征向量投影，建立特征向量投影变换模型，解算成像目标相对相机坐标系的旋转轴指向与转速，具有无需目标特征尺寸、普适性好的优点。
本发明解决其技术问题所采用的技术方案包括以下步骤：
1)选取与目标航天器本体固连的特征向量来代替目标航天器本身，以特征向量旋转状态代替目标旋转状态；选取不同时刻的目标图像，对图像中的特征向量进行提取，得到特征向量的投影矢量；所述的目标图像能够覆盖目标360度视角成像；
2)根据相机投影变换模型和欧拉轴/角旋转矩阵建立特征向量投影变换模型；
3)利用不同时刻的特征向量进行椭圆拟合，计算待求参数初值；
4)利用最小二乘平差方法解算成像目标相对相机坐标系的旋转轴指向与转速。
所述的步骤2)以相机的光学中心为坐标原点，x轴和y轴分别与CCD成像平面的两条垂直边平行，z轴与相机的光轴重合，建立相机坐标系；将相机坐标系平移至以目标航天器的旋转中心为原点，建立旋转中心坐标系；根据相似投影关系得到
其中，(xX,yX,zX)为物点P1在相机坐标系下的坐标，(xC,yC)为物点投影位置P2在像面坐标系下的物理单位坐标，f为相机的焦距；
设目标旋转中心为O，过旋转中心O的旋转轴在旋转中心坐标系中指向的天球坐标为(α,δ)，特征向量为AB，点A和B在旋转中心坐标系下的初始坐标分别为SA0和SB0，特征向量AB在旋转中心坐标系下的坐标SAB0＝SB0-SA0＝[x,y,z]T，转轴用单位向量表示为其中，α,δ为旋转轴指向的方位角和高度角，目标绕旋转轴e以w的角速度进行转动，经过时间t目标旋转的角度Φ＝wt；
由欧拉轴/角旋转矩阵可得，绕旋转轴e旋转Φ角的旋转矩阵为
则t时刻点A在旋转中心坐标系下的坐标SA＝RZJ T·SA0；
设t时刻旋转中心O在相机坐标系下的坐标为SO,X，t时刻点A在相机坐标系下的坐标SA,X＝SA+SO,X＝RZJ T·SA0+SO,X；t时刻点P在像平面上的投影坐标为
式中SA,X(3)、SB,X(3)分别表示点A、B在相机坐标系下的z轴坐标；设为常值μ，/>其中，/>
特征向量AB在像平面的投影坐标
令则有/>
t＝0时特征向量AB在像平面上的投影坐标
设待求参数X＝(α,δ,w,zμ)T，其中α,δ为旋转轴指向的方位角和高度角，w为旋转速度，zμ为构造参数，t时刻特征向量AB在像平面的投影坐标观测值为SAB,O＝(xAB,O,yAB,O)T；
建立误差方程误差方程线性化/>其中L为自由项，L＝SAB,O-SAB,C(X)。
所述的步骤3)将不同时刻的特征向量画在同一坐标系下，对这些向量的终点进行椭圆拟合，计算旋转轴指向天球坐标的初值α＝θ，其中，a、b分别为拟合椭圆的半长轴、半短轴，θ为椭圆短轴相对坐标系x轴的转角；用不同时刻特征向量间的夹角除以间隔时间得到转速w的初值，zμ的初值选为0。
所述的步骤4)将待求参数的初值作为估计值X，代入线性化误差方程计算出待求参数改正数△X，将改正后的估计值X+ΔX作为新的估计值X再代入线性化误差方程，反复迭代计算，直至△X小于设定精度要求，得到待求参数的最优估值。
本发明的有益效果是：
1)对不同旋转轴指向和转速均能准确解算，旋转轴指向偏差均在0.58°以内，转速误差均在0.0006°/s以内。
2)对任一成像相机和目标航天器都能适用，无需知道目标特征尺寸。
附图说明
图1是小孔成像原理图；
图2是相机坐标系与旋转中心坐标系示意图；
图3是特征向量椭圆拟合；
图4是本发明的算法流程图。
具体实施方式
下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。
本发明包括以下步骤：
1)选取与目标航天器本体固连的特征向量来代替目标航天器本身，以特征向量旋转状态代替目标旋转状态。选取6个以上的不同时刻目标图像，这些图像必须能够覆盖目标旋转360度，对这些图像中的特征向量进行提取，得到特征向量的投影矢量。
2)根据相机投影变换模型和欧拉轴/角旋转矩阵建立特征向量投影变换模型。
3)利用不同时刻的特征向量进行椭圆拟合，计算待求参数初值。
4)利用最小二乘平差方法解算成像目标相对相机坐标系的旋转轴指向与转速。
上述步骤的具体过程如下：
1)选取特征向量，对航天器图像中的特征向量进行提取。
选取与目标航天器本体固连的特征向量来代替目标航天器本身，以特征向量旋转状态代替目标旋转状态。以某卫星为例，其帆板已停止转动，与卫星本体固连，可以其南侧第一块太阳电池阵的对角线AB为特征向量。
对九张不同时刻该卫星图像中的特征向量进行提取，得到特征向量的投影如下表所示。
表1不同时刻特征向量AB在图像中的像素坐标
t(s) X(pixel) Y(pixel)
0 -242 -16
99.941 -197 -278
201.14 -30 -369
401.582 242 98
501.623 128 340
601.448 -92 311
701.406 -238 32
801.452 -194 -270
901.519 8 -351
2)建立目标旋转状态解算方程
如图1所示，以相机的光学中心为坐标原点，x轴和y轴分别与CCD成像平面的两条垂直边平行，z轴与相机的光轴重合，建立相机坐标系。将相机坐标系平移至以目标航天器的旋转中心为原点，建立旋转中心坐标系。相机坐标系和旋转中心坐标系如图2所示。根据相似投影关系可得
其中，(xX,yX,zX)为物点P1在相机坐标系下的坐标，(xC,yC)为物点投影位置P2在像面坐标系下的物理单位坐标，f为相机的焦距。
设目标旋转中心为O，过旋转中心O的旋转轴在旋转中心坐标系中指向的天球坐标为(α,δ)。特征向量为AB，点A和B在旋转中心坐标系下的初始坐标分别为SA0和SB0，特征向量AB在旋转中心坐标系下的坐标为
SAB0＝SB0-SA0＝[x,y,z]T (2)
转轴用单位向量表示为
其中，α,δ为旋转轴指向的方位角和高度角。目标绕旋转轴e以w的角速度进行转动，经过时间t目标旋转的角度为
Φ＝wt (4)
由欧拉轴/角旋转矩阵可得，绕旋转轴e旋转Φ角的旋转矩阵为
则t时刻点A在旋转中心坐标系下的坐标表示为
SA＝RZJ T·SA0 (1)
设t时刻旋转中心O在相机坐标系下的坐标为SO,X，t时刻点A在相机坐标系下的坐标为
SA,X＝SA+SO,X＝RZJ T·SA0+SO,X (2)
由式(1)可得，t时刻点P在像平面上的投影坐标为
公式中SA,X(3)、SB,X(3)分别表示点A、B在相机坐标系下的z轴坐标。由式(7)可得，t时刻点A和B在像平面上的投影坐标为
由于t时刻目标与相机的距离远大于相机的焦距，则
SA,X(3)≈SB,X(3)＞＞f (10)
故可设
因为目标与相机的距离随时间变化，故μ(t)为关于t的函数。但在较短时间内目标与相机的距离变化量远小于该距离，故可将μ(t)视为常值μ。
则式(8)(9)可写为
其中
特征向量AB在像平面的投影坐标计算公式为
SAB,C＝SB,C-SA,C＝Rμ·RZJ T·(SB0-SA0)＝Rμ·RZJ T·SAB0 (14)
将式(8)(13)代入，得
对式(15)，令
xμ＝μx
yμ＝μy (16)
zμ＝μz
则有
t＝0时，由式(3)(4)可得
代入式(20)可得
SAB,C|t＝0为t＝0时特征向量AB在像平面上的投影坐标，可由目标图像提取得到，故xμ,yμ分别为-242像素和-16像素。设待求参数X＝(α,δ,w,zμ)T，其中α,δ为旋转轴指向的方位角和高度角，w为旋转速度，zμ为构造参数，t时刻特征向量AB在像平面的投影坐标观测值为SAB,O＝(xAB,O,yAB,O)T，见表1。
建立误差方程
误差方程线性化
其中L为自由项，表达式为
L＝SAB,O-SAB,C(X) (22)
3)选取待求参数初值
将不同时刻的特征向量画在同一坐标系下，对这些向量的终点进行椭圆拟合，根据公式(23)、(24)计算旋转轴指向天球坐标α,δ的初值。
α＝θ (23)
其中，a、b分别为拟合椭圆的半长轴、半短轴，θ为椭圆短轴相对坐标系x轴的转角。用不同时刻特征向量间的夹角除以间隔时间得到转速w的初值，zμ的初值选为0。以上即为待求参数X的初值。
4)解算目标旋转轴指向和转速
将待求参数X的初值作为估计值X，代入式(21)中计算出待求参数改正数△X，将改正后的估计值X+ΔX作为新的估计值X再代入式(21)中，反复迭代计算，直至△X<ε(ε为精度要求)，即可得到待求参数X的最优估值(164.2199，40.8655，0.5031，276.5173)，其中包含旋转轴指向天球坐标α＝164.2°，δ＝40.9°和转速ω＝0.5°/s。

Claims (4)

1.一种基于特征向量的天基成像非合作目标旋转状态解算方法，其特征在于，包括以下步骤：
1)选取与目标航天器本体固连的特征向量来代替目标航天器本身，以特征向量旋转状态代替目标旋转状态；选取不同时刻的目标图像，对图像中的特征向量进行提取，得到特征向量的投影矢量；所述的目标图像能够覆盖目标360度视角成像；
2)根据相机投影变换模型和欧拉轴/角旋转矩阵建立特征向量投影变换模型；
3)利用不同时刻的特征向量进行椭圆拟合，计算待求参数初值；
4)利用最小二乘平差方法解算成像目标相对相机坐标系的旋转轴指向与转速。
2.根据权利要求1所述的基于特征向量的天基成像非合作目标旋转状态解算方法，其特征在于，所述的步骤2)以相机的光学中心为坐标原点，x轴和y轴分别与CCD成像平面的两条垂直边平行，z轴与相机的光轴重合，建立相机坐标系；将相机坐标系平移至以目标航天器的旋转中心为原点，建立旋转中心坐标系；根据相似投影关系得到
其中，(xX,yX,zX)为物点P1在相机坐标系下的坐标，(xC,yC)为物点投影位置P2在像面坐标系下的物理单位坐标，f为相机的焦距；
设目标旋转中心为O，过旋转中心O的旋转轴在旋转中心坐标系中指向的天球坐标为(α,δ)，特征向量为AB，点A和B在旋转中心坐标系下的初始坐标分别为SA0和SB0，特征向量AB在旋转中心坐标系下的坐标SAB0＝SB0-SA0＝[x,y,z]T，转轴用单位向量表示为其中，α,δ为旋转轴指向的方位角和高度角，目标绕旋转轴e以w的角速度进行转动，经过时间t目标旋转的角度Φ＝wt；
由欧拉轴/角旋转矩阵可得，绕旋转轴e旋转Φ角的旋转矩阵为
则t时刻点A在旋转中心坐标系下的坐标SA＝RZJ T·SA0；
设t时刻旋转中心O在相机坐标系下的坐标为SO,X，t时刻点A在相机坐标系下的坐标SA,X＝SA+SO,X＝RZJ T·SA0+SO,X；t时刻点P在像平面上的投影坐标为
式中SA,X(3)、SB,X(3)分别表示点A、B在相机坐标系下的z轴坐标；设为常值μ，/>其中，/>
特征向量AB在像平面的投影坐标
令则有/>
t＝0时特征向量AB在像平面上的投影坐标
设待求参数X＝(α,δ,w,zμ)T，其中α,δ为旋转轴指向的方位角和高度角，w为旋转速度，zμ为构造参数，t时刻特征向量AB在像平面的投影坐标观测值为SAB,O＝(xAB,O,yAB,O)T；
建立误差方程误差方程线性化/>其中L为自由项，L＝SAB,O-SAB,C(X)。
3.根据权利要求1所述的基于特征向量的天基成像非合作目标旋转状态解算方法，其特征在于，所述的步骤3)将不同时刻的特征向量画在同一坐标系下，对这些向量的终点进行椭圆拟合，计算旋转轴指向天球坐标的初值α＝θ，其中，a、b分别为拟合椭圆的半长轴、半短轴，θ为椭圆短轴相对坐标系x轴的转角；用不同时刻特征向量间的夹角除以间隔时间得到转速w的初值，zμ的初值选为0。
4.根据权利要求1所述的基于特征向量的天基成像非合作目标旋转状态解算方法，其特征在于，所述的步骤4)将待求参数的初值作为估计值X，代入线性化误差方程计算出待求参数改正数△X，将改正后的估计值X+ΔX作为新的估计值X再代入线性化误差方程，反复迭代计算，直至△X小于设定精度要求，得到待求参数的最优估值。
CN202011440247.7A 2020-12-07 2020-12-07 基于特征向量的天基成像非合作目标旋转状态解算方法 Active CN112559959B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011440247.7A CN112559959B (zh) 2020-12-07 2020-12-07 基于特征向量的天基成像非合作目标旋转状态解算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011440247.7A CN112559959B (zh) 2020-12-07 2020-12-07 基于特征向量的天基成像非合作目标旋转状态解算方法

Publications (2)

Publication Number Publication Date
CN112559959A CN112559959A (zh) 2021-03-26
CN112559959B true CN112559959B (zh) 2023-11-07

Family

ID=75061042

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011440247.7A Active CN112559959B (zh) 2020-12-07 2020-12-07 基于特征向量的天基成像非合作目标旋转状态解算方法

Country Status (1)

Country Link
CN (1) CN112559959B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113406629B (zh) * 2021-05-12 2023-07-25 北京理工大学 基于雷达长时间观测的天体目标转动估计与三维重建方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104200495A (zh) * 2014-09-25 2014-12-10 重庆信科设计有限公司 一种视频监控中的多目标跟踪方法
JP2016009391A (ja) * 2014-06-25 2016-01-18 Kddi株式会社 情報処理装置ならびにその特徴点選択方法、装置およびプログラム
CN110287873A (zh) * 2019-06-25 2019-09-27 清华大学深圳研究生院 基于深度神经网络的非合作目标位姿测量方法、系统及终端设备
CN110823214A (zh) * 2019-10-18 2020-02-21 西北工业大学 一种空间完全非合作目标相对位姿和惯量估计方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2016009391A (ja) * 2014-06-25 2016-01-18 Kddi株式会社 情報処理装置ならびにその特徴点選択方法、装置およびプログラム
CN104200495A (zh) * 2014-09-25 2014-12-10 重庆信科设计有限公司 一种视频监控中的多目标跟踪方法
CN110287873A (zh) * 2019-06-25 2019-09-27 清华大学深圳研究生院 基于深度神经网络的非合作目标位姿测量方法、系统及终端设备
CN110823214A (zh) * 2019-10-18 2020-02-21 西北工业大学 一种空间完全非合作目标相对位姿和惯量估计方法

Also Published As

Publication number Publication date
CN112559959A (zh) 2021-03-26

Similar Documents

Publication Publication Date Title
Peng et al. Virtual stereovision pose measurement of noncooperative space targets for a dual-arm space robot
CN109405835B (zh) 基于非合作目标直线与圆单目图像的相对位姿测量方法
CN106197425B (zh) 基于卫星姿态角的地面目标点位置的计算方法
CN107525492B (zh) 一种适用于敏捷对地观测卫星的偏流角仿真分析方法
CN106525001B (zh) 地球静止轨道遥感卫星相机视轴空间指向计算方法
CN111381256B (zh) 主动遥感卫星天线相位中心偏移误差计算的方法和系统
CN108663043B (zh) 基于单个相机辅助的分布式pos主子节点相对位姿测量方法
CN107564057B (zh) 顾及大气折光校正的高轨面阵光学卫星在轨几何标定方法
CN111102981B (zh) 一种基于ukf的高精度卫星相对导航方法
CN110030978B (zh) 一种全链路光学卫星几何成像模型构建方法及系统
CN110134134B (zh) 一种无人机悬停状态下的测风方法
CN112435301A (zh) 一种基于恒星轨迹的遥感相机在轨几何定标方法
Xu et al. A pose measurement method of a non-cooperative GEO spacecraft based on stereo vision
CN109599674B (zh) 一种基于解耦的相控阵天线稳定角跟踪方法
CN104406583A (zh) 双星敏感器联合确定载体姿态方法
CN112559959B (zh) 基于特征向量的天基成像非合作目标旋转状态解算方法
Barrows Videogrammetric model deformation measurement technique for wind tunnel applications
Liu et al. Hull deformation measurement for spacecraft TT&C ship by Photogrammetry
CN112179334A (zh) 基于两步Kalman滤波的星光导航方法及系统
CN107883925B (zh) 一种导航星座星间观测目标卫星图像模拟方法
CN113129377A (zh) 一种三维激光雷达快速鲁棒slam方法和装置
CN112577463B (zh) 姿态参数修正的航天器单目视觉测距方法
CN112378383B (zh) 基于圆和线特征非合作目标相对位姿双目视觉测量方法
CN116007657A (zh) 一种小视场相机与星敏测量基准同步的天文标定方法
CN112729305B (zh) 一种基于单飞行器导引头图像信息的多目标定位方法

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
GR01 Patent grant
GR01 Patent grant