一种基于二次反射的高精度卫星太阳光压摄动力建模方法

摘要

本发明公开了一种基于二次反射的高精度卫星太阳光压摄动力建模方法，采用模拟太阳光源照射卫星表面部件，研究光线入射后反射再次入射到表面部件后，表面部件的受到的光压摄动力，最终得到太阳光压摄动加速度模型，该建模方法能够建立复杂卫星包含二次反射的高精度太阳光压摄动模型，作为中高轨卫星最主要的非保守力，能够完善并提高高精度卫星动力学模型，进一步提高精密定轨、轨道预报精度。

说明书

技术领域

本发明属于飞行器设计领域，涉及一种基于二次反射的高精度卫星太阳光压摄动力建模方法。

背景技术

随着高精度卫星精密定轨与轨道预报精度要求提高，尽可能精确的卫星动力学模型成为工程需要。以导航卫星为例，光压摄动作为目前中高轨卫星最主要的误差源，其精确度直接影响卫星的精密定轨与轨道预报。太阳光压是光压摄动最主要的摄动源，因此，考虑二次反射的高精度卫星太阳光压摄动力建模的研究是必要的。

发明内容

有鉴于此，本发明提供了一种基于二次反射的高精度卫星太阳光压摄动力建模方法，能够完善并提高高精度卫星动力学模型，进一步提高精密定轨、轨道预报精度。

本发明的一种基于二次反射的卫星太阳光压摄动力建模方法，包括如下步骤：

步骤1、获取待建模卫星的表面三维模型，确定卫星的所有表面部件，将每个表面部件的表面划分成多个面元；

步骤2、用面积等于卫星外形最大包络的方形像元阵列作为模拟太阳光源，每个像元作为一个小光源并发出一条垂直于方形像元阵列的太阳光线，所述模拟太阳光源与卫星的距离为一个天文单位；

步骤3、针对所述卫星的各个表面部件，记录所述模拟太阳光源的每个像元发出的光线与表面部件最短距离的交点和相应的入射角，统计第i个表面部件上有所述交点的面元数量，即为初次照射第i个表面部件的有效被照射面元数量A_ieff1；i＝1,2,...N，其中N为表面部件数量；

根据有效被照射面元数量A_ieff1，得到第i个表面部件由模拟太阳光源初次照射所受的法向光压摄动力f_ni1以及切向光压摄动力f_si1；

步骤4、首先，针对各个表面部件，根据步骤3获得的所有入射光线的交点和入射角，继而确定入射光线在所在面元上的反射光线；

其次，统计第i个表面部件与所述反射光线有最短距离交点的面元个数A_ieff2，则第i个表面部件被二次照射中漫反射的有效面元数为A_ieff2·v(1-μ)，镜面反射的有效面元数为A_ieff2·vμ；

然后，根据A_ieff2·v(1-μ)，获得第i个表面部件在二次照射中由漫反射引起的法向光压摄动力f_ni2d；根据A_ieff2·vμ，获得第i个表面部件在二次照射中由镜面反射引起的法向光压摄动力f_ni2m和切向光压摄动力f_si2m；

最后，得到第i个部件所受的总的法向光压摄动力f_ni为：

f_ni＝f_ni1+f_ni2d+f_ni2m；

总的切向光压摄动力f_si为：

f_si＝f_si1+f_si2m；

步骤5、根据各个表面部件的安装位置和法线方向，将所有表面部件的法向光压摄动力f_ni和切向光压摄动力f_si分别分解到卫星本体坐标系的三个轴方向，并得到各轴方向的合力，最后得到各轴方向的光压摄动加速度；

步骤6、采用步骤2-5的方法，在一个卫星运动周期内，以设定角度为步长，计算设定时间内每个卫星运动周期不同时刻的太阳光压摄动加速度；针对计算结果，通过傅里叶多项式拟合得到太阳光压摄动加速度的数学模型。

较佳的，所述步骤6中的设定时间为半年。

较佳的，所述像元大小为1mm×1mm。

较佳的，所述面元大小为1mm×1mm。

较佳的，所述初次照射法向光压摄动力f_ni1为：

$-\tau ·\underset{j=1}{\overset{A_{ieff1}}{\Sigma }}\left(E_{sun}{\mathrm{cos\theta }}_{ij}/c\right)\left[\left(1+{\nu }_i{\mu }_i\right){\mathrm{cos\theta }}_{ij}+\frac{2}{3}{\nu }_i\left(1-{\mu }_i\right)\right]$

切向光压摄动力f_si1为：

$\tau \underset{j=1}{\overset{A_{ieff1}}{\Sigma }}\left(E_{sun}{\mathrm{cos\theta }}_{ij}/c\right)\left(1-{\nu }_i{\mu }_i\right){\mathrm{sin\theta }}_{ij}$

其中，j＝1,2,...,M，M为表面部件上面元数量，τ为日蚀因子，c为光速，E_sun＝1367W/m²为太阳辐射强度；θ_ij为第i个表面部件的第j个面元的光线入射角；ν_i表示第i个表面部件的反射率、μ_i第i个表面部件的镜面系数。

较佳的，二次照射中，第i个表面部件由漫反射引起的法向光压摄动力f_ni2d为：

$f_{ni2d}=-\tau \underset{j=1}{\overset{A_{ieff2}·v\left(1-\mu \right)}{\Sigma }}\left(E_{sun}{\mathrm{cos\theta }}_{ij}/c\right)\left[\left(1+{\nu }_i{\mu }_i\right){\mathrm{cos\theta }}_{ij}+\frac{2}{3}{\nu }_i\left(1-{\mu }_i\right)\right]$

第i个表面部件由镜面反射引起的法向光压摄动力f_ni2m和切向光压摄动力f_si2m为：

$f_{ni2m}=-\tau \underset{i=1}{\overset{A_{ieff2}·v\mu }{\Sigma }}\left(E_{sun}{\mathrm{cos\theta }}_{ij}/c\right)\left[\left(1+{\nu }_i{\mu }_i\right){\mathrm{cos\theta }}_{ij}+\frac{2}{3}{\nu }_i\left(1-{\mu }_i\right)\right]$

$f_{si2m}=\tau \underset{j=1}{\overset{A_{ieff2}·v\mu }{\Sigma }}\left(E_{sun}{\mathrm{cos\theta }}_{ij}/c\right)\left(1-{\nu }_i{\mu }_i\right){\mathrm{sin\theta }}_{ij}$

其中，j＝1,2,...,M，M为表面部件上面元数量，τ为日蚀因子，c为光速，E_sun＝1367W/m²为太阳辐射强度；θ_ij为第i个表面部件的第j个面元的光线入射角；ν_i表示第i个表面部件的反射率、μ_i第i个表面部件的镜面系数。

较佳的，所述步骤6中所述设定角度为10°。

本发明具有如下有益效果：

(1)本发明的建模方法能够建立复杂卫星包含二次反射的高精度太阳光压摄动模型，作为中高轨卫星最主要的非保守力，能够完善并提高高精度卫星动力学模型，进一步提高精密定轨、轨道预报精度。

附图说明

图1为卫星本体坐标系示意图；

图2为卫星轨道弧段示意图；

图3为太阳光线像元阵列示意图；

图4为本发明的建模方法流程图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

如图4所示，本发明的建模方法的具体步骤如下：

1)选取卫星本体坐标系作为太阳光压摄动建模的基准坐标系，如图1所示，以卫星的中心的坐标原点，沿太阳翼方向为Y轴，X轴垂直Y同时垂直卫星的一个表面，Z轴同时垂直于Y轴、X轴以及卫星的一个表面；

2)基于卫星本体坐标系，获取待建模研究卫星的表面三维模型，对卫星所有表面部件编号记为i，i为1，…，N，建立包括每个卫星表面部件的卫星表面参数数据库，每个表面部件的信息包括：部件的形状(例如矩形、圆柱、抛物面等)、部件相对本体系的安装位置(r_xi,r_yi,r_zi)^T、表面积A_i、反射率v_i、镜面系数μ_i、表面部件的外法线在本体系中的单位矢量P_i＝(x_i,y_i,z_i)^T；

3)卫星姿态控制模式和太阳矢量运动规律分析。设太阳在J2000惯性系中的矢量为卫星的轨道六根数中升交点经度Ω、轨道倾角i、升交点幅角u，J2000惯性系下的太阳矢量依次经过3-1-2三次坐标转换，可以转到轨道坐标系，旋转矩阵为：

$R_{oJ2000}=R_y(-u-\frac{\pi }{2})·R_x(i-\frac{\pi }{2})·R_z\left(\Omega \right)$

再设卫星按321转序姿态角为：滚动角俯仰角θ、偏航角ψ，则轨道坐标系与卫星本体坐标系之间的转换矩阵为：

则太阳矢量在基准本体系中的矢量为：

${\overrightarrow{r}}_{sb}=R_{bo}·R_{0J2000}·{\overrightarrow{r}}_s$

本体系中，沿太阳矢量方向的单位向量记为

4)根据卫星在轨运行参数，判断卫星的光照条件，卫星在轨所经历的不同光照弧段如图2所示。

设τ为日蚀因子，日照区τ＝1，半影区τ线性减小，本影区减为τ＝0。每个表面部件受的摄动力根据所处弧段乘对应的日蚀因子。

5)用面积等于卫星外形最大包络的方形像元阵列模拟太阳光源，像元阵列的指向垂直于太阳矢量，距离卫星本体系为一个天文单位，如图3所示，相邻像元阵列的间隔1mm。定义像元阵列坐标系为：坐标原点在像元阵列的中心点，Z轴沿太阳矢量方向，X轴在像元阵列内，平行于像元阵列某一边长，Z轴与X、Y轴符合右手旋转。设每个像元中心点在像元阵列坐标系中的坐标为P_j，k＝(j，k，0)，在每个像元中心添加一个太阳到卫星方向的单位矢量，作为模拟光线，记为R_j,k(j表示矢量在阵列中所处的行，k表示矢量在阵列中所处的列)。j、k的范围为-N，…，0，…，N，N的大小由卫星表面积的最大包络确定。

根据像元阵列坐标系的定义和本体系中太阳矢量的指向，本体坐标系与像元阵列坐标系的转换矩阵为：

L＝R_y(θ_y)·R_z(θ_z)

其中：${\theta }_z=arcsin\frac{y}{\sqrt{x^2+y^2}},{\theta }_y=arcsin\frac{x}{\sqrt{x^2+z^2}}.$

每个像元中心点在本体系中的坐标为：

$P_{j,k}^{\prime }=L^{-1}·P_{j,k}+{\overrightarrow{r}}_{sb}$

6)计算初次照射引起的太阳光压摄动力。由本体系中太阳指向、像元中心点坐标、以及部件数据库，建立光路跟踪函数，判断每个表面部件的实际受照射情况及有效受照面积。通过每个表面部件的安装位置，表面积，确定部件在本体系中的空间函数及边界。记某一像元中心点在本体系中的坐标为(x_si,y_si,z_si)，太阳矢量为(X_s,Y_s,Z_s)^T，则太阳光线的空间方程为：

$\frac{x-x_{si}}{X_s}=\frac{y-y_{si}}{Y_s}=\frac{z-z_{si}}{Z_s}$

计算并记录每条光线与卫星表面距离最短的交点(即为光线所照射到的点)和相应的入射角(光线和被照射面元法线的夹角)，并对该部件有效受照光线数量加1，对每一条照射光线对应有效被照射面元面积记为Δ＝1mm²。通过与卫星部件存在最短距离的光线数，确定初次照射时每个表面部件有效被照射面元数为A_ieff1。

对每一个受太阳光初次照射的部件，以i表示卫星表面部件代号，其所受法向光压摄动力f_ni1为：

$f_{ni1}=-\tau ·\underset{j=1}{\overset{A_{ieff1}}{\Sigma }}(E_{sun}{\mathrm{cos\theta }}_{ij}/c)[(1+{\nu }_i{\mu }_i){\mathrm{cos\theta }}_{ij}+\frac{2}{3}{\nu }_i(1-{\mu }_i)]$

切向光压摄动力f_si1为：

$f_{si1}=\tau \underset{j=1}{\overset{A_{ieff1}}{\Sigma }}(E_{sun}{\mathrm{cos\theta }}_{ij}/c)(1-{\nu }_i{\mu }_i){\mathrm{sin\theta }}_{ij}$

式中：E_sun＝1367W/m²为太阳辐射强度；θ_ij为第i个表面部件的第j个面元的光线入射角。

7)计算二次照射引起的太阳光压摄动力。由于卫星是个构型复杂的结构体，初次入射光线的反射光线会再次照射在某些表面部件上，形成二次太阳光压摄动。

根据卫星表面部件在空间的安装位置、外法向方向P_i＝(x_i,y_i,z_i)^T，确定固定平面部件面元法向与部件相同、固定曲面部件面元法向为星体坐标系原点指向面元中心，转动部件的面元法向由部件运动情况决定。综合考虑太阳光线在星本体坐标系中的矢量，确定反射光线的方向：漫反射光线沿面元法向，镜面反射光线沿入射光线的对称方向。

由漫反射和镜面反射光线的方向和光线出发点(即首次被照射面元点)，利用第6)步中光路跟踪函数法计算反射光线所照射到的表面面元，对每一条照射光线对应有效被照射面元面积记为Δ＝1mm²。通过反射光线与部件最短距离的判断，确定每个表面部件由二次照射的有效被照射面元数为A_ieff2。根据反射类型，被反射光线照射的有效被照射面元数为A_ieff2乘以相应的比例，漫反射的比例为v(1-μ)，镜面反射的比例为vμ。因此被二次照射中漫反射的有效面元数为A_ieff2·v(1-μ)，镜面反射的有效面元数为A_ieff2·vμ。

二次照射中，第i个表面部件由漫反射引起的法向光压摄动力f_ni2d为：

$f_{ni2d}=-\tau \underset{j=1}{\overset{A_{ieff2}·v\left(1-\mu \right)}{\Sigma }}\left(E_{sun}{\mathrm{cos\theta }}_{ij}/c\right)\left[\left(1+{\nu }_i{\mu }_i\right){\mathrm{cos\theta }}_{ij}+\frac{2}{3}{\nu }_i\left(1-{\mu }_i\right)\right]$

二次照射中，第i个表面部件由镜面反射引起的法向光压摄动力f_ni2m和切向光压摄动力f_si2m为：

$f_{ni2m}=-\tau \underset{i=1}{\overset{A_{ieff2}·v\mu }{\Sigma }}\left(E_{sun}{\mathrm{cos\theta }}_{ij}/c\right)\left[\left(1+{\nu }_i{\mu }_i\right){\mathrm{cos\theta }}_{ij}+\frac{2}{3}{\nu }_i\left(1-{\mu }_i\right)\right]$

$f_{si2m}=\tau \underset{j=1}{\overset{A_{ieff2}·v\mu }{\Sigma }}\left(E_{sun}{\mathrm{cos\theta }}_{ij}/c\right)\left(1-{\nu }_i{\mu }_i\right){\mathrm{sin\theta }}_{ij}$

所以，第i个部件所受法向光压摄动力f_ni为：

f_ni＝f_ni1+f_ni2d+f_ni2m

切向光压摄动力f_si为：

f_si＝f_si1+f_si2m

8)根据部件的安装位置、法线方向等，将所有表面部件的太阳光压摄动力(f_ni,f_si)分解为沿本体系三轴方向的三个分力，分别记为(F_xi,F_yi,F_zi)。

并对所有表面部件的地球反照光压摄动力矢量叠加，计算整星在本体系下的太阳光压摄动力和光压摄动加速度：

$F_X=\underset{i=1}{\overset{N}{\Sigma }}{F}_{xi};F_Y=\underset{i=1}{\overset{N}{\Sigma }}{F}_{yi};F_Z=\underset{i=1}{\overset{N}{\Sigma }}{F}_{zi};$

$a_X=\frac{F_X}{m};a_Y=\frac{F_Y}{m};a_Z=\frac{F_Z}{m}.$

9)步骤1)-8)是计算特定太阳矢量方向的太阳光压摄动加速度，在一个卫星运动周期内，以卫星在轨道内运动10°为步长，计算半年内每个卫星运动周期不同时刻的太阳光压摄动加速度。对计算结果，通过傅里叶多项式拟合得到太阳光压摄动加速度的数学模型。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。