GPS广播星历参数拟合的雅可比矩阵数值导数计算方法①

2012-07-23邢志成王解先

全球定位系统 2012年1期

邢志成，王解先

（同济大学土木工程学院测量与国土信息工程系，上海200092）

0 引言

用户利用GPS广播星历参数可以计算出GPS卫星的位置，继而解算用户的位置；反过来，也可以利用GPS卫星的位置，拟合出广播星历的参数［1－2]。GPS广播星历就是由地面控制中心在确定卫星轨道的基础上，以卫星的位置作为观测量拟合出广播星历参数再上传给卫星，然后，卫星再转发这些星历参数给用户［3]。

利用最小二乘原理拟合广播星历参数的过程中，偏导数矩阵的求解有解析法和数值导数法两种方式。前者需要推导位置矢量对待估参数的偏导数，形式复杂，若改变参数形式，偏导数的形式需要重新推导；与前者相比，数值导数计算方法形式简单，计算方便，便于计算机操作［4]。采用数值导数计算方法，利用IGS精密星历卫星坐标数据，根据最小二乘平差原理拟合出GPS广播星历的15个参数；并分析了计算参数的雅可比矩阵时，增量δ的选取对拟合精度和效率的影响，得出了一些结论。

1 GPS广播星历参数拟合算法

GPS广播星历参数以开普勒轨道根数为基础，加上表示轨道摄动调和系数以及随时间的变化率，一共有15个：6个开普勒轨道根数ω，Ω0，M）和9个摄动力参数（Δn，i，Ω，Cuc，Cus，Crc，Crs，Cic，Cis）.15个参数的具体意义参考文献［2].

1.1 由广播星历计算卫星位置

若要计算时刻t卫星的空间坐标，根据广播星历参数，按如下步骤计算：

1）求长半轴A

2）计算平角速度n0

3）计算从需要时刻到参考时刻的时间差tk

4）改正平角速度n

5）计算平近点角Mk

6）按下式迭代计算偏近点角Ek

7）由下式计算真近点角

8）计算纬度参数φk

9）周期改正项

10）计算改正后的纬度参数uk

11）计算改正后的向径rk

12）计算改正后的倾角ik

13）计算卫星在轨道平面内的坐标（xk，yk）

14）改正升交点的经度

15）最后计算卫星在地固坐标系中的坐标Rk＝（XkYkZk）

1.2 数值导数的概念及雅可比矩阵的确定

若y＝f（x1，x2，…，xi…），y到xi的导数的定义为［5]

当δ足够小时，这一导数可以近似写成数值导数的形式

结合1.1节，则一组卫星坐标Rk＝（XkYkZk）到Xi（i＝1，2，…15）的雅可比矩阵可表示为

若利用同一颗卫星的上下各5个历元共11组坐标，则可以得到式（4）的参数的雅可比矩阵B.

1.3 平差模型

对于给定的参考时刻t，设参数为Xi（i＝1，2，…15），由一个历元的卫星坐标Rk＝［XkYkZk]T，就可以根据式（1）建立如下的误差方程

式（6）中的R0代表R的近似值，X0代表X的近似值，F（X0）＝R0.

则（6）式可以写为

式（7）中L＝（R－R0）（根据最小二乘原理，利用间接平差求解，可得

通过式（5），对每个历元的卫星坐标可列出3个误差方程，若选取11个历元，那么一共有33个方程，对于15个参数，进行间接平差迭代求解。迭代的收敛条件为：｜σi＋1－σi｜／σi＜ξ，其中，ξ是预先给定的阈值，取1.0e－05.是第i次迭代的单位权方差［2]。

2 算例分析

使用2010年1月1日从0时45分到3时15分第02颗GPS卫星的共11组IGS精密星历数据（间隔15min）。15参数中六个轨道根数的初值如表2所示，其余的9个参数初值设为0.

表2 六个轨道根数初值

为了比较小量δi的取值对拟合精度和效率的影响，算例根据1.2节所述的方法，利用公式Δ＝f（x1，x2，…，xi＋δ…）－f（x1，x2…，xi…）得到的差值Δ所在的量级不同，采用12种方案确定增量δi的取值，分别为：Δ量级为千万米级、百万米级、十万米级、万米级、千米级；百米级、十米级、米级、分米级、厘米级、毫米级、微米级。当然，增量的取值对计算结果的影响和计算机的性能有很大关系，本算例是在32位操作系统下完成的。根据这个方法可以确定增量δi的取值范围。表3列出了方案8米级的增量δi的取值。