基于自适应测量噪声方差估计的星载GPS定轨方法与流程

本发明涉及一种星载gps定轨方法，尤其涉及一种基于自适应测量噪声方差估计的星载gps定轨方法。

背景技术：

随着全球定位系统(globalpositioningsystem,gps)的发展，为航天器导航定轨提供了一种全新的方式。与传统的地基定轨相比，星载gps能够提供更加精确、连续的导航、定轨服务，因此星载gps技术已逐渐成为低轨卫星定轨的一条有效途径。

目前，根据是否采用低轨卫星动力学信息，星载gps定轨主要分为以下3种：1)几何法，该方法直接利用星载gps接收机的伪距信息进行定位计算，其最大特点是不受动力学模型误差的影响，但对观测量误差及gps卫星几何分布敏感，只能确定观测时刻的卫星位置，不能对卫星轨道进行预报；2)动力法，该方法利用卫星动力学模型建立含参数的卫星运动方程，确定卫星理论轨迹，然后利用实测数据做轨道改进，能以较少的观测数据获得可靠的卫星轨道，并能对影响卫星运动的各种摄动力进行定量估计，对卫星长期运动趋势定量预报，但由于动力学模型的误差会随时间而增长，限制了其长弧段定轨的精度；3)约化动力法，该方法通过在低轨卫星的动力学模型上附加一个假想的力(过程噪声)，并在定轨过程中进行估计来平衡几何观测和动力学模型的贡献，解决了动力法对动力学模型误差非常敏感和几何法取决于相对几何与测量精度的问题。约化动力法定轨通常都可以采用卡尔曼滤波(kalmanfiltering,kf)进行数学描述，每次使用一个时刻的观测数据，序贯输出当前时刻的估计结果，十分适合于实时定轨，而且精度与事后处理相当。

kf算法只有在函数模型和随机模型精确已知的条件下才具有最小的均方误差，其中随机模型参数通常是根据经验选取的，需要良好的关于过程噪声方差q和测量噪声方差r的先验知识，这通常来自于大量的实验分析。然后，它们的值在整个估计过程中都被认为是保持不变的。这种不变性结构会导致kf不能获得最优的估计结果，有时甚至会导致滤波器发散，因为过程噪声和测量噪声往往取决于实际的应用环境和动态条件。为了能够在不同的环境下正常工作，q和r的值必须取得比较保守，以保持kf算法在最差情况下的稳定性，这会导致滤波结果变差。因此，需要寻找更加合理的确定随机模型参数的方法，这对于提高星载gps定轨精度具有十分重要的意义。

技术实现要素：

为了克服传统kf定轨方法在随机模型参数取值上的问题，本发明提出了一种基于自适应测量噪声方差估计(adaptivemeasurementnoisevarianceestimation,amnve)的星载gps定轨方法，该方法可以根据接收机载噪比(carriernoiseratio,c/n0)自适应调整r，从而使其更好地反映测量值的噪声统计特性，能够获得更高的定轨精度。

本发明解决技术问题所采用的技术方案如下：一种基于自适应测量噪声方差估计的星载gps定轨方法，该方法包括如下步骤：

(1)根据星载gps接收机跟踪环中的i、q两路信号，利用窄带宽带功率比值法估计信号载噪比c/n0；

(2)根据各通道跟踪环c/n0，计算l1、l2双频信号的伪距测量误差；

(3)组合l1、l2双频伪距测量值得到无电离层延迟的伪距测量值ρif及其误差σif；

(4)以ρif为观测量，利用ekf算法对卫星的轨道信息进行估计。

进一步的，所述利用窄带宽带功率比值法估计信号载噪比的具体过程如下：

(1.1)根据星载gps接收机跟踪环中的i、q支路信号，计算宽带功率和窄带功率，其中第k时刻的宽带功率wbpk为：

第k时刻的窄带功率nbpk为：

式中，m表示计算宽带和窄带功率的累加样点数，通常取1≤m≤20；ii表示第i时刻i(同相)支路积分清零输出；qi表示第i时刻q(正交)支路积分清零输出；

(1.2)进一步计算第k时刻的窄带和宽带功率比值：

(1.3)为了降低噪声误差，对n个时刻的np计算平均值μnp：

(1.4)求得gps接收机的信号载噪比c/n0：

其中，t表示相干积分时间。

进一步的，所述步骤(2)具体如下：

根据步骤(1)分别求得l1信号跟踪环的载噪比(c/n0)l1和l2信号跟踪环的载噪比(c/n0)l2，通过各自的信号载噪比，求得gps接收机l1信号跟踪环的测量噪声误差σl1和l2信号跟踪环的测量噪声误差σl2：

其中，tc表示码片长度，bl1、bl2分别表示l1、l2信号跟踪环的噪声带宽。

进一步的，所述步骤(3)具体如下：

由于伪距测量值受到电离层延迟的影响，需要对gpsl1、l2双频测量值ρl1、ρl2进行线性组合，从而得到校正后的伪距测量值ρif为：

ρif的测量噪声误差σif为：

其中，α＝(fl1/fl2)²，fl1、fl2分别表示gpsl1、l2信号所对应的载波频率。

进一步的，所述步骤(4)具体过程如下：

设状态向量xr表示位置向量，xv表示速度向量，表示接收机钟差；

(4.1)初始化：利用最小二乘法对gps伪距定位线性方程组的求解结果作为位置向量初始值和钟差初始值利用接收机多普勒频移测量值求得的瞬时速度作为速度向量的初始值

(4.2)预测：在前一时刻状态估计值的基础上预测当前时刻的状态值，具体过程为：

(4.2.1)状态预测：

其中，表示第k时刻位置和速度状态预测值，表示第k-1时刻位置和速度状态估计值，tk表示第k时刻的时间值，xrv表示卫星的位置和速度向量，表示第k时刻接收机钟差预测值，表示第k-1时刻接收机钟差估计值；

(4.2.2)状态预测方差矩阵计算：

其中，表示第k时刻的状态预测方差矩阵，φk表示第k时刻的状态转移矩阵，表示第k-1时刻的状态估计方差矩阵；q表示状态噪声方差矩阵；

(4.3)更新：利用当前第k时刻测得的m个伪距测量值来更新状态预测值，具体包括：

(4.3.1)滤波增益矩阵计算：

其中，kk表示第k时刻的滤波增益矩阵；hk表示第k时刻的观测矩阵；rk表示第k时刻的测量噪声方差矩阵；

根据伪距测量值噪声误差σif构建测量噪声方差矩阵rk：

其中，σif的下标k表示第k时刻，上标m表示第m个观测量；

(4.3.2)状态估计：

其中，表示第k时刻的卫星轨道估计值，表示第k时刻的卫星轨道预测值，h(.)表示非线性观测向量函数；

(4.3.3)状态估计方差矩阵计算：

其中，表示第k时刻的状态估计方差矩阵，i表示单位矩阵。

重复步骤(4.2)-步骤(4.3)可以估计得到任意第k时刻的卫星轨道信息

本发明的有益效果是：

(1)大大降低了对系统测量噪声先验信息的要求，而且在滤波过程中可以根据信号载噪比的变化情况自动调整测量噪声方差矩阵r，定轨精度高；

(2)采用伪距测量值作为观测量，直接在地心地固坐标系下进行卫星轨道确定，无需进行空间坐标系之间的转换，从而简化了计算量，而且适用于可见星少于4颗的情况；

(3)计算复杂度低，无需其它额外传感器辅助，结构简单，适用于计算资源受限的在轨应用。

附图说明

图1是本发明amnve信号处理结构框图；

图2是grace卫星与gps卫星prn07视线方向的距离变化动态图；

图3是消除电离层延迟后真实伪距误差与通过c/n0估计伪距误差曲线图；

图4是ekf算法与本发明amnve算法径向误差曲线图；

图5是ekf算法与本发明amnve算法切向误差曲线图；

图6是ekf算法与本发明amnve算法法向误差曲线图；

图7是ekf算法与本发明amnve算法3d定位误差曲线图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

如图1所示，本发明提供一种基于自适应测量噪声方差估计的星载gps定轨方法，具体包括如下步骤：

1)根据星载gps接收机跟踪环中的i、q两路信号，利用窄带宽带功率比值法估计信号载噪比c/n0，其基本原理是以两种不同噪声带宽下的信号功率的比值来计算c/n0，从而避开直接求信号功率和噪声功率的运算，具体过程如下：

1.1)根据星载gps接收机跟踪环中的i、q支路信号计算宽带功率和窄带功率。其中第k时刻的宽带(带宽为1/t)功率wbpk为：

第k时刻的窄带(带宽为1/mt)功率为：

式中，t表示相干积分时间，对于gpsl1c/a码信号，通常取t＝1ms；ii、qi分别表示第i时刻i、q支路的积分清零输出；m表示计算宽带和窄带功率的累加样点数，为了避免在累加过程中出现电文比特翻转，通常取1≤m≤20。

1.2)进一步计算第k时刻的窄带和宽带功率比值：

1.3)为了降低噪声误差，对n个时刻的np计算平均值μnp：

1.4)根据上式计算结果求得gps接收机该通道信号的c/n0(c/n0估计值的更新时间为kmt)：

2)根据各通道跟踪环c/n0，计算l1、l2双频信号的伪距测量误差：

gps接收机码环测量误差主要包括热噪声误差和动态应力误差，其中热噪声误差表示为：

n阶码环的动态应力误差表示为：

其中，tc表示码片长度，bl表示跟踪环路噪声带宽，ω0表示环路特征频率，r表示卫星与接收机之间的视线距离。

由图2可知，grace卫星与gps卫星prn07之间视线方向加加速度很小，即d³r/dt³≈0，这意味着对于三阶码环，其动态应力误差re≈0；而且，通过载波辅助码技术，可以去除码跟踪环中的所有动态。因此，码环测量误差可以忽略动态应力误差，只考虑热噪声误差。

由此可得，gpsl1、l2信号伪距测量误差σl1、σl2可以直接通过各自跟踪环的载噪比估计值(c/n0)l1、(c/n0)l2计算得到：

其中，tc表示码片长度，bl1、bl2分别表示l1、l2信号跟踪环的噪声带宽。

3)组合l1、l2双频伪距测量值得到无电离层延迟的伪距测量值ρif及其误差σif：

由于伪距测量值受到电离层延迟的影响，需要对gpsl1、l2双频测量值ρl1、ρl2进行线性组合，从而得到无电离层延迟的伪距测量值ρif为：

ρif的测量噪声误差σif为：

其中，α＝(fl1/fl2)²，fl1、fl2分别表示gpsl1、l2信号所对应的载波频率。

4)以ρif为观测量，利用ekf算法对卫星的轨道信息进行估计：

设状态向量xr表示位置向量，xv表示速度向量，表示接收机钟差，具体步骤如下：

4.1)初始化：利用最小二乘法对gps伪距定位线性方程组的求解结果作为位置向量初始值和钟差初始值利用接收机多普勒频移测量值求得的瞬时速度作为速度向量的初始值

4.2)预测：在前一时刻状态估计值的基础上预测当前时刻的状态值，具体过程为：

4.2.1)状态预测：

其中，表示第k-1时刻位置和速度状态估计值，表示第k时刻位置和速度状态预测值，tk表示第k时刻的时间值，积分运算采用四阶龙格-库塔方法，xrv表示卫星的位置和速度向量，表示第k-1时刻接收机钟差估计值，表示第k时刻接收机钟差预测值。

4.2.2)状态预测方差矩阵计算：

其中，表示第k-1时刻的状态估计方差矩阵；表示第k时刻的状态预测方差矩阵；φk表示第k时刻的状态转移矩阵；q表示状态噪声方差矩阵。

4.3)更新：利用当前第k时刻测得的m个伪距测量值来更新状态预测值，具体包括：

4.3.1)滤波增益矩阵计算：

其中，kk表示第k时刻的滤波增益矩阵；hk表示第k时刻的观测矩阵；rk表示第k时刻的测量噪声方差矩阵。由图3可知，直接通过c/n0估计的伪距噪声误差σif能够较为合理地反映真实伪距观测量的噪声统计特性，因此可以直接用来构建测量噪声方差矩阵：

其中，σif的下标k表示第k时刻，上标m表示第m个观测量；

4.3.2)状态估计：

其中，表示第k时刻的卫星轨道估计值；表示第k时刻的卫星轨道预测值；h(·)表示非线性观测向量函数。

4.3.3)状态估计方差矩阵计算：

其中，表示第k时刻的状态估计方差矩阵；i表示单位矩阵。

重复步骤(4.2)-步骤(4.3)可以估计得到任意第k时刻的卫星轨道信息在上述估计过程中，测量噪声方差矩阵r会根据gps接收机各通道信号c/n0的变化而自动调整，克服了传统固定设计ekf算法在先验信息不充分和动态环境中存在的不足，能够获得更高的定轨精度。

实施例：

以2008年1月5日-6日gracea卫星的level1b数据为例，该数据中包含了星载gps接收机采集的l1、l2双频伪距测量值和信噪比等数据。采用本发明提出的amnve方法和传统的ekf算法对gracea卫星的轨道进行实时估计，并将结果与事后精密星历进行对比得到其定轨精度。

由图4、5、6可知，传统的ekf算法在径向、切向和法向误差为0.93m、1.27m和0.83m，而本发明提出的amnve方法在三个方向的误差为0.78m、0.97m和0.68m，分别降低了16％、23％和18％；由图7可知，3d定位误差由1.78m下降至1.42m，降低了20％。在传统的ekf算法中，测量噪声方差r在整个滤波过程中保持不变，而实际测量值由于各种因素的影响，其噪声统计特性是一直在变化的，这会导致传统ekf算法性能降低；而本发明提出的amnve方法通过实时监测星载gps接收机各通道信号c/n0的变化来自适应调整r，使其正确反映观测量的噪声统计特性，因而可以获得更高的定轨精度。

同时，为了分析动力学模型对ekf算法与本发明amnve方法性能的影响，对4种不同阶数重力场模型下的定轨精度进行了比较，结果如表1所示。由表1可知，在20×20阶重力场模型情况下，本发明提出的amnve方法与传统的ekf算法性能接近，定轨精度只提升了0.09m，这是由于重力场模型比较精确，导致在估计过程中对测量值的依赖程度会降低，使得本发明的优势不明显，但重力场模型越精确，计算量也越大；随着重力模型精度降低，计算量下降，本发明amnve方法的优势越来越明显，比如在只考虑j2时，amnve方法将定轨精度提高了0.36m，这非常适合于在星载计算能力受限的环境中使用。

表1：