基于MATLAB的抗差多项式拟合方法在GNSS周跳探测中的应用研究

于佳慧 张世涛 张岩

摘  要： 为了提高GNSS载波相位的定位精度，本文采用历元间求差与抗差估计相结合的方法对过去常常使用的多项式进行改良，最终提出用抗差多项式拟合来探测并修复周跳。本文基于Matlab软件编制程序，并对GNSS载波相位周跳进行了试算。实验结果表明：此方法可以做到对GNSS周跳探测与修复功能，具有文件读写、卫星观测值提取、周跳探测与修复、绘制拟合的误差结果图形等功能。通过这种方法可以探测并修复1周以上的周跳。得出结论：通过抗差多项式拟合方法可以探测并修复GNSS周跳。

关键词： MATLAB;GNSS;周跳的探测与修复;抗差多项式

中图分类号： TP3    文献标识码： A    DOI：10.3969/j.issn.1003-6970.2019.08.041

本文著录格式：于佳慧，张世涛，张岩. 基于MATLAB的抗差多项式拟合方法在GNSS周跳探测中的应用研究[J]. 软件，2019，40（8）：175180

【Abstract】： In order to improve the positioning accuracy of GNSS carrier phase， this paper uses the method of inter-calendar difference and robust estimation to improve the traditional multinomial， and finally proposes a method of using robust multinomial fitting to detect and repair cycle slip. In this paper， based on Matlab software， the phase cycle jump of GNSS carrier is calculated. The experimental results show that this method can detect and repair GNSS cycle slip， and has the functions of file reading and writing， satellite observation data extraction， cycle slip detection and repair， drawing fitting error result figure and so on. Through this program， the cycle jump can be detected and repaired for more than 1 week. Draw Conclusion： the robust polynomial fitting method can be used in the detection and repair of GNSS cycle slip.

【Key words】： MATLAB; GNSS; Detection and repair; Robust Polynomials

0  引言

GNSS具有使用观测精度高、测量时间短、全天候、技术自动化程度高、观测范围大、能提供三维坐标测量成果等特点[1]。在GPS数据处理过程中，周跳的存在会使观测值中出现一个偏差，这会使观测值失真，若不修复周跳，剔除这些错误成果，会对测绘工作造成严重后果。因此，要想提高定位精度就必须能准确地探测周跳并且修复载波相位。目前可用于接收机的周跳探测方法主要有高次差法、电离层求差法、伪距和载波相位观测值组合、多項式拟合法来探测和修复周跳等。

总之，有很多方法可以探测和修复周跳，但不同的周跳探测与修复方法对提高计算效率及可靠性都有所不同。常规方法对载波相位观测值周跳探测和修复虽然都有其优点，但也存在其局限性，鉴于此，本文将以多项式拟合法为基础，先进行时间标准化，然后采用历元间求差与抗差估计相结合的方法对传统算法中存在的问题进行改进，并利用GNSS/GLONASS实测数据验证改进后的算法的有效性[2]。

1  MATLAB简介

1.1  什么是MATLAB

MATLAB软件是由美国Math Works公司推出的一种交互式、面向对象的程序设计语言，广泛用于工业界和学术界。不仅能很好的处理数值问题，同时还能解决符号问题，绘制图形。

1.2  MATLAB所具有的优势

MATLAB系统包括以下五个部分：

（1）MATLAB具有多种语言体系。在MATLAB可用的编程语言，对应于工具箱中六个目录，分别为ops（操作符与特殊字符）、lang（编程语言结构）、strfun（字符串操作）、iofun（文件的输入输出）、time fun（时间和日期）和data type（数据类型和结构）[3]。

（2） MATLAB具有完整工作界面。包括MATLAB编程和调试的环境，对工作空间中的变量进行管理及进行输入输出的数据[4]。

（3）MATLAB具有强大的图形处理能力。其中包括二、三维图形可视化，也包括定义图形外观的操作，最终可以为MATLAB应用程序创建整套的用户界面。

（4）数学函数库。包括初等函数;也包括更复杂的功能，如矩阵功能（矩阵求逆、求矩阵的特征值）、快速傅里叶变换等[3]。

（5）MATLAB应用程序接口（API）。提供接口程序使程序员可以编写C/C6、Fortran程序调用MATLAB作为计算引擎，并用MATLAB调用外部应用程序或者读写MATLAB文件[4]。MATLAB中的M文件的语法与其他的高级语言类似，但语法较其他一般的高级语言来说要简单，它只是一个ASCII码简单的文本文件，其语法编制与数学语言非常接近。它不仅是一种程序化的编程语言，也是一种解释性的编程语言，可逐行解释和运行程序，程序更易调试[4]。

2  抗差多项式拟合法

2.1  周跳概念和产生的原因

载波相位可以用在高精度卫星导航定位。完整且连续的周期通过使用多普勒计数完成，而载波相位在导航定位中主要是测量不完整的周期。信噪比低、信号质量差以及接收机接收到的质量差等方面都可能引起完整的周期部分突变，叫作整周跳变，简称为周跳（cycle slip）。

载波相位是由 、 、 这些部分组成的。虽然可以以极高的精度测定 ，但也要正确地测定 和 的情况下，才能确定载波相位值。在载波相位测量的实际观测值中，小于一周的部分称为观测时刻的观测值。卫星载波信号与接收机参考振荡信号形成的差频信号不足一周。接收机载波跟踪回路中的鉴相器测量可以实现对载波信号的测量。只有接收机的载波信号和参考振荡信号能在观测时段产生差频信号，才可以得到正确的观测值。整个星期的计数不一样。它是计数器从第一个观测时间到当前观测时间积累的差分频率信号中的整数频带数。如果由于外界因素，通过使用多普勒计数法的计数器在两个观测周期之间的某个时段停止计数，从而使整个周期的计数比实际值少n周，那么当计数器恢复正常运行时，所有载波相位波段将包含相同的误差值，该误差值小于整个周期的正常值（如图1）。这样一个系统性的误差，在一个周期内计数和不足一周的部分保持正确计数的情况，称为整周跳变，简称周跳[5]。

2.2  抗差估计求解拟合系数

传统的多项式拟合方法一般计算拟合系数是采用最小二乘法，但当观测值不服从正态分布的假设时，即m个拟合观测值出现粗差或周跳时，最小二乘法估计就不具有抗干扰性，错误的估计结果可能由单个观测值的偏差导致。而利用抗差估计法与利用最小二乘估计法不同，它不像最小二乘那样过分的追求有效性与无偏性等内部性质，抗差估计则更加重视估值的实际可靠性和抗差性[6]。从第一个历元开始，先一次取m个历元 和对应的载波相位观测值 代入式（1），通过抗差估计法解算其拟合系数。假设初始权矩阵为单位阵，即 ，之后进行平差计算。根据第一次平差后得到的改正数阵 和单位权中误差 ，利用下式（2）

对各相位观测值重新定权，如果 则认为m个历元观测值中不存在周跳，最后的平差结果与第一次平差结果相同。若有偏差，可利用权阵 在再次进行平差计算，直至 时停止迭代，取第j次迭代结果作为最终结果。这样可以确保参加拟合的m个相位观测值为不含周跳的值，從而确保得到的拟合系数更加准确和可靠。

2.3  利用抗差多项式探测和修复周跳

根据求得的拟合系数拟合，并带入m+1个历元值，来估计第m+1个历元的相位值，将利用抗差多项式拟合出来的拟合值与对应历元的实际观测作差，若差值 时，则认为 不含周跳，去掉本次拟合的第一个观测值，加入第m+1个历元观测值继续进行拟合;当差值 时，则认为 含有周跳，采用外推的整周计数去取代有周跳的实际观测值中的整周计数，但不足一周的部分仍然保持不动，然后继续上述过程，直到最后一个相位观测值为止[7]。

本部分首先对高次差法的原理和适用范围优缺点进行了介绍，其次介绍了利用传统多项式探测周跳原理，在此基础上进行时间标准化、历元间求差、重新定权等改进，从而建立抗差多项式模型，利用抗差多项式探测和修复周跳。

3  编制周跳探测与修复系统程序

3.1  程序的设计框架

研究利用多项式[8]拟合法周跳探测和修复方法，对传统的多项式拟合法周跳探测和修复方法进行改进，并确定多项式拟合的阶数。

在无周跳“干净”观测数据中加入周跳，通过差分技术对相邻历元之间的观测数据进行差分，对差分数据进行多项式建模拟合与预报，探测周跳;利用抗差多项式建模探测周跳，即编制抗差多项式程序设计框架（如图2）。

3.2  周跳探测程序算例分析

1. 解决多项式的基础问题

卫星到地面的距离对时间的四阶以上导数已趋于零，其变化规律是随机的，无法再用多项式进行拟合。故多项式拟合的阶数n取3～4阶，在本实验中n=3的条件下既能保证精度又保证了运算速度。

拟合窗口m选取时，当拟合窗口宽度越大时，外推值越准确，拟合中误差越小，但过大又会增加计算量。当拟合窗口宽度越小时，外推值精度越差。所以通过取不同的值进行实验，根据实验取m=10比较合适。

门槛系数k的取值是随不同的情况改变的。当k较大时，表示只有当观测值偏离外推值很大的时候才认为它是异常的;当k较小时，表示当观测值偏离外推值较小的时候就认为它是异常值了。一般取k在2到9之间的数。经过多次取值实验，结果表明门槛系数取k=4时效果最好。

2. 高次差法检验观测数据质量

算例基于2014年1月14日国内的GPS/GLONASS的实测数据，观测时间为10分钟，采样间隔为1秒。使用高阶差分的方法先对16时45分0秒到16时55分0秒GPS的10号卫星共600个L1载波的据进行周跳探测，对观测值求四阶差分的时间序列图（如图3），从（如图3）可以看出，实验结果中没有周跳[9]。

3. 多项式拟合法算例分析

本文选取上述无周跳的数据中从16时45分0秒到16时45分60秒的实测GPS数据进行多项式拟合实验，以第3秒（即第4个历元）到第12秒（即第13个历元）的数据作为拟合多项式的数据，外推出从第13秒到第60秒的数据，并设计实验方案：方案1，传统的多项式拟合法;方案2，历元间求差并采用最小二乘[10]平差的多项式拟合法;方案3，历元间求差并采用抗差估计的多项式拟合法。

在无周跳情况下（如图4），三种方案拟合能力的比较，前一段时间曲线变化比较平缓，较符合实际情况;随着观测历元数的增大，方案一与方案二曲线的波动越来越大，出现发散的情况，有可能造成周跳判断失误。方案三可明显看出采用抗差估计后每个历元的拟合中误差都有较大改善，GPS观测量拟合值与实际值之差在0.01周以内。

人为在第6秒（即第7个观测历元）上加上1周周跳（如图5）三种方案拟合效果的比较，前一段时间3条曲线变化比较平缓;随着观测历元数的增大，方案一与方案二拟合值与实际值之差越来越大，且方案一中拟合值与实际值之差最大，甚至达到了800周，这都是由于在参加拟合的m个观测值中存在周跳影响的，而方案三可明显看出采用抗差多项式拟合法后每个历元的拟合值与实际值之差都较小且在0.01周以内，这是因为方案三可对参加拟合的数据进行筛选，对存在周跳的数据不用事先利用高次差法判定后，予以剔除，这样可以确保参加拟合的m个观测值为“干净”的，得到的可视化[11]拟合系数也更加准确和可靠。 

人为在非参加拟合的多项式数据中第14秒（即第15个观测历元）上加上1周周跳（如图6）。

三种方案拟合能力的比较，前一段时间3条曲线变化比较平缓，较符合实际情况，并在第15个观测历元上发现周跳，探测结果（如表1）;随着观测历元数的增大，方案一与方案二拟合中误差越来越大，然而方案三的拟合误差在0.01周，这表明：利用抗差多项式法即可探测出周跳的位置和大小，又可达到修复周跳的目的。

人为在第14秒（即第15个观测历元）上加上1周周跳、在第20秒上加上5周周跳、在第30秒上加上10周周跳（如图7），三种方案拟合能力的比较，前一段时间3条曲线变化较符合实际情况，并在第14、20、30秒观测时刻发现周跳。随着观测历元数的增大，方案一与方案二拟合中误差越来越大，虽然有连续周跳却不影响其拟合精度，然而方案三的拟合误差在有连续周跳情况下仍然控制在0.01周内，这表明：利用抗差多项式法也可探测出连续周跳的位置和大小，见下（如表2）并对其周跳进行修复[12]。

4  结论

本文介绍了周跳产生的原因和周跳对定位精度的影响，在此基础上，以抗差多项式拟合法周跳探測与修复为研究对象，得出以下结论：采用抗差估计多项式拟合法，由于对传统多项式拟合进行了重新定权的改进，将含有周跳的拟合数据剔除，所以外推出来的拟合值与实际值之差会很小。在参加拟合的数据不存在周跳的情况下，显然采用抗差多项式拟合法探测出周跳的能力更好;由于参加拟合的数据不存在周跳，所以拟合值不受加入周跳的影响，而采用抗差多项式拟合法能将周跳修复。

采用抗差估计多项式拟合法对传统多项式拟合进行了改进，解决了传统多项式拟合中拟合值与实际观测值之差发散的问题，结合抗差估计保证了估值的实际抗差性和可靠性[13]。通过实测数据进行实验，结果验证了该方法探测与修复周跳的效果比传统的多项式拟合法和在传统多项式拟合基础上只进行历元间求差的方法更准确可靠。

参考文献

[1] 周忠谟， 易杰军， 周琪编著. GPS卫星测量原理与应用[M].湖北武汉. 测绘出版社， 1995： 90-91

[2] 李明， 高星伟， 徐爱功. 一种改进的周跳多项式拟合方法[J]. 测绘科学， 2008， 33（4）： 82-83.

[3] 孙晓， 王勇， 王宝山. MATLAB软件在测量平差中的应用[J]. 焦作工学院学（自然科学版）， 2002（5）.

[4] 蔡旭辉， 刘卫国， 蔡立艳. MATLAB基础与应用教程[M].人民邮电出版社， 2009： 1-2.

[5] 关昊. GPS载波相位周跳探测与修复方法的研究[D]. 辽宁工程技术大学， 2011.

[6] 王福丽， 成英燕， 韦铖， 王晓明. 利用抗差多项式拟合法探测修GNSS周跳[J]. 大地测量与地球动力学， 2013（3）： 129-132.

[7] 邢会敏， 杨福芹. 基于改进多项式拟合法的周跳探测与修复[J]. 科技信息， 2011， （22）： 108-109.

[8] 胡梦英， 贺祖国. 基于重开始共轭思想的改进多项式插值法[J]. 软件， 2015， 36（11）： 48-51

[9] BISNATH S B. Efficient automated cycle-slip correction of dual-frequency kinem atic GPS data[A]. Proceedings of ION GPS 2000. the l3th International Technical Meeting of The Institute of Navigation[C]， Salt Lake City， Utah， 2000， 145- 154.

[10] 聂敬云， 李春青， 李威威， 等. 关于遗传算法优化的最小二乘支持向量机在MBR仿真预测中的研究[J]. 软件， 2015， 36（5）： 40-44.

[11] 孙彦超， 章坚民. 基于MATLAB的DG选址定容可视化系统的设计与实现[J]. 软件， 2018， 39（4）： 142-147.

[12] 胡云康， 姜苏， 吴志荣， 等. 基于改进的纹理合成图像修复算法[J]. 软件， 2016， 37（4）： 60-63.

[13] 史小玲. 集团网网络可靠性的探讨及实例研究[J]. 软件， 2016， 37（01）： 117-119.