邵沛涵,王志斌,高成发,张 博,郝文世,高 平
(1.东南大学 交通学院,南京 211189;2.河北雄安京德高速公路有限公司,河北 保定 071799)
北斗卫星导航系统(BeiDou Navigation Satellite System,BDS)是我国自主研发的全球定位系统,结合我国国情形成了“三步走”发展战略[1]。2020年6月23日北斗三号全球卫星导航系统(BDS-3)的最后一颗全球组网卫星发射成功并于2020年7月31日正式宣布开通。对于MEO卫星,BDS-3不仅能够提供兼容BDS-2的频率(B1I和B3I),还增加了B1c和B2a两个信号[2],能够兼容GPS(L5)和Galileo(E5a),以便于实现系统间的兼容互操作。同时,BDS-3相对于BDS-2在卫星星座和信号设计方面有大幅改进并使用新的氢原子钟,多路径效应明显减弱[3]。因此,BDS-3的全面建成必将带来定位精度和稳定性的进一步提高。
与传统监测手段相比,GNSS测量具有高精度、全天候实时监测、工作效率高、人工需求小、测站间无需通视、可提供全球统一的三维地心坐标等特点。王世进等人对GPS/BDS的RTK定位算法进行研究分析[4],结果表明GPS/BDS的RTK定位精度能够达到cm级,且模糊度固定时间相对于单系统明显减少。隋春玲等人分析GPS/GLONASS/BDS组合的单历元单频短基线RTK定位精度[5],在NEU3个方向上的精度均达到mm级;王艺希等人[6]利用附带参数的卡尔曼滤波方法,计算得到31 km基线条件下,NE方向定位误差在1 cm以内,U方向误差在3 cm以内。综上,由于设计规范规定的分层厚度为亚米级,RTK能够满足精度需求,为数字化施工提供有利条件。
文中采用的方案是在压路机顶部安装GNSS接收机、在工程项目部安装GNSS基站(两者相距不超过15 km),使用附带参数的卡尔曼滤波模型进行RTK定位,实时采集地面点坐标,并和现场水准测量结果进行对比,分析RTK定位的精度和层厚计算的准确性。
文中以GPS和BDS组合定位为例,其双差伪距和载波相位观测方程[7-8]为:
(1)
式中:r和s分别代表参考卫星和差分卫星;b和m分别代表基准站和流动站;φ和P表示以m为单位的载波和伪距观测量;ρ表示卫星到接收机的距离;λ和N表示对应频率载波的波长和模糊度;T和I表示对流层和电离层延迟;ε和e表示伪距和相位的观测噪声等其他误差。
将上述载波观测方程写成误差方程并进行线性化,可得:
(2)
其中,
(3)
(4)
伪距误差方程线性化同理,定权时伪距和载波权重按照1∶1 000。GPS和BDS组合定位模型可在各系统内部选择参考卫星进行双差,再组合不同系统的双差观测值来求解基线向量。利用组合模型进行RTK解算可以将各系统的卫星数据都进行应用,并认为系统之间互不相关。
文中采用基于高度角的随机模型,即将观测值噪声σ表达成以卫星高度角E为变量的函数[9]。已有研究表明,高度角定权法比单位权法得到的数据精度有明显提高,且在卫星数量较少时尤为明显[10]。下面简述定权过程:
当考虑到电离层、对流层延迟影响与高度角的比例关系时[11],可以定义一种正弦函数模型。若选择卫星r作为观测历元ti的参考卫星,则卫星s的双差观测值的噪声可表示为:
(5)
式中:Er和Es分别为参考卫星r和卫星s的高度角。
若历元ti中共有n颗卫星,则会构成n-sys个观测方程,sys为使用的卫星系统个数,根据误差传播定律,则方差-协方差阵可表示为:
(6)
将线性化后的误差方程简化为:
V=BX-L.
(7)
由于文中采用的实测数据采样间隔为1.0 s,故对于动态RTK定位宜采用常加速度模型[12],则式中的X即状态向量可以描述为:
X=[xpos,ypos,zpos,xvel,yvel,zvel,xacc,yacc,zacc,
▽ΔN1,…,▽ΔNn-sys]T.
(8)
其中包含位置参数、速度参数、加速度参数和双差模糊度参数。位置参数的初值可以用单点定位的结果,速度参数的初值可以根据多普勒观测值计算得到,当多普勒观测值不足时可根据实际情况设置一个经验值,加速度初始值一般设置为零,模糊度参数初始值可由伪距和载波确定。
t(k)时刻卡尔曼滤波基本递推过程如下[9,13]:
(9)
(10)
(11)
(12)
(13)
系统噪声的方差-协方差阵Qk为:
Qk=
(14)
对于其初值,位置参数可赋予一个较大的方差,速度和加速度对应的方差较小,模糊度对应的方差较大。
状态转移矩阵Φk,k-1的设置如下:
(15)
通过附带参数的卡尔曼滤波模型得到模糊度参数的浮点解,然后使用LAMBDA算法[14]进行固定,将固定解回代到观测方程中即能得到基线向量的固定解。
层厚计算可分为实时处理和后处理,两者基本思路相同,区别是后处理时可以根据实际施工范围和最后的压实时间对所有坐标数据点进行约束,减少循环次数,从而提高计算效率。但考虑到施工现场的实际情况,一般压实后不会马上进行填土,实时处理时最后的压实高度可认为和后处理时选取的时间相同,因此在实际计算时,实时处理和后处理效率基本相同。
1)获取实时坐标数据。通过RTK计算得到当前时刻t(k)的三维坐标数据Pk(xk,yk,zk),并进行储存。
2)从初始时刻t(0)开始以每个历元的平面坐标为基准,建立纵向结构体Di(i=1,2,…,n),该结构体包括分层数据。设置一定的阈值,当其他历元的坐标落在该阈值以内时,认为是同一个纵向结构体。否则,新申请一个纵向结构体Di+1并将该点存入。
3)建立初始分层。判断最先进入纵向结构体的5个历元的高程坐标,剔除粗差后取平均值作为初始层高,建立初始分层结构体。
4)对每个新进入纵向结构体的三维坐标点Pk进行分层判断,如果高程方向差值在给定阈值以内(根据施工规范一般设置为h1=25 cm左右为宜),则将该点添加进这一分层结构体。否则,新申请一个分层结构体并将该点存入。需要注意的是,为避免粗差对分层造成较大影响,应设置高差的最大阈值(h2=35 cm),该值可根据现场填土情况取一个经验值,当超过该阈值时,认为不是填土造成的分层而是粗差,此时不应申请新的分层结构体。
5)重复上述步骤,逐历元判断。当存在两层及以上分层结构体时,即可计算得出层厚。
1)根据先验信息,设置施工区域的坐标范围,根据时间截取需要处理的数据,从而减少数据量提高后处理速度。
2)划分格网。地形起伏较大时格网边长可设置小一些。
3)循环所有数据点,并判断到格网中心点的距离,在阈值以内时将其加入到该格网结构体中。
4)根据时间最近原则生成所需数据,根据每个数据点到格网中心的距离取加权平均值作为格网的高程,具体算法如下:
(16)
5)生成三维地形数据。三维地形数据实时处理流程如图1所示。
图1 实时处理流程
测试时,在路基上事先布设若干点(见图2),相邻点间隔10 m左右,利用RTK测出点的平面坐标以便于后续放样。从施工现场的高等级水准点出发,采用二等水准的方式将高程引入路基,并测量出布设点的高程,每次压实后测量出真实高程,通过两次高程相减得出层厚,并和利用压路机实时数据计算的结果进行对比,以验证文中所提方法的准确性。
图2 RTK点位布设示意图
压路机采集数据坐标系为WGS84坐标,水准测得的高程为施工坐标系,可通过Bursa七参数模型实现两者的统一[15]。具体做法是首先采集3个及以上已知施工坐标的控制点的WGS84坐标,通过式(17)计算出七个参数,再利用七参数将压路机采集的数据按式(18)实时进行坐标转化。
(17)
(18)
文中使用的数据为“京德高速”现场实测数据,数据采集日期分别为2020年8月1日(年积日214)和2020年8月4日(年积日217),这两天分别进行两次填土和压实,如图3、图4所示。基站距离压路机施工现场4.5 km左右。
图3 施工现场示意图
图4 压路机工作过程示意图
3.2.1 BDS-3定位性能分析
文中在分析BDS-3的定位性能时,将原始观测数据导出,分析GPS、BDS-2和BDS-3不同系统组合情况下,RTK定位时的模糊度固定成功率和U方向上的定位误差(层厚计算对U方向的精度尤其重视)。图5、图6给出每个历元下各系统卫星数量和对应的PODP值,其中PDOP数值越小卫星分布程度越好,定位精度越高。
图5显示每个历元下BDS卫星数量已经远超过GPS卫星数量,同时BDS-3卫星数量基本在4颗以上,可以满足单BDS-3定位的需要,从图6可以看出,只有GPS卫星或BDS-2卫星的情况下卫星的几何分布不理想,同时也能看出即使BDS-2卫星数量比GPS卫星多的情况下有时PDOP值仍较大。但加入BDS-3卫星后整体PDOP值显著减小,可见BDS-3空间结构优于BDS-2,同时GPS+BDS-2+BDS-3多系统的融合进一步增强卫星的几何结构,使得PDOP值再次减小。
图5 各系统卫星数
图6 各系统PDOP统计
表1给出压路机在施工区域工作时,在GPS、BDS-2和BDS-3不同系统组合情况下基线模糊度固定成功的历元数和成功率,可以看出加入双系统固定成功率比单系统成功率要高,且加入BDS-3后模糊度固定成功率得到了极大的改善。
表1 模糊度固定情况比较
由于层厚计算最关键的是U(高程)方向的定位误差,在精度要求范围内,可认为二等水准测量结果为真值,因此文中比较GPS、BDS-2和BDS-3不同系统组合情况下计算得到的高程和二等水准测量结果的差值,如表2所示。从表中可以看出,使用BDS定位时,高程方向误差为cm级,能够满足层厚计算的精度要求。同时使用BDS-3定位相对于不使用在高程方向上的精度得到了明显的提升。
表2 定位结果比较
3.2.2 层厚计算
由上述分析可知,将GNSS-RTK技术应用于路基层厚计算中是可行的,通过压路机采集的离散数据,可将施工现场划分出若干格网(格网效果如图7所示),再通过内插和加权平均的方式,就可以计算得到施工现场任意点的高程坐标,两次压实后的高程相减就能得到任意点的层厚值。
图7 格网效果(m)
文中分析2020年8月1日(年积日214)和2020年8月4日(年积日217)两次填土后形成的新层的层厚,并和二等水准结果做比较(部分点如表3所示,整体结果如图8所示),结果显示,GNSS计算的层厚值(平均值40.2 cm)和水准测量的层厚值(平均值40.3 cm)整体偏差在1 cm以内,因此文中提出的方法能够准确地计算出层厚值。
图8 层厚整体情况
表3 层厚值比较 m
文中首先给出GNSS-RTK定位的函数模型和层厚计算的具体过程,分析比较GPS、BDS-2、BDS-3不同组合定位的精度,结果表明多系统定位在卫星空间结构和定位精度上均优于单系统,BDS-3相对于BDS-2的定位精度提升幅度较为明显。在层厚计算最关注的高程方向精度,多系统定位能够满足实时计算的精度要求。
文中讨论内容仅限于将GNSS应用于层厚计算,对于路基施工更关注的压实度暂未涉及,有关GNSS+传感器融合处理的内容仍需继续研究。此外,对于BDS-3多频,尤其是新增的几个频点的性能评估仍需持续分析。