一种基于并行计算的实时北斗相位小数偏差快速估计方法

1.本发明涉及卫星导航定位技术领域，尤其是一种基于并行计算的实时北斗相位小数偏差快速估计方法。


背景技术：

2.北斗卫星导航系统(bds，beidou navigation satellite system)是我国自主建设、独立运行的全球导航卫星系统(gnss，global navigation satellite system)，是国家重要空间基础设置。当前，基于北斗位置服务的产业年产值已高达四千亿元，并呈现稳定高速增长态势，无人驾驶、智慧城市建设、自然资源调查等均需要可靠的空间位置信息。
3.精密单点定位(ppp，precise point positioning)技术
1.可以在全球范围内获取地心参考框架下精确三维坐标，然而由于其模糊度参数吸收了伪距和相位硬件延迟，失去整周特性的ppp浮点解的精度和可靠性亟需进一步提高。
4.为实现ppp模糊度可靠固定，国内外学者先后提出相位小数偏差模型
2.、整数相位钟模型
3.、钟差解耦模型
4.等，并从服务端改正数类型和客户端模型自由度等方面对上述三种方法进行了对比。由于bds系统建设周期较短，且采用三种混合轨道类型卫星，包括地球静止轨道卫星(geo)、倾斜地球同步轨道卫星(igso)和中圆地球轨道卫星(meo)，卫星观测值质量、卫星姿态模型、卫星空间环境等与其余gnss系统存在一定差异，不顾及两者之间的差异降低了bds ppp定位性能。
5.同时，为了满足实时位置服务的时效性和长期稳定性，bds时变相位小数偏差估计的效率和数值稳定性显得尤为重要。近年来，众多学者采用引入外部矩阵库、非严格的网解模型、降低相位小数偏差更新时间等方法，提高相位小数偏差估计效率。然而，上述技术手段依赖于硬件平台，难以实现跨平台操作，降低了相位小数偏差的精度和可靠性，不利于bds ppp模糊度可靠固定，而严密的网解模型存在高维矩阵求逆和方差-协方差耗时等问题。
6.因此，可靠估计实时bds相位小数偏差实现可靠的bds ppp固定解，对于推动我国北斗提供优质的导航、定位和授时服务具有重要的实用价值和研究意义。


技术实现要素：

7.本发明所要解决的技术问题在于，提供一种基于并行计算的实时北斗相位小数偏差快速估计方法，能够实现实时bds相位小数偏差快速、稳健估计。
8.为解决上述技术问题，本发明提供一种基于并行计算的实时北斗相位小数偏差快速估计方法，包括如下步骤：
9.步骤1、数据流时间同步与有效性检核；
10.步骤2、根据步骤1获得的同步实时观测值并更新外部所需文件，多测站并行计算mw平均值和无电离层组合浮点模糊度；
11.步骤3、根据步骤2得到的多测站mw平均值，作为输入观测值，进行时变宽巷相位小
数偏差快速估计；
12.步骤4、根据步骤2得到的多测站无电离层组合浮点模糊度，作为输入观测值，进行时变窄巷相位小数偏差快速估计；
13.步骤5、线性转换各频相位小数偏差。
14.优选的，步骤1中，数据流时间同步与有效性检核具体包括如下步骤：
15.步骤11、获取本地时间，设置相位小数偏差更新率和数据流等待时间，以此作为逐历元获取同步实时数据流的时间参考；
16.步骤12、利用广播星历和实时轨道/钟差改正数，恢复实时卫星轨道和钟差数据，并监测其数据是否存在缺失和异常；
17.步骤13、检测更新多系统差分码偏差dcb文件、预报地球定向参数eop文件、天线文件、blq格式的海潮文件，作为后续数据预处理的外部输入数据；
18.步骤14、针对部分测站观测值缺失问题，特别是连续多个历元出现观测值数据中断现象，设置等待时间，在等待时间内该测站数据不参与后续相位小数偏差计算。
19.优选的，步骤2中，多测站并行计算浮点模糊度具体包括如下步骤：
20.步骤21、删除双频伪距和相位观测值缺失卫星，逐颗形成如下双频线性组合观测值：
[0021][0022]
式中，f1和f2分别表示bds第一频和第二频频率，p1/p2和l1/l2分别为bds第一频和第二频伪距观测值和相位观测值，pi为第一频率和第二频率的伪距观测值的差值，li为第一频率和第二频率的相位观测值的差值，由双频载波相位l1/l2和伪距观测值p1/p2构造的mw组合观测值，pc为双频无电离层伪距组合观测值，lc为双频无电离层相位组合观测值；
[0023]
步骤22、顾及geo/igso/meo三种轨道类型卫星的轨道高度和pi值阈值进行伪距粗差探测，利用mw组合和li组合探测周跳，考虑geo卫星的近乎静地特性，li线性组合探测周跳的电离层变化率阈值设置为meo/igso的一半，对超过阈值的观测数据进行剔除，并记录每个连续观测弧段的弧段数目及新旧弧段切换时间；
[0024]
步骤23、修正各种误差模型，包括对流层延迟、天线相位中心改正、相位缠绕、潮汐效应、卫星多路径空间误差，其中geo卫星姿态采用零偏模式，igso/meo卫星采用动偏模式，并采用移动滑动窗口计算逐历元mw平均值m
mw
及其中误差σ
mw
：
[0025]
[0026]
式中，k表示连续弧段内的窗口大小，最大移动窗口大小设置为120；
[0027]
步骤24、按照公式(3)设置geo/igso/meo三种轨道类型卫星观测值先验方差，由于接收机和不同轨道卫星环境不同，分别设置三类卫星的伪距和相位观测值的方差经验值σ0为3/0.3/0.3m和0.03/0.003/0.003m，其中，e表示卫星高度角；
[0028][0029]
步骤25、利用步骤23中修正各种模型误差后的pc和lc，采用如公式(4)中的无电离层组合ppp模型，利用igg-iii权函数模型，根据标准化残差判断并降低异常观测值的权重，从而获得可靠的无电离层组合浮点模糊度a
if
；
[0030][0031]
式中，dtr表示测站r的接收机钟差，zr和m分别表示测站天顶对流层湿延迟和湿延迟投影函数，ε
pc
和ε
lc
分别表示pc和lc观测值未建模误差；
[0032]
步骤26、当有多个测站的观测数据时，按照步骤21-步骤25的步骤，采用openmp并行处理，即可获得同步后的多测站mw平均值和无电离层组合浮点模糊度。
[0033]
优选的，步骤3中，时变宽巷相位小数偏差快速估计具体包括如下步骤：
[0034]
步骤31、以步骤2中获取的多测站mw平均值作为输入观测值，则宽巷浮点模糊度表达为：
[0035][0036]
式中，表示多测站mw平均值，λw表示宽巷波长，b
w,*
表示接收机宽巷相位小数偏差，表示卫星宽巷相位小数偏差，nw是对应的整周宽巷模糊度；
[0037]
假设共观测到q个接收机和p颗卫星，以宽巷整周模糊度和卫星、接收机宽巷相位小数偏差为估计参数，构建如下网解宽巷相位小数偏差函数模型：
[0038][0039]
式中，i表示单位矩阵，下标表示其维数，e
p
表示元素均为1的列向量，下标表示其行数，表示网解宽巷相位小数偏差的kalman滤波设计矩阵；
[0040]
步骤32、连续弧段的模糊度是一个不变的常数，其过程噪声应该为0，因此将连续弧段的模糊度参数设置为常数模型，由于接收机和卫星的小数偏差参数的时变特性，将接收机宽巷相位小数偏差和卫星接收机宽巷小数偏差分别设置为随机游走，其谱密度分别为1.0e-4
和1.0e-6
m2/s，利用状态转移矩阵的对角特性，进行kalman滤波时间更新的标量计算，即在上个历元的参数估计基础上，进行当前历元的状态预测，以实现快速的数值计算；
[0041]
步骤33、按照mw平滑值的中误差σ
mw
进行升序排列，以接收机、卫星、宽巷浮点模糊度构成无向连通图，引入kruskal最小生成树算法生成稳健的独立宽巷模糊度基准；
[0042]
步骤34、以宽巷模糊度基准值作为虚拟观测值，其虚拟方差设置为1.0e-9
，卫星相位小数偏差设置为0，即作为约束添加至宽巷相位小数偏差网解模型，利用kalman滤波估计每颗卫星相对于该基准约束下的宽巷相位小数偏差；
[0043]
步骤35、以每个测站每颗卫星mw平滑值及选择的宽巷模糊度基准值视为单个观测值，按照公式(7)实现kalman滤波标量计算，避免了标准滤波存在高维矩阵求逆耗时难题；并采用openmp并行计算实现滤波方差-协方差矩阵的快速更新，进一步加速网解相位小数偏差计算效率；
[0044][0045]
式中，为第k个历元的第i颗卫星观测值方差，和分别表示kalman滤波状态预测值和滤波值，和分别为状态对应的方差阵，i表示单位矩阵，k
k,i
表示kalman滤波增益矩阵，l
k,i
表示第k个历元的第i颗卫星观测值，a
k,i
为公式(6)中矩阵该测站第i颗卫星观测值对应的行向量；
[0046]
步骤36、以宽巷模糊度估计值偏差为0.15周、标准差为0.15周和模糊度固定为整数的概率1000为依据，若某颗卫星按步骤35进行滤波得到的宽巷模糊度满足上述要求，则将其固定为整数，遍历所有估计的宽巷模糊度，直至无新的宽巷模糊度可以被固定为止，输出此时的卫星宽巷相位小数偏差此步骤获取的固定的宽巷模糊度是步骤4的前提。
[0047]
优选的，步骤4中，时变窄巷相位小数偏差快速估计具体包括如下步骤：
[0048]
步骤41、以步骤3固定的宽巷整周模糊度为已知值，逐颗卫星修正步骤2中获取的无电离层浮点模糊度，如公式(8)所示：
[0049][0050]
式中，即为修正后的窄巷模糊度浮点值，λw和λn为bds宽巷和窄巷波长，以为输入观测值，以窄巷整周模糊度和卫星、接收机窄巷相位小数偏差为估计参数，构建的网解窄巷相位小数偏差函数模型的设计矩阵与公式(6)一样；
[0051]
步骤42、参考步骤3中步骤，使得窄巷相位小数偏差估计与宽巷相位小数偏差估计保持一致，便于程序扩展使用，并输出此时的卫星窄巷相位小数偏差b
1s
。
[0052]
优选的，步骤5中，线性转换各频相位小数偏差具体包括如下步骤：
[0053]
利用步骤3和步骤4获取的宽巷和窄巷相位小数偏差，依据公式(1)中mw组合和lc组合其对应频率的系数，按照公式(9)和公式(10)将其转向至各频相位小数偏差并标记其伪距观测值通道，使得用户可以采用灵活的函数模型实现广域单测站模糊度固定；
[0054][0055][0056]
式中，和分别为转化后的第一频和第二频相位小数偏差，即实时北斗相位小数偏差，dcb即为相应伪距观测值通道上的卫星伪距差分码延迟。
[0057]
本发明的有益效果为：本发明采用严格的数据质量控制策略，基于并行计算和单观测值迭代更新算法，实现了多测站网解时变相位小数偏差快速估计，既提高了实时北斗相位小数偏差估计的效率，又提高了估计的相位小数偏差的精度和可靠性。
附图说明
[0058]
图1为本发明的方法流程示意图。
[0059]
图2为本发明逐测站并行数据预处理与浮点模糊度计算的流程示意图。
[0060]
图3为本发明时变宽巷和窄巷小数偏差快速估计的流程示意图。
[0061]
图4为本发明独立模糊度基准选取示意图。
[0062]
图5为本发明利用openmp实现滤波方差-协方差快速更新的伪代码。
具体实施方式
[0063]
如图1所示，一种基于并行计算的实时北斗相位小数偏差快速估计方法，包括如下步骤：
[0064]
(一)数据流时间同步与有效性检核；
[0065]
该步骤具体包括：
[0066]
(1-1)获取本地时间，设置相位小数偏差更新率和数据流等待时间，以此作为逐历元获取同步实时数据流的时间参考；
[0067]
(1-2)利用广播星历和实时轨道/钟差改正数，恢复实时卫星轨道和钟差数据，并监测其数据是否存在缺失和异常；
[0068]
(1-3)检测更新多系统差分码偏差(dcb，differential code bias)文件，预报地球定向参数(eop，earth orientation parameters)文件，天线文件，blq格式的海潮文件，作为后续数据预处理的外部输入数据；
[0069]
(1-4)针对部分测站观测值缺失问题，特别是连续多个历元出现观测值数据中断现象，设置等待时间，在等待时间内该测站数据并不参与后续相位小数偏差计算。完成该步骤后，即可获得同步的实时观测值并更新外部所需文件，为步骤2的逐测站并行计算准备可靠的数据支持。
[0070]
(二)多测站并行计算浮点模糊度；
[0071]
逐测站并行数据预处理与浮点模糊度计算的流程如图2所示，该步骤具体包括：
[0072]
(2-1)删除双频伪距和相位观测值缺失卫星，逐颗卫星形成如下双频线性组合观测值：
[0073][0074]
式中，f1和f2分别表示bds第一频和第二频频率，p1/p2和l1/l2分别为bds第一频和第二频伪距观测值和相位观测值，pi为第一频率和第二频率的伪距观测值的差值，li为第一频率和第二频率的相位观测值的差值，mw(melbourne-wu
¨
bbena)观测值为双频载波相位l1/l2和伪距观测值p1/p2构造的组合观测值，pc为双频无电离层伪距组合观测值，lc为双频无电离层相位组合观测值。
[0075]
(2-2)顾及geo/igso/meo三种轨道类型卫星的轨道高度和pi值阈值(60m)进行伪距粗差探测，利用mw组合和li组合探测周跳，考虑geo卫星的近乎静地特性，li线性组合探测周跳的电离层变化率阈值设置为meo/igso的一半，对超过阈值的观测数据进行剔除，并记录每个连续观测弧段的弧段数目及新旧弧段切换时间；
[0076]
(2-3)修正各种误差模型，包括对流层延迟、天线相位中心改正、相位缠绕、潮汐效应、卫星多路径空间误差(此误差改正仅仅bds-2igso/meo卫星)等，其中geo卫星姿态采用零偏模式，igso/meo卫星采用动偏模式，并采用移动滑动窗口计算逐历元mw平均值m
mw
及其中误差σ
mw
：
[0077][0078]
式中，k表示连续弧段内的窗口大小(最大移动窗口大小设置为120)。
[0079]
(2-4)按照公式(3)设置geo/igso/meo三种轨道类型卫星观测值先验方差，由于接收机和不同轨道卫星环境不同，分别设置三类卫星的伪距和相位观测值的方差经验值σ0为3/0.3/0.3m和0.03/0.003/0.003m，其中，e表示卫星高度角。
[0080][0081]
(2-5)利用(2-3)步骤中修正各种模型误差后的pc和lc，采用如式(4)中的无电离层组合ppp模型，利用igg-iii权函数模型，根据标准化残差判断并降低异常观测值的权重，从而获得可靠的无电离层组合浮点模糊度a
if
。
[0082][0083]
式中，dtr表示测站r的接收机钟差，zr和m分别表示测站天顶对流层湿延迟和湿延迟投影函数，ε
pc
和ε
lc
分别表示pc和lc观测值未建模误差。
[0084]
(2-6)当有多个测站的观测数据时，按照(2-1)～(2-5)的步骤，采用openmp并行处理，即可获得同步后的多测站mw平均值和无电离层组合浮点模糊度a
if
，此步骤获得的mw平均值和无电离层组合浮点模糊度分别作为后续步骤(三)和步骤(四)的网解模型的输入观测值。
[0085]
(三)时变宽巷相位小数偏差快速估计；
[0086]
时变宽巷小数偏差快速估计的流程如图3所示，该步骤具体包括：
[0087]
(3-1)以步骤(二)中计算的多测站mw平均值作为输入观测值，则宽巷浮点模糊度可表达为：
[0088][0089]
式中，表示多测站mw平均值，λw表示宽巷波长，b
w,*
表示接收机宽巷相位小数偏差，表示卫星宽巷相位小数偏差，nw是对应的整周宽巷模糊度。
[0090]
假设共观测到q个接收机和p颗卫星，以宽巷整周模糊度和卫星、接收机宽巷相位小数偏差为估计参数，构建如下网解宽巷相位小数偏差函数模型：
[0091][0092]
式中，i表示单位矩阵，下标表示其维数，e
p
表示元素均为1的列向量，下标表示其行数，表示网解宽巷相位小数偏差的kalman滤波设计矩阵。
[0093]
(3-2)连续弧段的模糊度是一个不变的常数，其过程噪声应该为0，将连续弧段的模糊度参数设置为常数模型，由于接收机和卫星的小数偏差参数的时变特性，将接收机宽巷相位小数偏差和卫星接收机宽巷小数偏差分别设置为随机游走，其谱密度分别为1.0e-4
和1.0e-6
m2/s,利用状态转移矩阵的对角特性，进行kalman滤波时间更新的标量计算，即在上个历元的参数估计基础上，进行当前历元的状态预测，以实现快速的数值计算。
[0094]
(3-3)按照mw平滑值的中误差σ
mw
进行升序排列，以接收机、卫星、宽巷浮点模糊度构成无向连通图，引入kruskal最小生成树算法生成稳健的独立宽巷模糊度基准，图4为独立模糊度基准选取示意图，方框和圆圈分别表示测站和卫星，黑色加粗线条即表示该模糊度被选取为独立模糊度基准；
[0095]
(3-4)以宽巷模糊度基准值作为虚拟观测值，其虚拟方差设置为1.0e-9
,卫星相位小数偏差设置为0，即作为约束添加至宽巷相位小数偏差网解模型，利用kalman滤波估计每颗卫星相对于该基准约束下的宽巷相位小数偏差；
[0096]
(3-5)以每个测站每颗卫星mw平滑值及选择的宽巷模糊度基准值视为单个观测值，按照公式(7)实现kalman滤波标量计算，避免了标准滤波存在高维矩阵求逆耗时难题；并采用openmp并行计算实现滤波方差-协方差矩阵的快速更新，其伪代码如图5所示，以进
一步加速网解相位小数偏差计算效率；
[0097][0098]
式中，为第k个历元的第i颗卫星观测噪声的方差，和分别表示kalman滤波状态预测值和滤波值，和分别为状态对应的方差阵，i表示单位矩阵，k
k,i
表示kalman滤波增益矩阵，l
k,i
表示第k个历元的第i颗卫星观测值，a
k,i
为式(6)中矩阵该测站第i颗卫星观测值对应的行向量。
[0099]
(3-6)以宽巷模糊度估计值偏差(0.15周)、标准差(0.15周)和模糊度固定为整数的概率(1000)为依据，若某颗卫星按(3-5)进行滤波得到的宽巷模糊度满足上述要求，则将其固定为整数，遍历所有估计的宽巷模糊度，直至无新的宽巷模糊度可以被固定为止，输出此时的卫星宽巷相位小数偏差此步骤获取的固定的宽巷模糊度是步骤(四)的前提。
[0100]
(四)时变窄巷相位小数偏差快速估计；
[0101]
与上一步骤类似，时变窄巷小数偏差快速估计的流程与图3较为一致，主要区别在于输入数据，该步骤具体包括：
[0102]
(4-1)以步骤(三)固定的宽巷整周模糊度为已知值，逐颗卫星修正步骤(二)中获取的无电离层浮点模糊度，如式(8)所示：
[0103][0104]
式中，即为修正后的窄巷模糊度浮点值，λw和λn为bds宽巷和窄巷波长。以为输入观测值，以窄巷整周模糊度和卫星、接收机窄巷相位小数偏差为估计参数，构建的网解窄巷相位小数偏差函数模型的设计矩阵与式(6)一样。
[0105]
(4-2)～(4-6)参考(三)中步骤，使得窄巷相位小数偏差估计和宽巷相位小数偏差估计保持一致，便于程序扩展使用，并输出此时的卫星窄巷相位小数偏差
[0106]
(五)线性转换各频相位小数偏差；
[0107]
利用步骤(三)和步骤(四)获取的宽巷和窄巷相位小数偏差，依据(1)中mw组合和lc组合其对应频率的系数，按照式(9)和(10)将其转向至各频相位小数偏差并标记其伪距观测值通道，使得用户可以采用灵活的函数模型实现广域单测站模糊度固定；
[0108][0109][0110]
式中，和分别为转化后的第一频和第二频相位小数偏差，即实时北斗相位小数偏差，dcb即为相应伪距观测值通道上的卫星伪距差分码偏差。
[0111]
综上，本发明采用严格的数据质量控制策略，基于并行计算和单观测值迭代更新算法，实现了多测站网解时变相位小数偏差快速估计，既提高了实时北斗相位小数偏差估计的效率，又提高了估计的相位小数偏差的精度和可靠性。