星间频间伪距定位模型的精度分析及程序实现*

龚栎澎，杨久东，李 韧，梁 鹏

(华北理工大学 矿业工程学院，河北 唐山 063210)

伪距定位是通过利用GNSS接收机的伪距观测数据计算测站坐标的定位方法，现在已经发展成为比较成熟的技术，主要用于导航定位和差分定位中初始大概位置的定位[1]。随着精密单点定位(PPP)的逐渐深入研究，伪距单点定位方法及理论作为其基础方法及定位理论，对伪距单点定位方法的深入研究是极其重要的[2-3]。伪距定位已经不仅仅局限于传统的伪距定位，通过研究的不断深入，已经衍生出了多种伪距定位解算模型[4-5]。研究发现，采用最小二乘迭代法对伪距进行解算,其迭代次数少但定位精度较低[6-7]，所以采用卡尔曼( Kalman) 滤波法解算可使得精度有所提高，但迭代次数过多[8]；很多学者研究使用最小二乘法和Kalman滤波法组合使用的方法进行伪距解算，减少迭代次数，提高定位精度[9]。伪距单点定位方法的精度一般为米级，组合伪距观测值定位精度可高于单一伪距观测值精度[10]，GPS静态精密单点定位模型，点位精度可达毫米级[11]。BDS+GPS双系统组合伪距定位模型的定位精度可高于单一系统伪距定位精度[12-13]。

伪距单点定位方法及理论的研究虽然很多[14]，但关于伪距定位的数学模型的表述多不够详细，笔者通过对伪距定位的进一步研究，整理出加权最小二乘伪距定位和星间频间伪距定位模型，以Python3.8为工具进行编程实现，并利用实测数据验证该数学模型的可行性与正确性以及程序的稳定性，同时将两种伪距定位方法的定位结果进行比较，分析两种伪距定位方法的定位效果及精度。

1 伪距定位模型

笔者通过对伪距单点定位模型的研究，利用接收机中的GPS卫星系统的观测文件o文件中的伪距观测值和星历文件P文件，采用无电离层模型较短基线测量的实验方法，分别计算两种伪距单点定位模型中的测站坐标。

1.1 加权最小二乘伪距单点定位模型

若在历元时刻t，接收机接收到n颗GPS卫星，选取高度角最高的卫星r为参考卫星，其余的卫星记作卫星s，对于参考卫星r列出该卫星与地面测站点在L1频率的伪距观测方程[15-16]，如式(1)所示：

(1)

(2)

(3)

可得式(4)：

(4)

对未知参数进行权阵为高度角定权模型的加权最小二乘求解，求得坐标结果如式(5)所示：

(5)

将计算出的伪距定位结果与(x0,y0,z0)进行比较，从而判断各分量位置差值的大小，如果各分量差值大于0.01 m，则令x0=Xi，返回卫星位置继续进行迭代计算，如果小于0.01 m，则退出迭代计算，输出此历元的定位结果，继续下一历元的坐标计算。

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

(6)

1.2 星间频间伪距定位模型

文献中根据加权最小二乘传统伪距定位模型(星间单差伪距单点定位模型)提出了一种利用同一卫星系统中的卫星间及频率间的双差伪距定位模型[17]。由伪距定位原理可知，将式(2)～式(4)进行变形整理后可得式(7)：

(7)

根据参考卫星和其他卫星的观测方程进行星间作差，可消去式(7a)和式(7b)的公共项cVtR,如式(8)所示：

(8)

GPS卫星系统信号播发频率L1是调制基本频率10.23 MHz×154获得的频率，L2是调制基本频率10.23 MHz×120获得的频率，如式(9)所示：

(9)

(10)

依据公式进行最小二乘求解，因为星间频间伪距定位方法是依据L1、L2频率间作差进行计算定位，但是卫星的频率发射是不稳定的，存在一定的频率波动范围，在后期解算中需要进行一定的系数调整后再计算出测站的伪距定位坐标。在进行最小二乘解算中的权阵采用同一历元下卫星高度角定权模型，权阵为方差阵的逆，其方差公式如式(11)所示：

(11)

式中,σ为某颗卫星方差;θ为卫星高度角。根据误差传播定律可推导出式(12)：

D=(774+604)

(12)

式中,参考星为卫星高度角最大的卫星，其n-1余颗卫星的方差依据其高度角求得。再计算出权阵后通过式(4)～式(6)求出定位坐标结果以及均方根误差。

2 伪距定位的程序设计

在利用Python语言进行两种伪距定位模型的程序设计时，如图1所示，建立了不同的类来实现伪距定位坐标的计算，主要包括文件读取类、坐标转换类、时间统一类、测站类、伪距定位模型计算类、精度评价类等六大类。

3 定位坐标解算及精度分析

运用本文所述的两种伪距定位模型，利用2019年11月在河北省唐山市曹妃甸生态城控制网静态观测数据，选用C401、C410和GPS1共3个控制点，C401与C410间距5 493 m，C401与GPS1间距1 950 m，C410与GPS1间距4 348 m。同步观测时间为100 min，静态接收机型号为海星达IRTK2型双频接收机，且经过质检鉴定，静态数据文件为RINEX3.02文件，采样间隔1 s，高度截止角为10°，因观测时间较长，故选取观测历元中的300个历元参与模型计算。程序自动计算采用从1～300个历元计算两个伪距定位模型中的定位坐标，同时使用中海达HGO基线解算软件进行控制网的解算并进行严密二维约束平差后，最弱点中误差的最大值为1.19 mm，符合精度要求，相对伪距定位模型的米级定位结果，可以认为HGO解算的二维约束平差结果为真值，并将两种伪距定位模型的定位结果进行比较，如图(2)～图(10)所示。

图(2)～图(4)为测站C401的XYZ坐标值，图(5)～图(7)为测站C410的XYZ坐标值，图(8)～图(10)为测站GPS1的XYZ坐标值。从图中分析可知，三个测站的加权最小二乘伪距定位模型和星间频间伪距定位模型的每个历元解算的定位坐标都出现了不同程度的波动，符合伪距定位的解算规律。在两种定位方法模型解算的300个历元坐标中，加权最小二乘方法的图线波动较小，图形表现较为平稳，最大波动较差不足5 m，但是均在精密坐标图线下方，反观星间频间伪距定位模型解算的300历元定位坐标，波动较大，大部分图线波动在精密坐标图线的上方与下方，其平均值与精密坐标图线更为接近。

图1 程序设计流程图

图2 测站C401的X坐标值

图3 测站C401的Y坐标值

图4 测站C401的Z坐标值

图5 测站C410的X坐标值

图6 测站C410的Y坐标值

图7 测站C410的Z坐标值

图8 测站GPS1的X坐标值

图 9 测站GPS1的Y坐标值

图10 测站GPS1的Z坐标值

为了更加直观地分析两种伪距定位方法在3个测站的定位精度，绘制表1，其中，ΔX、ΔY、ΔZ为300历元的解算坐标平均值与精确坐标的差值，DX、DY、DZ为300历元解算坐标的变化范围。可以看出，3个测站在经过300个历元解算后，对比两种伪距定位方法的坐标平均值与精密平差的坐标的较差值发现，加权最小二乘伪距定位的X坐标精度要优于星间频间伪距定位模型，但是Y、Z的平均值坐标解算精度均是星间频间伪距定位的定位精度较高，加权最小二乘伪距定位的最小较差为0.863 7 m，最大较差为4.781 2 m，星间频间伪距定位的最小较差为0.451 1 m，最大较差为2.526 4 m，但在解算中的坐标变化范围波动较大，最大波动达到14.403 7 m，而加权最小二乘伪距定位最大波动仅为4.507 7 m。

表1 测站定位坐标误差及坐标变化范围/m

RMSE值反映观测值偏离真实有效数据值的波动幅度情况，这一值越小误差就越低，精度越高，伪距定位越准确。通过表2分析可以得出，用两种伪距定位方法在3个测站进行300个历元的解算，3个测站的RMSE值表现一致，均是加权最小二乘伪距定位的X坐标解算精度高于星间频间伪距定位，而星间频间伪距定位在3个测站的RMSEY、RMSEZ、RMSEXYZ的值均低于加权最小二乘伪距定位。综合分析，星间频间伪距定位方法的定位精度要高于加权最小二乘伪距定位，但是X坐标解算精度低于加权最小二乘伪距定位，主要原因可能在于选择的解算模型不够精确，在X定位计算中发生偏移，导致定位精度不高。

表2 测站RMSE值

4 结 论

本文对两种伪距定位的数学模型进行分析和研究，用Python程序实现了加权最小二乘伪距定位模型、星间频间伪距定位模型的计算，并利用静态实测数据进行分析，测站点的伪距定位精度在米级，加权最小二乘伪距定位方法计算的XYZ坐标误差在5 m内，星间频间伪距定位方法计算的XYZ坐标误差在3 m内，星间频间伪距定位的精度高于加权最小二乘伪距定位。经过大量的实测数据的计算检验，笔者编制的伪距定位程序是可靠的，其数学模型是可行的。由于实验数据区域较小，属于短基线解算，对流层及电离层误差影响较小，且仅采用了GPS单系统伪距观测值进行伪距计算，未加入其它卫星系统进行联合解算。若在后期加入其它卫星系统进行有效组合伪距定位，能进一步提高伪距定位精度及稳定性。