一种模糊度固定的长基线卫星编队GNSS相对定位方法

摘要

一种模糊度固定的长基线卫星编队GNSS相对定位方法，目的是提高模糊度固定成功率和相对定位结果的精度。技术方案是先收集和预处理输入数据，确定编队卫星绝对概略轨道；然后消除差分观测数据中的几何距离和钟差，估计单差相位模糊度浮点解和单差电离层延迟参数，采用双差变换得到双差宽巷模糊度浮点解和协方差矩阵，固定双差宽巷整周模糊度和双差窄巷整周模糊度，最后输出模糊度固定的相对定位结果。采用本发明避免了传统方法M-W组合中伪码、相位等权处理所导致的模糊度固定强烈依赖于观测精度较低的伪码的问题，提高了长基线卫星编队GNSS相对定位模糊度固定成功率和最终相对定位结果的精度，计算稳定，提高了相对定位结果的可靠性。

主权项

一种模糊度固定的长基线卫星编队GNSS相对定位方法，其特征在于包括以下步骤：第一步，收集和预处理输入数据：1.1从卫星平台收集观测数据：GNSS观测数据、卫星平台姿态数据；从网站收集辅助数据：GNSS卫星精密星历、钟差以及天线PCV信息、地球重力场、地球旋转参数、UTC时间跳秒数据、JPL太阳历、太阳辐射流量、地磁指数；1.2对收集的观测数据进行预处理：对GNSS观测数据的野值进行剔除，并探测相位周跳；第二步，确定编队卫星A、B的绝对轨道位置，输出编队卫星A、B的轨道力学模型参数y_A0、y_B0和每个时刻t的绝对轨道位置r_A(t)、r_B(t)，轨道位置包含x、y、z三个直角坐标分量；第三步，消除差分观测数据中的几何距离和钟差：3.1单差相对定位，方法是：单差是对两个接收机观测的同一颗GNSS卫星、同种类型观测数据做差，形成差分观测数据，采用t时刻的单差观测方程(1)消除GNSS卫星钟差，保留接收机相对钟差，$\begin{array}{c}{\mathrm{\Delta P}}_1^j\left(t\right)=P_{1,A}^j\left(t\right)-P_{1,B}^j\left(t\right)={\rho }_{AB}^j\left(t\right)+c·{\mathrm{\delta t}}_{AB}\left(t\right)+I_{AB}^j\left(t\right)+{ϵ}_{{\mathrm{\Delta P}}_1}^j\left(t\right)\\ {\mathrm{\Delta P}}_2^j\left(t\right)=P_{2,A}^j\left(t\right)-P_{2,B}^j\left(t\right)={\rho }_{AB}^j\left(t\right)+c·{\mathrm{\delta t}}_{AB}\left(t\right)+\frac{f_1^2}{f_2^2}·I_{AB}^j\left(t\right)+{ϵ}_{{\mathrm{\Delta P}}_2}^j\left(t\right)\\ {\mathrm{\Delta L}}_1^j\left(t\right)=L_{1,A}^j\left(t\right)-L_{1,B}^j\left(t\right)={\rho }_{AB}^j\left(t\right)+c·{\mathrm{\delta t}}_{AB}\left(t\right)-I_{AB}^j\left(t\right)+{\lambda }_1·{\mathrm{\Delta N}}_1^{\left(k\right)}+{ϵ}_{{\mathrm{\Delta L}}_1}^j\left(t\right)\\ {\mathrm{\Delta L}}_2^j\left(t\right)=L_{2,A}^j\left(t\right)-L_{2,B}^j\left(t\right)={\rho }_{AB}^j\left(t\right)+c·{\mathrm{\delta t}}_{AB}\left(t\right)-\frac{f_1^2}{f_2^2}·I_{AB}^j\left(t\right)+{\lambda }_2·{\mathrm{\Delta N}}_2^{\left(k\right)}+{ϵ}_{{\mathrm{\Delta L}}_2}^j\left(t\right)\end{array}---\left(1\right)$其中t表示观测时刻，上标j表示第j颗GNSS卫星，1≤j≤M，M是卫星的总数；下标1、2表示GNSS信号频率f₁、f₂；Δ表示单差；为在时刻t频率f₁的伪码差分观测数据；为在时刻t频率f₂的伪码差分观测数据；为在时刻t频率f₁的相位差分观测数据；为在时刻t频率f₂的相位数据差分观测数据；A、B表示两颗编队卫星；为编队卫星A在时刻t频率f₁的伪码观测数据；为编队卫星A在时刻t频率f₂的伪码观测数据；为编队卫星A在时刻t频率f₁的相位观测数据；为编队卫星A在时刻t频率f₂的相位观测数据；为编队卫星B在时刻t频率f₁的伪码观测数据；为编队卫星B在时刻t频率f₂的伪码观测数据；为编队卫星B在时刻t频率f₁的相位观测数据；为编队卫星B在时刻t频率f₂的相位观测数据；均来自第一步中经过预处理后的GNSS观测数据；为编队卫星A、B在时刻t与GNSS卫星j之间的单差几何距离；c表示光速，δt_AB为编队卫星A、B两个GNSS接收机之间的相对钟差；为在时刻t的单差电离层延迟；为第k个连续跟踪弧段f₁频率的单差相位模糊度，1≤k≤n_b，n_b为连续跟踪弧段总数；为第k个连续跟踪弧段f₂频率的单差相位模糊度；ε为观测误差；采用消电离层组合消除电离层延迟得到消电离层单差观测方程(2)，$\begin{array}{c}{\mathrm{\Delta P}}_{IF}^j\left(t\right)=\frac{f_1^2}{f_1^2-f_2^2}{\mathrm{\Delta P}}_1^j\left(t\right)-\frac{f_2^2}{f_1^2-f_2^2}{\mathrm{\Delta P}}_2^j\left(t\right)={\rho }_{AB}^j\left(t\right)+c·{\mathrm{\delta t}}_{AB}\left(t\right)+{ϵ}_{{\mathrm{\Delta P}}_{IF}}^j\left(t\right)\\ {\mathrm{\Delta L}}_{IF}^j\left(t\right)=\frac{f_1^2}{f_1^2-f_2^2}{\mathrm{\Delta L}}_1^j\left(t\right)-\frac{f_2^2}{f_1^2-f_2^2}{\mathrm{\Delta L}}_2^j\left(t\right)={\rho }_{AB}^j\left(t\right)+c·{\mathrm{\delta t}}_{AB}\left(t\right)+{\mathrm{\Delta A}}_{IF}^j+{ϵ}_{{\mathrm{\Delta L}}_{IF}}^j\left(t\right)\end{array}---\left(2\right)$其中“IF”表示消电离层组合；为时刻t的单差伪码消电离层组合数据；为时刻t的单差相位消电离层组合数据；为第k个连续跟踪弧段对应的单差相位消电离层模糊度；单差几何距离通过公式(3)计算，${\rho }_{AB}^j\left(t\right)=|r_B\left(t\right)-r_G^j\left(t\right)|-|r_A\left(t\right)-r_G^j\left(t\right)|---\left(3\right)$其中| |表示计算三维矢量的长度；r_A(t)、r_B(t)为编队卫星A、B在时刻t的绝对轨道位置；为GNSS卫星j在时刻t的位置，利用IGS提供的事后精密星历插值计获得；将在概略点处线性化展开，得到公式(4),${\rho }_{AB}^j\left(t\right)\approx {\rho }_{AB,0}^j\left(t\right)-e_B^j\left(t\right)·{\mathrm{\Delta r}}_B\left(t\right)---\left(4\right)$表示在时刻t从卫星B到GNSS卫星j的视线方向单位矢量；为的概略点，其初始值采用第二步得到的编队卫星A、B的绝对轨道位置计算得到；Δr_B(t)为卫星B在t时刻的轨道位置改进量，Δr_B(t)＝r_AC(y_B0+Δy_B,t)‑r_AC(y_B0,t)          (5)其中r_AC表示轨道力学积分函数，r_AC(y_B0,t)表示对轨道力学模型参数y_B0进行Adams‑Cowell多步积分得到的卫星B在时刻t的轨道位置；y_B0、Δy_B分别表示卫星B的轨道力学模型参数向量的初始值和改进量，包括初始时刻轨道位置和速度、太阳光压系数、大气阻力系数和经验加速度系数；y_B0采用第二步输出的轨道力学模型参数；Δy_B根据公式(6)采用最小二乘估计得到；将多个时刻的消电离层单差观测方程(2)整理成矩阵形式，得到单差观测改正值向量z_SD，z_SD＝H_SD·x_SD+e_SD               (6)其中“SD”表示单差；单差观测改正值向量$z_{SD}={[...,{\mathrm{\Delta P}}_{IF}^j\left(t\right)-{\rho }_{AB,0}^j\left(t\right),{\mathrm{\Delta L}}_{IF}^j\left(t\right)-{\rho }_{AB,0}^j\left(t\right),...]}^T,$[ ]^T表示矩阵转置；设计矩阵表示偏导数；误差向量$e_{SD}={[...,{ϵ}_{{\mathrm{\Delta P}}_{IF}}^j\left(t\right),{ϵ}_{{\mathrm{\Delta L}}_{IF}}^j\left(t\right),...]}^T;$待估参数向量$x_{SD}={[...,{\mathrm{\delta t}}_{AB}\left(t\right),...,{\mathrm{\Delta A}}_{IF}^{\left(k\right)},...,{\mathrm{\Delta y}}_B]}^T;$此时待估参数包括：接收机相对钟差δt_AB(t)、单差相位消电离层模糊度和轨道力学模型参数改进向量Δy_B；3.2采用加权最小二乘方法估计得到x_SD，$x_{SD}={[H_{SD}^T·W_{SD}·H_{SD}]}^{-1}·H_{SD}^T·W_{SD}·z_{SD}---\left(7\right)$其中单差观测权矩阵$W_{SD}=diag(...,1/2{\sigma }_P^2,1/2{\sigma }_L^2,...),$σ_L表示相位观测精度，σ_P表示伪码观测精度；3.3计算单差几何距离：卫星轨道位置r_B(t)通过Adams‑Cowell多步积分r_AC(y_B0+Δy_B,t)获得，再将r_B(t)代入公式(3)更新计算单差几何距离3.4消除差分观测数据中单差几何距离和钟差项δt_AB(t)：将3.3得到的单差几何距离和3.2得到的钟差项δt_AB(t)从公式(1)中的观测数据中扣除，得到多个时刻扣除几何距离和钟差后的单差观测方程(8)，$\begin{array}{c}{\mathrm{\Delta P}}_1^j\left(t\right)-{\rho }_{AB}^j\left(t\right)-c·{\mathrm{\delta t}}_{AB}\left(t\right)=I_{AB}^j\left(t\right)+{ϵ}_{{\mathrm{\Delta P}}_1}^j\left(t\right)\\ {\mathrm{\Delta P}}_2^j\left(t\right)-{\rho }_{AB}^j\left(t\right)-c·{\mathrm{\delta t}}_{AB}\left(t\right)=\frac{f_1^2}{f_2^2}·I_{AB}^j\left(t\right)+{ϵ}_{{\mathrm{\Delta P}}_2}^j\left(t\right)\\ {\mathrm{\Delta L}}_1^j\left(t\right)-{\rho }_{AB}^j\left(t\right)-c·{\mathrm{\delta t}}_{AB}\left(t\right)=-I_{AB}^j\left(t\right)+{\lambda }_1·{\mathrm{\Delta N}}_1^{\left(k\right)}+{ϵ}_{{\mathrm{\Delta L}}_1}^j\left(t\right)\\ {\mathrm{\Delta L}}_2^j\left(t\right)-{\rho }_{AB}^j\left(t\right)-c·{\mathrm{\delta t}}_{AB}\left(t\right)=-\frac{f_1^2}{f_2^2}·I_{AB}^j\left(t\right)+{\lambda }_2·{\mathrm{\Delta N}}_2^{\left(k\right)}+{ϵ}_{{\mathrm{\Delta L}}_2}^j\left(t\right)\end{array}---\left(8\right)$第四步，估计单差相位模糊度浮点解和单差电离层延迟参数：4.1设置伪码权值w_P与相位权值w_L之比，取w_P:w_L＝σ_L:σ_P，σ_L表示相位观测精度，σ_P表示伪码观测精度；4.2采用加权最小二乘方法估计单差相位模糊度浮点解和单差电离层延迟参数；在第k个连续跟踪弧段内，将多个时刻扣除几何距离和钟差的单差观测方程(8)整理成矩阵形式，z_k＝H_k·x_k+e_k          (9)其中单差观测值改正向量$z_k=\left[\begin{array}{c}·\\ ·\\ ·\\ \Delta P_1^j\left(t\right)-{\rho }_{AB}^j\left(t\right)-c·\delta t_{AB}\left(t\right)\\ \Delta P_2^j\left(t\right)-{\rho }_{AB}^j\left(t\right)-c·\delta t_{AB}\left(t\right)\\ \Delta L_1^j\left(t\right)-{\rho }_{AB}^j\left(t\right)-c·\delta t_{AB}\left(t\right)\\ \Delta L_2^j\left(t\right)-{\rho }_{AB}^j\left(t\right)-c·\delta t_{AB}\left(t\right)\\ ·\\ ·\\ ·\end{array}\right];$误差向量$e_k=\left[\begin{array}{c}·\\ ·\\ ·\\ {ϵ}_{{\mathrm{\Delta P}}_1}^j\left(t\right)\\ {ϵ}_{{\mathrm{\Delta P}}_2}^j\left(t\right)\\ {ϵ}_{{\mathrm{\Delta L}}_1}^j\left(t\right)\\ {ϵ}_{{\mathrm{\Delta L}}_2}^j\left(t\right)\\ ·\\ ·\\ ·\end{array}\right];$设计矩阵待估参数向量此时待估参数包括：单差电离层延迟，2个单差相位模糊度参数和得到x_k的加权最小二乘估计值为，${\hat{x}}_k={[H_k^T·W_k·H_k]}^{-1}·H_k^T·W_k·z_k---\left(10\right)$其中观测权矩阵$W_k=diag(...,w_P^2/w_L^2,w_P^2/w_L^2,1,1,...),$diag( )表示对角矩阵；输出估计得到第k个连续跟踪弧段的单差相位模糊度参数和单差电离层延迟参数第五步，采用双差变换得到双差宽巷模糊度浮点解和协方差矩阵：5.1将第四步得到的第k个连续跟踪弧段的单差相位模糊度参数相减，得到对应的单差宽巷模糊度则n_b个连续跟踪弧段对应的单差宽巷相位模糊度分别为定义浮点解向量及其协方差矩阵构造单差宽巷模糊度浮点解向量方差为此时$\mathrm{cov}\left({\hat{z}}_{W,SD}\right)=diag(...,{\sigma }_{W,SD}^2,...);$5.2根据双差和单差观测方程之间的关系，构造双差线性变换算子T_DD，得到双差宽巷模糊度浮点解为，${\hat{z}}_{W,DD}=T_{DD}·{\hat{z}}_{W,SD}---\left(11\right)$其中${\hat{z}}_{W,DD}={[...,▿\Delta {\hat{N}}_W^{\left(jk\right)},▿\Delta {\hat{N}}_W^{\left(jm\right)},...]}^T,$${\hat{z}}_{W,SD}={[...,\Delta {\hat{N}}_W^{\left(j\right)},\Delta {\hat{N}}_W^{\left(k\right)},\Delta {\hat{N}}_W^{\left(m\right)},...]}^T,$通过公式(12)计算得到双差宽巷模糊度协方差矩阵$\mathrm{cov}\left({\hat{z}}_{W,DD}\right)=T_{DD}·\mathrm{cov}\left({\hat{z}}_{W,SD}\right)·T_{DD}^T---\left(12\right)$第六步，固定双差宽巷整周模糊度：6.1采用整数最小二乘方法固定双差宽巷整周模糊度：将第五步得到的双差宽巷模糊度浮点解及其协方差矩阵作为输入，采用整数最小二乘方法固定双差宽巷整周模糊度，即搜索整数空间得到双差宽巷整周模糊度固定解使得公式(13)达到最小，$\underset{z_{W,DD}\in Z^n}{\min }\{{[z_{W,DD}-{\hat{z}}_{W,DD}]}^T·{\mathrm{cov}}^{-1}\left({\hat{z}}_{W,DD}\right)·[z_{W,DD}-{\hat{z}}_{W,DD}\left]\right\}---\left(13\right)$其中cov^‑1表示协方差矩阵的逆；min{ }表示求最小值；Zⁿ表示整数空间；采用最小二乘去相关算法LAMBDA进行搜索得到双差宽巷相位模糊度的整数解；6.2检验双差宽巷整周模糊度固定正确性：计算LAMBDA方法输出的次优整周模糊度对应的模糊度残差量R_S和最优整周模糊度对应的模糊度残差量R_B之比，若R_S/R_B≤k_S/B，则双差宽巷整周模糊度固定解z_W,DD被采用，转第七步；若R_S/R_B>k_S/B，则双差宽巷整周模糊度固定解z_W,DD固定不正确，双差宽巷模糊度仍采用其浮点解即令k_S/B为给定的门限值，转第七步；第七步，固定双差窄巷整周模糊度：7.1将第六步得到的双差宽巷整周模糊度固定解z_W,DD作为已知量从双差相位观测方程中扣除，双差观测值包括双差伪码和双差相位，是在单差观测方程的基础上再对不同GNSS卫星的单差观测数据作差，消除接收机相对钟差，即$▿{\mathrm{\Delta P}}_{IF}^{jk}\left(t\right)={\rho }_{AB}^{jk}\left(t\right)+{ϵ}_{▿{\mathrm{\Delta P}}_{IF}}^{jk}\left(t\right)$$▿{\mathrm{\Delta L}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}·▿{\mathrm{\Delta N}}_W^{\left(jk\right)}={\rho }_{AB}^{jk}\left(t\right)+\frac{c}{2(f_1+f_2)}·▿{\mathrm{\Delta N}}_N^{\left(jk\right)}+{ϵ}_{▿{\mathrm{\Delta L}}_{IF}}^{jk}\left(t\right)---\left(14\right)$其中表示双差；上标j、k表示第j、k颗GNSS卫星；表示GNSS卫星j、k之间的双差伪码消电离层观测值，表示GNSS卫星j、k之间的双差相位消电离层观测值，$▿{\mathrm{\Delta L}}_{IF}^{jk}\left(t\right)={\mathrm{\Delta L}}_{IF}^j\left(t\right)-{\mathrm{\Delta L}}_{IF}^k\left(t\right);$分别表示双差宽巷模糊度和双差窄巷模糊度；表示GNSS卫星j、k之间的双差几何距离，${\rho }_{AB}^{jk}\left(t\right)={\rho }_{AB}^j\left(t\right)-{\rho }_{AB}^k\left(t\right);$将多个时刻的双差相位观测方程(14)整理成矩阵形式，得到多个时刻双差观测方程(15)z_DD＝H_DD·x_DD+e_DD      (15)其中“DD”表示双差；z_DD为双差观测值改正向量，$z_{DD}={[...,▿{\mathrm{\Delta P}}_{IF}^{jk}\left(t\right)-P_{AB,0}^{jk}\left(t\right),▿{\mathrm{\Delta L}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}·▿N_W^{\left(jk\right)}-{\rho }_{AB,0}^{jk}\left(t\right),···]}^T;$设计矩阵误差向量$e_{DD}={[...,{ϵ}_{{\mathrm{\Delta P}}_{IF}}^{jk}\left(t\right),{ϵ}_{{\mathrm{\Delta L}}_{IF}}^{jk}\left(t\right),...]}^T;$待估参数向量$x_{DD}={[...,▿{\mathrm{\Delta N}}_N^{\left(jk\right)},...,{\mathrm{\Delta y}}_B]}^T;$7.2根据多个时刻双差观测方程(15)，采用最小二乘估计方法，估计得到窄巷相位模糊度浮点解向量及其协方差矩阵7.3采用整数最小二乘方法固定双差窄巷整周模糊度；将窄巷相位模糊度浮点解及其协方差矩阵作为输入，采用整数最小二乘方法固定双差窄巷整周模糊度，搜索整数空间得到双差窄巷整周模糊度固定解使得公式(16)达到最小，$\underset{z_{N,DD}\in Z^n}{\min }\{{[z_{N,DD}-{\hat{z}}_{N,DD}]}^T·{\mathrm{cov}}^{-1}\left({\hat{z}}_{N,DD}\right)·[z_{N,DD}-{\hat{z}}_{N,DD}\left]\right\}---\left(16\right)$7.4检验双差窄巷整周模糊度固定正确性：计算次优整周模糊度对应的模糊度残差量R_S和最优整周模糊度对应的模糊度残差量R_B之比，若R_S/R_B≤k_S/B，则双差窄巷整周模糊度固定解z_N,DD被采用，转第八步；若R_S/R_B>k_S/B，则双差窄巷整周模糊度固定解z_N,DD固定不正确，双差窄巷模糊度仍采用其浮点解即令转第八步；第八步，输出模糊度固定的相对定位结果：8.1将第七步得到的双差窄巷整周模糊度固定解作为已知量再一次从双差相位观测方程(14)中扣除，得到扣除模糊度后的双差观测方程(17)：$▿{\mathrm{\Delta P}}_{IF}^{jk}\left(t\right)={\rho }_{AB}^{jk}\left(t\right)+{ϵ}_{▿{\mathrm{\Delta P}}_{IF}}^{jk}\left(t\right)$$▿{\mathrm{\Delta L}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}·▿{\mathrm{\Delta N}}_W^{\left(jk\right)}-\frac{c}{2(f_1+f_2)}·▿{\mathrm{\Delta N}}_N^{\left(jk\right)}={\rho }_{AB}^{jk}\left(t\right)+{ϵ}_{▿{\mathrm{\Delta L}}_{IF}}^{jk}\left(t\right)---\left(17\right)$将多个时刻的扣除模糊度后的双差观测方程(17)整理成矩阵形式，得到多个时刻扣除模糊度后的双差观测方程(18)$z_{DD}^{*}=H_{DD}^{*}·{\mathrm{\Delta y}}_B+e_{DD}---\left(18\right)$其中双差观测值改正向量，$z_{DD}^{*}={[...,▿{\mathrm{\Delta P}}_{IF}^{jk}\left(t\right)-{\rho }_{AB,0}^{jk}\left(t\right),▿{\mathrm{\Delta L}}_{IF}^{jk}\left(t\right)-\frac{c}{2(f_1-f_2)}·▿{\mathrm{\Delta N}}_W^{\left(jk\right)}-\frac{c}{2(f_1-f_2)}·▿{\mathrm{\Delta N}}_N^{\left(jk\right)}-{\rho }_{AB,0}^{jk}\left(t\right),...]}^T;$设计矩阵待估参数向量Δy_B为轨道力学模型参数改进向量；8.2根据多个时刻扣除模糊度后的双差观测方程(18)，采用最小二乘估计方法，估计得到轨道力学模型参数改进向量Δy_B，卫星轨道位置r_B(t)通过轨道力学积分r_AC(y_B0+Δy_B,t)获得，输出最终模糊度固定后的编队相对轨道r_B(t)‑r_A(t)。