注：

1）本人在学习PPP过程中，对PPPH软件内所有源码进行了注释，相关理论进行了解析，并通过本文记录，由于是学习记录有些地方注释在了源码上，所以部分理论可能不够详细，请见谅。

2）内容来源包含相关论文及其他CSDN博客，如有侵权请联系修正。

3）一些理论理解可能不够深入，如内容存在问题欢迎交流和讨论。

 软件论文、操作手册、带有注释版本源码链接如下：

PPPH软件论文翻译

PPPH软件操作手册翻译

PPPH软件注释版

 

 

目录

 

PPPH软件解析

一、数据读取（data_hand）

二、数据预处理（preprocess）

1、elm_badclk函数：

2、decimation函数

3、cal_sat函数

4、elv_mask函数

5、cs_detect函数

6、clk_jmp2函数用于钟跳探测与修复

7、outlier函数异常值检测并剔除

8、smoothing函数：相位平滑伪距（标准Hatch均值方法，其中电离层延迟用相位组合替换）

三、误差模型（nmodel）

1、相对论效应

2、时间系统说明

3、对流层天顶延迟（saastamoinen模型）

4、天线误差改正

4.1卫星天线相位中心偏差PCO改正（sat_apc）

4.2接收机天线相位中心偏差PCO改正（rec_apc）

4.3接收机天线高误差

4.4天线相位缠绕改正

四、滤波求解（MGNSS_filter）

1、dtr_satlist函数

2、dtr_satno函数

3、dtr_sys函数

4、滤波初始、过程参数设置

5、卡尔曼滤波公式

6、PPPH的抗差自适应卡尔曼滤波

---

PPPH软件解析

 

• 一、数据读取（data_hand）

data_hand函数说明：

1、观测值文件（o文件）文件头存入data.inf，读取的dcb数据在读取伪距观测值时直接改正（仅改正了GPS），读取的双频观测数据存入data.obs且通过data.obs.st来判断该历元该卫星是否同时含有p1、p2、l1、l2数据。

2、可以读取四系统的精密星历数据，同时可以读取连续三天的精密星历文件，并取前一天的后五个历元和后一天的前五个历元进行合并，以便于后续插值。

3、支持四系统的钟差文件读取，若没有钟差文件可以利用精密星历的第四列数据赋值给钟差数据。

4、可以读取四系统双频接收机和卫星的天线数据NORTH / EAST / UP，只能进行PCO改正，并分别存入data.atx.sat.neu = sat_neu./1000，data.atx.rcv.neu = rcv_neu./1000。

二、数据预处理（preprocess）

preprocess函数说明：

1、elm_badclk函数：

        对读取的钟差数据进行检验：依次对每一列即每颗卫星进行检查（第 i 列数据既包含非 NaN 值，又包含 NaN 值或者值为 999999.999999时，执行相应的操作），对于含有NaN的历元卫星将data.obs.st对应位置赋值为0，后续不做处理。

2、decimation函数

        对所有历元进行判断，如果实际历元数大于或小于data.obs.st所记录的历元数，则将多于或少于的地方赋值为空，在读取文件体时已进行过该操作。

3、cal_sat函数

        根据精密星历文件、精密钟差文件计算卫星位置和卫星钟差（钟差和星历插值方法相同，使用9阶拉格朗日插值方法）拉格朗日插值方法   拉格朗日插值方法

• 根据第一、二频点伪距/光速获得信号近似传输时间tof，再用data.obs.ep（用于存储当前历元的天内秒）减去tof获取卫星发射信号的时间nep。
• 获取当前卫星的钟差数据data.clk(:,k);根据采样间隔、和nep插值获取当前历元卫星钟差dt。
• 若卫星钟差dt为NaN，则当前历当前卫星不做处理，data.obs.st(i,k) = 0。
• 根据卫星钟差修正卫星发射信号时刻，nep = nep - dt;。
• 插值获取卫星位置，若有前后两天精密星历数据则使用entrp_orbt，否则entrp。
• 得到卫星位置XYZ放在R数组，准备进行地球自转改正R = [ X ,Y, Z ]。
• 地球自转改正链接
• 卫星位置、速度和钟差保存在data.psat，信号传输时间保存在data.tofs。

4、elv_mask函数

用于判断卫星截止高度角：

• local函数根据卫星位置和接收机概略位置计算出卫星方位角az和高度角elev。
• xyz2plh将接收机概略位置xyz转换成blh再转换成enu最后计算卫星高度角方位角。

计算卫星高度角、方位角（XYZ转ENU，ENU转RAH）

卫星高度角和方位角的计算

RTKLIB——坐标系相互转换（ecef2pos，pos2ecef，ecef2enu，enu2ecef）

5、cs_detect函数

用于双频周跳探测与修复：

① 先通过[arc]=arc_dtr(data.obs);函数将每个卫星的连续段存储在一个cell数组中，每行代表一个段的开始和结束历元。

② 逐卫星进行周跳探测与修复。

③ MW和GF组合观测值如下：

④关于MW和GF组合：

GF组合周跳探测及后续钟跳探测论文：

Cycle Slip and Clock Jump Repair with Multi-Frequency Multi-Constellation GNSS data for Precise Point Positioning

6、clk_jmp2函数用于钟跳探测与修复

详解见《GNSS精密单点定位理论方法与研究》---张小红、李星星等P111，《多模GNSS融合精密单点定位理论与方法》---蔡昌盛P78。

7、outlier函数异常值检测并剔除

采用MW法进行检验！

最小二乘曲线拟合：在二维XOY平面上，曲线由多项式f（x）表示：

在ppph的outlier函数中横轴代表历元，纵轴为y值，即L

8、smoothing函数：相位平滑伪距（标准Hatch均值方法，其中电离层延迟用相位组合替换）

三、误差模型（nmodel）

E:\2024-7-17\多模卫星导航\误差改正_20231225.pptx

所有误差模型均存在model数组中：

model各列元素：

1、 年份	2、 年积日（一年中的第几天）	3、 日内秒data.obs.ep
4、 卫星编号	5、 卫星信号传输时间(天内秒）	6、 P1P2L1L2观测值
7、 误差总和	8、 卫星位置X	9、 卫星位置Y
10、卫星位置Z	11、卫星速度VX	12、卫星速度VY
13、卫星速度VZ	14、卫星接收机近似距离	15、卫星钟误差（m）
16、卫星天线相位中心改正PCO	17、接收机天线相位中心改正PCO	18、接收机天线高误差
19、狭义相对论误差（m）	20、相位缠绕误差	21、对流层干气总延迟
22、上一历元相位缠绕误差	23、广义相对论误差（m）	24、空
25、固体潮误差	26、卫星高度角	27、卫星方位角
28、对流层湿气映射函数	29、对流层湿气北梯度系数	30、对流层湿气东梯度系数
31、对流层干气天顶延迟

Xyz转经纬度原理：xyz转blh

1、相对论效应

rel_clk（model(t,19)）-狭义相对论，rpath（model(t,23)）-广义相对论，参考《多模GNSS融合精密单点定位理论与方法》---蔡昌盛P59。

2、时间系统说明

       cal2jul将年月日转换成儒略日期Julian Day，JD，指自公元4713年1月1日格林尼治子午线（UTC中午12：00）开始起算的累计天数，详情见书《卫星导航定位原理（白皮）》---黄丁发、张勤、张小红等P32。后为了简化得到简化儒略日MJD = JD – 2400000.5。

      此外，地球力学时TDT是相对地球质心的天体运动方程所采用的时间参数，基本单位是国际秒制SI，与原子时的尺度一致。国际天文联合会决定，于1977年1月1日原子时（TAI）0时与TDT的严格关系如下：

TDT = TAI +32.184″， 而TAI = GPST + 19″

3、对流层天顶延迟（saastamoinen模型）

标准大气层模型可以表示为：
­

4、天线误差改正

       GNSS观测值测量得是卫星至接收机之间的距离，对应卫星端为卫星天线相位中心。GNSS导航电文星历数据对应卫星位置为天线相位中心，而IGS提供的精密星历是基于卫星轨道力模型计算的，给出的是卫星质心的坐标。因此，当采用IGS精密星历进行数据处理时，需要考虑卫星质心与天线相位中心之间的差异。卫星天线相位中心与卫星质心之间的偏差称之为卫星天线偏差。同样在接收机端，GNSS观测值对应的是接收机的天线相位中心，而接收机的目标测量位置为天线参考点（Antenna Reference Point，ARP），数据处理时也需要考虑他们之间的差异。将接收机天线相位中心与接收机天线参考点之间的偏差称为接收机天线偏差。

       天线相位中心会随着卫星信号的高度角和方位角不同而不同，因而测量时会产生多个瞬时相位中心，因此通常根据瞬时相位中心的大致位置定义一个平均相位中心点。平均相位中心与卫星质心或与接收机天线参考点之间的偏差称为天线相位中心偏差（PCO-Phase Center Offset）。瞬时天线相位中心（对应单个观测值）与平均相位中心间的差异称为天线相位中心变化（PCV-Phase Center Variation）。对GNSS观测值进行天线改正时，必须同时考虑天线相位中心偏差PCO和天线相位中心变化PCV。
      如下图所示，设a为PCO矢量（天线参考点指向平均天线相位中心），ρ0 为卫星至接收机的单位向量，则PCO对相位观测值的影响可表示为：

ΔPCO=a·ρ0

      天线相位中心变化PCV对相位的影响则与方位角α、天顶距z和载波频率f 有关，可用下式表示为，而总的天线偏差是PCO和PCV的综合影响，将观测值改正至ARP点。

ΔPCV=ΔPCV(α,z,f)

       为了给用户提供相应的天线偏差改正方法，IGS相继推出了相对天线相位中心改正模型和绝对天线相位中心改正模型。两种莫小姐均是根据天线的类型、信号频率、卫星天顶距等，确定各种天线在对应不同频率、不同天顶距、方位角的天线偏差的坐标改正值或距离改正值，并将测定数据编列成表，提供给用户，用户则根据卫星类型或接收机接收信号的参数内插确定对应的天线偏差改正值。

4.1卫星天线相位中心偏差PCO改正（sat_apc）

天线相位中心偏差改正可用通过下式进行计算：

        式中ex  ey  ez 分别为星固坐标系轴在惯性坐标系中的单位矢量，Xphase 、Xmass 分别为惯性坐标系中以相位为中心的卫星坐标和以质量为中心的卫星坐标，Xoffset 为卫星天线相位中心在星固系的偏差PCO。

PCO改正值在天线文件  NORTH / EAST / UP 读取：

对于PCO改正，由于IGS提供的是基于星固系下的PCO值，因此需要将其转换至地固系。

        星固系原点位于卫星质心，z轴平行于天线指向地心，y轴平行于太阳帆板平面，x轴构成右手系。如图所示，太阳矢量从卫星质心指向太阳，卫星运动过程中z轴始终指向地球，y轴始终垂直于太阳矢量即垂直于太阳、地球、卫星组成的平面。太阳帆板可绕轴旋转以保证其垂直于太阳光线方向来最好地收集太阳能。

形式一：太阳角 β 定义为z轴和太阳单位矢量 nsum 的夹角，在定轨时需要用 β 计算太阳光压。星固坐标系坐标轴单位矢量表示为（ex   ey   ez），则太阳单位矢量可表示为nsum=( sinβ,0,cosβ)。
r = [X Y Z]T为卫星质心坐标，rsun=[Xsun Ysun  Zsun ]T为太阳坐标。

形式二：（ex   ey   ez）计算方式如下：

形式三：非矢量形式计算（ex   ey   ez）

4.2接收机天线相位中心偏差PCO改正（rec_apc）

        接收机天线相位中心矫正的过程与方法与卫星端类似（天线参考点位置 = 天线相位中心位置 - 天线相位中心偏差（PCO）；几何距离 = 观测距离 - 天线PCV + 其他改正值）。需要说明的是，IGS天线文件中给出的PCO值是在测站地平坐标系下的三个分量，因此需要将其转换到地心地固系，再进行PCO改正。

 

4.3接收机天线高误差

由观测文件头文件读取：      0.0280   0.0000  0.0000    ANTENNA: DELTA H/E/N

与4.2接收机天线相位中心偏差改正算法相同。

4.4天线相位缠绕改正

RTKLIB 中的天线相位缠绕误差修正

GPS从入门到放弃（二十三）、相位缠绕

      GNSS卫星信号采用右极化方式，这种极化方式的信号使得观测到的载波相位与卫星、接收机天线的朝向有关。接收机和卫星天线任何一方的旋转将使载波相位产生最大一周的变化，也就是在距离上一个波长的变化，这种效应称为天线旋转相位增加效应。对该效应改正通常也称为天线相位缠绕改正（wind-up）。

      接收机天线除非在动态定位中它在运动，通常是指向一个固定方向（北方向），但对于卫星天线，由于它的太阳能面板需要始终朝向太阳，因而它在运动过程中会慢慢旋转，这就造成卫星和接收机之间的几何距离发生变化。除此之外，在日食期间，卫星为了能重定向太阳能板朝向太阳，将快速旋转，在半小时之内，旋转将达到一周，在此期间，需要对相位数据进行改正或者将这段数据去除。
       该误差对于差分定位通常忽略，而对于PPP影响十分显著能达半个波长。通常在动态定位或导航中接收机天线才会发生旋转，天线相位缠绕的误差将会转移到接收机钟差解中，改正公式如下：

完整的相位缠绕改正Δϕ 是由其上个历元和本历元的小数部分共同组成：

理解：

四、滤波求解（MGNSS_filter）

1、dtr_satlist函数

   生成一个历元数（行）*1列的cell数组，每个cell数组（1*n）记录着当前历元可用卫星号。

2、dtr_satno函数

生成一个历元行 * 1列的数组，每行记录着当前历元共多少个可用卫星。

3、dtr_sys函数

根据选择的卫星系统以及是否估计对流层水平梯度误差输出每颗卫星的待估参数个数bp以及选择系统的指数sit，后续根据这个指数确定所选卫星系统。

待估参数个数bp：三个方向xyz偏差量（3个） + 接收机钟差*gnss系统数量（1*n个） + 接收机天顶对流层（1个） + 对流层水汽梯度（2个）

系统指数sit：

GPS：1	GLO：2	GPS+GLO：3	GPS+GAL：4
GPS+BDS：5	GLO+GAL：6	GLO+BDS：7	GPS+GLO+GAL：8
GPS+GLO+BDS：9	GPs+GAL+BDS：10	GLO+GAL+BDS：11	GPS+GLO+GAL+BDS：12

4、滤波初始、过程参数设置

5、卡尔曼滤波公式

6、PPPH的抗差自适应卡尔曼滤波

参考文献：Adaptive robust Kalman filtering for precise point positioning---Guo F, Zhang X (2014b)