一种伪距单点定位的数学模型研究及程序实现

郭平，周适，段太生，李学仕，王靠省

(1.中铁二局集团有限公司,四川 成都 610031,2.中铁二院工程集团有限责任公司，四川 成都 610031)

0 引 言

单点定位计算测站坐标已是比较成熟的技术．很多商业软件着重解决不同测站点的基线向量解算,采用伪距和载波相位观测值参与基线解算,笔者查阅文献资料,现存大量关于单点定位甚至精密单点定位(PPP)的论文,但其讲述的数学模型不够详细,通过论文依然较难整理出能直接应用于程序设计的一整套数学公式．笔者对单点定位的数学模型进行研究,提出一种能正确计算单点定位的数学模型,利用伪距观测值,通过C#编程实现计算测站单点定位坐标,并提出一些结论．

1 单点定位数学模型分析

单点定位可利用接收机实测获得原始数据,并转换为通用的RINEX格式．其中N文件为GPS卫星的广播星历文件,G文件为GLONASS卫星的广播星历文件,C文件为BDS卫星的广播星历文件,N、G、C文件也可整合形成为一个文件,即公共导航星历P文件[1],记载了各卫星系统的星历数据．通过星历文件可以计算卫星某一时刻的空间直角坐标XYZ, 再通过O文件各历元的伪距和相位观测值可以计算测站点的坐标．本文主要研究利用O文件的伪距观测值,在GPS卫星系统下,采用无电离层模型,计算测站单点定位的坐标．首先对利用一个历元的观测数据实现测站单点定位功能的数学模型进行详细介绍,在此基础上实现多历元数据共同解算测站单点定位的坐标．

若在历元时刻t,接收机接收到n颗GPS卫星,选取高度角最高的卫星作为参考卫星,记作卫星r,其余卫星记作卫星s．

对于参考卫星r,列出该卫星与地面测站点在L1频率的伪距观测方程[2],如式(1)所示．

(1)

(2)

式中：x0，y0，z0为测站初始坐标,初值可设为(0,0,0),或采用该历元所有GPS卫星的重心坐标在地面的投影值,初值不同对最终测站点的坐标计算没有任何影响,只是迭代次数和运算时间上的差异而已．xr，yr，zr为卫星r的空间直角坐标．

同理,可列出其余卫星s与地面测站点在L1频率下的伪距观测方程

(3)

需注意的是,1个历元共接收n颗GPS卫星,除了参考卫星r外,其余卫星数量为n-1,式(3)可列出n-1个方程．

利用其余卫星和参考卫星进行星间作差,即用式(3)减去式(1),可消去式(1)和式(3)的公共项C·dtr,即接收机钟差可消去,可得:

(4)

(5)

602)C·dTrs+(772-602)Trs+

(6)

将式(6)移项可得:

C·dTrs-(772-602)Trs

(7)

简写为

L=Bx+V.

(8)

对于1个历元接收到n颗GPS卫星,可列出n-1个误差方程,如式(8)所示．其中常数项矩阵L(n-1)×1,系数矩阵B(n-1)×3,未知数矩阵x3×1,误差项矩阵V(n-1)×1．根据最小二乘原理,可计算未知数x,如式(9)所示[3]．

x=(BTPB)-1BTPL.

(9)

此时,测站坐标可按式(10)计算．

(10)

由于测站初始坐标为(0,0,0),因此一次计算的结果还不准确,需进行迭代计算．迭代结束计算的判别指标可令式(9)计算的改正数x<1 cm．即改正数小于1 cm时程序停止迭代循环计算,笔者采用编制的单点定位程序运行时发现,迭代计算一般在5次以内,即可得到改正数x<1 cm的结果,大部分历元数据迭代次数在3～4次即可完成计算．

式(9)中的权阵P可通过方差阵D求逆运算得到．考虑到历元的每颗卫星高度角不同,高度角越高的卫星其方差应越小[4],考虑到式(7)随机误差所带的系数,将卫星不同频率伪距观测值的方差视作相同,采用式(11)确定其方差．

(11)

根据误差传播定律,不难推导方差阵D：

D=(774+604)

(12)

式中,r为高度角最高的参考卫星,其余n-1颗卫星的方差根据不同高度角求得．

通过以上分析可求得每个历元计算的单点定位坐标及精度信息．综合各历元数据计算的测站单点定位的坐标可通过法方程叠加的方法实现．具体分析如下:首先计算得到的每个历元的单点定位坐标,以前面两个历元的数据为例,未知数x可写为

(13)

现需要用前面两个历元的数据共同计算坐标,其误差方程式为

(14)

根据最小二乘原理,未知数矩阵x为

(15)

展开后可得式(16)

(16)

同理,若有k个历元共同参与解算,将k-1个历元计算的坐标作为初始值,通过第k-1个历元计算得到的坐标和累加的(BTPB)-1和BTPL,计算以k个历元数据共同计算的单点定位坐标:

(17)

采用k个历元数据共同计算测站坐标的单位权中误差σ0，

(18)

以k个历元共同计算测站点坐标的点位精度以及表征卫星分布状态的DOP值计算过程如下:首先计算空间直角坐标系下的协因素阵Qxx，可由式(19)[5]得出，

(19)

该协因素阵是在空间直角坐标系中给出的,为了便于估算测站的位置精度,常采用其在站心坐标系的表达形式．站心坐标系统中的协因素阵设为QB．

(20)

其中旋转矩阵H为[6]

(21)

式中,B、L分别为测站点对应的经纬度．

测站N、E、U方向上的点位精度MN、ME、MU可由式(22)计算得出．

(22)

平面位置精度因子(HDOP)值、高程精度因子(VDOP)值、空间位置精度因子(PDOP)值可由式(23)计算得出．

(23)

由于单独每个历元的数据都可解算一个测站点坐标．各历元数据可视作是独立观测值．可通过均方根误差(RMSE)指标判别各历元解算坐标与平均值的变化幅度．可通过式(24)求得坐标各分量的RMSE值和坐标的RMSE值．

(24)

RMSE的定义应是各历元坐标分量减去坐标真值,而非减去平均值(真值可采用精密星历计算的测站坐标,精度可达cm级,相对于广播星历m级的精度,可近似认为是真值)．若没有真值的情况下,可简化用平均值代替．

在计算坐标各分量RMSE值后,易于得到一种简易判别粗差的数据处理方法．即通过计算得到的RMSE值,判别各历元的坐标与平均值的吻合程度,若坐标分量较差超过2倍RMSE值,则不采用该历元的数据进行坐标计算的结果．

2 计算过程中需要用到的坐标转换模型

在进行单点定位计算时,需要用到的坐标转换有四个坐标转换模型．分别是:1)大地坐标BLH转换为XYZ;2)空间直角坐标XYZ转换为大地坐标BLH;3)空间直角坐标XYZ转换为站心坐标NEU;4)站心坐标NEU转换为空间直角坐标XYZ．

由BLH转换为XYZ的数学公式为

(25)

由XYZ转换成BLH的数学公式为[7]

(26)

式中,纬度B需要迭代计算[8]，可令初值

由某点的空间直角坐标XYZ和以测站点X0、Y0、Z0作为站心原点的站心坐标NEU的数学公式为

(27)

将式(27)方程两边乘以H-1,易求由站心坐标NEU转换为空间直角坐标XYZ的公式：

(28)

3 根据广播星历计算GPS卫星坐标的程序设计

由GPS广播星历可计算卫星在WGS-84系统下的空间直角坐标,星历文件每2小时发布一次,每颗卫星都有参考历元瞬间的开普勒轨道6参数、反映摄动力影响的9参数以及参考时刻参数,共计16个星历参数[9]．广播星历计算卫星的精度约为2 m,若需要卫星坐标达到厘米级精度,需下载精密星历产品,本文采用GPS导航N文件(广播星历)计算卫星坐标．广播星历的数据格式定义及卫星坐标计算公式可参见相关文献[2，6]．

在程序设计时,需建立各种类,采用面向对象的方法进行编程．首先应建立RINEX类用于读取N文件和O文件,RINEX类中定义两种方法,用于读取N文件和O文件,并将读取的数据进行存储．为使用方便,还需要定义卫星类(Satellite),星历类(Ephemeris),空间笛卡尔坐标类(Cartesian)．在星历类(Ephemeris)中定义N文件中所有星历参数的变量,在笛卡尔坐标类(Cartesian)中定义空间直角坐标XYZ的变量, 卫星类(Satellite)包含星历、接收机钟差改正等信息,在卫星类中编写一个Cartesian类型的函数,通过给定的时间参数,计算得到Cartesian类型的卫星空间直角坐标XYZ．

需注意的是,一天24小时共发布12个星历,可用一个变量表示星历号,程序通过给定的时间判断星历号,若O文件的历元时刻所对应的星历号在N文件中没有,则需要搜素与之相邻的星历,这样才能获取N文件中的星历参数用于计算GPS卫星坐标．

4 根据O、N文件计算测站单点定位的程序设计

通过对N文件的读取,获取各GPS卫星的星历参数,可计算任一时刻GPS卫星的空间直角坐标XYZ．通过对O文件的读取,获取各历元在L1、L2频率下的伪距观测值P1、P2[10].

在进行程序设计时,应建立如下类:数据类(data)、 椭球类(Ellipsoid)、大地坐标类(Geodetic)、站心坐标类(Local-Coordinates)、坐标转换类(Coordinate-Transform)、对流层改正类(Troposphere)、测站类(Station)．数据类(data)里定义如下变量:伪距观测值、GPS卫星的PRN号(1～32),GPS卫星在每个历元数据中的索引号．椭球类Ellipsoid定义椭球的参数变量,如椭球长半轴、第一和第二偏心率、扁率．Geodetic定义某点的大地坐标BLH．Local-Coordinates定义某点相对于测站点为原点的站心坐标NEU．坐标转换类Coordinate-Transform中定义四种静态类型的坐标转换的方法,如前文所述．Troposphere定义对流层改正模型的方法,本文采用Saastamonien对流层改正模型进行对流层改正．测站类Station中实例化之前定义的变量类型,以及定义式(9)中的各矩阵数组,通过之前所述的单点定位数学模型,实现利用伪距观测值进行单点定位计算,可计算测站接收机相位中心的三维坐标以及得到DOP值、点位精度等精度信息．最后根据接收机天线相位参考点ARP和天线垂高信息,可计算测站点的空间直角坐标X，Y，Z．程序按照如下流程图进行设计,可实现单历元单点定位计算, 在单历元计算得到定位坐标后,根据前述数学模型,可计算多历元数据共同作用下单点定位的三维坐标．

图1 根据O、N文件计算测站单点定位程序流程图

5 实测数据的坐标及精度分析

5.1 采用各历元单独计算的坐标和RMSE值分析

运用本文所述的数学模型编制单点定位程序,利用实测某静态RINEX数据,采样间隔10 s,共396个历元,高度截止角为15°．程序自动计算采用从1～396个历元每个历元单独计算的单点定位坐标．如图2～4所示．

图2 单独每个历元计算的X坐标

图3 单独每个历元计算的Y坐标

图4 单独每个历元计算的Z坐标

从图1～3所示,单独利用每个历元的数据计算的单点定位坐标是离散不连续的,没有规律可循,这是因为各历元之间的数据本身是独立的,各历元的坐标分量变化幅度在16 m以内,如表1所示．

表1 单独采用每个历元数据计算的坐标统计表

由表1可见,X坐标变化范围为10.81 m,Y坐标变化范围为15.56 m,Z坐标变化范围为11.47 m．各历元之间坐标分量变化范围在16 m以内,震荡稍大．这是因为仅采用了GPS伪距观测值,未使用相位观测值数据．若数学模型能加上相位观测值,则各历元单独计算的坐标值震荡会更小．

将各历元单独计算的坐标取平均值后与O文件中的近似坐标比较,X坐标较差为0.84 m,Y坐标较差为-1.78 m,Z坐标较差为0.07 m,各坐标分量较差均在2 m范围内,一方面,O文件中的近似坐标也并非准确坐标,是接收机自带程序计算的单点定位坐标,算法未知,可能使用了GPS和GLONASS的伪距观测值或者相位观测值．而本文仅使用了GPS卫星系统的伪距观测值．

5.2 采用多历元共同计算的坐标和精度分析

笔者根据前文所述的数学模型,采用法方程叠加的原理编制单点定位程序可计算多历元数据共同解算的坐标．从第2个历元开始,使用1～2个历元的数据共同解算,第3个历元则采用1～3个历元共同解算,直到第396个历元完成1～396个历元的数据共同解算．图5～7为多历元数据共同解算的XYZ坐标变化图．

图5 多历元共同解算的X坐标

图6 多历元共同解算的Y坐标

图7 多历元共同解算的Z坐标

由图4～6可知,采用多历元数据共同解算的单点定位坐标XYZ,其坐标值随着参与解算的历元数据增加,相邻多历元的坐标变化可认为是连续而非离散的,这是由式(17)中的法方程叠加决定的．相邻的多历元之间(如1～25和1～26)解算的坐标是很接近的,坐标变化不到0.5 m,但综合不同历元解算后的坐标,坐标震荡变化比较显著．这主要有以下三个原因:1)个别历元解算的精度较低影响了整体解算的精度,可通过粗差判别予以剔除;2)数学模型和个别未知参数估计有待进一步地优化和研究；3)本次解算只采用了GPS单系统的伪距观测值,未加入相位观测值及其它卫星系统的观测值,导致坐标变化范围较大．若采用多卫星系统,加入各卫星系统的相位观测值,采用多历元数据计算的坐标值变化范围将大大缩小．

在不考虑粗差影响的情况下,采用所有历元数据(1～396)解算的坐标是最终结果．其值与O文件中的近似坐标较差比较如表2所示．

表2 所有历元共同解算的坐标与O文件中近似坐标较差比较

由表2可知,所有历元数据共同解算的坐标与O文件近似坐标各分量较差在3 m以内,一方面O文件的近似坐标也不是准确的数值,只能作为一个参考．如果用所有历元共同解算的坐标与前文中单独各历元计算的坐标平均值比较,如表3所示,坐标各分量较差在3.5 m以内．

表3 共同解算与单独各历元解算的坐标平均值比较

在计算坐标的同时,点位精度Mx、My、Mz和DOP值也一并计算得到．如图8～9所示．

图8 多历元共同解算的点位精度Mx，My，Mz

图9 多历元共同解算的HDOP，VDOP，PDOP

采用所有历元数据(1～396)解算的点位精度和DOP值如表4所示．

表4 所有历元(1～396)解算的坐标和点位精度、DOP值m

由图8～9的点位精度和DOP值变化图可以看出,随着参与解算的历元数据增加,点位精度Mx，My，Mz的数值在降低,DOP值也降低,这符合法方程叠加的特点,参与解算的历元数据越多,精度越高．采用所有历元解算的点位精度在0.5 m以内,DOP值小于1,说明计算的结果较可靠．

6 结束语

本文对伪距单点定位的数学模型进行分析和研究,用C#程序实现伪距单点定位的计算,并采用实测数据作为分析,测站点单点定位的坐标精度在米级,本文中的伪距单点定位模型计算结果与O文件的近似坐标较为接近．利用双频消除电离层模型并非完全消除电离层影响[11],还存在高阶项的电离层误差,另外对流层延迟、接收机钟差、粗差探测等仍需要通过大量研究来进一步减小误差．笔者编制的单点定位程序计算的各历元单点定位的坐标,通过大量实测数据的演算,该数学模型是可行的．由于仅采用GPS单系统伪距观测值进行单点定位计算,未加入相位观测值及其它卫星系统,各历元计算的坐标仍有10余米的变化,震荡变化比较显著,还需要对数学模型进一步研究和优化．采用法方程叠加的原理计算多历元数据共同解算的单点定位坐标, 其计算结果与O文件的近似坐标较为接近,点位精度和DOP值可一并计算得到．若在计算多历元数据时,通过粗差探测的方法,将个别误差较大的历元观测值剔除,并加入相位观测值,加入其它卫星系统的观测数据,能进一步提高最终解算坐标的精度．