CN108279426B - 一种测控站至卫星星下点航路捷径的解析计算方法 - Google Patents

Abstract

本发明公开一种测控站至卫星星下点航路捷径的解析计算方法，包括：依据二体轨道理论对当前时间点的轨道参数进行变换，获取轨道六根数；将测控站的坐标转换为地心轨道升交点坐标系下的坐标；根据测控站在地心轨道升交点坐标系下的坐标计算伪航路点；根据伪航路点和基于中间轨道理论的解析算法进行外推计算得到卫星在伪航路点的轨道参数；将卫星在伪航路点的轨道参数返回第一步进行一次迭代计算得到真航路点；通过所述真航路点计算得到航路捷径值。本发明提供的计算方法，相较于现有技术，本发明能够在确保精度的前提下，克服现有积分算法计算时间过长的缺点，较为快速精确地计算得到航路捷径。

Description

一种测控站至卫星星下点航路捷径的解析计算方法

技术领域

本发明涉及航路捷径的计算。更具体地，涉及一种测控站至卫星星下点航路捷径的解析计算方法。

背景技术

地面测控站至卫星星下点轨迹的最小大地线长称为航路捷径。在航路捷径的计算过程中，一般是已知地面测控站的站址和当前时刻的卫星飞行参数，计算航路捷径。

现有的解决方法是采用积分方式，每隔一定时间步长，计算出卫星的星下点，然后计算测控站至该星下点的大地线长，对各星下点对应的大地线长进行比对，随着卫星由远及近再由近及远，大地线长具有由大变小，又由小变大的趋势，因此总可以找到最近的一点，此点即为航路点。但是采用积分算法时，需要对卫星的轨道参数进行多次外推，直至找到航路点，而且为保证计算精度，常常需要选择较小的计算步长，这导致计算时间较长而无法满足要求，尤其在大量计算任务时更为突出。因此有必要对提高航路捷径计算效率的算法进行研究。

发明内容

为达到上述目的，本发明采用下述技术方案：

本发明提供一种测控站至卫星星下点航路捷径的解析计算方法，包括：

依据二体轨道理论对当前时间点的轨道参数进行变换，获取轨道六根数；

将测控站的坐标转换为地心轨道升交点坐标系下的坐标；

根据测控站在地心轨道升交点坐标系下的坐标计算伪航路点；

根据伪航路点和基于中间轨道理论的解析算法进行外推计算得到卫星在伪航路点的轨道参数；

将卫星在伪航路点的轨道参数返回第一步进行一次迭代计算得到真航路点；

通过所述真航路点计算得到航路捷径值。

优选地，所述当前时间点的轨道参数包括当前时间点的位置矢量和速度矢量参数。

优选地，所述轨道六根数包括轨道倾角、升交点赤经、偏心率、半长轴、近地点幅角及平近点角。

优选地，所述轨道倾角采用轨道倾角公式计算，所述轨道倾角公式为：

所述升交点赤经采用升交点赤经公式计算，所述升交点赤经公式为：

所述偏心率采用偏心率公式计算，所述偏心率公式为：

所述半长轴采用半长轴公式计算，所述半长轴公式为：

所述近地点幅角采用近地点幅角公式计算，所述近地点幅角公式为：

所述平近点角采用平近点角公式计算，所述平近点角公式为：

M₀＝E-esinE；

其中，h_z为h的z轴分量，hx为h的x轴分量，hy为h的y轴分量；

其中，h的计算公式为：

其中，E为偏近点角，E的计算公式为：

其中，f为真近点角，f的计算公式为：f＝μ-ω，μ为升交点幅角，ω为近地点幅角，μ的计算公式为：

其中X为位置矢量的x轴分量，Y为位置矢量的y轴分量，Z为位置矢量的z轴分量。

优选地，所述测控站的坐标为84系坐标，所述转换测控站坐标的步骤包括：

将测控站的经纬高B0、L0、H0转换为84系直角坐标X0、Y0、Z0；

将X0、Y0、Z0转换为绕Z轴转第一角度得到X0.5、Y0.5、Z0.5，所述第一角度为升交点赤经减格林威治赤经；

将X0.5、Y0.5、Z0.5转换为绕X轴转第二角度得到X1、Y1、Z1，所述第二角度为轨道倾角。

优选地，所述计算伪航路点的步骤包括：

将测控站在地心轨道升交点坐标系下的坐标转换为大地纬度B1、大地经度L1以及大地高程H1；

确定测控站所在子午圈与卫星轨道面的交点为伪航路点。

优选地，所述大地纬度B1采用第一公式计算，所述第一公式为：

所述大地经度L1采用第二公式计算，所述第二公式为：

所述大地高程采用第三公式计算，所述第三公式为：

所述经度为L1，纬度为0的点为伪航路点。

优选地，所述计算得到航路捷径值包括计算地球椭球模型的大地线长的步骤。

本发明的有益效果如下：

本发明提供一种测控站至卫星星下点航路捷径的解析计算方法，依据二体轨道理论对当前时刻的轨道参数进行变换，得到伪航路点，从而通过外推算法和一次迭代计算进一步获得精确的航路点，相较于现有技术，本发明能够在确保精度的前提下，克服现有积分算法计算时间过长的缺点，较为快速精确地计算得到航路点。

附图说明

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

图1示出地心轨道升交点坐标系示意图。

图2示出本发明一个实施方式提供的一种测控站至卫星星下点航路捷径的解析计算方法流程图。

具体实施方式

为了更清楚地说明本发明，下面结合优选实施例和附图对本发明做进一步的说明。附图中相似的部件以相同的附图标记进行表示。本领域技术人员应当理解，下面所具体描述的内容是说明性的而非限制性的，不应以此限制本发明的保护范围。

定义：地心轨道升交点坐标系，原点为地球质心，Z轴垂直于轨道平面，X轴指向轨道升交点，Y轴满足右手坐标法则。如图1所示。

本发明提供一种测控站至卫星星下点航路捷径的解析计算方法，如图1所示，包括：

S1：依据二体轨道理论对当前时间点的轨道参数进行变换，获取轨道六根数。

具体的，利用二体轨道理论将当前时间点的位置速度参数进行处理，得到轨道半长轴、倾角、升交点赤经、近地点幅角和平近点角等六根数。同时，根据当前时间点计算出格林威治赤经。

输入参数为当前时刻的WGS84系下位置

速度

将其转换为J2000坐标系下的位置

速度

利用二体理论，将当前时刻的位置、速度转换为轨道六根数。

轨道倾角

升交点赤经

偏心率

半长轴

近地点幅角

f＝μ-ω (8)

平近点角M₀＝E-esinE (10)

其中，f为真近点角，μ为升交点幅角，ω为近地点幅角。

S2：将测控站的坐标转换为地心轨道升交点坐标系下的坐标。

具体的，首先将测控站的经纬高B0、L0、H0转换为84系直角坐标X0、Y0、Z0，然后进行两次坐标系转换。第一次转换为绕Z轴转一定角度，该角度为升交点赤经减格林威治赤经，第二次转换为绕X轴转一定角度，该角度为轨道倾角。经过坐标转换得到了地心轨道升交点坐标系下的直角坐标X1、Y1、Z1。

两次坐标转换：

R1＝L_x(i)L_z(θ)R0 (11)

R0为84系的直角坐标，L_z(θ)为绕Z轴进行的坐标转换，L_x(i)为绕X轴进行的坐标转换。R1为地心轨道升交点坐标系下的直角坐标。

其中i为轨道倾角。

其中θ＝Ω-α_g，Ω为轨道升交点赤经，α_g为当前时刻的格林威治赤经。

S3：根据测控站在地心轨道升交点坐标系下的坐标计算伪航路点。

上一步计算得到了测控站在地心轨道升交点坐标系下的直角坐标，将其转换为经纬高B1、L1、H1，测控站所在子午圈与卫星轨道面的交点即为航路点。因此，航路点在地心轨道升交点坐标系的经度为L1，纬度为0。由于卫星轨道面是基于二体轨道理论确定，故此航路点存在计算误差，称为“伪航路点”，该点在真实的航路点附近。

位置矢量R1：

地心距离r：r＝|R1| (15)

大地经度L1：

大地纬度B1：

大地高程H1：

经度为L1、纬度为0的点即为“伪航路点”。

S4：根据伪航路点和基于中间轨道理论的解析算法进行外推计算得到卫星在伪航路点的轨道参数。

卫星经过“伪航路点”时，在地心轨道升交点坐标系中具有同样的经度L1和纬度0，该时刻卫星的真近点角等于经度加上近地点幅角。结合卫星初始时刻的真近点角，可以得到卫星从初始位置至“伪航路点”上空的飞行时间。

基于中间轨道理论的解析外推算法为一种高精度轨道预报算法，可以精确计算得到卫星的轨道参数。

首先给出飞行时间的计算步骤：

卫星经过“伪航路点”时的真近点角：f1＝L1+ω (19)

根据式(9)和式(10)计算得到平近点角M1。

卫星在初始时刻的真近点角为f0，同理可以得到平近点角为M0。

飞行时间t_f计算方法：

根据转移时间，利用高精度解析预报算法计算得到飞行器在“伪航路点”上空的精确位置、速度。

其中，外推计算理论的步骤为：

1.将当前时刻J2000坐标系的矢量x(t_i)，转换为椭球坐标系

其中x(t_i)由x_i，y_i，z_i，

和

组成，并且

d＝(r_i ²-c²)+δ(2z_i+δ) (23)

r_i ²＝x_i ²+y_i ²+z_i ² (24)

D²＝(ρ_i ²+c²)(1-η_i ²) (28)

2.计算Jacobi前三个常数(α₁，α₃，α₂)

3.F(ρ)和G(η)的四次式

F(ρ)＝μ[c²p₀(1-S₀)+(ρ²+c²)(γ₀ρ2+2ρ-p₀)] (32)

F(ρ)＝μγ₁(γρ²+2ρ-p)(ρ²+2A₁ρ+B₁) (33)

G(η)＝μ[-p₀(1-S₀)+(1-η²)(p₀+2δη+c²γ₀η²)] (34)

G(η)＝μS₁p₀(S+2Pη-η²)(1+P₁η-Q₁η²) (35)

4.初始化六个积分系数(R₁，R₂，R₃，N₁，N₂，N₃)

计算R₁，R₂，R₃，需要A_k和W_k，其中A₀＝1，通过因数分解得到A₁及其它，W₀＝V₀＝W为真近地点角，W₁＝(W+eV₁)/p，V₁＝sinW。

N₂＝D₁[u+k₁T₂/2+k₁T₂/2+3k₁ ²T₄/8+5k₁ ³T₆/16] (40)

计算N₁，N₂，N₃，需要

C_1k，C_2k和T_k，其中

T₀＝u (46)

T₁＝1-cosu (47)

T_k＝[(k-1)T_k-2-cosusin^k-1u]/k,k＝2,...,6(48)

5.计算Jacobi后三个常数(β₁，β₂，β₃)

t_i+β₁＝R₁(ρ_i)+c²N₁(η_i) (49)

β₂＝-α₂R₂(ρ_i)+α₂N₂(η_i) (50)

β₃＝φ_i+c²α₃R₃(ρ_i)-α₃N₃(η_i) (51)

使用初始环境和已初始化的系数。设置：

β₁＝-T 时间

β₂＝ω 近地点幅角

β₃＝Ω 升交点赤经

6.使用Jacobi常数(α₁，α₃，α₂，β₁，β₂，β₃)替换运动方程，求解给定时间t_f时的ρ_f，η_f和φ_f

t_f+β₁＝R₁(ρ_f)+c2N₁(η_f) (52)

β₂＝-α₂R₂(ρ_f)+α₂N₂(η_f) (53)

β₃＝φ_i+c²α₃R₃(ρ_f)-α₃N₃(η_f) (54)

第一个运动方程是开普勒方程的通式。只要计算使用椭球坐标系，相应的值由下式计算：

其中

7.将时间t_f时椭球坐标系下的位置矢量X(t_f)，转换为J2000坐标系的矢量x(t_f)

其中

S5：将卫星在伪航路点的轨道参数返回第一步进行一次迭代计算得到真航路点。

上一步得到的轨道参数为卫星在航路点附近的精确轨道参数，返回第一步中完成数据更新，对第一步至第三步的内容进行二次计算，即可得到真实的航路点。

S6：通过所述真航路点计算得到航路捷径值。

将该航路捷径点从地心轨道升交点坐标系转换至84系中，得到84系下的经度和纬度，考虑地球椭球模型，计算出测控站至航路点的大地线长即为航路捷径值。

根据测控站和真航路点的经度、纬度，可以计算椭球模型下的大地线长，模型如下：

测控站的经纬高为LBH₀，航路点的经纬高为LBH₁。

LBH₀＝[L₀；B₀；H₀]，LBH₁＝[L₁；B₁；H₁]，经坐标转换分别计算表面点1、表面点2在WGS84系中的位置r₁、r₂，计算两者叉乘结果r_n＝r₁×r₂。

a)当|r_n|＜1时，执行以下步骤：

1)

如果|Angle₁|＞1，则令Angle₁＝sign(Angle₁)；

2)Angle_{1_2}＝acos(Angle₁)；

3)如果Angle_{1_2}＜π/2，执行：

——LatCore＝asin(Tmp₁)，

——

——R_{1_2}＝R₁*Angle_{1_2}。

否则执行：

如果r_z·r₁＝0且r_z·r₂＝0，其中r_z＝[0；0；1]，按执行步骤(1)计算点1与点2距离，否则执行步骤(2)：

(1)

如果|Angle₁|＞1，则令Angle₁＝sign(Angle₁)；

Angle_{1_2}＝acos(Angle₁)；R_{1_2}＝R_a*Angle_{1_2}。

(2)令N_i＝3600，按照步长π/N_i离散化[0,π]，得到角度序列Ang＝0:π/N_i:π后，计算得到序列rAng＝R_a·[1-e²·cos²(Ang)]，计算

得到两点距离

b)当|r_n|≥1时，r_x＝r_z×r_n，r_z＝[0；0；1]，执行以下步骤：

如果|r_x|＝0：

——

如果|Angle₁|＞1，则令Angle₁＝sign(Angle₁)；

——Angle_{1_2}＝acos(Angle₁)，R_{1_2}＝R_a*Angle_{1_2}。

否则：

——

如果|XA|＞1，则令XA＝sign(XA)；

——A_XA＝acos(XA)；

——

如果|XB|＞1，则令XB＝sign(XB)；

——A_XB＝acos(XB)；

——

如果|Angle₁|＞1，则令Angle₁＝sign(Angle₁)；

——Angle_{1_2}＝acos(Angle₁)；

——如果r₁(3)×r₂(3)≥0，执行步骤(1)，否则执行步骤(2)：

(1)A_min＝min(A_XA,A_XB)，A_max＝max(A_XA,A_XB)；

(2)如果(A_XA+Angle_{1_2})≥π，A_min＝A_XA，

A_max＝A_XA+Angle_{1_2}；否则A_min＝-A_XA，A_max＝A_XB；

——r_y＝r_n×r_x，Tmp₂＝r_y(3)/|r_y|，如果|Tmp₂|＞1，Tmp₂＝sign(Tmp₂)；

——LatCore＝asin(Tmp₂)，

——

——

——令N_i＝3600，按照步长π/N_i离散化[A_min,A_max]，得到角度序列Ang＝A_min:(A_max-A_min)/N_i:A_max后，计算得到序列

再计算得到两点距离

c)完成高程修正，输出点1与点2距离R_{1_2}：

H_min＝min(H₁,H₂)

dH＝|H₁-H₂|

本发明提供一种测控站至卫星星下点航路捷径的解析计算方法，依据二体轨道理论对当前时刻的轨道参数进行变换，得到伪航路点，从而通过外推算法和一次迭代计算进一步获得精确的航路点，相较于现有技术，本发明能够在确保精度的前提下，克服现有积分算法计算时间过长的缺点，较为快速精确地计算得到航路点。

以下以某测控站和空间飞行器为例，对本发明所提出方法的有益效果进行具体说明。

测控站的站址坐标见表1所示。

表1测控站的站址坐标

序号	经度(°)	纬度(°)	高程(m)
				1	46.5	100	0

卫星的当前轨道参数见表2所示。

表2卫星的轨道参数

分别采用积分算法和本算法进行了航路捷径的计算，与STK的计算结果进行了比对。见表3所示。

表3航路捷径计算精度比对

两种算法的计算时间见表4所示，积分算法中卫星星下点的时间间隔为1s。

表4计算时间比对

由表3、表4可知，本算法与积分算法的计算精度相当，但计算速度大大优于积分算法。

由上所述，可见该方法能够在保证计算精度的前提下，大幅提高航路捷径的计算速度。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定，对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (8)

1.一种测控站至卫星星下点航路捷径的解析计算方法，其特征在于，包括：

依据二体轨道理论对当前时间点的轨道参数进行变换，获取轨道六根数；

将测控站的坐标转换为地心轨道升交点坐标系下的坐标；

根据测控站在地心轨道升交点坐标系下的坐标计算伪航路点；

根据伪航路点和基于中间轨道理论的解析算法进行外推计算得到卫星在伪航路点的轨道参数；

将卫星在伪航路点的轨道参数返回第一步进行一次迭代计算得到真航路点；

通过所述真航路点计算得到航路捷径值。

2.根据权利要求1所述方法，其特征在于，所述当前时间点的轨道参数包括当前时间点的位置矢量和速度矢量参数。

3.根据权利要求1所述方法，其特征在于，所述轨道六根数包括轨道倾角、升交点赤经、偏心率、半长轴、近地点幅角及平近点角。

4.根据权利要求3所述方法，其特征在于，所述轨道倾角采用轨道倾角公式计算，所述轨道倾角公式为：

所述升交点赤经采用升交点赤经公式计算，所述升交点赤经公式为：

所述偏心率采用偏心率公式计算，所述偏心率公式为：

所述半长轴采用半长轴公式计算，所述半长轴公式为：

所述近地点幅角采用近地点幅角公式计算，所述近地点幅角公式为：

所述平近点角采用平近点角公式计算，所述平近点角公式为：

M₀＝E-esinE；

其中，h_z为h的z轴分量，hx为h的x轴分量，hy为h的y轴分量；

其中，h的计算公式为：

其中，E为偏近点角，E的计算公式为：

其中，f为真近点角，f的计算公式为：f＝μ-ω，μ为升交点幅角，ω为近地点幅角，μ的计算公式为：

其中X为位置矢量的x轴分量，Y为位置矢量的y轴分量，Z为位置矢量的z轴分量。

5.根据权利要求1所述方法，其特征在于，所述测控站的坐标为84系坐标，所述转换测控站坐标的步骤包括：

将测控站的经纬高B0、L0、H0转换为84系直角坐标X0、Y0、Z0；

将X0、Y0、Z0转换为绕Z轴转第一角度得到X0.5、Y0.5、Z0.5，所述第一角度为升交点赤经减格林威治赤经；

将X0.5、Y0.5、Z0.5转换为绕X轴转第二角度得到X1、Y1、Z1，所述第二角度为轨道倾角。

6.根据权利要求2所述方法，其特征在于，所述计算伪航路点的步骤包括：

将测控站在地心轨道升交点坐标系下的坐标转换为大地纬度、大地经度以及大地高程；

确定测控站所在子午圈与卫星轨道面的交点为伪航路点。

7.根据权利要求6所述方法，其特征在于，所述大地纬度采用第一公式计算，所述第一公式为：

所述大地经度采用第二公式计算，所述第二公式为：

所述大地高程采用第三公式计算，所述第三公式为：

所述经度为L1，纬度为0的点为伪航路点。

8.根据权利要求1所述方法，其特征在于，所述计算得到航路捷径值包括计算地球椭球模型的大地线长的步骤。

Images