一种考虑地球扁率的太阳光照系数确定方法

1.本发明属于天体物理及轨道动力学领域，具体涉及一种考虑地球扁率的太阳光照系数确定方法。


背景技术：

2.计算太阳光照系数是计算太阳对于卫星的光压摄动力的关键环节之一。太阳光压摄动力对于在轨航天器的控制，如姿态，轨道控制来说较为重要，尤其对于高轨航天器以及航天器的高精度控制来说，太阳光压摄动力的影响非常大。所以，计算太阳光照系数具有重要意义。然而，由于获得遮挡面积的解析解十分困难，同时也要计算大量内容。在现有技术中，并没有计算以卫星视角观测太阳被地球遮挡面积百分比的高精度计算方法，进而导致现有技术计算太阳光照系数的精度较低。


技术实现要素：

3.本发明的目的是为解决采用现有方法计算太阳光照系数的精度低的问题，而提出的一种考虑地球扁率的太阳光照系数确定方法。
4.本发明为解决上述技术问题所采取的技术方案是：
5.一种考虑地球扁率的太阳光照系数确定方法，所述方法具体包括以下步骤：
6.步骤s1、根据太阳和地球的几何关系计算地球到太阳方向单位向量e
sun
；
7.步骤s2、根据卫星和地球的几何关系以及e
sun
确定卫星所处的区域，再根据卫星所处的区域确定太阳光照系数；
8.卫星所处的区域为：卫星在考虑扁率的地球的本影区、半影区、全影区或光照区；
9.所述步骤s2中，根据卫星和地球的几何关系以及e
sun
确定卫星所处的区域，其具体过程为：
10.步骤s21、将地球到太阳方向单位向量e
sun
与卫星质心相对于考虑扁率的地球质心的位置向量r＝(r1，r2，r3)进行点乘，即s＝r
t
·esun
，其中，上角标t代表转置，r1，r2，r3为向量r中的元素；
11.若s＞0，则卫星在考虑扁率的地球的光照区，若s≤0，则继续执行步骤s22；
12.步骤s22、计算卫星到太阳方向单位向量e
ss
：
[0013][0014]
其中，r
sun
为太阳质心相对于考虑扁率的地球质心的位置向量；
[0015]
计算太阳、地球和卫星三者所在平面的法向量n：
[0016][0017]
计算卫星所在位置的视太阳半径a：
[0018][0019]
其中，r
sun
为太阳半径；
[0020]
卫星指向太阳质心方向向量和卫星指向太阳切点方向向量之间的旋转矩阵r为：
[0021][0022]
其中，n
×
为斜对称矩阵；
[0023][0024]
其中，n1、n2和n3为斜对称矩阵中的元素；
[0025]
计算卫星到太阳的两个切点s1和s3：
[0026]
s1＝r
t
*e
ss
ꢀꢀ
(7)
[0027]
s3＝r*e
ss
ꢀꢀ
(8)
[0028]
步骤s23、判断卫星与切点s1构成的直线以及卫星与切点s3构成的直线是否同地球相交；
[0029]
若卫星与切点s1构成的直线以及卫星与切点s3构成的直线均与地球相交，则卫星位于考虑扁率的地球的本影区，若卫星与切点s1构成的直线以及卫星与切点s3构成的直线均不与地球相交，则卫星位于考虑扁率的地球的全影区，否则卫星位于考虑扁率的地球的半影区。
[0030]
本发明的有益效果是：
[0031]
采用本发明方法可自动计算出高精度的卫星视角观测太阳被地球遮挡面积百分比，进而计算出高精度的太阳光照系数。将本发明方法用于计算太阳光压摄动力，可以发挥出高精度优势，为航天器高精度姿态控制和轨道控制提供技术支持，本发明方法的精度可以精确到小数点后五位。
[0032]
若后续需要计算不同条件下的遮挡太阳面积百分比时，直接修改输入的卫星相对于地球以及太阳相对于地球的坐标参数，即可获得计算结果，方法实用性强。
附图说明
[0033]
图1为双锥形地影数学模型示意图；
[0034]
图2为切点s1和s3几何意义示意图；
[0035]
图3为点s2的几何意义示意图；
[0036]
图4a)为地影可能情况示意图图一；
[0037]
图4b)为地影可能情况示意图图二；
[0038]
图4c)为地影可能情况示意图图三；
[0039]
图4d)为地影可能情况示意图图四；
[0040]
图4e)为地影可能情况示意图图五；
[0041]
图5为本发明方法的流程图；
[0042]
图6为采用本发明方法的计算结果与商业软件stk对比的结果图；
[0043]
图7为采用本发明方法的计算结果与商业软件stk对比的误差图。
具体实施方式
[0044]
具体实施方式一、结合图1、图2和图5说明本实施方式。本实施方式所述的一种考虑地球扁率的太阳光照系数确定方法，所述方法具体包括以下步骤：
[0045]
步骤s1、根据太阳和地球的几何关系计算地球到太阳方向单位向量e
sun
；
[0046]
步骤s2、根据卫星和地球的几何关系以及e
sun
确定卫星所处的区域，再根据卫星所处的区域确定太阳光照系数；
[0047]
卫星所处的区域为：卫星在考虑扁率的地球的本影区、半影区、全影区或光照区；
[0048]
所述步骤s2中，根据卫星和地球的几何关系以及e
sun
确定卫星所处的区域，其具体过程为：
[0049]
步骤s21、将地球到太阳方向单位向量e
sun
与卫星质心相对于考虑扁率的地球质心的位置向量r＝(r1，r2，r3)进行点乘，即s＝r
t
·esun
，其中，上角标t代表转置，r1，r2，r3为向量r中的元素；
[0050]
若s＞0，则卫星在考虑扁率的地球的光照区，若s≤0，则继续执行步骤s22；
[0051]
步骤s22、计算卫星到太阳方向单位向量e
ss
：
[0052][0053]
其中，r
sun
为太阳质心相对于考虑扁率的地球质心的位置向量；
[0054]
计算太阳、地球和卫星三者所在平面的法向量n：
[0055][0056]
计算卫星所在位置的视太阳半径a：
[0057][0058]
其中，r
sun
为太阳半径，r
sun
的取值为695700km；
[0059]
根据欧拉定理可得，卫星指向太阳质心方向向量和卫星指向太阳切点方向向量之间的旋转矩阵r为：
[0060][0061]
其中，n
×
为斜对称矩阵；
[0062][0063]
其中，n1、n2和n3为斜对称矩阵中的元素；
[0064]
计算卫星到太阳的两个切点s1和s3：
[0065]
s1＝r
t
*e
ss
ꢀꢀ
(7)
[0066]
s3＝r*e
ss
ꢀꢀ
(8)
[0067]
步骤s23、判断卫星与切点s1构成的直线以及卫星与切点s3构成的直线是否同地球
相交；
[0068]
若卫星与切点s1构成的直线以及卫星与切点s3构成的直线均与地球相交，则卫星位于考虑扁率的地球的本影区，若卫星与切点s1构成的直线以及卫星与切点s3构成的直线均不与地球相交，则卫星位于考虑扁率的地球的全影区，否则卫星位于考虑扁率的地球的半影区。
[0069]
具体实施方式二：本实施方式与具体实施方式一不同的是，所述步骤s1的具体过程为：
[0070]
步骤s11、在地球固连坐标系下，将太阳质心相对于考虑扁率的地球质心的位置向量表示为r
sun
＝(d1，d2，d3)；
[0071]
其中，d1，d2，d3为向量r
sun
中的元素；
[0072]
步骤s12、利用式(1)计算地球到太阳方向单位向量e
sun
：
[0073][0074]
其中，||r
sun
||表示r
sun
的2范数。
[0075]
其它步骤及参数与具体实施方式一相同。
[0076]
具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述卫星在考虑扁率的地球的光照区时，太阳光照系数为μ＝1。
[0077]
其它步骤及参数与具体实施方式一或二相同。
[0078]
具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述卫星位于考虑扁率的地球的本影区时，太阳光照系数为μ＝0。
[0079]
其它步骤及参数与具体实施方式一至三之一相同。
[0080]
具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述卫星位于考虑扁率的地球的全影区时，太阳光照系数为μ＝1。
[0081]
其它步骤及参数与具体实施方式一至四之一相同。
[0082]
具体实施方式六：本实施方式与具体实施方式一至五之一不同的是，所述判断卫星与切点s1构成的直线以及卫星与切点s3构成的直线是否同地球相交，其具体过程为：
[0083]
步骤s231、将切点s1＝(a1，a2，a3)与卫星构成的直线的方向向量表示为b＝(b1，b2，b3)；
[0084]
步骤s232、切点s1与卫星构成的直线的空间直线方程为：
[0085][0086]
在地球固连坐标系下，地球的方程为：
[0087][0088]
其中，p为地球极半径，p的取值为6356.752km，e为赤道半径，e的取值为6378.137km；
[0089]
联立式(10)和式(11)得：
[0090][0091]
计算是否有交点判据δ：
[0092]
δ＝b
2-4ac
ꢀꢀ
(12)
[0093]
若δ大于0，则切点s1与卫星构成的直线同地球相交，否则，切点s1与卫星构成的直线同地球不相交；
[0094]
步骤s233、同理，判断出切点s3与卫星构成的直线是否同地球相交。
[0095]
在直线同地球相交的情况下，计算两个交点坐标(x1，y1，z1)和(x2，y2，z2)，点(x1，y1，z1)距地心(0，0，0)的距离r1以及点(x2，y2，z2)距地心(0，0，0)的距离r2；
[0096]
两点距地心(0，0，0)的距离r1，r2为：
[0097][0098][0099][0100][0101]
并计算出两个交点的距离d＝|r
1-r2|，以及两点到原点的平均距离
[0102]
其它步骤及参数与具体实施方式一至五之一相同。
[0103]
具体实施方式七：结合图3说明本实施方式。本实施方式与具体实施方式一至六之一不同的是，所述卫星位于考虑扁率的地球的半影区时，太阳光照系数的确定方法为：
[0104]
步骤1、初始化迭代角迭代判据γ以及迭代中止阈值m；
[0105]
步骤2、计算迭代角的旋转矩阵a，再利用旋转矩阵a计算太阳上的点s2：
[0106]
旋转矩阵a的计算方法与旋转矩阵r相同；
[0107]
s2＝a*s1ꢀꢀ
(17)
[0108]
计算卫星与太阳上的点s2的连线是否与地球相交，若相交，则计算两个交点的距离d2以及两个交点距地心的平均距离e2，并将s1赋值为s2，将迭代判据赋值为γ＝d2；若不相交，则不做处理；
[0109]
步骤3、将迭代角进行二分，即将迭代角更新为
[0110]
步骤4，利用更新后的迭代角重复步骤2至步骤3的过程，直至迭代判据γ小于m时停止，将最后一次迭代获得的s2作为地球在投影面的切点，利用最后一次迭代获得的e2来执行步骤5；
[0111]
步骤5、根据e2计算地球切点的视半径卫星到地球方向矢量与卫星到太阳方向矢量的夹角
[0112]
步骤6、计算a、b、c构成的三角形的高，a在c上的投影，进而计算出阴影面积a：
[0113][0114][0115][0116]
其中，x
′
和y
′
为中间变量；
[0117]
地影可能情况的示意图如图4a)至图4e)所示；
[0118]
步骤7、根据阴影面积a计算出卫星位于考虑扁率的地球的半影区时的太阳光照系数μ：
[0119][0120]
其它步骤及参数与具体实施方式一至六之一相同。
[0121]
如图6所示为本发明方法的计算结果与商业软件stk对比的结果图，图7所示为本发明方法的计算结果与商业软件stk对比的误差图，通过图6和图7验证了本发明方法的有效性。
[0122]
本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。