导航GEO卫星太阳辐射压力矩计算方法

技术领域

本发明属于航天飞行器设计领域，涉及导航GEO卫星太阳辐射压力矩计 算方法。

背景技术

由于卫星导航具有覆盖区域大、作用距离远以及受自然灾害破坏小等优 点，因此卫星导航得到了广泛应用。由于GEO卫星能覆盖固定区域，因此我 国在区域星座中得到了应用。太阳辐射压力是导航GEO卫星定点以后长期工 作期间的主要干扰力，其作用力使卫星轨道偏心率呈长周期变化。当太阳辐 射压力相对应的压心和卫星质心不一致时，会产生太阳辐射压力力矩，这是 影响卫星姿态运动的最主要干扰力矩。太阳辐射压力力矩与卫星的形状、表 面材料的光学特性，部件之间的遮挡以及卫星位置等有关。

导航GEO卫星因其使命特点，采用动量轮控制卫星姿态，动量轮会不断 吸收干扰力矩，到一定程度达到转速极限就会饱和，从而丧失姿态控制能力。 因此需要定期对动量轮进行卸载。为了使得动量轮卸载周期满足用户要求， 卫星设计之初，需要准确知道太阳干扰力矩的大小。

从以上分析可知，有必要提供导航GEO卫星太阳辐射压力矩计算方法来 支持导航卫星总体设计，以支导航星座组网建设工作。

一代导航GEO卫星太阳辐射压计算时考虑了卫星的形状、表面材料的光 学特性，部件之间的遮挡等因素，但没有考虑S/L天线向北偏置的角度。本 算法不但考虑了卫星的形状、表面材料的光学特性，部件之间的遮挡等因素， 同时考虑了S/L天线向北偏置的角度(即天线反射面按右手系绕星体X轴转 过的角度)，也针对网孔状S/L天线建立模型，计算结果更接近实际。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提出导航GEO卫星太 阳辐射压力矩计算方法，该方法针对导航GEO卫星的形状、表面材料的光学 特性，考虑了部件之间的遮挡等因素并考虑了S/L天线向北偏置的角度，同 时针对网孔状S/L天线建立模型，对导航GEO卫星太阳辐射压力矩进行了计 算，保证了导航GEO卫星太阳辐射压力矩计算结果的精度。

为解决上述技术问题，本发明的技术方案，通过以下过程实现：

(1)获得输入条件，所述输入条件包括：卫星构型；卫星几何尺寸； 卫星表面特性参数；卫星质心位置；太阳翼位置、尺寸和表面特性参数；S/L 天线位置几何尺寸和表面特性参数；

(2)根据卫星构型、几何尺寸、卫星质心以及卫星表面特性参数计算 太阳对星本体辐射压力矩；星本体辐射压力矩包括星本体各个面板的辐射压 力矩；所述各面板包括西面板、东面板、对地板、背地板、南面板和北面板；

(3)根据卫星质心、太阳翼位置、太阳翼尺寸和表面特性参数计算太 阳对太阳翼辐射压力矩；

(4)根据卫星质心、S/L天线位置、S/L天线几何尺寸和表面特性参数 计算太阳对S/L天线辐射压力矩，具体过程如下：

对所述S/L天线进行微元划分；

判断所述S/L天线受东面板的遮挡：

$\left|\mathrm{AC}\right|=\frac{l_3-x_{\mathrm{ij}}}{\cos \delta \sin \alpha };$

y_a＝y_ij+|AC|sinδ；

z_a＝z_ij+|AC|cosδcosα，

若满足：-l₂＜y_a＜l₁，-l₅＜z_a＜l₆则存在遮挡，否则无遮挡，

式中：

|AC|是天线上某微元到东板上某微元的距离；

x_ij是天线上某微元X坐标；

y_ij是天线上某微元Y坐标；

z_ij是天线上某微元Z坐标；

y_a是东板上某微元的Y坐标；

z_a是东板上某微元的Z坐标；

力矩计算：

对于天线正面：

存在遮挡时，M_x8＝0，M_y8＝0，M_z8＝0；

无遮挡时，光压力矩：$M_{x8}=\underset{j=1}{\overset{n}{\Sigma }}\underset{i=1}{\overset{2j}{\Sigma }}{M}_{x(i,j)};$$M_{y8}=\underset{j=1}{\overset{n}{\Sigma }}\underset{i=1}{\overset{2j}{\Sigma }}{M}_{y(i,j)};$$M_{z8}=\underset{j=1}{\overset{n}{\Sigma }}\underset{i=1}{\overset{2j}{\Sigma }}{M}_{z(i,j)},$

式中：

M_x8、M_y8、M_z8分别是天线正面X、Y、Z轴太阳辐射压力矩；

M_x(i，j)、M_y(i，j)、M_z(i，j)分别是微元X、Y、Z轴太阳辐射压力矩，按下式计 算：

M_x(i，j)＝[(P+τcosβ)(y_ijcosβ_Ncosβ′+z_ijsinβ_Ncosβ′)-τ(y_ijcosδcosα-z_ijsinδ)]K_dΔs；

M_y(i，j)＝[-(z_ijsinβ′+x_ijcosβ′cosβ_N)(P+τcosβ)-τ(z_ijcosδsinα-x_ijcosδcosα)]K_dΔs；

M_z(i，j)＝[(P+τcosβ)(y_ijsinβ′-x_ijsinβ_Ncosβ′)-τ(x_ijsinδ-y_ijcosδsinα)]K_dΔs，

式中：

Δs-微元的面积；

A₁、A₂、A₃分别是微元的X、Y、Z坐标；

β′是天线平面与星体系yox平面的夹角；

β_N是天线向北偏置的角度；

K_d是网状天线面积比例因子；

对于天线背面：

天线背面三轴力矩分别为：

M_x9＝[-(P+τcosβ)(a₂cosβ′cosβ_N+a₃cosβ′sinβ_N)-τ(a₂cosδcosα-a₃sinδ)]K_dπR²；

M_y9＝πR²K_d[(P+τcosβ)(a₃sinβ′+a₁cosβ′cosβ_N)+τcosδ(a₁cosα-a₃sinα)]；

M_z9＝[(P+τcosβ)(a₁sinβ_Ncosβ′-a₂sinβ′)-τ(a₁sinδ-a₂cosδsinα)]K_dπR²；

(5)对所述星本体辐射压力矩、所述太阳翼辐射压力矩和所述天线辐 射压力矩进行傅立叶拟合得到太阳辐射压力矩的傅立叶拟合系数。

利用上述太阳辐射压力矩的傅立叶拟合系数进行卫星总体设计。

进一步的，在上述步骤(2)中，采用微元分割法计算所述东面板的辐 射压力矩，具体步骤为：

以所述东面板一角为坐标原点，对所述东面板进行微元划分；

根据所述微元确定所述东面板受所述S/L天线的遮挡情况；

分别计算所述东面板的有遮挡力矩和无遮挡的光压力矩。

进一步的，在上述步骤(5)中，所述傅立叶拟合以太阳在形体坐标系 中的位置为自变量，将所述太阳辐射压力矩表示成傅立叶级数形式，其中包 含长期项和以天为周期的项。

本发明与现有技术相比的有益效果是：

本发明分别星本体辐射压力矩、太阳翼辐射压力矩和S/L天线辐射压力 矩分别进行计算，并最终进行傅里叶拟合，从而使得获得的太阳辐射压力矩 的傅立叶拟合系数可更精确的指导卫星的总体设计。

本发明针对S/L天线的网状特点及于东面板相互遮挡的情况，采用微元划分 方法分别计算S/L天线的正反两面的力矩，且在整个步骤中考虑了S/L天线向北 偏置的角度，从而可保证获得更加精确的S/L天线的辐射压力矩。

进一步的，针对受S/L天线遮挡的东面板，采用微元划分方法对东面板进行 划分，从而可对受遮挡部分和不受遮挡部分分别计算力矩，因此，同样保证了结 果的精度。本发明考虑了S/L网状天线的独特构型，并建立了网孔状S/L天线的 辐射压力矩模型。

附图说明

图1为本发明流程图。

具体实施方式

参考图1所示流程图，本发明方法的具体步骤如下：

第一步获得输入条件。计算所需输入条件主要包括：卫星构型、几何尺 寸和卫星表面特性参数；卫星质心位置；太阳翼位置、尺寸和表面特性参数； S/L天线位置几何尺寸和表面特性参数，为计算太阳辐射压力矩做准备；

第二步根据卫星构型、几何尺寸、卫星质心以及卫星表面特性参数计算 太阳对星本体各个面板的辐射压力矩，星体东面板与S/L天线存在遮挡，为 了考虑遮挡的影响，采用微元分割法计算西面板辐射压力矩。

针对各面板的具体计算方法如下：

(1).西面板

X轴力矩：M_x1＝-τ(J₂cosδcosα-J₁sinδ)

Y轴力矩：M_y1＝-PJ₁-τl₄J₃cosδcosα

Z轴力矩：M_z1＝J₂P+τl₄J₃sinδ

式中：

M_x1、M_y1、M_z1分别是西面板X、Y、Z轴太阳辐射压力矩，X、Y、Z轴分 别为卫星本体坐标系的三个方向；

δ为太阳光矢量和星体系XOZ面的夹角；

$P=-k\cos \beta \{\frac{2\mathrm{\rho A}}{3}+\cos \beta [1+(1-A)\rho ]+\frac{2(1-\rho )\kappa }{3}\},$是法向调节常数；

τ＝kcosβ[1-ρ(1-B)]；

β-面元法向和太阳光方向间的夹角；

ρ-反射系数；

τ-切向系数；

A，B-法向、切向调节系数；

κ-温度调节系数；

l₄-卫星质心到西面板距离；

$J_1=(l_1+l_2)(l_6^5-l_5^2)/2,$计算常数

$J_2=(l_5+l_6)(l_1^2-l_2^2)/2,$计算常数；

J₃＝(l₁+l₂)(l₅+l₆)，计算常数；

l₁-卫星质心到南面板距离；

l₂-卫星质心到北面板距离；

l₅-卫星质心到背地板距离；

l₆-卫星质心到对地板距离；

α-太阳在星体在XOZ平面内投影与星体Z轴夹角。

(2).东面板

由于东面板和天线相互遮挡，因此采用微元分割法计算，具体方法如下：

A、微元划分

选东面板一角为坐标原点，该点的坐标为(l₃，-l₂，-l₅)

将Y方向分为m等份，Z方向分为n等份，每份的长度为：

$\mathrm{\Delta Y}=\frac{l_1+l_2}{m},$$\mathrm{\Delta Z}=\frac{l_5+l_6}{n}$

取i＝1，2，3，...，n，j＝1，2，3，...，m；

第(i，j)微元中心坐标为：(l₃，-l₂+(j-1/2)ΔY，-l₅+(i-1/2)ΔZ)

式中：

l₃是卫星质心到东面板距离；

ΔY是微元Y轴方向长度；

ΔZ是微元Z轴方向长度；

n是东面板Z方向等份数；

m是东面板Y方向等份数。

B、微元被天线遮挡判断

$l=\frac{(a_1-A_1)N_x+(a_2-A_2)N_y+(a_3-A_3)N_z}{N_x\cos \delta \sin \alpha +N_y\sin \delta +N_z\cos \delta \cos \alpha }$

x＝lcosδsinα+A₁

y＝lsinδ+A₂

z＝lcosδcosα+A₃

|A′O′|＝[(x-a₁)²+(y-a₂)²+(z-a₃)²]^1/2，

当|A′O′|≥R时不遮挡，|A′O′|＜R时遮挡，

式中：

l是比例因子；

(x，y，z)是天线反射面上任一点坐标；

(a₁，a₂，a₃)为天线中心点的坐标；

R是天线口径；

|A′O′|是天线反射面上任一点到天线中心点距离；

(N_x，N_y，N_z)为天线法向坐标。

C、力矩计算

存在遮挡时，M_x2＝0，M_y2＝0，M_z2＝0；

无遮挡时，东面板光压力矩：$M_{x2}=\underset{i=1}{\overset{n}{\Sigma }}\underset{j=1}{\overset{m}{\Sigma }}{M}_{x(i,j)},$$M_{y2}=\underset{i=1}{\overset{n}{\Sigma }}\underset{j=1}{\overset{m}{\Sigma }}{M}_{y(i,j)},$$M_{z2}=\underset{i=1}{\overset{n}{\Sigma }}\underset{j=1}{\overset{m}{\Sigma }}{M}_{z(i,j)};$

式中：

M_x2、M_y2、M_z2分别是东面板X、Y、Z轴太阳辐射压力矩；

M_x(i，j)、M_y(i，j)、M_z(i，j)分别是微元X、Y、Z轴太阳辐射压力矩，按下式计 算：

M_x(i，j)＝-τ(A₂cosδcosα-A₃sinδ)ΔS

M_y(i，j)＝(PA₃+τA₁cosδcosα)ΔS

M_z(i，j)＝(-PA₂-τA₁sinδ)ΔS

式中：

Δs-微元的面积；

A₁、A₂、A₃分别是微元的X、Y、Z坐标。

(3).对地板

X轴力矩：M_x3＝PW₁+τl₆W₃sinδ；

Y轴力矩：M_y3＝-PW₂-τl₆W₃cosδsinα；

Z轴力矩：M_z3＝τ(W₁cosδsinα-W₂sinδ)，

式中：

M_x3、M_y3、M_z3分别是对地板X、Y、Z轴太阳辐射压力矩；

$W_1=(l_3+l_4)(l_1^2-l_2^2)/2,$计算常数；

$W_2=(l_1+l_2)(l_3^2-l_4^2)/2,$计算常数；

W₃＝(l₃+l₄)(l₁+l₂)，计算常数。

(4).背地板

X轴力矩：M_x4＝-PW₁-τl₅W₃sinδ；

Y轴力矩：M_y4＝PW₂+τl₅W₃cosδsinα；

Z轴力矩：M_z4＝τ(W₁cosδsinα-W₂sinδ)，

式中：

M_x4、M_y4、M_z4分别是背地板X、Y、Z轴太阳辐射压力矩。

(5).南面板

X轴力矩：M_x5＝-PQ₁-τl₁Q₃cosδcosα；

Y轴力矩：M_y5＝-τ(Q₁cosδsinα-Q₂cosδcosα)；

Z轴力矩：M_z5＝PQ₂+τl₁Q₃cosδsinα，

式中：

M_x5、M_y5、M_z5分别是南面板X、Y、Z轴太阳辐射压力矩；

$Q_1=(l_3+l_4)(l_6^2-l_5^2)/2,$计算常数；

$Q_2=(l_5+l_6)(l_3^2-l_4^2)/2,$计算常数；

Q₃＝(l₃+l₄)(l₅+l₆)，计算常数。

(6).北面板

X轴力矩：M_x6＝PQ₁+τl₂Q₃cosδcosα；

Y轴力矩：M_y6＝-τ(Q₁cosδsinα-Q₂cosδcosα)；

Z轴力矩：M_z6＝-PQ₂-τl₂Q₃cosδsinα，

式中：

M_x6、M_y6、M_z6分别是北面板X、Y、Z轴太阳辐射压力矩。

第三步根据卫星质心、太阳翼位置、太阳翼尺寸和表面特性参数计算太 阳对太阳翼辐射压力矩。

X轴力矩：M_x7＝PI₁cosα+τI₃D₂sinδ；

Y轴力矩：M_y7＝I₃[P(D₂sinα-D₁cosα)]；

Z轴力矩：M_z7＝-PI₁sinα-τI₃D₁sinδ，

式中：

M_x7、M_y7、M_z7分别是太阳翼X、Y、Z轴太阳辐射压力矩；

$I_1=L\frac{G_1^2-{(H+G_1)}^2+{(H+G_2)}^2-G_2^2}{2},$计算常数；

I₃＝2HL，计算常数；

D₁、D₂分别是太阳翼转轴中心在星体XOZ面投影点的X、Z坐标；

L是太阳帆板宽；

H是太阳帆板长。

第四步根据卫星质心、S/L天线位置、S/L天线几何尺寸和表面特性参 数计算太阳对S/L天线辐射压力矩，采用微元分割法分别针对天线正面和反 面进行计算。

天线正面可能被东面板遮挡，用微元法求解，计算每一个微元力矩，然 后叠加。具体方法如下：

A、微元的划分

将圆分为n个同心圆，最里面的圆分为两个半圆，以后从里到外按倍数 增加。

B、遮挡判断

$\left|\mathrm{AC}\right|=\frac{l_3-x_{\mathrm{ij}}}{\cos \delta \sin \alpha };$

y_a＝y_ij+|AC|sinδ；

z_a＝z_ij+|AC|cosδcosα，

若满足：-l₂＜y_a＜l₁，-l₅＜z_a＜l₆则存在遮挡，否则无遮挡，

式中：

|AC|是天线上某微元到东板上某微元的距离；

x_ij是天线上某微元X坐标；

y_ij是天线上某微元Y坐标；

z_ij是天线上某微元Z坐标；

y_a是东板上某微元的Y坐标；

z_a是东板上某微元的Z坐标。

C、力矩计算

对于天线正面：

存在遮挡时，M_x8＝0，M_y8＝0，M_z8＝0；

无遮挡时，光压力矩：$M_{x8}=\underset{j=1}{\overset{n}{\Sigma }}\underset{i=1}{\overset{2j}{\Sigma }}{M}_{x(i,j)};$$M_{y8}=\underset{j=1}{\overset{n}{\Sigma }}\underset{i=1}{\overset{2j}{\Sigma }}{M}_{y(i,j)};$$M_{z8}=\underset{j=1}{\overset{n}{\Sigma }}\underset{i=1}{\overset{2j}{\Sigma }}{M}_{z(i,j)},$

式中：

M_x8、M_y8、M_z8分别是天线正面X、Y、Z轴太阳辐射压力矩；

M_x(i，j)、M_y(i，j)、M_z(i，j)分别是微元X、Y、Z轴太阳辐射压力矩，按下式计 算：

M_x(i，j)＝[(P+τcosβ)(y_ijcosβ_Ncosβ′+z_ijsinβ_Ncosβ′)-τ(y_ijcosδcosα-z_ijsinδ)]K_dΔs；

M_y(i，j)＝[-(z_ijsinβ′+x_ijcosβ′cosβ_N)(P+τcosβ)-τ(z_ijcosδsinα-x_ijcosδcosα)]K_dΔs；

M_z(i，j)＝[(P+τcosβ)(y_ijsinβ′-x_ijsinβ_Ncosβ′)-τ(x_ijsinδ-y_ijcosδsinα)]K_dΔs，

式中：

Δs-微元的面积；

A₁、A₂、A₃分别是微元的X、Y、Z坐标；

β′是天线平面与星体系yox平面的夹角；

β_N是天线向北偏置的角度；

K_d是网状天线面积比例因子。

对于天线背面：

天线背面三轴力矩分别为：

M_x9＝[-(P+τcosβ)(a₂cosβ′cosβ_N+a₃cosβ′sinβ_N)-τ(a₂cosδcosα-a₃sinδ)]K_dπR²；

M_y9＝πR²K_d[(P+τcosβ)(a₃sinβ′+a₁cosβ′cosβ_N)+τcosδ(a₁cosα-a₃sinα)]；

M_z9＝[(P+τcosβ)(a₁sinβ_N cosβ′-a₂sinβ′)-τ(a₁sinδ-a₂cosδsinα)]K_dπR²。

第五步根据前面计算结果得到最终太阳辐射压力矩，经过傅立叶拟合得 到太阳辐射压力矩傅立叶拟合系数。

从前面计算可得到卫星不同轨道位置的三轴干扰力矩为：

M_x＝M_x1+M_x2+M_x3+M_x4+M_x5+M_x6+M_x7+M_x8+M_x9

M_y＝M_y1+M_y2+M_y3+M_y4+M_y5+M_y6+M_y7+M_y8+M_y9

M_z＝M_z1+M_z2+M_z3+M_z4+M_z5+M_z6+M_z7+M_z8+M_z9

经过傅立叶拟合，卫星三轴太阳辐射压力矩的傅立叶级数形式如下：

$M_x=a_{x0}+\underset{k=1}{\overset{N}{\Sigma }}[a_{\mathrm{xk}}\cos \left(\mathrm{k\theta }\right)+b_{\mathrm{xk}}\sin \left(\mathrm{k\theta }\right)]$

$M_y=a_{y0}+\underset{k=1}{\overset{N}{\Sigma }}[a_{\mathrm{yk}}\cos \left(\mathrm{k\theta }\right)+b_{\mathrm{yk}}\sin \left(\mathrm{k\theta }\right)]$

$M_z=a_{z0}+\underset{k=1}{\overset{N}{\Sigma }}[a_{\mathrm{zk}}\cos \left(\mathrm{k\theta }\right)+b_{\mathrm{zk}}\sin \left(\mathrm{k\theta }\right)]$

式中：

N是傅立叶拟合阶数；

θ是太阳在ZOX平面内投影与Z轴的夹角，以Y轴为右手系度量，角度 从0～360；

a_x0是X轴力矩常值系数；

a_y0是Y轴力矩常值系数；

a_z0是Z轴力矩常值系数；

a_xk、b_xk、a_yk、b_yk、a_zk、b_zk是周期性力矩系数。

利用上述获得的结果制导卫星的总体设计，从而为获得更优的卫星总体 设计精度。我国某型号卫星用了本发明计算了导航GEO卫星太阳辐射压力 矩，经过卫星在轨飞行验证，导航GEO卫星太阳辐射压力矩计算结果满足实 际工程需求。

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知 技术。