《轨道力学讲义》——第四讲：轨道计算与预测

第四讲：轨道计算与预测

引言

在轨道力学的研究中，轨道计算与预测是将理论付诸实践的关键环节。当我们掌握了轨道运动的基本规律和数学描述后，下一步便是要能够准确地计算航天器在任意时刻的位置和速度，并对其未来的运动轨迹进行预测。这不仅是理论研究的延伸，更是航天任务规划与执行的基础。

想象一下，当一颗探测器被发射到火星，工程师们需要精确计算其数月后与火星的交会时间和位置；当一颗通信卫星需要保持在地球同步轨道上，控制中心需要随时掌握其轨道偏差并进行校正；当空间站需要与补给飞船对接，两个航天器之间的相对位置必须精确到厘米级——这些都依赖于轨道计算与预测技术。

本讲将深入探讨轨道计算的核心问题——开普勒方程的求解，以及轨道预测的各种方法和技术。我们将从基本的数学模型出发，逐步展开对复杂轨道问题的讨论，并通过实例说明这些理论和方法在实际航天工程中的应用。

一、开普勒方程求解

1.1 开普勒方程的由来与意义

在上一讲中，我们已经了解到，在二体问题中，天体沿椭圆轨道运动时，其位置与时间的关系可以通过开普勒方程来描述：

$$M=E-e\sin E$$

其中，$M$是平近点角（Mean Anomaly），代表时间的线性函数；$E$是偏近点角（Eccentric Anomaly），用于确定天体在轨道上的位置；$e$是轨道偏心率。

开普勒方程看似简单，却隐藏着深刻的物理内涵。它巧妙地连接了时间领域（平近点角$M$）与空间领域（偏近点角$E$），体现了开普勒第二定律中时间与扫过面积成正比的关系。当我们想知道航天器在某一时刻的位置时，我们通常已知轨道参数和时间（即$M$值），需要求解出$E$值，再通过轨道方程计算出航天器的位置坐标。

然而，开普勒方程的特殊性在于，方程中的未知数$E$同时出现在等式两边，且以超越函数$\sin E$的形式出现，这使得方程不存在解析解。换言之，我们无法直接写出$E$关于$M$的显式表达式。这一特性使得开普勒方程的求解成为轨道力学中一个经典的数学问题，也促使科学家们发展了多种数值和近似方法来解决这一问题。

1.2 牛顿迭代法求解开普勒方程

牛顿迭代法（Newton-Raphson Method）是求解非线性方程最常用的数值方法之一，特别适合于求解开普勒方程。其基本思想是：通过对函数进行线性近似，从一个初始猜测值开始，逐步逼近方程的真实解。

对于开普勒方程，我们可以定义函数：

$$f(E)=E-e\sin E-M$$

显然，当$f(E)=0$时，$E$就是开普勒方程的解。牛顿迭代法的迭代公式为：

$$E_{n+1}=E_n-\frac{f(E_n)}{f^{\prime }(E_n)}$$

对于开普勒方程，$f^{\prime }(E)=1-e\cos E$，因此迭代公式变为：

$$E_{n+1}=E_n-\frac{E_n-e\sin E_n-M}{1-e\cos E_n}$$

迭代过程可以概括为以下步骤：

1. 选择一个初始猜测值$E_0$。通常，对于小偏心率轨道，可以取$E_0=M$；对于高偏心率轨道，可以根据$M$的值选择更复杂的初始估计。
2. 利用迭代公式计算$E_1,E_2,...,E_n$，直到$|E_n-E_{n-1}|$小于预设的精度要求。

牛顿迭代法的收敛速度非常快，特别是当初始值选择得当时。对于大多数实际应用中的轨道（偏心率不太接近1），通常只需要3-5次迭代就能达到很高的精度。

以下是一个具体的例子：假设一颗卫星运行在偏心率$e=0.1$的轨道上，已知平近点角$M=0.5$弧度，求解偏近点角$E$。

初始猜测：$E_0=M=0.5$

第一次迭代：
$E_1=E_0-\frac{E_0-e\sin E_0-M}{1-e\cos E_0}=0.5-\frac{0.5-0.1\sin 0.5-0.5}{1-0.1\cos 0.5}=0.5+\frac{0.1\sin 0.5}{1-0.1\cos 0.5}\approx 0.5476$

第二次迭代：
$E_2=E_1-\frac{E_1-e\sin E_1-M}{1-e\cos E_1}\approx 0.5478$

可以看到，仅经过两次迭代，结果已经相当准确了。在实际的轨道计算中，通常会设定一个收敛条件，例如$|E_n-E_{n-1}|<10^{-8}$，以确保计算精度。

1.3 特殊情况下的解法与优化

尽管牛顿迭代法在大多数情况下都表现良好，但在某些特殊情况下可能会遇到问题：

1. 高偏心率轨道：当轨道偏心率接近1时（如彗星轨道），牛顿迭代法的收敛可能会变慢。这时，可以采用更复杂的初始值估计方法，或者使用其他改进的数值方法。

2. 接近近地点或远地点：当$M$接近0或$\pi$时，迭代可能需要更多步骤才能收敛。此时可以使用特殊的初始值策略。

3. 计算效率问题：在需要大量计算轨道的应用中（如航天器轨道的长期预报），开普勒方程求解的计算效率变得尤为重要。

为了解决这些问题，航天工程师们发展了多种技术：

• 改进的初始值估计：根据$M$和$e$的值范围，选择更接近真实解的初始猜测值。
• 变量代换技术：通过引入新变量，改变方程的数学特性，使迭代更加稳定。
• 查表法：对于特定轨道，预先计算并存储一系列$M$和$E$的对应值，在需要时通过插值快速获得近似解。
• 级数展开法：将解表示为无穷级数，通过截取有限项获得近似解。

1.4 级数展开法

对于偏心率较小的轨道，可以采用级数展开法直接求解开普勒方程。这种方法的基本思想是将$E$表示为$M$的级数：

$$E=M+\sum _{n=1}^{\infty }{C}_n(e)\sin (nM)$$

其中$C_n(e)$是偏心率$e$的函数。具体地，前几项为：

$$E=M+e\sin M+\frac{e^2}{2}\sin 2M+\frac{3e^3}{8}\sin 3M+\frac{e^4}{3}\sin 4M+...$$

在实际应用中，当偏心率较小（例如$e<0.1$）时，取前两三项就能获得相当高的精度。这种方法在计算效率上有明显优势，尤其适用于需要快速计算大量轨道点的情况。

例如，对于地球同步轨道的通信卫星，其偏心率通常很小（接近零），此时可以只取第一项近似：

$$E\approx M+e\sin M$$

这种简化大大减少了计算量，同时仍能保持较高精度。

二、轨道预报技术

2.1 轨道预报的基本概念

轨道预报是指根据已知的初始条件（位置和速度）和轨道力学模型，预测航天器在未来某一时刻的位置和速度的过程。它是航天任务规划与控制的基础，在卫星管理、深空探测、空间碎片监测等领域都有广泛应用。

轨道预报的精度取决于多个因素：

1. 初始条件的精度：任何测量都存在误差，航天器初始状态的测量误差会随时间传播和放大。
2. 力学模型的完备性：实际航天环境中存在各种摄动力，模型中包含的摄动力越全面，预报精度越高。
3. 计算方法的精度：不同的数值方法具有不同的截断误差和舍入误差特性。

根据所用方法的不同，轨道预报技术可分为解析预报方法和数值积分方法两大类。

2.2 解析预报方法

解析预报方法是基于轨道动力学方程的解析解或半解析解进行预报的方法。它的核心思想是将复杂的轨道动力学问题简化，找出关键变量之间的函数关系，从而避免繁重的数值计算。

2.2.1 开普勒轨道预报

最基本的解析预报方法是开普勒轨道预报，它建立在二体问题的精确解的基础上。其预报步骤如下：

1. 根据初始时刻的位置矢量${\vec{r}}_0$和速度矢量${\vec{v}}_0$，计算轨道的六个轨道根数：半长轴$a$、偏心率$e$、轨道倾角$i$、升交点赤经$\Omega$、近地点幅角$\omega$和初始时刻的真近点角${\nu }_0$（或偏近点角$E_0$、平近点角$M_0$）。

2. 计算预报时刻$t$对应的平近点角$M$：
   $$M=M_0+n(t-t_0)$$
   其中，$n=\sqrt{\frac{\mu }{a^3}}$是平均角速度（平均运动），$\mu$是中心天体的引力常数。

3. 求解开普勒方程，得到偏近点角$E$。

4. 计算真近点角$\nu$：
   $$\tan \frac{\nu }{2}=\sqrt{\frac{1+e}{1-e}}\tan \frac{E}{2}$$

5. 计算地心距$r$：
   $$r=\frac{a(1-e^2)}{1+e\cos \nu }$$

6. 在轨道平面坐标系中计算位置矢量，再通过坐标变换得到惯性坐标系中的位置矢量$\vec{r}$和速度矢量$\vec{v}$。

这种方法计算效率高，且对于短期预报（如地球轨道卫星的几天内预报）精度可以接受。然而，它忽略了各种摄动力的影响，因此长期预报精度会随时间迅速下降。

2.2.2 考虑摄动的解析预报

为了提高预报精度，可以在开普勒轨道预报的基础上，考虑主要摄动力的影响。常见的方法包括：

1. 平均根数理论：将轨道根数分解为长期项、周期项和短周期项，通过摄动方程计算各项的变化率，从而得到考虑摄动后的轨道根数演化。

2. SGP4/SDP4模型：这是一套广泛应用于双行根数（TLE）轨道预报的模型，能够较好地考虑地球非球形引力、大气阻力、太阳和月球引力等摄动影响。

3. 半解析理论：结合解析方法和数值方法的优点，如将短周期变化通过解析表达式处理，而长期演化通过数值积分求解。

例如，考虑地球扁率摄动时，轨道升交点赤经$\Omega$和近地点幅角$\omega$的变化率为：

$$\dot{\Omega }=-\frac{3}{2}{J}_2{\left(\frac{R_E}{p}\right)}^{2}n\cos i$$

$$\dot{\omega }=\frac{3}{4}{J}_2{\left(\frac{R_E}{p}\right)}^{2}n(5{\cos }^{2}i-1)$$

其中，$J_2$是地球引力场的带谐系数，$R_E$是地球半径，$p=a(1-e^2)$是半通径。

在预报时，可以近似认为这些变化率在短时间内保持恒定，从而计算考虑摄动后的轨道参数：

$$\Omega ={\Omega }_0+\dot{\Omega }(t-t_0)$$

$$\omega ={\omega }_0+\dot{\omega }(t-t_0)$$

解析预报方法的优点在于计算效率高、物理意义明确；缺点是难以全面考虑复杂的摄动环境，精度有限。因此，它主要用于对精度要求不高的初步轨道设计和分析，或作为数值预报的初值估计。

2.3 数值积分方法

数值积分方法是通过直接数值求解航天器运动微分方程来进行轨道预报的方法。它的基本思想是：将连续的时间区间离散化，在每个离散时间点上，根据当前状态和力学模型计算航天器的加速度，然后通过数值积分算法推进到下一时间点。

航天器在空间中的运动可以用如下微分方程描述：

$$\frac{d^2\vec{r}}{d{t}^2}=-\frac{\mu }{r^3}\vec{r}+{\vec{a}}_{pert}$$

其中，$\vec{r}$是航天器的位置矢量，$\mu$是中心天体的引力常数，${\vec{a}}_{pert}$表示各种摄动力产生的加速度。

为了数值求解，通常将二阶微分方程转化为一阶微分方程组：

$$\frac{d\vec{r}}{dt}=\vec{v}$$

$$\frac{d\vec{v}}{dt}=-\frac{\mu }{r^3}\vec{r}+{\vec{a}}_{pert}$$

然后采用数值积分算法求解这个方程组。常用的数值积分方法包括：

1. 欧拉法：最简单的数值积分方法，但精度较低，在轨道力学中很少使用。

2. 龙格-库塔法：包括经典的四阶龙格-库塔法（RK4）和变步长的龙格-库塔-费尔贝格法（RKF45）等，兼顾精度和效率，在轨道预报中应用广泛。

3. 多步法：如Adams-Bashforth法和Adams-Moulton法，适合于光滑问题，计算效率高。

4. 辛积分法：特别适合于哈密顿系统，能够保持系统的能量守恒特性，有利于长期轨道预报的稳定性。

5. 高精度数值积分器：如Bulirsch-Stoer方法和Gauss-Jackson方法等，用于高精度轨道预报。

以四阶龙格-库塔法（RK4）为例，其迭代公式为：

$${\vec{y}}_{n+1}={\vec{y}}_n+\frac{h}{6}({\vec{k}}_1+2{\vec{k}}_2+2{\vec{k}}_3+{\vec{k}}_4)$$

其中，$\vec{y}=[\vec{r},\vec{v}{]}^T$是状态向量，$h$是时间步长，而：

$${\vec{k}}_1=\vec{f}({\vec{y}}_n,t_n)$$

$${\vec{k}}_2=\vec{f}({\vec{y}}_n+\frac{h}{2}{\vec{k}}_1,t_n+\frac{h}{2})$$

$${\vec{k}}_3=\vec{f}({\vec{y}}_n+\frac{h}{2}{\vec{k}}_2,t_n+\frac{h}{2})$$

$${\vec{k}}_4=\vec{f}({\vec{y}}_n+h{\vec{k}}_3,t_n+h)$$

其中，$\vec{f}(\vec{y},t)$是微分方程右端的函数，即：

$$\vec{f}(\vec{y},t)=\left[\begin{array}{c}\vec{v}\\ -\frac{\mu }{r^3}\vec{r}+{\vec{a}}_{pert}(t)\end{array}\right]$$

数值积分方法的优点是可以考虑各种复杂的摄动力，精度高；缺点是计算量大，且长期积分可能累积误差。在实际应用中，需要根据任务要求、预报时间跨度和计算资源等因素选择合适的积分方法和步长策略。

三、轨道确定

3.1 轨道确定的基本问题

轨道确定是指根据对航天器的观测数据，确定其轨道参数的过程。与轨道预报相反，轨道确定是一个"反问题"：从已知的"结果"（观测数据）推导"原因"（轨道参数）。这一过程在航天器发射后的初始轨道确定、在轨航天器的轨道维护、空间目标监视等方面都具有重要应用。

轨道确定面临的主要挑战包括：

1. 观测数据的不完备性：单次观测通常只能获得航天器状态的部分信息（如角度、距离或视线方向的速度），无法直接得到完整的位置和速度信息。
2. 观测误差的影响：所有测量都存在随机误差和系统误差，需要通过统计方法处理。
3. 动力学模型的不确定性：实际轨道受到各种摄动力的影响，而这些摄动力的精确建模存在困难。

根据处理方法和应用场景的不同，轨道确定可分为初始轨道确定和轨道改进两大类。

3.2 初始轨道确定

初始轨道确定（Initial Orbit Determination, IOD）是指在缺乏先验信息的情况下，仅使用有限的观测数据确定航天器轨道的过程。它通常用于新发射的航天器、新发现的空间目标或失联后重新获得联系的航天器。

初始轨道确定的方法多种多样，以下介绍几种经典方法：

3.2.1 Gauss方法

Gauss方法是最经典的初始轨道确定方法之一，适用于三次角度观测数据。它的基本思想是利用开普勒运动的几何特性，将角度观测转化为位置向量的确定问题。

Gauss方法的计算步骤如下：

1. 根据三次观测的时间$t_1$、$t_2$、$t_3$和观测方向${\hat{\rho }}_1$、${\hat{\rho }}_2$、${\hat{\rho }}_3$，建立观测方程：
   $${\vec{r}}_i={\vec{R}}_i+{\rho }_i{\hat{\rho }}_i(i=1,2,3)$$
   其中，${\vec{r}}_i$是航天器的地心位置向量，${\vec{R}}_i$是观测站的地心位置向量，${\rho }_i$是未知的航天器到观测站的距离。

2. 利用面积关系建立方程：
   $${\vec{r}}_1({\vec{r}}_2\times {\vec{r}}_3)+{\vec{r}}_2({\vec{r}}_3\times {\vec{r}}_1)+{\vec{r}}_3({\vec{r}}_1\times {\vec{r}}_2)=\vec{0}$$
   这一方程反映了三个位置向量共面的条件，同时考虑了开普勒第二定律中时间与面积的关系。

3. 通过迭代求解得到${\rho }_1$、${\rho }_2$、${\rho }_3$，进而确定${\vec{r}}_1$、${\vec{r}}_2$、${\vec{r}}_3$。

4. 选择其中一对位置向量（如${\vec{r}}_1$和${\vec{r}}_2$），利用朗伯特问题的解法计算对应的速度向量${\vec{v}}_1$或${\vec{v}}_2$。

5. 根据得到的位置和速度矢量，计算轨道六根数。

Gauss方法的优点是仅需要角度观测数据，适用于光学观测；缺点是计算复杂，对观测时间间隔有要求（太短或太长都会导致精度下降）。

3.2.2 Laplace方法

Laplace方法同样适用于角度观测数据，但其思路与Gauss方法不同。它的核心是将位置矢量对时间的展开式与观测方程结合，建立关于距离的线性方程。

具体步骤包括：

1. 将航天器位置矢量在中心时刻$t_2$附近展开为Taylor级数：
   $$\vec{r}(t)={\vec{r}}_2+{\dot{\vec{r}}}_2(t-t_2)+\frac{1}{2}{\ddot{\vec{r}}}_2(t-t_2{)}^2+...$$

2. 利用多次观测数据，通过数值微分计算方向向量$\hat{\rho }$及其一阶、二阶导数。

3. 结合二体运动方程$\ddot{\vec{r}}=-\frac{\mu }{r^3}\vec{r}$，建立关于${\rho }_2$、${\dot{\rho }}_2$和${\ddot{\rho }}_2$的方程组。

4. 求解方程组得到${\rho }_2$，进而确定${\vec{r}}_2$和${\dot{\vec{r}}}_2$。

5. 计算轨道六根数。

Laplace方法适合于观测弧段较短的情况，但要求观测数据足够密集，以便准确计算方向矢量的导数。

3.2.3 朗伯特问题解法

朗伯特问题是指：已知两个位置矢量${\vec{r}}_1$和${\vec{r}}_2$及其对应的时间间隔$\Delta t=t_2-t_1$，求解连接这两点的轨道。它在轨道交会设计和初始轨道确定中都有重要应用。

朗伯特问题的经典解法包括：

1. Universal变量法：引入通用变量$\chi$，将椭圆、抛物线和双曲线轨道统一处理，通过迭代求解超越方程。

2. Battin方法：一种基于连分式展开的高效算法，具有良好的数值稳定性。

3. Gooding方法：结合二分法和Newton法的混合算法，适用于各种轨道类型和转移角。

朗伯特问题的求解通常需要迭代计算，但现代算法已经能够高效且稳定地处理各种情况，包括近180°转移等特殊情况。

3.3 轨道改进

轨道改进（Orbit Improvement）是在已有初始轨道估计的基础上，利用更多观测数据对轨道参数进行精化的过程。它通常采用最小二乘法或贝叶斯估计等统计方法，处理冗余观测数据，提高轨道确定的精度和可靠性。

3.3.1 最小二乘法

最小二乘法是轨道改进中最基本也是最常用的方法。其基本思想是：选择一组轨道参数，使得观测数据与理论预报值之间的平方和误差最小。

假设有$m$个观测量组成的观测向量$\vec{y}$，与$n$个轨道参数组成的状态向量$\vec{x}$（通常$m>n$），两者之间的关系可表示为：

$$\vec{y}=\vec{h}(\vec{x})+\vec{\epsilon }$$

其中，$\vec{h}$是非线性观测方程，$\vec{\epsilon }$是观测误差。

为了使用线性最小二乘法，需要将观测方程线性化：

$$\vec{y}-\vec{h}({\vec{x}}_0)\approx \mathbf{H}(\vec{x}-{\vec{x}}_0)+\vec{\epsilon }$$

其中，${\vec{x}}_0$是初始状态估计，$\mathbf{H}=\frac{\partial \vec{h}}{\partial \vec{x}}{|}_{\vec{x}={\vec{x}}_0}$是雅可比矩阵（或称灵敏度矩阵）。

定义残差向量$\vec{r}=\vec{y}-\vec{h}({\vec{x}}_0)$和状态修正量$\Delta \vec{x}=\vec{x}-{\vec{x}}_0$，最小二乘解为：

$$\Delta \vec{x}=({\mathbf{H}}^T\mathbf{WH}{)}^{-1}{\mathbf{H}}^T\mathbf{W}\vec{r}$$

其中，$\mathbf{W}$是权重矩阵，通常取为观测误差协方差矩阵的逆：$\mathbf{W}={\mathbf{R}}^{-1}$。

最小二乘法的实现通常采用迭代方式：

1. 从初始估计${\vec{x}}_0$开始。
2. 计算观测残差$\vec{r}$和雅可比矩阵$\mathbf{H}$。
3. 求解正规方程，得到状态修正量$\Delta \vec{x}$。
4. 更新状态估计：${\vec{x}}_0\leftarrow {\vec{x}}_0+\Delta \vec{x}$。
5. 检查收敛条件，如果未收敛则返回步骤2。

3.3.2 卡尔曼滤波

卡尔曼滤波是一种递推的最优估计算法，特别适合于处理连续的观测数据流。相比传统的批处理最小二乘法，卡尔曼滤波能够实时地融合新的观测信息，更新状态估计及其不确定性，非常适合于在轨航天器的实时轨道确定。

卡尔曼滤波的基本思想是将状态估计表示为概率分布，并通过预测-校正的循环不断更新。其典型实现包括以下步骤：

1. 时间更新（预测）：
   $${\vec{x}}_{k|k-1}=\vec{f}({\vec{x}}_{k-1|k-1})$$
   $${\mathbf{P}}_{k|k-1}={\mathbf{F}}_{k-1}{\mathbf{P}}_{k-1|k-1}{\mathbf{F}}_{k-1}^T+{\mathbf{Q}}_{k-1}$$
   其中，$\vec{f}$是状态转移函数，$\mathbf{F}$是其雅可比矩阵，$\mathbf{P}$是状态协方差矩阵，$\mathbf{Q}$是过程噪声协方差矩阵。

2. 测量更新（校正）：
   $${\mathbf{K}}_k={\mathbf{P}}_{k|k-1}{\mathbf{H}}_k^T({\mathbf{H}}_k{\mathbf{P}}_{k|k-1}{\mathbf{H}}_k^T+{\mathbf{R}}_k{)}^{-1}$$
   $${\vec{x}}_{k|k}={\vec{x}}_{k|k-1}+{\mathbf{K}}_k({\vec{y}}_k-\vec{h}({\vec{x}}_{k|k-1}))$$
   $${\mathbf{P}}_{k|k}=(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k){\mathbf{P}}_{k|k-1}$$
   其中，$\mathbf{K}$是卡尔曼增益，$\vec{h}$是观测函数，$\mathbf{H}$是其雅可比矩阵，$\mathbf{R}$是观测噪声协方差矩阵。

对于非线性系统，常用的变种包括扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）。在航天器轨道确定中，由于动力学模型和观测模型的非线性性，通常采用EKF或UKF。

3.3.3 批处理与连续处理的比较

轨道改进方法可分为批处理和连续处理两类：

• 批处理方法（如最小二乘法）一次性处理所有观测数据，适合于事后分析和高精度轨道确定。
• 连续处理方法（如卡尔曼滤波）逐步融合新的观测信息，适合于实时轨道确定和导航。

两类方法各有优缺点：

• 批处理方法能够充分利用所有观测数据之间的相关性，理论上可获得更优的估计结果，但计算量大，不适合实时应用。
• 连续处理方法计算效率高，能够实时估计轨道参数，但可能受到序贯处理数据顺序的影响。

在实际应用中，常根据任务需求选择合适的方法，或将两类方法结合使用，如采用短弧批处理+滤波的混合策略。

四、实际案例分析

4.1 GPS卫星轨道计算与预测

全球定位系统（GPS）是现代社会不可或缺的基础设施，其核心是一个由约30颗卫星组成的星座。GPS卫星的轨道计算与预测是系统精确运行的关键。

4.1.1 GPS广播星历

GPS系统采用广播星历（Broadcast Ephemeris）的方式向用户提供卫星轨道信息。广播星历是一组描述卫星轨道的参数，由地面控制段计算并上传至卫星，再由卫星播发给用户。广播星历采用改进的开普勒轨道模型，包括：

1. 轨道参考参数：参考时刻$t_{oe}$、半长轴$\sqrt{a}$、偏心率$e$、轨道倾角$i_0$、升交点赤经${\Omega }_0$、近地点幅角$\omega$和平近点角$M_0$。

2. 摄动改正项：轨道倾角变化率$\dot{i}$、升交点赤经变化率$\dot{\Omega }$、调和改正项（包括轨道倾角、半径和幅角的正弦和余弦周期项）等。

用户接收到广播星历后，按照以下步骤计算卫星在任意时刻$t$的位置：

1. 计算自参考时刻以来的时间$\Delta t=t-t_{oe}$。

2. 考虑摄动效应，计算修正后的轨道参数：
   $$M=M_0+n\Delta t$$
   $$i=i_0+\dot{i}\Delta t+\Delta i$$
   $$\Omega ={\Omega }_0+(\dot{\Omega }-{\omega }_e)\Delta t$$
   其中，${\omega }_e$是地球自转角速度，$\Delta i$是轨道倾角的调和改正项。

3. 求解开普勒方程，得到偏近点角$E$。

4. 计算卫星在轨道平面坐标系中的位置，并通过坐标变换得到地球固连坐标系中的位置。

广播星历的精度通常在几米量级，足以满足大多数导航定位应用的需求。

4.1.2 GPS精密星历

对于需要更高精度的科学应用（如大地测量、地壳运动监测等），国际GNSS服务组织（IGS）提供精密星历（Precise Ephemeris）服务。精密星历是基于全球分布的监测站网络观测数据，通过后处理计算得到的卫星轨道和钟差信息。

精密星历采用表格形式，给出一系列离散时刻的卫星位置和速度向量，用户通过插值计算任意时刻的卫星位置。常用的插值方法包括拉格朗日插值和切比雪夫多项式插值。

精密星历的精度可达厘米级，但由于需要收集和处理全球观测数据，通常有一定的延迟（从数小时到数天不等）。根据延迟时间的不同，IGS提供多种精密星历产品：

• 超快速产品（Ultra-rapid）：部分基于预测，延迟最短，精度在10cm左右。
• 快速产品（Rapid）：延迟约17小时，精度在5cm左右。
• 最终产品（Final）：延迟约13天，精度在3cm左右。

4.2 深空探测任务中的轨道计算

深空探测任务（如火星探测、小行星探测等）对轨道计算与预测提出了更高的要求。这类任务通常涉及长时间跨度、多种天体引力场和复杂的机动策略。

4.2.1 火星探测任务轨道设计

以火星探测任务为例，其轨道计算通常分为以下几个阶段：

1. 发射段：从地球表面到地球停泊轨道的轨道计算，主要考虑火箭的推力特性和地球引力场。

2. 地火转移段：从地球到火星的霍曼转移轨道（或其变种）计算，需要考虑发射窗口选择、能量优化等问题。

3. 中途修正：根据深空测量数据，计算并执行轨道修正机动，校正轨道偏差。

4. 捕获段：计算火星捕获所需的减速量，以及捕获轨道的设计。

5. 环绕段：火星轨道的长期预报和维持，需考虑火星重力场的非球形性、火卫一和火卫二的引力摄动等。

在火星探测任务中，轨道计算面临的主要技术挑战包括：

• 测量数据的延迟和稀疏：由于通信延迟（单程光行时间可达数分钟至数十分钟）和有限的测控资源，测量数据往往稀疏且延迟大。
• 多体引力场的复杂性：需要同时考虑太阳、地球、火星等多个天体的引力作用。
• 轨道预报的长时间跨度：从发射到到达火星，通常需要数月至数年，轨道预报误差会随时间累积。
• 自主导航的需求：在关键阶段（如火星捕获），可能需要航天器具备一定的自主导航能力，以应对地面控制延迟的问题。

为了应对这些挑战，深空探测任务的轨道计算通常采用高精度的数值积分方法，结合详细的摄动模型（包括太阳辐射压、多体引力、非球形引力场等），同时使用统计轨道确定方法处理不确定性。

4.2.2 嫦娥工程轨道计算案例

中国的嫦娥工程是深空探测轨道计算的成功案例。以嫦娥五号采样返回任务为例，其任务涉及地月转移、环月轨道、月面着陆、月面起飞、月地转移和再入返回等多个阶段，每个阶段都需要精确的轨道计算。

在嫦娥五号任务中，轨道计算的关键技术包括：

1. 多目标约束的轨道优化：在满足发射能力、通信视线、照明条件等多种约束的前提下，优化轨道设计。

2. 高精度月球引力场模型：利用前期嫦娥任务积累的数据，建立精确的月球引力场模型，提高环月轨道的预报精度。

3. 自主导航与地基测量相结合：在关键阶段（如月面着陆和起飞），结合光学导航、激光测距等自主导航手段与地基测量数据，提高轨道确定精度。

4. 多源数据融合的轨道确定：综合利用测距、测速、干涉测量等多种测量数据，通过滤波和平滑技术提高轨道确定精度。

嫦娥五号任务的成功实施，证明了中国在深空探测轨道计算领域已经具备了世界先进水平的技术能力。

五、总结与展望

5.1 轨道计算与预测的发展趋势

随着航天技术的不断进步和应用需求的日益多样化，轨道计算与预测技术也在持续发展。未来的主要发展趋势包括：

1. 高精度轨道动力学模型：更精确地描述各种摄动力（如高阶地球引力场、太阳辐射压、大气阻力等）的作用，提高长期轨道预报的精度。

2. 人工智能辅助的轨道计算：利用机器学习和深度学习技术，改进轨道预报模型，特别是在处理复杂非线性因素（如大气密度变化）方面。

3. 分布式和自主导航系统：减少对地面测控系统的依赖，增强航天器自主导航和轨道控制能力，特别是对于深空和小行星等复杂环境。

4. 实时大数据处理：随着空间物体数量的增加和观测数据的爆炸式增长，发展高效的大数据处理技术，实现对大量空间目标的实时轨道监测和预警。

5. 量子计算在轨道力学中的应用：探索量子算法解决轨道优化和轨道确定中的复杂计算问题，提高计算效率。

5.2 轨道计算与航天工程的关系

轨道计算与预测作为航天工程的基础技术，与航天任务的各个方面都有紧密联系：

1. 任务设计与规划：轨道计算为任务设计提供理论基础，决定发射窗口、能量需求和任务时序等关键参数。

2. 航天器设计：轨道环境（如辐射环境、热环境）直接影响航天器的设计参数，而轨道计算能够预测这些环境条件。

3. 测控与导航：轨道预报为测控系统提供天线指向预报，为导航系统提供参考轨道。

4. 任务运行与管理：实时轨道确定是评估任务状态、制定应急预案的基础。

5. 任务扩展与延长：精确的轨道维持策略能够延长卫星寿命，拓展任务能力。

随着商业航天的蓬勃发展和"太空交通"的日益繁忙，轨道计算与预测技术将在空间安全、碰撞避免、空间态势感知等领域发挥更加重要的作用。

5.3 学习建议与研究方向

对于希望深入学习轨道计算与预测的学生，建议：

1. 打牢基础：深入理解牛顿力学、数值分析和统计估计的基本原理，这是轨道力学研究的基石。

2. 掌握工具：学习使用STK、GMAT等轨道分析软件，以及MATLAB、Python等编程工具，提高实践能力。

3. 关注前沿：定期阅读相关学术期刊和会议论文，了解轨道力学研究的最新进展。

4. 参与实践：积极参与卫星设计竞赛、立方星项目等实践活动，将理论知识应用于实际问题。

对于研究生和研究人员，以下方向值得关注：

1. 高精度轨道动力学建模：包括复杂天体引力场建模、大气密度精确预报、太阳辐射压细致建模等。

2. 实时轨道确定算法改进：如鲁棒性滤波算法、粒子滤波在轨道确定中的应用等。

3. 轨道优化与控制：低推力长时间机动的轨道优化、燃料最优的轨道维持策略等。

4. 空间态势感知：大规模空间目标的轨道监测、预警和分析技术。

5. 行星防御：小行星轨道精确预报与偏转技术研究。

思考题

1. 轨道计算的误差来源：轨道计算中的误差主要来自哪些方面？如何减小或控制这些误差？考虑初始条件测量、力学模型、计算方法等方面。

2. 数值积分与解析方法的权衡：在实际航天任务中，如何选择适当的轨道预报方法（解析方法或数值积分方法）？请从计算效率、精度要求、预报时间跨度等方面进行分析。

3. 轨道确定的观测策略：对于一个新发现的近地小行星，如何设计最优的观测策略来确定其轨道？考虑观测频率、时间跨度和观测方式等因素。

4. 机器学习在轨道计算中的应用：人工智能和机器学习技术能否应用于轨道计算与预测？请举例说明可能的应用场景和面临的挑战。

5. 相对论效应的影响：在何种情况下需要考虑相对论效应对轨道计算的影响？请给出具体例子并估算这种影响的量级。

6. 高精度轨道预报的关键技术：针对导航卫星（如GPS、北斗）的精密轨道预报，讨论哪些关键技术和模型是必不可少的。分析这些卫星的轨道预报精度对定位精度的影响机制。

7. 多体问题中的轨道计算：在地月系统或太阳-地球-月三体系统中，讨论圆限制性三体问题模型的适用性和局限性。分析拉格朗日点附近的轨道特性及其计算方法。

8. 轨道协方差分析：轨道确定中的协方差矩阵具有什么物理意义？讨论如何利用协方差信息评估轨道预报的不确定性，以及如何通过观测规划来有效降低这种不确定性。

9. 空间碎片轨道预报：空间碎片的轨道预报面临哪些特殊挑战？分析大气阻力、太阳活动以及碎片特性不确定性对轨道预报的影响，并提出可能的改进方法。

10. 开普勒方程数值解法比较：评估不同的开普勒方程数值解法（如牛顿法、连分式方法、Danby方法等）在计算效率和精度方面的优缺点。特别分析高偏心率轨道情况下的计算挑战。

11. 轨道确定的可观测性分析：讨论不同类型观测数据（角度、距离、多普勒等）对轨道确定的贡献。分析轨道确定中的可观测性问题，特别是对于深空探测器或地球同步卫星等特殊情况。

12. 混合坐标系统在轨道动力学中的应用：分析等变分（equinoctial）轨道要素和平根数（mean elements）等特殊坐标系在轨道计算与预测中的优势。尤其讨论它们在处理奇异性和摄动建模方面的特点。

习题

1. 数值算法题：一颗卫星运行在偏心率为0.2的椭圆轨道上，已知平近点角M = 30°，求解该卫星的偏近点角E和真近点角ν。分别使用牛顿迭代法和级数展开法求解，并比较两种方法的计算精度与效率。

2. 轨道预报题：一颗地球卫星在J2000.0历元的开普勒轨道根数为：半长轴a = 8000 km，偏心率e = 0.1，轨道倾角i = 60°，升交点赤经Ω = 30°，近地点幅角ω = 45°，平近点角M = 0°。
   (a) 计算该卫星在历元时刻的位置矢量和速度矢量（ECI坐标系）。
   (b) 考虑J2摄动，预测卫星在30天后的轨道根数变化。

3. 行星际任务设计题：某深空探测器从地球出发执行火星任务，已知其离开地球后的轨道能量为C3 = 15 km²/s²（C3是轨道能量的两倍）。
   (a) 计算探测器的日心转移轨道参数。
   (b) 假设发射日期为2022年7月15日，计算探测器到达火星的大致日期。
   © 探测器到达火星时，需要怎样的减速量才能进入环火轨道（200 km × 10000 km）？

4. 轨道确定题：地面雷达站对一颗未知航天器进行了三次观测，得到以下数据：

   • t₁ = 2023-01-01 12:00:00 UTC, 方位角 = 120°, 仰角 = 45°, 距离 = 1000 km
   • t₂ = 2023-01-01 12:10:00 UTC, 方位角 = 150°, 仰角 = 60°, 距离 = 800 km
   • t₃ = 2023-01-01 12:20:00 UTC, 方位角 = 180°, 仰角 = 75°, 距离 = 700 km
     雷达站位于北纬30°，东经105°，海拔500米。
     (a) 使用Gauss或Laplace方法确定航天器的初始轨道。
     (b) 评估所确定轨道的不确定性。

5. 滤波算法设计题：设计一个卡尔曼滤波器，用于处理GPS导航卫星的轨道确定问题。假设有地面测站提供伪距观测数据，卫星受到地球非球形引力场（至少包含J2项）和太阳光压的影响。
   (a) 给出系统状态方程和观测方程。
   (b) 分析滤波器的性能，包括收敛速度和稳态精度。

6. 数值模拟题：使用四阶龙格-库塔法模拟一颗地球卫星在J2摄动下的轨道演化。卫星初始轨道为近圆轨道，半长轴8000 km，偏心率0.01，轨道倾角98°。
   (a) 选择合适的步长，模拟卫星一个月的轨道。
   (b) 分析升交点赤经和近地点幅角的长期变化趋势。

7. 轨道预报模型实现题：编写一个程序实现SGP4/SDP4轨道预报模型，使用TLE（两行根数）作为输入，计算卫星在未来24小时内的位置。
   (a) 以国际空间站为例，比较预报结果与NASA公布的星历的差异。
   (b) 分析模型中对大气阻力和地球非球形引力场的处理方法。

8. 测量数据处理题：某卫星发射后，地面测站获得了以下观测数据（假设为完美无误差观测）：

   • 北京站（北纬40°，东经116°）：t = 2023-06-01 08:00:00 UTC, 方位角 = 135°, 仰角 = 30°
   • 悉尼站（南纬34°，东经151°）：t = 2023-06-01 08:00:00 UTC, 方位角 = 225°, 仰角 = 45°
     (a) 利用这些同时观测的数据，确定卫星的位置矢量。
     (b) 如果再增加一次观测，如何确定卫星的完整轨道？

9. 批处理估计算法题：实现一个批处理最小二乘轨道确定算法，处理一组光学角度观测数据（赤经、赤纬）。
   (a) 讨论初始轨道估计的方法选择。
   (b) 分析迭代过程中的收敛特性和可能的数值问题。
   © 评估不同观测数据权重设置对结果的影响。

10. 特殊轨道积分方法题：设计一个Encke方法的轨道积分器，用于长期预报一颗高椭圆轨道卫星的轨道演化。
    (a) 与直接数值积分方法相比，分析Encke方法的优势和适用条件。
    (b) 选择合适的摄动模型，包括地球非球形引力场、月球和太阳引力以及太阳辐射压。
    © 讨论积分步长控制策略。

11. 综合摄动分析题：某低轨卫星的初始轨道为圆轨道，高度500公里，倾角51.6°。考虑地球非球形引力场（J2-J6项）、大气阻力（标准大气模型）和太阳光压：
    (a) 编写程序模拟卫星一年的轨道演化，分析各轨道要素的变化特征。
    (b) 评估不同摄动力对轨道长期变化的相对贡献。
    © 预测卫星的轨道寿命，并分析影响寿命预测准确性的关键因素。

12. 行星际轨道优化题：实现一个行星际轨道计算程序，计算从地球到火星的最优转移窗口。
    (a) 使用圆锥曲线（patched conic）近似方法，计算2025-2030年间的发射窗口和所需ΔV。
    (b) 考虑发射能量、飞行时间和到达速度等因素，确定最优发射日期。
    © 分析得到的转移轨道在行星位置误差和发射偏差条件下的敏感性。

参考文献

1. Vallado, D. A. (2013). Fundamentals of Astrodynamics and Applications (4th ed.). Microcosm Press.

2. Battin, R. H. (1999). An Introduction to the Mathematics and Methods of Astrodynamics. AIAA Education Series.

3. Montenbruck, O., & Gill, E. (2012). Satellite Orbits: Models, Methods and Applications. Springer.

4. Curtis, H. D. (2020). Orbital Mechanics for Engineering Students (4th ed.). Butterworth-Heinemann.

5. 黄勇，航天器轨道力学与控制，科学出版社，2019年。

6. 庞之浩，航天器轨道理论，国防工业出版社，2014年。

7. 崔平远，空间轨道力学，哈尔滨工业大学出版社，2012年。

8. 饶炜，航天器轨道确定原理，北京航空航天大学出版社，2010年。