BDS/GPS实时变形监测技术的软件实现

周明志，胡志刚，李鹏博

(武汉大学卫星导航定位技术研究中心，湖北 武汉 430079)



BDS/GPS实时变形监测技术的软件实现

周明志，胡志刚，李鹏博

(武汉大学卫星导航定位技术研究中心，湖北 武汉 430079)

Realization of BDS/GNSS Real-time Deformation Monitoring Software

ZHOU Mingzhi，HU Zhigang，LI Pengbo

摘要：阐述了BDS/GPS变形监测实时算法，并利用C#和C语言研发了BDS/GPS实时变形监测软件。该软件支持多种格式的数据采集、存储、解算和图形显示。精密导轨试验结果表明，5 min的BDS/GPS实时解在E、N、U 3个方向的重复性精度分别为1.0～1.3 mm、1.4～1.7 mm和2.7～3.5 mm，探测精度在平面和高程方向分别为0.2～0.4 mm和1.1 mm。

关键词：BDS；实时监测；Kalman滤波

以GPS为代表的卫星导航定位技术由于具有高精度、全天候、准实时的优点，在桥梁、尾矿、大坝、高层建筑等方面得到了广泛应用[1-4]。1998年清江隔河岩大坝建立了我国第一个GPS变形监测系统[5]，并在当年的抗洪抢险决策中发挥了巨大作用。我国北斗卫星导航系统于2012年12月27日完成区域阶段部署，目前有5颗GEO、5颗IGSO和4颗MEO卫星在轨运行，可为亚太大部分地区提供公开服务[6]。近年来围绕北斗系统的高精度导航定位应用研究正在逐步开展。已有研究表明，BDS/GPS联合定位精度比GPS单系统的精度提高20%～30%，尤其在动态定位模式其精度提升明显[7]。在GNSS变形监测方面，现代化的超高层建筑、大型桥梁等不但要求毫米级的定位精度，还需要尽量缩短观测周期来监测其受大风等外部策动力影响时建(构)筑物的摆动幅度和振动周期信息。这就对GNSS变形监测的精度、稳定性和实时性提出了较高的要求。

大坝、水库、滑坡体变形具有连续、缓慢、量级小的特点，常规监测中一般利用一段时间数据(数小时)作为一组，采用静态相对定位方式处理，但此方法比较滞后，实时性不强。为了提高数据处理的实时性，常采用动态单历元处理模式(OTF)[8]。OTF需要初始化几十分钟并且保证模糊度连续跟踪。李博峰、韩保民等从提高法方程病态性角度求得单历元稳定解[9]。余学祥从观测方程入手，提出直接求得变形量的方法[10]。上述两种方法虽然避免了载波周跳探测，但都需要外界提供相对准确的未知参数先验信息。

本文首先探讨了BDS/GPS多模数据融合应用于高精度GNSS形变监测的关键技术，包括BDS/GPS时空基准统一、实时探测周跳、相位中心改正、Kalman滤波参数估计；然后运用C#和C语言开发了BDS/GPS实时变形监测软件；最后基于精密导轨试验平台采用人工增加变形量的方法验证了算法的正确性和软件的解算精度。

一、BDS/GPS多模数据融合关键算法

1. BDS/GPS时空基准统一

由于广播星历拥有特定的坐标基准和时间基准，进行不同系统组合定位时，往往需要考虑系统间的时间基准和坐标基准差异。GPST和BDT都采用原子时作为时间基准，但是其起算点不同：BDT起始历元分别为UTC时2006年1月1日0时0分0秒，而GPST为UTC时1980年1月6日0时0分0秒。由于UTC闰秒的存在，BDT与GPST之间有14 s的差异。除此之外，两者之间还存在微小的同步误差τ，在相对定位中可以忽略不计。两者之间的转换关系为

GPST=BDT+14 s+τ

在多模接收机采集的数据中，需要考虑上式以换算正确的卫星发射时刻来计算正确的卫星位置。

另外，BDS广播星历采用CGCS2000框架，GPS广播星历采用WGS-84框架，两者之间存在一定差异。文献[11]指出，在坐标系的实现精度范围内，无须考虑坐标系差异引起的结果偏差。实际上，由于大坝等构筑物变形监测局限于几十千米范围以内，坐标系之间的差异对基线解的影响不超过1 mm，故本文未考虑系统间的坐标转换。

时空基准统一后的BDS/GPS双差数学模型为

(1)

线性化后函数模型为

(2)

BDS与GPS是两个独立系统，观测精度和空间分布不同导致对待估参数贡献不同，在组合定位中需要进行非等权处理。如表1所示，Odolinski通过大量的数据统计出BDS和GPS非差伪距噪声σp和非差载波噪声σΦ经验值[13]。因此，BDS/GPS组合定位随机模型为

(3)

表1　非差伪距载波噪声

2. 接收机天线相位中心改正

接收机天线相位中心一般表述为相对于天线参考点的天线相位中心偏差(antenna phase center offser，PCO)和相对于PCO的天线相位中心变化(phase center variations，PCVs)。研究表明，不同的频率其天线相位中心存在几毫米甚至几厘米的差异。不考虑相位中心改正项将会造成毫米级甚至厘米级的位置偏差[14-15]，特别是高程方向。然而，目前还没有官方发布北斗接收机天线相位中心精确的改正模型。在这种情况下建议所有参与变形监测的接收机天线都进行指北，从而在进行相对定位时可以消除天线相位中心变化(PCV)对基线解的影响，但PCO仍然无法消除。此时得到的基线是相对于北斗PCO和GPS PCO相应的频率组合PCO值而不是监测站标石。对于关心监测站相对位置的变形监测来说影响不大，但对于精密单点定位来说，如果要求计算测站标位置时，接收机天线相位中心改正是强制修改项。

3. 实时周跳处理

周跳的产生中断了高精度载波观测值跟踪弧段的连续性。1周的周跳将对载波产生约20 cm的偏差。可靠的周跳探测算法，尤其是小周跳探测在目前仍然具有一定的挑战性。准确可靠和能够实时周跳探测是高精度实时变形监测解决模糊度的关键之一。目前学者提出了许多方法，如高次差法、小波分析法等，但都是基于事后数据，不能应用于实时[16]。软件采用刘志赵提出的TECR(total electron contents rate)方法[17]：采用电离层变化率并结合宽巷组合可实现实时周跳探测。宽巷Nwl表达式为

(4)

(5)

进一步推导k历元TECR为

(6)

将电离层变化率推导k历元电离层总量(TEC)变化值与载波观测值推导出的电离层变化值作比较，可得

(7)

由式(5)—式(7)推导TECR探测周跳噪声为

(8)

(9)

4.Kalman滤波参数估计

卡尔曼滤波算法是利用观测向量来估计随时间不断变化的随机向量，是一种对动态监测数据进行处理的有效方法，具有提高数据处理效率、减少噪声、提高预测精度等优点，在变形监测中常用来实时估计监测站点位。其公式如下。

状态更新

(10)

测量更新

(11)

(12)

初始方差阵P0为对角矩阵，其中坐标参数方差由伪距单点定位求得，单差模糊度方差由式(12)通过误差传播定律确定。当Pk阵设置不正确将会导致滤波发散，对于上一步探测出有周跳的卫星，先验Pk阵不能反映其模糊度精度，需要对其相应的单差模糊度参数方差重新初始化。

(13)

再根据MLAMBDA[18]方法求双差模糊度整数解N，根据模糊度残差平和的次小值与最小值的比率ratio判断模糊度固定是否正确[19]。检验通过后固定解为

(14)

值得注意的是，当有新卫星进入视野时，ratio值偏低，致使检验无法通过。这是因为新参与解算的模糊度精度较低，若和其他模糊度一起在整数最小二乘估计下(IIS)估计会影响其他模糊度固定。Verhagen提出部分模糊度固定方法[20]：从精度最高的模糊度开始，依次寻找子集模糊度，使它们理论固定成功率不低于给定的数值p0，通常是99%。只固定这部分模糊度，其余模糊度保持不变。这种方法需要反复计算模糊度成功率，直至找到最佳子模糊度组合。此处参考Verhagen提出的模糊度精度越高，整数解越可靠的思想，设定只有卫星高度角大于25°其模糊才参与固定，25°以下模糊度只参与实数解。经过15°到25°的滤波收敛，其精度不断提高，理论成功固定率也提高，这时再参与LAMBDA解算[21]。

二、软件实现

软件采用客户端/服务端(C/S)模式编写，客户端和服务端通过TCP/IP协议进行数据交换。客户端负责结果图表显示和用户命令交换，采用C#编写。第1节所阐述的关键算法，以及测站管理、数据解码等核心功能通过服务端完成，采用C语言编写。服务端在后台连续运行，只要输入正确的账号和密码，客户端就可以获得服务端发送过来的结果数据，并进行图表显示，软件界面如图1所示。软件流程如图2所示，主要包括测站管理、数据解码、BDS/GPS实时定位算法、结果显示等功能。其中，测站管理是对测站进行添加删除等操作，数据解码是对实时二进制数据解码。虽然RTCM 3.x修订版开始支持BDS数据，message type为1121—1127。但现在市场上接收机大部分还是RTCM 2.0版本或不是最新版本的RTCM 3.0，并不包含BDS数据。监测中采用BINEX格式，它是瑞士unavco机构协调众多接收机厂商制定的标准二进制数据格式，支持BDS/GNSS实时数据编码，北斗观测值和导航星历Record ID分别是0x01-05和0x7f-05，可对解算的结果进行实时显示，包括星空图、信噪比、监测站位移等。

图1　软件界面

图2　软件流程

三、精密导轨试验验证

为了验证软件的实时监测精度，在国家光电中心楼顶对本文研发的变形监测软件进行了测试。测试楼顶有2个观测墩，都安置TRM59800天线，其中一个监测墩安置导轨，导轨每次移动大小可以通过千分尺准确量取，如图3所示，这样可以通过移动导轨来模拟变形量。观测时段模拟变形量见表2。在测试过程中，设Kalman滤波时间为5 min。当时刻为整5 min时，重新初始化滤波器；最后统计5 min滤波结果精度。这样既可以保证精度又可以保证探测灵敏度。

图3　试验环境

序号观测时段状态时段12015-04-0916:40:57—2015-04-1015:25:00 初始状态时段22015-04-1015:30:00—2015-04-1116:00:20 E方向移动10mm时段32015-04-1116:40:10—2015-04-1217:25:00 W方向移动10mmU方向挪动30mm时段42015-04-1217:25:25—2015-04-1317:20:25 E方向移动10mm时段52015-04-1317:20:25—2015-04-1415:20:20 W方向移动5mm

系统连接接收机IP后，自动对流动站和监测站传过来的BINEX格式BDS/GPS实时数据解码、同步、组基线、Kalman滤波、LAMBDA解算、图形显示并将结果保存到数据库中。从数据库中提取滤波5 min结果，其时间序列图和统计表分别见图4和表3。

图4　模拟形变E/N/U三维方向时间序列图

时段平均值/m均方根RMS/mmENUENU解个数1-3.925-0.479-0.2061.01.52.82752-3.915-0.479-0.2061.11.52.92863-3.925-0.479-0.2371.21.73.52964-3.915-0.478-0.2381.31.42.72885-3.920-0.479-0.2381.01.52.7233

注：每个时段起始截止观测时间不一致，网络传输可能存在丢包情况，因而每个时段解的个数不一样。

表4　各时段模拟形变探测精度　mm

从图4可以看出，当没有发生变形时，实时解稳定，当导轨移动，实时解有一个突变，可以很明显看出一个形变过程。从表3中可以看出，BDS/GPS实时解3个方向解的稳定性分别在1.0～1.3 mm、1.4～1.7 mm、2.7～3.5 mm之间，E方向解的稳定性高于N方向，这是BDS卫星几何结构决定的。从表4可以得到软件实时探测精度在平面和高程方向分别是0.2～0.4 mm、1.1 mm。

四、结束语

高精度实时变形监测软件对实时性要求很高，传统的最小二乘往往是基于数小时的数据，结果滞后严重，不能应用于实时。单历元解不稳定，需要额外提供相对准确的先验信息。Kalman滤波能够基于状态域和观测域进行状态参数实时更新，非常适合应用于实时变形监测。本文介绍了基于Kalman滤波的BDS/GPS定位的关键算法，并在实现中用单差模糊度替代双差模糊度，用部分模糊度替代全部模糊度进行固定。基于该算法，采用C和C#语言，笔者开发了BDS/GPS实时变形监测软件。通过精密导轨试验证明该软件BDS/GPS实时解稳定性为1.0～1.3 mm、1.4～1.7 mm和2.7～3.5 mm之间，探测精度在平面和高程方向分别为0.2～0.4 mm和1.1 mm，满足变形监测的需求。在高精度形变监测中，多路径是影响定位结果一个重要因素，在后续工作，需对其影响作进一步研究。

参考文献：

[1]刘超. 矿区 GPS 变形监测及其伪卫星增强技术 [D]. 徐州：中国矿业大学, 2011.

[2]李伟. GNSS 算法研究及其在变形监测中的应用 [D]. 成都：西南交通大学, 2013.

[3]任常. 超高层建筑 GPS 变形监测单历元算法研究[D]. 北京：清华大学, 2013.

[4]丁培江. 水库潜在滑坡区 GNSS 变形监测管理信息系统研究与开发[D]. 昆明：昆明理工大学, 2014.

[5]姜卫平, 刘经南. GPS技术在隔河岩大坝监测中的应用研究[J]. 武汉测绘科技大学学报, 1998,23(S0):20-22.

[6]Chinese Satellite Navigation System Management Office. Beidou Satellite Navigation System Development Report(2.1)[R]. Beijing：[s.n.], 2012.

[7]SHI C, ZHAO Q, HU Z, et al. Precise Relative Positioning Using Real Tracking Data from COMPASS GEO and IGSO Satellites[J]. GPS solutions, 2013, 17(1): 103-119.

[8]陈永奇,JAMES L. 单历元GPS变形监测数据处理方法的研究[J]. 武汉测绘科技大学学报, 1998，23(4): 45-49,84.

[9]李博峰, 沈云中. 附有约束条件的GPS模糊度快速解算[J]. 武汉大学学报(信息科学版), 2009,34(1): 117-121.

[10]余学祥, 徐绍铨, 吕伟才. GPS变形监测信息的单历元解算方法研究[J]. 测绘学报, 2002,31(2):123-128.

[11]魏子卿. 2000中国大地坐标系及其与WGS84的比较[J].大地测量与地球动力学, 2008, 28(5):1-5.

[12]卢祥. GPS/BDS短基线实时快速相对定位研究[D]. 北京：中国科学院大学，2014.

[13]ODOLINSKI R, TEUNISSEN P J G, ODIJK D. An Analysis of Combined COMPASS/BeiDou-2 and GPS Single-and Multiple-frequency RTK Positioning[C]∥Proceedings of the ION 2013 Pacific PNT Meeting. Honolulu, Hawaii：[s.n.], 2013: 69-90.

[14]ROTHACHER M, SCHAER S, MERVART L, et al. Determination of Antenna Phase Center Variations Using GPS Data[C]∥Proc.. IGS Workshop on Special Topics and New Directions. [S.l.]: IGS,1995: 205-220.

[15]MADER G L. GPS Antenna Calibration at the National Geodetic Survey[J]. GPS Solutions, 1999, 3(1): 50-58.

[16]方荣新, 施闯, 魏娜,等. GPS数据质量控制中实时周跳探测研究[J]. 武汉大学学报(信息科学版), 2009,34(9): 1094-1097.

[17]LIU Z. A New Automated Cycle Slip Detection and Repair Method for a Single Dual-frequency GPS Receiver[J]. Journal of Geodesy, 2011, 85(3): 171-183.

[18]CHANG X W, YANG X, ZHOU T. MLAMBDA: a Modified LAMBDA Method for Integer Least-squares Estimation[J]. Journal of Geodesy, 2005, 79(9): 552-565.

[19]TEUNISSEN P J G, VERHAGEN S. The GNSS Ambiguity Ratio-test Revisited: A Better Way of Using It[J]. Survey Review, 2009, 41(312): 138-151.

[20]VERHAGEN A A, TEUNISSEN P J G, VAN DER MAREL H, et al. GNSS Ambiguity Resolution: Which Subset to Fix?[C]∥IGNSS Symposium 2011, International Global Navigation Satellite Systems Society. Sydney, Australia: University of New South Wales, 2011.

[21]TAKASU T, YASUDA A. Kalman-filter-based Integer Ambiguity Resolution Strategy for Long-baseline RTK with Ionosphere and Troposphere Estimation[C]∥ION NTM. [S.l.]: ION, 2010: 161-171.

中图分类号：P228.4

文献标识码：B

文章编号：0494-0911(2016)02-0026-06

作者简介：周明志(1990—), 男，硕士，主要从事形变监测研究。E-mail: nanwang_luojia@foxmail.com

基金项目：国家863计划(2013AA122501)；国家自然科学基金(2014AA123101)；中央高校基本科研业务费专项资金(2042014kf0085)

收稿日期：2015-06-30

引文格式： 周明志，胡志刚，李鹏博. BDS/GPS实时变形监测技术的软件实现[J].测绘通报，2016(2)：26-31.DOI:10.13474/j.cnki.11-2246.2016.0042.