航天器编队飞行(3):笛卡尔坐标表述

1.建模

  本节中，相对位置用定义在主星上的Hill坐标系下的分量$\mathbf{\rho }=(x,y,z{)}^T$来表征
$${\mathbf{r}}_d={\mathbf{r}}_c+\mathbf{\rho }=\left(r_c+x\right){\hat{\mathbf{o}}}_r+y{\hat{\mathbf{o}}}_{\theta }+z{\hat{\mathbf{o}}}_h\text{\hspace{1em}(1)}$$

其中，$r_c$是主星的轨道半径长度。Hill系$\mathcal{O}$相对于惯性坐标系$\mathcal{N}$的角速度为$${\mathbf{\omega }}_{\mathcal{O}/N}=\dot{f}{\hat{\mathbf{o}}}_h\text{\hspace{1em}(2)}$$

其中,$f$是主星的真近点角。
相对惯性系求二阶导数，得
$$\begin{array}{rl}{\ddot{\mathbf{r}}}_d=& \left({\ddot{r}}_c+\ddot{x}-2\dot{y}\dot{f}-\ddot{f}y-{\dot{f}}^2\left(r_c+x\right)\right){\hat{\mathbf{o}}}_r\\ & +\left(\ddot{y}+2\dot{f}\left({\dot{r}}_c+\dot{x}\right)+\ddot{f}\left(r_c+x\right)-{\dot{f}}^{2}y\right){\hat{\mathbf{o}}}_{\theta }+\ddot{\mathbf{z}}{\hat{\mathbf{o}}}_h\end{array}\text{\hspace{1em}(3)}$$

考虑到轨道角动量的值$h=r_c^{2}f$是常数，则有
$$\dot{h}=0=2r_c{\dot{r}}_c\dot{f}+r_c^2\ddot{f}\text{\hspace{1em}(4)}$$

因此，真近点角的二阶导数可以被写为

$$\ddot{f}=-2\frac{{\dot{r}}_c}{r_c}\dot{f}\text{\hspace{1em}(5)}$$

求主星惯性位置${\mathbf{r}}_c=r_c{\hat{\mathbf{o}}}_r$相对惯性系的二阶导数，并考虑轨道运动方程，得$${\ddot{\mathbf{r}}}_c=\left({\ddot{r}}_c-r_c{\dot{f}}^2\right){\hat{\mathbf{o}}}_r=-\frac{\mu }{r_c^3}{\mathbf{r}}_c=-\frac{\mu }{r_c^2}{\hat{\mathbf{o}}}_r\text{\hspace{1em}(6)}$$

展开分量形式，得
$${\ddot{r}}_c=r_c{\dot{f}}^2-\frac{\mu }{r_c^2}=r_c{\dot{f}}^2\left(1-\frac{r_c}{p}\right)\text{\hspace{1em}(7)}$$

把(5)和(7)代入(3)式，得到从星的加速度矢量表达式$$\begin{array}{rl}{\ddot{\mathbf{r}}}_d=& \left(\ddot{x}-2\dot{f}\left(\dot{y}-y\frac{{\dot{r}}_c}{r_c}\right)-x{\dot{f}}^2-\frac{\mu }{r_c^2}\right){\hat{\mathbf{o}}}_r\\ & +\left(\ddot{y}+2\dot{f}\left(\dot{x}-x\frac{{\dot{r}}_c}{r_c}\right)-y{\dot{f}}^2\right){\hat{\mathbf{o}}}_{\theta }+\ddot{z}{\hat{\mathbf{o}}}_h\end{array}\text{\hspace{1em}(8)}$$

同时，从星也满足自己的轨道运动方程。$${\ddot{\mathbf{r}}}_d=-\frac{\mu }{r_d^3}{\mathbf{r}}_d=-\frac{\mu }{r_d^3}\left(\begin{array}{c}{r}_c+x\\ y\\ z\end{array}\right)\text{\hspace{1em}(9)}$$

其中，${r_d=\sqrt{(}{r}_c+x)}^2+y^2+z^2$。

令(8),(9)各分量相等，得到了无摄动情况下准确的非线性运动方程
$$\begin{array}{rl}\ddot{x}-2\dot{f}\left(\dot{y}-y\frac{{\dot{r}}_c}{r_c}\right)-x{\dot{f}}^2-\frac{\mu }{r_c^2}& =-\frac{\mu }{r_d^3}\left(r_c+x\right)\\ \ddot{y}+2\dot{f}\left(\dot{x}-x\frac{{\dot{r}}_c}{r_c}\right)-y{\dot{f}}^2& =-\frac{\mu }{r_d^3}y\\ \ddot{z}& =-\frac{\mu }{r_d^3}z\end{array}\text{\hspace{1em}(10)}$$

简化

  式(10)是不考虑摄动情况下的精确解，适用于任意大的轨道，也适用于椭圆轨道。
  如果和主星轨道半径长度$r_c$相比，相对轨道坐标$(x,y,z)$很小，那么从星的轨道半径长度可以近似为(通过去掉轨道坐标$x,y,z$与$r_c$比值的二阶项)$$r_d=r_c\sqrt{1+2\frac{x}{r_c}}+\frac{x^2+y^2+z^2}{r_c^2}\approx r_c\sqrt{1+2\frac{x}{r_c}}\text{\hspace{1em}(11)}$$

因此有下列近似(一阶)
$$\frac{\mu }{r_d^3}\approx \frac{\mu }{r_c^3}\left(1-3\frac{x}{r_c}\right)\text{\hspace{1em}(12)}$$

将(12)代入(9),去掉高阶项，得

$$-\frac{\mu }{r_d^3}\left(\begin{array}{c}{r}_c+x\\ y\\ z\end{array}\right)\approx -\frac{\mu }{r_c^3}{\left(1-3\frac{x}{r_c}\right)}^O\left(\begin{array}{c}{r}_c+x\\ y\\ z\end{array}\right)\approx -\frac{\mu }{r_c^3}\left(\begin{array}{c}{r}_c-2x\\ y\\ z\end{array}\right)\text{\hspace{1em}(12)}$$
将(12)代入(10),并利用
$$\frac{\mu }{r_c^3}=\frac{r_c}{p}{\dot{f}}^2=\frac{{\dot{f}}^2}{1+e\cos f}\text{\hspace{1em}(13)}$$
经过整理后可得
$$\begin{array}{rl}\ddot{x}-x{\dot{f}}^2\left(1+2\frac{r_c}{p}\right)-2\dot{f}\left(\dot{y}-y\frac{{\dot{r}}_c}{r_c}\right)& =0\\ \ddot{y}+2\dot{f}\left(\dot{x}-x\frac{{\dot{r}}_c}{r_c}\right)-y{\dot{f}}^2\left(1-\frac{r_c}{p}\right)& =0\\ \ddot{z}+\frac{r_c}{p}{\dot{f}}^{2}z& =0\end{array}\text{\hspace{1em}(14)}$$
式(14)就是假定$x,y,z$相比$r_c$很小，并据此后得到的相对运动方程。
利用等式$\theta =\omega +f$，以及式(5),则式(14)可进一步变为式(15)这一更常用的形式$$\begin{array}{rl}\ddot{x}-x\left({\dot{\theta }}^2+2\frac{\mu }{r_c^3}\right)-y\ddot{\theta }-2\dot{y}\dot{\theta }& =0\\ \ddot{y}+x\ddot{\theta }+2\dot{x}\dot{\theta }-y\left({\dot{\theta }}^2-\frac{\mu }{r_c^3}\right)& =0\\ \ddot{z}+\frac{\mu }{r_c^3}z& =0\end{array}\text{\hspace{1em}(15)}$$

进一步简化

  如果主星是圆轨道，则$e=0,p=r_c$。并且$r_c$是常数。同时对于圆轨道，平均轨道速率$n$等于真近点角速率$\dot{f}$。将以上条件代入式(14)，得
$$\begin{array}{rl}\ddot{x}-2n\dot{y}-3n^{2}x& =0\\ \ddot{y}+2n\dot{x}& =0\\ \ddot{z}+n^{2}z& =0\end{array}\text{\hspace{1em}(16)}$$

这就是CW方程。

无量纲化

  将函数无量纲化
$$u=\frac{x}{r_c}v=\frac{y}{r_c}w=\frac{z}{r_c}$$

自变量取为无量纲的真近点角$f$，将对时间求导改为对$f$求导，
$$({)}^{\prime }\equiv \frac{\mathrm{d}()}{\mathrm{d}f}$$

可得$$\begin{array}{rl}& \frac{\dot{x}}{r_c}=u^{\prime }\dot{f}+u\frac{{\dot{r}}_c}{r_c}\frac{\ddot{x}}{r_c}=u^{\prime \prime }{\dot{f}}^2+u{\dot{f}}^2\left(1-\frac{r_c}{p}\right)\\ & \frac{\dot{y}}{r_c}=v^{\prime }\dot{f}+v\frac{{\dot{r}}_c}{r_c}\frac{\ddot{y}}{r_c}=v^{\prime \prime }{\dot{f}}^2+v{\dot{f}}^2\left(1-\frac{r_c}{p}\right)\\ & \frac{\dot{z}}{r_c}=w^{\prime }\dot{f}+w\frac{{\dot{r}}_c}{r_c}\frac{\ddot{z}}{r_c}=w^{\prime \prime }{\dot{f}}^2+w{\dot{f}}^2\left(1-\frac{r_c}{p}\right)\end{array}$$

代入(17),得
$$\begin{array}{rl}{u}^{\prime \prime }-2v^{\prime }-\frac{3u}{1+e\cos f}& =0\\ v^{\prime \prime }+2u^{\prime }& =0\\ w^{\prime \prime }+w& =0\end{array}$$

该该式在形式上与CW很类似，但也仅仅是类似而已。

2. Hill坐标系下的相对轨道

  CW方程(16)的解为
$$\begin{array}{l}x(t)=A_0\cos (nt+\alpha )+x_{\mathrm{o}\mathrm{f}\mathrm{f}}\\ y(t)=-2A_0\sin (nt+\alpha )-\frac{3}{2}nt{x}_{\mathrm{o}\mathrm{f}\mathrm{f}}+y_{\mathrm{o}\mathrm{f}\mathrm{f}}\\ z(t)=B_0\cos (nt+\beta )\end{array}\text{\hspace{1em}(17)}$$

其中，$d={\dot{y}}_0+2n{x}_0$，$x_{\mathrm{o}\mathrm{f}\mathrm{f}}=\frac{2d}{n}$。
  要使轨道闭合，则式(17)第二项关于时间$t$的一次项系数必须为0，即$x_{\mathrm{o}\mathrm{f}\mathrm{f}}=0$，即$d=0$。
因此
$$d={\dot{y}}_0+2n{x}_0=0\text{\hspace{1em}(18)}$$

此时对应的闭合解为
$$\begin{array}{l}x(t)=A_0\cos (nt+\alpha )\\ y(t)=-2A_0\sin (nt+\alpha )+y_{\mathrm{o}\mathrm{f}\mathrm{f}}\\ z(t)=B_0\cos (nt+\beta )\end{array}\text{\hspace{1em}(17)}$$

  在式(17)中，线性化运动的解析解由六个不变量表示，即轨道面内以及垂直轨道面的运动幅度参数$A_0,B_0$,相位角$\alpha ,\beta$以及径向的偏置$x_{\mathrm{o}\mathrm{f}\mathrm{f}}$和轨迹方向的初始偏置$y_{\mathrm{o}\mathrm{f}\mathrm{f}}$。
$$\begin{array}{llllll}[A_1& \alpha & x_{\text{off }}& y_{\text{off }}& B_1& \beta {]}^T\end{array}\text{\hspace{1em}(18)}$$

  这六个无摄动情况下的不变量被称为相对轨道根数。由于本文对应的这些是线性化模型对应的不变量，因此也被称为线性化相对轨道根数。其优势在于能比较直观的看到相对运动的几何形状，因为其常数项，三角函数项以及随时间增长的长期项都很直观。
  这组相对轨道根数可由任意时刻的Hill系下的位置和速度坐标求得。
$$\begin{array}{l}\begin{array}{rl}{A}_0& =\frac{\sqrt{9n^2{x}^2+x^2+12nx\dot{y}+4{\dot{y}}^2}}{n}\\ \alpha & =\arctan \left(\frac{-\dot{x}}{-3nx-2y}\right)-nt\end{array}\\ x_{\mathrm{o}\mathrm{f}\mathrm{f}}=4x+2\frac{j}{n}\\ y_{\mathrm{o}\mathrm{f}\mathrm{f}}=-2\frac{\dot{x}}{n}+y+(6nx+3\dot{y})t\\ B_0=\frac{\sqrt{n^2{z}^2+{\dot{z}}^2}}{n}\\ \beta =\arctan \left(\frac{-\dot{z}}{nz}\right)-nt\end{array}\text{\hspace{1em}(19)}$$

  注意到如果式(17)中$A_0$或$B_0$为零时，分别导致式(19)中的$\alpha ,\beta$求解公式奇异，说明$\alpha ,\beta$不能确定。

另一种解的形式

  CW方程的另一种解的形式用初始的六个状态$\left(x_0,y_0,z_0,{\dot{x}}_0,{\dot{y}}_0,{\dot{z}}_0\right)$来表示。
$$\begin{array}{rl}x(t)=& (4-3\cos (nt))x_0+\frac{\sin (nt)}{n}{\dot{x}}_0+\frac{2}{n}(1-\cos (nt)){\dot{y}}_0\\ y(t)=& 6(\sin (nt)-nt)x_0+y_0-\frac{2}{n}(1-\cos (nt)){\dot{x}}_0\\ & +\left(\frac{4}{n}\sin (nt)-3t\right){\dot{y}}_0\\ z(t)=& \cos (nt)z_0+\frac{\sin (nt)}{n}{\dot{z}}_0\\ \dot{x}(t)=& 3n\sin (nt)x_0+\cos (nt){\dot{x}}_0+2\sin (nt){\dot{y}}_0\\ \dot{y}(t)=& -6n(1-\cos (nt))x_0-2\sin (nt){\dot{x}}_0+(4\cos (nt)-3){\dot{y}}_0\\ \dot{z}(t)=& -n\sin (nt)z_0+\cos (nt){\dot{z}}_0\end{array}\text{\hspace{1em}(19)}$$

  $\left(x_0,y_0,z_0,{\dot{x}}_0,{\dot{y}}_0,{\dot{z}}_0\right)$也是一组不变量。
  这一解的形式建立起了与初始位置和速度的直接关系。