精密单点定位（precise point positioning, PPP）技术自问世以来得到了极大的关注^[1-4]. 得益于国际全球导航卫星服务组织（International GNSS Service, IGS）产品精度的提升，PPP静态定位精度与差分GPS相当，动态定位精度可达厘米级. PPP定位基准通常基于国际地球参考框架（international terrestrial reference frame, ITRF），这使得其在地灾监测、智能交通、精准农业、地壳运动和海啸预警等领域具有独特优势^[5].

基于PPP技术的优势，国内外先后开发了一系列相关软件，比较著名的有：喷气推进实验室（Jet Propulsion Laboratory, JPL）的GIPSY OASIS （GNSS- inferred positioning system and orbit analysis simulation software）、新不伦瑞克大学（UNB）的GAPS （GNSS analysis and positioning software）、GMV （Grupo Mecánica del Vuelo Sistemas）的magicPPP （magic precise point positioning）、加拿大自然资源部（Natural Resource Canada, NRCan）的CSRS-PPP （canadian spatial reference system precise point positioning）、加泰罗尼亚理工大学（UPC）的gLab （GNSS-lab tool）、东京海洋大学的RTKLIB（real time kinematic library）和武汉大学的PRIDE （positioning racers to image & decipher the earth）等，其中GIPSY、GAPS、CSRS-PPP和magicPPP免费提供在线数据处理服务^[6-9]. 西南交通大学也一直致力于GPS数据处理软件的研发，其中，PLAOD （precise location and orbit determination）软件主要目的是精密定位、定轨和大气反演. 相比较于常规软件，PLAOD功能较为丰富，除提供精密静态坐标外，还可提供高时间分辨率的精密动态坐标、对流层天顶延迟和电离层天顶延迟，其中，在进行大气估计时可固定测站坐标以加快其收敛速度. 此外，在提供低轨卫星天线相位中心的前提条件下，PLAOD可以对其进行精密定轨. PLAOD采用C#语言编写，所有代码均是自主研发，在联网的情况下可根据观测文件自动下载对应的精密轨道（sp3）、精密钟差（clk）和码偏差（DCB）等相关文件. 本文将主要介绍PLAOD的基本理论和性能分析.

1.   PLAOD解算模型和数据处理策略

1.1   线性组合观测量方程

对于频率为j的卫星s至接收机r间的伪距$P_{r,j}^{\left( s \right)}$和相位 $\phi _{r,j}^{(s)}$ 观测量方程可用式（1）表示.

式中：$\rho _r^{(s)}$为空间几何距离；$c$为光速；${t_r}$为接收机钟差；${t^{(s)}}$为卫星钟差；$T_r^{(s)}$为信号路径上的对流层延迟；${\mu _j} = {({\lambda _1}/{\lambda _j})^2}$为电离层系数，${\lambda _j}$为频率j对应的波长；$I_{r,1}^{(s)}$为频率1上的斜路径电离层延迟；${d_{r,j}}$和$d_j^{(s)}$分别为频率为j的接收机端和卫星端伪距硬件延迟；$Z_{r,j}^{(s)}$为整周模糊度；${\delta _{r,j}}$和$\delta _j^{(s)}$分别为频率j的接收机端和卫星端相位硬件延迟；$\varepsilon _{r,j}^{(s)}$为伪距未模型化的误差；$\xi _{r,j}^{(s)}$为相位未模型化的误差.

利用式（1）构造消电离层（ionosphere-free, IF）观测方程，如式（2）所示.

式中：下标IF表示频率j为消电离层组合的频率.

由于相位观测量中存在模糊度，式（2）中的卫星钟差基准实际上来源于伪距，即IGS等机构提供的卫星钟差为$c{t^{(s)}} + d_{{\rm{IF}}}^{(s)}$，同理，能估计的接收机钟差为$c{t_r} + {d_{r,{\rm{IF}}}}$，这样就会在相位IF观测量中引入对应伪距硬件延迟，这些偏差将被吸收到模糊度参数中，会破坏其整周属性，故直接估计的模糊度参数为浮点解. 如要获取模糊度固定解，则须提供卫星端小数偏差或者耦合卫星钟等外部产品. 相关文献证明模糊度固定可一定程度改善PPP东方向坐标估计的准确度，同时加速待估参数的收敛速度，但效果均有限^[10]. 式（2）中的对流层延迟可以进一步分解为

式中：${T_{r,{\rm{D}}}}$为利用模型获取的先验对流层干延迟；$m_{r,{\rm{D}}}^{(s)}$和$m_{r,{\rm{W}}}^{(s)}$分别为对应卫星s的对流层干延迟和湿延迟投影函数；${T_{r,{\rm{W}}}}$为天顶总延迟（zenith total delay, ZTD）扣除${T_{r,{\rm{D}}}}$后的待求天顶湿延迟；${g_{r,{\rm{n}}}}$和${g_{r,{\rm{e}}}}$分别为待求量对流层南北和东西的梯度；${a^{(s)}}$和${e^{(s)}}$分别为卫星相对于测站的方位角和高度角；${c_0}$为经验参数.

要想获得真实可降水量，必须利用实际气象数据计算得到真实的天顶干延迟${T_{r,{\rm{D}}0}}$，然后利用ZTD扣除${T_{r,{\rm{D}}0}}$，通过转换公式最终获得可降水量.

利用式（1）的相位观测方程还可构造无几何距离（geometry free, GF）观测方程，如式（4）所示.

式中：j≠1；下标GF表示频率j为无几何距离组合频率；等式右边主要有电离层和模糊度参数，电离层参数可以进一步分解，如式（5）所示.

式中：$m_{{\rm{ION}}}^{(s)}$为电离层投影函数；${I_{r,0}}$为待求测站天顶电离层延迟；${g_{r,{\rm{B}}}}$和${g_{r,{\rm{L}}}}$分别为电离层纬度梯度和经度梯度；$B_p^{(s)}$和$L_p^{(s)}$分别为卫星s对应的电离层穿刺点p的纬度和经度；${B_r}$和${L_r}$分别为测站的纬度和经度.

1.2   数据处理策略

PLAOD主要数据处理流程如图1所示.

图  1  PLAOD数据处理流程

Figure  1.  Data processing flowchart of PLAOD

下载: 全尺寸图片 幻灯片

周跳探测和粗差剔除对于最终结果至关重要，其在PLAOD中采用逐历元处理模式，以便于后续拓展实时数据处理应用；不同于传统的周跳探测算法，这里采用MW （melbourne wubbena）组合或GF组合乘以$\sqrt {\sin {e^{(s)}}}$作为探测量以提高质量控制的可靠性. 即便是伪距观测量，粗差探测也必须进行，用以提高伪距单点定位结果的准确度. 伪距单点定位在PLAOD中不仅提供PPP的初始历元先验值，还提供后续历元的钟差先验值，这样可以避免因伪距和相位同时发生钟跳引起待估参数的跳变，而单一伪距或相位钟跳在现在的接收机内部即可实现相应处理.

图1中的各项精密误差改正参照如表1所示.

表  1  误差改正

Table  1.  Error correction

误差来源	改正参考		误差来源	改正参考
相对论效应	Neil Ashby [11]		DCB （需要时）	CODE
卫星天线相位中心	igs_yy.atx		重力延迟	Neil Ashby [11]
接收机天线相位中心	igs_yy.atx		海潮	IERS2010[12]
固体潮汐	IERS2010[12]		极潮	IERS2010[12]
天线相位转绕	Wu 等, 1993[13]

下载: 导出CSV 

| 显示表格

PLAOD的数据处理核心在于扩展卡尔曼滤波^[14]. 卡尔曼滤波不仅便于实时估计参数，同时可以提高数据处理效率，以单个测站间隔30 s采样24.0 h观测数据为例，PLAOD在CPU为I7、内存为16 GB的计算机上仅需10 s即可对其执行完毕. PLAOD在进行滤波时相应的参数处理策略见表2.

表  2  PLAOD参数处理方案

Table  2.  Parameter processing strategy of PLAOD

参数	处理方案
测站坐标	静态时过程噪声为 0，动态时先验值源自伪距单点定位
模糊度	过程噪声为 0
对流层	随机游走，过程噪声默认为 4 cm2/h，动态时顾及历元间高差引起的影响
对流层梯度	随机游走
电离层	随机游走
电离层梯度	随机游走
接收机钟差	先验值源自伪距单点定位
卫星钟差	20 min弧段二阶多项式拟合
卫星坐标	10 阶滑动拉格朗日插值

下载: 导出CSV 

| 显示表格

值得一提的是表2中关于卫星钟差的处理，早期的IGS钟差产品输出间隔通常为5 min，对于高采样率的观测文件直接采用线性内插方法来获取钟差可能会导致其精度损失较大，PLAOD中利用20 min弧段内的原始钟差数据构造二阶函数则可有效减少钟差拟合精度的损失.

常规软件中常用高度角的正弦来定义观测量的随机模型（测量噪声），PLAOD中的测量噪声矩阵依赖于对流层投影函数或者电离层投影函数，以便于低高度角观测量得到更合理的利用. 以定位测量噪声R为例，其对应的矩阵如式（6）所示.

式中：0.300表示该卫星伪距观测量在天顶方向的测量噪声为0.300 m； 0.003表示该卫星的相位观测量噪声为对应伪距的1/100.

2.   软件性能分析

2.1   静态定位

选择2020年中国区域1月1日—3日期间能下载的IGS站点单天观测文件作为测试数据，其采间隔为30 s，站点分布如表3所示.

表  3  测试站地理点位置

Table  3.  Geographic location of testing stations

站点	经度/（°）	纬度/（°）
CHAN（长春）	125.443	43.790
BJFS（北京房山）	115.892	39.609
BJNM（北京）	116.224	40.245
HKSL（香港小冷水）	113.928	22.372
HKWS（香港黄石）	114.335	22.434
JFNG（九峰）	114.491	30.516
LHAZ（拉萨）	91.104	29.657
TWTF（桃园）	121.164	24.954

下载: 导出CSV 

| 显示表格

基于IGS最终精密星历和卫星钟差，利用PLAOD解算上述观测文件后，所得结果是基于ITRF的时间序列. 通常解算坐标序列最终历元的结果能在一定程度上反映软件解算的准确性和稳定性，将最终历元坐标与IGS提供的SINEX （solution independent exchange）文件坐标作差，并将其转换为站心坐标系中北、东和高方向的偏差（d_N、d_E和d_U），其结果如表4所示. 此外，坐标收敛情况作为软件完好性统计的指标之一，表4罗列了每个测站收敛到10 cm准确度所需的观测时长，需要注意的是这里的时长并不一定代表首次收敛到10 cm所需观测时长，而强调的是后续结果不再出现超过10 cm的偏差，同时10 cm表示三维距离总偏差.

表  4  站点解算坐标偏差

Table  4.  Coordinate bias of solved stations

站点	年积日	dN /mm	dE /mm	dU /mm	收敛时间 /min
CHAN	001	1.9	−9.2	−11.3	32.0
	002	2.4	−5.7	−7.5	19.0
	003	2.0	−7.1	−9.2	43.0
BJFS	001	2.0	1.2	−1.8	30.0
	002	2.5	0.7	−0.5	21.0
	003	2.5	−0.2	2.8	23.0
BJNM	001	2.6	−3.6	−0.8	29.0
	002	4.5	0.4	1.4	50.0
	003	4.0	−4.5	4.3	46.5
HKSL	001	2.0	0.5	−3.9	30.5
	002	3.8	0.1	1.7	25.0
	003	4.4	−2.5	2.0	20.5
HKWS	001	1.0	1.5	2.2	30.0
	002	3.3	1.9	−0.1	52.5
	003	3.5	−0.3	1.9	104.5
JFNG	001	2.0	3.2	0.6	20.0
	002	2.8	4.4	9.0	13.0
	003	3.1	3.4	8.1	19.0
LHAZ	001	3.6	−1.5	−3.8	42.0
	002	4.1	1.5	−0.3	57.0
	003	3.5	−1.5	−4.4	44.0
TWTF	001	0.8	3.9	−9.3	26.5
	002	4.4	1.0	−5.4	21.5
	003	4.0	0.1	−6.3	19.0

下载: 导出CSV 

| 显示表格

通过表4可以得出：各测站每天最终历元的坐标单方向偏差均在毫米级，除了CHAN站1月1日的高程方向偏差，这表明PLAOD与IGS分析中心的误差模型吻合度较高. 通过对各测站的收敛情况进行检查发现它们的收敛时间均在1.0 h以内，除了HKWS站1月3日的结果，其中，JFNG站三天均保持在20.0 min以内，最短时为13.0 min；而HKWS站1月3日的收敛时间长达104.5 min，经进一步分析发现该站伪距噪声较大.

为进一步展示PLAOD解算测站坐标偏差序列整体情况，选取上述偏差最大的结果和收敛时间最长的结果，分别对应CHAN站 1月1日和HKWS站 1月3日偏差序列，如图2和图3所示，图中：dX、dY和dZ分别为X、Y和Z方向的偏差，历元间隔为采样间隔（如无特殊说明下同）.

图  2  CHAN站2020年1月1日静态坐标偏差序列

Figure  2.  Coordinate bias series of CHAN station on 1st Jan, 2020 in static mode

下载: 全尺寸图片 幻灯片

图  3  HKWS站2020年1月3日静态坐标偏差序列

Figure  3.  Coordinate bias series of HKWS station on 3rd Jan, 2020 in static mode

下载: 全尺寸图片 幻灯片

图2显示：CHAN站1月1日偏差序列在收敛后波动较小，其收敛后的前期结果优于后期结果，这与卫星的几何形状分布等因素有一定的关联，但从整体上看经近3.5 h收敛后的偏差持续保持在毫米级. 图3显示：HKWS站1月3日偏差序列收敛时间较长，但在近14.0 h收敛后波动一直保持较小，不超过5 mm. 上述结果均证实PLAOD软件具备较好的稳健性.

2.2   动态定位

相对于静态定位，动态定位的实践应用更广泛，但在动态定位过程中，待估位置参数个数有所增加. 为检验PLAOD的动态定位性能，首先，将上述CHAN站1月1日的静态数据模拟成动态模式进行解算，所得结果如图4所示.

图  4  CHAN站2020年1月1日动态坐标偏差序列

Figure  4.  Coordinate bias series of CHAN station on 1st Jan, 2020 in kinematic mode

下载: 全尺寸图片 幻灯片

如图4所示，CHAN站1月1日动态定位偏差序列在收敛后，单方向均在cm级，将偏差序列转换成站心坐标系后显示高程方向误差占主要部分. 尽管每历元均需解算位置参数，但图4显示动态定位收敛速度也较快，仅比静态定位时慢几分钟.

另外利用2011年3月11日日本台站MIZU的观测数据来进行测试，MIZU站距离宫城震中较近，并且由于地震原因导致该站观测数据在地震后不久便发生中断，其观测数据解算结果如图5所示.

图  5  2011年3月11日MIZU站动态坐标偏差序列

Figure  5.  Coordinate bias series of MIZU station in kinematic mode on 11th Mar, 2011

下载: 全尺寸图片 幻灯片

30 s分辨率的MIZU站坐标偏差序列显示，X方向最大震时位移将近2.5 m，Y方向最大震时位移超过1.5 m，Z方向最大震时位移超过1.0 m，这与相关文献的结果高度吻合^[15]. 此外，尽管数据观测发生中断，从图4仍可显示：X、Y和Z方向的位移在一定时间后趋于稳定，这表明PLAOD不仅可以反演弹性形变，也可以探测永久形变.

虽然地震时期，站点形变的速率不是很高，而生产实践中移动载体的速率千变万化，但PLAOD不受该因素的制约，可以应用到各种精密动态定位场景中.

2.3   对流层天顶延迟

作为PPP的副产品，GPS反演的ZTD是气象分析和预报的重要参考指标. 相比于其他软件，PLAOD的ZTD过程噪声设置值较大，这主要是用以应对恶劣天气而设定的，以暴雨事件为例，ZTD在1.0 h内变化完全可以达到2 cm. 以2020年8月12日BJFS站反演的ZTD为例，现给出PLAOD、CSRS-PPP和IGS分析中心JPL提供的静态ZTD解算结果：其中PLAOD和CSRS-PPP 30 s提供一次ZTD值，均只采用向前滤波法；JPL则是5 min提供一次ZTD值，而且采用了向后滤波法，即类似批处理.

图6显示：3条曲线吻合度较高，尤其是CSRS-PPP与PLAOD对应的时序变化趋势. 由于JPL的结果采用了向后滤波，无法判断其ZTD收敛速度，不过这可以为PLAOD和CSRS-PPP的ZTD收敛速度提供参照. 易见CSRS-PPP的ZTD收敛速度较PLAOD快，这与不同软件采用的先验ZTD模型有关：CSRS-PPP采用VMF （vienna mapping function）投影函数模型，该模型的优点是精度高，缺点是需要实时更新相应文件，即每6.0 h更新一次^[16]；PLAOD采用的是GPT （global pressure and temperature）和GMF （global mapping function）模型，GPT模型准确度较VMF略低，但该模型不需要额外提供相应文件^[17-18]. 此外，在各项误差模型均准确的前提下，投影函数是影响ZTD解算准确度的主要因素，上述结果验证了GMF与VMF的高度吻合性. 大量数据验算表明PLAOD所得静态ZTD转化成可降水量对应的精度通常在1～2 mm.

图  6  2020年8月12日BJFS 站ZTD序列

Figure  6.  ZTD series of BJFS station on 12th Aug, 2020

下载: 全尺寸图片 幻灯片

PLAOD在估计ZTD时可以将测站坐标固定，以减少参数个数从而加速参数估计的收敛速度，同时也可以利用动态模式估计ZTD. 静止测站只能获取地面有限点位上空的ZTD，而移动平台能反演空间更多点位上空的ZTD，但移动平台反演的ZTD精度如何有待验证. 这里以两个观测数据来进行测试：第一个观测数据即上述TWTF站2020年1月1日对应的观测量，将该静止观测量按动态模式进行解算；第二个观测数据则来自2008年8月1日直升机载观测量，采用的接收机为拓普康，但天线类型未知，故无法对其做接收机天线相位中心改正. 利用PLAOD解算TWTF站的ZTD时，采用了常规静止模式（S模式）、动态模式（K模式）和固定测站坐标模式（F模式），参考值同样为JPL提供的ZTD，具体结果如图7所示；PLAOD在解算动态测站包括机载观测量时会根据相邻历元间的高差判断是否自动增加对流层过程噪声，噪声增加的大小依赖于GPT模型，而CSRS-PPP中动态解算时不估计ZTD，直接利用VMF模型，两软件所得机载GPS ZTD如图8所示.

图  7  2020年1月1日TWTF站ZTD序列

Figure  7.  ZTD series of TWTF station on 1st Jan, 2020

下载: 全尺寸图片 幻灯片

图  8  2008年8月1日机载GPS ZTD序列

Figure  8.  ZTD series of airborne GPS on 1st Aug, 2008

下载: 全尺寸图片 幻灯片

从图7可以看出：PLAOD解算出的F模式ZTD和S模式ZTD与JPL的结果高度吻合，其中，F模式ZTD的收敛速度更快，而K模式ZTD的解算精度有所下降，其局部波动明显，但在整体范围内也可以达到1～2 cm的准确度，对应可降水量则为2～3 mm的准确度. 从图8可以看出：PLAOD的ZTD估计值与VMF的模型值之间具有一定的吻合性，ZTD的变化与直升机的高度有较大关联，如果加入接收机天线相位中心改正能进一步提高ZTD的估计精度. 算例结果意味着移动平台GPS接收机反演ZTD具有一定的可行性，这将会为ZTD的获取提供新的途径，尤其是弥补其空间分辨率，有利于水汽层析等研究.

2.4   电离层天顶延迟

前述内容中提及PLAOD具备电离层估计功能，由于电离层估计是建立在相位GF组合观测量的基础上，忽略相位硬件延迟改正，实际上所得的电离层中可能包含部分相位硬件延迟误差. 以上述TWTF站1月1日的第一波段频率天顶电离层延迟估计为例，利用IGS的IONEX （ionosphere exchange）文件值作为参考，所得结果如图9所示.

图  9  TWTF站2020年1月1日电离层天顶延迟序列

Figure  9.  Ionospheric zenith delay series of TWTF station on 1st Jan, 2020

下载: 全尺寸图片 幻灯片

图9显示：PLAOD解算的电离层天顶延迟与IONEX对应的电离层天顶延迟具有一定的吻合度，虽然在局部会有一定的差异，但整体趋势基本相同，在当地时间PM 1：00—2：00到达峰值（中国区域内地方时比GPS时快8.0 h，地方时PM 1：00对应图中历元为600处），在凌晨降至谷值. 实际上该结果是IGS跟踪站观测量反演电离层的典型结果，实验证明相比较利用非差非组合PPP技术估计电离层，该结果与IONEX的吻合度更高，这意味着利用PLAOD监测电离层的变化具有一定的可行性.

2.5   低轨卫星精密轨道确定

随着社会科技和应用需求的发展，地球周围上空发射了越来越多的低轨卫星，这些卫星的应用是建立在其精密轨道确定的基础之上. PLAOD具备低轨卫星精密轨道确定功能，采用的方式为几何法定轨，其本质为高动态PPP，所得结果也被称为运动学轨道，该方法有利于地球重力场模型反演. 以GRACE（gravity recovery and climate experiment）双星A和B为例，利用PLAOD对2016年1月3日对应的间隔10 s采样观测量进行全历元解算，与JPL提供的简化动力学轨道值作差，其结果如图10和图11所示.

图  10  2016年1月3日GRACE （A）卫星轨道偏差序列

Figure  10.  Orbit bias series of satellite GRACE (A) on 3rd Jan, 2016

下载: 全尺寸图片 幻灯片

图  11  2016年1月3日GRACE （B）卫星轨道偏差序列

Figure  11.  Orbit bias series of satellite GRACE (B) on 3rd Jan, 2016

下载: 全尺寸图片 幻灯片

从图10和图11可看出：低轨卫星轨道准确度收敛到厘米级所需时间较短，通常只需几分钟，这得益于低轨卫星和GPS卫星均在高速飞行，使得GPS卫星对应低轨卫星的几何形状分布变化较快，有利于提高待估参数的收敛速度. 从整体上看，图中大多数历元的X、Y和Z方向偏差均在厘米级，但部分历元存在重收敛或者跳变，这主要与所接收的GPS卫星颗数以及卫星几何分布形状有关，这种情况通过向后滤波可一定程度改善，但无法完全消除. 通常在可接收卫星颗数为4且对应几何精度因子值又较大时，就会导致解算的位置参数存在较大的偏差，这在低轨卫星精密运动学轨道确定时不可避免，因为低轨卫星天线不能像地面接收机天线持续朝向天顶方向，这也是为什么国外众多科研机构提供的运动学轨道不连续的原因. 上述GRACE卫星轨道解算结果表明：和地表运动测站一样，只要低轨卫星能够观测足够的GPS卫星，加上严密的数据预处理方法，几何法定轨也能获取厘米级的运动学轨道结果.

3.   结　论

1） PLAOD的长时间静态定位精度可达毫米级，可用于高精度控制点坐标解算和形变监测等领域；

2） 在观测条件尚可的情况下，PLAOD的动态定位精度可达厘米级，适用于精密轨迹追踪和动态形变反演等领域；

3） 作为PPP的副产品，PLAOD估计的ZTD可转换为精密水汽，用于气象监测等领域；

4） 同样作为PLAOD非差解算副产品的电离层延迟，与IGS提供的电离层产品具有较高的吻合度，可用于地球上空电子活动事件分析；

5） PLAOD可解算低轨卫星精密轨道，进而为地球重力场反演等提供前期数据.

为适应GNSS应用的发展，后期将会在PLAOD中加入模糊度固定和多系统多频数据处理功能，尤其是重点拓展北斗的数据处理和应用.

致谢：感谢IGS提供观测文件和精密产品；感谢NRCan提供CSRS-PPP服务作为参考；感谢GFZ提供GRACE相关数据；感谢审稿专家提出的宝贵意见；本文受到河北省高等学校科学技术研究项目（QN2020216）资助.