一种适用于椭圆轨道的半解析阴影区预报方法技术领域
本发明涉及一种适用于椭圆轨道的半解析阴影区预报方法，主要在以太阳能电推
进为动力的地球卫星和深空探测器上使用，用于在卫星上预测任意倾角和偏心率椭圆轨道
在一个轨道周期内的阴影区范围。
背景技术
对于太阳能电推进卫星，电推力器是实现卫星轨道转移、位置保持控制和角动量
卸载的执行机构，需要较长时间段连续、频繁的开机工作。然而太阳能电推进系统高度依赖
于太阳电池阵为其提供充足的电能，若卫星进入阴影区，则电推力器不能工作。因此，为保
证电推进卫星的星上自主控制能够正常施行，有必要在每一轨道周期内对阴影区的范围进
行快速预测。
现有的针对于任意倾角和偏心率的阴影区预报方法，主要依靠全解析计算、数值
迭代、数值积分等方法，当应用于星上时，具有以下不足之处：①全解析的阴影区预报方法
需要对四次方程式求解，涉及到对虚数方程的求解，且步骤复杂；②全数值方法计算量过
大，且在星上计算能力有限的情况下精度较低；③数值积分法对于椭圆轨道而言，需使用变
步长积分，且在近地点附近计算量很大。因此，需要对阴影区预报方法进行重新设计，满足
电推进卫星的星上自主计算需求。
发明内容
本发明所要解决的技术问题是：克服现有技术的不足，提供了一种适用于椭圆轨
道的半解析阴影区预报方法，本方法采用半解析半数值的方法，首先利用解析算式估算阴
影区的可能存在范围，之后对是否存在阴影区进行粗略判断，并对阴影区类型进行分类，最
后，利用数值方法根据在阴影区可能出现的位置内迭代求解阴影区存在范围。本发明中的
方法降低了阴影区预报方法的计算量，对星上计算能力要求不高，满足了电推进卫星进行
星上自主控制的任务需求。
本发明包括如下技术方案：
一种适用于椭圆轨道的半解析阴影区预报方法，该方法包括如下步骤：
(1)将太阳矢量从惯性坐标系转换到地心轨道坐标系；
(2)依据阴影区桶形模型，获得阴影区在椭圆轨道平面内的方向角、阴影区投影椭
圆半长轴；
(3)解析计算阴影区可能存在范围；
(4)根据阴影椭圆和航天器所处轨道间的几何关系对阴影区可能出现的位置进行
分类，并对阴影区的存在进行粗略判断，若阴影区存在范围属于a类，则进入步骤(5)，若阴
影区存在范围属于b类，则进入步骤(6)，若不存在阴影区，结束本次预报；
(5)获取a类阴影区范围，在阴影区可能存在范围区域内，直接利用数值方法求得
阴影区存在范围，进入步骤(7)；
(6)获取对b类阴影区范围，缩小阴影区存在范围，搜索在此范围内进入阴影区的
任意点，若在有限的迭代次数下不能搜到该点，则判定卫星在本轨道周期内不进入阴影，结
束本次预报；若找到阴影点，根据该点位置利用数值方法求得阴影区存在范围，进入步骤
(7)；
(7)将阴影区存在范围转换为由偏近点角E、平近点角M以及时间表示的参数。
所述(1)中获得太阳矢量在地心轨道坐标系中的投影具体为：
已知航天器轨道半长轴a，偏心率e，倾角i，升交点赤经Ω，近地点幅角ω，则太阳
矢量在地心轨道坐标系(如图2所示，原点o在中心天体质心，X轴指向近地点，Y轴沿轨道法
向，Z轴与X、Y轴呈右手系)中的投影为：
其中为太阳在中心天体惯性系下的三维单位矢量，
所述步骤(2)中获得阴影区在轨道平面内的方向角、阴影区投影椭圆半长轴，具体
为：
(2-1)阴影区在椭圆轨道平面内的方向角θu，具体由公式：
给出，其中，为太阳在中心天体轨道坐标系中的三维单位矢量，为的
第一个分量，为的第二个分量。
(2-2)阴影区投影椭圆半长轴au，具体由公式：
给出，其中Re为中心天体半径，为的第三个分量。
所述步骤(3)中获取阴影区可能存在范围，具体为：
(3-1)所述阴影区可能存在范围两侧的真近点角f1、f2为：
f1＝f(θu，a，e，Re)；
f2＝f(θu，a，e，-Re)；
则真近点角f1、f2与方向角θu的距离分别为
df1＝f1-θu；
df2＝f2-θu；
将df1、df2到[-π，π]；
其中，a为轨道半长轴，e为轨道偏心率，函数F＝f(θu，a，e，R)的具体方法为：
如果：θu＝0或θu＝π
则：L1＝-1/e；
L2＝a·(1-e2)/e；
L4＝a·(1-e2)；
L5＝R2e2+a2-2a2e2+a2e4；
Lx1＝L1+L2·(L3+L4)/L5；
Lx2＝L1-L2·(L3-L4)/L5；
Ly1＝(L3+L4)·R/L5；
Ly2＝-(L3-L4)·R/L5；
f＝[atan2(Ly1，Lx1)，atan2(Ly2，Lx2)]；
其中，L1、L2、L3、L4、L5、Lx1、Lx2、Ly1、Ly2为中间变量。如果：或
则：L1＝acos(R/(a(1-e2)-R·e))；
f＝[L1，2π-L1]；
如果：θu≠0、θu≠π、
k＝tanθu；
R2＝R/sinθu；
L1＝-R2·k/(a(1-e2))；
L2＝k(a(1-e2)-e·R2)/(a(1-e2))；
L4＝R2·k2(a(1-e2)-e·R2)；
L5＝k2(a(1-e2)-e·R2)2+a2(1-e2)2；
Lx1＝(L3+L4)/L5；
Lx2＝-(L3-L4)/L5；
Ly1＝L1+L2·(L3+L4)/L5；
Ly2＝L1-L2·(L3-L4)/L5；
f＝[atan2(Ly1，Lx1)，atan2(Ly2，Lx2)]；
其中，L1、L2、L3、L4、L5、Lx1、Lx2、Ly1、Ly2为中间变量。df＝|f-θu|；
将df化到[-π，π]；
如果：dF(1)≤df(2)
则：F＝f(1)；
如果：dF(1)＞df(2)
则：F＝f(2)；
(3-2)对阴影区可能存在范围两侧的真近点角进行排序：
如果：df1≤df2
阴影区可能存在范围ff＝[f1，f2]；
df＝[df1，df2]；
如果：df1＞df2
阴影区可能存在范围ff＝[f2，f1]；
df＝[df2，df1]；
所述步骤(4)中对阴影区可能出现的位置进行分类，具体为：
阴影区的分类依据，具体由步骤：
(4-1)如果：r(θu)≤au
则：阴影区类型为a类，如图3(a)所示；
(4-2)如果：在阴影区可能存在范围ff中存在f*，使r(f*)≤au
则：可能存在b类阴影区，如图3(b)所示；
(4-3)如果：在阴影区可能存在范围ff中不存在f*，使r(f*)≤au
则：不存在阴影区，退出本方法。
其中r(f)＝a(1-e2)/(1+ecosf)为任意真近点角所对应的轨道半径。
r(f*)＝a(1-e2)/(1+ecosf*)为真近点角为f*时所对应的轨道半径，
r(θu)＝a(1-e2)/(1+ecosθu)为真近点角为θu时所对应的轨道半径。
所述步骤(5)中获得a类阴影区范围，具体为：
(5-1)获取数值迭代次数iter，具体由公式：
给出，其中∈为容许阴影区的最大误差。
(5-2)对任意真近点角f，是否进入阴影区的判别方式，具体由公式：
给出。当则进入阴影区，当则不进入阴影区。
(5-3)用数值方法获得a类阴影区范围：根据(5-2)中的阴影区判别公式
用二分法在区间[ff(1),θu]和[θu,ff(2)]内寻找阴影区范围fu＝[fu(1),fu
(2)]，迭代次数分别为iter(1)和iter(2)。
所述步骤(6)中获得b类阴影区范围，具体为：
(6-1)在ff＝[f1，f2]中搜索一个进入阴影区的点，具体方法为：
如果：θu≥π
则：
其中i＝1,2,3,…,n，j＝1到2i-1，n为正整数。
阴影区的判别公式为
如果：
则：计算结束
fmid＝ftempi,j；
ff＝[θu，f2]；
df＝ff-fmid；
将df化到[-π，π]；
其中fmid为中间变量，随i，j更新。
如果i＝1,2,3,…,n，始终则本轨道周期内不存在阴影
区，计算结束。
如果：θu＜π
则：
其中i＝1,2,3,…,n，j＝1到2i-1，n为正整数。
阴影区的判别公式为
如果：
则：计算结束
fmid＝ftempi,j；
ff＝[f1，θu]；
df＝ff-fmid；
将df化到[-π，π]；
其中fmid为中间变量，随i，j更新。
如果i＝1,2,3,…,n，始终则本轨道周期内不存在阴影
区，计算结束。
(6-2)在ff中搜索b类阴影区范围，所述步骤为：
获取数值迭代次数
根据阴影区判别公式用二分法在区间[ff(1),fmid]和
[fmid，ff(2)]内寻找阴影区范围fu＝[fu(1),fu(2)]，迭代次数分别为iter(1)和
iter(2)。其中，ff(1)为ff的第一个分量，ff(2)为ff的第二个分量。
所述步骤(7)将阴影区存在范围转换为由偏近点角Eu、平近点角Mu或时间tu表示的
参数，具体为：
阴影区所在的偏近点角范围为
将Eu化到[0，2π]；
平近点角Mu＝Eu-e·sinEu；
ΔM＝Mu-M0；
将ΔM化到[0，2π]；
时间
tu＝Δt+t0；
其中M0为参考时刻平近点角，t0为参考时刻。ΔM、Δt为中间变量。
本发明与现有技术相比的有益效果是：
(1)本发明方法首先采用解析的方法事先求得阴影区的大致存在范围，避免利用
数值方法对卫星整轨遇阴影情况进行搜索，可显著减少计算量，在同等计算能力下提高阴
影范围的求解经度；
(2)本发明方法首先将阴影区的类型分为两类，数值计算阴影区的具体存在范围，
避免用全解析法所面临的虚数多项式求解问题，使方法更加简化；
(3)本发明采取的是几何分析的方法，避免进行数值积分，降低了阴影区预报的计
算量，提高了大偏心率椭圆轨道阴影区预测的计算精度；
(4)本发明采取的方法每个轨道周期仅需计算一次，降低了计算量。
附图说明
图1为本发明阴影区预报方法流程框图；
图2为地心轨道系示意图；
图3为阴影区类型；
图4仿真算例。
具体实施方式
下面结合附图对本发明的具体实施方式进行进一步的详细描述。
当中心天体与太阳距离较大时，可将太阳光看作是平行光源，则中心天体遮挡太
阳所形成的阴影为桶形模型。因此，中心天体造成的阴影在轨道平面上的投影为椭圆形，椭
圆形中心位于中心天体的几何中心。若将中心天体质心和几何中心看作重合，则求解航天
器轨道的阴影区范围即为求解一个中点在原点的椭圆和一个焦点在原点的椭圆的交点问
题。
图1为阴影区预报方法的流程框图。如图1所示，本发明所采用的适用于椭圆轨道
的半解析阴影区预报方法首先进行地影区粗判断，再利用解析方法界定阴影区的可能存在
范围，最后利用数值方法搜索阴影区的确切位置。具体步骤如下：
(1)计算太阳矢量在地心轨道坐标系中的投影。
已知航天器轨道半长轴a，偏心率e，倾角i，升交点赤经Ω，近地点幅角ω，则太阳
矢量在地心轨道坐标系(如图2所示，原点在中心天体质心，X轴指向近地点，Y轴沿轨道法
向，Z轴与X、Y轴呈右手系)中的投影为：
其中为太阳在中心天体惯性系下的三维单位矢量，
(2)计算阴影区在轨道平面内的方向角、阴影区投影椭圆半长轴。
其中Re为中心天体半径。
(3)计算阴影区可能存在范围。阴影区可能存在范围两侧的近地点为：
f1＝f(θu，a，e，Re)；
f2＝f(θu，a，e，-Re)；
df1＝f1-θu；
df2＝f2-θu；
将df1、df2到[-π，π]；
如果：df1≤df2
阴影区可能存在范围ff＝[f1，f2]；
df＝[df1，df2]；
其他：阴影区可能存在范围ff＝[f2，f1]；
df＝[df2，df1]；
结束
其中，函数F＝f(θu，a，e，R)的具体方法为
如果：θu＝0或θu＝π
则：L1＝-1/e；
L2＝a·(1-e2)/e；
L4＝a·(1-e2)；
L5＝R2e2+a2-2a2e2+a2e4；
Lx1＝L1+L2·(L3+L4)/L5；
Lx2＝L1-L2·(L3-L4)/L5；
Ly1＝(L3+L4)·R/L5；
Ly2＝-(L3-L4)·R/L5；
f＝[atan2(Ly1，Lx1)，atan2(Ly2，Lx2)]；
如果：或
则：L1＝acos(R/(a(1-e2)-R·e))；
f＝[L1，2π-L1]；
其他：
k＝tanθu；
R2＝R/sinθu；
L1＝-R2·k/(a(1-e2))；
L2＝k(a(1-e2)-e·R2)/(a(1-e2))；
L4＝R2·k2(a(1-e2)-e·R2)；
L5＝k2(a(1-e2)-e·R2)2+a2(1-e2)2；
Lx1＝(L3+L4)/L5；
Lx2＝-(L3-L4)/L5；
Ly1＝L1+L2·(L3+L4)/L5；
Ly2＝L1-L2·(L3-L4)/L5；
f＝[atan2(Ly1，Lx1)，atan2(Ly2，Lx2)]；
结束
df＝|f-θu|；
将df化到[-π，π]；
如果：df(1)≤df(2)
则：F＝f(1)；
其他：F＝f(2)；
结束
(4)根据阴影椭圆的大小和方向等参数对阴影区可能出现的位置进行分类。
航天器在任意近地点的轨道半径为r(f)＝a(1-e2)/(1+ecosf)
如果：r(θu)≤au
则：阴影区类型为a类，如图3(a)所示，转入步骤(5)；
如果：在阴影区可能存在范围ff中存在f*，使r(f*)≤au
则：可能存在b类阴影区，如图3(b)所示，转入步骤(6)；
其他：不存在阴影区，退出本方法。
结束
(5)计算a类阴影区范围。
计算数值迭代次数其中∈为容许阴影区计算的最
大误差。
对任意真近点角f，是否进入阴影区的判别公式为：
根据以上公式，用二分法在区间[ff(1),θu]和[θu,ff(2)]内寻找阴影区范围fu＝
[fu(1),fu(2)]，迭代次数分别为iter(1)和iter(2)(关于二分法的详细方法，详见《数学手
册》，北京：高等教育出版社，1979：103-104)。
(6)计算b类阴影区范围。
①在ff＝[f1，f2]中搜索一个进入阴影区的点：
设n为迭代次数。
如果：θu≥π
则：循环：i从1到n
循环：j从1到2i-1
如果：阴影区的判别公式
则：fmid＝ftemp；
ff＝[θu，f2]；
df＝ff-fmid；
将df化到[-π，π]；
结束i循环，结束j循环；
结束
结束
结束
如果：θu＜π
则：循环：i从1到n
循环：j从1到2i-1
如果：阴影区的判别公式
则：fmid＝ftemp；
ff＝[f1，θu]；
结束i循环，结束j循环；
结束
结束
结束
结束
df＝ff-fmid；
将df化到[-π，π]；
②在ff中搜索b类阴影区范围：
计算数值迭代次数
根据阴影区判别公式用二分法在区间[ff(1),fmid]和
[fmid,ff(2)]内寻找阴影区范围fu＝[fu(1),fu(2)]，迭代次数分别为iter(1)和
iter(2)(关于二分法的详细方法，详见《数学手册》，北京：高等教育出版社，1979：103-
104)。
(7)将阴影区存在范围转换为由偏近点角E、平近点角M或时间表示的参数。
将Eu化到[0，2π]；
平近点角Mu＝Eu-e·sinEu；
ΔM＝Mu-M0；
将ΔM化到[0，2π]；
时间
tu＝Δt+t0；
其中M0为参考时刻平近点角，t0为参考时刻。ΔM、Δt为中间变量。
(8)仿真结果
设航天器的轨道根数、太阳方向角如下表所示：
表1仿真条件
仿真结果如图4所示，算例I所求得的阴影类型为a类，阴影区范围为237.9°～
265.2°，算例I所求得的阴影类型为b类，阴影区范围为300.6°～313.2°。
经仿真验证，本发明方法能够实现对任意倾角、任意偏心率的椭圆轨道的阴影区
预测，计算量满足星载计算机的计算能力约束。
实施例
一种适用于椭圆轨道的半解析阴影区预报方法，该方法包括如下步骤：
(1)将太阳矢量从惯性坐标系转换到地心轨道坐标系；
(2)依据阴影区桶形模型，获得阴影区在椭圆轨道平面内的方向角、阴影区投影椭
圆半长轴；
(3)解析计算阴影区可能存在范围；
(4)根据阴影椭圆和航天器所处轨道间的几何关系对阴影区可能出现的位置进行
分类，并对阴影区的存在进行粗略判断，若阴影区存在范围属于a类，则进入步骤(5)，若阴
影区存在范围属于b类，则进入步骤(6)，若不存在阴影区，结束本次预报；
(5)获取a类阴影区范围，在阴影区可能存在范围区域内，直接利用数值方法求得
阴影区存在范围，进入步骤(7)；
(6)获取对b类阴影区范围，缩小阴影区存在范围，搜索在此范围内进入阴影区的
任意点，若在有限的迭代次数下不能搜到该点，则判定卫星在本轨道周期内不进入阴影，结
束本次预报；若找到阴影点，根据该点位置利用数值方法求得阴影区存在范围，进入步骤
(7)；
(7)将阴影区存在范围转换为由偏近点角E、平近点角M以及时间表示的参数。
所述(1)中计算太阳矢量在地心轨道坐标系中的投影具体为：
已知航天器近地点高度hp为200km，远地点高度ha为24371km，倾角i＝0°，升交点赤
经Ω＝0°，近地点幅角ω＝0°，则太阳矢量在地心轨道坐标系中的投影为：
所述步骤(2)中计算阴影区在轨道平面内的方向角、阴影区投影椭圆半长轴，具体
为：
(2-1)阴影区在椭圆轨道平面内的方向角θu，具体由公式：
给出。
(2-2)阴影区投影椭圆半长轴au，具体由公式：
给出。
所述步骤(3)中计算阴影区可能存在范围，具体为：
(3-1)所述阴影区可能存在范围两侧的真近点角f1、f2为：
f1＝f(θu，a，e，Re)＝222.2°；
f2＝f(θu，a，e，-Re)＝281.8°；
则真近点角f1、f2与方向角θu的距离分别为
df1＝f1-θu＝-17.8°；
df2＝f2-θu＝41.8°；
将df1、df2到[-π，π]；
(3-2)对阴影区可能存在范围两侧的真近点角进行排序：
df1≤df2；
阴影区可能存在范围ff＝[f1，f2]；
df＝[df1，df2]；
所述步骤(4)中对阴影区可能出现的位置进行分类，具体为：
阴影区的分类依据，具体由步骤：
r(θu)≤au；
阴影区类型为a类，如图3(a)所示；
所述步骤(5)中计算a类阴影区范围，具体为：
(5-1)计算数值迭代次数iter，设∈＝0.1°：
(5-2)对任意真近点角f，是否进入阴影区的判别方式，具体由公式：
给出。当则进入阴影区，当则不进入阴影区。
(5-3)用数值方法计算a类阴影区范围：根据(5-2)中的阴影区判别公式
用二分法在区间[ff(1),θu]和[u,f(2)]内寻找阴影区范围fu＝[u(1),fu(2)]
＝[279.0°265.2°]。
所述步骤(6)中计算b类阴影区范围，具体为：
本实施例不进入此分支。
所述步骤(7)将阴影区存在范围转换为由偏近点角Eu、平近点角Mu或时间tu表示的
参数，具体为：
阴影区所在的偏近点角范围为
将Eu化到[0，2]；
平近点角Mu＝Eu-e·sinEu＝[317.3°307.1°]；
利用本发明一种适用于椭圆轨道的半解析阴影区预报方法，可预估一个轨道周期
内阴影区的具体范围，使太阳能电推进卫星在获知阴影区位置后提前关机，避免阴影区电
推力器开机引起星上其他载荷发生能源不足的现象。
本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。