说明书一种基于抗差自适应Kalman滤波的多系统动态PPP解算方法
技术领域
本发明涉及一种多系统动态PPP解算方法，属于GNSS动态精密单点定位领域。
背景技术
卫星定位技术已成为现代最主要的导航定位方式，无论是在工程测量、生产生活还是军事应用中，均发挥着重要的作用。卫星定位的基本原理是测量出已知位置的卫星到用户接收机之间的距离，然后根据多颗卫星的瞬间位置采用距离后方交会的方法，确定待测点的位置坐标。依据卫星接收机设备的定位方式，GNSS定位技术可分为单点定位和差分定位。传统的卫星单点定位技术是指利用伪距码进行定位，即GNSS绝对定位，该技术能够采用一台接收机完成实时测量。在实际应用中，由于伪距码信号码长较大，其本身精度较低，并且通过结合广播星历计算的卫星位置和卫星钟差精度较差等因素，伪距单点定位难以达到较高的定位精度，一般用于精度要求较低，实时性要求较高的导航应用中。
利用载波相位观测值进行载波差分定位，能够达到较高的定位精度，但需要至少两台接收机进行同步观测，通过双差分模型才能实现厘米级的定位服务。随着接收机间距离的增加，误差相关性逐渐减小，到达一定距离将无法减弱误差影响。该定位技术不仅增加了作业的成本和复杂度，而且在很多应用场合也受到了限制。精密单点定位(Precise Point Positioning,PPP)的出现克服了差分定位的诸多缺陷。PPP由美国喷气推进实验室(JPL)的Zumberge等人于1997年提出，并在他们开发的数据处理软件GIPSY上给予实现。PPP具有厘米级的静态定位精度和分米级的动态定位精度，可彻底摆脱大范围、长距离测量对地面参考站的依赖，显著提高了作业效率，节约了用户成本。PPP近年来不仅是定位技术领域的研究热点，也在各种控制测量、工程测量、地形测量等工程应用领域得到了广泛的应用。
在动态PPP定位中，每个历元的接收机坐标在不断变化，定位环境复杂。利用多系统联合定位，能够大幅度增加可用卫星数目，增强卫星几何构型，进一步满足动态定位需求。对于做不规则运动的物体，其精确的函数模型和随机模型均难以得到，在传统Kalman滤波的基础上，采用抗差估计减弱观测残差较大的信息，并利用预测状态向量确定自适应因子，以得到可靠的动态定位结果。
发明内容
发明目的：为了克服现有技术中存在的不足，本发明提供一种基于抗差自适应Kalman滤波的多系统动态PPP解算方法，能够有效解决目前动态定位中，由于函数模型和随机模型不准确，以及卫星空间结构不稳定对结果导致定位结果较差的现象。
技术方案：为实现上述目的，本发明采用的技术方案为：一种基于抗差自适应Kalman滤波的多系统动态PPP解算方法，首先使用三系统广播星历计算各自卫星坐标和卫星钟差，对其进行时空基准统一；结合观测文件中的伪距观测值组合所得的无电离层伪距观测值，进行选权迭代伪距单点定位，反算历元时刻接收机概略坐标和各系统接收机钟差；从IGS等分析中心网站获取各系统的卫星精密星历和精密钟差产品，通过拉格朗日插值方法计算各系统卫星精密坐标和卫星精密钟差；随后进行各误差改正模型计算定位误差改正值，观测文件、卫星精密坐标和卫星精密钟差并结合观测文件进行观测数据的质量控制；最后建立多系统动态PPP定位模型，采用抗差自适应Kalman滤波解算出观测值信息的残差值和预测信息，确定观测权阵和自适应因子，实现高精度、高稳定性的动态PPP定位。
具体包括如下步骤：
步骤1，使用三系统广播星历计算各自卫星坐标和卫星钟差，对其进行时空基准统一，得到统一基准后的卫星坐标和卫星钟差；根据观测文件的P1(C1)和P2(C2)伪距观测值，进行无电离层组合得到无电离层组合伪距观测值；根据观测文件的L1和L2观测值，进行无电离层组合得到无电离层组合载波观测值；
步骤2，根据步骤1得到的无电离层组合伪距观测值以及统一基准后的卫星坐标和卫星钟差进行选权迭代伪距单点定位，得到接收机概略坐标和接收机钟差；
步骤3，通过网络从IGS等分析中心获取三系统卫星精密星历和卫星精密钟差，根据观测文件的历元时刻，进行拉格朗日插值得到相应历元时刻的卫星精密坐标和卫星精密钟差；
步骤4，根据观测文件提供的接收机信息、步骤2得到的接收机概略坐标和接收机钟差以及步骤3得到的卫星精密坐标和精密钟差，通过误差改正模型计算精密单点定位过程的各项误差改正值，并结合观测数据进行观测数据数据的质量控制；
步骤5，根据步骤1得到的无电离层组合伪距观测值和无电离层组合载波观测值，步骤2得到的接收机概略坐标和接收机钟差，步骤3得到的卫星精密坐标和卫 星精密钟差以及步骤4得到的计算所得各项误差改正值，结合精密单点定位模型，构建多系统动态精密单点定位方程；
步骤6，采用抗差自适应卡尔曼滤波解算步骤5得到的多系统动态精密单点定位方程，根据滤波过程中得到的观测信息和预测的状态信息，调整观测值的权值和计算自适应因子，实现高精度和高稳定性的动态定位。
所述步骤1中无电离层组合得到无电离层组合伪距观测值和无电离层组合载波观测值的公式：
PIF=f12·P1-f22·P2f12-f22=ρ+cli·(dtr-dts)+dtrop+ϵPIF]]>
LIF=f12·L1-f22·L2f12-f22=ρ+cli·(dtr-dts)+dtrop+f12·λ1·N1-f22·λ2·N2f12-f22+ϵLIF]]>
其中，PIF表示无电离层组合伪距观测值；f1,f2表示GPS/GLONASS/BDS系统的载波观测频率；P1,P2表示GPS/GLONASS/BDS系统的伪距观测值；若没有P1/P2观测值时，可采用C1/C2观测值进行DCB改正后使用；ρ是卫星和测站的距离；cli表示光速；dtr为表示接收机钟差；dts表示卫星钟差；dtrop是对流层延迟；LIF表示无电离层组合载波观测值；L1,L2表示GPS/GLONASS/BDS系统的载波相位观测值；λ1，λ2分别表示GPS/GLONASS/BDS系统的载波相位波长，N1，N2分别表示GPS/GLONASS/BDS系统的无电离层模糊度；分别表示伪距和载波无电离层组合观测噪声。
所述步骤2得到的卫星精密坐标和卫星精密钟差：
Pn(x)=Σi=1nli(x)·yi=Σi=0n(yi·Πj=0,j≠inx-xjxi-xj)]]>
式中，Pn(x)表示待插值时刻的卫星信息，即插值后每个历元的卫星轨道坐标和钟差值，x为插值时刻，xi(i＝0,1,…,n)为节点对应时刻；yi(i＝0,1,…,n)为节点时刻对应的卫星轨道坐标和钟差值；li(x)(i＝0,1,…,n)表示为n次插值的的基本差值多项式。
所述步骤5得到的多系统动态精密单点定位观测方程为：
PIFg=ρg+cli·dtg+dtropg+ϵPIFg]]>
ΦIFg=ρg+cli·dtg+dtropg+NIFg+ϵΦIF′g]]>
PIFr=ρr+cli·dtr+dtropr+ϵPIFr]]>
ΦIFr=ρr+cli·dtr+dtropr+NIFr+ϵΦIFr]]>
PIFc=ρc+cli·dtc+dtropg+ϵPIFc]]>
ΦIFc=ρc+cli·dtc+dtropg+NIFc+ϵΦIFc]]>
式中上标g、r、c分别表示GPS、GLONASS和BDS系统；cli是光速；PIF是误差改正后的无电离层组合伪距观测值；ΦIF是误差改正后的无电离层组合载波观测值；ρ是卫星和测站的距离(通过接收机概略坐标和卫星坐标求得)；dt是各系统接收机钟差；dtrop是对流层延迟；NIF是无电离层组合模糊度；εIF是多路径误差，观测值噪声和其他残留的误差。定位观测方程中卫星精密钟差同其余解算所得定位误差改正值已做消除。
所述步骤5中多系统动态精密单点定位观测方程线性化后，可表示为：
V＝AX-L
上式中，观测设计矩阵A可表示为：
A=ax1GPSay1GPSaz1GPS100M1GPS...........................ax2nGPSay2nGPSaz2nGPS100M2nGPS...ax1GLOay1GLOaz1GLO010M1GLO...........................ax2jGPSay2jGPSaz2jGPS010M2jGLO...ax1BDSay1BDSaz1BDS001M1BDS...........................ax2kBDSaz2kBDSaz2kBDS001M2kBDS...]]>
状态向量可表示为：
X=δxδyδzcli·dtgcli·dtrcli·dtcdtropNIFg1...NIFgnNIFr1...NIFrjNIFc1...NIFck]]>
式中：V表示观测值残差向量；A表示观测设计矩阵；X表示状态向量；L表示观测向量；上标g、r、c表示GPS、GLONASS、BDS三系统；n,j,k表示观测到的GPS，GLONASS，BDS卫星个数；A阵前三列为观测方程线性化后δx,δy,δz的系数；A阵第七列为参数dtrop的系数；A阵第七列后的其他列为模糊度参数的系数，系数设置为1；δx,δy,δz分别表示x,y,z方向卫地距残差；cli表示光速；dt表示三系统相对应的接收机钟差；dtrop表示对流层延迟；NIF表示无电离层整周模糊度。
所述步骤6中实现高精度和高稳定性的动态定位方法，包括以下步骤：
步骤6.1，根据PPP观测方程中得到的相关参数，建立动态Kalman滤波模型，并设置滤波参数的初始值；Kalman滤波的函数模型为：
Xk+1=ΦkXk+WkLk+1=Ak+1Xk+1+Vk+1]]>
式中：下标k表示历元序号；Xk表示tk时刻观测方程中的状态向量；Xk+1表示tk+1时 刻观测方程中的状态向量；Φk表示状态转移矩阵，一般情况下取单位矩阵；Wk表示状态模型输入噪声向量；Lk+1表示观测向量，分为伪距和载波观测值；Ak+1表示观测方程中的系数矩阵；Vk+1表示观测值残差向量。
Kalman滤波的统计模型为：
E(Wk)=0→,E(Vk+1)=0→Cov(Wk)=Qk,Cov(Vk+1)=Rk+1,Cov(Wk,Vk+1)=0]]>
式中：E(·)表示求取数学期望；Cov(·)表示求取协方差阵；Wk表示状态模型输入噪声向量；Vk+1表示观测值残差向量；Qk表示状态模型输入噪声向量的协方差阵；Rk+1表示观测值残差向量的协方差阵；
步骤6.2，根据卫星高度角及系统间的经验权值比设置初始随机模型，再利用观测值的残差向量信息，采用抗差估计模型重新确定权值：
p‾i=pi,|v~i|≤k0pik0|v~i|(k1-|v~i|k1-k0)2,k0<|v~i|≤k10,|v~i|>k1]]>
式中：下标i表示残差值个数；pi表示先验权值，由高度角定权公式得到；表示观测值的残差值，即观测值残差向量Vk中的元素；k0,k1表示抗差估计阈值，利用残差值的总体分布确定；表示抗差估计调整后的权值；
根据抗差估计后确定的权值，则可以得到第k个历元的状态参数抗差解为：
X~k=(AkTP‾kAk)-1AkTP‾kLk]]>
式中：下标k表示历元序号；上标T表示矩阵转置；表示抗差状态向量；表示抗差估计调整后的权阵；Lk表示观测向量；Ak表示观测方程中的系数矩阵。
步骤6.3，自适应Kalman滤波中，需要根据滤波过程中得到的当前观测信息和预测状态信息，确定自适应因子的取值为：
αk=1,|ΔX~k|≤c0c0|ΔX~k|(c1-|ΔX~k|c1-c0)2,c0<|ΔX~k|≤c10,|ΔX~k|>c1]]>
式中：
ΔX~k=||X~k-X‾k||/tr(ΣX‾k)]]>
式中：下标k表示历元序号；||·||表示求范数；表示抗差状态向量；表示状态预测信息向量；c0,c1表示自适应因子的取值阈值；表示状态预报向量协方差矩阵的迹；αk表示自适应因子；
步骤6.4，根据上面步骤得到的数据，更新观测信息：
Kk=1αkΣX‾kAkT(1αkAkΣX‾kAkT+P‾k-1)-1X^k=X‾k+Kk(Lk-AkX‾k)ΣX^k=(I-KkAk)ΣX‾k/αk]]>
式中：上标T表示矩阵转置；Kk表示增益矩阵；αk表示自适应因子；表示状态预报向量协方差矩阵；Ak表示观测方程中的系数矩阵；表示抗差估计调整后的权阵；表示状态估计向量；表示状态预报向量；Lk表示观测向量；表示状态估计向量协方差矩阵；I表示单位矩阵；
步骤6.5，利用步骤6.2和步骤6.3中确定的权值和自适应因子，在步骤6.4中更新观测信息，进行滤波解算，从而实现减弱粗差和动态噪声对状态参数估值的影响。
有益效果：本发明提供的一种基于抗差自适应Kalman滤波的多系统动态PPP解算方法，在可用卫星数量足够充足的基础上，充分利用观测值的残差信息重新确定随机模 型，剔除动态测试中存在粗差的卫星；依据Kalman滤波得到的预测状态量的信息，确定出自适应因子，减小误差较大的预测信息的影响。且采用伪距单点定位单历元解算接收机概略位置和钟差，不会在历元间传递和积累定位误差，有效避免了由于局部观测误差导致定位结果发散的现象。与常规的动态PPP解算方法相比，本专利所提出的方法可以大幅度增加可用卫星数，有效减弱粗差的影响，并改善动态定位中的动态模型噪声异常情况，最终达到高精度和高稳定性的动态定位结果。
附图说明
图1是基于抗差自适应Kalman滤波的多系统动态PPP解算方法流程图；
图2是实验所用天津某CORS基准站分布图；
图3是基准站1：BC站GPS、GLONASS、BDS系统可用卫星数；
图4是基准站2：NH站GPS、GLONASS、BDS系统可用卫星数；
图5是基准站1：BC站GPS、GLONASS、BDS三系统动态PPP定位结果；
图6是基准站2：NH站GPS、GLONASS、BDS三系统动态PPP定位结果；
图7是动态测试数据：基于抗差自适应Kalman滤波的多系统动态PPP平面结果；
图8是动态测试数据：基于抗差自适应Kalman滤波的多系统动态PPP高程结果；
图9是动态测试数据：运动过程中基于抗差自适应和一般Kalman滤波的多系统动态PPP高程结果。
具体实施方式
下面结合附图对本发明作更进一步的说明。
一种基于抗差自适应Kalman滤波的多系统动态PPP解算方法，如图1所示：首先使用三系统广播星历计算各自卫星坐标和卫星钟差，对其进行时空基准统一；结合观测文件中的伪距观测值组合所得的无电离层伪距观测值，进行选权迭代伪距单点定位，反算历元接收机概略坐标和各系统接收机钟差；从IGS等分析中心网站获取各系统的卫星精密星历和精密钟差产品，通过拉格朗日插值方法计算各系统卫星精密坐标和卫星精密钟差；随后进行各误差改正模型计算定位误差改正值，观测文件、卫星精密坐标和卫星精密钟差并结合观测数据进行观测数据的质量控制；最后建立多系统动态PPP定位模型，采用抗差自适应Kalman滤波解算出观测值信息的残差值和预测信息，确定观测权阵和自适应因子，实现高精度、高稳定性的多系统动态PPP定位。
具体包括以下步骤：
步骤1，使用三系统广播星历计算各自卫星坐标和卫星钟差，对其进行时空基准统一，得到统一基准后的卫星坐标和卫星钟差；根据观测文件的P1(C1)和P2(C2)伪距观测值，进行无电离层组合得到无电离层组合伪距观测值；根据观测文件的L1和L2观测值，进行无电离层组合得到无电离层组合载波观测值。
1)，每个导航定位系统都拥有各自的时间系统，通过各系统时间基准之间的差异，将各系统时间基准统一到GPST，从而消除各系统间的时间偏差。GPS、GLONASS、BDS的时间系统信息列于表1：
表1 GPS、GLONASS、BDS时间系统
结合GPST、GLONASST和BDT与UTC时间之间的关系，三者之间的时间转化关系可以表示为：
GPST＝GLONASST+1s×n-19s-3h (1.1)
GPST＝BDT+14s (1.2)
式中，GPST表示GPS时；GLONASST表示GLONASS时；n表示跳秒数；BDST表示北斗时。
由于各个系统的接收机钟差值也存在偏差，在进行选权迭代伪距单点定位时，需将各系统对应的接收机钟差作为待估参数进行解算；解算所得各系统接收机钟差代入后续精密单点定位滤波方程中作为初值参与解算。
对于GPS、GLONASS、BDS三个坐标系之间的转换，需要通过坐标原点的平移、坐标轴的旋转以及坐标轴刻度的单位的变换实现。两个任意的三维空间直角坐标系O-XYZ和O′-X′Y′Z′，当它们存在三个及其以上的已知点时，可以采用布尔莎七参数模型计算转换得到：
XYZ=δxδyδz+(1+m)1ϵz-ϵy-ϵz1ϵxϵy-ϵx1X′Y′Z′---(1.3)]]>
式中，X、Y、Z和X′、Y′、Z′分别表示两个三维空间直角坐标系下的三维坐标；δx、δy、δz分别表示两个坐标系原点之间的平移参数；εx、εy、εz表示坐标系旋转产生的旋转参数；m表示尺度变化参数。
通过长期的计算与测试，已有很多机构和学者对GLONASS所采用的PZ-90坐标系和GPS的WGS-84坐标系进行了统一转换研究，目前，世界公认精度最高的经验转换公式为：
XYZWGS-84=-0.47-0.51-1.56+(1+22×10-9)1-1.728×10-6-0.017×10-61.728×10-610.076×10-60.017×10-6-0.076×10-61UVWPZ-90---(1.4)]]>
式中，X、Y、Z表示WGS84坐标系下的三维坐标；U、V、W表示PZ-90坐标系下的三维坐标。
目前，对于BDS所采用的CGCS2000坐标系与WGS-84之间没有明确的转换公式，且因为两者之间的坐标系参数相差很小，因此BDS与GPS系统无需进行坐标转换操作。
2)，利用精密星历中的卫星轨道坐标和双频伪距观测值信息，通过“消电离层组合”模型，消除电离层延迟误差的一阶项，计算出米级精度的接收机概略位置和接收机钟差。
PIF=f12P1-f22P2f12-f22---(1.6)]]>
上式中，PIF表示无电离层组合伪距观测值；f1,f2表示GPS/GLONASS/BDS系统的载波观测频率；P1,P2表示GPS/GLONASS/BDS系统的伪距观测值；若没有P1/P2观测值时，可采用C1/C2观测值进行DCB改正后使用。在多系统伪距单点定位方程中，包括3个接收机位置参数、3个接收机钟差参数。
步骤2，根据步骤1得到的无电离层组合伪距观测值以及统一基准后的卫星坐标和卫星钟差进行选权迭代伪距单点定位，得到接收机概略坐标和接收机钟差。
由于精密单点定位时，只需采用一台接收机进行数据采集，无法像差分定位一样消除或减弱部分误差。因此，需精确处理各项误差改正值，根据接收机概略坐标和卫星精 确位置，将误差分为三类进行处理：与卫星相关、与信号传播相关以及与接收机相关。精密单点定位的传统模型如式(1.7)和(1.8)所示：
PIF=f12·P1-f22·P2f12-f22=ρ+cli·(dtr-dts)+dtrop+ϵPIF---(1.7)]]>
LIF=f12·L1-f22·L2f12-f22=ρ+cli·(dtr-dts)+dtrop+f12·λ1·N1-f22·λ2·N2f12-f22+ϵLIF---(1.8)]]>
上式中，其中，PIF表示无电离层组合伪距观测值；f1,f2表示GPS/GLONASS/BDS系统的载波观测频率；P1,P2表示GPS/GLONASS/BDS系统的伪距观测值；若没有P1/P2观测值时，可采用C1/C2观测值进行DCB改正后使用；ρ是卫星和测站的距离；cli表示光速；dtr为表示接收机钟差；dts表示卫星钟差；dtrop是对流层延迟；LIF表示无电离层组合载波观测值；L1,L2表示GPS/GLONASS/BDS系统的载波相位观测值；N1，N2分别表示GPS/GLONASS/BDS系统的无电离层模糊度；λ1,λ,2表示GPS/GLONASS/BDS系统的载波相位波长；分别表示伪距和载波的无电离层组合观测噪声。若观测到n颗卫星，则观测方程数为2n个，未知数有(7+n)个，包括3个位置参数、3个接收机钟差、1个对流层湿延迟误差以及n个模糊度参数。
步骤3，通过网络从IGS等分析中心获取三系统卫星精密星历和卫星精密钟差，根据观测文件的历元时间，进行拉格朗日插值得到相应时刻的卫星精密坐标和卫星精密钟差。
在IGS网站上下载精密星历和精密钟差文件，利用拉格朗日插值得到高采样率的卫星信息，如式(1.1)所示：
Pn(x)=Σi=0nli(x)·yi=Σi=0n(yi·Πj=0,j≠inx-xjxi-xj)---(1.5)]]>
式中，表示待插值时刻的卫星信息，即插值后每个历元的卫星轨道坐标和钟差值， 为插值时刻，为节点对应时刻；为节点时刻对应的卫星轨道坐标和钟差值；表示为n次插值的的基本差值多项式；Pn(x)表示待插值时刻的卫星信息，即插值后每个历元的卫星轨道坐标和钟差值。
采用拉格朗日插值多项式时，应选取合适的阶数，此外，拉格朗日插值法外推精度较低，故针对单天24小时观测文件，需结合前后两天的卫星精密星历和卫星精密钟差文件进行内插。
步骤4，根据观测文件提供的接收机信息、步骤2得到的接收机概略坐标和接收机钟差以及步骤3得到的卫星精密坐标和精密钟差，通过误差改正模型计算精密单点定位过程的各项误差改正值，并结合观测数据进行观测数据数据的质量控制。
步骤5，根据步骤1得到的无电离层组合伪距观测值和无电离层组合载波观测值，步骤2得到的接收机概略坐标和接收机钟差，步骤3得到的卫星精密坐标和卫星精密钟差以及步骤4得到的计算所得各项误差改正值，结合精密单点定位模型，构建多系统动态精密单点定位方程。
采用卫星精密星历和卫星精密钟差产品后，多系统动态精密单点定位观测方程可表示为：
PIFg=ρg+cli·dtg+dtropg+ϵPIFg]]>
ΦIFg=ρg+cli·dtg+dtropg+NIFg+ϵΦIF′g]]>
PIFr=ρr+cli·dtr+dtropr+ϵPIFr]]>
ΦIFr=ρr+cli·dtr+dtropr+NIFr+ϵΦIFr]]>
PIFc=ρc+cli·dtc+dtropg+ϵPIFc]]>
ΦIFc=ρc+cli·dtc+dtropg+NIFc+ϵΦIFc]]>
式中上标g、r、c分别表示GPS、GLONASS和BDS系统；cli是光速；PIF是误差改正后的无电离层组合伪距观测值；ΦIF是误差改正后的无电离层组合载波观测值；ρ是卫星和测站的距离；dt是各系统接收机钟差；dtrop是对流层延迟；NIF是无电离层组合模糊度；εIF是多路径误差，观测值噪声和其他残留的误差。其中ρ通过接收机概略坐标和卫星坐标求得；卫星精密钟差同其余定位误差改正值一并消除。
所述步骤6中多系统动态精密单点定位观测方程线性化后，可表示为：
V＝AX-L
上式中，观测设计矩阵A可表示为：
A=ax1GPSay1GPSaz1GPS100M1GPS...........................ax2nGPSay2nGPSaz2nGPS100M2nGPS...ax1GLOay1GLOaz1GLO010M1GLO...........................ax2jGPSay2jGPSaz2jGPS010M2jGLO...ax1BDSay1BDSaz1BDS001M1BDS...........................ax2kBDSaz2kBDSaz2kBDS001M2kBDS...]]>
状态向量可表示为：
X=δxδyδzcli·dtgcli·dtrcli·dtcdtropNIFg1...NIFgnNIFr1...NIFrjNIFc1...NIFck]]>
式中：V表示观测值残差向量；A表示观测设计矩阵；X表示状态向量；L表示观测向量；上标g、r、c表示GPS、GLONASS、BDS三系统；n,j,k表示观测到的GPS，GLONASS，BDS卫星个数；A阵前三列为观测方程线性化后δx,δy,δz的系数；A阵第七列为参数dtrop的系数；A阵第七列后的其他列为模糊度参数的系数，系数设置为1；δx,δy,δz分别表示x,y,z方向卫地距残差；cli表示光速；dt表示三系统相对应的接收机钟差；dtrop表示对流层延迟；NIF表示无电离层整周模糊度。
对于任意一颗卫星i的系数阵为：
A=aibici1Mi10aibici1Mi00---(1.11)]]>
式中：xi、yi、zi表示卫星i的坐标，Xr、Yr、Zr表示接收机r的坐标，表示卫星到接收机的距离；Mi表示天顶方向上的对流层湿延 迟投影函数；第一行是载波观测方程式系数，第二行是伪距观测方程式系数。
步骤7，采用抗差自适应卡尔曼滤波解算步骤6得到的多系统动态精密单点定位方程，根据滤波过程中得到的观测信息和预测的状态信息，调整观测值的权值和计算自适应因子，实现高精度和高稳定性的动态定位。
在精密单点定位时通常采用Kalman滤波进行解算，可靠的Kalman滤波要求准确的函数模型和随机模型。在实际生活中，一般物体难以确保做规则的运动，因而构建精确的函数模型十分困难；随机模型一般根据已有的先验信息确定，通常采用高度角的正弦值来计算，与实际情况也存在一定偏差。于是可利用当前的观测信息和预测的状态信息来减弱运动中部分误差，即采用抗差模型和自适应滤波。
所述步骤6中实现高精度和高稳定性的动态定位方法，包括以下步骤：
步骤6.1，根据PPP观测方程中得到的相关参数，建立动态Kalman滤波模型，并设置滤波参数的初始值；Kalman滤波的函数模型为：
Xk+1=ΦkXk+WkLk+1=Ak+1Xk+1+Vk+1]]>
式中：下标k表示历元序号；Xk表示tk时刻观测方程中的状态向量；Xk+1表示tk+1时刻观测方程中的状态向量；Φk表示状态转移矩阵，一般情况下取单位矩阵；Wk表示状态模型输入噪声向量，在多系统动态PPP中坐标对应数值取经验数值100；Lk+1表示观测向量，分为伪距和载波观测值；Ak+1表示观测方程中的系数矩阵；Vk+1表示观测值残差向量。
Kalman滤波的统计模型为：
E(Wk)=0→,E(Vk+1)=0→Cov(Wk)=Qk,Cov(Vk+1)=Rk+1,Cov(Wk,Vk+1)=0]]>
式中：E(·)表示求取数学期望；Cov(·)表示求取协方差阵；Wk表示状态模型输入噪声向量；Vk+1表示观测值残差向量；Qk表示状态模型输入噪声向量的协方差阵；Rk+1表示观测值残差向量的协方差阵；
步骤6.2，根据卫星高度角及系统间的经验权值比设置初始随机模型，再利用观测值的残差向量信息，采用抗差估计模型重新确定权值：
p‾i=pi,|v~i|≤k0pik0|v~i|(k1-|v~i|k1-k0)2,k0<|v~i|≤k10,|v~i|>k1]]>
式中：下标i表示残差值个数；pi表示先验权值，由高度角定权公式得到；表示观测值的残差值，即观测值残差向量Vk中的元素；k0,k1表示抗差估计阈值，利用残差值的总体分布确定；表示抗差估计调整后的权值；
根据抗差估计后确定的权值，则可以得到第k个历元的状态参数抗差解为：
X~k=(AkTP‾kAk)-1AkTP‾kLk]]>
式中：下标k表示历元序号；上标T表示矩阵转置；表示抗差状态向量；表示抗差估计调整后的权阵；Lk表示观测向量；Ak表示观测方程中的系数矩阵。
步骤6.3，自适应Kalman滤波中，需要根据滤波过程中得到的当前观测信息和预测状态信息，确定自适应因子的取值为：
αk=1,|ΔX~k|≤c0c0|ΔX~k|(c1-|ΔX~k|c1-c0)2,c0<|ΔX~k|≤c10,|ΔX~k|>c1]]>
式中：
ΔX~k=||X~k-X‾k||/tr(ΣX‾k)]]>
式中：下标k表示历元序号；||·||表示求范数；表示抗差状态向量；表示状态预测信息向量；c0,c1表示自适应因子的取值阈值；表示状态预报向量协方差矩阵的迹；αk表示自适应因子；
步骤6.4，根据上面步骤得到的数据，更新观测信息：
Kk=1αkΣX‾kAkT(1αkAkΣX‾kAkT+P‾k-1)-1X^k=X‾k+Kk(Lk-AkX‾k)ΣX^k=(I-KkAk)ΣX‾k/αk]]>
式中：上标T表示矩阵转置；Kk表示增益矩阵；αk表示自适应因子；表示状态预报向量协方差矩阵；Ak表示观测方程中的系数矩阵；表示抗差估计调整后的权阵；表示状态估计向量；表示状态预报向量；Lk表示观测向量；表示状态估计向量协方差矩阵；I表示单位矩阵；
步骤6.5，利用步骤6.2和步骤6.3中确定的权值和自适应因子，在步骤6.4中更新观测信息，进行滤波解算，从而实现减弱粗差和动态噪声对状态参数估值的影响。
以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。