CN106249256B - 基于粒子群优化算法的实时glonass相位偏差估计方法 - Google Patents

Abstract

本发明提出一种基于粒子群优化算法的实时GLONASS相位偏差估计方法，属于卫星定位系统和定位测量技术领域，利用IFB变化率粒子对所有单差模糊度参数进行改正，采用LAMBDA方法进行双差模糊度固定，将获得的RATIO值作为评价每个粒子的适应度，并使用粒子群优化算法对IFB变化率粒子的群体最优位置进行搜索，最终获得最优IFB变化率估值；解决了由于IFB和模糊度线性相关所带来的问题；实时在线对IFB变化率进行估计，避免由于观测条件变化导致已有的IFB变化率校正参数不可用现象，适合于实时动态定位；采用实时的GLONASS相位偏差估计方法实现实时模糊度固定，进一步拓展GLONASS的应用领域。

Description

基于粒子群优化算法的实时GLONASS相位偏差估计方法

技术领域

本发明属于卫星定位系统和定位测量技术领域，具体涉及一种基于粒子群优化算法的实时GLONASS相位偏差估计方法。

背景技术

2011年12月8日俄罗斯GLONASS系统控制中心宣布GLONASS已具有全运行能力。迄今为止，除了GPS以外，仅有GLONASS具有全球定位和授时能力，因此很多研究团队对GLONASS的精密定位方法进行研究。GLONASS采用频分多址技术，由于卫星所采用的频率不同，不同卫星在GLONASS接收机端的未校准相位延迟(uncalibratedphase delays，UPD)也不同，这个差异称为接收机端相位观测值的频率间偏差(inter-frequencybias，IFB)。尽管GLONASS卫星端的UPD可以通过站间差分消除，但对于具有不同频率的卫星之间组成的双差模糊度，由于受到接收机端的IFB影响，将不再具备整数特性，不能固定为整数。

相同品牌的GLONASS接收机间具有近似的IFB值，但是对于不同品牌的GLONASS接收机，IFB可能会不同(Wanninger L.Carrier-phase inter-frequency biases ofGLONASS receivers：Journal ofGeodesy，2012，86(2)：139-148)；为避免由于接收机间的IFB不同而造成模糊度无法固定的问题，需要使用相同品牌的接收机进行基线数据采集。但这个要求在实际野外数据采集过程中并不现实：由于GLONASS接收机生产厂家不断增加，使用不同品牌接收机共同作业的概率也在增加；此外，天线、天线连接线长度以及接收机的重启，也会对IFB产生影响(Chen J，Xiao P，Zhang Y，et al.GPS/GLONASS system biasestimation and application in GPS/GLONASS combined positioning：ChinaSatellite Navigation Conference (CSNC) 2013 Proceedings，2013：323-333)，这些情况均使得GLONASS模糊度固定变得更加复杂。

目前针对IFB特性的研究发现，IFB具有随频率线性变化的特征，相同类型接收机的IFB基本一致，两个频率(L1/L2)上IFB对应的几何距离偏差是一致的(Wanninger L，Wallstab-Freitag S.Combined processing of GPS，GLONASS，and SBAS code phase andcarrier phase measurements：Proceedings ofION GNSS，2007：866-875)。根据这些特征，提出多种GLONASS模糊度固定的方法。其主要思想是预先求出IFB参数，处理时作为改正数。然而，由于模糊度参数与IFB的线性相关性，这些方法在估计IFB时一般要求有相当长的观测数据和已知的测站坐标，在确定IFB时还要用同步的GPS观测数据来增加模糊度参数的稳定性，最终得到相应的IFB信息。也正是因为这个原因，这些方法的效率很低，精度也有待提高。同时，这样的方法也很难满足实时定位中对IFB的快速野外标定的需要，因为IFB可能会受到环境的变化以及接收机重新启动的影响而发生变化。为了实现IFB的快速估计，Tian等学者采用粒子滤波对IFB进行估计，进而实现GLONASS模糊度快速固定(Tian Y，Ge M，Neitzel F.Particle filter-based estimation of inter-frequency phase bias forreal-time GLONASS integer ambiguity resolution：Journal ofGeodesy，2015，89(11)：1145-1158)。这种方法的优点为(1)在没有增加待估参数数量情况下，采用搜索的方式实现GLONASS模糊度固定；(2)在没有IFB先验信息的情况下，采用较少历元的GLONASS数据便可实现模糊度固定。但该方法也存在一些不可忽视的缺点：已有研究结果表明，IFB变化率的取值范围非常小，通常只有给定的IFB变化率值距离真值在±4mm/FN范围内才能使GLONASS模糊度固定。因此在使用粒子滤波进行参数估计过程中，粒子的采样间隔一般设置为1mm/FN。由于最大的IFB变化率值小于0.10m/FN，因此粒子的取值区间设为[-0.10，0.10]m/FN。根据粒子采样间隔以及取值区间，可以得出在粒子集合中粒子的总数为200个。这说明在进行参数估计过程中，每一历元需要进行200次的法方程构建，解算以及模糊度固定，而其中只有少数的几次解算才是有效的，这表明该算法的效率并不高，增加了计算负担，并且当GLONASS卫星个数较多或GLONASS与其它GNSS系统进行组合定位时，进行上述解算需要更多的时间，很难满足实时模糊度固定的要求；此外，由于粒子滤波采用加权平均方法进行参数计算，可能会使解算出来的结果十分接近IFB变化率真值，但并不在±4mm/FN范围内，如果使用这样的估值对IFB变化率参数进行改正，这将导致模糊度并不能成功固定。

粒子群优化算法(PSO)是一种进化计算技术，该算法模拟鸟群觅食的行为，通过鸟之间的集体协作使得群体达到最优。与遗传算法(GA)类似，系统初始化为一组随机解，通过迭代搜寻最优值。它也是一种基于群体迭代搜寻最优解的优化工具，但它无需进行交叉和变异操作，而是通过粒子在解空间追随最优的粒子进行搜索。PSO以其结构简单、容易实现、精度高、收敛快等优点引起了学术界的重视，目前已广泛应用于函数优化、神经网络训练、模糊系统控制以及其他遗传算法的应用领域。

发明内容

针对现有技术的不足，本发明提出一种基于粒子群优化算法的实时GLONASS相位偏差估计方法，以达到解决由于IFB和模糊度线性相关所带来的问题，提高搜索IFB变化率的准确性，避免由于观测条件变化导致已有的IFB变化率校正参数不可用现象的目的。

一种基于粒子群优化算法的实时GLONASS相位偏差估计方法，包括以下步骤：

步骤1：获取基准站和移动站的GLONASS观测数据，包括GLONASS伪距和载波观测数据，并逐历元进行读取；

步骤2、获取基准站和移动站概略坐标，选择基准站和移动站的共视卫星，对位置坐标、对流层延迟、电离层延迟以及模糊度参数进行初始化或更新，对基准站和移动站双差观测方程进行线性化处理，采用kalman滤波获得单差模糊度参数向量及其方差协方差阵；

步骤3、在将以周为单位的单差模糊度参数投影为以周为单位的双差模糊度参数的过程中，利用IFB变化率粒子对所有单差模糊度参数进行改正，采用LAMBDA方法进行双差模糊度固定，将获得的RATIO值作为评价每个粒子的适应度，并使用粒子群优化算法对IFB变化率粒子的群体最优位置进行搜索，最终获得最优IFB变化率估值；

步骤4、将确定出来的接收机端相位观测值的IFB变化率估值作为改正数对单差模糊度参数进行改正，采用LAMBDA方法进行双差模糊度固定，即完成对GLONASS模糊度的实时固定。

步骤2所述的获取基准站和移动站概略坐标，选择基准站和移动站的共视卫星，对位置坐标、对流层延迟、电离层延迟以及模糊度参数进行初始化或更新，对基准站和移动站双差观测方程进行线性化处理，采用kalman滤波获得单差模糊度参数向量及其方差协方差阵；

具体步骤如下：

步骤2-1、采用伪距单点定位方法获得基准站和移动站的概率坐标；

步骤2-2、根据基准站和移动站的观测数据选择共视卫星；

步骤2-3、在忽略IFB影响的基础上，构建基准站和移动站的非差观测方程；

步骤2-4、判断当前是否为首历元，若是，则对包含位置、对流层延迟、电离层延迟以及模糊度参数的状态向量及其方差协方差阵进行初始化，否则根据上一历元的状态向量对当前历元的状态向量进行预测；

步骤2-5、在步骤2-3的基础上构建基准站和移动站双差观测方程；

短基线情况下的双差观测方程，具体公式如下：

其中，表示根据和求双差而获得的伪距双差观测值，表示卫星i与接收机a之间的伪距观测值；表示卫星j与接收机a之间的伪距观测值；表示卫星i与接收机b之间的伪距观测值；表示卫星j与接收机b之间的伪距观测值；表示根据 和获得的双差值；表示卫星i与接收机a之间的几何距离；表示卫星j与接收机a之间的几何距离；表示卫星i与接收机b之间的几何距离；表示卫星j与接收机b之间的几何距离；表示根据和获得的伪距双差观测值测量噪声；表示卫星i与接收机a之间伪距的测量噪声；表示卫星j与接收机a之间伪距的测量噪声；表示卫星i与接收机b之间伪距的测量噪声；表示卫星j与接收机b之间伪距的测量噪声；表示卫星i在接收机a和b之间的相位单差观测值；表示卫星j在接收机a和b之间的相位单差观测值；表示卫星i频段n的载波波长；表示卫星j频段n的载波波长；表示卫星i在接收机a和b之间的单差模糊度；表示卫星j在接收机a和b之间的单差模糊度；表示根据 和获得的载波相位双差测量值测量噪声；表示卫星i与接收机a之间的载波相位观测值的测量噪声；表示卫星j与接收机a之间的载波相位观测值的测量噪声；表示卫星i与接收机b之间的载波相位观测值的测量噪声；表示卫星j与接收机b之间的载波相位观测值的测量噪声；

步骤2-6、对公式(1)和公式(2)进行线性化，获得误差方程；

具体公式如下：

v＝HX-Z (3)

其中，v为误差项；H表示系数阵；X表示状态向量；Z为OMC值；

步骤2-7、采用kalman滤波对状态向量进行求解。

步骤2-3所述的在忽略IFB影响的基础上，构建基准站和移动站的非差观测方程；

具体公式如下：

其中，c表示光速；dtⁱ表示卫星i的钟差；dt_a表示接收机a的钟差；表示卫星i与接收机a之间的电离层延迟，表示卫星i与接收机a之间的对流层延迟；表示卫星i与接收机a之间的载波相位观测值；表示卫星i与接收机a之间的载波相位观测值所对应的非差模糊度。

步骤2-4所述的根据上一历元的状态向量对当前历元的状态向量进行预测；

预测公式如下：

其中，表示k时刻的状态向量预测值；F_k-1表示k-1时刻系统状态转移矩阵；表示k-1时刻的状态向量；P_k|k-1表示k时刻的状态向量的方差协方差阵预测值；P_k-1|k-1表示k-1时刻的状态向量的方差协方差阵；Q_k-1表示k-1时刻过程噪声的方差协方差阵。

步骤2-7所述的采用kalman滤波对状态向量进行求解；

具体公式如下：

其中，表示k时刻经过滤波计算所获得的状态向量；K_k表示k时刻的滤波增益矩阵；Z_k表示k时刻的OMC值；H_k表示k时刻的系数阵；P_k|k表示k时刻经过滤波计算所获得的状态向量的方差协方差阵；R_k表示k时刻的观测值的方差协方差阵。

步骤3所述的在将以周为单位的单差模糊度参数投影为以周为单位的双差模糊度参数的过程中，利用IFB变化率粒子对所有单差模糊度参数进行改正，采用LAMBDA方法进行双差模糊度固定，将获得的RATIO值作为评价每个粒子的适应度，并使用粒子群优化算法对IFB变化率粒子的群体最优位置进行搜索，最终获得最优IFB变化率估值；

具体步骤如下：

步骤3-1、确定由单差模糊度向双差模糊度投影的投影矩阵，将单差模糊度参数对应的方差协方差阵投影为双差模糊度对应的方差协方差阵；

步骤3-2、将IFB变化率作为粒子群优化算法的搜索目标，对粒子群参数进行设置，包括粒子个数、群体最大迭代次数、最大和最小惯性因数、学习因子、粒子位置变量的最大值和最小值、粒子速度变量的最大值和最小值以及迭代停止条件；

步骤3-3、在指定的范围内，对粒子群进行随机初始化，包括随机位置和速度；

步骤3-4、计算每个粒子的适应度；

步骤3-5、根据粒子适应度更新粒子个体的历史最优位置；

步骤3-6、根据所有粒子个体的历史最优位置更新粒子群体的历史最优位置；

步骤3-7、根据粒子个体的历史最优位置和粒子群体的历史最优位置更新粒子的位置和速度；

步骤3-8、判断更新后的位置或速度是否超出指定的范围，若是，则对该粒子的位置或速度重新进行初始化，并执行步骤3-9；否则直接执行步骤3-9；

步骤3-9、判断当前粒子群体的历史最优位置所对应的适应度是否大于给定阈值或者当前迭代次数是否大于群体最大迭代次数，若是，则当前粒子群体的历史最优位置作为IFB变化率的最优估值，否则，重复执行步骤3-4至步骤3-9，直到满足迭代停止条件。

步骤3-1所述的确定由单差模糊度向双差模糊度投影的投影矩阵，将单差模糊度参数对应的方差协方差阵投影为双差模糊度对应的方差协方差阵；

投影过程公式如下：

DPD^T＝P′ (8)

其中，P′表示投影后状态向量的方差协方差阵；P表示投影前状态向量的方差协方差阵；

投影矩阵如下：

其中，m为观测到的GLONASS卫星个数。

步骤3-4所述的计算每个粒子的适应度；

具体步骤如下：

步骤3-4-1、利用IFB变化率粒子对所有单差模糊度参数进行改正；

改正公式如下：

其中，表示改正后的卫星j在接收机a和b之间的单差模糊度参数；表示改正前的卫星j在接收机a和b之间的单差模糊度参数；k^j表示GLONASS卫星j所对应的频率编号；表示GLONASS卫星j在第n频段上的波长；x_i表示当前IFB变化率粒子的位置；

步骤3-4-2、利用投影矩阵，将改正后的单差模糊度参数投影为双差模糊度参数；

投影过程公式如下：

DX＝X′ (11)

其中，X表示投影前的状态向量，具体公式如下：

公式(12)中，x、y、z为位置参数，该参数在投影前后保持不变；为卫星i在接收机a和b之间的单差模糊度参数，i∈[1，m]，m为观测到的GLONASS卫星个数；

X′表示投影后的状态向量，具体公式如下：

公式(13)中，为投影后的和之间的双差模糊度参数，为卫星1在接收机a和b之间的单差模糊度参数，为卫星j在接收机a和b之间的单差模糊度参数，j∈[2，m]；

步骤3-4-3、在投影后的参数向量及方差协方差阵中提取出与双差模糊度有关的部分，采用LAMBDA方法进行双差模糊度固定，将获取到的RATIO值作为当前粒子的适应度。

本发明优点：

1、本发明采用搜索的方式对IFB进行估计，解决了由于IFB和模糊度线性相关所带来的问题；

2、本发明在不增加待估参数个数以及先验信息条件下，采用粒子群优化算法高效准确搜索出IFB变化率参数；

3、本发明可以实时在线对IFB变化率进行估计，可避免由于观测条件变化导致已有的IFB变化率校正参数不可用现象，适合于实时动态定位；

4、本发明采用实时的GLONASS相位偏差估计方法实现了GLONASS实时模糊度固定，可以进一步拓展GLONASS的应用领域。

附图说明

图1为本发明一种实施方式的基于粒子群优化算法的实时GLONASS相位偏差估计方法流程图；

图2为本发明一种实施方式的步骤2的具体流程图；

图3为本发明一种实施方式的步骤3的具体流程图；

图4为本发明一种实施方式的利用粒子群优化算法单历元估计IFB变化率最优值及所需搜索次数的时间序列示意图；

图5为本发明一种实施方式的利用粒子群优化算法单历元估计IFB变化率最优值所对应的模糊度固定解定位精度统计示意图；

图6为本发明一种实施方式的利用粒子群优化算法逐历元估计IFB变化率最优值及所需搜索次数的时间序列示意图；

图7为本发明一种实施方式的利用粒子群优化算法逐历元估计IFB变化率最优值所对应的模糊度固定解定位精度统计。

具体实施方式

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

本发明实施例中，基于粒子群优化算法的实时GLONASS相位偏差估计方法，如图1所示，包括以下步骤：

步骤1：获取基准站和移动站的GLONASS观测数据，包括GLONASS伪距和载波观测数据，并逐历元进行读取；

本发明实施例中，选用一条8.6km长的基线数据，该基线采用两种不同类型的接收机进行数据采集，两类接收机分别为TRIMBLE NETR9和JPS EGGDT，对应的天线分别为TRM57971.00和AOAD/M_T；基线数据采样间隔为30s，观测时长为11h20min；

步骤2、获取基准站和移动站概略坐标，选择基准站和移动站的共视卫星，对位置坐标、对流层延迟、电离层延迟以及模糊度参数进行初始化或更新，通过对基准站和移动站双差观测方程进行线性化处理的方式获得误差，采用kalman滤波获得单差模糊度参数向量及其方差协方差阵；具体流程图如图2所示，具体步骤如下：

步骤2-1、采用伪距单点定位(SPP)方法获得基准站和移动站的概率坐标(-2491490.2616，-4660803.2317，3559129.0005)和(-2493304.6796，-4655215.1032，3565497.5918)；

步骤2-2、根据基准站和移动站的观测数据选择共视卫星；

共视卫星的数量会随着时间变化而变化，在本实施中，共视卫星数量的变化范围为[5，9]颗；

步骤2-3、在忽略IFB影响的基础上，构建基准站和移动站的非差观测方程；

具体公式如下：

其中，c表示光速；dtⁱ表示卫星i的钟差；dt_a表示接收机a的钟差；表示卫星i与接收机a之间的电离层延迟，表示卫星i与接收机a之间的对流层延迟；表示卫星i与接收机a之间的载波相位观测值；表示卫星i与接收机a之间的载波相位观测值所对应的非差模糊度。

步骤2-4、判断当前是否为首历元，若是，则对包含位置、对流层延迟、电离层延迟以及模糊度参数的状态向量及其方差协方差阵进行初始化，否则根据上一历元的状态向量对当前历元的状态向量进行预测；预测公式如下：

其中，表示k时刻的状态向量预测值；F_k-1表示k-1时刻系统状态转移矩阵；表示k-1时刻的状态向量；P_k|k-1表示k时刻的状态向量的方差协方差阵预测值；P_k-1|k-1表示k-1时刻的状态向量的方差协方差阵；Q_k-1表示k-1时刻过程噪声的方差协方差阵。

步骤2-5、在步骤2-3的基础上构建基准站和移动站双差观测方程；

短基线情况下的双差观测方程，具体公式如下：

其中，表示根据和求双差而获得的伪距双差观测值，表示卫星i与接收机a之间的伪距观测值；表示卫星j与接收机a之间的伪距观测值；表示卫星i与接收机b之间的伪距观测值；表示卫星j与接收机b之间的伪距观测值；表示根据 和获得的双差值；表示卫星i与接收机a之间的几何距离；表示卫星j与接收机a之间的几何距离；表示卫星i与接收机b之间的几何距离；表示卫星j与接收机b之间的几何距离；表示根据和获得的伪距双差观测值测量噪声；表示卫星i与接收机a之间伪距的测量噪声；表示卫星j与接收机a之间伪距的测量噪声；表示卫星i与接收机b之间伪距的测量噪声；表示卫星j与接收机b之间伪距的测量噪声；表示卫星i在接收机a和b之间的相位单差观测值；表示卫星j在接收机a和b之间的相位单差观测值；表示卫星i频段n的载波波长；表示卫星j频段n的载波波长；表示卫星i在接收机a和b之间的单差模糊度；表示卫星j在接收机a和b之间的单差模糊度；表示根据 和获得的载波相位双差测量值测量噪声；表示卫星i与接收机a之间的载波相位观测值的测量噪声；表示卫星j与接收机a之间的载波相位观测值的测量噪声；表示卫星i与接收机b之间的载波相位观测值的测量噪声；表示卫星j与接收机b之间的载波相位观测值的测量噪声；

步骤2-6、对公式(1)和公式(2)进行线性化，获得误差方程；

具体公式如下：

v＝HX-Z (3)

其中，v为误差项；H表示系数阵；X表示状态向量；Z为OMC值；

步骤2-7、采用kalman滤波对状态向量进行求解。

具体公式如下：

其中，表示k时刻经过滤波计算所获得的状态向量；K_k表示k时刻的滤波增益矩阵；Z_k表示k时刻的OMC值；H_k表示k时刻的系数阵；P_k|k表示k时刻经过滤波计算所获得的状态向量的方差协方差阵；R_k表示k时刻的观测值的方差协方差阵。

步骤3、在将以周为单位的单差模糊度参数投影为以周为单位的双差模糊度参数的过程中，利用IFB变化率粒子对所有单差模糊度参数进行改正，采用LAMBDA方法进行双差模糊度固定，将获得的RATIO值作为评价每个粒子的适应度，并使用粒子群优化算法对IFB变化率粒子的群体最优位置进行搜索，最终获得最优IFB变化率估值；具体流程图如图3所示，具体步骤如下：

步骤3-1、确定由单差模糊度向双差模糊度投影的投影矩阵，将单差模糊度参数对应的方差协方差阵投影为双差模糊度对应的方差协方差阵；投影过程公式如下：

DPD^T＝P′ (8)

其中，P′表示投影后状态向量的方差协方差阵；P表示投影前状态向量的方差协方差阵；

投影矩阵如下：

其中，m为观测到的GLONASS卫星个数。本发明实施例中，m的取值范围为[5，9]；

步骤3-2、将IFB变化率作为粒子群优化算法的搜索目标，对粒子群参数进行设置；

本发明实施例中，设置粒子个数设为10；群体最大迭代次数设为10；最大和最小惯性因数分别设为0.2和1.2；学习因子c₁和c₂均设为2；粒子位置变量的最大值和最小值分为设为-0.1和0.1；粒子速度变量的最大值和最小值分别为0.03和0；迭代停止条件中gbest所对应的适应度阈值设为4.5；

步骤3-3、在[-0.1，0.1]范围内，对粒子群进行随机初始化，包括随机位置和速度；

步骤3-4、计算每个粒子的适应度；具体步骤如下：

步骤3-4-1、利用IFB变化率粒子对所有单差模糊度参数进行改正；

改正公式如下：

其中，表示改正后的卫星j在接收机a和b之间的单差模糊度参数；表示改正前的卫星j在接收机a和b之间的单差模糊度参数；k^j表示GLONASS卫星j所对应的频率编号；表示GLONASS卫星j在第n频段上的波长；x_i表示当前IFB变化率粒子的位置；

步骤3-4-2、利用投影矩阵，将改正后的单差模糊度参数投影为双差模糊度参数；

投影过程公式如下：

DX＝X′ (11)

其中，X表示投影前的状态向量，具体公式如下：

公式(12)中，x、y、z为位置参数，该参数在投影前后保持不变；为卫星i在接收机a和b之间的单差模糊度参数，i∈[1，m]，m为观测到的GLONASS卫星个数；

X′表示投影后的状态向量，具体公式如下：

公式(13)中，为投影后的和之间的双差模糊度参数，为卫星1在接收机a和b之间的单差模糊度参数，为卫星j在接收机a和b之间的单差模糊度参数，j∈[2，m]；

步骤3-4-3、在投影后的参数向量及方差协方差阵中提取出与双差模糊度有关的部分，采用LAMBDA方法进行双差模糊度固定，将获取到的RATIO值作为当前粒子的适应度。

步骤3-5、根据粒子适应度更新粒子个体的历史最优位置；

步骤3-6、根据所有粒子个体的历史最优位置更新粒子群体的历史最优位置；

步骤3-7、根据粒子个体的历史最优位置和粒子群体的历史最优位置更新粒子的位置和速度；公式如下：

式(14)中：

其中，是第i个粒子在第k次迭代时的位置；是第i个粒子在第k次迭代时的速度；是第i个粒子在第k+1次迭代时的位置；是第i个粒子在第k+1次迭代时的速度；w^k是第i个粒子在第k次迭代时的惯性因子；f_i ^k是第i个粒子在第k次迭代时的适应度，f_thres是为判断粒子适应度是否可信而设置的一个阈值，w_min和w_max分别是设定的最小和最大惯性因子；c₁和c₂是设定的学习因子；r₁和r₂为[0，1]区间的均匀分布随机函数；pbest_i是第i个粒子的历史最优位置；gbest是所有粒子的群体历史最优位置；

步骤3-8、判断更新后的位置是否超出[-0.1，0.1]范围或速度是否超出[0，0.03]范围，若是，则对该粒子的位置或速度重新进行初始化，并执行步骤3-9；否则直接执行步骤3-9；

步骤3-9、判断当前粒子群体的历史最优位置所对应的适应度是否大于给定阈值4.5或者当前迭代次数是否大于群体最大迭代次数10，若是，则当前粒子群体的历史最优位置作为IFB变化率的最优估值，否则，重复执行步骤3-4至步骤3-9，直到满足迭代停止条件。

步骤4、将确定出来的接收机端相位观测值的IFB变化率估值作为改正数对单差模糊度参数进行改正，采用LAMBDA方法进行双差模糊度固定，即完成对GLONASS模糊度的实时固定；具体如下：

步骤4-1、将确定出来的IFB变化率估值作为改正数对GLONASS单差模糊度参数进行改正；

步骤4-2、使用投影矩阵D，将改正后的单差模糊度参数及其方差协方差阵投影为双差模糊度参数及其方差协方差阵；

步骤4-3、采用LAMBDA方法进行双差模糊度固定；

步骤4-4、使用RATIO值对模糊度固定解进行验证，如果RATIO值大于3，那么说明当前历元模糊度固定成功，利用固定的双差模糊度值对其它非模糊度参数进行改正，获得当前历元的模糊度固定解，同时根据参考卫星的单差模糊度值以及固定的双差模糊度值对投影前的单差模糊度进行改正，为下一历元数据处理做准备；如果RATIO值小于3，说明当前历元无法模糊度固定，以浮点解作为当前历元的最终解；

步骤4-5、如果连续n(n≥5)个历元均为模糊度固定解，并且这n个历元所求得的n个IFB变化率估值十分接近，它们的标准差小于0.004，那么缩小下一历元粒子群的搜索范围，进一步提高粒子群搜索效率；

本发明实施例中，将本发明所述方法与基于粒子滤波的IFB变化率估计方法(TianY，Ge M，Neitzel F.Particle filter-based estimation of inter-frequency phasebias for real-time GLONASS integer ambiguity resolution：Journal ofGeodesy，2015，89(11)：1145-1158)进行比较，利用本发明所述方法采用单历元数据处理方式对观测数据进行处理；

如图4所示，在观测时段内，1360个历元中有98.6％的历元能够准确搜索出IFB变化率，只有极少历元无法准确确定出对应的IFB变化率估值，对于这些历元，其搜索次数已经达到最大搜索次数，并且对应的适应度也均小于3，说明无论采用何种方式进行搜索，均无法使这些历元的模糊度固定，出现这类现象的可能原因为这些历元的观测数据质量不佳，即使给定准确的IFB变化率改正值，采用单历元数据也无法进行模糊度固定，因而也无法准确搜索出IFB变化率。采用本发明所述方法进行IFB变化率搜索，每历元平均搜索次数为32次，而对于基于粒子滤波的IFB变化率估计方法，首先该方法无法实现单历元IFB变化率估计，通常需要一定数量的历元进行收敛；其次该方法的每历元平均搜索次数至少为200次。该实验中单历元模糊度固定成功率达到96.2％，为了进一步验证模糊度固定的可靠性，将模糊度固定解的定位精度进行分析。

如图5所示，将模糊度固定解的定位结果与其真值进行比较，其在东(E)、北(N)、高(U)三个方向上的最大坐标偏差分别为11.7、16.7和36.7mm，平均定位精度分别为0.3、0.1和0.5mm，中误差分别为4.3、4.8和11.1mm；根据以上的实验结果统计分析表明使用本发明所述方法可以高效准确的获得IFB变化率估值，同时也可以获得可靠的模糊度固定解。

单历元数据处理方式可以很好的检验算法的性能，但在实际数据处理过程中，通常采用滤波的方式对数据进行逐历元处理。利用本发明所述方法对上述观测数据进行逐历元数据处理；粒子群参数设置如单历元数据处理方式。

如图6所示，在观测时段内，所有历元都能够准确搜索出IFB变化率，其平均值为6.7mm/FN，标准差为2mm/FN；每历元平均搜索次数为9次；模糊度固定成功率为97.8％。相较于单历元数据处理方式，采用滤波方式的逐历元数据处理可以更好的发挥本发明所述方法的性能。需要说明的是：相对于使用最小二乘和基于粒子滤波的IFB变化率估计方法，使用本发明所述方法获取的IFB变化率并不够稳定，这是因为在使用粒子群优化算法进行IFB变化率最优值搜索过程中，需要同时顾及解算效率和精度，考虑到实时定位对解算效率的要求，我们需要的并不是绝对群体最优值，而是既能够满足模糊度固定所需基本条件又能够达到搜索次数最少效果的近似最优值，因此利用本发明所述方法获取到的IFB变化率估值虽然不是十分稳定，但可以达到较好的模糊度固定效果。

为了对本次实验中的模糊度固定解的定位精度进行分析，将模糊度固定解的定位结果与其真值进行比较，如图7所示，在E、N、U三个方向上的最大坐标偏差分别为11.7、14.4和31.4mm，平均定位精度分别为0.1、0.0和0.1mm，中误差分别为4.2、4.6和11.5mm，以上统计结果说明利用本专利所述方法获取的模糊度固定解十分可靠。

以上所述，仅为本发明中最基础的具体实施方式，但本发明的保护范围并不局限于此，任何本技术领域人士在本发明所揭露的技术范围内，可理解到的替换，都应涵盖在本发明的包含范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (8)

1.一种基于粒子群优化算法的实时GLONASS相位偏差估计方法，其特征在于，包括以下步骤：

步骤1：获取基准站和移动站的GLONASS观测数据，包括GLONASS伪距和载波观测数据，并逐历元进行读取；

步骤2、获取基准站和移动站概略坐标，选择基准站和移动站的共视卫星，对位置坐标、对流层延迟、电离层延迟以及模糊度参数进行初始化或更新，对基准站和移动站双差观测方程进行线性化处理，采用kalman滤波获得单差模糊度参数向量及其方差协方差阵；

步骤3、在将以周为单位的单差模糊度参数投影为以周为单位的双差模糊度参数的过程中，利用IFB变化率粒子对所有单差模糊度参数进行改正，采用LAMBDA方法进行双差模糊度固定，将获得的RATIO值作为评价每个粒子的适应度，并使用粒子群优化算法对IFB变化率粒子的群体最优位置进行搜索，最终获得最优IFB变化率估值；

步骤4、将确定出来的接收机端相位观测值的IFB变化率估值作为改正数对单差模糊度参数进行改正，采用LAMBDA方法进行双差模糊度固定，即完成对GLONASS模糊度的实时固定。

2.根据权利要求1所述的基于粒子群优化算法的实时GLONASS相位偏差估计方法，其特征在于，步骤2所述的获取基准站和移动站概略坐标，选择基准站和移动站的共视卫星，对位置坐标、对流层延迟、电离层延迟以及模糊度参数进行初始化或更新，对基准站和移动站双差观测方程进行线性化处理，采用kalman滤波获得单差模糊度参数向量及其方差协方差阵；

具体步骤如下：

步骤2-1、采用伪距单点定位方法获得基准站和移动站的概率坐标；

步骤2-2、根据基准站和移动站的观测数据选择共视卫星；

步骤2-3、在忽略IFB影响的基础上，构建基准站和移动站的非差观测方程；

步骤2-4、判断当前是否为首历元，若是，则对包含位置、对流层延迟、电离层延迟以及模糊度参数的状态向量及其方差协方差阵进行初始化，否则根据上一历元的状态向量对当前历元的状态向量进行预测；

步骤2-5、在步骤2-3的基础上构建基准站和移动站双差观测方程；

短基线情况下的双差观测方程，具体公式如下：

其中，表示根据和求双差而获得的伪距双差观测值，表示卫星i与接收机a之间的伪距观测值；表示卫星j与接收机a之间的伪距观测值；表示卫星i与接收机b之间的伪距观测值；表示卫星j与接收机b之间的伪距观测值；表示根据和获得的双差值；表示卫星i与接收机a之间的几何距离；表示卫星j与接收机a之间的几何距离；表示卫星i与接收机b之间的几何距离；表示卫星j与接收机b之间的几何距离；表示根据和获得的伪距双差观测值测量噪声；表示卫星i与接收机a之间伪距的测量噪声；表示卫星j与接收机a之间伪距的测量噪声；表示卫星i与接收机b之间伪距的测量噪声；表示卫星j与接收机b之间伪距的测量噪声；表示卫星i在接收机a和b之间的相位单差观测值；表示卫星j在接收机a和b之间的相位单差观测值；表示卫星i频段n的载波波长；表示卫星j频段n的载波波长；表示卫星i在接收机a和b之间的单差模糊度；表示卫星j在接收机a和b之间的单差模糊度；表示根据和获得的载波相位双差测量值测量噪声；表示卫星i与接收机a之间的载波相位观测值的测量噪声；表示卫星j与接收机a之间的载波相位观测值的测量噪声；表示卫星i与接收机b之间的载波相位观测值的测量噪声；表示卫星，与接收机b之间的载波相位观测值的测量噪声；

步骤2-6、对公式(1)和公式(2)进行线性化，获得误差方程；

具体公式如下：

v＝HX-Z (3)

其中，v为误差项；H表示系数阵；X表示状态向量；Z为OMC值；

步骤2-7、采用kalman滤波对状态向量进行求解。

3.根据权利要求2所述的基于粒子群优化算法的实时GLONASS相位偏差估计方法，其特征在于，步骤2-3所述的在忽略IFB影响的基础上，构建基准站和移动站的非差观测方程；

具体公式如下：

其中，c表示光速；dtⁱ表示卫星i的钟差；dt_a表示接收机a的钟差；表示卫星i与接收机a之间的电离层延迟，表示卫星i与接收机a之间的对流层延迟；表示卫星i与接收机a之间的载波相位观测值；表示卫星i与接收机a之间的载波相位观测值所对应的非差模糊度。

4.根据权利要求2所述的基于粒子群优化算法的实时GLONASS相位偏差估计方法，其特征在于，步骤2-4所述的根据上一历元的状态向量对当前历元的状态向量进行预测；

预测公式如下：

其中，表示k时刻的状态向量预测值；F_k-1表示k-1时刻系统状态转移矩阵；表示k-1时刻的状态向量；P_k|k-1表示k时刻的状态向量的方差协方差阵预测值；P_k-1|k-1表示k-1时刻的状态向量的方差协方差阵；Q_k-1表示k-1时刻过程噪声的方差协方差阵。

5.根据权利要求2所述的基于粒子群优化算法的实时GLONASS相位偏差估计方法，其特征在于，步骤2-7所述的采用kalman滤波对状态向量进行求解；

具体公式如下：

其中，表示k时刻经过滤波计算所获得的状态向量；K_k表示k时刻的滤波增益矩阵；Z_k表示k时刻的OMC值；H_k表示k时刻的系数阵；P_k|k表示k时刻经过滤波计算所获得的状态向量的方差协方差阵；R_k表示k时刻的观测值的方差协方差阵。

6.根据权利要求1所述的基于粒子群优化算法的实时GLONASS相位偏差估计方法，其特征在于，步骤3所述的在将以周为单位的单差模糊度参数投影为以周为单位的双差模糊度参数的过程中，利用IFB变化率粒子对所有单差模糊度参数进行改正，采用LAMBDA方法进行双差模糊度固定，将获得的RATIO值作为评价每个粒子的适应度，并使用粒子群优化算法对IFB变化率粒子的群体最优位置进行搜索，最终获得最优IFB变化率估值；

具体步骤如下：

步骤3-1、确定由单差模糊度向双差模糊度投影的投影矩阵，将单差模糊度参数对应的方差协方差阵投影为双差模糊度对应的方差协方差阵；

步骤3-2、将IFB变化率作为粒子群优化算法的搜索目标，对粒子群参数进行设置，包括粒子个数、群体最大迭代次数、最大和最小惯性因数、学习因子、粒子位置变量的最大值和最小值、粒子速度变量的最大值和最小值以及迭代停止条件；

步骤3-3、在指定的范围内，对粒子群进行随机初始化，包括随机位置和速度；

步骤3-4、计算每个粒子的适应度；

步骤3-5、根据粒子适应度更新粒子个体的历史最优位置；

步骤3-6、根据所有粒子个体的历史最优位置更新粒子群体的历史最优位置；

步骤3-7、根据粒子个体的历史最优位置和粒子群体的历史最优位置更新粒子的位置和速度；

步骤3-8、判断更新后的位置或速度是否超出指定的范围，若是，则对该粒子的位置或速度重新进行初始化，并执行步骤3-9；否则直接执行步骤3-9；

步骤3-9、判断当前粒子群体的历史最优位置所对应的适应度是否大于给定阈值或者当前迭代次数是否大于群体最大迭代次数，若是，则当前粒子群体的历史最优位置作为IFB变化率的最优估值，否则，重复执行步骤3-4至步骤3-9，直到满足迭代停止条件。

7.根据权利要求6所述的基于粒子群优化算法的实时GLONASS相位偏差估计方法，其特征在于，步骤3-1所述的确定由单差模糊度向双差模糊度投影的投影矩阵，将单差模糊度参数对应的方差协方差阵投影为双差模糊度对应的方差协方差阵；

投影过程公式如下：

DPD^T＝P′ (8)

其中，P′表示投影后状态向量的方差协方差阵；P表示投影前状态向量的方差协方差阵；

投影矩阵如下：

其中，m为观测到的GLONASS卫星个数。

8.根据权利要求7所述的基于粒子群优化算法的实时GLONASS相位偏差估计方法，其特征在于，步骤3-4所述的计算每个粒子的适应度；

具体步骤如下：

步骤3-4-1、利用IFB变化率粒子对所有单差模糊度参数进行改正；

改正公式如下：

其中，表示改正后的卫星j在接收机a和b之间的单差模糊度参数；表示改正前的卫星j在接收机a和b之间的单差模糊度参数；k^j表示GLONASS卫星j所对应的频率编号；表示GLONASS卫星j在第n频段上的波长；x_i表示当前IFB变化率粒子的位置；

步骤3-4-2、利用投影矩阵，将改正后的单差模糊度参数投影为双差模糊度参数；

投影过程公式如下：

DX＝X′ (11)

其中，X表示投影前的状态向量，具体公式如下：

公式(12)中，x、y、z为位置参数，该参数在投影前后保持不变；为卫星i在接收机a和b之间的单差模糊度参数，i∈[1，m]，m为观测到的GLONASS卫星个数；

X′表示投影后的状态向量，具体公式如下：

公式(13)中，为投影后的和之间的双差模糊度参数，为卫星1在接收机a和b之间的单差模糊度参数，为卫星j在接收机a和b之间的单差模糊度参数，j∈[2，m]；

步骤3-4-3、在投影后的参数向量及方差协方差阵中提取出与双差模糊度有关的部分，采用LAMBDA方法进行双差模糊度固定，将获取到的RATIO值作为当前粒子的适应度。