本发明涉及一种GNSS接收定位领域，尤其涉及一种GNSS多模单频RTK周跳探测方法及装置。
背景技术：
：GNSS全称GlobalNavigationSatelliteSystem，即全球卫星导航系统的统称，包括美国的GPS全球卫星导航系统，俄罗斯的GLONASS全球卫星导航系统、中国的BDS北斗全球卫星导航系统(在建)以及欧洲的Galileo全球卫星导航系统(在建)。GNSS的快速发展，卫星的可用数快速增加，多GNSS系统组合定位成为当前以及未来GNSS定位的首选，其精度和可靠性较之单系统而言都得到了提高。要实现GNSS载波相位高精度定位要解决两个关键问题，即整周模糊度固定和载波相位周跳探测问题，而周跳是否准确探测又影响着整周模糊度能否固定或者快速收敛，因此周跳的探测尤为重要。当前适用于RTK(Real-timeKinematic)定位的周跳探测方法按照其检测量主要可划分为两类，一类基于观测值域，一类基于残差域。基于观测值域进行周跳探测方法主要利用观测值在时序上是平滑的这一特性构造检测量，典型的方法有伪距相位组合法、高次差法、多项式拟合法、多普勒积分法、电离层残差法、TurboEdit法等。基于残差域探测周跳的方法主要考虑到周跳在定位模型中引入了至少分米级的粗差这一特点，常用的方法有动态三差法、站际历元二次差法、滤波法等。伪距相位组合法的探测精度依赖于伪距的观测精度，该方法只能探测大的周跳。高次差法的探测精度受到接收机时钟的稳定性、采样间隔和大气延迟等因素的影响，并且该方法要求参与差分的前几个历元的观测值不存在周跳，否则难以分辨当前历元是否发生了周跳。多项式拟合法的探测精度同高次差法，也要求参与拟合参数计算的前几个历元观测值不存在周跳。多普勒积分法的探测精度依赖于多普勒观测值的精度以及采样间隔，在实际应用中一般用该方法探测大周跳。电离层残差法的探测精度主要受电离层活跃度的影响，一般情况下可用其探测一周等小周跳，但是该方法需要双频数据，不适用于单频载波相位定位。TurboEdit方法可探测小周跳，但是需要双频观测数据，不适用于单频载波相位定位。动态三差法不依赖载体的运动状态，可以探测小周跳，但是由于采用星间差分，所以存在参考星的选取问题，参考星观测值不能发生周跳，否则会导致所有的观测值引入粗差。星间差分也造成三差观测量是相关的，这对模型观测值粗差的探测非常不利的，只能选用逐次剔除卫星的方式进行粗差探测。当多颗卫星观测值同时发生周跳时，其周跳探测的过程将变得非常复杂，计算量大。此外该方法需要至少含有4颗干净的卫星观测值才能探测周跳。站际历元二次差法是在三差法的基础上提出的，为了解决三差法主要缺陷，该方法仅仅采用了站际-历元二次差分观测值，即避免了星间差分。该方法不存在参考星的选取问题，且二次差观测值之间不相关，所以可直接采用经典的数据探测法进行粗差探测。该方法主要缺陷在于其要求至少含有4颗干净的卫星观测值才能探测周跳，并且该方法是针对单GNSS系统提出的。滤波法的探测精度主要依赖于其预推的观测值精度，而滤波器的预测精度依赖于其动力学模型是否与载体真实的运动状态相符，当载体机动比较强时，如果没有外在状态输入信息(如惯导信息)，滤波器的预测精度会下降，周跳探测的精度下降。技术实现要素：为了克服现有技术的不足，本发明的目的之一在于提供了一种GNSS多模单频RTK周跳探测方法，其要解决的是适用于多GNSS系统单频RTK周跳探测的问题。为了解决上述问题，本发明所采用的技术方案如下：一种GNSS多模单频RTK周跳探测方法，包括以下步骤：步骤1、与GPS卫星、GLONASS卫星、Galileo卫星和北斗卫星通信并获取相应数据；步骤2、根据公式(1)进行计算，得到各卫星的残差向量，VG=BGX+IGΔT-LGVR=BRX+IRΔT-LRVE=BEX+IEΔT-LEVC=BCX+ICΔT-LC---(1)]]>其中，V表示残差向量，B表示坐标参数对应的设计矩阵，I表示元素为1的列向量，L表示“观测值-计算值”向量，下标G、R、E、C依次表示GPS卫星、GLONASS卫星、Galileo卫星、北斗卫星，X表示[dx,dy,dz]T，ΔT表示站际历元相对钟差，dx＝xt2-xt1，dy＝xt2-xt1，dz＝zt2-zt1，t1和t2表示相邻的两个历元，(xtytzt)表示流动站在t时刻的坐标改正数；步骤3、根据VG、VR、VE、VC计算RMS值，若RMS值≥阈值EPS时，则判断存在卫星没有发生周跳，并执行步骤6，否则，执行步骤4；步骤4、根据公式(2)计算标准化残差vi‾=viδQvivi---(2)]]>其中，表示残差斜因数阵Qvv对角线上的第i个元素，Qvv＝P-1-B(BTPB)-1BT，P为观测方程权阵，δ0表示单位权中误差，vi表示VG、VR、VE、VC中其中一颗卫星的残差；步骤5、当时，将最大的对应的观测量剔除重新构建公式(1)进行平差计算，然后进入步骤4进行标准残差判断，直到所有的满足时，则判断存在载波观测值跳变的卫星已经全部被探测出来，执行步骤6，其中，uα/2为预设值；步骤6、将步骤5所判断到的周跳与利用多普勒积分法探测到的周跳进行融合。优选地，多普勒积分法的数据来源于步骤1。优选地，阈值EPS＝0.02m。优选地，uα/2＝2.0。优选地，所述观测方程权阵由卫星高度角定权法得到。本发明的另一个目的在于提供了一种应用于GNSS多模单频RTK周跳探测方法的装置，其要解决的是适用于多GNSS系统单频RTK周跳探测的问题。为了解决上述问题，本发明所采用的技术方案如下：一种应用于GNSS多模单频RTK周跳探测方法的装置，包括以下模块：信息获取模块：与GPS卫星、GLONASS卫星、Galileo卫星和北斗卫星通信并获取相应数据；第一计算处理模块：根据公式(1)进行计算，得到各卫星的残差向量，VG=BGX+IGΔT-LGVR=BRX+IRΔT-LRVE=BEX+IEΔT-LEVC=BCX+ICΔT-LC---(1)]]>其中，V表示残差向量，B表示坐标参数对应的设计矩阵，I表示元素为1的列向量，L表示“观测值-计算值”向量，下标G、R、E、C依次表示GPS卫星、GLONASS卫星、Galileo卫星、北斗卫星，X表示[dx,dy,dz]T，ΔT表示站际历元相对钟差，dx＝xt2-xt1，dy＝xt2-xt1，dz＝zt2-zt1，t1和t2表示相邻的两个历元，(xtytzt)表示流动站在t时刻的坐标改正数；第一判断模块：根据VG、VR、VE、VC计算RMS值，若RMS值≥阈值EPS时，则判断存在卫星没有发生周跳，并执行融合计算模块，否则，执行第二计算模块；第二计算处理模块：根据公式(2)计算标准化残差vi‾=viδQvivi---(2)]]>其中，表示残差斜因数阵Qvv对角线上的第i个元素，Qvv＝P-1-B(BTPB)-1BT，P为观测方程权阵，δ0表示单位权中误差，vi表示VG、VR、VE、VC中其中一颗卫星的残差；第二判断模块：当时，将最大的对应的观测量剔除重新构建公式(1)进行平差计算，然后进入步骤4进行标准残差判断，直到所有的满足时，则判断存在载波观测值跳变的卫星已经被全部探测出来，执行步骤6，其中，uα/2为预设值；融合计算模块：将第二判断模块所判断到的周跳与利用多普勒积分法探测到的周跳进行融合。优选地，多普勒积分法的数据来源于信息获取模块。优选地，阈值EPS＝0.02m。优选地，uα/2＝2.0。优选地，所述观测方程权阵由卫星高度角定权法得到。与现有技术相比，本发明的有益效果有：本发明提出了一种多模站际历元二次差法与多普勒积分法组合探测周跳的新方法。基于残差域探测周跳的多模站际历元二次差法，较之原来的站际历元二次差法，只设置了一个站际历元相对钟差参数，即在站际历元二次差观测值中忽略系统时间偏差影响，各GNSS的站际历元二次差观测值可认为处于同一时间系统中，参数个数减少，站际历元二次差法探测粗差需要的有效卫星数减少，并且所有的观测卫星都可以得到利用，由于不考虑站际历元二次差观测值中各系统的时间偏差，只有一个时间参数，所有的观测卫星等价于在同一系统中，即便单一系统只有1颗卫星，也可以有效利用和进行周跳探测。附图说明：图1为本发明的GNSS多模单频RTK周跳探测方法的流程图；图2为本发明的多模站际历元二次差法与多普勒积分法组合探测周跳流程示意图。具体实施方式下面，结合附图以及具体实施方式，对本发明做进一步描述：在描述本实施例之前，需要介绍任意一颗卫星的残差向量的公式推导：建立站际历元二次差方程，线性化后的载波相位站际单差观测方程可表示如下：式中：Δ表示站际单差算子；下标t表示历元，上标p表示卫星；表示载波相位观测值，λ为载波波长；ρ为由流动站坐标近似值计算的站星距离；δion，δtrop依次表示电离层，对流层延迟量；cδt表示钟差，N为模糊度参数；ε表示载波相位测量噪声。lx,tp=(Xt0-Xtp)/ρtp,ly,tp=(Yt0-Ytp)/ρtp,lz,tp=(Zt0-Ztp)/ρtp,]]>其中，Xt0Yt0Zt0]]>为流动站在t时刻的近似坐标，XtpYtpZtp]]>表示卫星p在t时刻的空间坐标；(xtytzt)表示流动站在t时刻的坐标改正数。对相邻两个历元t1，t2的站际单差观测方程进行二次差分计算，整理有：vp=(lx,t2pxt2-lx,t1pxt1)+(ly,t2pyt2-ly,t1pyt1)+(lz,t2pzt2-lz,t1pzt1)+ΔT-lp---(4)]]>式中ΔT＝Δcδtt2-Δcδtt1，对于短基线而言，站际差分消除了卫星钟差影响，极大削弱了卫星轨道误差影响，经站际历元二次差分后的电离层延迟、对流层延迟可认为基本消除。令dy＝xt2-xt1，dz＝zt2-zt1，则式(4)可表示为：vp=lx,t2pdx+ly,t2pdy+lz,t2pdz+ΔT-lp+γt1,t2p---(5)]]>由式(5)可知，若能将划入残差，上述方程只含有4个未知参数，方程求解较之式(4)减少3个未知参数，方程求解需要的卫星数减少。对式(5)中舍去γt1,t2项的舍入误差的进行分析。令dlxp=lx,t2p-lx,t1p,]]>则有dlxp=(Xt10-Xt1p)/ρt1p-(Xt20-Xt2p)/ρt2p---(6)]]>记dX0=X120-Xt10,dXp=Xt2p-Xt1p,dρp=ρt2p-ρt1p,]]>则式(6)可转化为dlxp=dρ·(Xt10-Xt1p)/(ρt1pρt2p)-(dX0-dXp)/ρt2p---(7)]]>同理可得：dlyp=dρ·(Yt10-Yt1p)/(ρt1pρt2p)-(dY0-dYp)/ρt2p---(8)]]>dlzp=dρ·(Zt10-Zt1p)/(ρt1pρt2p)-(dZ0-dZp)/ρt2p---(9)]]>定义同dY0，dZ0定义同dX0；dYp，dZp定义同dXp。由定义可知，(dX0dY0dZ0)的物理意义表示流动站从t1时刻到t2时刻坐标移动量的近似值。(dXpdYpdZp)的物理意义表示卫星p从t1时刻到t2时刻的坐标移动量。dρp表示流动站与卫星p从t1时刻到t2时刻站星几何距离的变化量。则有：γt1,t2=dlxpxt1+dlypyt1+dlzpzt1---(10)]]>记ΔRt1＝max(|xt1||yt1||zt1|)，其物理意义表示流动站在t1时刻的近似坐标的精度。|γt1,t2|≤(|dlxp|+|dlyp|+|dlzp|)·ΔRt2---(11)]]>将式(7)、(8)、(9)代入式(11)，整理有：|γt1,t2|≤|dρp/ρt2p|·(|Xt10-Xt1p|+|Yt10-Yt1p|+|Zt10-Zt1p|)/ρt1p·ΔRt2+(|dX0|+|dY0|+|dZ0|+|dXp|+|dYp|+|dZp|)/ρt2p·ΔRt2---(12)]]>定义：dSrov＝[(dX0)2+(dY0)2+(dZ0)2]1/2，dSp＝[(dXp)2+(dYp)2+(dZp)2]1/2。其物理意义依次是流动站和卫星从t1时刻到t2时刻移动的距离。由柯西不等式可知：由式(12)可得：|γt1,t2|≤3(|dρp|+dSrov+dSp)/ρt2p·ΔRt2---(13)]]>考虑到站星距离ρ一般大于20000km，站星距离变化率小于800m/s，GPS卫星运动速度小于3800m/s，当流动站运动速度小于250m/s时，项对站际历元间二次差方程的影响可由下式进行估计：|γt1,t2|≤4.2×10-4(t2-t1)ΔRt2(14)式中ΔRt2表示t2时刻流动站近似坐标的坐标分量精度。若先采用伪距差分定位计算流动站的近似坐标，一般可保证ΔRt2优于10m。在(t2-t1)分别取0.1s，1s，5s时，舍去项所造成的极限误差依次为0.42mm，4.2mm，21mm，远小于L1波长，实际情况中舍入误差远小于其极限误差，而在动态定位中，采样间隔一般比较小，因此，γt1,t2项的舍入误差对二次差分模型的影响可以忽略。由此，舍去γt1,t2后，式(5)可表示为：vp=lx,t2pdx+ly,t2pdy+lz,t2pdz+ΔT-lp---(15)]]>vp表示任意一颗卫星的残差向量，公式(15)即为本实施例的站际历元二次差方程。基于上述公式(15)结合图1所示，本实施例公开的一种GNSS多模单频RTK周跳探测方法，包括以下步骤：步骤1、与GPS卫星、GLONASS卫星、Galileo卫星和北斗卫星通信并获取相应数据。获取数据的方法为现有技术，在此不再赘述，获取流程如图2所示。步骤2、根据公式(1)进行计算，得到各卫星的残差向量，VG=BGX+IGΔT-LGVR=BRX+IRΔT-LRVE=BEX+IEΔT-LEVC=BCX+ICΔT-LC---(1)]]>其中，V表示残差向量，B表示坐标参数对应的设计矩阵，I表示元素为1的列向量，L表示“观测值-计算值”向量，下标G、R、E、C依次表示GPS卫星、GLONASS卫星、Galileo卫星、北斗卫星，X表示[dx,dy,dz]T，ΔT表示站际历元相对钟差。请结合公式(15)，以计算GPS卫星的残差向量为例进行说明，而且GPS卫星不仅仅是一颗，BGX即为多颗卫星观测方程项的组合形式，LG即为多颗卫星观测方程lp项的组合形式，其他类型的卫星也如此类推。步骤3、根据VG、VR、VE、VC计算RMS值，若RMS值≤阈值EPS时，则判断所有卫星没有发生周跳，并执行步骤6，否则，执行步骤4。站际历元二次差方程舍去|γt1,t2|项造成的误差在高采样率动态定位中远小于载波波长，所以可以采用平差后残差的RMS值作为探测周跳的指标，当RMS＞EPS时，认为存在卫星发生了周跳。借鉴文献4动态三差法阀值的取值方法，考虑到舍入误差远小于其极限误差，阀值EPS可以取0.02m。步骤4、根据公式(2)计算标准化残差vi‾=viδQvivi---(2)]]>其中，表示残差斜因数阵Qvv对角线上的第i个元素，Qvv＝P-1-B(BTPB)-1BT，P为观测方程权阵，δ0表示单位权中误差，vi表示VG、VR、VE、VC中其中一颗卫星的残差。步骤5、当时，将最大的对应的观测量剔除重新构建公式(1)进行平差计算，然后进入步骤4进行标准残差判断，直到所有的满足时，则判断存在载波观测值跳变的卫星已经全部被探测出来，执行步骤6，其中，uα/2为预设值。所述观测量是指站际历元二次差方程，即与最大的对应的卫星的站际历元二次差方程。上述步骤4、5的理论依据如下：Tiberius等人的研究表明GPS测量数据概率特性可以用正态概率密度函数描述[1]。记v，Δ依次表示残差和真误差，B为设计矩阵，P为观测方程权阵。基于最小二乘平差有[2-3]：-v＝QvvPΔ＝RΔ(16)式中Qvv＝P-1-B(BTPB)-1BT，R＝I-B(BTPB)-1BTP。式(16)表明v可由Δ线性表达，所以当Δ是偶然误差服从正态分布时，表示残差斜因数阵Qvv对角线上的第i个元素。标准化残差标准化残差由式(2)计算，式中δ0表示单位权中误差。vi‾=viδ0Qvivi---(2)]]>根据u检验法可知，当时，认为数据中存在粗差。由于最小二乘具有均摊作用，某一观测量中的粗差平差后会对其他观测方程残差造成影响，导致多个超限，此时只需将最大的对应的观测量剔除并进行粗差标记，然后重复“平差-检验-剔除粗差观测量”，直到所有的满足α一般取0.05，则uα/2为1.96，本实施例将其取整到2.0。步骤6、将步骤5所判断到的周跳与利用多普勒积分法探测到的周跳进行融合。多普勒积分法的周跳检测量构造公式为：式中，D表示多普勒观测值。当大于预设值时，则判断为发生了周跳。本实施例还公开了一种应用于GNSS多模单频RTK周跳探测方法的装置，包括以下模块：信息获取模块：与GPS卫星、GLONASS卫星、Galileo卫星和北斗卫星通信并获取相应数据；第一计算处理模块：根据公式(1)进行计算，得到各卫星的残差向量，VG=BGX+IGΔT-LGVR=BRX+IRΔT-LRVE=BEX+IEΔT-LEVC=BCX+ICΔT-LC---(1)]]>其中，V表示残差向量，B表示坐标参数对应的设计矩阵，I表示元素为1的列向量，L表示“观测值-计算值”向量，下标G、R、E、C依次表示GPS卫星、GLONASS卫星、Galileo卫星、北斗卫星，X表示[dx,dy,dz]T，ΔT表示站际历元相对钟差，dx＝xt2-xt1，dy＝xt2-xt1，dz＝zt2-zt1，t1和t2表示相邻的两个历元，(xtytzt)表示流动站在t时刻的坐标改正数；第一判断模块：根据VG、VR、VE、VC计算RMS值，若RMS值≥阈值EPS时，则则判断存在卫星没有发生周跳，并执行融合计算模块，否则，执行第二计算模块，其中EPS可以取0.02m。第二计算处理模块：根据公式(2)计算标准化残差vi‾=viδQvivi---(2)]]>其中，表示残差斜因数阵Qvv对角线上的第i个元素，Qvv＝P-1-B(BTPB)-1BT，P为观测方程权阵，δ0表示单位权中误差，vi表示VG、VR、VE、VC中其中一颗卫星的残差；第二判断模块：当时，将最大的对应的观测量剔除重新构建公式(1)进行平差计算，然后进入步骤4进行标准残差判断，直到所有的满足时，则判断存在载波观测值跳变的卫星已经被全部探测出来，执行步骤6，其中，uα/2为预设值，uα/2＝2.0；融合计算模块：将第二判断模块所判断到的周跳与利用多普勒积分法探测到的周跳进行融合。本发明所涉及的参考文献如下：[1]TiberiusC,BorreK.AreGPSdatanormallydistributed[M]//GeodesyBeyond2000.SpringerBerlinHeidelberg,2000:243-248。[2]於宗俦,李明峰.多维粗差的同时定位与定值[J].武汉测绘科技大学学报,1996,04:17-23。[3]宋力杰,杨元喜.均值漂移模型粗差探测法与LEGE法的比较[J].测绘学报,1999,04:295-300。[4]袁洪,万卫星,宁百齐,李静年.基于三差解检测与修复GPS载波相位周跳新方法[J].测绘学报,1998,03:3-8+87。对于本领域的技术人员来说，可根据以上描述的技术方案以及构思，做出其它各种相应的改变以及变形，而所有的这些改变以及变形都应该属于本发明权利要求的保护范围之内。当前第1页1 2 3