CN106570234B - 一种适用于椭圆轨道的半解析阴影区预报方法 - Google Patents

Abstract

一种适用于椭圆轨道的半解析阴影区预报方法，采用半解析半数值的方法，预测任意倾角椭圆轨道在一个轨道周期内的阴影区范围。首先，根据太阳与轨道的关系，利用解析方法得到阴影区的可能存在范围；之后，对是否存在阴影区进行粗略判断，并对阴影类型进行分类；最后，根据不同的阴影类型在有限的迭代次数内利用数值法求得阴影区存在范围。本发明在估算任意轨道阴影区范围的过程中尽可能使用解析方法代替数值积分和迭代，降低了阴影区预报的计算量，为实现以太阳能电推进为动力的地球卫星和深空探测器的星上自主阴影区预报提供了一种有效方法。

Description

一种适用于椭圆轨道的半解析阴影区预报方法

技术领域

本发明涉及一种适用于椭圆轨道的半解析阴影区预报方法，主要在以太阳能电推进为动力的地球卫星和深空探测器上使用，用于在卫星上预测任意倾角和偏心率椭圆轨道在一个轨道周期内的阴影区范围。

背景技术

对于太阳能电推进卫星，电推力器是实现卫星轨道转移、位置保持控制和角动量卸载的执行机构，需要较长时间段连续、频繁的开机工作。然而太阳能电推进***高度依赖于太阳电池阵为其提供充足的电能，若卫星进入阴影区，则电推力器不能工作。因此，为保证电推进卫星的星上自主控制能够正常施行，有必要在每一轨道周期内对阴影区的范围进行快速预测。

现有的针对于任意倾角和偏心率的阴影区预报方法，主要依靠全解析计算、数值迭代、数值积分等方法，当应用于星上时，具有以下不足之处：①全解析的阴影区预报方法需要对四次方程式求解，涉及到对虚数方程的求解，且步骤复杂；②全数值方法计算量过大，且在星上计算能力有限的情况下精度较低；③数值积分法对于椭圆轨道而言，需使用变步长积分，且在近地点附近计算量很大。因此，需要对阴影区预报方法进行重新设计，满足电推进卫星的星上自主计算需求。

发明内容

本发明所要解决的技术问题是：克服现有技术的不足，提供了一种适用于椭圆轨道的半解析阴影区预报方法，本方法采用半解析半数值的方法，首先利用解析算式估算阴影区的可能存在范围，之后对是否存在阴影区进行粗略判断，并对阴影区类型进行分类，最后，利用数值方法根据在阴影区可能出现的位置内迭代求解阴影区存在范围。本发明中的方法降低了阴影区预报方法的计算量，对星上计算能力要求不高，满足了电推进卫星进行星上自主控制的任务需求。

本发明包括如下技术方案：

一种适用于椭圆轨道的半解析阴影区预报方法，该方法包括如下步骤：

(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)阴影区投影椭圆半长轴a_u，具体由公式：

给出，其中R_e为中心天体半径，为的第三个分量。

所述步骤(3)中获取阴影区可能存在范围，具体为：

(3-1)所述阴影区可能存在范围两侧的真近点角f₁、f₂为：

f₁＝f(θ_u，a，e，R_e)；

f₂＝f(θ_u，a，e，-R_e)；

则真近点角f₁、f₂与方向角θ_u的距离分别为

df₁＝f₁-θ_u；

df₂＝f₂-θ_u；

将df₁、df₂到[-π，π]；

其中，a为轨道半长轴，e为轨道偏心率，函数F＝f(θ_u，a，e，R)的具体方法为：

如果：θ_u＝0或θ_u＝π

则：L1＝-1/e；

L2＝a·(1-e²)/e；

L4＝a·(1-e²)；

L5＝R²e²+a²-2a²e²+a²e⁴；

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-e²)-R·e))；

f＝[L1，2π-L1]；

如果：θ_u≠0、θ_u≠π、

k＝tanθ_u；

R2＝R/sinθ_u；

L1＝-R2·k/(a(1-e²))；

L2＝k(a(1-e²)-e·R2)/(a(1-e²))；

L4＝R2·k²(a(1-e²)-e·R2)；

L5＝k²(a(1-e²)-e·R2)²+a²(1-e²)²；

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)对阴影区可能存在范围两侧的真近点角进行排序：

如果：df₁≤df₂

阴影区可能存在范围ff＝[f₁，f₂]；

df＝[df₁，df₂]；

如果：df₁＞df₂

阴影区可能存在范围ff＝[f₂，f₁]；

df＝[df₂，df₁]；

所述步骤(4)中对阴影区可能出现的位置进行分类，具体为：

阴影区的分类依据，具体由步骤：

(4-1)如果：r(θ_u)≤a_u

则：阴影区类型为a类，如图3(a)所示；

(4-2)如果：在阴影区可能存在范围ff中存在f^*，使r(f^*)≤a_u

则：可能存在b类阴影区，如图3(b)所示；

(4-3)如果：在阴影区可能存在范围ff中不存在f^*，使r(f^*)≤a_u

则：不存在阴影区，退出本方法。

其中r(f)＝a(1-e²)/(1+ecosf)为任意真近点角所对应的轨道半径。

r(f^*)＝a(1-e²)/(1+ecosf^*)为真近点角为f^*时所对应的轨道半径，

r(θ_u)＝a(1-e²)/(1+ecosθ_u)为真近点角为θ_u时所对应的轨道半径。

所述步骤(5)中获得a类阴影区范围，具体为：

(5-1)获取数值迭代次数iter，具体由公式：

给出，其中∈为容许阴影区的最大误差。

(5-2)对任意真近点角f，是否进入阴影区的判别方式，具体由公式：

给出。当则进入阴影区，当则不进入阴影区。

(5-3)用数值方法获得a类阴影区范围：根据(5-2)中的阴影区判别公式用二分法在区间[ff(1),θ_u]和[θ_u,ff(2)]内寻找阴影区范围f_u＝[f_u(1),f_u(2)]，迭代次数分别为iter(1)和iter(2)。

所述步骤(6)中获得b类阴影区范围，具体为：

(6-1)在ff＝[f₁，f₂]中搜索一个进入阴影区的点，具体方法为：

如果：θ_u≥π

则：

其中i＝1,2,3,…,n，j＝1到2^i-1，n为正整数。

阴影区的判别公式为

如果：

则：计算结束

f_mid＝f_tempi,j；

ff＝[θ_u，f₂]；

df＝ff-f_mid；

将df化到[-π，π]；

其中f_mid为中间变量，随i，j更新。

如果i＝1,2,3,…,n，始终则本轨道周期内不存在阴影区，计算结束。

如果：θ_u＜π

则：

其中i＝1,2,3,…,n，j＝1到2^i-1，n为正整数。

阴影区的判别公式为

如果：

则：计算结束

f_mid＝f_tempi,j；

ff＝[f₁，θ_u]；

df＝ff-f_mid；

将df化到[-π，π]；

其中f_mid为中间变量，随i，j更新。

如果i＝1,2,3,…,n，始终则本轨道周期内不存在阴影区，计算结束。

(6-2)在ff中搜索b类阴影区范围，所述步骤为：

获取数值迭代次数

根据阴影区判别公式用二分法在区间[ff(1),f_mid]和

[f_mid，ff(2)]内寻找阴影区范围f_u＝[f_u(1),f_u(2)]，迭代次数分别为iter(1)和iter(2)。其中，ff(1)为ff的第一个分量，ff(2)为ff的第二个分量。

所述步骤(7)将阴影区存在范围转换为由偏近点角E_u、平近点角M_u或时间t_u表示的参数，具体为：

阴影区所在的偏近点角范围为

将E_u化到[0，2π]；

平近点角M_u＝E_u-e·sinE_u；

ΔM＝M_u-M₀；

将ΔM化到[0，2π]；

时间

t_u＝Δt+t₀；

其中M₀为参考时刻平近点角，t₀为参考时刻。ΔM、Δt为中间变量。

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

(1)本发明方法首先采用解析的方法事先求得阴影区的大致存在范围，避免利用数值方法对卫星整轨遇阴影情况进行搜索，可显著减少计算量，在同等计算能力下提高阴影范围的求解经度；

(2)本发明方法首先将阴影区的类型分为两类，数值计算阴影区的具体存在范围，避免用全解析法所面临的虚数多项式求解问题，使方法更加简化；

(3)本发明采取的是几何分析的方法，避免进行数值积分，降低了阴影区预报的计算量，提高了大偏心率椭圆轨道阴影区预测的计算精度；

(4)本发明采取的方法每个轨道周期仅需计算一次，降低了计算量。

附图说明

图1为本发明阴影区预报方法流程框图；

图2为地心轨道系示意图；

图3为阴影区类型；

图4仿真算例。

具体实施方式

下面结合附图对本发明的具体实施方式进行进一步的详细描述。

当中心天体与太阳距离较大时，可将太阳光看作是平行光源，则中心天体遮挡太阳所形成的阴影为桶形模型。因此，中心天体造成的阴影在轨道平面上的投影为椭圆形，椭圆形中心位于中心天体的几何中心。若将中心天体质心和几何中心看作重合，则求解航天器轨道的阴影区范围即为求解一个中点在原点的椭圆和一个焦点在原点的椭圆的交点问题。

图1为阴影区预报方法的流程框图。如图1所示，本发明所采用的适用于椭圆轨道的半解析阴影区预报方法首先进行地影区粗判断，再利用解析方法界定阴影区的可能存在范围，最后利用数值方法搜索阴影区的确切位置。具体步骤如下：

(1)计算太阳矢量在地心轨道坐标系中的投影。

已知航天器轨道半长轴a，偏心率e，倾角i，升交点赤经Ω，近地点幅角ω，则太阳矢量在地心轨道坐标系(如图2所示，原点在中心天体质心，X轴指向近地点，Y轴沿轨道法向，Z轴与X、Y轴呈右手系)中的投影为：

其中为太阳在中心天体惯性系下的三维单位矢量，

(2)计算阴影区在轨道平面内的方向角、阴影区投影椭圆半长轴。

其中R_e为中心天体半径。

(3)计算阴影区可能存在范围。阴影区可能存在范围两侧的近地点为：

f₁＝f(θ_u，a，e，R_e)；

f₂＝f(θ_u，a，e，-R_e)；

df₁＝f₁-θ_u；

df₂＝f₂-θ_u；

将df₁、df₂到[-π，π]；

如果：df₁≤df₂

阴影区可能存在范围ff＝[f₁，f₂]；

df＝[df₁，df₂]；

其他：阴影区可能存在范围ff＝[f₂，f₁]；

df＝[df₂，df₁]；

结束

其中，函数F＝f(θ_u，a，e，R)的具体方法为

如果：θ_u＝0或θ_u＝π

则：L1＝-1/e；

L2＝a·(1-e²)/e；

L4＝a·(1-e²)；

L5＝R²e²+a²-2a²e²+a²e⁴；

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-e²)-R·e))；

f＝[L1，2π-L1]；

其他：

k＝tanθ_u；

R2＝R/sinθ_u；

L1＝-R2·k/(a(1-e²))；

L2＝k(a(1-e²)-e·R2)/(a(1-e²))；

L4＝R2·k²(a(1-e²)-e·R2)；

L5＝k²(a(1-e²)-e·R2)²+a²(1-e²)²；

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-e²)/(1+ecosf)

如果：r(θ_u)≤a_u

则：阴影区类型为a类，如图3(a)所示，转入步骤(5)；

如果：在阴影区可能存在范围ff中存在f^*，使r(f^*)≤a_u

则：可能存在b类阴影区，如图3(b)所示，转入步骤(6)；

其他：不存在阴影区，退出本方法。

结束

(5)计算a类阴影区范围。

计算数值迭代次数其中∈为容许阴影区计算的最大误差。

对任意真近点角f，是否进入阴影区的判别公式为：

根据以上公式，用二分法在区间[ff(1),θ_u]和[θ_u,ff(2)]内寻找阴影区范围f_u＝[f_u(1),f_u(2)]，迭代次数分别为iter(1)和iter(2)(关于二分法的详细方法，详见《数学手册》，北京：高等教育出版社，1979：103-104)。

(6)计算b类阴影区范围。

①在ff＝[f₁，f₂]中搜索一个进入阴影区的点：

设n为迭代次数。

如果：θ_u≥π

则：循环：i从1到n

循环：j从1到2^i-1

如果：阴影区的判别公式

则：f_mid＝f_temp；

ff＝[θ_u，f₂]；

df＝ff-f_mid；

将df化到[-π，π]；

结束i循环，结束j循环；

结束

结束

结束

如果：θ_u＜π

则：循环：i从1到n

循环：j从1到2^i-1

如果：阴影区的判别公式

则：f_mid＝f_temp；

ff＝[f₁，θ_u]；

结束i循环，结束j循环；

结束

结束

结束

结束

df＝ff-f_mid；

将df化到[-π，π]；

②在ff中搜索b类阴影区范围：

计算数值迭代次数

根据阴影区判别公式用二分法在区间[ff(1),f_mid]和

[f_mid,ff(2)]内寻找阴影区范围f_u＝[f_u(1),f_u(2)]，迭代次数分别为iter(1)和iter(2)(关于二分法的详细方法，详见《数学手册》，北京：高等教育出版社，1979：103-104)。

(7)将阴影区存在范围转换为由偏近点角E、平近点角M或时间表示的参数。

将E_u化到[0，2π]；

平近点角M_u＝E_u-e·sinE_u；

ΔM＝M_u-M₀；

将ΔM化到[0，2π]；

时间

t_u＝Δt+t₀；

其中M₀为参考时刻平近点角，t₀为参考时刻。Δ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)中计算太阳矢量在地心轨道坐标系中的投影具体为：

已知航天器近地点高度h_p为200km，远地点高度h_a为24371km，倾角i＝0°，升交点赤经Ω＝0°，近地点幅角ω＝0°，则太阳矢量在地心轨道坐标系中的投影为：

所述步骤(2)中计算阴影区在轨道平面内的方向角、阴影区投影椭圆半长轴，具体为：

(2-1)阴影区在椭圆轨道平面内的方向角θ_u，具体由公式：

给出。

(2-2)阴影区投影椭圆半长轴a_u，具体由公式：

给出。

所述步骤(3)中计算阴影区可能存在范围，具体为：

(3-1)所述阴影区可能存在范围两侧的真近点角f₁、f₂为：

f₁＝f(θ_u，a，e，R_e)＝222.2°；

f₂＝f(θ_u，a，e，-R_e)＝281.8°；

则真近点角f₁、f₂与方向角θ_u的距离分别为

df₁＝f₁-θ_u＝-17.8°；

df₂＝f₂-θ_u＝41.8°；

将df₁、df₂到[-π，π]；

(3-2)对阴影区可能存在范围两侧的真近点角进行排序：

df₁≤df₂；

阴影区可能存在范围ff＝[f₁，f₂]；

df＝[df₁，df₂]；

所述步骤(4)中对阴影区可能出现的位置进行分类，具体为：

阴影区的分类依据，具体由步骤：

r(θ_u)≤a_u；

阴影区类型为a类，如图3(a)所示；

所述步骤(5)中计算a类阴影区范围，具体为：

(5-1)计算数值迭代次数iter，设∈＝0.1°：

(5-2)对任意真近点角f，是否进入阴影区的判别方式，具体由公式：

给出。当则进入阴影区，当则不进入阴影区。

(5-3)用数值方法计算a类阴影区范围：根据(5-2)中的阴影区判别公式用二分法在区间[ff(1),θ_u]和[_u,f(2)]内寻找阴影区范围f_u＝[_u(1),f_u(2)]＝[279.0°265.2°]。

所述步骤(6)中计算b类阴影区范围，具体为：

本实施例不进入此分支。

所述步骤(7)将阴影区存在范围转换为由偏近点角E_u、平近点角M_u或时间t_u表示的参数，具体为：

阴影区所在的偏近点角范围为

将E_u化到[0，2]；

平近点角M_u＝E_u-e·sinE_u＝[317.3°307.1°]；

利用本发明一种适用于椭圆轨道的半解析阴影区预报方法，可预估一个轨道周期内阴影区的具体范围，使太阳能电推进卫星在获知阴影区位置后提前关机，避免阴影区电推力器开机引起星上其他载荷发生能源不足的现象。

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

Claims (6)

1.一种适用于椭圆轨道的半解析阴影区预报方法，其特征在于该方法包括如下步骤：

(1)将太阳矢量从惯性坐标系转换到地心轨道坐标系；

(2)依据阴影区桶形模型，获得阴影区在椭圆轨道平面内的方向角、阴影区投影椭圆半长轴；

(3)解析计算阴影区可能存在范围；

(4)根据阴影区投影椭圆和航天器所处轨道间的几何关系对阴影区可能出现的位置进行分类，并对阴影区的存在进行粗略判断，若阴影区存在范围属于a类，则进入步骤(5)，若阴影区存在范围属于b类，则进入步骤(6)，若不存在阴影区，结束本次预报；

(5)获取a类阴影区范围，在阴影区可能存在范围区域内，直接利用数值方法求得阴影区存在范围，进入步骤(7)；

(6)获取b类阴影区范围，缩小阴影区存在范围，搜索在此范围内进入阴影区的任意点，若在有限的迭代次数下不能搜到该点，则判定卫星在本轨道周期内不进入阴影，结束本次预报；若找到阴影点，根据该阴影点位置利用数值方法求得阴影区存在范围，进入步骤(7)；

(7)将阴影区存在范围转换为由偏近点角E、平近点角M以及时间表示的参数。

2.根据权利要求1所述的一种适用于椭圆轨道的半解析阴影区预报方法，其特征在于：所述步骤(2)中获得阴影区在轨道平面内的方向角、阴影区投影椭圆半长轴，具体为：

(2-1)阴影区在椭圆轨道平面内的方向角θ_u，具体由公式：

给出，其中，为太阳在中心天体轨道坐标系中的三维单位矢量，为的第一个分量，为的第二个分量；

(2-2)阴影区投影椭圆半长轴a_u，具体由公式：

给出，其中R_e为中心天体半径，为的第三个分量。

3.根据权利要求2所述的一种适用于椭圆轨道的半解析阴影区预报方法，其特征在于：所述步骤(3)中解析计算阴影区可能存在范围，具体为：

(3-1)所述阴影区可能存在范围两侧的真近点角f₁、f₂为：

f₁＝f(θ_u，a，e，R_e)；

f₂＝f(θ_u，a，e，-R_e)；

则真近点角f₁、f₂与方向角θ_u的距离分别为

df₁＝f₁-θ_u；

df₂＝f₂-θ_u；

将df₁、df₂化到[-π，π]；

其中，a为轨道半长轴，e为轨道偏心率，函数F＝f(θ_u，a，e，R)的具体方法为：

如果：θ_u＝0或θ_u＝π

则：L1＝-1/e；

L2＝a·(1-e²)/e；

L4＝a·(1-e²)；

L5＝R²e²+a²-2a²e²+a²e⁴；

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-e²)-R·e))；

f＝[L1，2π-L1]；

如果：θ_u≠0、θ_u≠π、

k＝tanθ_u；

R2＝R/sinθ_u；

L1＝-R2·k/(a(1-e²))；

L2＝k(a(1-e²)-e·R2)/(a(1-e²))；

L4＝R2·k²(a(1-e²)-e·R2)；

L5＝k²(a(1-e²)-e·R2)²+a²(1一e²)²；

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)对阴影区可能存在范围两侧的真近点角进行排序：

如果：df₁≤df₂

阴影区可能存在范围ff＝[f₁，f₂]；

df＝[df₁，df₂]；

如果：df₁＞df₂

阴影区可能存在范围ff＝[f₂，f₁]；

df＝[df₂，df₁]。

4.根据权利要求3所述的一种适用于椭圆轨道的半解析阴影区预报方法，其特征在于：所述步骤(4)中对阴影区可能出现的位置进行分类，具体为：

(4-1)如果：r(θ_u)≤a_u

则：阴影区类型为a类；

(4-2)如果：在阴影区可能存在范围ff中存在f^*，使r(f^*)≤a_u

则：可能存在b类阴影区；

(4-3)如果：在阴影区可能存在范围ff中不存在f^*，使r(f^*)≤a_u

则：不存在阴影区，退出本方法；

其中r(f)＝a(1-e²)/(1+e cosf)为任意真近点角所对应的轨道半径；

r(f^*)＝a(1-e²)/(1+e cosf^*)为真近点角为f^*时所对应的轨道半径，

r(θ_u)＝a(1-e²)/(1+e cosθ_u)为真近点角为θ_u时所对应的轨道半径。

5.根据权利要求4所述的一种适用于椭圆轨道的半解析阴影区预报方法，其特征在于：所述步骤(5)中获取a类阴影区范围，具体为：

(5-1)获取数值迭代次数iter，具体由公式：

给出，其中∈为容许阴影区的最大误差；

(5-2)对任意真近点角f，是否进入阴影区的判别方式，具体由公式：

给出；当则进入阴影区，当则不进入阴影区；

(5-3)用数值方法获得a类阴影区范围：根据(5-2)中的阴影区判别公式用二分法在区间[ff(1)，θ_u]和[θ_u，ff(2)]内寻找阴影区范围f_u＝[f_u(1)，f_u(2)]，迭代次数分别为iter(1)和iter(2)。

6.根据权利要求5所述的一种适用于椭圆轨道的半解析阴影区预报方法，其特征在于：所述步骤(6)中获取b类阴影区范围，具体为：

(6-1)在ff＝[f₁，f₂]中搜索一个进入阴影区的点，具体方法为：

如果：θ_u≥π

则：

其中i＝1，2，3，...，n，j＝1到2^i-1，n为正整数；

阴影区的判别公式为

如果：

则：计算结束；

f_mid＝f_tempi，j；

ff＝[θ_u，f₂]；

df＝ff-f_mid；

将df化到[-π，π]；

其中f_mid为中间变量，随i，j更新；

如果i＝1，2，3，...，n，始终则本轨道周期内不存在阴影区，计算结束；

如果：θ_u＜π

则：

其中i＝1，2，3，...，n，j＝1到2^i-1，n为正整数；

阴影区的判别公式为

如果：

则：计算结束；

f_mid＝f_tempi，j；

ff＝[f₁，θ_u]；

df＝ff-f_mid；

将df化到[-π，π]；

其中f_mid为中间变量，随i，j更新；

如果i＝1，2，3，...，n，始终则本轨道周期内不存在阴影区，计算结束；

(6-2)在ff中搜索b类阴影区范围，所述步骤为：

获取数值迭代次数

根据阴影区判别公式用二分法在区间[ff(1)，f_mid]和[f_mid，ff(2)]内寻找阴影区范围f_u＝[f_u(1)，f_u(2)]，迭代次数分别为iter(1)和iter(2)；其中，ff(1)为ff的第一个分量，ff(2)为ff的第二个分量。