6.4定轨

  其本质是利用非线性最小二乘法来进行迭代求解。利用多个时间点的测量来求取轨道初值。因此重点就是要求各个时间点的测量对于轨道初值的导数。根据链式法则，这一导数又被分成两部分，即各时间点的测量值对于测量时间点的轨道的雅各比矩阵，以及测量时间点的轨道对于初始轨道的雅各比矩阵。前者是测量雅各比，后者即状态转移矩阵。即
$$H_k=\frac{\partial {\mathbf{h}}_k}{\partial {\mathbf{x}}_k}\Phi \left(t_k,t_0\right)\text{\hspace{1em}(1)}$$

对于轨道问题， 状态转移矩阵有解析表达式，但较为复杂。这里采用近似处理办法：

1. 当前估计值$x_0$利用轨道动力学外推得到$x_1$，$x_2$，$x_3$…$x_m$。它们都是$x_0$的函数。
2. $t_1$时刻，$F_1\equiv \frac{\partial \mathbf{f}}{\partial \mathbf{x}}{|}_{x_1}$,状态转移矩阵近似值为$\Phi \left(t_1,t_0\right)\approx exp(F_1(t_1-t_0))$。这种处理等同于将$F_1$在$t_0$到$t_1$的时间段内看作常数。
3. $t_2$时刻，$F_2\equiv \frac{\partial \mathbf{f}}{\partial \mathbf{x}}{|}_{x_2}$,状态转移矩阵近似值为$\Phi \left(t_2,t_1\right)\approx exp(F_2(t_2-t_1))$。这种处理等同于将$F_2$在$t_1$到$t_2$的时间段内看作常数。则根据状态转移矩阵性质，有$\Phi \left(t_2,t_0\right)=\Phi \left(t_2,t_1\right)\Phi \left(t_1,t_0\right)$。以此类推，获得所有时刻的状态转移矩阵$\Phi \left(t_m,t_0\right)$
4. 利用$x_1$，$x_2$，$x_3$…$x_m$计算各测量时刻的残差$e_1$,$e_2$…$e_m$以及利用公式(1)求取$H_1$,$H_2$…$H_m$。
5. 进行迭代更新，获得新的估计值$x_0$。重复步骤(1),直到某次更新前后的误差小于某个设定值。

流程图如下：

心得，评注

1.length(x)获得矩阵x的行数与列数中较大的那个。
2.y=zeros(m*n,1); y(:)=x;是把m*n维矩阵按列叠起来形成的列向量赋给y。y记得要初始化。