本发明属于卫星导航定位领域，涉及基于GNSS-RTK的定位方法。
背景技术：
GNSS：全球卫星导航系统，(Global Navigation Satellite System)，也称为全 球导航卫星系统，是能在地球表面或近地空间的任何地点为用户提供全天候的3 维坐标和速度以及时间信息的空基无线电导航定位系统。RTK(Real-time kinematic，实时动态)载波相位差分技术，是实时处理两个测量站载波相位观测 量的差分方法，将基准站采集的载波相位发给用户接收机，进行求差解算坐标。 随着人们对位置信息需求的不断增加，传统的高成本高精度定位模式已经无法满 足现实的需求，特别在建筑物、树木遮挡等复杂环境下定位时，存在精度、可靠 性差的问题。
技术实现要素：
为解决上述问题，本发明提出定位精度高的基于GNSS-RTK的定位方法。
基于GNSS-RTK的定位方法，包括以下步骤：
步骤1，根据解析的原始星地观测值进行数据质量控制；
步骤2，采用伪距观测值进行SPP定位，根据上一历元的RTD伪距残差和 当前历元SPP定位计算的伪距残差进行RTD定位；
步骤3，根据RTD定位结果作为初始值分别通过最小二乘、附加模糊度参数 的Kalman滤波计算浮点解；
步骤4，采用优化的部分模糊度固定策略将模糊度浮点解转换为固定解。
优选的，所述步骤1中根据解析的原始星地观测值进行数据质量控制包括以 下步骤：
通过观测数据预处理获取可靠的观测值；
利用数据探测法探测并改正粗差；
通过基于历元间二次差的周跳探测方法以及基于历元间差分的单站周跳探 测方法探测并剔除周跳。
优选的，所述通过观测数据预处理获取可靠的观测值包括：选取信噪比位于 设定范围内的观测值，通过伪距、载波精度信息和可用标识进行数据可用性判断， 并根据卫星连续跟踪、解算可用性选取观测数据。
优选的，所述利用数据探测法探测并改正粗差包括以下步骤：
构造线性化的伪距观测方程为：
V＝L-A·X
其中，V：残差向量；L：观测值减计算值；A：设计矩阵；X：参数向量；
构造标准化残差作为粗差检测量：
QVV＝P-1-A(ATPA)-1AT
其中，Vi：残差向量的元素；P：观测值的权阵；δ：单位权中误差；QViVi：残 差协因数矩阵QVV对角线元素；
如果标准化残差则存在粗差，标准化残差时，对权阵做 如下处理：
观测值减计算值L分为不包含粗差的l1和包含粗差向量S的l2，其维数分 别为n1和n2；则伪距观测方程变为：
设H＝[0E]T，E为n2×n2单位阵，计算粗差向量S并将其改正：
PSS＝HTPQVVPH。
优选的，所述通过基于历元间二次差的周跳探测方法以及基于历元间差 分的单站周跳探测方法探测并剔除周跳包括：先通过历元间二次差探测两周 以上的大周跳，然后通过基于历元间单差观测值的单站周跳探测方法探测与 剔除小周跳，所述通过历元间二次差探测两周以上的大周跳包括：
根据当前历元卫星的载波观测值与上一历元的载波观测值计算当 前历元的历元间载波观测值单差：
计算上一历元的历元间载波观测值单差：
计算当前历元的历元间载波观测值二次差：
当载波连续跟踪时，计算所有观测卫星的历元间载波观测值二次差
向量F中各元素集中分布，若某一元素超出数值 中心2周，则认为是大周跳。
优选的，所述步骤2中RTD定位包括以下步骤：
根据上一历元的RTD伪距残差和当前历元SPP定位计算的伪距残差计算得 到当前历元的观测值方差
其中，当前历元的观测值方差由以下模型进行计算：
其中，Vb：基准站SPP定位的伪距残差；Vr：流动站SPP定位的伪距残差，V0： 上一历元的RTD定位的伪距残差。
优选的，所述步骤3中根据RTD定位结果作为初始值分别通过最小二乘、 附加模糊度参数的Kalman滤波计算模糊度浮点解包括以下步骤:
构造双差观测方程：
其中，双差算子；Pi：伪距观测值；ρ：卫地距；I：电离层延迟；双差伪距观测噪声；φi：载波观测值；λi：载波波长；双差模糊度；
根据双差观测方程将流动站的位置、速度、加速度以及站间单差模糊度作为 待估参数，通过转换矩阵将站间单差模糊度转换成双差模糊度Kalman滤 波的参数矩阵和转换矩阵分别为：
Xk＝[x y z vx vy vz ax ay az △N1 … △Nn]T
其中，(x,y,z)：参考时刻的卫星位置；(vx,vy,vz)：参考时刻的卫星运 动速度；(ax,ay,az)：参考时刻的卫星运动加速度；
最小二乘法解算得到模糊度浮点解和验后残差σ，根据模糊度浮点解 和验后残差σ对Kalman滤波解算出的模糊度浮点解作如下变换：
优选的，所述步骤4中采用优化的部分模糊度固定策略将模糊度浮点解转换 为固定解，包括以下步骤：
根据浮点解方差对模糊度进行排序，模糊度固定按照浮点解精度顺序逐步进 行搜索；对双差模糊度参数和对应的方差－协方差矩阵进行降相关，得到相关性 较低的变换模糊度，降相关公式为：
其中，ZT：变换矩阵；模糊度浮点解；变换后的模糊度；模糊度 浮点解方差-协方差矩阵；变换后的模糊度浮点解方差-协方差矩阵，
矩阵元素全为整数，行列式的值为±1，变换后的搜索空间为：
其中，χ2：搜索空间；
搜索得到的整数解通过逆变换得到原始模糊度的整数解：
定义Ratio值对模糊度固定解进行检验：
其中，为次优解；为最优解；若满足上式则认为解算合格，否则保留 浮点解；
当连续若干个模糊度均固定，符合首次模糊度固定条件，保留模糊度固定解， 并设置固定解精度；对部分模糊度进行检验，若满足固定条件，将部分模糊度设 置为固定解。
通过使用本发明，可以实现以下效果：
(1)通过数据质量控制和残差处理，能够对原始数据进行有效改进，对数 据质量进行充分评估，为模糊度浮点解以及固定解提供可靠的原始观测数据和参 考依据，提高定位精度；
(2)通过最小二乘、附加模糊度参数的Kalman滤波计算模糊度浮点解，有 效提高收敛速度、防止滤波异常，提高浮点解的精度，进而缩短模糊度固定时间；
(3)充分利用解算信息，实现模糊度快速固定。
附图说明
下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1是本发明的整体流程示意图；
图2是本发明的具体流程示意图。
具体实施方式
以下结合附图，对本发明的技术方案作进一步的描述，但本发明并不限于这 些实施例。
一种基于GNSS-RTK的定位方法，其主要思想包括：
(1)数据预处理和质量控制。选取可靠的卫星观测值、发现和修复伪距粗差、 探测和剔除载波周跳，提供可靠的伪距、载波观测值，保证模糊度的连续性，降 低或者防止固定解失锁；
(2)实时更新随机模型。RTD(实时伪距相对定位)随机模型根据SPP(伪距 单点定位)残差和上一历元的验后残差确定，浮点解双差伪距随机模型根据RTD 残差确定，保证了随机模型的准确度，提高了定位结果的精度；
(3)最小二乘原理和Kalman滤波确定浮点解。由最小二乘解算的模糊度浮点 解对Kalman滤波解进行改进，有效提高浮点解精度，提高滤波收敛效率，为进 一步固定模糊度做准备；
(4)模糊度固定。将模糊度浮点解按方差-协方差进行排序，采用逐步固定、 部分模糊度固定和保留策略，有效提高模糊度固定效率，事后滤波提高了浮点解 的精度，快速可靠地实现高精度定位信息，满足位置需求。
结合附图1，本实施例包括以下步骤：
步骤1，根据解析的原始星地观测值进行数据质量控制；
步骤2，采用伪距观测值进行SPP定位，根据上一历元的RTD伪距残差和 当前历元SPP定位计算的伪距残差进行RTD定位；
步骤3，根据RTD定位结果作为初始值分别通过最小二乘、附加模糊度参数 的Kalman滤波计算浮点解；
步骤4，采用优化的部分模糊度固定策略将模糊度浮点解转换为固定解。
结合附图2，下面对各个步骤进行详细说明：
步骤1，根据解析的原始星地观测值进行数据质量控制。
基准站、流动站原始观测数据一般采用自定义协议或者标准协议进行编译， 通过对基准站原始观测数据和流动站原始观测数据解析得到观测值和卫星星历。 然而观测值和卫星星历存在误差较大值，同时一些观测值和卫星星历存在周跳等。
基准站差分数据一般采用RTCM标准协议进行编译。RTCM帧结构由同步码、 保留位、信息长度、数据信息和CRC组成，对于RTD需播发包含基准站位置信 息、GNSS各系统观测值的定义信息。流动站原始观测数据一般采用自定义的二 进制协议，需包含卫星观测数据以及数据精度指标。
卫星轨道信息通过历书和当地观测时间进行计算。对于CDMA(码分多址) 卫星定位系统，广播星历预报参数包括：
①时间参数：toe：星历参考时刻(秒)；toc：钟差参考时刻(秒)；IODE：星 历数据龄期；
②开普勒六参数：M0：toe时刻的平近点角(弧度)；e：轨道偏心率；长 半轴的平方根Ω0：按参考时间计算的升交点赤经(弧度)；i0：按参 考时间计算的轨道倾角(弧度)；ω：近地点角距(弧度)；
③轨道摄动九参数：△n：卫星平均运动速率与计算值之差(弧度)；升交 点赤经变化率(弧度/秒)；IDOT：轨道倾角变化率(弧度/秒)；Cuc：纬度幅 角的余弦调和改正项的振幅(弧度)；Cus：纬度幅角的正弦调和改正项的振幅 (弧度)；Crc：轨道半径的余弦调和改正项的振幅(弧度)；Crs：轨道半径的 正弦调和改正项的振幅(弧度)；Cic：轨道倾角的余弦调和改正项的振幅(弧 度)；Cis：轨道倾角的正弦调和改正项的振幅(弧度)。
对于FDMA(频分多址)卫星定位系统，广播星历预报参数包括：
①时间参数：toe：星历参考时刻(秒)；IODE：星历数据龄期；
②位置参数：(x,y,z)：参考时刻的卫星位置；(vx,vy,vz)：参考时刻的卫 星运动速度；(ax,ay,az)：参考时刻的卫星运动加速度。
具体的，步骤1中根据解析的原始星地观测值进行数据质量控制包括以 下步骤：
步骤1.1，观测数据预处理：通过观测数据预处理获取可靠的观测值。
首先根据观测数据的信噪比，选取信噪比位于合理范围的观测值；然后 根据数据可用标识，剔除不可用的数据；再利用数据输出的伪距、载波精度 信息，选取精度较高的数据；最后，根据卫星跟踪连续情况确定观测值是否 参与解算情况。
步骤1.2，伪距粗差探测与修复：利用数据探测法探测并改正粗差。
构造线性化的伪距观测方程为：
V＝L-A·X(1)
其中：V：残差向量(米)；L：观测值减计算值(米)；A：设计矩阵；X： 参数向量。
构造标准化残差作为粗差检测量：
QVV＝P-1-A(ATPA)-1AT(4)
其中，Vi：残差向量的元素；P：观测值的权阵；δ：单位权中误差；QViVi：残 差协因数矩阵QVV对角线元素。
如果标准化残差则存在粗差，标准化残差时，对权阵做 如下处理：
将公式(1)中L分为不包含粗差的l1和包含粗差向量S的l2，其维数分别为 n1和n2。
则(1)式变为：
设H＝[0E]T，E为n2×n2单位阵，计算粗差向量S并将其改正，
PSS＝HTPQVVPH(8)。
步骤1.3，载波周跳探测与剔除：先通过历元间二次差探测两周以上的大周 跳，然后通过基于历元间单差观测值的单站周跳探测方法探测与剔除小周跳。
具体的，通过历元间二次差探测两周以上的大周跳：
根据当前历元卫星的载波观测值与上一历元的载波观测值计算当前 历元的历元间载波观测值单差：
计算上一历元的历元间载波观测值单差：
计算当前历元的历元间载波观测值二次差：
当载波连续跟踪时，计算所有观测卫星的历元间载波观测值二次差 向量F中各元素集中分布，若某一元素超出数值 中心2周，则认为是大周跳，数值中心通过分布集中的元素求平均值得出。
当载体处于低速运动时(运动状态不固定)，载波观测值的二次差值取决于 运动状态，其数值与高次差法计算出的特征不一致，不服从于高斯白噪声分布， 但每一历元计算的二次差数值都集中分布在某个中心点附近，这是因为所有卫星 的观测值都和当前运动状态相关。根据该分布特点可确定载波是否发生周跳，具 体计算方法为：首先对载波观测值二次差值进行从小到大排序，选出该组数值的 中位数作为数值的中心点，计算所有数值到中心点的距离，超出2周的数值即认 为载波发生了大周跳。
具体的，通过基于历元间单差观测值的单站周跳探测方法探测与剔除小周跳：
构造历元间单差观测方程为：
其中，历元间单差观测值(周)；ρSp：历元t的卫地距(米)；ρSq：历 元t+1的卫地距(米)；△t：接收机钟变化量(秒)；c：光速(m/s)；λ:载 波波长(米)；△N：载波周跳值(周)；△ε：观测噪声(米)。
待估参数包括三个坐标参数和一个接收机参数，若验后单位权中误差在毫米 量级，则不存在周跳，否则采用抗差最小二乘算法方程，逐步进行迭代计算，直 至前后两次解的坐标差值符合限差要求，等价权函数选取如下：
其中，表示标准化残差；P：权阵模型，采用高度角定权。
若标准化残差认为可能存在小周跳，然后在(9)中每次 增加一颗卫星的周跳估值作为待估参数，其余异常卫星的权值置为零，循环计算 每颗卫星周跳值，若周跳取整值不等于零则认为存在小周跳，剔除小周跳，否则， 将异常卫星归为正常卫星。
步骤2，采用伪距观测值进行SPP定位，根据上一历元的RTD伪距残差和 当前历元SPP定位计算的伪距残差进行RTD定位包括以下步骤：
步骤2.1，构造伪距双差观测方程为：
其中，双差算子；Pi：伪距观测值(米)；ρ：卫地距(米)；I：电离层 延迟(米)；双差伪距观测噪声(米)。
根据获取的高质量观测值进行SPP定位解算，得到基准站和流动站SPP定 位计算的伪距残差。
步骤2.2，由上一历元的系统值和方差，通过Kalman滤波模型递推可得到当 前历元的状态估计，根据测量方程确定当前历元的系统值和方差。
其中，Kalman滤波模型为：
其中，I：单位矩阵；一步预测值；Pk,k-1：一步预测值方差-协方差矩阵； Tk-1：系统噪声驱动矩阵；Qk-1：系统噪声序列的方差阵；Kk：增益矩阵Lk：观 测值序列；待估参数；Hk：系数矩阵；Rk：测量噪声的方差阵；Xk：滤 波估值；Pk：滤波方差-协方差矩阵。
步骤2.3，流动站位置初始值根据滤波模型递推结果进行确定，通过Kalman 滤波模型实时更新伪距双差OMC(观测值减测量值)。设基准站和流动站SPP 定位的伪距残差为Vb和Vr，上一历元的RTD伪距残差为V0，当前历元的观测值 方差由以下模型进行计算：
在这一步骤中，首先进行SPP定位，若解算失败则不输出；若解算正确，则 进行RTD定位。根据上一历元的RTD伪距残差V0和当前历元SPP定位计算的伪 距残差Vb和Vr计算得到当前历元的观测值方差有利于提高RTD定位的精度 和稳定性。
步骤3，根据RTD定位结果作为初始值分别通过最小二乘解算、附加模糊度 参数的Kalman滤波计算模糊度浮点解包括以下步骤:
步骤3.1，构造双差观测方程：
其中，,双差算子；Pi：伪距观测值(米)；ρ：卫地距(米)；I：电离层 延迟(米)；双差伪距观测噪声(米)；φi：载波观测值(米)；λi：载波 波长(米)；双差模糊度(周)。
步骤3.2，根据双差观测方程将流动站的位置、速度、加速度以及站间单差 模糊度作为待估参数，通过转换矩阵将站间单差模糊度转换成双差模糊度 Kalman滤波的参数矩阵和转换矩阵分别为：
Xk＝[x y z vx vy vz ax ay az △N1 … △Nn]T (15)
步骤3.3，和伪距相对定位相同，浮点解Kalman滤波观测方程同样更新双差 观测方程OMC。为提高滤波的数值稳定性，控制滤波发散，Pk矩阵改写为：
步骤3.4，根据Kalman滤波解算出的模糊度浮点解可能因为滤波收敛特征会 出现较大误差，通过最小二乘解算出的模糊度浮点解和验后残差σ，由Kalman 滤波解算出的模糊度浮点解做如下变换以提高精度：
步骤4，采用优化的部分模糊度固定策略将模糊度浮点解转换为固定解，包 括以下步骤：
根据浮点解方差对模糊度进行排序，模糊度固定按照浮点解精度顺序逐步进 行搜索；对双差模糊度参数和对应的方差－协方差矩阵进行降相关，得到相关性 较低的变换模糊度，降相关公式为：
其中，ZT：变换矩阵；模糊度浮点解；变换后的模糊度；模糊度 浮点解方差-协方差矩阵；变换后的模糊度浮点解方差-协方差矩阵，
矩阵元素全为整数，行列式的值为±1，变换后的搜索空间为：
其中，χ2：搜索空间；
搜索得到的整数解通过逆变换得到原始模糊度的整数解：
定义Ratio值对模糊度固定解进行检验：
其中，为次优解；为最优解；若满足上式则认为解算正确，若不满足 上式则认为解算失败，保留浮点解，根据浮点解输出定位信息；
优化过程为：当连续若干个模糊度均固定，符合首次模糊度固定条件，保留 模糊度固定解，并设置固定解精度；对部分模糊度进行检验，若满足固定条件， 将部分模糊度设置为模糊度固定。该优化过程通过逐步固定并保留部分模糊度的 方法提高模糊度固定率，进一步提高系统在遮蔽、半遮蔽区域的定位的稳定性和 可靠性。
最后计算得到位置参数固定解，根据模糊度固定解、位置参数固定解输出定 位信息，位置参数固定解计算公式：
其中，位置参数浮点解；位置参数固定解；模糊度浮点解； 模糊度固定解；位置参数和模糊度参数协方差矩阵；模糊度 参数方差－协方差矩阵。
在本实施例中，对原始观测值进行数据预处理和质量控制。根据原始数据信 噪比、卫星高度角、连续跟踪情况等信息选取参与解算的卫星，对伪距粗差进行 探测和修复确保伪距数据质量，对载波观测值进行周跳探测与剔除，为高精度的 相对定位提供数据质量保障.
根据SPP和RTD解算结果对伪距数据质量进行实时评估。通过SPP计算的 伪距残差和上一历元的验后残差，确定RTD中双差伪距观测方程的随机模型， 根据RTD计算的双差伪距残差得出伪距残差，为浮点解的解算提供参考。
将站间单差模糊度作为参数通过Kalman滤波进行实时估计。Kalman滤波通 过对双差测量观测值进行实时更新提高浮点解精度，将单差模糊度作为参数进行 实时估计可以有效抵抗模糊度固定错误、避免参考卫星变化带来的复杂数据处理、 提高解算效率。
采用Kalman滤波和最小二乘原理确定浮点解。由于原始伪距观测值精度较 差，为减少收敛时间、防止滤波发散以及提高浮点解的精度，对Kalman滤波浮 点解进行合理控制，浮点解的结果根据Kalman滤波和最小二乘结果进行确定。
将浮点解转换成固定解。浮点解的精度越高，转换成固定解的成功率越高， 首先确定浮点解精度高的卫星观测值的模糊度，将其设定为固定值，通过Kalman 事后滤波提高剩余模糊度浮点解的精度，然后再逐步确定剩余的模糊度。
本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的 修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权 利要求书所定义的范围。