卫星编队飞行指若干个小卫星按特定的编队飞行, 各卫星协同工作, 形成一个分布式的编队卫星群。与大卫星相比, 编队卫星群成本低、可靠性强、灵活性高、编队自由, 能完成大卫星所不能完成的特殊任务。因此, 编队卫星技术在导航、通信、观测、侦查等军事和民用方面都扮演着重要角色, 具有广阔的应用前景^[1]。

编队的控制是编队卫星群运行的重要前提, 包括编队调整与编队保持, 编队控制的基础则是卫星间的相对运动动力学关系。1960年, Clohessy等^[2]最早开展卫星飞行的相对运动动力学关系研究, 并针对空间两临近飞行器交会问题, 提出Clohessy-Wiltshire方程(C-W方程)。C-W方程是线性化、常系数的微分方程组, 形式简单, 存在解析解^[1], 在早期的卫星编队飞行和卫星空间交会的动力学与控制中得到广泛应用^[3-4]。然而, 万有引力与卫星轨道半径之间的关系是非线性的, 因此描述卫星编队飞行的动力学方程是非线性方程^[5]。相对于卫星的轨道半径, 编队卫星群的卫星间距很小, 若将万有引力引起的非线性项做泰勒展开, 高阶项部分数量级相对较小, 其线性部分占有主导地位, 因此, 卫星编队飞行的动力学方程呈现弱非线性的特点。在研究卫星编队飞行动力学时, 若采用C-W模型而忽略非线性项的影响, 势必带来较大的误差, 其后果是消耗更多的燃料以维持卫星编队, 影响功能的充分发挥。对卫星编队飞行动力学方程的积分, 应当充分利用其弱非线性的特点, 首先将其做泰勒展开后写成半线性的形式(即方程包含线性部分和非线性部分), 然后用精细积分法^[6-7]求解线性部分, 用合适的积分方法求解高阶项的非线性部分。目的是精确计算数量级较大的线性部分, 而由于高阶非线性部分的数量级较小, 误差也比较小。

针对线性常系数微分方程, 经典常微分方程理论的主要困难是指数矩阵的计算。1994年, 钟万勰等^[6-7]针对指数矩阵的计算提出精细积分法(precise integration method, PIM), 用矩阵的加法代替矩阵乘法, 巧妙地避免了舍入误差。精细积分法有多方面的优势, 如可以采用大步长, 计算精度高, 稳定性好^[8]。另外, Ashi等^[9]比较了6种常用的指数矩阵的求解方法, 其中矩阵分解(matrix decomposition)方法无论在大步长还是小步长都能精确计算指数矩阵, 但当矩阵独立特征向量个数小于维数时, 此方法的应用比较困难。精细积分法则不受此限制, 无论步长大小都能精确计算指数矩阵, 是一种有效的计算方法。

针对半线性微分方程, Hersch^[10]于1958年首次提出指数积分的思想。Certaine^[11]于1960年用常数变异公式, 首次给出指数Adams-Moulton方法。Hochbruck等^[12]于2010年综述了指数积分方法的发展, 并系统地总结了指数积分法的计算格式。由于指数积分方法的思想是精确地求解微分方程的线性部分, 而方程的刚性和高振荡性通常表现在其线性部分, 所以在指数积分方法提出的初期, 主要用于解决刚性问题或高振荡问题^[13]。在指数积分方法的计算格式中, 指数矩阵往往难以精确计算, 尤其在维数较大的情况下。这限制了指数积分方法的进一步发展, 因此在指数积分方法提出后的很长时间都没有得到重视。21世纪初, 由于新的针对指数矩阵的计算方法的提出^[9], 指数积分方法得到快速发展。目前指数积分方法在各个领域得到广泛应用, 如延时问题^[14]、高振荡问题^[15]、抛物型方程的求解^[16]、李群算法^[17]以及Hamilton系统的保辛算法的构造^[18]等。针对此方法的研究重点由计算格式的构造逐渐转移到其计算性能的研究, 如稳定性^[19-20]、收敛性^[21]等。

本文在精细积分法与指数积分方法的基础上, 提出精细指数积分方法, 并将其用于弱非线性方程的求解, 为弱非线性方程的求解提供一种有效的计算方法, 为卫星编队飞行动力学仿真提供一种精确的算法。

1 卫星编队飞行动力学方程

卫星编队飞行动力学主要研究伴随卫星(follo-wer satellite)相对于参考卫星(reference satellite)的相对运动与受力之间的关系。如图 1所示, 假设参考卫星在圆形轨道上运行, 轨道半径为r=|r|, O₂X₂Y₂Z₂为轨道坐标系, 其中O₂Y₂轴的方向与地心到参考卫星的方向一致, O₂X₂方向为参考卫星运动的反方向, O₂Z₂方向由右手螺旋定则确定。在轨道坐标系中, 伴随卫星的位置矢量为r₂=[x, y, z]^T, 在不考虑摄动力和控制力的情况下, 伴随卫星相对于参考卫星运动的Lagrange函数为

$ \begin{array}{l} L = \frac{m}{2}[{(\dot x-\dot \theta y-\dot \theta r)^2} + {(\dot y + \dot \theta x)^2} + {{\dot z}^2}] + \\ \;\;\;\;\;\frac{{\mu m}}{{\sqrt {{x^2} + {{(y + r)}^2} + {z^2}} }}, \end{array} $

(1)

其中, μ为地球引力常数, m为伴随卫星的质量, $\dot \theta $为参考卫星的定常轨道角速度, 等号右边第一项为伴随卫星的动能, 第二项为势能。由Lagrange方程可导出伴随卫星相对于参考卫星的动力学方程:

$ \left\{ \begin{array}{c} \ddot x- {{\dot \theta }^2}x- 2\dot \theta \dot y\\ \ddot y- {{\dot \theta }^2}y + 2\dot \theta \dot x - r{{\dot \theta }^2}\\ {\ddot z} \end{array} \right\} = - \frac{\mu }{{{{[{x^2} + {{(y + r)}^2} + {z^2}]}^{\frac{3}{2}}}}}\left\{ \begin{array}{c} x\\ y + r\\ z \end{array} \right\}\; $

(2)

方程(2)的非线性项在于万有引力引起的方程右端的分母。卫星编队飞行的距离通常比参考卫星的轨道半径小得多, 即$\left| {{\boldsymbol{r}_2}} \right|=\sqrt {{x^2} + {y^2} + {z^2}} \ll r$, 故方程(2)右端的分母部分表现为弱非线性。为了简化方程, 寻求解析解, 一般可将其做泰勒展开并保留线性项, 可得

$ \left\{ \begin{array}{c} {\ddot x}\\ {\ddot y}\\ {\ddot z} \end{array} \right\} = \left\{ \begin{array}{c} 2\dot \theta \dot y\\ 3{{\dot \theta }^2}y-2\dot \theta \dot x\\ -{{\dot \theta }^2}z \end{array} \right\}。 $

(3)

由于方程的线性化会引入误差, 当误差积累到一定程度时需要对卫星施加控制, 从而消耗更多的燃料以维持编队, 甚至影响编队卫星群功能的实现。在求解卫星编队飞行的非线性方程时, 为了方便精细指数积分方法的使用, 将方程(2)写成半线性的形式:

$ \begin{array}{l} \left\{ \begin{array}{c} {\ddot x}\\ {\ddot y}\\ {\ddot z} \end{array} \right\} = \left\{ \begin{array}{c} 2\dot \theta \dot y\\ 3{{\dot \theta }^2}y- 2\dot \theta \dot x\\ - {{\dot \theta }^2}z \end{array} \right\} + \left\{ \begin{array}{c} {{\dot \theta }^2}x\\ - 2{{\dot \theta }^2}y + r{{\dot \theta }^2}\\ {{\dot \theta }^2}z \end{array} \right\} - \\ \;\;\;\;\;\;\;\;\;\frac{\mu }{{{{\left[{{x^2} + {{(y + r)}^2} + {z^2}} \right]}^{\frac{3}{2}}}}}\left\{ \begin{array}{c} x\\ y + r\\ z \end{array} \right\}。 \end{array} $

(4)

方程(4)右端的第一项为方程(2)的线性部分; 第二项和第三项合起来为方程(2)的高阶非线性部分, 由于方程表现为弱非线性, 该部分的数值比线性部分小若干个数量级。

2 精细积分法

对于卫星编队飞行的线性化动力学方程(式(3)), 可以利用精细积分法, 计算得到伴随卫星相对于参考卫星的近似运动规律。精细积分法的详细推导过程^[6-7]如下。

将方程(3)降阶成一阶微分方程组:

$ \boldsymbol{\dot u} = \boldsymbol{Lu}, $

(5)

其中, $\boldsymbol{u}={\left[{x, y, z, \dot x, \dot y, \dot z} \right]^{\rm{T}}}$, L为常数矩阵。将时间等步长离散, 步长为τ, 即t_k=kτ, u_k=u(t_k), k=0, 1, 2, …, 则有

$ {\boldsymbol{u}_{k + 1}} = \boldsymbol{T}{\boldsymbol{u}_k}, \;\boldsymbol{T} = \exp (\boldsymbol{L}\tau ) $

(6)

若能精确求解矩阵T, 则可以精确求解常系数齐次线性微分方程(5)。精细积分法的核心正在于矩阵T的精确求解, 求解思路是利用指数函数的加法原理, 将时间步长τ分成2ⁿ等分, 令$\Delta t=\frac{\tau }{{{2^n}}}$, 则有

$ \boldsymbol{T} = \exp (\boldsymbol{L}\tau ) = {[\exp (\boldsymbol{L}\Delta t)]^{{2^n}}}。 $

(7)

由于$\Delta t$是一个很小的时间段, 所以在exp (LΔt)的泰勒展开式中直接取其前5项, 即可达到足够的精度, 即

$ \left\{ \begin{array}{l} \exp (\boldsymbol{L}\Delta t) \approx \boldsymbol{I} + {\boldsymbol{T}_1}, \\ {\boldsymbol{T}_1}{\rm{ = }}\boldsymbol{L}\Delta t + \frac{1}{{2!}}{(\boldsymbol{L}\Delta t)^2} + \\ \;\;\;\;\;\frac{1}{{3!}}{(\boldsymbol{L}\Delta t)^3} + \frac{1}{{4!}}{(\boldsymbol{L}\Delta t)^4}, \end{array} \right. $

(8)

其中, I为单位矩阵。与I相比, T₁是一个很小的矩阵。为防止计算过程中“大数吃小数”, 要避免T₁与I直接相加。由方程(7)和(8)可知:

$ \begin{array}{c} \boldsymbol{T}{\bf{ = }}{(\boldsymbol{I}{\bf{ + }}{\boldsymbol{T}_1})^{{2^n}}}\\ = {[(\boldsymbol{I}{\bf{ + }}{\boldsymbol{T}_1})(\boldsymbol{I}{\bf{ + }}{\boldsymbol{T}_1})]^{{2^{n -1}}}}\\ = {(\boldsymbol{I} + 2{\boldsymbol{T}_1} + {\boldsymbol{T}_1}^2)^{{2^{n -1}}}}。 \end{array} $

(9)

若有矩阵序列:

$ {\boldsymbol{T}_{i + 1}} = 2{\boldsymbol{T}_i} + {\boldsymbol{T}_i}^2\;\;(i = 1, \;2, \;3, \;...), $

则根据式(9)有如下递推公式:

$ \begin{matrix} \boldsymbol{T}={{(\boldsymbol{I}+{{\boldsymbol{T}}_{1}})}^{{{2}^{n}}}}={{(\boldsymbol{I}+{{\boldsymbol{T}}_{2}})}^{{{2}^{n-1}}}}={{(\boldsymbol{I}+{{\boldsymbol{T}}_{3}})}^{{{2}^{n-2}}}} \\ =\cdots ={{(\boldsymbol{I}+{{\boldsymbol{T}}_{n+1}})}^{{{2}^{0}}}} 。 \\ \end{matrix} $

(10)

通过求此矩阵序列得到T_n+1, 最后与单位矩阵I相加, 便可计算矩阵T的值。在此过程中, 巧妙地运用指数函数的性质, 避免了严重的舍入误差, 所以这种算法有很高的精度, 一般取n=20可保证算法良好的性能。由于计算过程中在每个时间步长τ内插入2²⁰个点, 所以即使步长比较大, 也不影响算法的精确性。精细积分法求解一般的常系数齐次线性微分方程, 能够达到计算机意义上的精度。

3 精细指数积分方法

精细积分法只能处理卫星编队飞行的线性化动力学方程。为了更加精确地求解伴随卫星相对于参考卫星的动力学规律, 还需考虑其非线性因素。将半线性微分方程(4)降阶, 且写成如下形式:

$ \boldsymbol{\dot u} = \boldsymbol{Lu} + \boldsymbol{N}(\boldsymbol{u}, \;t), $

(11)

其中, L与方程(5)中的L相同, N(u, t)为方程的非线性部分, 当方程为弱非线性时有$\left\| {\boldsymbol{N} (\boldsymbol{u}, \; t)} \right\| \ll \left\| {\boldsymbol{Lu}} \right\|$。

对时间离散后, 方程(11)的精确解可写为如下形式:

$ {\boldsymbol{u}_{k + 1}} = \exp (\boldsymbol{L}\tau ){\boldsymbol{u}_k} + \int_{{t_k}}^{{t_{k + 1}}} {\exp \left[{\boldsymbol{L}({t_{k + 1}}-s)} \right]\boldsymbol{N}(\boldsymbol{u}(s), s)ds}, $

(12)

其中, 右端的积分项称为Duhamel积分。用不同的方法求解, 可得到不同的指数积分方法, 目前应用较广的是指数Runge-Kutta方法。s级的指数Runge-Kutta方法格式如下:

$ \left\{ \begin{array}{l} {\boldsymbol{u}_{k + 1}} = \exp (\boldsymbol{L}\tau ){\boldsymbol{u}_k} + \tau \sum\limits_{i = 1}^s {{\boldsymbol{b}_i}(\boldsymbol{L}\tau ){\boldsymbol{K}_i}}, \\ {\boldsymbol{K}_i} = \boldsymbol{N}\left( {\exp ({c_i}\boldsymbol{L}\tau ){\boldsymbol{u}_k} + \tau \sum\limits_{j = 1}^s {{\boldsymbol{a}_{ij}}(L\tau ){\boldsymbol{K}_i}}, {t_k} + {c_i}\tau } \right), \end{array} \right. $

(13)

其中, a_ij(Lτ)和b_i (Lτ)均为待定矩阵, 且为L和τ的函数, c_i为待定常数。当N (u, t)=0时, 指数Runge-Kutta方法退化为常系数齐次线性微分方程的精确解, 当L=0时, 此方法退化为经典的Runge-Kutta方法。

根据a_ij(Lτ)和b_i(Lτ)给出的方法不同, 指数Runge-Kutta方法可以分成不同的类型, 其中包括两种最常用的类型:一种是积分因子法(integration factor, IF), 也称为Lawson方法^[22], 这种方法的待定矩阵完全根据已有的Runge-Kutta方法确定, 构造和应用都比较简单; 另一种是指数时间微分方法(exponential time differencing, ETD), 此方法通过常数变异公式, 将非线性部分用代数插值多项式代替, 需要定义φ函数^{[12, 19]}。本文只对积分因子法做简单介绍。

积分因子法中的a_ij(Lτ)和b_i(Lτ)按如下规律^[19]给出:

$ \left\{ \begin{array}{l} {\boldsymbol{a}_{ij}}(\boldsymbol{L}\tau ) = {a_{ij}}\exp [({c_i}-{c_j})\boldsymbol{L}\tau], \\ {\boldsymbol{b}_i}(\boldsymbol{L}\tau ) = {b_i}\exp [(1-{c_i})\boldsymbol{L}\tau], \end{array} \right. $

(14)

其中, a_ij, b_i, c_i与经典Runge-Kutta方法的系数相同, 即给定一种经典的Runge-Kutta方法就有一种积分因子法与之对应, 且当N(u, t)=0时积分因子法退化为相应的经典的Runge-Kutta方法。例如, 经典四级四阶的Runge-Kutta方法所对应的积分因子法如下:

$ \left\{ \begin{array}{l} {\boldsymbol{u}_{k + 1}} = {\boldsymbol{T}_a}{\boldsymbol{u}_k} + \frac{\tau }{6}({\boldsymbol{T}_a}{\boldsymbol{K}_1} + 2{\boldsymbol{T}_b}{\boldsymbol{K}_2} + 2{\boldsymbol{T}_b}{\boldsymbol{K}_3} + {\boldsymbol{K}_4}), \\ {\boldsymbol{K}_1} = \boldsymbol{N}({\boldsymbol{u}_k}, \;{t_k}), \\ {\boldsymbol{K}_2} = \boldsymbol{N}\left( {{\boldsymbol{T}_b}{\boldsymbol{u}_k} + \frac{\tau }{2}{\boldsymbol{T}_b}{\boldsymbol{K}_1}, \;{t_k} + \frac{\tau }{2}} \right), \\ {\boldsymbol{K}_3} = \boldsymbol{N}\left( {{\boldsymbol{T}_b}{\boldsymbol{u}_k} + \frac{\tau }{2}{\boldsymbol{K}_2}, \;{t_k} + \frac{\tau }{2}} \right), \\ {\boldsymbol{K}_4} = \boldsymbol{N}({\boldsymbol{T}_a}{\boldsymbol{u}_k} + \tau {\boldsymbol{T}_b}{\boldsymbol{K}_3}, \;{t_k} + \tau ), \end{array} \right. $

(15)

其中, ${\boldsymbol{T}_a}=\exp (\boldsymbol{L}\tau), \; {\boldsymbol{T}_b}=\exp \left ({\frac{{\boldsymbol{L}\tau }}{2}} \right)$, 可由精细积分法精确计算。在其他文献(如文献[12, 19-20])中, 积分因子法的书面格式有所不同。本文将其写成式(15)的格式, 一方面为了精细积分法的应用, 另一方面为了编写程序的方便, 本质上是一致的。本文将精细积分法与指数积分法结合得到的方法称为精细指数积分方法。

4 卫星编队飞行数值仿真

在本节中, 先通过一个弱非线性的算例, 验证本文提出的积分因子法的有效性, 然后将此方法用于卫星编队飞行动力学方程的求解。

4.1 积分因子法有效性的验证

考虑以下弱非线性微分方程:

$ \left\{ \begin{array}{l} \dot x = \left( {-2 + a} \right)x-a{y^2}, \\ \dot y = ax-y - a{y^2}, \end{array} \right.\;\;\;\;其中, \;\;\;\left\{ \begin{array}{l} {x_0} = 1\;, \\ {y_0} = 1\;。 \end{array} \right. $

(16)

方程有解析解:

$ \left\{ \begin{array}{*{35}{l}} x={{\text{e}}^{-2t}}, \\ y={{\text{e}}^{-t}} 。 \\ \end{array} \right. $

(17)

当a取值较小时, 方程的非线性部分数量级比较小, 此处取a=0.1。分别采用经典的四级四阶Runge-Kutta方法(RK4_4)、八级六阶Runge-Kutta方法(RK8_6)^[23]、三级六阶隐式Runge-Kutta方法(RK3_6)^[23]、与RK4_4对应的积分因子法(IF4_4)、与RK3_6对应的积分因子法(IF3_6), 在不同的时间步长下求解方程(16), 结果如图 2所示。

图 2为以上各种方法相对误差绝对值的最大值与积分时间步长之间的关系。由图 2可知, RK4_4与IF4_4、RK3_6与IF3_6基本上平行, 说明积分因子法与对应的Runge-Kutta方法具有相同的阶数。在相同的步长时, IF4_4比RK4_4、IF3_6比RK3_6的计算误差都低4个数量级左右, 说明在弱非线性情况下, 指数Runge-Kutta方法比Runge-Kutta方法在精度上有比较明显的优势。将IF4_4与RK3_6进行比较, 在较大步长时, IF4_4方法比RK3_6更精确; 在中等步长时, 两种方法精度相当; 而在较小步长时, RK3_6更有优势。

为了分析方程非线性部分的数值大小与各种方法求解相对误差之间的关系, 取步长为$\tau=0.1$, a在0.01~1之间变化, 同样用上述5种方法求解方程(16), 结果如图 3所示。

由图 3可知, 对于普通的Runge-Kutta方法, 方程非线性的强弱对求解结果的相对误差影响不大, 然而, 精细指数积分方法求解结果的相对误差随非线性部分的数值增大而增加。方程的非线性越弱, 精细指数积分方法求解结果越准确。与同阶的Runge-Kutta方法相比, 精细指数积分方法在精度上有优势。

4.2 卫星编队调整动力学仿真

在卫星编队飞行的非线性动力学方程中, 对于Runge-Kutta方法, 方程(2)与方程(4)没有本质上的差别, 为了便于计算, 可直接积分方程(2);对于指数Runge-Kutta方法, 先将方程写成半线性形式(即方程(4)), 再进行计算。由于非线性动力学方程没有解析解, 所以用RK3_6和RK8_6两种高阶算法的计算结果作为参考, 主要比较RK4_4与IF4_4的计算结果。另外, 用精细积分法(PIM)求解线性化的动力学方程(方程(3)), 研究方程的线性化带来的误差。

假设在一次编队调整过程中, 参考卫星的轨道半径为8000 km, 伴随卫星相对于参考卫星的初始状态为

$ \begin{array}{l} {\boldsymbol{u}_0} = {\left[{x, y, z, \dot x, \dot y, \dot z} \right]^{\rm{T}}}\\ \;\;\;\; = {\left[{10000, 0, 0, -1, 0, 0} \right]^{\rm{T}}} \end{array} $

积分步长为τ=500 s, 经过4小时的运行, 轨迹图如图 4所示。

图 4中PIM即为精细积分法求解线性化方程的结果, 与其他曲线均有明显差距。由图 4可知, 经过4小时运行后, 线性化方程的计算误差为10² m量级。从图 4的局部放大图中可看到, IF4_4的计算轨迹与RK3_6和RK8_6的计算轨迹基本上重合, 而RK4_4的误差为10 m量级。对比IF4_4与RK4_4的计算结果可知, 在卫星编队飞行动力学方程的求解中, 精细指数积分方法能提供更好的计算精度。

4.3 卫星编队保持动力学仿真

当参考卫星的轨道半径为8000 km时, 伴随卫星相对于参考卫星的初始状态为

$ \begin{array}{l} {\boldsymbol{u}_0} = {\left[{x, y, z, \dot x, \dot y, \dot z} \right]^{\rm{T}}}\\ \;\;\;\; = {\left[{1000, 0, 0, 0.125\dot \theta, 0, 0} \right]^{\rm{T}}}, \end{array} $

即可使伴随卫星相对于参考卫星的运动轨道为周期轨道。分别用RK3_6, RK4_4和IF4_4方法求解动力学方程, 取积分步长为τ=600 s, 积分时间为100小时, 得到伴随卫星相对于参考卫星的运动轨迹图(图 5)。

由图 5可知, 六阶的RK3_6方法和四阶的IF4_4方法能精确地模拟周期轨道, RK4_4方法求解的周期轨道则出现较大的误差。由此可知, 在弱非线性问题的求解上, 与同阶的Runge-Kutta方法对比, 指数积分方法能提供较高的精度, 因此在求解时可以采用较大的步长。

5 结论

本文将精细积分法用于指数积分方法中的指数矩阵的求解, 得到精细指数积分方法, 并将其用于弱非线性方程的求解。针对弱非线性问题, 例如卫星编队飞行的动力学方程, 精细积分法不再适用。如果忽略其弱非线性的性质, 直接用Runge-Kutta方法求解, 则线性部分也会引入误差。

精细指数积分方法能精确求解方程的线性部分, 并用合适的算法求解非线性部分, 保证了算法的精度。与同阶的Runge-Kutta方法对比, 方程的非线性越弱, 精细指数积分方法在精度上越有优势。在卫星编队飞行的动力学仿真中, 精细指数积分方法取得了比同阶的Runge-Kutta方法更精确的结果。