CN112880669A - 一种航天器星光折射和单轴旋转调制惯性组合导航方法 - Google Patents

一种航天器星光折射和单轴旋转调制惯性组合导航方法 Download PDF

Info

Publication number
CN112880669A
CN112880669A CN202011469193.7A CN202011469193A CN112880669A CN 112880669 A CN112880669 A CN 112880669A CN 202011469193 A CN202011469193 A CN 202011469193A CN 112880669 A CN112880669 A CN 112880669A
Authority
CN
China
Prior art keywords
star
refraction
spacecraft
error
refracted
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.)
Granted
Application number
CN202011469193.7A
Other languages
English (en)
Other versions
CN112880669B (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.)
Beihang University
Original Assignee
Beihang University
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 Beihang University filed Critical Beihang University
Priority to CN202011469193.7A priority Critical patent/CN112880669B/zh
Publication of CN112880669A publication Critical patent/CN112880669A/zh
Application granted granted Critical
Publication of CN112880669B publication Critical patent/CN112880669B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/165—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/20—Instruments for performing navigational calculations
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Automation & Control Theory (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Navigation (AREA)

Abstract

本发明涉及一种航天器星光折射和单轴旋转调制惯性组合导航方法。首先根据捷联惯性导航***的误差方程建立航天器的状态模型，然后利用星敏感器获得星点像素坐标量测量并建立量测模型，使用UKF滤波估计航天器的位置、速度和姿态，对航天器的姿态误差和位置误差进行修正。本发明属于航天器自主导航领域，可为航天器提供高精度的位置及速度信息，对航天器自主导航具有重要的实际意义。

Description

一种航天器星光折射和单轴旋转调制惯性组合导航方法
技术领域
本发明属于航天器自主导航领域，涉及一种基于星点像素坐标量测量的航天器自主导航方法。
背景技术
航天器在导航、信息通讯以及军事国防等领域发挥着极其重要的作用。传统的航天器导航***需要来自地面测控站的帮助。但随着航天领域技术的加速发展，航天器的数量急剧增加，导致地面测控站的负荷快速增加。航天器实现自主运行进而实现自主导航不仅可以减轻地面测控站的压力，也可以增强航天器的在轨生存能力，提高运行的可靠性。但是由于惯性测量元件中常值偏差的影响，水平基准的精度较低，影响星光折射导航的定位精度，进而会降低组合导航***的导航精度，因此采用IMU旋转调制技术以提高捷联惯性导航***的导航精度。然而由于***初始误差的影响，单轴旋转调制式捷联惯性导航***的导航误差随时间积累，难以满足航天器对长时间、长距离和高精度自主导航任务的要求。天文自主导航技术在精度，自主性，可靠性以及抗干扰能力等很多方面具有优势，并且它的误差不受时间和距离的影响，可有效修正惯性导航的误差。
影响天文导航性能的一个比较重要的条件是地平测量精度，航天器的自主天文导航可分为直接敏感地平和利用星光折射间接敏感地平两种方法。直接敏感地平航天器自主天文导航原理简单且易于实现，但是地球敏感器测量精度低，是影响导航性能的主要原因。星光折射间接敏感地平方法只需要星敏感器本身就可以获得量测信息，且现阶段的星敏感器的测量精度要远高于地球敏感器，这可以使得导航精度有大幅的改善。因此，惯性导航和星光折射具有很好的优势互补特性，对于长时间、长距离的探测器自主导航任务，利用星光折射信息辅助惯性导航实现航天器的自主导航是一种可行的方法。
发明内容
本发明要解决的技术问题是：克服单独使用某种导航方法存在的缺点，为航天器提供一种将星光折射与单轴旋转调制惯性导航结合的自主导航方法，提供高精度的位置和速度信息。
本发明解决其技术问题所采用的技术方案为：一种航天器星光折射和单轴旋转调制惯性组合导航方法，包括以下步骤：
步骤1：基于单轴旋转调制捷联惯导***的误差方程建立航天器的状态方程；
步骤2：利用星敏感器获得折射星星点像素坐标量测量，根据量测量建立基于星点像素坐标的量测方程；
步骤3：由于步骤1得出的状态方程和步骤2的量测方程都为非线性，所以使用UKF滤波估计航天器的位置、速度和姿态，并修正航天器的姿态误差和惯性器件误差。
具体步骤如下：
1、以捷联惯导***的误差方程作为***状态方程：
单轴旋转调制式捷联惯导***的捷联解算流程与捷联惯导相同，单轴旋转调制式捷联惯导***的误差方程也与传统捷联惯导***一致。根据组合滤波原理，选取捷联惯性导航***的误差方程作为***的状态方程，用
Figure BDA0002835660670000021
表示状态量，根据惯性导航原理，其状态方程为：
Figure BDA0002835660670000031
式中φ＝[φE φN φU]T是姿态误差，
Figure BDA0002835660670000032
是速度误差，
Figure BDA0002835660670000033
是地理坐标系下各轴向的速度误差；δrn＝[δL δλ δH]T是位置误差，δL是经度误差δλ是纬度误差δH是高度误差，
Figure BDA0002835660670000034
是加速度计的输出在n系中的投影，
Figure BDA0002835660670000035
是加速度计测得的比力，
Figure BDA0002835660670000036
是IMU的姿态矩阵，通过姿态矩阵
Figure BDA0002835660670000037
和旋转矩阵
Figure BDA0002835660670000038
的乘积得到
Figure BDA0002835660670000039
式中
Figure BDA00028356606700000310
θ，ψ分别是俯仰角、横滚角和航向角，
Figure BDA00028356606700000311
由旋转角度决定，捷联惯导***是以zb轴为轴的单轴***，转台的转速用w表示，
Figure BDA00028356606700000312
表示为：
Figure BDA00028356606700000313
t是时间间隔，
Figure BDA00028356606700000314
是n系下wie的误差，
Figure BDA00028356606700000315
表示航天器在n系下的地球自转角速率，
Figure BDA00028356606700000316
是n系下wen的误差，
Figure BDA00028356606700000317
是n系相对e系的旋转角速率在n系中的表示，
Figure BDA00028356606700000318
是
Figure BDA00028356606700000319
的误差，
Figure BDA00028356606700000320
为n系相对惯性空间的转动角速度矢量在n系下的投影，
Figure BDA0002835660670000041
Figure BDA0002835660670000042
表示在p系中陀螺仪的漂移，
Figure BDA0002835660670000043
表示在p系中加速度计的偏置。
2、折射星点像素坐标量测量的获取
最原始的星光折射信息来自于星敏感器拍摄的折射星图，从中可以提取出折射星的像素坐标。在星敏感器观测的折射星图中，折射星分布在星图中靠近地球的部分，非折射星分布在星图中远离地球的部分。由于非折射星相互之间的几何位置未发生改变，因此通过传统的星图识别方法和质心提取方法，可以识别出星图中部分或者所有的非折射星，并得到从星敏坐标系到惯性坐标系的姿态矩阵
Figure BDA00028356606700000411
其中，上标i代表惯性坐标系，下标c代表星敏感器坐标系。
通过传统的星图识别方法和质心提取方法，星图中所有的非折射星可以识别出来，其像素坐标记为(u1,v1)，(u2,v2)，…，(un,vn)，(n为识别出的非折射星的个数)。通过星敏感器成像原理可以得到这些星点在星敏感器坐标系中的三维坐标
Figure BDA0002835660670000044
Figure BDA0002835660670000045
式中，f是焦距，
Figure BDA0002835660670000046
Fov是指视场大小，Nx和Ny是像素数。
根据星图识别方法，得到所有非折射星在惯性系中的坐标
Figure BDA0002835660670000047
Figure BDA0002835660670000048
假定从星敏感器坐标系到惯性坐标系的转换矩阵为
Figure BDA0002835660670000049
可以得到如下关系：
Figure BDA00028356606700000410
Figure BDA0002835660670000051
通过姿态矩阵
Figure BDA0002835660670000059
和星敏感器光轴的指向，结合标准导航星表，可以得到星敏感器所拍摄的视场内所有星未折射前的模拟非折射星图，逐个计算折射星图中每颗星与模拟非折射星图中全部星的欧式距离，若该距离的最小值大于一定阈值(根据模拟星图的位置精度和折射星的识别精度设定)，则认为其为折射星；反之，则为非折射星。将折射星折射前、后的星点坐标记做(ui,vi)和(uri,vri)(i＝1,2,…,nr,nr为折射星的个数)，利用其折射前的星点位置，通过对非折射模拟星图的识别，还可以得到其赤经、赤纬。
3、折射星像素坐标量测方程的建立
折射星像素坐标的量测模型的建立主要包括3个重要步骤：星光折射角的计算，惯性坐标系下折射星矢量的计算，以及星敏感器坐标系下的折射星矢量和折射星像素坐标的计算。
①利用航天器位置估计值
Figure BDA0002835660670000052
和星光矢量S解算出星光折射角的估计值
Figure BDA0002835660670000053
根据星光大气折射模型可以得到以下公式：
Figure BDA0002835660670000054
通过星光折射几何关系，折射视高度还可以表示为：
Figure BDA0002835660670000055
式中r＝[x y z]为状态量的航天器位置矢量，r为航天器位置矢量r的长度，u＝|r·S|＝r cosα，S为恒星单位矢量。Re为地球半径，α为由几何关系产生的一个极小量，可忽略。
结合公式(7)和(8)，得到：
Figure BDA0002835660670000056
通过解算公式(9)就可以得到星光折射角的估计值
Figure BDA0002835660670000057
②计算惯性坐标系下的折射星矢量估计值
Figure BDA0002835660670000058
Figure BDA0002835660670000061
式中：C为旋转矩阵，表达式为：
Figure BDA0002835660670000062
q1,q2,q3,q4为四元数，表达式分别为：
Figure BDA0002835660670000063
其中：
Figure BDA0002835660670000064
③计算星敏感器坐标下的折射星矢量
Figure BDA0002835660670000065
和折射星像素坐标的估计值
Figure BDA0002835660670000066
折射星矢量
Figure BDA0002835660670000067
在星敏感器坐标系中的折射矢量估计
Figure BDA0002835660670000068
为：
Figure BDA0002835660670000069
折射星像素坐标的估计值
Figure BDA00028356606700000610
为：
Figure BDA00028356606700000611
则基于折射星像素坐标或折射星矢量的量测模型可以简化为：
Figure BDA00028356606700000612
式中：v3为折射星像素坐标量测噪声。
4、进行UKF滤波获得航天器的位置速度估计
量测方程(15)是非线性的，***采用无迹卡尔曼滤波(Unscented KalmanFilter，UKF)来对非线性***的数据进行滤波解算，估计航天器的位置、速度和姿态，并修正航天器的姿态误差和惯性器件误差。
它的主要步骤为：在
Figure BDA00028356606700000613
附近选取一系列样本点，使这些样本点的均值和协方差分别为
Figure BDA00028356606700000614
和P(k|k)，设
Figure BDA00028356606700000615
为n×1向量，则可以推导出2n+1个样本点和它对应的权重：
Figure BDA0002835660670000071
Figure BDA0002835660670000072
Figure BDA0002835660670000073
式中，n是状态量X的维数，τ是刻度参数，当状态量噪声服从高斯分布时，通常选取n+τ＝3；
Figure BDA0002835660670000074
表示
Figure BDA0002835660670000075
的第i维列向量，wi是权值，代表的是第i个Sigma点。标准UKF算法如下。
①初始化
Figure BDA0002835660670000076
②计算采样点
Figure BDA0002835660670000077
③时间更新
χi,k|k-1＝f(χi,k-1) (19)
Figure BDA0002835660670000078
Figure BDA0002835660670000079
zi,k|k-1＝h(χk|k-1,k) (22)
Figure BDA00028356606700000710
④量测更新
Figure BDA00028356606700000711
Figure BDA00028356606700000712
Kk＝Pxy,kPyy,k -1 (26)
Figure BDA0002835660670000081
Pk＝Pk|k-1-KkPyy,kKk T (28)
对于线性***来说，UKF与EKF滤波精度基本相同；但对于非线性较强的***来说，UKF性能要优异许多。因此本发明选择UKF方法应用于之后的航天器自主导航***。
本发明与现有技术相比的优点在于：(1)充分利用两种导航信息，实现对航天器的高精度自主导航。(2)折射视高度和星光折射角是星光折射自主导航***中常用的两种量测量，但是无论是折射视高度还是星光折射角作为星光折射量测量，本质上都只能反映折射角的大小这一方面的折射信息，星光折射方向作为与航天器位置矢量相关的折射信息，对于提高航天器自主导航精度有重要影响，所以本文用星点像素坐标作为星光折射量测量，既可以提供折射角的大小，也可以提供星光折射的方向，这是提高星光折射自主导航精度的关键。(3)由于惯性测量元件中的常值偏差，导航误差随着时间积累，最终会导致导航误差发散，旋转调制式惯性导航***能够将短期内相对固定不变的惯性器件常值漂移误差调制成周期性变化的漂移误差，在惯性器件精度水平已定的情况下，采用旋转调制技术，能够有效地提高惯性导航***的导航精度。
附图说明
图1为本发明中航天器星光折射/单轴旋转调制惯性组合导航方法流程图；
图2为本发明中单轴旋转调制式惯性导航原理示意图；
图3为本发明中单轴旋转调制结构图。
具体实施方式
下面结合附图及实施例对本发明进行详细说明。
图1给出了航天器星光折射/单轴旋转调制惯性组合导航方法***流程图。具体实施过程：
1、以捷联惯导***的误差方程作为***状态方程：
单轴旋转调制式捷联惯导***的捷联解算流程与捷联惯导相同，单轴旋转调制式捷联惯导***的误差方程也与传统捷联惯导***一致。根据组合滤波原理，选取捷联惯性导航***的误差方程作为***的状态方程，用
Figure BDA0002835660670000091
表示状态量，根据惯性导航原理，其状态方程为：
Figure BDA0002835660670000092
式中φ＝[φE φN φU]T是姿态误差，
Figure BDA0002835660670000093
是速度误差，
Figure BDA0002835660670000094
是地理坐标系下各轴向的速度误差；δrn＝[δL δλ δH]T是位置误差，δL是经度误差δλ是纬度误差δH是高度误差，
Figure BDA0002835660670000095
是加速度计的输出在n系中的投影，
Figure BDA0002835660670000096
是加速度计测得的比力，
Figure BDA0002835660670000097
是IMU的姿态矩阵，通过姿态矩阵
Figure BDA0002835660670000098
和旋转矩阵
Figure BDA0002835660670000099
的乘积得到
Figure BDA00028356606700000910
式中
Figure BDA00028356606700000911
θ，ψ分别是俯仰角、横滚角和航向角，
Figure BDA00028356606700000912
由旋转角度决定，捷联惯导***是以zb轴为轴的单轴***，转台的转速用w表示，
Figure BDA00028356606700000913
表示为：
Figure BDA00028356606700000914
t是时间间隔，
Figure BDA00028356606700000915
是n系下wie的误差，
Figure BDA00028356606700000916
表示航天器在n系下的地球自转角速率，
Figure BDA00028356606700000917
是n系下wen的误差，
Figure BDA0002835660670000101
是n系相对e系的旋转角速率在n系中的表示，
Figure BDA0002835660670000102
是
Figure BDA0002835660670000103
的误差，
Figure BDA0002835660670000104
为n系相对惯性空间的转动角速度矢量在n系下的投影，
Figure BDA0002835660670000105
Figure BDA0002835660670000106
表示在p系中陀螺仪的漂移，
Figure BDA0002835660670000107
表示在p系中加速度计的偏置。
2、折射星像素坐标量测量的获取
最原始的星光折射信息来自于星敏感器拍摄的折射星图，从中可以提取出折射星的像素坐标。在星敏感器观测的折射星图中，折射星分布在星图中靠近地球的部分，非折射星分布在星图中远离地球的部分。由于非折射星相互之间的几何位置未发生改变，因此通过传统的星图识别方法和质心提取方法，可以识别出星图中部分或者所有的非折射星，并得到从星敏坐标系到惯性坐标系的姿态矩阵
Figure BDA00028356606700001012
其中，上标i代表惯性坐标系，下标c代表星敏感器坐标系。
通过传统的星图识别方法和质心提取方法，星图中所有的非折射星可以识别出来，其像素坐标记为(u1,v1)，(u2,v2)，…，(un,vn)，(n为识别出的非折射星的个数)。通过星敏感器成像原理可以得到这些星点在星敏感器坐标系中的三维坐标
Figure BDA0002835660670000108
Figure BDA0002835660670000109
式中，f是焦距，
Figure BDA00028356606700001010
Fov是指视场大小，Nx和Ny是像素数。
根据星图识别方法，可以得到所有非折射星在惯性系中的坐标
Figure BDA00028356606700001011
Figure BDA0002835660670000111
假定从星敏感器坐标系到惯性坐标系的转换矩阵为
Figure BDA0002835660670000112
可以得到如下关系
Figure BDA0002835660670000113
Figure BDA0002835660670000114
通过姿态矩阵
Figure BDA0002835660670000115
和星敏感器光轴的指向，结合标准导航星表，可以得到星敏感器所拍摄的视场内所有星未折射前的模拟非折射星图，逐个计算折射星图中每颗星与模拟非折射星图中全部星的欧式距离，若该距离的最小值大于一定阈值(根据模拟星图的位置精度和折射星的识别精度设定)，则认为其为折射星；反之，则为非折射星。将折射星折射前、后的星点坐标记做(ui,vi)和(uri,vri)(i＝1,2,…,nr,nr为折射星的个数)，利用其折射前的星点位置，通过对非折射模拟星图的识别，还可以得到其赤经、赤纬。
3、折射星像素坐标量测模型的建立
折射星像素坐标的量测模型的建立主要包括3个重要步骤：星光折射角的计算，惯性坐标系下折射星矢量的计算，以及星敏感器坐标系下的折射星矢量和折射星像素坐标的计算。
①利用航天器位置估计值
Figure BDA0002835660670000118
和星光矢量S解算出星光折射角的估计值
Figure BDA0002835660670000119
根据星光大气折射模型可以得到以下公式：
Figure BDA0002835660670000116
通过星光折射几何关系，折射视高度还可以表示为：
Figure BDA0002835660670000117
式中r＝[x y z]为状态量的航天器位置矢量，r为航天器位置矢量r的长度，u＝|r·S|＝r cosα，S为恒星单位矢量。Re为地球半径，α为由几何关系产生的一个极小量，可忽略。
结合公式(7)和(8)，可以得到：
Figure BDA0002835660670000121
通过解算公式(9)就可以得到星光折射角的估计值
Figure BDA0002835660670000122
②计算惯性坐标系下的折射星矢量估计值
Figure BDA0002835660670000123
Figure BDA0002835660670000124
式中：C为旋转矩阵，表达式为：
Figure BDA0002835660670000125
q1,q2,q3,q4为四元数，表达式分别为：
Figure BDA0002835660670000126
其中：
Figure BDA0002835660670000127
③计算星敏感器坐标下的折射星矢量
Figure BDA0002835660670000128
和折射星像素坐标的估计值
Figure BDA0002835660670000129
折射星矢量
Figure BDA00028356606700001210
在星敏感器坐标系中的折射矢量估计
Figure BDA00028356606700001211
为：
Figure BDA00028356606700001212
折射星像素坐标的估计值
Figure BDA00028356606700001213
为：
Figure BDA00028356606700001214
则基于折射星像素坐标或折射星矢量的量测模型可以简化为：
Figure BDA00028356606700001215
式中：v3为折射星像素坐标量测噪声。
4、进行UKF滤波获得航天器的位置速度估计
量测方程(15)是非线性的，***采用无迹卡尔曼滤波(Unscented KalmanFilter，UKF)来对非线性***的数据进行滤波解算，估计航天器的位置、速度和姿态，并修正航天器的姿态误差和惯性器件误差。
它的主要步骤为：在
Figure BDA0002835660670000131
附近选取一系列样本点，使这些样本点的均值和协方差分别为
Figure BDA0002835660670000132
和P(k|k)，设
Figure BDA0002835660670000133
为n×1向量，则可以推导出2n+1个样本点和它对应的权重：
Figure BDA0002835660670000134
Figure BDA0002835660670000135
Figure BDA0002835660670000136
式中，n是状态量X的维数，τ是刻度参数，当状态量噪声服从高斯分布时，通常选取n+τ＝3；
Figure BDA0002835660670000137
表示
Figure BDA0002835660670000138
的第i维列向量，wi是权值，代表的是第i个Sigma点。标准UKF算法如下。
①初始化
Figure BDA0002835660670000139
②计算采样点
Figure BDA00028356606700001310
③时间更新
χi,k|k-1＝f(χi,k-1) (19)
Figure BDA00028356606700001311
Figure BDA00028356606700001312
zi,k|k-1＝h(χk|k-1,k) (22)
Figure BDA00028356606700001313
④量测更新
Figure BDA0002835660670000141
Figure BDA0002835660670000142
Kk＝Pxy,kPyy,k -1 (26)
Figure BDA0002835660670000143
Pk＝Pk|k-1-KkPyy,kKk T (28)
对于线性***来说，UKF与EKF滤波精度基本相同；但对于非线性较强的***来说，UKF性能要优异许多。因此本文选择UKF方法应用于之后的航天器自主导航***。
图2给出了单轴旋转调制式惯性导航原理示意图。具体实施过程：
旋转调制式惯性导航***中，IMU置于旋转机构上进行旋转，因此需要把IMU的输出转换到航天器本体坐标系中，然后再转换到导航坐标系中进行捷联解算。利用陀螺仪和加速度计输出的角速度和加速度可以解算出航天器的位置、速度和姿态。具体解算过程如图2所示。
旋转调制技术对惯性导航***来说是一种误差补偿技术，其过程可以表示为：初始时刻航天器处于静止状态，即航天器的位置固定不变，在航天器内部，转位机构控制IMU绕竖直方向以一定的角速度开始进行转动，转速已知，则IMU坐标系和航天器本体坐标系之间的转换矩阵为：
Figure BDA0002835660670000144
式中，w为旋转轴的旋转角速度，t为转动时间。
所有的捷联解算都在导航坐标系中进行，因此需要把相关的量测量转换到导航系中，为了使旋转调制的效果简单易懂，设置航天器本体坐标系和导航坐标系重合，也就是
Figure BDA0002835660670000145
为单位矩阵，
Figure BDA0002835660670000146
则在t时刻惯性测量元件的常值偏差在导航坐标系可以表示为：
Figure BDA0002835660670000151
Figure BDA0002835660670000152
式中εn＝[εE εN εU]T和
Figure BDA0002835660670000153
分别代表导航坐标系中和IMU坐标系陀螺仪的常值漂移。
Figure BDA0002835660670000154
和
Figure BDA0002835660670000155
分别代表导航坐标系中和IMU坐标系中加速度计的常值偏置。
由公式(30)和(31)可见，在水平方向，惯性测量元件的常值偏差被调制成具有正弦或者余弦形式的变化曲线，因此在一个完整的转动周期，它们的均值是零，不会对***造成额外的偏差，不会影响组合导航***的精度。然而沿着旋转轴方向上的陀螺仪和加速度计的常值偏差没有变化，它会引起捷联惯导***的定位误差随着时间的累积而增加，综上所述，单轴旋转技术仅能调制垂直旋转轴方向的惯性器件常值偏差，而与旋转轴平行方向的常值偏差仍然按照原来的方式传播。为了使IMU常值偏差在三个方向上均有所调制，可采用敏感轴与旋转轴非重合的转位方案。
如图3所示，图中o-xbybzb代表本体坐标系，o-xpypzp代表IMU坐标系，θbp代表IMU和航天器之间的安装角。在单轴旋转调制过程中，旋转调制技术对陀螺仪和加速度计的偏差具有相同的效果，因此我们以陀螺仪的调制结果为例来进行分析，在起始时刻航天器的导航系和本体系重合，陀螺仪绕着zb轴以恒定的角速度进行连续正反转旋转，旋转矩阵可以表示为：
Figure BDA0002835660670000156
陀螺仪的漂移可以表示为：
Figure BDA0002835660670000161
Figure BDA0002835660670000162
式中，
Figure BDA0002835660670000163
代表一个正反转周期，从式(34)可知，当陀螺仪在x和z方向零偏相等时，且θbp为45°时，陀螺仪在三轴方向上的常值漂移在理论上可以被调制成零。

Claims (5)

1.一种航天器星光折射和单轴旋转调制惯性组合导航方法，其特征在于，包括以下步骤：
步骤1：基于单轴旋转调制捷联惯导***的误差方程建立航天器的状态方程；
步骤2：利用星敏感器获得折射星星点像素坐标量测量，根据量测量建立基于星点像素坐标的量测方程；
步骤3：由于步骤1得出的状态方程和步骤2的量测方程都为非线性，所以使用UKF滤波估计航天器的位置、速度和姿态，并修正航天器的姿态误差和惯性器件误差。
2.根据权利要求1所述的航天器星光折射和单轴旋转调制惯性组合导航方法，其特征在于：所述步骤1中的状态方程如下：
采用
Figure FDA0002835660660000011
表示状态量，状态方程为：
Figure FDA0002835660660000012
式中φ＝[φE φN φU]T是姿态误差，
Figure FDA0002835660660000013
是速度误差，
Figure FDA0002835660660000014
是地理坐标系下各轴向的速度误差；δrn＝[δL δλ δH]T是位置误差，δL是经度误差δλ是纬度误差δH是高度误差，
Figure FDA0002835660660000015
是加速度计的输出在n系中的投影，
Figure FDA0002835660660000016
是加速度计测得的比力，
Figure FDA0002835660660000017
是IMU的姿态矩阵，通过姿态矩阵
Figure FDA0002835660660000018
和旋转矩阵
Figure FDA0002835660660000019
的乘积得到
Figure FDA0002835660660000021
式中
Figure FDA0002835660660000022
θ,ψ分别是俯仰角、横滚角和航向角，
Figure FDA0002835660660000023
由旋转角度决定，捷联惯导***是以zb轴为轴的单轴***，转台的转速用w表示，
Figure FDA0002835660660000024
表示为：
Figure FDA0002835660660000025
t是时间间隔，
Figure FDA0002835660660000026
是n系下wie的误差，
Figure FDA0002835660660000027
表示航天器在n系下的地球自转角速率，
Figure FDA0002835660660000028
是n系下wen的误差，
Figure FDA0002835660660000029
是n系相对e系的旋转角速率在n系中的表示，
Figure FDA00028356606600000210
是
Figure FDA00028356606600000211
的误差，
Figure FDA00028356606600000212
为n系相对惯性空间的转动角速度矢量在n系下的投影，
Figure FDA00028356606600000213
Figure FDA00028356606600000214
表示在p系中陀螺仪的漂移，
Figure FDA00028356606600000215
表示在p系中加速度计的偏置。
3.根据权利要求1所述的航天器星光折射和单轴旋转调制惯性组合导航方法，其特征在于：步骤2中，利用星敏感器获得星点像素坐标量测量的具体如下：
选用视场为10°×10°的星敏感器捕获折射星和非折射星，最原始的星光折射信息来自于星敏感器拍摄的折射星图，从中提取出折射星的像素坐标；在得到的星图中，折射星分布在星图中靠近地球的部分，非折射星分布在星图中远离地球的部分，通过星图匹配和星图识别方法获取星图中的非折射星，得到姿态矩阵
Figure FDA0002835660660000031
其中，上标i代表惯性坐标系，下标c代表星敏感器坐标系；
通过星图识别方法和质心提取方法，星图中所有的非折射星识别出来，其像素坐标记为(u1,v1)，(u2,v2)，…，(un,vn)，n为识别出的非折射星的个数，通过星敏感器成像原理得到这些星点在星敏感器坐标系中的三维坐标
Figure FDA0002835660660000032
Figure FDA0002835660660000033
Figure FDA0002835660660000034
式中，f是焦距，
Figure FDA0002835660660000035
Fov是指视场大小，Nx和Ny是像素数；
根据星图识别方法，得到所有非折射星在惯性系中的坐标
Figure FDA0002835660660000036
Figure FDA0002835660660000037
从星敏感器坐标系到惯性坐标系的转换矩阵为
Figure FDA0002835660660000038
得到如下关系：
Figure FDA0002835660660000039
Figure FDA00028356606600000310
通过姿态矩阵
Figure FDA00028356606600000311
和星敏感器光轴的指向，结合标准导航星表，得到星敏感器所拍摄的视场内所有星未折射前的模拟非折射星图，逐个计算折射星图中每颗星与模拟非折射星图中全部星的欧式距离，若该距离的最小值大于设定阈值(根据模拟星图的位置精度和折射星的识别精度设定)，则认为其为折射星；反之，则为非折射星；将折射星折射前、后的星点坐标记做(ui,vi)和(uri,vri)，i＝1,2,…,nr,nr为折射星的个数。
4.根据权利要求1所述的航天器星光折射和单轴旋转调制惯性组合导航方法，其特征在于：所述步骤2中，折射星星点像素坐标量测方程的建立如下：
折射星星点像素坐标的量测方程的建立包括3个重要步骤：星光折射角的计算，惯性坐标系下折射星矢量的计算，以及星敏感器坐标系下的折射星矢量和折射星像素坐标的计算；
①利用航天器位置估计值
Figure FDA0002835660660000041
和星光矢量S解算出星光折射角的估计值
Figure FDA0002835660660000042
根据星光大气折射模型得到以下公式：
Figure FDA0002835660660000043
通过星光折射几何关系，折射视高度还表示为：
Figure FDA0002835660660000044
式中r＝[x y z]为状态量的航天器位置矢量，r为航天器位置矢量r的长度，u＝|r·S|＝r cosα，S为恒星单位矢量，Re为地球半径，α为由几何关系产生的一个极小量，可忽略；
结合公式(7)和(8)，得到：
Figure FDA0002835660660000045
通过解算公式(9)则得到星光折射角的估计值
Figure FDA0002835660660000046
②计算惯性坐标系下的折射星矢量估计值
Figure FDA0002835660660000047
Figure FDA0002835660660000048
式中：C为旋转矩阵，表达式为：
Figure FDA0002835660660000049
q1,q2,q3,q4为四元数，表达式分别为：
Figure FDA00028356606600000410
其中：
Figure FDA0002835660660000051
③计算星敏感器坐标下的折射星矢量
Figure FDA0002835660660000052
和折射星像素坐标的估计值
Figure FDA0002835660660000053
折射星矢量
Figure FDA0002835660660000054
在星敏感器坐标系中的折射矢量估计
Figure FDA0002835660660000055
为：
Figure FDA0002835660660000056
折射星像素坐标的估计值
Figure FDA0002835660660000057
为：
Figure FDA0002835660660000058
则基于折射星像素坐标或折射星矢量的量测模型简化为：
Figure FDA0002835660660000059
式中：v3为折射星像素坐标量测噪声。
5.根据权利要求1所述的航天器星光折射和单轴旋转调制惯性组合导航方法，其特征在于：所述步骤3：使用UKF滤波估计航天器的位置、速度和姿态估计的过程如下：
在
Figure FDA00028356606600000510
附近选取一系列样本点，使这些样本点的均值和协方差分别为
Figure FDA00028356606600000511
和P(k|k)，设
Figure FDA00028356606600000512
为n×1向量，则得出2n+1个样本点和它对应的权重：
Figure FDA00028356606600000513
式中，n是状态量X的维数，τ是刻度参数，当状态量噪声服从高斯分布时，通常选取n+τ＝3；
Figure FDA00028356606600000514
表示
Figure FDA00028356606600000515
的第i维列向量，wi是权值，代表的是第i个Sigma点，标准UKF算法如下：
①初始化
Figure FDA00028356606600000516
②计算采样点
Figure FDA0002835660660000061
③时间更新
χi,k|k-1＝f(χi,k-1) (19)
Figure FDA0002835660660000062
Figure FDA0002835660660000063
zi,k|k-1＝h(χk|k-1,k) (22)
Figure FDA0002835660660000064
④量测更新
Figure FDA0002835660660000065
Figure FDA0002835660660000066
Kk＝Pxy,kPyy,k -1 (26)
Figure FDA0002835660660000067
Pk＝Pk|k-1-KkPyy,kKk T (28)。
CN202011469193.7A 2020-12-14 2020-12-14 一种航天器星光折射和单轴旋转调制惯性组合导航方法 Active CN112880669B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011469193.7A CN112880669B (zh) 2020-12-14 2020-12-14 一种航天器星光折射和单轴旋转调制惯性组合导航方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011469193.7A CN112880669B (zh) 2020-12-14 2020-12-14 一种航天器星光折射和单轴旋转调制惯性组合导航方法

Publications (2)

Publication Number Publication Date
CN112880669A true CN112880669A (zh) 2021-06-01
CN112880669B CN112880669B (zh) 2024-01-16

Family

ID=76043334

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011469193.7A Active CN112880669B (zh) 2020-12-14 2020-12-14 一种航天器星光折射和单轴旋转调制惯性组合导航方法

Country Status (1)

Country Link
CN (1) CN112880669B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113551668A (zh) * 2021-07-21 2021-10-26 北京航空航天大学 一种航天器惯性/恒星星光矢量/星光折射组合导航方法
CN113916717A (zh) * 2021-10-11 2022-01-11 北京航空航天大学 一种基于低轨航天器掩星的平流层大气密度反演方法
CN113916217A (zh) * 2021-10-11 2022-01-11 北京航空航天大学 一种基于分区平流层大气折射模型的星光定位方法
CN115638796A (zh) * 2022-09-19 2023-01-24 北京控制工程研究所 一种基于折射星/非折射星信息融合与预测的快速星图识别方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103913169A (zh) * 2014-03-12 2014-07-09 哈尔滨工程大学 一种飞行器的捷联惯性/星光折射组合导航方法
CN103994763A (zh) * 2014-05-21 2014-08-20 北京航空航天大学 一种火星车的sins/cns深组合导航***及其实现方法
CN108871326A (zh) * 2018-07-09 2018-11-23 北京航空航天大学 一种单轴旋转调制惯性-天文深组合导航方法
CN111060097A (zh) * 2020-01-15 2020-04-24 东南大学 一种提高位置误差估计精度的惯性/天文组合导航方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103913169A (zh) * 2014-03-12 2014-07-09 哈尔滨工程大学 一种飞行器的捷联惯性/星光折射组合导航方法
CN103994763A (zh) * 2014-05-21 2014-08-20 北京航空航天大学 一种火星车的sins/cns深组合导航***及其实现方法
CN108871326A (zh) * 2018-07-09 2018-11-23 北京航空航天大学 一种单轴旋转调制惯性-天文深组合导航方法
CN111060097A (zh) * 2020-01-15 2020-04-24 东南大学 一种提高位置误差估计精度的惯性/天文组合导航方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
YANG SHUJIE ET AL.: "Stellar Refraction-Based SINS/CNS Integrated Navigation System for Aerospace Vehicles", 《JOURNAL OF AEROSPACE ENGINEERING》, vol. 29, no. 2, pages 1 - 11 *
宁晓琳等: "基于折射方向矢量的地球卫星星光折射导航新方法", 《飞控与探测》, vol. 3, no. 2, pages 8 - 16 *
宋伟等: "基于UKF的航天器多普勒/天文组合导航方法研究", 《载人航天》, no. 3, pages 20 - 28 *
杨淑洁等: "一种航天飞行器的INS/CNS自主导航方案", 《中国惯性技术学报》, vol. 22, no. 6, pages 728 - 733 *
王融;熊智;刘建业;: "星光折射定位辅助的惯性/星光全组合导航算法", 传感器与微***, no. 12, pages 140 - 142 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113551668A (zh) * 2021-07-21 2021-10-26 北京航空航天大学 一种航天器惯性/恒星星光矢量/星光折射组合导航方法
CN113551668B (zh) * 2021-07-21 2024-05-28 北京航空航天大学 一种航天器惯性/恒星星光矢量/星光折射组合导航方法
CN113916717A (zh) * 2021-10-11 2022-01-11 北京航空航天大学 一种基于低轨航天器掩星的平流层大气密度反演方法
CN113916217A (zh) * 2021-10-11 2022-01-11 北京航空航天大学 一种基于分区平流层大气折射模型的星光定位方法
CN113916217B (zh) * 2021-10-11 2023-06-16 北京航空航天大学 一种基于分区平流层大气折射模型的星光定位方法
CN113916717B (zh) * 2021-10-11 2023-08-11 北京航空航天大学 一种基于低轨航天器掩星的平流层大气密度反演方法
CN115638796A (zh) * 2022-09-19 2023-01-24 北京控制工程研究所 一种基于折射星/非折射星信息融合与预测的快速星图识别方法

Also Published As

Publication number Publication date
CN112880669B (zh) 2024-01-16

Similar Documents

Publication Publication Date Title
CN107655476B (zh) 基于多信息融合补偿的行人高精度足部导航方法
CN111947652B (zh) 一种适用于月球着陆器的惯性/视觉/天文/激光测距组合导航方法
CN110501024B (zh) 一种车载ins/激光雷达组合导航***的量测误差补偿方法
CN112880669B (zh) 一种航天器星光折射和单轴旋转调制惯性组合导航方法
CN112629538B (zh) 基于融合互补滤波和卡尔曼滤波的舰船水平姿态测量方法
CN111156994B (zh) 一种基于mems惯性组件的ins/dr&gnss松组合导航方法
CN109596018B (zh) 基于磁测滚转角速率信息的旋转弹飞行姿态高精度估计方法
CN101793523B (zh) 一种组合导航和光电探测一体化***
CN113551668B (zh) 一种航天器惯性/恒星星光矢量/星光折射组合导航方法
CN110954102B (zh) 用于机器人定位的磁力计辅助惯性导航***及方法
CN109596144B (zh) Gnss位置辅助sins行进间初始对准方法
CN111121766B (zh) 一种基于星光矢量的天文与惯性组合导航方法
CN111024070A (zh) 一种基于航向自观测的惯性足绑式行人定位方法
CN113063429B (zh) 一种自适应车载组合导航定位方法
CN112504275B (zh) 一种基于级联卡尔曼滤波算法的水面舰船水平姿态测量方法
CN111189442B (zh) 基于cepf的无人机多源导航信息状态预测方法
CN103076015A (zh) 一种基于全面最优校正的sins/cns组合导航***及其导航方法
CN111102993A (zh) 一种旋转调制型捷联惯导***晃动基座初始对准方法
Fu et al. Autonomous in-motion alignment for land vehicle strapdown inertial navigation system without the aid of external sensors
CN110296719B (zh) 一种在轨标定方法
CN103604428A (zh) 基于高精度水平基准的星敏感器定位方法
CN106352897B (zh) 一种基于单目视觉传感器的硅mems陀螺误差估计与校正方法
CN116105730A (zh) 基于合作目标卫星甚短弧观测的仅测角光学组合导航方法
CN112902956A (zh) 一种手持式gnss/mems-ins接收机航向初值获取方法、电子设备、存储介质
CN104501809B (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