马晓东1,2， 宋冰1,2， 张书毕1,2， 刘鑫1,2
(1.中国矿业大学 环境与测绘学院， 江苏 徐州 221116；
2.江苏省资源环境信息工程重点实验室， 江苏 徐州 221116)
摘要：在应用北斗单频接收机对矿区进行实时变形监测时，为提高北斗单频载波相位观测值预处理中周跳探测的精度，提出了基于Kalman滤波的单频北斗载波相位实时周跳探测方法，论述了粗差与周跳的识别方法。利用某矿区的北斗单频载波相位实测数据进行实验，结果表明，该方法能高效、灵敏、实时地探测出周跳发生的位置。
关键词：矿区变形监测； 周跳探测； 北斗； Kalman滤波
网络出版地址：http://www.cnki.net/kcms/detail/32.1627.TP.20160405.1132.014.html
北斗卫星定位系统是我国自主建设的全球卫星定位系统，到目前为止已经初步具备了向亚太地区提供实时定位、导航和授时服务的能力[1]。周跳探测是北斗载波相位数据预处理中的重要部分，是确保获得“干净”数据的重要保障。目前已有很多学者对周跳探测方面的问题进行了研究，并提出了很多可以实际应用的解决方案。其中比较普遍的是多项式拟合法[2]、高次差法[3]、电离层残差法[4-5]、小波滤波法[6]等，但是不同的方法适用情况不一样，也不够完善。多项式拟合法对于大周跳的探测效果较好，但其探测效果受到多项式类型、拟合次数以及截断误差等因素的影响；高次差法需要逐级做差，周跳探测无法实时进行，而且在探测的过程中放大了噪声，不适用于探测小周跳甚至近百周的周跳值；电离层残差法推导过程比较严密，解决了周跳解的多值性问题，适用于小周跳的探测，且是目前探测效果最好的方法之一，但需要利用多频载波相位观测数据，当2个载波相位观测值同时发生周跳且周跳之比与频率之比接近时，该方法无法探测周跳[7]，同时该方法只针对多频接收机。小波滤波法可以方便地探测到发生周跳历元的突变点，但周跳的具体幅度还需通过其他方式配合探测[8]。
由于多频接收机价格昂贵，所以在实际生产中普遍使用单频GNSS接收机。在应用北斗单频接收机对矿区重点形变区域进行实时变形监测时，为提高北斗单频载波相位观测值预处理中周跳探测的精度和效率，本文提出了一种基于Kalman滤波的单频北斗载波相位实时周跳探测方法，给出了相应的数学模型，论述了周跳与粗差的识别方法，并利用淮北市某矿区沉降变形监测中实测北斗单频载波相位数据对该方法进行了验证。实验结果表明，该方法能准确区分粗差与周跳，同时能较为精确地给出周跳大小，提高了北斗单频载波相位观测数据预处理的精度。
1.1 Kalman滤波模型
Kalman滤波的数学模型由状态方程和观测方程组成，对于线性系统，其离散模型为[9]
(1)
(2)
式中：Xk为tk时刻的状态向量；A为tk-1时刻到tk时刻的状态转移矩阵；Wk-1为tk-1时刻的预报模型噪声，其值取决于状态方程表述的好坏，表述得越精确取值越小[10]；B为tk-1时刻的动态噪声系数矩阵；φk为tk时刻的载波相位观测向量；C为tk时刻的状态向量系数矩阵；Vk为tk时刻的观测噪声，其方差为σ2，在静态载波相位测量中σ=λ/100[11]，即波长的分辨率。
用于周跳探测的Kalman滤波模型是建立在三阶多项式的基础上，记载波相位观测值的采样间隔为T，则
， 
。
状态向量Xk由φk及其对时间t的一阶导
k、二阶导φ″k、三阶导φ‴k组成，即
[12]。
根据最小二乘原理可得Kalman滤波递推步骤：
(1) 状态向量的一步预测值为
(3)
(2) 状态向量一步预测值方差矩阵为
(4)
式中: Q为状态噪声的协方差阵， Q=E(W2)。
(3) 滤波增益矩阵为
(5)
(4) 状态向量估计值为
(6)
(5) 状态向量估计值方差矩阵为
(7)
式中： I为单位矩阵。
1.2 滤波初始值确定
为了启动Kalman滤波程序，必须要有一个初始值 X0和 P0，它一般可由前N 点通过三阶多项式平滑取得，本文取N=12，但N并不是一个固定数值，N的选取应视观测数据质量情况而定，若N值选取不当，容易造成滤波发散，不能很好地探测周跳。
用N个连续的载波相位观测值进行三阶多项式平滑，求出初始值 X0和 P0：
(8)
式中： a(k)=[1 (k-N)T (k-N)2T2/2 (k-N)3T3/6]T； L为载波相位观测值向量。
X0的协方差矩阵为
(9)
在选取N个点进行三阶多项式平滑时，应确保这N个观测值中不含有周跳，这可以通过残差平方和来检验：
(10)
单位权中误差估计值为
(11)
检验时通常以2.5倍的中误差作为限差：
(12)
当式(12)成立时，便可开始进行滤波计算。
1.3 周跳探测具体实现
在进行滤波计算时，先计算得到每个历元的一步预测值 
k/k-1，估计值与实际观测值之间会存在一定的误差 V(k)= φ(k)- 
k/k-1，若其不满足预设限差要求，即
>u(u为预设的周跳检验阀值)，则表明该历元时刻的载波相位观测值含有周跳或粗差，针对此种情况需要做进一步处理。
当>u成立时，需进一步计算第k+1个历元时刻的载波相位观测值与其对应时刻的一步预测值之间的误差 V(k+1)= 
k+1/k-1，并判定 V(k+1)与检验阈值u′=2.5u间的关系，u′=2.5u是因为V(k+1)是第k+1个历元时刻的一步预测误差，在计算过程中第k个历元时刻的观测值并没有参加滤波计算。若
>u′成立，则表明在第k个历元时刻载波相位观测值发生了周跳，此时应重新进行初始化；若不成立，则表明在第k个历元时刻载波相位观测值中含有野值或者粗差，此时应令第k个历元时刻的增益矩阵Jk=0，即该时刻的观测值不参与滤波计算。
具体的周跳探测流程如图1所示。
图1 周跳探测流程
2.1 实验设计
实验数据来源：采集于淮北市某矿区，采样间隔为1 s，经过其他数据预处理软件检验为不含有周跳的“干净”的北斗单频载波相位观测数据。设计2组实验来验证本文所提出的方法在矿区变形监测数据预处理中的可行性以及周跳探测结果的精确度。实验1分别利用多项式拟合法和Kalman滤波法来检验实验数据的噪声水平，结果如图2所示。
(a) 多项式拟合法
(b) Kalman滤波法
图2 实测北斗单频载波相位观测数据无周跳噪声序列
实验2在上述一组干净数据中第100历元处加入大小为1的粗差，在第400，700，1 000，1 300历元处加入大小为2，3，4，5周的周跳后，分别利用多项式拟合法、Kalman滤波法进行周跳探测，结果如图3及表1所示。
2.2 实验结果分析
从图2可以看出，对于采样间隔为1 s的北斗单频实测数据，多项式拟合法和Kalman 滤波法探测出的无周跳噪声水平相差不多，相比之下Kalman滤波的效果要好一些。在实际观测时，受接收机所处环境的影响，使得原始数据含有较大观测噪声，Kalman滤波预报残差序列在[-0.497,0.461]周范围内波动，在一定程度上会影响该方法探测周跳的精度。
图3直观地反映出Kalman滤波可以更加灵敏、准确、实时地探测出周跳或粗差发生的位置，较好地满足实时变形监测的时效性。从表1中可以清晰地看到，本文所提出的适用于北斗单频载波相位周跳的Kalman滤波法能够准确地区分出滤波曲线发生波动时观测值的变动类型。总的来说，在观测噪声较大时，Kalman滤波依然能够灵敏、精确地探测出周跳发生的位置，准确地区分周跳与粗差，同时也能够较为精确地确定发生周跳的估值大小。
(a) 多项式拟合法
(b) Kalman滤波法
图3 北斗单频载波相位观测数据模拟周跳探测效果
表1 北斗单频载波相位Kalman滤波估计周跳值的结果
针对在矿区变形监测中如何提高北斗单频载波相位观测数据预处理的精度问题，提出了相应的Kalman滤波方法，并利用淮北市某矿区的北斗单频载波相位实测数据，从理论分析与实例验证2个方面证明了该方法的可行性及探测结果的可靠性，并得出以下结论：
(1) 提出的Kalman滤波法能较为准确、灵敏地探测和区分出北斗单频载波相位观测值中的粗差或周跳，能近似实时地探测出载波相位观测值中的周跳。
(2) 在启动Kalman滤波探测周跳前，需要N个干净的历元数据进行滤波初始化，同样在探测到发生周跳的位置后仍需要进行初始化，这就要求相邻2个发生周跳位置间的历元间隔不得小于N，因此，这在某种程度上限制了Kalman滤波法探测周跳的准实时性。
(3) 当观测数据本身带有较大噪声时，周跳探测结果与真实值存在一定误差。因此，如何提高该情况下的周跳探测精度，将是下一步的研究重点。
参考文献：
[1] 杨元喜,李金龙,王爱兵,等.北斗区域卫星导航系统基本导航定位性能初步评估[J].中国科学:地球科学,2014,44(1):72-81.
[2] 蔡诗响,张小红,李星星,等.一种基于多项式拟合的单频周跳探测改进方法[J].测绘信息与工程,2009,34(5):1-3.
[3] DAI Zhen. MATLAB software for GPS cycle-slip processing[J].GPS Solution,2012,16(2):267-272.
[4] CAI Changsheng, LIU Zhizhao, XIA Pengfei,et al. Cycle slip detection and repair for undifferenced GPS observations under high ionospheric activity[J].GPS Solution,2013,17(2):247-260.
[5] 王维,王解先,高俊强.GPS周跳探测的方法研究[J].武汉大学学报:信息科学版,2010,35(6):687-690.
[6] 蔡昌盛,高井祥.GPS周跳探测及修复的小波变换法[J].武汉大学学报:信息科学版,2007,32(1):39-42.
[7] BANVILLE S, LANGLEY R B. Mitigating the impact of ionospheric cycle slips in GNSS observations[J].Journal of Geodesy ,2013,87(2):179-193.
[8] 黄兵杰,柳林涛,高光星,等.基于小波变换的GPS精密单点定位中的周跳探测[J].武汉大学学报：信息科学版,2006,31(6):512-515.
[9] GENG Yanrui, WANG Jinling. Adaptive estimation of multiple fading factors in Kalman filter for navigation applications[J].GPS Solution, 2008,12(4):273-279.
[10] 宋迎春,朱建军,陈正阳.动态定位中测量噪声时间相关的Kalman滤波[J].测绘学报，2006,35(4):328-331.
[11] 李猛,文援兰,梁加红,等.基于卡尔曼滤波的周跳探测算法[J].测绘科学技术学报, 2010,27(6):407-411.
[12] 刘伟平,郝金明,汪平,等.Kalman滤波在周跳探测与修复中的应用[J].大地测量与地球动力学,2009,29(6):101-103.
MA Xiaodong1,2， SONG Bing1,2， ZHANG Shubi1,2， LIU Xin1,2
(1.School of Environment Science and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China; 2.Key Laboratory of Resources and Environment Information Engineering of Jiangsu Province, Xuzhou 221116, China)
Abstract：In application of Beidou single frequency receiver in real time mining deformation monitoring, in order to improve precision of cycle slip detection in measurement preconditioning of Beidou single frequency carrier phase, real time cycle slip detection of Beidou single frequency carrier phase based on Kalman filtering was proposed, and recognition method of gross errors and cycle slips was discussed. The experiment was carried out with the Beidou single frequency carrier phase data collected in a certain mining area, the results show that the method can efficiently, sensitively and real-timely detect the location of cycle slip.
Key words:mining deformation monitoring; cycle slip detection; Beidou; Kalman filtering
作者简介：马晓东(1992-)，男，江苏徐州人，硕士研究生，主要研究方向为GNSS数据处理，E-mail：maxdcumt@163.com。
基金项目：国家自然科学基金青年科学基金项目(41504032)；江苏省自然科学基金青年基金项目(BK20150175)。
收稿日期：2015-12-03；修回日期：2016-02-24；责任编辑：胡娴。
中图分类号：TD325.4
文献标志码：A 网络出版时间：2016-04-05 11:32
文章编号：1671-251X(2016)04-0058-04
DOI:10.13272/j.issn.1671-251x.2016.04.014
马晓东，宋冰，张书毕，等.基于Kalman滤波的单频北斗载波相位实时周跳探测[J].工矿自动化，2016,42(4)：58-61.