一种日心轨道航天器太阳光压力矩确定方法

【技术领域】

本发明针对日心轨道航天器，基于太阳光压力模型，提出一种确定太阳光压力矩的新方法，本发明属于航天器总体设计技术领域。

背景技术：

随着航天任务需求的不断提升，对航天器的姿态精度有了更高的要求。对于深空探测任务的高轨道航天器来说，大部分都通过太阳帆板来为星上载荷提供能源，由于太阳帆板的存在，太阳光压力矩通常不能忽略，将严重影响航天器姿态精度的提升，因此有必要精确获得航天器受到的太阳光压力矩。

确定太阳光压力矩的关键在于实时判断航天器的遮挡情况，目前已经有不少的学者对其进行了研究。其中刘墩等人将卫星体和太阳帆板组成的几何体进行投影变换后，判断其相互遮挡并采用凸多边形求交的算法求取太阳光压有效作用面积，从而计算太阳光压力矩，但其计算过程相对复杂。郭建等人将空间站视为n个矩形平面的集合，再将每个矩形划分为正方形微元，通过判断投影后的微元与各个矩形平面的距离来判断遮挡，但由于投影后微元形状不确定，精确度有待考验。如何既简单高效又准确的确定航天器所受太阳光压力矩的大小仍是一大难题。

因此，本发明提出了一套算法较为简单，通用性高的光照遮挡判断流程，给出了一种光压力矩的确定方法。通过本发明提出的方法，可以根据面矢量法消除航天器的光照背面，进而对光照面进行遮挡判断，得到有效照射单元，最终根据光压力矩数学模型确定航天器所受光压力矩，进而得到光压力矩的解析表达式。

技术实现要素：

本发明针对日心轨道航天器，基于太阳光压力模型，提出一种太阳光压力矩的确定方法，从而得到航天器所受的太阳光压力矩。

针对上述问题，本发明技术方案如下：

根据航天器运动特性建立投影坐标系，得到航天器在投影坐标系下的分量列阵，从而根据面矢量法消除航天器的光照背面，进而对光照面进行遮挡判断，得到有效照射单元，最终根据光压力矩数学模型确定航天器所受光压力矩，进而得到光压力矩的解析表达式。具体操作步骤为：

步骤1：建立投影坐标系并获得航天器在投影坐标系下的分量列阵。具体包含如下步骤：

步骤1.1：定义坐标系

a.日心惯性坐标系fe(oexeyeze)

日心惯性坐标系的原点oe固连在太阳中心处，oexe轴在航天器轨道平面内并且指向某个恒星，oeze轴垂直于轨道平面，oeye在轨道平面内，并且和oexe轴、oeze轴构成直角坐标系。

b.日心轨道坐标系fo(ooxoyozo)

日心轨道坐标系的原点oo固连在航天器质心，oozo轴指向太阳中心，ooxo轴位于日心轨道平面内，垂直于oozo轴，并指向航天器运动方向，ooyo轴与ooxo轴、oozo轴构成直角坐标系。该坐标系随着航天器的轨道运动以角速度ωo绕ooyo轴负向旋转，ωo即为航天器的轨道角速度。

c.本体坐标系fb(obxbybzb)

本体坐标系与航天器固连，原点ob位于航天器质心，obxb轴指向航天器运动方向，obzb轴垂直于飞行轨道平面指向飞行器下方，obyb轴与obxb轴、obzb轴构成直角坐标系。

d.帆板固连坐标系fak(oakxakyakzak)

帆板固连坐标系的原点oak位于帆板与中心体连接处的中心，三轴方向与中心本体坐标系一致，且帆板可以绕oakyak轴旋转。

步骤1.2：坐标系间旋转矩阵计算

定义日心轨道坐标系fo到本体坐标系fb之间的旋转矩阵为rbo，本体坐标系fb到日心轨道坐标系fo之间的旋转矩阵为rob，定义本体坐标系fb到帆板固连坐标系fak之间的旋转矩阵为rakb，帆板固连坐标系fak到本体坐标系fb之间的旋转矩阵为rbak。航天器的姿态采用本体系与日心轨道系的相对方位来描述，旋转顺序采用3-1-2，θ,ψ分别表示航天器的滚动角、俯仰角以及偏航角，定义帆板绕oakyak轴旋转的角度为βk。从而：

rob＝r′bo,rbak＝r′akb(3)

步骤1.3：航天器微元的坐标转换

选取日心轨道坐标系为投影坐标系，太阳方向矢量(本体坐标系原点指向太阳中心的矢量)在投影坐标系下的分量列阵^os为：

^os＝[001]^t(4)

将航天器各个部件的质量微元的坐标数据转换到投影坐标系中，分为中心体和帆板两部分，表达式分别为：

^orbi＝rob^brbi(5)

^orbj＝rob(^brbak+rbak^akrakj)(6)

其中，^brbi和^orbi分别表示本体坐标系原点ob到中心体质量微元的矢量在本体坐标系和投影坐标系下的分量列阵，^akrakj表示帆板固连坐标系原点oak到帆板质量微元的矢量在帆板固连坐标系下的分量列阵，^orbj表示本体坐标系原点ob到帆板质量微元的矢量在投影坐标系下的分量列阵，^brbak表示本体坐标系原点ob到帆板附件坐标系原点oak的矢量在本体坐标系下的分量列阵。

步骤2：对航天器进行遮挡判断，得到有效照射单元。具体步骤如下：

步骤2.1：消除光照背面，得到受到照射的表面

利用有限元分析软件，将航天器表面划分为若干个三角形面片，根据三角形面片顶点在投影坐标系下的分量列阵求解面片的法向量在投影坐标系下的分量列阵^on0，表达式为：

其中，a、b、c分别表示三角形面片的三个顶点，^ora、^orb、^orc分别表示投影坐标系原点到三角形面片三个顶点的矢量在投影坐标系下的分量列阵，||r||表示r的二范数。

面片的外法向量(指向航天器外侧的法向量)在投影坐标系下的分量列阵^on通过三角形面片上的点进行修正，表达式为：

^on＝sign(^orp·^on0)^on0(8)

其中，^orp表示投影坐标系原点到三角形面片上p点的矢量在投影坐标系下的分量列阵，sign(*)为符号函数。

通过函数h(^os·^on)确定受到照射的表面：

其中，^onz表示^on的z轴分量，h(^os·^on)＝1表示面片受到照射，h(^os·^on)＝0表示面片被遮挡。

步骤2.2：航天器部件间的遮挡判断

将所有受到光照的面片投影到投影面上，投影面方程为：

z＝h(10)

对光照面的投影图形进行微元划分，并假设微元位于面片的投影中，则存在正数λ1,λ2,λ3使得：

其中，o表示投影坐标系的原点，p′表示微元的中心，a′、b′、c′分别表示三角形面片三个顶点a、b、c在投影面上对应的点。

在投影坐标系中计算对应的λ1,λ2,λ3，即：

其中，^orp′＝[xpyph]^t表示投影坐标系原点o到微元中心p′的矢量在投影坐标系下的分量列阵，^ora′＝[xayah]^t、^orb′＝[xbybh]^t、^orc′＝[xcych]^t分别表示投影坐标系原点o到投影面上a′、b′、c′的矢量在投影坐标系下的分量列阵。

当min(λ1,λ2,λ3)≥0时，表示微元位于三角形面片的投影内部，反之，微元位于三角形面片的投影外部。

假设微元在三角形面片abc上对应的点为p，则p点在投影坐标系下的z轴分量可通过下式计算得到：

微元可能不止为一个面片上的点的投影，其中z轴分量最大所对应的面片为该微元的照射面，该面片上对应的微元判断为光照面，而其他面片上对应的微元判断为遮挡。

步骤3：确定航天器受到的太阳光压力矩。

面片上微元da上的太阳辐射压力在投影坐标系下的分量列阵为：

其中，p≈4.56×10^-6n/m²表示辐射压强，da表示照射微元的面积，γ表示面片外法向量和太阳方向之间的夹角，ρa表示吸收比例系数，ρs表示镜面反射系数，ρd表示漫反射系数，且满足ρa+ρs+ρd＝1。

航天器所受太阳光压力矩在本体坐标系下的分量列阵为：

tsrp＝rbo∑^ori^×df(15)

其中，^ori表示投影坐标系原点到面片上微元的矢量在投影坐标系下的分量列阵，^ori^×表示^ori的叉乘矩阵，定义为：

步骤4：建立航天器所受太阳光压力矩的解析表达式。具体步骤如下：

步骤4.1：本发明的光压力矩解析表达式基于如下假设

假设1：航天器的结构为中心本体加对称分布的两个太阳帆板，且中心本体为均匀对称的六面体，其形心与质心重合。

假设2：航天器的姿态角处于较小范围内，航天器的中心本体与帆板之间的相互遮挡很小，可以忽略不计。

步骤4.2：建立航天器中心本体所受太阳光压力矩的解析表达式

基于步骤4.1的假设，航天器中心本体的形心和质心重合，有：

^br1＝^bn1＝[100]^t,^br2＝^bn2＝[010]^t,^br3＝^bn3＝[001]^t(16)

其中，^brk(k＝1,2,3)表示本体坐标系原点到中心本体面片质心的矢量在本体坐标系下的分量列阵，^bnk(k＝1,2,3)表示面片外法向量在本体坐标系下的分量列阵。

太阳方向在本体坐标系下的分量列阵^bs可表示为：

航天器中心本体所受光压力矩在本体坐标系下的解析表达式为：

因此，基于步骤4.1的假设，可以不考虑中心本体的太阳光压力矩。

步骤4.3：建立航天器帆板所受太阳光压力矩的解析表达式

基于步骤4.1的假设，航天器两个帆板对称分布，有：

^br1＝-^br2(19)

^bnk＝rbak^aknk＝rbak[001]^t＝[sinβk0cosβk]^t(20)

其中，^brk(k＝1,2)为本体坐标系原点到帆板质心的矢量在本体坐标系下的分量列阵，^bnk(k＝1,2)表示帆板面片外法向量在本体坐标系下的分量列阵。

航天器帆板所受光压力矩在本体坐标系下的解析表达式为：

分析帆板所受太阳光压力矩的解析表达式可知，当两个帆板的转角相等，即β1＝β2时，帆板产生的太阳光压力矩为零，即ts_srp＝0，当两个帆板的转角不相等，即β1≠β2时，安装在本体系yb轴的帆板会产生绕xb和zb轴的力矩。

本发明针对日心轨道上的航天器，基于太阳光压力模型，提出一种太阳光压力矩的确定方法，本发明主要优势体现在：

1)本发明提出的光压力矩确定方法，利用凸集性质来判断航天器各个部分之间是否相互遮挡，原理简单易懂，且计算精度较高。

2)本发明提出的光压力矩确定方法具有普适性，几乎所有的日心轨道航天器均能采用此方法来计算其受到的太阳光压力矩。

3)本发明针对中心本体加两个对称分布帆板的航天器构型，在航天器姿态小角度的情况下，给出了太阳光压力矩的解析表达式，使用该表达式计算太阳光压力矩精度基本不会下降，并且能有效地减小计算时间。

【附图说明】

图1为航天器结构示意图。

图2为各坐标系示意图。

图3为微元受光照示意图。

图4为太阳光压力矩计算流程图。

【具体实施方式】

下面结合附图1-4所示，以某型号的日心轨道的航天器为例，具体说明本发明的实施流程。航天器的参数如下：

航天器由中心刚体和两个对称分布的帆板构成，且中心刚体和帆板均为均质六面体。中心刚体的尺寸为50×25×20cm³，帆板的尺寸为80×25×1cm³，两个太阳帆板安装点在中心体机械坐标系下的坐标分别为(0,±21,0)cm。航天器的轨道为圆轨道，位于日心黄道平面内，航天器到太阳中心的距离为一个天文单位。中心刚体和帆板的吸收比例系数为ρa＝0.75，镜面反射系数为ρs＝0.25，漫反射系数为ρd＝0。航天器姿态角为θ＝10°,ψ＝9°，两个帆板的转角分别为β1＝5°，β2＝8°。如图4所示，下面介绍具体实施过程：

1、建立投影坐标系并获得航天器在投影坐标系下的分量列阵。具体包含如下步骤：

1.1定义坐标系：按照步骤1.1定义日心惯性坐标系、日心轨道坐标系、本体坐标系和帆板固连坐标系。

1.2坐标系间旋转矩阵计算

根据航天器的姿态角θ,ψ，帆板转角β1,β2，按照步骤1.2计算日心轨道坐标系fo、本体坐标系fb和帆板固连坐标系fak之间的旋转矩阵，即日心轨道坐标系fo到本体坐标系fb之间的旋转矩阵为rbo，本体坐标系fb到日心轨道坐标系fo之间的旋转矩阵为rob，定义本体坐标系fb到帆板固连坐标系fak之间的旋转矩阵为rakb，帆板固连坐标系fak到本体坐标系fb之间的旋转矩阵为rbak，计算结果为

rob＝r′bo,rbak＝r′akb(24)

1.3航天器微元的坐标转换

选取日心轨道坐标系为投影坐标系，太阳方向矢量(本体坐标系原点指向太阳中心的矢量)在投影坐标系下的分量列阵^os为：

^os＝[001]^t(25)

根据航天器的结构和尺寸，得到各个部件的质量微元的坐标数据，并将其转换到投影坐标系中，分为中心体和帆板两部分，表达式分别为：

^orbi＝rob^brbi(26)

^orbj＝rob(^brbak+rbak^akrakj)(27)

其中，^brbi和^orbi分别表示本体坐标系原点ob到中心体质量微元的矢量在本体坐标系和投影坐标系下的分量列阵，^akrakj表示帆板固连坐标系原点oak到帆板质量微元的矢量在帆板固连坐标系下的分量列阵，^orbj表示本体坐标系原点ob到帆板质量微元的矢量在投影坐标系下的分量列阵，^brbak表示本体坐标系原点ob到帆板附件坐标系原点oak的矢量在本体坐标系下的分量列阵。

2、对航天器进行遮挡判断，得到有效照射单元。具体步骤如下：

2.1消除光照背面，得到受到照射的表面

利用有限元分析软件，将航天器表面划分为若干个三角形面片，根据三角形面片顶点在投影坐标系下的分量列阵求解面片的法向量在投影坐标系下的分量列阵^on0，表达式为：

其中，a、b、c分别表示三角形面片的三个顶点，^ora、^orb、^orc分别表示投影坐标系原点到三角形面片三个顶点的矢量在投影坐标系下的分量列阵，||r||表示r的二范数。

面片的外法向量(指向航天器外侧的法向量)在投影坐标系下的分量列阵^on通过三角形面片上的点进行修正，表达式为：

^on＝sign(^orp·^on0)^on0(29)

其中，^orp表示投影坐标系原点到三角形面片上p点的矢量在投影坐标系下的分量列阵，sign(*)为符号函数。

通过函数h(^os·^on)确定受到照射的表面：

其中，^onz表示^on的z轴分量，h(^os·^on)＝1表示面片受到照射，h(^os·^on)＝0表示面片被遮挡。

2.2航天器部件间的遮挡判断

将所有受到光照的面片投影到投影面上，投影面方程为：

z＝h＝1m(31)

对光照面的投影图形进行微元划分，并假设微元位于面片的投影中，则存在正数λ1,λ2,λ3使得：

其中，o表示投影坐标系的原点，p′表示微元的中心，a′、b′、c′分别表示三角形面片三个顶点a、b、c在投影面上对应的点。

在投影坐标系中计算对应的λ1,λ2,λ3，即：

其中，^orp′＝[xpyph]^t表示投影坐标系原点o到微元中心p′的矢量在投影坐标系下的分量列阵，^ora′＝[xayah]^t、^orb′＝[xbybh]^t、^orc′＝[xcych]^t分别表示投影坐标系原点o到投影面上a′、b′、c′的矢量在投影坐标系下的分量列阵。

保留满足不等式min(λ1,λ2,λ3)≥0的面片，其可能为微元受到光照的平面。

假设微元在三角形面片abc上对应的点为p，则p点在投影坐标系下的z轴分量可通过下式计算得到：

其中z轴分量最大所对应的面片为该微元的照射面，该面片上对应的微元判断为光照面，而其他面片上对应的微元判断为遮挡。

3、确定航天器受到的太阳光压力矩。

面片上微元da上的太阳辐射压力在投影坐标系下的分量列阵为：

其中，p≈4.56×10^-6n/m²表示辐射压强，da表示照射微元的面积，γ表示面片外法向量和太阳方向之间的夹角，ρa表示吸收比例系数，ρs表示镜面反射系数，ρd表示漫反射系数，且满足ρa+ρs+ρd＝1。

航天器所受太阳光压力矩在本体坐标系下的分量列阵为：

其中，^ori表示投影坐标系原点到面片上微元的矢量在投影坐标系下的分量列阵，^ori^×表示^ori的叉乘矩阵，定义为

4、建立航天器所受太阳光压力矩的解析表达式。具体步骤如下：

4.1本发明的光压力矩解析表达式基于如下假设

假设1：航天器的结构为中心本体加对称分布的两个太阳帆板，且中心本体为均匀对称的六面体，其形心与质心重合。

假设2：航天器的姿态角处于较小范围内，航天器的中心本体与帆板之间的相互遮挡很小，可以忽略不计。

4.2建立航天器中心本体所受太阳光压力矩的解析表达式

将航天器参数带入步骤4.2中的解析表达式中进行计算，航天器中心本体所受光压力矩在本体坐标系下的解析表达式及具体大小为：

4.3建立航天器帆板所受太阳光压力矩的解析表达式

将航天器参数带入步骤4.3中的解析表达式中进行计算，航天器帆板所受光压力矩在本体坐标系下的解析表达式及具体大小为：

经过上述步骤能够求解得到航天器所受太阳光压力矩的具体大小。对比步骤3和解析表达式计算得到的太阳光压力矩可知，使用解析表达式计算绕xb和zb轴的光压力矩时，误差在7％以内，而绕yb轴的光压力矩比绕xb和zb轴的光压力矩小4个数量级，解析表达式中将其视为零也是合理的，表明得到的太阳光压力矩解析表达式具有一定的精度。