CN101226561A - 用于微型航天器姿态轨道控制系统的微型仿真支持系统及工作方法 - Google Patents

Abstract

一种用于微型航天器姿态轨道控制系统的微型仿真支持系统及工作方法，属于航天器姿态轨道控制仿真领域。该系统包括接口输入电路、处理单元和输出接口电路；其工作方法是通过接口电路采样星载计算机的输入指令；根据执行机构的模型计算出控制力与控制力矩；计算出的扰动力与扰动力矩；结合控制对象的力学特性利用数值计算方法更新航天器的当前姿态与轨道；根据各星上传感器模型计算出传感器数据；通过接口电路输出，供星载计算机采集传感器信号，从而支持实现航天器姿态轨道控制系统的全闭环半物理仿真。体积、功耗小，便于携带和桌面联试；基于实时多任务操作系统开发软件，维护简单，移植性可靠性有保证。

Description

用于微型航天器姿态轨道控制系统的微型仿真支持系统及工作方法

技术领域：本发明属于航天器姿态轨道控制仿真技术领域。

背景技术：随着空间技术的飞速发展，对航天器的研制成本和速度提出了更高的要求；尤其对微小卫星的研制速度与成本的要求更加苛刻，所以应当尽可能采用新的研制手段和工具，提高卫星的制造速度，降低研制成本，并提高可靠性。

姿态轨道控制系统是航天器最关键的分系统之一，其任务就是按照飞行计划的要求，实现航天器在各个飞行阶段的姿态和轨道的测量，以及姿态和轨道的控制。加速姿态轨道控制分系统的研发、降低开发成本是提高航天器研制水平的一个重要的方向，对辅助开发的工具提出了更高的要求。目前相关的电测设备体积庞大，不方便桌面联试；相关的半物理仿真设备成本很高，联试复杂，并且受到一定条件的限制。

同时目前软件技术迅猛发展，很多设备都采用了软件代替部分硬件的方式，降低了仿真设备的成本，又可以通过软件升级的方式提高系统的性能，减少了再投资，提高了资源的利用率；同时随着半导体技术的进步，嵌入式计算机的成本不断降低，而性能却不断提高，在很多场合能完成较好的高性能计算；嵌入式多任务操作系统也为该类产品的开发提供了便利。这些新技术为开发价廉物美的仿真设备保证了可能性。

所以辅助航天器姿态轨道控制系统研发的微型仿真支持系统是航天控制仿真设备发展的重要方向。

发明内容：

为便于微型航天器姿态轨道控制系统的研制，本发明提供了用于微型航天器姿态轨道控制系统的微型仿真支持系统，该航天器微型仿真支持系统包括输入接口、处理单元与输出接口三部分，所述输入接口主要由AD信号转换电路、IO数字量输入电路组成，用于接收导航计算机的模拟量和数字量指令输出；所述输出接口主要由DA、串口构成，用于输出各传感器的仿真测量数据；所述处理单元包括处理器、内存与存储器，程序固化在非易失性存储器中，程序的运行在处理器和内存中完成，主要有以下处理步骤：

首先通过接口电路采样星载计算机的输入指令，然后根据执行机构的模型计算出控制力与控制力矩；同时维护一个时间系统计算出各种空间物理参数，并在此基础上计算出扰动力与扰动力矩；然后结合控制对象的力学特性利用数值计算方法更新航天器的当前姿态与轨道；进而根据各星上传感器模型计算出各种传感器的测量数据；最后通过接口电路输出各种传感器数据。

工作方法的具体步骤是：

(一)通过接口电路采样星载计算机的输入信号，将输入信号转换成物理量；

(二)计算出执行机构对航天器的控制力和控制力矩，包括磁力矩器的磁力矩、转速控制动量轮的控制力矩、喷气系统的推力，

1)磁力矩器的磁力矩，

仿真磁力矩器的输出磁矩：

${\stackrel{\→}{m}}_r={\stackrel{\→}{m}}_q+\stackrel{\→}{e}---\left(1\right)$

式中，表示磁矩误差，

表示期望的磁矩输出，

表示实际的磁矩输出，与当地地磁场作用，形成控制力矩为

$\stackrel{\→}{L}={\stackrel{\→}{m}}_r\×\stackrel{\→}{B}---\left(2\right)$

2)转速控制动量轮的控制力矩，

(a)计算动量轮的转速，动量轮转速的表达式为

$\mathrm{\ω}=\frac{{\mathrm{\ω}}_n^2}{s^2+2\mathrm{\δ}{\mathrm{\ω}}_{n}s+{\mathrm{\ω}}_n^2}{\mathrm{\ω}}_{\mathrm{in}}---\left(3\right)$

式中， $\frac{{\mathrm{\ω}}_n^2}{s^2+2\mathrm{\δ}{\mathrm{\ω}}_{n}s+{\mathrm{\ω}}_n^2}$ 为二阶响应环节，s表示微分算子，并且|sω|＜door，door表示角加速度的门限，ω_in为期望的输入，ω表示真实的动量轮转速；

(b)根据上述式(3)计算卫星所受的控制力矩

$L=-J\stackrel{\·}{\mathrm{\ω}}---\left(4\right)$

式中的J为动量轮转动惯量，

表示角加速度；

3)喷气系统的推力，

喷气系统推力的表达式为：

$\left\{\begin{array}{cc}F+\mathrm{\ϵ}& \mathrm{On}\\ 0& \mathrm{Off}\end{array}\right.---\left(5\right)$

式中，F是推力器的标称推力，ε为推力误差，On表示阀门开启，Off表示阀门关闭；

(三)计算空间物理参数，利用太阳月球运动模型计算出太阳月球的位置，利用大气模型计算出当地大气的密度，利用地磁模型计算出当地磁场矢量，

(四)计算环境摄动力和环境干扰力矩：利用上一时刻的姿态轨道参数、环境物理参数、摄动力和干扰力矩模型，计算出环境摄动力和环境干扰力矩：

1)地球非中心引力的计算，

采用下式计算地球非中心引力：

$V=\frac{\mathrm{\μ}}{r}\underset{n=2}{\overset{N}{\mathrm{\Σ}}}{\left(\frac{R_e}{r}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{P}_n^m\left(\sin \mathrm{\φ}\right)(c_{\mathrm{nm}}\cos \mathrm{m\λ}+s_{\mathrm{nm}}\sin \mathrm{m\λ})---\left(6\right)$

式中，μ为地球引力常数，r为地心距，R_e为地球半径，P_n ^m勒让德多项式，λ、φ为地心经、纬度，c_nm为带谐项和田谐项；

2)三体摄动力的计算，

根据第(三)步计算的太阳、月球在空间的方位，计算出太阳和月球对航天器的摄动，

太阳引起的摄动加速度是：

${\stackrel{\→}{a}}_s=-{\mathrm{\μ}}_s(\frac{{\stackrel{\→}{r}}_{\mathrm{sv}}}{r_{\mathrm{sv}}^3}-\frac{{\stackrel{\→}{r}}_{\mathrm{se}}}{r_{\mathrm{se}}^3})---\left(7\right)$

月球引起的摄动加速度是：

${\stackrel{\→}{a}}_m=-{\mathrm{\μ}}_m(\frac{{\stackrel{\→}{r}}_{\mathrm{mv}}}{r_{\mathrm{mv}}^3}-\frac{{\stackrel{\→}{r}}_{\mathrm{me}}}{r_{\mathrm{me}}^3})---\left(8\right)$

式中，μ_s和μ_m为日月引力常数，

和为日月至卫星的矢量，和

为日月至地球的矢量；

3)大气摄动力的计算，

通过软件设置面质比，计算大气摄动力的表达式为：

${\stackrel{\→}{F}}_a=\frac{c_d}{2}\left(\mathrm{\ρ}{V}_R^2{A}_{\mathrm{\ρ}}\right)\stackrel{\→}{v}---\left(9\right)$

式中，c_d为阻力系数；ρ为大气密度；V_R为卫星与大气之间的相对速度；A为迎风面面积；为来流的单位矢量；

4)光压摄动力的计算，

通过软件设置面质比，计算光压摄动力的表达式为：

$\left\{\begin{array}{c}{\stackrel{\→}{F}}_n=\frac{\mathrm{IS}}{c}\cos \mathrm{\φ}\{\frac{2}{3}\mathrm{\ρ\μ}+[1+\mathrm{\ρ}(1-\mathrm{\μ})\left]\cos \mathrm{\φ}\right\}\\ {\stackrel{\→}{F}}_{\mathrm{\τ}}=\frac{\mathrm{IS}}{c}\cos \mathrm{\φ}\sin \mathrm{\φ}(1-\mathrm{\ρ}+\mathrm{\ρ\μ})\end{array}\right.---\left(10\right)$

式中，

表示法向力，

表示切向力，φ为太阳入射角，I是太阳常数，c是光速，ρ表示反射系数，μ表示散射系数，S表示受力面积；

5)引力梯度力矩的计算，

通过下式计算引力梯度力矩：

$\stackrel{\→}{L}=\frac{3\mathrm{\μ}}{R^5}\tilde{R}I\stackrel{\→}{R}---\left(11\right)$

式中，μ为地球引力常数，

为星地矢量，

表示反对称矩阵；

6)气动力矩的计算，

利用式(9)的大气摄动力，结合摄动力作用力臂

求得气动干扰力矩

${\stackrel{\→}{T}}_g={\stackrel{\→}{F}}_a\×\stackrel{\→}{d}---\left(12\right)$

7)剩磁力矩的计算，

计算方法与式(2)类似，即：

${\stackrel{\→}{T}}_m={\stackrel{\→}{m}}_b\×\stackrel{\→}{B}---\left(12\right)$

式中，

为星体剩磁矩；

8)光压力矩的计算，

可利用式(10)的摄动力，结合对应的力臂，求得光压干扰力矩

${\stackrel{\→}{T}}_l={\stackrel{\→}{F}}_n\×{\stackrel{\→}{d}}_n+{\stackrel{\→}{F}}_{\mathrm{\τ}}\×{\stackrel{\→}{d}}_{\mathrm{\τ}}---\left(14\right)$

式中，

和

分布表示法向力与切向力的力臂；

(五)利用上述已计算好的控制力、控制力矩、环境摄动力和环境干扰力矩，结合航天器力学模型和数值计算方法更新航天器的姿态轨道参数，其中：

1)航天器的轨道力学模型是：

$a=({\stackrel{\→}{F}}_c+{\stackrel{\→}{F}}_O+{\stackrel{\→}{F}}_{\mathrm{\ϵ}})/m---\left(15\right)$

式中，a为航天器的加速度，

为卫星所受的控制力，为地球中心引力，

为各摄动力，m是航天器的质量；

2)航天器的姿态力学模型是：

$\stackrel{\stackrel{\·}{\→}}{h}+\stackrel{\→}{\mathrm{\ω}}\×\stackrel{\→}{h}={\stackrel{\→}{L}}_c+{\stackrel{\→}{L}}_{\mathrm{\ϵ}}---\left(16\right)$

$\stackrel{\·}{q}=0.5q\stackrel{\‾}{\mathrm{\ω}}---\left(17\right)$

式中，

表示航天器动量矩大小的导数，

为航天器动量矩，为控制力矩，

为干扰力矩，

为角速度，q为姿态四元数， $\stackrel{\‾}{\mathrm{\ω}}=\left[\begin{array}{cc}0& \stackrel{\→}{\mathrm{\ω}}\end{array}\right].$

轨道的数值方法和姿态的数值计算方法均采用RK；

(六)利用上述第(五)步骤更新后的航天器姿态轨道参数、第(三)步骤的太阳、月亮与当前时间等信息，结合传感器模型计算传感器应该输出的物理量大小：

1)陀螺角速度测量值计算，

${\mathrm{\ω}}_{\mathrm{gc}}=(I-\mathrm{\Λ})(I-\mathrm{\Δ})Q_b^g{\mathrm{\ω}}_b-b-\mathrm{\η}---\left(18\right)$

式中，Q_b ^g表示体坐标系到陀螺坐标系的转换矩阵，b为随机游走误差，Λ为标定因子误差，Δ为非正交安装误差；

2)磁强计磁场强度测量值计算，

${\stackrel{\→}{B}}_k={(I+D)}^{-1}({\mathrm{\Θ}}^T{A}_k{\stackrel{\→}{M}}_k+b+\mathrm{\ϵ})---\left(19\right)$

式中，

为磁强计输出，D表示刻度因子误差及非正交化矩阵，Θ为非对准矩阵，

为某参考坐标系下的地磁矢量，A_k为参考坐标系到磁强计安装坐标系的转换矩阵，b为磁强计的常值偏置，

由第(三)步计算的磁场强度获得；

3)模拟太阳敏感器输出电压计算

U＝U₀cosθ           (20)

式中，θ表示入射角，U₀表示硅片正射时模拟太阳敏感器的输出电压，U表示当前的电压输出；

4)GPS定位数据计算

x＝x_r+d+υ       (21)

$\stackrel{\·}{d}=-\frac{1}{\mathrm{\τ}}d+\mathrm{\ϵ}---\left(22\right)$

式中，x_r是卫星真实的位置，x为GPS输出的定位数据，d为一阶马尔科夫过程，υ为噪声，τ为马尔科夫过程的相关时间，ε为驱动噪声。GPS定位时需要初始化，程序中应仿真该初始化的时间；

(七)通过输出接口电路输出各种数据，根据各传感器的实际输出方式，如果是模拟器件，则要将输出的物理量转换成电压值；如果是串口输出的设备，则要按固定的格式打包；最终将所有数据通过DA，串口和IO输出接口输出。

本发明中，除接口电路与嵌入式计算机外，系统的其余功能采用算法和物理模型描述，并用软件实现，利于后期升级并节约成本；采用与星上传感器电性能一致的接口电路，实现传感器电性能级别的仿真；采用高性能嵌入式系统，并基于实时多任务操作系统开发软件开发。

本发明的优点是：

1、采用民用嵌入式系统开发仿真设备，成本、功耗低，体积小，适合桌面联试，携带方便；

2、除接口电路与嵌入式计算机外，其余都采用软件实现，具有改造升级容易的特点，扩展性强，也降低了整个系统的成本；

3、基于实时多任务操作系统开发，既提高了仿真设备的可靠性，降低了研制该设备的难度，并且具有较好的移植性；

4、仿真内容全面，能支持航天器姿态轨道控制系统的全过程开发，考核从信号采集到指令输出的大部分环节。

附图说明

图1是用于微型航天器姿态轨道控制系统的卫星仿真支持系统组成框图

图2是本发明的系统层次示意图

图3是本发明的软件流程图

具体实施方式：

下面结合附图进一步说明本发明。

微型仿真支持系统从构成的层次来分，主要有5个部分：控制指令输入、执行机构模型解算、星体动力学、传感器数据仿真以及传感器数据输出，如图1所示，这5个部分按时序依次计算，往复循环。微型仿真支持系统的输入层主要采用AD或者IO进行高频采样，捕捉星载计算机控制指令的变化；执行机构模型层、星体力学模型与传感器模型层主要就是调用各种算法与模型，实时计算航天器姿轨参数与传感器数据；输出层通过刷新DA、串口等输出，输出各种传感器的数据，供星载计算机采集使用。

图1中的“……”表示其余各种执行机构和传感器，整个系统除必备的硬件接口与计算机外，都有软件构成。应用软件运行在实时操作系统上，整个软件系统运行在嵌入式硬件上，并且支持关键数据的输出，便于监控。

微型仿真支持系统从硬件上看主要有输入接口电路、含有数据处理程序的处理单元和输出接口电路，其中输入接口电路主要由AD信号转换电路、IO数字量输入电路组成，用于接收导航计算机的模拟量和数字量指令输出；数据处理单元收集控制指令，通过各种模型计算出个传感器的输出，处理程序固化在非易失性的存储器中，计算机可由嵌入式微处理器等实现；最后将模拟式传感器数据通过DA输出，数字式传感器通过串口输出，1O数字量接口输出某些传感器上的同步信号。如图2所示。

应用软件运行在实时嵌入式操作系统上，应用软件的数据处理的基本原理如图3所示，程序通过多任务操作系统周期触发运行。其详细过程如下：

(1)周期性地采集各种控制指令的输入，通过软件周期扫描系统的硬件接口电路上的输入信号，并将其存储在内存中；然后，将输入的电信号转换成相应的物理量：将模拟电压转换成对应执行机构期望输出的物理量，根据IO数字量的电平高低判断对应执行机构的工作状态；

(2)计算控制力和控制力矩：利用各种执行机构的模型与期望的输入，计算出执行机构对航天器的控制力和控制力矩。典型的执行机构如下：

1)磁力矩器的磁力矩

磁力矩器的磁矩输出是在期望的输出上增加非线性等因素，可以较好地仿真磁力矩器的输出。

${\stackrel{\→}{m}}_r={\stackrel{\→}{m}}_q+\stackrel{\→}{e}---\left(1\right)$

式中，

表示磁矩误差，

表示期望的磁矩输出，

表示实际的磁矩输出。

与当地地磁场作用，形成控制力矩。

$\stackrel{\→}{L}={\stackrel{\→}{m}}_r\×\stackrel{\→}{B}---\left(2\right)$

2)转速控制的动量轮的控制力矩

转速控制动量轮的输入信号代表了动量轮期望的转速，动量轮对期望转速的相应可建模成二阶响应曲线，对其动量矩的微分，即可获得动量轮的控制力矩；一般而言，该类动量轮都会对角加速度限饱和，模型也应该考虑该因素，其角速度响应模型是：

$\mathrm{\ω}=\frac{a}{s^2+2\mathrm{\δ}{\mathrm{\ω}}_{n}s+{\mathrm{\ω}}_n^2}{\mathrm{\ω}}_{\mathrm{in}}---\left(3\right)$

式中， $\frac{{\mathrm{\ω}}_n^2}{s^2+2\mathrm{\δ}{\mathrm{\ω}}_{n}s+{\mathrm{\ω}}_n^2}$ 为二阶响应环节，s表示微分算子，|sω|＜door，door表示角加速度的门限，ω_in为期望的输入，ω表示真实的动量轮转速。

所以卫星所受的控制力矩为

$L=-J\stackrel{\·}{\mathrm{\ω}}---\left(4\right)$

式中的J为动量轮转动惯量，

表示角加速度。

3)喷气系统的推力

喷气系统通过喷嘴喷射出高速气流，从而系统获得一个较稳定的推力。

$\left\{\begin{array}{cc}F+\mathrm{\ϵ}& \mathrm{On}\\ 0& \mathrm{Off}\end{array}\right.---\left(5\right)$

式中，F是推力器的标称推力，ε为推力误差，On表示阀门开启，Off表示阀门关闭。

(3)维护一个时间系统，该时间作为系统的标准时间，并进行环境物理参数的计算：利用标准时间系统和太阳月球运动模型计算时刻太阳、月球的位置；利用大气模型计算出当地大气的密度，利用IGRF地磁模型计算出当地磁场矢量。

(4)计算环境摄动力和环境干扰力矩：利用上一时刻的姿态轨道参数、环境物理数、摄动力和干扰力矩模型，计算出环境摄动力和环境干扰力矩：

1)地球非中心引力的计算

地球的高阶摄动项是干扰近地卫星轨道的主要因素，引力模型有多种，主要表现在参数的略有不同，但都可以达到较高的精度，阶数的选择也可根据计算量和仿真任务进行参数选择，地球非中心的引力位函数是：

$V=\frac{\mathrm{\μ}}{r}\underset{n=2}{\overset{N}{\mathrm{\Σ}}}{\left(\frac{R_e}{r}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{P}_n^m\left(\sin \mathrm{\φ}\right)(c_{\mathrm{nm}}\cos \mathrm{m\λ}+s_{\mathrm{nm}}\sin \mathrm{m\λ})---\left(6\right)$

式中，μ为地球引力常数，r为地心距，R_e为地球半径，P_n ^m勒让德多项式，λ、φ为地心经、纬度，c_nm为带谐项和田谐项。

2)三体摄动力的计算

太阳和月球对航天器的运动会发生影响，通过第(3)步计算得到的太阳与月球的方位，进而计算出太阳和月球对航天器的摄动。

太阳引起的摄动加速度是：

${\stackrel{\→}{a}}_s=-{\mathrm{\μ}}_s(\frac{{\stackrel{\→}{r}}_{\mathrm{sv}}}{r_{\mathrm{sv}}^3}-\frac{{\stackrel{\→}{r}}_{\mathrm{se}}}{r_{\mathrm{se}}^3})---\left(7\right)$

月球引起的摄动加速度是：

${\stackrel{\→}{a}}_m=-{\mathrm{\μ}}_m(\frac{{\stackrel{\→}{r}}_{\mathrm{mv}}}{r_{\mathrm{mv}}^3}-\frac{{\stackrel{\→}{r}}_{\mathrm{me}}}{r_{\mathrm{me}}^3})---\left(8\right)$

式中，μ_s和μ_m为日月引力常数，

和

为日月至卫星的矢量，

和为日月至地球的矢量。

3)大气摄动力的计算

对于近地卫星，尤其是面质比大的卫星，大气阻尼因素是不能忽略的，大气阻尼最主要的影响是降低轨道高度。通过软件设置面质比，计算出大气摄动：

${\stackrel{\→}{F}}_a=\frac{c_d}{2}\left(\mathrm{\ρ}{V}_R^2{A}_{\mathrm{\ρ}}\right)\stackrel{\→}{v}---\left(9\right)$

式中，c_d为阻力系数；ρ为大气密度；V_R为卫星与大气之间的相对速度；A为迎风面面积；

为来流的单位矢量。

4)光压摄动力的计算

对于面质比大的卫星，太阳光压对卫星的轨道也能产生显著的影响。通过软件设置面质比，计算出光压摄动：

$\left\{\begin{array}{c}{\stackrel{\→}{F}}_n=\frac{\mathrm{IS}}{c}\cos \mathrm{\φ}\{\frac{2}{3}\mathrm{\ρ\μ}+[1+\mathrm{\ρ}(1-\mathrm{\μ})\left]\cos \mathrm{\φ}\right\}\\ {\stackrel{\→}{F}}_{\mathrm{\τ}}=\frac{\mathrm{IS}}{c}\cos \mathrm{\φ}\sin \mathrm{\φ}(1-\mathrm{\ρ}+\mathrm{\ρ\μ})\end{array}\right.---\left(10\right)$

式中，

表示法向力，

表示切向力，φ为太阳入射角，I是太阳常数，c是光速，ρ表示反射系数，μ表示散射系数，S为照射面积。

5)引力梯度力矩的计算

星体各部分质量在地球引力场中所受的引力略有差别，从而导致对星体的中心产生一个附加力矩。

$L=\frac{3\mathrm{\μ}}{R^5}\tilde{R}I\stackrel{\→}{R}---\left(11\right)$

式中，μ为地球引力常数，

为星地矢量，

表示反对称矩阵。

6)气动力矩的计算

大气对近地卫星产生阻尼作用，该阻尼力如果与星体质心不重合，也会形成对星体扰动的气动力矩。利用式(9)的大气摄动力，结合摄动力作用力臂

求得气动干扰力矩：

${\stackrel{\→}{T}}_g={\stackrel{\→}{F}}_a\×\stackrel{\→}{d}---\left(12\right)$

7)剩磁力矩的计算

卫星星体会残留一些磁性，星体的剩磁与当地地球磁场相互作用，形成对星体扰动的剩磁干扰力矩。计算方法与式(2)类似，即

${\stackrel{\→}{T}}_m={\stackrel{\→}{m}}_b\×\stackrel{\→}{B}---\left(13\right)$

式中，

为星体剩磁矩；

8)光压力矩的计算

光压的合力与质心不一致，也会形成对星体扰动的光压力矩，可利用式(10)的摄动力与对应的力臂，求得光压干扰力矩：

${\stackrel{\→}{T}}_l={\stackrel{\→}{F}}_n\×{\stackrel{\→}{d}}_n+{\stackrel{\→}{F}}_{\mathrm{\τ}}\×{\stackrel{\→}{d}}_{\mathrm{\τ}}---\left(14\right)$

式中，

和

分布表示法向力与切向力的力臂；

(5)利用计算好的控制力、控制力矩、环境摄动力和环境干扰力矩，结合航天器力学模型和数值计算方法更新航天器的姿态轨道参数。

航天器的轨道力学模型是：

$a=({\stackrel{\→}{F}}_c+{\stackrel{\→}{F}}_O+{\stackrel{\→}{F}}_{\mathrm{\ϵ}})/m---\left(16\right)$

式中，α为航天器的加速度，

为控制力，为地球中心引力，

为各摄动力，m是航天器的质量。

航天器的姿态力学模型是：

$\stackrel{\stackrel{\·}{\→}}{h}+\stackrel{\→}{\mathrm{\ω}}\×\stackrel{\→}{h}={\stackrel{\→}{L}}_c+{\stackrel{\→}{L}}_{\mathrm{\ϵ}}---\left(17\right)$

$\stackrel{\·}{q}=0.5q\stackrel{\‾}{\mathrm{\ω}}---\left(18\right)$

式中，

表示航天器动量矩大小的导数，

为航天器动量矩，

为控制力矩，

为干扰力矩，为角速度，q为姿态四元数， $\stackrel{\‾}{\mathrm{\ω}}=\left[\begin{array}{cc}0& \stackrel{\→}{\mathrm{\ω}}\end{array}\right].$

轨道的数值方法和姿态的数值计算方法均采用RK。

(6)利用更新后的航天器姿态轨道参数、太阳、月亮与当前时间等信息，结合传感器模型计算传感器应该输出的物理量大小：

1)陀螺输出值的计算

陀螺能测量相对于惯性空间的角速率，除真实的角速率信息外，陀螺输出还需考虑非正交误差、刻度因素误差，随机游走和高斯噪声等误差。

${\mathrm{\ω}}_{\mathrm{gc}}=(I-\mathrm{\Λ})(I-\mathrm{\Δ})Q_b^g{\mathrm{\ω}}_b-b-\mathrm{\η}---\left(19\right)$

式中，Q_b ^g表示体坐标系到陀螺坐标系的转换矩阵，b为随机游走误差，Λ为标定因子误差，Δ为非正交安装误差。

2)磁强计磁场强度测量值的计算

磁强计的测量精度主要受非正交误差、刻度因子误差，还有常值偏置、高斯噪声等制约，其中常值偏置主要来自于星体剩磁。

${\stackrel{\→}{B}}_k={(I+D)}^{-1}({\mathrm{\Θ}}^T{A}_k{\stackrel{\→}{M}}_k+b+\mathrm{\ϵ})---\left(20\right)$

其中，

为磁强计输出，D表示刻度因子误差及非正交化矩阵，Θ为非对准矩阵，

为某参考坐标系下的地磁矢量，A_k为参考坐标系到磁强计安装坐标系的转换矩阵，b为磁强计的常值偏置。

由第(3)步计算的磁场强度获得。

3)模拟太阳敏感器的电压值计算

模拟太阳敏感器的硅片受太阳光照射产生电路，转换后输出电压

U＝U₀cosθ      (20)

式中，θ表示入射角，U₀表示硅片正射时模拟太阳敏感器的输出电压，U表示实际输出电压。

4)GPS定位数据的输出

GPS能按固定的频率输出定位信息，由于空间GPS接收机工作于高动态环境，解算延迟已经不容忽略，一般为用户使用的方便，某时刻的定位结果会在下时刻输出，所以软件中需具备延迟功能；同时GPS的定位误差具有一定的相关性，采用一阶马尔科夫过程仿真该相关误差；并且采用姿态与轨道信息，判定GPS天线是否指向外太空，如果指向外太空，则GPS能正常工作，否则失效；GPS首次定位时需要初始化，程序中需要仿真该初始化的时间。GPS的定位模型为：

x＝x_r+d+υ               (21)

$\stackrel{\·}{d}=-\frac{1}{\mathrm{\τ}}d+\mathrm{\ϵ}---\left(22\right)$

式中，x_r是卫星真实的位置，x为GPS输出的定位数据，d为一阶马尔科夫过程，υ为噪声，τ为马尔科夫过程的相关时间，ε为驱动噪声。

(7)根据各传感器的实际输出方式，如果是模拟器件，则要将输出的物理量转换成电压值；如果是串口输出的设备，则要按固定的格式打包；最终将所有数据通过DA，串口和IO等输出接口输出。

利用上述系统的构架原理和方法，可构成一个用于微型航天器姿态轨道控制系统开发的仿真支持系统，为星载控制系统的开发提供便利；由于大部分功能由软件实现，很容易扩展新的执行机构与传感器模型。

本发明与传统的电测系统相比较，由于采用嵌入式系统开发，具有成本功耗低，体积质量小，便于携带与桌面联试的优点；同时接口电路采用与真实传感器一致的接口，算法模型具有很高的精度，保证了仿真的高保真；由于大部分功能由软件实现，便于升级；整个应用软件基于多任务实时操作系统开发，具有很高的可靠性与移植性。本发明适合用于航天器的姿态轨道控制系统的开发，特别是微型卫星的姿态轨道控制系统开发，具有良好的应用前景。

Claims (2)

1.一种用于微型航天器姿态轨道控制系统的微型仿真支持系统，其特征在于，包含输入接口电路、处理单元和输出接口电路，其中输入接口电路的输出连接于处理单元的输入，处理单元的输出连于输出接口电路，所述的输入接口电路包含AD转换电路和IO数字信号输入电路，用于接收导航计算机的模拟量和数字量指令的输出；所述处理单元包括处理器、内存与存储器，程序固化在非易失性存储器中，程序的运行在处理器和内存中完成；所述的输出接口电路包含DA转换电路、串口输出电路、数字量输出电路、数据监控输出电路，用于输出仿真测量数据与监控数据。

2.一种如权利要求1所述的用于微型航天器姿态轨道控制系统的微型仿真支持系统的工作方法，其特征在于，该工作方法包括如下步骤：

(一)通过接口电路采样星载计算机的输入信号，将输入信号转换成物理量；

(二)计算出执行机构对航天器的控制力和控制力矩，包括磁力矩器的磁力矩、转速控制动量轮的控制力矩、喷气系统的推力，

1)磁力矩器的磁力矩，

仿真磁力矩器的输出磁矩：

${\stackrel{\→}{m}}_r={\stackrel{\→}{m}}_q+\stackrel{\→}{e}---\left(1\right)$

式中，

表示磁矩误差，

表示期望的磁矩输出，

表示实际的磁矩输出，

与当地地磁场

作用，形成控制力矩为

$\stackrel{\→}{L}={\stackrel{\→}{m}}_r\×\stackrel{\→}{B}---\left(2\right)$

2)转速控制动量轮的控制力矩，

(a)计算动量轮的转速，动量轮转速的表达式为

$\mathrm{\ω}=\frac{{\mathrm{\ω}}_n^2}{s^2+2\mathrm{\δ}{\mathrm{\ω}}_{n}s+{\mathrm{\ω}}_n^2}{\mathrm{\ω}}_{\mathrm{in}}---\left(3\right)$

式中， $\frac{{\mathrm{\ω}}_n^2}{s^2+2\mathrm{\δ}{\mathrm{\ω}}_{n}s+{\mathrm{\ω}}_n^2}$ 为二阶响应环节，s表示微分算子，并且|sω|＜door，door表示角加速度的门限，ω_in为期望的输入，ω表示真实的动量轮转速；

(b)根据上述式(3)计算卫星所受的控制力矩：

$L=-J\stackrel{\·}{\mathrm{\ω}}---\left(4\right)$

式中的J为动量轮转动惯量，

表示角加速度；

3)喷气系统的推力，

喷气系统推力的表达式为：

$\left\{\begin{array}{cc}F+\mathrm{\ϵ}& \mathrm{On}\\ 0& \mathrm{Off}\end{array}\right.---\left(5\right)$

式中，F是推力器的标称推力，ε为推力误差，On表示阀门开启，Off表示阀门关闭；

(三)计算空间物理参数，利用太阳月球运动模型计算出太阳月球的位置，利用大气模型计算出当地大气的密度，利用地磁模型计算出当地磁场矢量，

(四)计算环境摄动力和环境干扰力矩：利用上一时刻的姿态轨道参数、环境物理参数、摄动力和干扰力矩模型，计算出环境摄动力和环境干扰力矩：

1)地球非中心引力的计算，

采用下式计算地球非中心引力：

$V=\frac{\mathrm{\μ}}{r}\underset{n=2}{\overset{N}{\mathrm{\Σ}}}{\left(\frac{R_e}{r}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\Σ}}}{P}_n^m\left(\sin \mathrm{\φ}\right)(c_{\mathrm{nm}}\cos \mathrm{m\λ}+s_{\mathrm{nm}}\sin \mathrm{m\λ})---\left(6\right)$

式中，μ为地球引力常数，r为地心距，R_e为地球半径，P_n ^m勒让德多项式，λ、φ为地心经、纬度，c_nm为带谐项和田谐项；

2)三体摄动力的计算，

根据第(三)步计算的太阳、月球在空间的方位，计算出太阳和月球对航天器的摄动，

太阳引起的摄动加速度是：

${\stackrel{\→}{a}}_s=-{\mathrm{\μ}}_s(\frac{{\stackrel{\→}{r}}_{\mathrm{sv}}}{r_{\mathrm{sv}}^3}-\frac{{\stackrel{\→}{r}}_{\mathrm{se}}}{r_{\mathrm{se}}^3})---\left(7\right)$

月球引起的摄动加速度是：

${\stackrel{\→}{a}}_m=-{\mathrm{\μ}}_m(\frac{{\stackrel{\→}{r}}_{\mathrm{mv}}}{r_{\mathrm{mv}}^3}-\frac{{\stackrel{\→}{r}}_{\mathrm{me}}}{r_{\mathrm{me}}^3})---\left(8\right)$

式中，μ_s和μ_m为日月引力常数，

和

为日月至卫星的矢量，

和

为日月至地球的矢量；

3)大气摄动力的计算，

通过软件设置面质比，计算大气摄动力的表达式为：

${\stackrel{\→}{F}}_a=\frac{c_d}{2}\left(\mathrm{\ρ}{V}_R^2{A}_{\mathrm{\ρ}}\right)\stackrel{\→}{v}---\left(9\right)$

式中，c_d为阻力系数；ρ为大气密度；V_R为卫星与大气之间的相对速度；A为迎风面面积；

为来流的单位矢量；

4)光压摄动力的计算，

通过软件设置面质比，计算光压摄动力的表达式为：

$\left\{\begin{array}{c}{\stackrel{\→}{F}}_n=\frac{\mathrm{IS}}{c}\cos \mathrm{\φ}\{\frac{2}{3}\mathrm{\ρ\μ}+[1+\mathrm{\ρ}(1-\mathrm{\μ})]\cos \mathrm{\φ}\}\\ {\stackrel{\→}{F}}_{\mathrm{\τ}}=\frac{\mathrm{IS}}{c}\cos \mathrm{\φ}\sin \mathrm{\φ}(1-\mathrm{\ρ}+\mathrm{\ρ\μ})\end{array}\right.----\left(10\right)$

式中，表示法向力，

表示切向力，φ为太阳入射角，I是太阳常数，c是光速，ρ表示反射系数，μ表示散射系数，S表示受力面积；

5)引力梯度力矩的计算，

通过下式计算引力梯度力矩：

$\stackrel{\→}{L}=\frac{3\mathrm{\μ}}{R^5}\widetilde{\mathrm{RI}}\stackrel{\→}{R}---\left(11\right)$

式中，μ为地球引力常数，

为星地矢量，

表示反对称矩阵；

6)气动力矩的计算，

利用式(9)的大气摄动力，结合摄动力作用力臂

求得气动干扰力矩

${\stackrel{\→}{T}}_g={\stackrel{\→}{F}}_a\×\stackrel{\→}{d}---\left(12\right)$

7)剩磁力矩的计算，

计算方法与式(2)类似，即：

${\stackrel{\→}{T}}_m={\stackrel{\→}{m}}_b\×\stackrel{\→}{B}---\left(13\right)$

式中，为星体剩磁矩；

8)光压力矩的计算，

利用式(10)的摄动力，结合对应的力臂，求得光压干扰力矩 ${\stackrel{\→}{T}}_l={\stackrel{\→}{F}}_n\×{\stackrel{\→}{d}}_n+{\stackrel{\→}{F}}_{\mathrm{\τ}}\×{\stackrel{\→}{d}}_{\mathrm{\τ}}$ (14)，式中，

和

分布表示法向力与切向力的力臂；

(五)利用上述已计算好的控制力、控制力矩、环境摄动力和环境干扰力矩，结合航天器力学模型和数值计算方法更新航天器的姿态轨道参数，其中：

1)航天器的轨道力学模型是：

$a=({\stackrel{\→}{F}}_c+{\stackrel{\→}{F}}_O+{\stackrel{\→}{F}}_{\mathrm{\ϵ}})/m---\left(15\right)$

式中，a为航天器的加速度，

为卫星所受的控制力，

为地球中心引力，

为各摄动力，m是航天器的质量；

2)航天器的姿态力学模型是：

$\stackrel{\stackrel{\·}{\→}}{h}+\stackrel{\·}{\mathrm{\ω}}\×\stackrel{\→}{h}={\stackrel{\→}{L}}_c+{\stackrel{\→}{L}}_{\mathrm{\ϵ}}---\left(16\right)$

$\stackrel{\·}{q}=0.5q\stackrel{\‾}{\mathrm{\ω}}---\left(17\right)$

式中，

表示航天器动量矩大小的导数，

为航天器动量矩，为控制力矩，

为干扰力矩，

为角速度，q为姿态四元数， $\stackrel{\‾}{\mathrm{\ω}}=\left[\begin{array}{cc}0& \stackrel{\→}{\mathrm{\ω}}\end{array}\right].$

轨道的数值方法和姿态的数值计算方法均采用RK；

(六)利用上述第(五)步骤更新后的航天器姿态轨道参数、第(三)步骤的太阳、月亮与当前时间信息，结合传感器模型计算传感器应该输出的物理量大小：

1)陀螺角速度测量值计算，

${\mathrm{\ω}}_{\mathrm{gc}}=(I-\mathrm{\Λ})(I-\mathrm{\Δ})Q_b^g{\mathrm{\ω}}_b-b-\mathrm{\η}---\left(18\right)$

式中，Q_b ^g表示体坐标系到陀螺坐标系的转换矩阵，b为随机游走误差，Λ为标定因子误差，Δ为非正交安装误差；

2)磁强计磁场强度测量值计算，

${\stackrel{\→}{B}}_k={(I+D)}^{-1}({\mathrm{\Θ}}^T{A}_k{\stackrel{\→}{M}}_k+b+\mathrm{\ϵ})---\left(19\right)$

式中，

为磁强计输出，D表示刻度因子误差及非正交化矩阵，Θ为非对准矩阵，

为某参考坐标系下的地磁矢量，A_k为参考坐标系到磁强计安装坐标系的转换矩阵，b为磁强计的常值偏置，

由第(三)步计算的磁场强度获得；

3)模拟太阳敏感器输出电压计算，

U＝U₀cosθ                 (20)

式中，θ表示入射角，U₀表示硅片正射时模拟太阳敏感器的输出电压，U表示当前的电压输出；

4)GPS定位数据计算，

x＝x_r+d+υ                      (21)

$\stackrel{\·}{d}=-\frac{1}{\mathrm{\τ}}d+\mathrm{\ϵ}---\left(22\right)$

式中，x_r是卫星真实的位置，x为GPS输出的定位数据，d为一阶马尔科夫过程，υ为噪声，τ为马尔科夫过程的相关时间，ε为驱动噪声。GPS定位时需要初始化，程序中应仿真该初始化的时间；

(七)通过输出接口电路输出各种数据，根据各传感器的实际输出方式，如果是模拟器件，则要将输出的物理量转换成电压值；如果是串口输出的设备，则要按固定的格式打包；最终将所有数据通过DA，串口和IO输出接口输出。

Images