matlab仿真|基于 CW 方程与解析解法的卫星状态解算方法

一、基于 CW 方程的卫星动力学原理

为描述己方卫星与敌方卫星之间的运动和相对运动，引入三个坐标系。在空间动力学研究中，惯性坐标系$N$，也被称作地心赤道惯性坐标系，是一个以地球质心为原点的参考框架，其中$X$轴指向春分点，$Z$轴垂直于赤道平面指向北极方向，而$Y$轴则按右手法则垂直于$X$轴和$Z$轴确定。参考卫星轨道坐标系（记为$C$）的原点同样位于参考卫星的质心，其$x$轴指向从地心到参考卫星的位置矢量，$y$轴在参考卫星轨道平面内并与$x$轴垂直，指向参考卫星运动的前方，$z$轴则垂直于轨道平面，与$x$轴和$y$轴构成右手坐标系。

基本假设包括：地球被视为均质的中心引力场球体，地球和卫星均被简化为质点，且敌方卫星不进行主动控制。这些假设有助于简化问题并使得动力学分析更为可行。

因为惯性坐标系下的位置速度信息都较大（1e7 米起）不利于项目展示或算法计算，因此在参考卫星坐标系中建立相对运动模型（1e5 米起），以方便展示与算法计算。简化后的 C-W 方程如下（简化条件需要自行阅读论文）：

$$\{\begin{array}{l}\ddot{x}-2n\dot{y}-3n^{2}x=a_x\\ \ddot{y}+2n\dot{x}=a_y\\ \ddot{z}+n^{2}z=a_z\end{array}$$

若我们将此参考卫星作为原点，此 C-W 方程则表示的是卫星在参考卫星坐标系下相对于参考卫星的位置，但因为参考卫星为原点因此此公式计算的位置、速度和加速度均为参考卫星坐标系下的绝对位置和状态。其中$n$为参考卫星的平均角速度计算方法如下：

$$n=\sqrt{\frac{\mu }{r^3}}$$

其中$\mu$为地球引力常数，$r$为参考卫星轨道半径，$a_x、a_y、a_z$为对卫星产生的控制加速度。

根据状态转移方程求解能够得到以下公式和矩阵（如有不对请参考代码）：

$$X_i\text{（}t\text{）}=\varphi \left(\tau \right)X_i\left(t_0\right)+\Psi \left(\tau \right)a_i$$

$$\tau =t-t_0$$

$$\varphi \left(\tau \right)=\left[\begin{array}{cccccc}4-3\cos \left(n\tau \right)& 0& 0& \frac{\sin \left(n\tau \right)}{n}& \frac{2}{n}\left(1-\cos \left(n\tau \right)\right)& 0\\ 6\left(\sin \left(n\tau \right)-n\tau \right)& 1& 0& \frac{2}{n}\left(\cos \left(nt\right)-1\right)& \frac{4\sin \left(n\tau \right)}{n}-3\tau & 0\\ 0& 0& \cos \left(n\tau \right)& 0& 0& \frac{\sin \left(n\tau \right)}{n}\\ 3n\sin \left(n\tau \right)& 0& 0& \cos \left(n\tau \right)& 2\sin \left(n\tau \right)& 0\\ 6n\left(\cos \left(n\tau \right)-1\right)& 0& 0& -2\sin \left(n\tau \right)& 4\cos \left(n\tau \right)-3& 0\\ 0& 0& -n\sin \left(n\tau \right)& 0& 0& \cos \left(n\tau \right)\end{array}\right]$$

$$\Psi \left(\tau \right)=\left[\begin{array}{ccc}\frac{1-\cos \left(n\tau \right)}{n^2}& \frac{2\left(n\tau -\sin \left(n\tau \right)\right)}{n^2}& 0\\ \frac{2\left(\sin \left(n\tau \right)-n\tau \right)}{n^2}& \frac{4\left(1-\cos \left(nt\right)\right)}{n^2}-\frac{3{\tau }^2}{2}& 0\\ 0& 0& \frac{1-\cos \left(n\tau \right)}{n^2}\\ \frac{\sin \left(n\tau \right)}{n}& \frac{2\left(1-\cos \left(n\tau \right)\right)}{n}& 0\\ \frac{2\left(\cos \left(nt\right)-1\right)}{n}& \frac{4\sin \left(n\tau \right)}{n}-3\tau & 0\\ 0& 0& \frac{\sin \left(n\tau \right)}{n}\end{array}\right]$$

结合以上的公式与矩阵，当我们确定$\tau$和这一段时间恒定不变的控制加速度就可以解算出$t$时刻的卫星的状态$X_i（t）$。注意此解算方法适合于此段时间内控制加速度不变的情况和时间较长的状态解算 100s 这样。

二、代码复现

本文以 matlab 函数仿真为例。

function [posOutput,velOutput] = posvelcalculateFunction(posInput,velInput,dtInput,accelerationInput)
%基于CW方程与解析解法的卫星位置解算函数
%输入：1.posInput：卫星0时刻位置信息，三行一列；
%      2.velInput：卫星0时刻速度信息，三行一列；
%      3.dtInput：恒定加速度应用时间，三行一列；
%      4.accelerationInput：施加在卫星上的恒定加速度，三行一列；
%输出：1.posOutput，velOutput在应用dtInput时长的accelerationInput后的卫星位置与速度信息，各为三行一列
    mu = 3.986005e14;%地球引力常数m^3/s^2
    R_spacestation0 = 42171e3;%参考卫星的轨道半径 单位m
    n = sqrt(mu/R_spacestation0^3);%参考卫星轨道平均角速度 单位rad/s
    n2=mu/R_spacestation0^3;
    nt=n*dtInput;
    snt=sin(nt);
    cnt=cos(nt);
    A=[ 4-3*cnt    ,0  ,0     ,snt/n        ,(2*(1-cnt))/n        ,0    ;...
        6*(snt-nt) ,1  ,0     ,(2*(cnt-1))/n,(4*snt/n)-(3*dtInput),0    ;...
        0          ,0  ,cnt   ,0            ,0                    ,snt/n;...
        3*n*snt    ,0  ,0     ,cnt          ,2*snt                ,0    ;...
        6*n*(cnt-1),0  ,0     ,-2*snt       ,4*cnt-3              ,0    ;...
        0          ,0  ,-n*snt,0            ,0                    ,cnt  ];
    B=[ (1-cnt)/n2   ,(2*(nt-snt))/n2             ,0         ;...
        2*(snt-nt)/n2,(4*(1-cnt)/n2-3*(dtInput^2)/2),0         ;...
        0            ,0                           ,(1-cnt)/n2;...
        snt/n        ,2*(1-cnt)/n                 ,0         ;...
        2*(cnt-1)/n  ,4*snt/n-3*dtInput           ,0         ;...
        0            ,0                           ,snt/n     ];
    posvelInput=[posInput;velInput];
    posvelOutput=A*posvelInput+B*accelerationInput;
    posOutput=posvelOutput(1:3);
    velOutput=posvelOutput(4:6);
end

当参数为$X,Y,Z=(0,40000,0);VX,VY,VZ=(n\ast 20000,0,0)$。控制加速度为 0，其中$n=\sqrt{\frac{\mu }{r^3}}$为参考卫星的平均角速度，我这边$r=42171e3$。能出现以下轨迹是没问题的。