CN103941593B

CN103941593B - 低轨卫星姿态仿真方法

Info

Publication number: CN103941593B
Application number: CN201410134688.2A
Authority: CN
Inventors: 唐新明; 张过; 黄文超; 方辰; 高小明; 岳庆兴
Original assignee: SATELLITE SURVEYING AND MAPPING APPLICATION CENTER NASG
Current assignee: Ministry Of Natural Resources Land Satellite Remote Sensing Application Center
Priority date: 2014-04-04
Filing date: 2014-04-04
Publication date: 2018-02-02
Anticipated expiration: 2034-04-04
Also published as: CN103941593A

Abstract

本发明属于摄影测量与计算机仿真领域，特别涉及低轨卫星的姿态仿真方法。本发明目的在于解决现有技术不足，提出了一种全新的低轨卫星姿态仿真方法，填补了国内在这方面的空白，为低轨卫星的预研工作提供了分析的依据。本发明的技术方案是：对低轨卫星姿态过程进行分析，分析其标称状态下的因素以及影响其测量精度的误差因素，然后通过对这些因素进行建模，在计算机上对低轨卫星的姿态过程进行仿真。

Description

低轨卫星姿态仿真方法

技术领域

本发明属于摄影测量与计算机仿真领域，特别涉及低轨卫星的姿态仿真方法。

背景技术

在40-50年代，由于计算机技术的限制，只有依靠物理仿真，在美国亚利桑那大学光学中心建立了世界上第一个航空航天遥感物理仿真系统。在地面实验室里利用人造光源提供各种辐亮度和各个光谱谱段的照明条件，布置了不同背景下的各种大小尺寸靶标和军用目标模型(包括飞机、坦克、火炮等)，可以模拟卫星在轨飞行情况下的环境条件以及目标的运动等，采用可控制位置和运动模式的相机对目标按照预定的程序进行照相，以验证卫星的设计参数和成像质量。

在60年代，美国发射了多颗地球环境探测卫星，获得了大量地表、大气和地球环境的数据，这些数据为仿真实验室提供了接近真实的模型。从60年代到90年代，美国多次发射地球地理环境探测、校验和测绘卫星，用于监视和补充数据资料，修正数学模型。

80年代末期，ES公司首先在美国GE公司13个部门中的八个部门中应用了该公司的集成设计软件iSIGHT。1995年，美国NASA资助的LaRC(Landley Research Center)公布了PATCOD集成设计软件平台。美国NASA所属JPL实验室的飞行系统测试平台(Flight SystemTested，FST)、Langley研究中心的SPASIM(Spacecrfat Simulation)、俄罗斯能源科学生产联合体(NPO Energiya)的综合仿真测试平台(KMC)以及德国VEGA信息技术公司开发的仿真卫星(Simulating Spacecraft)等是九十年代卫星仿真技术发展的综合反映。这些软件用于航天卫星(重点是对卫星平台等大系统)的设计和仿真。

目前国外计算机仿真技术发展很快，可以比较逼真地仿真出成像链路的特性，取得了一定成果，但是还是不能代替物理仿真。且在查到的仿真软件文献中，多数是应用成果介绍，很少见到详细的软件内容。因此发展国内的全链路仿真算法和系统来指导卫星的预研工作、从而减少物理仿真的成本是很有必要的，本发明就是涉及全链路仿真系统中比较重要的姿态仿真部分。

发明内容

本发明目的在于解决现有技术不足，提出了一种全新的低轨卫星姿态仿真方法，填补了国内在这方面的空白，为低轨卫星的预研工作提供了分析的依据。

本发明的技术方案是：对低轨卫星姿态过程进行分析，分析其标称状态下的因素以及影响其测量精度的误差因素，然后通过对这些因素进行建模，在计算机上对低轨卫星的姿态过程进行仿真。

本发明提供了一种低轨卫星标称姿态仿真方法，其特征在于，所述方法具体包括以下步骤：

步骤1.1，计算三轴欧拉角速度，具体的，根据三轴姿态稳定度的设计指标，按照模拟高斯白噪声的模式模拟一组符合三轴稳定度指标的三轴欧拉姿态角速度；

步骤1.2，计算三轴欧拉角，具体的，根据三轴欧拉姿态角速度和假定景中心时刻的三轴欧拉角，根据如下公式的积分获得三轴欧拉姿态角：

其中α_t-Δt为前一时刻的三轴欧拉姿态角，为前一时刻的三轴欧拉角速度，α_t(t)为当前时刻的三轴欧拉角；

步骤1.3，生成测姿坐标系到轨道坐标系的旋转矩阵，具体的，根据设定的旋转矩阵形成顺序以及三轴欧拉角，生成测姿坐标系到轨道坐标系的旋转矩阵

步骤1.4，生成测姿坐标系到J2000的旋转矩阵，根据以下公式，计算测姿坐标系带J2000坐标系的旋转矩阵以及相应的J2000坐标系下姿态角：

其中：

式中，X_s，Y_s，Z_s和分别为卫星质心在J2000坐标系中的位置和速度；

步骤1.5，叠加星敏测姿的高频可测颤振；颤振公式如下：

ω=2πf

式中，t为采样间隔，f为颤振的频率，A为颤振的幅度，为相位角，y为颤振大小；

步骤1.6，生成本体坐标系到测姿坐标系的旋转矩阵；具体的，根据设定的旋转矩阵形成顺序以及测姿仪器三轴欧拉安装角，形成本体坐标系到测姿坐标系的旋转矩阵

步骤1.7，生成相机坐标系到本体坐标系的旋转矩阵；具体的，根据设定的旋转矩阵形成顺序以及传感器三轴欧拉安装角，生成相机坐标系到本体坐标系的旋转矩阵

步骤1.8，生成相机坐标系到J2000坐标系的旋转矩阵；具体的，根据上述制作的测姿坐标系相对于J2000的旋转矩阵、本体坐标系到测姿坐标系的旋转矩阵、相机坐标系到本体坐标系的旋转矩阵，按如下式子构成相机坐标系相对J2000坐标系的旋转矩阵：

步骤1.9，生成姿态四元数。根据如下公式生成姿态四元组：

上述公式可以根据下面对应关系来计算：

步骤1.10，根据生成的姿态四元数对低轨卫星姿态进行仿真计算；

步骤1.11，在计算机显示器上显示仿真结果。

优选的，在所述步骤1.5中，所述颤振模型是许多个波函数的叠加，如下公式所示：

(其中当y<σ时，y=0)。

优选的，在所述步骤1.5中，在添加颤振的时候将上述测姿坐标系相对于J2000的旋转矩阵按照转序要求变化为三轴欧拉角，在三轴欧拉角上添加颤振，然后再将添加了颤振的三轴欧拉角转化为旋转矩阵。

本发明还提供了一种星敏定姿低轨卫星误差姿态生成方法，其特征在于，所述方法具体包括以下步骤：

步骤2.1，计算三轴欧拉角速度；根据三轴姿态稳定度的设计指标，按照模拟高斯白噪声的模式模拟一组符合三轴稳定度指标的三轴欧拉姿态角速度；

步骤2.2，计算三轴欧拉角；根据三轴欧拉姿态角速度和假定景中心时刻的三轴欧拉角，根据如下公式的积分获得三轴欧拉姿态角：

步骤2.3，生成星敏测姿坐标系到轨道坐标系旋转矩阵；具体的，根据设定的旋转矩阵形成顺序，根据三轴欧拉角，形成测姿坐标系到轨道坐标系的旋转矩阵

步骤2.4，生成测姿坐标系到J2000的旋转矩阵；根据以下公式，计算测姿坐标系到J2000坐标系的旋转矩阵：

其中：

式中，X_s，Y_s，Z_s和为卫星质心在J2000坐标系中的位置和速度；

步骤2.5，叠加星敏测姿的高频可测颤振；颤振公式如下：

ω=2πf

步骤2.6，叠加测姿系统误差和测姿随机误差；根据测姿系统随机误差的设计指标，按照模拟高斯白噪声的模式模拟一组符合测姿系统随机误差指标的三轴欧拉姿态角随机误差，然后再将随机误差及设定的系统误差直接叠加到J2000坐标系下的三轴欧拉角上；

步骤2.7，生成本体坐标系到测姿坐标系的旋转矩阵；具体的，根据设定的旋转矩阵形成顺序，根据测姿仪器三轴欧拉安装角以及安装角的系统误差，形成本体坐标系到测姿坐标系的旋转矩阵

步骤2.8，生成相机坐标系到本体坐标系的旋转矩阵；具体的，根据设定的旋转矩阵形成顺序，根据传感器三轴欧拉安装角以及安装角的系统误差，形成相机坐标系到本体坐标系的旋转矩阵

步骤2.9，生成相机坐标系到.J2000坐标系的旋转矩阵；根据上述制作的测姿坐标系相对于J2000的旋转矩阵、本体坐标系到测姿坐标系的旋转矩阵、相机坐标系到本体坐标系的旋转矩阵，按如下式子构成相机坐标系到J2000坐标系的旋转矩阵：

步骤2.10，叠加星敏测姿的高频不可测颤振以及相机的整体颤振；

步骤2.11，叠加姿态量化误差。

步骤2.12，生成姿态四元数。根据如下公式制作姿态四元组：

上述公式可以根据下面对应关系来计算：

步骤2.13，根据生成的姿态四元数对低轨卫星姿态进行仿真计算；

步骤2.14，在计算机显示器上显示仿真结果。

优选的，在所述步骤2.5中，所述颤振模型是许多个波函数的叠加，如下公式所示：

(其中当y<σ时，y=0)。

优选的，在所述步骤2.5，2.6和2.7中，在添加颤振的时候将上述测姿坐标系相对于J2000的旋转矩阵按照转序要求变化为三轴欧拉角，在三轴欧拉角上添加颤振，然后再将添加了颤振的三轴欧拉角转化为旋转矩阵。

本发明还提供了一种陀螺定姿或者陀螺星敏联合定姿的低轨卫星误差姿态生成方法，所述方法具体包括以下步骤：：

步骤3.1，计算三轴欧拉角速度；根据三轴姿态稳定度的设计指标，按照模拟高斯白噪声的模式模拟一组符合三轴稳定度指标的三轴欧拉姿态角速度；

步骤3.2，计算三轴欧拉角；根据三轴欧拉姿态角速度和假定景中心时刻的三轴欧拉角，根据如下公式的积分获得三轴欧拉姿态角：

其中α_t-Δt为前一时刻的三轴欧拉姿态角，为前一时刻的三轴欧拉角速度，α_t(t)为当前时刻的三轴欧拉角。

步骤3.3，生成陀螺坐标系到轨道坐标系旋转矩阵；具体的，根据设定的旋转矩阵形成顺序，根据三轴欧拉角，形成测姿坐标系到轨道坐标系的旋转矩阵

步骤3.4，生成测姿坐标系到J2000的旋转矩阵，根据以下公式，计算测姿坐标系到J2000坐标系的旋转矩阵：

其中：

式中，X_s，Y_s，Z_s和为卫星质心在J2000坐标系中的位置和速度。

步骤3.5，叠加星敏测姿的高频可测颤振，颤振公式如下：

ω=2πf

步骤3.6，将上述添加了颤振的三轴欧拉角从J2000坐标系转化到轨道坐标系下，转化方法类似步骤3.4；

步骤3.7，叠加测姿系统误差和测姿随机误差；根据测姿系统随机误差的设计指标，按照模拟高斯白噪声的模式模拟一组符合测姿系统随机误差指标的三轴欧拉姿态角随机误差，然后再将随机误差及设定的系统误差直接叠加到轨道坐标系下的三轴欧拉角上。

步骤3.8，生成测姿坐标系相对于J2000的旋转矩阵，根据设定的旋转矩阵形成顺序，根据三轴欧拉角，形成测姿坐标系相对于轨道坐标系的旋转矩阵根据以下公式，计算测姿坐标系相对于J2000坐标系的旋转矩阵：

步骤3.9，生成本体坐标系到测姿坐标系的旋转矩阵；根据设定的旋转矩阵形成顺序，根据测姿仪器三轴欧拉安装角以及安装角的系统误差，形成本体坐标系相对测姿坐标系的旋转矩阵

步骤3.10，生成相机坐标系到本体坐标系的旋转矩阵；根据设定的旋转矩阵形成顺序，根据传感器三轴欧拉安装角以及安装角的系统误差，形成相机坐标系相对本体坐标系的旋转矩阵

步骤3.11，生成相机坐标系到J2000坐标系的旋转矩阵；根据上述制作的测姿坐标系相对于J2000的旋转矩阵、本体坐标系到测姿坐标系的旋转矩阵、相机坐标系到本体坐标系的旋转矩阵，按如下式子构成相机坐标系到J2000坐标系的旋转矩阵：

步骤3.12，叠加星敏测姿的高频不可测颤振以及相机的整体颤振；

步骤3.13，叠加姿态量化误差；

步骤3.14，生成姿态四元数。根据如下公式制作姿态四元组：

上述公式可以根据下面对应关系来计算：

步骤3.15，根据生成的姿态四元数对低轨卫星姿态进行仿真计算；

步骤3.16，在计算机显示器上显示仿真结果。

优选的，在所述步骤3.5中，所述颤振模型是许多个波函数的叠加，如下公式所示：

(其中当y<σ时，y=0)。

优选的，在所述步骤3.5和3.12中，在添加颤振的时候将上述测姿坐标系相对于J2000的旋转矩阵按照转序要求变化为三轴欧拉角，在三轴欧拉角上添加颤振，然后再将添加了颤振的三轴欧拉角转化为旋转矩阵。

本发明提出了一种全新的低轨卫星姿态仿真方法，填补了国内在这方面的空白，为低轨卫星的预研工作提供了分析的依据。

附图说明

图1是本发明实施例一流程图；

图2是本发明实施例二流程图；

图3是本发明实施例三流程图。

具体实施方式

为了更好的理解本发明的技术方案，下面对本发明涉及到几个概念的进行介绍：

1)严密几何定位模型

卫星严密几何成像模型如下所示：

式中：表示卫星在WGS84坐标系下的位置矢量；m为比例系数；R_J20002WGS84为J2000坐标系到WGS84坐标系的变换矩阵；R_orbit2J2000为轨道坐标系到J2000坐标系的变换矩阵；R_star2orbit为测姿坐标系到轨道坐标系的变换矩阵；为本体坐标系到测姿坐标系的变换矩阵，它由测姿系统的安装确定；R_camera2body为传感器坐标系到本体坐标系的变换矩阵，由相机安装确定；为定位设备在本体坐标系下的偏移；为传感器安装到本体坐标系的偏移；为像素点对应的成像内方位元素。所有这些项的误差均会引起最终的定位精度。而在实际的卫星在轨运行中，由于测量设备的精度限制、成像平台的稳定性及成像环境等的变化，存在很多因素引起定位项的误差。如测轨误差会造成R_orbit2J2000、R_J20002WGS84的误差；测姿误差会造成R_star2orbit的误差；设备安装误差会造成 R_camera2body的误差；相机镜头畸变会造成的误差；这些误差项最终均降低了几何定位精度。在综合分析影响几何定位精度的各方面因素后，我们将影响定位的误差源分为：姿态误差、轨道误差、时间同步误差、内方位元素误差、高程误差和控制点数量及其分布影响。其中姿态误差是卫星定位误差一种比较重要的误差源。

2)姿态仿真考虑因素

姿态仿真需要考虑标称姿态和误差姿态，因此在考虑其影响因素的时候，分成两部分考虑。

在标称姿态中不含误差，为卫星在轨运行过程中测姿仪器能够测量得到的，主要有姿态稳定度和高频可测部分。姿态稳定度姿态角长周期变化量，姿态角速度统计量，单位为°／s，利用随机统计原理模拟成像时间段内的姿态角速率，进行叠加。高频可测部分指卫星在颤振过程中，如果颤振振幅能够被测姿仪器测得，则其不认为为误差，添加到标称姿态中。

在误差姿态中需要考虑各类误差，我们分为设备安装误差、测不准误差、测不出误差以及输不出误差。设备安装误差是由于相机、测姿仪器等安装角上天后环境改变而发生变化。测不准误差则是由于测姿仪器的物理限制而造成的测量误差。测不出误差则是由于测姿仪器敏感程度低于敏感点的姿态角无法测出以及高频颤振的不可测部分，高频颤振的不可测部分包括传感器的高频颤振以及测姿设备不可测的高频颤振部分。输不出误差则是由于星上存储设备具有一定存储瓶颈，对姿态测量数据的量化位数有限造成的误差。

3)高斯白噪声的生成

本项目高斯白噪声的模拟是利用随机数生成。目前，常用的随机数发生器大致可以分成两类：一类是软件控制的随机数发生器，由某种算法生成，称其为伪随机数发生器；另一类是硬件随机数发生器，如噪声随机数发生器，能生成真正的随机数据流，是当今最先进的要数量子随机数发生器，但价格十分昂贵。

本项目采用软件控制的伪随机数发生器，一次产生符合试验要求的若干个随机数，如在卫星位置误差对目标定位精度的影响试验中，一次产生n个符合高斯分布的伪随机数。为了使该随机数列服从正态分布，先利用乘同余法产生一组近似符合[0，1]分布的伪随机数，再以此伪随机数列为基础，将其标准化产生一组给定精度的随机数，该组随机数的均值为零，方差为给定测量精度的1±0.01倍，然后将其顺序地加入到以上各参数中即获得已知精度的轨道观测值。

4)姿态机动模型建模

假设姿态初始值为r₀、侧摆角速度为dr、TT表示起始时刻持续时间，r₁为侧摆最终数值。则我们可以构建姿态机动模型，如下所示：

r₁=r₀+dr*TT

其中用|r(t)-r₁|≤0.1°控制上式中的±，为高斯白噪声。在考虑姿态机动情况下，在制作三轴欧拉角的基础上，添加姿态机动的数值。

5)姿态稳定度模型建模

假设姿态稳定度为α_t-Δt为前一时刻的三轴欧拉姿态角，α_t(t)为当前时刻的三轴欧拉角，则我们可以根据三轴欧拉姿态角速度和假定景中心时刻的三轴欧拉角的物理关系建立如下所示的姿态稳定度积分模型：

这里需要说明的一点是，姿态稳定度有一个限值，当累加的姿态超过限值时可以被姿态测量设备检测出来，此时采用的策略是将累加后的值减去2倍的欧拉角速度累加值：

6)颤振模型建模

假设t为采样间隔、f为颤振的频率、A为颤振的幅度、为相位角、y为颤振大小。则我们可以得到单正弦波的颤振模型如公式所示

ω=2πf

在实际过程中，颤振模型应该是许多个波函数的叠加，如下所示：

(其中当y<σ时，y=0)

上式颤振叠加模型，相位角取0，对于采样中多个频率的颤振累加到相应的姿态上。

本发明针对的低轨卫星姿态仿真算法根据仿真要求分为如下三种情况：

(1)标称姿态仿真，根据姿态稳定度指标及颤振模型模拟轨道坐标系和本体坐标系、J2000坐标系和本体坐标系等之间三轴欧拉角，得到模拟的无误差状态下的姿态测量数据。

(2)星敏定姿误差姿态仿真，是根据姿态指向和测量误差指标模拟轨道坐标系和本体坐标系、J2000坐标系和本体坐标系等之间三轴欧拉角，得到模拟的星敏定姿状态下含有误差的姿态测量数据。

(3)陀螺定姿及星敏陀螺联合定姿误差姿态仿真，是根据姿态指向和测量误差指标模拟轨道坐标系和本体坐标系、J2000坐标系和本体坐标系等之间三轴欧拉角，得到模拟的陀螺定姿或星敏陀螺定姿状态下含有误差的姿态测量数据。

需要指出的是，本发明虽然提出了三种低轨卫星姿态仿真算法，但是这三种算法中的低轨卫星的标称和误差姿态仿真，都是基于姿态机动模型、姿态稳定性模型以及姿态颤振模型等建立，并根据各类姿态角之间的旋转关系来求取三轴欧拉角，所以，本发明的所提出的上述三种低轨卫星姿态仿真算法基于同一个发明构思，具备专利法规定的单一性的要求，可以合案申请。

在具体实施本发明所提出的低轨卫星姿态仿真算法时，即可以采用计算机硬件完成，也可以采用计算机软件完成，也可以采用计算机硬件和软件二者来共同完成，本领域技术人员可以根据实际需要来选择实现的方式，本发明对此不做限制，以下结合实施例和附图1～3提供详细说明本发明的技术方案。

实施例一

实施例一提供了一种低轨卫星标称姿态仿真方法，参见附图1，所述方法具体包括以下步骤：

步骤1.1，计算三轴欧拉角速度。根据三轴姿态稳定度的设计指标，按照模拟高斯白噪声的模式模拟一组符合三轴稳定度指标的三轴欧拉姿态角速度；

步骤1.2，计算三轴欧拉角。根据三轴欧拉姿态角速度和假定景中心时刻的三轴欧拉角，根据如下公式的积分获得三轴欧拉姿态角：

这里需要说明的一点是，姿态稳定度有一个限值，当累加的姿态超过限值时可以被姿态测量设备检测出来，此时采用的策略是将累加后的值减去2倍的欧拉角速度累加值，如下式所示：

步骤1.3，生成测姿坐标系到轨道坐标系的旋转矩阵。具体的，根据设定的旋转矩阵形成顺序以及三轴欧拉角，生成测姿坐标系到轨道坐标系的旋转矩阵

步骤1.4，生成测姿坐标系到J2000的旋转矩阵。根据以下公式，计算测姿坐标系带J2000坐标系的旋转矩阵以及相应的J2000坐标系下姿态角：

其中：

式中，X_s，Y_s，Z_s和分别为卫星质心在J2000坐标系中的位置和速度。

步骤1.5，叠加星敏测姿的高频可测颤振。颤振公式如下：

ω=2πf

式中，t为采样间隔，f为颤振的频率，A为颤振的幅度，为相位角，y为颤振大小，仿真过程中这些数值都可以从颤振文件中直接读取。颤振模型是许多个波函数的叠加，如下公式所示：

(其中当y<σ时，y=0)

上式表示本次仿真过程中所采用的星敏颤振叠加的具体方法，相位角取0，对于采样中多个频率的颤振累加到相应的姿态上。设姿态测量设备的测姿灵敏度为σ，当上述累加值超过σ时，说明此颤振可以被姿态测量设备探测出来，需要将此颤振添加到姿态上面。

在添加的时候需要将上述测姿坐标系相对于J2000的旋转矩阵按照转序要求变化为三轴欧拉角，在三轴欧拉角上添加颤振，然后再将添加了颤振的三轴欧拉角转化为旋转矩阵。

步骤1.6，生成本体坐标系到测姿坐标系的旋转矩阵。具体的，根据设定的旋转矩阵形成顺序以及测姿仪器三轴欧拉安装角，形成本体坐标系到测姿坐标系的旋转矩阵

步骤1.7，生成相机坐标系到本体坐标系的旋转矩阵。具体的，根据设定的旋转矩阵形成顺序以及传感器三轴欧拉安装角，生成相机坐标系到本体坐标系的旋转矩阵

步骤1.8，生成相机坐标系到J2000坐标系的旋转矩阵。具体的，根据上述制作的测姿坐标系相对于J2000的旋转矩阵、本体坐标系到测姿坐标系的旋转矩阵、相机坐标系到本体坐标系的旋转矩阵，按如下式子构成相机坐标系相对J2000坐标系的旋转矩阵：

步骤1.9，生成姿态四元数。根据如下公式生成姿态四元组：

上述公式可以根据下面对应关系来计算：

步骤1.11，在计算机显示器上显示仿真结果。

实施例二

实施例二提供了一种星敏定姿低轨卫星误差姿态生成方法，参见附图2，所述方法具体包括以下步骤：

步骤2.1，计算三轴欧拉角速度。根据三轴姿态稳定度的设计指标，按照模拟高斯白噪声的模式模拟一组符合三轴稳定度指标的三轴欧拉姿态角速度。

步骤2.2，计算三轴欧拉角。根据三轴欧拉姿态角速度和假定景中心时刻的三轴欧拉角，根据如下公式的积分获得三轴欧拉姿态角：

步骤2.3，生成星敏测姿坐标系到轨道坐标系旋转矩阵。具体的，根据设定的旋转矩阵形成顺序，根据三轴欧拉角，形成测姿坐标系到轨道坐标系的旋转矩阵

步骤2.4，生成测姿坐标系到J2000的旋转矩阵。根据以下公式，计算测姿坐标系到J2000坐标系的旋转矩阵：

其中：

步骤2.5，叠加星敏测姿的高频可测颤振。颤振公式如下：

ω=2πf

(其中当y<σ时，y=0)

上式表示本次仿真过程中所采用的星敏颤振叠加的具体方法，相位角取0，对于采样中多个频率的颤振累加到相应的姿态上。设姿态测量设备的测姿灵敏度为σ，当上述累加值超过σ时，说明此颤振可以被姿态测量设备探测出来，需要将此颤振添加到J2000坐标系的姿态上面。

步骤2.6，叠加测姿系统误差和测姿随机误差。根据测姿系统随机误差的设计指标，按照模拟高斯白噪声的模式模拟一组符合测姿系统随机误差指标的三轴欧拉姿态角随机误差，然后再将随机误差及设定的系统误差直接叠加到J2000坐标系下的三轴欧拉角上。

在添加的时候需要将上述测姿坐标系相对于J2000的旋转矩阵按照转序要求变化为三轴欧拉角，在三轴欧拉角上添加测姿随机误差和测姿系统误差，然后再将添加了测姿误差的三轴欧拉角转化为旋转矩阵。

步骤2.7，生成本体坐标系到测姿坐标系的旋转矩阵。具体的，根据设定的旋转矩阵形成顺序，根据测姿仪器三轴欧拉安装角以及安装角的系统误差，形成本体坐标系到测姿坐标系的旋转矩阵

步骤2.8，生成相机坐标系到本体坐标系的旋转矩阵。具体的，根据设定的旋转矩阵形成顺序，根据传感器三轴欧拉安装角以及安装角的系统误差，形成相机坐标系到本体坐标系的旋转矩阵

步骤2.9，生成相机坐标系到J2000坐标系的旋转矩阵。根据上述制作的测姿坐标系相对于J2000的旋转矩阵、本体坐标系到测姿坐标系的旋转矩阵、相机坐标系到本体坐标系的旋转矩阵，按如下式子构成相机坐标系到J2000坐标系的旋转矩阵：

步骤2.10，叠加星敏测姿的高频不可测颤振以及相机的整体颤振。仿照上述颤振的生成方法，生成星敏不可测部分的颤振以及相机的整体颤振。在添加的时候需要将上述测姿坐标系相对于J2000的旋转矩阵按照转序要求变化为三轴欧拉角，在三轴欧拉角上添加颤振，然后在将添加了颤振的三轴欧拉角转化为旋转矩阵。

步骤2.11，叠加姿态量化误差。由于星上存储设备具有一定存储瓶颈，对姿态测量数据的量化位数有限造成的误差。因此应该根据量化指标，将生成的姿态角数值按照这个指标进行一定的截断，得到带量化误差的姿态。

步骤2.12，生成姿态四元数。根据如下公式制作姿态四元组：

上述公式可以根据下面对应关系来计算：

步骤2.14，在计算机显示器上显示仿真结果。

实施例三

实施例三提供了一种陀螺定姿或者陀螺星敏联合定姿的低轨卫星误差姿态生成方法，参见附图3，所述方法具体包括以下步骤：：

步骤3.1，计算三轴欧拉角速度。根据三轴姿态稳定度的设计指标，按照模拟高斯白噪声的模式模拟一组符合三轴稳定度指标的三轴欧拉姿态角速度。

步骤3.2，计算三轴欧拉角。根据三轴欧拉姿态角速度和假定景中心时刻的三轴欧拉角，根据如下公式的积分获得三轴欧拉姿态角：

其中α_t-Δt为前一时刻的三轴欧拉姿态角，为前一时刻的三轴欧拉角速度，α_t(t)为当前时刻的三轴欧拉角。这里需要说明的一点是，姿态稳定度有一个限值，当累加的姿态超过限值时可以被姿态测量设备检测出来，此时采用的策略是将累加后的值减去2倍的欧拉角速度累加值：

步骤3.3，生成陀螺坐标系到轨道坐标系旋转矩阵。具体的，根据设定的旋转矩阵形成顺序，根据三轴欧拉角，形成测姿坐标系到轨道坐标系的旋转矩阵步骤3.4，生成测姿坐标系到J2000的旋转矩阵。根据以下公式，计算测姿坐标系到J2000坐标系的旋转矩阵：

其中：

步骤3.5，叠加星敏测姿的高频可测颤振。颤振公式如下：

ω=2πf

(其中当y<σ时，y=0)

步骤3.6，将上述添加了颤振的三轴欧拉角从J2000坐标系转化到轨道坐标系下，转化方法类似步骤3.4。

步骤3.7，叠加测姿系统误差和测姿随机误差。根据测姿系统随机误差的设计指标，按照模拟高斯白噪声的模式模拟一组符合测姿系统随机误差指标的三轴欧拉姿态角随机误差，然后再将随机误差及设定的系统误差直接叠加到轨道坐标系下的三轴欧拉角上。

步骤3.8，生成测姿坐标系相对于J2000的旋转矩阵。根据设定的旋转矩阵形成顺序，根据三轴欧拉角，形成测姿坐标系相对于轨道坐标系的旋转矩阵根据以下公式，计算测姿坐标系相对于J2000坐标系的旋转矩阵：

步骤3.9，生成本体坐标系到测姿坐标系的旋转矩阵。根据设定的旋转矩阵形成顺序，根据测姿仪器三轴欧拉安装角以及安装角的系统误差，形成本体坐标系相对测姿坐标系的旋转矩阵

步骤3.10，生成相机坐标系到本体坐标系的旋转矩阵。根据设定的旋转矩阵形成顺序，根据传感器三轴欧拉安装角以及安装角的系统误差，形成相机坐标系相对本体坐标系的旋转矩阵

步骤3.11，生成相机坐标系到J2000坐标系的旋转矩阵。根据上述制作的测姿坐标系相对于J2000的旋转矩阵、本体坐标系到测姿坐标系的旋转矩阵、相机坐标系到本体坐标系的旋转矩阵，按如下式子构成相机坐标系到J2000坐标系的旋转矩阵：

步骤3.12，叠加星敏测姿的高频不可测颤振以及相机的整体颤振。仿照上述颤振的生成方法，生成星敏不可测部分的颤振以及相机的整体颤振。在添加的时候需要将上述测姿坐标系相对于J2000的旋转矩阵按照转序要求变化为三轴欧拉角，在三轴欧拉角上添加颤振，然后在将添加了颤振的三轴欧拉角转化为旋转矩阵。

步骤3.13，叠加姿态量化误差。由于星上存储设备具有一定存储瓶颈，对姿态测量数据的量化位数有限造成的误差。因此应该根据量化指标，将生成的姿态角数值按照这个指标进行一定的截断，得到带量化误差的姿态。

步骤3.14，生成姿态四元数。根据如下公式制作姿态四元组：

上述公式可以根据下面对应关系来计算：

步骤3.16，在计算机显示器上显示仿真结果。

以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。

Claims

1.一种低轨卫星标称姿态仿真方法，其特征在于，所述方法具体包括以下步骤：

步骤1.1，计算三轴欧拉角速度，具体的，根据三轴姿态稳定度的设计指标，按照模拟高斯白噪声的模式模拟一组符合三轴稳定度指标的三轴欧拉角速度；

步骤1.2，计算三轴欧拉角，具体的，根据三轴欧拉角速度和假定景中心时刻的三轴欧拉角，根据如下公式的积分获得三轴欧拉角：

<mrow> <msub> <mi>&alpha;</mi> <mi>t</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>&Integral;</mo> <msub> <mover> <mi>&alpha;</mi> <mo>&CenterDot;</mo> </mover> <mrow> <mi>t</mi> <mo>-</mo> <mi>&Delta;</mi> <mi>t</mi> </mrow> </msub> <mi>d</mi> <mi>t</mi> <mo>+</mo> <msub> <mi>&alpha;</mi> <mrow> <mi>t</mi> <mo>-</mo> <mi>&Delta;</mi> <mi>t</mi> </mrow> </msub> </mrow>

其中α_t-Δt为前一时刻的三轴欧拉角，为前一时刻的三轴欧拉角速度，α_t(t)为当前时刻的三轴欧拉角；

步骤1.4，生成测姿坐标系到J2000的旋转矩阵，根据以下公式，计算测姿坐标系到J2000坐标系的旋转矩阵以及相应的J2000坐标系下姿态角：

其中：

<mrow> <msubsup> <mi>R</mi> <mrow> <mi>o</mi> <mi>r</mi> <mi>b</mi> <mi>i</mi> <mi>t</mi> </mrow> <mrow> <mi>J</mi> <mn>2000</mn> </mrow> </msubsup> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mrow> <mo>(</mo> <msub> <mi>X</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mi>X</mi> </msub> </mtd> <mtd> <msub> <mrow> <mo>(</mo> <msub> <mi>Y</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mi>X</mi> </msub> </mtd> <mtd> <msub> <mrow> <mo>(</mo> <msub> <mi>Z</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mi>X</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mrow> <mo>(</mo> <msub> <mi>X</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mi>Y</mi> </msub> </mtd> <mtd> <msub> <mrow> <mo>(</mo> <msub> <mi>Y</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mi>Y</mi> </msub> </mtd> <mtd> <msub> <mrow> <mo>(</mo> <msub> <mi>Z</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mi>Y</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mrow> <mo>(</mo> <msub> <mi>X</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mi>Z</mi> </msub> </mtd> <mtd> <msub> <mrow> <mo>(</mo> <msub> <mi>Y</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mi>Z</mi> </msub> </mtd> <mtd> <msub> <mrow> <mo>(</mo> <msub> <mi>Z</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mi>Z</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>;</mo> <msub> <mover> <mi>Z</mi> <mo>&RightArrow;</mo> </mover> <mn>2</mn> </msub> <mo>=</mo> <mfrac> <mrow> <mover> <mi>P</mi> <mo>&RightArrow;</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mo>|</mo> <mo>|</mo> <mover> <mi>P</mi> <mo>&RightArrow;</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>|</mo> </mrow> </mfrac> <mo>,</mo> <msub> <mover> <mi>X</mi> <mo>&RightArrow;</mo> </mover> <mn>2</mn> </msub> <mo>=</mo> <mfrac> <mrow> <mover> <mi>V</mi> <mo>&RightArrow;</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>&Lambda;</mi> <msub> <mover> <mi>Z</mi> <mo>&RightArrow;</mo> </mover> <mn>2</mn> </msub> </mrow> <mrow> <mo>|</mo> <mo>|</mo> <mover> <mi>V</mi> <mo>&RightArrow;</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>&Lambda;</mi> <msub> <mover> <mi>Z</mi> <mo>&RightArrow;</mo> </mover> <mn>2</mn> </msub> <mo>|</mo> <mo>|</mo> </mrow> </mfrac> <mo>,</mo> <msub> <mover> <mi>Y</mi> <mo>&RightArrow;</mo> </mover> <mn>2</mn> </msub> <mo>=</mo> <msub> <mover> <mi>Z</mi> <mo>&RightArrow;</mo> </mover> <mn>2</mn> </msub> <mi>&Lambda;</mi> <msub> <mover> <mi>X</mi> <mo>&RightArrow;</mo> </mover> <mn>2</mn> </msub> </mrow>

V(t)＝[Xv_s Yv_s Zv_s]^T，

式中，为测姿坐标系到J2000坐标系的旋转矩阵；为轨道坐标系到J2000坐标系的旋转矩阵；为测姿坐标系到轨道坐标系的旋转矩阵；[(Z₂)_X,(Z₂)_Y,(Z₂)_Z]^T为归一化位置向量的三个分量；为卫星质心的位置矢量，为卫星质心的速度矢量；X_s,Y_s,Z_s和分别为卫星质心在J2000坐标系中的位置和速度；

步骤1.5，叠加星敏测姿的高频可测颤振；颤振公式如下：

<mrow> <mi>T</mi> <mo>=</mo> <mfrac> <mrow> <mn>2</mn> <mi>&pi;</mi> </mrow> <mi>&omega;</mi> </mfrac> </mrow>

<mrow> <mi>f</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mi>T</mi> </mfrac> <mo>=</mo> <mfrac> <mi>&omega;</mi> <mrow> <mn>2</mn> <mi>&pi;</mi> </mrow> </mfrac> </mrow>

ω＝2πf

<mrow> <msubsup> <mi>R</mi> <mrow> <mi>c</mi> <mi>a</mi> <mi>m</mi> <mi>e</mi> <mi>r</mi> <mi>a</mi> </mrow> <mrow> <mi>J</mi> <mn>2000</mn> </mrow> </msubsup> <mo>=</mo> <msubsup> <mi>R</mi> <mrow> <mi>s</mi> <mi>t</mi> <mi>a</mi> <mi>r</mi> </mrow> <mrow> <mi>J</mi> <mn>2000</mn> </mrow> </msubsup> <mo>&CenterDot;</mo> <msubsup> <mi>R</mi> <mrow> <mi>b</mi> <mi>o</mi> <mi>d</mi> <mi>y</mi> </mrow> <mrow> <mi>s</mi> <mi>t</mi> <mi>a</mi> <mi>r</mi> </mrow> </msubsup> <mo>&CenterDot;</mo> <msubsup> <mi>R</mi> <mrow> <mi>c</mi> <mi>a</mi> <mi>m</mi> <mi>e</mi> <mi>r</mi> <mi>a</mi> </mrow> <mrow> <mi>b</mi> <mi>o</mi> <mi>d</mi> <mi>y</mi> </mrow> </msubsup> </mrow>

其中，

步骤1.9，生成姿态四元数；根据如下公式生成姿态四元组：

其中，q₁、q₂、q₃、q₄为姿态四元数的四个分量；

上述公式可以根据下面对应关系来计算：

其中，R为姿态四元数对应旋转矩阵；

步骤1.11，在计算机显示器上显示仿真结果。

2.根据权利要求1所述的低轨卫星标称姿态仿真方法，其特征在于，在所述步骤1.5中，所述颤振模型是许多个波函数的叠加，如下公式所示：

其中当y＜σ时，y＝0。

3.根据权利要求1或2所述的低轨卫星标称姿态仿真方法，其特征在于，在所述步骤1.5中，在添加颤振的时候将上述测姿坐标系相对于J2000的旋转矩阵按照转序要求变化为三轴欧拉角，在三轴欧拉角上添加颤振，然后再将添加了颤振的三轴欧拉角转化为旋转矩阵。

4.一种星敏定姿低轨卫星误差姿态生成方法，其特征在于，所述方法具体包括以下步骤：

步骤2.1，计算三轴欧拉角速度；根据三轴姿态稳定度的设计指标，按照模拟高斯白噪声的模式模拟一组符合三轴稳定度指标的三轴欧拉角速度；

步骤2.2，计算三轴欧拉角；根据三轴欧拉角速度和假定景中心时刻的三轴欧拉角，根据如下公式的积分获得三轴欧拉角：

其中：

V(t)＝[Xv_s Yv_s Zv_s]^T，

步骤2.5，叠加星敏测姿的高频可测颤振；颤振公式如下：

ω＝2πf

步骤2.6，叠加测姿系统误差和测姿随机误差；根据测姿系统随机误差的设计指标，按照模拟高斯白噪声的模式模拟一组符合测姿系统随机误差指标的三轴欧拉角随机误差，然后再将随机误差及设定的系统误差直接叠加到J2000坐标系下的三轴欧拉角上；

步骤2.9，生成相机坐标系到J2000坐标系的旋转矩阵；根据上述制作的测姿坐标系相对于J2000的旋转矩阵、本体坐标系到测姿坐标系的旋转矩阵、相机坐标系到本体坐标系的旋转矩阵，按如下式子构成相机坐标系到J2000坐标系的旋转矩阵：

其中，

步骤2.11，叠加姿态量化误差；

步骤2.12，生成姿态四元数；根据如下公式制作姿态四元组：

其中，q₁、q₂、q₃、q₄为姿态四元数的四个分量；

上述公式可以根据下面对应关系来计算：

其中，R为姿态四元数对应旋转矩阵；

步骤2.14，在计算机显示器上显示仿真结果。

5.根据权利要求4所述的星敏定姿低轨卫星误差姿态生成方法，其特征在于，在所述步骤2.5中，所述颤振模型是许多个波函数的叠加，如下公式所示：

其中当y＜σ时，y＝0。

6.根据权利要求4或5所述的星敏定姿低轨卫星误差姿态生成方法，其特征在于，在所述步骤2.5，2.6和2.7中，在添加颤振的时候将上述测姿坐标系相对于J2000的旋转矩阵按照转序要求变化为三轴欧拉角，在三轴欧拉角上添加颤振，然后再将添加了颤振的三轴欧拉角转化为旋转矩阵。

7.一种陀螺定姿或者陀螺星敏联合定姿的低轨卫星误差姿态生成方法，所述方法具体包括以下步骤：

步骤3.1，计算三轴欧拉角速度；根据三轴姿态稳定度的设计指标，按照模拟高斯白噪声的模式模拟一组符合三轴稳定度指标的三轴欧拉角速度；

步骤3.2，计算三轴欧拉角；根据三轴欧拉角速度和假定景中心时刻的三轴欧拉角，根据如下公式的积分获得三轴欧拉角：

其中：

V(t)＝[Xv_s Yv_s Zv_s]^T，

步骤3.5，叠加星敏测姿的高频可测颤振，颤振公式如下：

ω＝2πf

步骤3.7，叠加测姿系统误差和测姿随机误差；根据测姿系统随机误差的设计指标，按照模拟高斯白噪声的模式模拟一组符合测姿系统随机误差指标的三轴欧拉角随机误差，然后再将随机误差及设定的系统误差直接叠加到轨道坐标系下的三轴欧拉角上；

其中，

步骤3.13，叠加姿态量化误差；

步骤3.14，生成姿态四元数；根据如下公式制作姿态四元组：

其中，q₁、q₂、q₃、q₄为姿态四元数的四个分量；

上述公式可以根据下面对应关系来计算：

其中，R为姿态四元数对应旋转矩阵；

步骤3.16，在计算机显示器上显示仿真结果。

8.根据权利要求7所述的陀螺定姿或者陀螺星敏联合定姿的低轨卫星误差姿态生成方法，其特征在于，在所述步骤3.5中，所述颤振模型是许多个波函数的叠加，如下公式所示：

其中当y＜σ时，y＝0。

9.根据权利要求7或8所述的陀螺定姿或者陀螺星敏联合定姿的低轨卫星误差姿态生成方法，其特征在于，在所述步骤3.5和3.12中，在添加颤振的时候将上述测姿坐标系相对于J2000的旋转矩阵按照转序要求变化为三轴欧拉角，在三轴欧拉角上添加颤振，然后再将添加了颤振的三轴欧拉角转化为旋转矩阵。