卫星轨道动力学基本理论

卫星轨道动力学基本理论

1 卫星运动的状态参数$x(t)$

卫星轨道动力学的基本理论基础包括万有引力和开普勒运动定律，如果仅考虑卫星与地球的二体运动，可以直接使用开普勒定律来描述，在二体问题基础上，摄动力被定义为除主导引力外所有导致轨道偏离的扰动因素，按作用机制可分为^[1]：

• 引力摄动：包含N体（如太阳/月球）引力、地球非球形引力及潮汐效应。
• 非引力摄动：涵盖太阳辐射压、大气阻力、地球反照辐射和红外辐射压力等非接触作用。

当考虑各种摄动力影响时的运动时，需要利用牛顿运动定律描述运动方程，卫星在惯性系下的运动满足：
$$\begin{array}{c}m\ddot{\vec{r}}=\vec{F}={\vec{F}}_g+{\vec{F}}_{ng}+{\vec{F}}_{emp}\end{array}$$

$m$为卫星质量；$\ddot{\vec{r}}$为卫星位置的二阶导数，即卫星运动的加速度向量；$\vec{F}$ 为卫星所受所有引力的合力向量；${\vec{F}}_g$为卫星所受的保守力的合力，包括地球与N体引力、地球非球形引力、潮汐引力和相对论效应摄动力；${\vec{F}}_{ng}$ 为卫星所受的非保守力合力，包括太阳光压力、大气阻力和地球反照辐射摄动力；${\vec{F}}_{emp}$ 为经验加速度力，用于吸收未模型化的摄动力，提高定轨精度。

对于求解二阶微分方程较为困难，在实际的定轨过程中将其化为两个一阶微分方程。因此将卫星运动方程表示为一阶微分方程的形式：

$$\begin{array}{c}\{\begin{array}{l}\dot{\vec{r}}=\vec{v}\\ \dot{\vec{v}}=\frac{\vec{F}}{m}={\vec{a}}_g+{\vec{a}}_{ng}+{\vec{a}}_{emp}\end{array}\end{array}$$

令$\vec{x}=(\vec{r},\dot{\vec{r}},p)$，其中$p$为必要的动力学参数如太阳光压参数，大气阻力系数等，维度为$n_p$。这些参数在整个定轨过程中一般为常数，所以$\dot{p}$一般为0，因此可以令$\vec{F}=(\dot{\vec{r}},\ddot{\vec{r}},0)$，因此可以得到$\dot{\vec{x}}=\frac{\mathrm{d}x}{\mathrm{d}t}=F(\vec{x},t)$。

当已知卫星$t_0$时刻的初始状态${\vec{x}}_0=({\vec{r}}_0,{\dot{\vec{r}}}_0,p_0)$，则卫星的运动方程为：

$$\begin{array}{c}\{\begin{array}{l}\dot{\vec{x}}(t)=F(\vec{x},t),t\ge t_0\\ \vec{x}(t_0)=x_0\end{array}\end{array}$$

这是一个含有初值的常微分方程问题，因此若已知某一参考时刻$t_0$卫星的状态参数$x_0$，并根据卫星关于时间的合力向量，便可以获得后续任意时刻$t$卫星的运动状态参数，具体计算公式为：

$$\begin{array}{c}\vec{x}(t)=\vec{x}(t_0)+{\int }_{t_0}^{t}F(\vec{x},t)dt\end{array}$$

2 状态转移矩阵$\Phi (t,t_0)$

假设在时刻$t$，卫星真实的运动状态参数为${\vec{x}}^{\ast }(t)$，状态预报参数为$\vec{x}(t)$，因此卫星运动参数在$t$时刻的预报误差为：

$$\begin{array}{c}\vec{e}=\vec{x}(t)-{\vec{x}}^{\ast }(t)\end{array}$$

对上式左右两边进行求导并泰勒展开可以得到：

$$\begin{array}{c}\begin{array}{rl}\dot{\vec{e}}(t)& =\dot{\vec{x}}(t)-{\dot{\vec{x}}}^{\ast }(t)\\ & =F^{\ast }(\vec{x},t)+\frac{\partial F(\vec{x},t)}{\partial x^{\ast }}(\vec{x}(t)-{\vec{x}}^{\ast }(t))-F^{\ast }(\vec{x},t)\\ & =\frac{\partial F(\vec{x},t)}{\partial x^{\ast }}(\vec{x}(t)-{\vec{x}}^{\ast }(t))\\ & =A(t)\vec{e}(t)\end{array}\end{array}$$

在上式中，定义了$A(t)=\frac{\partial F(\vec{x},t)}{\partial x^{\ast }}$其表示在$t$时刻，卫星所受所有作用力产生的加速度对当前时刻卫星状态参数的偏导数，其维度为$(6+n_p)\ast (6+n_p)$。其具体形式如下：

$$\begin{array}{c}A(t)=\frac{\partial F(\vec{x},t)}{\partial x^{\ast }}=\left[\begin{array}{ccc}\frac{\partial \dot{\vec{r}}}{\partial \vec{r}}& \frac{\partial \dot{\vec{r}}}{\partial \dot{\vec{r}}}& \frac{\partial \dot{\vec{r}}}{\partial p}\\ \frac{\partial \ddot{\vec{r}}}{\partial \vec{r}}& \frac{\partial \ddot{\vec{r}}}{\partial \dot{\vec{r}}}& \frac{\partial \ddot{\vec{r}}}{\partial p}\\ \frac{\partial \dot{p}}{\partial \vec{r}}& \frac{\partial \dot{p}}{\partial \dot{\vec{r}}}& \frac{\partial \dot{p}}{\partial p}\end{array}\right]=\left[\begin{array}{ccc}0& I& 0\\ \frac{\partial \ddot{\vec{x}}}{\partial \vec{x}}& \frac{\partial \ddot{\vec{x}}}{\partial \dot{\vec{x}}}& \frac{\partial \ddot{\vec{x}}}{\partial p}\\ 0& 0& 0\end{array}\right]\end{array}$$

上式中，虽然位置$\vec{x}$关于时间的导数为速度$\dot{\vec{x}}$，但因为(\vec{x}) 和 (\dot{\vec{x}}) 在状态向量中被当作独立变量来处理，所以(\frac{\partial \dot{\vec{x}}}{\partial \vec{x}} = 0) 是由状态空间的数学定义所决定的。他们之间真正的物理耦合（即位置如何影响运动）体现在加速度项 (\ddot{\vec{x}}) 对位置 (\vec{x}) 的偏导数 (\frac{\partial \ddot{\vec{x}}}{\partial \vec{x}}) 上，而这个项通常是非常复杂的，包含了所有的引力摄动力模型。所以$A(t)$第二行一般不为0。

同时容知速度$\dot{\vec{x}}$与$p$无关，所以$\frac{\partial \dot{\vec{x}}}{\partial p}=0$，另外因为参数$p$为常数，所以$A(t)$第三行为0。

当$A(t)$是时变时，线性微分方程$\dot{\vec{e}}(t)=A(t)\vec{e}(t)$的通解的一般形式为:

$$\begin{array}{c}\vec{e}(t)=\Phi (t,t_0)\vec{e}(t_0)\end{array}$$

其中，$\vec{e}(t_0)$为参考时刻$t_0$卫星状态参数的误差，$\vec{e}(t)$为参考时刻$t$卫星状态参数的误差，$\Phi (t,t_0)$为从参考时刻$t_0$到参考时刻$t$的状态偏差向量的状态转移矩阵，其维度为$(6+n_p)\ast (6+n_p)$，同时$\Phi (t,t_0)$需要满足：

$$\begin{array}{c}\frac{\partial }{\partial t}\Phi (t,t_0)=A(t)\Phi (t,t_0),\Phi (t_0,t_0)=I\end{array}$$

从状态$x(t_0)$转移到状态$x(t)$的角度思考$\Phi (t,t_0)$，其具体形式为：

$$\begin{array}{c}\Phi (t,t_0)=\left[\begin{array}{ccc}\frac{\partial \vec{r}}{\partial {\vec{r}}_0}& \frac{\partial \vec{r}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \vec{r}}{\partial p}\\ \frac{\partial \dot{\vec{r}}}{\partial {\vec{r}}_0}& \frac{\partial \dot{\vec{r}}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \dot{\vec{r}}}{\partial p}\\ \frac{\partial p}{\partial {\vec{r}}_0}& \frac{\partial p}{\partial {\dot{\vec{r}}}_0}& \frac{\partial p}{\partial p}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial \vec{r}}{\partial {\vec{r}}_0}& \frac{\partial \vec{r}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \vec{r}}{\partial p}\\ \frac{\partial \dot{\vec{r}}}{\partial {\vec{r}}_0}& \frac{\partial \dot{\vec{r}}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \dot{\vec{r}}}{\partial p}\\ 0& 0& I\end{array}\right]\end{array}$$
当$A(t)$是常数矩阵时（即$A(t)=A$），$\Phi (t,t_0)=e^{A(t-t_0)}$。

对公式8进行求导并利用公式6和8，可以得到：

$$\begin{array}{c}\begin{array}{rlr}\dot{\Phi }(t,t_0)\vec{e}(t_0)& =\dot{\vec{e}}(t)& \\ & =A(t)\vec{e}(t)& \\ & =A(t)\Phi (t,t_0)\vec{e}(t_0)& \end{array}\end{array}$$

因此可得：

$$\begin{array}{c}\dot{\Phi }(t,t_0)=A(t)\Phi (t,t_0)\end{array}$$

因为在定轨过程中动力学参数为常数，且加速度对位置和速度的偏导数可以进一步简化，其大致情况为：

• 加速度对位置的偏导数
  加速度对位置的偏导数$\frac{\partial \ddot{\vec{r}}}{\partial \vec{r}}$是重力位（引力位）$U(\vec{r})=-\frac{GM}{r}$的二阶导数矩阵，即重力梯度张量$\nabla (\nabla U)$（也称为Marussi张量）。其中加速度与重力位满足以下关系：
  $$\begin{array}{c}\nabla U=\nabla (-\frac{GM}{r})=GM\cdot (\frac{1}{r^2}\cdot \frac{\vec{r}}{r})=\frac{GM}{r^3}\vec{r}=-\ddot{\vec{r}}\end{array}$$
  所以
  $$\frac{\partial \ddot{\vec{r}}}{\partial \vec{r}}=\nabla \ddot{\vec{r}}=\nabla (-\nabla U)=-{\nabla }^{2}U$$

• 加速度对速度的偏导数
  加速度对速度的偏导数$\frac{\partial \ddot{\vec{r}}}{\partial \dot{\vec{r}}}$与所受力的形式有关：

  • 如果只考虑保守力，即做功与路径无关的力，因为加速度等于重力位梯度的负数，只与位置有关，如地球引力、第三体引力、潮汐力，位置与速度属于独立的变量，所以加速度与速度无关，因此 $\frac{\partial \vec{a}}{\partial \dot{\vec{r}}}=0$。
  • 如果考虑非保守的耗散力（如大气阻力），即做功与路径有关的力，那么加速度显式地依赖于速度，比如大气阻力，因此 $\frac{\partial \vec{a}}{\partial \dot{\vec{r}}}\ne 0$。同时需要注意太阳辐射压导致的加速度对速度的偏导数等于0，地球辐射压导致的加速度对速度的偏导数约等于0。

当卫星轨道高度较高时，可以忽略大气阻力的影响，所以可以进一步化简$A(t)$和$\Phi (t,t_0)$的表达：

$$\begin{array}{c}A(t)=\left[\begin{array}{cc}\frac{\partial \dot{\vec{r}}}{\partial \vec{r}}& \frac{\partial \dot{\vec{r}}}{\partial \dot{\vec{r}}}\\ \frac{\partial \ddot{\vec{r}}}{\partial \vec{r}}& \frac{\partial \ddot{\vec{r}}}{\partial \dot{\vec{r}}}\end{array}\right]=\left[\begin{array}{cc}0& I\\ \frac{\partial \ddot{\vec{x}}}{\partial \vec{x}}& 0\end{array}\right]\end{array}$$

$$\begin{array}{c}\Phi (t,t_0)=\left[\begin{array}{ccc}\frac{\partial \vec{r}}{\partial {\vec{r}}_0}& \frac{\partial \vec{r}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \vec{r}}{\partial p}\\ \frac{\partial \dot{\vec{r}}}{\partial {\vec{r}}_0}& \frac{\partial \dot{\vec{r}}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \dot{\vec{r}}}{\partial p}\end{array}\right]\end{array}$$

因此：
$$\begin{array}{c}\begin{array}{rl}\dot{\Phi }(t,t_0)& =A(t)\Phi (t,t_0)+\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& \frac{\partial \ddot{\vec{r}}}{\partial p}\end{array}\right]\\ & =\left[\begin{array}{cc}0& I\\ \frac{\partial \ddot{\vec{r}}}{\partial \vec{r}}& 0\end{array}\right]\left[\begin{array}{ccc}\frac{\partial \vec{r}}{\partial {\vec{r}}_0}& \frac{\partial \vec{r}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \vec{r}}{\partial p}\\ \frac{\partial \dot{\vec{r}}}{\partial {\vec{r}}_0}& \frac{\partial \dot{\vec{r}}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \dot{\vec{r}}}{\partial p}\end{array}\right]+\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& \frac{\partial \ddot{\vec{r}}}{\partial p}\end{array}\right]\\ & =\left[\begin{array}{ccc}\frac{\partial \dot{\vec{r}}}{\partial {\vec{r}}_0}& \frac{\partial \dot{\vec{r}}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \dot{\vec{r}}}{\partial p}\\ \frac{\partial \ddot{\vec{r}}}{\partial \vec{r}}\frac{\partial \vec{r}}{\partial {\vec{r}}_0}& \frac{\partial \ddot{\vec{r}}}{\partial \vec{r}}\frac{\partial \vec{r}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \ddot{\vec{r}}}{\partial \vec{r}}\frac{\partial \vec{r}}{\partial p}+\frac{\partial \ddot{\vec{r}}}{\partial p}\end{array}\right]\end{array}\end{array}$$

其中$A(t)$的维度是$6\ast 6$，$\Phi (t,t_0)$和$\dot{\Phi }(t,t_0)$的维度都为$6\ast (6+n_p)$。

3 求解状态参数

变微分方程的具体形式非常复杂,采用求解析解的方法去解具体表达式是不可能的,也是不必要的。如果已知参考时刻t0卫星状态参数，并且卫星运动的加速度对卫星状态参数的偏导数可以准确计算得到，通过数值积分就可以计算后续任意时刻$t$卫星位置和速度对 t0 时刻卫星初始状态参数的偏导数，而卫星运动的加速度对卫星状态参数的偏导数也是可以根据摄动力模型得到，因此利用数值方法称为公式15表示的变微分方程的最佳选择^[2]。

前两节完成了对状态参数以及转移矩阵微分方程的建立，可以进一步地构造状态向量$\vec{X}$与右函数$F(\vec{X},t)$:

$$\begin{array}{c}\vec{X}={\left[\begin{array}{cccccccc}\vec{r}& \dot{\vec{r}}& \frac{\partial \vec{r}}{\partial {\vec{r}}_0}& \frac{\partial \vec{r}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \vec{r}}{\partial p}& \frac{\partial \dot{\vec{r}}}{\partial {\vec{r}}_0}& \frac{\partial \dot{\vec{r}}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \dot{\vec{r}}}{\partial p}\end{array}\right]}^T\end{array}$$

$$\begin{array}{c}\begin{array}{rl}& F(\vec{X},t)=\frac{\mathrm{d}\vec{X}}{\mathrm{d}t}\\ & ={\left[\begin{array}{cccccccc}\dot{\vec{r}}& \ddot{\vec{r}}& \frac{\partial \dot{\vec{r}}}{\partial {\vec{r}}_0}& \frac{\partial \dot{\vec{r}}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \dot{\vec{r}}}{\partial p}& \frac{\partial \ddot{\vec{r}}}{\partial \vec{r}}\frac{\partial \vec{r}}{\partial {\vec{r}}_0}& \frac{\partial \ddot{\vec{r}}}{\partial \vec{r}}\frac{\partial \vec{r}}{\partial {\dot{\vec{r}}}_0}& \frac{\partial \ddot{\vec{r}}}{\partial \vec{r}}\frac{\partial \vec{r}}{\partial p}+\frac{\partial \ddot{\vec{r}}}{\partial p}\end{array}\right]}^T\end{array}\end{array}$$

状态向量$\vec{X}$的初值为：

$$\begin{array}{c}\begin{array}{rl}{\vec{X}}_0& ={\left[\begin{array}{cccccccc}{\vec{r}}_0& {\dot{\vec{r}}}_0& \frac{\partial {\vec{r}}_0}{\partial {\vec{r}}_0}& \frac{\partial {\vec{r}}_0}{\partial {\dot{\vec{r}}}_0}& \frac{\partial {\vec{r}}_0}{\partial p}& \frac{\partial {\dot{\vec{r}}}_0}{\partial {\vec{r}}_0}& \frac{\partial {\dot{\vec{r}}}_0}{\partial {\dot{\vec{r}}}_0}& \frac{\partial {\dot{\vec{r}}}_0}{\partial p}\end{array}\right]}^T\\ & ={\left[\begin{array}{cccccccc}{\vec{r}}_0& {\dot{\vec{r}}}_0& I& 0& 0& 0& I& 0\end{array}\right]}^T\end{array}\end{array}$$

因此便可以利用数值积分的方法，计算出任意时刻$t$卫星的位置和速度。

4 关于转移矩阵的思考

$\Phi (t,t_0)$虽然定义为状态误差的转移矩阵，但我们分析其形式，结合公式5和公式8，可以得到：
$$\begin{array}{c}\begin{array}{rl}\vec{x}(t)-{\vec{x}}^{\ast }(t)& =\Phi (t,t_0)[\vec{x}(t_0)-{\vec{x}}^{\ast }(t_0)]\\ \vec{x}(t)& ={\vec{x}}^{\ast }(t)+\Phi (t,t_0)[\vec{x}(t_0)-{\vec{x}}^{\ast }(t_0)]\\ \vec{x}(t)& =\Phi (t,t_0)\vec{x}(t_0)\end{array}\end{array}$$

可以发现状态参数的转移矩阵和误差的转移矩阵是一致的，这样在进行数值求解的过程中，可以只计算状态转移矩阵，然后利用状态转移矩阵和前一时刻的状态参数得到下一时刻的状态参数。同时这一性质在卡尔曼滤波算法中也至关重要。

引用：

1. 百度百科.摄动力[在线].https://baike.baidu.com/item/%E6%91%84%E5%8A%A8%E5%8A%9B/8735869.

2. [1]程军龙.北斗三号星间链路集中式自主定轨关键技术研究[D].武汉大学,2019.DOI:10.27379/d.cnki.gwhdu.2019.001592.